(毕业论文)《控制系统计算机仿真》上机实验报告

1、兰 州 理 工 大 学控制系统计算机仿真上机实验报告院系： 电信学院 班级： 10级自动化五班 姓名： 学号： 时间： 2013 年 5 月 12 日电气工程与信息工程学院2-2.用matlab语言求下列系统的状态方程、传递函数、零极点增益和部分分式形式的模型参数，并分别写出其相应的数学模型表达式：（1）（2） 解：（1）传递函数转化为状态方程： 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状态方程为： 传递函数转换为零极点增益： 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)零极点增益方程为：传递函数转换为部分分时形式： 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） 解：状态方程转换为传递函数为： 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传递函数为：状态方程转换成零极点： 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）转换成部分分式形式： 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. 用殴拉法求下列系统的输出响应在上，时的数值解。，要求保留4位小数，并将结果与真解比较。解：(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的数值解，并与前两题结果比较。解：(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.设典型闭环结构控制系统如下图所示，当阶跃输入幅值时，用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. 下图中，若各环节的传递函数已知为：但；重新列写联接矩阵和非零元素矩阵，将程序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重求上题输出的数据与曲线，并与四阶龙格库塔法比较精度。解：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.求下图非线性系统的输出响应y(t)，并与无非线性环节情况进行比较。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);