控制系统仿真实验报告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(