《计算机仿真技术》试题(含完整答案).docx

精品文档

一、数值计算，编程完成以下各题（共20 分，每小题 5 分）

1、脉冲宽度为 d ，周期为 T 的矩形脉冲的傅里叶级数如下式描述：

f ( )d

2

sin( n d / T )

cos( 2 n ) [1

n d / T

T n 1

当 n 150， d T 1 4， 1 / 2 1 / 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 ( x) (sin 5x) 2 e0.05x25x5 cos1.5x 1.5 x 5.5x5在区间[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('(',num2str(xmin),',',num2str(min(y)),')');

text(xmin,min(y),str);

Ylabel('f(x)')

经过运行后得到的图像截图如下：

运行后的最小值点x

m in =4.6，

f

m in= -8337.8625

3、画出函数 f ( x)cos 2 x e 0 .3 x 2.5 x 在[1，3]区间的图形，并用编程求解该非线性方程 f ( x )0 的一个根，设初始点为x02.

解：

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('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 72

x 5 z 23

yz30

解： %在新建中建立函数文件fun2_4.m

function 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.m

x0=[1 0.5 -1];

fsolve(@fun2_4,x0)

运行后结果为：

ans =-1.3229 3.2264 -0.9298

即是 x=-1.3229 y=3.2264z=-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 2rad/sec）

解：由传递函数知，该传递函数是将其用零极点描述法描述的，将其化为用传递函数表

1.08s29.72s6

G(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;

lsim(sys,ramp,t3);%绘制出系统的斜坡信号响应曲线

figure(4)

w=10^(-2):10^2;

bode(sys,w);% 绘制出系统的Bode图

fig(1)系统的阶跃信号响应曲线fig(2)系统的脉冲信号响应曲线fig(3)系统的斜坡信号响应曲线

fig(4)系统的 Bode 图

三、曲线拟合（15 分）

已知某型号液力变矩器原始特性参数，要求用多项式拟合的方法编程完成以下各小题：

1）用二阶多项式拟合出K ( i ) 曲线；用三阶多项式拟合出( i ) 曲线；用三阶多项式拟合出B (i ) 曲线。

2）用不同的颜色和不同的线型，将K ( i )的原始特性参数数据点和二阶拟合曲线绘制

在同一张图形中；将 (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,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')

grid

xlabel(' 转速比 i')

ylabel(' 变矩比 K')

title('二阶多项式拟合k曲线 ')

%

pause

figure(2)

pf2=polyfit(x,y2,3)

plot(x,px2,'b')

grid

xlabel(' 转速比 i')

ylabel(' 效率 \eta')

title('三阶多项式拟合\eta曲线')

%

pause

figure(3)

pf3=polyfit(x,y3,4)

px3=polyval(pf3,x)

plot(x,px3,'-r')

grid

xlabel(' 转速比 i')

ylabel(' 泵轮转矩系数\lambdaB')

title('四阶多项式拟合\lambdaB 曲线 ' )

%

figure(4)

pf1=polyfit(x,y1,2)

px1=polyval(pf1,x)

plot(x,y1,'or',x,px1,'k')

grid

xlabel(' 转速比 i')

ylabel(' 变矩比 K')

title('二阶多项式拟合k曲线 ')

Legend(' 原始数据 ','拟合曲线')

%将的原始特性参数数据点和二阶拟合曲线绘制在同一张图形中pause

figure(5)

pf2=polyfit(x,y2,3)

plot(x,y2,'*m',x,px2,'b')

grid

xlabel(' 转速比 i')

ylabel(' 效率 \eta')

title('三阶多项式拟合\eta曲线')

Legend(' 原始数据 ','拟合曲线',0)

%将的原始特性参数数据点和三阶拟合曲线绘制在同一张图形中pause

figure(6)

pf3=polyfit(x,y3,4)

px3=polyval(pf3,x)

plot(x,y3,'pk',x,px3,'-r')

grid

xlabel(' 转速比 i')

ylabel(' 泵轮转矩系数\lambdaB')

title('四阶多项式拟合\lambdaB 曲线 ' )

Legend(' 原始数据 ','拟合曲线',0)

%将的原始特性参数数据点和四阶拟合曲线绘制在同一张图形中

y1=poly2str(pf1,'x')%K ( i )

曲线的二阶拟合公式

y2=poly2str(pf2,'x')%( i )

曲线的三阶拟合公式

y3=poly2str(pf3,'x')% B ( i )

曲线的四阶拟合公式

运行后的结果如下：

运行后的二阶，三阶，四阶拟合曲线函数为：

y1 = 0.01325 x^2 - 1.8035 x + 2.491

y2 =-0.12713 x^3 - 1.6598 x^2 + 2.4499 x + 0.0025474

y3 =106.7407 x^4 - 199.9852 x^3 + 95.8404 x^2 - 8.7272 x + 26.9754

四、微分方程求解。（ 25分）

自己选择确定一个三阶微分方程，自己设置初始条件，用ode45 方法求微分方程的解。

要求:（例如： d 3 y(t )2 d 2 y(t)4 dy(t )8y(t ) 1 ， y(0)0， dy (0)1，

dt 3dt 2dt dt

d 2 y( 0)

）

dt20

1）仿真时间 t=30秒

2）结果绘制在一张图中，包括y t 曲线，一阶 y t 曲线，二阶 y t 曲线，三阶 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

y (2 ) y y (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ˊ ˊˊ-t','y-ˊˊˊt')

title(' 微分方程的解 ')

xlabel('时间 ')

ylabel('输出 ')

运行程序后输出图形如下：

五、 PID 设计（25分）

自己选定一个控制系统，（例如：某单位负反馈系统的开环传递函数为

400

G( s)），设计一个PID 控制器，使系统响应满足较快的上升时间和s( s230s 200)

过渡过程时间、较小的超调量、静态误差尽可能小。方法要求：用Ziegler—— Nichols 方法对三个参数K p、 K i、 K d进行整定，并比较PID 控制前后的性能，性能的比较要求编

程实现（用未加PID 控制的系统闭环传递函数阶跃响应与加PID 控制后的闭环传递函数的

阶跃响应进行比较）

解：

1)分析：用 Ziegler —— Nichols 方法是一种经验方法，关键是首先通过根

轨迹图找出 Km和ωm，然后利用经验公式求增益，微分，积分时间常数。

程序：

ng=400;dg=[1 30 200 0];

rlocus(ng,dg); %画根轨迹图

axis([-30 1 -20 20]);grid

[km,pole]=rlocfind(ng,dg)

wm=imag(pole(2))

kp=0.6*km

kd=kp*pi/(4*wm)

ki=kp*wm/pi

nk=[kd kp ki],dk=[1 0]

pause

nd=conv(nk,ng),dd=conv(dk,dg)

[n1,d1]=feedback(ng,dg,1,1)

[n2,d2]=feedback(nd,dd,1,1);%加PID后的闭环传函

figure

step(n1,d1,2)

grid

hold on

pause

step(n2,d2,2)

hold off

在程序中，首先使用 rlocus及rlocfind命令求出系统穿越增益Km=12.2961和穿越频率ωm=13.0220rad/s ，然后使用 Z— N 方程求出参数。

selected_point =-0.4325 +12.9814i

kp =7.3777 kd =0.4450 ki =30.5807

为采用 PID 控制前后的系统闭环阶跃响应情况比较。

图 6-1 系统的根轨迹图

图 6-2 PID 控制前后的系统闭环阶跃响应

三参数 K P ， K i ， K d 的整定

利用系统的等幅振荡曲线的 Ziegler —— Nichols 方法

控制类

控制器的控制参数

型

Kp K i K d

P 0.5Km ∞

0 PI 0.45Km 0.54Km/Tm 0

PID

0.6Km

1.2Km/Tm

0.072Km/T

d

2) PID 控制系统的开环传函为：

K p K D S

K I

G (s)

S

因为式中具有积分项，故如果 G(s) 是 n 型系统，加 PID 控制后系统变

为 n+1 型，可由下式根据给定的稳态误差指标确定参数

K 。

i

lim n K I G(s)

1

s o s

e

ss ,

G( s)

400

s( s

2

30 s 200)

是个 I 型系统，由于系统的开环传递函数中有

因为

积分项，故为 II 型系统，假定单位斜坡输入稳态误差 e ss

0.1

，则可以计算出

Ki 。即：

K 2

sK i G ( s )

s 0

2 K i

1 10 K i 5

0 . 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;

ejtheta=cos(thetar)+j*sin(thetar);

eqn=(ejtheta/g)+j*(ki/wgc);

x=imag(eqn);

r=real(eqn);

kp=r

kd=x/wgc

if ki~=0

dk=[1 0];nk=[kd kp ki];

else dk=1;nk=[kd kp];

end

pause

nd=conv(nk,ng),dd=conv(dk,dg)

[n1,d1]=feedback(ng,dg,1,1)

[n2,d2]=feedback(nd,dd,1,1) %加PID控制后的闭环系统传递函数pause

[g1m,p1m,wpc1,wgc1]=margin(ng,dg)

[g2m,p2m,wpc2,wgc2]=margin(nd,dd)

%幅值裕度，相角裕度，相频曲线穿越-180 °时的频率，截止频率

w=logspace(-1,2,200);

pause

figure

bode(ng,dg,w)

grid

hold on

bode(nd,dd,w)

hold off

figure

step(n1,d1,5)

grid

hold on

pause

step(n2,d2,5)

hold off

可以得到 :

K P 2.0204K d 0.5281

p2m =80.0044wgc2 =4.0004 （即：系统相角裕量 PM=80°，增益穿越频率n =4rad/s ）

图 6-3系统Bode图

图 6-4 闭环系统的阶跃响应