(1)线系统结构分析与分解及标准型.docx

实验三十七 线性系统结构分析与分解及标准型实验类型：验证 难度系数：0.3 实验性质：必做 课内学时：0 课外学时：2 分组人数：2开课方式：在课外完成在MATLAB平台上完成实验。实验目的：掌握线性系统状态空间标准型、解及其模型转换。实验设备与软件：1、MATLAB数值分析软件实验原理：1、标准型变换、矩阵Jordan型变换阵、特征值(1)标准型变换 命令格式 csys=canon(sys,type)说明：type制定规范型的形式，包括两种选项：model(模态规范型)、companion(伴随规范型,友矩阵型,能控II型)。(2)矩阵的Jordan规范型 命令格式 V J=Jordan(A)说明：V特征向量，J是Jordan型 (3)求矩阵特征值和特征向量 命令格式 V J=eig(A) cv= eig(A) 说明：V特征向量，J是Jordan型;cv是特征值列向量2、状态模型的相似变换： 命令格式 sysb=ss2ss(sys,T) 说明：sys是状态空间模型，T是非奇异变换阵的逆阵 传递函数模型与状态空间模型之间的相互转换：命令格式 A,B,C,D=tf2ss(num,den)num,den=ss2tf(A,B,C,D,iu) 说明：tf2ss:传递函数状态空间; ss2tf状态空间传递函数; iu是第iu个输入有效 zpk模型与状态空间模型之间的相互转换：命令格式：A,B,C,D=zp2ss(z,p,k)z,p,k=ss2zp(A,B,C,D,iu) 说明：zp2ss :zpk模型状态空间模型; ss2zp状态空间模型zpk模型3、线性定常系统的可控制与可观性及结构分解关于可控制与可观性及结构分解的理论内容，请看教材学习。(1)可控性和可观性-一般采用能控性/能观性矩阵类别(适用于离散或连续的情况)状态可控性和输出可控性子函数如下：function str =pdctrb(A,B)Qc=ctrb(A,B);r=rank(Qc);l=length(A);if r=l str=系统是状态完全可控的!else str=系统不是状态完全可控的!end实际上，rank(Qc)为系统的状态可控性指数，即系统中可控的状态的数目。function str=pdctrbo(A,B)%输出可控性Co=ctrb(A,B);m=size(C,1);%返回行数Qyc=C*Co,D;Tm=rank(Qyc);if m=Tm str=系统输出是完全可控的!else str=系统输出不是完全可控的!end实际上，rank(Qyc) 为系统的输出可控性指数，即系统中可控的输出的数目。状态可观性判别子函数代码如下：function str=pdobsv(A,C)Qo=obsv (A,C);r=rank(Qo);l= size(A,1);if r=l str=系统是状态完全可观的!else str=系统不是状态完全可观的!end实际上，rank(Qo)为系统的状态可观性指数，即系统中可观测的状态的数目。有了上述子函数，在Matlab中可以直接调用这些子函数来判断可控性和可观性。(2)可控性和可观性Gram矩阵可由下面的函数求得 W=gram(sys,type)-sys是系统的状态空间模型，type可以是c或o。通过判断W的正定性判定其可控性或可观性。由于从W的数学表达式上看，Wc和Wo是对称的半正定矩阵，它们分别满足下面的Lyapunov方程(这种方程在稳定性分析中还将提到)、所以系统必须稳定才能得到Gram矩阵。Wc矩阵中的值对应于输入信号对相应状态的贡献：第i个元素越大，则说明输入信号对第i个输入状态的贡献越大。Wo矩阵中的值对应于每个状态对系统输出的贡献：第i个元素越大，则说明第i状态对系统输出的贡献越大。(3)结构分解a.在matlab中调用ctrbf()函数对系统按能控性分解Abar,Bbar,Cbar,T,K=ctrbf(A,B,C)Abar,Bbar,Cbar,T,K=ctrbf(A,B,C,TOL)说明：K是可控的状态个数，TOL为误差，这里的变换是这样令的：x=T-1xbar，所以有Abar=TAT-1，Bbar=TB，Cbar=CT-1。显然这与教材中的令法是不同的，且T的选取方式也不同，得到的分解形式也不同，要加以区别。这里的能观性分解形式为：，b.在matlab中调用ctrbf()函数对系统按能观性分解Abar,Bbar,Cbar,T,K=obsvf(A,B,C)Abar,Bbar,Cbar,T,K= obsvf(A,B,C,TOL)说明：K是可观的状态个数，TOL为误差，这里的变换是这样令的：x=T-1xbar，所以有Abar=TAT-1，Bbar=TB，Cbar=CT-1。显然这与教材中的令法是不同的，且T的选取方式也不同，得到的分解形式也不同，要加以区别。这里的能观性分解形式为：，c.按能控、能观性分解-Kalman分解Kalman分解是先按能控性进行分解，然后对两部分状态再按能观性进行分解，最后得到能控能观部分、能控不能观部分、不能控能观部分、不能控不能观部分。在Matlab中没有提供对应的函数对系统进行直接的Kalman分解，这里我们给出一个函数实现分解：function Gk,T,K=kalmdec(G)G=ss(G);A=G.a;B=G.b;C=G.c;Ac,Bc,Cc,Tc,Kc=ctrbf(A,B,C);nc=rank(ctrb(A,B);n=length(A);ic=n-nc+1:n;Ao1,Bo1,Co1,To1,Ko1=obsv(Ac(ic,ic),Bc(ic),Cc(ic);if ncj Tea(i,j)=cp(i-(j-1);endendendT=Ttt*Tea;Abar=inv(T)*A*T;Bar=inv(T)*B;Cbar=C*T,Dar=D;Gss=ss(Abar,Bbar,Cbar,Dbar)(3)能观标准型-限于SISO系统若系统能观，则可转换成能观标准I型和II型。转换成能观标准I型代码:A=;B=;C=;D=;Gs=ss(A,B,C,D);Tinv=obsv (Gs.a,Gs.b);T=inv(Tinv);Abar=inv(T)*A*T;Bar=inv(T)*B;Cbar=C*T,Dar=D;Gss=ss(Abar,Bbar,Cbar,Dbar)转换成能观标准II型代码A=;B=;C=;D=;Gs=ss(A,B,C,D);Tt=obsv (Gs.a,Gs.c);Ttt=flipud(Tt);cp=poly(Gs.a);n= length(Gs.a);Tea=eye(n)for i=2:n for j=1:(n-1) if ij Tea(i,j)=cp(i-(j-1);endendendTea=Tea;T= Tea *Ttt;Abar=inv(T)*A*T;Bar=inv(T)*B;Cbar=C*T,Dar=D;Gss=ss(Abar,Bbar,Cbar,Dbar)实验内容：1、已知线性系统(1) 判断其状态可控性、可观性和传递函数的关系，并加以说明分析。(2) 对系统分别按能控性分解、能观性分解以及能控能观性分解。2、在Matlab中建立并运行如下的.m代码，回答下面的问题。num=1 2 3;den=conv(1 6 25,1 12 35);G=tf(num,den)Gs=ss(G)V J=jordan(Gs.a) %求特征向量和Gs.a的Jordan标准型Gss=ss2ss(Gs,inv(V) %Jordan型系统Gsm=canon(Gs,model) %模态型系统Gsf=canon(Gs,companion) %能控标准II型系统(1)给出无分号行的运行结果，并比较几个状态方程。(2)在什么情况下,cannon得到的是对角型系统?请举例说明。(3)将原理中给出的能控标准与能观标准型转换代码写成子函数的形式，并通过调用你所编写的子函数将.m文件中给出的模型变换成能控标准I、II型和能观标准I、II型，并从结果说明能控与能观标准型间关系。要求上述工作的在一个MATLAB程序.m文件中完成。并给出实验结果和相关说明。实验报告要求简要的说明实验原理(少于100字)；简明扼要的总结实验内容(少于100字)；编制.mdl和.m文件，并给出运行结果。报告格式请按实验报告模板编写。1、(1)function str = pdctrb(A,B)A=-6 -0.625 0.75;8 0 0;0 2 0;B=1;0;0;Qc=ctrb(A,B);r=rank(Qc);l=length(A);if r=l str=系统状态是完全可控的!;else str=系统不是状态完全可控的！;end pdctrbans =系统状态是完全可控的!function str = pdctrbo(A,B)A=-6 -0.625 0.75;8 0 0;0 2 0;B=1;0;0;C=1 -0.25 0.0625;Co=ctrb(A,B);m=size(C,1);Qyc=C*Co,1;Tm=rank(Qyc);if m=Tm str=系统输出是完全可控的!;else str=系统不是输出完全可控的！;end pdctrboans =系统输出是完全可控的!function str = pdobsv(A,C)A=-6 -0.625 0.75;8 0 0;0 2 0;C=1 -0.25 0.0625;Qo=obsv(A,C);r=rank(Qo);l=size(A,1);if r=l str=系统状态是完全可观的!;else str=系统不是状态完全可观的！;end pdobsvans =系统不是状态完全可观的！与传递函数的关系：A=-6 -0.625 0.75;8 0 0;0 2 0;B=1;0;0;C=1 -0.25 0.0625;D=0;G=ss(A,B,C,D);Gs=tf(G)Gs = s2 - 2 s + 1 - s3 + 6 s2 + 5 s - 12 num=1 -2 1; den=1 6 5 -12;z,p,k=tf2zp(num,den);Gss=zpk(z,p,k)Gss = (s-1)2 - (s+4) (s+3) (s-1)系统能控能观的充要条件是传递函数中没有零极点对消现象。(2)按能控性分解： A=-6 -0.625 0.75;8 0 0;0 2 0;B=1;0;0;C=1 -0.25 0.0625; Abar ,Bbar Cbar,T,K=ctrbf(A,B,C)Abar = 0 -2.0000 0 0 0 8.0000 -0.7500 -0.6250 -6.0000Bbar = 0 0 1Cbar = -0.0625 -0.2500 1.0000T = 0 0 -1 0 1 0 1 0 0K = 1 1 1按能观性分解： Abar ,Bbar Cbar,T,K=obsvf(A,B,C)Abar = 1.0000 -2.2772 2.7134 -0.0000 0.3407 -5.5238 0.0000 2.6251 -7.3407Bbar = 0.0558 -0.2432 0.9684Cbar = 0.0000 -0.0000 1.0327T = 0.0558 0.4465 0.8930 -0.2432 -0.8614 0.4459 0.9684 -0.2421 0.0605K = 1 1 0按能控能观性分解：程序：function Gk,T,K =kalmdec(G)A=-6 -0.625 0.75;8 0 0;0 2 0;B=1;0;0;C=1 -0.25 0.0625;D=0;G=ss(A,B,C,D)D=G.d;A=G.a;B=G.b;C=G.c;Ac,Bc,Cc,Tc,Kc=ctrbf(A,B,C);nc=rank(ctrb(A,B);n=length(A);ic=n-nc+1:n;Ao1,Bo1,Co1,To1,Ko1=obsvf(Ac(ic,ic),Bc(ic),Cc(ic);%if nc kalmdecGk = a = x1 x2 x3 x1 1 -2.277 2.713 x2 -5.551e-16 0.3407 -5.524 x3 -1.11e-16 2.625 -7.341 b = u1 x1 0.05581 x2 -0.2432 x3 0.9684 c = x1 x2 x3 y1 8.327e-17 2.776e-17 1.033 d = u1 y1 0 Continuous-time state-space model.2、(1)num=1 2 3;den=conv(1 6 25,1 12 35);G=tf(num,den)Gs=ss(G)V J=jordan(Gs.a)Gss=ss2ss(Gs,inv(V)Gsm=canon(Gs,model)Gsf=canon(Gs,companion)运行Untitled6.m Untitled6G = s2 + 2 s + 3 - s4 + 18 s3 + 132 s2 + 510 s + 875 Continuous-time transfer function.Gs = a = x1 x2 x3 x4 x1 -18 -8.25 -3.984 -3.418 x2 16 0 0 0 x3 0 8 0 0 x4 0 0 2 0 b = u1 x1 0.25 x2 0 x3 0 x4 0 c = x1 x2 x3 x4 y1 0 0.25 0.0625 0.04688 d = u1 y1 0 Continuous-time state-space model.V = -0.4883 -1.3398 0.4570 - 0.1719i 0.4570 + 0.1719i 1.5625 3.0625 -0.4375 + 1.5000i -0.4375 - 1.5000i -2.5000 -3.5000 -1.5000 - 2.0000i -1.5000 + 2.0000i 1.0000 1.0000 1.0000 1.0000 J = -5.0000 0 0 0 0 -7.0000 0 0 0 0 -3.0000 - 4.0000i 0 0 0 0 -3.0000 + 4.0000iGss = a = x1 x2 x3 x1 -5+1.04e-15i -2.13e-14+2.19e-15i 2.43e-14+2.27e-15i x2 -1.78e-15-4.62e-16i -7-1.06e-15i -5.7e-15+2.51e-15i x3 -1.38e-15+1.78e-15i 4.01e-15+2.22e-15i -3-4i x4 -1.73e-15+4.44e-16i 7.54e-15+4.44e-16i -4.88e-15+1.11e-15i x4 x1 2.01e-14-4.83e-15i x2 -4.96e-15-2.82e-15i x3 0+1.78e-15i x4 -3+4i b = u1 x1 1.6-2.22e-16i x2 -1+1.67e-16i x3 -0.3-0.1i x4 -0.3+0.1i c = x1 x2 x3 x4 y1 0.281-1.77e-16i 0.594-2.78e-16i -0.156+0.25i -0.156-0.25i d = u1 y1 0 Continuous-time state-space model.Gsm = a = x1 x2 x3 x4 x1 -3 4 0 0 x2 -4 -3 0 0 x3 0 0 -7 0 x4 0 0 0 -5 b = u1 x1 0.9492 x2 -0.1931 x3 4.942 x4 5.042 c = x1 x2 x3 x4 y1 0.1699 0.09054 -0.1201 0.08925 d = u1 y1 0 Continuous-time state-space model.Gsf = a = x1 x2 x3 x4 x1 0 0 0 -875 x2 1 0 0 -510 x3 0 1 0 -132 x4 0 0 1 -18 b = u1 x1 1 x2 0 x3 0 x4 0 c = x1 x2 x3 x4 y1 0 1 -16 159 d = u1 y1 0 Continuous-time state-space model.(2) 当系统的 | sI - A |= 0 得到的特征多项式的解无重根时，则cannon得到的是对角型系统。例如： A=0 1 0;0 0 1;-6 -11 -6;B=0 ;0; 1;C=1,0,0;D=0;E=eig(A)sys=ss(A,B,C,D);Gt,P=canon(sys,model)E = -1.0000 -2.0000 -3.0000Gt = a = x1 x2 x3 x1 -3 0 0 x2 0 -2 0 x3 0 0 -1 b = u1 x1 -3.881 x2 -4.899 x3 1.436 c = x1 x2 x3 y1 -0.1288 0.2041 0.3482 d = u1 y1 0 Continuous-time state-space model.P = -7.7621 -11.6431 -3.8810 -14.6969 -19.5959 -4.89908.6168 7.1807 1.4361(3)能控标准II型程序：function Gss2 = NK2(A,B,C,D)d6;A=Gs.a;B=Gs.b;C=Gs.c;D=Gs.d;Gs2=ss(A,B,C,D);T=ctrb(Gs.a,Gs.b);Abar=inv(T)*A*T;Bbar=inv(T)*B;Cbar=C*T;Dbar=D;Gss2=ss(Abar,Bbar,Cbar,Dbar)end结果： NK2Gss2 = a = x1 x2 x3 x4 x1 0 0 0 -875 x2 1 0 0 -510 x3 0 1 0 -132 x4 0 0 1 -18 b = u1 x1 1 x2 0 x3 0 x4 0 c = x1 x2 x3 x4 y1 0 1 -16 159 d = u1 y1 0 Continuous-time state-space model.能控标准I型：程序：function Gss1 = NK1(A,B,C,D)d6;A=Gs.a;B=Gs.b;C=Gs.c;D=Gs.d;Gs1=ss(A,B,C,D);Tt=ctrb(A,B);Ttt=fliplr(Tt);cp=poly(A);n=length(A);Tea=eye(n);for i=2:n for j=1:(n-1) if ij Tea(i,j)=cp(i-(j-1); end endendT=Ttt*Tea;Abar=inv(T)*A*T;Bbar=inv(T)*B;Cbar=C*T;Dbar=D;Gss1=ss(Abar,Bbar,Cbar,Dbar) end结果： NK1Gss1 = a = x1 x2 x3 x4 x1 2.842e-14 1 0 0 x2 -2.842e-14 0 1 0 x3 -1.023e-12 0 0 1 x4 -875 -510 -132 -18 b = u1 x1 0 x2 0 x3 0 x4 1 c = x1 x2 x3 x4 y1 3 2 1 0 d = u1 y1 0 Continuous-time state-space model.能观I型：程序：function Gss3 = NG1(A,B,C,D)d6;A=Gs.a;B=Gs.b;C=Gs.c;D=Gs.d;Gs3=ss(A,B,C,D);Tinv=obsv(Gs.a,Gs.c);T=inv(Tinv);Abar=inv(T)*A*T;Bbar=inv(T)*B;Cbar=C*T;Dbar=D;Gss3=ss(Abar,Bbar,Cbar,Dbar)end结果： NG1Gss3 = a = x1 x2 x3 x4 x1 -1.176e-14 1 -4.677e-16 -3.431e-17 x2 1.99e-13 7.105e-14 1 7.772e-16 x3 -3.183e-12 -9.095e-13 -1.137e-13 1 x4 -875 -510 -132 -1