自动控制上机作业

自动控制上机大作业班级：学号：姓名：1.1设质量——阻尼——弹簧系统的微分运动方程为)()()()(M22tftKxdttdxBdttxd式中，x(t)为位移输出信号，f(t)为输入的力信号。质量为M=1kg，粘性摩擦系数为B=5N/m⋅s−1，弹簧的弹性系数为K=20N/m。当t=0时，施加外力f(t)=30N，试问系统何时达到稳定？并画出该机械系统位移、速度随时间变化的曲线以及速度与位移的关系曲线。提示：龙格-库塔法求解微分方程数值解的函数：odel13(),调用方式：[T,Y]=ODE113(ODEFUN,TSPAN,Y0,OPTIONS)。其中ODEFUN为用户自定义的系统微分方程的描述，本题中可使用xt4odefile.m文件定义的函数；TSPAN表示计算开始和结束的时间；Y0表示微分方程的初始条件；OPTION为计算精度的可选参数，由odese()函数设置。odel13()函数只接受一阶微分方程的形式，使用时需要先将高阶方程化为若干个一阶微分方程；绘图函数：plot()，subplot();程序：ft=30;M=1;B=5;K=20;%系统参数tspan=[020];%设置仿真开始和结束时间x0=[0,0];%系统初始值，零初始条件options=odeset('AbsTol',[1e-6;1e-6]);%设置仿真计算精度[t,x]=ode113('xt4odefile',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.4f\n',x(i,1));disp(result);result=sprintf('速度v=%8.6f\n',x(i,2));disp(result);result=sprintf('加速度a=%9.6f\n',a(i));disp(result);result=sprintf('时间t=%4.2f\n',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('速度（米/秒）');title('时间-速度曲线');grid;subplot(1,3,3),plot(d,v);%绘制位移-速度曲线xlabel('位移（米）');ylabel('速度（米/秒）');title('位移-速度曲线');grid;其中xt4odefile.m文件为：functionxt=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.461.2假设控制系统的传递函数为)6106()752(232sssss试求其零点、极点和增益，并进行部分分式展开。提示：传递函数描述：tf(),调用方式：SYS=TF(NUM,DEN)。求取零点和极点的函数：tf2zp(),调用方式：[Z,P,K]=TF2ZP(NUM,DEN)传递函数的部分分式展开：residue(),调用方式：[R,P,K]=RESIDUE(B,A)clcsys=tf([257],[16106])disp('零点极点分别为Z,P')[Z,P,K]=tf2zp([257],[16106])B=[257];A=[16106];disp('部分分式展开')[R,P,K]=RESIDUE(B,A)计算结果：Transferfunction:2s^2+5s+7----------------------s^3+6s^2+10s+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考虑由下式表示的高阶系统811.12183223.116811.1218s6.32232342sssss，试求取系统的单位阶跃响应，并计算系统的上升时间、峰值时间、超调量和调整时间（2%误差带）。提示：阶跃响应函数：step(),调用方式：[Y,T]=STEP(SYS)。法一clcsys=tf([6.32231812.811],[1611.32231812.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.32231812.811],[1611.32231812.811])[y,t]=step(sys);ltiview(sys)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.2clcnum=[20,20,10];dun=[1,11,10,0];sys=tf(num,dun)nyquist(sys)结果如下：2.3clcnum=[20002000];den=[114.54072000];sys=tf(num,den)nichols(sys)v=[-270-90-4040];axis(v)ngrid结果：2.4clcnum=[20002000];den=[114.54072000];sys=tf(num,den)figure(1)bode(sys)[gmpmwgwp]=margin(sys)结果如下：2.5clcnum=[1];den=[0.51.510];sys=tf(num,den)sys2=feedback(sys,1)bode(sys2)[gmpmwgwp]=margin(sys)结果如下：3.1A程序如下：clcnum=[123];den=[1331];[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=[01;1-2];B=[0;1];C=[13];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=[01;-1-2];B1=[1;0];C1=[13];D1=1;A2=[01;-1-3];B2=[0;1];C2=[14];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）程序如下：clcA=[0-2;1-3];t=0.2;eAt=expm(A*t)结果如下：（2）结果如下：clcA=[0-2;1-3];B=[2;0];C=[03];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=[-31;1-3];B=[11;11];C=[11;1-1];D=[00;00];N=size(A);n=N(1);[num,den]=ss2tf(A,B,C,D,2);disp('可控矩阵')S=ctrb(A,B)%计算可控矩阵Sf=rank(S)%通过rank命令求可控矩阵if(f==n)%判断系统的能空性disp('系统是可控的')elsedisp('系统是不可控的')enddisp('')disp('可观测矩阵')V=obsv(A,C)%计算可观测性矩阵Vm=rank(V)%求可观测性矩阵的秩if(f==n)%判断系统的能观测性disp('系统是可观测的')elsedisp('系统是不可观测的')end结果如下：3.5程序如下：clcA=[01;-2-3];B=[0;1];C=[20];D=0;P_S=[-1-2];%系统的配置极点k=acker(A,B,P_S);%计算系统的反馈增益向量kP_O=[-3-3];%观测器的期望配置极点h=(acker(A',C',P_O))';%计算观测器输出反馈阵A1=[A,-B*k;h*CA-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控制算法被控对象为2400()50Gsss，PID控制参数为：kp=8，ki=0.10，kd=10。下面程序是输入为正弦信号，采样信号是1ms，控制器出的曲线以及误差曲线程序：ts=0.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;fork=1:1:1000time(k)=k*ts;kp=8;ki=0.1;kd=10;%SineSignalrin(k)=0.5*sin(2*pi*k*ts);du(k)=kp*x(1)+kd*x(2)+ki*x(3);%PIDControlleru(k)=u_1+du(k);%Restrictingtheoutputofcontroller%Linearmodelyout(k)=-den(2)*y_1-den(3)*y_2+num(2)*u_1+num(3)*u_2;error(k)=rin(k)-yout(k);%Returnofparametersu_2=u_1;u_1=u(k);y_2=y_1;y_1=yout(k);x(1)=error(k)-error_1;%CalculatingPx(2)=error(k)-2*error_1+error_2;%CalculatingDx(3)=error(k);%CalculatingIerror_2=error_1;error_1=error(k);endfigure(1);plot(time,rin,'b',time,yout,'r'),g