控制系统仿真实验报告

控制系统仿真实验报告
班级：测控1402班
姓名：王玮
学号：07
2018年01月
实验一经典的连续系统仿真建模方法
一实验目的:
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）中的仿真结果有何区别。

三程序代码:
龙格库塔:
%RK4文件
clc
close
H=[,]';u=; 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=;
u=;
Qd=;
A=2;
a1=;
a2=;
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=;h=1; U=;h=1;
U=;h=1;
2 步长h对仿真结果的影响:
U=;h=5; U=;h=20;
U=;h=39 U=;h=50
由以上结果知,仿真步长越大,仿真结果越不稳定。

采用ode45算法程序如下:
function dH=liu(t,H)
k=;
u=;
Qd=;
A=2;
a1=;
a2=;
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],[ ]);
plot(t,H(:,1),['r','+'],t,H(:,2),['g','*'])
u= u=
u=
用ode45与用龙格库塔法仿真结果基本一致。

2编写四阶 Runge_Kutta 公式的计算程序，对线性状态方程（18）式进行仿真：%RK4文件
clc
clear
close
x=[,]';u=; 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=;
a1=;a2=;
A=2;k=;
dx=zeros(2,1);
dx(1)=1/A*(k*u-x(1)/(2*sqrt/a1)+Qd);
dx(2)=1/A*(x(1)/(2*sqrt/a1)-x(2)/(2*sqrt/a2));
1 阀值u对仿真结果的影响:
U=；h=1； U=；h=1；
U=；h=1；
2 步长h对仿真结果的影响:
U=；h=5； U=；h=20；
U=；h=35； U=；h=50；
当步长为50时仿真结果开始不稳定，仿真步长越大,仿真结果越不稳定。

采用ode45算法程序如下:
function dx=liu2(t,x)%ÏßÐÔ
Qd=;
u=;
a1=;a2=;
A=2;k=;
dx=zeros(2,1);
dx(1)=1/A*(k*u-x(1)/(2*sqrt/a1)+Qd);
dx(2)=1/A*(x(1)/(2*sqrt/a1)-x(2)/(2*sqrt/a2));在命令窗口运行以下程序:
>> [t,x]=ode45('liu2',[1 200],[ ]);
>> plot(t,x(:,1),['r','+'],t,x(:,2),['g','*']) U=； u=；
U=
用ode45与用龙格库塔法仿真结果基本一致。

四思考题：
1讨论仿真步长对稳定性和仿真精度的影响
仿真步长越短，稳定性越高，精确度越高。

2你是怎样实现阀位增大和减小的，对于非线性模型和线性模型方法一样吗通过改变u来改变阀的开度，线性系统和非线性系统方法一样。

实验二面向结构图的仿真
第一部分线性系统仿真
一实验目的
1．掌握理解控制系统闭环仿真技术。

2．掌握理解面向结构图的离散相似法的原理和程序结构。

3．掌握 MATLAB 中C2D 函数的用法，掌握双线性变换的原理。

二实验内容
根据上面的各式，编写仿真程序，实现无扰动时给定值阶跃仿真实验
1. 取K P = ，T i = 85 s T = 10s，ΔH2 S =H2set_ percent = 80, ΔQ d = 0，tend = 700，进行仿真实验，绘制响应曲线。

clc
clear all
A=2;
ku=;
H10=;
alpha12 = sqrt(H10); alpha2 = sqrt(H20);
R12=2*sqrt(H10)/alpha12; R2=2*sqrt(H20)/alpha2; H1SpanLo=0;
H2SpanLo=0;
H1SpanHi=;
H2SpanHi=;
Kp=;
Ti=85;
R12*A
R12
ad = 1/(A*R12);
a1 = 1/(A*R12);
a2 = 1/(A*R2);
Kc=Kp/Ti;
Kd = 1/A;
K1 = ku/A;
K2 = 1/(A*R12);
uc(1)=0;ud(1)=0;u1(1)=0;u2(1)=0;
xc(1)=0;xd(1)=0;x1(1)=0;x2(1)=0;
yd(1)=0;yc(1)=0;y1(1)=0;y2(1)=0;
nCounter = 70;
T=10;
k=1;
deltaQd=0;
H20_percent=(H20-H2SpanLo)/(H2SpanHi-H2SpanLo)*100; H2=80;
tend = nCounter*T;
for t=T:T:tend
k=k+1;
uc(k)= (H2 - (y2(k-1)+H20-H2SpanLo)/(H2SpanHi-H2SpanLo)*100)/100;
ud(k)=deltaQd;
u1(k)=yc(k-1);
u2(k)=y1(k-1);
xc(k) = xc(k-1) + Kc*T*uc(k-1);
yc(k)=xc(k)+bc*Kc*uc(k);
xd(k) = exp(-ad*T)*xd(k-1) + Kd/ad*(1-exp(-ad*T))*ud(k);yd(k)=xd(k); x1(k) = exp(-a1*T)*x1(k-1) + K1/a1*(1-exp(-a1*T))*u1(k);y1(k)=x1(k); x2(k) = exp(-a2*T)*x2(k-1) + K2/a2*(1-exp(-a2*T))*u2(k);y2(k)=x2(k); end
Hlevel(:,1)=(y1+H10-H1SpanLo)/(H1SpanHi-H1SpanLo)*100;
Hlevel(:,2)=(y2+H20-H2SpanLo)/(H2SpanHi-H2SpanLo)*100;
yc=(yc+*100;
y2sp=H2*ones(size(y1'));
yv=yc;
textPositionH1=max(Hlevel(:,1));
textPositionH2=max(Hlevel(:,2));
H2Steady=Hlevel(size(Hlevel(:,1),1),1)*ones(size(y1')); xmax=max(0:T:tend);
xmin=0;
ymax=110;
ymin=50;
scrsz = get(0,'ScreenSize');
gca=figure('Position',[5 10 scrsz(3)-10 scrsz(4)-90]);
%gca=figure('Position',[5 10 scrsz(3)/2 scrsz(4)/])
set(gca,'Color','w');
plot(0:T:tend,Hlevel(:,1),'r','LineWidth',2)
hold on
plot(0:T:tend,Hlevel(:,2),'b','LineWidth',2)
hold on
plot(0:T:tend,yv,'k','LineWidth',2)
hold on
plot(0:T:tend,y2sp,'g','LineWidth',2)
hold on
plot(0:T:tend,H2Steady,'y','LineWidth',2)
line([tend/2 tend/2+27],[(ymax-ymin)/2+ymin-(ymax-ymin)/10
(ymax-ymin)/2+ymin-(ymax-ymin)/10],'Color','r','LineWidth',6)
text(tend/2+27,(ymax-ymin)/2+ymin-(ymax-ymin)/10,' 第一个水箱的液位
H1','FontSize',16)
line([tend/2 tend/2+27],[(ymax-ymin)/2+ymin-(ymax-ymin)/6
(ymax-ymin)/2+ymin-(ymax-ymin)/6],'Color','b','LineWidth',6)
text(tend/2+27,(ymax-ymin)/2+ymin-(ymax-ymin)/6,' 第二个水箱的液位
H2','FontSize',16)
line([tend/2 tend/2+27],[(ymax-ymin)/2+ymin-(ymax-ymin)/
(ymax-ymin)/2+ymin-(ymax-ymin)/],'Color','g','LineWidth',6)
text(tend/2+27,(ymax-ymin)/2+ymin-(ymax-ymin)/,' 第二个水箱的液位给定值
','FontSize',16)
line([tend/2 tend/2+27],[(ymax-ymin)/2+ymin-(ymax-ymin)/
(ymax-ymin)/2+ymin-(ymax-ymin)/],'Color','k','LineWidth',6)
text(tend/2+27,(ymax-ymin)/2+ymin-(ymax-ymin)/,'阀位变化情况
','FontSize',16)
axis([xmin xmax ymin ymax]);
text(tend/5,ymax+,' 实验二不考虑阀位饱和特性时的控制效果','FontSize',22) grid
2. 用 MATLAB 求出从输入到输出的传递函数，并将其用c2d 函数，利用双线性变换法转换为离散模型，再用dstep()函数求离散模型的阶跃响应，阶跃幅值为3。

clc
clear all
A=2;
ku=;
H10=;
H20=;
alpha12 = sqrt(H10);
alpha2 = sqrt(H20);
R12=2*sqrt(H10)/alpha12;
R2=2*sqrt(H20)/alpha2;
H1SpanLo=0;
H2SpanLo=0;
H1SpanHi=;
H2SpanHi=;
Kp=;
Ti=85;
R12*A
R12
ad = 1/(A*R12); a1 = 1/(A*R12); a2 = 1/(A*R2); Kc=Kp/Ti;
bc=Ti;
Kd = 1/A;
K1 = ku/A;
K2 = 1/(A*R12);
numc=[Kc*bc,Kc];% 用 MATLAB 求出从输入到输出的传递函数，
denc=[1];
num1=[K1];
den1=[1,a1];
num2=[K2];
den2=[1,a2];
gc=tf(numc,denc);
g1=tf(num1,den1);
g2=tf(num2,den2);
Sysq=gc*g1*g2;
SysG=feedback(Sysq,1);
gg=c2d(SysG,10,’tustin’);% 用c2d 函数，利用双线性变换法转
换为离散模型
dstep(3*{1},{1});%用dstep()函数求离散模型的阶跃响应，阶跃幅值为3结果
0.060.08
0.1
0.12
0.14
0.160.18
0.2
0.22
0.24
Step Response
Time (seconds)A m p l i t u d e
三 思考题
在未考虑调节阀饱和特性时，讨论一下两个水箱液位的变化情况，工业上是否允许讨论阀位的变化情况，工业上是否能实现
在一开始阀位大开，H1，H2液位上升迅速，很快就达到预期值。

但显然不能在工业上实现。

阀位有其本身的最大最小的限制，在仿真中出现的超过100%的情况在现实生活中不可能出现，因此这一部分对应的控制效果也是无效的。

第二部分 含有非线性环节的控制系统仿真
一 实验目的
4． 掌握理解控制系统闭环仿真技术。

5． 掌握理解面向结构图的离散相似法的原理和程序结构。

6． 掌握理含有非线性环节的控制系统的仿真方法。

二实验内容
根据上面的各式，编写仿真程序，实现无扰动时给定值阶跃仿真实验
1. 取K P = ，T i = 85 s T = 10s，ΔH2 S =H2set_ percent = 80, ΔQ d = 0，tend = 700，进行仿真实验，绘制响应曲线。

clc
clear all
A=2;
ku=;
H10=;
H20=;
alpha12 = sqrt(H10);
alpha2 = sqrt(H20);
R12=2*sqrt(H10)/alpha12;
R2=2*sqrt(H20)/alpha2;
H1SpanLo=0;
H2SpanLo=0;
H1SpanHi=;
H2SpanHi=;
Ti=*100;
%Kp=;
%Ti=999999;
ad = 1/(A*R12);
a1 = 1/(A*R12);
a2 = 1/(A*R2);
Kc=Kp/Ti;
bc=Ti;
Kd = 1/A;
K1 = ku/A;
K2 = 1/(A*R12);
uc(1)=0;uv(1)=0;ud(1)=0;u1(1)=0;u2(1)=0; xc(1)=0;xv(1)=0;xd(1)=0;x1(1)=0;x2(1)=0; yc(1)=0;yv(1)=0;yd(1)=0;y1(1)=0;y2(1)=0;
nCounter = 70;
k=1;
deltaQd=0;
c=;
H20_percent=(H20-H2SpanLo)/(H2SpanHi-H2SpanLo)*100;
H2set_percent=80;
tend = nCounter*T;
for t=T:T:tend
k=k+1;
uc(k)= (H2set_percent - (y2(k-1)+H20-H2SpanLo)/(H2SpanHi-H2SpanLo)*100)/100; uv(k)=yc(k-1);
ud(k)=deltaQd;
if uv(k)>c
yv(k)=c;
end
if uv(k)<-c
yv(k)=0;
end
if uv(k)<=c & uv(k)>=-c
yv(k)=uv(k);
end
u1(k)=yv(k);
u2(k)=y1(k-1);
xc(k) = xc(k-1) + Kc*T*uc(k-1);
yc(k)=xc(k)+bc*Kc*uc(k);
xd(k) = exp(-ad*T)*xd(k-1) + Kd/ad*(1-exp(-ad*T))*ud(k);yd(k)=xd(k); x1(k) = exp(-a1*T)*x1(k-1) + K1/a1*(1-exp(-a1*T))*u1(k);y1(k)=x1(k); x2(k) = exp(-a2*T)*x2(k-1) + K2/a2*(1-exp(-a2*T))*u2(k);y2(k)=x2(k); end
Hlevel(:,1)=(y1+H10-H1SpanLo)/(H1SpanHi-H1SpanLo)*100;
Hlevel(:,2)=(y2+H20-H2SpanLo)/(H2SpanHi-H2SpanLo)*100;
yv=(yv+*100;
y2sp=H2set_percent*ones(size(y1'));
textPositionH1=max(Hlevel(:,1));
textPositionH2=max(Hlevel(:,2));
H2Steady=Hlevel(size(Hlevel(:,1),1),1)*ones(size(y1'));
xmax=max(0:T:tend);
xmin=0;
ymax=110;
ymin=50;
scrsz = get(0,'ScreenSize');
gca=figure('Position',[5 10 scrsz(3)-10 scrsz(4)-90]) %gca=figure('Position',[5 10 scrsz(3)/2 scrsz(4)/])
set(gca,'Color','w');
plot(0:T:tend,Hlevel(:,1),'r','LineWidth',2)
hold on
plot(0:T:tend,Hlevel(:,2),'b','LineWidth',2)
hold on
plot(0:T:tend,yv,'k','LineWidth',2)
hold on
plot(0:T:tend,y2sp,'g','LineWidth',2)
hold on
plot(0:T:tend,H2Steady,'y','LineWidth',2)
line([tend/2 tend/2+27],[(ymax-ymin)/2+ymin-(ymax-ymin)/10
(ymax-ymin)/2+ymin-(ymax-ymin)/10],'Color','r','LineWidth',6)
text(tend/2+27,(ymax-ymin)/2+ymin-(ymax-ymin)/10,' 第一个水箱的液位
H1','FontSize',16)
line([tend/2 tend/2+27],[(ymax-ymin)/2+ymin-(ymax-ymin)/6
(ymax-ymin)/2+ymin-(ymax-ymin)/6],'Color','b','LineWidth',6)
text(tend/2+27,(ymax-ymin)/2+ymin-(ymax-ymin)/6,' 第二个水箱的液位
H2','FontSize',16)
line([tend/2 tend/2+27],[(ymax-ymin)/2+ymin-(ymax-ymin)/
(ymax-ymin)/2+ymin-(ymax-ymin)/],'Color','g','LineWidth',6)
text(tend/2+27,(ymax-ymin)/2+ymin-(ymax-ymin)/,' 第二个水箱的液位给定值','FontSize',16)
line([tend/2 tend/2+27],[(ymax-ymin)/2+ymin-(ymax-ymin)/
(ymax-ymin)/2+ymin-(ymax-ymin)/],'Color','k','LineWidth',6)
text(tend/2+27,(ymax-ymin)/2+ymin-(ymax-ymin)/,'阀位变化情况',16)
axis([xmin xmax ymin ymax]);
text(tend/5,ymax+,'实验三考虑阀位饱和特性时的控制效果','FontSize',22) grid
三思考题
与实验三相比，考虑调节阀饱和特性前后，响应有何不同
H1 H2的液位在考虑饱和特性之后，响应曲线比不考虑的时候略微平缓一些。

实验三采样系统的仿真
一实验目的
1．掌握理解数字控制系统的仿真技术。

2．掌握理解增量式 PID 数字控制器的实现方法。

二实验内容
根据上面的各式，编写仿真程序。

取K P = ， T i = 30 s, T d = ，T = 10s，H2set _ percent = 80，tend = 700 ,进行仿真实验，绘制仿真曲线。

clc
clear all
A=2;
ku=;
H10=;
H20=;
alpha12 = sqrt(H10); alpha2 = sqrt(H20);
R12=2*sqrt(H10)/alpha12; R2=2*sqrt(H20)/alpha2; H1SpanLo=0;
H2SpanLo=0;
H1SpanHi=;
H2SpanHi=;
Kp=;
Ti=30;
Td=10;
ad = 1/(A*R12);
a1 = 1/(A*R12);
a2 = 1/(A*R2);
Kc=Kp/Ti;
bc=Ti;
Kd = 1/A;
K1 = ku/A;
K2 = 1/(A*R12);
beta1=1/(a1*a2);
beta3=beta1*a1/(a2-a1); beta2=-beta1-beta3;
u(1)=0;u(2)=0;u(3)=0; y(1)=0;y(2)=0;y(3)=0; nCounter = 70;
T=10;
k=2;
deltuU=0;
e(1)=0;e(2)=0;e(3)=0; uc(1)=0;
H20_percent=(H20-H2SpanLo)/(H2SpanHi-H2SpanLo)*100;
H2set_percent=80;
tend = nCounter*T;
for t=2*T:T:tend
k=k+1;
e(3)=(H2set_percent - (y(k-1)+H20-H2SpanLo)/(H2SpanHi-H2SpanLo)*100)/100; deltaU=Kp*(e(3)-e(2))+Kp*T/Ti*e(3)+Kp*Td/T*[e(3)-2*e(2)+e(1)];
e(1)=e(2);
e(2)=e(3);
u(k)=u(k-1)+deltaU;
y(k)=(exp(-a1*T)+exp(-a2*T))*y(k-1) -
exp(-(a1+a2)*T)*y(k-2)+K1*K2*(beta1+beta2+beta3)*u(k)- ...
K1*K2*(beta1*(exp(-a1*T)+exp(-a2*T))+beta2*(1+exp(-a2*T)) ...
+beta3*(1+exp(-a1*T)))*u(k-1)+K1*K2*(beta1*exp(-(a1+a2)*T)+beta2*exp(-a2*T)+ .. .
beta3*exp(-a1*T))*u(k-2) ;
y(k-2)=y(k-1);
y(k-1)=y(k);
u(k-2)=u(k-1);
u(k-1)=u(k);
end
y=(y+H20-H2SpanLo)/(H2SpanHi-H2SpanLo)*100;
y2sp=H2set_percent*ones(size(y'));
u=(u+*100;
Hlevel(:,1)=y;
Hlevel(:,2)=y;
y1=y;
y2=y;
textPositionH1=max(Hlevel(:,1));
textPositionH2=max(Hlevel(:,2));
H2Steady=Hlevel(size(Hlevel(:,1),1),1)*ones(size(y1')); xmax=max(0:T:tend);
xmin=0;
ymax=90;
ymin=50;
scrsz = get(0,'ScreenSize');
gca=figure('Position',[5 10 scrsz(3)-10 scrsz(4)-90])
%gca=figure('Position',[5 10 scrsz(3)/2 scrsz(4)/])
set(gca,'Color','w');
plot(0:T:tend,Hlevel(:,1),'b','LineWidth',2)
hold on
plot(0:T:tend,y2sp,'k','LineWidth',2)
line([tend/2 tend/2+27],[(ymax-ymin)/2+ymin-(ymax-ymin)/6
(ymax-ymin)/2+ymin-(ymax-ymin)/6],'Color','b','LineWidth',6)
text(tend/2+27,(ymax-ymin)/2+ymin-(ymax-ymin)/6,' 第二个水箱的液位
H2','FontSize',16)
line([tend/2 tend/2+27],[(ymax-ymin)/2+ymin-(ymax-ymin)/
(ymax-ymin)/2+ymin-(ymax-ymin)/],'Color','k','LineWidth',6)
text(tend/2+27,(ymax-ymin)/2+ymin-(ymax-ymin)/,' 第二个水箱的液位给定值','FontSize',16)
axis([xmin xmax ymin ymax]);
text(tend/5,ymax+,' 实验三 PID数字控制器控制效果','FontSize',22) grid
结果： 0
100200300400500600700
5055
6065
70
75
80
85
90
第2个水箱的液位H2
第2个水箱的液位给定值
实验三 PID 数字控制器控制效果
三 实验内容
（1） 针对实验二的第一部分内容，利用Simulink 建立该系统，取K P = ， T i = 85 s , ΔH2 s =0， ΔQd = ，进行仿真实验，观察响应曲线。

0100200300400500600700
5052
54
56
58
60
62
64
（2） 针对实验二的第二部分内容，利用Simulink 建立该系统，取K P = ，T i = 85 s , T = 10s ，ΔH2 s =H 2set _ percent = 80 ， ΔQd =0 ，tend = 700，进行仿真实验，观察
响应曲线。

四思考题
讨论增量式PID算法优点
1、算式中不需要累加。

控制增量Δu(k)的确定仅与最近3次的采样值有关，容易通过加权处理获得比较好的控制效果；
2、计算机每次只输出控制增量，即对应执行机构位置的变化量，故机器发生故障时影响范围小、不会严重影响生产过程；
3、手动—自动切换时冲击小。

当控制从手动向自动切换时，可以作到无扰动切换。