(畢業(yè)論文)《控制系統(tǒng)計算機仿真》上機實驗報告

1、蘭 州 理 工 大 學控制系統(tǒng)計算機仿真上機實驗報告院系： 電信學院 班級： 10級自動化五班 姓名： 學號： 時間： 2013 年 5 月 12 日電氣工程與信息工程學院2-2.用matlab語言求下列系統(tǒng)的狀態(tài)方程、傳遞函數(shù)、零極點增益和部分分式形式的模型參數(shù)，并分別寫出其相應的數(shù)學模型表達式：（1）（2） 解：（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

2、 1 0b = 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 = -2.7306 + 2.8531

3、i -2.7306 - 2.8531i -1.5388 p = -4.0000 -3.0000 -2.0000 -1.0000k = 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 = 4.0000 -6.0000 2.0000 1.0000p = -4.0000

4、 -3.0000 -2.0000 -1.0000h = g3=residue(r,p,h)g3 =1.0000 7.0000 24.0000 24.0000部分分式形式為：（2） 解：狀態(tài)方程轉(zhuǎn)換為傳遞函數(shù)為： 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,-0.25,-0.75a = 2.2500 -5.0000 -1.2500 -0.5000 2.2500 -4.2500 -1.2500 -0.5000 0.2500 -0.5000 -1.2500 -1.0000 1.2500 -1.750

5、0 -0.2500 -0.7500b=4,2,2,0b = 4 2 2 0c=0,2,0,2c = 0 2 0 2d=0d = 0num,den=ss2tf(a,b,c,d)num = 0 4.0000 14.0000 19.7500 13.0000den = 1.0000 4.0000 5.8125 4.1562 1.5937g1=tf(num,den) transfer function: 4 s3 + 14 s2 + 19.75 s + 13-s4 + 4 s3 + 5.812 s2 + 4.156 s + 1.594傳遞函數(shù)為：狀態(tài)方程轉(zhuǎn)換成零極點： a=2.25,-5,-1.25,-

6、0.5;2.25,-4.25,-1.25,-0.5;0.25,-0.5,-1.25,-1;1.25,-1.75,-0.25,-0.75a = 2.2500 -5.0000 -1.2500 -0.5000 2.2500 -4.2500 -1.2500 -0.5000 0.2500 -0.5000 -1.2500 -1.0000 1.2500 -1.7500 -0.2500 -0.7500b=4,2,2,0b = 4 2 2 0c=0,2,0,2c = 0 2 0 2d=0d = 0z,p,k=ss2zp(a,b,c,d)z = -0.8835 + 1.0463i -0.8835 - 1.0463

7、i -1.7331 p = -0.4318 + 0.6803i -0.4318 - 0.6803i -1.6364 -1.5000 k = 4.0000g2=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 0.0880 - 2.5548i 0.0880 + 2.5548ip = -1.6364 -1.5000 -0.4

8、318 + 0.6803i -0.4318 - 0.6803ih = g3=residue(r,p,h)g3 =4.0000 14.0000 19.7500 13.0000部分分式形式的方程為：2-3. 用毆拉法求下列系統(tǒ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 0.9000 0.8100 0.7290 0.6561 0.5905 0.5314 0.4783 0.4305 0.38740.3487(2). h=0.1;disp(y=);fo

9、r t=0:h:1y=exp(-t);disp(y);endy= 1 0.9048 0.8187 0.7408 0.6703 0.6065 0.5488 0.4966 0.4493 0.40660.3679 比較歐拉方法求解與真值的差別歐拉10.90000.81000.72900.65610.59050.53140.47830.43050.38740.3487真值10.90480.81870.74080.67030.60650.54880.49660.44930.40660.3679誤差0-0.0048-0.00070.01180.01420.01600.01740.0180.0188-0.0

10、192-0.0192顯然誤差與h2為同階無窮小，歐拉法具有一階計算精度，精度較低，但算法簡單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 0.9048 0.8187 0.7408 0.6703 0.6065 0.5488 0.4966 0.4493 0.40660.3679(2) 比較這幾種方法：對于四階龍格-庫塔方法真值1

11、0.90480.81870.7408 0.6703 0.6065 0.54880.4966 0.4493 0.40660.3679龍庫10.90480.81870.7408 0.6703 0.6065 0.54880.4966 0.44930.40660.3679誤差00000000000 顯然四階龍格-庫塔法求解精度很高，基本接近真值。三種方法比較可以得到精度（四階） 精度（二階） 精度（歐拉）。3-2.設典型閉環(huán)結(jié)構(gòu)控制系統(tǒng)如下圖所示，當階躍輸入幅值時，用sp3_1.m求取輸出y(t)的響應。y(t)r(t)_2 解： a=0.016 0.864 3.27 3.42 1;b=30 25;x

12、0=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=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

13、+2*k2+2*k3+k4)/6;y=y,c*x;t=t,t(i)+h;endplot(t,y)3-5. 下圖中，若各環(huán)節(jié)的傳遞函數(shù)已知為：但；重新列寫聯(lián)接矩陣和非零元素矩陣，將程序sp3_2.m完善后，應用sp3_2.m求輸出的響應曲線。解：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

14、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.01;l1=1;t0=0;tf=10;nout=7;a=diag(p(:,1);b=

15、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*k3+k4)/6; end y=y,y(nout); t=t,t(i)+h*l1;endp

16、lot(t,y)3-7.用離散相似法仿真程序sp3_4.m重求上題輸出的數(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 0 0 0 1 0 -0.212 0 0;

17、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;fid(i)=0; if(d(i)=0);

18、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(t); for k=1:n-1 for l=

19、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:l1; k1=ab*y+b1*y0 k2=ab*(y+h*k1/2)+b1*y0 k3=ab*(y+h*

20、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,yplot(t,y,*r)3-8.求下圖非線性系統(tǒng)的輸出響應y(t)，并與無非線性環(huán)節(jié)情況進行比較。y(t)e(t)r(t)=105-5p=0.1 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;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;

21、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,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);