(1)線性系統(tǒng)結(jié)構(gòu)分析與分解及標準型

1、實驗三十七 線性系統(tǒng)結(jié)構(gòu)分析與分解及標準型實驗類型：驗證 難度系數(shù)：0.3 實驗性質(zhì)：必做 課內(nèi)學時：0 課外學時：2 分組人數(shù)：2開課方式：在課外完成在MATLAB平臺上完成實驗。實驗目的：掌握線性系統(tǒng)狀態(tài)空間標準型、解及其模型轉(zhuǎn)換。實驗設(shè)備與軟件：1、MATLAB數(shù)值分析軟件實驗原理：1、標準型變換、矩陣Jordan型變換陣、特征值(1)標準型變換 命令格式 csys=canon(sys,type)說明：type制定規(guī)范型的形式，包括兩種選項：model(模態(tài)規(guī)范型)、companion(伴隨規(guī)范型,友矩陣型,能控II型)。(2)矩陣的Jordan規(guī)范型 命令格式 V J=Jordan(A

2、)說明：V特征向量，J是Jordan型 (3)求矩陣特征值和特征向量 命令格式 V J=eig(A) cv= eig(A) 說明：V特征向量，J是Jordan型;cv是特征值列向量2、狀態(tài)模型的相似變換： 命令格式 sysb=ss2ss(sys,T) 說明：sys是狀態(tài)空間模型，T是非奇異變換陣的逆陣² 傳遞函數(shù)模型與狀態(tài)空間模型之間的相互轉(zhuǎn)換：命令格式 A,B,C,D=tf2ss(num,den)num,den=ss2tf(A,B,C,D,iu) 說明：tf2ss:傳遞函數(shù)à狀態(tài)空間; ss2tf狀態(tài)空間à傳遞函數(shù); iu是第iu個輸入有效² zpk模

3、型與狀態(tài)空間模型之間的相互轉(zhuǎn)換：命令格式：A,B,C,D=zp2ss(z,p,k)z,p,k=ss2zp(A,B,C,D,iu) 說明：zp2ss :zpk模型à狀態(tài)空間模型; ss2zp狀態(tài)空間模型àzpk模型3、線性定常系統(tǒng)的可控制與可觀性及結(jié)構(gòu)分解關(guān)于可控制與可觀性及結(jié)構(gòu)分解的理論內(nèi)容，請看教材學習。(1)可控性和可觀性-一般采用能控性/能觀性矩陣類別(適用于離散或連續(xù)的情況)狀態(tài)可控性和輸出可控性子函數(shù)如下：function str =pdctrb(A,B)Qc=ctrb(A,B);r=rank(Qc);l=length(A);if r=l str=系統(tǒng)是狀態(tài)完全可

4、控的!else str=系統(tǒng)不是狀態(tài)完全可控的!end實際上，rank(Qc)為系統(tǒng)的狀態(tài)可控性指數(shù)，即系統(tǒng)中可控的狀態(tài)的數(shù)目。function str=pdctrbo(A,B)%輸出可控性Co=ctrb(A,B);m=size(C,1);%返回行數(shù)Qyc=C*Co,D;Tm=rank(Qyc);if m=Tm str=系統(tǒng)輸出是完全可控的!else str=系統(tǒng)輸出不是完全可控的!end實際上，rank(Qyc) 為系統(tǒng)的輸出可控性指數(shù)，即系統(tǒng)中可控的輸出的數(shù)目。狀態(tài)可觀性判別子函數(shù)代碼如下：function str=pdobsv(A,C)Qo=obsv (A,C);r=rank(Qo);l

5、= size(A,1);if r=l str=系統(tǒng)是狀態(tài)完全可觀的!else str=系統(tǒng)不是狀態(tài)完全可觀的!end實際上，rank(Qo)為系統(tǒng)的狀態(tài)可觀性指數(shù)，即系統(tǒng)中可觀測的狀態(tài)的數(shù)目。有了上述子函數(shù)，在Matlab中可以直接調(diào)用這些子函數(shù)來判斷可控性和可觀性。(2)可控性和可觀性Gram矩陣可由下面的函數(shù)求得 W=gram(sys,type)-sys是系統(tǒng)的狀態(tài)空間模型，type可以是c或o。通過判斷W的正定性判定其可控性或可觀性。由于從W的數(shù)學表達式上看，Wc和Wo是對稱的半正定矩陣，它們分別滿足下面的Lyapunov方程(這種方程在穩(wěn)定性分析中還將提到)、所以系統(tǒng)必須穩(wěn)定才能得到G

6、ram矩陣。Wc矩陣中的值對應于輸入信號對相應狀態(tài)的貢獻：第i個元素越大，則說明輸入信號對第i個輸入狀態(tài)的貢獻越大。Wo矩陣中的值對應于每個狀態(tài)對系統(tǒng)輸出的貢獻：第i個元素越大，則說明第i狀態(tài)對系統(tǒng)輸出的貢獻越大。(3)結(jié)構(gòu)分解a.在matlab中調(diào)用ctrbf()函數(shù)對系統(tǒng)按能控性分解Abar,Bbar,Cbar,T,K=ctrbf(A,B,C)Abar,Bbar,Cbar,T,K=ctrbf(A,B,C,TOL)說明：K是可控的狀態(tài)個數(shù)，TOL為誤差，這里的變換是這樣令的：x=T-1xbar，所以有Abar=TAT-1，Bbar=TB，Cbar=CT-1。顯然這與教材中的令法是不同的，且T

7、的選取方式也不同，得到的分解形式也不同，要加以區(qū)別。這里的能觀性分解形式為：，b.在matlab中調(diào)用ctrbf()函數(shù)對系統(tǒng)按能觀性分解Abar,Bbar,Cbar,T,K=obsvf(A,B,C)Abar,Bbar,Cbar,T,K= obsvf(A,B,C,TOL)說明：K是可觀的狀態(tài)個數(shù)，TOL為誤差，這里的變換是這樣令的：x=T-1xbar，所以有Abar=TAT-1，Bbar=TB，Cbar=CT-1。顯然這與教材中的令法是不同的，且T的選取方式也不同，得到的分解形式也不同，要加以區(qū)別。這里的能觀性分解形式為：，c.按能控、能觀性分解-Kalman分解Kalman分解是先按能控性進

8、行分解，然后對兩部分狀態(tài)再按能觀性進行分解，最后得到能控能觀部分、能控不能觀部分、不能控能觀部分、不能控不能觀部分。在Matlab中沒有提供對應的函數(shù)對系統(tǒng)進行直接的Kalman分解，這里我們給出一個函數(shù)實現(xiàn)分解：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 nc<ninc=1:n-nc;A

9、o2,Bo2,Co2,To2,Ko2=obsvf(Ac(inc,inc),Bc(inc),Cc(inc);endm1,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

10、,Ck,G.d);此函數(shù)分解的形式與教材中的是順序不一樣(變換的令法也不一樣)：這里是不能控不能觀、不能控能觀、能控不能觀、能控能觀。這一點要注意。4、定常線性系統(tǒng)的標準型(轉(zhuǎn)換限于SISO系統(tǒng))(1)Jordan標準型化成Jordan標準形是一種并聯(lián)分解的策略，即將傳遞函數(shù)展開成部分分式和的形式。這種分解的原理在課堂上已講的很清楚，請復習相關(guān)內(nèi)容。這里就如何用Matlab求取此標準型做說明。a.互異根的情況下，代碼可以采用如下的形式，出可以有采用Jordan指令：num=;den=;G=tf num,den;r,p,k=residue(num,den);A=diag(p);B=ones(le

11、ngth(r),1);C=r;%為了得到更準確的結(jié)果，可以采用C=(rat(r)得到有理分式的形式D=k;Gss=SS(A,B,C,D)b.有重根的情況 我們直接可以采用Jordan指令。num=;den=;Gtf=tf(num,den)Gs=ss(G)V J=jordan(Gs.a) %V是特征向量，由特征向量組成非奇異變換矩陣, x=Tz,T=VGss=ss2ss(Gs,inv(V) %Jordan型系統(tǒng) (2)能控標準型-限于SISO系統(tǒng)若系統(tǒng)能控，則可轉(zhuǎn)換成能控標準I型和II型。轉(zhuǎn)換成能控標準II型代碼：A=;B=;C=;D=;Gs=ss(A,B,C,D);T=ctrb(Gs.a,Gs

12、.b) % x =Tz，T=B,AB,A2B, A n-1B， 能控標準II型Abar=inv(T)*A*T;Bar=inv(T)*B;Cbar=C*T,Dar=D;Gss=ss(Abar,Bbar,Cbar,Dbar)轉(zhuǎn)換成能控標準I型代碼：A=;B=;C=;D=;Gs=ss(A,B,C,D);Tt=ctrb(Gs.a,Gs.b);Ttt=fliplr(Tt);cp=poly(Gs.a);n= length(Gs.a);Tea=eye(n)for i=2:n for j=1:(n-1) if i>j Tea(i,j)=cp(i-(j-1);endendendT=Ttt*Tea;Abar

13、=inv(T)*A*T;Bar=inv(T)*B;Cbar=C*T,Dar=D;Gss=ss(Abar,Bbar,Cbar,Dbar)(3)能觀標準型-限于SISO系統(tǒng)若系統(tǒng)能觀，則可轉(zhuǎn)換成能觀標準I型和II型。轉(zhuǎn)換成能觀標準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)轉(zhuǎn)換成能觀標準II型代碼A=;B=;C=;D=;Gs=ss(A,B,C,D);Tt=obsv (Gs

14、.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 i>j 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)實驗內(nèi)容：1、已知線性系統(tǒng)(1) 判斷其狀態(tài)可控性、可觀性和傳遞函數(shù)的關(guān)系，并加以說明分析。(2) 對系統(tǒng)分別按能控性分解、能觀性分解以及能控能觀性分解。2、在Ma

15、tlab中建立并運行如下的.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型系統(tǒng)Gsm=canon(Gs,'model') %模態(tài)型系統(tǒng)Gsf=canon(Gs,'companion') %能控標準II型系統(tǒng)(1)給出無分號行的運行結(jié)果，并比較幾個狀態(tài)方程。(2)在什么情況下,cannon得到的是對角型系統(tǒng)?請舉例說明。(3)將原理中給出的能控標

16、準與能觀標準型轉(zhuǎn)換代碼寫成子函數(shù)的形式，并通過調(diào)用你所編寫的子函數(shù)將.m文件中給出的模型變換成能控標準I、II型和能觀標準I、II型，并從結(jié)果說明能控與能觀標準型間關(guān)系。要求上述工作的在一個MATLAB程序.m文件中完成。并給出實驗結(jié)果和相關(guān)說明。實驗報告要求簡要的說明實驗原理(少于100字)；簡明扼要的總結(jié)實驗內(nèi)容(少于100字)；編制.mdl和.m文件，并給出運行結(jié)果。報告格式請按實驗報告模板編寫。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=len

17、gth(A);if r=l str='系統(tǒng)狀態(tài)是完全可控的!'else str='系統(tǒng)不是狀態(tài)完全可控的！'end>> pdctrbans =系統(tǒng)狀態(tài)是完全可控的!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='系統(tǒng)輸出是完全可控的!'else str='系統(tǒng)不是輸出完全可控的！'en

18、d>> pdctrboans =系統(tǒng)輸出是完全可控的!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='系統(tǒng)狀態(tài)是完全可觀的!'else str='系統(tǒng)不是狀態(tài)完全可觀的！'end>> pdobsvans =系統(tǒng)不是狀態(tài)完全可觀的！與傳遞函數(shù)的關(guān)系：A=-6 -0.625 0.75;8 0 0;0 2 0;B=1;0;0;C=1 -0.25 0.062

19、5;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)系統(tǒng)能控能觀的充要條件是傳遞函數(shù)中沒有零極點對消現(xiàn)象。(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 Cb

20、ar,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 =

21、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

22、:n;Ao1,Bo1,Co1,To1,Ko1=obsvf(Ac(ic,ic),Bc(ic),Cc(ic);%if nc<n inc=1:n-nc; Ao2,Bo2,Co2,To2,Ko2=obsvf(Ac(inc,inc),Bc(inc),Cc(inc);%endm1,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

23、,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);end結(jié)果：>> 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 Continuo

24、us-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 functio

25、n.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.5

26、000i -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.3

27、8e-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

28、.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 x

29、1 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) 當系統(tǒng)的 | sI - A |= 0 得到的特征多項式的解無重根時，則cannon得到的是對角型系統(tǒng)。例如：>> 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=

30、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)能控

31、標準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結(jié)果：>> 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

32、= 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 i>j Tea(i,j)=cp(i-(j-1); end endendT=Ttt*Tea;Aba

33、r=inv(T)*A*T;Bbar=inv(T)*B;Cbar=C*T;Dbar=D;Gss1=ss(Abar,Bbar,Cbar,Dbar) end結(jié)果：>> 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結(jié)果：>