1、.wd.wd53/53.wd?控制工程根基經典控制局部?的MATLAB分析第一章MATLAB的 根本使用1-1MATLAB語言簡介MATLAB是一種高級矩陣語言，它由Math Works公司于1984年正式推出，它的 根本處理對象是矩陣，即使是一個標量純數，MATLAB也認為它是只有一個元素的矩陣。隨著MATLAB的開展，特別是它所包含的大量工具箱應用程序集的集結，使MATLAB已經成為帶有獨特數據構造、輸入輸出、流程控制語句和函數、并且面向對象的高級語言。MATLAB語言被稱為一種“演算紙式的科學計算語言，它在數值計算、符號運算、數據處理、自動控制、信號處理、神經網絡、優化計算、模糊邏輯、系

2、統辨識、小波分析、圖象處理、統計分析、甚至于金融財會等廣闊領域有著十分廣泛的用途。MATLAB語言在工程計算與分析方面具有無可比較的優異性能。它集計算、數據可視化和程序設計于一體，并能將問題和解決方案以使用者所熟悉的數學符號或圖形表示出來。MATLAB語言和C語言的關系與C語言和匯編語言的關系類似。例如當我們需要求一個矩陣的特征值時,在MATLAB下只需由幾個字符組成的一條指令即可得出結果，而不必去考慮用什么算法以及若何實現這些算法等低級問題，也不必深入了解相應算法的具體內容。就象在C語言下不必象匯編語言中去探究乘法是若何實現的，而只需要采用乘積的結果就可以了。MATLAB語言還有一個巨大的優

3、點是其高度的可靠性。例如對于一個病態矩陣的處理，MATLAB不會得出錯誤的結果，而用C或其它高級語言編寫出來的程序可能會得出錯誤的結果。這是因為MATLAB函數集及其工具箱都是由一些在該領域卓有研究成果，造詣很深的權威學者經過反復比較所得出來的最優方法，而且經過多年的實踐檢驗被證明是正確可靠的。1-2MATLAB的工作窗口下面以MATLAB6.1為例介紹。從實用的角度MATLAB的工作窗口包括命令窗口、M文件編輯器窗口、圖形編輯窗口、數學函數庫、應用程序接口及在線窗口。下面首先介紹MATLAB的命令窗口及M文件編輯器。一、命令窗口啟動MATLAB之后，屏幕上自動出現命令窗口MATLAB，它是M

4、ATLAB提供應用戶的操作界面，用戶可以在命令窗口內提示符“之后有的MATLAB版本命令窗口沒有提示符鍵入MATLAB命令，回車即獲得該命令的答案。命令窗口內有File、Edit、View、Web、Window、Help等菜單條。二、M文件編輯窗口M文件是MATLAB語言所特有的文件。用戶可以在M文件編輯窗口內，編寫一段程序，調試，運行并存盤，所保存的用戶程序即是用戶自己的M文件。MATLAB工具箱中大量的應用程序也是以M文件的形式出現的，這些M文件可以翻開來閱讀，甚至修改，但應注意，不可改開工具箱中的M文件！1進入M文件窗口有兩種方法命令窗口 FileNew M-File；命令窗口點擊“Fi

5、le字樣下面的圖標。M文件編輯窗口的標記是“Untitled無標題的。當用戶編寫的程序要存盤時，Untitled作為默認文件名提供應用戶，自然，用戶可以，也應當自己命名。假設用戶不自己命名，那么MATLAB會對Untitled進展編號。2M文件的執行：返回命令窗口，在當前目錄Current Directory內選擇所要運行的M文件的目錄，在命令窗口提示符“后，直接鍵入文件名不加后綴即可運行。注意：1機器默認路徑為一級子目錄MATLAB6 p1work；2MATLAB 6.1以前的版本，運行M文件的方法稍有不同，它必須在File菜單下，翻開“Run Script子菜單，鍵入需要運行的文件路徑及名

6、稱再回車，在這種情況下，work作為根目錄對待，不出現在M文件的路徑之中。本講義的參考程序都是在M文件窗口下編制的。三、在線幫助窗口在命令窗口中鍵入Help(空格) 函數名，可以立即獲得該函數的使用方法。1-3 MATLAB最 根本的矩陣操作作為命令窗口及M文件編輯器的應用實例，介紹幾個最 根本的矩陣運算命令。一、矩陣的輸入在方括號內依次按行鍵入矩陣元素，在一行內的各元素之間用空格或逗號分開，每行之間用分號分開。例如，在命令窗內輸入A=2 2 3;4 5 4;7 8 9 (注意：方括號，分號為矩陣行標記)B=1,3,5;6,-4,2;3,5,1 (逗號與空格功能一樣)A223B135 4546

7、-42 789 351同理：輸入A1246得到行矢量，輸入A22；4；6得到列矢量，于是，當輸入C=A；A1有C=123456789246A1作為矩陣C的最后一行，C和A相比，增加了一行。二、矩陣的轉置矩陣A的轉置用A表示，顯然，A1與A2互為轉置，即A1會得到以2,4，6為元素的列矢量。思考一下輸入C1 AA2 C2 AA1有什么結果而輸入A；A1有無意義三、矩陣的四那么運算1矩陣的加減法：當兩個矩陣維數一樣時可以直接進展“或“運算。如 D1AB，D2AB2 矩陣的乘法：當矩陣A，B維數相容時C3AB：普通意義下的矩陣相乘C4A .B：矩陣A與B的對應元素相乘顯然，ABBA一般情況，而A .

8、BB .A。A .B稱為數列型乘法，它要求參加運算的矩陣或數列具有一樣的行列數，這是MATLAB語言中的一種特殊運算，它在今后求取函數值等運算時是很重要的。實際上，前面所述的矩陣加、減法就是一種數列型運算。3 矩陣的除法D4AB：表示A-1B或inv(A)*B，即A的逆矩陣左乘矩陣B。D5B/A：表示BA-1或Binv(A)，即A的逆矩陣右乘B。D6A .B：表示B的每一個元素被A的對應元素除。D7A ./B：表示A的每一個元素被B的對應元素除。顯然，A .B與A ./B的各對應元素互為倒數。讀者可以思考一下，D6 .D7等于什么D8inv(A)：A的逆矩陣。翻開M文件編輯窗口，將上述命令依次

9、鍵入，得到fanli001如下：參考程序fanli001:矩陣的四那么運算A=2 2 3;4 5 4;7 8 9 % 三階矩陣輸入B=1,3,5;6,-4,2;3,5,1 % 三階矩陣輸入A1=2 4 6 %行向量A2=2;4;6 %列向量C=A;A1 %矩陣A增加一行C1=A A2%矩陣A增加一列C2=A A1 %矩陣A增加一列D1=A+B %矩陣相加D2=A-B%矩陣相減C3=A*B%矩陣與矩陣相乘C4=A.*B%矩陣的對應元素相乘D3=AB%A的逆左乘BD4=B/A%A的逆右乘BD6=A.B%B的各元素被A的對應元素除D7=A./B%A的各元素被B的對應元素除D8=inv(A) %A的逆

10、矩陣語句后面的%為語句說明符。MATLAB中矩陣運算的其它主要命令可通過在線幫助獲得。1-4MATLAB的符號運算操作一、進入符號運算功能在命令窗口鍵入syms x y z t 此后，即可以使用x，y，z，t等作自變量定義函數。syms x y z t real 規定所定義的變量為實型。二、代數方程求解使用命令solve可以求解代數方程，如求下例方程組的解，命令為x,y=solve(x2+x*y+y-3=0,x2-4*x-2*y+3=0)程序見范例程序fanli002。參考程序fanli002:代數方程求解syms x y%進入符號運算功能f1=x2+x*y+y-3%函數f1f2=x2-4*x

11、-2*y+3 %函數f2x,y=solve(x2+x*y+y-3=0,x2-4*x-2*y+3=0)%求解方程組f1a=simplify(subs(f1) %用求解出的x,y檢驗方程1f2a=simplify(subs(f2) %用求解出的x,y檢驗方程2x=double(x) %將符號變量轉換成浮點數y=double(y) %將符號變量轉換成浮點數f1=subs(f1) %用浮點數x,y檢驗方程1f2=subs(f2) %用浮點數x,y檢驗方程2由solve求出的根是根的符號表達形式，是準確解。命令simplify(f)表示化簡，subs(f1)表示將求解出x，y代回f1中；double是將

12、符號變量轉換成浮點數，是準確解x，y的近似值。這從程序運行f1a=0，f2a=0，f10，f20可以確認。三、符號矩陣運算符號矩陣可以和數值矩陣一樣進展運算，例如：求矩陣特征值eig，求矩陣的逆inv等命令都支持符號運算。設計算其特征值eig(A)，逆矩陣inv(A)，程序見fanli003。參考程序fanli003:符號矩陣的特征值syms t real %定義為實型變量A=sin(t) -cos(t);cos(t) sin(t) %定義矩陣AB1=eig(A) %求矩陣A的特征值B1=simple(B1) %化簡A的特征值表達式B2=inv(A) %求矩陣A的逆矩陣B2=simple(B2

13、) %化簡A的逆矩陣表達式C1=A*B2 %檢驗A的逆矩陣C1=simple(C1) %C1為單位矩陣注意函數的輸入方法，自變量用圓括號括起來。四、微積分運算設函數，那么MATLAB中微積分運算命令為fiff(f)：求函數f對自變量x的一階導數；diff(f，2)：求函數f對自變量x的二階導數；int(f)：求函數f的不定積分例1.1，設試計算其一階，二階導數，積分運算，并作出函數圖象，見范例fanli004。參考程序fanli004:函數的微分與積分syms xf=1/(5+4*cos(x)ezplot(f) %函數f的曲線f1=diff(f) %函數f的一階導數figure,ezplot(

14、f1) %函數f一階導數的曲線f2=diff(f,2) %函數f的二階導數figure,ezplot(f2) %函數f二階導數的曲線g=int(int(f2) %函數f的二階導數f2的二重積分figure,ezplot(g) %函數f2二重積分的曲線e=f-g %二階導數的二重積分與原函數的差e=simple(e)figure,ezplot(e)程序中ezplot(f)：作函數的圖形，x的取值范圍默認值為-2 x 2fanli004的函數曲線fanli004的一階導函數曲線fanli004的二階導函數曲線fanli004的二階導函數的二重積分曲線ezplot(f)是一個很有用的作圖命令，它的其

15、它應用形式，請查在線幫助。細心的讀者會發現，一個函數求二階導數后再對二階導數進展二重積分，其結果與原函數相差一個常數。相當于縱坐標發生平移。第二章系統的時域特性2-1傳遞函數一、傳遞函數的兩種形式傳遞函數通常表達成s的有理分式形式及零極點增益形式。設傳遞函數1有理分式形式分別將分子、分母中s多項式的系數按降冪排列成行矢量，缺項的系數用0補齊。上述函數可表示為num1=2 1%注意：方括號，同一行的各元素間留空格或逗號。den1=1 2 2 1syss1=tf(num1,den1)運行后，返回傳遞函數的形式。這種形式不能直接進展符號運算！2零極點增益形式Z,P,K = tf2zp(num1,de

16、n1)sys2 = zpk(Z,P,K)返回零、極點、增益表達式，其Z，P分別將零點和極點表示成列向量，假設無零點或極點用 空矩陣代替。運行得到的Z= -0.5極點P= -1，-0.5j0.866增益K = 2指令zp2tf(Z，P，K)將零極點增益變換成有理分式形式，見程序fanli005。參考程序fanli005:傳遞函數的有理分式及零極點增益模型num1=2 1 % 傳遞函數的分子系數向量den1=1 2 2 1 % 傳遞函數的分母系數向量sys1=tf(num1,den1) % 傳遞函數的有理分式模型Z,P,K=tf2zp(num1,den1) % 有理分式模型轉換成零極點增益模型nu

17、m2,den2=zp2tf(Z,P,K) % 零極點增益模型轉換成有理分式模型sys2=zpk(Z,P,K) % 傳遞函數的零極點增益模型A1,B1,C1,D1=tf2ss(num1,den1) % 有理分式模型轉換成狀態空間模型A2,B2,C2,D2=zp2ss(Z,P,K)% 零極點及增益模型轉換成狀態空間模型num1,den1=ss2tf(A1,B1,C1,D1)% 狀態空間模型轉換成有理分式模型Z,P,K=ss2zp(A2,B2,C2,D2) % 狀態空間模型轉換成零極點增益模型程序中，命令tf2ss，zp2ss及ss2tf，ss2zp是狀態空間模型與有理分式及零、極點、增益模型之間的

18、相互轉換。二、傳遞函數框圖的處理用框圖可以方便地表示傳遞函數的并聯，串聯及反響。為簡潔，僅以有理分式模型為例。G1G1G1+G21 并聯sysp = parallel(sys1,sys2)num,den=parallel(num1,den1,num2,den2)2 串聯G1(s)G2(s)G1(s)G2(s)syss = series(sys1,sys2)nums,dens = series(num1,den1,num2,den2)3 反響G1(s)G2(s)G1(s)G2(s)1+G1(s)G2(s) G3(s)G3(s)sysc=feedback(syss,sys3,1) %默認值-1nu

19、mc, denc = feedback(nums, dens, num3, den3)4 單位反響G1(s)G2(s)G1(s)G2(s)1+G1(s)G2(s) sysd = feedback(syss, 1)numd, dend = feedback(nums, dens, 1, 1)%(單位反響)上面給出了同一指令的兩種形式，相當于兩套平行指令。對于零極點增益形式，書寫稍復雜一些，可先用zpk轉換成系統形式，或用zp2tf轉折換成有理分式形式后再進展框圖化簡操作。三、簡單函數的拉普拉斯變換在MATLAB的符號功能中，可以對簡單函數進展拉普拉斯正、逆變換。拉氏正變換：laplace(f(t

20、)拉氏逆變換：ilaplace(L(s)其中為原函數，為象函數。命令格式參見fanli007。參考程序fanli007:拉普拉斯變換syms s t w a b cf1=sqrt(b-a)2+w2)/w*exp(-a*t)*sin(w*t+atan(w/(b-a)%原函數f1L1=laplace(f1) %f1的拉氏變換象函數L1=simple(L1) %化簡f2=ilaplace(L1)%L1的拉氏逆變換f2=simple(f2) %化簡在MATLAB中使用laplace及ilaplace命令時，要注意象、原函數的符號，特別是對初相不等于零的振蕩系統，運行結果常常同手冊上的結果相差一個符號，

21、這要注意函數表達式成立的條件。保險的方法是再使用拉氏變換的初值定理確定象、原函數的符號。2-2系統時域特性曲線在MATLAB中，當傳遞函數時，可以方便地求出系統的單位脈沖響應、單位階躍響應等曲線。一、系統的單位階躍響應stepstep有以下幾種格式step(sys)：直接作出sys的單位階躍響應曲線。其中sys = tf(num, den)或sys = zpk(z, p, k)，MATLAB自動決定響應時間。step(sys, t)設定響應時間的單位階躍響應。t可以設定為最大響應時間t = t終值(秒)，也可以設置為一個向量t = 0 : t: t終值注意冒號的使用。它產生一個從0到t終值的行

22、矢量，元素之間的間隔為t。step(sys1, sys2, , sysn)在同一幅圖上畫出幾個系統的單位階躍響應。y, t = step(sys)；命令輸出對應時刻t的各個單位階躍響應值，不畫圖。語句后的分號控制數據的屏幕顯示。如果要查看機器計算了多少個數據，可以使用命令size(y)得出的結果也說明數據作為列矢量的行數。要將計算出的y, t作成曲線，使用一般的作圖命令plot(t, y)plot后面跟的兩個參數橫坐標在前，縱坐標在后。參考程序見fanli008：參考程序fanli008:系統的單位階躍響應num1=4 2den1=2 8 14 11 4sys1=tf(num1,den1) %

23、系統G1(s)num2=2 1den2=1 4 6 7 3sys2=tf(num2,den2) %系統G2(s)y1,t1=step(sys1); %系統G1(s)的單位階躍響應數據y2,t2=step(sys2); %系統G2(s)的單位階躍響應數據step(sys1,sys2) %系統G1(s)、G2(s)的單位階躍響曲線figure,step(sys1,sys2,20) %系統G1(s)、G2(s)在自選時間20秒內的單位階躍響曲線figure,plot(t1,y1) %系統G1(s)的單位階躍響應曲線figure,plot(t2,y2) %系統G2(s)的單位階躍響應曲線fanli00

24、8:step(sys1,sys2,t)單位階躍曲線二、系統的單位脈沖響應impulse命令格式與單位階躍響應step的命令格式完全一樣，只需將語句中的step用impulse代替即可。針對同樣的系統，其單位脈沖響應的參考程序見fanli009。參考程序fanli009:系統的單位脈沖響應num1=4 2den1=2 8 14 11 4sys1=tf(num1,den1)num2=2 1den2=1 4 6 7 3sys2=tf(num2,den2)y1,t1=impulse(sys1); %系統G1(s)的單位脈沖響應數據y2,t2=impulse(sys2); %系統G2(s)的單位脈沖響應

25、數據impulse(sys1,sys2) %系統G1(s)、G2(s)的單位脈沖響應曲線figure,impulse(sys1,sys2,20)%系統G1(s)、G2(s)在自選時間20秒內的單位脈沖響應曲線figure,plot(t1,y1) %系統G1(s)的單位脈沖響應曲線figure,plot(t2,y2) %系統G2(s)的單位脈沖響應曲線hold on,step(sys2) %系統G2(s)的單位階躍和單位脈沖響應曲線fanli009:impulse(sys1,sys2)單位脈沖響應曲線程序的最后一句hold on是當前圖形保護模式。當要將新圖形作在當前圖形上時，必須使用hold

26、on。而figure的含意是另開一個新的圖形窗口，如果不用figure或hold on，那么新的圖形會占用原圖形窗口，始終只保存一個最新的圖形窗口。三、一階系統及二階系統的時域特性一階系統及二階系統是最 根本也是最重要的系統，高階系統總可以視為由假設干個一階和或二階系統組合構成。1 一階系統設增益為1影響系統特性的參數是其時間常數T，T越大，系統慣性越大，響應越慢。參考程序fanli010給出了T0.4, 1.2, 2.0, 2.8, 3.6, 4,4六條單位階躍響應曲線。參考程序fanli010:一階系統的單位階躍響應曲線num=1;i=1;for del=0.1:0.2:1.1 %一階系統

27、時間常數遞增間隔den=4*del 1; %一階系統分母向量step(tf(num,den) %一階系統單位階躍響應曲線hold on, %不同時間常數的一階系統單位階躍響應曲線簇i=i+1;end同理，可以作出對應的單位脈沖響應曲線，參考程序fanli011。參考程序fanli011:一階系統的單位脈沖響應曲線num=1;i=1;fanli010一階系統時間常數對單位階躍響應的影響fanli011一階系統時間常數對單位脈沖響應的影響for del=0.1:0.2:1.1den=4*del 1;impulse(tf(num,den),10) %一階系統單位階脈沖應曲線hold on, %不同時

28、間常數的一階系統單位脈沖響應曲線簇i=i+1;end注意MATLAB中for語句的構造。讀者可以改變不同的增益，看看圖形有何變化。2 二階系統設01設二階系統為二階系統的特征參數為固有頻率及阻尼比。當增大，系統振動頻率加快，振蕩加??；而隨著減小，系統振蕩加劇，振蕩峰鋒利。參考程序fanli012示出了當0.5, =1, 2, 3, 4, 5 rad/s時的間接階躍曲線簇。參考程序fanli012: 不同固有頻率的二階系統的單位階躍響應曲線=0.5i=1;for del=1:1:5; % 二階系統固有頻率遞增間隔num=del2; % 二階系統傳遞函數分子系數向量den=1 del del2;

29、% 不同固有頻率的二階系統分母系數向量step(tf(num,den),6) %二階系統單位階躍響應曲線hold on, %不同固有頻率的二階系統單位階躍響應曲線簇i=i+1;endfanli012二階系統固有頻率對單位階躍響應的影響參考程序fanli013示出了同一二階系統的單位脈沖響應曲線簇。參考程序fanli013: 不同固有頻率的二階系統的單位脈沖響應曲線=0.5i=1;for del=1:1:5;num=del2;den=1 del del2; % 不同固有頻率的二階系統分母系數向量impulse(tf(num,den),6) %二階系統單位脈沖響應曲線hold on, %不同固有頻

30、率的二階系統單位脈沖響應曲線簇i=i+1;endfanli013二階系統固有頻率對單位脈沖響應的影響參考程序fanlio14示出了當1，0.1，0.3，0.5，0.7，0.9的二階系統的單位階躍響應曲線簇。參考程序fanli014: 不同阻尼比的二階系統的單位階躍響應曲線=1i=1;for del=0.1:0.2:0.9; % 二階系統阻尼比遞增間隔 num=1; den=1 2*del 1; % 不同阻尼比的二階系統分母系數向量 step(tf(num,den),30) hold on, %不同阻尼比的二階系統單位階躍響應曲線簇 i=i+1;endfamli014二階系統阻尼比對單位階躍響應

31、的影響參考程序fanli015示出了同一二階系統當取0.1，0.3，0.5，0.7，0.9的單位脈沖響應曲線簇。參考程序fanli015: 不同阻尼比的二階系統的單位脈沖響應曲線=1i=1;for del=0.1:0.2:0.9; num=1; den=1 2*del 1;impulse(tf(num,den),30) hold on, %不同阻尼比的二階系統單位脈沖響應曲線簇 i=i+1;endfanli015二階系統阻尼比對單位脈沖響應的影響四、帶延時環節系統的典型響應設具有純延時環節的傳遞函數為計算這種系統的單位階躍響應不能使用一般方法。首先應使用Pade法對延時環節進展近似展開，num

32、T, dent = pade(T, 5)表示使用5階Pade級數得出的分子、分母系數的向量，然后將系統視為慣性環節與延遲環節的串聯，再求其階躍脈沖響應。設K = 2.5，T1，程序見fanli016。參考程序fanli016: 帶延時環節系統的單位階躍響應K=2.5;T=1; T1=0.5;numt,dent=pade(T,5);%5階Pade法近似延時環節syst=tf(numt,dent); %延時環節的近似傳遞函數num1=K;den1=T1 1;sys1=tf(num1,den1);%慣性環節sys=series(sys1,syst);% 帶延時環節的慣性系統step(sys,6),g

33、rid % 帶延時環節慣性系統的單位階躍響應figure,impulse(sys,6),grid% 帶延時環節慣性系統的單位脈沖響應fanli016帶延遲的一階系統的單位階躍響應fanli016帶延遲的一階系統的單位脈沖響應從響應曲線上可見，具有延時環節的慣性系統當加上單位階躍或單位脈沖輸入前后，其響應初始有振蕩特性，脈沖響應更為明顯。這種振蕩現象使得對這種系統進展校正需要使用特殊的方法。在作圖命令后，加上格線命令grid，機器自動給曲線圖加上格線。2-3響應曲線的動態分析一、MATLAB圖形的編輯在MATLAB中，用于編輯圖形的命令很多，因為可視化正是MATLAB的一種強大而優越的性能。從工

34、程實用角度，僅在圖形窗口進展編輯已經夠用了。在圖形窗口頂部的第二排依次排有以下命令按鈕：Edit Plot：圖形編輯。點擊進入圖形編輯功能。點擊選擇編輯對象。例如，曲線，標題，縱、橫坐標的說明等，可對選中項進展編輯。Insert Text：插入文本。點擊可在圖區內插入文字，并可選擇字體及大小等。Insert Arrow：插入箭頭。點擊可在圖區內畫箭頭，并可對該箭頭進展編輯。Insert Line()：插入直線。點擊可在圖區內插入直線，并可對該直線進展編輯，改變線寬及顏色。Zoom In：曲線放大。激活該按鈕，可以逐次放大曲線，以觀察曲線某些局部的細節，這對觀察變化劇烈的曲線局部很有幫助。Zoo

35、m Out：曲線縮小。激活該按鈕，可以逐次從放大狀態返回到原狀態。Rotate 3D：三維旋轉。激活該按鈕，在圖區內按住鼠標左鍵，拖曳鼠標，原二維曲線可三維空間內任意旋轉。此外，還有平面旋轉、三維動畫、前后放大、平移、飛逸等。例如，將曲線加粗，以利作圖曲線默認寬度為0.5。按Edit Plot選擇曲線雙擊左鍵 出現線寬選擇框line width 選擇線寬 OK。二、對曲線進展數據分析不激活Edit Plot按鈕圖形編輯不使能，讀取曲線上任一點的坐標值：鼠標指向選擇點，單擊左鍵，顯示曲線編號，橫坐標及縱坐標值；研究曲線擬合情況：按Tools Basic Fitting，出現Plot fitt圖形

36、擬合窗口，選擇需要擬合的曲線，點擊擬合方式，即可對所選圖形進展曲線擬合。MATLAB提供了樣條插值、保形插值shape-preserving interpolant及直到10階的多項式擬合。中選擇多項式擬合時，機器會給出曲線及擬合方程。在所給出的擬合方式中，三次樣條插值Cubic spline interpolant及保形插值shape-preserving interpolant效果最好。第三章系統的頻率特性系統的頻率特性就是，這是一個復變函數。3-1MATLAB中的復數及復變量在MATLAB中復數x的輸入方法：x = 2 + 3*i或2 + 3*j也可以s = 2 + 3i或s = 3 +

37、 2j在符號運算功能下，復變函數表示為設a，b都是實數：syms a b realx = ab*I或ab*j*不能省略在MATLAB中與復數有關的命令有abs(x)：求實數的絕對值及復數的模，支持符號運算；angle(x)：求復數的相角，單位為弧度，不支持符號運算。conj(x)：求復數的共軛復數，支持符號運算；imag(x)：求復數的虛部，支持符號運算；real(x)：求復數的實部，支持符號運算。復數的運算見參考程序fanli017。參考程序fanli017:復數運算含符號功能syms a b c d real%符號均為實變量X1=2-3i X2=a+c/b*jX3=c/a-2*d*jS1=

38、X2*X3/X1 %復數的乘除S2=simple(abs(S1) %復數的模S3=simple(real(S1) %復數的實部S4=simple(imag(S1) %復數的虛部S5=simple(conj(S1) %復數的共軛S6=simple(expand(S1*S5) %復數與其共軛的積S7=simple(S32+S42) %復數實部、虛部的平方和S8=simple(S6-S7) %檢驗復數運算的正確性fanli018給出了系統的頻率特性。參考程序fanli018:系統頻率特性的解析表達syms w real %定義頻率w 為實型變量g=50*(0.6*w*i+1)/(i*w)2/(4*w

39、*j+1) %系統頻率特性g1=simple(conj(g) %系統頻率特性的復共軛gr=simple(real(g) %系統實頻特性gi=simple(imag(g) %系統虛頻特性ga=simple(abs(g) %系統幅頻特性模gb=simple(sqrt(gr2+gi2) %系統幅頻特性的另一計算方法gc=simple(ga-gb) %檢驗系統幅頻特性的兩種計算方法3-2頻率特性的Nyquist圖將系統傳遞函數中的復變數s用純虛數j為角頻率，rad/s代替，即得到系統的頻率特性，又稱為諧波傳遞函數，有三種表示方法：式中：實頻特性：虛頻特性：幅頻特性：相頻特性在研究控制系統的特性時，采用

40、圖示方法比采用解析方法直觀、簡單，便于對動態特性進展更深入的研究。以下討論設系統傳遞函數為sys。一、Nyquist圖Nyquist 圖的有關命令nyquist(sys)直接返回系統sys的Nyquist圖，機器自己確定頻率范圍通常為-+nyquist(sys, w)返回在規定頻率rad/s范圍內的Nyquist圖，頻率范圍的格式為w = wmin, wmax (0wmin，那么意味著選擇了超前校正；假設選g，校正步驟為選0.1g為滯后環節零點轉角頻率，即選校正后系統的幅值穿越頻率為待校正系統的相位穿越頻率g，由G(s)的Bode圖可求出g時的幅值，設為AdB。這一幅值，應當由起始于-20dB

41、線，終止于0dB線的超前環節所補償。于是，該超前環節的零點轉角頻率應為于是，參考程序清單見fanli030。參考程序fanli030：滯后-超前校正程序ng=1;dg=0.5 1.5 1 0; % 校正前系統kc=10;dpm=45; %增益及相位裕度k_k=0.1; %滯后幅值穿越頻率調整系數w=logspace(-2,2,200); %計算頻率范圍bode(tf(kc*ng,dg),w),grid %校正前系統Bode圖mu,pu=bode(tf(kc*ng,dg),w); %校正前系統Bode響應mu_db=20*log10(mu); %幅值的分貝表示wg1=spline(mu_db,w

42、,0) %校正前系統的幅值穿越頻率wp2=spline(pu,w,-180) %校正前系統的相位穿越頻率wp3=k_k*wp2 %滯后環節零點轉角頻率t2_z=1/wp3 %滯后環節零點時間常數t2_p=10*t2_z %滯后環節極點時間常數nk2=t2_z 1 %滯后環節分子dk2=t2_p 1 %滯后環節分母ng1,dg1=series(kc*ng,dg,nk2,dk2) %含滯后環節的開環傳遞函數a1_db=spline(w,mu_db,wp2) %校前系統相位穿越頻率處的幅值wt=0.1*wp2/10(-a1_db/20) %超前環節的零點轉角頻率t1_z=1/wt %超前環節的零點時

43、間常數t1_p=0.1*t1_z %超前環節的極點時間常數nk1=t1_z 1 %超前環節的分子dk1=t1_p 1 %超前環節的分母nkc,dkc=series(nk2,dk2,nk1,dk1); %滯后超前校正環節sys1=series(tf(kc*ng,dg),tf(nkc,dkc) %滯后超前校正后的開環系統sysc=feedback(sys1,1) %滯后超前校正后的閉壞系統t=0:0.02:4;figure,step(tf(10,0.5 1.5 1 10),sysc,t),grid %校正前、后系統的單位階躍響應figure,bode(tf(kc*ng,dg),tf(ng1,dg1

44、),sys1,tf(nkc,dkc),0.001,1000),grid %校正前、后系統、滯后校正后系統及校正環節的Bode圖程序中k_k為滯后環節零點轉角頻率的調整系數，根據被校正系統不同，調整范圍超過一個量級。該例得出的滯后超前校正環節為對于穩定系統也可以采用滯后超前校正提高系統性能。fanli030滯后超前校正前后系統的單位階躍響應fanli030：滯后-超前校正前后系統及校正環節的Bode圖方法2：介紹兩次使用滯后超前兩用公式進展滯后超前校正的程序。滯后校正：設待校正系統G(s)的幅值穿越頻率為c，第一次進展滯后校正，選擇滯后校正后的系統幅值穿越頻率為c1=k1c，k11。使用公式求出

45、滯后校正環節G1(s)。超前校正：以滯后校正的系統G1(s)G(s)作為待校正系統，選擇超前校正后系統幅值穿越頻率c2=k2c，k21。使用同一組公式進展超前校正，得出超前校正環節G2(s)。校正后系統為G2(s) G1(s) G(s)。仿真說明，k1，k2的調節范圍較寬，超過一個數量級。例4.7 對系統進展滯后超前校正。ess10dB。參考程序見fanli031。參考程序fanli031:使用“兩用公式的滯后-超前校正ng=400;dg=1 200 200 0; % 校正前系統kc=5;dpm=45; %增益及相位裕度k1_k=0.01;k2_k=10; %滯后及超前幅值穿越頻率調整系數w=

46、logspace(-2,2,200);%計算頻率范圍bode(tf(kc*ng,dg),w),grid %校正前系統Bode圖mu,pu=bode(tf(kc*ng,dg),w); %校正前系統Bode響應mu_db=20*log10(mu); %幅值的分貝表示wc=spline(mu_db,w,0) %校正前系統的幅值穿越頻率wp1=spline(pu,w,-180) %校正前系統的相位穿越頻率pm1=spline(w,pu,wc*k1_k) %滯后環節零點轉角頻率處的相角wc1=wc*k1_k %滯后環節零點轉角頻率ngv1=polyval(ng,j*wc1);dgv1=polyval(d

47、g,j*wc1); %校正前系統在滯后環節零點轉角頻率處的函數值g1=ngv1/dgv1 %校正前系統在滯后環節零點轉角頻率處的頻率特性thetag1=angle(g1) %校正前系統在滯后環節零點轉角頻率處的相位弧度thetag_d1=thetag1*180/pi % 同前條相位度mg1=abs(g1) %校正前系統在滯后環節零點轉角頻率處的幅值mg_db1=20*log10(mg1) %同前條幅值分貝dpm_rad1=dpm*pi/180; %相位增益裕度弧度t_z1=(1+kc*mg1*cos(dpm_rad1-thetag1)/.(-wc1*kc*mg1*sin(dpm_rad1-th

48、etag1);t_p1=(cos(dpm_rad1-thetag1)+kc*mg1)/. (wc1*sin(dpm_rad1-thetag1);%滯后環節零點、極點時間常數計算式nk1=t_z1,1 %滯后環節分子dk1=t_p1,1 %滯后環節分母sysk1=series(tf(kc*ng,dg),tf(nk1,dk1) %含滯后環節的開環傳遞函數ng1,dg1=series(kc*ng,dg,nk1,dk1) %含滯后環節的開環傳遞函數mu2,pu2=bode(tf(ng1,dg1),w);%滯后校正之后系統的Bode響應mu_db2=20*log10(mu2); %滯后校正之后系統的幅值

49、分貝wcc=spline(mu_db2,w,0)% 滯后校正之后系統的幅值穿越頻率wpc=spline(pu2,w,-180)% % 滯后校正之后系統的相位穿越頻率pm2=spline(w,pu2,wc*k2_k) % 滯后校正之后的系統在期望幅值穿越頻率處的相位wc2=wc*k2_k %校正后系統的期望幅值穿越頻率ngv2=polyval(ng1,j*wc2);dgv2=polyval(dg1,j*wc2); %超前校正之前，系統在期望幅值穿越頻率處的函數值g2=ngv2/dgv2 %超前校正之前，系統在期望幅值穿越頻率處的頻率特性thetag2=angle(g2) %該點的相角弧度thet

50、ag_d2=thetag2*180/pi %該點的相角度mg2=abs(g2) %該點的幅值mg_db2=20*log10(mg2) %該點的幅值分貝dpm_rad2=dpm*pi/180;相位裕度t_z2=(1+mg2*cos(dpm_rad2-thetag2)/.(-wc2*mg2*sin(dpm_rad2-thetag2);t_p2=(cos(dpm_rad2-thetag2)+mg2)/. (wc2*sin(dpm_rad2-thetag2);%超前環節零點、極點時間常數計算式nk2=t_z2,1 %超前環節分子dk2=t_p2,1 %超前環節分母nk,dk=series(nk1,dk

51、1,nk2,dk2) %滯后超前環節的傳遞函數sysk=series(sysk1,tf(nk2,dk2)% 滯后超前校正后系統的開環傳遞函數sysc=feedback(sysk,1) %校正后系統的閉環傳遞函數t=0:0.02:10;figure,step(sysc,tf(kc*400,1 200 200 kc*400),t),grid %校正前后系統的單位階躍響應figure,bode(tf(kc*ng,dg),tf(nk,dk),sysk),grid%校正前、后系統及校正環節的Bode圖fanli031滯后超前校正前后系統的單位階躍響應(兩用公式)fanli031：滯后-超前校正前后系統及

52、校正環節的Bode圖(兩用公式)由仿真結果看，校正后快速性與超調均有明顯改善，唯嫌缺乏的是穩態誤差偏大，這可通過適當提高增益獲得解決。校正環節為由校正環節Bode可知，超前校正局部的相位補償已接近90。給出的校正參考程序，只要鍵入不同待校正系統及性能指標，運行即可獲得校正環節。不過，這些校正環節還僅是理論上的結果，假設需實際應用，尚須對其中的參數進展微調。1基于頻率法的串聯超前校正超前校正裝置的主要作用是通過其相位超前效應來改變頻率響應曲線的形狀，產生足夠大的相位超前角，以補償原來系統中元件造成的過大的相位滯后。因此校正時應使校正裝置的最大超前相位角出現在校正后系統的開環截止頻率處。例5-1：

53、單位反響系統的開環傳遞函數為，試確定串聯校正裝置的特性，使系統滿足在斜坡函數作用下系統的穩態誤差小于0.1，相角裕度。解：根據系統靜態精度的要求，選擇開環增益取，求原系統的相角裕度。num0=12; den0=2,1,0; w=0.1:1000;gm1,pm1,wcg1,wcp1=margin(num0,den0);mag1,phase1=bode(num0,den0,w);gm1,pm1,wcg1,wcp1margin(num0,den0) %計算系統的相角裕度和幅值裕度，并繪制出Bode圖grid;ans =Inf 11.6548 Inf 2.4240圖5-1 原系統的Bode圖由結果可知

54、，原系統相角裕度，不滿足指標要求，系統的Bode圖如圖5-1所示。考慮采用串聯超前校正裝置，以增加系統的相角裕度。確定串聯裝置所需要增加的超前相位角及求得的校正裝置參數。e=5; r=45; r0=pm1;phic=(r-r0+e)*pi/180;alpha=(1+sin(phic)/(1-sin(phic);將校正裝置的最大超前角處的頻率作為校正后系統的剪切頻率。那么有：即原系統幅頻特性幅值等于時的頻率，選為。根據=，求出校正裝置的參數。即。 il,ii=min(abs(mag1-1/sqrt(alpha);wc=w( ii);T=1/(wc*sqrt(alpha); numc=alpha*

55、T,1;denc=T,1; num,den=series(num0,den0,numc,denc); %原系統與校正裝置串聯gm,pm,wcg,wcp=margin(num,den);%返回系統新的相角裕度和幅值裕度 printsys(numc,denc) %顯示校正裝置的傳遞函數 disp(校正之后的系統開環傳遞函數為:); printsys(num,den) %顯示系統新的傳遞函數 mag2,phase2=bode(numc,denc,w);%計算指定頻率內校正裝置的相角范圍和幅值范圍 mag,phase=bode(num,den,w); %計算指定頻率內系統新的相角范圍和幅值范圍 sub

56、plot(2,1,1);semilogx(w,20*log10(mag),w,20*log10(mag1),-,w,20*log10(mag2),-.);grid; ylabel(幅值(db); title(-Go,-Gc,GoGc);subplot(2,1,2); semilogx(w,phase,w,phase1,-,w,phase2,-,w,(w-180-w),:);grid; ylabel(相位(0); xlabel(頻率(rad/sec);title(校正前：幅值裕量=,num2str(20*log10(gm1),db,相位裕量=,num2str(pm1),0;校正后：幅值裕量=,n

57、um2str(20*log10(gm),db,相位裕量=,num2str(pm),0);圖5-2 系統校正前后的傳遞函數及Bode圖2基于頻率法的串聯滯后校正滯后校正裝置將給系統帶來滯后相角。引入滯后裝置的真正目的不是為了提供一個滯后相角，而是要使系統增益適當衰減，以便提高系統的穩態精度。滯后校正的設計主要是利用它的高頻衰減作用，降低系統的截止頻率，以便能使得系統獲得充分的相位裕量。例5-2：單位反響系統的開環傳遞函數為，試確定串聯校正裝置的特性，使校正后系統的靜態速度誤差系數等于30/s，相角裕度，幅值裕量不小于10dB，截止頻率不小于2.3rad/s。解：根據系統靜態精度的要求，選擇開環增

58、益利用MATLAB繪制原系統的bode圖和相應的穩定裕度。num0=30; den0=conv(1,0,conv(0.1,1,0.2,1); w=logspace(-1,1.2);gm1,pm1,wcg1,wcp1=margin(num0,den0);mag1,phase1=bode(num0,den0,w);gm1,pm1,wcg1,wcp1margin(num0,den0) 圖5-3 原系統的Bode圖grid;ans = 0.5000 -17.2390 7.0711 9.7714由結果可知，原系統不穩定，且截止頻率遠大于要求值。系統的Bode圖如圖5-3所示，考慮采用串聯超前校正無法滿足要求，應選用滯后校正裝置。根據對相位裕量的要求，選擇相角為處的頻率作為校正后系統的截止頻率。確定原系統在新處的幅值衰減到0dB時所需的衰減量為。一般取校正裝置的轉折頻率分別為和。e=10; r=40; r0=pm1;phi=(-180+r+e); il,ii=min(abs(phase1-phi);wc=w( ii);