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

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

一实验目的:

1 了解和掌握利用仿真技术对控制系统进行分析的原理和步骤。

2 掌握机理分析建模方法。

3 深入理解阶常微分方程组数值积分解法的原理和程序结构，学习用Matlab编写

数值积分法仿真程序。

4 掌握和理解四阶Runge-Kutta法，加深理解仿真步长与算法稳定性的关系。

二实验原理:

1非线性模型仿真

三实验内容:

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）中的仿真结果有何区别。

四程序代码:

龙格库塔:

%RK4文件

clc

close

H=[1.2,1.4]';u=0.55; h=1;

TT=[];

XX=[];

for i=1:h:200

k1=f(H,u);

k2=f(H+h*k1/2,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 on

plot(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))-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;

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))-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文件

clc

x=[1.2,1.4]';u=0.5; h=5;

TT=[];

XX=[];

for i=1:h:200

k1=f2(x,u);

k2=f2(x+h*k1/2,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 on

plot(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)/(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；