《计算机仿真技术》报告

《计算机仿真技术》实验报告

实验一 数字仿真方法验证

一、实验目的

1．掌握基于数值积分法的系统仿真、了解各仿真参数的影响； 2．掌握基于离散相似法的系统仿真、了解各仿真参数的影响； 3．掌握SIMULINK 动态仿真；

4．熟悉MATLAB 语言及应用环境。

二、实验环境

网络计算机系统，MATLAB 语言环境

三、实验内容、要求

（一）试将示例1的问题改为调用ode45函数求解，并比较结果。

示例1：设方程如下，取步长 h =0.1。上机用如下程序可求出数值解。

调用ode45函数求解： 1）建立一阶微分方程组 du=u-2*t/u

2）建立描述微分方程组的函数m 文件 function du=sy11vdp(t,u) du=u-2*t/u

3）调用解题器指令ode45求解y

[t,u]=ode45('sy11vdp',[0 1],1) plot(t,u,'r-'); xlabel('t'); ylabel('u'); 结果对比：euler 法：

t=1,u=1.7848; RK 法：

t=1,u=1.7321; ode45求解:

t=1,u=1.7321;

[]

1,

01

)0(2∈?????=-=t u u t u dt du

ode45求解t-u 图：

00.10.20.30.4

0.50.60.70.80.91

1

1.11.21.31.41.51.61.71.8t

u

（二）试用四阶RK 法编程求解下列微分方程初值问题。仿真时间2s ，取步长h=0.1。

?????=-=1

)0(2y t y dt dy 四阶RK 法程序：

clear t=2; h=0.1; n=t/h; t0=0; y0=1;

y(1)=y0; t(1)=t0;

for i=0:n-1 k1=y0-t0^2;

k2=(y0+h*k1/2)-(t0+h/2)^2; k3=(y0+h*k2/2)-(t0+h/2)^2 k4=(y0+h*k3)-(t0+h)^2;

y1=y0+h*(k1+2*k2+2*k3+k4)/6; t1=t0+h; y0=y1; t0=t1;

y(i+2)=y1; t(i+2)=t1;

end y t

plot(t,y,'r'); 结果：t=2,y=2.6109

0.5

1

1.5

2

2.5

11.21.41.61.822.22.42.62.83：

（三）试求示例3分别在周期为5s 的方波信号和脉冲信号下的响应，仿真时间20s ，采样周期Ts=0.1。 已知系统的状态空间模型为

u X X X X ??

?

???+????

??????

??--=??????0107814

.07814.05572.02121

[]??

?

???=214493.69691

.1X X Y

程序：

clear

% Create system model

A=[-0.5572 -0.7814 ;0.7814 0]; B=[1;0];

C=[1.9691 6.4493]; D=0;

sys=ss(A,B,C,D);

% impulse response of the system subplot(211)

impulse(sys,[0:0.1:20])

% square response of the system

[u,t] = gensig('square',5,20,0.1); subplot(212) lsim(sys,u,t)

02468101214161820

-2

024

6Impulse Response

Time (sec)

A m p l i t u d e

02468101214161820

5

10

Linear Simulation Results

Time (sec)

A m p l i t u d e

（四）某系统框图如图所示，试用SIMULINK 进行仿真，并比较在无饱和非线性环节下系统仿真结果。

模型：

2Out1

1Out2

20s +12s +20s 32Transfer Fcn3

s+0.5s+0.1Transfer Fcn2

20s +12s +20s 32Transfer Fcn1

s+0.5s+0.1

Transfer Fcn

Step

Scope3

Scope2

Scope1

Scope

1.5Gain1

1.5Gain

0.5

1

2

3

4

56

7

8

9

10

00.10.20.30.40.50.60.70.80.91t

y

y1(saturation)y2

实验二 PID 控制器设计

一、实验目的

1．了解PID 控制原理，掌握相应PID 控制器设计仿真程序的应用； 2．掌握计算机辅助系统瞬态性能指标的计算； 3．掌握计算机辅助系统频率性能分析；

二、实验环境

网络计算机系统（采矿楼四楼测试实验室），MATLAB 语言环境

三、实验内容、要求

已知如图所示单位反馈系统

要求：

1．绘制系统的开环Nyquist 图和Bode 图，并判断该闭环系统是否稳定。 程序：Sy21.m

% Create system model

sys=tf([500 5000],[1 33 337 1775 4950 5000]); % Draw the Nyquist plot figure(1)

nyquist(sys) grid on

[Re,Im,w1]=nyquist(sys); % Draw the Bode plot figure(2) bode(sys)

[mag,phase,w2]=bode(sys); grid on

5000

49501775337335000

5002

3

4

5

++++++s s s s s s

% Create closed-loop feedback system model sysclose=feedback(sys,1);

% Quick access to zero-pole-gain data. [z,p,k]=zpkdata(sysclose,'v'); P

结果：

-1

-0.8-0.6-0.4-0.200.20.40.60.81

-1-0.8-0.6-0.4-0.2

00.20.40.60.8

10 dB

-20 dB

-10 dB

-6 dB

-4 dB

-2 dB

20 dB 10 dB 6 dB 4 dB

2 dB

Nyquist Diagram

Real Axis

I m a g i n a r y A x i

s

-200-150

-100

-50

M a g n i t u d e (d B )10

-2

10

-1

10

10

1

10

2

10

3

-360

-270-180-90

0P h a s e (d e g )

Bode Diagram

Frequency (rad/sec)

p = -19.9388 -5.1469 + 2.5108i -5.1469 - 2.5108i -1.3837 + 3.6577i -1.3837 - 3.6577i

无闭环右极点，系统稳定。

2．应用Ziegler —Nichols 方法设计P 控制器、PI 控制器和PID 控制器。 程序：sy22.m

% Design PID controler % Create System Model clf,clc

sys=tf([500 5000],[1 33 337 1775 4950 5000]); sysgroup=feedback(sys,1); %Design P-PI-PIDControler for i=1:3 type=i;

[sysc,Kp,Ti,Td]=pidmargin(sys,type); sysopen=sysc*sys;

sysclose=feedback(sysopen,1);

sysgroup=append(sysgroup,sysclose);

end

% pidmargin.m

function [sysc,Kp,Ti,Td]=pidmargin(sys,type)

margin(sys)

[Gm,Pm,Wg,Wc]=margin(sys);

Kcr=Gm;

Wcr=Wg;

Tcr=2*pi/Wcr;

switch type

case 1

disp('P Controler')

Kp=0.5*Kcr

Ti='No Design'

Td='No Design'

sysc=Kp;

case 2

disp('PI Controler')

Kp=0.4*Kcr

Ti=0.8*Tcr

Td='No Design'

sysc=Kp*(1+tf(1,[Ti,0]));

case 3

disp('PID Controler')

Kp=0.6*Kcr

Ti=0.5*Tcr

Td=0.12*Tcr

sysc=Kp*(1+tf(1,[Ti,0])+tf([Td,0],1));

end

结果：

Controller\ Kp Ti Td

P 1.7849 / /

PI 1.4279 1.0882 /

PID 2.1419 0.6801 0.1632

3．计算比较原系统与P控制系统、PI控制系统、PID控制系统的瞬态性能指标。

程序：sy23.m

% Design PID controler

% Create System Model

clf

clear

sys=tf([500 5000],[1 33 337 1775 4950 5000]);

sysgroup=feedback(sys,1);

%Design P-PI -PIDControler

for i=1:3

type=i;

[sysc,Kp,Ti,Td]=pidmargin(sys,type);

sysopen=sysc*sys;

sysclose=feedback(sysopen,1);

sysgroup=append(sysgroup,sysclose);

end

clf

for i=1:4

subplot(2,2,i)

step(sysgroup(i,i))

pidsys(i)=sysgroup(i,i); %%

% This program performs an analysis of step response of Design P-PI -PIDControler

[num,den]=tfdata(pidsys(i),'v'); %%

% Compute steady value

Finalvalue(i)=polyval(num,0)/polyval(den,0); %%

% Compute overshoot

[y,t]=step(pidsys(i)); %%

[Ymax,k]=max(y);

PeakTime(i)=t(k); %%

OvershootPercent(i)=100*(Ymax-Finalvalue)/Finalvalue;

% Compute rise time

n=1;

while y(n)<0.1*Finalvalue,n=n+1;end

m=1;

while y(m)<0.9*Finalvalue,m=m+1;end

RiseTime(i)=t(m)-t(n); %%

% Compute settling time

r=length(t);

while (y(r)>0.98*Finalvalue & y(r)<1.02*Finalvalue)

r=r-1;

end

SettlingTime(i)=t(r) ; %%

end

Finalvalue,PeakTime,OvershootPercent,RiseTime,SettlingTim e

结果：

01234

0.20.40.6

0.8sys

Time (sec)A m p l i t u d e

P -sys

Time (sec)A m p l i t u d e

02468

0.5

1

1.5 P I-sys

Time (sec)

A m p l i t u d e

P ID-sys

Time (sec)

A m p l i t u d e

02468

0.5

1

012345

0.5

1

1.5

System\ Finalvalue PeakTime, OvershootPercent RiseTime SettlingTime sys 0.5000 1.1572 22.9897 0.4788 2.9528 P_sys 0.6409 1.0266 59.9731 0.2933 7.9680 PI_sys 1.0000 1.1604 45.2370 0.3094 7.9678 PID_sys

1.0000

0.8954

45.7340

0.2132

4.9886