【正文】 Pm = ? ? Wcg = ? ? Wcp = ? ? 幅值裕量 =，相位裕量 =，即 ? ω＝ ，相位曲線穿越－ 180 線，幅值裕量為 Gm==。 ? F=freqresp(num,den,sqrt(1)*w) ? 結(jié)果顯示： 324 3 27 2 4 2 4()1 0 3 5 5 0 2 4s s sGss s s s? ? ??? ? ? ?? F = ? ? ? ? ? ? ? ? MATLAB也提供了直接求取 LTI系統(tǒng)（ LTI系統(tǒng)的定義參見(jiàn)第九章）在單個(gè)復(fù)數(shù)頻率點(diǎn)處的頻率響應(yīng)數(shù)據(jù)的函數(shù) evalfr()，其調(diào)用格式為： ? F= evalfr(sys,w) ? 式中， F為頻率響應(yīng)， w為給定的頻率向量。 ? Nyquist穩(wěn)定性判據(jù)可表示為：當(dāng) ω從－ ∞ → + ∞變化時(shí) ， Nyquist曲線 G(jω)H(jω)逆時(shí)針包圍 (1,j0)點(diǎn)的次數(shù) N， 等于系統(tǒng)開(kāi)環(huán)傳遞函數(shù)G(s)H(s)位于右半 s平面的極點(diǎn)數(shù) P， 即 N＝ P，則閉環(huán)系統(tǒng)穩(wěn)定 ， 否則閉環(huán)系統(tǒng)不穩(wěn)定 。 )()(23 ???? ssssHsG圖 516 例 5－ 20的 Nyquist曲線 ? 例 【 5－ 21】 已知系統(tǒng)的開(kāi)環(huán)傳遞函數(shù)為： ? 繪制系統(tǒng)的 Nyquist曲線，并判別閉環(huán)系統(tǒng)的穩(wěn)定性。 ? ? ? ? ? ?1 6 .7()0 .8 1 0 .2 5 1 0 .0 6 2 5 1sGss s s? ? ? ? G (s) 10 R (s) Y (s) 圖 5－ 18 例 5－ 22的圖 ? 解：先算出內(nèi)環(huán)傳遞函數(shù)， ? 然后以 G0(s)為開(kāi)環(huán)傳遞函數(shù)繪制出 Nyquist曲線，但這里不能直接采用奈氏判據(jù)，因?yàn)樵谇跋蛲ǖ郎嫌幸环糯笙禂?shù) k=10,因此，奈氏判據(jù)中的臨界點(diǎn)應(yīng)改成 (1/k,j0)點(diǎn)。[num1,den1]=zp2tf(z,p,10*k)。 解 MATLAB程序?yàn)? % num=1。當(dāng) w省略時(shí)，頻率點(diǎn)在 0~ 弧度之間自動(dòng)選取。1]。由此可見(jiàn)，帶有時(shí)間延遲的系統(tǒng)從某種意義上說(shuō)，相當(dāng)于在一個(gè)不帶時(shí)間延遲的傳遞函數(shù)模型后面串接一個(gè)純時(shí)間延遲環(huán)節(jié) 。 ? 解： MATLAB程序?yàn)椋? ? % ? num1=[1 1]。 ? numT=conv(n2,num1)。39。 ? [num,den]=invfreqs(mag.*exp(sqrt(1)*phase),w,m,n) ? 除了由頻率響應(yīng)數(shù)據(jù)來(lái)辨識(shí)系統(tǒng)模型外，還可以根據(jù)階躍響應(yīng)及脈沖響應(yīng)數(shù)據(jù)對(duì)系統(tǒng)的傳遞函數(shù)進(jìn)行辨識(shí)，其具體的做法是，首先將階躍響應(yīng)數(shù)據(jù)或脈沖響應(yīng)數(shù)據(jù)轉(zhuǎn)化為響應(yīng)的頻率響應(yīng)數(shù)據(jù)，然后再根據(jù)上面的方法來(lái)辨識(shí)原系統(tǒng)的模型。 ? F=freqresp(num,den,sqrt(1)*w) ? [num,den]=invfreqs(F,w,3,4),printsys(num,den) ? 結(jié)果顯示： ? num = ? ? den = ? ? ? num/den = ? ? 1 s^3 + 7 s^2 + 24 s + 24 ? ? s^4 + 10 s^3 + 35 s^2 + 50 s + 24 ? 利用 MATLAB的控制系統(tǒng)工具箱，不僅可以分析系統(tǒng)的能控性和能觀測(cè)性。 ? 。num=[1 7 24 24]。返回的 num和 den分別為辨識(shí)出傳遞函數(shù)分子和分母的系數(shù)向量，即系統(tǒng)的傳遞函數(shù)模型。) ? subplot(2,1,2)。*180/pi。一般情況下，取 pade近似的擬合階次為 3或 4就可以獲得相當(dāng)滿意的精度。由此可見(jiàn)，帶有時(shí)間延遲的系統(tǒng)從某種意義上說(shuō)，相當(dāng)于在一個(gè)不帶時(shí)間延遲的傳遞函數(shù)模型后面串接一個(gè)純時(shí)間延遲環(huán)節(jié) 。H=[2。 例如 ， dbode()函數(shù)的調(diào)用格式為： [mag,phase]=dbode(num,den,Ts,w) ? 或 [mag,phase]=dbode(G,H,C,D,Ts,iu,w) ? 式中， num,den和 G,H,C,D分別為離散系統(tǒng)傳遞函數(shù)和狀態(tài)方程模型的參數(shù)。但 Nichols圖的繪制方式和 Bode圖是不同的 ， 它可由以下命令繪制 ? plot(phase,20*log10(mag)) ? 當(dāng)然 ， Nichols圖也可采用與 Bode圖類似的簡(jiǎn)單命令來(lái)直接繪制 。nyquist(num,den)。由圖可知， ? Nyquist曲線按逆時(shí)針?lè)较虬鼑?(1,j0)點(diǎn) 1次，而開(kāi)環(huán)系統(tǒng)包含右半 s平面上的一個(gè)極點(diǎn)，所以以此構(gòu)成的閉環(huán)系統(tǒng)穩(wěn)定。den=[1 2 1 ]。F=evalfr(sys,w)。num=[1 7 24 24]。])。0。 利用 MATLAB控制系統(tǒng)工具箱提供的 margin( )函數(shù)可以求出系統(tǒng)的幅值裕量與相位裕量 ， 該函數(shù)的調(diào)用格式為 ? [Gm,Pm,Wcg,Wcp]=margin(num,den) ? 或 [Gm,Pm,Wcg,Wcp]=margin(A,B,C,D) ? 式中 Gm和 Pm分別為系統(tǒng)的幅值裕量和相位裕量 ， 而 Wcg和 Wcp分別為幅值裕量和相位裕量處相應(yīng)的頻率值 ， 即相位穿越頻率 和幅值穿越 。num=9。當(dāng)帶輸出變量引用函數(shù)時(shí)，可得到系統(tǒng) Bode圖相應(yīng)的幅值 mag，相位 phase及頻率點(diǎn)向量 ，其 ( , , , )A B C D???? 相互關(guān)系為 ? 相位以度為單位，幅值可轉(zhuǎn)換成以分貝為單位，即 1( ) ( )G j j?? ?? C I A B + D( ) ( ) , ( ) ( )m a g G j p h a s e G j? ? ? ?? ? ?( ) 2 0 * l o g 1 0 ( )m a g d B m a g?? 有 了 這 些 數(shù) 據(jù) 就 可 以 利 用 下 面 的MATLAB命令 ? subplot(2,1,1)。因此，也稱為幅頻和相頻曲線。 ? 執(zhí)行后可得如圖 5－ 13所示的仿真結(jié)果。（）函數(shù)可在連續(xù)系統(tǒng)的根軌跡或零極點(diǎn)圖上繪制出柵格線，柵格線由等阻尼系數(shù)和等自然頻率線構(gòu)成，阻尼系數(shù) 的步長(zhǎng)為 ；自然頻率 的步長(zhǎng)為 范圍為 。 )2)(1()()(???sssKsHsG 執(zhí)行以上程序，并移動(dòng)鼠標(biāo)到根軌跡與虛軸的交點(diǎn)處單擊鼠標(biāo)左鍵后可得如圖 512所示的根軌跡和如下結(jié)果： selected a point in the graphics window ? selected_point = ? ? K = ? ? poles = ? ? + ? 圖 512根軌跡圖 ? 由此可見(jiàn)根軌跡與虛軸交點(diǎn)處的增益K=6，這說(shuō)明當(dāng) K 6時(shí)系統(tǒng)穩(wěn)定，當(dāng) K 6時(shí)，系統(tǒng)不穩(wěn)定；利用 rlocfind( )函數(shù)也可找出根軌跡從實(shí)軸上的分離點(diǎn)處的增益 K =, 這說(shuō)明當(dāng) 0 K ，系統(tǒng)為單調(diào)衰減穩(wěn)定，當(dāng) K 6時(shí)系統(tǒng)為振蕩衰減穩(wěn)定的。 rlocus（ num,den）或 rlocus（ A,B,C,D）繪制系統(tǒng)根軌跡時(shí)，增益 K是自動(dòng)選取的， rlocus（ num,den,K）或 rlocus（ A,B,C,D,K）可利用指定的增益 K來(lái)繪制系統(tǒng)的根軌跡。 ? 關(guān)于控制系統(tǒng)的根軌跡分析，在 MATLAB控制系統(tǒng)工具箱中提供了幾個(gè)函數(shù)，如表 5－ 2所示。 ? u1=[ones(1,50),1*ones(1,50)]。title(39。 其他用法同 step( )函數(shù) 。title(39。dimpulse(A,B,C,D),title(39。[A,B,C,D]=c2dm(A1,B1,C1,D1,T,39。0 0 0]。D=zeros(2,2)。 1。 解： MATLAB程序?yàn)? % num=[2 ]。hold on ? for wn = w ? num=wn.^2。Step Response39。 2222)(nnnsssG???????解 MATLAB程序?yàn)? % wn=6。[num,den]=cloop(num0,den0)。 ? 解 Matlab窗口中執(zhí)行以下命令可得圖 52所示結(jié)果 。)。P=lyap(A,Q)。在這種情況下，利用李雅普諾夫第二法比較有效，尤其在系統(tǒng)含有非線性環(huán)節(jié)時(shí)更是如此。)。r=roots(P)。當(dāng)然判斷系統(tǒng)的穩(wěn)定性同樣可利用特征值來(lái)判斷。 unstable pole39。den0=[1 0 0 0 0 0]。)。 ? [例 51] 已知閉環(huán)系統(tǒng)的傳遞函數(shù)為 ? 判斷系統(tǒng)的穩(wěn)定性，并給出不穩(wěn)定極點(diǎn)。這樣做不但會(huì)大大提高運(yùn)算的效率，而且可以提高仿真的精度和可靠性。其實(shí)，對(duì)于各種線性系統(tǒng)模型在典型輸入信號(hào)作用下來(lái)說(shuō)，當(dāng)然沒(méi)有必要采用那些通用的算法來(lái)完成這種任務(wù)，而是應(yīng)該充分地利用線性系統(tǒng)的特點(diǎn)，采取更簡(jiǎn)單的方法來(lái)得到問(wèn)題的解。 ? ? 判斷一個(gè)線性系統(tǒng)穩(wěn)定性的一種最有效的方法是直接求出系統(tǒng)所有的極點(diǎn)，然后根據(jù)極點(diǎn)的分布情況來(lái)確定系統(tǒng)的穩(wěn)定性，對(duì)于極點(diǎn)的求取我們?cè)谏瞎?jié)中已作過(guò)介紹，下面舉例說(shuō)明其判斷方法。The Unstable Poles are:39。 解：則可利用下面的 MATLAB程序： % num0=[5 4 1 3 ]。,int2str(n1),39。 令特征多項(xiàng)式等于零，即得系統(tǒng)的特征方程 |sIA|=sn+a1sn1+…+ an1s+an＝ 0 的根稱為系統(tǒng)的特征值，即系統(tǒng)的閉環(huán)極點(diǎn)。 P=poly(A)。System is Stable39。end ? 在高階系統(tǒng)或者特征多項(xiàng)式中，當(dāng)某些系數(shù)不是數(shù)值時(shí)，利用求閉環(huán)極點(diǎn)或特征值的方法來(lái)判斷系統(tǒng)的穩(wěn)定性是比較困難的。Q=eye(size(A))。P0,正定，系統(tǒng)在原點(diǎn)處的平衡狀態(tài)是漸進(jìn)穩(wěn)定的 39。 ? 例 55 生成一個(gè)周期為 5秒 ， 持續(xù)時(shí)間為30秒 ， 采樣時(shí)間為 。den0=[1 8 36 40 0]。 S t e p R e s p o n s eT i m e ( s e c )Amplitude0 1 2 3 4 5 6 7 8 9 1000 . 20 . 40 . 60 . 811 . 21 . 4S y s t e m : s y sP e a k a m p l i t u d e : 1 . 0 3O v e r s h o o t ( % ) : 2 . 5 5A t t i m e ( s e c ) : 5 . 8圖 53 例 5－ 6的單位階躍響應(yīng)曲線 ? 例 57 對(duì)于典型二階系統(tǒng) ? ? 試?yán)L制出無(wú)阻尼自然振蕩頻率 ωn=6， 阻尼比 ζ分別為 ,… ,位階躍響應(yīng)曲線 。end title(39。 ? figure(1)。 例 59 已知二階離散系統(tǒng) 試求其單位階躍響應(yīng) 。 。0 2 0 2]。 。 ? T=。) ? subplot(2,2,2)。dinitial(A,B,C,D,x0) ? axis([0 6 ])。) 0 5 10 15 20 2500 . 511 . 52D i s c r e t e s t e p r e s p o n s eT i m e ( s e c )Amplitude0 5 10 15 0 . 200 . 20 . 40 . 6D i s c r e t e i m p u l s e r e s p o n s eT i m e ( s e c )Amplitude0 2 4 6012D i s c r e t e I n i t i a l R e s p o n s eT i m e ( s e c )Amplitude1 0 . 5 0 0 . 5 1 0 . 2 0 . 100 . 10 . 2D i s c r e t e P o l e Z e r