第3-1章-連續(xù)系統(tǒng)數(shù)值積分法仿真Matlab編程課件

第三章Matlab編程(數(shù)值積分法仿真)3.1連續(xù)系統(tǒng)數(shù)值積分法仿真編程思路目的：針對(duì)下面的系統(tǒng)編寫通用的數(shù)值積分法仿真程序(3.1-1)對(duì)這樣的系統(tǒng)進(jìn)行仿真，實(shí)際上涉及到控制的計(jì)算、狀態(tài)的數(shù)值積分計(jì)算和輸出的計(jì)算。3個(gè)函數(shù)g,f,h確定后，就可以完整地描述一個(gè)系統(tǒng)。給定初值：u0，y0,系統(tǒng)中的向量都采用列向量表示1第三章Matlab編程(數(shù)值積分法仿真)3.1連續(xù)系統(tǒng)數(shù)function[t,y]=w_DigiInteSimu(tstart,tstop,h,x0,u0,cnty,InteMethod,StateModel,OutputFile,ControlFile)%函數(shù)功能：用數(shù)值積分法仿真%輸入?yún)?shù)：tstart,tstop,h分別是起始時(shí)間、結(jié)束時(shí)間和仿真步長(zhǎng)，是標(biāo)量%x0,u0是狀態(tài)、輸入的初值，都是列向量%cnty是輸出變量的個(gè)數(shù)%InteMethod時(shí)數(shù)值積分方法，可以是'EUL','RK2','RK4'%StateModel是狀態(tài)模型的文件名%ControlFile是控制作用的文件名%OutputFile是系統(tǒng)輸出的文件名%輸出參數(shù)：t是仿真結(jié)果的時(shí)間序列%y是仿真結(jié)果系統(tǒng)的輸出序列(1)數(shù)值積分法仿真框架函數(shù)2function[t,y]=w_DigiInteSimu(function[t,y]=w_DigiInteSimu(tstart,tstop,h,x0,u0,cnty,InteMethod,StateModel,OutputFile,ControlFile)t=[tstart:h:tstop];%t數(shù)一個(gè)行序列cntt=size(t,2);%返回列數(shù)y=zeros(cnty,cntt);%構(gòu)造一個(gè)空矩陣，用來(lái)存儲(chǔ)結(jié)果y0=eval([OutputFile,'(tstart,x0,u0)']);%計(jì)算初始輸出y(:,1)=y0’;%將cury作為輸出的第1列curx=x0;%當(dāng)前一步的xcuru=u0;%當(dāng)前一步的ucury=y0;%當(dāng)前一步的yfori=1:1:cntt-1curu=eval([ControlFile,'(t(i),h,curx,curu,cury)']);%計(jì)算控制時(shí)傳遞的參數(shù)：當(dāng)前時(shí)間，步長(zhǎng)，當(dāng)前狀態(tài)和輸出curx=w_StepIntegral(t(i),h,curx,curu,InteMethod,StateModel);%單步積分運(yùn)算cury=eval([OutputFile,'(t(i),curx,curu)']);%計(jì)算輸出y(:,i+1)=cury‘;%將輸出加入到輸出序列里end3function[t,y]=w_DigiInteSimu(functionNewX=w_StepIntegral(curt,h,curx,curu,InteMethod,StateModel)%函數(shù)功能：?jiǎn)尾椒e分運(yùn)算,與模型方程有關(guān)%輸入?yún)?shù)：curt是當(dāng)前時(shí)間，h是數(shù)值積分步長(zhǎng)% curx,curu分別是當(dāng)前的狀態(tài)和控制向量%InteMethod是積分算法，字符串類型，可以取'EUL','RK2','RK4'%StateModel是狀態(tài)模型文件名稱，字符串類型%輸出參數(shù)：NewX是這一步計(jì)算的新的狀態(tài)向量(2)單步數(shù)值積分法函數(shù)單步數(shù)值積分函數(shù)只是對(duì)微分方程組StateModel進(jìn)行一步的計(jì)算，計(jì)算法由InteMethod參數(shù)指定，可以上歐拉法，RK2或RK4。4functionNewX=w_StepIntegral(cfunctionNewX=w_StepIntegral(curt,h,curx,curu,InteMethod,StateModel);ifInteMethod=='RK4'k1=eval([StateModel,'(curt,curx,curu)']);k2=eval([StateModel,'(curt+0.5*h,curx+0.5*h*k1,curu)']);k3=eval([StateModel,'(curt+0.5*h,curx+0.5*h*k2,curu)']);k4=eval([StateModel,'(curt+h,curx+h*k3,curu)']);NewX=curx+h*(k1+2*k2+2*k3+k4)/6;elseifInteMethod=='RK2'k1=eval([StateModel,'(curt,curx,curu)']);k2=eval([StateModel,'(curt+h,curx+h*k1,curu)']);NewX=curx+0.5*h*(k1+k2);else%歐拉法EULQk=eval([StateModel,'(curt,curx,curu)']);%eval用來(lái)執(zhí)行一個(gè)函數(shù)，傳遞函數(shù)名和函數(shù)的輸入?yún)?shù)NewX=curx+h*Qk;end5functionNewX=w_StepIntegral(c函數(shù)w_DigiInteSimu和w_StepIntegral構(gòu)造了一個(gè)數(shù)值積分法仿真的框架，并不涉及具體的系統(tǒng)。具體的系統(tǒng)由StateModel,ControlFile,OutputFile參個(gè)參數(shù)決定，實(shí)際上就是三個(gè)函數(shù)文件名，這三個(gè)函數(shù)輸入輸出參數(shù)必須遵循特定的格式，在準(zhǔn)備好由這3個(gè)函數(shù)描述的系統(tǒng)后，調(diào)用w_DigiInteSimu即可進(jìn)行仿真。還需要一個(gè)調(diào)用w_DigiInteSimu進(jìn)行仿真的程序，它指定模型文件，指定初始參數(shù)，并且對(duì)仿真結(jié)果繪圖。6函數(shù)w_DigiInteSimu和w_StepIntegra3.2算法穩(wěn)定性分析仿真編程針對(duì)下面的系統(tǒng)，求用解析法、歐拉法和RK4分別求解，計(jì)算歐拉法最大允許步長(zhǎng)，將步長(zhǎng)從0.1逐漸增大，比較三種解的效果。解：1）該系統(tǒng)是穩(wěn)定的，解析解為2）用歐拉法計(jì)算本例時(shí)，其步長(zhǎng)應(yīng)該滿足3）RK4步長(zhǎng)條件式是一個(gè)高階不等式，無(wú)法直接求解，只能用試探法確定RK4的步長(zhǎng)上限。(3.2-1)73.2算法穩(wěn)定性分析仿真編程針對(duì)下面的系統(tǒng)，求用解析法、歐我們使用3.1介紹的兩個(gè)基本函數(shù)，同時(shí)要針對(duì)系統(tǒng)模型編寫3個(gè)函數(shù)來(lái)描述該系統(tǒng)，最后編寫一個(gè)實(shí)現(xiàn)仿真初值設(shè)置、仿真求解和仿真結(jié)果顯示的函數(shù)。functionderX=w_SysStateEqu(curt,curx,curu)%functionderX=w_stateEqu(curt,curx,curu)%函數(shù)功能：在此函數(shù)里寫系統(tǒng)的狀態(tài)方程%輸入?yún)?shù)：curt當(dāng)前時(shí)間，curx和curu是當(dāng)前狀態(tài)和控制%輸出參數(shù)：derX是狀態(tài)的X的導(dǎo)數(shù)值derX(1)=-4*curx(1);（1）狀態(tài)模型函數(shù)w_SysStateEqu在該函數(shù)里描述(3.2-1)式的狀態(tài)微分方程8我們使用3.1介紹的兩個(gè)基本函數(shù)，同時(shí)要針對(duì)系統(tǒng)模型編寫3個(gè)functionOutputY=w_SysOutput(curt,curx,curu)%functionOutputY=w_SysOutput(curt,curx,curu)%函數(shù)功能：計(jì)算系統(tǒng)的輸出向量%輸入?yún)?shù)：curt是當(dāng)前時(shí)間，curx,curu分別是當(dāng)前的狀態(tài)和控制向量%輸出參數(shù)：OutputY是計(jì)算出來(lái)的輸出向量OutputY(1)=curx(1);%根據(jù)具體的模型寫輸出方程（2）系統(tǒng)輸出函數(shù)w_SysOutputfunctionControlU=w_SysControl(curt,h,curx,curu,cury)%ControlU=w_sysControl(curt,h,curx,curu,cury)%函數(shù)功能：計(jì)算系統(tǒng)的控制，如u=PID(t,curx,cury)%輸入?yún)?shù)：curt是當(dāng)前時(shí)間，h是數(shù)值積分步長(zhǎng)％ curx,curu,cury分別是當(dāng)前的狀態(tài)、控制和輸出向量%輸出參數(shù)：ControlU是計(jì)算出來(lái)的控制向量ControlU=curu;%根據(jù)實(shí)際系統(tǒng)寫控制（3）系統(tǒng)控制計(jì)算函數(shù)w_SysControl9functionOutputY=w_SysOutput(c（4）仿真組織函數(shù)functionsimu_Stability%函數(shù)功能：穩(wěn)定性分析仿真tstart=0;tstop=7;h=0.7;StateModel='w_SysStateEqu';ControlFile='w_SysControl';OutputFile='w_SysOutput';x0=[1];u0=[0];cnty=1;[t,y1]=w_DigiInteSimu(tstart,tstop,h,x0,u0,cnty,'EUL',StateModel,OutputFile,ControlFile);%用歐拉法求解[t,y4]=w_DigiInteSimu(tstart,tstop,h,x0,u0,cnty,'RK4',StateModel,OutputFile,ControlFile);%用RK4求解y5=exp(-4*t);%計(jì)算解析解plot(t,y1,'--m',t,y2,'-.r',t,y5,'-k');%繪制疊加曲線xlabel('時(shí)間'); ylabel('y');title('算法穩(wěn)定性分析仿真');該函數(shù)的主要任務(wù)是：1)設(shè)置仿真起始時(shí)間、結(jié)束時(shí)間、仿真步長(zhǎng)、初始的x和u2)設(shè)置模型的輸出方程、狀態(tài)方程和控制方程所在函數(shù)文件名3)調(diào)用w_DigiInteSimu進(jìn)行仿真計(jì)算4)將計(jì)算結(jié)果繪圖進(jìn)行比較。10（4）仿真組織函數(shù)functionsimu_Stabili修改simu_Stability函數(shù)里的參數(shù)h，就可以作出不同步長(zhǎng)時(shí)的仿真結(jié)果。下面是步長(zhǎng)從0.1不斷增大時(shí)的仿真結(jié)果h=0.1時(shí)的仿真曲線。1】RK4法的曲線與解析解重合，說(shuō)明RK4法非常準(zhǔn)確。2】歐拉法有較大誤差。11修改simu_Stability函數(shù)里的參數(shù)h，就可以作出不h=0.4時(shí)的仿真曲線。1】RK4法還是比較準(zhǔn)確。2】歐拉法誤差很大，出現(xiàn)衰減震蕩。3】?jī)煞N數(shù)值算法的步長(zhǎng)條件仍然滿足，算法仍然穩(wěn)定。12h=0.4時(shí)的仿真曲線。12h=0.6時(shí)的仿真曲線。歐拉法的步長(zhǎng)條件已經(jīng)不滿足，仿真解發(fā)散h=0.6時(shí)的仿真曲線。RK4法的步長(zhǎng)條件仍然滿足。但誤差增大13h=0.6時(shí)的仿真曲線。歐拉法的步長(zhǎng)條件已經(jīng)不滿足，仿真解發(fā)h=0.7時(shí)的仿真曲線。RK4法的步長(zhǎng)條件已經(jīng)不滿足，仿真結(jié)果發(fā)散。14h=0.7時(shí)的仿真曲線。143.3衛(wèi)星發(fā)射仿真衛(wèi)星發(fā)射運(yùn)動(dòng)方程G為重力系數(shù)系統(tǒng)是高階微分方程，不能直接用于仿真計(jì)算，需要轉(zhuǎn)化為狀態(tài)方程。定義狀態(tài)變量：運(yùn)動(dòng)方程用極坐標(biāo)表示，同時(shí)還需要用直角坐標(biāo)輸出。1.模型轉(zhuǎn)換153.3衛(wèi)星發(fā)射仿真衛(wèi)星發(fā)射運(yùn)動(dòng)方程G為重力系數(shù)系統(tǒng)是高階微得到系統(tǒng)仿真模型取X-Y平面坐標(biāo)作為輸出狀態(tài)初值為：v是初始發(fā)射速度，可分別取該系統(tǒng)沒(méi)有控制，主要是研究初始發(fā)射速度對(duì)軌道的影響。16得到系統(tǒng)仿真模型取X-Y平面坐標(biāo)作為輸出狀態(tài)初值為：v是初始2.模型描述函數(shù)構(gòu)造functionderX=sat_StateEqu(curt,curx,curu)G=401408;derX(1)=curx(3);derX(2)=curx(4);derX(3)=-G/(curx(1)*curx(1))+curx(1)*curx(4)*curx(4);derX(4)=-2*curx(3)*curx(4)/curx(1);derX=derX’;%轉(zhuǎn)換為列向量(1)狀態(tài)方程functionOutputY=sat_Output(curt,curx,curu)OutputY(1)=curx(1)*cos(curx(2));OutputY(2)=curx(1)*sin(curx(2));OutputY=OutputY’;%轉(zhuǎn)換為列向量(2)輸出方程(3)控制方程由于該系統(tǒng)沒(méi)有控制，因而采用與3.2節(jié)同樣的控制函數(shù)。172.模型描述函數(shù)構(gòu)造functionderX=sat_S3.仿真組織函數(shù)functionsimu_Satelite%函數(shù)功能：對(duì)衛(wèi)星發(fā)射軌道仿真tstart=0;StateModel='sat_StateEqu';ControlFile='w_SysControl';OutputFile='sat_Output'; u0=[0];cnty=2;h=200;tstop=100*h; v=10;x0=[6400,0,0,v/6400]‘;[t,y10]=w_DigiInteSimu(tstart,tstop,h,x0,u0,cnty,‘RK4’,StateModel,OutputFile,ControlFile);%用RK4求解h=200;tstop=50*h; v=9;x0=[6400,0,0,v/6400]’;[t,y9]=w_DigiInteSimu(tstart,tstop,h,x0,u0,cnty,‘RK4’,StateModel,OutputFile,ControlFile);%用RK4求解h=100;tstop=60*h; v=8;x0=[6400,0,0,v/6400]‘;[t,y8]=w_DigiInteSimu(tstart,tstop,h,x0,u0,cnty,‘RK4’,StateModel,OutputFile,ControlFile);%用RK4求解plot(y10(1,:),y10(2,:),'-',y9(1,:),y9(2,:),'--',y8(1,:),y8(2,:),':');%繪制疊加曲線xlabel('x'); ylabel('y'); title('衛(wèi)星發(fā)射軌道');183.仿真組織函數(shù)functionsimu_Satelit直角坐標(biāo)顯示的衛(wèi)星發(fā)射軌道：初始速度越大，規(guī)大越大。V不同時(shí)，需要的步長(zhǎng)和計(jì)算步數(shù)不一樣。19直角坐標(biāo)顯示的衛(wèi)星發(fā)射軌道：初始速度越大，規(guī)大越大。193.4線性系統(tǒng)仿真編程前面所介紹的仿真程序是針對(duì)(3.1-1)所表示的通用的非線性系統(tǒng)，針對(duì)如下所表示的線性狀態(tài)空間模型(3.4-1)我們雖然也可以用前面介紹的程序仿真，但比較麻煩，實(shí)際上系統(tǒng)由矩陣A、B、C、D就可以指定微分方程和輸出方程，而不必單獨(dú)用函數(shù)文件來(lái)表示。為此，我們單獨(dú)編寫針對(duì)系統(tǒng)(3.4-1)的線性系統(tǒng)數(shù)值積分法仿真程序。程序仍然由仿真框架和單步積分構(gòu)成。203.4線性系統(tǒng)仿真編程前面所介紹的仿真程序是針對(duì)(3.1-(1)線性系統(tǒng)數(shù)值積分仿真框架function[t,y]=w_LinearSimu(tstart,tstop,h,x0,u0,A,B,C,D,InteMethod,ControlFile)%函數(shù)功能：用數(shù)值積分法仿真對(duì)線性系統(tǒng)dx=AX+Bu,y=Cx+Du進(jìn)行仿真%輸入?yún)?shù)：tstart,tstop,h分別是起始時(shí)間、結(jié)束時(shí)間和仿真步長(zhǎng)，是標(biāo)量%x0,u0是狀態(tài)、輸入的初值，都是列向量%A,B,C,D是線性狀態(tài)空間模型的系數(shù)矩陣%InteMethod是數(shù)值積分方法，可以是'EUL','RK2','RK4'%ControlFile是控制作用函數(shù)，沒(méi)有控制算法時(shí)，可以省略%輸出參數(shù)：t是仿真結(jié)果的時(shí)間序列%y是仿真結(jié)果系統(tǒng)的輸出序列21(1)線性系統(tǒng)數(shù)值積分仿真框架function[t,y]=function[t,y]=w_LinearSimu(tstart,tstop,h,x0,u0,A,B,C,D,InteMethod,ControlFile)t=[tstart:h:tstop];%t數(shù)一個(gè)行序列 cntt=size(t,2);%返回列數(shù)cnty=size(C,1);%返回y的維數(shù)y=zeros(cnty,cntt);%構(gòu)造一個(gè)空矩陣，用來(lái)存儲(chǔ)結(jié)果y0=C*x0+D*u0;%eval([OutputFile,'(tstart,x0,u0)']);%計(jì)算初始輸出y(:,1)=y0;%將y0作為輸出的第1列curx=x0;%當(dāng)前一步的x curu=u0;%當(dāng)前一步的u cury=y0;%當(dāng)前一步的yfori=1:1:cntt-1curu=eval([ControlFile,'(t(i),h,curx,curu,cury,A,B,C,D)']);%計(jì)算控制時(shí)傳遞的參數(shù)：當(dāng)前時(shí)間，步長(zhǎng)，當(dāng)前狀態(tài)和輸出,模型系數(shù)矩陣curx=w_LinearStepIntegral(h,curx,curu,InteMethod,A,B);%單步積分運(yùn)算,針對(duì)線性模型cury=C*curx+D*curu;%計(jì)算輸出y(:,i+1)=cury;%將輸出加入到輸出序列里end22function[t,y]=w_LinearSimu(ts(2)線性系統(tǒng)單步積分函數(shù)functionNewX=w_LinearStepIntegral(h,curx,curu,InteMethod,A,B);%函數(shù)功能：?jiǎn)尾椒e分運(yùn)算,與模型方程有關(guān)%輸入?yún)?shù)：h是數(shù)值積分步長(zhǎng)，curx,curu分別是當(dāng)前的狀態(tài)和控制向量%InteMethod是積分算法，字符串類型，% 可以取'EUL','RK2','RK4',由于要用于判斷，字符串必須等長(zhǎng)%A,B是線性狀態(tài)模型的微分方程的系數(shù)矩陣%輸出參數(shù)：NewX是這一步計(jì)算的新的狀態(tài)向量Bu=B*curu;ifInteMethod=='RK4'k1=A*curx+Bu;k2=A*(curx+0.5*h*k1)+Bu;k3=A*(curx+0.5*h*k2)+Bu;k4=A*(curx+h*k3)+Bu;NewX=curx+h*(k1+2*k2+2*k3+k4)/6;elseifInteMethod=='RK2'k1=A*curx+Bu;%eval([StateModel,'(curt,curx,curu)']);k2=A*(curx+h*k1)+Bu;%eval([StateModel,'(curt+h,curx+h*k1,curu)']);NewX=curx+0.5*h*(k1+k2);else%歐拉法EULNewX=curx+h*(A*curx+Bu);end23(2)線性系統(tǒng)單步積分函數(shù)functionNewX=w_L3.5二階系統(tǒng)傳遞函數(shù)仿真1.模型轉(zhuǎn)換寫為狀態(tài)模型就是基本參數(shù)取值為針對(duì)本例，線性系統(tǒng)仿真程序進(jìn)行仿真。243.5二階系統(tǒng)傳遞函數(shù)仿真1.模型轉(zhuǎn)換寫為狀態(tài)模型就是基本2.仿真組織程序functionsimu_StateSpace%函數(shù)功能：對(duì)一個(gè)狀態(tài)空間模型進(jìn)行仿真tstart=0;tstop=20;h=0.1;B=[0;1];C=[1,0];D=[0];x0=[0;0];u0=[1];%單位階躍輸入dampratio=0.1;A=[0,1;-1,-2*dampratio];[t,y1]=w_LinearSimu(tstart,tstop,h,x0,u0,A,B,C,D,'RK4');dampratio=0.5;A(2,2)=-2*dampratio;[t,y5]=w_LinearSimu(tstart,tstop,h,x0,u0,A,B,C,D,'RK4');dampratio=1;A(2,2)=-2*dampratio;[t,y10]=w_LinearSimu(tstart,tstop,h,x0,u0,A,B,C,D,'RK4');dampratio=1.5;A(2,2)=-2*dampratio;[t,y15]=w_LinearSimu(tstart,tstop,h,x0,u0,A,B,C,D,'RK4');plot(t,y1,'-r',t,y5,':g',t,y10,'*y',t,y15,'.b');%繪制疊加曲線xlabel('x');ylabel('y');title('二階系統(tǒng)階躍輸入響應(yīng)');252.仿真組織程序functionsimu_StateSpa二階系統(tǒng)在階躍輸入作用下，不同阻尼比時(shí)的輸出曲線。26二階系統(tǒng)在階躍輸入作用下，不同阻尼比時(shí)的輸出曲線。263.6面向結(jié)構(gòu)圖仿真對(duì)下圖的系統(tǒng)進(jìn)行面向結(jié)構(gòu)圖的變換，并進(jìn)行仿真。G1(s)G2(s)G3(s)G4(s)β其中273.6面向結(jié)構(gòu)圖仿真對(duì)下圖的系統(tǒng)進(jìn)行面向結(jié)構(gòu)圖的變換，并進(jìn)行(1)模型處理每個(gè)環(huán)節(jié)用3個(gè)參數(shù)來(lái)描述，得到的矩陣28(1)模型處理每個(gè)環(huán)節(jié)用3個(gè)參數(shù)來(lái)描述，得到的矩陣28(2)仿真組織程序利用第2章編寫的結(jié)構(gòu)圖到狀態(tài)空間模型轉(zhuǎn)換的函數(shù)w_bd2ss進(jìn)行模型處理functionsimu_BlockDiagram%函數(shù)功能：面向結(jié)構(gòu)圖的仿真G=[2,1,1;0,1,3;5,1,6;0,1,1];W=[0,0,0,0;1,0,-2,0;0,1,0,1;0,0,0,0];W0=[1;0;0;1];[P,Q]=w_bd2ss(G,W,W0);%先通過(guò)面向結(jié)構(gòu)圖的模型變換得到狀態(tài)模型tstart=0;tstop=5;h=0.1;C=[1,0,0,0;0,1,0,0;0,0,1,0;0,0,0,1];%44個(gè)狀態(tài)都輸出D=[0]; x0=[0;0;0;0]; u0=[1];%單位階躍輸入[t,y]=w_LinearSimu(tstart,tstop,h,x0,u0,P,Q,C,D,'RK4');%size(y)%plot(t,y(1,:),'-r',t,y(2,:),':g',t,y(3,:),'*m',t,y(4,:),'.b');%繪制疊加曲線plot(t,y(1,:),'-r',t,y(3,:),'*m');%繪制疊加曲線xlabel('time'); ylabel('y'); title('面向結(jié)構(gòu)圖的仿真');29(2)仿真組織程序functionsimu_BlockDG=[2,1,1;0,1,3;5,1,6;0,1,1];只繪制y1和y3，可以看出兩個(gè)輸出的區(qū)別，它們都趨向于穩(wěn)定。30G=[2,1,1;0,1,3;5,1,6;0,1,1];只繪同時(shí)繪制4個(gè)輸出時(shí)，由于y2,y4變化很大，發(fā)散，所以y1和y3看起來(lái)好像重疊了。31同時(shí)繪制4個(gè)輸出時(shí)，由于y2,y4變化很大，發(fā)散，所以y1和附錄：plot函數(shù)參數(shù)取值意義Plot函數(shù)的一般調(diào)用是：plot(x1,y1,option1,x2,y2,option2,….)32附錄：plot函數(shù)參數(shù)取值意義32第三章Matlab編程(數(shù)值積分法仿真)3.1連續(xù)系統(tǒng)數(shù)值積分法仿真編程思路目的：針對(duì)下面的系統(tǒng)編寫通用的數(shù)值積分法仿真程序(3.1-1)對(duì)這樣的系統(tǒng)進(jìn)行仿真，實(shí)際上涉及到控制的計(jì)算、狀態(tài)的數(shù)值積分計(jì)算和輸出的計(jì)算。3個(gè)函數(shù)g,f,h確定后，就可以完整地描述一個(gè)系統(tǒng)。給定初值：u0，y0,系統(tǒng)中的向量都采用列向量表示33第三章Matlab編程(數(shù)值積分法仿真)3.1連續(xù)系統(tǒng)數(shù)function[t,y]=w_DigiInteSimu(tstart,tstop,h,x0,u0,cnty,InteMethod,StateModel,OutputFile,ControlFile)%函數(shù)功能：用數(shù)值積分法仿真%輸入?yún)?shù)：tstart,tstop,h分別是起始時(shí)間、結(jié)束時(shí)間和仿真步長(zhǎng)，是標(biāo)量%x0,u0是狀態(tài)、輸入的初值，都是列向量%cnty是輸出變量的個(gè)數(shù)%InteMethod時(shí)數(shù)值積分方法，可以是'EUL','RK2','RK4'%StateModel是狀態(tài)模型的文件名%ControlFile是控制作用的文件名%OutputFile是系統(tǒng)輸出的文件名%輸出參數(shù)：t是仿真結(jié)果的時(shí)間序列%y是仿真結(jié)果系統(tǒng)的輸出序列(1)數(shù)值積分法仿真框架函數(shù)34function[t,y]=w_DigiInteSimu(function[t,y]=w_DigiInteSimu(tstart,tstop,h,x0,u0,cnty,InteMethod,StateModel,OutputFile,ControlFile)t=[tstart:h:tstop];%t數(shù)一個(gè)行序列cntt=size(t,2);%返回列數(shù)y=zeros(cnty,cntt);%構(gòu)造一個(gè)空矩陣，用來(lái)存儲(chǔ)結(jié)果y0=eval([OutputFile,'(tstart,x0,u0)']);%計(jì)算初始輸出y(:,1)=y0’;%將cury作為輸出的第1列curx=x0;%當(dāng)前一步的xcuru=u0;%當(dāng)前一步的ucury=y0;%當(dāng)前一步的yfori=1:1:cntt-1curu=eval([ControlFile,'(t(i),h,curx,curu,cury)']);%計(jì)算控制時(shí)傳遞的參數(shù)：當(dāng)前時(shí)間，步長(zhǎng)，當(dāng)前狀態(tài)和輸出curx=w_StepIntegral(t(i),h,curx,curu,InteMethod,StateModel);%單步積分運(yùn)算cury=eval([OutputFile,'(t(i),curx,curu)']);%計(jì)算輸出y(:,i+1)=cury‘;%將輸出加入到輸出序列里end35function[t,y]=w_DigiInteSimu(functionNewX=w_StepIntegral(curt,h,curx,curu,InteMethod,StateModel)%函數(shù)功能：?jiǎn)尾椒e分運(yùn)算,與模型方程有關(guān)%輸入?yún)?shù)：curt是當(dāng)前時(shí)間，h是數(shù)值積分步長(zhǎng)% curx,curu分別是當(dāng)前的狀態(tài)和控制向量%InteMethod是積分算法，字符串類型，可以取'EUL','RK2','RK4'%StateModel是狀態(tài)模型文件名稱，字符串類型%輸出參數(shù)：NewX是這一步計(jì)算的新的狀態(tài)向量(2)單步數(shù)值積分法函數(shù)單步數(shù)值積分函數(shù)只是對(duì)微分方程組StateModel進(jìn)行一步的計(jì)算，計(jì)算法由InteMethod參數(shù)指定，可以上歐拉法，RK2或RK4。36functionNewX=w_StepIntegral(cfunctionNewX=w_StepIntegral(curt,h,curx,curu,InteMethod,StateModel);ifInteMethod=='RK4'k1=eval([StateModel,'(curt,curx,curu)']);k2=eval([StateModel,'(curt+0.5*h,curx+0.5*h*k1,curu)']);k3=eval([StateModel,'(curt+0.5*h,curx+0.5*h*k2,curu)']);k4=eval([StateModel,'(curt+h,curx+h*k3,curu)']);NewX=curx+h*(k1+2*k2+2*k3+k4)/6;elseifInteMethod=='RK2'k1=eval([StateModel,'(curt,curx,curu)']);k2=eval([StateModel,'(curt+h,curx+h*k1,curu)']);NewX=curx+0.5*h*(k1+k2);else%歐拉法EULQk=eval([StateModel,'(curt,curx,curu)']);%eval用來(lái)執(zhí)行一個(gè)函數(shù)，傳遞函數(shù)名和函數(shù)的輸入?yún)?shù)NewX=curx+h*Qk;end37functionNewX=w_StepIntegral(c函數(shù)w_DigiInteSimu和w_StepIntegral構(gòu)造了一個(gè)數(shù)值積分法仿真的框架，并不涉及具體的系統(tǒng)。具體的系統(tǒng)由StateModel,ControlFile,OutputFile參個(gè)參數(shù)決定，實(shí)際上就是三個(gè)函數(shù)文件名，這三個(gè)函數(shù)輸入輸出參數(shù)必須遵循特定的格式，在準(zhǔn)備好由這3個(gè)函數(shù)描述的系統(tǒng)后，調(diào)用w_DigiInteSimu即可進(jìn)行仿真。還需要一個(gè)調(diào)用w_DigiInteSimu進(jìn)行仿真的程序，它指定模型文件，指定初始參數(shù)，并且對(duì)仿真結(jié)果繪圖。38函數(shù)w_DigiInteSimu和w_StepIntegra3.2算法穩(wěn)定性分析仿真編程針對(duì)下面的系統(tǒng)，求用解析法、歐拉法和RK4分別求解，計(jì)算歐拉法最大允許步長(zhǎng)，將步長(zhǎng)從0.1逐漸增大，比較三種解的效果。解：1）該系統(tǒng)是穩(wěn)定的，解析解為2）用歐拉法計(jì)算本例時(shí)，其步長(zhǎng)應(yīng)該滿足3）RK4步長(zhǎng)條件式是一個(gè)高階不等式，無(wú)法直接求解，只能用試探法確定RK4的步長(zhǎng)上限。(3.2-1)393.2算法穩(wěn)定性分析仿真編程針對(duì)下面的系統(tǒng)，求用解析法、歐我們使用3.1介紹的兩個(gè)基本函數(shù)，同時(shí)要針對(duì)系統(tǒng)模型編寫3個(gè)函數(shù)來(lái)描述該系統(tǒng)，最后編寫一個(gè)實(shí)現(xiàn)仿真初值設(shè)置、仿真求解和仿真結(jié)果顯示的函數(shù)。functionderX=w_SysStateEqu(curt,curx,curu)%functionderX=w_stateEqu(curt,curx,curu)%函數(shù)功能：在此函數(shù)里寫系統(tǒng)的狀態(tài)方程%輸入?yún)?shù)：curt當(dāng)前時(shí)間，curx和curu是當(dāng)前狀態(tài)和控制%輸出參數(shù)：derX是狀態(tài)的X的導(dǎo)數(shù)值derX(1)=-4*curx(1);（1）狀態(tài)模型函數(shù)w_SysStateEqu在該函數(shù)里描述(3.2-1)式的狀態(tài)微分方程40我們使用3.1介紹的兩個(gè)基本函數(shù)，同時(shí)要針對(duì)系統(tǒng)模型編寫3個(gè)functionOutputY=w_SysOutput(curt,curx,curu)%functionOutputY=w_SysOutput(curt,curx,curu)%函數(shù)功能：計(jì)算系統(tǒng)的輸出向量%輸入?yún)?shù)：curt是當(dāng)前時(shí)間，curx,curu分別是當(dāng)前的狀態(tài)和控制向量%輸出參數(shù)：OutputY是計(jì)算出來(lái)的輸出向量OutputY(1)=curx(1);%根據(jù)具體的模型寫輸出方程（2）系統(tǒng)輸出函數(shù)w_SysOutputfunctionControlU=w_SysControl(curt,h,curx,curu,cury)%ControlU=w_sysControl(curt,h,curx,curu,cury)%函數(shù)功能：計(jì)算系統(tǒng)的控制，如u=PID(t,curx,cury)%輸入?yún)?shù)：curt是當(dāng)前時(shí)間，h是數(shù)值積分步長(zhǎng)％ curx,curu,cury分別是當(dāng)前的狀態(tài)、控制和輸出向量%輸出參數(shù)：ControlU是計(jì)算出來(lái)的控制向量ControlU=curu;%根據(jù)實(shí)際系統(tǒng)寫控制（3）系統(tǒng)控制計(jì)算函數(shù)w_SysControl41functionOutputY=w_SysOutput(c（4）仿真組織函數(shù)functionsimu_Stability%函數(shù)功能：穩(wěn)定性分析仿真tstart=0;tstop=7;h=0.7;StateModel='w_SysStateEqu';ControlFile='w_SysControl';OutputFile='w_SysOutput';x0=[1];u0=[0];cnty=1;[t,y1]=w_DigiInteSimu(tstart,tstop,h,x0,u0,cnty,'EUL',StateModel,OutputFile,ControlFile);%用歐拉法求解[t,y4]=w_DigiInteSimu(tstart,tstop,h,x0,u0,cnty,'RK4',StateModel,OutputFile,ControlFile);%用RK4求解y5=exp(-4*t);%計(jì)算解析解plot(t,y1,'--m',t,y2,'-.r',t,y5,'-k');%繪制疊加曲線xlabel('時(shí)間'); ylabel('y');title('算法穩(wěn)定性分析仿真');該函數(shù)的主要任務(wù)是：1)設(shè)置仿真起始時(shí)間、結(jié)束時(shí)間、仿真步長(zhǎng)、初始的x和u2)設(shè)置模型的輸出方程、狀態(tài)方程和控制方程所在函數(shù)文件名3)調(diào)用w_DigiInteSimu進(jìn)行仿真計(jì)算4)將計(jì)算結(jié)果繪圖進(jìn)行比較。42（4）仿真組織函數(shù)functionsimu_Stabili修改simu_Stability函數(shù)里的參數(shù)h，就可以作出不同步長(zhǎng)時(shí)的仿真結(jié)果。下面是步長(zhǎng)從0.1不斷增大時(shí)的仿真結(jié)果h=0.1時(shí)的仿真曲線。1】RK4法的曲線與解析解重合，說(shuō)明RK4法非常準(zhǔn)確。2】歐拉法有較大誤差。43修改simu_Stability函數(shù)里的參數(shù)h，就可以作出不h=0.4時(shí)的仿真曲線。1】RK4法還是比較準(zhǔn)確。2】歐拉法誤差很大，出現(xiàn)衰減震蕩。3】?jī)煞N數(shù)值算法的步長(zhǎng)條件仍然滿足，算法仍然穩(wěn)定。44h=0.4時(shí)的仿真曲線。12h=0.6時(shí)的仿真曲線。歐拉法的步長(zhǎng)條件已經(jīng)不滿足，仿真解發(fā)散h=0.6時(shí)的仿真曲線。RK4法的步長(zhǎng)條件仍然滿足。但誤差增大45h=0.6時(shí)的仿真曲線。歐拉法的步長(zhǎng)條件已經(jīng)不滿足，仿真解發(fā)h=0.7時(shí)的仿真曲線。RK4法的步長(zhǎng)條件已經(jīng)不滿足，仿真結(jié)果發(fā)散。46h=0.7時(shí)的仿真曲線。143.3衛(wèi)星發(fā)射仿真衛(wèi)星發(fā)射運(yùn)動(dòng)方程G為重力系數(shù)系統(tǒng)是高階微分方程，不能直接用于仿真計(jì)算，需要轉(zhuǎn)化為狀態(tài)方程。定義狀態(tài)變量：運(yùn)動(dòng)方程用極坐標(biāo)表示，同時(shí)還需要用直角坐標(biāo)輸出。1.模型轉(zhuǎn)換473.3衛(wèi)星發(fā)射仿真衛(wèi)星發(fā)射運(yùn)動(dòng)方程G為重力系數(shù)系統(tǒng)是高階微得到系統(tǒng)仿真模型取X-Y平面坐標(biāo)作為輸出狀態(tài)初值為：v是初始發(fā)射速度，可分別取該系統(tǒng)沒(méi)有控制，主要是研究初始發(fā)射速度對(duì)軌道的影響。48得到系統(tǒng)仿真模型取X-Y平面坐標(biāo)作為輸出狀態(tài)初值為：v是初始2.模型描述函數(shù)構(gòu)造functionderX=sat_StateEqu(curt,curx,curu)G=401408;derX(1)=curx(3);derX(2)=curx(4);derX(3)=-G/(curx(1)*curx(1))+curx(1)*curx(4)*curx(4);derX(4)=-2*curx(3)*curx(4)/curx(1);derX=derX’;%轉(zhuǎn)換為列向量(1)狀態(tài)方程functionOutputY=sat_Output(curt,curx,curu)OutputY(1)=curx(1)*cos(curx(2));OutputY(2)=curx(1)*sin(curx(2));OutputY=OutputY’;%轉(zhuǎn)換為列向量(2)輸出方程(3)控制方程由于該系統(tǒng)沒(méi)有控制，因而采用與3.2節(jié)同樣的控制函數(shù)。492.模型描述函數(shù)構(gòu)造functionderX=sat_S3.仿真組織函數(shù)functionsimu_Satelite%函數(shù)功能：對(duì)衛(wèi)星發(fā)射軌道仿真tstart=0;StateModel='sat_StateEqu';ControlFile='w_SysControl';OutputFile='sat_Output'; u0=[0];cnty=2;h=200;tstop=100*h; v=10;x0=[6400,0,0,v/6400]‘;[t,y10]=w_DigiInteSimu(tstart,tstop,h,x0,u0,cnty,‘RK4’,StateModel,OutputFile,ControlFile);%用RK4求解h=200;tstop=50*h; v=9;x0=[6400,0,0,v/6400]’;[t,y9]=w_DigiInteSimu(tstart,tstop,h,x0,u0,cnty,‘RK4’,StateModel,OutputFile,ControlFile);%用RK4求解h=100;tstop=60*h; v=8;x0=[6400,0,0,v/6400]‘;[t,y8]=w_DigiInteSimu(tstart,tstop,h,x0,u0,cnty,‘RK4’,StateModel,OutputFile,ControlFile);%用RK4求解plot(y10(1,:),y10(2,:),'-',y9(1,:),y9(2,:),'--',y8(1,:),y8(2,:),':');%繪制疊加曲線xlabel('x'); ylabel('y'); title('衛(wèi)星發(fā)射軌道');503.仿真組織函數(shù)functionsimu_Satelit直角坐標(biāo)顯示的衛(wèi)星發(fā)射軌道：初始速度越大，規(guī)大越大。V不同時(shí)，需要的步長(zhǎng)和計(jì)算步數(shù)不一樣。51直角坐標(biāo)顯示的衛(wèi)星發(fā)射軌道：初始速度越大，規(guī)大越大。193.4線性系統(tǒng)仿真編程前面所介紹的仿真程序是針對(duì)(3.1-1)所表示的通用的非線性系統(tǒng)，針對(duì)如下所表示的線性狀態(tài)空間模型(3.4-1)我們雖然也可以用前面介紹的程序仿真，但比較麻煩，實(shí)際上系統(tǒng)由矩陣A、B、C、D就可以指定微分方程和輸出方程，而不必單獨(dú)用函數(shù)文件來(lái)表示。為此，我們單獨(dú)編寫針對(duì)系統(tǒng)(3.4-1)的線性系統(tǒng)數(shù)值積分法仿真程序。程序仍然由仿真框架和單步積分構(gòu)成。523.4線性系統(tǒng)仿真編程前面所介紹的仿真程序是針對(duì)(3.1-(1)線性系統(tǒng)數(shù)值積分仿真框架function[t,y]=w_LinearSimu(tstart,tstop,h,x0,u0,A,B,C,D,InteMethod,ControlFile)%函數(shù)功能：用數(shù)值積分法仿真對(duì)線性系統(tǒng)dx=AX+Bu,y=Cx+Du進(jìn)行仿真%輸入?yún)?shù)：tstart,tstop,h分別是起始時(shí)間、結(jié)束時(shí)間和仿真步長(zhǎng)，是標(biāo)量%x0,u0是狀態(tài)、輸入的初值，都是列向量%A,B,C,D是線性狀態(tài)空間模型的系數(shù)矩陣%InteMethod是數(shù)值積分方法，可以是'EUL','RK2','RK4'%ControlFile是控制作用函數(shù)，沒(méi)有控制算法時(shí)，可以省略%輸出參數(shù)：t是仿真結(jié)果的時(shí)間序列%y是仿真結(jié)果系統(tǒng)的輸出序列53(1)線性系統(tǒng)數(shù)值積分仿真框架function[t,y]=function[t,y]=w_LinearSimu(tstart,tstop,h,x0,u0,A,B,C,D,InteMethod,ControlFile)t=[tstart:h:tstop];%t數(shù)一個(gè)行序列 cntt=size(t,2);%返回列數(shù)cnty=size(C,1);%返回y的維數(shù)y=zeros(cnty,cntt);%構(gòu)造一個(gè)空矩陣，用來(lái)存儲(chǔ)結(jié)果y0=C*x0+D*u0;%eval([OutputFile,'(tstart,x0,u0)']);%計(jì)算初始輸出y(:,1)=y0;%將y0作為輸出的第1列curx=x0;%當(dāng)前一步的x curu=u0;%當(dāng)前一步的u cury=y0;%當(dāng)前一步的yfori=1:1:cntt-1curu=eval([ControlFile,'(t(i),h,curx,curu,cury,A,B,C,D)']);%計(jì)算控制時(shí)傳遞的參數(shù)：當(dāng)前時(shí)間，步長(zhǎng)，當(dāng)前狀態(tài)和輸出,模型系數(shù)矩陣curx=w_LinearStepIntegral(h,curx,curu,InteMethod,A,B);%單步積分運(yùn)算,針對(duì)線性模型cury=C*curx+D*curu;%計(jì)算輸出y(:,i+1)=cury;%將輸出加入到輸出序列里end54function[t,y]=w_LinearSimu(ts(2)線性系統(tǒng)單步積分函數(shù)functionNewX=w_LinearStepIntegral(h,curx,curu,InteMethod,A,B);%函數(shù)功能：?jiǎn)尾椒e分運(yùn)算,與模型方程有關(guān)%輸入?yún)?shù)：h是數(shù)值積分步長(zhǎng)，curx,curu分別是當(dāng)前的狀態(tài)和控制向量%InteMethod是積分算法，字符串類型，% 可以取'EUL','RK2','RK4',由于要用于判斷，字符串必須等長(zhǎng)%A,B是線性狀態(tài)模型的微分方程的系數(shù)矩陣%輸出參數(shù)：NewX是這一步計(jì)算的新的狀態(tài)向量Bu=B*curu;ifInteMethod=='RK4'k1=A*curx+Bu;k2=A*(curx+0.5*h*k1)+Bu;k3=A*(curx+0.5*h*k2)+Bu;k4=A*(curx+h*k3)+Bu;NewX=curx+h*(k1+2*k2+2*k3+k4)/6;elseifInteMethod=='RK2'k1=A*curx+Bu;%eval([StateModel,'(curt,curx,curu