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

實(shí)驗(yàn)二利用MATLAB進(jìn)行時域分析實(shí)驗(yàn)二利用MATLAB進(jìn)行時域分析實(shí)驗(yàn)二利用MATLAB進(jìn)行時域分析實(shí)驗(yàn)二利用MATLAB進(jìn)行時域分析編制僅供參考審核批準(zhǔn)生效日期地址：電話：傳真:郵編:實(shí)驗(yàn)二利用MATLAB進(jìn)行時域分析本實(shí)驗(yàn)內(nèi)容包含以下三個部分：基于MATLAB的線性系統(tǒng)穩(wěn)定性分析、基于MATLAB的線性系統(tǒng)動態(tài)性能分析、和MATALB進(jìn)行控制系統(tǒng)時域分析的一些其它實(shí)例。一、基于MATLAB的線性系統(tǒng)穩(wěn)定性分析線性系統(tǒng)穩(wěn)定的充要條件是系統(tǒng)的特征根均位于S平面的左半部分。系統(tǒng)的零極點(diǎn)模型可以直接被用來判斷系統(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)例已知系統(tǒng)的特征多項式為：特征方程的解可由下面的MATLAB命令得出。>>p=[1,0,3,2,1,1];v=roots(p)結(jié)果顯示：v=+-+-利用多項式求根函數(shù)roots(),可以很方便的求出系統(tǒng)的零點(diǎn)和極點(diǎn)，然后根據(jù)零極點(diǎn)分析系統(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()是互為逆運(yùn)算的。（3）多項式求值在MATLAB中通過函數(shù)polyval()可以求得多項式在給定點(diǎn)的值，該函數(shù)的調(diào)用格式為：polyval(p,v)對于上例中的p值，求取多項式在x點(diǎn)的值，可輸入如下命令：>>p=[1,0,3,2,1,1];x=1polyval(p,x)結(jié)果顯示x=1ans=8(4)部分分式展開考慮下列傳遞函數(shù)：式中，但是和中某些量可能為零。MATLAB函數(shù)可將展開成部分分式，直接求出展開式中的留數(shù)、極點(diǎn)和余項。該函數(shù)的調(diào)用格式為：則的部分分式展開由下式給出：式中,,…,,為極點(diǎn)，,,…,為各極點(diǎn)的留數(shù)，為余項。例設(shè)傳遞函數(shù)為：該傳遞函數(shù)的部分分式展開由以下命令獲得：>>num=[2,5,3,6];den=[1,6,11,6];[r,p,k]=residue(num,den)命令窗口中顯示如下結(jié)果r=p=k=2中留數(shù)為列向量r，極點(diǎn)為列向量p，余項為行向量k。由此可得出部分分式展開式：該函數(shù)也可以逆向調(diào)用，把部分分式展開轉(zhuǎn)變回多項式之比的形式，命令格式為：[num,den]=residue(r,p,k)對上例有：>>[num,den]=residue(r,p,k)結(jié)果顯示num=den=應(yīng)當(dāng)指出，如果p(j)=p(j+1)=…=p(j+m-1),則極點(diǎn)p(j)是一個m重極點(diǎn)。在這種情況下，部分分式展開式將包括下列諸項：例設(shè)傳遞函數(shù)為：則部分分式展開由以下命令獲得：>>v=[-1,-1,-1]num=[0,1,2,3];den=poly(v);[r,p,k]=residue(num,den)結(jié)果顯示v=-1-1-1r=p=k=[]其中由poly()命令將分母化為標(biāo)準(zhǔn)降冪排列多項式系數(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={1}注意：式19中要求的G必須是零極點(diǎn)模型對象，且出現(xiàn)了矩陣的點(diǎn)運(yùn)算“.”和大括號｛｝表示的矩陣元素，詳細(xì)內(nèi)容參閱后面章節(jié)。例已知傳遞函數(shù)為：輸入如下命令：num=[,,];den=[1,,22,,0];G=tf(num,den);G1=zpk(G);Z=tzero(G)P={1}結(jié)果顯示Z=-7-2P=0+-（6）零極點(diǎn)分布圖在MATLAB中，可利用pzmap()函數(shù)繪制連續(xù)系統(tǒng)的零、極點(diǎn)圖，從而分析系統(tǒng)的穩(wěn)定性，該函數(shù)調(diào)用格式為：pzmap(num,den)例給定傳遞函數(shù):利用下列命令可自動打開一個圖形窗口，顯示該系統(tǒng)的零、極點(diǎn)分布圖，如圖3-1所示。>>num=[3,2,5,4,6];den=[1,3,4,2,7,2];pzmap(num,den)title(1Pole-ZeroMap1)%圖形標(biāo)題。二.系統(tǒng)動態(tài)特性分析。控制系統(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ù)為：對函數(shù)c(s)進(jìn)行部分分式展開，我們可以用num,[den,0]來表示c(s)的分子和分母。例給定系統(tǒng)的傳遞函數(shù)：用以下命令對進(jìn)行部分分式展開。>>num,[den,0]num=[1,7,24,24]den=[1,10,35,50,24][r,p,k]=residue(num,[den,0])輸出結(jié)果為num=172424den=110355024r=p=0k=[]輸出函數(shù)c(s)為：C（s）=num/（den*s）num=172424ans=1103550240拉氏變換得：L=laplace((s^3+7*s+24*s+24)/(s^4+10*s^3+35*s^2+50*s+24),s,t)（2）單位階躍響應(yīng)的求法：控制系統(tǒng)工具箱中給出了一個函數(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)在單位階躍輸入條件下的動態(tài)響應(yīng)圖，同時給出穩(wěn)態(tài)值。對于式(b)和(d)，t為圖像顯示的時間長度，是用戶指定的時間向量。式(a)和(c)的顯示時間由系統(tǒng)根據(jù)輸出曲線的形狀自行設(shè)定。如果需要將輸出結(jié)果返回到MATLAB工作空間中，則采用以下調(diào)用格式：c=step(G)此時，屏上不會顯示響應(yīng)曲線，必須利用plot()命令去查看響應(yīng)曲線。plot可以根據(jù)兩個或多個給定的向量繪制二維圖形，詳細(xì)介紹可以查閱后面的章節(jié)。例已知傳遞函數(shù)為：利用以下MATLAB命令可得階躍響應(yīng)曲線如圖3-2所示。圖3-2MATLAB繪制的響應(yīng)曲線>>num=[0,0,25];den=[1,4,25];step(num,den)grid%繪制網(wǎng)格線title(1Unit-StepResponseofG(s)=25/(s^2+4s+25)1)%圖像標(biāo)題我們還可以用下面的語句來得出階躍響應(yīng)曲線>>G=tf([0,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=dcgain(G)%求取穩(wěn)態(tài)值。系統(tǒng)顯示的圖形類似于上一個例子，在命令窗口中顯示了如下結(jié)果Css=1（3）求階躍響應(yīng)的性能指標(biāo)MATLAB提供了強(qiáng)大的繪圖計算功能，可以用多種方法求取系統(tǒng)的動態(tài)響應(yīng)指標(biāo)。我們首先介紹一種最簡單的方法――游動鼠標(biāo)法。對于例16，在程序運(yùn)行完畢后，用鼠標(biāo)左鍵點(diǎn)擊時域響應(yīng)圖線任意一點(diǎn)，系統(tǒng)會自動跳出一個小方框，小方框顯示了這一點(diǎn)的橫坐標(biāo)（時間）和縱坐標(biāo)（幅值）。按住鼠標(biāo)左鍵在曲線上移動，可以找到曲線幅值最大的一點(diǎn)――即曲線最大峰值，此時小方框中顯示的時間就是此二階系統(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)輸出量，若將輸出量返回到變量y中，可以調(diào)用如下格式[y,t]=step(G)該函數(shù)還同時返回了自動生成的時間變量t，對返回的這一對變量y和t的值進(jìn)行計算，可以得到時域性能指標(biāo)。=1\*GB3①峰值時間(timetopeak)可由以下命令獲得：[Y,k]=max(y);timetopeak=t(k)應(yīng)用取最大值函數(shù)max()求出y的峰值及相應(yīng)的時間，并存于變量Y和k中。然后在變量t中取出峰值時間，并將它賦給變量timetopeak。=2\*GB3②最大(百分比)超調(diào)量(percentovershoot)可由以下命令得到：C=dcgain(G);[Y,k]=max(y);percentovershoot=100*(Y-C)/Cdcgain()函數(shù)用于求取系統(tǒng)的終值，將終值賦給變量C，然后依據(jù)超調(diào)量的定義，由Y和C計算出百分比超調(diào)量。=3\*GB3③上升時間(risetime)可利用MATLAB中控制語句編制M文件來獲得。首先簡單介紹一下循環(huán)語句while的使用。while循環(huán)語句的一般格式為：while<循環(huán)判斷語句>循環(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)控制變量的命令。要求出上升時間，可以用while語句編寫以下程序得到：C=dcgain(G);n=1;whiley(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;whiley(n)<*Cn=n+1;endm=1;whiley(n)<*Cm=m+1;endrisetime=t(m)-t(n)=4\*GB3④調(diào)節(jié)時間(setllingtime)可由while語句編程得到：C=dcgain(G);i=length(t);while(y(i)>*C)&(y(i)<*C)i=i-1;endsetllingtime=t(i)用向量長度函數(shù)length()可求得t序列的長度，將其設(shè)定為變量i的上限值。例已知二階系統(tǒng)傳遞函數(shù)為:利用下面的程序可得到階躍響應(yīng)如圖3-3及性能指標(biāo)數(shù)據(jù)。G(S)=3/(S^2+2S+10)>>num=[0,0,3];den=[1,2,10];step(num,den)gridtitle(1Unit-StepResponseofG(s)=3/(s^2+2s+10)1)或者G=tf([0,0,3],[1,2,10]);t=0::5;c=step(G,t);plot(t,c)>>Css=dcgain(G)Css=G=tf([0,0,3],[1,2,10]);>>[y,t]=step(G)y=0