實驗二 利用MATLAB進行時域分析

實驗二利用MATLAB進行時域分析實驗二利用MATLAB進行時域分析實驗二利用MATLAB進行時域分析實驗二利用MATLAB進行時域分析編制僅供參考審核批準生效日期地址：電話：傳真:郵編:實驗二利用MATLAB進行時域分析本實驗內(nèi)容包含以下三個部分：基于MATLAB的線性系統(tǒng)穩(wěn)定性分析、基于MATLAB的線性系統(tǒng)動態(tài)性能分析、和MATALB進行控制系統(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在復數(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)的零點和極點，然后根據(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ù)polyval()可以求得多項式在給定點的值，該函數(shù)的調(diào)用格式為：polyval(p,v)對于上例中的p值，求取多項式在x點的值，可輸入如下命令：>>p=[1,0,3,2,1,1];x=1polyval(p,x)結(jié)果顯示x=1ans=8(4)部分分式展開考慮下列傳遞函數(shù)：式中，但是和中某些量可能為零。MATLAB函數(shù)可將展開成部分分式，直接求出展開式中的留數(shù)、極點和余項。該函數(shù)的調(diào)用格式為：則的部分分式展開由下式給出：式中,,…,,為極點，,,…,為各極點的留數(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，極點為列向量p，余項為行向量k。由此可得出部分分式展開式：該函數(shù)也可以逆向調(diào)用，把部分分式展開轉(zhuǎn)變回多項式之比的形式，命令格式為：[num,den]=residue(r,p,k)對上例有：>>[num,den]=residue(r,p,k)結(jié)果顯示num=den=應當指出，如果p(j)=p(j+1)=…=p(j+m-1),則極點p(j)是一個m重極點。在這種情況下，部分分式展開式將包括下列諸項：例設(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()命令將分母化為標準降冪排列多項式系數(shù)向量den,k=[]為空矩陣。由上可得展開式為：（5）由傳遞函數(shù)求零點和極點。在MATLAB控制系統(tǒng)工具箱中，給出了由傳遞函數(shù)對象G求出系統(tǒng)零點和極點的函數(shù)，其調(diào)用格式分別為：Z=tzero(G)P={1}注意：式19中要求的G必須是零極點模型對象，且出現(xiàn)了矩陣的點運算“.”和大括號｛｝表示的矩陣元素，詳細內(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）零極點分布圖在MATLAB中，可利用pzmap()函數(shù)繪制連續(xù)系統(tǒng)的零、極點圖，從而分析系統(tǒng)的穩(wěn)定性，該函數(shù)調(diào)用格式為：pzmap(num,den)例給定傳遞函數(shù):利用下列命令可自動打開一個圖形窗口，顯示該系統(tǒng)的零、極點分布圖，如圖3-1所示。>>num=[3,2,5,4,6];den=[1,3,4,2,7,2];pzmap(num,den)title(1Pole-ZeroMap1)%圖形標題。二.系統(tǒng)動態(tài)特性分析?？刂葡到y(tǒng)軟件包提供了控制系統(tǒng)工程需要的基本的時域與頻域分析工具函數(shù)。連續(xù)時間系統(tǒng)分析函數(shù)Impulse脈沖響應Step階躍響應Lsim任意輸入的仿真Bode波特圖Nyquist奈奎斯特圖Lyap李雅普諾夫方程Gram可控性與可觀性（1）時域響應解析算法――部分分式展開法用拉氏變換法求系統(tǒng)的單位階躍響應，可直接得出輸出c(t)隨時間t變化的規(guī)律，對于高階系統(tǒng)，輸出的拉氏變換象函數(shù)為：對函數(shù)c(s)進行部分分式展開，我們可以用num,[den,0]來表示c(s)的分子和分母。例給定系統(tǒng)的傳遞函數(shù)：用以下命令對進行部分分式展開。>>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）單位階躍響應的求法：控制系統(tǒng)工具箱中給出了一個函數(shù)step()來直接求取線性系統(tǒ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)響應圖，同時給出穩(wěn)態(tài)值。對于式(b)和(d)，t為圖像顯示的時間長度，是用戶指定的時間向量。式(a)和(c)的顯示時間由系統(tǒng)根據(jù)輸出曲線的形狀自行設(shè)定。如果需要將輸出結(jié)果返回到MATLAB工作空間中，則采用以下調(diào)用格式：c=step(G)此時，屏上不會顯示響應曲線，必須利用plot()命令去查看響應曲線。plot可以根據(jù)兩個或多個給定的向量繪制二維圖形，詳細介紹可以查閱后面的章節(jié)。例已知傳遞函數(shù)為：利用以下MATLAB命令可得階躍響應曲線如圖3-2所示。圖3-2MATLAB繪制的響應曲線>>num=[0,0,25];den=[1,4,25];step(num,den)grid%繪制網(wǎng)格線title(1Unit-StepResponseofG(s)=25/(s^2+4s+25)1)%圖像標題我們還可以用下面的語句來得出階躍響應曲線>>G=tf([0,0,25],[1,4,25]);t=0::5;%從0到5每隔取一個值。c=step(G,t);%動態(tài)響應的幅值賦給變量cplot(t,c)%繪二維圖形，橫坐標取t,縱坐標取c。Css=dcgain(G)%求取穩(wěn)態(tài)值。系統(tǒng)顯示的圖形類似于上一個例子，在命令窗口中顯示了如下結(jié)果Css=1（3）求階躍響應的性能指標MATLAB提供了強大的繪圖計算功能，可以用多種方法求取系統(tǒng)的動態(tài)響應指標。我們首先介紹一種最簡單的方法――游動鼠標法。對于例16，在程序運行完畢后，用鼠標左鍵點擊時域響應圖線任意一點，系統(tǒng)會自動跳出一個小方框，小方框顯示了這一點的橫坐標（時間）和縱坐標（幅值）。按住鼠標左鍵在曲線上移動，可以找到曲線幅值最大的一點――即曲線最大峰值，此時小方框中顯示的時間就是此二階系統(tǒng)的峰值時間，根據(jù)觀察到的穩(wěn)態(tài)值和峰值可以計算出系統(tǒng)的超調(diào)量。系統(tǒng)的上升時間和穩(wěn)態(tài)響應時間可以依此類推。這種方法簡單易用，但同時應注意它不適用于用plot()命令畫出的圖形。另一種比較常用的方法就是用編程方式求取時域響應的各項性能指標。與上一段介紹的游動鼠標法相比，編程方法稍微復雜，但通過下面的學習，讀者可以掌握一定的編程技巧，能夠?qū)⒖刂圃碇R和編程方法相結(jié)合，自己編寫一些程序，獲取一些較為復雜的性能指標。通過前面的學習，我們已經(jīng)可以用階躍響應函數(shù)step()獲得系統(tǒng)輸出量，若將輸出量返回到變量y中，可以調(diào)用如下格式[y,t]=step(G)該函數(shù)還同時返回了自動生成的時間變量t，對返回的這一對變量y和t的值進行計算，可以得到時域性能指標。=1\*GB3①峰值時間(timetopeak)可由以下命令獲得：[Y,k]=max(y);timetopeak=t(k)應用取最大值函數(shù)max()求出y的峰值及相應的時間，并存于變量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)判斷語句為某種形式的邏輯判斷表達式。當表達式的邏輯值為真時，就執(zhí)行循環(huán)體內(nèi)的語句；當表達式的邏輯值為假時，就退出當前的循環(huán)體。如果循環(huán)判斷語句為矩陣時，當且僅當所有的矩陣元素非零時，邏輯表達式的值為真。為避免循環(huán)語句陷入死循環(huán)，在語句內(nèi)必須有可以自動修改循環(huán)控制變量的命令。要求出上升時間，可以用while語句編寫以下程序得到：C=dcgain(G);n=1;whiley(n)<Cn=n+1;endrisetime=t(n)在階躍輸入條件下，y的值由零逐漸增大，當以上循環(huán)滿足y=C時，退出循環(huán)，此時對應的時刻，即為上升時間。對于輸出無超調(diào)的系統(tǒ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ù)為:利用下面的程序可得到階躍響應如圖3-3及性能指標數(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