蒙特卡罗方法介绍和其建模应用视频课程0331专题培训课件

主要内容蒙特卡洛方法应用实例22009-B眼科病床安排应用3蒙特卡洛方法介绍12023/7/26南京信息工程大学蒙特卡洛方法介绍蒙特卡洛的优缺点及其适用范围2随机数的生成3蒙特卡洛方法概述12023/7/26南京信息工程大学MonteCarlo的起源设总计投了M个点，落入阴影部分N个，则不规则图形的面积为蒙特卡罗方法数学建模的解有两类精确解近似解自然现象有两类确定性现象不确定性现象随机现象模糊现象1随机投点试验求近似解引例---计算面积在一定条件下必然发生现象2023/7/26南京信息工程大学MonteCarlo的起源MonteCarlo方法：又称随机模拟方法，对研究的系统进行随机观察抽样,通过对样本值的统计分析，求得所研究系统的某些参数不同于一般数值计算方法，它是以概率统计理论为基础的一种方法，由于能够比较逼真地描述事物的特点及物理实验过程，解决一些数值方法难以解决的问题，因而该方法的应用领域日趋广泛它是在上世纪四十年代中期为了适应当时原子能事业的发展而发展起来，“曼哈顿计划”主持人之一、数学家：冯·诺伊曼用驰名世界的赌城—摩纳哥最大的城市MonteCarlo—来命名这种方法Monte-CarloMonacoJohnVonNeumann(1903-1957)

2023/7/26南京信息工程大学MonteCarlo方法的应用物理：核物理，热力学与统计物理，粒子输运问题等数学：多重积分、解微分方程、非线性方程组求解等工程领域：真空技术，水力学，激光技术等经济学领域：期权定价、项目管理、投资风险决策等其他领域：化学、医学，生物，生产管理、系统科学、公用事业等方面，随着科学技术的发展，其应用范围将更加广泛。2023/7/26南京信息工程大学MonteCarlo方法的基本思想

法国数学家ComtedeBuffon2023/7/26南京信息工程大学MonteCarlo方法的基本思想

每次投针试验实际上变成从两个均匀分布的随机变量中抽样，从而针线相交的概率为:

2023/7/26南京信息工程大学计算机模拟程序实现：functionpiguji=buffon(llength1,llength2,mm)%llength1,llength2分别表示1/2线的宽度和针的长度%mm是随机实验次数frq=0;xrandnum=unifrnd(0,llength1,1,mm);theta=unifrnd(0,pi,1,mm);forii=1:mmif(xrandnum(1,ii)<=(llength2*sin(theta(1,ii))))frq=frq+1;endendpiguji=(2*llength2/llength1)/(frq/mm)endbuffon(1,.6,1000)piguji=3.14622023/7/26南京信息工程大学建立统计模型，主要特征参量方面要与实际问题或系统相一致，问题的解对应于模型中随机变量的概率分布或其某些数字特征根据模型中各个随机变量的分布，在计算机上产生随机数，实现一次模拟过程所需的足够数量的随机数，进而进行随机模拟实验根据概率模型的特点和随机变量的分布特性，设计和选取合适的抽样方法，并对每个随机变量进行抽样（包括直接抽样、分层抽样、相关抽样、重要抽样等）按照所建立模型进行仿真试验、计算，求出问题的随机解统计分析模拟试验结果，给出问题的估计以及其精度估计。必要时，还应改进模型以降低估计方差和减少试验费用，提高模拟计算的效率。蒙特卡洛方法计算机模拟基本步骤2023/7/26南京信息工程大学蒙特卡洛模拟的理论基础大数定律---贝努里（Bernoulli）大数定律中心极限定理2023/7/26南京信息工程大学蒙特卡洛模拟的误差分析由中心极限定理可知：这表明，不等式近似地以概率1成立。上式也表明，收敛到的阶为O(n-1/2)。通常，蒙特卡罗方法的误差ε定义为：2023/7/26南京信息工程大学蒙特卡洛方法优缺点及其适用范围MonteCarlo方法及其程序结构简单产生随机数，通过大量简单重复抽样和简单计算计算相应的值收敛速度与问题维数无关MonteCarlo方法的收敛速度为O(n-1/2)，与一般数值方法相比很慢。因此，用MonteCarlo方法不能解决精确度要求很高的问题MonteCarlo方法误差只与标准差和样本容量n有关，而与样本所在空间无关，即MonteCarlo方法的收敛速度与问题维数无关，而其他数值方法则不然。MonteCarlo方法的适用性强MonteCarlo方法对多维问题的适用性在解题时受问题条件限制的影响较小例如：要计算s维空间中的任一区域Ds上的积分2023/7/26南京信息工程大学蒙特卡洛方法优缺点及其适用范围优点能够比较逼真地描述具有随机性质的事物的特点及物理实验过程受几何条件限制小收敛速度与问题的维数无关误差容易确定程序结构简单，易于实现缺点收敛速度慢。误差具有概率性。进行模拟的前提是各输入变量是相互独立的蒙特卡罗方法适应于随机性问题，对于非随机性问题，必须将它转化为随机问题2023/7/26南京信息工程大学Matlab中的随机数生成函数rand(m,n):(0,1)均匀分布unifrnd(a,b,m,n):(a,b)均匀分布randn(m,n):标准正态normrnd(mu,sigma,m,n):

正态分布exprnd(theta,m,n):指数分布

poissrnd(lamda,m,n):

泊松分布binornd(n,p,m,n):二项分布>>hist(poissrnd(3,1,1000),50)2023/7/26南京信息工程大学Matlab中的随机数生成函数random('name',A1,A2,A3,M,N)

name的取值可以是'norm'or'Normal''unif'or'Uniform''poiss'or'Poisson''beta'or'Beta''exp'or'Exponential''gam'or'Gamma''geo'or'Geometric''unid'or'DiscreteUniform'......2023/7/26南京信息工程大学一般分布随机数产生方法基本方法有如下三种：逆变换法复合抽样方法筛选法定理设随机变量U服从(0,1)上的均匀分布，则X=F-1(U)的分布函数为F(x)因此，要产生来自F(x)的随机数，只要先产生来自U(0,1)的随机数，然后计算F-1(u)即可。其步骤为

2023/7/26南京信息工程大学连续分布对于连续型分布，如果分布函数F(x)的反函数F－1(x)能够解析表示，则直接抽样方法是:2023/7/26南京信息工程大学指数分布为连续型分布，其一般形式如下：其分布函数为：

则(1)由U(0,1)抽取u

因为1-u也是(0,1)上均匀随机数，可将上式简化为例.指数分布2023/7/26南京信息工程大学Randnum=（-2）*log(unifrnd(0,1,1,1000));cdfplot(Randnum)例.产生指数分布的随机数2023/7/26南京信息工程大学例：病人视网膜疾病术后观察时间一般为[5,15]天，其人数分布如直方图所示，求满足该分布的术后观察天数序列。离散分布天数概率分布函数值区间53/1013/101162/1015/101277/10112/1013816/10128/1014911/10139/10151018/10157/10161115/10172/10171211/10183/1018139/10192/1019145/10197/10110154/1011112023/7/26南京信息工程大学离散分布x=[5:15];X=[];y=[3,2,7,16,11,18,15,11,9,5,4];yy=y(1);fori=2:length(y)yy=[yy,y(i)+yy(i-1)];endyy=yy/yy(end);fori=1:10000r=rand(1);d=size(find(r>=yy),2)+1;X(i)=x(d);Endhist(X,x)南京信息工程大学2023/7/26蒙特卡洛方法应用实例2023/7/26南京信息工程大学定积分的MC计算事实上，不少的统计问题最后都归结为定积分的近似计算问题!相对于其它方法，用MC方法比一般的数值方法有优点，主要体现在它的误差与维数m无关！下面考虑一个简单的定积分为了说明问题，这里介绍两种求的简单的MC方法。2023/7/26南京信息工程大学

方法简述：设a，b有限，0<f(x)<M，={(x,y):axb，0yM},并设(X,Y)是在上均匀分布的二维随机变量，其联合密度函数为则易见是中y=f(x)曲线下方的面积假设我们向中进行随机投点，则点落在y=f(x)下方的概率p,随机投点法2023/7/26南京信息工程大学若我们进行了n次投点，其中n0次点落入y=f(x)曲线下方，则用频率n0/n来估计概率p。即那么我们可以得到的一个估计2023/7/26南京信息工程大学注1随机投点法的思想简单明了，且每n次投点结果服从二项分布，故，其中注2可证是的无偏估计。若用估计的标准差来衡量其精度，则估计的精度的阶为。2023/7/26南京信息工程大学求解定积分的算例例7

计算定积分事实上，其精确解为用随机投点法求解：liti27(0,4,4,1000000)result=7.2336注：增加样本数目，可提高计算精度，但计算时间也会提高。functionresult=liti27(a,b,m,mm)%a是积分的下限%b是积分的上限%m是函数的上界%mm是随机实验次数frq=0;xrandnum=unifrnd(a,b,1,mm);yrandnum=unifrnd(0,m,1,mm);forii=1:mmif(cos(xrandnum(1,ii))+2>=yrandnum(1,ii))frq=frq+1;endendresult=frq*m*(b-a)/mm2023/7/26南京信息工程大学蒙特卡罗求解积分例

计算两条抛物线y=x2,x=y2所围面积.在正方形[0,1]×[0,1]区域投入2000个均匀随机点则随机点落入抛物线所围区域的概率为所求面积与正方形面积之比functionS=area1(N)ifnargin==0,N=2000;endX=rand(N,1);Y=rand(N,1);II=find(Y<=sqrt(X)&Y>=X.^2);m=length(II);S=m/N;x1=0:0.01:1;x2=1:-0.01:0;y1=sqrt(x1);y2=x2.^2;fill([x1,x2],[y1,y2],'c')南京信息工程大学2023/7/26蒙特卡罗求解积分南京信息工程大学例：冰淇淋锥的体积计算&冰淇淋锥含于体积=8的正方体22N个点均匀分布于正方体中,锥体中占有m个,则锥体与正方体体积之比近似为m:N2023/7/26蒙特卡罗求解积分南京信息工程大学functiondata=icecream(L)ifnargin==0,L=7;endN=10000;fork=1:LP=rand(N,3);x=2*P(:,1)-1;y=2*P(:,2)-1;z=2*P(:,3);R2=x.^2+y.^2;R=sqrt(R2);II=find(z>=R&z<=1+sqrt(1-R2));m=length(II);q(k)=8*m/N;enddata=[q;q-pi];蒙特卡罗方法计算冰淇淋锥的体积2023/7/26设g(x)是(a,b)上的一个密度函数，改写基本原理：对积分其中，X是服从g(x)的随机变量．可见，积分可以表示为X的函数的期望。由矩法，若有n个来自g(x)的观测值x1,…,xn，则可给出的一个矩估计:样本平均值法2023/7/26南京信息工程大学特别地，若a，b有限，可取g(x)为[a,b]上均匀分布．此时,设x1,…,xn是来自U(a,b)的随机数，则的一个估计为:具体步骤为:注可证是的无偏估计。一般而言，样本均值法要比随机投点法更有效。样本平均值法2023/7/26南京信息工程大学例9

计算定积分事实上，其精确解为样本平均值法求解：liti29(0,4,1000)result=7.1854liti29(0,4,10000)result=7.2153litti29(0,4,100000)result=7.2419注:增加样本数目，可提高计算精度，但计算时间也会提高。

求解定积分的算例functionresult=liti29(a,b,mm)%a是积分的下限%b是积分的上限%积分函数cos(x)+2%mm是随机实验次数sum=0;xrandnum=unifrnd(a,b,1,mm);forii=1:mmsum=sum+cos(xrandnum(1,ii))+2;endresult=sum*(b-a)/mm2023/7/26南京信息工程大学用蒙特卡洛法解非线性规划问题南京信息工程大学2023/7/26用蒙特卡洛法解非线性规划问题南京信息工程大学基本假设试验点的第j个分量xj服从[aj,bj]内的均匀分布．符号假设求解过程先产生一个随机数作为初始试验点,以后则将上一个试验点的第j个分量随机产生，其它分量不变而产生一新的试验点．这样，每产生一个新试验点只需一个新的随机数分量．当K>MAXK或P>MAXP时停止迭代．2023/7/26用蒙特卡洛法解非线性规划问题南京信息工程大学初始化:给定MAXK,MAXP;k=0,p=0,Q:大整数xj=aj+R(bj-aj)j=1,2,…,nj=0j=j+1,p=p+1P>MAXP?YNxj=aj+R(bj-aj)gi(X)≥0?i=1,2…nYNj<n?f(X)≥Q?YNX*=X,Q=f(X)k=k+1k>MAXK?YN输出X,Q,停止YN2023/7/26用蒙特卡洛法解非线性规划问题%mylp.mfunctionz=mylp(x)%目标函数z=2*x(1)^2+x(2)^2-x(1)*x(2)-8*x(1)-3*x(2);%转化为求最小值问题

在Matlab软件包中编程，共需三个Ｍ－文件:randlp.m,mylp.m,

lpconst.m.主程序为randlp.m.南京信息工程大学2023/7/26%randlp.m

function[sol,r1,r2]=randlp(a,b,n)%随机模拟解非线性规划debug=1;a=0;%试验点下界b=10;%试验点上界n=1000;%试验点个数r1=unifrnd(a,b,n,1);%n1阶的［a,b］均匀分布随机数矩阵r2=unifrnd(a,b,n,1);sol=[r1(1)r2(1)];z0=inf;fori=1:nx1=r1(i);x2=r2(i);lpc=lpconst([x1x2])；

iflpc==1z=mylp([x1x2])；

ifz<z0z0=z；sol=[x1x2]；

end

endend2023/7/26南京信息工程大学复杂概率模拟例:在我方某前沿防守地域，敌人以一个炮排（含两门火炮）为单位对我方进行干扰和破坏．为躲避我方打击，敌方对其阵地进行了伪装并经常变换射击地点．经过长期观察发现，我方指挥所对敌方目标的指示有50％是准确的，而我方火力单位，在指示正确时，有1/3的概率能毁伤敌人一门火炮，有1/6的概率能全部消灭敌人．现在希望能用某种方式把我方将要对敌人实施的1次打击结果显现出来，利用频率稳定性，确定有效射击(毁伤一门炮或全部消灭)的概率.2023/7/26南京信息工程大学复杂概率模拟分析：这是一个复杂概率问题，可以通过理论计算得到相应的概率.为了直观地显示我方射击的过程，现采用模拟的方式。1.问题分析需要模拟出以下两件事：[1]观察所对目标的指示正确与否模拟试验有两种结果，每一种结果出现的概率都是1/2．因此，可用投掷一枚硬币的方式予以确定,当硬币出现正面时为指示正确，反之为不正确．

2023/7/26南京信息工程大学复杂概率模拟[2]当指示正确时，我方火力单位的射击结果情况模拟试验有三种结果：毁伤一门火炮的可能性为1/3(即2/6)毁伤两门的可能性为1/6没能毁伤敌火炮的可能性为1/2(即3/6)这时可用投掷骰子的方法来确定：如果出现的是１、２、３三个点：则认为没能击中敌人；如果出现的是４、５点：则认为毁伤敌人一门火炮；若出现的是６点：则认为毁伤敌人两门火炮．2023/7/26南京信息工程大学复杂概率模拟2.符号假设i：要模拟的打击次数；k1：没击中敌人火炮的射击总数；k2：击中敌人一门火炮的射击总数；k3：击中敌人两门火炮的射击总数．E：有效射击(毁伤一门炮或两门炮)的概率2023/7/26南京信息工程大学复杂概率模拟3.模拟框图初始化:i=0,k1=0,k2=0,k3=0i=i+1骰子点数?k1=k1+1k2=k2+1k3=k3+1k1=k1+1i＜mm?E=(k2+k3)/mm停止硬币正面?YNNY1，2，34，562023/7/26南京信息工程大学复杂概率模拟functionliti26(p,mm)efreq=zeros(1,mm);randnum1=binornd(1,p,1,mm);randnum2=unidrnd(6,1,mm);k1=0;k2=0;k3=0;fori=1:mmifrandnum1(i)==0k1=k1+1;elseifrandnum2(i)<=3k1=k1+1;elseifrandnum2(i)==6k3=k3+1;elsek2=k2+1;endendefreq(i)=(k2+k3)/i;endnum=1:mm;plot(num,efreq)2023/7/26南京信息工程大学复杂概率模拟liti26(0.5,2000)liti26(0.5,20000)2023/7/26南京信息工程大学复杂概率模拟5.理论计算

模拟结果与理论计算近似一致，能更加真实地表达实际战斗动态过程．

2023/7/26南京信息工程大学报童的诀窍2023/7/26南京信息工程大学报童的诀窍问题报童售报：a(零售价)

>b(购进价)

>c(退回价)售出一份赚a-b；退回一份赔b-c

每天购进多少份可使收入最大？分析购进太多卖不完退回赔钱购进太少不够销售赚钱少应根据需求确定购进量每天需求量是随机的优化问题的目标函数应是长期的日平均收入每天收入是随机的存在一个合适的购进量等于每天收入的期望2023/7/26南京信息工程大学报童的诀窍建模

设每天购进n份，日平均收入为G(n)调查需求量的随机规律——每天需求量为r的概率f(r),r=0,1,2…准备求n使G(n)最大

已知售出一份赚a-b；退回一份赔b-c2023/7/26南京信息工程大学报童的诀窍求解将r视为连续变量2023/7/26南京信息工程大学报童的诀窍结果解释nP1P2取n使

a-b~售出一份赚的钱

b-c~退回一份赔的钱0rp2023/7/26南京信息工程大学报童的诀窍

然而，我们仍然难以应用，且不直观。鉴于报童问题的随机性特点，用仿真是非常有效的方法。为此，我们举一个具体例子，设a=0.5,b=0.3,c=0.15,且每天的报纸需求服从P(50)的分布，模拟见程序。由模拟可知每天宜订购51份报纸。2023/7/26南京信息工程大学系统的可靠性计算问题一个元件(或系统)能正常工作的概率称为元件(或系统)的可靠性系统由元件组成,常见的元件连接方式：串联并联12212023/7/26南京信息工程大学例13

设两系统都是由

4个元件组成,每个元件正常工作的概率为p,每个元件是否正常工作相互独立.两系统的连接方式如下图所示，比较两系统的可靠性.A1A2B2B1S1:2023/7/26南京信息工程大学A1A2B2B1例14

S2:2023/7/26南京信息工程大学例13'14'设两系统都是由

4个元件组成,每个元件的寿命服从平均寿命为θa1,θa2,θb1,θb2的指数分布，每个元件是否正常工作相互独立.两系统的连接方式如下图所示，求两系统寿命大于T=100的概率.functionRguji=liti213(t,thetaa1,thetaa2,thetab1,thetab2,mm)%t是要求系统生存的寿命%thetaa1是元件A1的数学期望%thetaa2是元件A2的数学期望%thetab1是元件B1的数学期望%thetab2是元件B2的数学期望%mm是随机实验次数frq=0;randnuma1=exprnd(thetaa1,1,mm);randnuma2=exprnd(thetaa2,1,mm);randnumb1=exprnd(thetab1,1,mm);randnumb2=exprnd(thetab2,1,mm);forii=1:mmif(randnuma1(1,ii)>t)&(randnuma2(1,ii)>t)pass1=1;elsepass1=0;endif(randnumb1(1,ii)>t)&(randnumb2(1,ii)>t)pass2=1;elsepass2=0;endif(pass1+pass2)>=1frq=frq+1;endend,Rguji=frq/mm系统1：2023/7/26南京信息工程大学functionRguji=liti214(t,thetaa1,thetaa2,thetab1,thetab2,mm)%t是要求系统生存的寿命%thetaa1是元件A1的数学期望%thetaa2是元件A2的数学期望%thetab1是元件B1的数学期望%thetab2是元件B2的数学期望%mm是随机实验次数frq=0;randnuma1=exprnd(thetaa1,1,mm);randnuma2=exprnd(thetaa2,1,mm);randnumb1=exprnd(thetab1,1,mm);randnumb2=exprnd(thetab2,1,mm);forii=1:mmif(randnuma1(1,ii)>t)|(randnumb1(1,ii)>t)pass1=1;elsepass1=0;endif(randnuma2(1,ii)>t)|(randnumb2(1,ii)>t)pass2=1;elsepass2=0;endif(pass1*pass2)==1frq=frq+1;endendRguji=frq/mm系统2：2023/7/26南京信息工程大学坎雷渔业公司问题2023/7/26南京信息工程大学0.坎雷渔业公司问题

克林特坎雷经营着Massachusetts一家拥有50条鳕鱼捕捉船的渔业公司，每个工作日，渔船早上离港，中午作业完毕，每次每条船能捕鱼3500单位。有许多港口都可以停靠并出售鳕鱼。每个港口每条的价格是不确定的，并且变化很大；而且港口之间价格也不一样，另外，每个港口的需求量是有限的，如果一条船比别的船晚到一个港口，那么它的鱼就卖不出去，要倒进海洋中。2023/7/26南京信息工程大学1.坎雷渔业公司问题简化简化问题假设渔业公司只有一条船，每次出海的成本为10,000美元，每次出海捕鱼3500单位，两个港口[格洛斯特，岩石港]可以停靠：格洛斯特是鳕鱼的集散地，价格一直稳定在每单位3.25美元，需求几乎是无限的岩石港比较小，价格较高但波动比较大；岩石港的价格服从均值为3.65标准差为0.20的正态分布,需求量服从表1的离散分布；并且我们假设两个港口之间的价格，需求量之间是相互独立的，每天渔船只能在一个港口停靠并出售它的鳕鱼。而岩石港每天的价格事先并不知道。坎雷想挣得尽可能大的利润，哪一个港口停靠更好？日需求单位01,0002,0003,0004,0005,0006,000概率0.020.030.050.080.330.290.2表1岩石港鳕鱼日需求分布表2023/7/26南京信息工程大学1.坎雷渔业公司问题简化渔船在格洛斯特港停靠的利润G为：但是，停靠在岩石港的利润计算出P没这简单，因为价格和需求量都是不确定的，每天的利润是一个随机变量，为了决定选择哪个港口，下面的问题将是很有帮助的：(a)使用岩石港日利润的概率分布大概是什么形状？(b)使用岩石港利润高于使用格洛斯特港利润的概率是多少？(c)使用岩石港亏本的概率是多少？(d)使用岩石港日利润的期望值是多少？(e)使用岩石港日利润的标准差是多少？2023/7/26南京信息工程大学2.坎雷渔业公司初步分析

定义两个随机变量：PR=岩石港的鳕鱼价格：PR~N(3.65,0.202)D=停靠岩石港坎雷面临的需求量：D的分布如表1记F为停靠岩石港的日利润，那么有：上面5个问题更简洁的表达为：(a)F的概率密度函数是什么形状？(b)P（F>1375）是多少？(c)P（F<0）是多少？(d)F的期望值是多少？(e)F的标准差是多少？注：F是两个随机变量乘积的函数，它的分布用前面的方面不易求得，我们必须运用新的方法2023/7/26南京信息工程大学3.坎雷渔业公司的随机模拟模型模拟实验很简单，首先我们做出表2：第二列[岩石港的需求]：服从表1所示的离散分布的随机数第三列为第二列与3,500中的小者第四列[岩石港的鳕鱼价格]：均值为3.65，标准差为0.2正态分布的随机数最后一列：日利润数据，利用前面公式计算可得日期岩石港需求D(lbs.)Min(D,3,500)(lbs.)岩石港价格PR($/lb.)日利润F($)12…10000表2坎雷渔业公司电脑模拟表2023/7/26南京信息工程大学4.样本数据（即随机数）生成产生服从离散分布D的样本数据（即随机数）：将[01]进行分割，每个区间对应离散分布的可能值，并且每个区间的长度就是该概率值用随机数发生器产生服从[01]间均匀分布的随机数对于每一个[01]的随机数，赋予它所位于区间的长度对应的离散分布所取的值产生服从连续分布PR的样本数据：利用随机数发生器产生一系列服从[0,1]均匀分布的随机数对于第1步中产生的随机数u，计算x值使得分布函数值等于u，即利用F(x)=u求出x,那么y就是我们要的随机数将上述结果填入表2中，我们就得到了分析样本数据2023/7/26南京信息工程大学5.样本数据分析利用模拟结果中的数据，可以回答前面的问题：(a)F的概率密度是什么形状？给出10000个样本观测值的经验分布函数和直方图(b)P(F>1375)的估计值是多少？计算出F大于1375的频率,它即为该估计值我们估计停靠岩石港利润更高的概率，这个结果将支持我们是否选择岩石港(c)P(F<0)的估计值是多少？计算出F小于0频率，它为所求估计值(d)F的期望的估计值是多少？样本均值即为F的期望的估计值2023/7/26南京信息工程大学历年赛题随机问题和排队问题在数模竞赛中多次出现过。如：1997A题：零件的参数设计2009B题：眼科病床的合理安排2013A题：车道被占用对城市道路通行能力的影响2014年美赛题目：除非超车否则靠右行驶2023年7月26日排队问题随机模拟2023/7/26南京信息工程大学排队现象到达的顾客要求服务内容服务机构1.不能运转的机器修理修理技工2.修理技工领取修配零件发放修配零件的管理员3.病人诊断或动手术医生(或包括手术台)4.电话呼唤通话交换台5.文件稿打字打字员6.提货单提取存货仓库管理员7.到达机场上空的飞机降落跑道8.驶入港口的货船装(卸)货装(卸)货码头(泊位)9.上游河水进入水库放水，调整水位水闸管理员10.进入我方阵地的敌机我方高射炮进行射击我方高射炮在日常生活中，人们会遇到各种各样的排队问题：2023/7/26南京信息工程大学排队论排队论(QueueingTheory)：通过研究各种服务系统等待现象中的概率特征，从而解决服务系统最优设计与最优控制的一种理论，是运筹学的一个重要分支由于顾客到达和服务时间的随机性，排队过程通常是一个随机过程，排队论又称“随机服务系统理论”排队论发源于上世纪初，1909年丹麦的哥本哈根电话公司A.K.埃尔浪(Erlang)在热力学统计平衡概念的启发下予以解决了美国贝尔电话公司遇到的通话线路与电话用户呼叫的数量关系问题2023/7/26南京信息工程大学排队系统排队服务过程---三个基本部分顾客输入过程排队结构与排队规则服务机构与服务规则2023/7/26南京信息工程大学排队系统顾客输入过程：顾客源(总体)：有限/无限;顾客到达方式：逐个/逐批;(仅研究逐个情形)顾客到达间隔：随机型/确定型;顾客前后到达是否独立：相互独立/相互关联；输入过程是否平稳：平稳/非平稳；(仅研究平稳性)2023/7/26南京信息工程大学排队系统排队结构与排队规则：顾客排队方式：等待制/损失制/混合制排队系统容量：有限制/无限制;排队队列数目:单列/多列;是否中途退出:允许/禁止;是否列间转移:允许/禁止;(仅研究禁止退出和转移的情形)2023/7/26南京信息工程大学排队系统等待时间+服务时间无限排队损失制排队混合制排队系统排队有限排队队长有限等待时间有限逗留时间有限等待制系统2023/7/26南京信息工程大学排队系统服务机构服务台(员)数目：单个/多个服务台(员)排列形式：并列/串列/混合服务台(员)服务方式：逐个/逐批服务时间分布：随机型/确定型服务时间分布是否平稳:平稳/非平稳112c

…12c

…12c

…服务规则：先到先服务(FCFS)；后到先服务(LCFS)；随机服务(RS)；优先服务2023/7/26南京信息工程大学排队模型1.排队模型表示方法：D.G.Kendall(1953)表示法：X/Y/ZX顾客到达间隔时间分布Y服务台（员）服务时间分布Z服务台（员）个数（单个或多个并列）国际排队论标准化会议(1971)表示法：X/Y/Z/A/B/CA系统容量限制B顾客源（总体）数目C服务规则（FCFS，LCFS等）约定：如略去后三项，即指X/Y/Z/∞/∞/FCFS的情形2023/7/26南京信息工程大学排队模型2.到达间隔和服务时间典型分布：M：负指数分布(M是Markov的字头，因为负指数分布具有无记忆性，即Markov性)D：确定型(deterministic)Ek：k阶爱尔朗(erlang)分布GI：一般相互独立(generalindependent)的时间间隔的分布G：一般(general)服务时间的分布2023/7/26南京信息工程大学排队模型3.排队模型示例：M/M/S/∞

：表示输入过程是Poisson流,服务时间服从负指数分布,系统有S个服务台平行服务,系统容量为无穷的等待制排队系统.M/G/1/∞：表示输入过程是Poisson流，顾客所需的服务时间为独立、服从一般概率分布，系统中只有一个服务台，容量为无穷的等待制系统.GI/M/1/∞：表示输入过程为顾客独立到达且相继到达的间隔时间服从一船概率分布，服务时间是相互独立、服从负指数分布，系统中只有一个服务台，容量为无穷的等待制系统D/M/S/K：表示相继到达的间隔时间独立、服从定长分布、服务时间相互独立、服从负指数分布，系统中有S个服务台平行服务，容量为K的混合制系统.2023/7/26南京信息工程大学系统参数---运行状态参数系统状态N(t)：指排队系统在时刻t时的全部顾客数N(t)，包括“排队顾客数”和“正被服务顾客数”系统状态的可能值如下：系统容量无限制，N(t)=0，1，2，…系统容量为N时,N(t)=0,1,2,…,N服务台个数为c/损失制,N(t)=0,1,2,…,c系统处于这些状态的概率一般是随时间t变化的，所以在时刻t、系统状态为n的概率可以用Pn(t)表示。2023/7/26南京信息工程大学系统参数---运行状态参数系统状态概率：瞬态概率Pn(t)：表示时刻系统状态N(t)=n的概率;稳态概率Pn：一般，排队系统运行了一定长的时间后，系统状态的概率分布不再随时间t变化，即初始时刻（t=0）系统状态的概率分布(Pn(0)，n≥0）的影响将消失求稳态概率Pn时，不需要求t→∞时Pn(t)的极限，而只需令导数P'(t)=0即可。2023/7/26南京信息工程大学系统参数---运行指标参数系统运行指标参数---评价排队系统的优劣队长与等待队长：队长:系统中的顾客数（n），期望值为：等待队长:系统中排队等待服务的顾客数，期望值为：2023/7/26南京信息工程大学系统参数---运行指标参数系统运行指标参数---评价排队系统的优劣逗留时间与等待时间逗留时间:指一个顾客在系统中全部停留时间；期望值，记为Ws等待时间:指一个顾客在系统中排队等待时间；期望值，记为Wq［逗留时间］＝［等待时间］＋［服务时间］其他相关指标:忙期:指从顾客到达空闲服务机构起到服务机构再次空闲的时间长度忙期服务量：指一个忙期内系统平均完成服务的顾客数损失率:指顾客到达排队系统，未接受服务而离去的概率服务强度：=/c

2023/7/26南京信息工程大学排队系统时间参数分布规律到达间隔的分布和服务时间的分布

经验分布泊松流负指数分布爱尔朗分布2023/7/26南京信息工程大学经验分布例1.大连港大港区载货500吨以上船舶到达(不包括定期到达的船舶)逐日记录：将数据整理成船舶到达数的分布表：可以计算出：平均到达率=到达总数/总天数=1271/365=3.48(艘/天)船舶到达数n频数频率(%)0120.0331430.1182640.1753740.2034710.1955490.1346260.0717190.052840.011920.00510以上10.003合计3651.0002023/7/26南京信息工程大学经验分布实际中测定相继到达时间间隔的方法以τi表示第i号顾客到达的时刻，以si表示对它的服务时间，这样可算出相继到达的间隔时间ti(ti=τi+1-τi)和排队等待时间wi，它们的关系如下：相继到达时间间隔等待时间：2023/7/26南京信息工程大学经验分布例2某服务机构是单服务台，先到先服务，对41个顾客记录到达时刻τ和服务时间s(单位为分钟)，在表中以第1号顾客到达时刻为0，对所有顾客的全部服务时间为127分钟。将原始记录整理成下表。到达间隔/分钟次数162103846536272819110以上1合计40服务时间/分钟次数1102103745546271819以上1合计412023/7/26南京信息工程大学经验分布例2平均间隔时间=142/40=3.55(分钟/人)平均到达率=41/142=0.28(人/分钟)平均服务时间=127/41=3.12(分钟/人)平均服务率=41/127=0.32(人/分钟)2023/7/26南京信息工程大学泊松流设表示在时间区间内到达的顾客数，令表示在时间区间内有个顾客到达的概率，即：当合于下列三个条件时，我们说顾客的到达形成泊松流：

1.无后效性(Markov性):在不相重叠的时间区间内顾客到达数是相互独立2.平稳性:对充分小的，在时间区间内有1个顾客到达的概率与t无关，而与区间长度成正比，即

是常数，表示单位时内有一个顾客到达的概率，称为概率强度。3.普通性：对于充分小的，在时间区间内有2个或2个以上顾客到达的概率极小，以至于可以忽略，即：2023/7/26南京信息工程大学泊松流表示长为t的时间区间内到达n个顾客的概率，随机变量服从泊松分布，它的数学期望和方差分别是：2023/7/26南京信息工程大学负指数分布负指数分布的定义：如果随机变量T的概率密度是：则称T服从负指数分布。它的分布函数是：数学期望为方差为标准差为2023/7/26南京信息工程大学负指数分布负指数分布的性质（无记忆性或马尔柯夫性）(1)由条件概率公式容易证明该性质说明：一个顾客到来所需的时间与过去顾客到达的时间无关，也就是说该顾客是随机地到达。(2)当输入过程是泊松流时，那么顾客相继到达的间隔时间T必然服从负指数分布。这是因为对于泊松流，在区间内至少有1个顾客到达的概率是又可表示为

结合负指数分布函数即得所证。2023/7/26南京信息工程大学令，则ρ称为服务强度。负指数分布相继到达时间间隔为负指数分布与到达为泊松流的等价性顾客相继到达的间隔时间独立且为同负指数分布(密度函数为)与输入过程为泊松流(参数为λ)是等价的。对于泊松流，λ表示单位时间平均到达的顾客数，1/λ表示相继顾客到达平均间隔时间。服务时间v的分布对顾客的服务时间也就是在忙期相继离开系统的两顾客的间隔时间，有时也服从负指数分布。设它的分布函数和密度分别是其中表示单位时间能被服务完成的顾客数，称为平均服务率，表示对顾客的平均服务时间。2023/7/26南京信息工程大学爱尔朗分布爱尔朗分布的定义：设是k个相互独立的随机变量，服从相同参数的负指数分布，那么的概率密度是：称T服从k阶爱尔朗分布，其均值和方差分别为例子：串列的k个服务台，每台服务时间相互独立，服从相同的负指数分布(参数kμ)，那么一顾客走完这k个服务台总共所需要服务时间就服从k阶爱尔朗分布。2023/7/26南京信息工程大学爱尔朗分布爱尔朗分布的意义当k=1时，爱尔朗分布化为负指数分布，可看成是一种完全随机的分布；当k增大时，爱尔朗分布的图形逐渐变为对称的当k≥30时爱尔朗分布近似于正态分布；k→∞时,Var[T]→0，这时爱尔朗分布化为确定型分布。一般k阶爱尔朗分布可看成完全随机与完全确定的中间型，能对现实世界提供更为广泛的适应性。2023/7/26南京信息工程大学Little（利特尔）公式对一般排队系统，性能指标之间均有下式成立：其中有效到达率为：2023/7/26南京信息工程大学蒙特卡洛排队模拟[1]系统的假设：（1）顾客源是无穷的；（2）排队的长度没有限制；（3）到达系统的顾客按先后顺序依次进入服务，“先到先服务”。在某商店有一个售货员，顾客陆续来到，售货员逐个地接待顾客．当到来的顾客较多时，一部分顾客便须排队等待，被接待后的顾客便离开商店．设：1．顾客到来间隔时间服从参数为10的指数分布．2．对顾客的服务时间服从［4，15］上的均匀分布．3．排队按先到先服务规则，队可以无限制长．假定一个工作日为8小时，时间以分钟为单位。[1]模拟一个工作日内完成服务的个数及顾客平均等待时间t．[2]模拟100个工作日，求出平均每日完成服务的个数及每日顾客的平均等待时间。2023/7/26南京信息工程大学蒙特卡洛排队模拟c1b1c3c4c5c2e1b2e2b3e3b4e4b5t建立排队模型：ci=ci-1+xibi=max(ci,ei-1)ei=bi+yiwi=bi-ciw：总等待时间；ci：第i个顾客的到达时刻；bi：第i个顾客开始服务时刻；ei：第i个顾客服务结束时刻；xi：第i-1个顾客与第i个顾客之间到达的间隔时间；yi：对第i个顾客的服务时间。2023/7/26南京信息工程大学蒙特卡洛排队模拟初始化：令i=1，ei-1=0，w=0产生间隔时间随机数xi~参数为10的指数分布ci=xi,bi=xi

产生服务时间随机数yi~[4,15]的均匀分布ei=bi+yi累计等待时间：w=w+bi-ci准备下一次服务：i=i+1产生间隔时间随机数xi~参数为10的指数分布ci=ci-1+xi

确定开始服务时间：bi=max(ci,ei-1)bi＞480?YNi=i-1,t=w/i输出结果：完成服务个数：m=i

平均等待时间：t停止[2]模拟100日ToMatlab(simu2)[1]模拟一日ToMatlab(simu1)2023/7/26南京信息工程大学蒙特卡洛排队模拟function[i,t]=simu1i=1;e=0;w=0;x(i)=random('exp',10);c(i)=x(i);b(i)=x(i);while(b(i)<=480)y(i)=11*rand(1)+4;e(i)=b(i)+y(i);w=w+b(i)-c(i);i=i+1;x(i)=random('exp',10);c(i)=c(i-1)+x(i);b(i)=max(c(i),e(i-1));endi=i-1;t=w/i;fprintf('一天完成服务人数：%f位\n',i);fprintf('每位顾客平均等待时间：%f分钟\n',t);functionsimu2M(1)=0;T(1)=0;fori=1:100[M(i),T(i)]=simu1;endmean_M=mean(M);mean_T=mean(T);fprintf(‘一百天完成服务人数：%f位\n',mean_M);fprintf(‘一百天每位顾客平均等待时间：%f分钟\n',mean_T);2023/7/26南京信息工程大学眼科病床的合理安排4全国大学生数学建模竞赛2009年B题2023/7/26南京信息工程大学一、赛题医院就医排队是大家都非常熟悉的现象，它以这样或那样的形式出现在我们面前，例如，患者到门诊就诊、到收费处划价、到药房取药、到注射室打针、等待住院等，往往需要排队等待接受某种服务我们考虑某医院眼科病床的合理安排的数学建模问题。该医院眼科门诊每天开放，住院部共有病床79张。该医院眼科手术主要分四大类：白内障、视网膜疾病、青光眼和外伤。附录中给出了2008年7月13日至2008年9月11日这段时间里各类病人的情况。白内障手术较简单，而且没有急症。目前该院是每周一、三做白内障手术，此类病人的术前准备时间只需1、2天。做两只眼的病人比做一只眼的要多一些，大约占到60%。如果要做双眼是周一先做一只，周三再做另一只。2023/7/26南京信息工程大学外伤疾病通常属于急症，病床有空时立即安排住院，住院后第二天便会安排手术。其他眼科疾病比较复杂，有各种不同情况，但大致住院以后2-3天内就可以接受手术，主要是术后的观察时间较长。这类疾病手术时间可根据需要安排，一般不安排在周一、周三。由于急症数量较少，建模时这些眼科疾病可不考虑急症。该医院眼科手术条件比较充分，在考虑病床安排时可不考虑手术条件的限制，但考虑到手术医生的安排问题，通常情况下白内障手术与其他眼科手术（急症除外）不安排在同一天做。当前该住院部对全体非急症病人是按照FCFS（Firstcome,Firstserve）规则安排住院，但等待住院病人队列却越来越长，医院方面希望你们能通过数学建模来帮助解决该住院部的病床合理安排问题，以提高对医院资源的有效利用。2023/7/26南京信息工程大学问题一：试分析确定合理的评价指标体系，用以评价该问题的病床安排模型的优劣。问题二：试就该住院部当前的情况，建立合理的病床安排模型，以根据已知的第二天拟出院病人数来确定第二天应该安排哪些病人住院。并对你们的模型利用问题一中的指标体系作出评价。问题三：作为病人，自然希望尽早知道自己大约何时能住院。能否根据当时住院病人及等待住院病人的统计情况，在病人门诊时即告知其大致入住时间区间。问题四：若该住院部周六、周日不安排手术，请你们重新回答问题二，医院的手术时间安排是否应作出相应调整?问题五：有人从便于管理的角度提出建议，在一般情形下，医院病床安排可采取使各类病人占用病床的比例大致固定的方案，试就此方案，建立使得所有病人在系统内的平均逗留时间（含等待入院及住院时间）最短的病床比例分配模型。2023/7/26南京信息工程大学【附录】2008-07-13到2008-09-11的病人信息序号类型门诊时间入院时间第一次手术时间第二次手术时间出院时间1外伤2008-7-132008-7-142008-7-15/2008-7-192视网膜疾病2008-7-132008-7-252008-7-27/2008-8-83白内障2008-7-132008-7-252008-7-28/2008-7-314视网膜疾病2008-7-132008-7-252008-7-27/2008-8-45青光眼2008-7-132008-7-252008-7-27/2008-8-56视网膜疾病2008-7-132008-7-262008-7-29/2008-8-117白内障(双眼)2008-7-132008-7-262008-7-282008-7-302008-8-28视网膜疾病2008-7-142008-7-262008-7-29/2008-8-69白内障(双眼)2008-7-142008-7-262008-7-282008-7-302008-8-110白内障2008-7-142008-7-262008-7-28/2008-7-30…………………348外伤2008-9-32008-9-42008-9-5/2008-9-10349外伤2008-9-42008-9-52008-9-6/2008-9-112023/7/26南京信息工程大学【附录】2008-07-13到2008-09-11的病人信息序号类型门诊时间入院时间第一次手术时间第二次手术时间出院时间1视网膜疾病2008-8-152008-8-292008-8-31//2视网膜疾病2008-8-162008-8-292008-8-31//3白内障(双眼)2008-8-192008-9-12008-9-82008-9-10/4………………579外伤2008-9-92008-9-102008-9-11/1白内障(双眼)2008-8-30////2视网膜疾病2008-8-30////…………………101视网膜疾病2008-9-11////102视网膜疾病2008-9-11////2023/7/26南京信息工程大学建模准备---数据分析数据分析：门诊情况等待入院等待手术术后恢复各类就诊病人占总就诊病人比例各类病人每天出院情况

原方案分析：排队机制：先到先服务，不存在插队情况，对各类病人都公平术前等待时间长，入院等待时间长。方案不合理的地方。理论分析认为：每天就诊人数与出院人数基本相等，可以做到队长不变的理想情况。2023/7/26南京信息工程大学建模准备---数据分析白内障病人入院后1-2天即可做手术，但是统计该指标得白内障单眼病人需要等待2.38天，白内障双眼病人需要等待3.63天，远大于1-2天。关键信息：术前住院时间过长，病床有效利用率不高2023/7/26南京信息工程大学建模准备---分布检验分布检验：病人到达人数：服从Poisson分布（画频数（率）直方图，分布检验，分布参数提取）术后住院时间：服从正态分布or

Г分布等（估计5类病人住院时间的均值与方差，分布检验）2023/7/26南京信息工程大学病人到达人数分布拟合优度法检验（卡方检验或F检验）

2023/7/26南京信息工程大学术后恢复时间分布2023/7/26南京信息工程大学第一问问题一：试分析确定合理的评价指标体系，用以评价该问题的病床安排模型的优劣。评价指标可分两类：效率指标和公平性指标

效率指标——平均术前住院时间，或病床有效利用率。也可以是等待入院时间等。公平性指标——“插队人数比例”。（容易被忽略也可以有不同理解）这两个指标可以有各种不同的定义，其合理性是评分依据，注意：不必要引入过多实质上相同的指标原方案的唯一准则，自然有其合理性。2023/7/26南京信息工程大学第二问问题二：试就该住院部当前的情况，建立合理的病床安排模型，以根据已知的第二天拟出院病人数来确定第二天应该安排哪些病人住院。并对你们的模型利用问题一中的指标体系作出评价。本问主要考核能否给出一个相对合理的病床安排模型，主要目标为：提高病床有效利用率以及提高公平度。就提高病床有效利用率而言，病人术后住院时间是一个不可优化的量，所以只能在术前等待时间上作文章。经对问题的分析可知：对白内障病人的入院时间加以限制成为提高效率的必然选择。2023/7/26南京信息工程大学本问主要解决方法是仿真方法，大致可分为“先仿真，再优化”与“边仿真，边优化”两类：前者是先确定若干种住院规则，然后根据仿真统计结果选出较优规则；后者是先确定一个优化原则，然后在仿真时，对每一个排队病人按照该优化原则决定住院先后。显然后者要更好一些。2023/7/26南京信息工程大学优化的入院规则：规则一：按照病人的可入院时间，在确定一周中某天哪个病人入院时，病人采用FCFS服务制（白内障除外）；规则二：优先排急症病人；规则三：白内障病人的术前住院时间为1—2天；视网膜疾病、青光眼两类病人术前住院时间为2-3天；外伤手术有空床即安排入院，第二天手术规则四：选择一周为一周期；规则五：周六、周日除了外伤病人以外，优先安排白内障（双眼）病人入院根据以上准则及原来的方案通过计算机模拟给出一段时间的模拟，分别计算评价指标，说明方案的优劣。各类病人的平均住院时间缩短1-2天，公平性下降5%左右（模拟值）。2023/7/26南京信息工程大学仿真注意事项如何产生就诊人的队列（计算插队时用）如何处理出院分布如何根据出院情况及入院准则安排入院计算总等待人数计算插队人数计算床位有效利用率2023/7/26南京信息工程大学第三问问题三：作为病人，自然希望尽早知道自己大约何时能住院。能否根据当时住院病人及等待住院病人的统计情况，在病人门诊时即告知其大致入住时间区间。此问希望学生给出一个满足一定置信度（例如：90%）的预约住院时间区间，区间长度越短越好。一种自然的想法是通过同类病人术后住院时间的概率分布从理论上得到这一区间，但这样做的一个困难是已处于术后住院状态的该类病人的继续住院时间不服从同一分布，从而将该类病人（含已住院与未住院）的预计住院时间求和后的随机变量的分布不知道。2023/7/26南京信息工程大学预约住院时间置信区间设当前时刻为T0，当前排队人数为P，预计住院时刻为T，病人每日出院人数的统计平均值为α，则设一个已出院病人实际住院时刻为T1，通过仿真统计一段时间内所有病人的根据90%的置信度确定两个阈值从而得到当前病人的预计住院时间区间为。

2023/7/26南京信息工程大学预约住院时间置信区间现有等待人数条件下：2023/7/26南京信息工程大学第四问问题四：若该住院部周六、周日不安排手术，请你们重新回答问题二，医院的手术时间安排是否应作出相应调整?若仍采用“一三方案”，效率较低，通过分析可以发现主要原因是对视网膜与青光眼病人而言，会造成病床使用效率降低。通过有限种方案的仿真计算比较可知，采用“二四方案”或“三五方案”可使病床使用效率有所提高。前者效率＋公平总体效果较好，后者效率较高，但公平性较差一些。模拟与计算：模拟一、三，二、四方案，三、五方案计算三种方案的两个评价指标值，说明方案的优劣结果：一、三方案较差，三、五方案效率最高，公平性不好，二、四方案二者兼而有之。2023/7/26南京信息工程大学第五问问题五：有人从便于管理的角度提出建议，在一般情形下，医院病床安排可采取使各类病人占用病床的比例大致固定的方案，试就此方案，建立使得所有病人在系统内的平均逗留时间（含等待入院及住院时间）最短的病床比例分配模型。主要有三种模型：仿真计算模型：床位分配只有有限种组合情形，可以通过穷举仿真方法得到各种组合的评价指标统计值，再比较得到最佳组合方案。此方案计算量较大，且模型通用性有一定局限。服务强度平衡模型：当各分类系统的服务强度相等时，效果最佳。可以通过建立条件极值模型，利用拉格朗日方法证明这一结论。排队论近似模型：通过经验公式将M/G/K系统近似为M/M/K系统，然后利用排队论的现成结论写出优化模型。2023/7/26南京信息工程大学近似模型的参数确定较困难（如各类病人的入院等待时间）。需要通过仿真实现。最好用仿真模型比如：（M/M/k)2023/7/26南京信息工程大学综合评述数据检验是本问题中必须做的，但被许多参赛队所忽略，从而意外成为区分点之一。公平性指标被许多人忽略，反映出对问题本质认识不到位。效率指标也可以适当精简。优化模型的多样性是本题目最大的亮点，涌现许多意料之外的解法。2023/7/26南京信息工程大学重要的数据变量T病人排队数组C门诊时间B入院时间W空等开始时间O手术开始时间E术后出院时间蒙特卡洛仿真2023/7/26南京信息工程大学例：建立评价指标定义指标1：平均术前住院时间M1M2=∑(Oi-Wi)/∑i定义指标2：公平性M2M2=“延期入院”人数/总人数总指标：

M=x1*M1+x2*M2参数x1,x2按要求设定。2023/7/26南京信息工程大学例：建立评价指标定义指标1：入院等待不满