天津大學(xué)自動化2010級控制系統(tǒng)設(shè)計與仿真實(shí)驗(yàn)報告

1、TIANJIN UNIVERSITY控制系統(tǒng)設(shè)計與仿真實(shí)驗(yàn)報告姓名：葉林奇 班級：自動化2班 年級：2010級 學(xué)號：3010203109第一次上機(jī)實(shí)驗(yàn)任務(wù)1、熟悉matlab軟件的運(yùn)行環(huán)境，包括命令窗體，workspace等，熟悉繪圖命令。2、采用四階龍格庫塔法求如下二階系統(tǒng)的在幅值為1脈寬為1刺激下響應(yīng)的數(shù)值解。，程序：1） function r=u(t)r=(t>=0)&(t<=1)2） dt=0.01;T=5;s=0.5;w=10;A=0 1;-w2 -2*s*wB=0;w2X=0;0;Y=X;for i=dt:dt:T K1=A*X+B*u(i); K2=A*(X

2、+K1*dt/2)+B*u(i); K3=A*(X+K2*dt/2)+B*u(i); K4=A*(X+K2*dt)+B*u(i); K=(K1+2*K2+2*K3+K4)/6; X=X+K*dt; Y=Y X;endplot(0:dt:T,Y(1,:),'r.')結(jié)果：3、采用四階龍格庫塔法求高階系統(tǒng)階單位躍響應(yīng)曲線的數(shù)值解。，程序：dt=0.01;N=50/dt;U=1;s=0.5;w=10;T=5;A=0 1 0;0 0 1;-w2/T -(w2*T+2*s*w)/T -(2*s*w*T+1)/T;B=0;0;w2/T;X=0;0;0;Y=X;for i=1:N K1=A*

3、X+B*U; K2=A*(X+K1*dt/2)+B*U; K3=A*(X+K2*dt/2)+B*U; K4=A*(X+K2*dt)+B*U; K=(K1+2*K2+2*K3+K4)/6; X=X+K*dt; Y=Y X;endy=Y(1,:);t=0:dt:50;plot(t,y)結(jié)果：4、 自學(xué)OED45指令用法，并求解題2中二階系統(tǒng)的單位階躍響應(yīng)。程序：1) function r=u1(x);r=(x>=0)&(x<=1)2) function dx=myderv(t,x);s=0.5;w=10;dx=x(2);-w2*x(1)-2*s*w*x(2)+w2*u1(t)3

4、) x0=0;0;t y=ode45(myderv,0,5,x0);plot(t,y(:,1)結(jié)果：第二次上機(jī)任務(wù)1、試用simulink方法解微分方程，并封裝模塊，輸出為。得到各狀態(tài)變量的時間序列，以及相平面上的吸引子。參數(shù)入口為的值以及的初值。（其中，以及初值分別為） 提示：模塊輸入是輸出量的微分。simulink:封裝模塊：程序：t x y=sim('task2_1')figure(1)plot(t,y)figure(2)plot3(y(:,1),y(:,2),y(:,3)結(jié)果：x1、x2、x3時間序列：x1、x2、x3在相平面構(gòu)成的吸引子：2、用simulink搭建PI

5、控制器的控制回路，被控對象傳遞函數(shù)：，分別分析 （1）、比例系數(shù)由小到大以及積分時間由小到大對階躍響應(yīng)曲線的影響。 （2）、控制器輸出有飽和以及反饋有時滯情況下，階躍響應(yīng)曲線的變化。 （3）、主控制回路傳遞函數(shù)為：，副回路為：，主回路采用PI控制器，副回路采用P控制器，分析控制系統(tǒng)對主回路以及副回路的階躍擾動的抑制。注：PI控制器表達(dá)式為，串級控制如圖所示。(1)simulink：程序：Ti=1;figure(1)for Kp=0.5:0.5:5;t x y=sim('PID')hold onplot(t,y)endKp=1;figure(2)for Ti=0.5:0.5:8;

6、t x y=sim('PID')hold onplot(t,y)end結(jié)果：Kp增大結(jié)論：隨著Kp增大，比例作用增強(qiáng)，系統(tǒng)快速性提高Ti增大結(jié)論：隨著Ti增大，積分作用減弱，系統(tǒng)超調(diào)量減小，同時快速性減弱（2）simulink:結(jié)果：輸出有飽和反饋有時滯原響應(yīng)曲線結(jié)論：當(dāng)輸出有飽和時，穩(wěn)態(tài)誤差增大；當(dāng)反饋有時滯時，階躍響應(yīng)振蕩加劇。(3)simulink：結(jié)果：副回路加入干擾主回路加入干擾t=100時加入干擾結(jié)論：系統(tǒng)對副回路擾動具有更強(qiáng)的抑制作用。3、編寫S函數(shù)模塊，實(shí)現(xiàn)兩路正弦信號的疊加，正弦信號相位差為60度。S函數(shù)：functionsys,x0,str,ts=syst(

7、t,x,u,flag)switch flag, case 0 sys,x0,str,ts=mdlInitializeSizes(); case 3 sys=mdlOutputs(u); case 1,2,4,9,sys=; otherwise error('ErrFLag=',num2str(flag);endfunction sys,x0,str,ts=mdlInitializeSizes()sizes=simsizes;sizes.NumContStates=0;sizes.NumDiscStates=0;sizes.NumOutputs=1;sizes.NumInputs

8、=2;sizes.DirFeedthrough=1;sizes.NumSampleTimes=1;sys=simsizes(sizes);x0=;str= ;ts=-1 0;function sys=mdlOutputs(u)sys=u(1)+u(2);simulink：結(jié)果：4、熟悉SimPowerSystem模塊，試用Thyristor Converter模塊以及Synchronized 6-pulse generator模塊實(shí)現(xiàn)三相電整流。（自學(xué)內(nèi)容）simulink：(三相橋式整流，阻感負(fù)載)結(jié)果：=0 =30=60 =90 5、利用觸發(fā)子系統(tǒng)構(gòu)建一零階保持器，實(shí)現(xiàn)對正弦信號的采樣，并

9、比較不同采樣周期下的采樣波形。Simulink：結(jié)果：6、若被控對象傳遞函數(shù)為,控制器為，試用simulink搭建一單位反饋控制系統(tǒng)，分析采樣周期T對系統(tǒng)單位階躍響應(yīng)的影響。simulink：程序：T=0.01;for i=1:3T=10*T; sim('task2_6s')plot(tout,yout,'r')hold onend結(jié)果：T=10T=1T=0.1從上圖可以看出，隨著采樣周期T增大，階躍響應(yīng)的振蕩減小，超調(diào)量減少；快速性降低，上升時間變長。7、設(shè)一單位反饋控制系統(tǒng)，控制器采用PI控制，Kp=200,Ki=10, 控制器飽和非線性寬度為2，受控對象為

10、時變模型，由微分方程給出，如下： 求系統(tǒng)單位階躍響應(yīng)，并分析不同Kp取值對響應(yīng)曲線的影響。Simulink:程序：Kp=200;Ki=10;t x y=sim('task2_7s');plot(t,y)結(jié)果：改變Kp：Ki=10;for Kp=2 20 200t x y=sim('task2_7s');plot(t,y)hold onendKp=20Kp=200Kp=2結(jié)論：Kp增大，系統(tǒng)由不穩(wěn)定變?yōu)榉€(wěn)定第三次上機(jī)任務(wù)1、熟悉控制系統(tǒng)各個模型表示方法的命令以及它們之間的相互轉(zhuǎn)化。（展開形式，零極點(diǎn)形式，狀態(tài)空間形式以及部分分式形式。）2、試用至少三種方法，判斷一

11、下系統(tǒng)的穩(wěn)定性:（1） 程序：num=1 2 3 1;den=1 5 2 1 1;G=tf(num,den)%*·方法1*z,p,k=tf2zp(num,den)%*·方法2*r=roots(den)%*·方法3*pzmap(num,den)結(jié)論：系統(tǒng)有兩個不穩(wěn)定極點(diǎn)(實(shí)部>0)（2）程序：A=1 3;5 2;%*方法1*P1=poly(A);r1=roots(P1)%*方法2*r1=eig(A)%求矩陣A的全部特征值%*方法*Q=eye(size(A);P=lyap(A,Q)det=det(P)結(jié)論：系統(tǒng)有1個不穩(wěn)定極點(diǎn)(實(shí)部>0)。用李雅普諾夫第二

12、法也可看出系統(tǒng)不穩(wěn)定，因所得矩陣P不是正定。3、試產(chǎn)生一周期為5秒，時長為30秒，最大值為1，最小值為0的三角波；得到如下一階系統(tǒng)在三角波輸入下的時間響應(yīng)曲線。 程序：x = 0:0.01:30;y1=mod(x,2.5)/2.5;y2=1-mod(x,2.5)/2.5;y=0for i=0:5:30y=y+y1.*(x>=i)&(x<=i+2.5)+y2.*(x>=i+2.5)&(x<i+5);endnum=1den=2 1lsim(num,den,y,x)結(jié)果：5、對如下二階系統(tǒng)做時域分析，得到阻尼比在01之間變化的時候，階躍響應(yīng)的上升時間，調(diào)節(jié)時間

13、，峰值時間，超調(diào)量以及衰減比（第一個峰值與穩(wěn)態(tài)值之差與第二個峰值與穩(wěn)態(tài)值之差的比）其中。 程序：t=0:0.005:10;w=5;num=w2;tr=0;ts=0;tp=0;pos=0;n=0;for s=0.2 0.4 0.6 0.8 den=1 2*s*w w2; y x t=step(num,den,t); mp,dtp=max(y); %峰值mp和對應(yīng)的時間點(diǎn) tp=tp;t(dtp); %峰值時間tp yss=dcgain(num,den); %穩(wěn)態(tài)數(shù)值yss pos=pos;100*(mp-yss)/yss; %超調(diào)量pos dtr=min(find(abs(y-yss)<0

14、.01); tr=tr;t(dtr); %上升時間tr dts=max(find(abs(y-0.95*yss)<0.01|abs(y-1.05*yss)<0.01); ts=ts;t(dts); %調(diào)節(jié)時間ts peak=findpeaks(y); n=n;(peak(1)-yss)/(peak(2)-yss); %衰減比n step(num,den,t) hold onendtable=tr ts tp pos n結(jié)果：阻尼比增大衰減比超調(diào)量峰值時間調(diào)節(jié)時間上升時間0.80.60.40.2阻尼比6、已知開環(huán)傳遞函數(shù)如下，1）試用根軌跡方法得到其臨界穩(wěn)定增益。2）若k=10,試用

15、伯德圖方法，判斷其穩(wěn)定性。程序：num=1;den=conv(2,1,conv(1,1,0.1,1); figure(1)rlocus(num,den);K,poles=rlocfind(num,den)num=10;figure(2)margin(num,den)結(jié)果：（1）根軌跡：臨界穩(wěn)定增益：（2）結(jié)論：穿越頻率對應(yīng)的相角在-180度線上方，故k=10時，系統(tǒng)穩(wěn)定7、已知系統(tǒng)開環(huán)傳遞函數(shù)如下試設(shè)計一超前校正環(huán)節(jié)，使得超調(diào)量為20%，調(diào)節(jié)時間為1s。系統(tǒng)單位斜坡穩(wěn)態(tài)響應(yīng)誤差為10%。并作出校正前后后的系統(tǒng)單位階躍響應(yīng)時域曲線加以比較。（1）求校正器增益 Kc在斜坡信息作用下，系統(tǒng)的穩(wěn)態(tài)誤差

16、ess1/Kv= 1/2Kc =0.1,可得：Kv=5（2）校驗(yàn)原系統(tǒng)的階躍響應(yīng)超調(diào)量是否滿足要求t=0:0.02:50;n1=10; d1=conv(conv(1 0,0.5 1),0.1 1);sope=tf(n1,d1)sys=feedback(sope,1);y,t=step(sys,t);mp,dtp=max(y); %峰值mp和對應(yīng)的時間點(diǎn)L=length(t); %數(shù)據(jù)長度yss=y(L) %穩(wěn)態(tài)數(shù)值ysspos=100*(mp-yss)/yss %超調(diào)量posdts=max(find(abs(y-0.98*yss)<0.01);ts=t(dts) %調(diào)節(jié)時間ts結(jié)果：超調(diào)

17、量=85.62%>20%調(diào)節(jié)時間=35.22s>1s不滿足要求 （3）確定期望極點(diǎn)位置pos=0.2; %期望超調(diào)量ts=1; %期望調(diào)節(jié)時間s=(log(1/pos)2/(pi)2+(log(1/pos)2)(1/2)%阻尼比wn=4.5/(s*ts)%無阻尼自然震蕩頻率p=1 2*s*wn wn2;roots(p)%求系統(tǒng)的主導(dǎo)極點(diǎn)結(jié)果：(4)求校正器的傳遞函數(shù)kc=5;s_1=-4.5000 + 8.7839i;nk1=2;dk1=conv(conv(1 0,0.5 1),0.1 1);ngv=polyval(nk1,s_1);dgv=polyval(dk1,s_1);g=ngv/dgv;zetag=angle(g);mg=abs(g);ms=abs(s_1);zetas=angle(s_1);tz=(sin(zetas)-kc*mg*sin(zetag-zetas)/(kc*mg*ms*sin(zetag);tp=-(kc*mg*sin(zetas)+sin(zetag+zetas)/(ms*sin(zetag);nk=tz,1;dk=tp,1