1、實驗一  面向微分方程的數值積分法仿真一、實驗目的1掌握數值積分法的基本概念、原理及應用；2用龍格-庫塔法解算微分方程，增加編寫仿真程序的能力；3分析數值積分算法的計算步長與計算精度、速度、穩定性的關系；4. 對數值算法中的“病態問題”進行研究。二、實驗內容1、已知系統微分方程及初值條件取步長h=0.1 ，試分別用歐拉方程法和RK4法求t=2h時的y值，并將求得的值與解析解比較(將三個解繪于同一坐標中，且用數值進行比較)，說明造成差異的原因。（編程完成；選用MATLAB ode函數完成。） 程序如下：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')分析:紅線為用歐拉法得到的結果，綠線為用四階龍格庫塔法得到的結果，藍線為根據解析方程得到的結果。其差異原因主要有兩個：（1）二者的方法不同，歐拉法是根據一階微分方程計算得到的，龍格庫塔法是根據四階微分方程得到的；（2）由于步長取為0.1，所以得到的圖像與解析解之間存在差異，若將步長取小，則得到的解將更靠近解析解。  2、已知系統的傳遞函數為在單位階躍輸入下，系統響應的解析解為 試分別用歐拉方程法和RK4法對系統進行仿真（

4、編程完成）：1）比較兩種數值積分解與解析解的逼近程度；（繪圖）2）改變步長，分析步長對數值解精度的影響；3）不斷加大步長，分析計算穩定性的變化。程序如下： 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') 當h=0.01時，仿真圖如下： 可見當步長很小時，數值積分解逼近解析解，其數值解得的精度也很高。 當h=0.1時，仿真圖如下： 可見隨著步長的增加，其數值積分解與解析解的逼近程度減小，精度降低。 當h=0.3時，仿真圖如下： 可見其數值積分解已變得不穩定了，即步長不能取得太大3、求下圖所示系統的階躍響應y(t)的數值解。（v=1,k=1,開始時間t0=0,結束時間tf=10,h=0.25）分析k、v對系統響應的影響。（編程用RK4求解；Simulink）(修正：G(s)表達式的分母中，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); %形成能控標準型A陣 B=zeros(1,n-1),1' %形成輸入陣B(n維列向量)m1=length(b); %分子系數向量維數m+1 C=fliplr(b),zeros(1,n-m1);

8、 %形成輸出陣C(n維行向量) Ab=A-B*C*V; %形成閉環系數陣AbX=X0'y=0;t=T0; %賦初值，準備開始遞推運算 N=round(Tf-T0)/h; %確定輸出點數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; %求狀態 y=y,C*X; %求輸出并以向量形式保存 t=t,t(i)+h; %輸出對應時刻以向量形式保存endt',y' %輸

9、出數據形式結果plot(t,y) %輸出曲線形式結果 k=0.5 K=1，v=1 K=1，V=0.5 K=1,v=5時， K=2，v=1時分析：當k取值增大，v值不變時，系統輸出的波頭增多，而且也變陡，穩態精度降低，當k增加到一定程度時系統便發散了（即不穩定了）。當v值增大，k值不變時，波頭也是變多變陡，當v值增大到一定程度時系統便不穩定了。4、已知系統狀態狀態方程為（修正：x的系數矩陣中最后一個元素為-40，而非40）采用RK4法，步長分別取h=0.01、0.04、0.06，求系統的零輸入響應，并繪圖分析各狀態變量的響應狀態及產生的原因。（提示：病態系統）程序如下：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時 H=0.02時 H=0.04時分析：如圖，當h=0.01時，在t=0.2s以后系統輸出便趨于平穩，當取h=0.02時，系統輸出振蕩劇烈，趨于穩定的時間也變長，當取h=0.04

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

12、與計算步長h有關，不同的數值方法對h都有不同的穩定性限制范圍，且與被仿真對象的時間常數有關。一般所選步長與系統最小時間常數有以下數量級關系： 而多步法、隱式算法有較好的數值穩定性，在對穩定性較注重時，應予以優先選擇。 實驗二第一部分  面向結構圖的數值積分法仿真一、實驗目的加深理解連續系統面向結構圖仿真的原理及特點，進一步掌握數字積分法解算微分方程的方法，增加編寫仿真程序的能力。二、實驗內容1、用面向方框圖的數字仿真方法對下列系統進行仿真。    2、求解下圖所示系統在f=-1(t)階躍擾動作用下第、第環節的動態過程。分別用面向框圖的

13、數值積分法（RK4法）、MATLAB中有關系統建模的命令和Simulink三種方法求解。三實驗程序及結果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);結果：第二部分  面向結構圖的離散相似法仿真一、實驗目的   1掌握離散

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

19、。三實驗程序及結果1.(1)Simulink法結果：(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)結果：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)實驗三 采樣控制系統的數字仿真一、實驗目的  1、掌握采樣控制系統數字仿真的一般方法；2、掌握采樣周期與計算步長的關系與確定方法。二、實驗內容1、已知采樣系統結構如圖所示，分別用程序設計方法和Simulink法求系統的輸出響應

25、。（取仿真時間:Tf=10；采樣周期:T=0.2；計算步長：h=0.01）(注：第一個框圖中，分子分母中z的指數均為-1；第二個框圖中e的指數為-Ts)2、設某數字控制系統如圖所示，采樣周期為T=0.1s，初始狀態。(注：第二個框圖中e的指數為-Ts) （1）數字控制器為PI調節器 ，其中（2）數字控制器為PID調節器，其中要求： 1）用程序設計法、MATLAB控制工具箱時域響應分析函數、Simulink法求系統在單位階躍輸入信號作用下的輸出響應；  2）比較兩種控制器作用下系統的響應，由此得出微分調節的作用。三，程序及結果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);結果：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、實驗四  數字仿真技術的綜合應用一、 實驗目的使學生在系統地學習了計算機仿真的基本概念和方法后，對計算機仿真的全過程有一個清楚的概念和認識，能夠將所學的知識應用到實際系統中去。二、實驗內容    選用經典控制理論或現代控制理論中某一種你所熟知的控制器設計方法，為下列某一系統（三選一）進行控制器設計：  （三）直流電動機調速系統 1 系統物理模型、參數及要求1.1 系統物理模型本組設計的被控對象為前一課程設計直流拖動自動控制系統所采用的直流電機，其物理模型如下圖1所示。圖1 雙閉環控制電流調速系統的特點是電機的轉速和電流分別由兩個獨立的調節器分別

38、控制，且轉速調節器的輸出就是電流調節器的給定，因此電流環能夠跟隨轉速的偏差調節電機電樞的電流。當轉速低于給定轉速時，轉速調節器的積分作用使輸出增加，即電流給定上升，并通過電流環調節使電機電樞電流增大，從而使電機獲得加速轉矩，電機轉速上升。當實際轉速高于給定轉速時，轉速調節器的輸出減小，即給定電流減小，并通過電流環調節使電機電樞電流下降，電機將因為電磁轉矩減小而減速。在當轉速調節器飽和輸出達到限幅值時，電流環即以最大電流限制實現電機的加速，使電機的啟動時間最短。雙閉環調速系統的原理框圖如圖2所示。圖21.2 系統參數電機型號：DJ15minrad=n /1600 額定參數： ， ， ， ， 電樞

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

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

41、251;ùêëé+úûùêëéúûùêëé=úûùêëé21212110012.10118.0-3.833358.32-xxyxxxx&&u由式帶入數據進行計算化簡可得：其中u代表給定電壓值，代表系統內部電流變化，代表系統輸出轉速變化，y=x表示系統的轉速輸出?？傻茫篈=-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) %求系統階次qc=ctrb(A,B),nc=rank(qc) %求能控性判別矩陣及其秩if n=nc,disp('系統是可控的!'), %判斷系統的能控性else disp('系統是不可控的!'),end運行結果為：>> 原系統的秩為：n = 2qc = 1.0e+003 * 0.0011 -0.0365 0 9.3333nc = 2系統是可控

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

44、的!2.2.2 系統穩定性分析(1) 利用MATLAB程序分析系統穩定性是系統正常工作的首要條件。只要系統的狀態矩陣A的特征根全部具有負實部，系統就是狀態穩定的。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('系統的零極點模型為:')z,p,kif Flags=1 disp('系統不穩定');else disp('系統是穩定的');end運行結果為：>> 系統的零極點模型為:z = Empty matrix: 0-by-1p = -16.2900+26.7949i -16.2900-26.7949ik = 9.3333e+00系統是穩定的(2) 利用零極點圖分析MATLAB程序如下:A=-32.58 -0.118;8333.33 0;B=1.12;0;C