控制系统仿真实验报告0717013819.docx
《控制系统仿真实验报告0717013819.docx》由会员分享，可在线阅读，更多相关《控制系统仿真实验报告0717013819.docx（29页珍藏版）》请在冰点文库上搜索。
控制系统仿真实验报告0717013819
控制系统仿真实验报告
班级：
测控1402班
姓名：
王玮
学号：
1405040207
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=[1.2,1.4]';u=0.55;h=1;
TT=[];
XX=[];
fori=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=[TTi];
XX=[XXH];
end;
holdon
plot（TT,XX（1,:
）,'--',TT,XX（2,:
））;
xlabel（'time'）
ylabel（'H'）
gtext（'H1'）
gtext（'H2'）
holdon
水箱模型:
functiondH=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=39U=0.5;h=50
由以上结果知,仿真步长越大,仿真结果越不稳定。
采用ode45算法程序如下:
functiondH=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*u+Qd-a1*sqrt（H
（1）））;
dH
（2）=1/A*（a1*sqrt（H
（1））-a2*sqrt（H
（2）））;
在命令窗口运行以下程序:
[t,H]=ode45（'liu',[1200],[1.21.1]）;
plot（t,H（:
1）,['r','+'],t,H（:
2）,['g','*']）
u=0.45u=0.5
u=0.55
用
ode45
与用龙格库塔法仿真结果基本一致。
2编写四阶
Runge_Kutta
公式的计算程序，对线性状态方程（
18）式进行仿真：
%RK4文件
clc
clear
close
x=[1.2,1.4]';u=0.5;h=5;
TT=[];
XX=[];
fori=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=[TTi];
XX=[XXx];
end;
holdon
plot（TT,XX（1,:
）,'--',TT,XX（2,:
））;
xlabel（'time'）
ylabel（'x'）
gtext（'x1'）
gtext（'x2'）
holdon
线性函数：
functiondx=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；
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）
%?
?
D?
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',[1200],[1.21.4]）;
>>plot（t,x（:
1）,['r','+'],t,x（:
2）,['g','*']）
U=0.45；u=0.5；
U=0.55
用ode45与用龙格库塔法仿真结果基本一致。
四思考题：
1讨论仿真步长对稳定性和仿真精度的影响仿真步长越短，稳定性越高，精确度越高。
2你是怎样实现阀位增大和减小的，对于非线性模型和线性模型方法一样吗？
通过改变u来改变阀的开度，线性系统和非线性系统方法一样。
实验二面向结构图的仿真
第一部分线性系统仿真
一实验目的
1．掌握理解控制系统闭环仿真技术。
2．掌握理解面向结构图的离散相似法的原理和程序结构。
3．掌握MATLAB中C2D函数的用法，掌握双线性变换的原理。
二实验内容
根据上面的各式，编写仿真程序，实现无扰动时给定值阶跃仿真实验
1.取KP=1.78，Ti=85sT=10s，H2S=H2set_percent=80,Qd=0，
tend=700，进行仿真实验，绘制响应曲线。
clc
clearall
A=2;
ku=0.1/0.5;
H10=1.5;
H20=1.4;
alpha12=0.25/sqrt（H10）;
alpha2=0.25/sqrt（H20）;
R12=2*sqrt（H10）/alpha12;
R2=2*sqrt（H20）/alpha2;
H1SpanLo=0;
H2SpanLo=0;
H1SpanHi=2.52;
H2SpanHi=2.52;
Kp=1.78;
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）;
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;
fort=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+0.5）*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',[510scrsz（3）-10scrsz（4）-90]）;
%gca=figure（'Position',[510scrsz（3）/2scrsz（4）/1.5]）
set（gca,'Color','w'）;
plot（0:
T:
tend,Hlevel（:
1）,'r','LineWidth',2）
holdon
plot（0:
T:
tend,Hlevel（:
2）,'b','LineWidth',2）
holdon
plot（0:
T:
tend,yv,'k','LineWidth',2）
holdon
plot（0:
T:
tend,y2sp,'g','LineWidth',2）
holdon
plot（0:
T:
tend,H2Steady,'y','LineWidth',2）
line（[tend/2tend/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/2tend/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/2tend/2+27],[（ymax-ymin）/2+ymin-（ymax-ymin）/4.2
（ymax-ymin）/2+ymin-（ymax-ymin）/4.2],'Color','g','LineWidth',6）
text（tend/2+27,（ymax-ymin）/2+ymin-（ymax-ymin）/4.2,'第二个水箱的液位给定值
','FontSize',16）
line（[tend/2tend/2+27],[（ymax-ymin）/2+ymin-（ymax-ymin）/3.2
（ymax-ymin）/2+ymin-（ymax-ymin）/3.2],'Color','k','LineWidth',6）
text（tend/2+27,（ymax-ymin）/2+ymin-（ymax-ymin）/3.2,'阀位变化情况
','FontSize',16）
axis（[xminxmaxyminymax]）;
text（tend/5,ymax+1.5,'实验二不考虑阀位饱和特性时的控制效果
','FontSize',22）grid
110
100
90
80
70
60
实验二不考虑阀位饱和特性时的控制效果
第1个水箱的液位H1
第2个水箱的液位H2
第2个水箱的液位给定值
阀位变化情况
50
0100200300400500600700
2.用MATLAB求出从输入到输出的传递函数，并将其用c2d函数，利用双线性变换法转
换为离散模型，再用dstep（）函数求离散模型的阶跃响应，阶跃幅值为3。
clc
clearall
A=2;
ku=0.1/0.5;
H10=1.5;
H20=1.4;
alpha12=0.25/sqrt（H10）;
alpha2=0.25/sqrt（H20）;
R12=2*sqrt（H10）/alpha12;
R2=2*sqrt（H20）/alpha2;
H1SpanLo=0;
H2SpanLo=0;
H1SpanHi=2.52;
H2SpanHi=2.52;
Kp=1.78;
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*gg.num{1},gg.den{1}）;%用dstep（）函数求离散模型的阶跃响应，阶跃幅值为3
结果
StepResponse
0.24
0.22
0.2
0.18
e
d0.16
u
til
p
m0.14
A
0.12
0.1
0.08
0.06
0510152025
Time（seconds）
三思考题
在未考虑调节阀饱和特性时，讨论一下两个水箱液位的变化情况，工业上是否允许？
讨论阀位的变化情况，工业上是否能实现？
在一开始阀位大开，H1，H2液位上升迅速，很快就达到预期值。
但显然不能在工业上实
现。
阀位有其本身的最大最小的限制，在仿真中出现的超过100%的情况在现实生活中不可能出现，因此这一部分对应的控制效果也是无效的。
第二部分含有非线性环节的控制系统仿真
一实验目的
4．掌握理解控制系统闭环仿真技术。
5．掌握理解面向结构图的离散相似法的原理和程序结构。
6．掌握理含有非线性环节的控制系统的仿真方法。
二实验内容
根据上面的各式，编写仿真程序，实现无扰动时给定值阶跃仿真实验
1.取
KP
=1.78，
Ti
=85
sT
=10
s
，
H2S
=2
_
percent
=80,
Qd
=0，
Hset
tend=700，进行仿真实验，绘制响应曲线。
clc
clearall
A=2;
ku=0.1/0.5;
H10=1.5;
H20=1.4;
alpha12=0.25/sqrt（H10）;
alpha2=0.25/sqrt（H20）;
R12=2*sqrt（H10）/alpha12;
R2=2*sqrt（H20）/alpha2;
H1SpanLo=0;
H2SpanLo=0;
H1SpanHi=2.52;
H2SpanHi=2.52;
Kp=3.91/2.2;;
Ti=0.85*100;
%Kp=3.21;
%Ti=99999999999999;
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;
T=10;
k=1;
deltaQd=0;
c=0.5;
H20_percent=（H20-H2SpanLo）/（H2SpanHi-H2SpanLo）*100;
H2set_percent=80;
tend=nCounter*T;
fort=T:
T:
tend
k=k+1;
uc（k）=（H2set_percent
uv（k）=yc（k-1）;
ud（k）=deltaQd;
ifuv（k）>c
yv（k）=c;
-（y2（k-1）+H20-H2SpanLo）/（H2SpanHi-H2SpanLo）*100）/100;
end
ifuv（k）<-c
yv（k）=0;
end
ifuv（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+0.5）*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'
[510scrsz（3）-10scrsz（4）-90]）
%gca=figure（'Position',[510scrsz（3）/2scrsz（4）/1.5]）
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/2tend/2+27],[（ymax-ymin）/2+ymin-（ymax-ymin）/10
（ymax-ymin）/2+ymin-（ymax-ymin）/10],
'Color'
'r'
'LineWidth',6）
text（
升级会员 