线性系统结构分析与分解及标准型

1广西大学实验报告纸实验内容：线性系统结构分析与分解及标准型【实验时间】2014年5月28日【实验地点】宿舍(课外)【实验目的】掌握线性系统状态空间标准型、解及其模型转换。【实验设备与软件】1、MATLAB数值分析软件【实验原理】1、标准型变换、矩阵Jordan型变换阵、特征值。其格式分别为csys=canon（sys，’type’）；[VJ]=Jordan（A）；[VJ]=eig（A）、cv=eig（A）。2、状态模型的相似变换：命令格式sysb=ss2ss（sys，T）、传递函数模型与状态空间模型之间的相互转换：命令格式[A,B,C,D]=tf2ss（num，den）[num,den]=ss2tf（A,B,C,D,iu）、zpk模型与状态空间模型之间的互相转换：命令格式[A,B,C,D]=zp2ss（z,p,k）[z,p,k]=ss2zp（A,B,C,D,iu）3、线性定常系统的可控制性与可观性及结构分解（1）可控性和可观性——一般采用能控性或能观性矩阵类别状态可观性判别子函数代码如下：Functionstr=pdobsv(A,C)Qo=obsv(A,C);r=rank(Qo);1=size(A,1);ifr==1Str=’系统是状态完全可观的！’elsestr=’系统不是状态完全可观的！’end（2）可控性和可观性Gram矩阵可由下面函数求得W=gram（sys，type）---sys是系统的状态空间模型，type可以是’c’和’o’。（3）结构分解a.在matlab中调用ctrbf（）函数对系统按能控性分解[Abar,Bbar,Cbar,T,K]=ctrbf(A,B,C)[Abar,Bbar,Cbar,T,K]=ctrbf(A,B,C,TOL)b.在matlab中调用ctrbf（）函数对系统按能观性分解[Abar,Bbar,Cbar,T,K]=obsvf(A,B,C)[Abar,Bbar,Cbar,T,K]=obsvf(A,B,C,TOL)c.按能控、能观性分解--Kalman分解4、定常线性系统的标准型（1）Jordan标准型2a.互异根的情况下，代码可以采用如下的形式，出可以有采用Jordan指令。num=[...];den=[...];G=tf[num,den];[r,p,k]=residue(num,den);A=diag(p);B=ones(length(r),1);C=k;Gss=SS(A,B,C,D)b.有重根的情况num=[...];den=[...];Gtf=tf(num,den)Gs=ss(G)[VJ]=jordan(Gs.a)Gss=ss2ss(Gs,inv(V))（2）能控标准型-限于SISO系统转换成能控标准Ⅱ型代码：A=[...];B=[...];C=[...];D=[...];Gs=ss(A,B,C,D);T=ctrb(Gs.a,Gs.b)Abar=inv(T)*A*T；Bar=inv(T)*B;Cbar=C*T,Dar=D;Gss=ss(Abar,Bbar,Cbar,Dbar)（3）能观标准型--限于SISO系统转换成能观标准Ⅰ型代码：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)【实验内容、方法、过程与分析】1.实验内容1.已知线性系统uxyux0625.025.0100102000875.0625.06x（1）判断其状态可控性、可观性和传递函数的关系，并加以说明分析。程序代码如下：子函数如下：functionstr=pdctrb(a,b)Qc=ctrb(a,b);r=rank(Qc);l=length(a);ifr==lstr='系统是状态完全可控的！';else3str='系统是状态不完全可控的！';endfunctionstr=pdobsv(a,c)Qo=obsv(a,c);r=rank(Qo);l=size(a,1);ifr==lstr='系统是状态完全可观的！';elsestr='系统是状态不完全可观的！';end主函数为：a=[-6-0.6250.75;800;020];b=[100]';c=[1-0.250.0625];d=1;G=ss(a,b,c,d);Gz=zpk(G)[y,t,x]=step(G,36);plot(t,x,'r',t,y,'b');gridon;axis([036010]);title('系统状态和单位阶跃输出曲线');xlabel('t(s)');ylabel('y(t)/x(t)');str=pdctrb(a,b)str=pdobsv(a,c)运行结果：Zero/pole/gain:(s+6.236)(s+1.764)(s-1)-------------------------(s+4)(s+3)(s-1)str=系统是状态完全可控的！str=系统是状态不完全可观的！说明分析：由运行结果表明，传递函数发生了s=1的零极点对消。因此，系统是不完全可观的。（2）对系统分别按能控性分解、能观性分解以及能控能观性分解。能控性分解：[abar,bbar,cbar,T,K]=ctrbf(a,b,c)4abar=02.0000000-8.0000-0.75000.6250-6.0000bbar=001cbar=-0.06250.25001.0000T=00-10-10100K=111能观性分解：[abar,bbar,cbar,T,K]=obsvf(a,b,c)abar=1.00002.2772-2.71340.00000.3407-5.5238-0.00002.6251-7.3407bbar=0.05580.2432-0.9684cbar=0.00000.0000-1.0327T=0.05580.44650.89300.24320.8614-0.4459-0.96840.2421-0.0605K=110能控能观性分解：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));ifncninc=1:n-nc;5[Ao2,Bo2,Co2,To2,Ko2]=obsvf(Ac(inc,inc),Bc(inc),Cc(inc));end[m1,n1]=size(To1);[m2,n2]=size(To2);To=[To2,zeros(m2,n1);zeros(m1,n2),To1];T=To*Tc;n1=rank(obsv(Ac(ic,ic),Cc(ic)));n2=rank(obsv(Ac(inc,inc),Cc(inc)));K=[zeros(1,n-nc-n2),ones(1,n2),...,2*ones(1,nc-n1),3*ones(1,n1)];Ak=T*A*inv(T);Bk=T*B;Ck=C*inv(T);Gk=ss(Ak,Bk,Ck,G.d);2.在MATLAB建立并运行如下的.m代码，回答下面的问题。num=[123];den=conv([1625],[11235]);G=tf(num,den)Gs=ss(G)[VJ]=jordan(Gs.a)Gss=ss2ss(Gs,inv(V))Gsm=canon(Gs,'model')Gsf=canon(Gs,'companion')运行结果：Transferfunction:s^2+2s+3------------------------------------s^4+18s^3+132s^2+510s+875a=x1x2x3x4x1-18-8.25-3.984-3.418x216000x30800x40020b=u1x10.25x20x30x40c=x1x2x3x4y100.250.06250.04688d=6u1y10Continuous-timemodel.V=-0.4883-1.33980.4570+0.1719i0.4570-0.1719i1.56253.0625-0.4375-1.5000i-0.4375+1.5000i-2.5000-3.5000-1.5000+2.0000i-1.5000-2.0000i1.00001.00001.00001.0000J=-5.00000000-7.00000000-3.0000+4.0000i0000-3.0000-4.0000ia=x1x2x3x4x1-5+8.4e-015i-1.07e-014+8.55e-015i1.9e-014+8.38e-015i2.19e-014+1.19e-014ix2-3.55e-015-4.3e-015i-7-4.43e-015i-5.7e-015-6.06e-015i-4.07e-015+1.57e-016ix34.97e-016-3.55e-015i1.81e-015-3.55e-015i-3+4i4.44e-016-4.22e-015ix43.91e-016-1.78e-015i4.41e-015-1.78e-015i-4.88e-015-3.11e-015i-3-4ib=u1x11.6+2.22e-016ix2-1-1.67e-016ix3-0.3+0.1ix4-0.3-0.1ic=x1x2x3x4y10.281+1.25e-016i0.594+1.94e-016i-0.156-0.25i-0.156+0.25id=u1y10Continuous-timemodel.a=x1x2x3x4x1-3400x2-4-300x300-70x4000-5b=u1x10.9492x2-0.1931x34.942x45.0427c=x1x2x3x4y10.16990.09054-0.12010.08925d=u1y10Continuous-timemodel.a=x1x2x3x4x1000-875x2100-510x3010-132x4001-18b=u1x11x20x30x40c=x1x2x3x4y101-16159d=u1y10（2）canon得到的对角型系统为：num=[21];den=[17148];G=tf(num,den);Gs=ss(G);Gsm=canon(Gs,'model')运行结果：a=x1x2x3x1-400x20-20x300-1b=u1x1-3.83x2-4.899x31.915c=x1x2x38y10.3046-0.3062-0.1741d=u1y10说明：传递函数中只含有各个相异的实数极点，可以转换为对角型系统。（3）将原理中给出的能控标准型与能观标准型转换代码写成子函数的形式，并通过调用你所编写的子函数将.m文件中给出的模型变换成能控标准Ⅰ、Ⅱ型和能观标准Ⅰ、Ⅱ型，并从结果说明能控与能观标准型间关系。因为系统可控不能观，故只能能转化成能控型。转换成能控标准Ⅱ型代码：A=[-6-0.6250.75;800;020];B=[100]';C=[1-0.250.0625];D=1;Gs=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;Gss=ss(Abar,Bbar,Cbar,Dbar)运行结果：T=1-63108-480016Cbar=1-844a=x1x2x3x10012x210-5x301-6b=u1x11x20x30c=x1x2x3y11-844d=u1y11转换成