实验一经典的连续系统仿真建模方法
一实验目的
1．了解和掌握利用仿真技术对控制系统进行分析的原理和步骤。
2．掌握机理分析建模方法。
3．深入理解一阶常微分方程组数值积分解法的原理和程序结构，学习用Matlab 编写
数值积分法仿真程序。
4．掌握和理解四阶Runge-Kutta 法，加深理解仿真步长与算法稳定性的关系。
二实验内容
1. 编写四阶Runge_Kutta 公式的计算程序，对非线性模型（3）式进行仿真。（1）将阀位u 增大10％和减小10％，观察响应曲线的形状;
u=0.45时的图像：
u=0.55
01002003004005006007008009001000 1.25
1.3
1.35
1.4
1.45
1.5
开大或关小阀位之后，稳态值会相应的从原液位上升或下降，这是符合实际的。
（2） 研究仿真步长对稳定性的影响，仿真步长取多大时RK4 算法变得不稳定？ 由（1）可知，当步长为40时，仿真结果是稳定的 当步长为80时的图像
01002003004005006007008009001000
1.35
1.4
1.45
1.5
1.55
1.6
1.65
1.7
1
2
3
4
5
6
7
h （1,1）的数值稳定，但是并不是实际求得的稳态值。h （1,2）的值显然发散。 进一步取小步长，取hstep=42时，图像出现偏差，但是稳态值不变
Hstep=65时，图像偏差明显
0200400600800100012001400160018002000
-140
-120-100-80-60-40-200
20020040060080010001200
1.35
1.41.451.51.551.61.651.7
1.75
而hsetp=65.7时，图像就发散了
（3） 利用 MATLAB 中的ode45()函数进行求解，比较与（1）中的仿真结果有何区别。 Ode45调用的函数：
function [dh]=daoshu(t,x)
020040060080010001200140016001800
0.5
1
1.5
2
2.5
3
020040060080010001200140016001800
-25
-20
-15
-10
-5
5
u(1)=0.4;%此处可以修改阀位
dh=zeros(2,1);
u(2)=0.15;
A=2;
ku=0.1/0.5;
alpha12 = 0.25/sqrt(1.5);
alpha2 = 0.25/sqrt(1.4);
dh(1,1)=(ku*(u(1))+u(2)-alpha12*sqrt(x(1,1)))/A;
dh(2,1)=(alpha12*sqrt(x(1,1))-alpha2*sqrt(x(2,1)))/A;
end
在主程序中添加ode45的算法：
[T,Y]=ode45('daoshu',[0,1000],[1.5,1.4]);
figure(1)
plot(T,Y)
hold on
plot([0:hStep:nCounter*hStep]',Hlevel)
grid
将ode45与编写的龙格库塔算法画到同一个坐标系中（点表示ode45）：
01002003004005006007008009001000
可以发现，ode45更快达到稳定值，实际上，在缩小了编写算法中的步长后，两种算法的曲线基本重合，说明ode45的精度很高：
01002003004005006007008009001000 2. 编写四阶Runge_Kutta 公式的计算程序，对线性状态方程（18）式进行仿真
（1）将阀位增大10％和减小10％，观察响应曲线的形状;
增大10％，令u(1)=0.05
令u(1)=-0.05
（2） 研究仿真步长对稳定性的影响，仿真步长取多大时RK4 算法变得不稳定？
方法同上，大约在hstep=61时就不稳定了
（3） 阀位增大10％和减小10％，利用MATLAB 中的ode45()函数进行求解阶跃响
050100150200250
1.35
1.4
1.45
1.5
1.55
1.6
1.65
1.7
050100150200250
1.25
1.3
1.35
1.4
1.45
1.5
减小10%（实线为编写算法，点为龙格库塔法）：
050100150200250增大10%：
050100150200250程序实现：
调用微分函数：
function [dh] = dx( t,x )
%DX 被ode45调用的微分方程函数
u=zeros(2,1);
u(1)=0.05; %此处可以修改阀位
u(2)=0;
A=2;
ku=0.1/0.5;
alpha12 = 0.25/sqrt(1.5);
alpha2 = 0.25/sqrt(1.4);
R12=2*sqrt(1.5)/alpha12;
R2=2*sqrt(1.4)/alpha2;
AMTRIX = [-1/(A*R12) 0;1/(A*R12) -1/(A*R2)];
BMATRIX = [ku/A 1/A;0 0];
dh=AMTRIX*x+BMATRIX*u;
end
主程序中添加：
[T,Y]=ode45('dx',[0,250],[0,0]);
Y(:,1)=Y(:,1)+1.5
Y(:,2)=Y(:,2)+1.4
figure(3)
plot(T,Y,'.')
hold on
plot([0:hStep:nCounter*hStep]',Hlevel)
grid
三实验报告
实验完成后，要写出实验报告，内容包括：
1．实验步骤及说明；
2．实验所用的仿真程序清单，以及程序结构的简单说明；
3．实验结果曲线及分析，稳态值是多少；
四思考题
1．讨论仿真步长对稳定性和仿真精度的影响。
一般来说，仿真步长越长，系统越倾向于不稳定，仿真步长越短，系统越稳定，但是需要的仿真时间越长。
2．你是怎样实现阀位增大和减小10%的？对于非线性模型和线性模型方法一样吗？
线性模型与非线性模型实际改变方式不一样。
线性模型中，由于采用的是增量化的方程，所以在修改时只需要修改增量，即0.05或-
0.05。
非线性模型中，采用的绝对量的方程，所以在修改时应该修改绝对量，即改为0.45或
0.55。
五实验总结
通过本次实验，熟悉了对matlab软件的使用，加深了对微分方程数值积分解法的理
