计算机模拟版本3.ppt

,单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,数学建模培训,计算机模拟,主讲：何仁斌,概述,计算机模拟是用计算机对,系统,的结构和行为进行动态演示，以评价或预测系统的行为效果，为决策者提供信息的一种方法。它是解决较复杂的实际问题的一条有效途径。,计算机模拟也可以说是用,计算机程序,直接建立真实系统的模型，并通过计算了解系统随时间变化的行为或特征。,应用领域：航空、机电、冶金、社会经济、交通运输、生态系统等。,计算机模拟分为,连续系统,模拟和,离散系统,模拟。,状态随时间连续变化的系统称为连续系统。通常该系统的模型一般可以用微分方程的形式表达，通过一些物理机理推导出来。模拟结果往往是近似的。,例如，追逐问题，浓度问题。,一、连续系统,(x,1,(t),x,n,(t,),1,、追逐问题,AB,BC,CD,DA,状态变量,A(,x,(,t,),y,(,t,),假设：,v=,1,m/s,追逐的方向是对准追逐目标,A(0,10),B(10,10),D(0,0),C(0,10),试确定四人追逐的运动轨迹。,A(,x,(,t+,t,),y,(,t,+,t,),模拟方法：,1.,建立平面直角坐标系；,2.,以时间间隔,t,进行采样，并且计算各个,时刻下的状态：,A:(,x,1,(,t,),y,1,(,t,)(,x,1,(,t,+,t,),y,1,(,t,+,t,),B:(,x,2,(,t,),y,2,(,t,)(,x,2,(,t,+,t,),y,2,(,t,+,t,),(,x,1,(,t,+,t,),y,1,(,t,+,t,),（由几何原理）,(,x,1,(,t,)+,v,t,cos(,),y,1,(,t,)+,v,t,sin(,),x,1,y,1,x,2,y,2,在,t,时刻下,3.,选取足够小的,t,，模拟到任意两人的距,离小于,v,t,为止。,4.,连接四人在各时刻下的位置，就得到所求的,运动轨迹。,Matlab,(simu2.m),v=1;,dt,=0.05;,d=20;%x=x(1),x(2),x(3),x(4),x(5),x(6),x(7),x(8),x=0 0 0 10,10,10,10,0;,x(9)=x(1);,x(10)=x(2);,hold on,axis(,equal,);,axis(0 10 0 10);,while,(d,0.1),for,i=1:2:7,d=sqrt(x(i)-x(i+2)2+(x(i+1)-x(i+3)2);,x(i,)=,x(i)+v,*,dt,*(x(i+2)-x(i)/d;,x(i+1)=x(i+1)+v*,dt,*(x(i+3)-x(i+1)/d;,plot(x(i),x(i+1),.,),end,x(9)=x(1);x(10)=x(2);,end,hold on,A,B,C,D,某军的一导弹基地发现正北方向,120km,处海面上有敌艇一艘以,90km/h,的速度向正东方向行驶。该基地立即发射导弹跟踪追击敌艇，导弹速度为,450km/h,，自动导航系统使导弹在任一时刻都能对准敌艇。试问导弹在何时何处击中敌艇？,2,、导弹跟踪问题,P,(,x,(,t,),y,(,t,),O,x,M,(,v,e,t,H,),y,A,(0,H,),正东方向,正北方向,H,=120,（千米）,km/h,km/h,敌,导,(,x,k,(,t,),y,k,(,t,),simu1.m,stest1.m,huatu.m,微分方程建模,计算机模拟,消去,t,：,微分方程建模,1,）,2,）,x=dsolve(D2x*(120-y)/sqrt(Dx)2+1)=90/450,Dx(0)=0,x(0)=0,y),微分方程求解,simplify(x,),ans,=(45*120(1/5)*(y-120)(4/5)-1800*(-1)(4/5)+(-1)(3/5)*202500(1/5)*(y-120)(6/5)/(72*(-(-1)(3/5)(1/2),-(45*120(1/5)*(y-120)(4/5)-1800*(-1)(4/5)+(-1)(3/5)*202500(1/5)*(y-120)(6/5)/(72*(-(-1)(3/5)(1/2),y=120,eval(x,),ans,=,25.0000+0.0000i,-25.0000-0.0000i,如何决定导弹位置,P,2,(x,2,y,2,)?,计算机模拟,1),当,t=0,时，导弹位置,O(0,0);,敌艇位置,A(0,120),；,2),当,t=,时，导弹位置,P,1,(x,1,y,1,);,敌艇位置,M,1,(90,120),；,3),当,t=2,时，导弹位置,P,2,(x,2,y,2,);,敌艇位置,M,2,(902,120),；,P,(,x,(,t,),y,(,t,),O,x,M,(,v,e,t,H,),y,A,(0,H,),正东方向,正北方向,H,=120,（千米）,km/h,km/h,敌,导,(,x,k,(,t,),y,k,(,t,),function,num,y_j,L,T,=simu1(x,y,t,eps),k=0;,while,k1000 (simu1.m),p=90*k*,t-x,;,q=120-y;,d1=p/(p2+q2)0.5;,d2=q/(p2+q2)0.5;,x=x+450*t*d1;,y=y+450*t*d2;,if,(,abs(q,)0.1)%,终止条件,(),for,k=1:280,p=90*k*t-x(1);,q=120-x(2);,d1=p/(p2+q2)0.5;,d2=q/(p2+q2)0.5;,x(1)=x(1)+450*t*d1;,x(2)=x(2)+450*t*d2;,x1(1)=90*k*t;,x1(2)=120;,h1=line(,color,0,0.2,0.4,linewidth,2);,h2=line(,color,0,0.6,0.9,linewidth,3);,set(h1,xdata,x1(1),ydata,x1(2);,set(h2,xdata,x(1),ydata,x(2);,end,end,hold on,（,2,）如果当基地发射导弹的同时，敌艇立即由仪器发觉。假定敌艇为高速快艇，它即刻以,135,千米,/,小时的速度与导弹方向成垂直的方向逃逸，问导弹何时何地击中敌艇,？,不停地改变逃跑策略，运动轨迹如何？,修正敌艇速度：,90 135,千米,/,小时,思考：,二、离散系统,离散系统,是指系统状态只在有限的时间点或可数的时间点上发生变化的系统。并假设离散系统状态的变化是在一个时间点上瞬间完成。模型一般用流程图或网络来表示。其间可能涉及到随机事件等。,关键：模拟步骤、数据收集、模拟时钟,在我方某前沿防守地域，敌人以一个炮排（含两门火炮）为单位对我方进行干扰和破坏为躲避我方打击，敌方对其阵地进行了伪装并经常变换射击地点,经过长期观察发现，我方指挥所对敌方目标的指示有,50,是准确的，而我方火力单位，在指示正确时，有,1/3,的射击效果能毁伤敌人一门火炮，有,1/6,的射击效果能全部消灭敌人,现在希望能用某种方式把我方将要对敌人实施的,20,次打击结果显现出来，确定有效射击的比率及毁伤敌方火炮的平均值。,分析,：,这是一个概率问题，可以通过理论计算得到相应的概率和期望值,.,但这样只能给出作战行动的最终静态结果，而显示不出作战行动的动态过程,.,为了能显示我方,20,次射击的过程，现采用模拟的方式。,射击问题,需要模拟出以下两件事：,1.,问题分析,2,当指示正确时，我方火力单位的射击结果情况,1,观察所对目标的指示正确与否,模拟试验有两种结果，每一种结果出现的概率都是,1/2,因此，,可用投掷一枚硬币的方式予以确定,，当硬币出现正面时为指示正确，反之为不正确,模拟试验有三种结果：毁伤一门火炮的可能性为,1/3(,即,2/6),，毁伤两门的可能性为,1/6,，没能毁伤敌火炮的可能性为,1/2(,即,3/6),这时,可用投掷骰子的方法来确定,：,如果出现的是、三个点：则认为没能击中敌人；,如果出现的是、点：则认为毁伤敌人一门火炮；,若出现的是点：则认为毁伤敌人两门火炮,2.,符号假设,i,：要模拟的打击次数；,k,1,：没击中敌人火炮的射击总数；,k,2,：击中敌人一门火炮的射击总数；,k,3,：击中敌人两门火炮的射击总数,E,：有效射击比率；,E,1,：,20,次射击平均每次毁伤敌人的火炮数,3.,模拟框图,初始化,:i=0,k,1,=0,k,2,=0,k,3,=0,i=i+1,骰子,点数,?,k,1,=k,1,+1,k,2,=k,2,+1,k,3,=k,3,+1,k,1,=k,1,+1,i,20?,E=(k,2,+k,3,)/20 E,1,=0*k,1,/20+1*k,2,/20+2*k,3,/20,停止,硬币正面,?,Y,N,N,Y,1,，,2,，,3,4,，,5,6,4.,模拟结果,5.,理论计算,6.,结果比较,虽然模拟结果与理论计算不完全一致，但它却能更加真实地表达实际战斗动态过程,蒙特卡罗方法,(,Monte Carlo method,),蒙特卡罗方法,(,Monte Carlo method,),，或称计算机随机模拟方法，是一种基于“,随机数,”的计算方法。源于美国在第,二,次世界大战研制原子弹的“曼哈顿计划”,，,该计划的主持人之一数学家,冯诺伊曼,用驰名世界的,赌城,摩纳哥的Monte Carlo来命名这种方法，为它蒙上了一层神秘色彩。,蒙特卡罗,方法的,基本思想,很早以前就被人们所发现和利用。早在17世纪，人们就知道,用事件发生的“频率”来决定事件的“概率”,。1,9,世纪人们用,蒲丰,投针,的方法来,计算,圆周率,，上,世纪40年代电子计算机的出现，特别是近年来,高速电子计算机,的出现，使得用数学方法在计算机上大量、快速地模拟这样的试验成为可能。,蒲丰投针实验,近似计算圆周率,蒲丰投针实验,：,法国科学家蒲丰,(,Buffon,),在,1777,年提,出的蒲丰投针实验是早期几何概率一个,非常著名的例子。,蒲丰投针实验,的重要,性并非是为了求得比其它方法更精确的,值，而是它,开创了,使用,随机数处理确定性数学问题,的先河，是用,偶然性方法去解决确定性计算,的前导，,由此,可以,领略到从“概率土壤”上开出的一朵瑰丽的鲜花,蒙特卡罗,方法,(,MC),蒲丰投针实验,可归结为下面的数学问题,：平面上画有距离为,a,的一些平行线，向平面上任意投一根长为,l(l,a),的针，假设针落在任意位置的可能性相同，试求,针与平行线相交的概率,P(,从而求,),蒲丰投针实验,近似,计算圆周率,蒲丰投针实验,：,如右图所示，以,M,表示针落下后的中点，,以,x,表示,M,到最近一条,平行线的距离，以,表示针与此线的交角：,针落地的所有可能结果满足：,其样本空间视作矩形区域,面积是,:,针与平行线相交的条件：,它是样本空间,子集,A,，面积是：,syms,l,phi;int(l/2*sin(phi),phi,0,pi);%,ans,=,l,因此，针与平行线相交的概率为：,从而有：特别当 时,蒲丰投针实验,近似,计算圆周率,蒲丰投针实验,近似,计算圆周率,蒲丰投针实验,的计算机模拟：,format long;%,设置,15,位显示精度,a=1;,l=0.6;,%两平行线间的宽度,和,针长,figure;axis(0,pi,0,a/2);%,初始化绘图板,set(gca,nextplot,a,dd,);,%,初始化绘图方式为叠加,counter=0;,n=,2010,;,%,初始化计数器和设定,投针次数,x=unifrnd(0,a/2,1,n);,p,hi,=unifrnd(0,pi,1,n);,%,样本空间,for i=1:n,if x(i)l*sin(,p,hi,(i)/2,%满足此条件表示针与线的相交,plot(phi(i),x(i),r.),;,frame(i,)=,getframe,;%,描点并取帧,counter=counter+1;,%,统计针与线相交的次数,end,end,fren=counter/n;,pihat=2*l/(a*fren),%,用,频率,近似计算,蒲丰投针实验,近似,计算圆周率,蒲丰投针实验,的计算机模拟：,意大利数学家拉泽里尼得到了准确到,6,位小数的,值，不过他的实验因为,太准确,而受到了质疑,蒲丰投针实验计算圆周率,蒙特卡罗投点法是,蒲丰投针实验,的推广：,在一个边长为,a,的正方形内随机投点，该点落在此正方形的内切圆中的概率应为该内切圆与正方形的面积比值，即,n=10000;a=2;m=0;,for i=1:n,x=rand(1)*a;y=rand(1)*a;,if(x-a/2)2+(y-a/2)2=(a/2)2),m=m+1;,end,end,disp,(,近似为,:,num2str(4*,m/n,);,%,以下为向量化程序,x=rand(1,10000);y=rand(1,10000);,num2str(length(find(x-1).2+(y-1).2=1)*4/10000),x,y,o,(a/2,a/2,),一些,Monte Carlo,思想,的例子,例,1,(,相遇问题,),甲、乙两船在,24,小时内独立地随机到达码头,.,设两船到达码头时刻都服从,0,24,上的均匀分布，如果甲船到达码头后停留,2,小时,乙船到达码头后停留,1,小时,.,问两船相遇的概率有多大？（可用,几何概率，此处略,）,分析：如果知道甲、乙两船到达的时刻,x,和,y,，两船能相遇的条件就是：,两船到达的时刻,x,和,y,可以随机生成，生成一组到达时刻，可以确定是否能相遇。,如果重复很多次，统计相遇的比例就可近似为相遇的概率。,这也是,Monte Carlo,思想,.,两船到达码头时刻服从,0,24,上的均匀分布,，,甲船停留,2,小时,乙船停留,1,小时,相遇概率？,方法一,利用软件产生一组,0,，,24,上的随机数,x,i:,10.4 5.3 0.3 9.0 20.1 6.8 7.8,5.4,分别代表甲船到达码头时间；,再产生一组,0,，,24,上的随机数,y,i,:,6.0 3.4 8.1 17.5 0.8 13.4,14.0,分别,代表,乙船到达码头时间；,把（,x,i,，,y,i,）,(i,1,，,2,),当作一天中甲船乙船到达的时间，统计出能相遇的频数，计算出相遇的频率做为相遇的概率。,可以使用软件（,excel,、,C,等实现）,两船到达码头时刻服从,0,24,上的均匀分布,，,甲船停留,2,小时,乙船停留,1,小时,相遇概率？,方法二,clc,;N=1000;c=0;,%,模拟次数N，相遇次数c清零,for i=1:N,%重复N次到达时间,x=unifrnd(0,24,1,1);,%,甲船到达时间x(随机数,),y=unifrnd(0,24,1,1);,%,乙船到达时间y(随机数,),if(x,=y&y=x+2)|(y=x&x=y+1),c=c+1;,%如果能相遇，则计数器加,1,end,end,P=c/N,%,显示相遇的概率近似值,注意：每次运行的结果一般都不一样,%,计算机模拟程序,（报童诀窍的简化版）报童每天清晨从报社购进报纸零售，晚上退回没有卖掉的报纸,.,若每份报纸的购进价为,b=0.75,元，售出价为,a=1,元，退回价为,c=0.6,元,.,每天需求量,X,是离散型随机变量，其分布为,例,2,X,500,510,520,P,0.34,0.36,0.30,问：如果报童每天购进报纸为,n=510,份，每天的平均利润是多少？,方法二：如果我们知道每天的需求量，可直接计算利润。而每天需求量可以按分布生成（随机模拟思想）,方法一：概率方法（略）,售出价,a=1,，购进价,b=0.75,，退回价,c=0.6,，,购进数量,n=510,份,需求量,X,500,510,520,概率,P,0.34,0.36,0.30,按照需求量的分布规律，随机生成,N=20,个数据,:,510 520 500 510 520 500,500,500,500,510 500,500,500,520 510 520 510,510,代表,20,天的需求量，计算出报童在这,20,天的总利润和平均利润，用,平均利润,来近似报童的平均收入。,这也是,Monte Carlo,方法,.,问,题,如何按分布规律产生随机数据？,随机数据很多时，如何编程？,报童诀窍的简化版,需求量,X,500,510,520,概率,P,0.34,0.36,0.30,问题,1,：如何产生以下分布规律的随机数据？,注,:,rand(m,n,),可以生成,0,1,上均匀分布随机数,把,0,1,分成长度为,0.34,、,0.36,、,0.30,的三个区间,0,0.34,、,(0.34,0.70,、,(0.70,1,用,rand(1,1),产生,1,个,0,1,上均匀分布随机数，如该数在,0,0.34,、,(0.34,0.70,或,(0.70,1,内，相当于该天的需求量相应为,500,、,510,和,52,0,重复多次就可以若干天的需求量了,生成,N=20,天的需求量的,matlab,代码可以为,报童诀窍的简化版,生成,N=20,天的需求量的,matlab,代码（定义函数）,function,y=randfun1(N),y=zeros(1,N);,for,i=1:N,t=rand(1,1);,if,t=0.34 X=500;,elseif,tn),y=n*(a-b);,else,y=x*(a-b)-(n-x)*(b-c),end,模拟程序代码,N=1000;,x=randfun1(N);,y=0;,for i=1:N,y=y+fun2(x(i);,end,y/N,报童诀窍的简化版,a=1;b=0.75;c=0.6;n=510,;n=510;N=1000,;,y=0;,f,or,i=1:N,t=rand(1,1);,if,(t=0.34)x=500;,elseif,(tn)y=,y+n,*(a-b);,else,y=,y+x,*(a-,b)-(n-x,)*(,b-c,);,end,End,fprintf(1,平均利润=%.3f,y/N);,也可完整模拟程序：,根据随机数,t,，计算需求量,x,值,根据需求量,x,，计算利润并累加到,y,中。,显示平均利润,售出价,a,购进价,b,退回价,c,购进数n,总利润y,模拟天数,N,报童诀窍的简化版,粒子分离器的某关键参数（记作,y,）由,7,个零件的参数（记作,x,1,，,x,2,，,，,x,7,）决定，,均值,=,标定值 标准差,=,容差等级*标定值,3,关键参数,y,的目标值是,1.50,，当偏离为,0.1,但未达到,0.3,时，产品为次品，损失为,1000,元；当偏离达到,0.3,时，产品为废品，损失为,9000,元。由于工艺原因，,7,个零件参数可以看着是正态随机变量，在后面的标定值及容差等级情况下，求产品的平均损失？,容差等级,A=1%B=5%,C=15%,零件参数设计（,1997A,）,例,3,均值,=,标定值 标准差,=,容差等级*标定值,3,容差等级,A=1%B=5%,C=15%,零件损失,Q,是一个随机变量，平均损失就是期望,E(Q),由于,y,的表达式很复杂，要想计算,y,的分布和上述概率很困难，我们必须寻找较为有效的近似方法。,模拟：通过产生指定分布的随机数，来代表,7,个零件的参数值，计算,y,值，确定损失大小。多做几次可得均值,零件参数,%,编写计算,y,值的函数：,y=funY(x1,x2,x3,x4,x5,x6,x7)%,编写计算损失费函数,:Q=,funQ(y,),D1=0.1,0.3,0.1,0.1,1.5,16,0.75%,标定值,D2=5,5,5,15,15,5,5./100%,容差等级,a=D1b=D1.*D2./3,x1=normrnd(a(1),b(1),1,1);,x7=normrnd(a(7),b(7),1,1),y=funY(x1,x2,x3,x4,x5,x6,x7)Q=,funQ(y,）,产生一组零件参数，相当制作一个产品。重复,N,1000,次（即生产,N,个产品），求出损失费用的平均,.(,代码略,),按指定分布产生随机数，作为零件参数,计算该产品的关键参数,y,值和损失的大小,零件期望,a,零件标准差,b,零件参数,用蒙特卡洛法解非线性规划问题,基本假设,试验点的第,j,个分量,x,j,服从,a,j,b,j,内的均匀分布,符号假设,求解过程,先产生一个随机数作为初始试验点,以后则将上一个试验点的第,j,个分量随机产生，其它分量不变而产生一新的试验点这样，每产生一个新试验点只需一个新的随机数分量当,KMAXK,或,PMAXP,时停止迭代,框 图,初始化,:,给定,MAXK,MAXP;k=0,p=0,Q:,大整数,x,j,=,a,j,+R(b,j,-a,j,)j=1,2,n,j=0,j=j+1,p=p+1,PMAXP?,Y,N,x,j,=,a,j,+R(b,j,-a,j,),g,i,(X)0?,i=1,2n,Y,N,jMAXK?,Y,N,输出,X,Q,停止,Y,N,在,Matlab,软件包中编程，共需三个文件,:,randlp.m,mylp.m,lpconst.m,.,主程序为,randlp.m,.,%,mylp.m,function,z=,mylp(x,),%,目标函数,z=2*x(1)2+x(2)2-x(1)*x(2)-8*x(1)-3*x(2);,%,转化为求最小值问题,%,randlp.m,function,sol,r1,r2=,randlp(a,b,n,),%,随机模拟解非线性规划,debug=1;,a=0;,%,试验点下界,b=10;,%,试验点上界,n=1000;,%,试验点个数,r1=unifrnd(a,b,n,1);,%n,1,阶的,a,b,均匀分布随机数矩阵,r2=unifrnd(a,b,n,1);,sol=r1(1)r2(1);,z0=,inf,;,for,i=1:n,x1=r1(i);,x2=r2(i);,lpc,=lpconst(x1 x2),；,if,lpc,=1,z=mylp(x1 x2),；,if,zz0,z0=z,；,sol=x1 x2,；,end,end,end,To,Matlab(randlp,),返回,随机数的产生,1,、均匀随机数,(,均匀分布,U0,1),rand(),2,、产生其他随机数的方法,逆变换法、舍选法、近似抽样法等。,设概率分布函数,F(x,),是严格单调增的，,F,的反函数记作,F,-1,。先产生,UU(0,1),再取,X,F,-1,即为所求。,如指数分布，分布函数,F(x,)=1-exp(-),可得,3,、,(,非,),常见分布随机数如何产生？（离散）,经验分布函数法,:,设一组实际数据，将它们分组整理形成频率图或表格。,x0,X,f,3.5 5.5 8,0.3 0.5 0.2,3.5 3.5,5.5,）,5.5,8)8,构造经验分布函数,X,F,0 0.3 0.8 1,由均匀随机数,R0,1,，决定,X,的抽样值。,当,0R,0.3,，抽样值,x,i,（,0,3.5,；,x,i,=3.5,当,0.3R,0.8,，抽样值,x,i,（,3.5,5.5,；,x,i,=5.5,当,R0.8,，抽样值,x,i,（,5.5,8,；,x,i,=8,r=rand;,if,0r&r0.3,y=3.5;,elseif,0.3=r&r0.8,y=5.5;,else,y=8;,end,y%,产生一个随机数,Matlab,程序：,r=rand(10);,lisanrnd1.m,y=;,for,i=1:10,if,0,r(i)&r(i,)0.3,y(i,)=3.5;,elseif,0.3=,r(i)&r(i,)0,为常数，则称,X,服从,参数,为 的,指数分布,。,指数分布的期望值为,排队服务系统中顾客到达率为常数时的到达间隔、故障率为常数时零件的寿命都服从指数分布。,指数分布在排队论、可靠性分析中有广泛应用。,注意：,Matlab,中，产生,参数,为 的指数分布的命令为,exprnd,(),例,顾客到达某商店的间隔时间服从参数为,0.1,的指数分布,指数分布的均值为,1/0.1=10,。,指两个顾客到达商店的平均间隔时间是,10,个单位时间,.,即平均,10,个单位时间到达,1,个顾客,.,顾客到达的间隔时间可用,exprnd(10),模拟。,设离散型随机变量,X,的所有可能取值为,0,1,2,且取各个值的概率为,其中,0,为常数，则称,X,服从参数为 的,帕松分布,。,帕松分布在排队系统、产品检验、天文、物理等领域有广泛应用。,帕松分布的期望值为,如相继两个事件出现的间隔时间服从参数为 的指数分布，则在单位时间间隔内事件出现的次数服从参数为 的泊松分布即单位时间内该事件出现,k,次的概率为：,反之亦然。,指数分布与帕松分布的关系：,(1),指两个顾客到达商店的平均间隔时间是,10,个单位时间,.,即平均,10,个单位时间到达,1,个顾客,.,(2),指一个单位时间内平均到达,0.1,个顾客,例,(1),顾客到达某商店的间隔时间服从参数为,0.1,的指数分布,(2),该商店在单位时间内到达的顾客数服从参数为,0.1,的帕松分布,排队问题,机械故障等候修理 飞机跑道,日常生活中经常遇到的排队问题：,自选商场收款台 医院里病人等候就诊,输入情况：顾客到达时间和服务时间。,系统状态：排队等候的顾客数目（队长）,L(t,),服务员是否在工作或服务效率等；,简图,：,第二顾客接受,服务时间,s,2,x,5,0,x,1,x,2,x,3,x,4,y,1,y,2,y,3,y,4,y,5,D,2,关系,：,系统在什么条件下处于空闲状态？,(,y,i,x,i+1,),排队系统中，顾客到达时刻数据如何收集？对每个顾客的服务时间如何？,X,：,x,1,，,x,2,，,，,x,n,第一个顾客到达的时刻,第二个顾客到达的时刻,计算机遵循某种规则进行随机抽样。,S:s,1,x,2,x,n,排队服务系统的模拟,1,、模拟目的,2,、系统假设与输入,3,、系统的状态,4,、初始条件和终止条件,5,、模拟过程的运行,6,、系统的性能指标,7,、与分析模型的比较,8,、其它排队系统,1,、模拟目的,如果知道了顾客到达的时间和服务时间的统计规律（大量实际数据或概率分布），怎样通过模拟得到衡量系统的性能指标，并与分析模型的结果进行比较。,2,、系统假设与输入,首行作如下基本假设：,（,1,）顾客源是无穷的,（,2,）排队的长度没有限制,（,3,）到达系统的顾客按先后顺序进入服务,假设由概率分布产生了如下的数据：,k,1,2,3,4,5,6,7,8,9,i,k,0,1,3,4,1,3,1,8,2,s,k,2,3,1,3,4,1,3,4,1,其中,i,k,表示第,k,位顾客与第,k-1,位顾客到达的时间间隔，,s,k,是第,k,位顾客的服务时间，下面就以这组数据作为系统的输入。,3,、系统的状态,队长：排队等候的顾客数目。记作,L(t,).,服务员的状态用,S(t,),表示，当服务员工作时，令,S(t,)=1,，空闲时令,S(t,)=0.,事件：引起系统状态,L(t,),和,S(t,),改变的行为。,“,顾客到达,”,和,“,顾客离开,”,设第,k,位顾客到达和离开时刻分别为,a,k,和,d,k,，则有,4,、初始条件和终止条件,不妨假设,t=0,是第,1,位顾客到达的时刻，此,时服务员处于空闲状态。模拟的终止时刻是,T,（给定值），或模拟直到第,n,个（给定值）顾,客进入服务的时间。,5,、模拟过程,进行模拟时，可以事先产生一批数据,i,k,s,k,再计算,a,k,d,k,，然后让时钟,t,按,a,k,d,k,从小到大的顺序跳动；也可在时钟,t,跳到某一事件发生时，才产生一个所需要的,i,或,s,。这里采用后一种方法。,设当前时钟为,t,，队长,L(t,),记作,WL,，服务员状态,S(t,),记作,SS,，,t,以后下一个,“,顾客到达,”,事件的发生时刻记作,AT,t,以后下一个,“,顾客离开,”,事件的发生时刻记作,DT.,AT0,WL=WL+1,令,SS,1,产生,s,置,DT=,t+s,产生,i,置,AT,t+i,t=T,令,SS,0,WL,WL,1,产生,s,置,DT,t+s,停止,置,DT,999,是,否,是,否,是,否,否,是,st,=6;it=5;WL=0;SS=0;AT=0;DT=99999;T=10;t=0;ttt=0;www=0;sss=0;,while(t=T),if AT0,WL=WL-1;s=exprnd(1/st,1,1);DT=,t+s,;,else,SS=0;DT=99999;,end,end,ttt,=,ttt,t,;www=,www,WL,;,sss,=,sss,SS,;,end,ttt(2:length(ttt),www(2:length(ttt),sss(2:length(ttt),6,、系统的性能指标,平均队长,平均待时间,服务利用率,7,、与分析模型的比较,对于一般的排队服务系统，很难建立分析模型并求出结果，即使是单服务员系统，也只有在顾客到达间隔和服务时间服从某种特殊的概率分布时，才可以在稳态情况下由分析模型得到上述三个指标的解析形式的结果。,分析模型基本假设：,1,、到达间隔服从参数,的指数分布；,2,、服务时间服从参数为,的指数分布；,3,、排队规则是,“,先到先服务,”,，顾客数目和排队长度都无限制。,排队论中将这样的单服务员系统记作,M/M/1,。,服务强度：,稳态（,1,）,8,、其它的排队服务系统,有两个（或多个）服务员,串联服务系统，每个顾客接受多次服务,从排队顾客中选服务时间最短的最先服务,队的长度有限制,带有优化目标的排队服务系统,库存太多造成浪费或资金积压，库存太少不能满足需求也造成损失，需要进行决策：何时进货，进多少，使得平均费用最少，而收益最大。,库存系统,方案甲：按前一天的销售量作为当天的生产量；,方案乙：按前二天的平均销售量作为当天的生产量；,假定市场对该产品的每天,需要量,是一个随机变量，但从以往的统计分析得知它服从正态分布,N,(135,22.4,2,）,.,例,1,某企业生产易变质的产品。（如蛋糕）当天生产的产品必须当天售出，否则就会变质。该产品单位成本为,2.5,元，单位产品售价为,5,元。企业为避免存货过多而造成损失，拟从以下两种货存方案中选出一个较优的方案,.(,如何决定当天的生产数量？,),模拟的基本思路：,需要获得市场对该产品的需要量的样本,值；,(,产生抽样值,),寻找货存量与需要量之间的关系；,按照两种不同方案（货存量,需要量）计算出经,T,天后企业的利润值,(,累计,),；,比较大小，从中选出一个较优的方案；,方案甲：货存量,Y1,分析两种状态下,X,Y1,和,XY1,的利润值,L1,置初始状态,获取产品的需求量,X,X=,normrnd,（,135,，,22.5,）,N=1,方案乙：货存量,Y2,分析两种状态下,X,Y2,和,XY2,的利润值,L2,累计利润值,TL1,累计利润值,TL2,N=N+1,判断：,NT(T,天,),比较计算：,maxTL1,，,TL2,输出最大利润值及方案,MATLAB,编制,的,程序,：,kucun.m。,运行程序：,kutest.m,初始状态：,T=100;S1=136;S21=126;S22=148;,TL1,TL2=kucun(T,S1,S21,S22),观察计算结果,如果多数情况有：,TL1TL2,则,认为方案乙较优.,思考：,1,）一般模拟结果波动性较大，如何减少这种波动？,2,）修改上述程序。,可靠性问题,一设备上有三个相同的轴承，每个轴承正常工作寿命为随机变量，概率分布如下：,寿命,/h,1000,1100,1200,1300,1400,1500,1600,1700,1800,1900,概率,0.10,0.13,0.25,0.13,0.09,0.12,0.02,0.06,0.05,0.05,有轴承损坏,设备停止工作,检修工准备开始更换部件，称为一个延迟时间，它也是随机变量，分布如下：,延迟时间,/min,5,10,15,概 率,0.6,0.3,0.1,主要费用：,1,、设备停工损失费：,5,元,/,分钟；,2,、检修工人的工时费：,12,元,/,小时；,3,、轴承的成本费：,16,元,/,个,更换轴承所需要的时间：,一个 两个 三个,20 30 40,（,min,）,问题：,现在有两种方案：,方案一：损坏一个轴承只更换一个轴承；,方案二：一旦有轴承损坏就全部更换；,试通过计算机模拟对以上两种方案做出评价。,随机数怎样产生？,模拟时选用时间步长法还是,事件步长法？,关于随机数的产生见,Lrnd.m,（零件寿命）,Yrnd.m,（延迟时间）,方案一的数学模型：,kekao1.m,目标函数,min,c,=,U,i,/T,其中,损失费,工时费,成本费,t,i,表示延迟时间,设事件发生时刻,T,i,，它是由轴承工作寿命,L,i,、延迟时间和修理时间迭加产生。总的运行时间为,T,（,10000,）。,例：（手工模拟）,事件类型,发生时刻,延迟时间,A,1400 h,5 min,B,1500 h,15 min,C,1500 h,15 min,方案一的初始事件表,t=0t=1400,；,cost=(5+20)5+12 1/3+16=145,；,产生下一个,A,事件发生时刻,1400+25/60+1000,（,L,i,）,=2400,小时,25,分钟,损失费,工时费,成本费,事件类型,发生时刻,延迟时间,B,1500 h,15 min,C,1500 h,15 min,A,2400 h25min,5 min,第一次刷新后的事件表,t=1400 t=1500,；,cost=145+(15+30)5+12 1/2+16 2,=408,；,产生下一个,B,事件发生时刻,1500+45/60+1200,（,L,i,）,=2700,小时,45,分钟,c1=5;c2=12;c3=16;(kekao1.m),g1=20/60;g2=30/60;g3=40/60;,Lrnd,;,%Random function,Lrnd,gets a random,life(three,),yrnd,;,%Random function,yrnd,gets a random,delate(one,),lm=,sort(l,);,U=;t=;,for,j=1:50,if,lm(1)lm(2),U(j,)=c1*(y+g1*60)+c2*g1+c3;,elseif,lm(1)=lm(2)lm(3),U(j,)=c1*(y+g2*60)+c2*g2+2*c3;,else,U(j,)=c1*(y+g3*60)+c2*g3+3*c3;,end,Lrnd1;,%(one),t(j,)=lm(1)+y/60+l1;,L=lm(2),lm(3),t(j);,lm=,sort(L,);,yrnd,;,end,U;t,;,zU,=,sum(U,)C=zU/t(50),计算结果：,ZU=7850,（元），,C=0.3354,（元,/,小时）。,方案二的数学模型：,目标函数,min,c,=,U,i,/T,其中,方案二的情况比较单一。,Kekao3.m,模拟框图如下：,数据初始化,T100000,产生下一个事件发生时刻,根据不同方案确定更换策略,产生延迟时间,计算当前更换费用,累积费用,产生新轴承的寿命时间,确定下一个事件发生的时刻,系统时刻跳转到下一个事件发生时刻,输出总费用,停止,否,是,c1=5;c2=12;c3=16;(kekao3.m),g1=20/60;g2=30/60;g3=40/60;,Lrnd,;,%(three),yrn