2018年12月27日星期四 第八章 状态空间分析的MATLAB 实现

xt(t ) = Tx(t )，或 xt(n ) = Tx(n ) 2018年12月27日星期四 变换后系统为 1. 系统实现的基本概念 1.1 坐标变换 设LTI系统的状态方程为： 或 其中 ARnn, B Rnm, C Rln, D Rlm. 或 其中： At = TAT –1, Bt = TB, Ct = CT –1, Dt = D 设T Rnn 为非奇异矩阵，则可定义如下坐标变换 xt(t ) = Tx(t )，或 xt(n ) = Tx(n ) 坐标变换不改变系统的传递函数

1.2 系统的能控性、能控标准型实现及系统的能控与不能控分解 若系统完全能控，则必然存在一个坐标变换T, 使得 2018年12月27日星期四 1.2 系统的能控性、能控标准型实现及系统的能控与不能控分解 若系统完全能控，则必然存在一个坐标变换T, 使得 设LTI系统状态空间模型由矩阵(A, B, C , D )描述 Gram 矩阵判据：系统完全能控的充分必要条件是存在时刻t1 > 0, 使得Gram矩阵 非奇异。 det(sI A)=sn + a1sn – 1 +  + an – 1s + an 当rankCo < n 时，存在坐标变换T，对系统进行能控和不能控分解 秩判据：连续LTI系统完全能控的充分必要条件为系统的可控性矩阵 的秩为n , 即rank Co = n。 则(Ac , Bc )构成系统的能控子空间。

1.3 系统的能观性、能观标准型实现及系统的能观与不能观分解 若系统完全能观，则必然存在一个坐标变换T, 使得 2018年12月27日星期四 1.3 系统的能观性、能观标准型实现及系统的能观与不能观分解 若系统完全能观，则必然存在一个坐标变换T, 使得 设LTI系统状态空间模型由矩阵(A, B, C , D )描述 Gram 矩阵判据：系统完全能观的充分必要条件是存在时刻t1 > 0, 使得Gram矩阵 非奇异。 det(sI -A)=sn + a1sn – 1 +  + an – 1s + an 当rankOb < n 时，存在坐标变换T，对系统进行能控可不能控分解 秩判据：连续LTI系统完全能控的充分必要条件为系统的可控性矩阵 的秩为n , 即rankOb = n。 则(Ao , Co )构成系统的能观子空间。

Jordan标准型是通过一个相似变换T，使系统的A 矩阵变为： 2018年12月27日星期四 1.4系统的Jordan标准型实现 Jordan标准型是通过一个相似变换T，使系统的A 矩阵变为： 称为矩阵A的模态矩阵(model matrix ), 而Ji 称为系统的Jordan子矩阵，且有： i 为矩阵A的特征值。 1.5系统的可控可观分解及最小实现 通过相似变换，可将系统(A, B, C, D )等效地变换为下述规范形式： 该分解称为Kalman分解，而可控可观子空间称为原系统的最小实现。

连续LTI系统(A, B, C, D) 的可控可观Gramian矩阵为： 2018年12月27日星期四 1.6 系统的平衡实现 连续LTI系统(A, B, C, D) 的可控可观Gramian矩阵为： 且满足Lyapunov方程： 从而存在一个坐标变换将系统变换为(Ab , Bb , Cb , Db ), 且满足如下Lyapunov方程： 为对角阵。 (Ab , Bb , Cb , Db )称为系统的平衡实现。 1.7 系统的平衡降阶实现 设系统平衡实现中的可控可观Gram矩阵为＝diag(1, 2), 将1, 2的对角元从大到小排列有： 1＝diag(1 ,  ,r ), 2＝diag(r +1 ,  ,n ) 且r >> r +1 , 从而相应地将系统分解为： 截掉对应对应2的子系统，则降阶系统可写成 称为系统的平衡降阶实现。

2.状态空间实现函数 状态空间实现函数列表 函数名称 功 能 ss2ss 坐标变换 canon 状态空间的正则实现 ctrb 可控矩阵计算 2018年12月27日星期四 2.状态空间实现函数 状态空间实现函数列表 函数名称 功 能 ss2ss 坐标变换 canon 状态空间的正则实现 ctrb 可控矩阵计算 ctrbf 系统的可控与不可控分解 gram 求系统的可控与可观gramian矩阵 obsv 可观矩阵计算 obsvf 系统的可观与不可观分解 ssbal 状态空间的对角平衡实现 balreal 状态空间的平衡实现 minreal 状态空间的最小实现 modred 模型降阶

调用格式：sysT=ss2ss(sys, T) 其中：T为坐标变换矩阵 sys为LTI对象 sysT为变换后的LTI对象 2018年12月27日星期四 2.1 ss2ss 功能： 坐标变换 调用格式：sysT=ss2ss(sys, T) 其中：T为坐标变换矩阵 sys为LTI对象 sysT为变换后的LTI对象 例1.已知系统的状态空间表示为 >>ch8eg1 a = 0 1 -2 -3 b = 1 c = 6 0 d = >> 设坐标变换矩阵为 求变换后系统状态方程

变换后的LTI对象 >>[t,v]=eig(a)； >> t=inv(t)； 2018年12月27日星期四 变换后的LTI对象 >>[t,v]=eig(a)； >> t=inv(t)； >> ss2ss(sysT,t) a = x1 x2 x1 -1 0 x2 0 -2 b = u1 x1 1.414 x2 2.236 c = x1 x2 y1 4.243 -2.683 d = y1 0 Continuous-time model. >> >> sysT a = x1 x2 x1 0 1 x2 -2 -3 b = u1 x1 0 x2 1 c = y1 6 0 d = y1 0 Continuous-time model. >>

格式：[csys,T]=canon(sys, ‘type’) 其中： sys: 为LTI对象; ‘type’: 可取字符串为: 2018年12月27日星期四 2.2 canon 功能：状态空间的正则实现 格式：[csys,T]=canon(sys, ‘type’) 其中： sys: 为LTI对象; ‘type’: 可取字符串为: ‘modal’－计算模态矩阵正则 实现, ‘companion’－计算伴随阵正 则实现, csys: 为变换后的LTI对象; T: 为坐标变换矩阵，当给出的 LTI对象不是ss对象时，T为 空矩阵。 注意: 1. 模态实现要求A可对角化。 2. 伴随实现要求关于第一输 入可控。 例2.考虑如下线性系统 计算系统的模态正则实现以及伴随正则实现 >> a=[-4 -1.5 -0.5 -0.125;… 4 0 0 0;0 2 0 0;0 0 1 0]; >> b=[8;0;0;0]; >> c=[-8.074 -2.897 0.2116 -0.4306]; >> d=0.7538; >> sys=ss(a,b,c,d); >> [msys, mT]=canon(sys,'modal'); >> [csys, cT]=canon(sys,'companion');

>> msys.a >> csys.a ans = ans = -1.0001 0.0001 0 0 2018年12月27日星期四 >> msys.a ans = -1.0001 0.0001 0 0 -0.0001 -1.0001 0 0 0 0 -0.9999 0.0001 0 0 -0.0001 -0.9999 >> mT mT = 1.0e+011 * -7.1374 -5.3530 -2.6765 -0.8922 -7.1379 -5.3530 -2.6763 -0.8920 -7.1370 -5.3528 -2.6764 -0.8921 7.1383 5.3542 2.6773 0.8925 >> csys.a ans = 0 0 0 -1 1 0 0 -4 0 1 0 -6 0 0 1 -4 >> cT cT = 0.1250 0.1250 0.0938 0.0625 0 0.0313 0.0625 0.0938 0 0 0.0156 0.0625 0 0 0 0.0156

例3.计算例2系统的可控性矩阵 2.3 ctrb 功能：可控矩阵计算 格式：Co=ctrb(A,B) Co=ctrb(sys) 其中： 2018年12月27日星期四 例3.计算例2系统的可控性矩阵 2.3 ctrb 功能：可控矩阵计算 格式：Co=ctrb(A,B) Co=ctrb(sys) 其中： A,B: 为状态方程矩阵； sys: 为LTI对象; Co : 为可控性矩阵。 说明：该函数返回状态方程： (d/dt) x = Ax +Bu 的可控性矩阵 Co = [B AB A2B … An–1B ] >> a=[-4 -1.5 -0.5 -0.125;… 4 0 0 0;0 2 0 0;0 0 1 0]; >> b=[8;0;0;0]; >> Co=ctrb(a,b) Co = 8 -32 80 -160 0 32 -128 320 0 0 64 -256 0 0 0 64 >> [b a*b a^2*b a^3*b] ans = 0 0 64 -256 0 0 0 64 >>rank(ctrb(a,b)) 4

格式: [Ab,Bb,Cb,T,k]=ctrbf(A,B,C) [Ab,Bb,Cb,T,k]=ctrbf(A,B,C,tol) 其中： 2018年12月27日星期四 2.4 ctrbf 功能：系统的可控与不可控分解 格式: [Ab,Bb,Cb,T,k]=ctrbf(A,B,C) [Ab,Bb,Cb,T,k]=ctrbf(A,B,C,tol) 其中： A,B,C: 为系统矩阵； tol: 误差容限，缺省为： 10*n*norm(a,1)*eps Ab,Bb,Cb: 为分解后系统矩阵; T : 为分解的变换矩阵； k: 长度为n的向量，其分 量之和为可控矩阵的秩 说明：该函数返回(A,B,C )的可控 与不可控分解： 例3.将如下系统进行可与不可控分解 >> a=[1 1 ;4 -2]; >> b=[1 -1;1 -1]; >> c=[1 0;0 1]; >> [Ab,Bb,Cb,T,k]=ctrbf(a,b,c) Ab = -3.0000 0.0000 3.0000 2.0000 Bb = 0 0 -1.4142 1.4142 Cb = -0.7071 -0.7071 0.7071 -0.7071

‘c’ :字符, 表示计算可控gramian阵; ‘o’ :字符, 表示计算可观gramian阵; Wc: 返回的可控gramian阵; 2018年12月27日星期四 T = -0.7071 0.7071 -0.7071 -0.7071 k = 1 0 2.5 gram 功能：计算可控与可观gramian阵 格式: Wc=gram(sys,’c’) Wo=gram(sys,’o’) 其中： sys: 状态空间LTI对象； ‘c’ :字符, 表示计算可控gramian阵; ‘o’ :字符, 表示计算可观gramian阵; Wc: 返回的可控gramian阵; Wo: 返回的可控gramian阵. 说明：通过求解Lyapunov方程的解 i) 连续时间系统： AWc + WcAT + BBT = 0 ATWo + WoA + CTC = 0 ii) 离散时间系统： AWcAT – Wc + BBT = 0 ATWoA – Wo + CTC = 0 分解后系统为： 分解表明该系统有一个状态不可控。

例4.求系统的可控可观gramian阵 说明系统稳定. >> Wc=gram(sys,'c') Wc = 0.0313 0 0 2018年12月27日星期四 例4.求系统的可控可观gramian阵 说明系统稳定. >> Wc=gram(sys,'c') Wc = 0.0313 0 0 0 0 0 >> Wo=gram(sys,'o') Wo = 0.0098 0.3132 0.0153 0.3132 16.9927 1.3466 0.0153 1.3466 0.2940 >>A*Wc+Wc*A'+B*B' ans = 0 0 0 >> A=[-1 64.9 -2.44;0 -0.8 0.2;… 0 -5 -0.8]; >> B=[0.25;0;0]; >> C=[0.14 0.03 -0.08]; >> D=0; >> sys=ss(A,B,C,D); >> eig(sys.a) ans = -1.0000 -0.8000 + 1.0000i -0.8000 - 1.0000i

2.6 obsv 功能：可观矩阵计算 格式: Ob=obsv(A,C) Ob=obsv(sys) 其中： 2018年12月27日星期四 2.6 obsv 功能：可观矩阵计算 格式: Ob=obsv(A,C) Ob=obsv(sys) 其中： A,C: 为系统状态及输出方程矩阵； sys: 状态空间LTI对象； Ob: 系统可观性矩阵。 说明： 若rank Ob= n，则系统可观测。 例5.求系统的可观性阵 >> A=[-1 64.9 -2.44;0 -0.8 0.2;… 0 -5 -0.8]; >> B=[0.25;0;0]; >> C=[0.14 0.03 -0.08]; >> D=0; >> sys=ss(A,B,C,D);

格式: [Ab,Bb,Cb,T,k]=obsvf(A,B,C) [Ab,Bb,Cb,T,k]=obsvf(A,B,C,tol) 其中： 2018年12月27日星期四 >> Ob=obsv(A,C) Ob = 0.1400 0.0300 -0.0800 -0.1400 9.4620 -0.2716 0.1400 -15.2976 2.4513 >> Ob=obsv(sys) >> rank(Ob) ans = 3 2.7 obsvf 功能：系统的可观与不可观分解 格式: [Ab,Bb,Cb,T,k]=obsvf(A,B,C) [Ab,Bb,Cb,T,k]=obsvf(A,B,C,tol) 其中： A,B,C: 为系统矩阵； tol: 误差容限，缺省为： 10*n*norm(a,1)*eps Ab,Bb,Cb: 为分解后系统矩阵; T : 为分解的变换矩阵； k: 长度为n的向量，其分 量之和为可观矩阵的秩 说明：该函数返回(A,B,C )的可观 与不可观分解： 说明系统可观测。

例6.求系统的可观与不可观性分解 变换后的系统为： 分解表明系统有一个状态不可观测。 T = 0.7071 0.7071 2018年12月27日星期四 例6.求系统的可观与不可观性分解 T = 0.7071 0.7071 -0.7071 0.7071 k = 1 0 >> A=[ 1 4;1 -2];B=eye(2);C=[1 -1]; >> [Ab,Bb,Cb,T,k]=obsvf(A,B,C) Ab = 2.0000 -3.0000 0.0000 -3.0000 Bb = 0.7071 0.7071 -0.7071 0.7071 Cb = 0 -1.4142 变换后的系统为： 分解表明系统有一个状态不可观测。

格式: [sysb,T]=ssbal(sys) [sysb,T]=ssbal(sys,condT) 其中： 2018年12月27日星期四 2.8 ssbal 功能：状态空间的对角平衡实现 格式: [sysb,T]=ssbal(sys) [sysb,T]=ssbal(sys,condT) 其中： sys： 为状态空间表示的LTI对象; sysb: 返回的状态空间平衡实现； T : 为对角形变换矩阵； condT : 为指定T的条件数上界,缺 省情况下为inf。 说明：计算对角形坐标变换矩阵T,使得： 的行和列有近似相等的范数。 2.9 balreal 功能：状态空间的平衡实现 格式: sysb =balreal(sys) [sysb,g,T,Ti]=balreal(sys) 其中： sys： 为状态空间表示的LTI对象; sysb: 返回的状态空间平衡实现； g: 平衡实现gramian阵对角元 排成的向量; T : 变换矩阵； Ti : 为T的逆. 说明：该平衡实现主要用于模型降阶, 该函数可同时用于连续或离散系统, 对于传递函数模型, 则首先进行模型转换.

例7.求下述系统的平衡实现 b = u1 x1 0.125 x2 0.5 x3 32 c = x1 x2 x3 2018年12月27日星期四 b = u1 x1 0.125 x2 0.5 x3 32 c = x1 x2 x3 y1 0.8 20 3.125 d = y1 0 Continuous-time model. T = 0.1250 0 0 0 0.5000 0 0 0 32.0000 例7.求下述系统的平衡实现 >> a=[1 1e4 1e2;0 1e2 1e5;10 1 0]; >> b=[1;1;1]; >> c=[0.1 10 1e2]; >> sys=ss(a,b,c,0); >> [ssy,T]=ssbal(sys) a = x1 x2 x3 x1 1 2500 0.3906 x2 0 100 1563 x3 2560 64 0

gramian矩阵表明后面两个状态可以删除. 例8. 求下述系统的平衡实现, 并对 系统进行降阶. 2018年12月27日星期四 gramian矩阵表明后面两个状态可以删除. 例8. 求下述系统的平衡实现, 并对 系统进行降阶. >> sysr=modred(sysb,[2 3],'del') %删除后两个状态 >> zpk(sysr) zero/pole/gain: 1.0001 -------- (s+4.97) >> bode(sys,'-',sysr,'x') >> grid >> sys=zpk([-10 -20.01],… [-5 -9.9 -20.1],1) Zero/pole/gain: (s+10) (s+20.01) ---------------------- (s+5) (s+9.9) (s+20.1) >> [sysb,g]=balreal(sys); >> g' ans = 0.1006 0.0001 0.0000

格式: sysr=minreal(sys) sysr=minreal(sys,tol) 其中： sys： 为LTI对象; 2018年12月27日星期四 2.10 minreal 功能：状态空间的最小实现 格式: sysr=minreal(sys) sysr=minreal(sys,tol) 其中： sys： 为LTI对象; tol: 为传递函数零极相消的容许 误差; sysr: 返回最小实现. 说明： 删除状态空间模型中不可控和不可观的状态; 或者对消零极点模型中相同的零极点对, 输出系统具有最小的阶次. 例9. 求下述闭环系统的最小实现. – >> g=zpk([],1, 1)； >> h=tf([2 1],[1 0])； >> cloop=inv(1+g*h)*g zero/pole/gain: s (s-1) -------------------- (s-1) (s^2 + s + 1) >> cloop=minreal(cloop) s -------------- (s^2 + s + 1)

格式: rsys=modred(sys,elim) rsys=modred(sys,elim,’mdc’) 2018年12月27日星期四 2.11 modred 功能：模型降阶 格式: rsys=modred(sys,elim) rsys=modred(sys,elim,’mdc’) rsys=modred(sys,elim,’del’) 其中： sys： 为LTI对象; elim: 给出需要删除的状态索引; ‘mdc’: 降阶后系统与原系统有相同 的直流增义(DC gian); ‘del’: 简单删除需要删除的状态, 不保证降阶后系统与原系统 有相匹配的直流增义, 但在 频域性能与原系统更接近. 说明： 用于降低系统的维数, 常与balreal 函数连用. 例10. 系统的传递函数为: 降阶前先要进行降阶实现. >> h=tf([1 11 36 26],… [1 14.6 74.96 153.7 99.65]); >> [hb,g]=balreal(h); >> g' ans = 0.1394 0.0095 0.0006 0.0000 结果表明gramian阵对角元后面三个元素相对较小, 因此可以将其对应的状态删除. 同时比较两种方法.

>> hmdc=modred(hb,2:4,'mdc'); 2018年12月27日星期四 >> hmdc=modred(hb,2:4,'mdc'); >> hdel=modred(hb,2:4,'del'); >> bode(h,'-',hmdc,'x',hdel,'*') >> step(h,'-', hmdc,'-.', hdel,'--')

3.状态空间表达式及状态方程的解 与传递函数的关系 对状态方程进行拉普拉斯变换 3.1 状态空间表达式的建立 状态空间方程: 2018年12月27日星期四 3.状态空间表达式及状态方程的解 与传递函数的关系 对状态方程进行拉普拉斯变换 3.1 状态空间表达式的建立 状态空间方程: 假定 x (0) = 0, 则有 从而: G(s) = C(sI –A)–1B + D 且记为 状态空间方框图

输入的导数出现在微分方程系统中的状态空间表示. 2018年12月27日星期四 输入的导数出现在微分方程系统中的状态空间表示. 考虑如下动力学方程: 以及 可得状态空间方程 方法1 令

2018年12月27日星期四 方法2 令 例11. 考虑如下质量-弹簧-阻尼系统. 假定t < 0时小车静止, u 是小车的位移也是系统的输入, 在t =0时刻, 小车以恒定速度移动, 即du/dt = constant, 质量块的位移为系统的输出. 建立系统的状态空间方程. 可得状态空间方程 系统的动力学方程为:

取: m =10kg, b = 20 N-s/m, k = 100 N/m, 在斜坡输入du/dt =1m/s下观察系统响应 2018年12月27日星期四 令: >> t=0:0.01:4; >> A=[0 1;-10 -2]; >>B=[2;6]; >>C=[1 0]; >>D=0; >> sys=ss(A,B,C,D); >> u=t; >> lsim(sys,u,t); >>grid; >> title('Unit-Ramp Response… (Method 1)') >> xlabel('t') >> ylabel('Output y and Unit-Ramp… Input u') >> text(0.85,0.25,'y') >> text(0.15,0.8,'u') 由方法1可得: 取: m =10kg, b = 20 N-s/m, k = 100 N/m, 在斜坡输入du/dt =1m/s下观察系统响应

2018年12月27日星期四 类似地, 用方法2: 令: 可得方程 带入数据后:

>> sys=ss(A,B,C,D); >> u=t; >> lsim(sys,u,t); 2018年12月27日星期四 >> A=[0 1;-10 -2]; >> B=[0;1]; >> C=[10 2]; >> D=0; >> sys=ss(A,B,C,D); >> u=t; >> lsim(sys,u,t); >> title('Unit-Ramp Response… (Method 1)') >> xlabel('t') >> ylabel('Output y and Unit-Ramp… Input u') >> text(0.85,0.25,'y') >> text(0.15,0.8,'u') >> grid

其中w1, w2 为外界扰动, u为控制作用. 设m1 = m2 =1, 其动力学方程为: 2018年12月27日星期四 则写成状态空间的表示为: 例12. 考虑如下由弹簧连接的两个无 摩擦小车系统: 则系统矩阵为 其中w1, w2 为外界扰动, u为控制作用. 设m1 = m2 =1, 其动力学方程为:

设k = 0.5~2.0, 如取k = 1.25, 建立ss对象 系统需要控制器镇定. >> step(sys) 2018年12月27日星期四 设k = 0.5~2.0, 如取k = 1.25, 建立ss对象 >> step(sys) >> a=[0 0 1 0;0 0 0 1;-1.25 1.25 0 0; 1.25 –1.25 0 0]; >> b=[0 0;0 0 ;1 0;-1 0]; >> b=[0 0;0 0 ;1 0;0 -1]; >> c=[0 1 0 0];d=0 >> sys=ss(a,b,c,d) >> eig(sys.a) ans = 0.0000 0.0000 + 1.5811i 0.0000 - 1.5811i -0.0000 开环系统有两个纯虚特征值和两个零特征值, 系统不稳定, 阶跃响应如图所示: 系统需要控制器镇定.

说明该系统是可控可观的, 设控制器状态方程为: 2018年12月27日星期四 >> rank(ctrb(sys)) ans = 4 >> rank(obsv(sys)) >> Ac=[0 -0.7195 1 0; 0 -2.9732 0 1; -2.5133 4.8548 -1.7287 -0.9616; 1.0063 -5.4097 -0.0081 0.0304]; >>Bc=[0.720;2.973;-3.37;4.419]; >>Cc=[-1.506 0.494 -1.738 -0.932]; >> Dc=0; >> sysc=ss(Ac, Bc, Cc,Dc); >> csys=feedback(sys,sysc,1,1,+1); >>step(csys) 说明该系统是可控可观的, 设控制器状态方程为: 其中

2018年12月27日星期四 系统的simulink模型如图 其中

提取闭环系统线性模型 >> [A B,C,D]=linmod2('Two_m_one_sp_system'); 2018年12月27日星期四 提取闭环系统线性模型 >> [A B,C,D]=linmod2('Two_m_one_sp_system'); >> clsys=ss(A,B,C,D); >> step(clsys)