控制系统仿真实验报告0717013819.docx

1、控制系统仿真实验报告0717013819控制系统仿真实验报告班级：测控 1402 班姓名：王玮学号： 14050402072018 年 01 月实验一经典的连续系统仿真建模方法一 实验目的 :1 了解和掌握利用仿真技术对控制系统进行分析的原理和步骤。2 掌握机理分析建模方法。3深入理解阶常微分方程组数值积分解法的原理和程序结构，学习用Matlab 编写数值积分法仿真程序。4掌握和理解四阶 Runge-Kutta法，加深理解仿真步长与算法稳定性的关系。二 实验内容 :1. 编写四阶 Runge_Kutta 公式的计算程序，对非线性模型（3）式进行仿真。（1） 将阀位 u 增大 10和减小 10，

2、观察响应曲线的形状 ;（2） 研究仿真步长对稳定性的影响，仿真步长取多大时 RK4 算法变得不稳定？（3） 利用 MATLAB 中的 ode45() 函数进行求解，比较与（ 1）中的仿真结果有何区别。2. 编写四阶 Runge_Kutta 公式的计算程序，对线性状态方程（18）式进行仿真（1） 将阀位增大 10和减小 10，观察响应曲线的形状 ;（2） 研究仿真步长对稳定性的影响，仿真步长取多大时 RK4 算法变得不稳定？（ 4） 阀位增大 10和减小 10，利用 MATLAB中的 ode45() 函数进行求解阶跃响应，比较与（ 1）中的仿真结果有何区别。三 程序代码 :龙格库塔 :%RK4文

3、件clccloseH=1.2,1.4;u=0.55; h=1;TT=;XX=;for i=1:h:200k1=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 onplot(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

4、=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=

5、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,1 200,1.2 1.1);plot(t,H(:,1),r,+,t,H(:,2),g,*)u=0.45 u=0.5u=0.55用ode45与用龙格库塔法仿真结果基本一致。2 编写四阶Runge_Kutta公式的计算程序，对线性状态方程（18）式进行仿真：%RK4文件cl

6、cclearclosex=1.2,1.4;u=0.5; h=5;TT=;XX=;for i=1:h:200k1=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 onplot(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

7、;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 算法程序如下 :, 仿真结果越不稳定。fu

8、nctiondx=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,1 200,1.2 1.4); plot(t,x(:,1),r,+,t,x(:,2),g,*)U=0.45 ； u=0.5 ；U=0.55用 ode45 与用龙格库塔法仿真结果基本一致。

9、四 思考题：1 讨论仿真步长对稳定性和仿真精度的影响仿真步长越短，稳定性越高，精确度越高。2 你是怎样实现阀位增大和减小的，对于非线性模型和线性模型方法一样吗？通过改变 u 来改变阀的开度，线性系统和非线性系统方法一样。实验二 面向结构图的仿真第一部分 线性系统仿真一 实验目的1 掌握理解控制系统闭环仿真技术。2 掌握理解面向结构图的离散相似法的原理和程序结构。3 掌握 MATLAB 中 C2D 函数的用法，掌握双线性变换的原理。二 实验内容根据上面的各式，编写仿真程序，实现无扰动时给定值阶跃仿真实验1. 取 K P = 1.78 ， T i = 85 sT = 10 s， H2 S =H2s

10、et _ percent = 80, Q d = 0 ，tend = 700 ，进行仿真实验，绘制响应曲线。clcclear allA=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*AR12ad = 1/(A*R12);a1 = 1/(A*R12);a2 = 1

11、/(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;for t=T:T:tendk=k+1;uc(k)= (H2 - (y2(k-

12、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

13、(k-1) + K2/a2*(1-exp(-a2*T)*u2(k);y2(k)=x2(k);endHlevel(:,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

14、)*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 onplot(0:T:tend,Hlevel(:,2), b , LineWidt

15、h ,2)hold onplot(0:T:tend,yv, k , LineWidth ,2)hold onplot(0:T:tend,y2sp, g , LineWidth ,2)hold onplot(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)/

16、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)/4.2(ymax-ymin)/2+ymin-(ymax-y

17、min)/4.2, Color , g , LineWidth ,6)text(tend/2+27,(ymax-ymin)/2+ymin-(ymax-ymin)/4.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 ,

18、16)axis(xmin xmax ymin ymax);text(tend/5,ymax+1.5, 实验二 不考虑阀位饱和特性时的控制效果 , FontSize ,22) grid11010090807060实验二 不考虑阀位饱和特性时的控制效果第1个水箱的液位 H1第2个水箱的液位 H2第2个水箱的液位给定值阀位变化情况500 100 200 300 400 500 600 7002. 用 MATLAB 求出从输入到输出的传递函数，并将其用 c2d 函数，利用双线性变换法转换为离散模型，再用 dstep() 函数求离散模型的阶跃响应，阶跃幅值为 3。clcclear allA=2;ku=0

19、.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*AR12ad = 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;% 用

20、 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.num1,gg.den1);% 用 dstep() 函数求离散模型的阶跃响应，阶跃幅值为 3结果Step Response0.240.220.20.18ed 0.16uti

21、lpm 0.14A0.120.10.080.060 5 10 15 20 25Time (seconds)三 思考题在未考虑调节阀饱和特性时， 讨论一下两个水箱液位的变化情况， 工业上是否允许？讨论阀位的变化情况，工业上是否能实现？在一开始阀位大开， H1， H2 液位上升迅速，很快就达到预期值。但显然不能在工业上实现。阀位有其本身的最大最小的限制，在仿真中出现的超过 100%的情况在现实生活中不可能出现，因此这一部分对应的控制效果也是无效的。第二部分 含有非线性环节的控制系统仿真一 实验目的4 掌握理解控制系统闭环仿真技术。5 掌握理解面向结构图的离散相似法的原理和程序结构。6 掌握理含有非

22、线性环节的控制系统的仿真方法。二 实验内容根据上面的各式，编写仿真程序，实现无扰动时给定值阶跃仿真实验1. 取K P= 1.78 ，T i= 85sT= 10s，H2 S= 2_percent= 80,Q d= 0 ，H settend = 700 ，进行仿真实验，绘制响应曲线。clcclear allA=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;

23、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;nCounte

24、r = 70;T=10;k=1;deltaQd=0;c=0.5;H20_percent=(H20-H2SpanLo)/(H2SpanHi-H2SpanLo)*100;H2set_percent=80;tend = nCounter*T;for t=T:T:tendk=k+1;uc(k)= (H2set_percentuv(k)=yc(k-1);ud(k)=deltaQd;if uv(k)cyv(k)=c;- (y2(k-1)+H20-H2SpanLo)/(H2SpanHi-H2SpanLo)*100)/100;endif uv(k)-cyv(k)=0;endif uv(k)=-cyv(k)=u

25、v(k);endu1(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);endHleve

26、l(:,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=11

27、0;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)holdonplot(0:T:tend,Hlevel(:,2),b, LineWidth,2)holdonplot(0:T:tend,yv,k, LineWidth,2)holdonplot(0:T:tend,y2sp,g , LineWidth,2)holdonplot(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(