基于MATLAB的数字模拟仿真..

基于MATLAB的数字模拟仿真

摘要：本文阐述了计算机模拟仿真在解决实际问题时的重要性，并较为系统的介绍了使用计算机仿真的原理及方法。对于计算机模拟仿真的三大类方法：蒙特卡罗法、连续系统模拟和离散事件系统模拟，在本文中均给出了与之对应的实例及基于MATLAB模拟仿真的相关程序，并通过实例深入的分析了计算机模拟解决实际问题的优势及不足。

关键词：计算机模拟；仿真原理；数学模型；蒙特卡罗法；连续系统模拟；离散事件系统模拟

在实际问题中，我们通常会面对一些带随机因素的复杂系统，用分析方法建模常常需要作许多简化假设，这样进行处理过后的模型与我们面临的实际问题可能相差很远，以致求解得到答案根本无法应用，这时，计算机模拟几乎成为唯一的选择。本文通过对计算机模拟仿真进行系统地介绍，寻求利用模拟仿真来解决问题的一般方法，并深入探讨了这些方法的长处和不足。我们定义一些具有特定的功能、相互之间以一定的规律联系的对象所组成的总体为一个系统，模拟就是利用物理的、数学的模型以系统为问题解决对象，来类比、模仿现实系统及其演变过程，以寻求过程规律的一种方法。模拟的基本思想是建立一个实验的模型，这个模型包含所研究系统的主要特点，这样做的目的就是通过对这个实验模型的运行，获得所要研究系统的必要信息。另外，系统的运行离不开算法，仿真算法是将系统模型转换成仿真模型的一类算法，在数字仿真模型中起核心和关键作用。

1、所谓计算机仿真

计算机仿真是利用计算机对一个实际系统的结构和行为进行动态演示，以评价或预测该系统的行为效果。它是解决较复杂的实际问题的一条有效途径。针对一个确定的系统，根据运行的相似原理，利用计算机来逼真模仿研究对象（研究对象可以是真实的系统，也可以是设想中的系统），计算机仿真是将研究对象进行数学描述，建模编程，且在计算机中运行实现。

对比于物理模拟通常花费较大、周期较长，且在物理模型上改变系统结构和系数都较困难的诸多缺陷，计算机模拟不怕破坏、易修改、可重用，有更强的系统适应能力。但是计算机模拟也有缺陷，比如受限于系统建模技术，即系统数学模型不易建立、程序调试复杂等。

计算机仿真可以用于研制产品或设计系统的全过程中，包括方案论证、技术指标确定、设计分析、生产制造、试验测试、维护训练、故障处理等各个阶段。

2、计算机仿真的目的

对于一个系统，是否选择进行计算机模拟的问题，基于判断计算机模拟与非计算机模拟方法孰优孰劣的问题。归纳以下运用计算机模拟的情况：

（1）在一个实际系统还没有建立起来之前，要对系统的行为或结果进行分析研究时，计算机仿真是一种行之有效的方法。

（2）在有些真实系统上做实验会影响系统的正常运行，这时进行计算机模拟就是为了避免给实际系统带来不必要的损失。如在生产中任意改变工艺参数可能会导致废品，在经济活动中随意将一个决策付诸行动可能会引起经济混乱。

（3）当人是系统的一部分时，他的行为往往会影响实验的效果，这时运用系统进行仿真研究，就是为了排除人的主观因素的影响。

（4）在实际系统上做实验时，由于系统误差和偶然误差的存在，对实验结果的好坏我们很难作出正确的判断，这时运用计算机模拟，就可以保证每次操作的条件相同，排除误差。（5）有些系统一旦建立起来之后就无法复原，利用计算机模拟可重复性的这一优势，可以获得显著的经济效益。例如要投资建立一个大型企业，要分析它建成之后的经济效益和社会效益，不能用建立起来试试看的办法，因为建成后就无法回到原来的状态了。

3、计算机仿真的分类

计算机模拟分为动态模拟和静态模拟，数值分析中的蒙特卡罗法就是典型的静态模拟，动态模拟可以分为连续系统模拟和离散事件系统模拟。连续系统模拟研究系统的状态随时间连续变化的情况，在解决实际问题时，一般要建立微分方程模型，先确定系统的连续状态变化量，然后将它在时间上进行适当的离散量化处理，并由此模拟系统的运行状态。而离散事件系统模拟讨论的是系统状态只在一些离散时间点上，由于随机时间的推进而发生变化，其问题解决模型一般用流程图或网络来表示。

注：连续系统模拟中虽然有时也考虑一些随机作用，但经常还是把它当作确定性问题去考虑，而在离散事件系统模拟中事件的出现和系统状态变量几乎总都是随机的。但是，值得注意的一点是，在这两种系统模拟中时间具有重要作用，我们都是要观察系统在时间过程中的变化。

4、计算机仿真的原理

事实告诉我们，现实世界充满不确定性，我们所研究的现实对象往往难以摆脱随机因素的影响。要使我们的数学模型能够较真实地刻画实际对象，必须面对这个现实。概率论是用数学的思想和方法处理和研究随机现象的一个有效的工具。但有时它还难以用来处理复杂系统中的随机性。而我们运用计算机来模拟随机现象的方法基于随机数，它经常应用于复杂系统的动态仿真的研究当中。仿真模型是处理复杂系统中随机性的计算机模型,也是使用计算机研究和解决实际问题的一条重要途径。

对随机现象进行模拟，实质上要给出随机变量的模拟，也就是说利用计算机随机地产生一系列数值(称为随机数)，它们的出现要服从一定的概率分布。目前，经常使用的是按照

在实际应用中用哪种随机数生成法，要针对具体的系统做出与之相适应的选择。

（1）当研究对象视为大量相互独立的随机变量之和，且其中每一种变量对总和的影响都很小时，可以认为该对象服从正态分布。

（2）考试分数的偏差、射击命中点与目标的偏差、人的身高、体重等，都可近似看成服从正态分布。

（3）排队服务系统中顾客到达率为常数时的到达间隔、故障率为常数时零件的寿命都服从指数分布。指数分布在排队论、可靠性分析中有广泛应用。

（4）涉及到排队系统、产品检验、天文、物理等领域时可用到泊松分布。

5、计算机仿真的方法

5.1解决计算机模拟的一般步骤：

进行计算机仿真一般要进过四个步骤：系统分析，模型构造，运行与改进和输出结果[]

1。

图中A 表示系统分析，B 表示模型构造，C 表示运行与改进，D 表示输出结果

第一步：系统分析，明确目标。就是要明确问题和提出总体方案。首先要把被仿真系统的内容表达清楚，弄清仿真的目的，然后选择描述这些目标的主要环节和状态变量，明确定义所研究问题的范围、边界和初始条件，并充分估计初始条件对系统主要性能的影响。

第二步：模型构造、采集数据。包括建立模型、收集数据、编写程序、程序验证和模型确认等。建立模型就是选择合适的仿真方法，如时间步长法、事件表法等，确定系统的初始状态， 设计整个系统的仿真流程。最后选择合适的通用语言或仿真语言编写、调试程序。

第三步：模型的运行与改进。首先确定一些具体的运行方案，如初始条件、参数、步长、重复次数等，然后输入数据，运行程序，直到符合实际系统的要求及精度为止。

第四步：模型输出、统计分析。包括提供文件的清单，记录重要的中间结果，输出格式要有利于用户了解整个仿真过程，分析和使用仿真结果。

D 图4-1：仿真步骤流程

5.2针对不同的模拟方法举例

5.2.1基于随机数来进行计算机模拟的蒙特卡罗法（Random simulation ）

蒙特卡罗法属于实验数学的一个分支，它利用随机数进行统计试验，以求得的统计特征值（如均值、概率等）作为待解决问题的数值解。而在这一过程中所作的统计试验称为蒙特卡罗法。蒙特卡罗法的基本思想原理是首先建立一个概率模型，使所求问题的解正好是该模型的参数或其他有关的特征量。然后通过模拟、统计，即多次随机抽样实验，统计出某事件发生的百分比。只要实验次数很大，该百分比便近似于事件发生的概率。而这实际上也就是概率的统计定义。

5.2.1.1蒙特卡罗法的原理

这里引用资料[]2，根据车比雪夫定理，设1X ,2X ,, ,n X 是相独立的随机变量序列，它们服从相同的分布，且有有限的数学期望a 和方差 ，则1X ,2X ,, ,n X 的算术平均值当时∞→n 按概率1收敛于a ，即对于任意0>ε有：

111lim =?

??

???<-∑=∞→εn i i n a x n P 由中心极限定理得到：

?∑∞+∞--=∞→=???

?

??????????≤-dt e a x n P t n i i n 212

211lim πεσ 即当n 很大时，

σ

a x n n

i i -∑=1

1近似服从标准正态分布。

5.2.1.2蒙特卡罗法的实例应用

这里应用蒙特卡洛法解决最古老但又不失经典的例子——浦丰(Buffon)投针问题。

问题重述：若干条相互平行线间的距离为d ，针的长度为l ，将针投出后其与平行线相交的概率为P ,求出其概率。

模型建立：通过计算可知针与平行线相交的概率为：πd l 2=P

其中d 为平行线间距离，l 为针的长度，其中d l <，将针投n 次，若有m 次与平行线相交，则可用n m 作为P 的近似值。由此可计算出π的近似值为：m d nl 2≈π 注：要使结果更加准确就需使n 尽量大一些

编程求解：设x 是一随机变量，它服从区间[]d 0，上的均匀分布。同理，?是服从区间[]π，0上的均匀分布。按照某种抽样法，产生随机变量的可能值，例如进行n 次抽样，得到样本值()i i x ?，，其中n 2,1i ，， =，统计出满足不等式：()d lsin x i i ≥+?的次数

()n

10000=n 时，

模拟得到：概率3185.0=P ，圆周率1397.3=π 这与真实值已经很接近了。

5.2.2基于离散时间点上的离散型随机变量的模拟(Discrete system )

离散系统是指系统状态只在有限的时间点或可数的时间点上有随机事件发生的系统，在这些点上发生的离散事件改变了系统的状态，并假设系统状态的变化是在该时间点上瞬时完成的，即离散事件是瞬时出现的。 5.2.2.1离散型随机变量模拟的原理

设离散随机变量X 的可能值为{}i x ，相应的概率为i i p x X P ==)(，,,2,1 =i 设

()

00=p

，()

∑==n

i i n p p

1

，,,2,1 =n 将)(n p 作为分点，把区间()1,0分为一系列小区间

(

)

()()

n n p p ,1-。设Y 是()1,0上均匀随机变量，则有：

()()()

n n p Y p P ≤<-1=()()1--n n p p =n p =()n x x P =， ,2,1=n

因此可以用随机变量Y 落在内的情况来模拟离散的随机变量X 的取值情况。

5.2.2.2离散型随机变量模拟的实例应用

在离散型随机变量模拟所解决的问题中，排队系统是一类很常见又非常重要的问题，且这类问题广泛存在于银行、医院、超市等现实生活中的领域。由一个或多个服务台和一些等待服务的顾客组成的离散时间系统称为排队系统。这里列举一个简单的排队系统问题，并从中了解离散型随机变量模拟的思想方法。

问题重述：一个理发店有两位服务员A 和B ,顾客随机地到达该理发店，每分钟有一个顾客到达和没有顾客到达的概率均是21，其中%60的顾客理发仅用5分钟,另外%40的顾

客用8分钟。试对前60分钟的情况进行仿真，模拟出完成服务的个数及顾客平均等待时 间t 。

模型建立：本题属于多服务台C M M //随即服务模型。

假设客源是无穷的，且队长不受限制，顾客接受先来后到的服务规则。设顾客总的等待时间为w ，第i 个顾客的到达时刻为i c ，第i 个顾客开始接受服务的时刻为i b ，第i 个顾客服务结束的时刻为i e ，第1-i 个顾客与第i 个顾客之间的时间间隔为i x ，第i 个顾客的服务时间为i y 。这样便容易得到如下关系：

???

??=+=+=---)

,,max(211i i i i

i

i i i i i e e c b y b e x c c

5.2.3基于状态随时间连续变化的连续型随机变量的模拟(Continuity system )

连续系统模拟研究系统的状态随时间连续变化的情况。连续系统虽然研究的是确定性问题，但在实际处理中，是将连续状态变量在时间上进行离散化处理，并由此模拟系统的运行状态，即其状态由于随机时间的推进而发生变化。因此，连续系统的计算机模拟只能是近似的，只要这种近似达到一定的精度，一般就可以满足要求。 5.2.3.1连续型随机变量模拟的原理

处理连续型随机变量模拟问题有多种方法，其中反函数法[]

3是最常用的。其方法核心是通过求概率分布的反函数产生随机数。众所周知随机变量X 的概率分布函数()x F 为定义在

()1,0 区间的单调递增函数，设Y 为区间()1,0的均匀随机变量, 令()Y x F =,只要求出反函

数()Y F

x 1

-=,x 即为具有概率分布函数的随机机数。由概率论的理论可以证明X 和

()Y F 1-相同的概率分布。

若X 的概率密度为()x ?，由()X F Y =,()()?∞

-=

=x

dx x X F Y ?是区间()1,0上均匀分布

的随机变量，如果给定区间()1,0上均匀分布的随机数i r , 则具有给定分布的随机数i x 可由方程中解出：()()?

==i

X i i dx x x F r 0

?

5.2.3.2连续型随机变量模拟的实例应用

下面介绍利用连续型随机变量模拟解决追击性问题。

问题重述：正方形ABCD 的4个顶点各有1人。在某一时刻，4人同时出发以匀速m 1v =，按顺时针方向追逐下一人，如果他们始终保持对准目标，则最终按螺旋状曲线于中心点O 。试求出这种情况下每个人的行进轨迹。

模型建立：建立平面直角坐标系：()11y x ，A ,()22y x ，B ，()33y x ，C ，()44y x ，D ，取时间间隔为t ?，计算每一点在各个时刻的坐标。设某点在t 时刻的坐标为: ()i i y x ，， 则在t t ?+时刻的坐为: ()ααsin ,cos t v y t v x i i ?+?+ 其中 d x x i i -=

+1cos α d

y y i i -=+1sin α ()()2

121i i i i y y x x d -+-=++ 编程求解：设置初始化条件为()0,0A ，()0,10B ，()10,0C ，()10,10D ，取足够小的ε,当ε

述建立的模型，用MATLAB7.0编程（程序见附录三）。

图4-2：追击问题仿真分解图

将整个仿真过程分解(如图4-2)，我们可以直观地看到整个追击过程。为了更直观的看到整个追击过程，本文提供了利用MATLAB7.0编制的仿真动画程序(见附录四)。通过仿真动画,追击的过程可以一目了然地表现出来。

6、计算机仿真的效果分析

在前面已经提到多种在处理实际问题时可以用到计算机仿真的情况，这里将看看仿真效果到底如何。

我们知道随机现象存在于现实的各种事件中，对于一特定事件，如果我们可以求出理论值固然是件好事，但事实上并不是所有事件都可以如我们所愿，然而即便能够求出其理论上的情形，也与现实情况存在不小的出入，这时仿真便显示出其强大的适应能力。计算机仿真可以很好的顾及理论值和现实中的随机情况。

（1）对于可以进行理论计算的情况 这里列举一个简单的例子来加以说明：

观察所对目标的指示真确与否。以投币为例，该模拟实验有两种结果，且每种情况出现的概率都相等。对其100次投币情况进行模拟，得到以下结果：

表6-1：模拟值与理论值对比

从表中我们可以看出，虽然模拟值与理论值不完全一样，但它却更能真实地表现实际情况。

（2）对于无法进行理论计算的情况 当一个系统过于复杂，以至于我们无法进行理论计算时，计算机仿真可以很好地表现出系统的整个动态过程。

以2009年全国大学生数学建模竞赛B 题为例： 医院就医排队是大家都非常熟悉的现象，它以这样或那样的形式出现在我们面前，例如，患者到门诊就诊、到收费处划价、到药房取药、到注射室打针、等待住院等，往往需要排队等待接受某种服务。题目要求解决目前队列越来越长的问题。

对于这种情况，理论计算是相当困难的，通过计算机仿真，这个问题却可以得到很好的解决（仿真结果如图[]4），通过图像可以直观地看到队列的长度随天数的变化趋势，相对于理论计算有显著的优越性。

天数

对应天数的队列长度

队长随天数的变化趋势

图6-2:队长随时间的变化趋势

7、计算机仿真与数学建模

数学建模是利用数学方法解决实际问题的一种实践。即通过抽象、简化、假设、引进变量等处理过程后，将实际问题用数学表达式，通过建立起数学模型，运用数学方法及计算机

技术去解决实际问题。在建立起数学模型后，接下来的工作就是检验模型的正确性、及适应能力，而往往这一步是相当困难的。然而，计算机仿真由于其的易修改、可重用且强大的适应力等优越性，可使这一过程简化并且与实际情况吻合的很好。

图7-1：建模—仿真流程图

另外，计算机仿真在数学建模竞赛中有广泛的应用。在每年举行的数学建模大赛中，据不完全统计：1992 A、施肥效果分析，1997 A、零件的参数设计，1999 A、自动化车床管理， 2001 B、公交车调度，评卷策略，2009 B、眼科病床安排等题目均可用计算机仿真来进行很好的解答。计算机仿真往往是解决数学建模竞赛题目的一把利器。

结束语

计算机仿真作为一种解决问题的方法，有优势，当然也有它的不足之处。但是总的来看，它凭借自身的优势已成为解决数学建模问题最好的方法之一。根据具体的实际问题，建立模型，正确又科学地使用计算机仿真才能将问题解决好。

8、参考文献

[1] 颜薇娜. 基于蒙特卡洛模拟的商业银行排队问题研究[J]. 技术经济与管理研究. 2009年第1期.

[2] 莫兴德. 计算机仿真建模的几种方法.广西大学学报(自然科学版)[J]. 第28卷增刊2003年10月.

[3] 莫兴德. 计算机仿真. PPT.

[4] 辜继明，赵闪，余明明. 眼科病床安排的数学模型. 2009年9月14日. 第14页.

[5] 宋来中，王志明. 数学建模与实验[M]. 科学出版社.2005年8月.

9、附录

附录一：浦丰(Buffon)投针问题程序

n=10000;l=1;a=2;

m=0;

for k=1:n

x=unifrnd(0,a);

e=unifrnd(0,pi);

if x+l*sin(e)>=a

m=m+1;

end

end

p=m/n

pai=2*n*l/(m*a)

附录二：

附录三：追击问题仿真程序

v=1;

dt=0.05;

x=[0 0 10 10];

y=[0 10 10 0];

for i=1:4

plot(x(i),y(i),'.'),hold on

end

d=20;

subplot(2,2,4);

while(d>0.1)

x(5)=x(1);y(5)=y(1);

for i=1:4

d=sqrt((x(i+1)-x(i))^2+(y(i+1)-y(i))^2); x(i)=x(i)+v*dt*(x(i+1)-x(i))/d;

y(i)=y(i)+v*dt*(y(i+1)-y(i))/d;

plot(x(i),y(i),'.'),

hold on

end

end

附录四：追击问题的仿真动画程序

v=1;

dt=0.05;

x=[0 0 10 10];

y=[0 10 10 0];

for i=1:4

plot(x(i),y(i),'.'),hold on

end

d=20;

axis equal,

M=moviein(16);

while(d>0.1)

x(5)=x(1);y(5)=y(1);

for i=1:4

d=sqrt((x(i+1)-x(i))^2+(y(i+1)-y(i))^2); x(i)=x(i)+v*dt*(x(i+1)-x(i))/d;

y(i)=y(i)+v*dt*(y(i+1)-y(i))/d;

plot(x(i),y(i),'.'),

M(:,i)=getframe;

end

hold on

end

movie(M,30)