西电电院自动控制技术上机报告.doc

自动控制技术上机实验报告班级：021215 学号：021214 姓名：时域分析程序源代码：close all;clear all;ft = 30;M=1;B=5;K=20; %系统参数t0=0;tfinal = 5;tspan = t0 tfinal; %设置仿真开始和结束时间x0 = 0,0; %系统初始值，零初始条件options = odeset(AbsTol,1e-6;1e-6); %设置仿真计算精度t,x = ode113(xt4odefile,tspan,x0,options);%微分方程求解，计算位移x(:,1)和速度x(:,2) 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 disp(稳态时系统的位移、速度和加速度及对应的时间分别为： );%显示计算结果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,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(速度（米/秒）);title(时间-速度曲线);grid; subplot(1,3,3),plot(d,v); %绘制位移-速度曲线xlabel(位移（米）);ylabel(速度（米/秒）);title(位移-速度曲线);grid; 运行结果： EX11稳态时系统的位移、速度和加速度及对应的时间分别为： 位移 d=1.5000速度 v=-0.000086加速度 a=-0.000084时间 t=4.46源程序代码：close all;clear all; num = 0,2,5,7;den = 1,6,10,6;z1,p1,k1 = tf2zp(num,den) r2,p2,k2 = residue(num,den) 运行结果： EX12z1 = -1.2500 + 1.3919i -1.2500 - 1.3919ip1 = -3.7693 + 0.0000i -1.1154 + 0.5897i -1.1154 - 0.5897ik1 = 2r2 = 2.2417 + 0.0000i -0.1208 - 1.0004i -0.1208 + 1.0004ip2 = -3.7693 + 0.0000i -1.1154 + 0.5897i -1.1154 - 0.5897ik2 = 源程序代码：clca=6.3223 18 12.811;b=1 6 1.322 18 12.811;sys=tf(a,b);t=0:.005:35;step(sys)title(系统的单位阶跃响应)y,t=step(a,b,t);r=1;while(y(r)0.98&y(s) EX14sys1 = s + 1 - s3 + 5 s2 + 6 s Continuous-time transfer function.r = -2.0505 + 4.3225i -2.0505 - 4.3225i -0.8989 + 0.0000isys2 = 1 - s + 0.8989 Continuous-time transfer function.sys = s + 1 - s4 + 5.899 s3 + 10.49 s2 + 5.393 s Continuous-time transfer function.频域分析：源程序代码：clcnum=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)运行结果：sys = 0.01 s2 + 0.0001 s + 0.01 - 0.25 s4 + 0.01 s3 + s2 Continuous-time transfer function.sys2 = 0.01 s2 + 0.0001 s + 0.01 - 0.25 s4 + 0.01 s3 + 1.01 s2 + 0.0001 s + 0.01 Continuous-time transfer function.源程序代码：close all;clear all; num = 0 20 20 10; %开环传递函数分子den = conv(1 1 0,1 10); %开环传递函数的分母nyquist(num,den)%v = -2 3 -3;axis(-2 2 -3 3)grid x = -pi:0.01:pi;plot(x,sin(x), grid on运行结果：源程序代码：close all;clear all; num = 2000,2000; %开环传递函数的分子den = conv(1 0.5 0,1 14 400); %开环传递函数的分母nichols(num,den) %绘制nichols图v = -270 -90 -40 40;axis(v)ngrid %标出nichols图线运行结果：源程序代码：num=0 2000 2000;den=conv(1 0.5 0,1 14 400);h=tf(num,den);gm,pm,wg,wc=margin(h);gm,pm,wg,wc运行结果： EX24gm = 2.7493pm = 73.3527wg = 19.8244wc =5.3477源程序代码：clcnum=1;den=0.5 1.5 1 0;sys=tf(num,den)sys2=feedback(sys,1)bode(sys2)gm pm wg wp=margin(sys)运行结果： - 0.5 s3 + 1.5 s2 + s + 1 Continuous-time transfer function.gm = 3.0000pm = 32.6133wg = 1.4142wp =0.7494现代控制理论源程序代码：close all;clear all; %3.1_Anum = 1 2 3; %传递函数分子多项式的系数den = 1 3 3 1; %传递函数分母多项式的系数A,B,C,D = tf2ss(num,den) %3.1_Bz = -1;-3; %传递函数的零点p = 0;-2;-4;-6; %传递函数的极点k = 4;A,B,C,D = zp2ss(z,p,k) %3.1_CA = 0,1;-1,-2;B = 0;1;C = 1,3;D = 1;num,den = ss2tf(A,B,C,D)printsys(num,den,s)z,p,k = ss2zp(A,B,C,D) 运行结果： EX31A = -3 -3 -1 1 0 0 0 1 0B = 1 0 0C = 1 2 3D = 0A = -10.0000 -4.8990 0 0 4.8990 0 0 0 -6.0000 -4.2866 -2.0000 0 0 0 1.0000 0B = 1 0 1 0C = 0 0 0 4D = 0num = 1.0000 5.0000 2.0000den = 1 2 1 num/den = s2 + 5 s + 2 - s2 + 2 s + 1z = -0.4384 -4.5616p = -1 -1k = 1源程序代码：close all;clear all; A1 = 0,1;-1,-2;B1 = 0;1;C1 = 1,3;D1 = 1;A2 = 0,1;-1,-3;B2 = 0;1;C2 = 1,4;D2 = 0;A,B,C,D = series(A1,B1,C1,D1,A2,B2,C2,D2)A,B,C,D = parallel(A1,B1,C1,D1,A2,B2,C2,D2)A,B,C,D = feedback(A1,B1,C1,D1,A2,B2,C2,D2)A,B,C,D = feedback(A1,B1,C1,D1,A2,B2,C2,D2,+1)运行结果： EX32A = 0 1 0 0 -1 -3 1 3 0 0 0 1 0 0 -1 -2B = 0 1 0 1C = 1 4 0 0D = 0A = 0 1 0 0 -1 -2 0 0 0 0 0 1 0 0 -1 -3B = 0 1 0 1C = 1 3 1 4D = 1A = 0 1 0 0 -1 -2 -1 -4 0 0 0 1 1 3 -2 -7B = 0 1 0 1C = 1 3 -1 -4D = 1A = 0 1 0 0 -1 -2 1 4 0 0 0 1 1 3 0 1B = 0 1 0 1C = 1 3 1 4D = 1源程序代码：clc;close all;clear all; A = 0 -2;1 -3;t = 0.2;Phi = expm(A*t); %求状态转移矩阵B = 2;0;C = 0 3;D = 0;x0 = 1 1;t = 0 0.2;u = 0*t;y,x = lsim(A,B,C,D,u,t,x0) %求系统响应运行结果：y = 3.0000 2.0110x = 1.0000 1.0000 0.6703 0.6703源程序代码：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)f=rank(S)if(f=n) disp(系统是可控的)else disp(系统是不可控的)enddisp()disp(可观测矩阵：)V=obsv(A,C)m=rank(V)if(f=n) disp(系统时可观测的)else disp(系统是不可观测的)end 运行结果：可控矩阵：S = 1 1 -2 -2 1 1 -2 -2f = 1系统是不可控的可观测矩阵：V = 1 1 1 -1 -2 -2 -4 4m = 2系统是不可观测的源程序代码：clcA=0 1;-2 -3;B=0;1;C=2 0;D=0;P_S=-1 -2;k=acker(A,B,P_S);P_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)运行结果：校正设计源程序代码：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;for k=1:1:1000 time(k)=k*ts; kp=8;ki=0.1;kd=10; rin(k)=0.5*sin(2*pi*k*ts); du(k)=kp*x(1)+kd*x(2)+ki*x(3); u(k)=u_1+du(k); yout(k)=-den(2)*y_1-den(3)*y_2+num(2)*u_1+num(3)*u_2; error(k)=rin(k)-yout(k); u_2=u_1;u_1=u(k); y_2=y_1;y_1=yout(k); x(1)=error(k)-error_1; x(2)=error(k)-2*error_1+error_2; x(3)=error(k); error_2=error_1; error_1=error(k);endfigure(1);plot(time,rin,b,time,yout,r),grid ongtext(rinrightarrow)gtext(leftarrowyout)title(系统输出曲线)xlabel(time(s),ylabel(rin,yout);figure(2);plot(time,error,r),grid ontitle(误差曲线)xlabel(time(s);ylabel(error);注：输入信号为正弦，采样信号为1ms。运行结果：源程序代码：clc;clear all;ts=20;sys=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=0;y_1=0;y_2=0;y_3=0;error_1=0;ei=0;for k=1:1:200 time(k)=k*ts; rin(k)=40; kp=0.80; ki=0.005; kd=3.0; yout(k)=-den(2)*y_1+num(2)*u_5; error(k)=rin(k)-yout(k); M=2; if M=1 if abs(error(k)=30&abs(errpr(k)=20&abs(error(k)=10&abs(error(k)=110 u(k)=110; end if u(k)=um u(k)=um; end if u(k)=num if error(k)0 alpha=0; else alpha=1; end elseif u(k)0 alpha=1; else alpha=0; end else alpha=1; end elseif M=2 alpha=1; end 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); x(2)=(error(k)-error_1)/ts; x(3)=x(3)+alpha*error(k)*ts; xi(k)=x(3);endfigure(1);subplot(3,1,1);plot(time,rin,b,time,yout,r);xlabel(time(s);ylabel(Controller?output);subplot(3,1,3);plot(time,xi,r);xlabel(time(s);ylabel(Integration);运行结果：M=2时，无抗积分饱和：M=1时，有抗积分饱和：源程序代码：clc;clear all;ts=0.001;sys=tf(5.235e005,1,87.35,1.047e004,0);dsys=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);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; yout(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; D(k)=0.50*rands(1); yyout(k)=yout(k)+D(k); filty(k)=-den1(2)*f_1+num1(1)*(yyout(k)+yy_1); error(k)=rin(k)-filty(k); if abs(error(k)=0.20 ei=ei+error(k)*ts; else ei=0; end kp=0.50;ki=0.10;kd=0.020; u(k)=kp*error(k)+ki*ei+kd*(error(k)-error_1)/ts; M=2; if M=1 u(k)=u(k); elseif M=2 if abs(error(k)=10 u(k)=10; end if u(k)=-10 u(k)=-10; end rin_1=rin(k); u_3=u_2;u_2=u_1;u_1=u(k); y_3=y_2;y_2=y_1;y_1=yout(k); f_1=filty(k); yy_1=yyout(k); error_2=error_1; error_1=error(k);endfigure(1);subplot(2,1,1);plot(time,rin,r,time,filty,b);xlabel(time(s);ylabel(rin,yout);subplot(2,1,2);plot(time,u,r);xlabel(time(s);ylabel(u);figure(2);plot(time,D,r);xlabel(time(s);ylabel(Disturbance?signal);运行结果：M=1时，不带死区结果：M=2时，带死区结果：