数学建模论文(蒙特卡罗的多服务台和单服务台排队系统)

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

专业：

姓名：

学号:

指导老师：

利用Monte Carlo方法模拟单服务台排队系统和多服务台排队系统

摘要

蒙特卡罗方法（Monte Carlo）又称统计模拟法随机抽样技术，是一种随机模拟方法，以概率和统计理论方法为基础的一种计算方法，是使用随机数（或更常见的伪随机数）来解决很多计算问题的方法。将所求解的问题同一定的概率模型相联系，用电子计算机实现统计模拟或抽样，以获得问题的近似解。本文通过两个具体的服务机构为例，分别说明如何利用蒙特卡洛方法模拟单服务台排队系统和多服务台排队系统。

单服务台排队系统（排队模型之港口系统）：通过排队论和蒙特卡洛方法解决了生产系统的效率问题，通过对工具到达时间和服务时间的计算机拟合，将基本模型确定在//1

M M排队模型，通过对此基本模型的分析和改进，在概率论相关理论的基础之上使用计算机模拟仿真（蒙特卡洛法）对生产系统的整个运行过程进行模拟，得出最后的结论。

多服务台排队系统（开水供应模型）：为了解决水房打水时的拥挤问题。根据相关数据和假设推导，最终建立了多服务窗排队M/G/n模型，用极大似然估计和排队论等方法对其进行了求解，并用Matlab软件对数据进行了处理和绘图。用灵敏度分析对结果进行了验证。本模型比较完美地解决了水房排队拥挤问题，而且经过简单的修改，它可以用于很多类似的排队问题。

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

一、问题重述

港口排队系统：一个带有船只卸货设备的小港口，任何时间仅能为一艘船只卸货。船只进港是为了卸货，响铃两艘船到达的时间间隔在15分钟到145分钟变化。一艘船只卸货的时间有所卸货物的类型决定，在15分钟到90分钟之间变化。

开水供应系统：学院开水房的供水时间有限，水房面积有限，水管易受水垢堵塞。根据调查数据可知：通畅时几乎无人排队，堵塞时水房十分拥挤。由此可以看出水房设计存在问题，我们可以把开水房看成是一个随即服务系统，应用排队论的方法对系统运行状态做定量的描述。

二、基本假设

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

若一艘船只的等待时间是从到达到开始卸货的时间，每艘船只的平均等待时间和最长等待时间是多少？

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

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

开水供应系统：

假设Ⅰ、顾客流满足参数为λ的Poisson分布，其中λ为单位时间到达的顾客平均数。每个顾客所需的服务时间相互独立，顾客流是无限的，在观测期间平稳。

假设Ⅱ、排队方式为单一队列的等候制，先到先服务。虽然水房内有多个服务台，每个服务台都有自己的队列，但同时顾客总是自由转移到最短的队列上，不可能出现有顾客排队而服务器空闲的情况。本文最后对两种排队方式的比较也表明这一假设是合理的。

假设Ⅲ、水房共有20个并联的服务台（水龙头），设每个服务台的服务时间服从某个相同的分布，t和σ分别是服务时间的均值和均方差，γ=σ/ t为偏离系数。由于锅炉及输水管容量的限制，使t依赖于正在进行服务的水龙头个数m，设此时平均服务时间t(m)。且存在一临界值当m<= m0 时，t(m)为常数

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) ，是研究系统随机聚散现象和随机服务系统工作过程的数学理论和方法，又称随机服务系统理论，为运筹学的一个分支。本题研究的是生产系统的效率问题，可以将磨损的工具认为顾客，将打磨机当做服务系统。

//1

M M：较为经典的一种排队论模式，按照前面的Kendall记号定义，前面的M代表顾客(工具)到达时间服从泊松分布，后面的M则表示服务时间服从

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

排队论研究的基本问题

1.排队系统的统计推断:即判断一个给定的排队系统符合于哪种模型，以便根据排队理论进行研究。

2.系统性态问题:即研究各种排队系统的概率规律性，主要研究队长分布、等待时间分布和忙期分布等统计指标,包括了瞬态和稳态两种情形。

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

为了得到一些合理的答案，利用计算器或可编程计算器来模拟港口的活动。假定相邻两艘船到达的时间间隔和每艘船只卸货的时间区间分布，加入两艘船到达的时间间隔可以是15到145之间的任何数，且这个区间内的任何整数等可能的出现。再给出模拟这个系统的一般算法之间，考虑有5艘传至的假象情况。

对每艘船只有以下数据：

相邻两艘船到达的时间间隔20 30 15 120 25

卸货时间55 45 60 75 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分钟。

五、模型的建立和求解

港口排队系统：对于问题中存在的服务系统，建立排队论模型，在仅能为一艘船通过是一个标准的//1M G 模型：

所谓//1M G 模型，就是输入过程为泊松流时，服务时间为任意的条件之下的，服务机器只有一个得时候。对于//1M G 模型，服务时间T 的分布式一般的，（但是要求期望值()E T 和()Var T 方差都存在），其他条件和标准的//1M M 型相同。为了达到稳态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)

单服务台单队系统

…

…

船只到达

进入队列

服务台

接受服务

船只离去

开水供应模型：由假设Ⅱ、Ⅲ可知，若 n时,则 n 个服务台是相互独立，服从相同分布，即是一个 M/G/n 型排队模型。如果则相当于服务台之间可以相互帮助的服务系统，平均服务时间 t 为正在服务的服务台数 m 的函数。考虑一简单情形：当 m 时，t(m)=;当＜ m ≤n 时，t(m)=，此模

框

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

即n>时的系统实际相当于

M/G/的排队模型。

首先得求出临界服务台数 ,设水龙头及输出管直径分别为;水的流速为v,从而由的含义知：

（1-2）

即。由实际估测，=6.5cm,=1.3cm.于是>20=n,因此现有的水房

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

; （1-3）

; （1-4）

; （1-5）

（1-6）

L=；（1-7）

Wq=L/. （1-8）另外公式中要求ρ<1,否则系统永远不能到稳定状态，排队的人越来越多，即队长将趋于无穷大。

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

ρ=0.8224<1代入解出:

＝0.292, L =14.97,=0.134,W=0.134,=0.945

根据假设Ⅳ，水垢的积累与时间成线性递增变化，f(x)=kT。随着水垢的积累，服务时间相应增加。那么处于水房通畅和爆满这两个极端状态之间的水房运营情况又如何呢?下面的模型当=12.10时, >1,水房爆满，进一步分析以了解拥挤情况，拥挤原因以及缓解的办法。

六、模型的检验与评价

港口排队系统：

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

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

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

一艘船的平均等待时间

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

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

以上就是对港口问题的具体分析，其实港口问题还可以从船只的排队角度出发，我们还可以对多个港口通行做相应的模拟试验，让船主尽量减少等待时间且港口卸货设备的利用率达到最高，从而是港口的主人获得更大的利润。从排队角度来解决问题，可以使问题的广度增加，选秘书问题就是一个很典型的例子，可以从排队角度解决，如果用我在文章中应用的方法来解决也是可以的，

这仅仅是一个港口的小问题，甚至可以说是一个非常简单的问题，但是已经让我感觉到了数学的美，在老师的引导下慢慢接近一种抽象的美，在写论文的这几天中，数据的整理和分析是最值得享受的时刻，在Excel里输入自己的数据，是一种忐忑的感觉，因为在那么多的数据面前，我真的不知道将会发生什么，拟合的过程就更是有意思了，一次一次的尝试，一次一次的比较，在这个过程中，如果有一点点的进步都会让我兴奋，数学建模在生活中处处存在，如果真的能够掌

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

开水供应系统：

一、灵敏度分析：

由公式（1-3）、（1-4）、（1-5）和（1-8）知，直接影响系统各运行指标λ，，t,其中λ为不可控的参数，在分析中可以看成不变。

的参数是n,γ

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

n=20,λ=2.17,γ=0.5281均是不变的。用Matlab绘出图1.如下：

图1

由上图可以看出，当服务时间t>8时，随服务时间的增加，系统的排队顾客迅速增长，而服务台的空闲率Po快速下降。由此可见，该水房在大部分时间不拥挤，服务台利用率较小，与实际观察相符合。

接着，我们继续讨论服务强度ρ与Lq、Po之间的关系。已知γ.作出图形变化趋势图2，如下：

t=7.56,5281

=

.0

图2

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

由上图可知，当服务台很少时（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

图

4

从上图可知，最优服务台数在=7.58

与=12.10两种情况下分别为：n1 =18,

=28。

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

如果假定系统中 n 个队列间没有顾客转移，则每个队列平均到达率为λ / n ，从而成为 n 个 M/G/1 系统，平均服务时间 t 不变，λ=0.1，仍利用原公式计算得到多队时系统的运行指标。得到表 4，如下：

表4

从上表可以看出，单队时等待队长、等待时间都比多队时低，而服务台的利用率都比多队时高因而具有明显优越性，因此建立一个大水房明显优于建立多个小水房。

参考文献：

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

（2）Jerry Banks，John S.Carson，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

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 服从什么规律。

定理：设(){},0N t t ≥是具有参数λ的泊松过程，即(){}(){}

,0,1,2,

,0,,1!

n

t n t P N t n e n t T n n λλ-==

=>≥是对应的时间间隔序列，则随机变量()

0,1,2,

,0n

T n t =>是独立同分布的，且服从均

值为1λ-的负指数分布，即()-t

e

t 0

0 t 0f t λλ?≥?=?

。

证明 因为1T 是Possion 过程中第一个顾客到达的时间，所以时间{}1

T t ≥等价于[)0,t 内没有顾客到达。故{}(){}

()0

1

00!

t t t P T t P N t e e λλλ--≥===

=，进而可得

{}{}111t P T t P T t e λ-<=-≤=

所以1T 是服从均值为1λ-的负指数分布。

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

{}[){}[){}()()(){}

()(){}(){}

2112,, 000 t P T t T s P t t s T s P t t s Possion P N t s N s P N t N Possion e P T t λ-≥==+==+=+-==-===≥在内没有顾客到达在内没有顾客到达过程的独立性过程的平稳增量性质

即{}{}2

211t

P T

t P T t e λ-<=-≥=-，故2T 也是服从均值为1λ-的负指数分布。

2、对于任意的1n ≥和1,,0n t s s ≥有

{}()(){}()(){}11221-111-1,,

,000t n n n n n P T t T s T s T s P N t s s N s s P N t N e λ--≥====++

+-+

+==-==

即 {}t

n

1e P T t λ-<=-，所以对任一()1n

T n ≥，它都服从均值为1λ-的负指数分布。证毕。

开水供应系统： 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);

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;

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);

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;