控制系統(tǒng)數(shù)字仿真實(shí)驗(yàn)指導(dǎo)書(shū)

1、淮 陰 工 學(xué) 院控制系統(tǒng)數(shù)字仿真實(shí)驗(yàn)指導(dǎo)書(shū)編者：王建國(guó)適用學(xué)院：電子與電氣工程學(xué)院電子與電氣工程學(xué)院2015年 6 月 1 日實(shí)驗(yàn)一面向一階微分方程組的連續(xù)系統(tǒng)數(shù)字仿真實(shí)驗(yàn)一、實(shí)驗(yàn)?zāi)康?掌握以一階微分方程組形式描述的連續(xù)系統(tǒng)的數(shù)字仿真方法和步驟。2. 了解如何編程實(shí)現(xiàn)對(duì)一階微分方程組模型的仿真研究。二、實(shí)驗(yàn)內(nèi)容單位反饋控制系統(tǒng)的結(jié)構(gòu)如圖A.1所示。圖A.1 系統(tǒng)結(jié)構(gòu)1根據(jù)系統(tǒng)閉環(huán)傳遞函數(shù)，求出對(duì)應(yīng)的狀態(tài)空間模型，并利用MATLAB中的ODE解函數(shù)編程，對(duì)系統(tǒng)進(jìn)行仿真研究，并繪制對(duì)應(yīng)的單位階躍響應(yīng)曲線(xiàn)。2考慮微分方程模型試編程對(duì)其進(jìn)行仿真（0t100），并用不同的色彩、線(xiàn)型繪制仿真結(jié)果曲線(xiàn)。三

2、、預(yù)習(xí)要求1. 復(fù)習(xí)ODE解函數(shù)使用方法和調(diào)用格式。2復(fù)習(xí)一階微分方程組仿真程序的編制方法。四、實(shí)驗(yàn)報(bào)告要求1. 整理各種實(shí)驗(yàn)條件下的數(shù)據(jù)和曲線(xiàn)。2. 寫(xiě)出實(shí)驗(yàn)報(bào)告。五、教案1利用MATLAB編程求出圖A.1所示的系統(tǒng)狀態(tài)空間模型的程序（文件名為Appl1_1_1.m）如下：% 這是實(shí)驗(yàn)一中第一個(gè)系統(tǒng)求取狀態(tài)空間模型的程序clear% 輸入三個(gè)環(huán)節(jié)的系統(tǒng)傳遞函數(shù)sys1=tf(1,0.5,1);sys2=tf(1,0.2,1);sys3=tf(1,1,0);% 求系統(tǒng)的閉環(huán)傳遞函數(shù)sys=sys1*sys2*sys3;sys_cl_tf=feedback(sys, -1);% 求系統(tǒng)的狀態(tài)空間

3、模型sys_cl_ss=ss(sys_cl_tf)運(yùn)行Appl1_1_1.m后，得到系統(tǒng)的狀態(tài)空間模型為：據(jù)此，編制仿真程序（文件名為Appl1_1_1_fang.m）如下：% 這是實(shí)驗(yàn)一中第一個(gè)系統(tǒng)的仿真程序clearx0=0,0,0;tspan=0,10;t,x=ode45('Appl_1_1_func', tspan, x0);plot(t, 1.25*x(:,3);xlabel('t');ylabel('y(t)');grid;title('實(shí)驗(yàn)一 第一個(gè)系統(tǒng)的單位階躍響應(yīng)');其中，Appl_1_1_func.m程序如下

4、：% 這是實(shí)驗(yàn)一中第一個(gè)系統(tǒng)的對(duì)應(yīng)的數(shù)學(xué)模型function xdot=Appl_1_1_func(t, x)xdot1(1)=-7*x(1)-2.5*x(2)-1.25*x(3)+1;xdot1(2)=4*x(1);xdot1(3)=2*x(2);xdot=xdot1'運(yùn)行Appl1_1_1_fang.m后，得到系統(tǒng)的單位階躍響應(yīng)曲線(xiàn)如圖A.1_1所示。圖A.1_1 實(shí)驗(yàn)一第一個(gè)系統(tǒng)的單位階躍響應(yīng)2. 利用MATLAB編程編制求解系統(tǒng)的仿真程序（文件名為Appl1_1_2_fang.m）如下：% 這是實(shí)驗(yàn)一中第二個(gè)系統(tǒng)的仿真程序clearx0=10,10;tspan=0,100;t,

5、x=ode45('Appl_1_2_func', tspan, x0);plot(t, x(:,1), 'k', 'LineWidth',2);holdplot(t, x(:,2), 'r-.', 'LineWidth',2);legend('x(1)', 'x(2)');xlabel('t');ylabel('狀態(tài)向量');grid;title('實(shí)驗(yàn)一 第二個(gè)系統(tǒng)的仿真輸出');其中，Appl_1_2_func.m程序如下：% 這是

6、實(shí)驗(yàn)一中第二個(gè)系統(tǒng)的對(duì)應(yīng)的數(shù)學(xué)模型function xdot=Appl_1_2_func(t, x)xdot1(1)=- x(1) *(x(1)*x(1)+x(2)*x(2)+x(2);xdot1(2)=-x(1)-x(2) * (x(1)*x(1)+x(2)*x(2);xdot=xdot1'運(yùn)行Appl1_1_2_fang.m后，得到系統(tǒng)的仿真曲線(xiàn)如圖A.1_2所示。圖A.1_2 實(shí)驗(yàn)一第二個(gè)系統(tǒng)的仿真曲線(xiàn)實(shí)驗(yàn)2 面向系統(tǒng)結(jié)構(gòu)圖的連續(xù)系統(tǒng)數(shù)字仿真實(shí)驗(yàn) 一、實(shí)驗(yàn)?zāi)康模?）掌握以系統(tǒng)結(jié)構(gòu)圖形式描述的連續(xù)系統(tǒng)的數(shù)字仿真方法和步驟。（2）初步了解如何用仿真方法來(lái)分析系統(tǒng)的動(dòng)態(tài)性能。（3）了解

7、不同的數(shù)值積分算法與仿真計(jì)算精度之間的關(guān)系。（4）學(xué)會(huì)一種初步尋求合理仿真步長(zhǎng)的方法。（5）了解RK4法計(jì)算穩(wěn)定性和步長(zhǎng)的關(guān)系。二、實(shí)驗(yàn)內(nèi)容單位反饋控制系統(tǒng)的結(jié)構(gòu)如附錄圖1所示，其開(kāi)環(huán)傳遞函數(shù)和閉環(huán)傳遞函數(shù)分別為在零初始條件下閉環(huán)系統(tǒng)單位階躍響應(yīng)的標(biāo)準(zhǔn)解為1 按實(shí)驗(yàn)?zāi)康摹⒁蠛鸵阎獥l件，建立系統(tǒng)的Simulink模型。2. 按下面的經(jīng)驗(yàn)公式（1）初選仿真步長(zhǎng) （為此，需要學(xué)習(xí)MATLAB中求取開(kāi)環(huán)剪切頻率 的函數(shù)margin的使用方法）。 （1） 為系統(tǒng)中反應(yīng)最快的小閉環(huán)的開(kāi)環(huán)剪切頻率。 用margin函數(shù)編寫(xiě)如下程序求開(kāi)環(huán)剪切頻率：clearnum=1;den=0.1,0.7,1,0;G0=

8、tf(num,den);margin(G0)結(jié)果如圖：由此我們得到。根據(jù)經(jīng)驗(yàn)公式可得，仿真步長(zhǎng)的取值范圍為：。3. 選擇RK4法，運(yùn)行仿真模型，適當(dāng)調(diào)整步長(zhǎng)和仿真起止時(shí)間，以得到比較完整的過(guò)渡過(guò)程，觀(guān)察紀(jì)錄過(guò)渡過(guò)程的數(shù)據(jù)。根據(jù)實(shí)驗(yàn)內(nèi)容2結(jié)果得到的仿真步長(zhǎng)的取值范圍：，結(jié)合實(shí)驗(yàn)內(nèi)容1中建立的Simulink仿真模型，設(shè)定仿真參數(shù)時(shí)，我們選取定步長(zhǎng)的ode4（即RK4）算法，仿真時(shí)間0.0到10.0，對(duì)比在仿真步長(zhǎng)在之間變化時(shí)（選取h=0.06,0.10,0.14,0.18,0.22五組數(shù)據(jù)），步長(zhǎng)對(duì)精度的影響。誤差評(píng)價(jià)指標(biāo)公式的程序：t(1)=0, yAE(1)=0;h=0.06(0.10,0.

9、14,0.18,0.22);M=10/h;For k=1:M t(k+1)= t(k)+h;yAE(k+1)=1-*exp(*t(k+1)- *exp(*t(k+1) .*sin( ;endyAE= yAEMAX_ERR=max(abs(yAE-y2);MEAN_ERR=sum(abs(yAE -y2)/size(yAE,1);ERR_NUM=MAX_ERR MEAN_ERR;求解得到的RK4算法步長(zhǎng)與誤差關(guān)系如下：仿真步長(zhǎng)3.7394×10-63.4533×10-51.5475×10-45.2613×10-40.00142.4459×10-7

10、2.1378×10-69.3895×10-62.9336×10-57.4717×10-5上表中數(shù)據(jù)是我們按照RK4法步長(zhǎng)選擇經(jīng)驗(yàn)公式，選取的仿真步長(zhǎng)實(shí)驗(yàn)結(jié)果，為了比較在更小步長(zhǎng)下RK4法精度變化，我們補(bǔ)做在=0.01，0.03,0.05時(shí)仿真，并比較精度，結(jié)果如下：仿真步長(zhǎng)0.010.030.05680×10-52.0467×10-71.7180×10-6130×10-51.3853×10-81.1378×10-7通過(guò)以上兩表比較，得出數(shù)值積分算法在舍入誤差和截?cái)嗾`差二者的影響下，仿真總誤差與步

11、長(zhǎng)的關(guān)系并不是單調(diào)函數(shù)關(guān)系，而是一個(gè)具有極值的函數(shù)關(guān)系。同時(shí)可知，根據(jù)RK4法步長(zhǎng)選取的經(jīng)驗(yàn)公式選擇步長(zhǎng)，能讓我們較為準(zhǔn)確的找到最佳步長(zhǎng)。4改變I環(huán)節(jié)的放大倍數(shù)（將1變?yōu)?），不斷加大步長(zhǎng)觀(guān)測(cè)記錄計(jì)算穩(wěn)定性的變化。(1)h=0.01(2)h=0.1(3)h=0.25(4)h=0.4(5)h=0.5 經(jīng)過(guò)以上五張圖的對(duì)比得知，當(dāng)仿真步長(zhǎng)過(guò)大時(shí)，數(shù)值積分算法的解會(huì)發(fā)生錯(cuò)誤。因此針對(duì)數(shù)值積分算法穩(wěn)定性的要求，必須合理選擇步長(zhǎng)。5. 在相同的條件下，選擇歐拉法和RK2法，再讓仿真模型運(yùn)行，觀(guān)察紀(jì)錄過(guò)渡過(guò)程的數(shù)據(jù)。我們選取=0.03,0.13,0.23三個(gè)步長(zhǎng)，分別利用歐拉法、RK2法、RK4法進(jìn)行仿真

12、，比較在相同步長(zhǎng)情況下，各種方法的精度評(píng)價(jià)指標(biāo)MAX_ERR、MEAN_ERR，實(shí)驗(yàn)結(jié)果如下：評(píng)價(jià)指標(biāo)算法h=0.03h=0.13h=0.23歐拉法0.01540.07400.1473RK2法2.0535×10-40.00400.0138RK4法2.0467×10-71.0773×10-40.0017歐拉法0.00200.00970.0197RK2法2.7789×10-55.4921×10-40.0017RK4法6.9507×10-93.4069×10-64.6734×10-5由上表可以看出，在相同步長(zhǎng)下仿真，RK

13、4法精度明顯高于RK2法和歐拉法。三實(shí)驗(yàn)總結(jié) 本實(shí)驗(yàn)主要研究的是數(shù)值積分中RK4法的步長(zhǎng)選取、精度問(wèn)題和穩(wěn)定性。從實(shí)驗(yàn)結(jié)果中我們得出，步長(zhǎng)的選取對(duì)數(shù)值積分仿真的精度有著很大的影響。我們一般采用經(jīng)驗(yàn)公式或選取仿真步長(zhǎng)，以保證小誤差仿真。如果步長(zhǎng)選擇不合理，會(huì)造成誤差的積累越來(lái)越大，使計(jì)算出現(xiàn)不穩(wěn)定情況，得出錯(cuò)誤的結(jié)果。實(shí)驗(yàn)3 連續(xù)系統(tǒng)離散相似法的數(shù)字仿真實(shí)驗(yàn)一、實(shí)驗(yàn)?zāi)康模?）掌握以系統(tǒng)結(jié)構(gòu)圖形式描述的連續(xù)系時(shí)域離散相似法的數(shù)字仿真方法和步驟。（2）學(xué)會(huì)利用時(shí)域離散相似法分析線(xiàn)性和非線(xiàn)性控制系統(tǒng)的動(dòng)態(tài)性能以及典型非線(xiàn)性環(huán)節(jié)對(duì)控制系統(tǒng)動(dòng)態(tài)性能的影響。二、實(shí)驗(yàn)內(nèi)容含死區(qū)環(huán)節(jié)的非線(xiàn)性控制系統(tǒng)的結(jié)構(gòu)圖如課

14、本所示。（1）具體見(jiàn)課本；（2）求出圖A.2中傳遞函數(shù)對(duì)應(yīng)的狀態(tài)空間模型，并在該模型前加入虛擬的采樣開(kāi)關(guān)和零階保持器，得到離散化狀態(tài)空間模型。程序如下：clearsys1=tf(0.8 1,1 0)；sys2=tf(10,0.1 1)；sys3=tf(1,1 0)；sys= sys1* sys2*sys3;sysc=ss (sys);T=0.05sysd=c2d(sysc,T);（3）在c=0，0.1，0.5，1.0情況下，利用時(shí)域離散相似法編程完成對(duì)該系統(tǒng)的仿真研究。Ad=sysd.a;Bd=sysd.b;Cd=sysd.c;Dd=sysd.d;X=0;0;0;yt=0;tt=0;R=1;C

15、=0(0.1，0.5，1.0)M=30/T;for k=1:M yt1=0; if yt(end)>C; yt1=yt(end)-C; elseif yt(end)<-C; yt1=yt(end)+C; end U=R-yt1; X=Ad*X+Bd*U; Y=Cd*X+Dd*U; yt=yt,Y; tt=tt,k*T;endplot(tt,yt,'k'),grid三、預(yù)習(xí)要求 具體見(jiàn)課本。四、實(shí)驗(yàn)報(bào)告要求 具體見(jiàn)課本。實(shí)驗(yàn)四連續(xù)系統(tǒng)的快速數(shù)字仿真實(shí)驗(yàn)一、實(shí)驗(yàn)?zāi)康?熟悉快速數(shù)字仿真中雙線(xiàn)性變換法、增廣矩陣法和時(shí)域矩陣法的原理及步驟。2. 學(xué)會(huì)編程實(shí)現(xiàn)雙線(xiàn)性變換法、增廣

16、矩陣法和時(shí)域矩陣法，對(duì)控制系統(tǒng)進(jìn)行仿真研究。二、實(shí)驗(yàn)內(nèi)容1已知系統(tǒng)狀態(tài)方程及初值條件 系統(tǒng)響應(yīng)的解析解為試分別采用增廣矩陣法和RK4法對(duì)系統(tǒng)進(jìn)行仿真，并比較兩種方法所得數(shù)值解的計(jì)算精度、所耗費(fèi)的計(jì)算時(shí)間以及計(jì)算穩(wěn)定性。（1）編寫(xiě)程序如下：clearA=-1000.25 999.75 0.5; 999.75 -1000.25 0.5; 0 0 0; t(1)=0; x1(1)=1; % 置解析解x1的初值x2(1)=-1;hh=0.05; M=round(10/hh); % 求解析解for n=1:M tt=t(n)+hh; t(n+1)=tt; x1(n+1)=-exp(-0.5*tt)+ex

17、p(-2000*tt)+1; x2(n+1)=-exp(-0.5*tt)-exp(-2000*tt)+1;end% 采用增廣矩陣法求解T=0.003125; m=4; t_frog=zeros(1,M+5); x=zeros(3, M+5); t_frog(1)=0; x(:,1)=1; -1 ; 1; % 置狀態(tài)向量的初始值star_time_frog=clock; F=expm(A*T); for k=0:m % 完成蛙跳式計(jì)算t_frog(k+2)= 2(k)*T;x(:, k+2)=F*x(:,1);if k<mF=F*F; end end qT=2(m)*T; for k=1:

18、M-1 t_frog(k+6)= t_frog(k+5)+ qT ;x(:, k+6)=F* x(:, k+5);endpass_time_frog=etime(clock, star_time_frog) % 采用RK4法求解h=0.001; % 置RK4法步長(zhǎng)h=0.001st_rk4=zeros(1, 50*M+1); y=zeros(3, 50*M+1); t_rk4(1)=0; y(:,1)=1; -1 ; 1; % 置狀態(tài)向量的初始值star_time_rk4=clock; for k=1:50*Mt_rk4(k+1)=k*h; % 自變量時(shí)間K1=A*y(:,k);K2=A*(y

19、(:,k)+h* K1/2);K3=A*(y(:,k)+h* K2/2);K4=A*(y(:,k)+h* K3);y(:, k+1)=y(:, k)+h*( K1+2* K2+2* K3+K4)/6; endpass_time_rk4=etime(clock, star_time_rk4) xlswrite('file.xlsx',x')（2）計(jì)算精度比較：下面三張圖分別從excel表格中截取了解析解、RK4法和增廣矩陣法最后一頁(yè)的值進(jìn)行比較。三者均是不斷增大。由上圖，增廣矩陣法的計(jì)算精度比較高。（3）耗費(fèi)計(jì)算時(shí)間比較>> y4_1pass_time_fro

20、g = 0.0010pass_time_rk4 =0.0870由以上的程序運(yùn)行結(jié)果，增廣矩陣法速度明顯快于RK4法。（4）計(jì)算穩(wěn)定性增廣矩陣法是一種絕對(duì)穩(wěn)定的算法，允許選擇大采樣周期T。 RK4法是條件穩(wěn)定算法，步長(zhǎng)h的大小除了與算法本身有關(guān)外，還與方程本身的性質(zhì)有關(guān)。對(duì)于實(shí)際系統(tǒng)，由于其特征值不一定為實(shí)數(shù)。一般而言，步長(zhǎng)h必須滿(mǎn)足一定的條件，對(duì)于RK4法，步長(zhǎng)的選擇應(yīng)滿(mǎn)足：2已知控制系統(tǒng)的閉環(huán)傳遞函數(shù)為試采用雙線(xiàn)性變換法（）對(duì)系統(tǒng)進(jìn)行仿真，并與RK4法（）所得結(jié)果（作為標(biāo)準(zhǔn)解）進(jìn)行計(jì)算精度比較。編寫(xiě)程序如下：clear;T=0.5;n1=8;d1=1 4 8;sys1=tf(n1,d1);n

21、2=8;d2=1 4 8;sys2=tf(n2,d2);n3=4;d3=1 4;sys3=tf(n3,d3);sysd=sys1*sys2*sys3;G=c2d(sysd,T,'tustin');A=G.den1;L=length(A);B=G.num1;B=B/A(1);A=A(2:L)/A(1);R=zeros(L,1);Y=zeros(L-1,1);M=10/T;yt=0; tt=0;for k=1:M R(k)=1; R=R(k);R(1:L-1); yk=-A*Y+B*R; Y=yk;Y(1:L-2); yt=yt,yk; tt=tt,k*T; endplot(tt,

22、yt,'k*');hold on;sim('y4');plot(tout,y4,'k')ttRK4法增廣矩陣法00.0239843453605750.1200000000000000.5000000000000000.2641402369901490.41599999999999910.6736723136533440.7855999999999971.5000000000000000.9593505296240661.03471999999999421.044923198775741.0895999999999932.5000000000000

23、001.031322876999861.04129279999999331.007222706515090.9969894399999943.5000000000000000.9978893524483020.98741478399999440.9977727275212970.9946405887999944.5000000000000000.9993374528614051.00062679039999451.000068663152341.0016094904319945.5000000000000001.000137623044091.00058249256959461.0000494

24、84790660.9998793423257546.5000000000000000.9999998674351800.99980972487474670.9999922312849910.9999441673912267.5000000000000000.9999966842249181.00001928305245680.9999997736062331.0000210749898048.5000000000000001.000000407983231.00000473927007091.000000206933870.9999973080145069.5000000000000001.0

25、00000025992390.999997793096632根據(jù)上表，將增廣矩陣法的結(jié)果與RK4法進(jìn)行比較，發(fā)現(xiàn)t越大，計(jì)算的精度越高實(shí)驗(yàn)五采樣控制系統(tǒng)的數(shù)字仿真實(shí)驗(yàn)一、實(shí)驗(yàn)?zāi)康?掌握采樣控制系統(tǒng)數(shù)字仿真的特點(diǎn)。2. 了解數(shù)字控制器對(duì)系統(tǒng)動(dòng)態(tài)性能的影響。3學(xué)會(huì)編制雙重循環(huán)法的仿真程序。二、實(shí)驗(yàn)內(nèi)容某單位反饋控制系統(tǒng)中，被控對(duì)象的傳遞函數(shù)為在單位速度信號(hào)輸入下，選擇零階保持器在時(shí)設(shè)計(jì)的“最少拍”控制器為1. 按實(shí)驗(yàn)?zāi)康?、要求和已知條件，建立系統(tǒng)的Simulink模型，并且編制雙重循環(huán)法的仿真程序。 先用如下程序?qū)⑦B續(xù)系統(tǒng)轉(zhuǎn)換為離散系統(tǒng)c2d函數(shù)，將系統(tǒng)離散化。clear;h=0.01; A=0

26、0;1 -1;B=10;0;G,H=c2d(A,B,h);運(yùn)行程序后得到故連續(xù)系統(tǒng)被控對(duì)象的等價(jià)離散化狀態(tài)方程為： 而題中控制系統(tǒng)的差分方程為： 上式中 雙重循環(huán)法仿真程序：clearG=1 0;0.01 0.99;n=2;m=3;H=0.1;0.0005;c=0 1; h=0.01;T=1; N=round(T/h); TF=15; M=round(TF/T); x=0;0;u=zeros(n,1);e=zeros(m,1);r=zeros(M,1);A=0.282 0.718;B=0.543 -0.471 0.0999;t=0;xt=x;for a=1:M for b=1:N x=G*x+

27、H*u(1); xt=xt,x; t=t,b*h+(a-1)*T; end y=c*x; r(a)=a; ea=r(a)-y; e=ea;e(1:m-1); ua=A*u+B*e; u=ua;u(1:n-1);endx1t=xt(1,:);x2t=xt(2,:);plot(t,x1t,':k',t,x2t,'k');legend('x1', 'x2');grid;xlabel('tame(s)');2.分別運(yùn)行雙重循環(huán)法的仿真程序和Simulink模型，觀(guān)察紀(jì)錄系統(tǒng)動(dòng)態(tài)性能的變化。（1）雙重循環(huán)法仿真程序（2）simulink模型 從雙重循環(huán)法和Simulink模型仿真的結(jié)果可以看出二者基本相同。系統(tǒng)的輸出開(kāi)始會(huì)有振蕩，隨時(shí)間的增加，振幅減小，最后趨于穩(wěn)定。由此我們可以看出，雙重循環(huán)法這種按連續(xù)系統(tǒng)離散相似算法建立連續(xù)部分各個(gè)環(huán)節(jié)仿真模型的方法，不僅仿真結(jié)果準(zhǔn)確，而且在得到系統(tǒng)輸出值的同時(shí)，也得到了被控對(duì)象內(nèi)部狀態(tài)變量的相應(yīng)，實(shí)用性強(qiáng)，值得采用。3若數(shù)字控制器改為“無(wú)紋波”控制器試對(duì)兩種情況下的輸出進(jìn)行比較。設(shè)定仿真時(shí)間為15S ，分別得出結(jié)果：由此可以看出，仿真中，在仿真快結(jié)束時(shí)系統(tǒng)會(huì)出現(xiàn)不穩(wěn)定