第五讲蒙特卡洛方法排队论讲稿

《第五讲蒙特卡洛方法排队论讲稿》由会员分享，可在线阅读，更多相关《第五讲蒙特卡洛方法排队论讲稿（66页珍藏版）》请在装配图网上搜索。

1、计量数模提高班专题五蒙特卡罗仿真与排队论模型蒙特卡罗仿真与排队论模型柴中林柴中林 2011/4/16计量数模提高班专题五 数学模型的解有两种数学模型的解有两种精确解精确解近似解近似解当然，对解的近似程度以及求解的复杂程度做一定的讨论对建当然，对解的近似程度以及求解的复杂程度做一定的讨论对建模来讲都是有益的。模来讲都是有益的。求解问题，人们总希望得到精确解。但是在很多情况下，精确解求解问题，人们总希望得到精确解。但是在很多情况下，精确解是求不出或很难求出的。在此情况下，求得问题的近似解就是必是求不出或很难求出的。在此情况下，求得问题的近似解就是必然的了。此外，从应用的角度讲，一定程度的近似解就够

2、了。然的了。此外，从应用的角度讲，一定程度的近似解就够了。引言引言计量数模提高班专题五排队论是重要的一类随机模型，而蒙特拉罗方法则是基于随机排队论是重要的一类随机模型，而蒙特拉罗方法则是基于随机理论的一种重要的仿真模拟方法，它们都与不确定现象相关。理论的一种重要的仿真模拟方法，它们都与不确定现象相关。引言引言 自然现象有两类自然现象有两类确定性现象确定性现象不确定性现象不确定性现象随机现象随机现象模糊现象模糊现象计量数模提高班专题五2、能用蒙特卡罗方法编程求解问题；、能用蒙特卡罗方法编程求解问题；1、了解蒙特卡罗方法的原理和适用范围；、了解蒙特卡罗方法的原理和适用范围；3、了解排队问题的特点、

3、基本类型和理论；、了解排队问题的特点、基本类型和理论；4、能对简单的排队问题编程求解。、能对简单的排队问题编程求解。本专题的学习目的本专题的学习目的计量数模提高班专题五一、蒙特卡罗方法一、蒙特卡罗方法简介简介 蒙特卡罗(Monte Carlo，美国一赌城的名称)方法又称统计模拟法、随机抽样技术，是一种以概率和统计理论方法为基础的基于随机模拟的数值计算方法。它将所求解的问题同一定的概率模型相联系，用计算机和随机数实现统计模拟或抽样，再根据统计理论获得问题的近似解。计量数模提高班专题五蒙特卡罗方法的概率论依据：1 设A表示一随机事件，P(A)，fn(A)分别表示A发生的概率和频率，则当n很大时有P

4、(A)fn(A)。2 设X是一随机变量（总体）,E(X)=,D(X)=2分别是X 的期望和方差，X1,X2,Xn,是来自总体X的一个样本，则 分别是 和2的估计。3 该方法也可以估计参数，如（）中含参数，利用的估计式就可估计出参数的值。niiXXn122)(11niiXXn11模拟得到模拟得到计量数模提高班专题五随机数的产生rand产生一个服从U(0,1)分布的随机数rand(m,n)产生mn个服从U(0,1)分布的随机数exprnd()产生一个服从e()分布的随机数poissrnd()产生一个服从P()分布的随机数binornd(n,p)产生一个服从B(n,p)分布的随机数normrnd(,

5、2)产生一个服从N(,2)分布的随机数unifrnd(a,b)产生一个服从U(a,b)分布的随机数其他一些随机变量可利用U(0,1)分布通过适当数学方法得到（参见下面例子）计量数模提高班专题五仿真与模拟的目的和原理 仿真和模拟可以说是针对同一内容的不仿真和模拟可以说是针对同一内容的不同角度的看法描述，当需要对某一问题观同角度的看法描述，当需要对某一问题观察研究而相应的观察和实验时间和成本花察研究而相应的观察和实验时间和成本花费太高时，可以考虑用一个模型代替原型，费太高时，可以考虑用一个模型代替原型，用模型的研究达到原型的研究的目的（以用模型的研究达到原型的研究的目的（以节约时间和成本），这就是

6、仿真，其在计节约时间和成本），这就是仿真，其在计算机上等的实现过程也称为模拟。算机上等的实现过程也称为模拟。计量数模提高班专题五例1：中子穿过原子反应堆屏障问题模拟 原子反应堆外的铅屏障是用来屏障原子反应堆中逸出的中子的，以免给人类造成危害。经试验和分析，可做以下简单假设：每一个进入屏障的中子在撞到铅原子前行进的距离为D，然后这个中子以随机方向弹回来，并且在它的下一次撞击中又行进距离D。假设屏障厚度为3D，每一个中子能经受10次撞击，问进入的中子能以多大的比例穿透铅屏障？二、仿真例子与分析计量数模提高班专题五 分析：该问题显然难以用概率论解决，用蒙特卡罗方法却很容易。为方便，不妨做进一步假设：

7、1 中子反弹回反应堆后认为不能穿过屏障。2 与铅原子相撞后任意方向等可能反弹。3 中子撞击十次后毁灭。画图如右，模拟流程图如下x轴y轴中子屏障计量数模提高班专题五中子撞击铅屏模拟流程图 初始化系统状态 产生一个新中子的初试方向和运行终点 中子回到 反应堆了吗求频率,结束YN 碰撞，产生新方向和运行终点 模拟次数 到了吗NY 中子出了 铅屏了吗 碰撞次数到了吗N 频数增加YYN对复杂的对复杂的对象编程，对象编程，画一个流画一个流程图是很程图是很有价值的有价值的计量数模提高班专题五程序为n=10000;%test numberm=0;%frequencyfor i=1:n theta=unifrn

8、d(-pi,pi);%initinal angle x=cos(theta);%only x needed for j=1:10 theta=2*pi*rand;%new angle hitted x=x+cos(theta);if x3 m=m+1;break;end end endfn=m/n%rate计量数模提高班专题五例2：计算定积分 分析：这个积分应该有精确解，因为原函数的缘故这个积分不易求得，精确解难以得到，故求一个近似解是必然的选择。可以用其他方法求近似解，这里用蒙特卡罗方法。用蒙特卡罗方法离不开随机变量。当问题本身具有随机性时，随机变量的选取与这个随机问题应当一致（如上例）；而

9、当问题本身不具有随机性时（本例），就要引入随机变量，将确定性问题转化为不确定问题，以求得问题的解。dxxI)exp(302计量数模提高班专题五 根据积分区域，引入随机变量 X,且XU(0,3),记其密度含函数为(x)，又记f(x)=exp(-x2),且在0,3上记Y=F(x)=f(x)/(x)，则 模拟结果为0.8704,软件算得结果为0.8862.YYEXFEdxxxFdxxxxfdxxfI)()()()()()()()(303030计算重积分原理相同，且更有价值计量数模提高班专题五例3：用蒙特卡罗法求的（近似）值 求的值已有多种方法，而且要多精确都可以。蒙特卡罗方法求的值效果并不好，这里主

10、要介绍原理。本问题相当于用蒙特卡罗法求一个参数的近似值。计量数模提高班专题五 可以根据圆面积及其积分如上构造模型，也可如下构造：设XU（-1，1），YU（-1，1），且相互独立，则（X，Y）在如右图所示的正方形内服从均匀分布。今考虑事件X2+Y21，该事件就相当于随机变量（X，Y）落在圆周内。利用均匀分布的特征容易得到 ()1()()14nsP AP AfAs阴正，即从而 4fn(A)，程序如下x轴y轴11-1-1计量数模提高班专题五 n=10000;%test number m=0;for i=1:n x=unifrnd(-1,1);%generate the rand number of

11、x y=unifrnd(-1,1);%generate the rand number of y if x2+y20)0)，称，称X X服从参数为服从参数为的泊松分布，的泊松分布，若在上式中引入时间参数若在上式中引入时间参数t t，即令，即令tt代替代替，则有：，则有：在概率论中，我们曾学过泊松分布，设随机变在概率论中，我们曾学过泊松分布，设随机变量为量为X，则有：，则有：!neP Xnn n=0,1,2,（1）与时间有关的随机变量的概率与时间有关的随机变量的概率，是一个，是一个随机过程随机过程，即即泊松过程泊松过程。t0，n=0,1,2,（2）计量数模提高班专题五)()(,1221ntNtN

12、PttPn（t2t1，n0）若设若设N(t)N(t)表示在时间区间表示在时间区间0,t)0,t)内到达的顾客数内到达的顾客数(t0),P(t0),Pn n(t(t1 1,t,t2 2)表示在时间区间表示在时间区间tt1 1,t,t2 2)(t)(t2 2tt1 1)内有内有n(0)n(0)个顾客到达的概率。即：个顾客到达的概率。即：在一定的假设条件下在一定的假设条件下 顾客的到达过程就是顾客的到达过程就是一个泊松过程。一个泊松过程。当当P Pn n(t(t1 1,t,t2 2)符合下述三个条件时，顾客到达过程符合下述三个条件时，顾客到达过程就是泊松过程就是泊松过程(顾客到达形成普阿松流顾客到达

13、形成普阿松流)。计量数模提高班专题五 无后效性：无后效性：各区间的到达相互独立各区间的到达相互独立,即即MarkovMarkov性。性。.t0 t1 t2 tn-1 tn|)(|)(11112211)()(,.,)(,)(nnnnxtxnxtxxtxxtxnntxPntxP 也就是说过程在也就是说过程在t+tt+t所处的状态与所处的状态与t t以前所处的状以前所处的状态无态无关。关。平稳性：平稳性：即对于足够小的即对于足够小的tt，有：，有：)()(tttttP ，1普阿松流具有如下特性：普阿松流具有如下特性：在在t,t+tt,t+t内有一个顾客到达的概率与内有一个顾客到达的概率与t t无关无

14、关,而与而与tt成正比。成正比。计量数模提高班专题五 普通性：普通性：对充分小的对充分小的t,t,在时间区间（在时间区间（t,t+tt）内有内有2 2个或个或2 2个以上顾客到达的概率是一高阶无穷小个以上顾客到达的概率是一高阶无穷小.由此知，在由此知，在(t，t+t)t)区间内没有顾客到达的概率为：区间内没有顾客到达的概率为：)(1),(0tottttP 令令t1 1=0,t=0,t2 2=t,=t,则则P(tP(t1 1,t,t2 2)=P)=Pn n(0,t)=P(0,t)=Pn n(t)(t)0 0 是常数，它是常数，它表示单位时间到达的顾客数，称表示单位时间到达的顾客数，称为概率强度。

15、为概率强度。2)(),(nntotttP 即即 P P0 0+P+P1 1+P+P22=1=1计量数模提高班专题五 其概率密度函数为：其概率密度函数为：tTTedtdF)t(f t0t0 当输入过程是泊松流时，我们研究两顾客相继到当输入过程是泊松流时，我们研究两顾客相继到达的时间间隔的概率分布。达的时间间隔的概率分布。设设T T为时间间隔，分布函数为为时间间隔，分布函数为F FT T（t t），则：），则：F FT T（t t）=PTt=PTt 此概率等价于在此概率等价于在00，t t）区间内至少有）区间内至少有1 1个顾客到个顾客到达的概率。达的概率。tTetPtF 1)(1)(0 t0t0

16、tetP)(0 没有顾客到达的概率没有顾客到达的概率为：为：(由（由（5）式而来）式而来)间隔：间隔：间隔：间隔：间隔间隔 对分布函对分布函数微分数微分计量数模提高班专题五 表示单位时间内顾客平均到达数。表示单位时间内顾客平均到达数。1/表示顾客到达的平均间隔时间。表示顾客到达的平均间隔时间。可以证明，间隔时间可以证明，间隔时间T T独立且服从负指数分布与独立且服从负指数分布与顾客到达形成泊松流是等价的。顾客到达形成泊松流是等价的。对顾客的服务时间对顾客的服务时间：系统处于忙期时系统处于忙期时两顾客相继离两顾客相继离开系统的时间间隔开系统的时间间隔，一般地也服从负指数分布，设：，一般地也服从负

17、指数分布，设：即即T服从负指数分布，它的期望及方差为：服从负指数分布，它的期望及方差为：1TE21 TVar 接受服务，然后离开接受服务，然后离开服务时间的分布：服务时间的分布：计量数模提高班专题五其中：其中：表示单位时间内能被服务的顾客数，即平均表示单位时间内能被服务的顾客数，即平均 服务率。服务率。1/1/表示一个顾客的平均服务时间。表示一个顾客的平均服务时间。tetF 1)(tetf )(，则，则 令令 ，则，则称为服务强度称为服务强度。一般的，要设一般的，要设1（否则队列会无限，永远达不到稳态）计量数模提高班专题五 对排队模型，在给定输入和服务条件下，主要对排队模型，在给定输入和服务条

18、件下，主要研究系统的下述运行指标：研究系统的下述运行指标：(1)(1)系统的系统的平均队长平均队长LsLs(期望值期望值)和和平均队列长平均队列长LqLq(期望值期望值)；(2)(2)系统中系统中顾客平均逗留时间顾客平均逗留时间WsWs与队列中与队列中平均平均等待时间等待时间WqWq；本节只研究本节只研究M/M/1M/M/1模型模型.4.5 M/M/1 4.5 M/M/1 模型研究模型研究计量数模提高班专题五（M/M/1/FCFS）系统中有系统中有n n个顾客个顾客 在任意时刻在任意时刻t t，状态为，状态为n n的概率的概率P Pn n(t(t)（瞬态概率），（瞬态概率），它决定了系统的运行

19、特征。它决定了系统的运行特征。已知顾客到达服从参数为已知顾客到达服从参数为的泊松过程，服务时的泊松过程，服务时间服从参数为间服从参数为的负指数分布。现仍然通过研究区间的负指数分布。现仍然通过研究区间 t，t+tt）的变化来求解。在时刻）的变化来求解。在时刻t+tt，系统中有，系统中有n n个顾客不外乎有下列四种情况（个顾客不外乎有下列四种情况（t，t+tt）内到达）内到达或离开或离开2 2个以上个以上没列入）。没列入）。？计量数模提高班专题五区区间间(t t,t+t t)情情况况 时时刻刻t t的的顾顾客客 到到达达 离离去去 时时刻刻t+t t的的顾顾客客 (t t,t+t t)的的概概率率

20、 0 0,t+t t 的的概概率率 A n n 1-t+O(t)1-t+O(t)Pn(1-t+O(t)（1-t+O(t)）B n+1 n 1-t+O(t)t+O(t)Pn+1(1-t+O(t)（t+O(t)）C n-1 n t+O(t)1-t+O(t)Pn-1(t+O(t)（1-t+O(t)）D n n t+O(t)t+O(t)Pn(t+O(t)（t+O(t)）由于这四种情况是互不相容的，所以由于这四种情况是互不相容的，所以Pn(t+t)t)应应是这四项之和，则有：是这四项之和，则有：tttPtttPtttPttPnnnn)1)()()1)(1)()(1)()1()(1tOtttPn 所有的高

21、阶无所有的高阶无穷小和并穷小和并计量数模提高班专题五t)t(P)t(P)t(O)t(Ott)(t(Pnnn 11)t(O)t(O)t(Pt)t(P)t(O)t(Pnnn 111)t(Ot)t(Pt)t(P)tt)(t(Pnnn 111t)t(O)t(P)()t(P)t(Pt)t(P)tt(Pnnnnn 11 令令t0t0，得关于，得关于P Pn n(t)(t)的微分差分方程：的微分差分方程：)()()()()(11tPtPtPdttdPnnnn(1)当当n=0时，只有表中的（时，只有表中的（A）、（）、（B）两种情况，）两种情况，因为在较小的因为在较小的tt内不可能发生（内不可能发生（D D）

22、（到达后即离）（到达后即离去），若发生可将去），若发生可将tt取小即可。取小即可。计量数模提高班专题五)t()t)(t(P)t)(t(P)tt(P 11100)t(P)t(Pdt)t(dP100 (2)(2)这种系统状态这种系统状态(n)随时间变化的过程就是生灭过随时间变化的过程就是生灭过程（程（Birth and Death Process）。）。方程方程(1)、(2)解是瞬态解解是瞬态解,无法应用。无法应用。它 对 时它 对 时间的导数为间的导数为0，所以由，所以由(1)、(2)两式得：两式得：稳态时，稳态时，Pn(t)与时间无关与时间无关,可以写成可以写成Pn,011 nnnP)(PP0

23、10 PP(3)(3)(4)(4)计量数模提高班专题五 由此可得该排队系统的状态转移图：由此可得该排队系统的状态转移图：由（由（4）得：）得：001PPP 其中其中服务强度服务强度 将其代入（将其代入（3）式并令）式并令n=1,2,(也可从状态转移也可从状态转移图中看出状态平衡方程图中看出状态平衡方程)得：得：011 nnnP)(PP010 PP(3)(3)(4)(4)关于关于Pn的差的差分方程分方程 n-1n-1 n n n+1n+1 2 2 0 0 1 1 状态转移图状态转移图计量数模提高班专题五0120 P)(PP n=1n=10020 P)(PP 0202021PP)(P)(P n=2

24、n=20231 P)(PP00230 P)(PP 0303022231PP)(P)(P 计量数模提高班专题五 以此类推以此类推，当，当n=n时，时，00)(PPPnnn (5)(5)1 10 nnP 及概率性质知：及概率性质知：111000 PPnn(数列的极限为数列的极限为 )11 10Pnn)(P 1(6)(6)否则排队无限远否则排队无限远 系统系统稳态稳态概率概率系统的运行指标系统的运行指标计量数模提高班专题五 (1)系统中的队长系统中的队长Ls（平均队长）（平均队长）nnnnsnPnL 001.)(n.)()()(n 11312132.nn.nn 1433223322 132.n(01

25、)1 Ls 即即:(7)(7)期望期望计量数模提高班专题五(2)队列中等待的平均顾客数队列中等待的平均顾客数Lq nnnnqnPnL 111)1()1(nnnn)(n 1111 12sL(8)(8)(3)顾客在系统中的平均逗留时间顾客在系统中的平均逗留时间Ws 顾客在系统中的逗留时间是随机变量，可以证顾客在系统中的逗留时间是随机变量，可以证明，它服从参数为明，它服从参数为-的负指数分布，分布函数的负指数分布，分布函数 11nn计量数模提高班专题五 和密度函数为：和密度函数为：w)(e)w(F 1w)(e)()w(f （w00）)WsLs(1wEWs (4)(4)顾客在队列中的平均等待时间顾客在

26、队列中的平均等待时间W Wq q 111WsWq)WqLq(等待时等待时间间 顾客在队列中的平均等待时间应为顾客在队列中的平均等待时间应为W Ws s减去平均服减去平均服务时间。务时间。考虑考虑L LS S与与W WS S的关的关系系计量数模提高班专题五 LsLsLLqs 12WsLs WqLq 四个指标的关系为四个指标的关系为(Little Little 公式公式)：计量数模提高班专题五知识点评：从上述分析可以看出：排队论是有用的，同时又非常难。根据顾客的到来，排队方式和服务方式可以存在很多排队模型。最简单的M/M/1模型就这么复杂，何况其他模型。事实上能进行理论概率研究的排队模拟是不多的，

27、而这些模型多少还有些理想化的成分。对于生活中许多一般化的排队问题，是求不出其概率分布，得不到其平均队长、等待时间等标志性参数的，怎么办？这就需要仿真和模拟，而蒙特卡罗方法是求解这类问题的有效方法计量数模提高班专题五 此外，对于给定的排队问题，若给出了顾客到达和服务时间的分布类型，可直接模拟。若给出的是一些顾客到达和服务时间的数据，则要根据数据特点判断可能的分布类型，再进行分布类型的检验，以得到顾客到达和服务时间的分布类型及相关参数，再进行模拟。计量数模提高班专题五例4：M/M/1模型及其模拟流程图 M/M/1排队模型即是顾客到来的时间间隔和服务时间都服从负指数分布，只有一个服务台（员），队长无

28、限（来了就排队）。我们先分析排队和服务的特点和规则，再给出模拟流程图，而后用适当参数计算顾客的平均等待时间和平均逗留时间。计量数模提高班专题五排队过程分析 因为 Matlab需要，队伍中至少要有两个人（可虚拟，将来会来）。假设有一个顾客刚开始服务，在它服务期间，顾客仍会来。当他服务结束后，若没有等待顾客，时间推到下一个顾客到来，等待时间不增加。若有等待顾客，则第一个顾客离开队列，开始服务，所有等待的顾客的等待时间加到总的等待时间中去。模拟结束后，将总的等待时间除服务人数得平均等待时间。计量数模提高班专题五M/M/1模型流程图 初始化系统状态 第一个人接受服务 有人等吗 计算相关结果 结束 时间

29、推到该人服务结束 同时他离开系统Y 时间推到下一个顾客到来 服务人数 到了吗N 服务中会有 人来吗 产生新到达顾客NNYY计量数模提高班专题五 当模拟 100次得 Wq=9.0459,Ws=12.5612,而当模拟 10000次得 Wq=10.1294,Ws=13.1731,而从理论上应该是Ws=3/4/(1/3-1/4)=9,Wq=1/(1/3-1/4)=12，它们还是接近的。它们还是接近的。计量数模提高班专题五五、作业五、作业 1 两人约定下午1-3点在某处见面，设两人的到达分别服从分布U(1,3)和N(2,1)，且一人最多等另一人15分钟，求两人见面的概率。2 编程模拟顾客到达服从负指数分布（=10），服务时间分别服从正态分布（均值为9，方差为4）和定长时间（T=9）的排队过程，求顾客平均等待时间和逗留时间。3（选做）一列火车大约在下午1点离开A站，规律如下：离站时间 13：00 13：05 13：10 概率 0.7 0.2 0.1火车从A到B途中所需的平均时间为30分，有2分钟的标准差。如果你要赶的是这趟火车的下一站B，而你到达的时间分布为 时间 13：28 13：30 13：32 13：34 概率 0.3 0.4 0.2 0.1问你能赶上这列火车的概率是多少？4 用蒙特卡罗方法求积分8.13.0/sinxdxxI计量数模提高班专题五Thank you!