第5章用MATLAB進(jìn)行控制系統(tǒng)頻域分析

1、第5章 用MATLAB進(jìn)行控制系統(tǒng)頻域分析一、基于MATLAB的線性系統(tǒng)的頻域分析基本知識（）頻率特性函數(shù)。設(shè)線性系統(tǒng)傳遞函數(shù)為：則頻率特性函數(shù)為：由下面的MATLAB語句可直接求出G(jw)。i=sqrt(-1)%求取-1的平方根 GW=polyval(num，i*w)./polyval(den，i*w) 其中（num，den）為系統(tǒng)的傳遞函數(shù)模型。而w為頻率點(diǎn)構(gòu)成的向量，點(diǎn)右除（./）運(yùn)算符表示操作元素點(diǎn)對點(diǎn)的運(yùn)算。從數(shù)值運(yùn)算的角度來看，上述算法在系統(tǒng)的極點(diǎn)附近精度不會很理想，甚至出現(xiàn)無窮大值，運(yùn)算結(jié)果是一系列復(fù)數(shù)返回到變量GW中。（）用MATLAB作奈魁斯特圖。控制系統(tǒng)工具箱中提供了一個

2、MATLAB函數(shù)nyquist( )，該函數(shù)可以用來直接求解Nyquist陣列或繪制奈氏圖。當(dāng)命令中不包含左端返回變量時，nyquist（）函數(shù)僅在屏幕上產(chǎn)生奈氏圖，命令調(diào)用格式為：nyquist(num,den) nyquist(num,den,w) 或者nyquist(G) nyquist(G,w) 該命令將畫出下列開環(huán)系統(tǒng)傳遞函數(shù)的奈氏曲線： 如果用戶給出頻率向量w,則w包含了要分析的以弧度/秒表示的諸頻率點(diǎn)。在這些頻率點(diǎn)上，將對系統(tǒng)的頻率響應(yīng)進(jìn)行計算，若沒有指定的w向量，則該函數(shù)自動選擇頻率向量進(jìn)行計算。w包含了用戶要分析的以弧度/秒表示的諸頻率點(diǎn),MATLAB會自動計算這些點(diǎn)的頻率響

3、應(yīng)。當(dāng)命令中包含了左端的返回變量時，即:re,im,w=nyquist(G) 或 re,im,w=nyquist(G,w) 函數(shù)運(yùn)行后不在屏幕上產(chǎn)生圖形，而是將計算結(jié)果返回到矩陣re、im和w中。矩陣re和im分別表示頻率響應(yīng)的實部和虛部，它們都是由向量w中指定的頻率點(diǎn)計算得到的。在運(yùn)行結(jié)果中，w數(shù)列的每一個值分別對應(yīng)re、im數(shù)列的每一個值。例5.1 考慮二階典型環(huán)節(jié): 試?yán)肕ATLAB畫出奈氏圖。利用下面的命令，可以得出系統(tǒng)的奈氏圖，如圖5-1所示。 num=0,0,1;den=1,0.8,1;nyquist(num,den)% 設(shè)置坐標(biāo)顯示范圍v=-2,2,-2,2;axis(v)gr

4、id圖5-1 二階環(huán)節(jié)奈氏圖title(Nyquist Plot of G(s)=1/(s2+0.8s+1) )（3）用MATLAB作伯德圖控制系統(tǒng)工具箱里提供的bode()函數(shù)可以直接求取、繪制給定線性系統(tǒng)的伯德圖。當(dāng)命令不包含左端返回變量時，函數(shù)運(yùn)行后會在屏幕上直接畫出伯德圖。如果命令表達(dá)式的左端含有返回變量，bode()函數(shù)計算出的幅值和相角將返回到相應(yīng)的矩陣中，這時屏幕上不顯示頻率響應(yīng)圖。命令的調(diào)用格式為： mag,phase,w=bode(num,den) mag,phase,w=bode(num,den,w) 或 mag,phase,w=bode(G) mag,phase,w=bo

5、de(G,w) 矩陣mag、phase包含系統(tǒng)頻率響應(yīng)的幅值和相角，這些幅值和相角是在用戶指定的頻率點(diǎn)上計算得到的。用戶如果不指定頻率w,MATLAB會自動產(chǎn)生w向量，并根據(jù)w向量上各點(diǎn)計算幅值和相角。這時的相角是以度來表示的，幅值為增益值，在畫伯德圖時要轉(zhuǎn)換成分貝值，因為分貝是作幅頻圖時常用單位?？梢杂梢韵旅畎逊缔D(zhuǎn)變成分貝：magdb=20log10(mag) 繪圖時的橫坐標(biāo)是以對數(shù)分度的。為了指定頻率的范圍，可采用以下命令格式：logspace(d1，d2) 或logspace(d1，d2，n) 公式logspace(d1，d2) 是在指定頻率范圍內(nèi)按對數(shù)距離分成50等分的，即在兩個十

6、進(jìn)制數(shù)和之間產(chǎn)生一個由50個點(diǎn)組成的分量，向量中的點(diǎn)數(shù)50是一個默認(rèn)值。例如要在弧度/秒與弧度/秒之間的頻區(qū)畫伯德圖，則輸入命令時，,在此頻區(qū)自動按對數(shù)距離等分成50個頻率點(diǎn)，返回到工作空間中，即w=logspace(-1，2)要對計算點(diǎn)數(shù)進(jìn)行人工設(shè)定，則采用公式logspace(d1，d2，n)。例如，要在與之間產(chǎn)生100個對數(shù)等分點(diǎn)，可輸入以下命令：w=logspace(0，3，100)在畫伯德圖時，利用以上各式產(chǎn)生的頻率向量w，可以很方便地畫出希望頻率的伯德圖。由于伯德圖是半對數(shù)坐標(biāo)圖且幅頻圖和相頻圖要同時在一個繪圖窗口中繪制，因此，要用到半對數(shù)坐標(biāo)繪圖函數(shù)和子圖命令。1）對數(shù)坐標(biāo)繪圖函

7、數(shù)利用工作空間中的向量x，y繪圖，要調(diào)用plot函數(shù)，若要繪制對數(shù)或半對數(shù)坐標(biāo)圖，只需要用相應(yīng)函數(shù)名取代plot即可，其余參數(shù)應(yīng)用與plot完全一致。命令公式有：semilogx(x，y，s) 上式表示只對x軸進(jìn)行對數(shù)變換，y軸仍為線性坐標(biāo)。semilogy(x，y，s) 上式是y軸取對數(shù)變換的半對數(shù)坐標(biāo)圖。Loglog(x，y，s) 上式是全對數(shù)坐標(biāo)圖，即x軸和y 軸均取對數(shù)變換。2）子圖命令MATLAB允許將一個圖形窗口分成多個子窗口，分別顯示多個圖形，這就要用到subplot()函數(shù)，其調(diào)用格式為：subplot(m，n，k)該函數(shù)將把一個圖形窗口分割成mn個子繪圖區(qū)域，m為行數(shù)，n為列

8、數(shù)，用戶可以通過參數(shù)k調(diào)用各子繪圖區(qū)域進(jìn)行操作，子圖區(qū)域編號為按行從左至右編號。對一個子圖進(jìn)行的圖形設(shè)置不會影響到其它子圖，而且允許各子圖具有不同的坐標(biāo)系。例如，subplot(4，3，6)則表示將窗口分割成43個部分。在第6部分上繪制圖形。 MATLAB最多允許99的分割。例5.2給定單位負(fù)反饋系統(tǒng)的開環(huán)傳遞函數(shù)為： 試畫出伯德圖。利用以下MATLAB程序，可以直接在屏幕上繪出伯德圖如圖5-2。 num=10*1,1;den=1,7,0;bode(num,den)gridtitle(Bode Diagram of G(s)=10*(s+1)/s(s+7) )該程序繪圖時的頻率范圍是自動確定的

9、，從0.01弧度/秒到30弧度/秒，且幅值取分貝值，軸取對數(shù)，圖形分成2個子圖，均是自動完成的。圖5-2自動產(chǎn)生頻率點(diǎn)畫出的伯德圖如果希望顯示的頻率范圍窄一點(diǎn)，則程序修改為： num=10*1,1;den=1,7,0;w=logspace(-1,2,50); % 從0.1至100，取50個點(diǎn)。mag, phase, w=bode(num, den, w);magdB=20*log10(mag) % 增益值轉(zhuǎn)化為分貝值。% 第一個圖畫伯德圖幅頻部分。subplot(2,1,1);semilogx(w,magdB, -r) % 用紅線畫gridtitle(Bode Diagram of G(s)=

10、 10*(s+1)/s(s+7) )xlabel(Frequency(rad/s)ylabel(Gain(dB)% 第二個圖畫伯德圖相頻部分。subplot(2,1,2); semilogx(w,phase, -r); gridxlabel(Frequency(rad/s)ylabel(Phase(deg) ) 圖5-3 用戶指定的頻率點(diǎn)畫出的伯德圖修改程序后畫出的伯德圖如5-3所示。（4）. 用MATLEB求取穩(wěn)定裕量同前面介紹的求時域響應(yīng)性能指標(biāo)類似，由MATLAB里bode()函數(shù)繪制的伯德圖也可以采用游動鼠標(biāo)法求取系統(tǒng)的幅值裕量和相位裕量。例如，我們可以在圖20的幅頻曲線上按住鼠標(biāo)左鍵

11、游動鼠標(biāo)，找出縱坐標(biāo)（Magnitude）趨近于零的點(diǎn)，從提示框圖中讀出其頻率約為7.25dB。然后在相頻曲線上用同樣的方法找到橫坐標(biāo)（Frequence）最接近7.25dB的點(diǎn)，可讀出其相角為-53.9度，由此可得，此系統(tǒng)的相角裕量為126.1度。幅值裕量的計算方法與此類似。 此外，控制系統(tǒng)工具箱中提供了margin()函數(shù)來求取給定線性系統(tǒng)幅值裕量和相位裕量，該函數(shù)可以由下面格式來調(diào)用：Gm, Pm, Wcg, Wcp=margin(G); (67)可以看出，幅值裕量與相位裕量可以由LTI對象G求出，返回的變量對（Gm, Wcg）為幅值裕量的值與相應(yīng)的相角穿越頻率，而（Pm, Wcp）則為

12、相位裕量的值與相應(yīng)的幅值穿越頻率。若得出的裕量為無窮大，則其值為Inf，這時相應(yīng)的頻率值為NaN(表示非數(shù)值)，Inf和NaN均為MATLAB軟件保留的常數(shù)。如果已知系統(tǒng)的頻率響應(yīng)數(shù)據(jù)，我們還可以由下面的格式調(diào)用此函數(shù)。Gm, Pm, Wcg, Wcp=margin(mag, phase, w);其中(mag, phase, w)分別為頻率響應(yīng)的幅值、相位與頻率向量。例5.3 已知三階系統(tǒng)開環(huán)傳遞函數(shù)為：利用下面的MATLAB程序，畫出系統(tǒng)的奈氏圖，求出相應(yīng)的幅值裕量和相位裕量，并求出閉環(huán)單位階躍響應(yīng)曲線。 G=tf(3.5,1,2,3,2);subplot(1,2,1);% 第一個圖為奈氏圖

13、nyquist(G);gridxlabel(Real Axis)ylabel(Imag Axis)% 第二個圖為時域響應(yīng)圖Gm,Pm,Wcg,Wcp=margin(G)G_c=feedback(G,1);subplot(1,2,2);step(G_c)gridxlabel(Time(secs) )ylabel(Amplitude) 顯示結(jié)果為：ans=1.1429 1.1578 1.7321 1.6542圖5-4三階系統(tǒng)的奈氏圖和階躍響應(yīng)圖畫出的圖形如圖5-4 所示。由奈氏曲線可以看出，奈氏曲線并不包圍（-1，j0）點(diǎn)，故閉環(huán)系統(tǒng)是穩(wěn)定的。由于幅值裕量雖然大于1，但很接近1，故奈氏曲線與實軸的

14、交點(diǎn)離臨界點(diǎn)（-1，j0）很近，且相位裕量也只有7.1578o，所以系統(tǒng)盡管穩(wěn)定，但其性能不會太好。觀察閉環(huán)階躍響應(yīng)圖，可以看到波形有較強(qiáng)的振蕩。如果系統(tǒng)的相角裕量45o,我們一般稱該系統(tǒng)有較好的相角裕量。例5.4 考慮一個新的系統(tǒng)模型，開環(huán)傳遞函數(shù)為：由下面MATLAB語句可直接求出系統(tǒng)的幅值裕量和相位裕量: G=tf(100*conv(1,5,1,5), conv(1,1,1,1,9);Gm, Pm, Wcg, Wcp=margin(G)結(jié)果顯示 Gm =Pm = Inf85.4365Wcg =Wcp =NaN100.3285再輸入命令 G_c=feedback(G,1); step(G_

15、c)gridxlalel(Time(sec) )ylalel(Amplitude)圖5-5較理想的系統(tǒng)響應(yīng)可以看出，該系統(tǒng)有無窮大幅值裕量，且相角裕量高達(dá)85.4365o。所以系統(tǒng)的閉環(huán)響應(yīng)是較理想的，閉環(huán)響應(yīng)圖如圖5-5.（5）.時間延遲系統(tǒng)的頻域響應(yīng)(1) 時間延遲系統(tǒng)的傳遞函數(shù)模型帶有延遲環(huán)節(jié)e-Ts的系統(tǒng)不具有有理函數(shù)的標(biāo)準(zhǔn)形式，在MATLAB中，建立這類系統(tǒng)的模型。要由一個屬性設(shè)置函數(shù)set()來實現(xiàn)。該函數(shù)的調(diào)用格式為：set(H, 屬性名, 屬性值)其中H為圖形元素的句柄（handle）。在MATLAB中，當(dāng)對圖形元素作進(jìn)一步操作時，只需對該句柄進(jìn)行操作即可。例如以下調(diào)用格式h=

16、plot(x,y)G=tf(num,den)Plot()函數(shù)將返回一個句柄h，tf()函數(shù)返回一個句柄G，要想改變句柄h所對應(yīng)曲線的顏色，則可以調(diào)用下面命令：Set(h,color,1,0,0);即對“color”參數(shù)進(jìn)行賦值，將曲線變成紅色（由1,0,0決定）同樣，要想對G句柄所對應(yīng)模型的延遲時間Td進(jìn)行修改，則可調(diào)用下面命令Set(G, Td,T)其中T為延遲時間。由此修改后，模型G就已具有時間延遲特性。(2) 時間延遲系統(tǒng)的頻域響應(yīng)含有一個延遲環(huán)節(jié)的系統(tǒng)，其開環(huán)頻域響應(yīng)為 可見，該系統(tǒng)的幅頻特性不變，只加大了相位滯后。例5.5 考慮系統(tǒng)的開環(huán)模型為: 當(dāng)T=1時，我們可以由下面的MATL

17、AB命令繪出系統(tǒng)的奈氏圖，如圖5-6所示,此系統(tǒng)對應(yīng)的時域響應(yīng)圖為5-7。 G=tf(1,1, 1); T=1;w=0, logspace(-3, 1, 100), logspace(1,2,200);set(G,Td, T); % 延遲1秒。nyquist(G,w)gridfigure% 建立一個新的繪圖窗口step(G)圖5-6 時間延遲系統(tǒng)奈氏圖圖5-7 時間延遲系統(tǒng)的階躍響應(yīng)例5.6 控制系統(tǒng)如圖5-8所示，試分析其閉環(huán)特性。圖5-8 下面是該系統(tǒng)的演示過程，帶符號語句為注釋語句。 %清文字屏 %clc %清圖形屏 %clf %系統(tǒng)1描述 a1=1 3 0; b1=2; %系統(tǒng)2描述

18、a2=1 2 2; b2=1 2; %兩系統(tǒng)串聯(lián)后的傳遞函數(shù)： b12,a12=series(b1,a1,b2,a2); %串聯(lián)后頻率響應(yīng)（波特圖） figure(1) bode(b12,a12);pause %clf %求振幅及相位裕度 figure(2) m,p,w=bode(b12,a12) %求開環(huán)振幅Gm和相位裕度Pm及其對應(yīng)頻率 margin(m,p,w); Gm,Pm,Wg,Wp=margin(m,p,w);pause %畫nichols曲線，先用ngrid函數(shù)畫坐標(biāo)網(wǎng) figure(3) ngrid(new); nichols(b12,a12);pause %clf %畫根軌跡

19、，先消除畫nichols曲線時凍結(jié)的坐標(biāo)系， %用axis命令恢復(fù)自動定標(biāo)。 figure(4) axis; rlocus(b12,a12); %用鼠標(biāo)器找適當(dāng)?shù)臉O點(diǎn) k1,p1= rlocfind(b12,a12);pause %將環(huán)路閉合 b1=k1*b1; b121=k1*b12; %再看閉環(huán)振幅及相位裕度 figure(5) m,p,w=bode(b121,a12); margin(m,p,w); Gm,Pm,Wg,Wp=margin(m,p,w);pause %求閉環(huán)傳遞函數(shù)， bb，ab指單位反饋情況， % bb1，ab1指將系統(tǒng)2置于反饋通道上： bb,ab=cloop(b121

20、,a12);pause bb1,ab1=feedback(b1,a1,b2,a2,-1);pause %求閉環(huán)頻率響應(yīng) figure(6) bode(bb1,ab1);pause %clg; %閉環(huán)過渡過程 figure(7) step(bb1,ab1);pause %clg figure(8) impulse(bb1,ab1);pause %clg; %再求閉環(huán)零一極點(diǎn)，作為校核： zb,pb,kb=tf2zp(bb,ab);pause printsys(zb,pb,kb); format echo off %演示結(jié)束運(yùn)行結(jié)果：開環(huán)系統(tǒng)幅值裕量及相位裕量：Gm=3.5435,Pm=55.02

21、2,Wg=1.6159,Wp=0.66898在根軌跡上選擇一點(diǎn)：selected_point = -4.0158 + 0.0309i此時閉環(huán)系統(tǒng)幅值裕量及相位裕量：Gm= 0.34677,Pm= -29.162,Wg= 1.6159,Wp= 2.5192 根據(jù)選定根軌跡上一點(diǎn)的位置可計算出相應(yīng)的閉環(huán)零點(diǎn)，閉環(huán)極點(diǎn)及閉環(huán)增益值。二、附錄 基于MATLAB的線性系統(tǒng)的頻域分析的其它實例例5.7 已知單位負(fù)反饋系統(tǒng)前向通道的傳遞函數(shù)為： 試?yán)L制出Bode圖并計算系統(tǒng)的頻域性能指標(biāo)。程序及結(jié)果如下：num=0 0 2 8 12 8 2;den=1 5 10 10 5 1 0;sys=tf(num,de

22、n);mag,phase,w=bode(sys);gm,pm,wcp,wcg=margin(mag,phase,w);margin(mag,phase,w)kg=20*log(gm)結(jié)果：gm=332.17, pm=38.692, wg=25.847, wc=1.249kg=50.4272 例5.8 已知一帶延遲因子的系統(tǒng)開環(huán)傳函為：，試求其有理傳遞函數(shù)的頻率響應(yīng)，同時在同一坐標(biāo)中繪制以Pade近似延遲因子式系統(tǒng)的Bode圖，并求此時系統(tǒng)的頻域性能指標(biāo)。程序及結(jié)果如下：n1=10 50;d1=conv(1 3,1 3);s1=tf(n1,d1);w=logspace(-1,2);mag1,phase1=bode(n1,d1,w);np,dp=pade(0.8,4);%s1=tf(n1,d1)*tf(np,dp);numt=conv(n1,np);dent=conv(d1,dp);s2=