1、控制系統(tǒng)的仿真第1頁，共100頁，2022年，5月20日，15點11分，星期三1本章教學目的及要求熟悉MATLAB的基本應用掌握利用MATLAB建立數(shù)學模型的方法熟悉控制系統(tǒng)的時域和頻域分析掌握MATLAB的仿真編程應用第7章控制系統(tǒng)的MATLAB仿真第2頁，共100頁，2022年，5月20日，15點11分，星期三2 7.1 控制系統(tǒng)的模型表示 MATLAB 提供了數(shù)學模型的建立函數(shù)和各模型之間的轉(zhuǎn)換功能函數(shù)，可以分別采用傳遞函數(shù)、零極點增益、狀態(tài)空間以及動態(tài)結(jié)構圖等4種數(shù)學模型來表示控制系統(tǒng)，前3種是用數(shù)學表達式描述的系統(tǒng)模型，每種模型都有連續(xù)系統(tǒng)及離散系統(tǒng)兩種類別的表示；而動態(tài)結(jié)構圖是基于

2、傳遞函數(shù)的圖形化形式，是采用MATLAB中提供的SIMULINK結(jié)構圖來實現(xiàn)的。MATLAB中使用的數(shù)學模型之間的轉(zhuǎn)換也很方便，使得采用MATLAB編制的程序更加簡單、精煉而高效。第7章控制系統(tǒng)的MATLAB仿真第3頁，共100頁，2022年，5月20日，15點11分，星期三3 7.1.1 系統(tǒng)的傳遞函數(shù)模型表示 傳遞函數(shù)模型通常表示線性定常時不變系統(tǒng)（LTI），可以是連續(xù)的時間系統(tǒng)，也可以是離散的時間系統(tǒng)。 對于離散的時間系統(tǒng)，其脈沖傳遞函數(shù)可表示為：第7章控制系統(tǒng)的MATLAB仿真 不論是連續(xù)的還是離散的時間系統(tǒng)，其傳遞函數(shù)的分子/分母多項式均按s或z的降冪來排列。在MATLAB中可直接采

3、用分子/分母多項式系數(shù)構成的兩個向量num與den來表示系統(tǒng)，即：第4頁，共100頁，2022年，5月20日，15點11分，星期三4 在MATLAB中，可用函數(shù)命令tf()來建立控制系統(tǒng)的傳遞函數(shù)模型，其調(diào)用格式和功能分別為：（1）sys= tf(num, den)；（2）sys= tf(num, den, Ts) ；。（3）sys= tf(M) ；（4）tfsys= tf(sys)； 第7章控制系統(tǒng)的MATLAB仿真第5頁，共100頁，2022年，5月20日，15點11分，星期三5第7章7.1.2 零極點增益模型 當連續(xù)系統(tǒng)的傳遞函數(shù)表達式采用系統(tǒng)增益、系統(tǒng)零點與系統(tǒng)極點來表示時，稱之為系統(tǒng)

4、零極點增益模型。系統(tǒng)零極點增益模型是傳遞函數(shù)模型的一種特殊形式。 離散系統(tǒng)的傳遞函數(shù)零極點增益模型： 控制系統(tǒng)的MATLAB仿真第6頁，共100頁，2022年，5月20日，15點11分，星期三6第7章 在MATLAB里，連續(xù)與離散系統(tǒng)都可直接用向量z、p、k構成的矢量組z, p, k來表示系統(tǒng)，即： 控制系統(tǒng)的MATLAB仿真第7頁，共100頁，2022年，5月20日，15點11分，星期三7第7章 在MATLAB中，可用函數(shù)命令zpk()來建立控制系統(tǒng)的零極點增益模型，其調(diào)用格式和功能分別為：（1）sys= zpk(num, den) ；（2）sys= zpk(num,den, Ts)； （3

5、）sys= zpk(M)； （4）tfsys= zpk(sys)；控制系統(tǒng)的MATLAB仿真第8頁，共100頁，2022年，5月20日，15點11分，星期三8 7.1.3 狀態(tài)空間模型 控制系統(tǒng)在主要工作區(qū)域內(nèi)的一定條件下可近似為線性時不變（LTI）模型，連續(xù)LTI對象系統(tǒng)總是能用一階微分方程組來表示，寫成矩陣形式即為狀態(tài)空間模型：第7章控制系統(tǒng)的MATLAB仿真其中： 式（a）系統(tǒng)的狀態(tài)方程，是由n個一階微分方程組成的微分方程組；式（b）系統(tǒng)的輸出方程，是由1個線性代數(shù)方程組成的；第9頁，共100頁，2022年，5月20日，15點11分，星期三9第7章離散系統(tǒng)的狀態(tài)空間模型可表示為：控制系統(tǒng)

6、的MATLAB仿真式中：U系統(tǒng)的控制輸入向量；X系統(tǒng)的狀態(tài)向量；Y系統(tǒng)的輸出向量；k特定時刻的采樣點；A狀態(tài)矩陣，由控制對象的參數(shù)決定；B控制矩陣；C輸出矩陣；D直接傳輸矩陣。第10頁，共100頁，2022年，5月20日，15點11分，星期三10 MATLAB中的函數(shù)ss()可用來建立控制系統(tǒng)的狀態(tài)空間模型，或者將傳遞函數(shù)模型與零極點增益模型轉(zhuǎn)換為系統(tǒng)狀態(tài)空間模型。ss()函數(shù)的調(diào)用格式為：（1）sys= ss(a, b, c, d)；（2）sys= ss(a, b, c, d, Ts)；（3）sys= ss(d)；該函數(shù)等價于sys= ss( , , ,d)。（4）sys_ss= ss(sy

7、s)；第7章控制系統(tǒng)的MATLAB仿真第11頁，共100頁，2022年，5月20日，15點11分，星期三11第7章【例7.1】已知某系統(tǒng)的狀態(tài)空間表達式為：控制系統(tǒng)的MATLAB仿真試采用MATLAB語言求出該系統(tǒng)的狀態(tài)空間模型。第12頁，共100頁，2022年，5月20日，15點11分，星期三12解：采用狀態(tài)空間模型表示時，可在MATLAB命令窗口中輸入以下命令： A=1 0 0 0;0 1 0 0;0 0 1 0;-1 -5 0 -2; B=0;0;0;1; C=3 2 1 0;D=0; sys=ss(A,B,C,D)上述指令執(zhí)行后可得指定系統(tǒng)的狀態(tài)空間模型為：a = x1 x2 x3 x

8、4 x1 1 0 0 0 x2 0 1 0 0 x3 0 0 1 0 x4 -1 -5 0 -2第7章控制系統(tǒng)的MATLAB仿真第13頁，共100頁，2022年，5月20日，15點11分，星期三13第7章 b = u1 x1 0 x2 0 x3 0 x4 1 c = x1 x2 x3 x4 y1 3 2 1 0 d = u1 y1 0 Continuous-time model. 控制系統(tǒng)的MATLAB仿真第14頁，共100頁，2022年，5月20日，15點11分，星期三14 7.1.4 系統(tǒng)不同模型間的相互轉(zhuǎn)換 在MATLAB 6.1的信號處理工具箱與控制系統(tǒng)工具箱中，提供了傳遞函數(shù)模型、零

9、極點增益模型與狀態(tài)空間模型之間轉(zhuǎn)換的函數(shù)：ss2tf()，ss2zp()，tf2ss()，tf2zp()，zp2ss()，zp2tf()。 這些函數(shù)之間的轉(zhuǎn)換功能如表7-1所示。 第7章控制系統(tǒng)的MATLAB仿真第15頁，共100頁，2022年，5月20日，15點11分，星期三15第7章控制系統(tǒng)的MATLAB仿真表7-1 數(shù)學模型之間的轉(zhuǎn)換函數(shù)及其功能函 數(shù) 名函 數(shù) 功 能ss2tf將系統(tǒng)狀態(tài)空間模型轉(zhuǎn)換為傳遞函數(shù)模型ss2zp將系統(tǒng)狀態(tài)空間模型轉(zhuǎn)換為零極點增益模型tf2ss將系統(tǒng)傳遞函數(shù)模型轉(zhuǎn)換為狀態(tài)空間模型tf2zp將系統(tǒng)傳遞函數(shù)模型轉(zhuǎn)換為零極點增益模型zp2ss將系統(tǒng)零極點增益模型轉(zhuǎn)換

10、為狀態(tài)空間模型zp2tf將系統(tǒng)零極點增益模型轉(zhuǎn)換為傳遞函數(shù)模型第16頁，共100頁，2022年，5月20日，15點11分，星期三16【例7.2】 已知某系統(tǒng)的傳遞函數(shù)為：第7章控制系統(tǒng)的MATLAB仿真 試用MATLAB語言求出該系統(tǒng)的傳遞函數(shù)模型、狀態(tài)空間模型和零極點增益模型。解：（1）求系統(tǒng)的傳遞函數(shù)模型在MATLAB命令窗口輸入以下命令：num=12 24 12 20;den=2 4 6 2 2;sys=tf(num, den)第17頁，共100頁，2022年，5月20日，15點11分，星期三17第7章執(zhí)行以上語句后可得系統(tǒng)的傳遞函數(shù)模型為：12 s3 + 24 s2 + 12 s +

11、20-2 s4 + 4 s3 + 6 s2 + 2 s + 2（2）求系統(tǒng)的狀態(tài)空間模型該系統(tǒng)的狀態(tài)空間模型可以通過MATLAB的模型轉(zhuǎn)換函數(shù)來完成。在MATLAB命令窗口輸入以下命令：a,b,c,d=tf2ss(num,den);sys=ss(a,b,c,d) 控制系統(tǒng)的MATLAB仿真第18頁，共100頁，2022年，5月20日，15點11分，星期三18執(zhí)行完上述語句后，可得系統(tǒng)的狀態(tài)空間模型的狀態(tài)矩陣a，系統(tǒng)控制矩陣b，系統(tǒng)輸出矩陣c，系統(tǒng)直接傳輸矩陣d分別為：a = x1 x2 x3 x4 x1 -2 -3 -1 -1 x2 1 0 0 0 x3 0 1 0 0 x4 0 0 1 0第

12、7章控制系統(tǒng)的MATLAB仿真第19頁，共100頁，2022年，5月20日，15點11分，星期三19第7章 b = u1 x1 1 x2 0 x3 0 x4 0 c = x1 x2 x3 x4 y1 6 12 6 10 d = u1 y1 0 Continuous-time model.控制系統(tǒng)的MATLAB仿真第20頁，共100頁，2022年，5月20日，15點11分，星期三20由以上數(shù)據(jù)可寫出系統(tǒng)的狀態(tài)空間模型為：第7章控制系統(tǒng)的MATLAB仿真第21頁，共100頁，2022年，5月20日，15點11分，星期三21第7章（3）求系統(tǒng)的零極點增益模型該系統(tǒng)的零極點增益模型也可以通過MATLA

13、B的模型轉(zhuǎn)換函數(shù)來完成。在MATLAB命令窗口輸入以下命令： z, p, k=tf2zp(num, den); sys= zpk(z, p, k)執(zhí)行以上語句后可得系統(tǒng)的零極點增益模型為：Zero/pole/gain: 6 (s+1.929) (s2 + 0.07058s + 0.8638)-(s2 + 0.08663s + 0.413) (s2 + 1.913 s + 2.421)控制系統(tǒng)的MATLAB仿真第22頁，共100頁，2022年，5月20日，15點11分，星期三227.2 環(huán)節(jié)方框圖模型的化簡 7.2.1 環(huán)節(jié)串聯(lián)連接的化簡 多個環(huán)節(jié)串聯(lián)的連接形式是控制系統(tǒng)最基本的組成結(jié)構形式之一

14、?？刂葡到y(tǒng)的環(huán)節(jié)串聯(lián)及其化簡就是模塊方框圖模型的串聯(lián)及其化簡。可以用MATLAB的函數(shù)命令series()將串聯(lián)模塊進行等效變換。 使用series()函數(shù)命令不必做多項式的乘除運算即可實現(xiàn)兩個環(huán)節(jié)傳遞函數(shù)的串聯(lián)連接。如果令sys1= tf(num1,den1)，sys2= tf(num2,den2)，其命令格式為：sys= series(sysl,sys2)第7章控制系統(tǒng)的MATLAB仿真第23頁，共100頁，2022年，5月20日，15點11分，星期三23如果已知兩個環(huán)節(jié)的狀態(tài)空間模型矩陣組分別為：（a1,b1,c1,d1）與（a2,b2,c2,d2），則求兩個環(huán)節(jié)串聯(lián)連接等效系統(tǒng)狀態(tài)空間

15、模型a,b,c,d矩陣組的命令格式為：a,b,c,d= series(a1,bl,cl,dl,a2,b2,c2,d2) series()函數(shù)命令還可以將多個環(huán)節(jié)按兩兩串聯(lián)的形式多次遞歸調(diào)用加以連接，進行等效化簡。 sys= series(sysl,sys2)命令可以用命令sys=sys1*sys2*sysn取代，不僅省掉“series()”字符，且可以實現(xiàn)多個環(huán)節(jié)的串聯(lián)等效傳遞函數(shù)的求取。第7章控制系統(tǒng)的MATLAB仿真第24頁，共100頁，2022年，5月20日，15點11分，星期三24第7章【例7.3】 已知雙閉環(huán)調(diào)速系統(tǒng)電流環(huán)內(nèi)的前向通道3個模塊傳遞函數(shù)分別為：控制系統(tǒng)的MATLAB仿真

16、試求串聯(lián)連接的等效傳遞函數(shù)及其等效狀態(tài)空間模型。第25頁，共100頁，2022年，5月20日，15點11分，星期三25解：（1）根據(jù)MATLAB程序設計的基本方法和函數(shù)命令series()，可以編寫出MATLAB程序如下：n1=0.0128 1;d1=0.04 0;sys1=tf(n1,d1);n2=30;d2=0.00167 1;sys2=tf(n2,d2);n3=2.5;d3=0.0128 1;sys3=tf(n3,d3);sys=sys1*sys2*sys3s1=ss(sys1);s2=ss(sys2);s3=ss(sys3);sys12=series(s1,s2);sys123=ser

17、ies(sys12,s3)第7章控制系統(tǒng)的MATLAB仿真第26頁，共100頁，2022年，5月20日，15點11分，星期三26第7章控制系統(tǒng)的MATLAB仿真（2）在MATLAB命令窗口輸入程序名，程序運行后得到如下電流環(huán)內(nèi)前向通道的等效傳遞函數(shù)及等效狀態(tài)空間模型：Transfer function: 0.96 s + 75-8.55e-007 s3 + 0.0005788 s2 + 0.04 s a = x1 x2 x3 x1 -78.13 2246 0 x2 0 -598.8 800 x3 0 0 0 第27頁，共100頁，2022年，5月20日，15點11分，星期三27 b = u1

18、x1 0 x2 40.96 x3 4 c = x1 x2 x3 y1 12.21 0 0 d = u1 y1 0Continuous-time model.第7章控制系統(tǒng)的MATLAB仿真第28頁，共100頁，2022年，5月20日，15點11分，星期三28第7章 7.2.2 環(huán)節(jié)并聯(lián)連結(jié)的化簡 環(huán)節(jié)并聯(lián)是指多個環(huán)節(jié)的輸入信號相同，所有環(huán)節(jié)輸出的代數(shù)和為其總輸出。 采用parallel()函數(shù)命令可以等效化簡兩個環(huán)節(jié)的并聯(lián)連接。parallel()函數(shù)命令調(diào)用格式為：num，den= parallel（numl,denl,num2,den2）該命令已由命令sys=sysl+sys2+sysn所

19、取代，省掉了“parallel()”字符，且可以實現(xiàn)多個環(huán)節(jié)的并聯(lián)等效處理。parallel()函數(shù)命令調(diào)用格式還有：a,b,c,d= parallel(a1,bl,cl,dl,a2,b2,c2,d2) 控制系統(tǒng)的MATLAB仿真第29頁，共100頁，2022年，5月20日，15點11分，星期三29第7章【例7.4】 已知兩個環(huán)節(jié)的傳遞函數(shù)分別為： 試求兩環(huán)節(jié)并聯(lián)連接等效傳遞函數(shù)的num與den向量及等效的狀態(tài)空間模型?？刂葡到y(tǒng)的MATLAB仿真第30頁，共100頁，2022年，5月20日，15點11分，星期三30第7章解：（1）采用函數(shù)命令parallel()，給出MATLAB程序如下：nu

20、m1=1;den1=1 2;sys1=tf(num1,den1);num2=2 1;den2=1 1 2;sys2=tf(num2,den2);s1=ss(sys1);s2=ss(sys2);sys=sys1+sys2sys12=parallel(s1,s2)（2）在MATLAB命令窗口運行程序后得到等效傳遞函數(shù)及等效狀態(tài)空間模型為：Transfer function: 3 s2 + 6 s + 4-s3 + 3 s2 + 4 s + 4 控制系統(tǒng)的MATLAB仿真第31頁，共100頁，2022年，5月20日，15點11分，星期三31第7章a = x1 x2 x3 x1 -2 0 0 x2 0

21、 -1 -1 x3 0 2 0b = u1 x1 1 x2 2 x3 0c = x1 x2 x3 y1 1 1 0.25d = u1 y1 0控制系統(tǒng)的MATLAB仿真第32頁，共100頁，2022年，5月20日，15點11分，星期三32 7.2.3 環(huán)節(jié)反饋連結(jié)的化簡 利用MATLAB中的feedback()函數(shù)命令可將兩個環(huán)節(jié)按反饋形式進行連接后求其等效傳遞函數(shù)。feedback()函數(shù)命令形式：（1）G=feedback(G1,G2,sign)； （2）a,b,c,d=feedback(a1,b1,c1,d1,a2,b2,c2,d2,sign)（3）num,den=feedback(nu

22、m1,den1,num2,den2,sign)；第7章控制系統(tǒng)的MATLAB仿真第33頁，共100頁，2022年，5月20日，15點11分，星期三33第7章 【例7.5】已知晶閘管-直流電機單閉環(huán)調(diào)速系統(tǒng)的動態(tài)結(jié)構圖如圖7-4所示，求該閉環(huán)系統(tǒng)的傳遞函數(shù)?？刂葡到y(tǒng)的MATLAB仿真圖7-4 直流單閉環(huán)調(diào)速系統(tǒng)第34頁，共100頁，2022年，5月20日，15點11分，星期三34第7章解：（1）給出MATLAB程序如下：n1=1;d1=0.017 1;s1=tf(n1,d1);n2=1;d2=0.075 0;s2=tf(n2,d2);sys1=feedback(s1*s2,1)n3=0.049

23、1;d3=0.088 0;s3=tf(n3,d3);n4=44;d4=0.00167 1;s4=tf(n4,d4);n5=1;d5=0.1925;s5=tf(n5,d5);n6=0.01178;d6=1;s6=tf(n6,d6);sysq=sys1*s3*s4*s5;sys=feedback(sysq,s6)控制系統(tǒng)的MATLAB仿真第35頁，共100頁，2022年，5月20日，15點11分，星期三35第7章 （2）程序運行結(jié)果為：Transfer function: 1-0.001275 s2 + 0.075 s + 1 Transfer function: 2.156 s + 44-3.6

24、07e-008 s4 + 2.372e-005 s3 + 0.001299 s2+ 0.04234 s + 0.5183 控制系統(tǒng)的MATLAB仿真第36頁，共100頁，2022年，5月20日，15點11分，星期三36第7章由以上運算數(shù)據(jù)可得單閉環(huán)系統(tǒng)的小閉環(huán)的傳遞函數(shù)為：單閉環(huán)系統(tǒng)的閉環(huán)的傳遞函數(shù)（略去分母的項）為：控制系統(tǒng)的MATLAB仿真第37頁，共100頁，2022年，5月20日，15點11分，星期三37第7章7.3 控制系統(tǒng)時域分析的MATLAB實現(xiàn)7.3.1 概述1系統(tǒng)時域響應仿真的主要問題時域分析法是根據(jù)系統(tǒng)微分方程用拉普拉斯變換求解系統(tǒng)的動態(tài)響應。典型的動態(tài)過程響應有單位階躍響

25、應、單位斜坡響應、單位加速度響應與單位脈沖響應等。時域分析的另外一個目的是求解響應的性能指標。通常將控制系統(tǒng)跟蹤或復現(xiàn)階躍輸入信號響應的指標作為系統(tǒng)控制性能的指標。階躍響應的一般性能指標有：峰值時間tp、超調(diào)量%、調(diào)節(jié)時間ts、及穩(wěn)態(tài)誤差ess。 控制系統(tǒng)的MATLAB仿真第38頁，共100頁，2022年，5月20日，15點11分，星期三38第7章控制系統(tǒng)的MATLAB仿真 2時域分析的MATLAB實現(xiàn) 時域響應的MATLAB仿真方法有兩種，一種是在MATLAB的函數(shù)指令方式下進行時域仿真；另一種是在SIMULINK窗口菜單操作方式下進行時域仿真。 對于連續(xù)系統(tǒng)，可利用MATLAB提供的求取單

26、位階躍響應函數(shù)step()、單位脈沖響應函數(shù)impulse()、零輸入響應函數(shù)initial()等處理系統(tǒng)的響應。 對于離散系統(tǒng)，可利用MATLAB提供的單位階躍響應函數(shù)dstep()、單位脈沖響應函數(shù)dimpulse()、零輸入響應函數(shù)dinitial()等來處理。 第39頁，共100頁，2022年，5月20日，15點11分，星期三39第7章 7.3.2 連續(xù)系統(tǒng)的單位階躍響應 利用step()函數(shù)可求出連續(xù)系統(tǒng)的單位階躍響應，該命令的調(diào)用格式和功能如下：（1）step(sys)； （2）step(sys，t)；（3）step(sysl,sys2,，sysn)；（4）y，t，x=step(s

27、ys)； 求出了系統(tǒng)的階躍響應曲線，要分析系統(tǒng)的穩(wěn)定性能，還需要計算系統(tǒng)的性能指標。包括峰值時間tp、超調(diào)量%、調(diào)節(jié)時間ts、及穩(wěn)態(tài)誤差ess等。這些參數(shù)可以直接從圖上讀取，也可通過MABLAB的M函數(shù)來求取。 控制系統(tǒng)的MATLAB仿真第40頁，共100頁，2022年，5月20日，15點11分，星期三40第7章下面介紹利用M函數(shù)來計算性能指標。定義該指標函數(shù)為terg函數(shù)，其程序語句如下：function mp,tp,ess,b1,b2,sigma,n,pusi,T,f=targ(y,t) mp,tf=max(y); ct=length(t);mp; tp=t(tf);tm=max(tf);

28、yss=(ct);q=1;m=q-1;while m3,控制系統(tǒng)的MATLAB仿真第41頁，共100頁，2022年，5月20日，15點11分，星期三41第7章 for a=(tm/100):0.01:tm j=0:a:tm; for i=1:length(j); if (y(i+1)-y(i)0 m=m+1; pm(m)=y(i); tp(m)=t(i); end end endendyss=y(ct);ess=1-yss;控制系統(tǒng)的MATLAB仿真第42頁，共100頁，2022年，5月20日，15點11分，星期三42第7章 b1=pm(1)-yss;b2=pm(2)-yss;sigma=10

29、0*b1/yss;n=b1/b2;pusi=(b1-b2)/b1;T=(tp(2)-tp(1);f=1/Ttp=tp(1); 上述程序中：mp為曲線的峰值；tp為峰值時間；ess 為階躍響應的余差；b1為第一峰值的超調(diào)量；b2為第二峰值的超調(diào)量；sigma為階躍響應的超調(diào)量；n為階躍響應的衰減比；pusi為階躍響應的衰減率；T為衰減振蕩周期；f振蕩頻率?？刂葡到y(tǒng)的MATLAB仿真第43頁，共100頁，2022年，5月20日，15點11分，星期三43第7章7.3.3 連續(xù)系統(tǒng)的單位脈沖響應求連續(xù)系統(tǒng)單位脈沖響應的函數(shù)為impulse()，該函數(shù)命令的調(diào)用格式和功能如下： （1）impulse(s

30、ys)；（2）impulse(sys，t)；（3）impulse(sysl,sys2,，sysN)； （4）y，t，x=impulse(sys)； 【例7.6】已知單位負反饋系統(tǒng)前向通道的傳遞函數(shù)為： ，試繪出系統(tǒng)的單位階躍響應和脈沖響應曲線。并計算系統(tǒng)的階躍響應性能指標?？刂葡到y(tǒng)的MATLAB仿真第44頁，共100頁，2022年，5月20日，15點11分，星期三44第7章解：（1）求階躍響應曲線在命令窗口輸入：num=4; den=1 1 0;sys=tf(num,den);closys=feedback(sys,1);% 求閉環(huán)傳遞函數(shù)step(closys)可得如圖7-5所示的階躍響應曲

31、線。控制系統(tǒng)的MATLAB仿真第45頁，共100頁，2022年，5月20日，15點11分，星期三45第7章控制系統(tǒng)的MATLAB仿真圖7-5 系統(tǒng)的階躍響應曲線第46頁，共100頁，2022年，5月20日，15點11分，星期三46第7章 （2）求脈沖響應曲線在命令窗口繼續(xù)輸入：impulse(closys) ；可得如圖7-6所示的脈沖響應曲線。控制系統(tǒng)的MATLAB仿真圖7-6 系統(tǒng)的脈沖響應曲線第47頁，共100頁，2022年，5月20日，15點11分，星期三47第7章 （3）求系統(tǒng)的階躍響應性能指標 首先在MATLAB程序編輯窗口編輯前面給出的targ函數(shù)，并保存函數(shù)名為targ.m。然后

32、繼續(xù)在MATLAB命令窗口輸入以下語句：y,t=step(slosys);mp,tp,ess,b1,b2,sigma,n,pusi,T,f=targ(y,t)運行完上述語句后，可得以下性能指標：mp = 1.4441 ；tp = 1.6062ess = -0.0017 ；b1 = 0.4425b2 = 0.0857 ；sigma = 44.1728n = 5.1651 ；pusi = 0.8064T = 3.2125 ；f = 0.3113控制系統(tǒng)的MATLAB仿真第48頁，共100頁，2022年，5月20日，15點11分，星期三48第7章7.3.4 連續(xù)系統(tǒng)的零輸入響應求連續(xù)系統(tǒng)零輸入響應的

33、函數(shù)為initial()，該函數(shù)命令的格式如下：（1）y，t，x= initial(sys)；（2）initial(sys,x0)； （3）initial(sys,x0, t)； （4）initial(sysl,sys2,，sysN,x0)?？刂葡到y(tǒng)的MATLAB仿真第49頁，共100頁，2022年，5月20日，15點11分，星期三49第7章 7.4 控制系統(tǒng)頻域分析的MATLAB實現(xiàn)7.4.1 概述1控制系統(tǒng)頻域響應仿真的主要問題 對于頻域分析，MATLAB仿真就是用先進的計算機仿真技術來精確繪制Bode圖、Nyquist曲線圖、Nichols曲線圖等三種曲線，并計算系統(tǒng)的頻域性能指標：剪切

34、頻率（也叫截止頻率）、- 穿越頻率、相角穩(wěn)定裕度、幅值穩(wěn)定裕度，以便研究系統(tǒng)控制過程的穩(wěn)定性、快速性及穩(wěn)態(tài)精度等性能。控制系統(tǒng)的MATLAB仿真第50頁，共100頁，2022年，5月20日，15點11分，星期三50第7章 2MATLAB中頻域分析的實現(xiàn)方法 頻域分析在經(jīng)典控制里占有重要的地位。借助于MATLAB軟件，運行它所提供的頻率分析函數(shù)，能夠方便、簡單、快捷地繪制Bode圖、Nyquist曲線和Nichols圖，并計算出頻域性能指標，還可以借助于這些曲線對系統(tǒng)進行分析。 應用MATLAB提供的上述頻域仿真函數(shù)與其他函數(shù)命令、語句，編制成MATLAB程序，這種在MATLAB的指令方式下進行

35、頻域仿真同樣是最常用的實現(xiàn)方法?？刂葡到y(tǒng)的MATLAB仿真第51頁，共100頁，2022年，5月20日，15點11分，星期三51第7章 7.4.2 求連續(xù)系統(tǒng)的Bode圖1求連續(xù)系統(tǒng)Bode圖的函數(shù) 在MATLAB中求連續(xù)系統(tǒng)Bode圖的函數(shù)為bode()，該函數(shù)命令的調(diào)用格式和功能為：（1）mag,phase,w=bode(sys)； （2）bode(sys)；（3）bode(sys,w)； （4）bode(sysl,sys2,sysn)或bode(sysl,sys2,，sysn,w)； 控制系統(tǒng)的MATLAB仿真第52頁，共100頁，2022年，5月20日，15點11分，星期三52第7章

36、2求離散系統(tǒng)Bode圖的函數(shù)dbode()求離散系統(tǒng)Bode圖的函數(shù)為dbode()，該函數(shù)命令的調(diào)用格式為：（1）mag,phase,w=dbode(a,b,c,d,Ts)（2）mag,phase,w=dbode(a,b,c,d,Ts,iu) （3）mag,phase,w=dbode(a,b,c,d,Ts,iu,w)或 mag,phase,w=dbode(num,den,Ts,w) （4）mag,phase,w=dbode(num,den,Ts) 控制系統(tǒng)的MATLAB仿真第53頁，共100頁，2022年，5月20日，15點11分，星期三53第7章 7.4.3 繪制系統(tǒng)Nyquist曲線的函

37、數(shù) 計算或繪制系統(tǒng)Nyquist曲線的函數(shù)為nyquist()，該函數(shù)命令的調(diào)用格式為：（1）re,im,w=nyquist(sys)（2）nyquist(sys)（3）nyquist(sys,w)（4）nyquist(sysl,sys2,sysn) （5）nyquist(sysl,sys2,sysn,w)控制系統(tǒng)的MATLAB仿真第54頁，共100頁，2022年，5月20日，15點11分，星期三54第7章 如果系統(tǒng)為離散系統(tǒng)，則相應的nyquist函數(shù)為 dnyquist()，函數(shù)命令調(diào)用為：（1）re,im,w= dnyquist(a,b,c,d,Ts)（2）re,im,w= dnyqui

38、st(a,b,c,d,Ts,iu)（3）re,im,w= dnyquist(a,b,c,d,Ts,iu,w)（4）re,im,w= dnyquist(num,den,Ts)（5）re,im,w= dnyquist(num,den,Ts,w)控制系統(tǒng)的MATLAB仿真第55頁，共100頁，2022年，5月20日，15點11分，星期三55第7章 7.4.4 求連續(xù)系統(tǒng)Nichols曲線的函數(shù) 求連續(xù)系統(tǒng)Nichols曲線的函數(shù)為nichols()，該函數(shù)命令的調(diào)用格式為：（1）mag,phase,w= nichols(sys)（2）nichols(sys)（3）nichols(sys,w)（4）n

39、ichols(sysl,sys2,sysn)（5）nichols(sysl,sys2,，sysn,w)控制系統(tǒng)的MATLAB仿真第56頁，共100頁，2022年，5月20日，15點11分，星期三56第7章控制系統(tǒng)的MATLAB仿真 7.4.5 求系統(tǒng)幅值裕度與相位裕度 求系統(tǒng)幅值裕度與相位裕度的函數(shù)為margin()，該函數(shù)命令的調(diào)用格式為：（1）Gm,Pm,Wcp,Wcg= margin(sys)； （2）Gm,Pm,Wcp,Wcg= margin(mag,phase,w)（3）margin(sys)； （4）margin(mag,phase,w)； 對于離散時間系統(tǒng)，可以先用函數(shù)dbode

40、()計算頻率響應，然后再調(diào)用margin()函數(shù)，即用margin(mag,phase,w)函數(shù)命令，再求系統(tǒng)性能指標或繪制Bode圖。 第57頁，共100頁，2022年，5月20日，15點11分，星期三57第7章控制系統(tǒng)的MATLAB仿真 7.5 利用MATLAB語言進行編程仿真的應用7.5.1 面向傳遞函數(shù)的線性系統(tǒng)仿真 對線性連續(xù)系統(tǒng)進行單輸入、單輸出的仿真，可以將用戶輸入的系統(tǒng)傳遞函數(shù)模型轉(zhuǎn)化為仿真計算模型，再應用數(shù)值積分法進行處理。采用四階龍格庫塔法保證系統(tǒng)仿真過程中具備一定的精度和性能指標要求，并可實現(xiàn)重復運行，便于研究參數(shù)的變化對系統(tǒng)動態(tài)性能的影響，且運行過程直觀、形象，修改參數(shù)

41、容易。第58頁，共100頁，2022年，5月20日，15點11分，星期三58第7章控制系統(tǒng)的MATLAB仿真 1仿真系統(tǒng)的典型結(jié)構框圖如圖7-11所示的典型閉環(huán)系統(tǒng)結(jié)構，有：圖7-11 典型閉環(huán)系統(tǒng)結(jié)構圖 第59頁，共100頁，2022年，5月20日，15點11分，星期三59第7章控制系統(tǒng)的MATLAB仿真 2仿真程序框圖與實現(xiàn)（1）程序框圖的描述 如圖7-12所示。（2）程序設計的實現(xiàn) 用MATLAB語言編寫的面向傳遞函數(shù)的線性系統(tǒng)仿真程序定義為CSS1(Control System Simulation 1) 。第60頁，共100頁，2022年，5月20日，15點11分，星期三60第7章圖

42、7-12 面向傳遞函數(shù)的線性系統(tǒng)仿真框圖第61頁，共100頁，2022年，5月20日，15點11分，星期三61第7章控制系統(tǒng)的MATLAB仿真設源程序名為CSS1.m，保存為M文件。1）輸入數(shù)據(jù)a=;% n+1維分母系數(shù)向量c=;% m+1維分子系數(shù)向量 X0=;% 狀態(tài)向量初值V=;% 反饋系數(shù)n=;% 系統(tǒng)階次T0=;% 起始時間Tf=;% 終止時間h=;% 計算步長R=r;% 階躍輸入函數(shù)幅值 第62頁，共100頁，2022年，5月20日，15點11分，星期三62第7章控制系統(tǒng)的MATLAB仿真 2）形成開、閉環(huán)系數(shù)陣c=c/a(1); % 首一化處理a=a/a(1);A=a(2:n+1

43、); A=rot90(rot90(eye(n-1,n);-fliplr(A);% 形成能控標準形B=zeros(1,n-1),1;% 形成輸入陣Bm1=length(c);% 分子系數(shù)向量維數(shù)m+1C=fliplr(c),zeros(1,n-m1);% 形成輸出陣CAb=A-B*C*V;% 形成閉環(huán)系數(shù)陣AbX=X0;y=0;t=T0;% 設初值，準備遞推運算第63頁，共100頁，2022年，5月20日，15點11分，星期三63第7章控制系統(tǒng)的MATLAB仿真 3）運算求解N=round(Tf-T0)/h;% 確定輸出點數(shù)for i=1:N% 四階龍格-庫塔法K1=Ab*X+B*R;K2=Ab

44、*(X+h*K1/2)+B*R;K3=Ab*(X+h*K2/2)+B*R;K4=Ab*(X+h*K3)+B*R;% 求各次斜率KX=X+h*(K1+2*K2+2*K3+K4)/6;% 求狀態(tài)y=y,C*X;% 求輸出并以向量形式保存 t=t,t(i)+h;% 輸出對應時刻以向量形式保存end第64頁，共100頁，2022年，5月20日，15點11分，星期三64第7章控制系統(tǒng)的MATLAB仿真 4）輸出結(jié)果t,y % 輸出數(shù)據(jù)形式結(jié)果plot(t,y)% 輸出曲線形式結(jié)果 以上即為采用MATLAB語言實現(xiàn)面向典型閉環(huán)系統(tǒng)數(shù)字仿真的程序，可編輯為CSS1.m文件存儲起來。使用時，只要進入到MATL

45、AB命令窗口，按要求輸入系統(tǒng)參數(shù)和運行參數(shù)，再調(diào)用該文件，即可得到系統(tǒng)的仿真結(jié)果。第65頁，共100頁，2022年，5月20日，15點11分，星期三65第7章控制系統(tǒng)的MATLAB仿真 該程序中用到了幾個MATLAB特殊功能函數(shù)，簡要說明如下：eye(m,n)產(chǎn)生mn型I陣，即主對角元素均為1，其余元素為零的長方陣length( )求取括號內(nèi)向量的維數(shù)zeros(m,n)產(chǎn)生mn型零陣rot90( )將括號內(nèi)矩陣左旋90fliplr( )將括號內(nèi)矩陣或向量左右翻轉(zhuǎn)round( )對括號內(nèi)運算取整plot(t,y)打印以t為橫軸，y為縱軸的y(t)曲線 第66頁，共100頁，2022年，5月20

46、日，15點11分，星期三66第7章控制系統(tǒng)的MATLAB仿真 7.5.2 面向結(jié)構圖的線性系統(tǒng)仿真1仿真程序框圖與程序?qū)崿F(xiàn)（1）程序框圖的描述 當控制系統(tǒng)的數(shù)學模型以動態(tài)結(jié)構圖的形式給出時，可以用以下的方法處理，如圖7-15所示。（2）程序設計的實現(xiàn) 采用MATLAB語言編程，程序名為CSS2.M，保存為M文件。第67頁，共100頁，2022年，5月20日，15點11分，星期三67圖7-15 面向結(jié)構圖的線性系統(tǒng)的仿真框圖第7章第68頁，共100頁，2022年，5月20日，15點11分，星期三68第7章控制系統(tǒng)的MATLAB仿真輸入系統(tǒng)原始數(shù)據(jù)：P % 輸入各環(huán)節(jié)參數(shù)WIJ=% 輸入連接陣非零

47、元素n=;% 環(huán)節(jié)個數(shù)（系統(tǒng)階次）Y0=;% 階躍輸入幅值Yt0=% 各環(huán)節(jié)輸出初值h=;% 計算步長L1=;%每隔L1個點打印一次 T0=;% 起始時間 Tf=;% 終止時間nout=Nout;% 輸出環(huán)節(jié)編號 第69頁，共100頁，2022年，5月20日，15點11分，星期三69第7章控制系統(tǒng)的MATLAB仿真程序處理開始：A=diag(P(:,1);B=diag(P(:,2);% 求開環(huán)系統(tǒng)系數(shù)陣C=diag(P(:,3);D=diag(P(:,4);m=length(WIJ(:,1);% 求非零元素個數(shù)W0=zeros(n,1);W=zeros(n,n);% 建立初始W、W0陣for

48、k=1:mif (WIJ(k,2)=0);W0(WIJ(k,1)=WIJ(k,3);% 求W0陣else W(WIJ(k,1),WIJ(k,2)=WIJ(k,3);% 求W陣endend第70頁，共100頁，2022年，5月20日，15點11分，星期三70第7章控制系統(tǒng)的MATLAB仿真 Q=B-D*W;Qn=inv(Q);% 求Q和Q逆陣R=C*W-A;V1=C*W0;% 求R,V1陣Ab=Qn*R;b1=Qn*V1;% 形成閉環(huán)系數(shù)陣Y=Yt0; y=Y(nout);t=T0;% 置初值N=round(Tf-T0)/(h*L1);% 總輸出點數(shù)for i=1:N; % 每循環(huán)一次，輸出一點

49、數(shù)據(jù)for j=1:L1;% 每輸出點之間計算L1次K1=Ab*Y+b1*Y0; K2=Ab*(Y+h*K1/2)+b1*Y0;K3=Ab*(Y+h*K2/2)+b1*Y0;K4=Ab*(Y+h*K3)+b1*Y0;Y=Y+h*(K1+2*K2+2*K3+K4)/6;% 龍格庫塔法求解方程end第71頁，共100頁，2022年，5月20日，15點11分，星期三71第7章控制系統(tǒng)的MATLAB仿真y=y,Y(nout);% 保存輸出環(huán)節(jié)動態(tài)響應值t=t,t(i)+h*L1;% 保存時間向量endt,y% 輸出響應數(shù)據(jù)plot(t,y)% 輸出響應圖形 程序中用到MATLAB特殊功能函數(shù)diag(

50、 )，其作用是產(chǎn)生以括號內(nèi)向量元素為主對角元素的對角陣；P(:,I)用于取P陣中第I列的所有元素作為單獨列向量。第72頁，共100頁，2022年，5月20日，15點11分，星期三72第7章控制系統(tǒng)的MATLAB仿真7.5.3 線性系統(tǒng)的離散相似法仿真1仿真框圖描述 離散相似法仿真線性系統(tǒng)的程序框圖如圖7-18所示。2仿真程序設計 本程序定義為CSS3.M，輸入系統(tǒng)各環(huán)節(jié)參數(shù)陣P，連接陣WIJ，運行參數(shù)n,h,L1,T0,Tf，階躍輸入幅值Y0，輸出環(huán)節(jié)編號nout等，方法同CSS2程序。第73頁，共100頁，2022年，5月20日，15點11分，星期三73圖7-18 離散相似法仿真線性系統(tǒng)程序

51、框圖第7章第74頁，共100頁，2022年，5月20日，15點11分，星期三74第7章控制系統(tǒng)的MATLAB仿真CSS3仿真程序設計如下：A=diag(P(:,1);B=diag(P(:,2);% 求開環(huán)系統(tǒng)系數(shù)陣C=diag(P(:,3);D=diag(P(:,4);m=length(WIJ(:,1);% 求非零元素個數(shù)W0=zeros(n,1);W=zeros(n,n);% 建立初始W、W0陣for k=1:mif (WIJ(k,2)=0);W0(WIJ(k,1)=WIJ(k,3);% 求W0陣else W(WIJ(k,1),WIJ(k,2)=WIJ(k,3);% 求W陣endend第75

52、頁，共100頁，2022年，5月20日，15點11分，星期三75第7章控制系統(tǒng)的MATLAB仿真for i=1:nif (A(i,i)=0);% 求積分或積分比例環(huán)節(jié)各系數(shù)FI(i)=1;FIM(i)=h*C(i,i)/B(i,i);FIJ(i)=h*h*C(i,i)/B(i,i)/2;FIC(i)=1;FID(i)=0;% 求積分環(huán)節(jié)系數(shù)if (D(i,i)=0);% 若為積分比例環(huán)節(jié)，修正系數(shù)FID(i)=D(i,i)/B(i,i);elseendelse% 求慣性或慣性比例環(huán)節(jié)各系數(shù)FI(i)=exp(-h*A(i,i)/B(i,i);FIM(i)=(1-FI(i)*C(i,i)/A(i

53、,i);FIJ(i)=h*C(i,i)/A(i,i)-FIM(i)*B(i,i)/A(i,i);第76頁，共100頁，2022年，5月20日，15點11分，星期三76第7章控制系統(tǒng)的MATLAB仿真FIC(i)=1;FID(i)=0;% 求慣性環(huán)節(jié)系數(shù)if (D(i,i)=0);% 若為慣性比例環(huán)節(jié)，修正系數(shù)FIC(i)=C(i,i)/D(i,i)-A(i,i)/B(i,i);FID(i)=D(i,i)/B(i,i);elseendendendY=zeros(n,1);X=Y;y=0;Uk=zeros(n,1);Ub=Uk;t=T0;% 置初值 N=round(Tf-T0)/(h*L1);%

54、計算總輸出點數(shù) 第77頁，共100頁，2022年，5月20日，15點11分，星期三77第7章控制系統(tǒng)的MATLAB仿真for k=1:Nfor l=1:L1Ub=Uk; % 保存前一次輸入值Uk=W*Y+W0*Y0; % 求當前UkUdot=(Uk-Ub)/h; % 求當前Uk的導數(shù)Uf=2*Uk-Ub; % 求下一步的輸入UfX=FI.*X+FIM.*Uk+FIJ.*Udot;% 求狀態(tài)向量Y=FIC.*X+FID.*Uf;% 求輸出endy=y,Y(nout);t=t,t(k)+h*L1;endt,yplot(t,y)第78頁，共100頁，2022年，5月20日，15點11分，星期三78第

55、7章控制系統(tǒng)的MATLAB仿真該程序用到了幾個MATLAB語言中的特殊運算符： 冒號符，用于生成行向量；= = 邏輯等號符，MATLAB中的規(guī)定寫法； = 邏輯不等號符，MATLAB中的規(guī)定寫法； * 點乘運算符，MATLAB中的規(guī)定寫法，表示兩個同維矩陣中各相應位置上的元素相乘運算。第79頁，共100頁，2022年，5月20日，15點11分，星期三79第7章控制系統(tǒng)的MATLAB仿真7.5.4 非線性系統(tǒng)離散相似法仿真1典型非線性環(huán)節(jié)的求解（1）飽和非線性環(huán)節(jié)：采用MATLAB語言來編制M函數(shù)如下，函數(shù)名為satu( )：function Uc=satu(Ur,S1)if(abs(Ur)=S

56、1)if(Ur0)Uc=S1;else Uc=-S1;endelse Uc=Ur;end 第80頁，共100頁，2022年，5月20日，15點11分，星期三80第7章控制系統(tǒng)的MATLAB仿真（2）死區(qū)非線性：采用MATLAB語言來編制M函數(shù)如下，函數(shù)名為dead( )：function Uc=dead(Ur,S1)if (abs(Ur)=S1)if(Ur0)Uc=Ur-S1;else Uc=Ur+S1;endelse Uc=0;end第81頁，共100頁，2022年，5月20日，15點11分，星期三81第7章控制系統(tǒng)的MATLAB仿真（3）滯環(huán)非線性：M函數(shù)名為backlash ( )：fu

57、nction Uc,Ubb=backlash(Urb,Ur,Ucb,S1)if (UrUrb)if(Ur-S1)=Ucb)Uc=Ur-S1；else Uc=Ucb；else if(UrUrb) if(Ur-S1)=Ucb)Uc=Ur+S1；else Uc=Ucb；else Uc=Ucb；endendUbb=Ur； 第82頁，共100頁，2022年，5月20日，15點11分，星期三82第7章控制系統(tǒng)的MATLAB仿真非線性系統(tǒng)的離散相似法仿真程序命名為css4.m。根據(jù)上述分析可以畫出CSS4程序框圖如7-21所示。CSS4程序設計：輸入數(shù)據(jù)：輸入系統(tǒng)各環(huán)節(jié)參數(shù)陣P，連接陣WIJ，運行參數(shù)n,h

58、,L1,T0,Tf,階躍輸入幅值Y0，輸出環(huán)節(jié)編號nout(方法同CSS3)，以及非線性標志向量Z，參數(shù)向量S。程序開始：求開環(huán)各系數(shù)陣A,B,C,D 及連接關系矩陣 W,W0 (方法同CSS3)；求取各環(huán)節(jié)離散系數(shù)FI(i),FIM(i),FIJ(i),FIC(i),FID(i)。第83頁，共100頁，2022年，5月20日，15點11分，星期三83第7章圖7-21 離散相似法仿真非線性系統(tǒng)程序框圖第84頁，共100頁，2022年，5月20日，15點11分，星期三84第7章控制系統(tǒng)的MATLAB仿真求解運行部分如下：Y=zeros(n,1);X=Y;y=0; Uk=zeros(n,1);Ub

59、b=Uk;t=T0;N=round(Tf-T0)/(h*L1);for k=1:Nfor l=1:L1Ub=Uk; Uk=W*Y+W0*Y0;for i=1:nif(Z(i)=0)if (Z(i)=1)Uk(i)=satu(Uk(i),S(i);end 第85頁，共100頁，2022年，5月20日，15點11分，星期三85第7章控制系統(tǒng)的MATLAB仿真if (Z(i)=2) Uk(i)=dead(Uk(i),S(i);endif (Z(i)=3) Uk(i),Ubb(i)=backlash(Ubb(i),Uk(i),Ub(i),S(i);endend end Udot=(Uk-Ub)/h;U

60、f=2*Uk-Ub;X=FI.*X+FIM.*Uk+FIJ.*Udot; Yb=Y;Y=FIC.*X+FID.*Uf;第86頁，共100頁，2022年，5月20日，15點11分，星期三86第7章控制系統(tǒng)的MATLAB仿真for i=1:nif(Z(i)=0)if(Z(i)=4)Y(i)=satu(Y(i),S(i);endif(Z(i)=5)Y(i)=dead(Y(i),S(i);endif (Z(i)=6)Y(i),Ubb(i)=backlash(Ubb(i),Y(i),Yb(i),S(i);第87頁，共100頁，2022年，5月20日，15點11分，星期三87第7章控制系統(tǒng)的MATLAB仿