系統(tǒng)仿真報(bào)告一

1、.實(shí)驗(yàn)一經(jīng)典的連續(xù)系統(tǒng)仿真建模方法一 實(shí)驗(yàn)?zāi)康? 了解和掌握利用仿真技術(shù)對(duì)控制系統(tǒng)進(jìn)行分析的原理和步驟。2 掌握機(jī)理分析建模方法。3 深入理解一階常微分方程組數(shù)值積分解法的原理和程序結(jié)構(gòu)，學(xué)習(xí)用Matlab編寫數(shù)值積分法仿真程序。4 掌握和理解四階Runge-Kutta法，加深理解仿真步長(zhǎng)與算法穩(wěn)定性的關(guān)系。二 實(shí)驗(yàn)內(nèi)容1. 編寫四階 Runge_Kutta 公式的計(jì)算程序，對(duì)非線性模型（ 3）式進(jìn)行仿真。（1） 將閥位 u 增大 10 和減小 10 ，觀察響應(yīng)曲線的形狀;（2） 研究仿真步長(zhǎng)對(duì)穩(wěn)定性的影響，仿真步長(zhǎng)取多大時(shí)RK4算法變得不穩(wěn)定？（3） 利用MATLAB中的 ode45() 函

2、數(shù)進(jìn)行求解，比較與（1）中的仿真結(jié)果有何區(qū)別。2. 編寫四階 Runge_Kutta 公式的計(jì)算程序，對(duì)線性狀態(tài)方程（ 18 ）式進(jìn)行仿真（ 1）將閥位增大 10 和減小 10 ，觀察響應(yīng)曲線的形狀 ;（2） 研究仿真步長(zhǎng)對(duì)穩(wěn)定性的影響，仿真步長(zhǎng)取多大時(shí)RK4算法變得不穩(wěn)定？（ 4） 閥位增大 10 和減小 10 ，利用 MATLAB 中的 ode45() 函數(shù)進(jìn)行求解階躍響應(yīng)，比較與（ 1）中的仿真結(jié)果有何區(qū)別。三 實(shí)驗(yàn)代碼及結(jié)果編寫四階Runge_Kutta公式的計(jì)算程序，對(duì)非線性模型（3）式進(jìn)行仿真。function dY=f(Y,u)%f 函數(shù)k=0.2;Qd=0.15;A=2;a1=

3、0.20412;a2=0.21129;dY=zeros(2,1);dY(1)=1/A*(k*u+Qd-a1*sqrt(Y(1)dY(2)=1/A*(a1*sqrt(Y(1)-a2*sqrt(Y(2)%RK4 文件clccloseY=1.2,1.4'u=0.45; Y=0.5;TT=;XX=;for i=1:Y:100k1=f(Y,u);k2=f(Y+Y*k1/2,u);k3=f(Y+Y*k2/2,u);k4=f(Y+Y*k3,u);.Y=Y+Y*(k1+2*k2+2*k3+k4)/6;TT=TT i;XX=XX Y;end;plot(TT,XX,'-.')xlabel(

4、'x')ylabel('Y')仿真曲線如下1.551.5H11.451.4HH21.351.31.25020406080100120140160180200timeu=0.5穩(wěn)態(tài)值Y1=1.2 ，Y2=1.4 。 u 增大 10 和減小10之仿真曲線如下1.41.38H11.361.341.32H 1.31.28H21.261.241.221.2020406080100120140160180200timeu=0.451.651.61.551.51.45H1.41.351.31.250u=0.552.62.42.22H1.81.61.40.H1H22040608

5、0100120140160180200timeH2H120406080100120140160180200time.步長(zhǎng)為60 ， u=0.5 仿真曲線可見步長(zhǎng)越大越不穩(wěn)定采用 ode45 算法程序如下functiondY=f(Y,u)k=0.2;u=0.5;Qd=0.15;A=2;a1=0.20412;a2=0.21129;dY=zeros(2,1);dY(1)=1/A*(k*u+Qd-a1*sqrt(Y(1)dY(2)=1/A*(a1*sqrt(Y(1)-a2*sqrt(Y(2)T,Y = ode45('f',1,200,1.2,1.4);%在命令窗口運(yùn)行以下程序plot(

6、T,Y(:,1),'-',T,Y(:,2),'-')1.41.38H11.361.341.32H1.31.28H21.261.241.221.2020406080100120140160180200time仿真曲線與四階 Runge_Kutta公式的計(jì)算一致編寫四階Runge_Kutta公式的計(jì)算程序，對(duì)線性狀態(tài)方程（18 ）式進(jìn)行仿真functiondY=f1(Y,u) %f1函數(shù)k=0.2;Qd=0.00001;A=2;a1=0.20412;a2=0.21129;R1=2*sqrt(1.5)/a1;.R2=2*sqrt(1.4)/a2;dY=zeros(2,

7、1);dY(1)=k/A*u+1/A*Qd-1/(A*R1)*Y(1)dY(2)=1/(A*R1)*Y(1)-1/(A*R2)*Y(2)%RK4clccloseY=0.001,0.001'u=0.00001;Y=0.1;TT=;XX=;fori=1:Y:200k1=f1(Y,u);k2=f1(Y+Y*k1/2,u);k3=f1(Y+Y*k2/2,u);k4=f1(Y+Y*k3,u);Y=Y+Y*(k1+2*k2+2*k3+k4)/6;TT=TT i;XX=XX Y;end ;Yoldonplot(TT,XX(1,:),'-' ,TT,XX(2,:);xlabel(

8、9;time')ylabel('Y')gtext('Y1')gtext('Y2')Yoldon仿真曲線如下.-4x 10109876H2HH154321020406080100120140160180200timeu=0.00001u 增大 10 和減小10 之仿真曲線如下-4x 1010987H26HH154321020406080100120140160180200timeu=0.000011.-4x 101098H27H16H54321020406080100120140160180200timeu=0.000009當(dāng)步長(zhǎng)40 時(shí)，

9、曲線如下，可見步長(zhǎng)越大越不穩(wěn)定-4x 103H1H2H21020406080100120140160180time采用 ode45 函數(shù)求解程序如下functiondY=f1(Y,u).k=0.2;Qd=0.00001;A=2;a1=0.20412;a2=0.21129;R1=2*sqrt(1.5)/a1;R2=2*sqrt(1.4)/a2;dY=zeros(2,1);dY(1)=k/A*u+1/A*Qd-1/(A*R1)*Y(1)dY(2)=1/(A*R1)*Y(1)-1/(A*R2)*Y(2)T,Y = ode45('f',1,200,0.00001,0.00001);plot(T,Y(:,1),'-',T,Y(:,2),'-')-4x 10109876H2HH1543210204060801001201