数模竞赛培训之计算机模拟

数学建模培训计算机模拟主讲：何仁斌概述

计算机模拟是用计算机对系统的结构和行为进行动态演示，以评价或预测系统的行为效果，为决策者提供信息的一种方法。它是解决较复杂的实际问题的一条有效途径。计算机模拟也可以说是用计算机程序直接建立真实系统的模型，并通过计算了解系统随时间变化的行为或特征。应用领域：航空、机电、冶金、社会经济、交通运输、生态系统等。计算机模拟分为连续系统模拟和离散系统模拟。

状态随时间连续变化的系统称为连续系统。通常该系统的模型一般可以用微分方程的形式表达，通过一些物理机理推导出来。模拟结果往往是近似的。例如，追逐问题，浓度问题。一、连续系统(x1(t),…,xn(t))1、追逐问题A→BB→CC→DD→A状态变量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:(x1(t),y1(t))→(x1(t+△t),y1(t+△t))B:(x2(t),y2(t))→(x2(t+△t),y2(t+△t))(x1(t+△t),y1(t+△t))（由几何原理）≈(x1(t)+v△tcos(θ),y1(t)+v△tsin(θ)θx1,y1x2,y2在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=[0001010

10

100];x(9)=x(1);x(10)=x(2);holdonaxis('equal');axis([010010]);while(d>0.1)

fori=1:2:7d=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);endholdonABCD

某军的一导弹基地发现正北方向120km处海面上有敌艇一艘以90km/h的速度向正东方向行驶。该基地立即发射导弹跟踪追击敌艇，导弹速度为450km/h，自动导航系统使导弹在任一时刻都能对准敌艇。试问导弹在何时何处击中敌艇？2、导弹跟踪问题P(x(t),y(t))OxM(vet,H)yA(0,H)正东方向正北方向H=120（千米）km/hkm/h敌导(xk(t),yk(t))simu1.mstest1.mhuatu.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如何决定导弹位置

P2(x2,y2)?计算机模拟1)当t=0时，导弹位置O(0,0);敌艇位置A(0,120)；2)当t=τ时，导弹位置P1(x1,y1);敌艇位置M1(90τ,120)；3)当t=2τ时，导弹位置P2(x2,y2);敌艇位置M2(90×2τ,120)；function[num,y_j,L,T]=simu1(x,y,t,eps)k=0;whilek<1000(simu1.m)p=90*k*t-x;q=120-y;d1=p/(p^2+q^2)^0.5;d2=q/(p^2+q^2)^0.5;x=x+450*t*d1;y=y+450*t*d2;

if(abs(q)<eps)

break;

end

k=k+1;endnum=k;L=x;y_j=y;T=L/90;(stest1.m)t=0.1;x=0;y=0;[num,y_j,L,T]=simu1(x,y,t,eps)计算结果：Num=2;y_j=131.2155;L=22.675（km）;T=0.2519（h）导弹的位置(x(t),y(t))能否用图形描述它们的追逐运动过程？敌艇的位置(水平方向(90*k*t,120)x=[0,0];(huatu.m)axis('equal');axis([0,140,0,140]);gridonholdont=0.001;%步长while(abs(x(2)-120)>0.1)%终止条件(<)

fork=1:280p=90*k*t-x(1);q=120-x(2);d1=p/(p^2+q^2)^0.5;d2=q/(p^2+q^2)^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));

endendholdon（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：要模拟的打击次数；k1：没击中敌人火炮的射击总数；k2：击中敌人一门火炮的射击总数；k3：击中敌人两门火炮的射击总数．E：有效射击比率；E1：20次射击平均每次毁伤敌人的火炮数．3.模拟框图初始化:i=0,k1=0,k2=0,k3=0i=i+1骰子点数?k1=k1+1k2=k2+1k3=k3+1k1=k1+1i＜20?E=(k2+k3)/20E1=0*k1/20+1*k2/20+2*k3/20停止硬币正面?YNNY1，2，34，564.模拟结果5.理论计算6.结果比较

虽然模拟结果与理论计算不完全一致，但它却能更加真实地表达实际战斗动态过程．蒙特卡罗方法(MonteCarlomethod)

蒙特卡罗方法(MonteCarlomethod)，或称计算机随机模拟方法，是一种基于“随机数”的计算方法。源于美国在第二次世界大战研制原子弹的“曼哈顿计划”，该计划的主持人之一数学家冯·诺伊曼用驰名世界的赌城—摩纳哥的MonteCarlo—来命名这种方法，为它蒙上了一层神秘色彩。蒙特卡罗方法的基本思想很早以前就被人们所发现和利用。早在17世纪，人们就知道用事件发生的“频率”来决定事件的“概率”。19世纪人们用蒲丰投针的方法来计算圆周率π，上世纪40年代电子计算机的出现，特别是近年来高速电子计算机的出现，使得用数学方法在计算机上大量、快速地模拟这样的试验成为可能。蒲丰投针实验近似计算圆周率π蒲丰投针实验：法国科学家蒲丰(Buffon)在1777年提出的蒲丰投针实验是早期几何概率一个非常著名的例子。蒲丰投针实验的重要性并非是为了求得比其它方法更精确的π值，而是它开创了使用随机数处理确定性数学问题的先河，是用偶然性方法去解决确定性计算的前导，由此可以领略到从“概率土壤”上开出的一朵瑰丽的鲜花——蒙特卡罗方法(MC)蒲丰投针实验可归结为下面的数学问题：平面上画有距离为a的一些平行线，向平面上任意投一根长为l(l<a)的针，假设针落在任意位置的可能性相同，试求针与平行线相交的概率P(从而求π)蒲丰投针实验近似计算圆周率π蒲丰投针实验：如右图所示，以M表示针落下后的中点，以x表示M到最近一条平行线的距离，以φ表示针与此线的交角：针落地的所有可能结果满足：其样本空间视作矩形区域Ω,面积是:针与平行线相交的条件：它是样本空间Ω子集A，面积是：syms

lphi;int('l/2*sin(phi)',phi,0,pi);%ans=l因此，针与平行线相交的概率为：从而有：特别当时蒲丰投针实验近似计算圆周率π蒲丰投针实验近似计算圆周率π蒲丰投针实验的计算机模拟：formatlong;%设置15位显示精度a=1;

l=0.6;

%两平行线间的宽度和针长figure;axis([0,pi,0,a/2]);%初始化绘图板set(gca,'nextplot','add');%初始化绘图方式为叠加counter=0;

n=2010;

%初始化计数器和设定投针次数x=unifrnd(0,a/2,1,n);phi=unifrnd(0,pi,1,n);

%样本空间Ωfori=1:nifx(i)<l*sin(phi(i))/2

%满足此条件表示针与线的相交

plot(phi(i),x(i),'r.');frame(i)=getframe;%描点并取帧

counter=counter+1;%统计针与线相交的次数endendfren=counter/n;

pihat=2*l/(a*fren)

%用频率近似计算π%movie(frame,1)%播放帧动画1次蒲丰投针实验近似计算圆周率π蒲丰投针实验的计算机模拟：意大利数学家拉泽里尼得到了准确到6位小数的π值，不过他的实验因为太准确而受到了质疑蒲丰投针实验计算圆周率π蒙特卡罗投点法是蒲丰投针实验的推广：在一个边长为a的正方形内随机投点，该点落在此正方形的内切圆中的概率应为该内切圆与正方形的面积比值，即n=10000;a=2;m=0;fori=1:nx=rand(1)*a;y=rand(1)*a;if((x-a/2)^2+(y-a/2)^2<=(a/2)^2)m=m+1;endenddisp([‘π近似为:',num2str(4*m/n)]);%以下为向量化程序x=rand(1,n);y=rand(1,n);num2str(length(find((x-1).^2+(y-1).^2<=1))*4/n)xyo(a/2,a/2)一些MonteCarlo思想的例子[例1](相遇问题)甲、乙两船在24小时内独立地随机到达码头.设两船到达码头时刻都服从[0,24]上的均匀分布，如果甲船到达码头后停留2小时,乙船到达码头后停留1小时.问两船相遇的概率有多大？（可用几何概率，此处略）分析：如果知道甲、乙两船到达的时刻x和y，两船能相遇的条件就是：两船到达的时刻x和y可以随机生成，生成一组到达时刻，可以确定是否能相遇。如果重复很多次，统计相遇的比例就可近似为相遇的概率。这也是MonteCarlo思想.两船到达码头时刻服从[0,24]上的均匀分布，甲船停留2小时,乙船停留1小时,相遇概率？方法一利用软件产生一组[0，24]上的随机数xi:

10.45.30.39.020.16.87.8……5.4

分别代表甲船到达码头时间；再产生一组[0，24]上的随机数yi:

6.03.48.117.50.813.4……14.0

分别代表乙船到达码头时间；把（xi，yi）(i＝1，2…)当作一天中甲船乙船到达的时间，统计出能相遇的频数，计算出相遇的频率做为相遇的概率。可以使用软件（excel、C等实现）两船到达码头时刻服从[0,24]上的均匀分布，甲船停留2小时,乙船停留1小时,相遇概率？方法二clc;N=1000;c=0;

%模拟次数N，相遇次数c清零fori=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

endendP=c/N

%显示相遇的概率近似值注意：每次运行的结果一般都不一样%计算机模拟程序（报童诀窍的简化版）报童每天清晨从报社购进报纸零售，晚上退回没有卖掉的报纸.若每份报纸的购进价为b=0.75元，售出价为a=1元，退回价为c=0.6元.每天需求量X是离散型随机变量，其分布为[例2]X500510520P0.340.360.30问：如果报童每天购进报纸为n=510份，每天的平均利润是多少？方法二：如果我们知道每天的需求量，可直接计算利润。而每天需求量可以按分布生成（随机模拟思想）方法一：概率方法（略）售出价a=1，购进价b=0.75，退回价c=0.6，购进数量n=510份需求量X500510520概率P0.340.360.30按照需求量的分布规律，随机生成N=20个数据:510520500510520500500

500

500510500500

500520510520510510代表20天的需求量，计算出报童在这20天的总利润和平均利润，用平均利润来近似报童的平均收入。这也是MonteCarlo方法.问题如何按分布规律产生随机数据？随机数据很多时，如何编程？报童诀窍的简化版需求量X500510520概率P0.340.360.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和520重复多次就可以若干天的需求量了生成N=20天的需求量的matlab代码可以为报童诀窍的简化版生成N=20天的需求量的matlab代码（定义函数）functiony=randfun1(N)y=zeros(1,N);fori=1:Nt=rand(1,1);

ift<=0.34

X=500;

elseift<=0.70

X=510;

else

X=520;

end

y(i)=X;end根据随机数t的范围，确定需求量X值，并保存到数组的相应位置中（关键部分）产生1行N列的全0矩阵，目的分配好矩阵大小，可以省略文件名为randfun1.m报童诀窍的简化版问题2：如何从需求量计算利润？售出价a=1，购进价b=0.75，退回价c=0.6，购进数n=510函数文件名fun2.m

functiony=fun2(x)a=1;%售出价

b=0.75;%购进价

c=0.6;%退回价

n=510;%购进数量if(x>n)

y=n*(a-b);

else

y=x*(a-b)-(n-x)*(b-c)

end模拟程序代码N=1000;x=randfun1(N);y=0;fori=1:Ny=y+fun2(x(i));endy/N报童诀窍的简化版a=1;b=0.75;c=0.6;n=510;

n=510;N=1000;y=0;fori=1:Nt=rand(1,1);

if(t<=0.34)x=500;

elseif(t<=0.70)x=510;

elsex=520;

end

if(x>n)y=y+n*(a-b);

elsey=y+x*(a-b)-(n-x)*(b-c);

endEndfprintf(1,'平均利润=%.3f',y/N);也可完整模拟程序：根据随机数t，计算需求量x值根据需求量x，计算利润并累加到y中。显示平均利润售出价a

购进价b

退回价c

购进数n

总利润y

模拟天数N报童诀窍的简化版粒子分离器的某关键参数（记作y）由7个零件的参数（记作x1，x2，…，x7）决定，均值=标定值标准差=容差等级*标定值÷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=D1

b=D1.*D2./3x1=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个分量xj服从[aj,bj]内的均匀分布．符号假设求解过程先产生一个随机数作为初始试验点,以后则将上一个试验点的第j个分量随机产生，其它分量不变而产生一新的试验点．这样，每产生一个新试验点只需一个新的随机数分量．当K>MAXK或P>MAXP时停止迭代．框图初始化:给定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,停止YN

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

mylp.m,

lpconst.m.主程序为randlp.m.%mylp.mfunctionz=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;fori=1:nx1=r1(i);x2=r2(i);

lpc=lpconst([x1x2])；

if

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

ifz<z0z0=z；

sol=[x1x2]；

end

endendToMatlab(randlp)返回随机数的产生1、均匀随机数(均匀分布U[0,1])rand()2、产生其他随机数的方法逆变换法、舍选法、近似抽样法等。

设概率分布函数F(x)是严格单调增的，F的反函数记作F-1。先产生U~U(0,1),再取X＝F-1即为所求。如指数分布，分布函数F(x)=1-exp(-λ),可得

3、(非)常见分布随机数如何产生？（离散）经验分布函数法:①设一组实际数据，将它们分组整理形成频率图或表格。x≥0Xf3.55.580.30.50.2

＜3.5[3.5,5.5）[5.5,8)≥8构造经验分布函数XF00.30.81②由均匀随机数R∈[0,1]，决定X的抽样值。当0<R＜0.3，抽样值xi∈（0,3.5]；xi=3.5

当0.3≤R＜0.8，抽样值xi∈（3.5,5.5]；xi=5.5

当R≥0.8，抽样值xi∈（5.5,8]；xi=8r=rand;if0<r&r<0.3y=3.5;elseif

0.3<=r&r<0.8y=5.5;else

y=8;endy%产生一个随机数Matlab程序：r=rand(10);lisanrnd1.my=[];fori=1:10if0<r(i)&r(i)<0.3

y(i)=3.5;

elseif

0.3<=r(i)&r(i)<0.8

y(i)=5.5;else

y(i)=8;endendy%产生随机变量数组。

计算结果y=Columns1through75.55.53.55.53.55.55.5Columns8through108.05.55.5MATLAB的实现常见分布的随机变量获取抽样值1、U[0,1]R=rand(m,n)

2、U[a,b]R=unifrnd(a,b,m,n)3、E()

R=exprnd(1/

,m,n)4、N(0,1)R=randn(m,n)5、N(mu,sigma)R=normrnd(mu,sigma,m,n)6、B(N,p)R=binornd(N,p,m,n)7、P(

)R=poissrnd(

,m,n)以上语句都是产生m

n

的随机矩阵.(hist(x))产生模拟随机数的计算机命令2．产生m

n阶［０，１］均匀分布的随机数矩阵：rand(m,n)

产生一个［０，１］均匀分布的随机数：rand1．产生m

n阶［a，b］均匀分布U（a，b）的随机数矩阵：

unifrnd(a,b,m,n)

产生一个［a，b］均匀分布的随机数：unifrnd(a,b)

当只知道一个随机变量取值在（a，b）内，但不知道（也没理由假设）它在何处取值的概率大，在何处取值的概率小，就只好用U（a，b）来模拟它。当研究对象视为大量相互独立的随机变量之和，且其中每一种变量对总和的影响都很小时，可以认为该对象服从正态分布。机械加工得到的零件尺寸的偏差、射击命中点与目标的偏差、各种测量误差、人的身高、体重等，都可近似看成服从正态分布。若连续型随机变量X的概率密度函数为其中>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)

服务员是否在工作或服务效率等；简图：第二顾客接受服务时间s2x50x1x2x3x4y1y2y3y4y5D2关系：系统在什么条件下处于空闲状态？(yi<xi+1)

排队系统中，顾客到达时刻数据如何收集？对每个顾客的服务时间如何？X：x1，x2，…，xn第一个顾客到达的时刻第二个顾客到达的时刻计算机遵循某种规则进行随机抽样。S:s1,x2,…,xn排队服务系统的模拟1、模拟目的2、系统假设与输入3、系统的状态4、初始条件和终止条件5、模拟过程的运行6、系统的性能指标7、与分析模型的比较8、其它排队系统1、模拟目的如果知道了顾客到达的时间和服务时间的统计规律（大量实际数据或概率分布），怎样通过模拟得到衡量系统的性能指标，并与分析模型的结果进行比较。2、系统假设与输入首行作如下基本假设：（1）顾客源是无穷的（2）排队的长度没有限制（3）到达系统的顾客按先后顺序进入服务假设由概率分布产生了如下的数据：k123456789…ik013413182…sk231341341…其中ik表示第k位顾客与第k-1位顾客到达的时间间隔，sk是第k位顾客的服务时间，下面就以这组数据作为系统的输入。3、系统的状态队长：排队等候的顾客数目。记作L(t).服务员的状态用S(t)表示，当服务员工作时，令S(t)=1，空闲时令S(t)=0.事件：引起系统状态L(t)和S(t)改变的行为。

“顾客到达”和“顾客离开”设第k位顾客到达和离开时刻分别为ak和dk，则有

4、初始条件和终止条件不妨假设t=0是第1位顾客到达的时刻，此时服务员处于空闲状态。模拟的终止时刻是T（给定值），或模拟直到第n个（给定值）顾客进入服务的时间。5、模拟过程

进行模拟时，可以事先产生一批数据ik,sk,再计算ak,dk，然后让时钟t按ak,dk从小到大的顺序跳动；也可在时钟t跳到某一事件发生时，才产生一个所需要的i或s。这里采用后一种方法。设当前时钟为t，队长L(t)记作WL，服务员状态S(t)记作SS，t以后下一个“顾客到达”事件的发生时刻记作AT,t以后下一个“顾客离开”事件的发生时刻记作DT.AT<DTWL=0,SS=0,AT=0,DT=999置t=ATSS=0置t=DTWL>0WL=WL+1令SS＝1产生s置DT=t+s产生i置AT＝t+it>=T令SS＝0WL＝WL－1产生s置DT＝t+s停止置DT＝999是否是否是否否是6、系统的性能指标平均队长平均待时间服务利用率7、与分析模型的比较对于一般的排队服务系统，很难建立分析模型并求出结果，即使是单服务员系统，也只有在顾客到达间隔和服务时间服从某种特殊的概率分布时，才可以在稳态情况下由分析模型得到上述三个指标的解析形式的结果。分析模型基本假设：1、到达间隔服从参数λ的指数分布；2、服务时间服从参数为μ的指数分布；3、排队规则是“先到先服务”，顾客数目和排队长度都无限制。

排队论中将这样的单服务员系统记作M/M/1。

服务强度：稳态（ρ<1）8、其它的排队服务系统有两个（或多个）服务员串联服务系统，每个顾客接受多次服务从排队顾客中选服务时间最短的最先服务队的长度有限制带有优化目标的排队服务系统

库存太多造成浪费或资金积压，库存太少不能满足需求也造成损失，需要进行决策：何时进货，进多少，使得平均费用最少，而收益最大。库存系统方案甲：按前一天的销售量作为当天的生产量；方案乙：按前二天的平均销售量作为当天的生产量；假定市场对该产品的每天需要量是一个随机变量，但从以往的统计分析得知它服从正态分布N(135,22.42）.例1某企业生产易变质的产品。（如蛋糕）当天生产的产品必须当天售出，否则就会变质。该产品单位成本为2.5元，单位产品售价为5元。企业为避免存货过多而造成损失，拟从以下两种货存方案中选出一个较优的方案.(如何决定当天的生产数量？)模拟的基本思路：需要获得市场对该产品的需要量的样本值；(产生抽样值)②寻找货存量与需要量之间的关系；③按照两种不同方案（货存量≥{≤}需要量）计算出经T天后企业的利润值(累计)；④比较大小，从中选出一个较优的方案；方案甲：货存量Y1分析两种状态下X＞Y1和X≤Y1的利润值L1置初始状态获取产品的需求量XX=normrnd（135，22.5）N=1方案乙：货存量Y2分析两种状态下X＞Y2和X≤Y2的利润值L2累计利润值TL1累计利润值TL2N=N+1判断：N≤T(T天)比较计算：max{TL1，TL2}输出最大利润值及方案MATLAB编制的程序：

kucun.m。运行程序：kutest.m初始状态：T=100;S1=136;S21=126;S22=148;[TL1,TL2]=kucun(T,S1,S21,S22)观察计算结果如果多数情况有：TL1<TL2,则认为方案乙较优.思考：1）一般模拟结果波动性较大，如何减少这种波动？2）修改上述程序。可靠性问题

一设备上有三个相同的轴承，每个轴承正常工作寿命为随机变量，概率分布如下：寿命/h1000110012001300140015001600170018001900概率0.100.130.250.130.090.120.020.060.050.05有轴承损坏→设备停止工作→检修工准备开始更换部件，称为一个延迟时间，它也是随机变量，分布如下：延迟时间/min51015概率0.60.30.1主要费用：1、设备停工损失费：5元/分钟；2、检修工人的工时费：12元/小时；3、轴承的成本费：16元/个更换轴承所需要的时间：一个两个三个

203040（min）问题：现在有两种方案：方案一：损坏一个轴承只更换一个轴承；方案二：一旦有轴承损坏就全部更换；试通过计算机模拟对以上两种方案做出评价。①随机数怎样产生？②模拟时选用时间步长法还是事件步长法？关于随机数的产生见Lrnd.m（零件寿命）

Yrnd.m（延迟时间）

方案一的数学模型：kekao1.m目标函数minc=∑Ui

/T其中损失费工时费成本费ti

表示延迟时间

设事件发生时刻Ti，它是由轴承工作寿命Li、延迟时间和修理时间迭加产生。总的运行时间为T（10000）。例：（手工模拟）事件类型发生时刻延迟时间A1400h5minB1500h15minC1500h15min方案一的初始事件表t=0→t=1400；

cost=(5+20)×5+12×1/3+16=145；产生下一个A事件发生时刻1400+25/60+1000（Li）=2400小时25分钟损失费工时费成本费事件类型发生时刻延迟时间B1500h15minC1500h15minA2400h25min5min第一次刷新后的事件表t=1400→t=1500；

cost=145+(15+30)×5+12×1/2+16×2=408；③产生下一个B事件发生时刻1500+45/60+1200（Li）=2700小时45分钟c1=5;c2=12;c3=16;(kekao1