自动控制上机作业

1、自动控制上机大作业班级：学号：姓名：1.1设质量阻尼弹簧系统的微分运动方程为 式中，x(t)为位移输出信号，f(t)为输入的力信号。质量为M=1kg，粘性摩擦系数为B = 5N / m s1，弹簧的弹性系数为K=20N/m。当t=0 时，施加外力f(t)=30N，试问系统何时达到稳定？并画出该机械系统位移、速度随时间变化的曲线以及速度与位移的关系曲线。提示：龙格-库塔法求解微分方程数值解的函数：odel13(),调用方式：T,Y = ODE113(ODEFUN,TSPAN,Y0,OPTIONS)。其中ODEFUN 为用户自定义的系统微分方程的描述，本题中可使用xt4odefile.m 文件定义

2、的函数；TSPAN 表示计算开始和结束的时间；Y0 表示微分方程的初始条件；OPTION 为计算精度的可选参数，由odese()函数设置。odel13()函数只接受一阶微分方程的形式，使用时需要先将高阶方程化为若干个一阶微分方程；绘图函数：plot()，subplot();程序：ft = 30;M=1;B=5;K=20; %系统参数tspan = 0 20; %设置仿真开始和结束时间x0 = 0,0; %系统初始值，零初始条件options = odeset('AbsTol',1e-6;1e-6); %设置仿真计算精度 t,x = ode113('xt4odefile&

3、#39;,tspan,x0,options);%微分方程求解。 a = 1/M*(ft-B*x(:,2)-K*x(:,1); %计算加速度 i = 1; while (abs(a(i)>0.0001|(abs(x(i,2)>0.0001) i = i+1; end %显示计算结果result = sprintf('位移 d=%6.4fn',x(i,1);disp(result);result = sprintf('速度 v=%8.6fn',x(i,2);disp(result);result = sprintf('加速度 a=%9.6fn&#

4、39;,a(i);disp(result);result = sprintf('时间 t=%4.2fn',t(i);disp(result);d = x(:,1); subplot(1,3,1),plot(t,d); %绘制时间-位移曲线xlabel('时间（秒）');ylabel('位移（米）');title('时间-位移曲线');grid; v = x(:,2);subplot(1,3,2),plot(t,v); %绘制时间-速度曲线xlabel('时间（秒）');ylabel('速度（米/秒）'

5、;);title('时间-速度曲线');grid; subplot(1,3,3),plot(d,v); %绘制位移-速度曲线xlabel('位移（米）');ylabel('速度（米/秒）');title('位移-速度曲线');grid;其中xt4odefile.m文件为：function xt = odefileC(t,x);ft = 30;M=1;B=5;K=20;xt = x(2);1/M*(ft-B*x(2)-K*x(1);计算结果：位移 d=1.5000速度 v=-0.000086加速度 a=-0.000084时间 t=4

6、.461.2假设控制系统的传递函数为试求其零点、极点和增益，并进行部分分式展开。提示：传递函数描述：tf(), 调用方式：SYS = TF(NUM,DEN)。求取零点和极点的函数：tf2zp(), 调用方式：Z,P,K = TF2ZP(NUM,DEN)传递函数的部分分式展开：residue(), 调用方式：R,P,K = RESIDUE(B,A)clcsys=tf(2 5 7,1 6 10 6)disp('零点极点分别为Z,P')Z,P,K=tf2zp(2 5 7,1 6 10 6)B=2 5 7;A=1 6 10 6;disp('部分分式展开')R,P,K =

7、 RESIDUE(B,A)计算结果：Transfer function: 2 s2 + 5 s + 7-s3 + 6 s2 + 10 s + 6 零点极点分别为Z,PZ = -1.2500 + 1.3919i -1.2500 - 1.3919iP = -3.7693 -1.1154 + 0.5897i -1.1154 - 0.5897iK = 2部分分式展开R = 2.2417 -0.1208 - 1.0004i -0.1208 + 1.0004iP = -3.7693 -1.1154 + 0.5897i -1.1154 - 0.5897iK = 1.3 考虑由下式表示的高阶系统 ，试求取系统

8、的单位阶跃响应，并计算系统的上升时间、峰值时间、超调量和调整时间（2%误差带）。提示：阶跃响应函数：step(), 调用方式：Y,T = STEP(SYS)。法一clcsys=tf(6.3223 18 12.811,1 6 11.3223 18 12.811 )y,t=step(sys);mp=max(y);tp=spline(y,t,mp)%峰值时间cs=length(t);yss=y(cs)%稳态ct=(mp-yss)/yss%超调量结果：法二clcsys=tf(6.3223 18 12.811,1 6 11.3223 18 12.811)y,t=step(sys);ltiview(sys

9、)1.4clcnum=1,1;dun=1,5,6,0;rlocus(num,dun)当K=20.575时，clcnum=20.575,20.575;dun=1,5,6,0;sys=tf(num,dun)ltiview(sys)2.1clcnum=0.01,0.0001,0.01;dun=0.25,0.01,1,0,0;sys=tf(num,dun)figure(1)bode(sys)figure(2)sys2=feedback(sys,1)%加单位反馈bode(sys2)结果：系统总的波特图2.2 clcnum=20,20,10;dun=1,11,10,0;sys=tf(num,dun)nyq

10、uist(sys)结果如下：2.3 clcnum=2000 2000;den=1 14.5 407 200 0;sys=tf(num,den)nichols(sys)v = -270 -90 -40 40;axis(v)ngrid 结果：2.4clcnum=2000 2000;den=1 14.5 407 200 0;sys=tf(num,den)figure(1)bode(sys)gm pm wg wp=margin(sys)结果如下：2.5clcnum=1;den=0.5 1.5 1 0;sys=tf(num,den)sys2=feedback(sys,1)bode(sys2)gm pm

11、wg wp=margin(sys)结果如下：3.1A程序如下：clcnum=1 2 3;den=1 3 3 1;A,B,C,D=tf2ss(num,den)结果如下：B程序如下：clcZ=-1 -3;P=-2 -4 -6;K=4;A,B,C,D=zp2ss(Z,P,K)结果如下：C程序如下：clcA=0 1;1 -2;B=0;1;C=1 3;D=1;num,den=ss2tf(A,B,C,D);tf(num,den)z,p,k=ss2zp(A,B,C,D);zpk(z,p,k)结果如下：3.2程序如下：clcA1=0 1;-1 -2;B1=1;0;C1=1 3;D1=1;A2=0 1;-1 -

12、3;B2=0;1;C2=1 4;D2=0;num1,den1=ss2tf(A1,B1,C1,D1);sys1=tf(num1,den1)num2,den2=ss2tf(A2,B2,C2,D2);sys2=tf(num2,den2)sysc=series(sys1,sys2)sysb=parallel(sys1,sys2)numc1,denc1=feedback(num1,den1,num2,den2,-1);sysf=tf(numc1,denc1)numc2,denc2=feedback(num1,den1,num2,den2,1);sysf=tf(numc2,denc2)结果如下：3.3（1

13、）程序如下：clcA=0 -2;1 -3;t=0.2;eAt=expm(A*t)结果如下：（2）结果如下：clcA=0 -2;1 -3;B=2;0;C=0 3;D=0;t=0:0.01:0.2;a=size(t);u=zeros(a);x0=1;1;lsim(A,B,C,D,u,t,x0)title('系统的响应')结果如下：3.4程序如下：clcA=-3 1;1 -3;B=1 1;1 1;C=1 1;1 -1;D=0 0;0 0;N=size(A);n=N(1);num,den=ss2tf(A,B,C,D,2);disp('可控矩阵')S=ctrb(A,B)%

14、计算可控矩阵Sf=rank(S)%通过rank命令求可控矩阵if (f=n) %判断系统的能空性 disp('系统是可控的')else disp('系统是不可控的')enddisp(' ')disp('可观测矩阵')V=obsv(A,C)%计算可观测性矩阵Vm=rank(V)%求可观测性矩阵的秩if (f=n) %判断系统的能观测性 disp('系统是可观测的')else disp('系统是不可观测的')end结果如下：3.5程序如下：clcA=0 1;-2 -3;B=0;1;C=2 0;D=0;P

15、_S=-1 -2;%系统的配置极点k=acker(A,B,P_S);%计算系统的反馈增益向量kP_O=-3 -3;%观测器的期望配置极点h=(acker(A',C',P_O)'%计算观测器输出反馈阵A1=A,-B*k;h*C A-B*k-h*C;B1=B;B;C1=C,zeros(1,2);D1=0;sys=ss(A1,B1,C1,D1)%建立复合系统动态模型tf(sys)结果如下：校正设计4.1增量式PID 控制算法被控对象为，PID 控制参数为：kp=8，ki=0.10，kd=10。下面程序是输入为正弦信号，采样信号是1ms，控制器出的曲线以及误差曲线程序：ts=0

16、.001;sys=tf(400,1,50,0);dsys=c2d(sys,ts,'z');num,den=tfdata(dsys,'v'); u_1=0.0;u_2=0.0;y_1=0.0;y_2=0.0;x=0,0,0'error_1=0;error_2=0;for k=1:1:1000time(k)=k*ts; kp=8;ki=0.1;kd=10; %Sine Signal rin(k)=0.5*sin(2*pi*k*ts); du(k)=kp*x(1)+kd*x(2)+ki*x(3); %PID Controlleru(k)=u_1+du(k);%

17、Restricting the output of controller %Linear modelyout(k)=-den(2)*y_1-den(3)*y_2+num(2)*u_1+num(3)*u_2; error(k)=rin(k)-yout(k); %Return of parametersu_2=u_1;u_1=u(k);y_2=y_1;y_1=yout(k); x(1)=error(k)-error_1; %Calculating Px(2)=error(k)-2*error_1+error_2; %Calculating Dx(3)=error(k); %Calculating

18、I error_2=error_1;error_1=error(k);endfigure(1);plot(time,rin,'b',time,yout,'r') ,grid ongtext('rinrightarrow ')gtext('leftarrow yout')title('系统输出曲线')xlabel('time(s)'),ylabel('rin,yout'); figure(2);plot(time,error,'r') ,grid on title(&#

19、39;误差曲线')xlabel('time(s)');ylabel('error');结果如下：4.2积分分离PID 控制算法被控对象为，采样时间为20s，延迟时间为4 个采样时间。解：离散积分分离控制算法表达式：程序如下：clear all;close all;ts=20; %delay plantsys=tf(1,60,1,'inputdelay',80);dsys=c2d(sys,ts,'z');num,den=tfdata(dsys,'v');u_1=0;u_2=0;u_3=0;u_4=0;u_5=

20、0;y_1=0;y_2=0;y_3=0;error_1=0;ei=0;for k=1:1:200time(k)=k*ts;rin(k)=40;kp=0.80;ki=0.005;kd=3.0;%Delay?plant?yout(k)=-den(2)*y_1+num(2)*u_5;%I?separation?error(k)=rin(k)-yout(k);M=2;if M=1 %Using?integration?separation? if abs(error(k)>=30&abs(error(k)<=40 beta=0.3; elseif abs(error(k)>=

21、20&abs(error(k)<=30 beta=0.6; elseif abs(error(k)>=10&abs(error(k)<=20 beta=0.9; else beta=1.0; endelseif M=2 beta=1.0; %Not?using?integration?separation?endei=ei+error(k)*ts; %积分作用u(k)=kp*error(k)+kd*(error(k)-error_1)/ts+beta*ki*ei;if u(k)>=110 u(k)=110;endif u(k)<=-110 u(k)

22、=-110;endu_5=u_4;u_4=u_3;u_3=u_2;u_2=u_1;u_1=u(k);y_3=y_2;y_2=y_1;y_1=yout(k);error_1=error(k);endfigure(1);plot(time,rin,'b',time,yout,'r');title('加积分作用效果')xlabel('time(s)');ylabel('rin,yout');figure(2);plot(time,u,'r');xlabel('time(s)');ylabe

23、l('u');当M=2时，是不加积分作用的效果，结果如下：当M=1时，是加积分作用的效果，结果如下：4.3抗积分饱和PID 控制算法被控对象为，采样时间为1ms，取指令信号rin（k）=30。M=1 和M=2 时仿真结果分别是怎样的？原理：积分饱和现象是在系统存在一个方向的偏差，PID控制器的输出由于积分作用的不断加大而加大，从而导致执行器达到极限位置，如果控制器的输出继续增大，阀门开度却不能再增大，此时计算机输出超出了执行器正常操作的范围而进入了饱和区。一旦系统出现反方向的偏差，控制器的输出逐渐退出饱和区，但若进入的饱和区很长，那么退出饱和区的时间也就很长，在饱和区这段时间执

24、行器仍停留在极限位置而不能随偏差的反向立即做出相应的改变，就好像系统失去了控制一样，造成控制性能恶化，这就是积分饱和。作为防止积分饱和方法之一就是抗积分饱和法。思路是在计算U(k)时，首先判断上一时刻控制器的输出量U(k-1)是否已经超出限制范围若 U(k-1)>= Umax,则只累加负偏差，正偏差则去掉积分作用； 若 U(k-1)<=-Umax,则只累加正偏差，负偏差则去掉积分作用；程序如下：clear allclose allts=0.001;sys=tf(5.235e005,1,87.35,1.047e004,0);dsys=c2d(sys,ts,'z

25、9;);num,den=tfdata(dsys,'v');u_1=0.0;u_2=0.0;u_3=0.0;y_1=0;y_2=0;y_3=0;x=0,0,0'error_1=0;um=6;kp=0.85;ki=9.0;kd=0.0;rin=30; %step signalfor k=1:1:800 time(k)=k*ts; u(k)=kp*x(1)+kd*x(2)+ki*x(3); %?PID?Controller if u(k)>=um u(k)=um; end if u(k)<=-um u(k)=-um; end %Linear?model? yout

26、(k)=-den(2)*y_1-den(3)*y_2-den(4)*y_3+num(2)*u_1+num(3)*u_2+num(4)*u_3; error(k)=rin-yout(k); M=2; if M=1 %Using?intergration?sturation? if u(k)>=um if error(k)>0 alpha=0; else alpha=1; end elseif u(k)<=-um if error(k)>0 alpha=1; else alpha=0; end else alpha=1; end elseif M=2 % Not?using

27、?intergration?sturation alpha=1; end %Return?of?PID?parameters? u_3=u_2;u_2=u_1;u_1=u(k); y_3=y_2;y_2=y_1;y_1=yout(k); error_1=error(k); x(1)=error(k); %?Calculating?P? x(2)=(error(k)-error_1)/ts; %?Calculating?D? x(3)=x(3)+alpha*error(k)*ts; %?Calculating?I? xi(k)=x(3);endfigure(1);subplot(3,1,1);p

28、lot(time,rin,'b',time,yout,'r');xlabel('time(s)');ylabel('Position?tracking');subplot(3,1,2);plot(time,u,'r');xlabel('time(s)');ylabel('Controller?output');subplot(3,1,3);plot(time,xi,'r');xlabel('time(s)');ylabel('Integrati

29、on'); 结果如下： 结果如下：当M=2时，为无抗积分饱和，结果如图：当M=1时，为有抗积分饱和，结果如下4.4带死区的PID 控制算法被控对象为，采样时间为1ms，对象输出上有一个幅值为0.5的正态分布的随机干扰信号。取=0.20，死区参数e0=0.10，采用低通滤波器对输出信号进行滤波，滤波器为：。取M=2 时仿真结果分别是怎样的？原理：为了避免控制作用过于频繁，消除由于频繁动作所引起的震荡，可以采用带死区的PID控制算法控制，计算式为：程序：clear allclose allts=0.001;sys=tf(5.235e005,1,87.35,1.047e004,0);dsys

30、=c2d(sys,ts,'z');num,den=tfdata(dsys,'v');u_1=0;u_2=0;u_3=0;u_4=0;u_5=0;y_1=0;y_2=0;y_3=0;yy_1=0;error_1=0;error_2=0;ei=0;sys1=tf(1,0.04,1); %Low?Freq?Signal?Filter?dsys1=c2d(sys1,ts,'tucsin');num1,den1=tfdata(dsys1,'v');f_1=0;for k=1:1:2000 time(k)=k*ts; rin(k)=1; %Step?Signal %Linear?model? yout(k)=-den(2)*y_1-den(3)*y_2-den(4)*y_