自動控制上機作業(yè)

1、自動控制上機大作業(yè)班級：學號：姓名：1.1設質量阻尼彈簧系統的微分運動方程為 式中，x(t)為位移輸出信號，f(t)為輸入的力信號。質量為M=1kg，粘性摩擦系數為B = 5N / m s1，彈簧的彈性系數為K=20N/m。當t=0 時，施加外力f(t)=30N，試問系統何時達到穩(wěn)定？并畫出該機械系統位移、速度隨時間變化的曲線以及速度與位移的關系曲線。提示：龍格-庫塔法求解微分方程數值解的函數：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,tspan,x0,option

3、s);%微分方程求解。 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,a(i);disp(result);result = sprintf(時間 t=%4.2fn,t(i);di

4、sp(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;其中xt4ode

5、file.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.加速度 a=-0.時間 t=4.461.2假設控制系統的傳遞函數為試求其零點、極點和增益，并進行部分分式展開。提示：傳遞函數描述：tf(), 調用方式：SYS = TF(NUM,DEN)。求取零點和極點的函數：tf2zp(), 調用方式：Z,P,K = TF2ZP(NUM,DEN)傳遞函數的部分分式展開：residue(), 調用方式：R,P,K = RESIDUE(

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

7、 -0.1208 - 1.0004i -0.1208 + 1.0004iP = -3.7693 -1.1154 + 0.5897i -1.1154 - 0.5897iK = 1.3 考慮由下式表示的高階系統 ，試求取系統的單位階躍響應，并計算系統的上升時間、峰值時間、超調量和調整時間（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);

8、yss=y(cs)%穩(wěn)態(tài)ct=(mp-yss)/yss%超調量結果：法二clcsys=tf(6.3223 18 12.811,1 6 11.3223 18 12.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(s

9、ys)figure(2)sys2=feedback(sys,1)%加單位反饋bode(sys2)結果：系統總的波特圖2.2 clcnum=20,20,10;dun=1,11,10,0;sys=tf(num,dun)nyquist(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(s

10、ys)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 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(

11、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 -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,den

12、1,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=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

13、;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)%計算可控矩陣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程序如下：cl

14、cA=0 1;-2 -3;B=0;1;C=2 0;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*C A-B*k-h*C;B1=B;B;C1=C,zeros(1,2);D1=0;sys=ss(A1,B1,C1,D1)%建立復合系統動態(tài)模型tf(sys)結果如下：校正設計4.1增量式PID 控制算法被控對象為，PID 控制參數為：kp=8，ki=0.10，kd=10。下面程序是輸入為正弦信號，采樣信號是1ms，控制器出

15、的曲線以及誤差曲線程序：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: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);%Restrict

16、ing 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 I error_

17、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(誤差曲線)xlabel(time(s);ylabel(error);結果如下：4.2積分分離PID 控制算法被控對象為，采樣時間為20s，延遲時間為4 個采樣時間。解：離散積分分離控

18、制算法表達式：程序如下：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=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?separati

19、on?error(k)=rin(k)-yout(k);M=2;if M=1 %Using?integration?separation? if abs(error(k)=30&abs(error(k)=20&abs(error(k)=10&abs(error(k)=110 u(k)=110;endif u(k)= Umax,則只累加負偏差，正偏差則去掉積分作用；若 U(k-1)=um u(k)=um; end if u(k)=um if error(k)0 alpha=0; else alpha=1; end elseif u(k)0 alpha=1; else alpha=0; end el

20、se alpha=1; end elseif M=2 % Not?using?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

21、(k)=x(3);endfigure(1);subplot(3,1,1);plot(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(Integration); 結果如下： 結果如下：當M=2時，為無抗積分飽和，結果如圖：當M=1時，為有抗積分飽和，結果如下4.4帶死區(qū)的

22、PID 控制算法被控對象為，采樣時間為1ms，對象輸出上有一個幅值為0.5的正態(tài)分布的隨機干擾信號。取=0.20，死區(qū)參數e0=0.10，采用低通濾波器對輸出信號進行濾波，濾波器為：。取M=2 時仿真結果分別是怎樣的？原理：為了避免控制作用過于頻繁，消除由于頻繁動作所引起的震蕩，可以采用帶死區(qū)的PID控制算法控制，計算式為：程序：clear allclose allts=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

23、=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_3+num(2)*u_1+num(3)*u_2+num(4)*u_3; D(k)=0.50*rands(1); %Disturbance?signal? y