表示系统在控制过程结束后，要求系统的终端状态x (tf )应达到某些要求，终端时刻tf 可以固定，也可以自由，视最优控制问题的性质而定.

Presentation on theme: "表示系统在控制过程结束后，要求系统的终端状态x (tf )应达到某些要求，终端时刻tf 可以固定，也可以自由，视最优控制问题的性质而定."— Presentation transcript:

1 表示系统在控制过程结束后，要求系统的终端状态x (tf )应达到某些要求，终端时刻tf 可以固定，也可以自由，视最优控制问题的性质而定.
2018年9月17日星期一 (2)末值型性能指标 1.3 最优控制问题的数学分类法 (1) 拉格朗日问题 已知 表示系统在控制过程结束后，要求系统的终端状态x (tf )应达到某些要求，终端时刻tf 可以固定，也可以自由，视最优控制问题的性质而定. 求最优控制u (t ), 使得以下性能指标 (3) 复合型性能指标 为极小。 (2) 迈耶尔问题 已知 表示对控制过程及控制过程结束后的终端状态均有要求，是最一般形式的性能指标。 但是终端状态的某些分量x (tf )与终端时刻tf 没有给定，求最优控制u (t) 使得以下性能指标为极小。

2 求最优控制u (t )使得以下性能指标为极小。
2018年9月17日星期一 (3) 波尔扎问题 已知 2. 连续系统 LQ 最优控制的MATLAB 实现 给定连续LTI系统状态方程 求最优控制u (t )使得以下性能指标为极小。 其中 A R nn , B  R nr , x  R n, u R r 。 系统的性能指标为： 1.4 线性二次型最优控制 系统：状态空间表示的线性系统； 指标：状态与控制输入的二次型函 数； 问题：在约束条件下寻求控制输入 使得二次型目标函数达到极 小。 其中:Q R nn , R  R rr , R=RT > 0. 假定条件 (1) Q > 0, Q = QT , (A, B)完全可控； 或 (2) Q  0, Q = QT , (A, B)完全可控; (A , C )完全可观测，DDT = Q , D 为任意矩阵。

3 的唯一解，同时还反回闭环特征值L以及增义矩阵G = BTX。 存在最优状态反馈矩阵： K = R–1BT X 与唯一的最优控制律：
2018年9月17日星期一 MATLAB命令care 功能：求解连续时间系统的代数 Riccati方程. 格式1: [X,L,G]=care(A,B,Q) 计算Riccati方程 ATX + XA –XBBTX + Q = 0 的唯一解，同时还反回闭环特征值L以及增义矩阵G = BTX。 存在最优状态反馈矩阵： K = R–1BT X 与唯一的最优控制律： u*(t )=–Kx (t ) = –R–1BT Xx (t ) 以及最优性能指标： 其中X 是下列Riccati方程的唯一解 XA + ATX –XB R–1BTX + Q = 0 格式2: [X,L,G]=care(A,B,Q,R,S,E) 求解更一般的Riccati方程 ATXE + ETXA –(ETXB +S)R–1(BTXE + ST ) + Q = 0 当R, S, E 缺省时， R ＝I, S = 0, E = I , 系统在返回解X 的同时还返回增义矩阵G =R–1(BTXE +ST ), 以及闭环特征值L。 闭环系统： 是渐近稳定的，其解为最优轨线 x*(t )。

4 例 A 和 A–B R–1BTX 的特征值为： 需要求解求Riccati方程 ATX + XA –XB R–1BTX + CTC = 0
2018年9月17日星期一 例 l = g = A 和 A–B R–1BTX 的特征值为： 需要求解求Riccati方程 ATX + XA –XB R–1BTX + CTC = 0 >> [eig(a) eig(a-b*g)] ans = >> a=[-3 2;1 1]; >> b=[0;1]; >> c=[1 -1]; >> r=3; >> [x,l,g]=care(a,b,c'*c,r) x = 从而向量 l 给出了闭环系统的特征值。

5 说明：[K,S,e]=lqr(A,B,Q,R,N)计算连续时间系统的最优反馈增义矩阵K，使系统 MATLAB命令lqr2
2018年9月17日星期一 MATLAB命令lqr 功能：连续系统的LQ调节器设计 格式: [K,S,e]=lqr(A,B,Q,R) [K,S,e]=lqr(A,B,Q,R,N) 说明：[K,S,e]=lqr(A,B,Q,R,N)计算连续时间系统的最优反馈增义矩阵K，使系统 MATLAB命令lqr2 功能：连续系统的LQ调节器设计 格式: [K,S,e]=lqr2(A,B,Q,R,N) 说明：[K,S,e]=lqr2(A,B,Q,R,N)计算连续时间系统的最优反馈增义矩阵K，使系统 采用反馈控制律 u = –Kx, 使性能指标函数 采用反馈控制律 u = –Kx, 使性能指标函数 最小。同时返回代数Riccati方程 最小。同时返回代数Riccati方程 ATS + SA –(SB +N)R–1(BTS + NT ) + Q = 0 的解S 及闭环特征值e。N缺省时为0。这里K = R–1(BTS+NT ). ATS + SA –(SB +N)R–1(BTS + NT ) + Q = 0 的解S 及闭环特征值e。N缺省时为0。这里K = R–1(BTS+NT ).

6 ATS A–S –(ATSB +N)(BTS B+R)–1 (BTSA + NT ) + Q = 0
2018年9月17日星期一 MATLAB命令dlqr 功能：离散系统的LQ调节器设计 格式: [K,S,e]=dlqr(A,B,Q,R) [K,S,e]=dlqr(A,B,Q,R,N) 说明：[K,S,e]=dlqr(A,B,Q,R,N)计算离散时间系统的最优反馈增义矩阵K，使系统 MATLAB命令lqry 功能：连续系统的LQ调节器设计 格式: [K,S,e]=lqry(sys,Q,R) [K,S,e]=lqry(sys,Q,R,N) [K,S,e]=lqry(A,B,Q,R,N) 说明：[K,S,e]=lqry(A,B,Q,R,N)计算连续时间系统的最优反馈增义矩阵K，使系统 采用反馈控制律 u [n]= –Kx[n], 使性能指标函数 采用反馈控制律 u = –Kx, 使性能指标函数 最小。同时返回代数Riccati方程 ATS A–S –(ATSB +N)(BTS B+R)–1 (BTSA + NT ) + Q = 0 的解S 及闭环特征值e。N缺省时为0. 这里K = (BTSB+R)–1(BTSA+NT ). 最小。同时返回代数Riccati方程的解S 及闭环特征值e。 N缺省时为0.

7 试计算最优状态反馈矩阵K , 使得J 最小并对其闭环系统进行单位阶跃给定响应仿真。 由结构图有系统的控制信号：
2018年9月17日星期一 例 已知系统动态系统方程为 式中反馈增义矩阵K ： K = [k1 k2 k3] 系统性能指标为： 系统结构图如图所示 试计算最优状态反馈矩阵K , 使得J 最小并对其闭环系统进行单位阶跃给定响应仿真。 由结构图有系统的控制信号： u = – k1(x1 – r) – (k2x2+k3x3) = k1r – (k1x1+ k2x2+k3x3) = k1r – Kx

8 首先求出最优反馈矩阵K , 由结构图可得系统的状态方程为：
2018年9月17日星期一 首先求出最优反馈矩阵K , 由结构图可得系统的状态方程为： >> l8201 K = S = e = i i

9 试计算最优状态反馈矩阵K , 使得J 最小并对其闭环系统进行单位阶跃给定响应仿真。
2018年9月17日星期一 例 续上题。要求采用输出作为性能指标， >> l8202 K = 性能指标为 试计算最优状态反馈矩阵K , 使得J 最小并对其闭环系统进行单位阶跃给定响应仿真。

10 利用linmode2(‘t8204’) 可以得到该系统的状态方程具有上面例子的形式。
2018年9月17日星期一 例. 已知可控直流电源供电给直流电机的系统结构图如图所示。欲对系统进行最优状态反馈与输出反馈控制，试分别计算状态指标及输出指标下反馈增义矩阵，并对其闭环系统进行阶跃响应仿真。 利用状态指标时： 利用输出指标时： Qy =1000; R = 1 利用linmode2(‘t8204’) 可以得到该系统的状态方程具有上面例子的形式。

11 运行MATLAB程序求出Kx 和Ky 如下 超调量稍大，峰值时间较短，有一次振荡 超调量较小，超调后衰减。 >> l8203
2018年9月17日星期一 运行MATLAB程序求出Kx 和Ky 如下 >> l8203 Kx = Ky = 1.0e+003 * 超调量稍大，峰值时间较短，有一次振荡 超调量较小，超调后衰减。

12 3. 离散系统稳态线性二次型最优控制的MATLAB实现 则有最优反馈矩阵为:
2018年9月17日星期一 3. 离散系统稳态线性二次型最优控制的MATLAB实现 则有最优反馈矩阵为: 给定离散线性系统 x (k +1) = Ax (k)+Bu (k), x (0) = x0 (k = 0, 1, … , N – 1) 式中: x (k)R n, u (k) R r, A R nn, B R nr. 与之对应的最优控制序列: 以及最优性能指标： 系统性能指标为: P (k) > 0 是下列Riccati差分方程的解 令N , 系统最优控制解成为稳态解，性能指标改为： Q R nn 为正定或半正定, R R rr 为正定, S R nn 为正定或半正定。

13 ATPA –P – ATPB (BTP B + R)–1BTPA + Q = 0 闭环系统的状态方程为：
2018年9月17日星期一 K (k)变成常数增义矩阵K : K = [R + BTPB]–1BTPA P (k)也变成常数矩阵P : P = Q + AT [P –1+BR–1BT]–1A 易证明： [P –1+BR–1BT ]–1 = P –PB(BTPB +R)–1BTP 从而得代数Riccati方程： ATPA –P – ATPB (BTP B + R)–1BTPA + Q = 0 闭环系统的状态方程为： 最优性能指标仍为： 最优控制序列为:

14 v(k ) = r (k) –y(k) + v(k – 1) y(k ) = cx(k ) 由v(k ) 有：
2018年9月17日星期一 例. 已知伺服系统动态结构如图 由结构图有： x(k + 1) = ax(k) + bu(k) u(k) = k1v(k) – k2x(k) v(k ) = r (k) –y(k) + v(k – 1) y(k ) = cx(k ) 由v(k ) 有： v(k ＋1) = r (k ＋1) –y(k ＋1) + v(k ) = －cax(k)+v(k ) – cbu(k ) + r(k + 1) 且： x(k +1) = ax (k )+bu (k ) y(k ) = c x(k ) a = 0.5, b = 1, c = 1, d = 0 试计算稳态最优反馈增益矩阵并求系统闭环后的单位阶跃给定响应。

15 由x (k +1)与v (k +1)写出矩阵： 由u (k) = k1v (k) –k2x(k) 有
2018年9月17日星期一 由x (k +1)与v (k +1)写出矩阵： 由u (k) = k1v (k) –k2x(k) 有 ue (k) = k1ve (k) –k2xe(k) 再令 x1(k) ＝xe (k), x2(k) = ve(k) , w(k) = ue (k) , 以上矩阵可写成： 对于系统稳态：k = , 有： 令： xe (k) = x (k) – x (), ve (k) = v (k) – v (), ue (k) = u (k) – u (), r (k ) = r () 可得： 当系统有单位阶跃输入给定时， r (k ＋1) = r () = r 并考虑到 u(k) = k1v(k) – k2x(k)

16 2018年9月17日星期一 原问题可整理为： 设定性指标为： 式中Q、R 选择为： 程序运行后的结果如下：

17 例. 考虑如图所示倒立摆系统， 是摆杆与垂直线夹角，假定 很小从而sin  , cos 1, 同时假定 也很小，从而 。已知
2018年9月17日星期一 >> l8301 Kx = 例. 考虑如图所示倒立摆系统， 是摆杆与垂直线夹角，假定 很小从而sin  , cos 1, 同时假定 也很小，从而 。已知 M = 2kg, m = 0.1kg, l = 0.5m

18 2018年9月17日星期一 应用牛顿定律得 考虑m 绕P点的旋转运动 或 该式简化为 进一步有

19 2018年9月17日星期一 利用 得 或写成向量矩阵形式 经过简单的变量代换后： 令 得状态方程： 带入数据后：

20 利用MATLAB命令c2d转化为离散模型 从而得离散模型 x(k +1) = Gx(k) ＋Hu(k)
2018年9月17日星期一 利用MATLAB命令c2d转化为离散模型 从而得离散模型 x(k +1) = Gx(k) ＋Hu(k) y (k) = Cx(k) + Du(k) >> A=[ ]; >> B=[0;-1;0;0.5]; >> [G,H]=c2d(A,B,0.1) G = H = 0.0025 0.0501

21 ue(k ) = –Kxe(k) + K1ve(k)
2018年9月17日星期一 现在考虑伺服系统的二次最优控制 现定义 利用前面结论，令 xe (k ) = x(k ) – x() ve (k ) = v(k ) – v() ue(k ) = –Kxe(k) + K1ve(k)

22 2018年9月17日星期一 我们得到离散系统 针对下列二次指标 选择 MATLAB程序如下

23 2018年9月17日星期一 >> l8301a Kx =

24 首先调用函数linmod2() 将结构图模型转换为状态空间模型：
2018年9月17日星期一 例. 如图所示的可控整流装置供电给直流电机传动系统。试对T=0.1s设计数字最优控制系统(计算稳态最优反馈增益矩阵) 并求闭环离散系统的阶跃响应。 首先调用函数linmod2() 将结构图模型转换为状态空间模型： 调用函数c2d( )对系统按T=0.1s进行采样离散化，得模型： 并得数字控制系统：

25 二次性能指标为： >> l8302 Kx = 0.1036 0.1452 6.5705 -0.0050
2018年9月17日星期一 二次性能指标为： >> l8302 Kx =

26 Kalman滤波器就是带有噪声情况下的最优观测器
2018年9月17日星期一 即有： 4. 最优观测器的MATLAB实现 Kalman滤波器就是带有噪声情况下的最优观测器 4.1 LTI系统的Kalman滤波 给定连续时间系统方程为： 式中x Rn, u  Rr, y Rm, A Rnn, BRnr, G Rnr, C Rmn, D Rmr, w(t )为随机噪声干扰输入，是零均值的r 维白噪声过程；v(t )为随机量测噪声，是零均值的m 维白噪声过程。 w(t )和v(t )两噪声过程平稳且互不相关。 令状态向量与其估计值的之差为： 则在{C, A}完全可观测的假定下，使估计误差平方和的期望值最小，即有

27 [Kest,L,P]=kalman(sys,Q,R,N)
2018年9月17日星期一 MATLAB命令kalman 功能：系统的Kalman滤波器设计 格式: [Kest,L,P]=kalman(sys,Q,R,N) [Kest,L,P,M,Z]=kalman(sys,Q,R,N) %仅用于离散系统(参考帮助文件) [Kest,L,P]=kalman(sys,Q,R,N, sensors,known) 说明：kalman函数用于带有系统噪声和测量噪声的连续时间或离散时间系统设计一个Kalman最优状态估计器。 其最优估计器为： 式中 其中P 0是下列Riccati方程的解： 可以证明： 4.2 Kalman滤波的MATLAB实现 sys为状态空间模型 , 当模型有两个输入时

28 N 为可选项，对应模型噪声与量测噪声相关项； Kest为Kalman滤波器的状态估计器,其方程如下： MATLAB命令lqe
2018年9月17日星期一 Q 为模型噪声的协方差矩阵； R 为量测噪声的协方差矩阵； N 为可选项，对应模型噪声与量测噪声相关项； Kest为Kalman滤波器的状态估计器,其方程如下： MATLAB命令lqe 功能：连续系统Kalman估计器设计 格式: [L,P,E]=lqe(A,G,C,Q,R,N) 说明：lqe返回系统 的观测器增义矩阵L； Q为模型噪声的协方差矩阵; R 为量测噪声的协方差矩阵; N 为可选项，对应模型噪声与量测噪声相关项； P 为对应Riccati方程的解； E 为估计器的闭环特征值。 L为滤波器的增义矩阵；P 为对应Riccati方程的解； sensors用于指定可测输出； known 用于指定已知输入，其余输入为噪声输入(假定已知输入和噪声输入混在一起，系统输出非完全可测)。

29 2018年9月17日星期一 例 已知系统的状态方程与量测方程分别为： 其中，w (t), v(t )都是零均值白噪声,且w (t ) =1e-1, v(t ) = 1e-3。试计算Kalman滤波器的增义矩阵与估计协的方差。 >> l8401 L = 4.2544 9.0499 P =

30 2018年9月17日星期一 例. 已知可控直流电源供电给直流电机的系统结构如图所示，试对系统进行Kalman滤波器的设计并对Kalman滤波后的系统闭环进行阶跃响应仿真.

31 根据题意给出如下MATLAB程序 >> l8402 L = 7.3619 31.9083 0.1031 P =
2018年9月17日星期一 >> l8402 L = 7.3619 0.1031 P = 根据题意给出如下MATLAB程序

32 5. 线性二次型Guass最优控制的MATLAB实现
2018年9月17日星期一 5. 线性二次型Guass最优控制的MATLAB实现 考虑系统随机输入噪声与随机量测噪声的线性二次型最优控制叫做线性(L)二次型(Q)高斯(Gauss)最优控制，也称为LQG控制。 5.1 LQG最优控制求解 给定LTI系统方程： 假定随机干扰噪声w (t )是零均值r 维白噪声过程,随机量测噪声v (t )是零均值m 维白噪声过程。两过程均为平稳的且互不相关。

33 其中P > 0为下列Riccati方程的唯一解： PA + ATP –PBR–1BTP + Q = 0
2018年9月17日星期一 系统的性能指标为： LQG最优控制是两个方面的问题: 一是二次型调节器问题，二是最优估计器问题。 (1)最优调节器： K = R –1BT P 其中P > 0为下列Riccati方程的唯一解： PA + ATP –PBR–1BTP + Q = 0

34 AP0 + P0AT +GQ0GT–P0CTR0–1CP0 = 0
2018年9月17日星期一 (2)Kalman滤波器 L = P0CTR0–1 其中P0为下列Riccati方程的解： AP0 + P0AT +GQ0GT–P0CTR0–1CP0 = 0 其中：Q0, R0 分别为模型噪声及量测噪声协方差阵。 将两问题的解综合在一起便得LQG最优控制问题的解。

35 功能：根据Kalman估计器增益和状态反馈增益建立 LQG调节器。 格式：rlqg=lqgreg(kest,k)
2018年9月17日星期一 5.2 LQG最优控制的MATLAB实现 MATLAB命令 lqgreg 功能：根据Kalman估计器增益和状态反馈增益建立 LQG调节器。 格式：rlqg=lqgreg(kest,k) rlqg=lqgreg(kest,k,controls) 其中：kest 为Kalman滤波器的状态估计器, 可由命 令[Kest,L,P]=kalman(sys,Q,R,N)建立； k 为最优状态反馈增义，可由lqr 或lqry建立。 rlqg 为连接kalman滤波器与最优状态反馈矩 阵后形成的LQG调节器。 controls 为LQG的额外控制通道编号。

36 说明：lqgreg函数根据kalman滤波器和最优状态反馈 增益建立LQG调节器，调节器方程为
2018年9月17日星期一 说明：lqgreg函数根据kalman滤波器和最优状态反馈 增益建立LQG调节器，调节器方程为

37 功能：根据估计器增益L和状态反馈增益K建立来调 节器。 格式：rsys=reg(sys,K,L)
2018年9月17日星期一 MATLAB命令 reg 功能：根据估计器增益L和状态反馈增益K建立来调 节器。 格式：rsys=reg(sys,K,L) rsys=reg(sys,K,L,sensors,known,controls) 其中：sys： 为受控对象状态空间模型; K: 状态反馈增益矩阵; L: 状态估计器增益矩阵； rsys: 返回的调节器； sensor: 传感器通道号向量； known: 已知外部信号通道号向量； control: 状态反馈增益通道号向量.

38 说明: 增益矩阵K、L通常可利用极点配置或二次型最 优以及LQG方法求得。 调节器模型为
2018年9月17日星期一 说明: 增益矩阵K、L通常可利用极点配置或二次型最 优以及LQG方法求得。 调节器模型为

39 例. 已知可控直流电源供电给直流电机的系统结构图如图所示。对系统进行LQG最优控制设计。
2018年9月17日星期一 例. 已知可控直流电源供电给直流电机的系统结构图如图所示。对系统进行LQG最优控制设计。 系统性能指标为： 设系统模型为 其中 G = [ ]T, 且

40 2018年9月17日星期一 MATLAB程序如下： 运行后有 >> l8501

41 控制系统工具箱其他常用命令简介 1.系统生成命令 rss、drss、ord2 rss( ) 生成随机稳定的连续时间状态空间模型
2018年9月17日星期一 控制系统工具箱其他常用命令简介 1.系统生成命令 rss、drss、ord2 rss( ) 生成随机稳定的连续时间状态空间模型 sys= rss(n) 生成n阶SISO 状态空间模型. sys= rss(n,p) 生成n阶单输入p个输出的状态空间模型. sys= rss(n,p,m) 生成n阶m输入p个输出的状态空间模型. drss( ) 生成随机稳定的离散时间状态空间模型 sys= drss(n) 生成n阶SISO 状态空间模型. sys= drss(n,p) 生成n阶单输入p个输出的状态空间模型. sys= drss(n,p,m) 生成n阶m输入p个输出的状态空间模型. 在所有生成的离散时间状态空间模型中，采样时间均未指定

42 Continuous-time model. >> >> eig(sys.a) ans = -0.6488
2018年9月17日星期一 >> sys=rss(2,2,3) a = x x2 x x b = u u u3 x x c = y y d = y y Continuous-time model. >> >> eig(sys.a) ans = >> impulse(sys)

43 Sampling time: unspecified Discrete-time model. >> eig(sys.a)
2018年9月17日星期一 >> sys=drss(2,2,3) a = x x2 x x b = u u u3 x x c = y y d = u u u3 y y Sampling time: unspecified Discrete-time model. >> eig(sys.a) ans = i i >> step(sys)

44 说明：[A, B, C, D]=ord2(wn,z)生成由(A,B,C,D)描述的二
2018年9月17日星期一 ord2( ) 生成连续的二阶系统 格式：[A, B, C, D]=ord2(wn,z) [num, den]=ord2(wn,z) 说明：[A, B, C, D]=ord2(wn,z)生成由(A,B,C,D)描述的二 阶状态空间系统，其固有频率为wn( 即n ), 阻尼系 数为z(即 )传递函数为： [num, den]=ord2(wn,z)生成传递函数描述的二阶系 统的分子分母多项式num和den

45 例. 生成一个阻尼系数为= 0.4, 固有频率为n=2.4的二阶 系统。并画出系统的Bode图。
2018年9月17日星期一 例. 生成一个阻尼系数为= 0.4, 固有频率为n=2.4的二阶 系统。并画出系统的Bode图。 >> margin(sys) >> grid >> [num,den]=ord2(2.4,0.4) num = 1 den = >> sys=tf(num,den) Transfer function: s^ s

46 2.模型测试命令 eig、pole、damp、zero,
2018年9月17日星期一 2.模型测试命令 eig、pole、damp、zero, eig( ) 计算矩阵的特征值和特征向量 格式：E=eig(X) [V,D]=eig(X) 说明：E=eig(X)返回由矩阵X的特征值构成的向量E. [V,D]=eig(X)返回对角矩阵D和满矩阵V，使得 X*V=V*D, D的主对角线上元为特征值。 pole( ) 计算LTI系统的极点 格式：p=pole(sys) 说明：返回由矩阵sys.a的特征值构成的向量p. 多重极点 是数值敏感的，计算精度不能很高。

47 该命令生成一个3阶系统, 检验系统模态。 >> eig(g.a) >> g=rss(3,2,2)； ans =
2018年9月17日星期一 >> eig(g.a) ans = i i >> [V,D]=eig(g.a) V = i i i i D = i i >> g=rss(3,2,2)； 该命令生成一个3阶系统, 检验系统模态。 >> pole(g) ans = i i

48 damp( ) 计算LTI系统固有频率、阻尼系数和极点 格式：[wn,z]=damp(sys) [wn,z,p]=damp(sys)
2018年9月17日星期一 damp( ) 计算LTI系统固有频率、阻尼系数和极点 格式：[wn,z]=damp(sys) [wn,z,p]=damp(sys) 说明：返回LTI系统sys的自然频率wn, 阻尼因子z 和系 统极点构成的向量p. 当不带有返回参数时，命令 给出返回数据列表。对于采样时间为Ts 的离散系 统，该命令通过令z = exp(sTs ), 对等价的连续时 间系统进行计算，当不指定采样时间时，wn和z 为空的。

49 例. 计算并显示下列连续系统的自然频率，阻尼因子
2018年9月17日星期一 例. 计算并显示下列连续系统的自然频率，阻尼因子 >> H=tf([2 5 1],[1 2 3]) Transfer function: 2 s^2 + 5 s + 1 s^2 + 2 s + 3 >> damp(H) Eigenvalue Damping Freq. (rad/s) -1.00e e+000i e e+000 -1.00e e+000i e e+000 >> [wn,z,p]=damp(H) wn = 1.7321 z = 0.5774 p = i i

50 [z,gain]=zero(sys) %用于 SISO模型 说明：计算SISO系统的零点或MIMO系统的传输零点， 并返回零点构成的向量z.
2018年9月17日星期一 zero( ) 计算LTI系统传输零点 格式：z=zero(sys) [z,gain]=zero(sys) %用于 SISO模型 说明：计算SISO系统的零点或MIMO系统的传输零点， 并返回零点构成的向量z. 对于SISO系统， [z,gain] = zero(sys) 还返回零 极点模型意义下的增益gain. MIMO系统的传输零点是使得 发生降秩的C

51 例.计算下列系统的传输零点和极点。 >> s=tf('s'); >> pole(G)
2018年9月17日星期一 例.计算下列系统的传输零点和极点。 >> s=tf('s'); >> G=[1,2*s+1;1,s^2+4*s+2]/(s^2+3*s+2)； >> zero(G) ans = >> pole(G) ans = -2 -1

52 gensig 是用于命令lsim( )仿真的周期信号发生器。 调用格式1： [u,t] = gensig(type,tau)
2018年9月17日星期一 3.模型仿真命令 3.1 gensig、lsim gensig 是用于命令lsim( )仿真的周期信号发生器。 调用格式1： [u,t] = gensig(type,tau) 其中： u为标量信号 type = ‘sin’ 正弦波 tyep = ‘square’ 方波 type = ‘pulse’ 周期冲击 tau 为指定周期 t 信号u的采样时刻 所有的信号高度为一个单位。

53 调用格式2： [u,t] = gensig(type,tau,Tf,Ts) 进一步指定持续时间Tf和采样周期Ts
2018年9月17日星期一 调用格式2： [u,t] = gensig(type,tau,Tf,Ts) 进一步指定持续时间Tf和采样周期Ts 例. 生成一个周期为5秒，持续时间为30秒，采样周期为 0.1秒的方波 >> [u,t]=gensig('square',5,30,0.1); >> plot(t,u); >> axis([ ]) >>

54 lsim(sys1,’PlotStyle1’,…,sysN,’PlotStyleN’,u,t)
2018年9月17日星期一 lsim 对系统的任意输入进行仿真 格式: lsim(sys,u,t) lsim(sys,u,t,x0) lsim(sys1,sys2,…,sysN,u,t) lsim(sys1,’PlotStyle1’,…,sysN,’PlotStyleN’,u,t) [y,t,x]=lsim(sys,u,t,x0) 其中: u、t 是时间信号u (t ), 对于多输入系统信号维数要 与输入通道数相匹配； x0 是系统的初始状态； sys1,…,sysN 是同时被仿真的多个LTI对象； ‘PlotStyle1’,…, ’PlotStyleN’ 是plot命令支持的各种 属性标识符; 当命令返回仿真数据时，系统不产生图形。

55 其中：方波周期为4秒，持续时间为10秒，采样周期为0.1秒
2018年9月17日星期一 例. 计算下述系统的方波响应 其中：方波周期为4秒，持续时间为10秒，采样周期为0.1秒 >> [u,t]=gensig('square',4,10,0.1); >> h=tf({[2 5 1];[1 -1]},{[1 2 3];[1 1 5]}); >> lsim(h,u,t) >>

56 继续上例，计算系统在0～30秒对指数信号u (t ) = e –0.3t 的响应。
2018年9月17日星期一 继续上例，计算系统在0～30秒对指数信号u (t ) = e –0.3t 的响应。 >> t=0:0.1:10; >> u=exp(-0.3*t); >> h=tf({[2 5 1];[1 -1]},… {[1 2 3];[1 1 5]}); >> lsim(h,u,t)

57 考虑两个输入通道的系统的仿真：首先建立一个两输入系统.
2018年9月17日星期一 考虑两个输入通道的系统的仿真：首先建立一个两输入系统. >> sys=rss(2,2,2) a = x x2 x x b = u u2 x x c = y y d = u u2 y y >> u=[(exp(-0.3*t).*sin(3*t))',(exp(-0.1*t).*cos(2*t))']; >> lsim(sys,u,t)

58 Initial( ) 计算状态空间模型的初始条件响应 格式：initial(sys,x0) initial(sys,x0,t)
2018年9月17日星期一 3.2 initial Initial( ) 计算状态空间模型的初始条件响应 格式：initial(sys,x0) initial(sys,x0,t) initial(sys1,…,sysN, x0) initial(sys1,’PlotStyle1’,…,sysN,’PlotSyleN’,x0) [y,t,x]=initial(sys,x0) 说明：1)画出零输入的连续以及离散状态空间模型sys (sys1, …, sysN)对初始条件x0的时间响应曲线. initial(sys,x0,t)可以用两种方式指定仿真时间t t=T仿真终止时间 t=0:dt(仿真步长):T仿真终止时间 2) 对于离散系统，仿真步长要与采样周期相匹配。

59 例. 画出下列状态空间模型的初值响应 3) initial(sys1,…,sysN, x0) 将多个系统的仿真结果画在一张 图上.
2018年9月17日星期一 3) initial(sys1,…,sysN, x0) 将多个系统的仿真结果画在一张 图上. 4) [y,t,x]=initial(sys,x0) 返回仿真数据，不产生图形，其中 y 是输出响应数据向量，x是状态轨迹向量。 例. 画出下列状态空间模型的初值响应

60 >> sys=ss(a,[],c,[]) >> initial(sys,x0)
2018年9月17日星期一 >> a=[ , ; , 0]; >> c=[1.9691,6.4493]; >> x0=[1;0]; >> sys=ss(a,[],c,[]) >> initial(sys,x0) >> [y,t,x]=initial(sys,x0); >> subplot(2,1,1); >> plot(t,x(:,1),t,x(:,2)) >> subplot(2,1,2); >> plot(x(:,1),x(:,2))

61 evalfr() 计算LTI模型在单个(复)频率点的响应 格式：frsp=evalfr(sys,f)
2018年9月17日星期一 4.模型频域响应命令 4.1 evalfr、freqresp evalfr() 计算LTI模型在单个(复)频率点的响应 格式：frsp=evalfr(sys,f) 说明：计算LTI模型sys (tf,ss,zpk模型)的传递函数在频率点f的 值：H(f ) = D + C (f I –A)–1B freqresp() 计算LTI模型在频率格点(1, 2, … , n)的响应 格式：H=freqresp(sys,w) 说明：计算LTI模型sys在由向量w指定的频率格点的响应阵列 H, H是一个3D数据阵列, 如H (i,j,k)为传递函数H(s)的第 j个输入到第i 个输出在 w(k)= k 点的值Hi j (k )。频率 为实数, 单位为:弧度/秒. 对采样时间为Ts 的离散系统， 取z = exp(jkTs )。

62 例. 计算下列系统在频率点  =1, 10, 100的频率响应. >> H=freqresp(P,w) H(:,:,1) =
2018年9月17日星期一 例. 计算下列系统在频率点  =1, 10, 100的频率响应. >> H=freqresp(P,w) H(:,:,1) = i i H(:,:,2) = i i H(:,:,3) = i i >> w=[ ]; >> s=tf('s'); >> P=[0,1/(s+1);(s-1)/(s+2),1]

63 sigma() 计算LTI模型奇异值的频率响应 格式：sigma(sys) sigma(sys,w)
2018年9月17日星期一 4.2 sigma sigma() 计算LTI模型奇异值的频率响应 格式：sigma(sys) sigma(sys,w) sigma(sys,w,type) % 有相同输入输出个数 sigma(sys1,…,sysN) sigma(sys1,…,sysN,w) sigma(sys1,…,sysN,w,type) % 有相同输入输出个数 sigma(sys1,’PlotStyle1’,…,sysN,’PlotStyleN’) [sv,w]=sigma(sys) [sv]=sigma(sys,w)

64 说明：1)计算LTI模型sys的奇异值的频率响应，它是SISO系 统 Bode图的推广。对于H(s), 在s =j点的奇异值为
2018年9月17日星期一 说明：1)计算LTI模型sys的奇异值的频率响应，它是SISO系 统 Bode图的推广。对于H(s), 在s =j点的奇异值为 H(j)* H(j) 的特征值的平方根。最大奇异值对系统有重要影响。 2) 对于采样周期为Ts 的离散系统，取H (exp(jTs ))。 3) 对于输入输出个数相等的系统，频率响应矩阵H(j) 是方阵 type=1, 计算频率响应H(j)的逆H(j)–1的奇异值; type=2,计算频率响应I +H(j) 的奇异值; type=3,计算频率响应I + H(j)–1的奇异值. 4) 可以将多个系统的奇异值响应画在同一图上。 5) 当带有返回响应数据时，不产生图形。 6) w为用户指定的频率点向量 (radians/second)。

65 例. 画出系统H(s) 和I +H(s)的奇异值响应。
2018年9月17日星期一 例. 画出系统H(s) 和I +H(s)的奇异值响应。 >> s=tf('s') >> H=[0,(3*s)/(s^2+s+10);… (s+1)/(s+5),(2)/(s+6)]; >> subplot(211) >> sigma(H) >> grid >> subplot(212) >> sigma(H,[],2)