实验一 经典的连续系统仿真建模方法
一 实验目的
1． 了解和掌握利用仿真技术对控制系统进行分析的原理和步骤。2． 掌握机理分析建模方法。
3． 深入理解一阶常微分方程组数值积分解法的原理和程序结构，学习用Matlab 编写
数值积分法仿真程序。
4． 掌握和理解四阶 Runge-Kutta 法，加深理解仿真步长与算法稳定性的关系。
二 实验内容
1. 编写四阶 Runge_Kutta 公式的计算程序，对非线性模型（3）式进行仿真。
（1） 将阀位u 增大10％和减小10％，观察响应曲线的形状;（2） 研究仿真步长对稳定性的影响，仿真步长取多大时RK4 算法变得不稳定？
（3） 利用 MATLAB 中的ode45()函数进行求解，比较与（1）中的仿真结果有何区别。
2. 编写四阶 Runge_Kutta 公式的计算程序，对线性状态方程（18）式进行仿真
（1） 将阀位增大10％和减小10％，观察响应曲线的形状;
（2） 研究仿真步长对稳定性的影响，仿真步长取多大时RK4 算法变得不稳定？
（4） 阀位增大10％和减小10％，利用MATLAB 中的ode45()函数进行求解阶跃响
应，比较与（1）中的仿真结果有何区别。三 实验代码及结果
编写四阶 Runge_Kutta 公式的计算程序，对非线性模型（3）式进行仿真。
function dY=f(Y,u) %f函数k=0.2;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)))%RK4文件 clcclose
Y=[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('x')ylabel('Y')仿真曲线如下
1.551.5H11.451.4HH21.351.31.25020406080100time120140160180200u=0.5稳态值Y1=1.2，Y2=1.4。 u 增大10％和减小10％之仿真曲线如下1.41.38H11.361.341.321.31.281.261.241.221.2020406080100time120140160180200H2Hu=0.451.651.61.551.51.45H2H1H1.41.351.31.25020406080100time120140160180200u=0.552.62.42.2H22H1.81.6H11.4020406080100time120140160180200步长为60，u=0.5仿真曲线可见步长越大越不稳定采用ode45算法程序如下function dY=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]);%在命令窗口运行以下程序plot(T,Y(:,1),'-',T,Y(:,2),'--')
1.41.38H11.361.341.321.31.281.261.241.221.2020406080100time120140160180200H2仿真曲线与四阶 Runge_Kutta 公式的计算一致H编写四阶 Runge_Kutta 公式的计算程序，对线性状态方程（18）式进行仿真function dY=f1(Y,u) %f1函数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)%RK4 clcclose
Y=[0.001,0.001]';u=0.00001;Y=0.1;TT=[];XX=[];for i=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;Yold on
plot(TT,XX(1,:),'--',TT,XX(2,:));xlabel('time')
ylabel('Y')gtext('Y1')gtext('Y2')Yold on仿真曲线如下-4109876x 10H2H1H54321020406080100time120140160180200u=0.00001u 增大10％和减小10％之仿真曲线如下109876x 10-4H2H1H54321020406080100time120140160180200u=0.000011109876H1x 10-4H2H54321020406080100time120140160180200u=0.000009当步长40时，曲线如下，可见步长越大越不稳定x 10-43H1H2H21020406080time100120140160180采用ode45函数求解程序如下function dY=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),'--')109876H2H1x 10-4H54321020406080100time120140160180200仿真曲线与四阶 Runge_Kutta 公式的计算一致五思考题1． 讨论仿真步长对稳定性和仿真精度的影响。答：步长越大越不稳定，但仿真速度快2． 你是怎样实现阀位增大和减小10%的？对于非线性模型和线性模型方法一样吗？答：u*1.1和u*0.9，两种方法一样
