《控制系統(tǒng)計算機(jī)仿真》上機(jī)實驗報告

1、蘭 州 理 工 大 學(xué)?控制系統(tǒng)計算機(jī)仿真?上機(jī)實驗報告院系： 電信學(xué)院 班級： 10級自動化五班 姓名： 學(xué)號： 時間： 2021 年 5 月 12 日電氣工程與信息工程學(xué)院2-2.用MATLAB語言求以下系統(tǒng)的狀態(tài)方程、傳遞函數(shù)、零極點增益和局部分式形式的模型參數(shù)，并分別寫出其相應(yīng)的數(shù)學(xué)模型表達(dá)式：12 解：1傳遞函數(shù)轉(zhuǎn)化為狀態(tài)方程： num=1,7,24,24num = 1 7 24 24den=1,10,35,50,24den = 1 10 35 50 24A,B,C,D=tf2ss(num,den)A = -10 -35 -50 -24 1 0 0 0 0 1 0 0 0 0 1 0

2、B = 1 0 0 0C = 1 7 24 24D = 0G2=ss(A,B,C,D) a = x1 x2 x3 x4 x1 -10 -35 -50 -24 x2 1 0 0 0 x3 0 1 0 0 x4 0 0 1 0 b = u1 x1 1 x2 0 x3 0 x4 0 c = x1 x2 x3 x4 y1 1 7 24 24 d = u1 y1 0狀態(tài)方程為： 傳遞函數(shù)轉(zhuǎn)換為零極點增益：>> num=7,24,24num = 7 24 24den=10,35,50,24den = 10 35 50 24Z,P,K=tf2zp(num,den)Z = -1.5388 P =K

3、 = 1 G1=zpk(Z,P,K) Zero/pole/gain:(s+1.539) (s2 + 5.461s + 15.6)- (s+4) (s+3) (s+2) (s+1)零極點增益方程為：傳遞函數(shù)轉(zhuǎn)換為局部分時形式：>> num=7,24,24num = 7 24 24den=10,35,50,24 den = 10 35 50 24 R,P,H=residue(num,den)R =P =H = G3=residue(R,P,H)G3 =局部分式形式為：2 解：狀態(tài)方程轉(zhuǎn)換為傳遞函數(shù)為：>> A=2.25,-5,-1.25,-0.5;2.25,-4.25,-1

4、.25,-0.5;0.25,-0.5,-1.25,-1;1.25,-1.75,-0.25,-0.75A = 1B=4,2,2,0'B = 4 2 2 0C=0,2,0,2C = 0 2 0 2D=0D = 0num,den=ss2tf(A,B,C,D)num =den =G1=tf(num,den) Transfer function: 4 s3 + 14 s2 + 19.75 s + 13-傳遞函數(shù)為：狀態(tài)方程轉(zhuǎn)換成零極點：>> A=2.25,-5,-1.25,-0.5;2.25,-4.25,-1.25,-0.5;0.25,-0.5,-1.25,-1;1.25,-1.75

5、,-0.25,-0.75A = 0B=4,2,2,0'B = 4 2 2 0C=0,2,0,2C = 0 2 0 2D=0D = 0Z,P,K=ss2zp(A,B,C,D)Z = -1.7331 P = -1.6364 -1.5000 K =G2=zpk(Z,P,K) Zero/pole/gain:4 (s+1.733) (s2 + 1.767s + 1.875)-(s+1.5) (s+1.636) (s2 + 0.8636s + 0.6493)零極點增益方程為：3轉(zhuǎn)換成局部分式形式：>> R,P,H=residue(num,den)R = -2.4618 6.2857 P

6、 = -1.6364 -1.5000 -0.4318 H = G3=residue(R,P,H)G3 =局部分式形式的方程為：2-3. 用毆拉法求以下系統(tǒng)的輸出響應(yīng)在上，時的數(shù)值解。，要求保存4位小數(shù)，并將結(jié)果與真解比擬。解：(1). >> h=0.1;disp('y=');y=1;for t=0:h:1m=y;disp(y);y=m-m*h;endy= 1 (2). >> h=0.1;disp('y=');for t=0:h:1y=exp(-t);disp(y);endy= 1 比擬歐拉方法求解與真值的差異歐拉1真值1誤差0顯然誤差與h

7、2為同階無窮小，歐拉法具有一階計算精度，精度較低，但算法簡單2-5. 用四階龍格庫塔法求解2-3的數(shù)值解，并與前兩題結(jié)果比擬。解：(1)>> h=0.1;disp('y=');y=1;for t=0:h:1disp(y);K1=-y;K2=-(y+K1*h/2);K3=-(y+K2*h/2);K4=-(y+K3*h);y=y+(K1+2*K2+2*K3+K4)*h/6;endy= 1 (2) 比擬這幾種方法：對于四階龍格-庫塔方法真值10.4966 龍庫10.4066誤差00000000000 顯然四階龍格-庫塔法求解精度很高，根本接近真值。三種方法比擬可以得到精度

8、四階 精度二階 精度歐拉。3-2.設(shè)典型閉環(huán)結(jié)構(gòu)控制系統(tǒng)如以下圖所示，當(dāng)階躍輸入幅值時，用sp3_1.m求取輸出y(t)的響應(yīng)。y(t)r(t)_2 解：>> a=0.016 0.864 3.27 3.42 1;b=30 25;X0=0 0 0 0;V=2;n=4;T0=0;Tf=10;h=0.01;R=20;b=b/a(1);a=a/a(1);A=a(2: n+1);A=rot90(rot90(eye(n-1,n);-fliplr(A);B=zeros(1,n-1),1'm1=length(b);C=fliplr(b),zeros(1,n-m1);Ab=A-B*C*V;X

9、=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/2)+B*R;X=X+h*(K1+2*K2+2*K3+K4)/6;y=y,C*X;t=t,t(i)+h;endplot(t,y)3-5. 以下圖中，假設(shè)各環(huán)節(jié)的傳遞函數(shù)為：但；重新列寫聯(lián)接矩陣和非零元素矩陣的響應(yīng)曲線。解：7 70 0;1 0.15 0.21 0;0 1 130 0;1 0.01 0.1 0;1 0.01 0.0044 0;W=0 0 0 0 0 0 0

10、0 0;1 0 0 0 0 0 0 0 -1 ; 0 1 0 0 0 0 0 0 0 ; 0 0 1 0 0 0 0 -1 0;0 0 0 1 0 0 0 0 0 ;0 0 0 0 1 0 0 0 0 ;0 0 0 0 0 1 0 0 0 ;0 0 0 0 0 1 0 0 0 ;0 0 0 0 0 0 1 0 0 ;W0=1;0;0;0;0;0;0;0;0;WIJ=1 0 1;2 1 1;2 9 -1;3 2 1;4 3 1 ;4 8 -1;5 4 1;6 5 1;6 10 -1;7 6 1 ;8 6 1 ;9 7 1;n=9;Y0=1;Yt0=0 0 0 0 0 0 0 0 0 ;h=0.0

11、1;L1=1;T0=0;Tf=10;nout=7;A=diag(P(:,1);B=diag(P(:,2);C=diag(P(:,3);D=diag(P(:,4); 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*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*

12、K3+K4)/6; end y=y,Y(nout); t=t,t(i)+h*L1;endplot(t,y)3-7.的數(shù)據(jù)與曲線，并與四階龍格庫塔法比擬精度。解：P=1 0.01 1 0;0 0.085 1 0.17;1 0.01 1 0;0 0.051 1 0.15;1 0.0067 70 0;1 0.15 0.21 0;0 1 130 0;1 0.01 0.1 0;1 0.01 0.0044 0; W=0 0 0 0 0 0 0 0 0;1 0 0 0 0 0 0 0 -1; 0 1 0 0 0 0 0 0 0; 0 0 1 0 0 0 0 -1 0;0 0 0 1 0 0 0 0 0 ;0

13、 0 0 0 1 0 -0.212 0 0;0 0 0 0 0 1 0 0 0;0 0 0 0 0 1 0 0 0;0 0 0 0 0 0 1 0 0 ;W0=1;0;0;0;0;0;0;0;0;n=9;Y0=1;Yt0=0 0 0 0 0 0 0 0 0;h=0.01;L1=10;T0=0;Tf=20;nout=7;A=diag(P(:,1);B=diag(P(:,2);C=diag(P(:,3);D=diag(P(:,4);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;

14、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); FIC(i)=C(i)/D(i)-A(i)/B(i); FID(i)=D(i)/B(i); else end end end Y=zeros(n,1);X=Y;y=0;Uk=zeros(9,1);Ub=Uk; t=T0:h*L1:Tf;N=length(

15、t); for k=1:N-1 for l=1:L1 Ub=Uk; Uk=W*Y+W0*Y0; Udot=(Uk-Ub)/h; Uf=2*Uk-Ub; X=FI'.*X+FIM'.*Uk+FIJ'.*Udot; Y=FIC'.*X+FID'.*Uf; end y=y,Y(nout); end plot(t,y)hold onQ=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:L

16、1; K1=Ab*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;endt',y'plot(t,y,'*r')3-8.求以下圖非線性系統(tǒng)的輸出響應(yīng)y(t)，并與無非線性環(huán)節(jié)情況進(jìn)行比擬。y(t)e(t)r(t)=105-5 1 0.5 1;0 1 20 0;2 1 1 0;10 1 1 0;WIJ=1 0 1;1 4 -1;2 1 1;3 2 1;

17、4 3 1 ;Z=0 0 0 0;S=0 0 0 0;h=0.01;L1=25;n=4;T0=0Tf=20;nout=4;Y0=10;sp3_8; plot(t,y,'r') hold onZ=4 0 0 0;S=5 0 0 0;sp3_8;plot(t,y,'-')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:mif (WIJ(k,2)=0); W0(WIJ(k,1)=WIJ(k,3);else W(WIJ(k,1),WIJ(k,

18、2)=WIJ(k,3);end;end;for i=1:nif(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);elseendelseFI(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);FIC(i)=C(i)/D(i)-A(i)/B(i);FID(i)=D(i)/B(i);elseendendendY=zeros(n,1);X=Y;y=0;Uk=zeros(n,1);Ubb=Uk;t=T0:h*L1:Tf;