1、2019/6/13,Lxy, China Jiliang Universty,1,数学建模专题三,Monte Carlo 模拟,2019/6/13,Lxy, China Jiliang Universty,2,内容提纲,1.引言 2.Monte Carlo模拟基本思想 3.随机数生成函数 4.应用实例举例 5.排队论模拟 6.Monte Carlo模拟求解规划问题,2019/6/13,Lxy, China Jiliang Universty,3,引言(Introduction),Monte Carlo方法：,蒙特卡罗方法，又称随机模拟方法，属于计算数学的一个分支，它是在上世纪四十年代中期为了适
2、应当时原子能事业的发展而发展起来的。 亦称统计模拟方法，statistical simulation method 利用随机数进行数值模拟的方法,Monte Carlo名字的由来：,Monte Carlo是摩纳哥（monaco)的首都，该城以赌博闻名,Nicholas Metropolis (1915-1999),Monte-Carlo, Monaco,2019/6/13,Lxy, China Jiliang Universty,4,二十世纪四十年代中期，由于科学技术的发展和电子计算机的发明，蒙特卡罗方法作为一种独立的方法被提出来，并首先在核武器的试验与研制中得到了应用。但其基本思想并非新颖，
3、人们在生产实践和科学试验中就已发现，并加以利用。,一个简单实例,Monte Carlo方法的基本思想,2019/6/13,Lxy, China Jiliang Universty,5,举例,Buffon投针实验,1768年，法国数学家Comte de Buffon利用投针实验估计的值,2019/6/13,Lxy, China Jiliang Universty,6,Solution,The positioning of the needle relative to nearby lines can be described with a random vector which has comp
4、onents:,The random vector is uniformly distributed on the region 0,d)0,). Accordingly, it has probability density function 1/d.,The probability that the needle will cross one of the lines is given by the integral,2019/6/13,Lxy, China Jiliang Universty,7,例1. 蒲丰投针问题,利用关系式：求出值其中为投计次数，n为针与平行线相交次数。这就是古典概
5、率论中著名的蒲丰氏问题。,2019/6/13,Lxy, China Jiliang Universty,8,一些人进行了实验，其结果列于下表 ：,2019/6/13,Lxy, China Jiliang Universty,9,基本思想,由上面的例子可以看出，当所求问题的解是某个事件的概率，或者是某个随机变量的数学期望，或者是与之有关的量时，通过某种试验的方法，得出该事件发生的频率，再通过它得到问题的解。这就是蒙特卡罗方法的基本思想。,蒙特卡罗方法的基本思想虽然早已被人们提出，却很少被使用。直到电子计算机出现后，使得人们可以通过电子计算机来模拟巨大数目的随机试验过程，使得蒙特卡罗方法得以广泛地
6、应用，在现代化的科学技术中发挥应有的作用。,2019/6/13,Lxy, China Jiliang Universty,10,例1 在我方某前沿防守地域，敌人以一个炮排（含两门火炮）为单位对我方进行干扰和破坏为躲避我方打击，敌方对其阵地进行了伪装并经常变换射击地点,经过长期观察发现，我方指挥所对敌方目标的指示有50是准确的，而我方火力单位，在指示正确时，有1/3的射击效果能毁伤敌人一门火炮，有1/6的射击效果能全部毁伤敌人火炮,现在希望能用某种方式把我方将要对敌人实施的20次打击结果显现出来，确定有效射击的比率及毁伤敌方火炮的平均值。,分析： 这是一个概率问题，可以通过理论计算得到相应的概率
7、和期望值.但这样只能给出作战行动的最终静态结果，而显示不出作战行动的动态过程.,为了能显示我方20次射击的过程，现采用模拟的方式。,举例,2019/6/13,Lxy, China Jiliang Universty,11,需要模拟出以下两件事：,2 当指示正确时，我方火力单位的射击结果情况,1 观察所对目标的指示正确与否,模拟试验有两种结果，每一种结果出现的概率都是1/2,因此，可用投掷一枚硬币的方式予以确定，当硬币出现正面时为指示正确，反之为不正确,模拟试验有三种结果：毁伤一门火炮的可能性为1/3(即2/6)，毁伤两门的可能性为1/6，没能毁伤敌火炮的可能性为1/2(即3/6),这时可用投掷
8、骰子的方法来确定： 如果出现的是、三个点：则认为没能击中敌人； 如果出现的是、点：则认为毁伤敌人一门火炮； 若出现的是点：则认为毁伤敌人两门火炮,问题分析,2019/6/13,Lxy, China Jiliang Universty,12,i：要模拟的打击次数； k1：没击中敌人火炮的射击总数； k2：击中敌人一门火炮的射击总数； k3：击中敌人两门火炮的射击总数； E：有效射击比率； E1：20次射击平均每次毁伤敌人的火炮数,符号说明,2019/6/13,Lxy, China Jiliang Universty,13,模拟框图,2019/6/13,Lxy, China Jiliang Uni
9、versty,14,模拟结果,2019/6/13,Lxy, China Jiliang Universty,15,2019/6/13,Lxy, China Jiliang Universty,16,理论计算,2019/6/13,Lxy, China Jiliang Universty,17,结果比较,虽然模拟结果与理论计算不完全一致，但它却能更加真实地表达实际战斗动态过程,用蒙特卡洛方法进行计算机模拟的步骤：,1 设计一个逻辑框图，即模拟模型 2 根据流程图编写程序，模拟随机现象可通过具有各种概率分布的模拟随机数来模拟随机现象 3 分析模拟结果，计算所需要结果.,2019/6/13,Lxy,
10、China Jiliang Universty,18,注：rand(n)=rand(n,n),Matlab 中的随机数生成函数,2019/6/13,Lxy, China Jiliang Universty,19,name 的取值可以是,norm or Normal unif or Uniform poiss or Poisson beta or Beta exp or Exponential gam or Gamma geo or Geometric unid or Discrete Uniform. .,random(name,A1,A2,A3,M,N),Matlab 中的随机数生成函数,2
11、019/6/13,Lxy, China Jiliang Universty,20,fix(x) : 截尾取整，直接将小数部分舍去 floor(x) : 不超过 x 的最大整数 ceil(x) : 不小于 x 的最小整数 round(x) : 四舍五入取整,Matlab中的取整函数,2019/6/13,Lxy, China Jiliang Universty,21,随机投掷均匀硬币，验证国徽朝上与朝下的概率是否都是 1/2,n=10000; % 给定试验次数 m=0; for i=1:nx=randperm(2)-1;y=x(1);if y=0 % 0 表示国徽朝上，1 表示国徽朝下m=m+1;
12、end end fprintf(国徽朝上的频率为：%fn,m/n);,小实例一：投掷硬币,2019/6/13,Lxy, China Jiliang Universty,22,随机投掷骰子，验证各点出现的概率是否为 1/6,n=10000; m1=0; m2=0; m3=0; m4=0; m5=0;m6=0; for i=1:nx=randperm(6); y=x(1);switch ycase 1, m1=m1+1;case 2, m2=m2+1;case 3, m3=m3+1;case 4, m4=m4+1;case 5, m5=m5+1;otherwise, m6=m6+1;end end
13、 . % 输出结果,小实例二：投掷骰子,2019/6/13,Lxy, China Jiliang Universty,23,用蒙特卡罗 ( Monte Carlo ) 投点法计算 的值,n=100000; a=2; m=0; for i=1:nx=rand(1)*a/2;y=rand(1)*a/2;if ( x2 + y2 = (a/2)2 )m=m+1;end end fprintf(计算出来的 pi 为：%fn,4*m/n);,小实例三：蒙特卡罗投点法,2019/6/13,Lxy, China Jiliang Universty,24,在画有许多间距为 d 的等距平行线的白纸上，随机投掷一
14、根长为 l ( l d ) 的均匀直针，求针与平行线相交的概率，并计算的 值,小实例四：蒲丰投针实验,2019/6/13,Lxy, China Jiliang Universty,25,n=100000; l=0.5; d=1; m=0; for i=1:nalpha=rand(1)*pi;y=rand(1)*d/2;if y=l/2*sin(alpha)m=m+1;end end fprintf(针与平行线相交的频率为：%fn,m/n); fprintf(计算出来的 pi 为：%fn,2*n*l/(m*d);,小实例四源程序,2019/6/13,Lxy, China Jiliang Univ
15、ersty,26,设某班有 m 个学生，则该班至少有两人同一天生日的概率是多少？,小实例五：生日问题,2019/6/13,Lxy, China Jiliang Universty,27,n=1000; p=0; m=50; % 设该班的人数为 50 for t=1:na=; q=0;for k=1:mb=randperm(365);a=a,b(1);endc=unique(a);if length(a)=length(c)p=p+1;end end fprintf(至少两人同一天生日的概率为：%fn,p/n);,试验五源程序,2019/6/13,Lxy, China Jiliang Unive
16、rsty,28,clear; m = 50; p1= 1:365; p2= 1:365-m, 365*ones(1,m); p = p1./p2; p = 1- prod(p); fprintf(至少两人同一天生日的概率为：%fn,p);,小实例五的理论值计算,2019/6/13,Lxy, China Jiliang Universty,29,排队问题随机模拟,排队论主要研究随机服务系统的工作过程。,在排队系统中，服务对象的到达时间和服务时间都是随机的。排队论通过对随机服务现象的统计研究，找出反映这些随机现象平均特性的规律指标，如排队的长度、等待的时间及服务利用率。,2019/6/13,Lxy
17、, China Jiliang Universty,30,1 系统的假设： （1） 顾客源是无穷的； （2） 排队的长度没有限制； （ 3）到达系统的顾客按先后顺序依次进入服务， “先到先服务”。,在某商店有一个售货员，顾客陆续来到，售货员逐个地接待顾客当到来的顾客较多时，一部分顾客便须排队等待，被接待后的顾客便离开商店设： 1顾客到来间隔时间服从参数为0.1的指数分布 对顾客的服务时间服从，上的均匀分布 排队按先到先服务规则，队长无限制,假定一个工作日为8小时，时间以分钟为单位。 1模拟一个工作日内完成服务的个数及顾客平均等待时间t 2模拟100个工作日，求出平均每日完成服务的个数及每日顾客
18、的平均等待时间。,单服务员的排队模型模拟,2019/6/13,Lxy, China Jiliang Universty,31,w：总等待时间； ci：第i个顾客的到达时刻； bi：第i个顾客开始服务时刻； ei：第i个顾客服务结束时刻； xi：第i-1个顾客与第i个顾客之间到达的间隔时间； yi：对第i个顾客的服务时间。,符号说明,2019/6/13,Lxy, China Jiliang Universty,32,c1,b1,c3,c4,c5,c2,e1,b2,e2,b3,e3,b4,e4,b5,ci=ci-1+ xi ei=bi+yi bi=max(ci,ei-1),t,思路分析,2019/
19、6/13,Lxy, China Jiliang Universty,33,1模拟一日 To Matlab(simu1),2模拟100日 To Matlab(simu2),流程框图,2019/6/13,Lxy, China Jiliang Universty,34,用蒙特卡洛法解非线性规划问题,2019/6/13,Lxy, China Jiliang Universty,35,基本假设,试验点的第j个分量xj服从aj ,bj内的均匀分布,符号假设,求解过程,先产生一个随机数作为初始试验点,以后则将上一个试验点的第j个分量随机产生，其它分量不变而产生一新的试验点这样，每产生一个新试验点只需一个新的
20、随机数分量当KMAXK或PMAXP时停止迭代,2019/6/13,Lxy, China Jiliang Universty,36,2019/6/13,Lxy, China Jiliang Universty,37,在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); %转化为求最小值问题,2019/6/13,Lxy, China Jiliang Universty,38,%
21、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; for i=1:nx1=r1(i);x2=r2(i);lpc=lpconst(x1 x2)；if lpc=1z=mylp(x1 x2)；if zz0z0=z；sol=x1 x2；endend end,To Matlab(randlp),返回,2019/6/13,Lxy, China Jiliang Universty,39,Thanks!,
