2024年4月30日发(作者：)

课程名称：数学建模与数学实验

学 院：

专 业：

姓 名：

学 号:

指导老师：

1 / 22

利用

Monte Carlo

方法模拟单服务台排队系统和多服务台排队系统

摘 要

蒙特卡罗方法（Monte Carlo）又称统计模拟法随机抽样技术，是一种随机

模拟方法，以概率和统计理论方法为基础的一种计算方法，是使用随机数（或更

常见的伪随机数）来解决很多计算问题的方法。将所求解的问题同一定的概率模

型相联系，用电子计算机实现统计模拟或抽样，以获得问题的近似解。本文通过

两个具体的服务机构为例，分别说明如何利用蒙特卡洛方法模拟单服务台排队系

统和多服务台排队系统。

单服务台排队系统（排队模型之港口系统）：通过排队论和蒙特卡洛方法解

决了生产系统的效率问题，通过对工具到达时间和服务时间的计算机拟合，将基

本模型确定在

M/M/1

排队模型，通过对此基本模型的分析和改进，在概率论相

关理论的基础之上使用计算机模拟仿真（蒙特卡洛法）对生产系统的整个运行过

程进行模拟，得出最后的结论。

多服务台排队系统（开水供应模型）：为了解决水房打水时的拥挤问题。根

据相关数据和假设推导，最终建立了多服务窗排队M/G/n模型，用极大似然估计

和排队论等方法对其进行了求解，并用Matlab软件对数据进行了处理和绘图。

用灵敏度分析对结果进行了验证。本模型比较完美地解决了水房排队拥挤问题，

而且经过简单的修改，它可以用于很多类似的排队问题。

关键词：蒙特卡洛方法，排队论，拟合优度，泊松流，灵敏度分析。

2 / 22

一、问题重述

港口排队系统：一个带有船只卸货设备的小港口，任何时间仅能为一艘船只

卸货。船只进港是为了卸货，响铃两艘船到达的时间间隔在15分钟到145分钟

变化。一艘船只卸货的时间有所卸货物的类型决定，在15分钟到90分钟之间变

化。

开水供应系统：学院开水房的供水时间有限，水房面积有限，水管易受水垢

堵塞。根据调查数据可知：通畅时几乎无人排队，堵塞时水房十分拥挤。由此可

以看出水房设计存在问题，我们可以把开水房看成是一个随即服务系统，应用排

队论的方法对系统运行状态做定量的描述。

二、基本假设

港口排队系统：通过对问题的重述，那么，每艘船只在港口的平均时间和最

长时间是多少？

若一艘船只的等待时间是从到达到开始卸货的时间，每艘船只的平均等待时

间和最长等待时间是多少？

卸货设备空闲时间的百分比是多少？

船只排队最长的长度是多少？

开水供应系统：

假设Ⅰ、顾客流满足参数为



的Poisson分布，其中



为单位时间到达的顾

客平均数。每个顾客所需的服务时间相互独立，顾客流是无限的，在观测期间平

稳。

假设Ⅱ、排队方式为单一队列的等候制，先到先服务。虽然水房内有多个服

务台，每个服务台都有自己的队列，但同时顾客总是自由转移到最短的队列上，

不可能出现有顾客排队而服务器空闲的情况。本文最后对两种排队方式的比较也

表明这一假设是合理的。

假设Ⅲ、水房共有20个并联的服务台（水龙头），设每个服务台的服务时

间服从某个相同的分布，t和σ分别是服务时间的均值和均方差，γ=σ/ t为偏

离系数。由于锅炉及输水管容量的限制，使t依赖于正在进行服务的水龙头个数

m，设此时平均服务时间t(m)。且存在一临界值 当m<= m0 时，t(m)为常数

3 / 22

t0;m>m0时，管道中的水便分给 m 个龙头流出，从而 t(m)> t0,且 t(m)是 m 的

单增函数。

假设Ⅳ、污垢的积累与时间成线性变化，设为f(x)=kT(k>0,表示污垢积累速

率；T为距上次清理污垢时间间隔。

假设Ⅴ、单位时间为 10 秒。

显然，假设Ⅱ、Ⅲ、Ⅳ都是合理的，对假设 Ⅰ进行拟合优度检验，得出假设Ⅰ

也是合理的。

三、符号约定

开水供应系统用到的符号和参数：

L ——系统内顾客数的期望值；

Lq——系统内排队顾客数的数学期望；

W ——顾客在系统内的平均逗留时间；

Wq——顾客排队等待时间的期望；

P0——系统内有服务台空闲的概率；

ρ=t /n ——系统的服务强度（即用水龙头的程度）；

n ——水龙头的个数。



——Wq的上限值



——Po的上限值

四、问题分析

港口排队系统：

排队论：排队论(Queuing Theory) ，是研究系统随机聚散现象和随机服务

系统工作过程的数学理论和方法，又称随机服务系统理论，为运筹学的一个分支。

本题研究的是生产系统的效率问题，可以将磨损的工具认为顾客，将打磨机当做

服务系统。

M/M/1

：较为经典的一种排队论模式，按照前面的Kendall记号定义，

前面的M代表顾客(工具)到达时间服从泊松分布，后面的M则表示服务时间服从

4 / 22

负指数分布，1为仅有一个打磨机。

排队论研究的基本问题

1.排队系统的统计推断:即判断一个给定的排队系统符合于哪种模型，以便

根据排队理论进行研究。

2.系统性态问题:即研究各种排队系统的概率规律性，主要研究队长分布、

等待时间分布和忙期分布等统计指标,包括了瞬态和稳态两种情形。

3.最优化问题：即包括最优设计(静态优化)，最优运营（动态优化）。

为了得到一些合理的答案，利用计算器或可编程计算器来模拟港口的活动。

假定相邻两艘船到达的时间间隔和每艘船只卸货的时间区间分布，加入两艘船到

达的时间间隔可以是15到145之间的任何数，且这个区间内的任何整数等可能

的出现。再给出模拟这个系统的一般算法之间，考虑有5艘传至的假象情况。

对每艘船只有以下数据：

相邻两艘船到达的时间间隔

卸货时间

船1

20

55

船2

30

45

船3

15

60

船4

120

75

船5

25

80

因为船1在时钟于t=0分钟计时开始后20分钟到达，所以港口卸货设备在

开始时空空闲了20分钟。船1立即开始卸货，卸货用时55分，其间，船2在时

钟开始计时后t=20+30=50分中到达。在船1与t=20+55=75分钟卸货完毕之前，

船2不能开始卸货，这意味着船2在卸货前必须等待75-50=25分钟。

在船2开始卸货之前，船2于t=50+15=65分钟到达，因为船2在t=75分钟

开始卸货，并且卸货需45分钟，所以在船2与t=75+45=120分钟卸货完毕之前，

船3不能开始卸货。这样，船3必须等待120分钟。

船4在t=65+120=185分钟之前没有到达，因此船3已经在t=120+60=180分

钟卸货完毕，港口卸货设备空闲185-180=5分钟，并且，船4到达后立即卸货。

最后，在船4于t=185+75=260分钟卸货完毕之前，船5在t=185+25=210到达，

于是船5在开始卸货前等待260-210=50分钟。

5 / 22

五、模型的建立和求解

港口排队系统：对于问题中存在的服务系统，建立排队论模型，在仅能为一

艘船通过是一个标准的

M/G/1

模型：

所谓

M/G/1

模型，就是输入过程为泊松流时，服务时间为任意的条件之下

的，服务机器只有一个得时候。对于

M/G/1

模型，服务时间T的分布式一般的，

（但是要求期望值

E(T)

和

Var(T)

方差都存在），其他条件和标准的

M/M/1

型相

同。为了达到稳态



1

还是必要的，其中有







E(T)

。

进入队列

船只到达

服务台

…

接受服务

…

船只离去

单服务台单队系统

单服务员的排队模型设：

（1） 船只到来间隔时间服从参数为0.1的指数分布．

（2） 对船只的服务时间服从［４，１５］上的均匀分布．

（3） 排队按先到先服务规则，队长无限制．

系统的假设：

（1） 船只源是无穷的；

（2） 排队的长度没有限制；

（3） 到达系统的船只按先后顺序依次进入服务， 即“先到先服务”。

符号说明

w：总等待时间；c

i

：第i个顾客的到达时刻；b

i

：第i个顾客开始服务

时刻；e

i

：第i个顾客服务结束时刻；x

i

：第i-1个顾客与第i个顾客之间到达

的间隔时间；y

i

：对第i个顾客的服务时间

c

i

=c

i-1

+ x

i

e

i

=b

i

+y

i

b

i

=max(c

i

,e

i-1

)

6 / 22

模拟

框图

初始化：令i=1，e

i-1

=0，w=0

产生间隔时间随机数x

i

~参数为0.1的指数分布

c

i

=x

i

, b

i

=x

i

Y

b

i

＞480

N

确定开始服务时间：b

i

=max(c

i

,e

i-1

)

产生间隔时间随机数x

i

~参数为0.1的指数分布

c

i

=c

i-1

+ x

i

累计等待时间：w=w+b

i

-c

i

准备下一次服务：i=i+1

产生服务时间随机数y

i

~[4,15]的均匀分布

e

i

=b

i

+y

i

i=i-1,t=w/i

输出结果：完成服务个数：m=i

平均等待时间：t

停止

开水供应模型：由假设Ⅱ、Ⅲ可知，若 n时,则 n 个服务台是相互独立，

则相当于服务台之间服从相同分布，即是一个 M/G/n 型排队模型。如果

可以相互帮助的服务系统，平均服务时间 t 为正在服务的服务台数 m 的函数。

考虑一简单情形：当 m时，t(m)=;当

＜

m

≤

n

时，

t

(

m

)=，此

7 / 22

时 m个服务员以的速率进行服务，但总的服务速率总是

即n>

M/G/的排队模型。

时的系统实际相当于

首先得求出临界服务台数

速为v,从而由的含义知：

,设水龙头及输出管直径分别为;水的流

（1-2）

即。由实际估测，=6.5

cm

,

=1.3

cm

.于是>20=n,因此现有的水房

系统服从M/G/20的排队模型。

; （1-3）

; （1-4）

; （1-5）

L=

Wq=L/

（1-6）

； （1-7）

. （1-8）

另外公式中要求<1,否则系统永远不能到稳定状态，排队的人越来越多，即

队长将趋于无穷大。

对水房系统，=2.17,n=20, 当管道通畅时，=7.58,

=0.8224<1代入解出:

8 / 22

=3.45,

＝0.292,

L

=14.97,=0.134,W=0.134,=0.945

可见此时水房内为由15人，而水龙头有20个，面积有10平方米，几乎无人排

队，不会产生拥挤现象。这与我观察的实际情况相符。

根据假设Ⅳ ，水垢的积累与时间成线性递增变化，f(x)=kT。随着水垢的积

累，服务时间相应增加。那么处于水房通畅和爆满这两个极端状态之间的水房运

营情况又如何呢?下面的模型当=12.10时,>1,水房爆满，进一步分析以了解拥

挤情况，拥挤原因以及缓解的办法。

六、模型的检验与评价

港口排队系统：

表1 100艘船港口和系统的模拟结果

一艘船呆在港口97

的平均时间

一艘船呆在港口174

的最长时间

一艘船的平均等23

待时间

一艘船的最长等99

待时间

9 / 22

79 78 81 85 99

121 111 141 140 159

8 5 9 12 24

46 33 64 68 93

卸货设备空闲时0.067

间的百分比

0.079 0.093 0.07 0.069 0.028

上图为一艘船呆在港口的平均时间

上图为一艘船呆在港口的最长时间

一艘船的平均等待时间

10 / 22

上图为一艘船的最长等待时间

上图为一艘船的最长等待时间

以上就是对港口问题的具体分析，其实港口问题还可以从船只的排队角度出

发，我们还可以对多个港口通行做相应的模拟试验，让船主尽量减少等待时间且

港口卸货设备的利用率达到最高，从而是港口的主人获得更大的利润。从排队角

度来解决问题，可以使问题的广度增加，选秘书问题就是一个很典型的例子，可

以从排队角度解决，如果用我在文章中应用的方法来解决也是可以的，

这仅仅是一个港口的小问题，甚至可以说是一个非常简单的问题，但是已经让

我感觉到了数学的美，在老师的引导下慢慢接近一种抽象的美，在写论文的这几

天中，数据的整理和分析是最值得享受的时刻，在Excel里输入自己的数据，是

一种忐忑的感觉，因为在那么多的数据面前，我真的不知道将会发生什么，拟合

的过程就更是有意思了，一次一次的尝试，一次一次的比较，在这个过程中，如

果有一点点的进步都会让我兴奋，数学建模在生活中处处存在，如果真的能够掌

11 / 22

握这个本领，生活一定会变得简单而精彩！

开水供应系统：

一、灵敏度分析：

由公式（1-3）、（1-4）、（1-5）和（1-8）知，直接影响系统各运行指标

的参数是n,



，



，t,其中



为不可控的参数，在分析中可以看成不变。

首先，我们讨论Lq、Po和服务时间t之间的关系。已知服务台数

n=20,



=2.17,



=0.5281均是不变的。用Matlab绘出图1.如下：

图1

由上图可以看出，当服务时间t>8时，随服务时间的增加，系统的排队顾客

迅速增长，而服务台的空闲率Po快速下降。由此可见，该水房在大部分时间不拥

挤，服务台利用率较小，与实际观察相符合。

接着，我们继续讨论服务强度



与Lq、Po之间的关系。已知

t=7.56,



0.5281

.作出图形变化趋势图2，如下：

图2

服务台数目n对Lq、Po的影响

12 / 22



2.17

t7.58

m025

由上图可知，当服务台很少时（n

足够多时，增加服务台数量，对于缩短平均等待队长效果不明显。

二、系统的最优化：

上面我们只讨论了服务时间t、服务台数量n与系统内排队顾客数学期望Lq、

系统内服务台的空闲概率Po之间的关系，但是对于固定的m0,存在Po和Lq之

间的合理分配问题，顾客流大时，Po较小，Lq较大。由于没有给出Lq和Po的

相关数据，不能找到Lq和Po的最优解。现在的另外一种方法是:在两种互相矛

盾的度量（平均等待时间

定上限值和。

现在讨论系统服务台空闲率Po、顾客排队等待时间Wq和服务台数n之间的

关系。假设=0.6,=0.4,服务时间=7.58,=12.10得到图3、图4，如下：

图3

和服务台空闲概率)之间取折中值，即对和规

13 / 22

图4

从上图可知，最优服务台数在=7.58与=12.10两种情况下分别为：n1

=18,=28。

因为顾客到达率在不同时段内是不一样的，所以我们还得继续讨论如何安排

顾客流问题，

如果假定系统中 n 个队列间没有顾客转移，则每个队列平均到达率为 /

n

，

从而成为 n 个 M/G/1 系统，平均服务时间 t 不变，=0.1，仍利用原公式计算

得到多队时系统的运行指标。得到表 4，如下：

表4

运行指标

服务时间 排队方式

t

1

=7.58

L

q

L

W

q

P

0

单队

多对

0.294

1.43×20

32.73

35.35

16.97

48.4

54.37

799.0

0.134

14.3

15.3

353.5

0.945

0.998

0.089

0.307

t

2

=12.10

单队

多对

（每个子系（整个子系

统）

14 / 22

统）

从上表可以看出，单队时等待队长、等待时间都比多队时低，而服务台的利

用率都比多队时高因而具有明显优越性，因此建立一个大水房明显优于建立多个

小水房。

参考文献：

（1）《运筹学》教材编写组编. 运筹学. 北京：清华大学出版社，2008

（2）Jerry Banks，John ，Barry L Nelson 等著. 离散事件系统仿真.

北京：机械工业出版社，2007

（3）《排队论模型与蒙特卡洛仿真》

（4）茆诗松 周纪芗. 概率论与数理统计. 北京：中国统计出版社. 2007

（5）张德丰 等. MATLAB概率与数理统计分析. 北京：机械工业出版社. 2010

（6）姜启源 谢金星. 数学建模案例选集. 北京：高等教育出版社. 2006

（7）杨启帆.数学建模案例集. 北京：高等教育出版社. 2006

附录：

港口排队模型：

编程如下：

clear

cs=100;

for j=1:cs

w(j)=0;

i=1;

x(i)=exprnd(10);

c(i)=x(i);

b(i)=x(i);

while b(i)<=480

15 / 22

y(i)=unifrnd(4,15);

e(i)=b(i)+y(i);

w(j)=w(j)+b(i)-c(i);

i=i+1;

x(i)=exprnd(10);

c(i)=c(i-1)+x(i);

b(i)=max(c(i),e(i-1));

end

i=i-1;

t(j)=w(j)/i;

m(j)=i;

end

pt=0;

pm=0;

for j=1:cs

pt=pt+t(j);

pm=pm+m(j);

end

pt=pt/cs

pm=pm/cs

附录二

排队论中一个感兴趣的问题时，当输入过程是Possion流时，顾客相继到达

的间隔时间T服从什么规律。

定理：设



N



t



,t0



是具有参数



的泊松过程，即

P



N



t



n



是对应的时间间隔序列，则随机变量

T



n0,1,2,

n





t





n!

n

e





t

,n0,1,2,,t0,



T

n

,n1



,t0



是独立同分布的，且服从均

16 / 22





e

值为



1

的负指数分布，即

f



t









t0





0 t0

-



t

。

1

证明 因为

T

1

是Possion过程中第一个顾客到达的时间，所以时间



Tt



等价于



0,t



内没有顾客到达。故

P



Tt



P



N



t



0



1





t





0!

0

e





t

e





t

，进而可得

P



T

1

t



1P



T

1

t



e





t

所以

T

1

是服从均值为



1

的负指数分布。

1、利用Possion过程的独立、平稳增量性质，得

P



T

2

tT

1

s



P



在



t,ts



内没有顾客到达T

1

s



P



在



t,ts



内没有顾客到达





Possion过程的独立性



P



N



ts



N



s



0



P



N



t



N



0



0





Possion过程的平稳增量性质



e





t

P



T

2

t



即

P



T

2

t



1P



T

2

t



1e





t

，故

T

2

也是服从均值为



1

的负指数分布。

2、对于任意的

n1

和

t,s

1

,s

n

0

有

P



T

n

tT

1

s

1

,T

2

s

2

,,T

n1

s

n-1



P



N



ts

1

s

n



N



s

1

s

n-1



0



P



N



t



N



0



0



e





t

即

P



T

n

t



1e





t

，所以对任一

T



n1



，它都服从均值为



1

的负指数分布。证毕。

n

开水供应系统：

MATLAB程序：

1: 顾客到达率 的极大似然估计程序：

x0=zeros(1,66);

x1=ones(1,132);

x2=2.*ones(1,131);

x3=3.*ones(1,110);

x4=4.*ones(1,50);

x5=5.*ones(1,22);

x6=6.*ones(1,10);

x7=7.*ones(1,4);

x8=8.*ones(1,3);

17 / 22

x=[x0,x1,x2,x3,x4,x5,x6,x7,x8];

[Lambdahat,Lambdaci]=poissfit(x,0.05)

结果为 ：Lambdahat =2.1705 Lambdaci =[ 2.0448，2.2961] 即： =2.17

2: 在管道畅通时，服务时间的均值和样本标准差程序：

clear all;

x=[30 35 35 40 40 40 45 45 50 50 55 60 60 60 ...

65 65 65 65 65 65 65 65 65 70 70 70 70 ...

75 75 75 75 75 80 80 80 85 85 85 85 85 95 95 95 95 105 105 ...

125 125 155 245];

x1=mean(x)

x2=std(x)

结果为 ：x1 =75.8000， x2 =34.4840

3: 在管道堵塞时，服务时间的均值和样本标准差程序：

clear all;

x=[30 30 40 45 45 45 55 55 55 65 65 70 70 70 75 75 ...

75 75 80 85 90 95 95 95 95 100 105 105 110 110 ...

125 130 130 130 135 135 140 140 145 145 155 155 ...

160 175 185 185 190 200 205 205 215 240 255 ...

265 300];

x3=mean(x)

x4=std(x)

结果为 ：x3 = 120.9091， x4 = 63.8581

4：Lq、Po和服务时间t之间的关系的程序：

t=0:0.1:10;

c=20; s1=1; s2=0; b=0.5281;x=2.17*t;

p=x./c ; %p为服务强度。x为(



t)

d=1;

for m=1:19

d=m*d;

18 / 22

s1=x.^2/d+s1;

end

c1=20*d;

s2=x.^2/c1/(1-p);

s=s1+s2;

p0=1./s ;

x1=x.^20; %(



t^22)

x2=p0.*x1 ;

n=1/c1./(1-p);

P0=(1-n.*x2); %P0的表达式；

j1=1./(1-p);

j2=p.*j1;

b1=(1+b^2)/2;

Lq=j2.*n.*x2*b1;

plot(x,P0),axis([0,14,0,1])

hold on;

plot(x,Lq,'*'),axis([0,14,0,1.3])

legend('Lq和服务时间t的关系','P0和服务时间t的关系',0)

5：t1=7.58时系统服务台空闲率Po、顾客排队等待时间Wq和服务台数n的关系的

程序：

clear all;

x=16.445;

a=2.17;

c=16:1:25 ;

s1=1;

s2=0;

b=0.5281;

for i=1:10

p(i)=x/c(i);

19 / 22

end

for i=1:10

s1(i)=1;

d(i)=1;

for m=1:c(i)-1

d(i)=m.*d(i);

s1(i)=x.^m/d(i)+s1(i);

end

end

for i=1:10

c1(i)=c(i)*d(i);

f(i)=x^c(i);

s2(i)=f(i)/c1(i)/(1-p(i));

s(i)=s1(i)+s2(i);

p0(i)=1/s(i);

x1(i)=x^c(i);

x2(i)=p0(i)*x1(i);

n(i)=1/c1(i)/(1-p(i));

P0(i)=(1-n(i)*x2(i)) ;

j1(i)=1/(1-p(i)); j2(i)=p(i)*j1(i);

b1=(1+b^2)/2;

Lq(i)=j2(i)*n(i)*x2(i)*b1;

Wq(i)=Lq(i)/a;

end

plotyy(c,P0,c,Wq,'stem','plot')

6：t2=12.10时系统服务台空闲率Po、顾客排队等待时间Wq和服务台数n的关系

的程序：

x=12.10*2.17;

a=2.17;

20 / 22

c=26:1:35

s1=1;

s2=0;

b=0.5128;

for i=1:10

p(i)=x/c(i);

end

for i=1:10

s1(i)=1;

d(i)=1;

for m=1:c(i)-1

d(i)=m.*d(i);

s1(i)=x.^m/d(i)+s1(i)

end

end

for i=1:10

c1(i)=c(i)*d(i);

f(i)=x^c(i);

s2(i)=f(i)/c1(i)/(1-p(i));

s(i)=s1(i)+s2(i);

p0(i)=1/s(i);

x1(i)=x^c(i);

x2(i)=p0(i)*x1(i);

n(i)=1/c1(i)/(1-p(i))

P0(i)=(1-n(i)*x2(i))

j1(i)=1/(1-p(i));

j2(i)=p(i)*j1(i);

b1=(1+b^2)/2;

Lq(i)=j2(i)*n(i)*x2(i)*b1;

21 / 22

Wq(i)=Lq(i)/a;

end

plotyy(c,P0,c,Wq,'stem','plot')

grid on

22 / 22