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

Presentation on theme: "2018年12月27日星期四 第八章 状态空间分析的MATLAB 实现."— Presentation transcript:

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

2 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 ) 坐标变换不改变系统的传递函数

3 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 )构成系统的能控子空间。

4 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 )构成系统的能观子空间。

5 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分解，而可控可观子空间称为原系统的最小实现。

6 连续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的子系统，则降阶系统可写成 称为系统的平衡降阶实现。

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

8 调用格式：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 = b = 1 c = d = >> 设坐标变换矩阵为 求变换后系统状态方程

9 变换后的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 x x b = u1 x x c = x x2 y d = y1 0 Continuous-time model. >> >> sysT a = x1 x2 x x b = u1 x1 0 x2 1 c = y d = y1 0 Continuous-time model. >>

10 格式：[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=[ ;… ; ; ]; >> b=[8;0;0;0]; >> c=[ ]; >> d=0.7538; >> sys=ss(a,b,c,d); >> [msys, mT]=canon(sys,'modal'); >> [csys, cT]=canon(sys,'companion');

11 >> msys.a >> csys.a ans = ans = -1.0001 0.0001 0 0
2018年12月27日星期四 >> msys.a ans = >> mT mT = 1.0e+011 * >> csys.a ans = >> cT cT =

12 例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=[ ;… ; ; ]; >> b=[8;0;0;0]; >> Co=ctrb(a,b) Co = >> [b a*b a^2*b a^3*b] ans = >>rank(ctrb(a,b)) 4

13 格式: [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 = Bb = Cb =

14 ‘c’ :字符, 表示计算可控gramian阵; ‘o’ :字符, 表示计算可观gramian阵; Wc: 返回的可控gramian阵;
2018年12月27日星期四 T = k = 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 分解后系统为： 分解表明该系统有一个状态不可控。

15 例4.求系统的可控可观gramian阵 说明系统稳定. >> Wc=gram(sys,'c') Wc = 0.0313 0 0
2018年12月27日星期四 例4.求系统的可控可观gramian阵 说明系统稳定. >> Wc=gram(sys,'c') Wc = >> Wo=gram(sys,'o') Wo = >>A*Wc+Wc*A'+B*B' ans = >> A=[ ; ;… ]; >> B=[0.25;0;0]; >> C=[ ]; >> D=0; >> sys=ss(A,B,C,D); >> eig(sys.a) ans = i i

16 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=[ ; ;… ]; >> B=[0.25;0;0]; >> C=[ ]; >> D=0; >> sys=ss(A,B,C,D);

17 格式: [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 = >> 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 )的可观 与不可观分解： 说明系统可观测。

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

19 格式: [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的逆. 说明：该平衡实现主要用于模型降阶, 该函数可同时用于连续或离散系统, 对于传递函数模型, 则首先进行模型转换.

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

21 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([ ],… [ ],1) Zero/pole/gain: (s+10) (s+20.01) (s+5) (s+9.9) (s+20.1) >> [sysb,g]=balreal(sys); >> g' ans =

22 格式: 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)

23 格式: 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([ ],… [ ]); >> [hb,g]=balreal(h); >> g' ans = 结果表明gramian阵对角元后面三个元素相对较小, 因此可以将其对应的状态删除. 同时比较两种方法.

24 >> 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,'--')

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

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

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

28 取: 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下观察系统响应

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

30 >> 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

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

32 设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=[ ; ; ; 1.25 – ]; >> b=[0 0;0 0 ;1 0;-1 0]; >> b=[0 0;0 0 ;1 0;0 -1]; >> c=[ ];d=0 >> sys=ss(a,b,c,d) >> eig(sys.a) ans = 0.0000 i i 开环系统有两个纯虚特征值和两个零特征值, 系统不稳定, 阶跃响应如图所示: 系统需要控制器镇定.

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

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

35 提取闭环系统线性模型 >> [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)