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

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

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

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

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

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

7、s)；第7章控制系統的MATLAB仿真第11頁，共100頁，2022年，5月20日，15點11分，星期三11第7章【例7.1】已知某系統的狀態空間表達式為：控制系統的MATLAB仿真試采用MATLAB語言求出該系統的狀態空間模型。第12頁，共100頁，2022年，5月20日，15點11分，星期三12解：采用狀態空間模型表示時，可在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)上述指令執行后可得指定系統的狀態空間模型為：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章控制系統的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. 控制系統的MATLAB仿真第14頁，共100頁，2022年，5月20日，15點11分，星期三14 7.1.4 系統不同模型間的相互轉換 在MATLAB 6.1的信號處理工具箱與控制系統工具箱中，提供了傳遞函數模型、零

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

10、為狀態空間模型zp2tf將系統零極點增益模型轉換為傳遞函數模型第16頁，共100頁，2022年，5月20日，15點11分，星期三16【例7.2】 已知某系統的傳遞函數為：第7章控制系統的MATLAB仿真 試用MATLAB語言求出該系統的傳遞函數模型、狀態空間模型和零極點增益模型。解：（1）求系統的傳遞函數模型在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章執行以上語句后可得系統的傳遞函數模型為：12 s3 + 24 s2 + 12 s +

11、20-2 s4 + 4 s3 + 6 s2 + 2 s + 2（2）求系統的狀態空間模型該系統的狀態空間模型可以通過MATLAB的模型轉換函數來完成。在MATLAB命令窗口輸入以下命令：a,b,c,d=tf2ss(num,den);sys=ss(a,b,c,d) 控制系統的MATLAB仿真第18頁，共100頁，2022年，5月20日，15點11分，星期三18執行完上述語句后，可得系統的狀態空間模型的狀態矩陣a，系統控制矩陣b，系統輸出矩陣c，系統直接傳輸矩陣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章控制系統的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.控制系統的MATLAB仿真第20頁，共100頁，2022年，5月20日，15點11分，星期三20由以上數據可寫出系統的狀態空間模型為：第7章控制系統的MATLAB仿真第21頁，共100頁，2022年，5月20日，15點11分，星期三21第7章（3）求系統的零極點增益模型該系統的零極點增益模型也可以通過MATLA

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

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

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

16、試求串聯連接的等效傳遞函數及其等效狀態空間模型。第25頁，共100頁，2022年，5月20日，15點11分，星期三25解：（1）根據MATLAB程序設計的基本方法和函數命令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章控制系統的MATLAB仿真第26頁，共100頁，2022年，5月20日，15點11分，星期三26第7章控制系統的MATLAB仿真（2）在MATLAB命令窗口輸入程序名，程序運行后得到如下電流環內前向通道的等效傳遞函數及等效狀態空間模型：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章控制系統的MATLAB仿真第28頁，共100頁，2022年，5月20日，15點11分，星期三28第7章 7.2.2 環節并聯連結的化簡 環節并聯是指多個環節的輸入信號相同，所有環節輸出的代數和為其總輸出。 采用parallel()函數命令可以等效化簡兩個環節的并聯連接。parallel()函數命令調用格式為：num，den= parallel（numl,denl,num2,den2）該命令已由命令sys=sysl+sys2+sysn所

19、取代，省掉了“parallel()”字符，且可以實現多個環節的并聯等效處理。parallel()函數命令調用格式還有：a,b,c,d= parallel(a1,bl,cl,dl,a2,b2,c2,d2) 控制系統的MATLAB仿真第29頁，共100頁，2022年，5月20日，15點11分，星期三29第7章【例7.4】 已知兩個環節的傳遞函數分別為： 試求兩環節并聯連接等效傳遞函數的num與den向量及等效的狀態空間模型。控制系統的MATLAB仿真第30頁，共100頁，2022年，5月20日，15點11分，星期三30第7章解：（1）采用函數命令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命令窗口運行程序后得到等效傳遞函數及等效狀態空間模型為：Transfer function: 3 s2 + 6 s + 4-s3 + 3 s2 + 4 s + 4 控制系統的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控制系統的MATLAB仿真第32頁，共100頁，2022年，5月20日，15點11分，星期三32 7.2.3 環節反饋連結的化簡 利用MATLAB中的feedback()函數命令可將兩個環節按反饋形式進行連接后求其等效傳遞函數。feedback()函數命令形式：（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章控制系統的MATLAB仿真第33頁，共100頁，2022年，5月20日，15點11分，星期三33第7章 【例7.5】已知晶閘管-直流電機單閉環調速系統的動態結構圖如圖7-4所示，求該閉環系統的傳遞函數。控制系統的MATLAB仿真圖7-4 直流單閉環調速系統第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)控制系統的MATLAB仿真第35頁，共100頁，2022年，5月20日，15點11分，星期三35第7章 （2）程序運行結果為：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 控制系統的MATLAB仿真第36頁，共100頁，2022年，5月20日，15點11分，星期三36第7章由以上運算數據可得單閉環系統的小閉環的傳遞函數為：單閉環系統的閉環的傳遞函數（略去分母的項）為：控制系統的MATLAB仿真第37頁，共100頁，2022年，5月20日，15點11分，星期三37第7章7.3 控制系統時域分析的MATLAB實現7.3.1 概述1系統時域響應仿真的主要問題時域分析法是根據系統微分方程用拉普拉斯變換求解系統的動態響應。典型的動態過程響應有單位階躍響

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

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

27、ys)； 求出了系統的階躍響應曲線，要分析系統的穩定性能，還需要計算系統的性能指標。包括峰值時間tp、超調量%、調節時間ts、及穩態誤差ess等。這些參數可以直接從圖上讀取，也可通過MABLAB的M函數來求取。 控制系統的MATLAB仿真第40頁，共100頁，2022年，5月20日，15點11分，星期三40第7章下面介紹利用M函數來計算性能指標。定義該指標函數為terg函數，其程序語句如下：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,控制系統的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;控制系統的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為第一峰值的超調量；b2為第二峰值的超調量；sigma為階躍響應的超調量；n為階躍響應的衰減比；pusi為階躍響應的衰減率；T為衰減振蕩周期；f振蕩頻率。控制系統的MATLAB仿真第43頁，共100頁，2022年，5月20日，15點11分，星期三43第7章7.3.3 連續系統的單位脈沖響應求連續系統單位脈沖響應的函數為impulse()，該函數命令的調用格式和功能如下： （1）impulse(s

30、ys)；（2）impulse(sys，t)；（3）impulse(sysl,sys2,，sysN)； （4）y，t，x=impulse(sys)； 【例7.6】已知單位負反饋系統前向通道的傳遞函數為： ，試繪出系統的單位階躍響應和脈沖響應曲線。并計算系統的階躍響應性能指標。控制系統的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);% 求閉環傳遞函數step(closys)可得如圖7-5所示的階躍響應曲

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

32、繼續在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控制系統的MATLAB仿真第48頁，共100頁，2022年，5月20日，15點11分，星期三48第7章7.3.4 連續系統的零輸入響應求連續系統零輸入響應的

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

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

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

36、2求離散系統Bode圖的函數dbode()求離散系統Bode圖的函數為dbode()，該函數命令的調用格式為：（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) 控制系統的MATLAB仿真第53頁，共100頁，2022年，5月20日，15點11分，星期三53第7章 7.4.3 繪制系統Nyquist曲線的函

37、數 計算或繪制系統Nyquist曲線的函數為nyquist()，該函數命令的調用格式為：（1）re,im,w=nyquist(sys)（2）nyquist(sys)（3）nyquist(sys,w)（4）nyquist(sysl,sys2,sysn) （5）nyquist(sysl,sys2,sysn,w)控制系統的MATLAB仿真第54頁，共100頁，2022年，5月20日，15點11分，星期三54第7章 如果系統為離散系統，則相應的nyquist函數為 dnyquist()，函數命令調用為：（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)控制系統的MATLAB仿真第55頁，共100頁，2022年，5月20日，15點11分，星期三55第7章 7.4.4 求連續系統Nichols曲線的函數 求連續系統Nichols曲線的函數為nichols()，該函數命令的調用格式為：（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)控制系統的MATLAB仿真第56頁，共100頁，2022年，5月20日，15點11分，星期三56第7章控制系統的MATLAB仿真 7.4.5 求系統幅值裕度與相位裕度 求系統幅值裕度與相位裕度的函數為margin()，該函數命令的調用格式為：（1）Gm,Pm,Wcp,Wcg= margin(sys)； （2）Gm,Pm,Wcp,Wcg= margin(mag,phase,w)（3）margin(sys)； （4）margin(mag,phase,w)； 對于離散時間系統，可以先用函數dbode

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

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

42、7-12 面向傳遞函數的線性系統仿真框圖第61頁，共100頁，2022年，5月20日，15點11分，星期三61第7章控制系統的MATLAB仿真設源程序名為CSS1.m，保存為M文件。1）輸入數據a=;% n+1維分母系數向量c=;% m+1維分子系數向量 X0=;% 狀態向量初值V=;% 反饋系數n=;% 系統階次T0=;% 起始時間Tf=;% 終止時間h=;% 計算步長R=r;% 階躍輸入函數幅值 第62頁，共100頁，2022年，5月20日，15點11分，星期三62第7章控制系統的MATLAB仿真 2）形成開、閉環系數陣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);% 分子系數向量維數m+1C=fliplr(c),zeros(1,n-m1);% 形成輸出陣CAb=A-B*C*V;% 形成閉環系數陣AbX=X0;y=0;t=T0;% 設初值，準備遞推運算第63頁，共100頁，2022年，5月20日，15點11分，星期三63第7章控制系統的MATLAB仿真 3）運算求解N=round(Tf-T0)/h;% 確定輸出點數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;% 求狀態y=y,C*X;% 求輸出并以向量形式保存 t=t,t(i)+h;% 輸出對應時刻以向量形式保存end第64頁，共100頁，2022年，5月20日，15點11分，星期三64第7章控制系統的MATLAB仿真 4）輸出結果t,y % 輸出數據形式結果plot(t,y)% 輸出曲線形式結果 以上即為采用MATLAB語言實現面向典型閉環系統數字仿真的程序，可編輯為CSS1.m文件存儲起來。使用時，只要進入到MATL

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

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

47、元素n=;% 環節個數（系統階次）Y0=;% 階躍輸入幅值Yt0=% 各環節輸出初值h=;% 計算步長L1=;%每隔L1個點打印一次 T0=;% 起始時間 Tf=;% 終止時間nout=Nout;% 輸出環節編號 第69頁，共100頁，2022年，5月20日，15點11分，星期三69第7章控制系統的MATLAB仿真程序處理開始：A=diag(P(:,1);B=diag(P(:,2);% 求開環系統系數陣C=diag(P(:,3);D=diag(P(:,4);m=length(WIJ(:,1);% 求非零元素個數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章控制系統的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;% 形成閉環系數陣Y=Yt0; y=Y(nout);t=T0;% 置初值N=round(Tf-T0)/(h*L1);% 總輸出點數for i=1:N; % 每循環一次，輸出一點

49、數據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章控制系統的MATLAB仿真y=y,Y(nout);% 保存輸出環節動態響應值t=t,t(i)+h*L1;% 保存時間向量endt,y% 輸出響應數據plot(t,y)% 輸出響應圖形 程序中用到MATLAB特殊功能函數diag(

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

51、框圖第7章第74頁，共100頁，2022年，5月20日，15點11分，星期三74第7章控制系統的MATLAB仿真CSS3仿真程序設計如下：A=diag(P(:,1);B=diag(P(:,2);% 求開環系統系數陣C=diag(P(:,3);D=diag(P(:,4);m=length(WIJ(:,1);% 求非零元素個數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章控制系統的MATLAB仿真for i=1:nif (A(i,i)=0);% 求積分或積分比例環節各系數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;% 求積分環節系數if (D(i,i)=0);% 若為積分比例環節，修正系數FID(i)=D(i,i)/B(i,i);elseendelse% 求慣性或慣性比例環節各系數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章控制系統的MATLAB仿真FIC(i)=1;FID(i)=0;% 求慣性環節系數if (D(i,i)=0);% 若為慣性比例環節，修正系數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、計算總輸出點數 第77頁，共100頁，2022年，5月20日，15點11分，星期三77第7章控制系統的MATLAB仿真for k=1:Nfor l=1:L1Ub=Uk; % 保存前一次輸入值Uk=W*Y+W0*Y0; % 求當前UkUdot=(Uk-Ub)/h; % 求當前Uk的導數Uf=2*Uk-Ub; % 求下一步的輸入UfX=FI.*X+FIM.*Uk+FIJ.*Udot;% 求狀態向量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章控制系統的MATLAB仿真該程序用到了幾個MATLAB語言中的特殊運算符： 冒號符，用于生成行向量；= = 邏輯等號符，MATLAB中的規定寫法； = 邏輯不等號符，MATLAB中的規定寫法； * 點乘運算符，MATLAB中的規定寫法，表示兩個同維矩陣中各相應位置上的元素相乘運算。第79頁，共100頁，2022年，5月20日，15點11分，星期三79第7章控制系統的MATLAB仿真7.5.4 非線性系統離散相似法仿真1典型非線性環節的求解（1）飽和非線性環節：采用MATLAB語言來編制M函數如下，函數名為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章控制系統的MATLAB仿真（2）死區非線性：采用MATLAB語言來編制M函數如下，函數名為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章控制系統的MATLAB仿真（3）滯環非線性：M函數名為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章控制系統的MATLAB仿真非線性系統的離散相似法仿真程序命名為css4.m。根據上述分析可以畫出CSS4程序框圖如7-21所示。CSS4程序設計：輸入數據：輸入系統各環節參數陣P，連接陣WIJ，運行參數n,h

58、,L1,T0,Tf,階躍輸入幅值Y0，輸出環節編號nout(方法同CSS3)，以及非線性標志向量Z，參數向量S。程序開始：求開環各系數陣A,B,C,D 及連接關系矩陣 W,W0 (方法同CSS3)；求取各環節離散系數FI(i),FIM(i),FIJ(i),FIC(i),FID(i)。第83頁，共100頁，2022年，5月20日，15點11分，星期三83第7章圖7-21 離散相似法仿真非線性系統程序框圖第84頁，共100頁，2022年，5月20日，15點11分，星期三84第7章控制系統的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章控制系統的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章控制系統的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章控制系統的MATLAB仿