1、计算机仿真技术试题含完整答案docx精品文档一、数值计算，编程完成以下各题（共 20 分，每小题 5 分）1、脉冲宽度为 d ，周期为 T 的矩形脉冲的傅里叶级数如下式描述：f ( )d2sin( n d / T )cos( 2 n )1nd / TTn 1当 n 150， d T1 4，1 / 21 / 2 ，绘制出函数f ( ) 的图形。解：syms n t;f=(sin(n*pi/4)/(n*pi/4)*cos(2*pi*n*t);s=symsum(f,n,1,150);y=(1+2*s)/4;x=-0.5:0.01:0.5;Y=subs(y,t,x);plot(x,Y)2、画出函数 f
2、 ( x) (sin 5x) 2 e0.05x2 5x5 cos1.5x 1.5 x 5.5 x5 在区间 3,5 的图形，求出该函数在区间 3, 5 中的最小值点 x m in 和函数的最小值 f m in .解：程序如下x=3:0.05:5;y=(sin(5*x).2).*exp(0.05*x.2)-5*(x.5).*cos(1.5*x)+1.5*abs(x+5.5)+x.2.5;mix_where=find(y=min(y);xmin=x(mix_where);hold on;plot(x,y);plot(xmin,min(y),go,linewidth,5);str=strcat(,n
3、um2str(xmin),num2str(min(y),);text(xmin,min(y),str);Xlabel(x).精品文档Ylabel(f(x)经过运行后得到的图像截图如下：运行后的最小值点 x m in =4.6 ， f m in = -8337.86253、画出函数 f ( x)cos 2 x e 0 .3 x2.5 x 在 1 ，3 区间的图形， 并用编程求解该非线性方程 f ( x )0 的一个根，设初始点为 x0 2 .解：x=1:0.02:3;x0=2;y=(x)(cos(x).2).*exp(-0.3*x)-2.5*abs(x);fplot(y,1,3);Xlabel(
4、x)Ylabel(f(x)X1=fzero(cos(x).2).*exp(-0.3*x)-2.5*abs(x),x0)运行后求得该方程的一个根为 z=0.3256 。4、已知非线性方程组如下，编程求方程组的解，设初始点为1 0.5 -1.x 2x 72x5 z 23yz 3 0.精品文档解： %在新建中建立函数文件 fun2_4.mfunction f=fun2_4(x)f=x(1).2+x(1)*sqrt(7)+2;x(1)+5*x(3).2-3;x(2).*x(3)+3;%非线性方程组求解主程序 fxxfcz.mx0=1 0.5 -1;fsolve(fun2_4,x0)运行后结果为：ans
5、 =-1.3229 3.2264 -0.9298即是 x=-1.3229 y=3.2264 z=-0.9298 .二、控制系统仿真（ 15 分）.精品文档某控制系统的开环传递函数为：6(1.5s1)(0.12s1)G(S)1)(0.05 s，要求：编制一个完整s(6s1)的程序完成以下各小题的要求，所绘制的图形分别定义为四张图。1）绘制出系统的阶跃信号响应曲线（响应时间为0 30s ）2）绘制出系统的脉冲信号响应曲线（响应时间为0 20 s）3）绘制出系统的斜坡信号响应曲线（响应时间为0 10 s）4）绘制出系统的 Bode 图（要求频率范围为10 2 10 2 rad/sec ）解：由传递函
6、数知， 该传递函数是将其用零极点描述法描述的， 将其化为用传递函数表1.08s29.72s6G(S)6.05 s2s ，所以 num=0 1.08 9.72 6,den=0.3 6.05 1 0 。述的形式为：0.3s3%用传递函数编程求解num=0 1.08 9.72 6;den=0.3 6.05 1 0;sys=tf(num,den);t1=0:0.1:30;figure(1)step(sys) % 绘制出系统的阶跃信号响应曲线t2=0:0.1:20;figure(2)impulse(sys) % 绘制出系统的脉冲信号响应曲线t3=0:0.1:10;figure(3)ramp=t3;lsi
7、m(sys,ramp,t3);% 绘制出系统的斜坡信号响应曲线figure(4)w=10(-2):102;bode(sys,w);% 绘制出系统的 Bode图.精品文档fig(1) 系统的阶跃信号响应曲线fig(2) 系统的脉冲信号响应曲线fig(3) 系统的斜坡信号响应曲线.精品文档fig(4) 系统的 Bode 图三、曲线拟合（ 15 分）已知某型号液力变矩器原始特性参数，要求用多项式拟合的方法编程完成以下各小题：1）用二阶多项式拟合出 K ( i ) 曲线；用三阶多项式拟合出 ( i ) 曲线；用三阶多项式拟合出 B (i ) 曲线。2）用不同的颜色和不同的线型， 将 K ( i ) 的
8、原始特性参数数据点和二阶拟合曲线绘制在同一张图形中； 将 (i ) 的原始特性参数数据点和三阶拟合曲线绘制在同一张图形中；将B ( i ) 的原始特性参数数据点和四阶拟合曲线绘制在同一张图形中。3）运行程序， 写出 K ( i ) 曲线的二阶拟合公式、( i ) 曲线的三阶拟合公式和B (i )曲线的四阶拟合公式。.精品文档解：% 曲线拟合（ Curve fitting ）disp(Input Data-i; Output Data-k(i),eta(i),lambdaB(i):) x=0.065,0.098,0.147,0.187,0.243,0.295,0.344,0.398,0.448,
9、0.499;y1=2.37,2.32,2.23,2.15,2.05,1.96,1.87,1.78,1.69,1.59;y2=0.154,0.227,0.327,0.403,0.497,0.576,0.644,0.707,0.757,0.795;y3=26.775,26.845,27.147,27.549,28.052,28.389,28.645,28.756,28.645,28.243;figure(1)pf1=polyfit(x,y1,2)px1=polyval(pf1,x)plot(x,px1,k)gridxlabel( 转速比 i)ylabel( 变矩比 K)title( 二阶多项式拟合
10、 k曲线 )%pausefigure(2)pf2=polyfit(x,y2,3)px2=polyval(pf2,x).精品文档plot(x,px2,b)gridxlabel( 转速比 i)ylabel( 效率 eta)title( 三阶多项式拟合 eta 曲线 )%pausefigure(3)pf3=polyfit(x,y3,4)px3=polyval(pf3,x)plot(x,px3,-r)gridxlabel( 转速比 i)ylabel( 泵轮转矩系数 lambdaB)title( 四阶多项式拟合 lambdaB 曲线 )%figure(4)pf1=polyfit(x,y1,2)px1=p
11、olyval(pf1,x)plot(x,y1,or,x,px1,k)gridxlabel( 转速比 i)ylabel( 变矩比 K)title( 二阶多项式拟合 k曲线 )Legend( 原始数据 , 拟合曲线 )%将的原始特性参数数据点和二阶拟合曲线绘制在同一张图形中pausefigure(5)pf2=polyfit(x,y2,3)px2=polyval(pf2,x).精品文档plot(x,y2,*m,x,px2,b)gridxlabel( 转速比 i)ylabel( 效率 eta)title( 三阶多项式拟合 eta 曲线 )Legend( 原始数据 , 拟合曲线 ,0)%将的原始特性参数
12、数据点和三阶拟合曲线绘制在同一张图形中pausefigure(6)pf3=polyfit(x,y3,4)px3=polyval(pf3,x)plot(x,y3,pk,x,px3,-r)gridxlabel( 转速比 i)ylabel( 泵轮转矩系数 lambdaB)title( 四阶多项式拟合 lambdaB 曲线 )Legend( 原始数据 , 拟合曲线 ,0)%将的原始特性参数数据点和四阶拟合曲线绘制在同一张图形中y1=poly2str(pf1,x)%K ( i ) 曲线的二阶拟合公式y2=poly2str(pf2,x)%( i ) 曲线的三阶拟合公式y3=poly2str(pf3,x)%
13、B ( i ) 曲线的四阶拟合公式运行后的结果如下：运行后的二阶，三阶，四阶拟合曲线函数为：y1 = 0.01325 x2 - 1.8035 x + 2.491y2 =-0.12713 x3 - 1.6598 x2 + 2.4499 x + 0.0025474y3 =106.7407 x4 - 199.9852 x3 + 95.8404 x2 - 8.7272 x + 26.9754.精品文档四、微分方程求解。（ 25分）自己选择确定一个三阶微分方程，自己设置初始条件，用ode45 方法求微分方程的解。要求 : （例如： d 3 y(t )2 d 2 y(t)4 dy(t )8y(t )1 ，
14、 y(0)0 ， dy (0)1，dt 3dt 2dtdtd 2 y( 0)dt20 ）1）仿真时间 t=30秒2）结果绘制在一张图中， 包括 yt 曲线，一阶 yt 曲线，二阶 yt 曲线，三阶 y t曲线3）用图例命令分别说明四条曲线为“ y t ”，“ y t ”，“ y t ” ，“ y t ”4）定义横坐标为“时间”，纵坐标为“输出”，图形标题名称为“微分方程的解”解：系统方程为d 3 y(t)2 d 2 y(t )4 dy(t )8y(t ) 1, 这是一个单变量三阶常微dt 3dt 2dt分方程。将上式写成一个一阶方程组的形式，这是函数ode45 调用规定的格式。令： y(1)y
15、y (2 ) yy (1)y(3) y y(1) y( 2)y(1)y(2)y(2)y(3)1 8 y(1) 2 y(3) 4y(2)y(3)函数文件程序：function ydot=myfun1(t,y)ydot=y(2);y(3);1-8*y(1)-2*y(3)-4*y(2);主文件程序：.精品文档t=0 30;y0=0;1;0;tt,yy=ode45(myfun1,t,y0);y=(1-yy(:,3)-2*yy(:,2)-4*yy(:,1)/8;plot(tt,y,r,tt,yy(:,1),k,tt,yy(:,2),-g,tt,yy(:,3),-.b);legend(y-t,y -t,y
16、 -t,y -t)title( 微分方程的解 )xlabel(时间 )ylabel(输出 )运行程序后输出图形如下：五、 PID 设计（ 25 分）自己选定一个控制系统，（例如：某单位负反馈系统的开环传递函数为400G( s) ），设计一个 PID 控制器，使系统响应满足较快的上升时间和s( s2 30s 200)过渡过程时间、较小的超调量、静态误差尽可能小。方法要求：用 Ziegler Nichols 方法对三个参数 K p 、 K i 、 K d 进行整定，并比较 PID 控制前后的性能，性能的比较要求编.精品文档程实现（用未加 PID 控制的系统闭环传递函数阶跃响应与加 PID 控制后的
17、闭环传递函数的阶跃响应进行比较）解：1)分析：用 Ziegler Nichols 方法是一种经验方法，关键是首先通过根轨迹图找出 Km和 m，然后利用经验公式求增益，微分，积分时间常数。程序：ng=400;dg=1 30 200 0;rlocus(ng,dg); % 画根轨迹图axis(-30 1 -20 20);gridkm,pole=rlocfind(ng,dg)wm=imag(pole(2)kp=0.6*kmkd=kp*pi/(4*wm)ki=kp*wm/pink=kd kp ki,dk=1 0pausend=conv(nk,ng),dd=conv(dk,dg)n1,d1=feedbac
18、k(ng,dg,1,1)n2,d2=feedback(nd,dd,1,1);% 加PID后的闭环传函figurestep(n1,d1,2)gridhold onpausestep(n2,d2,2)hold off在程序中，首先使用 rlocus 及 rlocfind 命令求出系统穿越增益 Km=12.2961和穿越频率 m=13.0220rad/s ，然后使用 Z N 方程求出参数。.精品文档selected_point =-0.4325 +12.9814ikp =7.3777 kd =0.4450 ki =30.5807为采用 PID 控制前后的系统闭环阶跃响应情况比较。图 6-1 系统的根
19、轨迹图图 6-2 PID 控制前后的系统闭环阶跃响应.精品文档三参数 KP， Ki ， Kd 的整定利用系统的等幅振荡曲线的 Ziegler Nichols 方法控制类控制器的控制参数型KpKiKdP0.5Km0PI0.45Km0.54Km/Tm0PID0.6Km1.2Km/Tm0.072Km/Td2)PID 控制系统的开环传函为：K p K D S K I G (s)S因为式中具有积分项，故如果 G(s) 是 n 型系统，加 PID 控制后系统变为 n+1 型，可由下式根据给定的稳态误差指标确定参数K 。ilim nK I G(s)1s o sess ,G( s)400s( s230 s 2
20、00) 是个 I 型系统，由于系统的开环传递函数中有因为积分项，故为 II型系统，假定单位斜坡输入稳态误差ess0.1，则可以计算出Ki 。即：K 2 sK i G ( s ) s 0 2 K i110 K i 50 . 1已知系统性能指标为：系统相角裕量 PM=80，增益穿越频率 n =4rad/s ，故利用这两个参数来求 Kp，Kd。程序如下：ng=400;dg=1 30 200 0;.精品文档ki=5;wgc=4;pm=80;ngv=polyval(ng,j*wgc);dgv=polyval(dg,j*wgc);g=ngv/dgv;thetar=(pm-180)*pi/180;ejthe
21、ta=cos(thetar)+j*sin(thetar);eqn=(ejtheta/g)+j*(ki/wgc);x=imag(eqn);r=real(eqn);kp=rkd=x/wgcif ki=0dk=1 0;nk=kd kp ki;else dk=1;nk=kd kp;endpausend=conv(nk,ng),dd=conv(dk,dg)n1,d1=feedback(ng,dg,1,1)n2,d2=feedback(nd,dd,1,1) % 加PID控制后的闭环系统传递函数pauseg1m,p1m,wpc1,wgc1=margin(ng,dg)g2m,p2m,wpc2,wgc2=margin(nd,dd)%幅值裕度，相角裕度，相频曲线穿越 -180 时的频率，截止频率w=logspace(-1,2,200);pausefigurebode(ng,dg,w)gridhold on.精品文档bode(nd,dd,w)hold offfigurestep(n1,d1,5)gridhold onpausestep(n2,d2,5)hold off可以得到 :K P 2.0204Kd 0.5281p2m =80.0044wgc2 =4.0004 （即：系统相角裕量 PM=80，增益穿越频率n =4rad/s ）图 6-3 系统 Bode 图.精品文档图 6-4 闭环系统的阶跃响应.
