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

1、.自动控制技术上机实验报告班级：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 =

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

3、%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('速度

4、（米/秒）');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 =

5、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.5897i

6、k2 = 源程序代码：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)<1.001) r=r+1;end;rise_time=(r-1)*0.005 %求系统上升时间ymax,tp=max(y);peak_time=(tp-1)*0.005k=dcgain(sys);max_overshoot=(ymax-k)/ks=7001;while(y(s)>0.98&

7、y(s)<1.02) s=s-1;endsetting_time=(s-1)*0.005运行结果：rise_time = 0.7050peak_time = 34.6100max_overshoot = 2.6819e+06setting_time =35源程序代码：num=1 1;den=1 5 6 0;sys1=tf(num,den)subplot(3,1,1)step(sys1)subplot(3,1,2)rlocus(num,den)title('系统的根轨迹曲线')r=rlocus(num,den,20.575)a=1;b=1 0.8989;sys2=tf(a,

8、b)sys=series(sys1,sys2) %串联subplot(3,1,3)step(sys)运行结果：rise_time = 0.7050peak_time = 34.6100max_overshoot = 2.6819e+06setting_time = 35>> 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 Continuo

9、us-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 Co

10、ntinuous-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;

11、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运行结果：

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

13、%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

14、 = -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;B

15、1 = 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

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

17、 = 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('系统是可控的')

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

19、_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:1

20、000 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);

21、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('误差曲线

22、')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; ri

23、n(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)<=40 beta=0.3; elseif abs(error(k)>=20&abs(error(k)<=30 beta=0.6; elseif abs(error(k)>=10&abs(error(k)<=20 beta=0.9; else beta=1.0; e

24、nd elseif M=2 beta=1.0; end ei=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; end if u(k)<=-110 u(k)=-110; end u_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,

25、'r');title('不加积分作用效果')xlabel('time(s)');ylabel('rin,yout');figure(2);plot(time,u,'r');xlabel('time(s)');ylabel('u');运行结果：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,

26、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;for k=1:1:800 time(k)=k*ts; u(k)=kp*x(1)+kd*x(2)+ki*x(3); if u(k)>=um u(k)=um; end if u(k)<=-um y(k)=-um; end yout(k)=-den(2)*y_1-den(3)*y_2-den(4)*y_3+num(2)*u_1+num(

27、3)*u_2+num(4)*u_3; error(k)=rin-yout(k); M=1; if M=1 if u(k)>=num 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 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)=err

28、or(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('Integra

29、tion');运行结果：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