【正文】 [mag,phase]=dbode(G,H,C,D,Ts,iu,w) ? 式中， num,den和 G,H,C,D分別為離散系統(tǒng)傳遞函數(shù)和狀態(tài)方程模型的參數(shù)。H=[2。由此可見(jiàn)，帶有時(shí)間延遲的系統(tǒng)從某種意義上說(shuō)，相當(dāng)于在一個(gè)不帶時(shí)間延遲的傳遞函數(shù)模型后面串接一個(gè)純時(shí)間延遲環(huán)節(jié) 。一般情況下，取 pade近似的擬合階次為 3或 4就可以獲得相當(dāng)滿意的精度。*180/pi。) ? subplot(2,1,2)。返回的 num和 den分別為辨識(shí)出傳遞函數(shù)分子和分母的系數(shù)向量，即系統(tǒng)的傳遞函數(shù)模型。num=[1 7 24 24]。而系統(tǒng)的能觀測(cè)性是指系統(tǒng)狀態(tài)的變化能否由系統(tǒng)的輸出反映出來(lái)。 ?00( ) ( ) ( )fTj t j tG j g t e dt g t e dt??? ? ??????? 【 例 5－ 26】 利用例 5－ 18在給定頻率范圍上求得的頻率響應(yīng)值 F，辨識(shí)系統(tǒng)的傳遞函數(shù)模型。這種獲得數(shù)學(xué)模型的過(guò)程稱為系統(tǒng)辨識(shí)。 ? semilogx(w,20*log10(m1),w,20*log10(m2),39。 ? [m1,p1]=bode(num1,den1,w)。 A, B, C, D1()Gs1()Gs 1()Gs??11( ) ( ) , ( ) ( )G jw G jw G jw G jw ??? ? ? ? ?se??se??? 純時(shí)間延遲環(huán)節(jié) 常采用 近似方法，它是 1892年由法國(guó)數(shù)學(xué)家 提出的一種著名的有理近似方法，其表達(dá)式為： ? 式中， 稱為 近似系數(shù)。grid ? 如圖 5－ 21所示。0 1 1。 plot(phase,20*log10(mag)) ngrid 執(zhí)行后可得如圖 520所示的 Nichols圖 。由圖可知， Nyquist曲線不包圍 (,j0)點(diǎn)，而開(kāi)環(huán)系統(tǒng)的三個(gè)極點(diǎn)均位于左半 s平面，因此這個(gè)系統(tǒng)是穩(wěn)定的。 ? [num1,den1]=zp2tf(z1,p1,k1)。 ? [num,den]=zp2tf(z,p,k)。 ? 例 526 已知系統(tǒng)的開(kāi)環(huán)傳遞函數(shù)為 ? 繪制 Nyquist圖 ， 并判斷系統(tǒng)的穩(wěn)定性 。den=[1 10 35 50 24]。 例 【 5－ 18】 已知系統(tǒng)的傳遞函數(shù)為 求在頻率為 ~10之間的頻率響應(yīng)。 ? disp([‘幅值裕量 =’,num2str(20*log10(Gm)),‘dB,’,‘相位裕量=’,num2str(Pm),‘度 。 54]。 利用 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 o M a pR e a l A x i sImaginary Axis圖 58 例 5－ 11的曲線 ?6. 任意輸入函數(shù)的響應(yīng) ? 連續(xù)系統(tǒng)對(duì)任意輸入函數(shù)的響應(yīng)可利用 MATLAB的函數(shù) lsim( )求取 ， 其調(diào)用格式為 ? [y,x]=lsim(num,den, u,t) ? 或 [y,x]=lsim(A,B,C,D,iu,u,t) ? 其中 u為給定輸入序列構(gòu)成的矩陣 ， 它的每列對(duì)應(yīng)一個(gè)輸入 ， 每行對(duì)應(yīng)一個(gè)新的時(shí)間點(diǎn) ， 其行數(shù)與時(shí)間 t的長(zhǎng)度相等 。 ? lsim(num,den,u,t)。den=[1 ]。一般地，將這一參數(shù)選作開(kāi)環(huán)系統(tǒng)的增益 K，而在無(wú)零極點(diǎn)對(duì)消時(shí)，閉環(huán)系統(tǒng)特征方程的根就是閉環(huán)傳遞函數(shù)的極點(diǎn)。 0)()(1 ?? sHsG0)( )(1 ?? sd en sn u mK+G ( s )H ( s )圖7 1 2圖 5－ 11 ? rlocus函數(shù)既適用于連續(xù)系統(tǒng)，也適用于離散系統(tǒng)。 ? rlocus(num,den), ? [K,poles]=rlocfind(num,den)。在已有的圖形上繪制柵格線（’ ） 先清除圖形屏幕，然后繪制出柵格線，并設(shè)置成 使后續(xù)的繪圖命令能繪制在柵格上；（ ） 可指定阻尼系數(shù) 和自然頻率并且在繪制柵格線之前清除圖形窗口。sgrid。 （ Bode圖） 設(shè)系統(tǒng)的開(kāi)環(huán)傳遞函數(shù)為： 10111...()...mmmnnnb s b s bGss a s a??? ? ??? ? ?則系統(tǒng)的頻率響應(yīng)為： 系統(tǒng)的伯德圖就是幅值 與相位 分別對(duì)角頻率 進(jìn)行繪圖。 ? bode函數(shù)本身可以通過(guò)輸入元素的個(gè)數(shù)來(lái)自動(dòng)地識(shí)別給出的是傳遞函數(shù)模型還是狀態(tài)方程模型。 figure(1)。 3 0 2 0 1 0010Magnitude (dB)100101 1 8 0 1 3 5 9 0 4 50Phase (