實驗二利用MATLAB進(jìn)行時域分析

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

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

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

4、(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為各極點的留數(shù)，k(s)為余項。例設(shè)傳遞函數(shù)為:G(s)2s3 5s2 3s 6s3 6s211s 6該傳遞函數(shù)的部分分式展開由以下命令獲得: nu m=2,5,3,6;den=1,6,11,6;r,p,k=residue( nu m,de n)命令窗口中顯示如下結(jié)果中留數(shù)為列向量r，極點為列向量p,余項為行向量 k。由此可得出部分分式展開式:64G(s)32該函數(shù)也可以逆向調(diào)用,把部分分式展開轉(zhuǎn)變回多項式器之比的形式，命令格式為:nu m,de

5、 n=residue(r,p,k)對上例有： nu m,de n=residue(r,p,k)結(jié)果顯示num =den =應(yīng)當(dāng)指出，如果p(j)=p(j+1)= -=p(j+m-1),則極點p(j)是一個m重極點。在這種情況下，部分分式展開式將包括下列諸項:r(j)s p( j)r(j 1)r(j m 1)s p(j) 2p(j)例設(shè)傳遞函數(shù)為:G(s)s2 2s 3(s 1)3s2 2s 3s3 3s2 3s 1則部分分式展開由以下命令獲得: v=-1,-1,-1num=0,1,2,3;den=poly(v);r,p,k=residue( nu m,de n)結(jié)果顯示v =-1-1 -1其中

6、由poly()命令將分母化為標(biāo)準(zhǔn)降幕排列多項式系數(shù)向量由上可得展開式為：den, k=為空矩陣。G(s)1 02s 1 (s 1)2(s 1)3(5) 由傳遞函數(shù)求零點和極點。在MATLAB控制系統(tǒng)工具箱中，給出了由傳遞函數(shù)對象G求出系統(tǒng)零點和極點的函數(shù)，其調(diào)用格式分別為：Z=tzero(G)P=1注意：式19中要求的G必須是零極點模型對象，且出現(xiàn)了矩陣的點運算“”和大括號訂表示的矩陣元素，詳細(xì)內(nèi)容參閱后面章節(jié)。例已知傳遞函數(shù)為：2、6.8s61.2s95.2G(s) 3s47.5s3 22s2 19.5s輸入如下命令：num=,;den=1,22,0;G=tf( nu m,de n);G1=

7、zpk(G);Z=tzero(G)P=1結(jié)果顯示Z =-7-2(6) 零極點分布圖在MATLAB中，可利用pzmap()函數(shù)繪制連續(xù)系統(tǒng)的零、極點圖，從而分析系統(tǒng)的穩(wěn)定性，該函數(shù)調(diào)用格式為：pzmap( nu m,de n)例給定傳遞函數(shù)：3s4 2s3 5s2 4s 6s5 3s4 4s3 2s2 7s 2利用下列命令可自動打開一個圖形窗口，顯示該系統(tǒng)的零、極點分布圖，如圖3-1所示。 nu m=3,2,5,4,6;den=1,3,4,2,7,2;pzmap (nu m,de n)title( 1 Pole-Zero Map1) % 圖形標(biāo)題。1.5cfxa wanraay-1.5-2-IL

8、X-一-rr-1.5-1-0.5Real Axis0Pole-Zero Map5 oo5 o-0.5系統(tǒng)動態(tài)特性分析?？刂葡到y(tǒng)軟件包提供了控制系統(tǒng)工程需要的基本的時域與頻域分析工具函數(shù)。連續(xù)時間系統(tǒng)分析函數(shù)Impulse脈沖響應(yīng)Step階躍響應(yīng)Lsim任意輸入的仿真Bode波特圖Nyquist奈奎斯特圖Lyap李雅普諾夫方程Gram可控性與可觀性(1)時域響應(yīng)解析算法部分分式展開法用拉氏變換法求系統(tǒng)的單位階躍響應(yīng)，可直接得出輸出c(t)隨時間t變化的規(guī)律，對于高階系統(tǒng)，輸出的拉氏變換象函數(shù)為：1num 1C(s) G(s)sden s對函數(shù)c(s)進(jìn)行部分分式展開，我們可以用num,den,0

9、來表示c(s)的分子和分母。例給定系統(tǒng)的傳遞函數(shù)：、s3 7s2 24s 24G(s) 3s410s335s250s 24用以下命令對9進(jìn)行部分分式展開。s nu m,de n,0num=1,7,24,24den=1,10,35,50,24r,p,k=residue( nu m,de n, 0)輸出結(jié)果為num =24241103550r =den =240k =輸出函數(shù)c(s)為：C( s) =num/ ( den*s )num =1 72424ans =24 01103550拉氏變換得：L=laplace(sA3+7*s+24*s+24)/(sA4+10*sA3+35*sA2+50*s+2

10、4),s,t)(2 )單位階躍響應(yīng)的求法:控制系統(tǒng)工具箱中給出了一個函數(shù)step()來直接求取線性系統(tǒng)的階躍響應(yīng)，如果已知傳遞函數(shù)為:G(s)numden則該函數(shù)可有以下幾種調(diào)用格式:step (nu m,de n)step (nu m,de n,t)或(b)step(G)(c)step(G,t)(d)該函數(shù)將繪制出系統(tǒng)在單位階躍輸入條件下的動態(tài)響應(yīng)圖，同時給出穩(wěn)態(tài)值。對于式(b)和(d)，t為圖像顯示的時間長度，是用戶指定的時間向量。式(a)和(c)的顯示時間由系統(tǒng)根據(jù)輸出曲線的形狀自行設(shè)定。如果需要將輸出結(jié)果返回到MATLAB工作空間中，則采用以下調(diào)用格式:c=step(G)此時，屏上不會

11、顯示響應(yīng)曲線，必須利用plot()命令去查看響應(yīng)曲線。plot可以根據(jù)兩個或多個給定的向量繪制二維圖形，詳細(xì)介紹可以查閱后面的章節(jié)。例已知傳遞函數(shù)為:G(s)25s2 4s 251 .40.3050.4UllbSt即住訊皿梵 GT G(S23CSss-25)UU.51.522.5I me sec0圖3-2 MATLAB繪制的響應(yīng)曲線num=0,0,25;den=1,4,25;step( nu m,de n)grid %繪制網(wǎng)格線title( 1Unit-Step Response of G(s)=25(s2+4s+25) 1) % 圖像標(biāo)題我們還可以用下面的語句來得出階躍響應(yīng)曲線 G=tf(0

12、,0,25,1,4,25);t=0:5;%從0到5每隔取一個值。c=step(G,t); %動態(tài)響應(yīng)的幅值賦給變量cplot(t,c) %繪二維圖形，橫坐標(biāo)取t,縱坐標(biāo)取C。Css=dcgai n(G) % 求取穩(wěn)態(tài)值。系統(tǒng)顯示的圖形類似于上一個例子，在命令窗口中顯示了如下結(jié)果Css=1(3 )求階躍響應(yīng)的性能指標(biāo)MATLAB提供了強大的繪圖計算功能，可以用多種方法求取系統(tǒng)的動態(tài)響應(yīng)指標(biāo)。我們首先介紹一種最簡單的方法一一游動鼠標(biāo)法。對于例16,在程序運行完畢后，用鼠標(biāo)左鍵 點擊時域響應(yīng)圖線任意一點， 系統(tǒng)會自動跳出一個小方框， 小方框顯示了這一點的橫坐標(biāo) (時 間)和縱坐標(biāo)(幅值)。按住鼠標(biāo)左

13、鍵在曲線上移動，可以找到曲線幅值最大的一點一一即曲線最大峰值， 此時小方框中顯示的時間就是此二階系統(tǒng)的峰值時間， 根據(jù)觀察到的穩(wěn)態(tài)值 和峰值可以計算出系統(tǒng)的超調(diào)量。 系統(tǒng)的上升時間和穩(wěn)態(tài)響應(yīng)時間可以依此類推。 這種方法 簡單易用，但同時應(yīng)注意它不適用于用plot()命令畫出的圖形。另一種比較常用的方法就是用編程方式求取時域響應(yīng)的各項性能指標(biāo)。與上一段介紹的游動鼠標(biāo)法相比，編程方法稍微復(fù)雜，但通過下面的學(xué)習(xí)，讀者可以掌握一定的編程技巧， 能夠?qū)⒖刂圃碇R和編程方法相結(jié)合， 自己編寫一些程序， 獲取一些較為復(fù)雜的性能指標(biāo)。通過前面的學(xué)習(xí)， 我們已經(jīng)可以用階躍響應(yīng)函數(shù) step( )獲得系統(tǒng)輸出量

14、， 若將輸出量返 回到變量 y 中，可以調(diào)用如下格式y(tǒng),t=step(G)該函數(shù)還同時返回了自動生成的時間變量t，對返回的這一對變量 y和t的值進(jìn)行計算，可以得到時域性能指標(biāo)。 峰值時間 (timetopeak) 可由以下命令獲得：Y,k=max(y);timetopeak=t(k)應(yīng)用取最大值函數(shù) max()求出y的峰值及相應(yīng)的時間，并存于變量Y和k中。然后在變量 t 中取出峰值時間，并將它賦給變量 timetopeak 。 最大 (百分比 )超調(diào)量 (percentovershoot) 可由以下命令得到：C=dcgain(G);Y,k=max(y);percentovershoot=100

15、*(Y-C)/Cdcgain()函數(shù)用于求取系統(tǒng)的終值，將終值賦給變量C,然后依據(jù)超調(diào)量的定義，由Y 和C計算出百分比超調(diào)量。 上升時間(risetime)可利用MATLAB中控制語句編制 M文件來獲得。首先簡單介紹 一下循環(huán)語句 while 的使用。while 循環(huán)語句的一般格式為：while循環(huán)體end其中 ,循環(huán)判斷語句為某種形式的邏輯判斷表達(dá)式。當(dāng)表達(dá)式的邏輯值為真時， 就執(zhí)行循環(huán)體內(nèi)的語句； 當(dāng)表達(dá)式的邏輯值為假時， 就退出 當(dāng)前的循環(huán)體。 如果循環(huán)判斷語句為矩陣時， 當(dāng)且僅當(dāng)所有的矩陣元素非零時， 邏輯表達(dá)式 的值為真。為避免循環(huán)語句陷入死循環(huán)， 在語句內(nèi)必須有可以自動修改循環(huán)控制

16、變量的命令。要求出上升時間，可以用 while 語句編寫以下程序得到：C=dcgain(G);n=1;while y(n)Cn=n+1;endrisetime=t(n)在階躍輸入條件下， y 的值由零逐漸增大，當(dāng)以上循環(huán)滿足 y=C 時，退出循環(huán)，此時對 應(yīng)的時刻，即為上升時間。對于輸出無超調(diào)的系統(tǒng)響應(yīng)，上升時間定義為輸出從穩(wěn)態(tài)值的10%上升到 90%所需時間，則計算程序如下：C=dcgain(G);n=1;while y(n)*Cn=n+1;endm=1;while y(n)*C)&(y(i) 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); timetop