自控系統(tǒng)數(shù)字仿真實(shí)驗(yàn)

1、實(shí)驗(yàn)一  面向微分方程的數(shù)值積分法仿真一、實(shí)驗(yàn)?zāi)康?掌握數(shù)值積分法的基本概念、原理及應(yīng)用；2用龍格-庫塔法解算微分方程，增加編寫仿真程序的能力；3分析數(shù)值積分算法的計(jì)算步長與計(jì)算精度、速度、穩(wěn)定性的關(guān)系；4. 對數(shù)值算法中的“病態(tài)問題”進(jìn)行研究。二、實(shí)驗(yàn)內(nèi)容1、已知系統(tǒng)微分方程及初值條件取步長h=0.1 ，試分別用歐拉方程法和RK4法求t=2h時(shí)的y值，并將求得的值與解析解比較(將三個(gè)解繪于同一坐標(biāo)中，且用數(shù)值進(jìn)行比較)，說明造成差異的原因。（編程完成；選用MATLAB ode函數(shù)完成。） 程序如下：t0=0;tf=2;h=0.1;y1=1;y2=1;y

2、3=1;t1=0;t2=0;t3=0;n=round(tf-t0)/h;for i=1:n %歐拉法 y1(i+1)=y1(i)+h*(2*h+y1(i); t1=t1,t1(i)+h;endfor i=1:n %四階龍格-庫塔法 k1=y2(i)+t2(i); k2=y2(i)+h*k1/2+t2(i)+h/2; k3=y2(i)+h*k2/2+t2(i)+h/2; k4=y2(i)+h*k3+t2(i)+h; y2(i+1)=y2(i)+h*(k1+2*k2+2*k3+k4)/6; t2=t2,t2(i)+h;endfor i=1:n y3(i+1)=2*exp(t3(i)-t3(i)-1

3、; t3=t3,t3(i)+h;endplot(t1,y1,'r',t2,y2,'g',t3,y3,'k')分析:紅線為用歐拉法得到的結(jié)果，綠線為用四階龍格庫塔法得到的結(jié)果，藍(lán)線為根據(jù)解析方程得到的結(jié)果。其差異原因主要有兩個(gè)：（1）二者的方法不同，歐拉法是根據(jù)一階微分方程計(jì)算得到的，龍格庫塔法是根據(jù)四階微分方程得到的；（2）由于步長取為0.1，所以得到的圖像與解析解之間存在差異，若將步長取小，則得到的解將更靠近解析解。  2、已知系統(tǒng)的傳遞函數(shù)為在單位階躍輸入下，系統(tǒng)響應(yīng)的解析解為 試分別用歐拉方程法和RK4法對系統(tǒng)進(jìn)行仿真（

4、編程完成）：1）比較兩種數(shù)值積分解與解析解的逼近程度；（繪圖）2）改變步長，分析步長對數(shù)值解精度的影響；3）不斷加大步長，分析計(jì)算穩(wěn)定性的變化。程序如下： a=0 1 0; 0 0 1; -22.06 -27 -10;b=0;0;1; c=40.6 0 0;X1=0;0;0; t=0;Y1=0;X=0; u=1; Y2=0; Y3=0;X2=0;0;0; x=0; h=0.1;t0=0; tf=2; t1=0; t2=0; t3=0;N=(tf-t0)/h;for i=1:N k1=a*X1+b; k2=b+a*(h*k1/2+X1); k3=b+a*(h*k2/2+X1); k4=b+a*(

5、h*k3+X1); X1=X1+h*(k1+2*k2+2*k3+k4)/6; Y1=Y1,c*X1; t1=t1,t1(i)+h;endfor i=1:N x=X2(:,i)+h*(a*X2(:,i)+b*u); y=c*x; X2=X2,x; Y2=Y2,y; t2=t2,t2(i)+h;endfor i=1:N y=1.84-4.95*i*exp(-1.88*i)-1.5*exp(-1.88*i)-0.34*exp(-6.24*i); Y3=Y3,y; t3=t3,t3(i)+h;endplot(t1,Y1,'r',t2,Y2,'g',t3,Y3,'

6、b') 當(dāng)h=0.01時(shí)，仿真圖如下： 可見當(dāng)步長很小時(shí)，數(shù)值積分解逼近解析解，其數(shù)值解得的精度也很高。 當(dāng)h=0.1時(shí)，仿真圖如下： 可見隨著步長的增加，其數(shù)值積分解與解析解的逼近程度減小，精度降低。 當(dāng)h=0.3時(shí)，仿真圖如下： 可見其數(shù)值積分解已變得不穩(wěn)定了，即步長不能取得太大3、求下圖所示系統(tǒng)的階躍響應(yīng)y(t)的數(shù)值解。（v=1,k=1,開始時(shí)間t0=0,結(jié)束時(shí)間tf=10,h=0.25）分析k、v對系統(tǒng)響應(yīng)的影響。（編程用RK4求解；Simulink）(修正：G(s)表達(dá)式的分母中，0.25s+1改為（0.25s+1）的二次方)    程序如下： k=1;

7、 a=conv(1 0 0,conv(0.25 1,0.25 1) b=2*k k; X0=0 0 0 0 ;v=1; n0=4; tf=10; h0=0.25; r=1; V=v; n=n0; T0=0; Tf=tf; h=h0; R=r; b=b/a(1); a=a/a(1); A=a(2: n+1); %首一化處理 A=rot90(rot90(eye(n-1,n); -fliplr(A); %形成能控標(biāo)準(zhǔn)型A陣 B=zeros(1,n-1),1' %形成輸入陣B(n維列向量)m1=length(b); %分子系數(shù)向量維數(shù)m+1 C=fliplr(b),zeros(1,n-m1);

8、 %形成輸出陣C(n維行向量) Ab=A-B*C*V; %形成閉環(huán)系數(shù)陣AbX=X0'y=0;t=T0; %賦初值，準(zhǔn)備開始遞推運(yùn)算 N=round(Tf-T0)/h; %確定輸出點(diǎn)數(shù)for i=1:N %四階龍格-庫塔法 k1=Ab*X+B*R; k2=Ab*(X+h*k1/2)+B*R; k3=Ab*(X+h*k2/2)+B*R; k4=Ab*(X+h*k3)+B*R; %求各階次斜率 X=X+h*(k1+2*k2+2*k3+k4)/6; %求狀態(tài) y=y,C*X; %求輸出并以向量形式保存 t=t,t(i)+h; %輸出對應(yīng)時(shí)刻以向量形式保存endt',y' %輸

9、出數(shù)據(jù)形式結(jié)果plot(t,y) %輸出曲線形式結(jié)果 k=0.5 K=1，v=1 K=1，V=0.5 K=1,v=5時(shí)， K=2，v=1時(shí)分析：當(dāng)k取值增大，v值不變時(shí)，系統(tǒng)輸出的波頭增多，而且也變陡，穩(wěn)態(tài)精度降低，當(dāng)k增加到一定程度時(shí)系統(tǒng)便發(fā)散了（即不穩(wěn)定了）。當(dāng)v值增大，k值不變時(shí)，波頭也是變多變陡，當(dāng)v值增大到一定程度時(shí)系統(tǒng)便不穩(wěn)定了。4、已知系統(tǒng)狀態(tài)狀態(tài)方程為（修正：x的系數(shù)矩陣中最后一個(gè)元素為-40，而非40）采用RK4法，步長分別取h=0.01、0.04、0.06，求系統(tǒng)的零輸入響應(yīng)，并繪圖分析各狀態(tài)變量的響應(yīng)狀態(tài)及產(chǎn)生的原因。（提示：病態(tài)系統(tǒng)）程序如下：a=-21 19 -20;

10、19 -21 20;40 -40 -40;x=1;0;-1;X=x;t=0;t0=0;tf=2;h=0.01;n=round(tf-t0)/h;for i=1:n x=X(:,i)+h*a*X(:,i); X=X,x; t=t,t(i)+h;endb=X(1,:);c=X(2,:);d=X(3,:);plot(t,b,'r',t,c,'g',t,d,'b') h=0.01時(shí) H=0.02時(shí) H=0.04時(shí)分析：如圖，當(dāng)h=0.01時(shí)，在t=0.2s以后系統(tǒng)輸出便趨于平穩(wěn)，當(dāng)取h=0.02時(shí)，系統(tǒng)輸出振蕩劇烈，趨于穩(wěn)定的時(shí)間也變長，當(dāng)取h=0.04

11、后，系統(tǒng)輸出呈發(fā)散振蕩形式。當(dāng)h=0.06后系統(tǒng)仍然是發(fā)散的，即當(dāng)h的取值改變時(shí)，原先穩(wěn)定的系統(tǒng)變得不穩(wěn)定了，這便是病態(tài)系統(tǒng)。但是這個(gè)結(jié)果與書上的不同。三、思考題1、在進(jìn)行仿真計(jì)算時(shí)，是否選用的數(shù)值積分法的階次越高越好？ 不是，因?yàn)殡A次越高，雖然精度高了，但其計(jì)算的復(fù)雜程度也越高，有時(shí)甚至無法完成計(jì)算。 2、選用數(shù)值積分法進(jìn)行仿真的原則。(1) 精度 仿真結(jié)果的精度主要受三項(xiàng)誤差影響： 截?cái)嗾`差：由算法本身的精度階次所決定； 舍入誤差：由計(jì)算機(jī)字長決定； 累積誤差：由以上兩項(xiàng)誤差隨計(jì)算時(shí)間長短累積情況決定。(2) 計(jì)算速度 取決于所用數(shù)值方法和計(jì)算步長。(3) 穩(wěn)定性 數(shù)值穩(wěn)定性主要

12、與計(jì)算步長h有關(guān)，不同的數(shù)值方法對h都有不同的穩(wěn)定性限制范圍，且與被仿真對象的時(shí)間常數(shù)有關(guān)。一般所選步長與系統(tǒng)最小時(shí)間常數(shù)有以下數(shù)量級(jí)關(guān)系： 而多步法、隱式算法有較好的數(shù)值穩(wěn)定性，在對穩(wěn)定性較注重時(shí)，應(yīng)予以優(yōu)先選擇。 實(shí)驗(yàn)二第一部分  面向結(jié)構(gòu)圖的數(shù)值積分法仿真一、實(shí)驗(yàn)?zāi)康募由罾斫膺B續(xù)系統(tǒng)面向結(jié)構(gòu)圖仿真的原理及特點(diǎn)，進(jìn)一步掌握數(shù)字積分法解算微分方程的方法，增加編寫仿真程序的能力。二、實(shí)驗(yàn)內(nèi)容1、用面向方框圖的數(shù)字仿真方法對下列系統(tǒng)進(jìn)行仿真。    2、求解下圖所示系統(tǒng)在f=-1(t)階躍擾動(dòng)作用下第、第環(huán)節(jié)的動(dòng)態(tài)過程。分別用面向框圖的

13、數(shù)值積分法（RK4法）、MATLAB中有關(guān)系統(tǒng)建模的命令和Simulink三種方法求解。三實(shí)驗(yàn)程序及結(jié)果1.P=0,0.07,1,0.14;1,0.012,1,0;0,0.05,1,0.15;10,1,1,0;1,0.01,0,0.0008;WIJ=1,0,1;1,4,-1;2,1,1;3,2,1;3,5,-1;4,3,1;5,4,1n=5;Y0=1;Yt0=0 0 0 0 0;h=0.01;t=0;L1=5;T0=0;Tf=2;nout=4;A=diag(P(:,1);B=diag(P(:,2);C=diag(P(:,3);D=diag(P(:,4);m=length(WIJ(:,1);W0

14、=zeros(n,1);W=zeros(n,n);for k=1:m if (WIJ(k,2)=0) W0(WIJ(k,1)=WIJ(k,3); else W(WIJ(k,1),WIJ(k,2)=WIJ(k,3); end;end;Q=B-D*W;Qn=inv(Q);R=C*W-A;V1=C*W0;Ab=Qn*R;b1=Qn*V1;Y=Yt0'y=Y(nout);t=T0;N=round(Tf-T0)/(h*L1);for i=1:N; for j=1:L1; k1=Ab*Y+b1*1; k2=Ab*(Y+h*k1/2)+b1*1; k3=Ab*(Y+h*k2/2)+b1*1; k4=

15、Ab*(Y+h*k3)+b1*1; Y=Y+h*(k1+2*k2+2*k3+k4)/6; end; y=y,Y(nout); t=t,t(i)+L1*h; end; plot(t,y)2.P=1 0.0149 1 0; 1 0.00254 26.6667 0; 1 0.391 0 0.4199; 0 0.248 1 0; 1 2.88 1 0;WIJ=1 5 1; 2 1 1; 3 4 1; 4 2 1; 4 3 -1; 5 0 1; 5 4 -1;n=5;Y0=-1;Yt0=0 0 0 0 0;h=0.005;L1=10;T0=0;Tf=10;nout=5;A=diag(P(:,1);B=d

16、iag(P(:,2);C=diag(P(:,3);D=diag(P(:,4);m=length(WIJ(:,1);W0=zeros(n,1);W=zeros(n,n);for k=1:m if (WIJ(k,2)=0);W0(WIJ(k,1)=WIJ(k,3); else W(WIJ(k,1),WIJ(k,2)=WIJ(k,3); end;end;Q=B-D*W;Qn=inv(Q);R=C*W-A;V1=C*W0;Ab=Qn*R;b1=Qn*V1;Y=Yt0'y=Y(nout);t=T0;N=round(Tf-T0)/(h*L1);for i=1:N; for j=1:L1; K1=A

17、b*Y+b1*Y0; K2=Ab*(Y+h*K1/2)+b1*Y0; K3=Ab*(Y+h*K2/2)+b1*Y0; K4=Ab*(Y+h*K3)+b1*Y0; Y=Y+h*(K1+2*K2+2*K3+K4)/6; end;y=y,Y(nout);t=t,t(i)+h*L1;end;t',y'plot(t,y)grid on;xlabel('t/s');ylabel('Amplitude');axis(0 10 -0.08 0.02);結(jié)果：第二部分  面向結(jié)構(gòu)圖的離散相似法仿真一、實(shí)驗(yàn)?zāi)康?#160;  1掌握離散

18、相似法的基本概念和原理；2典型環(huán)節(jié)離散系數(shù)的求取及差分方程表示；3. 掌握非線性系統(tǒng)的數(shù)字仿真方法。二、實(shí)驗(yàn)內(nèi)容1、已知控制系統(tǒng)結(jié)構(gòu)圖如圖所示，設(shè)輸入階躍函數(shù)幅值Y0=10，滯環(huán)非線性參數(shù)s=1(滯環(huán)寬度），請用離散相似法編程和Simulink法對系統(tǒng)進(jìn)行如下分析：1）不考慮非線性環(huán)節(jié)影響時(shí)，求解y(t)的階躍響應(yīng)；2）考慮非線性環(huán)節(jié)影響，其余參數(shù)不變，求解y(t)并與線性情況所得結(jié)果進(jìn)行比較；3）改變的滯環(huán)非線性參數(shù)s，分析該非線性對系統(tǒng)的影響。2、系統(tǒng)結(jié)構(gòu)圖如圖所示，先理論分析該系統(tǒng)是否會(huì)產(chǎn)生自激振蕩，若會(huì)求出振蕩的振幅和頻率，并用離散相似法編程和Simulink法驗(yàn)證分析的結(jié)果

19、。三實(shí)驗(yàn)程序及結(jié)果1.(1)Simulink法結(jié)果：(2)MATLAB法P=1 10 5 10;1 0.5 1 0;1 0.1 1 0;0 1 1 0;WIJ=1 0 1;2 1 1;3 2 1;4 3 1;1 4 -1;X0=0 0 0 0;Z=0 0 0 6;S=0 0 0 1; h=0.01;L1=25;n=4;T0=0;Tf=10;nout=4;Y0=10;A=(P(:,1);B=(P(:,2);C=(P(:,3);D=(P(:,4);m=length(WIJ(:,1);W0=zeros(n,1);W=zeros(n,n);for k=1:m if (WIJ(k,2)=0);W0(WI

20、J(k,1)=WIJ(k,3); else W(WIJ(k,1),WIJ(k,2)=WIJ(k,3); endendfor i=1:n if(A(i)=0); FI(i)=1; FIM(i)=h*C(i)/B(i); FIJ(i)=h*h*C(i)/B(i)/2; FIC(i)=1;FID(i)=0; if(D(i)=0); FID(i)=D(i)/B(i); else end else FI(i)=exp(-h*A(i)/B(i); FIM(i)=(1-FI(i)*C(i)/A(i); FIJ(i)=h*C(i)/A(i)-FIM(i)*B(i)/A(i); FIC(i)=1;FID(i)=

21、0; if(D(i)=0); FIM=(1-FI(i)*D(i)/A(i); FIJ(i)=h*D(i)/A(i)-FIM(i)*B(i)/A(i); FIC(i)=C(i)/D(i)-A(i)/B(i); FID(i)=D(i)/B(i); else end endend Y=zeros(n,1);X=Y;y=0;Uk=zeros(n,1);Ubb=Uk;t=T0:h*L1:Tf;N=length(t);for k=1:N-1 for l=1:L1 Ub=Uk; Uk=W*Y+W0*Y0; for i=1:n if(Z(i)=0) if(Z(i)=1) Uk(i)=satu(Uk(i),S(

22、i); end if(Z(i)=2) Uk(i)=dead(Uk(i),S(i); end if(Z(i)=3) Uk(i),Ubb(i)=backlash(Ubb(i),Uk(i),Ub(i),S(i); end end end Udot=(Uk-Ub)/h; Uf=2*Uk-Ub; X=FI'.*X+FIM'.*Uk+FIJ'.*Udot; Yb=Y; Y=FIC'.*X+FID'.*Uf; for i=1:n if(Z(i)=0) if(Z(i)=4) Y(i)=satu(Y(i),S(i); end if(Z(i)=5) Y(i)=dead(Y(

23、i),S(i); end if(Z(i)=6) Y(i),Ubb(i)=backlash(Ubb(i),Y(i),Yb(i),S(i); end end end end y=y,Y(nout);endplot(t,y)結(jié)果：S4=1S4=52. (1)Simulink法： (2)Matlab法P=0,1,1,0;1,1,1,0;2,1,1,0;WIJ=1 0 1;1 3 -1;2 1 1;3 2 1;n=3;Y0=10;X0=0 0 0 0;Yt0=0 0 0 0;h=0.01;L1=25;t0=0;tf=10;nout=3;Z=0 0 6;S=0 0 1;A=(P(:,1);B=(P(:,2

24、);C=(P(:,3);D=(P(:,4);m=length(WIJ(:,1);W0=zeros(n,1);W=zeros(n,n);for k=1:m if (WIJ(k,2)=0);W0(WIJ(k,1)=WIJ(k,3); else W(WIJ(k,1),WIJ(k,2)=WIJ(k,3); end;end;sp3_3sp3_4plot(t,y)實(shí)驗(yàn)三 采樣控制系統(tǒng)的數(shù)字仿真一、實(shí)驗(yàn)?zāi)康?#160; 1、掌握采樣控制系統(tǒng)數(shù)字仿真的一般方法；2、掌握采樣周期與計(jì)算步長的關(guān)系與確定方法。二、實(shí)驗(yàn)內(nèi)容1、已知采樣系統(tǒng)結(jié)構(gòu)如圖所示，分別用程序設(shè)計(jì)方法和Simulink法求系統(tǒng)的輸出響應(yīng)

25、。（取仿真時(shí)間:Tf=10；采樣周期:T=0.2；計(jì)算步長：h=0.01）(注：第一個(gè)框圖中，分子分母中z的指數(shù)均為-1；第二個(gè)框圖中e的指數(shù)為-Ts)2、設(shè)某數(shù)字控制系統(tǒng)如圖所示，采樣周期為T=0.1s，初始狀態(tài)。(注：第二個(gè)框圖中e的指數(shù)為-Ts) （1）數(shù)字控制器為PI調(diào)節(jié)器 ，其中（2）數(shù)字控制器為PID調(diào)節(jié)器，其中要求： 1）用程序設(shè)計(jì)法、MATLAB控制工具箱時(shí)域響應(yīng)分析函數(shù)、Simulink法求系統(tǒng)在單位階躍輸入信號(hào)作用下的輸出響應(yīng)；  2）比較兩種控制器作用下系統(tǒng)的響應(yīng)，由此得出微分調(diào)節(jié)的作用。三，程序及結(jié)果1.（1）matlab法：G=2.72

26、-1;F=0.717;P=0 1 1 0; 1 1 1 0;WIJ=1 0 1; 2 1 1;n=2;Y0=1;X0=0 0; h=0.01;T0=0;Tf=10;Ts=0.2;nout=2;A=P(:,1);B=P(:,2);C=P(:,3);D=P(:,4);m=length(WIJ(:,1);W0=zeros(n,1);W=zeros(n,n); for k=1:m if(WIJ(k,2)=0);W0(WIJ(k,1)=WIJ(k,3); else W(WIJ(k,1),WIJ(k,2)=WIJ(k,3); endend for i=1:n if(A(i)=0) FI(i)=1; FIM

27、(i)=h*(C(i)/B(i); FIJ(i)=h*h*C(i)/B(i)/2; FIC(i)=1;FID(i)=0; if(D(i)=0) FID(i)=D(i)/B(i); else end else FI(i)=exp(-h)*A(i)/B(i); FIM(i)=(1-FI(i)*(C(i)/A(i); FIJ(i)=h*C(i)/A(i)-FIM(i)*B(i)/A(i); FIC(i)=1;FID(i)=0; if(D(i)=0) FIM(i)=(1-FI(i)*(D(i)/A(i); FIJ(i)=h*D(i)/A(i)-FIM(i)*B(i)*A(i); FIC(i)=C(i)

28、/D(i)-A(i)/B(i); FID(i)=D(i)/B(i); else end endendY=zeros(n,1);X=Y;y=0;Uk=zeros(n,1);Ub=Uk;U=0;uk=0;ek=0;E=zeros(n,1);x2=0;t=T0:h:Tf;t0=T0:h:Ts;N=length(t0);t1=T0:Ts:Tf;M=length(t1);for k=1:M-1; U=uk; ek=Y0-x2; E=ek;E(1); uk=-F*U+G*E; for l=1:N-1 Ub=Uk; Uk=W*Y+W0*uk; Udot=(Uk-Ub)/h; Uf=2*Uk-Ub; X=FI

29、'.*X+FIM'.*Uk+FIJ'.*Udot; Yb=Y; Y=FIC'.*X+FID'.*Uf; y(N-1)*(k-1)+l),1)=Y(nout); end x2=Y(nout); endplot(1:1000,y);結(jié)果：2.（1）G=15.2 -15;F=-1;P=1 1 1 0; 1 4 1 0;WIJ=1 0 1; 2 1 1;n=2;Y0=1;X0=0 0; h=0.01;L1=25;T0=0;Tf=10;Ts=0.1;nout=2;A=P(:,1);B=P(:,2);C=P(:,3);D=P(:,4);m=length(WIJ(:,

30、1);W0=zeros(n,1);W=zeros(n,n); for k=1:m if(WIJ(k,2)=0);W0(WIJ(k,1)=WIJ(k,3); else W(WIJ(k,1),WIJ(k,2)=WIJ(k,3); endend for i=1:n if(A(i)=0) FI(i)=1; FIM(i)=h*(C(i)/B(i); FIJ(i)=h*h*C(i)/B(i)/2; FIC(i)=1;FID(i)=0; if(D(i)=0) FID(i)=D(i)/B(i); else end else FI(i)=exp(-h)*A(i)/B(i); FIM(i)=(1-FI(i)*(C

31、(i)/A(i); FIJ(i)=h*C(i)/A(i)-FIM(i)*B(i)/A(i); FIC(i)=1;FID(i)=0; if(D(i)=0) FIM(i)=(1-FI(i)*(D(i)/A(i); FIJ(i)=h*D(i)/A(i)-FIM(i)*B(i)*A(i); FIC(i)=C(i)/D(i)-A(i)/B(i); FID(i)=D(i)/B(i); else end endendY=zeros(n,1);X=Y;y=0;Uk=zeros(n,1);Ub=Uk;U=0;uk=0;ek=0;E=zeros(2,1);x2=0;t=T0:h:Tf;t0=T0:h:Ts;N=l

32、ength(t0);t1=T0:Ts:Tf;M=length(t1);for k=1:M-1; U=uk; ek=Y0-x2; E=ek;E(1); uk=-F*U+G*E; for l=1:N-1 Ub=Uk; Uk=W*Y+W0*uk; Udot=(Uk-Ub)/h; Uf=2*Uk-Ub; X=FI'.*X+FIM'.*Uk+FIJ'.*Udot; Yb=Y; Y=FIC'.*X+FID'.*Uf; y(N-1)*(k-1)+l),1)=Y(nout); end x2=Y(nout); endt=0:0.01:9.99;plot(t,y,'

33、-');grid on;xlabel('t/s');ylabel('Amplitude');hold on（2）matlab法G=55.2 -95 40;F=-1;P=1 1 1 0; 1 4 1 0;WIJ=1 0 1; 2 1 1;n=2;Y0=1;X0=0 0; h=0.01;L1=25;T0=0;Tf=10;Ts=0.1;nout=2;A=P(:,1);B=P(:,2);C=P(:,3);D=P(:,4);m=length(WIJ(:,1);W0=zeros(n,1);W=zeros(n,n); for k=1:m if(WIJ(k,2)=0);

34、W0(WIJ(k,1)=WIJ(k,3); else W(WIJ(k,1),WIJ(k,2)=WIJ(k,3); endend for i=1:n if(A(i)=0) FI(i)=1; FIM(i)=h*(C(i)/B(i); FIJ(i)=h*h*C(i)/B(i)/2; FIC(i)=1;FID(i)=0; if(D(i)=0) FID(i)=D(i)/B(i); else end else FI(i)=exp(-h)*A(i)/B(i); FIM(i)=(1-FI(i)*(C(i)/A(i); FIJ(i)=h*C(i)/A(i)-FIM(i)*B(i)/A(i); FIC(i)=1;

35、FID(i)=0; if(D(i)=0) FIM(i)=(1-FI(i)*(D(i)/A(i); FIJ(i)=h*D(i)/A(i)-FIM(i)*B(i)*A(i); FIC(i)=C(i)/D(i)-A(i)/B(i); FID(i)=D(i)/B(i); else end endendY=zeros(n,1);X=Y;y=0;Uk=zeros(n,1);Ub=Uk;U=0;uk=0;ek=0;E=zeros(3,1);x2=0;t=T0:h:Tf;t0=T0:h:Ts;N=length(t0);t1=T0:Ts:Tf;M=length(t1);for k=1:M-1; U=uk; ek

36、=Y0-x2; E=ek;E(1:2); uk=-F*U+G*E; for l=1:N-1 Ub=Uk; Uk=W*Y+W0*uk; Udot=(Uk-Ub)/h; Uf=2*Uk-Ub; X=FI'.*X+FIM'.*Uk+FIJ'.*Udot; Yb=Y; Y=FIC'.*X+FID'.*Uf; y(N-1)*(k-1)+l),1)=Y(nout); end x2=Y(nout); endt=0:0.01:9.99;plot(t,y);grid on;xlabel('t/s');ylabel('Amplitude');

37、實(shí)驗(yàn)四  數(shù)字仿真技術(shù)的綜合應(yīng)用一、 實(shí)驗(yàn)?zāi)康氖箤W(xué)生在系統(tǒng)地學(xué)習(xí)了計(jì)算機(jī)仿真的基本概念和方法后，對計(jì)算機(jī)仿真的全過程有一個(gè)清楚的概念和認(rèn)識(shí)，能夠?qū)⑺鶎W(xué)的知識(shí)應(yīng)用到實(shí)際系統(tǒng)中去。二、實(shí)驗(yàn)內(nèi)容    選用經(jīng)典控制理論或現(xiàn)代控制理論中某一種你所熟知的控制器設(shè)計(jì)方法，為下列某一系統(tǒng)（三選一）進(jìn)行控制器設(shè)計(jì)：  （三）直流電動(dòng)機(jī)調(diào)速系統(tǒng) 1 系統(tǒng)物理模型、參數(shù)及要求1.1 系統(tǒng)物理模型本組設(shè)計(jì)的被控對象為前一課程設(shè)計(jì)直流拖動(dòng)自動(dòng)控制系統(tǒng)所采用的直流電機(jī)，其物理模型如下圖1所示。圖1 雙閉環(huán)控制電流調(diào)速系統(tǒng)的特點(diǎn)是電機(jī)的轉(zhuǎn)速和電流分別由兩個(gè)獨(dú)立的調(diào)節(jié)器分別

38、控制，且轉(zhuǎn)速調(diào)節(jié)器的輸出就是電流調(diào)節(jié)器的給定，因此電流環(huán)能夠跟隨轉(zhuǎn)速的偏差調(diào)節(jié)電機(jī)電樞的電流。當(dāng)轉(zhuǎn)速低于給定轉(zhuǎn)速時(shí)，轉(zhuǎn)速調(diào)節(jié)器的積分作用使輸出增加，即電流給定上升，并通過電流環(huán)調(diào)節(jié)使電機(jī)電樞電流增大，從而使電機(jī)獲得加速轉(zhuǎn)矩，電機(jī)轉(zhuǎn)速上升。當(dāng)實(shí)際轉(zhuǎn)速高于給定轉(zhuǎn)速時(shí)，轉(zhuǎn)速調(diào)節(jié)器的輸出減小，即給定電流減小，并通過電流環(huán)調(diào)節(jié)使電機(jī)電樞電流下降，電機(jī)將因?yàn)殡姶呸D(zhuǎn)矩減小而減速。在當(dāng)轉(zhuǎn)速調(diào)節(jié)器飽和輸出達(dá)到限幅值時(shí)，電流環(huán)即以最大電流限制實(shí)現(xiàn)電機(jī)的加速，使電機(jī)的啟動(dòng)時(shí)間最短。雙閉環(huán)調(diào)速系統(tǒng)的原理框圖如圖2所示。圖21.2 系統(tǒng)參數(shù)電機(jī)型號(hào)：DJ15minrad=n /1600 額定參數(shù)： ， ， ， ， 電樞

39、電阻：R=29，電樞電感：L=0.89mH電機(jī)飛輪慣量：GD=0.045N/m電樞回路電磁時(shí)間常數(shù)：T=0.03s, 系統(tǒng)的機(jī)電時(shí)間常數(shù)：T=0.034s電動(dòng)機(jī)電勢時(shí)間常數(shù)：C=0.105， 轉(zhuǎn)矩常數(shù)：C=1.0Vmin/r電流反饋系數(shù)：=5V/A, 轉(zhuǎn)速反饋系數(shù)：=0.003V/(rpm) 1.3 要求設(shè)計(jì)狀態(tài)反饋控制器，使得系統(tǒng)的單位階躍響應(yīng)性能指標(biāo)為：（1）調(diào)節(jié)時(shí)間小于2秒； （2）系統(tǒng)的超調(diào)量小于5%； （3）穩(wěn)態(tài)誤差小于1%。2 系統(tǒng)模型的建立2.1 模型抽象直流電機(jī)轉(zhuǎn)矩和電樞電流的關(guān)系為： 電樞旋轉(zhuǎn)產(chǎn)生反電動(dòng)勢e與旋轉(zhuǎn)運(yùn)動(dòng)角速度的關(guān)系為： 根據(jù)牛頓第二定律列寫運(yùn)動(dòng)平衡方程式為： 其

40、中b為電機(jī)摩擦系數(shù)，此處忽略不計(jì)。根據(jù)回路電壓法列寫電機(jī)電樞回路方程式為： 由于： ，可得： ， 其中，m為一個(gè)旋轉(zhuǎn)體上的一個(gè)質(zhì)點(diǎn)的質(zhì)量，質(zhì)量m為該質(zhì)量的重量G和重力加速度g之比，R和D分別為旋轉(zhuǎn)體的半徑和直徑，綜合上兩式可得： 從而可以得到電機(jī)電樞回路電壓平衡和電機(jī)運(yùn)動(dòng)平衡的一組微分方程式 其中，摩擦系數(shù) =b/9.55,此處忽略不計(jì)。設(shè)系統(tǒng)的狀態(tài)變量為： ，以輸入電壓u為輸入，轉(zhuǎn)速n為輸出。建立系統(tǒng)狀態(tài)空間表達(dá)式為： ïïîïïíìúûùêëé=ú&#

41、251;ùêëé+úûùêëéúûùêëé=úûùêëé21212110012.10118.0-3.833358.32-xxyxxxx&&u由式帶入數(shù)據(jù)進(jìn)行計(jì)算化簡可得：其中u代表給定電壓值，代表系統(tǒng)內(nèi)部電流變化，代表系統(tǒng)輸出轉(zhuǎn)速變化，y=x表示系統(tǒng)的轉(zhuǎn)速輸出。可得：A=-32.58 -0.118;8333.33 0; B=1.12;0; C=0 1; D=0。

42、2.2 所建模型的性能分析2.2.1 能空能觀測性分析(1) 能空性分析MATLAB程序如下： A=-32.58 -0.118;8333.33 0;B=1.12;0;C=0 1;D=0;n=length(A) %求系統(tǒng)階次qc=ctrb(A,B),nc=rank(qc) %求能控性判別矩陣及其秩if n=nc,disp('系統(tǒng)是可控的!'), %判斷系統(tǒng)的能控性else disp('系統(tǒng)是不可控的!'),end運(yùn)行結(jié)果為：>> 原系統(tǒng)的秩為：n = 2qc = 1.0e+003 * 0.0011 -0.0365 0 9.3333nc = 2系統(tǒng)是可控

43、的! (2) 能觀測性分析MATLAB程序如下：A=-32.58 -0.118;8333.33 0;B=1.12;0;C=0 1;D=0;disp('原系統(tǒng)的秩為：');n=length(A) %求系統(tǒng)階次qo=obsv(A,C),no=rank(qo) %求系統(tǒng)能觀測性判別矩陣及其秩if n=no,disp('系統(tǒng)是可觀測的!'), %判斷系統(tǒng)的能觀測性else disp('系統(tǒng)是不可觀測的!'),end運(yùn)行結(jié)果為：>>原系統(tǒng)的秩為：n = 2qo = 1.0e+003 * 0 0.0010 8.3333 0no = 2系統(tǒng)是可觀測

44、的!2.2.2 系統(tǒng)穩(wěn)定性分析(1) 利用MATLAB程序分析系統(tǒng)穩(wěn)定性是系統(tǒng)正常工作的首要條件。只要系統(tǒng)的狀態(tài)矩陣A的特征根全部具有負(fù)實(shí)部，系統(tǒng)就是狀態(tài)穩(wěn)定的。MATLAB程序如下:A=-32.58 -0.118;8333.33 0; B=1.12;0; C=0 1; D=0;z,p,k=ss2zp(A,B,C,D,1); Flags=0; n=length(A);for i=1:n; if real(p(i)>0; Flags=1; endenddisp('系統(tǒng)的零極點(diǎn)模型為:')z,p,kif Flags=1 disp('系統(tǒng)不穩(wěn)定');else disp('系統(tǒng)是穩(wěn)定的');end運(yùn)行結(jié)果為：>> 系統(tǒng)的零極點(diǎn)模型為:z = Empty matrix: 0-by-1p = -16.2900+26.7949i -16.2900-26.7949ik = 9.3333e+00系統(tǒng)是穩(wěn)定的(2) 利用零極點(diǎn)圖分析MATLAB程序如下:A=-32.58 -0.118;8333.33 0;B=1.12;0;C