控制工程基礎(chǔ)(經(jīng)典控制部分)MATLAB分析.doc

控制工程基礎(chǔ)（經(jīng)典控制部分）的MATLAB分析第一章MATLAB的基本使用1-1MATLAB語言簡介MATLAB是一種高級矩陣語言，它由Math Works公司于1984年正式推出，它的基本處理對象是矩陣，即使是一個標量純數(shù)，MATLAB也認為它是只有一個元素的矩陣。隨著MATLAB的發(fā)展，特別是它所包含的大量工具箱（應(yīng)用程序集）的集結(jié)，使MATLAB已經(jīng)成為帶有獨特數(shù)據(jù)結(jié)構(gòu)、輸入輸出、流程控制語句和函數(shù)、并且面向?qū)ο蟮母呒壵Z言。MATLAB語言被稱為一種“演算紙式的科學(xué)計算語言”，它在數(shù)值計算、符號運算、數(shù)據(jù)處理、自動控制、信號處理、神經(jīng)網(wǎng)絡(luò)、優(yōu)化計算、模糊邏輯、系統(tǒng)辨識、小波分析、圖象處理、統(tǒng)計分析、甚至于金融財會等廣大領(lǐng)域有著十分廣泛的用途。MATLAB語言在工程計算與分析方面具有無可比擬的優(yōu)異性能。它集計算、數(shù)據(jù)可視化和程序設(shè)計于一體，并能將問題和解決方案以使用者所熟悉的數(shù)學(xué)符號或圖形表示出來。MATLAB語言和C語言的關(guān)系與C語言和匯編語言的關(guān)系類似。例如當我們需要求一個矩陣的特征值時,在MATLAB下只需由幾個字符組成的一條指令即可得出結(jié)果，而不必去考慮用什么算法以及如何實現(xiàn)這些算法等低級問題，也不必深入了解相應(yīng)算法的具體內(nèi)容。就象在C語言下不必象匯編語言中去探究乘法是怎樣實現(xiàn)的，而只需要采用乘積的結(jié)果就可以了。MATLAB語言還有一個巨大的優(yōu)點是其高度的可靠性。例如對于一個病態(tài)矩陣的處理，MATLAB不會得出錯誤的結(jié)果，而用C或其它高級語言編寫出來的程序可能會得出錯誤的結(jié)果。這是因為MATLAB函數(shù)集及其工具箱都是由一些在該領(lǐng)域卓有研究成果，造詣很深的權(quán)威學(xué)者經(jīng)過反復(fù)比較所得出來的最優(yōu)方法，而且經(jīng)過多年的實踐檢驗被證明是正確可靠的。1-2MATLAB的工作窗口下面以MATLAB6.1為例介紹。從實用的角度MATLAB的工作窗口包括命令窗口、M文件編輯器窗口、圖形編輯窗口、數(shù)學(xué)函數(shù)庫、應(yīng)用程序接口及在線窗口。下面首先介紹MATLAB的命令窗口及M文件編輯器。一、命令窗口啟動MATLAB之后，屏幕上自動出現(xiàn)命令窗口MATLAB，它是MATLAB提供給用戶的操作界面，用戶可以在命令窗口內(nèi)提示符“”之后（有的MATLAB版本命令窗口沒有提示符）鍵入MATLAB命令，回車即獲得該命令的答案。命令窗口內(nèi)有File、Edit、View、Web、Window、Help等菜單條。二、M文件編輯窗口M文件是MATLAB語言所特有的文件。用戶可以在M文件編輯窗口內(nèi)，編寫一段程序，調(diào)試，運行并存盤，所保存的用戶程序即是用戶自己的M文件。MATLAB工具箱中大量的應(yīng)用程序也是以M文件的形式出現(xiàn)的，這些M文件可以打開來閱讀，甚至修改，但應(yīng)注意，不可改動工具箱中的M文件！1進入M文件窗口有兩種方法1) 命令窗口 File New M-File；2) 命令窗口 點擊“File”字樣下面的圖標。M文件編輯窗口的標記是“Untitled”（無標題的）。當用戶編寫的程序要存盤時，Untitled作為默認文件名提供給用戶，自然，用戶可以，也應(yīng)當自己命名。若用戶不自己命名，則MATLAB會對Untitled進行編號。2M文件的執(zhí)行：返回命令窗口，在當前目錄（Current Directory）內(nèi)選擇所要運行的M文件的目錄，在命令窗口提示符“”后，直接鍵入文件名（不加后綴）即可運行。注意：（1）機器默認路徑為一級子目錄MATLAB6 p1work；（2）MATLAB 6.1以前的版本，運行M文件的方法稍有不同，它必須在File菜單下，打開“Run Script”子菜單，鍵入需要運行的文件路徑及名稱再回車，在這種情況下，work作為根目錄對待，不出現(xiàn)在M文件的路徑之中。本講義的參考程序都是在M文件窗口下編制的。三、在線幫助窗口在命令窗口中鍵入Help(空格) 函數(shù)名，可以立即獲得該函數(shù)的使用方法。1-3 MATLAB最基本的矩陣操作作為命令窗口及M文件編輯器的應(yīng)用實例，介紹幾個最基本的矩陣運算命令。一、矩陣的輸入在方括號內(nèi)依次按行鍵入矩陣元素，在一行內(nèi)的各元素之間用空格或逗號分開，每行之間用分號分開。例如，在命令窗內(nèi)輸入A=2 2 3;4 5 4;7 8 9 (注意：方括號，分號為矩陣行標記)B=1,3,5;6,-4,2;3,5,1 (逗號與空格功能相同)A223B135 4546-42 789 351同理：輸入A1246得到行矢量，輸入A22；4；6得到列矢量，于是，當輸入C=A；A1有C=123456789246A1作為矩陣C的最后一行，C和A相比，增加了一行。二、矩陣的轉(zhuǎn)置矩陣A的轉(zhuǎn)置用A表示，顯然，A1與A2互為轉(zhuǎn)置，即A1會得到以2,4，6為元素的列矢量。思考一下輸入C1 AA2 C2 AA1有什么結(jié)果？而輸入A；A1有無意義？三、矩陣的四則運算1矩陣的加減法：當兩個矩陣維數(shù)相同時可以直接進行“”或“”運算。如 D1AB，D2AB2 矩陣的乘法：當矩陣A，B維數(shù)相容時C3AB：普通意義下的矩陣相乘C4A .B：矩陣A與B的對應(yīng)元素相乘顯然，ABBA（一般情況），而A .BB .A。A .B稱為數(shù)列型乘法，它要求參加運算的矩陣或數(shù)列具有相同的行列數(shù)，這是MATLAB語言中的一種特殊運算，它在今后求取函數(shù)值等運算時是很重要的。實際上，前面所述的矩陣加、減法就是一種數(shù)列型運算。3 矩陣的除法D4AB：表示A-1B或inv(A)*B，即A的逆矩陣左乘矩陣B。D5B/A：表示BA-1或Binv(A)，即A的逆矩陣右乘B。D6A .B：表示B的每一個元素被A的對應(yīng)元素除。D7A ./B：表示A的每一個元素被B的對應(yīng)元素除。顯然，A .B與A ./B的各對應(yīng)元素互為倒數(shù)。讀者可以思考一下，D6 .D7等于什么？D8inv(A)：A的逆矩陣。打開M文件編輯窗口，將上述命令依次鍵入，得到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%矩陣的對應(yīng)元素相乘D3=AB%A的逆左乘BD4=B/A%A的逆右乘BD6=A.B%B的各元素被A的對應(yīng)元素除D7=A./B%A的各元素被B的對應(yīng)元素除D8=inv(A) %A的逆矩陣語句后面的%為語句說明符。MATLAB中矩陣運算的其它主要命令可通過在線幫助獲得。1-4MATLAB的符號運算操作一、進入符號運算功能在命令窗口鍵入syms x y z t 此后，即可以使用x，y，z，t等作自變量定義函數(shù)。syms x y z t real 規(guī)定所定義的變量為實型。二、代數(shù)方程求解使用命令solve 可以求解代數(shù)方程，如求下例方程組的解，命令為x,y=solve(x2+x*y+y-3=0,x2-4*x-2*y+3=0)程序見范例程序fanli002。參考程序fanli002:代數(shù)方程求解syms x y%進入符號運算功能f1=x2+x*y+y-3%函數(shù)f1f2=x2-4*x-2*y+3 %函數(shù)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) %將符號變量轉(zhuǎn)換成浮點數(shù)y=double(y) %將符號變量轉(zhuǎn)換成浮點數(shù)f1=subs(f1) %用浮點數(shù)x,y檢驗方程1f2=subs(f2) %用浮點數(shù)x,y檢驗方程2由solve求出的根是根的符號表達形式，是準確解。命令simplify(f)表示化簡，subs(f1)表示將求解出x，y代回f1中；double是將符號變量轉(zhuǎn)換成浮點數(shù)，是準確解x，y的近似值。這從程序運行f1a=0，f2a=0，f10，f20可以確認。三、符號矩陣運算符號矩陣可以和數(shù)值矩陣一樣進行運算，例如：求矩陣特征值eig，求矩陣的逆inv等命令都支持符號運算。設(shè)計算其特征值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) %化簡A的逆矩陣表達式C1=A*B2 %檢驗A的逆矩陣C1=simple(C1) %C1為單位矩陣注意函數(shù)的輸入方法，自變量用圓括號括起來。四、微積分運算設(shè)函數(shù)，則MATLAB中微積分運算命令為fiff(f)：求函數(shù)f對自變量x的一階導(dǎo)數(shù)；diff(f，2)：求函數(shù)f對自變量x的二階導(dǎo)數(shù)；int(f)：求函數(shù)f的不定積分例1.1，設(shè)試計算其一階，二階導(dǎo)數(shù)，積分運算，并作出函數(shù)圖象，見范例fanli004。參考程序fanli004:函數(shù)的微分與積分syms xf=1/(5+4*cos(x)ezplot(f) %函數(shù)f的曲線f1=diff(f) %函數(shù)f的一階導(dǎo)數(shù)figure,ezplot(f1) %函數(shù)f一階導(dǎo)數(shù)的曲線f2=diff(f,2) %函數(shù)f的二階導(dǎo)數(shù)figure,ezplot(f2) %函數(shù)f二階導(dǎo)數(shù)的曲線g=int(int(f2) %函數(shù)f的二階導(dǎo)數(shù)f2的二重積分figure,ezplot(g) %函數(shù)f2二重積分的曲線e=f-g %二階導(dǎo)數(shù)的二重積分與原函數(shù)的差e=simple(e)figure,ezplot(e)程序中ezplot(f)：作函數(shù)的圖形，x的取值范圍默認值為-2 x 2fanli004的函數(shù)曲線fanli004的一階導(dǎo)函數(shù)曲線fanli004的二階導(dǎo)函數(shù)曲線fanli004的二階導(dǎo)函數(shù)的二重積分曲線ezplot(f)是一個很有用的作圖命令，它的其它應(yīng)用形式，請查在線幫助。細心的讀者會發(fā)現(xiàn)，一個函數(shù)求二階導(dǎo)數(shù)后再對二階導(dǎo)數(shù)進行二重積分，其結(jié)果與原函數(shù)相差一個常數(shù)。相當于縱坐標發(fā)生平移。第二章系統(tǒng)的時域特性2-1傳遞函數(shù)一、傳遞函數(shù)的兩種形式傳遞函數(shù)通常表達成s的有理分式形式及零極點增益形式。設(shè)傳遞函數(shù)1有理分式形式分別將分子、分母中s多項式的系數(shù)按降冪排列成行矢量，缺項的系數(shù)用0補齊。上述函數(shù)可表示為num1=2 1 %（注意：方括號，同一行的各元素間留空格或逗號）。den1=1 2 2 1syss1=tf(num1,den1)運行后，返回傳遞函數(shù)的形式。這種形式不能直接進行符號運算！2零極點增益形式Z,P,K = tf2zp(num1,den1)sys2 = zpk(Z,P,K)返回零、極點、增益表達式，其Z，P分別將零點和極點表示成列向量，若無零點或極點用 （空矩陣）代替。運行得到的點 Z = -0.5極點P= -1，-0.5j0.866增益K = 2指令zp2tf(Z，P，K)將零極點增益變換成有理分式形式，見程序fanli005。參考程序fanli005:傳遞函數(shù)的有理分式及零極點增益模型num1=2 1 % 傳遞函數(shù)的分子系數(shù)向量den1=1 2 2 1 % 傳遞函數(shù)的分母系數(shù)向量sys1=tf(num1,den1) % 傳遞函數(shù)的有理分式模型Z,P,K=tf2zp(num1,den1) % 有理分式模型轉(zhuǎn)換成零極點增益模型num2,den2=zp2tf(Z,P,K) % 零極點增益模型轉(zhuǎn)換成有理分式模型sys2=zpk(Z,P,K) % 傳遞函數(shù)的零極點增益模型A1,B1,C1,D1=tf2ss(num1,den1) % 有理分式模型轉(zhuǎn)換成狀態(tài)空間模型A2,B2,C2,D2=zp2ss(Z,P,K)% 零極點及增益模型轉(zhuǎn)換成狀態(tài)空間模型num1,den1=ss2tf(A1,B1,C1,D1)% 狀態(tài)空間模型轉(zhuǎn)換成有理分式模型Z,P,K=ss2zp(A2,B2,C2,D2) % 狀態(tài)空間模型轉(zhuǎn)換成零極點增益模型程序中，命令tf2ss，zp2ss及ss2tf，ss2zp是狀態(tài)空間模型與有理分式及零、極點、增益模型之間的相互轉(zhuǎn)換。二、傳遞函數(shù)框圖的處理用框圖可以方便地表示傳遞函數(shù)的并聯(lián)，串聯(lián)及反饋。為簡潔，僅以有理分式模型為例。G1G1G1+G21 并聯(lián)sysp = parallel(sys1,sys2)num,den=parallel(num1,den1,num2,den2)2 串聯(lián)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) %默認值（-1）numc, 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) %(單位反饋)上面給出了同一指令的兩種形式，相當于兩套平行指令。對于零極點增益形式，書寫稍復(fù)雜一些，可先用zpk轉(zhuǎn)換成系統(tǒng)形式，或用zp2tf轉(zhuǎn)折換成有理分式形式后再進行框圖化簡操作。三、簡單函數(shù)的拉普拉斯變換在MATLAB的符號功能中，可以對簡單函數(shù)進行拉普拉斯正、逆變換。拉氏正變換：laplace(f(t)拉氏逆變換：ilaplace(L(s)其中為原函數(shù)，為象函數(shù)。命令格式參見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)%原函數(shù)f1L1=laplace(f1) %f1的拉氏變換（象函數(shù)）L1=simple(L1) %化簡f2=ilaplace(L1)%L1的拉氏逆變換f2=simple(f2) %化簡在MATLAB中使用laplace及ilaplace命令時，要注意象、原函數(shù)的符號，特別是對初相不等于零的振蕩系統(tǒng)，運行結(jié)果常常同手冊上的結(jié)果相差一個符號，這要注意函數(shù)表達式成立的條件。保險的辦法是再使用拉氏變換的初值定理確定象、原函數(shù)的符號。2-2系統(tǒng)時域特性曲線在MATLAB中，當傳遞函數(shù)已知時，可以方便地求出系統(tǒng)的單位脈沖響應(yīng)、單位階躍響應(yīng)等曲線。一、系統(tǒng)的單位階躍響應(yīng)stepstep有以下幾種格式step(sys)：直接作出sys的單位階躍響應(yīng)曲線。其中sys = tf(num, den) 或sys = zpk(z, p, k)，MATLAB自動決定響應(yīng)時間。step(sys, t)設(shè)定響應(yīng)時間的單位階躍響應(yīng)。t可以設(shè)定為最大響應(yīng)時間 t = t終值(秒)，也可以設(shè)置為一個向量t = 0 : t : t終值注意冒號的使用。它產(chǎn)生一個從0到t終值的行矢量，元素之間的間隔為t。step(sys1, sys2, , sysn)在同一幅圖上畫出幾個系統(tǒng)的單位階躍響應(yīng)。y, t = step(sys)；命令輸出對應(yīng)時刻t的各個單位階躍響應(yīng)值，不畫圖。語句后的分號控制數(shù)據(jù)的屏幕顯示。如果要查看機器計算了多少個數(shù)據(jù)，可以使用命令size(y)得出的結(jié)果也表明數(shù)據(jù)作為列矢量的行數(shù)。要將計算出的y, t作成曲線，使用一般的作圖命令plot(t, y)plot后面跟的兩個參數(shù)橫坐標在前，縱坐標在后。參考程序見fanli008：參考程序fanli008:系統(tǒng)的單位階躍響應(yīng)num1=4 2den1=2 8 14 11 4sys1=tf(num1,den1) %系統(tǒng)G1(s)num2=2 1den2=1 4 6 7 3sys2=tf(num2,den2) %系統(tǒng)G2(s)y1,t1=step(sys1); %系統(tǒng)G1(s)的單位階躍響應(yīng)數(shù)據(jù)y2,t2=step(sys2); %系統(tǒng)G2(s)的單位階躍響應(yīng)數(shù)據(jù)step(sys1,sys2) %系統(tǒng)G1(s)、G2(s)的單位階躍響曲線figure,step(sys1,sys2,20) %系統(tǒng)G1(s)、G2(s)在自選時間（20秒）內(nèi)的單位階躍響曲線figure,plot(t1,y1) %系統(tǒng)G1(s)的單位階躍響應(yīng)曲線figure,plot(t2,y2) %系統(tǒng)G2(s)的單位階躍響應(yīng)曲線fanli008:step(sys1,sys2,t)單位階躍曲線二、系統(tǒng)的單位脈沖響應(yīng)impulse命令格式與單位階躍響應(yīng)step的命令格式完全相同，只需將語句中的step用impulse代替即可。針對同樣的系統(tǒng)，其單位脈沖響應(yīng)的參考程序見fanli009。參考程序fanli009:系統(tǒng)的單位脈沖響應(yīng)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); %系統(tǒng)G1(s)的單位脈沖響應(yīng)數(shù)據(jù)y2,t2=impulse(sys2); %系統(tǒng)G2(s)的單位脈沖響應(yīng)數(shù)據(jù)impulse(sys1,sys2) %系統(tǒng)G1(s)、G2(s)的單位脈沖響應(yīng)曲線figure,impulse(sys1,sys2,20)%系統(tǒng)G1(s)、G2(s)在自選時間（20秒）內(nèi)的單位脈沖響應(yīng)曲線figure,plot(t1,y1) %系統(tǒng)G1(s)的單位脈沖響應(yīng)曲線figure,plot(t2,y2) %系統(tǒng)G2(s)的單位脈沖響應(yīng)曲線hold on,step(sys2) %系統(tǒng)G2(s)的單位階躍和單位脈沖響應(yīng)曲線fanli009:impulse(sys1,sys2)單位脈沖響應(yīng)曲線程序的最后一句hold on是當前圖形保護模式。當要將新圖形作在當前圖形上時，必須使用hold on。而figure的含意是另開一個新的圖形窗口，如果不用figure或hold on，則新的圖形會占用原圖形窗口，始終只保留一個最新的圖形窗口。三、一階系統(tǒng)及二階系統(tǒng)的時域特性一階系統(tǒng)及二階系統(tǒng)是最基本也是最重要的系統(tǒng)，高階系統(tǒng)總可以視為由若干個一階和（或）二階系統(tǒng)組合構(gòu)成。1 一階系統(tǒng)（設(shè)增益為1）影響系統(tǒng)特性的參數(shù)是其時間常數(shù)T，T越大，系統(tǒng)慣性越大，響應(yīng)越慢。參考程序fanli010給出了T0.4, 1.2, 2.0, 2.8, 3.6, 4,4六條單位階躍響應(yīng)曲線。參考程序fanli010:一階系統(tǒng)的單位階躍響應(yīng)曲線num=1;i=1;for del=0.1:0.2:1.1 %一階系統(tǒng)時間常數(shù)遞增間隔den=4*del 1; %一階系統(tǒng)分母向量step(tf(num,den) %一階系統(tǒng)單位階躍響應(yīng)曲線hold on, %不同時間常數(shù)的一階系統(tǒng)單位階躍響應(yīng)曲線簇i=i+1;end同理，可以作出對應(yīng)的單位脈沖響應(yīng)曲線，參考程序fanli011。參考程序fanli011:一階系統(tǒng)的單位脈沖響應(yīng)曲線num=1;i=1;fanli010一階系統(tǒng)時間常數(shù)對單位階躍響應(yīng)的影響fanli011一階系統(tǒng)時間常數(shù)對單位脈沖響應(yīng)的影響for del=0.1:0.2:1.1den=4*del 1;impulse(tf(num,den),10) %一階系統(tǒng)單位階脈沖應(yīng)曲線hold on, %不同時間常數(shù)的一階系統(tǒng)單位脈沖響應(yīng)曲線簇 i=i+1;end注意MATLAB中for語句的結(jié)構(gòu)。讀者可以改變不同的增益，看看圖形有何變化。2 二階系統(tǒng)（設(shè)01）設(shè)二階系統(tǒng)為二階系統(tǒng)的特征參數(shù)為固有頻率及阻尼比。當增大，系統(tǒng)振動頻率加快，振蕩加劇；而隨著減小，系統(tǒng)振蕩加劇，振蕩峰尖銳。參考程序fanli012示出了當0.5, =1, 2, 3, 4, 5 rad/s時的間接階躍曲線簇。參考程序fanli012: 不同固有頻率的二階系統(tǒng)的單位階躍響應(yīng)曲線（=0.5）i=1;for del=1:1:5; % 二階系統(tǒng)固有頻率遞增間隔num=del2; % 二階系統(tǒng)傳遞函數(shù)分子系數(shù)向量den=1 del del2; % 不同固有頻率的二階系統(tǒng)分母系數(shù)向量step(tf(num,den),6) %二階系統(tǒng)單位階躍響應(yīng)曲線hold on, %不同固有頻率的二階系統(tǒng)單位階躍響應(yīng)曲線簇i=i+1;endfanli012二階系統(tǒng)固有頻率對單位階躍響應(yīng)的影響參考程序fanli013示出了同一二階系統(tǒng)的單位脈沖響應(yīng)曲線簇。參考程序fanli013: 不同固有頻率的二階系統(tǒng)的單位脈沖響應(yīng)曲線（=0.5）i=1;for del=1:1:5;num=del2;den=1 del del2; % 不同固有頻率的二階系統(tǒng)分母系數(shù)向量impulse(tf(num,den),6) %二階系統(tǒng)單位脈沖響應(yīng)曲線hold on, %不同固有頻率的二階系統(tǒng)單位脈沖響應(yīng)曲線簇i=i+1;endfanli013二階系統(tǒng)固有頻率對單位脈沖響應(yīng)的影響參考程序fanlio14示出了當1，0.1，0.3，0.5，0.7，0.9的二階系統(tǒng)的單位階躍響應(yīng)曲線簇。參考程序fanli014: 不同阻尼比的二階系統(tǒng)的單位階躍響應(yīng)曲線（=1）i=1;for del=0.1:0.2:0.9; % 二階系統(tǒng)阻尼比遞增間隔 num=1; den=1 2*del 1; % 不同阻尼比的二階系統(tǒng)分母系數(shù)向量 step(tf(num,den),30) hold on, %不同阻尼比的二階系統(tǒng)單位階躍響應(yīng)曲線簇 i=i+1;endfamli014二階系統(tǒng)阻尼比對單位階躍響應(yīng)的影響參考程序fanli015示出了同一二階系統(tǒng)當取0.1，0.3，0.5，0.7，0.9的單位脈沖響應(yīng)曲線簇。參考程序fanli015: 不同阻尼比的二階系統(tǒng)的單位脈沖響應(yīng)曲線（=1）i=1;for del=0.1:0.2:0.9; num=1; den=1 2*del 1; impulse(tf(num,den),30) hold on, %不同阻尼比的二階系統(tǒng)單位脈沖響應(yīng)曲線簇 i=i+1;endfanli015二階系統(tǒng)阻尼比對單位脈沖響應(yīng)的影響四、帶延時環(huán)節(jié)系統(tǒng)的典型響應(yīng)設(shè)具有純延時環(huán)節(jié)的傳遞函數(shù)為計算這種系統(tǒng)的單位階躍響應(yīng)不能使用一般方法。首先應(yīng)使用Pade法對延時環(huán)節(jié)進行近似展開，numT, dent = pade(T, 5) （表示使用5階Pade級數(shù)）得出的分子、分母系數(shù)的向量，然后將系統(tǒng)視為慣性環(huán)節(jié)與延遲環(huán)節(jié)的串聯(lián)，再求其階躍脈沖響應(yīng)。設(shè)K = 2.5，T1，程序見fanli016。參考程序fanli016: 帶延時環(huán)節(jié)系統(tǒng)的單位階躍響應(yīng)K=2.5;T=1; T1=0.5;numt,dent=pade(T,5);%5階Pade法近似延時環(huán)節(jié)syst=tf(numt,dent); %延時環(huán)節(jié)的近似傳遞函數(shù)num1=K;den1=T1 1;sys1=tf(num1,den1);%慣性環(huán)節(jié)sys=series(sys1,syst);% 帶延時環(huán)節(jié)的慣性系統(tǒng)step(sys,6),grid % 帶延時環(huán)節(jié)慣性系統(tǒng)的單位階躍響應(yīng)figure,impulse(sys,6),grid% 帶延時環(huán)節(jié)慣性系統(tǒng)的單位脈沖響應(yīng)fanli016帶延遲的一階系統(tǒng)的單位階躍響應(yīng)fanli016帶延遲的一階系統(tǒng)的單位脈沖響應(yīng)從響應(yīng)曲線上可見，具有延時環(huán)節(jié)的慣性系統(tǒng)當加上單位階躍或單位脈沖輸入前后，其響應(yīng)初始有振蕩特性，脈沖響應(yīng)更為明顯。這種振蕩現(xiàn)象使得對這種系統(tǒng)進行校正需要使用特殊的方法。在作圖命令后，加上格線命令grid，機器自動給曲線圖加上格線。2-3響應(yīng)曲線的動態(tài)分析一、MATLAB圖形的編輯在MATLAB中，用于編輯圖形的命令很多，因為可視化正是MATLAB的一種強大而優(yōu)越的性能。從工程實用角度，僅在圖形窗口進行編輯已經(jīng)夠用了。在圖形窗口頂部的第二排依次排有下列命令按鈕：Edit Plot（）：圖形編輯。點擊進入圖形編輯功能。點擊選擇編輯對象。例如，曲線，標題，縱、橫坐標的說明等，可對選中項進行編輯。Insert Text（）：插入文本。點擊可在圖區(qū)內(nèi)插入文字，并可選擇字體及大小等。Insert Arrow（）：插入箭頭。點擊可在圖區(qū)內(nèi)畫箭頭，并可對該箭頭進行編輯。Insert Line()：插入直線。點擊可在圖區(qū)內(nèi)插入直線，并可對該直線進行編輯，改變線寬及顏色。Zoom In（）：曲線放大。激活該按鈕，可以逐次放大曲線，以觀察曲線某些部分的細節(jié)，這對觀察變化劇烈的曲線部分很有幫助。Zoom Out（）：曲線縮小。激活該按鈕，可以逐次從放大狀態(tài)返回到原狀態(tài)。Rotate 3D（）：三維旋轉(zhuǎn)。激活該按鈕，在圖區(qū)內(nèi)按住鼠標左鍵，拖曳鼠標，原二維曲線可三維空間內(nèi)任意旋轉(zhuǎn)。此外，還有平面旋轉(zhuǎn)、三維動畫、前后放大、平移、飛逸等。例如，將曲線加粗，以利作圖（曲線默認寬度為0.5）。按Edit Plot 選擇曲線雙擊左鍵 出現(xiàn)線寬選擇框（line width） 選擇線寬 OK。二、對曲線進行數(shù)據(jù)分析不激活Edit Plot按鈕（圖形編輯不使能）， 讀取曲線上任一點的坐標值：鼠標指向選擇點，單擊左鍵，顯示曲線編號，橫坐標及縱坐標值； 研究曲線擬合情況：按Tools Basic Fitting，出現(xiàn)Plot fitt（圖形擬合）窗口，選擇需要擬合的曲線，點擊擬合方式，即可對所選圖形進行曲線擬合。MATLAB提供了樣條插值、保形插值（shape-preserving interpolant）及直到10階的多項式擬合。當選擇多項式擬合時，機器會給出曲線及擬合方程。在所給出的擬合方式中，三次樣條插值（Cubic spline interpolant）及保形插值（shape-preserving interpolant）效果最好。第三章 系統(tǒng)的頻率特性系統(tǒng)的頻率特性就是，這是一個復(fù)變函數(shù)。3-1MATLAB中的復(fù)數(shù)及復(fù)變量在MATLAB中復(fù)數(shù)x的輸入方法：x = 2 + 3*i或2 + 3*j也可以s = 2 + 3i或s = 3 + 2j在符號運算功能下，復(fù)變函數(shù)表示為（設(shè)a，b都是實數(shù)）：syms a b realx = ab*I或ab*j（*不能省略）在MATLAB中與復(fù)數(shù)有關(guān)的命令有abs(x)：求實數(shù)的絕對值及復(fù)數(shù)的模，支持符號運算；angle(x)：求復(fù)數(shù)的相角，單位為弧度，不支持符號運算。conj(x)：求復(fù)數(shù)的共軛復(fù)數(shù)，支持符號運算；imag(x)：求復(fù)數(shù)的虛部，支持符號運算；real(x)：求復(fù)數(shù)的實部，支持符號運算。復(fù)數(shù)的運算見參考程序fanli017。參考程序fanli017:復(fù)數(shù)運算（含符號功能）syms a b c d real%符號均為實變量X1=2-3i X2=a+c/b*jX3=c/a-2*d*jS1=X2*X3/X1 %復(fù)數(shù)的乘除S2=simple(abs(S1) %復(fù)數(shù)的模S3=simple(real(S1) %復(fù)數(shù)的實部S4=simple(imag(S1) %復(fù)數(shù)的虛部S5=simple(conj(S1) %復(fù)數(shù)的共軛S6=simple(expand(S1*S5) %復(fù)數(shù)與其共軛的積S7=simple(S32+S42) %復(fù)數(shù)實部、虛部的平方和S8=simple(S6-S7) %檢驗復(fù)數(shù)運算的正確性fanli018給出了系統(tǒng)的頻率特性。參考程序fanli018:系統(tǒng)頻率特性的解析表達syms w real %定義頻率w 為實型變量g=50*(0.6*w*i+1)/(i*w)2/(4*w*j+1) %系統(tǒng)頻率特性g1=simple(conj(g) %系統(tǒng)頻率特性的復(fù)共軛gr=simple(real(g) %系統(tǒng)實頻特性gi=simple(imag(g) %系統(tǒng)虛頻特性ga=simple(abs(g) %系統(tǒng)幅頻特性（模）gb=simple(sqrt(gr2+gi2) %系統(tǒng)幅頻特性的另一計算方法gc=simple(ga-gb) %檢驗系統(tǒng)幅頻特性的兩種計算方法3-2頻率特性的Nyquist圖將系統(tǒng)傳遞函數(shù)中的復(fù)變數(shù)s用純虛數(shù)j（為角頻率，rad/s）代替，即得到系統(tǒng)的頻率特性，又稱為諧波傳遞函數(shù)，有三種表示方法：式中：實頻特性：虛頻特性：幅頻特性：相頻特性在研究控制系統(tǒng)的特性時，采用圖示方法比采用解析方法直觀、簡單，便于對動態(tài)特性進行更深入的研究。以下討論設(shè)系統(tǒng)傳遞函數(shù)為sys。一、Nyquist圖1 Nyquist 圖的有關(guān)命令nyquist(sys)直接返回系統(tǒng)sys的Nyquist圖，機器自己確定頻率范圍（通常為- +）nyquist(sys, w)返回在規(guī)定頻率（rad/s）范圍內(nèi)的Nyquist圖，頻率范圍的格式為w = wmin, wmax (0wminwmin)nyquist(sys1,sys2,sysn)nyquist(sys1,sys2,sysn,w)在同一幅度圖上畫出幾個系統(tǒng)的Nyquist圖，可以由機器確定頻率范圍，也可以自選頻率范圍。re,im,w=nyquist(sys)分別返回系統(tǒng)sys對應(yīng)頻率w的實部（re）及虛部（im）值。讀取實部，虛部及對應(yīng)頻率點的方法：re(1,1,k)im(1,1,k)w(k) (rad/s)其中，k為機器計算Nyquist響應(yīng)的序列編號，必須是正整數(shù)而且應(yīng)當在計算范圍內(nèi)?？梢允褂妹顂ize(w)查看。k值隨問題的簡、繁而變化。這是一個非常有用的命令，當然也可以采用直接在Nyquist圖上用鼠標讀取的方法。re,im = nyquist(sys,w)返回在某一頻率w（rad/s）或某一頻率范圍w = wmin, wmax 系統(tǒng)的實部及虛部值。當選定為頻率范圍時，調(diào)出實部，虛部及對應(yīng)頻率值的方法同前一指令。例3.1 分析系統(tǒng)的頻率特性。該傳遞函數(shù)是一個零極點及有理分式的混合形式。首先將子系統(tǒng)變換成零極點增益形式z1，p1，k1，再將系統(tǒng)變換成零極點增益模型（z，p，k）再作Nyquist圖。使用re，im，wnyquist(sys)語句，可計算出若干頻率點所對應(yīng)的實部、虛部值，而且在w(25)=0.8802rad/sw(26)=0.8951rad/s之間，曲線將從第三象限穿過負實軸進入第二象限，在w(0.8802，0.8951)之間使用語句re，im=nyquist(sys，w)0.8802，w，0.8951進行更細致的計算，在需要的精度內(nèi)求出穿越負實軸的參數(shù)。參考程序見fanli019。參考程序fanli019:用Nyquist圖分析系統(tǒng)頻率特性num=5;den=1 0.2 1;%子系統(tǒng)的傳遞函數(shù)Z1,P1,K1=tf2zp(num,den) %子系統(tǒng)的極點Z=-2.5;-2 %系統(tǒng)的零點P=P1;1/1.2;-1/1.4 %系統(tǒng)的極點K=num %系統(tǒng)的增益sys=zpk(Z,P,K) %系統(tǒng)的零、極點、增益模型nyquist(sys) %系統(tǒng)的Nyquist圖re,im,w=nyquist(sys) %計算系統(tǒng)若干頻率點所對應(yīng)的實部與虛部值w(25) %w(25)=0.8802,曲線處于第三象限靠近負實軸的頻率點w(26) %w(26)=0.8951,曲線處于第二象限靠近負實軸的頻率點re1,im1=nyquist(sys,0.8802:0.00005:0.8951) %在（0.8802, 0.8951）(rad/s)范圍內(nèi)，以0.00005rad/s為間隔計算實部及虛部值fanli019的Nyquist圖2 帶延時環(huán)節(jié)系統(tǒng)的Nyquist圖對帶有延時環(huán)節(jié)的系統(tǒng)，由于不能表示成傳遞函數(shù)的有理分式，所以MATLAB中沒有直接命令可用。但由于Nyquist圖實際上是以頻率為參變量，橫坐標為實頻特性，縱坐標為虛頻特性的直角坐標系圖。所以，可以通過計算實頻特性與虛頻特性，再使用plot(x，y)作出Nyquist圖。例3.2 求作具有延時環(huán)節(jié)的系統(tǒng)的Nyquist圖。程序見fanli020。參考程序fanli020:帶延時環(huán)節(jié)系統(tǒng)的Nyquist圖w=0:0.01:30;%設(shè)定頻率變化范圍gjw=5*exp(-1.5*j*w)./(1+j*2*w);%計算頻率特性x=real(gjw); %實頻特性y=imag(gjw); %虛頻特性plot(x,y),grid %帶延時環(huán)節(jié)系統(tǒng)的Nyquist圖fanli020帶延時環(huán)節(jié)的慣性系統(tǒng)的Nyquist圖注意，程序中兩項相除時采用了對應(yīng)元素相除的方法，格式為：分子./分母。3-3系統(tǒng)的Bode圖Bode圖是系統(tǒng)的對數(shù)頻率特性曲線，由兩幅圖構(gòu)成：對數(shù)幅頻特性圖和對數(shù)相頻特性圖?？v坐標以線性分度，幅值，以分貝表示；以度表示。兩張圖共用橫軸。橫坐標為頻率，采用對數(shù)分度（），但習(xí)慣上仍標真數(shù)。因而，橫坐標是按十倍頻程均勻分度的。對于寫成典型環(huán)節(jié)形式的傳遞函數(shù)，其對應(yīng)頻率特性為幅頻特性相頻特性一、Bode圖的有關(guān)命令Bode圖的命令格式與Nyquist圖的命令格式十分相似，其作用對象仍然是由有理分式或零極點模型構(gòu)成的系統(tǒng)傳遞函數(shù)sys。Bode圖的有關(guān)命令如下：Bode(sys) 直接返回bode圖bode(sys, w) 在設(shè)定頻率范圍w=wmin, wmax或w=wmin: w:wmax內(nèi)的Bode圖。Bode(sys1，sys2，sysN)Bode(sys1，sys2，sysN，w)在同一幅圖上作幾個系統(tǒng)的Bode圖。mag，phase，wbode(sys)其中，w表示一系列計算點的頻率（rad/s）phase：對應(yīng)計算頻率點的相角（度）mag：對應(yīng)計算頻率點的幅值（非分貝）若要用分貝表示mag(db)=20*log10(mag)讀取mag，phase，w的方法與Nyquistg命令中讀取實部、虛部及頻率的方法相同。mag，phasebode(sys，w)返回某一頻率點w，或某一頻率范圍w =wmin, wmax或w = wmin: w:wmax的幅值及相角。當選它為頻率范圍時，讀取方法相同。例3.3 作系統(tǒng)的Bode圖。注意該系統(tǒng)是一個II型系統(tǒng)。程序見fanli021。參考程序fanli021:型系統(tǒng)的Bode圖num=1 0.1 7.5;den=1 0.12 9 0 0;sys=tf(num,den)bode(sys,1,5),gridmag,phase,w=bode(sys); %計算系統(tǒng)若干頻率點所對應(yīng)的幅值與相角mag,phase=bode(sys,2.8479:0.001:2.8719) %計算小頻率范圍的幅值與相角fanli021的Bode圖程序的最后一句，用以計算相頻特性曲線峰點所對應(yīng)的相角，頻率及幅值。例3.4對比Nyquist圖及Bode圖。fanli022系統(tǒng)的Nyquist圖fanli022系統(tǒng)的Bode圖參考程序fanli022:系統(tǒng)Nyquist 圖與Bode圖的對比num=50;den=0.01 0.14 1;Z1,P1,K1=tf2zp(num,den) %子系統(tǒng)的極點Z=-2.5;-1/0.2P=P1;1/0.9;-1/0.3K=numsys=zpk(Z,P,K) %系統(tǒng)的零極點增益模型bode(sys),grid %系統(tǒng)的Bode圖figure,nyquist(sys) %系統(tǒng)的Nyquist圖mag,phase,w=bode(sys);%系統(tǒng)幅值及相角的計算二、典型環(huán)節(jié)的Nyquist圖和Bode圖根據(jù)典型環(huán)節(jié)的傳遞函數(shù)，使用前述的Nyquist圖及Bode圖的相關(guān)命令，很容易作出它們的這兩種頻率特性圖。讀者應(yīng)重視典型環(huán)節(jié)Nyquist圖及Bode圖的作法及走向，以及它們所反映的系統(tǒng)動態(tài)特性，研究在原有系統(tǒng)增加或減少一個典型環(huán)節(jié)時對系統(tǒng)的特性曲線所引起的變化以及這種變化的物理及工程含義。例3.5：對比及的Nyquist圖及Bode圖，特別注意Bode圖中的相角變化，程序見fanli023。參考程序fanli023:系統(tǒng)Nyquist 圖與Bode圖的對比num1=1 -2;num2=-1 2;den=1 0.1 4;sys1=tf(num1,den) %系統(tǒng)1sys2=tf(num2,den) %系統(tǒng)2bode(sys1),grid %系統(tǒng)1的Bode圖figure,bode(sys2),grid %系統(tǒng)2的Bode圖figure,nyquist(sys1),grid %系統(tǒng)1的Nyquist圖figure,nyquist(sys2),grid %系統(tǒng)2的Nyquist圖fanli023:系統(tǒng)1的Bode圖fanli023:系統(tǒng)2的Bode圖fanli023:系統(tǒng)1的Nyquist圖fanli023:系統(tǒng)2的Nyquist圖Nyquist圖中的虛線是等幅值圓。第四章系統(tǒng)的校正4-1系統(tǒng)的穩(wěn)定性分析系統(tǒng)穩(wěn)定的基本條件是其極點具有負實部。判定系統(tǒng)穩(wěn)定性的主要方法有1. 若已知系統(tǒng)的特征多項式 可以用求解方程的辦法，求出特征根，常用如下指令：solve(f(s)返回特征根的值。Solve命令必須在符號運算功能下運行。2. 設(shè)特征多項式為其降排列的系數(shù)向量為CC1C2Cn-1 則指令C1roots(C)返回特征根的數(shù)值。若使用命令poly(C1)則返回原多項式系數(shù)向量，即roots與poly互為逆操作。Roots命令不需要符號功能支持。若給出的是開環(huán)傳遞函數(shù)，需要判斷對應(yīng)的單位負反饋系統(tǒng)的穩(wěn)定性，必須用feedback(sys，1)命令求出對應(yīng)的閉環(huán)系統(tǒng)，再求閉環(huán)特征根。2 若已知傳遞函數(shù)。可以用tf2zp轉(zhuǎn)換成零、極點模型判斷其穩(wěn)定性。3 根據(jù)Nyquist穩(wěn)定性判據(jù)，作出開環(huán)系統(tǒng)的Nyquist圖，通過該系統(tǒng)圍繞(-1，j0)點的情況判斷相應(yīng)的單位負反饋閉環(huán)系統(tǒng)的穩(wěn)定性。4 Bode判據(jù)，通過開環(huán)系統(tǒng)的Bode圖判斷相應(yīng)的單位負反饋閉環(huán)系統(tǒng)穩(wěn)定性。對于復(fù)雜系統(tǒng)，可以綜合使用上述幾種方法，還可用step(sys)命令檢驗系統(tǒng)的單位階躍響應(yīng)是否收斂來映證判定結(jié)果的正確性。例4.1在系統(tǒng)的Nyquist圖上，疊加上單位圓，以便使用Nyquist判據(jù)。設(shè)開環(huán)系統(tǒng)作出其Nyquist圖并疊加上單位圓圖。由圖可見，該開環(huán)系統(tǒng)的Nyquist圖圍繞（1，j0）點逆時針及順時針各一次，因而相應(yīng)的閉環(huán)系統(tǒng)穩(wěn)定。程序見fanli024。參考程序fanli024: Nyquist圖與單位圓圖w2=linspace(0