實(shí)驗(yàn)二 利用MATLAB進(jìn)行時(shí)域分析參考模板

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

2、 已知系統(tǒng)的特征多項(xiàng)式為： 特征方程的解可由下面的MATLAB命令得出。>> p=1,0,3,2,1,1; v=roots(p)結(jié)果顯示：v = 0.3202 + 1.7042i 0.3202 - 1.7042i -0.7209 0.0402 + 0.6780i 0.0402 - 0.6780i利用多項(xiàng)式求根函數(shù)roots(),可以很方便的求出系統(tǒng)的零點(diǎn)和極點(diǎn)，然后根據(jù)零極點(diǎn)分析系統(tǒng)穩(wěn)定性和其它性能。（2）由根創(chuàng)建多項(xiàng)式如果已知多項(xiàng)式的因式分解式或特征根，可由MATLAB函數(shù)poly()直接得出特征多項(xiàng)式系數(shù)向量，其調(diào)用格式為：p=poly(v) 如上例中：v=0.3202+1.7

3、042i;0.3202-1.7042i; -0.7209;0.0402+0.6780i; 0.0402-0.6780i; >> p=poly(v)結(jié)果顯示p =1 / 19 1.0000 0.0001 3.0000 2.0001 0.9998 0.9999由此可見，函數(shù)roots()與函數(shù)poly()是互為逆運(yùn)算的。（3）多項(xiàng)式求值在MATLAB 中通過函數(shù)polyval()可以求得多項(xiàng)式在給定點(diǎn)的值，該函數(shù)的調(diào)用格式為： polyval(p,v) 對于上例中的p值，求取多項(xiàng)式在x點(diǎn)的值，可輸入如下命令：>> p=1,0,3,2,1,1;x=1 polyval(p,x)

4、結(jié)果顯示x = 1ans = 8 (4)部分分式展開 考慮下列傳遞函數(shù)： 式中，但是和中某些量可能為零。MATLAB函數(shù)可將展開成部分分式，直接求出展開式中的留數(shù)、極點(diǎn)和余項(xiàng)。該函數(shù)的調(diào)用格式為： 則的部分分式展開由下式給出： 式中, ,為極點(diǎn)，, 為各極點(diǎn)的留數(shù)，為余項(xiàng)。例3.2 設(shè)傳遞函數(shù)為： 該傳遞函數(shù)的部分分式展開由以下命令獲得：>> num=2,5,3,6; den=1,6,11,6; r,p,k=residue(num,den)命令窗口中顯示如下結(jié)果r = -6.0000 -4.0000 3.0000p = -3.0000 -2.0000 -1.0000k = 2中留數(shù)

5、為列向量r，極點(diǎn)為列向量p，余項(xiàng)為行向量k。由此可得出部分分式展開式：該函數(shù)也可以逆向調(diào)用，把部分分式展開轉(zhuǎn)變回多項(xiàng)式之比的形式，命令格式為： num,den=residue(r,p,k) 對上例有： >> num,den=residue(r,p,k)結(jié)果顯示num = 2.0000 5.0000 3.0000 6.0000den = 1.0000 6.0000 11.0000 6.0000應(yīng)當(dāng)指出，如果p(j)=p(j+1)=p(j+m-1),則極點(diǎn)p(j)是一個(gè)m重極點(diǎn)。在這種情況下，部分分式展開式將包括下列諸項(xiàng)： 例3.3 設(shè)傳遞函數(shù)為： 則部分分式展開由以下命令獲得：>

6、;> v=-1,-1,-1num=0,1,2,3; den=poly(v);r,p,k=residue(num,den)結(jié)果顯示v = -1 -1 -1r = 1.0000 0.0000 2.0000p = -1.0000 -1.0000 -1.0000k = 其中由poly()命令將分母化為標(biāo)準(zhǔn)降冪排列多項(xiàng)式系數(shù)向量den, k=為空矩陣。由上可得展開式為： （5）由傳遞函數(shù)求零點(diǎn)和極點(diǎn)。在MATLAB控制系統(tǒng)工具箱中，給出了由傳遞函數(shù)對象G求出系統(tǒng)零點(diǎn)和極點(diǎn)的函數(shù)，其調(diào)用格式分別為：Z=tzero(G) P=G.P1 注意：式19中要求的G必須是零極點(diǎn)模型對象，且出現(xiàn)了矩陣的點(diǎn)運(yùn)算“

7、.”和大括號(hào)表示的矩陣元素，詳細(xì)內(nèi)容參閱后面章節(jié)。例3.4 已知傳遞函數(shù)為： 輸入如下命令：num=6.8,61.2,95.2;den=1,7.5,22,19.5,0;G=tf(num,den);G1=zpk(G);Z=tzero(G)P=G1.P1結(jié)果顯示Z = -7 -2P = 0 -3.0000 + 2.0000i -3.0000 - 2.0000i -1.5000 （6）零極點(diǎn)分布圖 在MATLAB中，可利用pzmap()函數(shù)繪制連續(xù)系統(tǒng)的零、極點(diǎn)圖，從而分析系統(tǒng)的穩(wěn)定性，該函數(shù)調(diào)用格式為：pzmap(num,den) 例 3.5 給定傳遞函數(shù): 利用下列命令可自動(dòng)打開一個(gè)圖形窗口，顯

8、示該系統(tǒng)的零、極點(diǎn)分布圖，如圖3-1所示。>> num=3,2,5,4,6; den=1,3,4,2,7,2;pzmap(num,den)title(¹Pole-Zero Map¹)% 圖形標(biāo)題。二. 系統(tǒng)動(dòng)態(tài)特性分析?？刂葡到y(tǒng)軟件包提供了控制系統(tǒng)工程需要的基本的時(shí)域與頻域分析工具函數(shù)。連續(xù)時(shí)間系統(tǒng)分析函數(shù)Impulse脈沖響應(yīng)Step階躍響應(yīng)Lsim任意輸入的仿真Bode波特圖Nyquist奈奎斯特圖Lyap李雅普諾夫方程Gram可控性與可觀性（1）時(shí)域響應(yīng)解析算法部分分式展開法用拉氏變換法求系統(tǒng)的單位階躍響應(yīng)，可直接得出輸出c(t)隨時(shí)間t變化的規(guī)律，對于高

9、階系統(tǒng)，輸出的拉氏變換象函數(shù)為： 對函數(shù)c(s)進(jìn)行部分分式展開，我們可以用num,den,0來表示c(s)的分子和分母。例 3.6 給定系統(tǒng)的傳遞函數(shù)： 用以下命令對進(jìn)行部分分式展開。>> num,den,0num=1,7,24,24den=1,10,35,50,24r,p,k=residue(num,den,0)輸出結(jié)果為num = 1 7 24 24den = 1 10 35 50 24r = -1.0000 2.0000 -1.0000 -1.0000 1.0000p = -4.0000 -3.0000 -2.0000 -1.0000 0k = 輸出函數(shù)c(s)為： C（s

10、）=num/（den*s） num = 1 7 24 24ans = 1 10 35 50 24 0拉氏變換得：L=laplace(s3+7*s+24*s+24)/(s4+10*s3+35*s2+50*s+24),s,t)（2）單位階躍響應(yīng)的求法：控制系統(tǒng)工具箱中給出了一個(gè)函數(shù)step()來直接求取線性系統(tǒng)的階躍響應(yīng)，如果已知傳遞函數(shù)為：則該函數(shù)可有以下幾種調(diào)用格式：step(num,den) (a)step(num,den,t) (b)或 step(G) (c)step(G,t) (d)該函數(shù)將繪制出系統(tǒng)在單位階躍輸入條件下的動(dòng)態(tài)響應(yīng)圖，同時(shí)給出穩(wěn)態(tài)值。對于式(b)和(d)，t為圖像顯示的時(shí)

11、間長度，是用戶指定的時(shí)間向量。式(a)和(c)的顯示時(shí)間由系統(tǒng)根據(jù)輸出曲線的形狀自行設(shè)定。如果需要將輸出結(jié)果返回到MATLAB工作空間中，則采用以下調(diào)用格式： c=step(G) 此時(shí)，屏上不會(huì)顯示響應(yīng)曲線，必須利用plot()命令去查看響應(yīng)曲線。plot 可以根據(jù)兩個(gè)或多個(gè)給定的向量繪制二維圖形，詳細(xì)介紹可以查閱后面的章節(jié)。例3.7 已知傳遞函數(shù)為： 利用以下MATLAB命令可得階躍響應(yīng)曲線如圖3-2所示。圖3-2 MATLAB繪制的響應(yīng)曲線>> num=0,0,25; den=1,4,25;step(num,den)grid % 繪制網(wǎng)格線title(¹Unit-St

12、ep Response of G(s)=25/(s2+4s+25) ¹) % 圖像標(biāo)題我們還可以用下面的語句來得出階躍響應(yīng)曲線>> G=tf(0,0,25,1,4,25); t=0:0.1:5; % 從0到5每隔0.1取一個(gè)值。 c=step(G,t);% 動(dòng)態(tài)響應(yīng)的幅值賦給變量c plot(t,c) % 繪二維圖形，橫坐標(biāo)取t,縱坐標(biāo)取c。 Css=dcgain(G) % 求取穩(wěn)態(tài)值。系統(tǒng)顯示的圖形類似于上一個(gè)例子，在命令窗口中顯示了如下結(jié)果 Css= 1（3）求階躍響應(yīng)的性能指標(biāo)MATLAB提供了強(qiáng)大的繪圖計(jì)算功能，可以用多種方法求取系統(tǒng)的動(dòng)態(tài)響應(yīng)指標(biāo)。我們首先介紹一

13、種最簡單的方法游動(dòng)鼠標(biāo)法。對于例16，在程序運(yùn)行完畢后，用鼠標(biāo)左鍵點(diǎn)擊時(shí)域響應(yīng)圖線任意一點(diǎn)，系統(tǒng)會(huì)自動(dòng)跳出一個(gè)小方框，小方框顯示了這一點(diǎn)的橫坐標(biāo)（時(shí)間）和縱坐標(biāo)（幅值）。按住鼠標(biāo)左鍵在曲線上移動(dòng)，可以找到曲線幅值最大的一點(diǎn)即曲線最大峰值，此時(shí)小方框中顯示的時(shí)間就是此二階系統(tǒng)的峰值時(shí)間，根據(jù)觀察到的穩(wěn)態(tài)值和峰值可以計(jì)算出系統(tǒng)的超調(diào)量。系統(tǒng)的上升時(shí)間和穩(wěn)態(tài)響應(yīng)時(shí)間可以依此類推。這種方法簡單易用，但同時(shí)應(yīng)注意它不適用于用plot()命令畫出的圖形。另一種比較常用的方法就是用編程方式求取時(shí)域響應(yīng)的各項(xiàng)性能指標(biāo)。與上一段介紹的游動(dòng)鼠標(biāo)法相比，編程方法稍微復(fù)雜，但通過下面的學(xué)習(xí)，讀者可以掌握一定的編程技巧

14、，能夠?qū)⒖刂圃碇R(shí)和編程方法相結(jié)合，自己編寫一些程序，獲取一些較為復(fù)雜的性能指標(biāo)。通過前面的學(xué)習(xí)，我們已經(jīng)可以用階躍響應(yīng)函數(shù)step( )獲得系統(tǒng)輸出量，若將輸出量返回到變量y中，可以調(diào)用如下格式 y,t=step(G) 該函數(shù)還同時(shí)返回了自動(dòng)生成的時(shí)間變量t，對返回的這一對變量y和t的值進(jìn)行計(jì)算，可以得到時(shí)域性能指標(biāo)。 峰值時(shí)間(timetopeak)可由以下命令獲得： Y,k=max(y); timetopeak=t(k) 應(yīng)用取最大值函數(shù)max()求出y的峰值及相應(yīng)的時(shí)間，并存于變量Y和k中。然后在變量t中取出峰值時(shí)間，并將它賦給變量timetopeak。 最大(百分比)超調(diào)量(per

15、centovershoot)可由以下命令得到：C=dcgain(G);Y,k=max(y); percentovershoot=100*(Y-C)/C dcgain( )函數(shù)用于求取系統(tǒng)的終值，將終值賦給變量C，然后依據(jù)超調(diào)量的定義，由Y和C計(jì)算出百分比超調(diào)量。 上升時(shí)間(risetime)可利用MATLAB中控制語句編制M文件來獲得。首先簡單介紹一下循環(huán)語句while的使用。while循環(huán)語句的一般格式為：while<循環(huán)判斷語句> 循環(huán)體end其中,循環(huán)判斷語句為某種形式的邏輯判斷表達(dá)式。當(dāng)表達(dá)式的邏輯值為真時(shí)，就執(zhí)行循環(huán)體內(nèi)的語句；當(dāng)表達(dá)式的邏輯值為假時(shí)，就退出當(dāng)前的循環(huán)體。

16、如果循環(huán)判斷語句為矩陣時(shí)，當(dāng)且僅當(dāng)所有的矩陣元素非零時(shí)，邏輯表達(dá)式的值為真。為避免循環(huán)語句陷入死循環(huán)，在語句內(nèi)必須有可以自動(dòng)修改循環(huán)控制變量的命令。要求出上升時(shí)間，可以用while語句編寫以下程序得到：C=dcgain(G);n=1;while y(n)<C n=n+1;endrisetime=t(n)在階躍輸入條件下，y 的值由零逐漸增大，當(dāng)以上循環(huán)滿足y=C時(shí)，退出循環(huán)，此時(shí)對應(yīng)的時(shí)刻，即為上升時(shí)間。對于輸出無超調(diào)的系統(tǒng)響應(yīng)，上升時(shí)間定義為輸出從穩(wěn)態(tài)值的10%上升到90%所需時(shí)間，則計(jì)算程序如下：C=dcgain(G);n=1; while y(n)<0.1*C n=n+1;

17、endm=1; while y(n)<0.9*C m=m+1; endrisetime=t(m)-t(n) 調(diào)節(jié)時(shí)間(setllingtime)可由while語句編程得到：C=dcgain(G);i=length(t); while(y(i)>0.98*C)&(y(i)<1.02*C) i=i-1;endsetllingtime=t(i)用向量長度函數(shù)length( )可求得t序列的長度，將其設(shè)定為變量i的上限值。例 3.8 已知二階系統(tǒng)傳遞函數(shù)為: 利用下面的stepanalysis.m程序可得到階躍響應(yīng)如圖 3-3及性能指標(biāo)數(shù)據(jù)。 G(S)=3/(S2+2S+10

18、)>> num=0,0,3; den=1,2,10;step(num,den)grid title(¹Unit-Step Response of G(s)=3/(s2+2s+10) ¹)或者G=tf(0,0,3,1,2,10);t=0:0.1:5;c=step(G,t);plot(t,c)>> Css=dcgain(G)Css = 0.3000G=tf(0,0,3,1,2,10);>> y,t=step(G)y = 0 0.0044 0.0168 0.0361 0.0610 0.0902 0.1223 0.1562 0.1907 0.22

19、48 0.2574 0.2879 0.3155 0.3398 0.3604 0.3771 0.3898 0.3987 0.4037 0.4053 0.4036 0.3992 0.3923 0.3835 0.3733 0.3620 0.3500 0.3379 0.3260 0.3146 0.3039 0.2943 0.2858 0.2786 0.2728 0.2683 0.2653 0.2636 0.2631 0.2637 0.2653 0.2677 0.2708 0.2744 0.2784 0.2826 0.2868 0.2910 0.2950 0.2987 0.3021 0.3051 0.3

20、076 0.3096 0.3112 0.3122 0.3128 0.3130 0.3127 0.3122 0.3113 0.3102 0.3089 0.3075 0.3061 0.3046 0.3031 0.3017 0.3004 0.2992 0.2982 0.2973 0.2966 0.2961 0.2957 0.2955 0.2955 0.2955 0.2957 0.2960 0.2964 0.2969 0.2974 0.2979 0.2984 0.2989 0.2994 0.2999 0.3003 0.3006 0.3010 0.3012 0.3014 0.3015 0.3016 0.

21、3016 0.3016 0.3015 0.3014 0.3012 0.3011 0.3009 0.3007 0.3006 0.3004 0.3002 0.3000 0.2999 0.2998t = 0 0.0552 0.1104 0.1656 0.2209 0.2761 0.3313 0.3865 0.4417 0.4969 0.5521 0.6074 0.6626 0.7178 0.7730 0.8282 0.8834 0.9386 0.9939 1.0491 1.1043 1.1595 1.2147 1.2699 1.3252 1.3804 1.4356 1.4908 1.5460 1.6012 1.6564 1.7117 1.7669 1.8221 1.8773 1.9325 1.9877 2.0429 2.0982 2.1534 2.2086 2.2638