系統(tǒng)仿真技術實驗報告

1、系統(tǒng)仿真技術實驗報告學 院：信息科學與工程學院專 業(yè)：電氣工程及其自動化班 級：XXXXX 姓 名：XXXXX 學 號：XXXXXXXXXXX 目 錄一、 實驗一（MATLAB中矩陣與多項式的基本運算）.41.1 實驗題目.41.2 實驗擴展.131.3 實驗心得.20二、 實驗二（MATLAB繪圖命令）.212.1 實驗題目.212.2 實驗擴展.262.3 實驗心得.29三、 實驗三（MATLAB程序設計）.303.1 實驗題目.303.2 實驗擴展.373.3 實驗心得.38四、 實驗四（MATLAB的符號計算與SIMULINK的使用）.394.1 實驗題目.394.2 實驗擴展.484

2、.3 實驗心得.49五、實驗五（MATLAB在控制系統(tǒng)分析中的應用）.505.1 實驗題目.505.2 實驗心得.71六、 實驗六（連續(xù)系統(tǒng)數字仿真的基本算法）.726.1實驗題目.72七、總結.78 實驗一：MATLAB中矩陣與多項式的基本運算1.1 實驗題目1.1.1 實驗任務： 1了解MATLAB命令窗口和程序文件的調用。2熟悉如下MATLAB的基本運算： 矩陣的產生、數據的輸入、相關元素的顯示； 矩陣的加法、乘法、左除、右除； 特殊矩陣：單位矩陣、“1”矩陣、“0”矩陣、對角陣、隨機矩陣的產生和運算； 多項式的運算：多項式求根、多項式之間的乘除。1.1.2 基本命令訓練：1eye(m)

3、2one(n)、ones(m,n) 3zeros(m,n)實驗結果：4rand(m,n) 5diag(v)分析：rand(m,n):產生m*n的隨機陣diag（v）：以v的元素作為對角陣的元素產生對角陣實驗結果：6AB 、A/B、 inv(A)*B 、B*inv(A)分析：X=AB是AX=B的解，且AB=A(-1)BX=A/B是XB=A的解，且A/B=AB(-1)Inv(A)*B:A的逆矩陣與左乘BB*inv(A)：B左乘A的逆矩陣實驗結果：7roots(p)8Poly分析：roots(p):多項式p的根poly(C):用C的元素為根構造多項式實驗結果：9conv 、deconv分析：conv

4、:進行多項式的乘法deconv:進行多項式的除法實驗結果：10A*B 與 A.*B的區(qū)別分析：A*B：進行矩陣運算A.*B：進行數組運算實驗結果：11who與whos的使用分析：who:列出當前工作空間中所有變量名字whos:列出當前工作空間中所有變量，以及它們的名字、尺寸、所占字節(jié)數、屬性等信息。實驗結果：12disp、size(a)、length(a)的使用分析：disp:屏幕輸出函數size(a):求矩陣的大小(獲取數組的行數和列數)length(a):數組長度(即行數或列數中較大者)實驗結果：12 實驗擴展1.2.1 向量及其運算1）、向量產生實驗結果：2）、向量基本運算分析：dot:

5、維數相同的兩個向量的點積，結果是一個標量cross:兩個向量的叉積就是一個過兩個相交向量的交點且垂直于兩個向量所在平面的向量，結果是一個矢量實驗結果：1.2.2 矩陣及其運算1）、矩陣產生通過冒號可以截取指定矩陣中的部分2）、矩陣基本運算det(A):返回矩陣A的行列式的值eig(A):返回矩陣A的特征值inv(A):求A的逆矩陣fliplr(A):矩陣A關于中軸的反轉矩陣(左右反轉)flipud(A):上下反轉rot90(A):矩陣A逆時針旋轉90度diag(A):抽取主對角線元素構成向量rem(A,x):A矩陣除以模數x后的余數1.2.3 數組及其運算實驗結果：1.2.4 多項式及其運算分

6、析：polyval(P,X):計算多項式P在X點的值（若X是向量或矩陣，則采用數組運算（點運算）polyvalm(P,X):以方陣X為自變量，計算多項式的值，采用矩陣運算1.3、 實驗心得MATLAB中數值運算的部分較為繁多，但和之前學過的C語言等也有相同的部分。但MATLAB的語言更為精簡，比較方便。MATLAB是一款工程性質十分強的軟件，只有大量的應用實踐，才能夠對其中的函數運用自如，才能理解其中的內涵。實驗二 MATLAB繪圖命令2.1 實驗題目21.1實驗任務熟悉MATLAB基本繪圖命令，掌握如下繪圖方法：1坐標系的選擇、圖形的繪制；2圖形注解（題目、標號、說明、分格線）的加入；3圖形

7、線型、符號、顏色的選取。2.1.2 基本命令訓練1plot 2loglog 3semilogx 4semilogy5polar 6title 7xlabel 8ylabel9text 10grid 11bar 12stairs13contour實驗結果：2.2 實驗擴展 2.2.1 圖形的修飾 1. 曲線顏色: b(blue),c(cyan),g(green),k(black),m(magenta),r(red),w(white),y(yellow)2.曲線線型：-（實線，默認）， ：（點線）， -.(點畫線), -(虛線) 3.數據點形：.(實點), +(十字形), 。(圓圈) , *(星號

8、), x(叉號), s(正方形), d(棱形), h(六角形), p(五角形), (右三角), (左三角) 實驗舉例：2.2.2 三維圖形的繪制仿真舉例：t=0:pi/30:10*pi; %設定時間范圍plot3(cos(t),sin(t),t,-b,Linewidth,4)%繪制三維圖形，并作修飾grid on %加網格axis square %坐標為方形figure(2) %新建圖形窗口2stem3(cos(t),sin(t),t,-.k) %繪制三維火柴桿圖實驗結果：2.3 實驗心得用戶應用MATLAB，只需指定繪圖方式，并提供充足的繪圖數據，可以得到所需的圖形，也可對圖形加以修飾。二維

9、圖形的繪制基本函數為plot。較為復雜的三維圖形分為三維曲線，三維網線圖和三維曲面圖。對符號函數或是無法求出函數參數之間顯示關系的隱函數，可以使用符號函數繪制圖形的函數。實驗三 MATLAB程序設計3.1 實驗題目3.1.1實驗任務1熟悉MATLAB程序設計的方法和思路；2掌握循環(huán)、分支語句的編寫，學會使用look for、help命令。31.2 程序訓練1計算11000之內的斐波那契亞數列 f=1,1;i=1;while f(i)+f(i+1)1000 f(i+2)=f(i)+f(i+1); i=i+1;endf,i實驗結果：2 m=3;n=4;for i=1:m for j=1:n a(i

10、,j)=1/(i+j-1); endendformat rata實驗結果：3 m=3;n=4;for i=1:m for j=1:n a(i,j)=1/(i+j-1);endenda實驗結果：程序2和程序3的區(qū)別：format rat 小數分數表示。4 x=input(請輸入x的值:); if x=10 y=cos(x+1)+sqrt(x*x+1); else y=x*sqrt(x+sqrt(x); end y實驗結果：5去掉多項式或數列開頭的零項p=0 0 0 1 3 0 2 0 0 9;for i=1:length(p),if p(1)=0,p=p(2:length(p); end;end

11、;p實驗結果：6 建立MATLAB的函數文件，程序代碼如下，以文件名ex2_4.m存盤function f=ffibno(n)%ffibno 計算斐波那契亞數列的函數文件%n可取任意自然數%程序如下f=1,1;i=1;while f(i)+f(i+1)n f(i+2)=f(i)+f(i+1); i=i+1;end輸入完畢后在MATLAB的命令窗口輸入ex2_4(200)，得到運行結果。在MATLAB的命令窗口輸入lookfor ffibno，得到結果：ex2_4.m: %ffibno 計算斐波那契亞數列的函數文件在MATLAB的命令窗口輸入help ex2_4，得到結果： ffibno 計算斐

12、波那契亞數列的函數文件 n可取任意自然數 程序如下實驗結果：3.1.3 程序設計用一個MATLAB語言編寫一個程序：輸入一個自然數，判斷它是否是素數，如果是，輸出“It is one prime”，如果不是，輸出“It is not one prime.”。要求通過調用子函數實現。最好能具有如下功能：設計較好的人機對話界面，程序中含有提示性的輸入輸出語句。能實現循環(huán)操作，由操作者輸入相關命令來控制是否繼續(xù)進行素數的判斷。如果操作者希望停止這種判斷，則可以退出程序。如果所輸入的自然數是一個合數，除了給出其不是素數的結論外，還應給出至少一種其因數分解形式。例：輸入 6， 因為6不是素數。則程序中除

13、了有“It is not one prime”的結論外，還應有：“6=2*3”的說明。實驗結果：3.2 實驗擴展3.2.1 選擇結構舉例：g=input(please input grade:,s); %由用戶交互輸入成績等級switch(g) %判斷用戶輸入caseA,adisp(85100);caseB,bdisp(7084);caseC,cdisp(6069);caseD,ddisp(=period./2);%看rem函數功能lsim(num,den,u,t);實驗結果： 7. 已知開環(huán)系統(tǒng)傳遞函數，繪制系統(tǒng)的根軌跡，并分析其穩(wěn)定性 %程序如下： num=1 2;den1=1 4 3;d

14、en=conv(den1,den1);figure(1)rlocus(num,den)k,p= rlocfind(num,den) figure(2)k=55;num1=k*1 2;den=1 4 3;den1=conv(den,den);num,den=cloop(num1,den1,-1);impulse(num,den)title(impulse response (k=55) )figure(3)k=56;num1=k*1 2;den=1 4 3;den1=conv(den,den);num,den=cloop(num1,den1,-1);impulse(num,den)title(i

15、mpulse response(k=56)實驗結果： 8. 作如下系統(tǒng)的bode圖 %程序如下： n=1 , 1 ; d=1 , 4 , 11 , 7 ; bode(n , d)實驗結果：9. 系統(tǒng)傳函如下 求有理傳函的頻率響應，然后在同一張圖上繪出以四階伯德近似表示的系統(tǒng)頻率響應 %程序如下： num=1;den=conv(1 2,conv(1 2,1 2); w=logspace(-1,2); t=0.5;m1,p1=bode(num,den,2);p1=p1-t*w*180/pi;n2,d2=pade(t,4);numt=conv(n2,num);dent=(conv(den,d2);m

16、2,p2=bode(numt,dent,w);subplot(2,1,1);semilogx(w,20*log10(m1),w,20*log10(m2),g-);grid on ; title(bode plot);xlabel(frequency);ylabel(gain);subplot(2,1,2);semilogx(w,p1,w,p2,g-);grid on;xlabel(frequency);ylabel(phase);實驗結果：10. 已知系統(tǒng)模型為 求它的幅值裕度和相角裕度 %程序如下： n=3.5; d=1 2 3 2; Gm,Pm,Wcg,Wcp=margin(n,d)實驗結

17、果：11. 二階系統(tǒng)為：令wn=1,分別作出=2 , 1 , 0.707 , 0.5時的nyquist曲線。 %程序如下：n=1 ; d1=1 , 4 , 1 ; d2=1 , 2 , 1 ; d3=1 , 1.414 , 1; d4=1,1,1; nyquist(n,d1) ; hold on nyquist(n,d2) ; nyquist(n,d3) ; nyquist(n,d4) ;實驗結果：12. 已知系統(tǒng)的開環(huán)傳遞函數為 繪制系統(tǒng)的Nyqusit圖，并討論系統(tǒng)的穩(wěn)定性 %程序如下： G=tf(1000,conv(1,3,2,1,5);nyquist(G);axis(square)實驗

18、結果：13. 分別由w的自動變量和人工變量作下列系統(tǒng)的nyquist曲線： %程序如下： n=1 ; d=1 , 1 ,0 ;nyquist(n ,d) ; %自動變量n=1 ; d=1 , 1 ,0 ; w=0.5 : 0.1 : 3 ;nyquist(n , d , w) ; %人工變量實驗結果：14. 一多環(huán)系統(tǒng)，其結構圖如下，使用Nyquist頻率曲線判斷系統(tǒng)的穩(wěn)定性。 %程序如下： k1=16.7/0.0125;z1=0;p1=-1.25 -4 -16;num1,den1=zp2tf(z1,p1,k1);num,den=cloop(num1,den1);z,p,k=tf2zp(num

19、,den);p figure(1)nyquist(num,den)figure(2)num2,den2=cloop(num,den);impulse(num2,den2);實驗結果：15. 已知系統(tǒng)為：作該系統(tǒng)的nichols曲線。 %程序如下： n=1 ; d=1 , 1 , 0 ; ngrid(new) ; nichols(n , d) ;實驗結果：16. 已知系統(tǒng)的開環(huán)傳遞函數為：當k=2時，分別作nichols曲線和波特圖。 %程序如下：num=1;den=conv(conv(1 0,1 1),0.5 1);subplot(1,2,1);nichols(num,den);grid; %

20、 nichols曲線subplot(1,2,2);g=tf(num,den);bode(feedback(g,1,-1);grid; %波特圖實驗結果：17. 系統(tǒng)的開環(huán)傳遞函數為： 分別確定k=2和k=10時閉環(huán)系統(tǒng)的穩(wěn)定性。%程序如下：d1=1 , 3 , 2 , 0 ; n1=2 ;nc1 , dc1=cloop(n1 , d1 ,-1) ;roots(dc1)d2=d1 ; n2=10 ;nc2 , dc2=cloop(n2 , d2,-1) ; roots(dc2)實驗結果：18. 系統(tǒng)的狀態(tài)方程為： 試確定系統(tǒng)的穩(wěn)定性。 %程序如下： a=-4,-3,0 ; 1,0,0 ; 0,1

21、,0 ; b=1;0;0 ; c=0,1,2 ; d=0 ;eig(a) %求特征根rank(ctrb(a,b)實驗結果：5.2 實驗心得通過本次試驗，熟悉了解了MATLAB中SIMULINK的應用，特別是在自動控制系統(tǒng)建模中，它非常的簡便易操作，相比輸入語句的建模來說，它形象而快速，不僅便于理解和分析，而且大大提高了仿真速度。同時，也鞏固了MATLAB的畫圖功能，同樣為之后對控制系統(tǒng)的仿真打下基礎，對于相關傳遞函數圖形的繪制有很大的幫助。實驗六 連續(xù)系統(tǒng)數字仿真的基本算法6.1 實驗題目6.1.1 實驗任務1理解歐拉法和龍格-庫塔法的基本思想；2理解數值積分算法的計算精度、速度、穩(wěn)定性與步長

22、的關系。6.1.2 程序練習1. 取h=0.2，試分別用歐拉法、RK2法和RK4法求解微分方程的數值解，并比較計算精度。 注：解析解： %程序如下cleart(1)=0; % 置自變量初值y(1)=1; y_euler(1)=1; y_rk2(1)=1; y_rk4(1)=1; % 置解析解和數值解的初值h=0.2; % 步長% 求解析解for k=1:5 t(k+1)=t(k)+h; y(k+1)=sqrt(1+2*t(k+1);end% 利用歐拉法求解for k=1:5 y_euler(k+1)=y_euler(k)+h*(y_euler(k)-2*t(k)/y_euler(k);end%

23、 利用RK2法求解for k=1:5 k1=y_rk2(k)-2*t(k)/y_rk2(k); k2=(y_rk2(k)+h*k1)-2*(t(k)+h)/(y_rk2(k)+h*k1); y_rk2(k+1)=y_rk2(k)+h*(k1+k2)/2;end% 利用RK4法求解for k=1:5 k1=y_rk4(k)-2*t(k)/y_rk4(k); k2=(y_rk4(k)+h*k1/2)-2*(t(k)+h/2)/(y_rk4(k)+h*k1/2); k3=(y_rk4(k)+h*k2/2)-2*(t(k)+h/2)/(y_rk4(k)+h*k2/2); k4=(y_rk4(k)+h*

24、k3)-2*(t(k)+h)/(y_rk4(k)+h*k3); y_rk4(k+1)=y_rk4(k)+h*(k1+2*k2+2*k3+k4)/6;end% 輸出結果disp( 時間 解析解 歐拉法 RK2法 RK4法)yt=t, y, y_euler, y_rk2, y_rk4;disp(yt)實驗結果：2. 考慮如下二階系統(tǒng)： 在上的數字仿真解（已知：，），并將不同步長下的仿真結果與解析解進行精度比較。 說明：已知該微分方程的解析解分別為： 采用RK4法進行計算，選擇狀態(tài)變量： 則有如下狀態(tài)空間模型及初值條件 采用RK4法進行計算。程序如下：clearh=input(請輸入步長h=); %

25、 輸入步長M=round(10/h); % 置總計算步數t(1)=0; % 置自變量初值y_0(1)=100; y_05(1)=100; % 置解析解的初始值(y_0和y_05分別對應于為R=0和R=0.5)x1(1)=100; x2(1)=0; % 置狀態(tài)向量初值y_rk4_0(1)=x1(1); y_rk4_05(1)=x1(1); % 置數值解的初值 % 求解析解for k=1:M t(k+1)=t(k)+h; y_0(k+1)=100*cos(t(k+1); y_05(k+1)=100*exp(-t(k+1)/2).*cos(sqrt(3)/2*t(k+1)+100*sqrt(3)/3

26、*exp(-t(k+1)/2).*sin(sqrt(3)/2*t(k+1);end% 利用RK4法求解% R=0for k=1:M k11=x2(k); k12=-x1(k); k21=x2(k)+h*k12/2; k22=-(x1(k)+h*k11/2); k31=x2(k)+h*k22/2; k32=-(x1(k)+h*k21/2); k41=x2(k)+h*k32; k42=-(x1(k)+h*k31); x1(k+1)=x1(k)+h*(k11+2*k21+2*k31+k41)/6; x2(k+1)=x2(k)+h*(k12+2*k22+2*k32+k42)/6; y_rk4_0(k+

27、1)=x1(k+1);end% R=0.5for k=1:M k11=x2(k); k12=-x1(k)-x2(k); k21=x2(k)+h*k12/2; k22=-(x1(k)+h*k11/2)-(x2(k)+h*k12/2); k31=x2(k)+h*k22/2; k32=-(x1(k)+h*k21/2)-(x2(k)+h*k22/2); k41=x2(k)+h*k32; k42=-(x1(k)+h*k31)-(x2(k)+h*k32); x1(k+1)=x1(k)+h*(k11+2*k21+2*k31+k41)/6; x2(k+1)=x2(k)+h*(k12+2*k22+2*k32+k42)/6; y_rk4_05(k+1)=x1(k+1);end% 求出誤差最