數(shù)字仿真報告

1、哈爾濱工程大學(xué)實 驗 報 告 實 驗 名 稱： 控制系統(tǒng)數(shù)字仿真_班 級：_ _學(xué) 號：_ _姓 名：_ _指 導(dǎo) 教 師：_ _實驗室名稱：_ 實驗一：數(shù)值積分法的應(yīng)用 1.實驗?zāi)康氖煜?shù)值積分法在系統(tǒng)仿真中的應(yīng)用，會靈活使用歐拉法，龍格庫塔法等數(shù)值仿真方法解決實際問題，了解步長h的大小對數(shù)值積分法仿真精確度的影響。2.實驗器材電腦一臺，Matlab 7.0軟件3.實驗內(nèi)容分別用歐拉法、二階龍格庫塔及四階龍格庫塔法計算系統(tǒng)：在階躍函數(shù)下的過渡過程。1.選擇相同的步距h0.005，試比較計算結(jié)果。2.選擇不同的步距：歐拉法h0.0001，二階龍格庫塔法h0.005，四階龍格庫塔法h0.05，試

2、比較計算結(jié)果。4.實驗分析 目標(biāo)傳遞函數(shù) ，要求其在階躍函數(shù)下的過渡過程，需要求閉環(huán)傳遞函數(shù)GS(S)，如下： 將其轉(zhuǎn)換為可控標(biāo)準(zhǔn)型的狀態(tài)空間方程為： 針對本系統(tǒng)，需要利用歐拉法，二階龍格庫塔及四階龍格庫塔法的矩陣形式，其中：歐拉法 ，化為矩陣形式如下： 同理：二階-龍格庫塔法矩陣形式為：四階-龍格庫塔法矩陣形式為：5.程序流程圖 經(jīng)過上述分析，可得出程序流程圖：6.實驗結(jié)果選擇相同的步距h0.005，Matlab軟件Command Window內(nèi)的仿真程序為：>>仿真時間Tf=1歐拉法步長h1=0.005二階龍格庫塔法步長h2=0.005四階龍格庫塔法步長h3=0.005 得到的

3、仿真圖形如下： 其中圖例中的matlab為直接仿真結(jié)果輸出，紅色曲線均為利用歐拉法，二階龍格庫塔及四階龍格庫塔法的仿真結(jié)果。由圖可以看出，歐拉法誤差較大，二階龍格庫塔及四階龍格庫塔法仿真結(jié)果接近目標(biāo)函數(shù)實際圖像。選擇不同的步距：歐拉法h0.0001，二階龍格庫塔法h0.005，四階龍格庫塔法h0.05。Matlab軟件Command Window內(nèi)的仿真程序為：>>仿真時間Tf=1歐拉法步長h1=0.0001二階龍格庫塔法步長h2=0.005四階龍格庫塔法步長h3=0.05仿真圖像為：可發(fā)現(xiàn)在不同的步長下，歐拉法的精度很高，二階龍格庫塔比四階龍格庫塔法的仿真效果更好。這可以說明：在

4、步長無限制的情況下，任何仿真方法都可以達(dá)到目標(biāo)精度，但相對應(yīng)犧牲了運行速度，加大了工作量。7.程序clc;clear;clf;num0=(conv(100,5 1);dem0=(conv(10,1,conv(1,1,0.15,1);num1=1;dem1=1;num den=feedback(num0,dem0,num1,dem1);sys1=tf(num,den);%Transfer function:% 500 s + 100%-%1.5 s3 + 11.65 s2 + 511.1 s + 101A,B,C,D=tf2ss(num,den);Tf=input('仿真時間Tf=

5、9;);u=1;%歐拉法h1=input('歐拉法步長h1=');x1=zeros(length(A),1);y1=0;t1=0;for i1=1:Tf/h1 x1=x1+h1*(A*x1+B*u); y1=y1;C*x1;t1=t1;t1(i1)+h1;endsubplot(3,1,1);hold onplot(t1,y1,'r');step(sys1);axis(0,1,0,1.6);title('歐拉法');xlabel('時間（s）');ylabel('輸出量y');legend('歐拉法'

6、,'matlab');hold off;grid%二階龍格庫塔法h2=input('二階龍格庫塔法步長h2=');x2=zeros(length(A),1);y2=0;t2=0;for i2=1:Tf/h2 k1=A*x2+B*u; k2=A*(x2+h2*k1)+B*u; x2=x2+0.5*h2*(k1+k2); y2=y2;C*x2;t2=t2;t2(i2)+h2;endsubplot(3,1,2);hold onplot(t2,y2,'r');step(sys1);axis(0,1,0,1.6);title('二階龍格庫塔法

7、9;);xlabel('時間（s）');ylabel('輸出量y');legend('二龍','matlab');hold off;grid%四階龍格庫塔法h3=input('四階龍格庫塔法步長h3=');x3=zeros(length(A),1);y3=0;t3=0;for i3=1:Tf/h3 k1=A*x3+B*u; k2=A*(x3+h3*k1/2)+B*u; k3=A*(x3+h3*k2/2)+B*u; k4=A*(x3+h3*k3)+B*u; x3=x3+h3*(k1+2*k2+2*k3+k4)/6;

8、y3=y3;C*x3;t3=t3;t3(i3)+h3;endsubplot(3,1,3);hold onplot(t3,y3,'r');step(sys1);axis(0,1,0,1.6);title('四階龍格庫塔法');xlabel('時間（s）');ylabel('輸出量y');legend('四龍','matlab');hold off;gridY1=;Y2=;Y3=;for i=1:10:200Y1=Y1;y1(i);Y2=Y2;y2(i);Y3=Y2;y2(i);end實驗二：控制系統(tǒng)優(yōu)

9、化設(shè)計與仿真1. 實驗?zāi)康氖煜渭冃畏ǖ脑砑盎緦?yōu)過程，了解改進的單純形法基本原理。會使用改進的單純尋優(yōu)法解決實際問題。2.實驗器材電腦一臺，Matlab 7.0軟件3.實驗內(nèi)容利用單純形法尋優(yōu)程序，并在計算機上對下圖所示的系統(tǒng)中的控制器參數(shù)T1,K1進行尋優(yōu)。指標(biāo)函數(shù)取 Q=。初值選擇為，初始單純形的尺寸為0.05。4.實驗分析 本實驗利用了改進單純形法對系統(tǒng)進行了尋優(yōu)。改進單純形法的基本思想是：給定初始點(0)和步長a，產(chǎn)生初始單純形S0，通過反射、擴張、收縮和緊縮等一系列動作將單純形翻滾、變形，從而產(chǎn)生一系列的單純形S1，S2，S3，,逐漸向極小值點靠攏。當(dāng)滿足精度指標(biāo)時，迭代停止，

10、取當(dāng)前單純形的“最好點”作為極小點的近似。改進單純形法的迭代規(guī)則: 假設(shè)當(dāng)前單純形為Sk ,對組成單純形的(m+1)個頂點，記L為“最好點”, H為“最壞點”, G為“次壞點”,即 首先，計算當(dāng)前單純形的(m+1)個頂點中去掉最壞點H后的形心： 判別是否滿足終止條件，即計算 如果error<(為給定的精度指標(biāo)),則停止迭代，取當(dāng)前單純形的“最好點”L作為所極小值點*的近似。否則，計算“最壞點”H關(guān)于形心的反射點R ，其中： 。根據(jù)反射點函數(shù)值與目標(biāo)函數(shù)值相比較之后得出的不同結(jié)論，進行相對應(yīng)的擴張，收縮操作，經(jīng)過多次循環(huán)之后，得到最優(yōu)點。 5.程序流程圖根據(jù)改進單純性法的尋優(yōu)過程，分析反射

11、點函數(shù)值與各目標(biāo)函數(shù)值的大小，進行相對應(yīng)操作，得出程序流程圖如下：6.實驗結(jié)果根據(jù)改進的單純尋優(yōu)法的尋優(yōu)過程，初始值選擇為，選擇擴張因子=2，收縮因子=0.5，對初始值進行擴張和收縮。其中表1給出了尋優(yōu)過程中最好點及最壞點的參數(shù)變化。表1 尋優(yōu)過程中最好點與最壞點的變化尋優(yōu)次數(shù)最壞點k1最壞點T1最好點k1最好點T1備注00.300095.00000.350095.0000最初點10.350095.0500.350095.0000反射點20.450094.97500.350095.0000擴張點30.450094.92500.450094.9750反射點40.500094.92500.4500

12、94.9250收縮點 50.425094.96250.450094.9250收縮點60.456294.93440.450094.9250收縮點 70.437594.94380.450094.9250擴張點80.437594.91880.437594.9438反射點9 0.431394.93440.437594.9188收縮點 100.440694.92970.437594.9188收縮點11 0.440694.90740.440694.9297反射點 120.441494.91050.440694.9047收縮點 130.439194.91170.4406 94.9074擴張點 140.438

13、394.89020.439194.9117擴張點 150.440694.90740.438394.8902最優(yōu)點 最優(yōu)點可計算出誤差s= 8.0781e-006<e=0.00001,滿足尋優(yōu)條件，則可繪出最壞點與最好點k1,T1的變化圖像如下： 同樣，通過運行編寫的matlab程序，可給出尋優(yōu)前后動態(tài)響應(yīng)的比較圖如下：如圖可以看出，采用改進后的單純尋優(yōu)法尋優(yōu)以后，動態(tài)過程的調(diào)節(jié)時間有所減少，但超調(diào)量與振蕩次數(shù)明顯增大，但經(jīng)過一段時間的不穩(wěn)定后，控制過程趨于穩(wěn)定?？梢钥闯鲞@次尋優(yōu)的結(jié)果并不是很令人滿意?？筛淖償U大因子與收縮因子的大小或重新選取k1,T1來改變這種情況。7.程序控制主m文件c

14、lear;clc;%初始化global yy;global kk;N=2;h=0.05;gama=2;beita=0.5;k_max=25;%最大搜索次數(shù)k=0;%當(dāng)前搜索次數(shù)disp('初始點');k1=0.3;T1=95;alpha=zeros(N+1,N+1);alpha(1,1)=k1;alpha(2,1)=T1;e=eye(N,N);alpha_h=zeros(2,1);%最壞點alpha_l=zeros(2,1);%最好點alpha_g=zeros(2,1);%次壞點alpha_r=zeros(2,1);%反射點alpha_e=zeros(2,1);%擴張點alph

15、a_c=zeros(2,1);%壓縮點%2.計算單純形的初始點for i=1:N for j=2:N+1 alpha(i,j)=alpha(i,1)+h*e(i,j-1); endend%3.計算每個初始點對應(yīng)的函數(shù)值for i=1:N+1 alpha(3,i)=Q(alpha(1,i),alpha(2,i);endalphafor i=1:N+1for j=i+1:N+1 if alpha(3,i)<alpha(3,j) l_k1=alpha(1,i); l_T1=alpha(2,i); l_Q=alpha(3,i); alpha(1,i)=alpha(1,j); alpha(2,i)

16、=alpha(2,j); alpha(3,i)=alpha(3,j); alpha(1,j)=l_k1; alpha(2,j)=l_T1; alpha(3,j)=l_Q; endendendfor i=1:N alpha_h(i,1)=alpha(i,1); alpha_g(i,1)=alpha(i,2); alpha_l(i,1)=alpha(i,3);endc_h=alpha(3,1);c_g=alpha(3,2);c_l=alpha(3,3);e=0.00001;aH=;aL=;while k<=k_max if c_h<c_g%重新找出最壞點、最好點和次壞點 l_Q=c_h

17、; l_s=alpha_h; c_h=c_g; alpha_h=alpha_g; c_g=l_Q; alpha_g=l_s; end if c_g<c_l l_Q=c_g; l_s=alpha_g; c_g=c_l; alpha_g=alpha_l; c_l=l_Q; alpha_l=l_s; end %確定形心 alpha_x=(alpha_g+alpha_l)/2; c_x=Q(alpha_x(1,1),alpha_x(2,1); s=sqrt(c_g-c_x).2+(c_l-c_x).2)/3); aH=aH,alpha_h' aL=aL,alpha_l' if s

18、>=e alpha_r=2*alpha_x-alpha_h; %計算反射點 c_r=Q(alpha_r(1,1),alpha_r(2,1); if c_r>=c_l if c_r>c_g if c_r<c_h alpha_c=alpha_x+beita*(alpha_r-alpha_x); c_c=Q(alpha_c(1,1),alpha_c(2,1); else alpha_c=alpha_x+beita*(alpha_r-alpha_x); c_c=Q(alpha_c(1,1),alpha_c(2,1); end if c_c>c_h alpha_g=alph

19、a_l+(alpha_g-alpha_l)/2; alpha_h=alpha_l+(alpha_h-alpha_l)/2; else alpha_h=alpha_c; c_h=c_c; disp('收縮'); alpha_h alpha_l end else alpha_h=alpha_r; c_h=c_r; disp('反射點'); alpha_h alpha_l end else alpha_e=alpha_x+gama*(alpha_r-alpha_x); c_e=Q(alpha_e(1,1),alpha_e(2,1); if c_e<c_l alp

20、ha_h=alpha_e; c_h=c_e; disp('擴張'); alpha_h alpha_l else alpha_h=alpha_r; c_h=c_r; disp('反射點'); alpha_h alpha_l end end k=k+1; else break; endend disp('最優(yōu)點'); alpha_h alpha_l Q1(k1,T1);figure(1);Q(alpha_l(1,1),alpha_l(2,1);plot(kk,yy,'b','LineWidth',2);legend(&

21、#39;優(yōu)化前','優(yōu)化后');hold off;grid on;title('優(yōu)化前后輸出曲線');ah_K=;ah_T=;al_T=;al_K=;for j=1:2:21 ah_K=ah_K,aH(j); ah_T=ah_T,aH(j+1); al_K=al_K,aL(j); al_T=al_T,aL(j+1);endfigure(2);grid on;subplot(2,2,1);plot(ah_K,'-rs','LineWidth',2);title('最壞點K');grid on;subplot(

22、2,2,2);plot(ah_T,'-rs','LineWidth',2);title('最壞點T');grid on;subplot(2,2,3);plot(al_K,'-rs','LineWidth',2);title('最好點K');grid on;subplot(2,2,4);plot(al_T,'-rs','LineWidth',2);title('最好點T');grid on;子m文件-Q函數(shù)function Q=Q(k1,T1)num=10*k1*T1 0;den=2*T1 2+T1 1 0;a b c d=tf2ss(num,den);T=0.005;h=T;N=5000;x=0;0;0;R=1*ones(1,N);yy(1,1)=0;f