控制系统仿真实验报告

控制系统仿真实验报告班级：测控1402班姓名：王玮学号：14050402072018年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文件clccloseH=[1.2,1.4]';u=0.55;h=1;TT=[];XX=[];fori=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=[TTi];XX=[XXH];end;holdonplot(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.5u=0.55用ode45与用龙格库塔法仿真结果基本一致。2编写四阶Runge_Kutta公式的计算程序，对线性状态方程（18）式进行仿真：%RK4文件clcclearclosex=[1.2,1.4]';u=0.5;h=5;TT=[];XX=[];fori=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=[TTi];XX=[XXx];end;holdonplot(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算法程序如下:functiondx=liu2(t,x)%ÏßÐÔ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，进行仿真实验，绘制响应曲线。clcclearallA=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/(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:tendk=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);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)*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)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/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([xminx