蘭州理工大學(xué)線性系統(tǒng)理論期末MATLAB大作業(yè)

1、線性系統(tǒng)理論的Matlab實踐1、在造紙流程中，投料箱應(yīng)該把紙漿流變成2cm的射流，并均勻噴灑在網(wǎng)狀傳送帶上。為此，要精確控制噴射速度和傳送速度之間的比例關(guān)系。投料箱內(nèi)的壓力是需要控制的主要變量，它決定了紙漿的噴射速度。投料箱內(nèi)的總壓力是紙漿液壓和另外灌注的氣壓之和。由壓力控制的投料箱是個耦合系統(tǒng)，因此，我們很難用手工方法保證紙張的質(zhì)量。在特定的工作點上，將投料箱線性化，可以得到下面的狀態(tài)空間模型： x = -0.8+0.02-0.020 x+0.0510.0010 u y =x1 , x2 其中，系統(tǒng)的狀態(tài)變量x1=液面高度，x2=壓力，系統(tǒng)的控制變量u1=紙漿流量u2=氣壓閥門的開啟量。在

2、上述條件下，試設(shè)計合適的狀態(tài)變量反饋控制器，使系統(tǒng)具有實特征根，且有一個根大于5解：下面是對此設(shè)計的MATLAB程序?qū)崿F(xiàn)：>> A=-0.8 0.02;-0.02 0;>> B=0.05 1;0.001,0;>> r=rank(ctrb(A,B)r =2>> C=1 1;>> P=1 6;>> K=place(A,B,P)K = 1.0e+003 * -0.0200 -6.0000 -0.0008 0.30002、描述恒速制導(dǎo)導(dǎo)彈的運動方程為： x = 01000-0.1-0.50000.500000 010000.510

3、00 x + 01000 u y= 0 0 0 1 0 x (a) 運用ctrb函數(shù)計算系統(tǒng)的能控型矩陣，并驗證系統(tǒng)是不可控的；(b) 計算從u到Y(jié)的傳遞函數(shù)，并消去傳遞函數(shù)中的分子和分母公因式，由此可以得到能控的狀態(tài)空間模型。在消去了公因子之后，請用tf2ss函數(shù)確定新的狀態(tài)變量模型；(c) 證明(b)中得到的狀態(tài)變量模型是能控的；(d) 說明恒速制導(dǎo)導(dǎo)彈是否穩(wěn)定？(e) 討論狀態(tài)變量模型的能控性和復(fù)雜性的關(guān)系（假設(shè)用狀態(tài)變量的數(shù)目來度量復(fù)雜性）解程序如下:clearA=input('請輸入系統(tǒng)矩陣:');B=input('請輸入輸入矩陣:');C=input

4、('請輸入輸出矩陣:');Qc1=ctrb(A,B)N1=size(A);n1=N1(1) %判斷狀態(tài)方程維數(shù)rc1=rank(Qc1)if rc1=n1 disp('系統(tǒng)可控')elseif rc1<n1 disp('系統(tǒng)不可控')endsyms s I=eye(n1);Q=inv(s*I-A);sys=collect(C*Q*B) %求解原狀態(tài)方程的頻域傳遞函數(shù)并化簡num=500 250 50;den=1 0 0;A1 B1 C1 D1=tf2ss(num,den)Qc2=ctrb(A1,B1)N2=size(A1);n2=N2(1)

5、 %判斷狀態(tài)方程維數(shù)rc2=rank(Qc2)if rc2=n2 disp('系統(tǒng)可控')elseif rc2<n2 disp('系統(tǒng)不可控')endd1=eig(A)'d2=eig(A1)'flag1=0;flag2=0;for i=1:n1 if real(d1(i)>0 flag1=1; endendif flag1=1 disp('原系統(tǒng)不穩(wěn)定')else disp('原系統(tǒng)穩(wěn)定')endfor j=1:n2 if real(d2(j)>0 flag2=1; endendif flag2=

6、1 disp('新系統(tǒng)不穩(wěn)定')else disp('新系統(tǒng)穩(wěn)定')end運行結(jié)果：請輸入系統(tǒng)矩陣:0 1 0 0 0;-0.1 -0.5 0 0 0;0.5 0 0 0 0;0 0 10 0 0;0.5 1 0 0 0請輸入輸入矩陣:0;1;0;0;0請輸入輸出矩陣:0 0 0 1 0Qc1 = 0 1.0000 -0.5000 0.1500 -0.0250 1.0000 -0.5000 0.1500 -0.0250 -0.0025 0 0 0.5000 -0.2500 0.0750 0 0 0 5.0000 -2.5000 0 1.0000 0 -0.100

7、0 0.0500n1 = 5rc1 = 4系統(tǒng)不可控 sys = 50/s2/(10*s2+5*s+1) A1 = 0 0 1 0B1 = 1 0C1 = 250 50D1 = 500Qc2 = 1 0 0 1n2 = 2rc2 = 2系統(tǒng)可控d1 = 0 0 0 -0.2500 - 0.1936i -0.2500 + 0.1936id2 = 0 0原系統(tǒng)穩(wěn)定新系統(tǒng)穩(wěn)定分析：由上述分析結(jié)果可知原系統(tǒng)和新系統(tǒng)均穩(wěn)定，而實際上由系統(tǒng)的極點可知，原系統(tǒng)是穩(wěn)定的，新系統(tǒng)實際上處于臨界穩(wěn)定狀態(tài)也可認(rèn)為是不穩(wěn)定的；若以狀態(tài)變量的數(shù)目來度量復(fù)雜性，可知系統(tǒng)的完全可控性與復(fù)雜性存在類似反比的關(guān)系，及復(fù)雜性越高

8、系統(tǒng)完全可控的難度越大，復(fù)雜性越低系統(tǒng)完全可控的難度越低。3、垂直起降的飛機的線性化模型為：x=Ax+B1u1+B2u2其中A = -0.03660.02710.0188-0.04550.0482-0.01000.0214-4.02080.10020.3681-0.70701.42000000B1 = 0.44223.5466-5.52000 ， B2 = 0.1761-7.59224.49000系統(tǒng)的狀態(tài)變量為水平速度（節(jié)）、垂直速度（節(jié)）、傾斜率（度/秒）和傾斜角（度）；系統(tǒng)的控制輸入為和，其中用于控制垂直運動，用于控制水平運動。(a) 計算系統(tǒng)矩陣的特征值，并由此判斷系統(tǒng)是否穩(wěn)定；(b)

9、 利用poly函數(shù)確定的特征多項式，計算特征根，并與(a)中得到的特征根相比較；(c) 當(dāng)只有發(fā)揮作用時，系統(tǒng)能控嗎？當(dāng)只有發(fā)揮作用時，結(jié)果又如何？請比較解釋你的結(jié)論。解：矩陣A的特征值由下列方式實現(xiàn)：>> A=-0.0366 0.0271 0.0188 -0.4555;0.0482 -1.0100 0.0024 -4.0208;0.1002 0.3681 -0.7070 1.4200;0 0 1 0;>> Lambda=eig(A)Lambda = 0.2758 + 0.2576i 0.2758 - 0.2576i -0.2325 -2.0727 由上可知，系統(tǒng)有兩個

10、負(fù)實根-0.2325和-2.0727，兩個實部為正的共軛復(fù)根0.2758 + 0.2576i和0.2758 - 0.2576i，而要使系統(tǒng)漸進穩(wěn)定，所有的特征根必須都具有負(fù)實部，所以此系統(tǒng)為不穩(wěn)定系統(tǒng)。由poly函數(shù)確定系統(tǒng)特征多項式的實現(xiàn)如下：>> A=-0.0366 0.0271 0.0188 -0.4555;0.0482 -1.0100 0.0024 -4.0208;0.1002 0.3681 -0.7070 1.4200;0 0 1 0;>> b=poly(A)b = 1.0000 1.7536 -0.6472 0.0625 0.0686上述b的值為系統(tǒng)特征多項

11、式的系數(shù)，則系統(tǒng)特征多項式為a(s)= s4+1.7536s3-0.6472s2+0.0625s+0.0686，計算此特征多項式的根有如下實現(xiàn)過程：>> b=1 1.7536 -0.6472 0.0625 0.0686;>> roots(b)ans = -2.0727 0.2758 + 0.2576i 0.2758 - 0.2576i -0.2324 由特征多項式所得特征根為兩個負(fù)實根 -2.0727和 -0.2324，兩個具有正實部的共軛復(fù)根0.2758 + 0.2576i和0.2758 - 0.2576i，將其與(a)中所得特征根比較如下：(a): -2.0727

12、-0.2325 0.2758 + 0.2576i 0.2758 - 0.2576i (b): -2.0727 -0.2324 0.2758 + 0.2576i 0.2758 - 0.2576i 可以看出，(a)和(b)所得的系統(tǒng)特征值只有一個負(fù)實根不相同之外其他的特征根都相同，而不相同的兩個負(fù)實根-0.2325和-0.2324只相差0.0001，相差甚微，僅僅是這么小的差距對分析系統(tǒng)性能并沒有很大的影響，完全可以忽略。(c)當(dāng)只有u1作用于系統(tǒng)時，即輸入矩陣只有B1矩陣，由系統(tǒng)矩陣A和輸入矩陣B1確定的能控性判別矩陣可知系統(tǒng)的能控性，對此有如下實現(xiàn)過程：>> A=-0.0366 0

13、.0271 0.0188 -0.4555;0.0482 -1.0100 0.0024 -4.0208;0.1002 0.3681 -0.7070 1.4200;0 0 1 0;>> B1=0.4422;3.5466;-5.5200;0;>> Co1=ctrb(A,B1)Co1 = 0.4422 -0.0238 2.5171 -2.0270 3.5466 -3.5740 25.8160 -47.1028 -5.5200 5.2525 -12.8699 26.3126 0 -5.5200 5.2525 -12.8699>> r=rank(Co1)r =4系統(tǒng)狀態(tài)

14、維數(shù)n=4，又有r=n=4，所以只有u1單獨作用時系統(tǒng)完全能控。當(dāng)只有u2作用于系統(tǒng)時，相應(yīng)的輸入矩陣只有B2，可以由系統(tǒng)矩陣A和輸入矩陣B2確定的能控性判別矩陣判斷系統(tǒng)的能控性，具體的實現(xiàn)如下過程：>> B2=0.1761;-7.5922;4.4900;0;>> Co2=ctrb(A,B2)Co2 = 0.1761 -0.1278 -1.9441 2.3338 -7.5922 7.6874 -25.8381 49.9646 4.4900 -5.9515 13.4004 -27.6310 0 4.4900 -5.9515 13.4004>> r=rank(C

15、o2)r =4同樣有r=n=4，所以只有u2單獨作用于系統(tǒng)時系統(tǒng)也是完全可控的。4、為了探究月球背面（遠離地球的一面）的奧秘，人們付出了不懈的努力。例如，在地球-太陽-月球系統(tǒng)中，人們希望通信衛(wèi)星能定點在不受月球遮擋的軌道上，并為此開展了廣泛的論證研究工作。圖中給出了預(yù)期衛(wèi)星軌道的示意圖，從地球上看上去，衛(wèi)星軌道的光影恰似環(huán)繞月球的外層光暈，因此這種軌道又稱為光暈軌道。軌道控制的目的是，使通信衛(wèi)星在地球可見的光暈軌道上運行，從而保證通信鏈路的暢通，所需的通信鏈路包括從地球到衛(wèi)星和從衛(wèi)星到月球背面共兩段線路。衛(wèi)星繞定點位置運動時，經(jīng)過標(biāo)準(zhǔn)化和線性化的漂移運動方程為：其中，狀態(tài)變量是衛(wèi)星在三個方向

16、上的位置和速度漂移，輸入分別是軌控發(fā)動機在、和方向上產(chǎn)生的加速度。(a) 衛(wèi)星的定點位置是否穩(wěn)定？(b) 如果只有發(fā)揮作用，衛(wèi)星是否能控？(c) 如果只有發(fā)揮作用，衛(wèi)星是否能控？(d) 如果只有發(fā)揮作用，衛(wèi)星是否能控？(e) 如果能夠測得方向的位置漂移，請確定由到該位置漂移量的傳遞函數(shù)。（提示：可以令觀測輸出為）(f) 用tf2ss函數(shù)，計算(e)中得到的傳遞函數(shù)的狀態(tài)變量模型，并驗證該軌跡子系統(tǒng)是能控系統(tǒng)；(g) 采用狀態(tài)反饋，設(shè)計合適的反饋控制器，使(f)中得到的系統(tǒng)的閉環(huán)極點為和。解：(a)由系統(tǒng)矩陣A的特征值可以判斷系統(tǒng)的穩(wěn)定性，如下可知系統(tǒng)含有正實數(shù)特征根，故系統(tǒng)是不穩(wěn)定的>&

17、gt; A=0 0 0 1 0 0;0 0 0 0 1 0;0 0 0 0 0 1;7.3809 0 0 0 2 0;0 -2.1904 0 -2 0 0;0 0 -3.1904 0 0 0>> Lambda=eig(A)Lambda = -2.1587 2.1587 0.0000 + 1.8626i 0.0000 - 1.8626i 0 + 1.7862i 0 - 1.7862i(b) 系統(tǒng)只有u1發(fā)揮作用時，利用系統(tǒng)矩陣A和輸入矩陣B1得到的能控性判別矩陣Co1判斷系統(tǒng)的能控性，有如下實現(xiàn)過程：>> A=0 0 0 1 0 0;0 0 0 0 1 0;0 0 0 0

18、 0 1;7.3809 0 0 0 2 0;0 -2.1904 0 -2 0 0;0 0 -3.1904 0 0 0;>> B1=0;0;0;1;0;0;>> Co1=ctrb(A,B1)Co1 = 0 1.0000 0 3.3809 0 20.1921 0 0 -2.0000 0 -2.3810 0 0 0 0 0 0 0 1.0000 0 3.3809 0 20.1921 0 0 -2.0000 0 -2.3810 0 -35.1688 0 0 0 0 0 0>> r1=rank(Co1)r1 =4系統(tǒng)狀態(tài)變量的維數(shù)n=5，由于r1=4<n=5，所

19、以只有u1發(fā)揮作用時，衛(wèi)星定點系統(tǒng)不完全能控。(c)當(dāng)系統(tǒng)只有u2發(fā)揮作用時，由系統(tǒng)矩陣A和輸入矩陣B2構(gòu)成能控性判別矩陣Co2，判斷系統(tǒng)能控性的實現(xiàn)過程如下：>> A=0 0 0 1 0 0;0 0 0 0 1 0;0 0 0 0 0 1;7.3809 0 0 0 2 0;0 -2.1904 0 -2 0 0;0 0 -3.1904 0 0 0;>> B2=0;0;0;0;1;0;>> Co2=ctrb(A,B2)Co2 = 0 0 2.0000 0 2.3810 0 0 1.0000 0 -6.1904 0 8.7975 0 0 0 0 0 0 0 2.

20、0000 0 2.3810 0 35.1688 1.0000 0 -6.1904 0 8.7975 0 0 0 0 0 0 0>> r2=rank(Co2) r2=4同樣有r2=4<n=6，故只有u2作用時，衛(wèi)星定點系統(tǒng)也是不完全能控的。(d)當(dāng)系統(tǒng)只有u3發(fā)揮作用時，由系統(tǒng)矩陣A和輸入矩陣B3構(gòu)成能控性判別矩陣Co3，判斷系統(tǒng)能控性的實現(xiàn)過程如下：>> A=0 0 0 1 0 0;0 0 0 0 1 0;0 0 0 0 0 1;7.3809 0 0 0 2 0;0 -2.1904 0 -2 0 0;0 0 -3.1904 0 0 0;>> B3=0;

21、0;0;0;0;1;>> Co3=ctrb(A,B3)Co3 = 0 0 0 0 0 0 0 0 0 0 0 0 0 1.0000 0 -3.1904 0 10.1787 0 0 0 0 0 0 0 0 0 0 0 0 1.0000 0 -3.1904 0 10.1787 0>> r3=rank(Co3)r3 =2可知r3=2<n=6,所以當(dāng)只有u3作用于系統(tǒng)時，衛(wèi)星定點系統(tǒng)也是不可控的。(e) 以U2為輸入到該位置漂移量求得系統(tǒng)的傳遞函數(shù)MATLAB實現(xiàn)如下：>> syms s>> A=0 0 0 1 0 0;0 0 0 0 1 0;0

22、0 0 0 0 1;7.3809 0 0 0 2 0;0 -2.1904 0 -2 0 0;0 0 -3.1904 0 0 0；>> B=0;0;0;0;0;1;>> C=0 1 0 0 0 0;>> I=eye(6)；>> F=inv(s*I-A)；>> G=simple(simple(C*F*B)；最后得到系統(tǒng)傳遞函數(shù)為：G =-(6250000*s2 - 46130625)/(7440625*s2 6250000*s4 + 101044521)(f)由上述所得傳遞函數(shù)獲取狀態(tài)變量模型有如下實現(xiàn)：>> num=-625

23、0000 0 46130625;>> den=-6250000 0 7440625 0 101044521;>> G=tf(num,den)>> sys=ss(G)a = x1 x2 x3 x4x1 0 0.5952 0 2.021x2 2 0 0 0x3 0 2 0 0x4 0 0 2 0b = u1 x1 1 x2 0 x3 0 x4 0c = x1 x2 x3 x4 y1 0 0.5 0 -0.9226d = u1 y1 0 上述系統(tǒng)的能控性判別矩陣及其秩如下：QC = 1.0000 0 1.1904 0 0 2.0000 0 2.3808 0 0 4

24、.0000 0 0 0 0 8.0000>> r=rank(QC)r = 4r=4=n故系統(tǒng)是完全能控的(g)由于是單輸入系統(tǒng)，故在極點配置的過程中對于重極點的情形會在實現(xiàn)時程序出錯，故將給定的兩個重極點-10中的一個取為與其很接近的-10.0001，這樣的取法對系統(tǒng)的性能影響甚微，故具備可行性，對此的極點配置的狀態(tài)反饋矩陣有如下實現(xiàn)：>> A=0 0.5952 0 2.021;2 0 0 0;0 2 0 0;0 0 2 0;>> B=1;0;0;0;>> C=0 0.5 0 -0.9266;>> P=-1+j -1-j -10 -1

25、0.0001;>> K=place(A,B,P)K = 22.0001 71.5958 60.0005 27.0212則有u2=Kx所得控制輸入為：u2 = 22.0001 71.5958 60.0005 27.0212x5、在8.2風(fēng)力機的一階模型中，采用漿距角控制風(fēng)力機的轉(zhuǎn)速，風(fēng)速的變化視為擾動，設(shè)計風(fēng)力機轉(zhuǎn)速的閉環(huán)PI控制，使轉(zhuǎn)速恒定。解：要設(shè)計風(fēng)力機轉(zhuǎn)速的閉環(huán)PI控制，首先繪制出帶有PI控制器的風(fēng)力機控制系統(tǒng)，如圖所示：帶有PI控制器的風(fēng)力機控制器系統(tǒng)框圖圖中PI控制器的模型為K1+K2/s，從而可以求出該系統(tǒng)閉環(huán)系統(tǒng)的特征方程為S2+0.3397K1*S+0.3397K2

26、=0要使得該風(fēng)力機穩(wěn)定運行，則需讓特征方程的根具有負(fù)實部。設(shè)計程序如下：syms s k1 k2;>> s=solve(s2+0.339*k1*s+0.339*k2) s = - (339*k1)/2000 - (114921*k12)/1000000 - (339*k2)/250)(1/2)/2 (114921*k12)/1000000 - (339*k2)/250)(1/2)/2 - (339*k1)/2000使其為負(fù)數(shù)，從而可得出K1<=180，K2<=86,得到的PI控制器就能滿足要求。6、在8.2風(fēng)力機的三階模型中，采用漿距角控制風(fēng)力機的轉(zhuǎn)速，風(fēng)速的變化視為擾

27、動，電磁轉(zhuǎn)矩視為常數(shù)，采用狀態(tài)反饋和極點配置算法，設(shè)計風(fēng)力機轉(zhuǎn)速的閉環(huán)控制系統(tǒng)。解：即確定狀態(tài)反饋矩陣K，從而使得系統(tǒng)矩陣的特征根為期望特征根，不妨設(shè)期望的系統(tǒng)矩陣特征根為-2，-3，-4；設(shè)計程序如下：A=-0.8468 -8.1598e-008 0.5071;867637000 0 -867637000;1.1636e+004 0.0019 -1.1636e+004;>> B=-0.0262;0;0;>> Qc=ctrb(A,B)Qc = 1.0e+011 * -0.0000 0.0000 -0.0000 0 -0.0002 2.6453 0 -0.0000 0.0

28、000>> Rc=rank(Qc)Rc = 3>> P=-2 -3 -4;>> K=place(A,B,P)K = 1.0e+005 *4.4381 0.0000 -4.4380按照所求K所構(gòu)成的狀態(tài)反饋矩陣從而形成的狀態(tài)反饋能夠滿足要求。7、問題同上，設(shè)計風(fēng)力機轉(zhuǎn)速的LQR控制器。解：以漿距角為控制輸入的三階模型如下：=x+v要對系統(tǒng)進行狀態(tài)反饋極點配置需先判斷系統(tǒng)能控性，能控性判別矩陣為：Qc = 1.0e+011 * -0.0000 0.0000 -0.0000 0 -0.0002 2.6453 0 -0.0000 0.0000可以看出rankQc=3

29、，系統(tǒng)完全能控，則可以進行任意極點的配置。如下計算系統(tǒng)的特征值：>> A=-0.8468 -8.1598e-008 0.5071;867637000 0 -867637000;1.1636e+004 0.0019 -1.1636e+004;>> Lambda=eig(A)Lambda = 1.0e+004 * -1.1493 -0.0000 -0.0143為了滿足控制要求，為系統(tǒng)取定幾個位于左半S平面的幾個閉環(huán)極點，首先取離虛軸較近的一對共軛主導(dǎo)極點，其余的極點取到離虛軸的距離為主導(dǎo)極點距離虛軸3-6倍，則系統(tǒng)性能主要取決于主導(dǎo)極點，而受其他非主導(dǎo)極點的影響不大，基此取

30、定如下三個期望的系統(tǒng)閉環(huán)極點：S1,2=-2±4j，S3=-10對期望極點進行配置算法如下：A=-0.8468 -8.1598e-008 0.5071;867637000 0 -867637000;1.1636e+004 0.0019 -1.1636e+004;>> B=-0.0262;0;0;>> P=-2+4j -2-4j -10;>> K=place(A,B,P)K = 1.0e+005 * 4.4362 0.0000 -4.4361狀態(tài)反饋之后閉環(huán)系統(tǒng)的方框圖如下所示：X=AX+B1+B2V U -K8、風(fēng)力機三階模型中表示系統(tǒng)軸上的扭轉(zhuǎn)彈

31、力第二個狀態(tài)是不容易測量的狀態(tài)變量，設(shè)計關(guān)于該狀態(tài)的觀測器，并構(gòu)成狀態(tài)反饋系統(tǒng)。解：本題目要設(shè)計風(fēng)力機三階模型的狀態(tài)反饋及其狀態(tài)觀測器，并且為降維觀測器，由于輸出矩陣C的秩為1，故而降維觀測器的維數(shù)為3-1=2；不妨設(shè)狀態(tài)反饋的期望極點為-3，-4，-5，狀態(tài)觀測器的期望極點為-7，-8則相應(yīng)的設(shè)計程序如下：設(shè)計程序：>> A=-0.8468 -8.1598e-008 0.5071;867637000 0 -867637000;1.1636e+004 0.0019 -1.1636e+004;>> B=-0.0262;0;0;>> C=0 0 1;>> Qc=ctrb(A,B)Qc = 1.0e+011 * -0.0000 0.0000 -0.0000 0 -0.0002 2.6453 0 -0.0000 0.0000>> Rc=rank(Qc)Rc = 3>> Qo=obs