控制系统仿真
实验报告
姓名：大葱哥 学号：
班级：测控1202
2015.10.28
实验一 经典的连续系统仿真建模方法
一 实验目的
1． 了解和掌握利用仿真技术对控制系统进行分析的原理和步骤。 2． 掌握机理分析建模方法。
3． 深入理解一阶常微分方程组数值积分解法的原理和程序结构，学习用 Matlab 编写 数值积分法仿真程序。
4． 掌握和理解四阶 Runge-Kutta 法，加深理解仿真步长与算法稳定性的关系。
二 实验内容
1. 编写四阶 Runge_Kutta 公式的计算程序， 对非线性模型 （3） 式进行仿真。 仿真程序: 主程序：
clc clear all u=zeros(2,1); u(1)=0.5%稳态
%u(1)=0.55;%阀位增大%10 %u(1)=0.45;%阀位减小%10 u(2)=0.15;
h=zeros(1,2); h(1,1)=1.5; h(1,2)=1.4; hStep = 10; Hlevel = h; nCounter = 25;
for t=0:hStep:(nCounter-1)*hStep
h = Z06_SystemSimulation_Lab01_Nonlinear_RK4(hStep,t,h,u); Hlevel=[Hlevel;h]; end figure(1)
plot([0:hStep:nCounter*hStep]',Hlevel) grid
四阶龙格库塔算法：
function RK4_Result =
Z06_SystemSimulation_Lab01_Nonlinear_RK4(h,t0,x0,u0) K1=Z06_SystemSimulation_Lab01_Linear_dxCompute(t0,x0,u0);
K2=Z06_SystemSimulation_Lab01_Linear_dxCompute(t0+h/2,x0+h*K1/2,u0); K3=Z06_SystemSimulation_Lab01_Linear_dxCompute(t0+h/2,x0+h*K2/2,u0); K4=Z06_SystemSimulation_Lab01_Linear_dxCompute(t0+h,x0+h*K3,u0); RK4_Result = x0 + h*(K1+2*K2+2*K3+K4)/6;
（1） 将阀位 u 增大 10％和减小 10％，观察响应曲线的形状;
稳态仿真曲线：
增大10%：
减小10%：
（2） 研究仿真步长对稳定性的影响，仿真步长取多大时 RK4 算法变得不稳定？
仿真步长的选取范围可以通过理论计算获知，这里我们通过仿真得到当步长大于63时，曲线发散，系统不再稳定，发散曲线如下图：
（3） 利用 MATLAB 中的 ode45()函数进行求解， 比较与 （1） 中的仿真结果有何区别。 2. 编写四阶 Runge_Kutta 公式的计算程序，对线性状态方程（18）式进行 仿真 仿真程序
clc clear all close all u=zeros(2,1); u(1)=0.0;%稳态
% u(1)=0.05;%阀位增大%10 % u(1)=-0.05;%阀位减小%10 u(2)=0;
h=zeros(1,2); h(1,1)=0; h(1,2)=0; hStep = 10; Hlevel = h; nCounter = 25;
for t=0:hStep:(nCounter-1)*hStep
h = Z06_SystemSimulation_Lab01_Linear_RK4(hStep,t,h,u); Hlevel=[Hlevel;h]; end
Hlevel(:,1) = Hlevel(:,1)+1.5; Hlevel(:,2) = Hlevel(:,2)+1.4; figure(3)
plot([0:hStep:nCounter*hStep]',Hlevel) grid hold on
（1） 将阀位增大 10％和减小 10％，观察响应曲线的形状;
稳态仿真曲线：
阀位增大%10：
阀位减小%10：
（2） 研究仿真步长对稳定性的影响，仿真步长取多大时 RK4 算法变得不稳定？
当步长大于62时，算法变得不再稳定
（3） 阀位增大 10％和减小 10％，利用 MATLAB 中的 ode45()函数进行求解阶跃响应，
比较与（1）中的仿真结果有何区别。
三 思考题
1． 讨论仿真步长对稳定性和仿真精度的影响。
仿真步长越长，系统越倾向于不稳定，仿真步长越短，系统越稳定，但是需要的仿真时间越长。
2． 你是怎样实现阀位增大和减小10%的？对于非线性模型和线性模型方法一样吗？ 通过改变u(1)的初值可以实现对阀位的增大和减小，由于线性模型法是对变量的变化量进行计算的，而非线性模型法是直接对变量进行计算，所以两种方法实现对阀位的增大和减小方面存在一些区别，不完全一致。
四 实验总结
通过本次实验我学会了Matlab的基本使用方法，对该软件有了一定程度的认识，同时对仿真步长对系统的稳定性和仿真精度的关系有了一定的理解，对一阶常微分方程组数值积分解法的原理和程序结构有了进一步的掌握，深入了解了Runge-Kutta 法，对MATLAB软件的仿真有了一定程度的认识，为以后的学习奠定了基础。
实验二 面向结构图的仿真 第一部分 线性系统仿真
一 实验目的
1． 掌握理解控制系统闭环仿真技术。
2． 掌握理解面向结构图的离散相似法的原理和程序结构。
3．掌握 MATLAB 中 C2D 函数的用法，掌握双线性变换的原理。
二 实验内容
1.编写仿真程序，实现无扰动时给定值阶跃仿真实验
仿真程序：
%面向结构图的离散相似法,线性模型 clc clear all
%%%%%%%%%%%%%%%%%%%%%定义参数%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% A=2; %水箱横截面积 ku=0.1/0.5; %阀门流量系数 H10=1.5; %水箱1的平衡高度 H20=1.4; %水箱2的平衡高度
alpha12 = 0.25/sqrt(H10); %水箱1流向水箱2流量系数 alpha2 = 0.25/sqrt(H20); %水箱2流出水流量系数 R12=2*sqrt(H10)/alpha12; %线性化数学模型中的参数 R2=2*sqrt(H20)/alpha2; %线性化数学模型中的参数 H1SpanLo=0; %水箱1量程下限 H2SpanLo=0; %水箱2量程下限 H1SpanHi=2.52; %水箱1量程上限 H2SpanHi=2.52; %水箱2量程上限 Kp=1.78; %PI控制器比例系数
Ti=85; %PI控制器积分时间常数 R12*A R12
Kc=Kp/Ti; %式(4)中参数 bc=Ti; %式(4)中参数
Kd = 1/A; %式(5)中参数 ad = 1/(A*R12); %式(5)中参数
a1 = 1/(A*R12); %式(6)中参数 K1 = ku/A; %式(6)中参数
a2 = 1/(A*R2); %式(7)中参数
K2 = 1/(A*R12); %式(7)中参数
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%定义初值%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% uc(1)=0; %Gc环节(8)的输入初值 xc(1)=0; %Gc环节(8)的状态初值 yc(1)=0; %Gc环节(8)的输出初值
ud(1)=0; %Gd环节(9)的输入初值 xd(1)=0; %Gd环节(9)的状态初值 yd(1)=0; %Gd环节(9)的输出初值
u1(1)=0; %G1环节(10)的输入初值 x1(1)=0; %G1环节(10)的状态初值 y1(1)=0; %G1环节(10)的输出初值
u2(1)=0; %G2环节(11)的输入初值 x2(1)=0; %G2环节(11)的状态初值
y2(1)=0; %G2环节(11)的输出初值
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 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; %系统(8)的输入百分比 ud(k)=deltaQd; %系统(9)的输入 u1(k)=yc(k-1); %系统(10)的输入 u2(k)=y1(k-1); %系统(11)的输入
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',[5 10 scrsz(3)-10 scrsz(4)-90]); %gca=figure('Position',[5 10 scrsz(3)/2 scrsz(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/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,' 第1个水箱的液位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,' 第2个水箱的液位H2','FontSize',16)
line([tend/2 tend/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,' 第2个水箱的液位给定值','FontSize',16)
line([tend/2 tend/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([xmin xmax ymin ymax]);
text(tend/5,ymax+1.5,' 实验二 不考虑阀位饱和特性时的控制效果','FontSize',22) grid
%hold off
仿真实验曲线：
由控制的效果可以看出，阀位的变化出现了超过100%的情况，这是控制系统中不允许出现的，说明控制功能需要改进。
三 思考题
在未考虑调节阀饱和特性时， 讨论一下两个水箱液位的变化情况， 工业上是否允许？
讨论阀位的变化情况，工业上是否能实现？
答：在未考虑调节阀饱和特性时，两个水箱的液位在阀位一开始大开时，H1，H2液位上升迅速，很快就达到预期值。但是其控制不能在工业上实现，没有实际可行度。除此之外，阀位有其本身的最大最小的限制，在仿真中出现的超过100%的情况在现实生活中不可能出现，因此这一部分对应的控制效果在工业上不能实现。
四 实验总结
通过本次实验，我对控制系统闭环仿真技术有了进一步的认识，掌握了面向结构图的离散相似法的原理和程序结构，了解了双线性变换的原理，收获颇丰。
第二部分 含有非线性环节的控制系统仿真
一 实验目的
4． 掌握理解控制系统闭环仿真技术。
5． 掌握理解面向结构图的离散相似法的原理和程序结构。 6． 掌握理含有非线性环节的控制系统的仿真方法。
二 实验内容
根据上面的各式，编写仿真程序，实现无扰动时给定值阶跃仿真实验
仿真程序：
%面向结构图的离散相似法,非线性模型 clc
clear all 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=82; 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+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',[5 10 scrsz(3)-10 scrsz(4)-90]) %gca=figure('Position',[5 10 scrsz(3)/2 scrsz(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/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,' 第1个水箱的液位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,' 第2个水箱的液位H2','FontSize',16)
line([tend/2 tend/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,' 第2个水箱的液位给定值','FontSize',16)
line([tend/2 tend/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([xmin xmax ymin ymax]);
text(tend/5,ymax+1.5,' 实验三 考虑阀位饱和特性时的控制效果','FontSize',22) grid
%hold off 仿真实验曲线：
三 思考题
与实验三相比，考虑调节阀饱和特性前后，响应有何不同？
答：与实验三相比，调节饱和特性阀前后调进行对比，H1 、H2的液位在考虑饱和特性之后，响应曲线比不考虑的时候略微平缓一些，阀位的调节没有出现超过100%的情况。
四 实验总结
通过本次实验，我认识了面向结构图的离散相似法的原理和程序结构，掌握了理含有非线性环节的控制系统的仿真方法。
实验三 采样系统的仿真
一 实验目的
1． 掌握理解数字控制系统的仿真技术。
2． 掌握理解增量式 PID 数字控制器的实现方法。
二 实验内容 1
根据上面的各式，编写仿真程序。
仿真程序： clc
clear all close all; 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.89; 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(k)=(H2set_percent - (y(k-1)+H20-H2SpanLo)/(H2SpanHi-H2SpanLo)*100)/100; deltaU=Kp*(e(k)-e(k-1))+Kp*T/Ti*e(k)+Kp*Td/T*[e(k)-2*e(k-1)+e(k-2)]; e(k-2)=e(k-1); e(k-1)=e(k);
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+0.5)*100;
% 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)/1.5]) set(gca,'Color','w');
% plot(0:T:tend,Hlevel(:,1),'g','LineWidth',2) plot(0:T:tend,y,'b') figure
plot(0:T:tend,e,'b')
% 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,' 第2个水箱的液位H2','FontSize',16)
% line([tend/2 tend/2+27],[(ymax-ymin)/2+ymin-(ymax-ymin)/4.2 (ymax-ymin)/2+ymin-(ymax-ymin)/4.2],'Color','k','LineWidth',6)
% text(tend/2+27,(ymax-ymin)/2+ymin-(ymax-ymin)/4.2,' 第2个水箱的液位给定值','FontSize',16)
% axis([xmin xmax ymin ymax]); % text(tend/5,ymax+1.5,' 实验三 PID数字控制器控制效果','FontSize',22) % Grid
%hold off
仿真曲线：
被控变量H的变化趋势：
控制量的变化趋势：
三 实验内容 2
搭建系统如图所示：
运行仿真：
结构图：
仿真曲线H2：
仿真曲线H1：
对比实验二中的编程实现结果，仿真曲线基本一致。
四 思考题
讨论增量式PID算法优点
答：①算式中不需要累加。控制增量Δu(k)的确定仅与最近3次的采样值有关，容易通过加权处理获得比较好的控制效果；
②计算机每次只输出控制增量，即对应执行机构位置的变化量，故机器发生故障时影响范围小、不会严重影响生产过程；
③手动—自动切换时冲击小。当控制从手动向自动切换时，可以作到无扰动切换。
五 实验总结
通过本次实验，我对采样控制有了更深一步的认识，基本掌握了理解数字控制系统的仿真技术，学会了增量式 PID 数字控制器的实现方法。最后，感谢老师的悉心指导。
因篇幅问题不能全部显示，请点此查看更多更全内容