数学建模-计算机模拟技术

模拟技术

河南师范大学数学学院

1模拟基础

2Monte-Carlo模拟

3模拟模型案例

1模拟基础

1.0模拟的背景、思路

应用领域：

•第二次世界大战期间，J.V.Neumann等人将进行的“中

子扩散”的科研项目取名为“Monte-Carlo”

・运输系统模拟

•摩天大楼安全疏散系统模拟

•国民经济发展模拟

•人口增长系统模拟

•供水系统模拟

•管理系统模拟

•雷达系统模拟

・战争系统模拟

模拟思路：

・“模拟”一对系统抽象建模

・“试验”一根据模型设计算法，编程进行反复试验

・“估计”一根据试验数据

・“收集”一根据试验结果作出判断

1.1模拟的基本知识

1.1.1模拟的概念及作用

•现实系统的数学或逻辑模型可能十分复

杂，例如大多数具有随机因素的复杂系

统，其中的一些随机性因素很难用准确的

数学公式表述，从而也无法对整个系统采

用解析法求解。模拟是处理这类实际问题

的有力工具。

•模拟通常借助于计算机进行。

计算机模拟：在已经建立的数学、逻辑模

型的基础之上，通过计算机试验，对一个

系统按照一定的决策原则或作业规则，由

一个状态变换为另一个状态的行为进行描

述和分析。

模拟的作用：

•对于很难用解析方法加以处理的问题，模

拟是一种有效的技术；

•对建模过程中的假设进行鉴定，对理论研

究的结论加以检验；

•对不同的实现方案进行多次模拟，按照既

定的目标函数对不同方案进行比较，从中

选择最优方案。

1.1.2模拟的分类

通常，模拟时间是模拟的主要自变量。

设计正确的模拟时间推进机理：模拟过程中

应根据系统的特性正确推进模拟时间，使

系统中各要素与发生的事件保持同步。

推进模拟时间的基本方法：

■下次事件法：将模拟时间由一个事件发生

的时间点推进到紧接着的下一次事件发生

的时间点。

•固定时间步长法：模拟时间每次均以相等

的固定步长向前推进，每到达一个新的模

拟时间点需检查相应时间段内是否发生了

事件。需根据实际问题合理设置模拟时间

发生改变的步长。

根据模拟过程中因变量的变化情况进行分类：

1）离散型模拟：因变量在与事件时间有关的具体模

拟时间点呈离散性变化。大多数系统（如排队服务

系统）可采用离散型模拟。

时间推进方法：一般采用下次事件法_

应当重点对系统状态可能发生改变的事件进行描述，

并确迫这些事件之间的逻辑关系。

/K

事件时间时间

图1离散型模拟

排队系统通常采用离散型模拟模型。其中，发生系

统状态变化的事件有两个：一是有顾客到达；二

是服务员完成服务。将最近发生上述两种事件之

一的时刻设置为下次事件发生点。

2）连续型模拟：因变量随时间的改变呈连续性变化。

在大多数计算机模拟过程中，按固定的步长推进

模拟时间。

通常需建立一系列的由系统状态变量组成的

状态方程组，以描述状态变量与模拟时间的关系。

3）混合型模拟：因变量随时间的推移而作连续性的

变化并具有离散性的突变，如库存控制系统。

1.1.3模拟的方式

终态模拟：在规定的时间T内进行模拟运行，时间达

至打时，模拟终止。其性能指标明显取决于系统的

初始状态。

稳态模拟：随着模拟时间的推移，系统的性能逐渐

趋于平稳。其目的是研究非终态系统长期运行条

件下的稳态性能，模拟时间的长短取决于能否获

得系统性能的优良估计（可由模拟输出的精度确

定）。

1.1.4模拟的一般步骤

•明确问题，建立模型。正确描述待研究问题，明

确规定模拟的目的和任务，确定衡量系统性能或

模拟输出结果的目标函数，然后根据系统的结构

及作业规则，分析系统各状态变量之间的关系，

以此为基础建立所研究的系统模型。

•收集和整理数据资料。模拟技术的正确运用，往

往由大量的输入数据作依靠。在随机模拟中，应

认真分析具体收集到的随机性数据资料，确定系

统中随机性因素的概率分布特性，以此为依据产

生模拟过程所必需的抽样数据O

•编制程序，模拟运行。

均值、方差、置信区间等），灵敏性分析，选择

最优方案。

注：模拟结果的统计分析模拟的输出结果是分布特

征未知的随机变量，每次运行的结果仅仅是对该

随机变量所有观察值总体的一次抽样，对总体的

代表性很差，虽然可以增加模拟运行的时间从而

增加抽样次数，但这些数据总是由一个“种子”经

过一定的算法而获得的伪随机序列，它们是自相

关的，并不能构成统计上独立的随机样本。

1.2随机模拟案例：赶上火车的概率

【问题】如图，一列火车从A站开往B站，某人每天

赶往B站上这趟火车.

/某人

火车运行方向

AB

思考：请研究他能否赶上这趟火车。

他已了解到：

1）火车从A站到B站的运行时间是均值为30分

钟，标准差为2分钟的随机变量；

2）火车在下午大约1点离开A站，离开时刻

的频率分布如下：

出发时刻午后1:00午后1:05午后1:10

频率0.70.20.1

他到达B站的时刻的频率分布为

时刻午后1:28午后1:30午后1:32午后1:34

频率0.30.40.20.1

他能否及时赶上火车?

明确问题：他能及时赶上火车的概率是多少?

建模方向（思路）：

i）分析法：用概率统计知识建立分析模型，求解

析解。（思考）

ii）模拟法：用概率统计知识建立模型，通过模拟

求近似解。

即先建立模拟模型，然后通过计算机模拟得到

问题的近似解。在同样条件下多次试验，计算

他能及时赶上火车的频率。

问题分析：能及时赶上火车的充要条件是:

心盟+5

即Ky

其中

T；—火车从A站出发的时刻;

—火车的运行时间；

是什么变

丁3—他到达B站的时刻。量？如何

模拟?

基本假设：

i)假设TjT2,T3都是相互独立的随机变量;

ii)将午后1时记为t=0,设火车运行时间T2

服从正态分布：丫2〜N(30,22)o

建立模型：为了简化计算，将下午1点记为初始时

刻。得到随机变量%和T3的分布律如下：

火车出发时刻明和人到达B站时刻T3的分布

律分别为：

（分）0510

P⑴0.70.20.1

T3（分）28303234

P⑴0.30.40.20.1

能及时赶上火车的概率p=P{T3〈TI+T2}

如果r为在(0,1)均匀分布的随机数，为了模拟随

机变量%和T3,可以通过如下方法：

0,0<r<0.7,

5,0.7<r<0.9,

10,0.9<r<l.则。和可分别

用来模拟随机变

28,0<r<0,3,

量和T3。

30,0.3<r<0.7,

t=<

332,0.7<r<0.9,

34,0.9<r<1.0.

模拟算法设计

输入：赶火车次数（天数）

输出：赶上火车的频率

两种不同风格的算法描述

主要变量说明：

n模拟次数

k临时变量，存储当前累计模拟次数

count存储赶上火车的次数

算法1（分步骤描述）

第1步输入模拟次数n

第2步k=l,count=0

第3步当k〈二n,执行Step4-9,否则执行第H步

第4步生成均匀分布随机数赋给r

第5步由r及公式确定T1模拟火车出发时刻

第6步生成均匀分布随机数赋给r；

第7步由r及公式确定T3模拟人达到时刻

第8步生成正态分布随机数T2模拟火车运行时间

第9步IFT1+T2>T3,count=count+l,END

第10步执行第3步

第11步输出赶上火车频率p二count/n

算法2（伪代码描述）

D初始化：

输入模拟次数n;

count=0;

打）模拟口次

fori=lton,

模拟随机变量L，T2,T3,

分别赋给t3;

iftl+t2>t3,

count二count+1

endif

endfor

app_prb=count/n;

n=input（1输入模拟次数：，）;

count=0;

fori=l:n,

模拟rtl=rand;%模拟随机变量tl（火车从A站出发的时刻）

ifrtl<0.7

程序Tl=0;

elseifrtl>=0.7&rtl<0.9

Tl=5;

elseTl=10•end

T2=30+randn*2;%模均随机菱量t2（火车的运行时间）

%模拟施机变量七3（彳也到达B站的时亥J）

rt3=rand;

ifrt3<0.3

T3=28;

elseifrt3>=0.3&rt3<0.7

T3=30;

elseifrt3>=0.7&rt3<0.9

T3=32;

else

T3=34;

end

ifT3<T1+T2,count=count+l;end

end%for

prob=count/n

模拟结果：每次模拟1000次或5000次

序号模拟次数近似概率p

110000.6280

210000.6920

310000.6530

450000.6490

550000.6260

650000.6288

系统模拟注意事项：

•一次模拟结果毫无意义！

•模拟是试验性的,是思维结果的验证。

•必须进行足够多次的模拟，并对结果进行统计分

析。

系统模拟特点：系统模拟是研究系统，特别是动态

系统的重要方法，对于：

•结构复杂的系统；

•很难用解析方法求出变量关系的系统；

•内部机理不明的“黑箱”系统；

•为验证用其他方法建立的模型及结果,应是较好

的选择。

1.3随机变量的建模

•利用理论分布，基于对问题的实际的、合理的假

设，选择适当的理论分布模拟随机变量

优点是给出了各种理论结果出现的概率，便于进行数学分

析和处理。但此方法仅限于十分简单的情况，当问题越复

杂，数学处理变得越困难，并且丢失了试验数据的信息。

・基于实际数据的频率做近似模拟

优点在于完全与观察数据相符，并且随实际问题的复杂程

度增大不会产生更大的困难，仅增大工作量而已。缺点是

不便于进行数学分析，不得不依赖于模拟得到的统计结果。

•应用常将两种模拟方法结合起来使用

1.4均匀分布随机变量模拟

1.4.1平方取中法

。1.从一个4位数.“开始，称为种子；

。2.将它平方得到一个8位数(必要时前面加0):

93.取中间的4位数作为下一个随机数。

。按照如上方式，就能得到0—9999随机出现的整

数。因而可以用来模拟(0,1)上均匀分布的随机

数。

。缺点：会退化为0或几个数的循环，所以现在一般不

采用这种方法

。思考：如何的断随机数生成算法的好坏？

9functionsimrandl

%平方取中法:会退化为0

x0=6532161%2040,2041不好

%653217出现循环

n=6;。/on位整数

count=0;disp(sprintf('n))

whilex0〜=0%fori=l:100,

t=xOA2;t=mod(t,10A(n+n/2));

t=fix(t/l0A(n/2));disp(sprintf(,%10d,,t))

x0=t;

count=count+1;

end

1.4.2线性同余法

。给定4个正整数即b,c,4

。产生规则：

一+可以取很大,如。

i=(axn+6)modc,cc=231

。特点(问题)：

。1.序列不断重复(循环周期至多，・个数)；

。2.序列各数之间缺少统计独立性

。3.序列是个确定性

•两组在实践中效果较好的参数设置：

531

n九+1=(5rz:n)mod(2-31)

31

1Tl+1=(8%九)mod(2-1)

线性同余法程序：产生o—c整数

9产生n个随机整数

。functionr=randint(x,a,b,c,n)

%输入x为种子

%线性同余法

fori=l:n,

x=mod((a*x+b)fc);

r(i)=x;

end

9调用例子

。x=7;a=lrb=7,c=10;n=20;

r=randint(x,a,b,c,n)

。x=7;a=l,b=7,c=231;n=20;

r=randint(x,a,b,c,n)

。思考题：利用线性同余法实现如下功能

9产生(0,1)上的随机数；

。产生U(a,b)均匀分布的随机数；

。产生0,1，一n的随机整数；

。产生从J到k(jWk)上的随机整数。

。思考解答

(设randint为线性同余法随机数生成函数)

°(1)randint/c;

°(2)a+(b-a)*randint/c;

•(3)INT((n+l)*randint/c)

9(4)j+INT((k-j+l)*randint/c)

1.5Matlab随机模拟函数：

1.5.1常见分布随机变量的模拟

❷1.产生mxn阶［0,1］均匀分布的随机数矩

阵：rand(m,n)

9思考题：请问如何用rand产生均匀分布U(a,b)的

随机数？

•2.产生mxn阶［a,b］均匀分布U(a,b)的随

机数矩阵：unifrnd(a,b,m,n)

•思考题：在Matlab中用rand和unifrnd分别产生10个均

匀分布U(0,5)的随机数。

O3.产生mxn阶均值为〃，标准差为。的正态分布的

随机数矩阵：normrnd(出胆m,n)

0

1

72只G

9randn(m,n)产生标准正态分布的随机数

。当研究对象视为大量相互独立的随机变量之和，且其

中每一种变量对总和的影响都很小时，可以认为该对

象服从正态分布。

。机械加工得到的零件尺寸的偏差、射击命中点与目标

的偏差、各种测量误差、人的身高、体重等，都可近

似看成服从正态分布。

。4.产生mxn阶期望值为〃的指数分布的随机数矩

阵：exprnd("m,n)

❷若连续型随机变量X的概率密度函数为

Xe~xtx>0

/(1)=0x<Q

。其中入＞0为常数，则称X服从参数为入的指数分

布。

。指数分布的期望值为1/入

9排队服务系统中顾客到达率为常数时的到达间

隔、故障率为常数时零件的寿命都服从小指数分

布。

。指数分布在排队论、可靠性分析中有广泛应用。

•注意：Matlab中，产生参数为X的指数分布的命

令为exprnd(l")

顾客到达某商店的间隔时间服从参数为0」的指数分布

指数分布的均值为1/0.1=10。指两个顾客到达商店的

平均间隔时间是10个单位时间.即平均10个单位时间到

达1个顾客.顾客到达的间隔时间可用exprnd(lO)模

拟。

•5.产生mxn阶参数为人的泊松分布的随机数矩

阵：poissrnd(A,m,n)

9设离散型随机变量X的所有可能取值为0,1,2,・..

，且取各个值的概率为

Xke~x

P(X=k)=一卜厂,A：=0,1,2,--,

•其中入＞0为常数，则称X服从参数为人的泊松分

布。

。泊松分布的期望值为入

•泊松分布在排队系统、产品检验、天文、物理等

领域有广泛应用。

O指数分布与泊松分布的关系：

如相继两个事件出现的间隔时间服从参数为1的指数

分布，则在单位时间间隔内事件出现的次数服从参数

为人的泊松分布.即单位时间内该事件出现k次的概率

为：

P(X=A-)=---1—.k=0.1,2,…，

反之亦然。

(1)顾客到达某商店的间隔时间服从参数为0.1的指数

分布今(2)该商店在单位时间内到达的顾客数服从参

数为0」的泊松分布

❷说明：

。(1)指两个顾客到达商店的平均间隔时间是10个单位时

间.即平均10个单位时间到达1个顾客.

。(2)指一个单位时间内平均到达0.1个顾客

敌坦克分队对我方阵地实施突袭,其到达规律服从泊

松分布,平均每分钟到达4辆.

(1)模拟敌坦克在3分钟内到达目标区的数量，以

及在第1,2.3分钟内各到达几辆坦克.

(2)模拟在3分钟内每辆敌坦克的到达时刻。

。求解思路：先明确随机变量及其分布

g(1)每分钟到达车辆数：4用poissrnd(4)进行模

拟。A=4

。(2)两辆车到达的平均间隔时间0.25分钟.

一坦克到达的间隔时间应服从参数为4的负指数分

布，用exprnd(1/4)模拟。

。练习：如何用Matlab编程实现上述问题的模拟？

1.5.2其它随机变量的模拟

X01

概率0.50.5

。模拟方法：

0<r<0.5,

I={I:r〜U(0,l).

。0.5<r<1.

。实现：

。u=rand

ifu<0.5,

X=0;

else

X=l;

X012|3

概率p0.30.20.40.1

。模拟方法:

0,0<r<0.3,

1,0.3<r<0.5,

r〜U(0,l).

2.0.5<r<0.9,

3.0.9<r<1.

•u=rand

ifu<0.3,

X=0;

elseifu<0.5%0.3+0.2

X=l;

elseifu<0.9%0.3+0.2+0.4

X=2;

else

X=3;

end

2Monte-Carlo模拟

2.1Monte-Carlo原理

o蒙特卡罗(MonteCarlo)方法，也称为随机模拟

(randomsimulation)

o基本思想：为了解决数学、物理、工程技术等方

面的问题

。1.首先建立一个概率模型或随机过程，使它的参数等

于问题的解；

。2.然后通过对模型或过程的观察或抽样试验来计算所

求参数的统计特征，最后给出所求解的近似值。

2.2蒙特卡罗法应用

2.2.1求解非线性规划

O对于非线性规划问题：

min/(N),weEn

fgi(X)>0,z=1,2,•••,m.

s，1cij<Xj<bj,j=1,2,•••,n

。用蒙特卡罗法求解的基本思想是：在估计的区域

{（/匕4）|勺€1,2,...,n}

内随机取若干实验点，然后从试验点中找出可行

点，再从可行点中选择最小点.

。基本假设

9试验点的第J个分量叼服从向,%］内的均匀分布.

。符号假设

。P：试验点总数；MAXP：最大试验点总数；

。K:可行点总数；MAXK：最大可行点数；

eX*：迭代产生的最优占・

。Q：迭代产生的最小殖其初始值为计算机所能

表示的最大数.

。算法思路

。先产生一个随机数作为初始试验点，以后则将上一个试

验点的第j个分量随机产生，其它分量不变而产生一新

的试验点.这样，每产生一个新试验点只需一个新的

随机数分量.当K>MAXK或P>MAXP时停止迭代.

框图

例子1

maxz———J72+x^X2+8<z：］+3x2

3叼+/2=10

N1>0

(N2>0

9在Matlab软件包中编程，一般需要设计三个M-文

件：

。主程序；

❷目标函数程序；

。约束条件测试程序.

例子2（飞行管理模型的蒙特卡罗法求解程

序）

在MATLAB软件包中编程，共需3个M文件:randlp.m,mylp.m,

Ipconst.m.主程序为randlp.m.

%mylp.m

functionz=mylp(x)%目标函数

z=2*x(1)A2+x(2)A2-x(1)*x(2)-8*x(1)-3*x(2);%转化为求最小值问

题

%Ipconst.m

functionlpc=lpconst(x)%约束条件

If3*x(1)+x(2)-10<=0.5&3*x(1)+x(2)-10>=-0.5%约束条件的误差

为±0.5

lpc=1;

else

lpc=0;

end

%randlp.m

function[sol,r1,r2]=randlp(a,b,n)%随机模拟解非线性规划

debug=l;

a=0;%试验点下界

b=10;%试验点上界

n=1000;%试验点个数

rl=unifrnd(a,b,n,1);%n/邛介的［a,b］均匀分布随机数矩阵

r2=unifrnd(a,b,n,1);

sol=[rl(1)r2(1)];

z0=inf;

fori=l:n

xl=rl(i);

x2=r2(i);

lpc=lpconst([xlx2])；

iflioc==l

z=mylp([xlx2])；

ifz<z0

zO=z；

sol=[xlx2]；

end

end

end

2.2.2估算圆周率

•(1)采用确定性模型估算圆周率

。问题：请建立确定性模型估算因周率.

❷利用定积分的几何意义：估算曲边梯形面积

•(2)采用蒙特卡罗法估算圆周率

9平行四边形ABCD的面积为22=4

0扇的面积S=7fl2=7T

。现在模拟产生在正方形ABCD中均匀分布的点n个，如

果这n个点中有m个点在该圆内，则圆的面积与正方

形ABCD的面积之比可近似为m/n；

。即菅7詈=万七等

9模拟程序：估算圆周率：

。n=input('请输入产生点的个数：');m=0;

fori=l:n,

x=-l+2*rand;y=-l+2*rand;

ifx"2+y"2<=l,m=m+l;end

end

mypi=4*m/n

。更简洁实现：

9n=10000;xy=-l+2*rand(2,n);

x=xy(l,:);y=xy(2,:);

mypi=4*sum(x.A2+y.'2j=l)/n

。下面是一组运行结果：

模拟次数〃7T的近似

100003.1028

100003.1395

1000003.1418

1000003.1403

9小结

。运用蒙特卡罗法求解问题要先建立概率模型或者随机

过程。

。思考

9请建立更精确的估计圆周率的数学模型。

2.2.3估算定积分

例子3

求定积分二;N2dH=3

。Uo

图3.1:直线〃=2,直线］=1与曲线y=X2

1)频率法

9分析：，/")dN

•条件:

。考虑/(02。的情形

。构造一个矩形包含了曲边梯形，矩形区域高“满足：

d>maxf(x]

a<x<b

。利用大数定理，事件发生频率依概率收敛与事件发生

的概率。

9具体思路：

9在矩形区域内产生”(足够大)个点

。落在曲边梯形内的点为m个，则

Xdx

faf()〜772

P=------------------------------------

Sn

其中矩形区域面积S=(b—Q)d

。则所求定积分

b

f(x)dx七—(6—Q)d

n

。n=l0A6;%产生随机点的个数

a=0;%积分下限

b=l;%积分上限

maxfx=l;%函数最大值(估计)

d=maxfx+l;%确定一个d(作为大矩形的高)

m=0;

fori=l:n,

x=a+rand*(b-a);%随机产生x坐标

y=d*rand;%随机产生y坐标

ify<=xr,％在曲线下方

m=m+l;%计算增力口

end

end

s=m/n*d*(b-a)

O运行程序:

/x2dx70.33321600000000

Jo

9思考题

❷（1）模型中的参数（如矩形区域的高d）对结果有影响

吗？

。（2）这种方法对于函数值可正可负的情形该怎么应用

呢？

2)平均值法

•平均值法是利用随机变量的平均值(数学期望)

来计算定积分1

I=(3.1)

Ja

e考虑一系列在(0,1)上均匀分布的独立随机变

量：则｛/(&)｝，，=1,2,・♦•,n是

一系列相互独立同分布的随机变量，有

6[/(8)]=-j—/f3日=-—

b—a八b—a

。由强大数定律知以概率为1成立

段！£/(金)=£(3.2)

82=1

。因此当八足够大时，可得近似公式

bn

(3.3)

Zi=l

。由此得到计算/的平均值法；

•(1)产生在(0,1)上均匀分布随机数71r2,

•(2)令%=a+(6—a加)=1,2,…，71；

n

。(3)计算Mg/(%)作为/的估计值。

•分析式(3.1)〜皆；，3.3),可发现本质上平均值法是用样

本平均值作总体数学期望的估计2。

。%program:sim_ini_mean

%蒙特卡罗模拟法求定积分(平均值法)

n=input('输入产生随机数个数：')；

ifisempty(n)|(n<10000),

n=10000;

end

s=0;%存储函数值之和

a=0;%积分下限

b=l;%积分上限

fori=l:n,

x=a+(b-a)*rand;%rand(l)

funvalue=x2;%计算函数值

s=s+funvalue;%结合函数计算函数值，并累计到变

量s

end

I=(b-a)*s/n

。运行结果：

sim-int_mean

e输入产生随机数个数：1。0

I=

0.33795894941753

练习1

有30个电子部件D1，。2,••・，。30，它们的使用情况

如下:一开始D1工作;若D1损坏则D,2立即开始工

作;若D.2损坏则D3立即开始工作,其余类推。

没D的寿命X/i=l,2,…,30）服从参数

为人=0.1（小时）的指数分布且相互独立。

令T为30个器件的总使用时间,试用M—C（蒙特卡

罗）模拟法求解P[T>30｝和E（T）a。

3模拟模型案例

例1：某港口提供有足够的泊位供船舶停靠，但现在

仅有一个可供装卸的泊位，船舶先到则先进行装

卸，如果船舶得不到及时装卸而造成的滞期费为

每小时100元。现要弄清该系统的性能，重点考

察船舶进入该港后等待装卸的滞留时间以及卸位

（即装卸用的泊位）的利用率，从而进行经济效

益分析。

首先，对进入该港口的100艘船进行实际调查，记录其活动情况，得到

这100艘船到达港口的时间间隔和装卸时间的分布情况的频数和累积频

率分布。

表1100艘船到达港口的时间间隔频数表

到达间隔56789101112131415161718

（小时）

到达船舶136791011111197654

数（艘）

表2100艘船装卸时间频数表

装卸时间（小时）910111213141516

装卸船舶数（艘）2022191610832

表3船舶到达港口的时间间隔累积分布表

到达间隔56789101112131415161718

累积频率0.010.040.100.170.260.360.470.580.690.780.850.910.961.00

表4船舶装卸时间累积分布表

装卸时间910111213141516

累积频率0.200.420.610.770.870.950.981.00

设：一年内到达港口的船舶数为N,到港船舶滞留

总时间为T1,船舶平均滞留时间为T2,港口支付

滞留费为Y,装卸泊位空闲时间装卸泊位利用率S。

Stepl:模拟船舶到达港口的时间间隔的分布P1,

Step2:模拟船舶装卸所需时间的分布P2,

Setp3:编写程序。

为了比较准确地反映系统的性能，我们采用稳态

模拟方式，可以考察系统运行360天的情况。假

定港口每天24小时连续工作，因此模拟时间T设

置为8640小时。考虑到船舶到达港口的事件是影

响整个系统状态的主要因素，可以将模拟事件设

置为该事件的发生时刻。

所用变量说明：

t:一艘船到达港口的时间；

tl:前一艘船驶离港口的时间；

td：两艘船到达港口的时间间隔；

ts：当前船舶装卸所需时间；

tw：t时刻所有已到达港口的船舶等待装卸时

间总和；

tf：t时刻装卸位已空闲时间总和；

表5某港口对船舶服务的模拟结果

模

拟到达船到港船舶滞留船舶平均滞港口支付滞装卸泊位空装卸泊位

次

数

舶数总时间（小时）留时间（小时）留费阮）闲时间（小时）利用率

1736748010.1674800035396%

272243696.0543690051594%

3737910912.3690190042695%

473657587.8257580038196%

5741750110.1275010034196%

673968409.2668400042295%

经过6次模拟计算，装卸泊位的平均利用率为

95.3%,但到达的船舶的却因得不到及时的服务

而造成平均每艘滞留9.295小时，每年到港船舶总

滞留6843小时，因此港口每年约需支付68万元的

滞留费，这显然是一笔不小的开支。为此，可修

改模拟模型，考虑增设一个装卸泊位，并重新进

行模拟运行，将这些数据提供给决策者一确定是

否需要再新建一个装卸泊位。

练习：某公共汽车站每隔30分钟到达一辆汽车，但

可能有［0,司分钟的误差，此误差大小与前一辆汽

车的运行无关。汽车最多容纳50名旅客，到达该

汽车站时车内旅客人数服从［20,50］的均匀分布。

到站下车的旅客人数服从［3,7］的均匀分布，每

名旅客下车的时间服从口，7］秒的均匀分布。旅

客按每30分钟到达12个人的泊松分布到达汽车

站，单列排队等车，先到先上，如果某位旅客未

能上车，他将不再等候。旅客上车时间服从

[4,12]秒的均匀分布。上下车的规则是：先下后

上，逐个上车，逐个下车。假设每天共发25

辆汽车，现在要求模拟一天汽车的运行情况，了

解一天中在站内等候汽车的总人数、能上车及不

能上车的人数、旅客排队时间分布情况、不能上

车人数的分布情况等。

例2：露天矿用电铲采掘矿石，然后用卡车运送到卸

场。假设共有m辆卡车装运，有n台(nvm)电铲

同时采掘供n辆卡车同时装卸；卸场有s个卸位

(s<=m),可供s辆卡车同时卸料。装运过程以

班为单位，每个班开始工作时，m辆中的n辆卡车

由n台电铲装车，其余m—n辆卡车排队等候，每

辆卡车装满矿石后驶离装料处，排于队列首位的

待装卡车驶至空闲电铲前并调转车头接收装料，

而刚才已装满的卡车驶至卸场时，若卸位空闲则

调头卸完矿石并重新驶会采掘场，排在待装卡车

的队尾。

为了充分发挥每台电铲、每辆卡车的效率，提高

班产量，就需要确定电铲、卡车及卸位之间的一

个适当匹配数量。

假定我们总共要模拟20个班运行的情况。已知每

个班连续工作6小时，并且每个班开始工作时的初

始状态均一致，即全部卡车都在采掘场等待装料。

经过对该矿开采情况进行收集与处理，得到以下有

关数据：

(1)装车时间〜N(1.32,0.27*0.27)；

(2)重载运行时间近似一个常数4分钟；

(3)卡车调头时间近似一个常数0.67分钟；

(4)卸车时间〜e(1/0.74)；

(5)卡车空载运行时间近似一个常数3分钟；

(6)每辆卡车装载量〜N(22.5,0.83*0.83);

记C[i](i=1〜m)为第i辆卡车。每辆卡车在运行过程

中均要经过装车、重载运行、卸车和空载运行四

个阶段(这里将卸车时和返回时接受装料的调头

分别纳入卸车和空载运行)，因此，可为每辆卡

车C[i]定义一个状态C[i].s,代表卡车所处的各个

不同阶段(i=l〜m)。可以假定每辆卡车开始时

均在采掘场，第一辆卡车开始装料后就开始影响

整个系统中各辆卡车的状态，随着时间的推进，

各卡车进入各自的阶段，所以又为每辆卡车C[i]设

置了一个进入某一阶段的时间(i=1〜m)。

这样，就可将模拟时间推进到所有卡车进入各自

阶段时刻的最小值，相应地安排该阶段的服务，

服务结束后进入下一阶段并重新设置状态值。

为简单计，这里仅考虑了在只有一个电铲和

一个卸位(BPn=s=1)的情况下，应该配

备多少辆卡车m才能尽量提高电铲、卡车的

效率以及班产量。

表6不同卡车数的模拟结果表

卡车数每班电铲平每班每辆卡车每班每辆卡车电铲效率每辆卡车平均班产

（辆）均空闲时间装车等待时间卸车等待时间平均效率量（吨）

（分钟）（分钟）（分钟）

5138.77.110.761.47%95.05%3658.07

6101.810.816.271.72%92.50%4265.52

773.215.526.379.66%88.41%4739.19

851.220.838.685.78%83.49%5114.29

9