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

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

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

3、2one(n)、ones(m,n) 3zeros(m,n)實驗結(jié)果：4rand(m,n) 5diag(v)分析：rand(m,n):產(chǎn)生m*n的隨機陣diag（v）：以v的元素作為對角陣的元素產(chǎn)生對角陣實驗結(jié)果：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的逆矩陣實驗結(jié)果：7roots(p)8Poly分析：roots(p):多項式p的根poly(C):用C的元素為根構(gòu)造多項式實驗結(jié)果：9conv 、deconv分析：conv

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

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

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

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

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

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

10、,j)=1/(i+j-1); endendformat rata實驗結(jié)果：3 m=3;n=4;for i=1:m for j=1:n a(i,j)=1/(i+j-1);endenda實驗結(jié)果：程序2和程序3的區(qū)別：format rat 小數(shù)分?jǐn)?shù)表示。4 x=input(請輸入x的值:); if x=10 y=cos(x+1)+sqrt(x*x+1); else y=x*sqrt(x+sqrt(x); end y實驗結(jié)果：5去掉多項式或數(shù)列開頭的零項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實驗結(jié)果：6 建立MATLAB的函數(shù)文件，程序代碼如下，以文件名ex2_4.m存盤function f=ffibno(n)%ffibno 計算斐波那契亞數(shù)列的函數(shù)文件%n可取任意自然數(shù)%程序如下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)，得到運行結(jié)果。在MATLAB的命令窗口輸入lookfor ffibno，得到結(jié)果：ex2_4.m: %ffibno 計算斐波那契亞數(shù)列的函數(shù)文件在MATLAB的命令窗口輸入help ex2_4，得到結(jié)果： ffibno 計算斐

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

13、了有“It is not one prime”的結(jié)論外，還應(yīng)有：“6=2*3”的說明。實驗結(jié)果：3.2 實驗擴展3.2.1 選擇結(jié)構(gòu)舉例：g=input(please input grade:,s); %由用戶交互輸入成績等級switch(g) %判斷用戶輸入caseA,adisp(85100);caseB,bdisp(7084);caseC,cdisp(6069);caseD,ddisp(=period./2);%看rem函數(shù)功能lsim(num,den,u,t);實驗結(jié)果： 7. 已知開環(huán)系統(tǒng)傳遞函數(shù)，繪制系統(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)實驗結(jié)果： 8. 作如下系統(tǒng)的bode圖 %程序如下： n=1 , 1 ; d=1 , 4 , 11 , 7 ; bode(n , d)實驗結(jié)果：9. 系統(tǒng)傳函如下 求有理傳函的頻率響應(yīng)，然后在同一張圖上繪出以四階伯德近似表示的系統(tǒng)頻率響應(yī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);實驗結(jié)果：10. 已知系統(tǒng)模型為 求它的幅值裕度和相角裕度 %程序如下： n=3.5; d=1 2 3 2; Gm,Pm,Wcg,Wcp=margin(n,d)實驗結(jié)

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) ;實驗結(jié)果：12. 已知系統(tǒng)的開環(huán)傳遞函數(shù)為 繪制系統(tǒng)的Nyqusit圖，并討論系統(tǒng)的穩(wěn)定性 %程序如下： G=tf(1000,conv(1,3,2,1,5);nyquist(G);axis(square)實驗

18、結(jié)果：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) ; %人工變量實驗結(jié)果：14. 一多環(huán)系統(tǒng)，其結(jié)構(gòu)圖如下，使用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);實驗結(jié)果：15. 已知系統(tǒng)為：作該系統(tǒng)的nichols曲線。 %程序如下： n=1 ; d=1 , 1 , 0 ; ngrid(new) ; nichols(n , d) ;實驗結(jié)果：16. 已知系統(tǒng)的開環(huán)傳遞函數(shù)為：當(dāng)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; %波特圖實驗結(jié)果：17. 系統(tǒng)的開環(huán)傳遞函數(shù)為： 分別確定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)實驗結(jié)果：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)實驗結(jié)果：5.2 實驗心得通過本次試驗，熟悉了解了MATLAB中SIMULINK的應(yīng)用，特別是在自動控制系統(tǒng)建模中，它非常的簡便易操作，相比輸入語句的建模來說，它形象而快速，不僅便于理解和分析，而且大大提高了仿真速度。同時，也鞏固了MATLAB的畫圖功能，同樣為之后對控制系統(tǒng)的仿真打下基礎(chǔ)，對于相關(guān)傳遞函數(shù)圖形的繪制有很大的幫助。實驗六 連續(xù)系統(tǒng)數(shù)字仿真的基本算法6.1 實驗題目6.1.1 實驗任務(wù)1理解歐拉法和龍格-庫塔法的基本思想；2理解數(shù)值積分算法的計算精度、速度、穩(wěn)定性與步長

22、的關(guān)系。6.1.2 程序練習(xí)1. 取h=0.2，試分別用歐拉法、RK2法和RK4法求解微分方程的數(shù)值解，并比較計算精度。 注：解析解： %程序如下cleart(1)=0; % 置自變量初值y(1)=1; y_euler(1)=1; y_rk2(1)=1; y_rk4(1)=1; % 置解析解和數(shù)值解的初值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% 輸出結(jié)果disp( 時間 解析解 歐拉法 RK2法 RK4法)yt=t, y, y_euler, y_rk2, y_rk4;disp(yt)實驗結(jié)果：2. 考慮如下二階系統(tǒng)： 在上的數(shù)字仿真解（已知：，），并將不同步長下的仿真結(jié)果與解析解進行精度比較。 說明：已知該微分方程的解析解分別為： 采用RK4法進行計算，選擇狀態(tài)變量： 則有如下狀態(tài)空間模型及初值條件 采用RK4法進行計算。程序如下：clearh=input(請輸入步長h=); %

25、 輸入步長M=round(10/h); % 置總計算步數(shù)t(1)=0; % 置自變量初值y_0(1)=100; y_05(1)=100; % 置解析解的初始值(y_0和y_05分別對應(yīng)于為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); % 置數(shù)值解的初值 % 求解析解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% 求出誤差最