1、實驗二 利用 MATLAB 進行時域分析本實驗內容包含以下三個部分：基于MATLAB的線性系統穩定性分析、基于MATLAB的線性系統動態性能分析、和MATALB進行控制系統時域分析的一些其它實例。一、 基于MATLAB的線性系統穩定性分析線性系統穩定的充要條件是系統的特征根均位于S平面的左半部分。系統的零極點模型可以直接被用來判斷系統的穩定性。另外，MATLAB語言中提供了有關多項式的操作函數，也可以用于系統的分析和計算。(1) 直接求特征多項式的根設p為特征多項式的系數向量，則MATLAB函數roots()可以直接求出方程 p=0在復數范圍內的解v,該函數的調用格式為：v=roots(p)例

2、 已知系統的特征多項式為：5 3 2x 3x 2x x 1特征方程的解可由下面的MATLAB命令得出。>> p=1,0,3,2,1,1;v=roots(p)結果顯示：v =+利用多項式求根函數roots(), 可以很方便的求出系統的零點和極點，然后根據零極點分析系統穩定性和其它性能。(2) 由根創建多項式如果已知多項式的因式分解式或特征根，可由MATLAB函數poly()直接得出特征多項式系數向量，其調用格式為：p=poly(v)如上例中：v=+; ;+;>> p=poly(v)結果顯示p =由此可見，函數roots()與函數poly()是互為逆運算的。(3) 多項式求

3、值在MATLAB中通過函數polyval()可以求得多項式在給定點的值，該函數的調用格式為:polyval(p,v)對于上例中的p值，求取多項式在x點的值，可輸入如下命令：>> p=1,0,3,2,1,1;x=1polyval(p,x)結果顯示x =1ans =(4) 部分分式展開考慮下列傳遞函數：M (s)num b°sn biSn 1bnN(s)den a0sn a1 sn 1an式中a。0，但是ai和bj中某些量可能為零。MATLAB函數可將器展開成部分分式，直接求出展開式中的留數'極點和余項。該函數的調用格式為:r,p,k residue( nu m,de

4、 n)則總的部分分式展開由下式給出:N(s)k(s) s P(n)M (s)r(1)r(2)応 s p(1) s p(2)式中 p(1)P1,p(2)P2,，P5)Pn，為極點,r(1)r1 ,r(2) Q,r(n)rn為各極點的留數，k(s)為余項。例設傳遞函數為:G(s)2s3 5s2 3s 6s3 6s211s 6該傳遞函數的部分分式展開由以下命令獲得:>> nu m=2,5,3,6;den=1,6,11,6;r,p,k=residue( nu m,de n)命令窗口中顯示如下結果中留數為列向量r，極點為列向量p,余項為行向量 k。由此可得出部分分式展開式:64G(s)32該

5、函數也可以逆向調用,把部分分式展開轉變回多項式器之比的形式，命令格式為:nu m,de n=residue(r,p,k)對上例有：>> nu m,de n=residue(r,p,k)結果顯示num =den =應當指出，如果p(j)=p(j+1)= -=p(j+m-1),則極點p(j)是一個m重極點。在這種情況下，部分分式展開式將包括下列諸項:r(j)s p( j)r(j 1)s p(j) 2r(j m 1)p(j)例設傳遞函數為:G(s)s2 2s 3(s 1)3s2 2s 3s3 3s2 3s 1則部分分式展開由以下命令獲得:>> v=-1,-1,-1num=0,

6、1,2,3;den=poly(v);r,p,k=residue( nu m,de n)結果顯示v =-1-1 -1其中由poly()命令將分母化為標準降幕排列多項式系數向量由上可得展開式為：den, k=為空矩陣。G(s)1 02s 1 (s 1)2(s 1)3(5) 由傳遞函數求零點和極點。在MATLAB控制系統工具箱中，給出了由傳遞函數對象G求出系統零點和極點的函數，其調用格式分別為：Z=tzero(G)P=1注意：式19中要求的G必須是零極點模型對象，且出現了矩陣的點運算“”和大括號訂表示的矩陣元素，詳細內容參閱后面章節。例已知傳遞函數為：2、6.8s61.2s95.2G(s) 3s47

7、.5s3 22s2 19.5s輸入如下命令：num=,;den=1,22,0;G=tf( nu m,de n);G1=zpk(G);Z=tzero(G)P=1結果顯示Z =-7-2(6) 零極點分布圖在MATLAB中，可利用pzmap()函數繪制連續系統的零、極點圖，從而分析系統的穩定性，該函數調用格式為：pzmap( nu m,de n)例給定傳遞函數：3s4 2s3 5s2 4s 6s5 3s4 4s3 2s2 7s 2利用下列命令可自動打開一個圖形窗口，顯示該系統的零、極點分布圖，如圖3-1所示。>> nu m=3,2,5,4,6;den=1,3,4,2,7,2;pzmap

8、(nu m,de n)title( 1 Pole-Zero Map1) % 圖形標題。1.5cfxa wanraay-1.5-2-ILX-一-rr-1.5-1-0.5Real Axis0Pole-Zero Map5 oo5 o-0.5系統動態特性分析。控制系統軟件包提供了控制系統工程需要的基本的時域與頻域分析工具函數。連續時間系統分析函數Impulse脈沖響應Step階躍響應Lsim任意輸入的仿真Bode波特圖Nyquist奈奎斯特圖Lyap李雅普諾夫方程Gram可控性與可觀性(1)時域響應解析算法部分分式展開法用拉氏變換法求系統的單位階躍響應，可直接得出輸出c(t)隨時間t變化的規律，對于高

9、階系統，輸出的拉氏變換象函數為：1num 1C(s) G(s)sden s對函數c(s)進行部分分式展開，我們可以用num,den,0來表示c(s)的分子和分母。例給定系統的傳遞函數：、s3 7s2 24s 24G(s) 3s410s335s250s 24用以下命令對9®進行部分分式展開。s>> nu m,de n,0num=1,7,24,24den=1,10,35,50,24r,p,k=residue( nu m,de n, 0)輸出結果為num =24241103550r =den =240k =輸出函數c(s)為：C( s) =num/ ( den*s )num =

10、1 72424ans =24 01103550拉氏變換得：L=laplace(sA3+7*s+24*s+24)/(sA4+10*sA3+35*sA2+50*s+24),s,t)(2 )單位階躍響應的求法:控制系統工具箱中給出了一個函數step()來直接求取線性系統的階躍響應，如果已知傳遞函數為:G(s)numden則該函數可有以下幾種調用格式:step (nu m,de n)step (nu m,de n,t)或(b)step(G)(c)step(G,t)(d)該函數將繪制出系統在單位階躍輸入條件下的動態響應圖，同時給出穩態值。對于式(b)和(d)，t為圖像顯示的時間長度，是用戶指定的時間向量

11、。式(a)和(c)的顯示時間由系統根據輸出曲線的形狀自行設定。如果需要將輸出結果返回到MATLAB工作空間中，則采用以下調用格式:c=step(G)此時，屏上不會顯示響應曲線，必須利用plot()命令去查看響應曲線。plot可以根據兩個或多個給定的向量繪制二維圖形，詳細介紹可以查閱后面的章節。例已知傳遞函數為:G(s)25s2 4s 251 .40.3050.4UllbSt即住訊皿梵 GT G(S>23CSss-25)UU.51.522.5I me sec0圖3-2 MATLAB繪制的響應曲線>>num=0,0,25;den=1,4,25;step( nu m,de n)gr

12、id %繪制網格線title( 1Unit-Step Response of G(s)=25(s2+4s+25) 1) % 圖像標題我們還可以用下面的語句來得出階躍響應曲線>> G=tf(0,0,25,1,4,25);t=0:5;%從0到5每隔取一個值。c=step(G,t); %動態響應的幅值賦給變量cplot(t,c) %繪二維圖形，橫坐標取t,縱坐標取C。Css=dcgai n(G) % 求取穩態值。系統顯示的圖形類似于上一個例子，在命令窗口中顯示了如下結果Css=1(3 )求階躍響應的性能指標MATLAB提供了強大的繪圖計算功能，可以用多種方法求取系統的動態響應指標。我們首

13、先介紹一種最簡單的方法一一游動鼠標法。對于例16,在程序運行完畢后，用鼠標左鍵 點擊時域響應圖線任意一點， 系統會自動跳出一個小方框， 小方框顯示了這一點的橫坐標 (時 間)和縱坐標(幅值)。按住鼠標左鍵在曲線上移動，可以找到曲線幅值最大的一點一一即曲線最大峰值， 此時小方框中顯示的時間就是此二階系統的峰值時間， 根據觀察到的穩態值 和峰值可以計算出系統的超調量。 系統的上升時間和穩態響應時間可以依此類推。 這種方法 簡單易用，但同時應注意它不適用于用plot()命令畫出的圖形。另一種比較常用的方法就是用編程方式求取時域響應的各項性能指標。與上一段介紹的游動鼠標法相比，編程方法稍微復雜，但通過

14、下面的學習，讀者可以掌握一定的編程技巧， 能夠將控制原理知識和編程方法相結合， 自己編寫一些程序， 獲取一些較為復雜的性能指標。通過前面的學習， 我們已經可以用階躍響應函數step( )獲得系統輸出量，若將輸出量返回到變量 y 中，可以調用如下格式y,t=step(G)該函數還同時返回了自動生成的時間變量 t ，對返回的這一對變量 y 和 t 的值進行計算， 可以得到時域性能指標。 峰值時間 (timetopeak) 可由以下命令獲得：Y,k=max(y);timetopeak=t(k)應用取最大值函數 max()求出y的峰值及相應的時間，并存于變量Y和k中。然后在變量 t 中取出峰值時間，并

15、將它賦給變量 timetopeak 。 最大 (百分比 )超調量 (percentovershoot) 可由以下命令得到：C=dcgain(G);Y,k=max(y);percentovershoot=100*(Y-C)/Cdcgain()函數用于求取系統的終值，將終值賦給變量C,然后依據超調量的定義，由Y 和C計算出百分比超調量。 上升時間(risetime)可利用MATLAB中控制語句編制 M文件來獲得。首先簡單介紹 一下循環語句 while 的使用。while 循環語句的一般格式為：while< 循環判斷語句 >循環體end其中 ,循環判斷語句為某種形式的邏輯判斷表達式。當表

16、達式的邏輯值為真時， 就執行循環體內的語句； 當表達式的邏輯值為假時， 就退出 當前的循環體。 如果循環判斷語句為矩陣時， 當且僅當所有的矩陣元素非零時， 邏輯表達式 的值為真。為避免循環語句陷入死循環， 在語句內必須有可以自動修改循環控制變量的命令。要求出上升時間，可以用 while 語句編寫以下程序得到：C=dcgain(G);n=1;while y(n)<Cn=n+1;endrisetime=t(n)在階躍輸入條件下， y 的值由零逐漸增大，當以上循環滿足 y=C 時，退出循環，此時對 應的時刻，即為上升時間。對于輸出無超調的系統響應，上升時間定義為輸出從穩態值的10%上升到 90

17、%所需時間，則計算程序如下：C=dcgain(G);n=1;while y(n)<*Cn=n+1;endm=1;while y(n)<*Cm=m+1;endrisetime=t(m)-t (n) 調節時間(setllingtime)可由while語句編程得到：C=dcgai n( G);i=le ngth(t);while(y(i)>*C)&(y(i)<*C)i=i-1;endsetlli ngtime=t(i)i的上限值。用向量長度函數length()可求得t序列的長度，將其設定為變量例已知二階系統傳遞函數為：G(s)3(s 13i)(s 1 3i)利用下面的

18、程序可得到階躍響應如圖3-3及性能指標數據。G(S)=3SA2+2S+10)>> nu m=0,0,3;den=1,2,10;step( nu m,de n)gridtitle( 1Un it-Step Respo nse of G(s)=3/(sA2+2s+10) 1)或者G=tf(0,0,3,1,2,10);t=0:5;c=step(G,t);plot(t,c)>> Css=dcga in(G)Css =G=tf(0,0,3,1,2,10);>> y,t=step(G)y =0>> Y,k=max(y);>> timetopeak=t(k)timetopeak =>> C=dcgain(G);>> Y,k=max(