数学建模之计算机仿真-PPT.ppt

数学建模之China Undergraduate Mathematical Contest in Modeling 一个问题l我们做一个实验:把一个硬币掷一万次,统计两个面出现的次数。这样做很简单但却需要大师时间，有没有一咱较快的办法把这个实验完成呢？掷硬币仿真流程图概述l 计算机科学技术的迅猛发展，给许多学科带来了巨大的影响计算机不但使问题的求解变得更加方便、快捷和精确，而且使得解决实际问题的领域更加广泛计算机适合于解决那些规模大、难以解析化以及不确定的数学模型例如对于一些带随机因素的复杂系统，用分析方法建模常常需要作许多简化假设，与面临的实际问题可能相差甚远，以致解答根本无法应用，这时仿真几乎成为人们的唯一选择在历届的美国和中国大学生的数学建模竞赛（MCM）中，学生们经常用到计算机仿真方法去求解、检验等计算机仿真(computer simulation)是建模过程中较为重要的一类方法 计算机仿真可以解决以下5类问题:l(1)难以用数学公式表示的系统,或者没有建立和求解的有效方法.l(2)虽然可以用解析的方法解决问题,但数学的分析与计算过于复杂,这时计算机仿真可能提供简单可行的求解方法.l(3)希望能在较短的时间内观察到系统发展的全过程,以估计某些参数对系统行为的影响.l(4)难以在实际环境中进行试验和观察时,计算机仿真是唯一可行的方法,如太空飞行的研究.l(5)需要对系统或过程进行长期运行比较,从大量方案中寻找最优方案.计算机仿真案例1模型建立：由于本题要求使从搅拌中心到各个工地运输混凝土的总的吨公里数最少，所以，该问题的目标函数是 求解方法：1、高数中的方法2、数值计算方法3、计算机仿真：离散化，遍历！计算机仿真案例2n 例2（赶火车过程仿真）一列火车从A站经过B站开往C站，某人每天赶往B站乘这趟火车。已知火车从A站到B站的运行时间是均值为30min、标准差为2min的正态随机变量。火车大约在下午1点离开站。火车离开时刻的频率分布和这个人到达站时刻的频率分布如下表所示。问他能赶上火车的概率有多大？出发时刻1:00 1:05 1:10到达时刻1:28 1:30 1:32 1:34频率0.7 0.2 0.1频率0.3 0.4 0.2 0.1u仿真过程：u1、生成火车的发车时间、运行时间，从而达得到其到达B站的时间。u2、生成此人达到B站的时间。u3、如果此人到达站的时间早于火车到达时间，则算赶上火车一次。u4、将上述过程重复一万次，统计赶上火车的频率作为所求概率。分析：这个问题用概率论的方法求解十分困难，它涉及此人到达时刻、火车离开A站的时刻、火车运行时间几个随机变量。我们可以用计算机仿真的方法来解决。概述l 计算机仿真在计算机中运行实现,不怕破坏,易修改,可重用,安全经济,不受外界条件和场地空间的限制.l 仿真分为静态仿真(static simulation)和动态仿真(dynamic simulation).数值积分中的蒙特卡洛方法是典型的静态仿真动态仿真又分为连续系统仿真和离散系统仿真连续系统是指状态变量随着时间连续变化的系统，例如传染病的检测与预报系统.离散系统是指系统状态变量只在有限的时间点或可数的时间点上发生变化的系统,例如排队系统.概述l 仿真系统，必须设置一个仿真时钟(simulate clock)，它能将时间从一个时刻向另一个时刻进行推进，并且能随时反映系统时间的当前值其中，模拟时间推进方式有两种:时间步长法(均匀间隔时间推进法,连续系统常用)和事件步长法(下次事件推进法,离散系统常用)主要内容l 一:准备知识:随机数的产生l 二:随机变量的模拟l 三:连续系统的模拟-时间步长法l 四:离散系统的模拟-事件步长法l 五：蒙特卡洛方法一:准备知识:随机数的产生l 由于仿真研究的实际系统要受到多种随机因素的作用和影响,在仿真过程中必须处理大量的随机因素.要解决此问题的前提是确定随机变量的类型和选择合适的随机数产生的方法.l 对随机现象进行模拟,实质是要给出随机变量的模拟,也就是说要利用计算机随机产生一系列数值,使它们服从一定的概率分布服从一定的概率分布,称这些数值为随机数.l 最基本,最常用的是(0,1)区间内均匀分布的随机数.其他分布的随机数均可利用它来产生.1:产生模拟随机数的计算机命令l 在MATLAB中,可以直接产生满足各种分布的随机数,命令如下:l 常见的分布函数 MATLAB语句 l 均匀分布U0,1 R=rand(m,n)l 均匀分布Ua,b R=unifrnd(a,b,m,n)l 指数分布 E()R=exprnd(,m,n)l 正态分布N(mu,sigma)R=normrnd(mu,sigma,m,n)l 标准正态分布N(0,1)R=randn(m,n)l 二项分布B(n,p)R=binornd(n,p,m,n1)l 泊松分布 P()R=poissrnd(,m,n)l 以上语句均产生m n 的矩阵2:案例分析l 例1:unifrnd(2,3)l unifrnd(1,32,1,4)l normrnd(1,2)l normrnd(1,2,2,3)l rand(2,3)l randn(2,3)l ans=2.8132l ans=1.3057 5.3056 7.2857 7.1604l ans=0.2527l ans=l 2.7429 0.0219 2.7759l 2.2756 0.0992-0.9560l ans=l 0.6038 0.1988 0.7468l 0.2722 0.0153 0.4451l ans=l-0.0945-1.3089-0.2440l-0.2141 0.8248-0.1778l 2:案例分析l 例2:敌空战部队对我方港口进行空袭,其到达规律服从泊松分布,平均每分钟到达4架飞机.l(1)模拟敌机在3分钟内到达目标区域的数量,以及在第1,2,3分钟内各到达几架飞机;l(2)模拟在3分钟内每架飞机的到达时刻.分析:(1)n1=poissrnd(4),n2=poissrnd(4),n3=poissrnd(4),n=n1+n2+n3(2)由排队论知识,敌机到达规律服从泊松分布等价于敌机到达港口的间隔时间服从参数为1/4的指数分布,故可由指数分布模拟每架飞机的到达时刻.注：如果单位时间发生的次数（如到达的人数）服从参数为r的泊松分布,则任连续发生的两次时间的间隔时间序列服从参数为r的指数分布!泊松分布的期望是，根据到泊松分布和指数分布的关系，可以推出指数分布的期望是1/。2:案例分析l clearl t=0;l j=0;%到达的飞机数 l while t3l j=j+1l t=t+exprnd(1/4)l end二:随机变量的模拟l 在很多实际问题中，我们需要模拟服从一定分布的随机变量,来进行计算和预测。l 利用均匀分布的随机数可以产生具有任意分布的随机变量的样本,从而可以对随机变量的取值情况进行模拟.l 1 连续型随机变量的模拟l 具有给定分布的连续型随机变量可以利用在区间(0,1)上均匀分布的随机数来模拟,最常用的方法是逆变换法.l 结论:若随机变量Y有连续的分布函数F(y),l 则Z与Y有相同的分布.1 连续型随机变量的模拟 一般说来，具有给定分布的连续型随机变量可以利用在区间(0,1)上均匀分布的随机数来模拟最常用的方法是反函数法由概率论的理论可以证明，若随机变量Y有连续的分布函数F(y)，而X是区间(0,1)上均匀分布的随机变量，令，则Z与Y有相同的分布由此，若已知Y的概率密度为，由 可得反函数法l 如果给定区间(0,1)上均匀分布的随机数,则具有给定分布Y的随机数 可由方程 l 解出.l 例:模拟服从参数为 的指数分布时,由 可得 注：指数分布的密度函数2 离散型随机变量的模拟l 设随机变量X的分布律为:有相同的发生的概率.因此我们可以用随机变量R落在小区间内的情况来模拟离散的随机变量X的取值情况.因此我们可以用随机变量R落在小区间内的情况来模拟离散的随机变量X的取值情况具体执行的过程是：产生一个(0,1)上均匀分布的随机数 r(简称随机数)，若p(n-1)r p(n)则理解为发生事件(X=xn)于是就可以模拟随机变量的取值情况2 离散型随机变量的模拟l 例 3:随机变量 表示每分钟到达银行柜台的顾客数.X的分布列见下表,试模拟10分钟内顾客到达柜台的情况.l 表1 10分钟内顾客到达柜台的情况l Xk 0 1 2l pk 0.4 0.3 0.3l 分析:因为每分钟到达柜台的人数是随机的,所以可用计算机随机生成一组(0,1)的数据,由X的概率分布情况,可认为随机数在(0,0.4)范围内时没有顾客光顾,在0.4,0.7)时,有一个顾客光顾,在0.7,1)时,有两个顾客光顾.l 从而有MATLAB程序:2 离散型随机变量的模拟l r=rand(1,10);l for i=1:10;l if r(i)0.4l n(i)=0;l elseif 0.4=r(i)&r(i)0.7l n(i)=1;l else n(i)=2;l end;l endl rl n三:连续系统的模拟-时间步长法l 对连续系统的计算机模拟是近似地获取系统状态在一些离散时刻点上的数值在一定假设条件下，利用数学运算模拟系统的运行过程连续系统模型一般是微分方程，它在数值模拟中最基本的算法是数值积分算法例如有一系统可用微分方程来描述:已知输出量y的初始条件，现在要求出输出量y随时间变化的过程y(t)。1 理论介绍l 最直观的想法是：首先将时间离散化，令l，称为第k步的计算步距（一般是等间距的），然后按以下算法计算状态变量在各时刻 上的近似值:l 其中初始点 按照这种作法即可求出整个的曲线这种最简单的数值积分算法称为欧拉法除此之外，还有其他一些算法 1 理论介绍l 因此，连续系统模拟方法是：首先确定系统的连续状态变量，然后将它在时间上进行离散化处理，并由此模拟系统的运行状态模拟过程分为许多相等的时间间隔,时间步长的长度可以根据实际问题分别取为秒,分,小时,天等.仿真时钟按时间步长等距推进,每次推进都要扫描系统中所有活动,按预定计划和目标进行分析,计算,记录系统状态的变化,直到满足某个终止条件为止.计算机仿真举例例4：追击问题 我辑私舰雷达发现距cc里处里处有一艘走私船正以匀速度速度aa沿直线行驶,辑私舰立即以最大的速度速度 bb追赶,若用雷达进行跟踪,保持船的瞬时速度方向始终指向走私船,试求辑私舰追逐路线和追上的时间,并且给出缉私舰和走私船的路线图。选取坐标系如图:假设缉私艇初始位置在点(c,0),走私船初始位置在点(0,0),走私船的行驶方向为y方向.Oxy(c,0)D(x,y).R(0,at)设缉私艇为动点D,走私船为动点R.在时刻t,缉私艇的位置是D(x,y),走私船的位置是R(0,at).取时间间隔(步长)为,则在时刻 t+,D的位置是,追赶方向(D的运动方向)为DR.用方向余弦表示为:(*)计算机仿真举例算法:赋初值:初始时刻 t0,时间步长,速度a,b,初始位置c循 环:由公式(*)计算动点D 在各时刻的坐标 D 计算动点R在各时刻的坐标R.终 止:当D,R 的距离小于给定值 时终止仿真算法：第二步：计算动点缉私艇D在时刻 时的坐标 计算走私船R在时刻 时的坐标 第一步：设置时间步长,速度a,b及初始位置第三步：计算缉私艇与走私船这两个动点之间的距离：根据事先给定的距离，判断缉私艇是否已经追上了走私船，从而判断退出循环还是让时间产生一个步长，返回到第二步继续进入下一次循环；第四步：当从上述循环退出后，由点列 和可分别绘制成两条曲线即为缉私艇和走私船走过的轨迹曲线。l试求缉私舰追逐路线和追上的时间。l取c=3千米，a=0.4千米/分钟，b=0.8千米/分钟计算机仿真举例将上述算法中求出的点D用直线段连接起来，就得到追赶曲线。循环终止时的即为追赶时间。可以看出，步长选取得越小，所求的曲线越精确。例5:追逐问题 1问题提出 假设正方形ABCD的四个顶点处各站一人在某一时刻，四人同时以匀速v沿顺时针方向追逐下一个人，并且在任意时刻他们始终保持追逐的方向是对准追逐目标，例如，A追逐B，任意时刻A始终向着B追可以证明四人的运动轨迹将按螺旋曲线状汇合于中心O 怎样证明呢？有两种证明方法一是分别求出四人的运动轨迹曲线解析式，求证四条曲线在某时刻相交于一点另一方法则是用计算机模拟将四人的运动轨迹直观地表示在图形上.2 应用举例2建立模型及模拟方法模拟步骤：1）建立平面直角坐标系l 2）以时间间隔t进行采样，在每一时t计算每个人在下一时t+t时的坐标l 3）不妨设甲的追逐对象是乙，在时间t时，甲2建立模型及模拟方法l 4)选取足够小的，模拟到 时为止(即已追到）5）连接四人在各时刻的位置，就得到所求的轨迹l 根据以上模拟步骤，编出MATLAB程序如下：l%取v=1,t=12,A,B,C,D点的坐标分另为（0，10），（10，10），（10，0），(0,0)l v=1;l dt=0.05;l d=20;l x=0 0 0 10 10 10 10 0;%DABC点x,y坐标l x(9)=x(1);l x(10)=x(2);%D点l holdl axis(equal)%各坐标轴刻度增量相同l axis(0 10 0 10);%坐标轴范围3 MATLAB实现l for k=1:2:7l plot(x(k),x(k+1),.)l endl while(d0.1)l for i=1:2:7%循环求每个点的下一个时间点的坐标l d=sqrt(x(i)-x(i+2)2+(x(i+1)-x(i+3)2);%与追击点的距离l x(i)=x(i)+v*dt*(x(i+2)-x(i)/d;%下一个时间点的x坐标l x(i+1)=x(i+1)+v*dt*(x(i+3)-x(i+1)/d;%下一个时间点的y坐标l plot(x(i),x(i+1),.)l endl x(9)=x(1);x(10)=x(2);l endl hold3 MATLAB实现3 MATLAB实现例6 水池含盐量问题l 某水池有2000m3水,其中含盐2kg,以6m3/min的速率向水池内注入含盐为0.5kg/m3的盐水,同时又以4m3/min的速率从水池流出搅拌均匀的盐水.试用计算机仿真该水池内盐水的变化过程,并每隔10min计算水池中水的体积,含盐量,含盐率.欲使池中盐水含盐率达到0.2kg/m3,需经过多长时间?l 分析:这是一个连续系统,首先要将系统离散化,在一些离散点上进行考察,这些离散点的间隔就是时间步长.可取步长为1min,即隔1min考察一次系统的状态,并相应地记录和分析.在注入和流出的作用下,池中水的体积与含盐量,含盐率均随时间变化,初始时刻含盐率为0.001kg/m3,以后每分钟注入含盐率为0.5kg/m3的水6m3,流出混合均匀的盐水为4m3,当池中水的含盐率达到0.2kg/m3时,仿真过程结束.例6 水池含盐量问题l 记T时刻的体积为w(m3)，水的含盐量为s(kg)，水的含盐率为r=s/w(kg/m3)，每隔1min池水的动态变化过程如下：每分钟水的体积增加6-4=2(m3)；每分钟向池内注入盐60.5=3(kg)；每分钟向池外流出盐4r(kg)；每分钟池内增加盐3-4r(kg).l 本例还可以用微分方程建立数学模型，并求出它的解析解，这个解析解就是问题的精确解，有兴趣的读者可以按照这个思路求出该问题的精确解，考察相应时刻精确解与仿真解的差异，还可以进一步调整仿真过程的时间步长，通过与精确解的比较来研究时间步长的大小对仿真度的影响。MATLAB实现l clearl h=1;%时间步长为1l s0=2;%初始含盐2kgl w0=2000;%初始水池有水2000m3l r0=s0/w0;%初始浓度l s(1)=s0+0.5*6*h-4*h*r0;%一分钟后的含盐量l w(1)=w0+2*h;%一分钟后水池中的盐水体积l r(1)=s(1)/w(1);%一分钟后的浓度 l t(1)=h;l y(1)=(2000000+3000000*h+3000*h2+h3)/(1000+h)2;l for i=2:200l t(i)=i*h;l s(i)=s(i-1)+0.5*6*h-4*h*r(i-1);%第i步后的含盐量l w(i)=w(i-1)+2*h;%第i步后的盐水体积l r(i)=s(i)/w(i);%第i步后的盐水浓度l y(i)=(2000000+3000000*t(i)+3000*t(i)2+t(i)3)/(1000+t(i)2;l m=floor(i/10);MATLAB实现l if i/10-m0.2%若第i步后的盐水浓度大于0.2l t02=i*h;l r02=r(i);l break l endl endl t02,r02l 10*tm,sm,rm%表示逆l subplot(1,2,1),plot(t,s,blue);l hold onl subplot(1,2,2),plot(t,y,red);应用举例-面积计算例7:面积计算 求由曲线x2+y2=16,x2/36+y2=1,以及(x-2)2+(y+1)2=9所围成图形的面积。用Matlab语句作出图形得区域图 t=0:0.01:2*pi;x=sin(t);y=cos(t);plot(4*x,4*y,6*x,y,2+3*x,3*y-1)axis(-6,6,-6,6)应用举例-面积计算应用举例-面积计算此区域形状复杂，理论分析困难，可以用计算机仿真实现。将可能的区域等分，考察每个小区域是否在此区域中，将在此区域中的小面积相加即可其仿真图如下：应用举例-面积计算Matlab程序为x=-2:0.01:6;y=-2:0.01:2;s=0;h=0.01;for i=1:800 for j=1:400 xx=-2+i*h;yy=-2+j*h;if xx2+yy2=16 if xx2/36+yy2=1 if(xx-2)2+(yy+1)2=9 s=s+h2;end end end endends运行后给出面积的值 8.8310。四:离散系统的模拟-事件步长法 l 离散系统(discrete system)是指系统状态只在有限的时间点或可数的时间点上有随机事件驱动的系统例如排队系统（queue system），显然状态量的变化只是在离散的随机时间点上发生假设离散系统状态的变化是在一个时间点上瞬间完成的 l 常用的是事件步长法（下次事件推进法）其过程是：置模拟时钟的初值为0，跳到第一个事件发生的时刻，计算系统的状态，产生未来事件并加入到队列中去，跳到下一事件，计算系统状态，重复这一过程直到满足某个终止条件为止例7 收款台前的排队过程l 假设:(1)顾客到达收款台是随机的,平均时间间隔为0.5min,即间隔时间服从lamda=2的指数分布.l(2)对不同的顾客收款和装袋的时间服从正态分布N(1,1/3).l 试模拟20位顾客到收款台前的排队情况,我们关心的问题是每个顾客的平均等待时间,队长及服务员的工作效率例7 收款台前的排队过程l 分析:单服务台结构的排队系统有两类原发事件即到来和离去,顾客到来的后继事件是顾客接受服务,顾客离去的后继事件是服务台寻找服务,这四类事件各自的子程序框图如图1所示.l 假设:t(i)为第i位顾客到达时刻;t2(i)为第i位顾客受到服务的时间(随机变量);T(i)为第i位顾客离去时刻.l 将第i位顾客到达作为事件发生:t(i+1)-t(i)=r(i)(随机变量);平衡关系:当t(i+1)T(i)时,T(i+1)=t(i+1)+t2(i+1);否则,T(i+1)=T(i)+t2(i+1).图1 到来事件子程序l 系统人数+1产生下一个顾客到来时刻调用接受服务事件的程序图 2 离去事件子程序 系统人数-1置服务台”闲”已服务人数+1调用寻找服务事件子程序图3:接受服务事件子程序 服务台空置服务台”忙”产生服务结束时刻登记到事件表排队人数+1否是图4:寻找服务事件子程序 排队人数1置服务台”忙”产生服务结束时刻排队人数-1排队人数+1否是五：蒙特卡洛方法l 在用传统方法难以解决的问题中，有很大一部分可以用概率模型进行描述由于这类模型含有不确定的随机因素，分析起来通常比确定性的模型困难有的模型难以作定量分析，得不到解析的结果，或者是虽有解析结果，但计算代价太大以至不能使用在这种情况下，可以考虑采用Monte Carlo方法。Monte Carlo方法是计算机模拟的基础，它的名字来源于世界著名的赌城摩纳哥的蒙特卡洛，其历史起源于1777年法国科学家蒲丰提出的一种计算圆周的方法随机投针法，即著名的蒲丰投针问题。l 18世纪，法国数学家蒲丰和勒可莱尔提出的“投针问题”，记载于蒲丰1777年出版的著作中：“在平面上画有一组间距为a的平行线，将一根长度为L（L=a/2）的针任意掷在这个平面上，求此针与平行线中任一条相交的概率。”蒲丰本人证明了，这个概率是：l p=2L/（a）为圆周率利用这个公式可以用概率的方法得到圆周率的近似值。下面是一些资料试验者时间 投掷次数 相交次数圆周率估计值Wolf 1850年 5000 2532 3.1596Smith 1855年 3204 1218.5 3.1554C.De Morgan1860年 600 382.5 3.137Fox 1884年 1030 489 3.1595Lazzerini 1901年 3408 1808 3.1415929Reina 1925年 2520 859 3.1795五：蒙特卡洛方法l Monte Carlo方法的基本思想是首先建立一个概率模型，使所求问题的解正好是该模型的参数或其他有关的特征量然后通过模拟一统计试验，即多次随机抽样试验（确定m和n），统计出某事件发生的百分比只要试验次数很大，该百分比便近似于事件发生的概率这实际上就是概率的统计定义利用建立的概率模型，求出要估计的参数蒙特卡洛方法属于试验数学的一个分支五：蒙特卡洛方法l 蒙特卡洛方法适用范围很广泛，它既能求解确定性的问题，也能求解随机性的问题以及科学研究中的理论问题例如利用蒙特卡洛方法可以近似地计算定积分，即产生数值积分问题蒙特卡洛法求圆周率l clearl n=50000l X=rand(n,1);l Y=rand(n,1);l k=0;l for i=1:n;l if X(i)2+Y(i)2=1l k=k+1;l endl endl 4*k/n%1/4*PI*r2=k/n蒲丰投针求圆周率的M程序l function pi_value=pinpi(k,a,l)l%k投针次数l%a线距l%l针长（必须小于等于a）l%pi_value返回pi值l m=0;l if la;l fprintf(error：针长必须小于等于%dn,a);l fprintf(请重新调用函数pinpi(k,d,l)n);l pi_value=0;l else l for i=1:kl a*rand(1)=l*sin(pi*rand(1)%判断投针下沿与最近平行线的距离与投针垂直高度的关系，距离小于等于高度就相交l if m=m+1;l end l end l p=m/k;l pi_value=2*l/(a*p);l fprintf(投针法求得pi=%dn,pi_value);l endl format long gl pinpi(100000,4,3)l 投针法求得pi=3.143797e+000l ans=l 3.14379728795087任意曲边梯形面积的近似计算l 一个古老的问题：用一堆石头测量一个水塘的面积应该怎样做呢？测量方法如下：假定水塘位于一块面积已知的矩形农田之中如图82所示随机地向这块农田扔石头使得它们都落在农田内被扔到农田中的石头可能溅上了水，也可能没有溅上水，估计被“溅上水的”石头量占总的石头量的百分比试想如何利用这估计的百分比去近似计算该水塘面积？任意曲边梯形面积的近似计算任意曲边梯形面积的近似计算l 结合图82中的图形(1)分析，只要已知各种参数及函数（a，b，H，f(x)），有以下两种方法可近似计算水塘面积 l 1随机投点法 1）赋初值：试验次数n=0，成功次数m=0；规定投点试验的总次数N；l 2）随机选择m个数对xi,yi,1im，其中axib,0yiH，置 nnl；l 3）判断nN，若是，转4，否则停止计算；l 4）判断条件(表示一块溅水的石头)是否成立，若成立则置m=m+1，转2，否则转2；l 5）计算水塘面积的近似值S=H(b-a)m/N 任意曲边梯形面积的近似计算l 2平均值估计法l 1）产生a,b区间的均匀随机数x,i=1,2,N l 2)计算f(xi），i=1,2,Nl 3）计算l 该方法的特点是估计函数f(x)在a,b上的平均值，面积近似等于该平均值乘以(b-a)例9 库存问题l 在物资的供应过程中,由于到货和销售不可能做到同步同量,故总要保持一定的库存储备.如果库存过多,就会造成积压浪费以及保管费的上升;如果库存过少,会造成缺货.如何选择库存和订货策略,就是一个需要研究的问题.库存问题有多种类型,一般都比较复杂,下面讨论一种简单的情形.l 某电动车行的仓库管理人员采取一种简单的订货策略,当库存降低到P辆电动车时就向厂家例9 库存问题l 订货,每次订货Q辆,如果某一天的需求量超过了库存量,商店就有销售损失和信誉损失,但如果库存量过多,就会导致资金积压和保管费增加.若现在有如下表所示的两种库存策略,试比较选择一种策略以使总费用最少.l 表 订货方案 方案 重新订货点P/辆 重新订货量Q/辆 方案1 125 150 方案2 150 250例9 库存问题l 这个问题的已知条件是:l(1)从发出订货到收到货物需隔3天.l(2)每辆电动车保管费为0.50元/天,每辆电动车的缺货损失为1.60元/天,每次的订货费为75元.l(3)每天电动车需求量是0到99之间均匀分布的随机数.l(4)原始库存为110辆,假设第一天没有发出订货.例9 库存问题l 分析:这一问题用解析法讨论比较麻烦,但用计算机按天仿真仓库货物的变动情况却很方便.我们以30天为例,依次对这两种方案进行仿真,最后比较各方案的总费用,从而就可以作出决策.l 计算机仿真时的工作流程是早上到货,全天销售,晚上订货,输入一些常数和初始数据后,以一天为时间步长进行仿真.首先检查这一天是否为预订到货日期,如果是,则原有库存量加Q,并把预定到货量清为0;如果不是,则库存量不变.接着用计算机中的随机函数仿真随机需求量,若库存量大于需求量,则新的库存量减去需求量;反之,则新库存量变为0,并且要在总费用上加缺货损失,然后检查实际例9 库存问题l 库存量加上预定到货量是否小于重新订货点P,如果是,则需要重新订货,这时就加一次订货费.如此重复运行30天,即可得所需费用总值.由此比较这两种方案的总费用,即得最好方案.l clearl days=30;l P=125,150;%重新订货点P，单位辆l Q=150,250;%重新订货量Q，单位辆l cost=0,0;%花费初值l arrivalinterval=2;%到货时间间隔为2天l storagefee=0.5;lossfee=1.6;bookfee=75;%各种费用l storage0=110;booknumber=0;arrivedate=0;l nr=rand(days,1);%随机产生30天的随机数l for i=1:2%两种方案分别仿真l%下面是第一天的情况 l storage(1)=storage0;l n=round(99*nr(1);%30天中每天的需求量,0-99之间的均匀分布随机数l sale=n;l remain=storage(1)-n;%计算库存l if remain=P(i);%第i种方案，库存少于重新订货量，就需要订货l booknumber=Q(i);%订货数量l arrivedate=4;%到货日期l orderfee=bookfee;%订单费用l elsel orderfee=0;%未出现缺货情况，订单费用为0l endl storage(1)=remain;%更新货存l cost(i)=cost(i)+remain*storagefee+orderfee;%总费用累加，第一天计算结束l for j=2:days%第2天到第30天仿真l dh=j;l if abs(dh-arrivedate)0.01%订单到货l storage(j)=storage(j-1)+booknumber;%货存增加l booknumber=0;%已到货订货数量归0l arrivedate=j;%到货日期l elsel storage(j)=storage(j-1);%未到货，当前库存与前一天相等l end%结束判断是否订单到货l n=round(99*nr(j);%第j天随机需求量nl if storage(j)=n%库存足够l sale=n;%成交量l remain=storage(j)-n;%库存量l shortagenumber=0;%缺货数为0l else%否则，库存不够的处理l sale=storage(j);%只能销售已有货存的量，有缺货l remain=0;%库存卖完l shortagenumber=n-storage(j);%缺货量l endl storage(j)=remain;%更新库存l if remain+booknumber=P(i);%库存量加上预定到货量小于重新订货点P，重新订货l booknumber=Q(i);%订货数l arrivedate=dh+arrivalinterval;%到货日期l orderfee=bookfee;%订单费用l elsel orderfee=0;%无需重新订货，未发生费用l endl cost(i)=cost(i)+remain*storagefee+shortagenumber*lossfee+orderfee;%总费用增加l end;l mincost=min(cost);%求两个方案30天总费用最小值l endl costl mincost例10 赶火车过程仿真l 一列火车从A站经B站开往C站,某人每天赶往B站乘这趟火车.已知火车从A站到B站运行时间为均值30min,标准差为2min的正态随机变量.火车大约在下午1点离开A站,离开时刻的频率分布见表1.这个人到达B站时的频率分布表见表2.用计算机仿真火车开出,火车到达B站,这个人到达B站的情况,并给出能赶上火车的仿真结果.表1 离开时刻的频率分布出发时刻T 1:00 1:05 1:10频率 0.7 0.2 0.1例10 赶火车过程仿真l 表2 人到达B站时的频率分布l 到达时刻T 1:28 1:30 1:32 1:34l 频率:0.3 0.4 0.2 0.1l 分析:引入以下变量:T1为火车从A站开出的时刻;T2为火车从A站运行到B站所需要的时间;T3为此人到达B站的时刻.则有表3.l 表3 T1,T2的频率分布l T1,/min 0 5 10 T2/min 28 30 32 34l p 0.7 0.2 0.1 0.3 0.4 0.2 0.1例10 赶火车过程仿真l 显然,这位旅客要赶上火车的条件是T3 T1+T2,可以通过计算机模拟出这三个时间,再检验是否满足T3 T1+T2.如果满足,即能够赶上火车,若不然,则不能赶上火车.火车运行时间的仿真程序x=randn(10000,1);%实验10000次for i=1:10000%每次实验 y(i)=30+2*x(i);%A站运行到B站所需要的时间end开车时间的仿真程序l s1=0;s2=0;s3=0;l x=rand(10000,1);l for i=1:10000l if x(i)0.9l s3=s3+1;l endl endl s1/10000,1-s1/10000-s3/10000,s3/10000%三种开车时间出现的比率 人到达时刻的仿真程序l s1=0;s2=0;s3=0;s4=0;l x=rand(10000,1);l for i=1:10000l if x(i)0.3l s1=s1+1;l elseif x(i)0.7l s2=s2+1;l elsel if x(i)0.9l s3=s3+1;l else l s4=s4+1;l endl endl endl s1/10000,s2/10000,s3/10000,s4/10000%人四个时间点到达车站的比率赶上火车的仿真程序l s=0;l x1=rand(10000,1);l x2=rand(10000,1);l x3=randn(10000,1);l for i=1:10000l if x1(i)0.7%下午1点准时发车l T1=0;l elseif x1(i)0.9%延迟5分钟发车l T1=5;l elsel T1=10;%延迟10分钟发车l endl T2=30+2*x3(i);%A站运行到B站所需要的时间l if x2(i)0.3%此人到达B站的时刻，下午1：28l T3=28;l elseif x2(i)0.7%此人到达B站的时刻，下午1：30l T3=30;l elseif x2(i)0.9%此人下午1：32到达l T3=32;l else l T3=34;%此人下午1：34到达l endl if T3T1+T2%能赶上火车l s=s+1;%累加成功赶上火车的次数l endl endl s/10000%成功赶上火车的概率仿真方法评论l 计算机仿真的精度受到多方面因素的影响,比较难以控制.因此,在能用解析方法求解时,通常不用计算机仿真.l 在科技发展一日千里的今天,计算机仿真方法的应用仍然是非常广泛并且不可替代.展望将来,计算机方法将在人工智能领域,军事领域,医学领域,社会领域有着广泛的应用前景.