经典的连续系统仿真建模方法实验报告

1、实验一经典 的连续系统仿真建模方法一 实验目的:1 了解和掌握利用仿真技术对控制系统进行分析的原理和步骤。2 掌握机理分析建模方法。3 深入理解阶常微分方程组数值积分解法的原理和程序结构，学习用matlab编写 数值积分法仿真程序。4 掌握和理解四阶runge-kutta法，加深理解仿真步长与算法稳定性的关系。二 实验原理:1非线性模型仿真三 实验内容:1. 编写四阶 runge_kutta 公式的计算程序，对非线性模型（3）式进行仿真。（1） 将阀位u 增大10和减小10，观察响应曲线的形状;（2） 研究仿真步长对稳定性的影响，仿真步长取多大时rk4 算法变得不稳定？（3） 利用 matla

2、b 中的ode45()函数进行求解，比较与（1）中的仿真结果有何区别。2. 编写四阶 runge_kutta 公式的计算程序，对线性状态方程（18）式进行仿真（1） 将阀位增大10和减小10，观察响应曲线的形状;（2） 研究仿真步长对稳定性的影响，仿真步长取多大时rk4 算法变得不稳定？（4） 阀位增大10和减小10，利用matlab 中的ode45()函数进行求解阶跃响应，比较与（1）中的仿真结果有何区别。四 程序代码:龙格库塔:%rk4文件 clccloseh=1.2,1.4;u=0.55; h=1;tt=;xx=;for i=1:h:200k1=f(h,u);k2=f(h+h*k1/2,

3、u);k3=f(h+h*k2/2,u);k4=f(h+h*k3,u);h=h+h*(k1+2*k2+2*k3+k4)/6;tt=tt i;xx=xx h;end;hold onplot(tt,xx(1,:),-,tt,xx(2,:);xlabel(time)ylabel(h)gtext(h1)gtext(h2)hold on水箱模型:function dh=f(h,u)k=0.2;u=0.5;qd=0.15;a=2;a1=0.20412;a2=0.21129;dh=zeros(2,1);dh(1)=1/a*(k*u+qd-a1*sqrt(h(1);dh(2)=1/a*(a1*sqrt(h(1)

4、-a2*sqrt(h(2);三 实验结果:2编写四阶 runge_kutta 公式的计算程序，对线性状态方程（18）式进行仿真：1 阀值u对仿真结果的影响u=0.45;h=1; u=0.5;h=1; u=0.55;h=1;2 步长h对仿真结果的影响:u=0.5;h=5; u=0.5;h=20; u=0.5;h=39 u=0.5;h=50 由以上结果知,仿真步长越大,仿真结果越不稳定。采用ode45算法程序如下:function dh=liu(t,h)k=0.2;u=0.45;qd=0.15;a=2;a1=0.20412;a2=0.21129;dh=zeros(2,1);dh(1)=1/a*(k

5、*u+qd-a1*sqrt(h(1);dh(2)=1/a*(a1*sqrt(h(1)-a2*sqrt(h(2);在命令窗口运行以下程序:t,h=ode45(liu,1 200,1.2 1.1);plot(t,h(:,1),r,+,t,h(:,2),g,*)u=0.45 u=0.5 u=0.55用ode45与用龙格库塔法仿真结果基本一致。2编写四阶 runge_kutta 公式的计算程序，对线性状态方程（18）式进行仿真：%rk4文件 clcclearclosex=1.2,1.4;u=0.5; h=5;tt=;xx=;for i=1:h:200k1=f2(x,u);k2=f2(x+h*k1/2,

6、u);k3=f2(x+h*k2/2,u);k4=f2(x+h*k3,u);x=x+h*(k1+2*k2+2*k3+k4)/6;tt=tt i;xx=xx x;end;hold onplot(tt,xx(1,:),-,tt,xx(2,:);xlabel(time)ylabel(x)gtext(x1)gtext(x2)hold on线性函数：function dx=f2(x,u)%线性qd=0.1;a1=0.20412;a2=0.21129;a=2;k=0.2;dx=zeros(2,1);dx(1)=1/a*(k*u-x(1)/(2*sqrt(1.5)/a1)+qd);dx(2)=1/a*(x(1

7、)/(2*sqrt(1.5)/a1)-x(2)/(2*sqrt(1.4)/a2);1 阀值u对仿真结果的影响:u=0.45；h=1； u=0.5；h=1； u=0.55；h=1；2 步长h对仿真结果的影响:u=0.5；h=5； u=0.5；h=20； u=0.5；h=35； u=0.5；h=50； 当步长为50时仿真结果开始不稳定，仿真步长越大,仿真结果越不稳定。采用ode45算法程序如下:function dx=liu2(t,x)%qd=0.1;u=0.45;a1=0.20412;a2=0.21129;a=2;k=0.2;dx=zeros(2,1);dx(1)=1/a*(k*u-x(1)/(2*sqrt(1.5)/a1)+qd);dx(2)=1/a*(x(1)/(2*sqrt(1.5)/a1)-x(2)/(2*sqrt(1.4)/a2);在命令窗口运行以下程序: t,x=ode45(liu2,1 200,1.2 1.4); plot(t,x(:,1),r