蒙特卡罗方法教材课件

1、数学建模专题之蒙特卡罗方法主讲：李培峦时间：2014-08-032022/7/291内容提纲1.引言2.Monte Carlo模拟基本思想3.随机数生成函数4.应用实例举例5.随机游走模拟6.机械零件的可靠度计算7.排队论模拟8.Monte Carlo方法求解规划问题9.Monte Carlo方法预测搜救区域2022/7/292引言(Introduction)Monte Carlo方法：蒙特卡罗方法，又称随机模拟方法，属于计算数学的一个分支，它是在上世纪四十年代中期为了适应当时原子能事业的发展而发展起来的。亦称统计模拟方法，statistical simulation method 利用随机数

2、进行数值模拟的方法Monte Carlo名字的由来：Monte Carlo是摩纳哥（monaco)最大的城市，该城以赌博闻名。Monte-Carlo, Monaco2022/7/293引言(Introduction)Monte Carlo方法的起源：二十世纪四十年代中期，由于科学技术的发展和电子计算机的发明，蒙特卡罗方法作为一种独立的方法被提出来，并首先在核武器的试验与研制中得到了应用(中子的链锁反应)。但其基本思想并非新颖，可追溯到18世纪后半叶的蒲丰（Buffon）随机投针实验（1777年）。Nicholas Metropolis (1915-1999)John Von Neumann (

3、1903-1957) 2022/7/294引言(Introduction)Monte Carlo方法的应用： 物理：核物理，热力学与统计物理，粒子输运问题等 数学：多重积分、解微分方程、非线性方程组求解等 工程领域：真空技术，水力学，激光技术等 经济学领域：期权定价、项目管理、投资风险决策等 其他领域：化学、医学，生物，生产管理、系统科学、公用事业等方面。随着科学技术的发展，其应用范围将更加广泛。2022/7/295Monte Carlo方法的基本思想Buffon投针实验1768年，法国数学家Comte de Buffon利用投针实验估计的值dL2022/7/296SolutionThe po

4、sitioning of the needle relative to nearby lines can be described with a random vector which has components: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

5、the integral2022/7/297例1. 蒲丰投针问题利用关系式：求出值 其中 为投计次数，n 为针与平行线相交次数。这就是古典概率论中著名的蒲丰氏问题。2022/7/298一些人进行了实验，其结果列于下表 ：实验者年份投计次数的实验值沃尔弗(Wolf)185050003.1596斯密思(Smith)185532043.1553福克斯(Fox)189411203.1419拉查里尼(Lazzarini)1901340837/299基本思想 MS基本思想：为了求解数学、物理、工程技术或随机服务系统等方面的问题，首先构造一个模型（概率统计模型），使所求问题的解正好

6、是该模型的参数或特征量或有关量，然后通过模拟（统计试验），给出模型参数或特征量的估计值，最后得出所求问题的近似解。 MS特点： 1.方法新颖、应用面广、实用性强； 2.随机模拟方面的算法简单，但计算量大； 3.模拟结果具有随机性，精度较低； 4.模拟结果的收敛过程服从概率规律。2022/7/2910例1在我方某前沿防守地域，敌人以一个炮排（含两门火炮）为单位对我方进行干扰和破坏为躲避我方打击，敌方对其阵地进行了伪装并经常变换射击地点 经过长期观察发现，我方指挥所对敌方目标的指示有50是准确的，而我方火力单位，在指示正确时，有1/3的射击效果能毁伤敌人一门火炮，有1/6的射击效果能全部毁伤敌人火

7、炮 现在希望能用某种方式把我方将要对敌人实施的20次打击结果显现出来，确定有效射击的比率及毁伤敌方火炮的平均值。分析：这是一个概率问题，可以通过理论计算得到相应的概率和期望值.但这样只能给出作战行动的最终静态结果，而显示不出作战行动的动态过程. 为了能显示我方20次射击的过程，现采用模拟的方式。举例2022/7/2911 需要模拟出以下两件事： 2 当指示正确时，我方火力单位的射击结果情况1 观察所对目标的指示正确与否模拟试验有两种结果，每一种结果出现的概率都是1/2 因此，可用投掷一枚硬币的方式予以确定，当硬币出现正面时为指示正确，反之为不正确 模拟试验有三种结果：毁伤一门火炮的可能性为1/

8、3(即2/6)，毁伤两门的可能性为1/6，没能毁伤敌火炮的可能性为1/2(即3/6) 这时可用投掷骰子的方法来确定：如果出现的是、三个点：则认为没能击中敌人；如果出现的是、点：则认为毁伤敌人一门火炮；若出现的是点：则认为毁伤敌人两门火炮问题分析2022/7/2912i：要模拟的打击次数； k1：没击中敌人火炮的射击总数； k2：击中敌人一门火炮的射击总数；k3：击中敌人两门火炮的射击总数；E：有效射击比率； E1：20次射击平均每次毁伤敌人的火炮数符号说明2022/7/2913模拟框图初始化:i=0,k1=0,k2=0,k3=0i=i+1骰子点数?k1=k1+1k2=k2+1k3=k3+1k1

9、=k1+1i20?E=(k2+k3)/20 E1=(0*k1+1*k2+2*k3)/20停止硬币正面?YNNY1，2，34，562022/7/2914模拟结果2022/7/2915从以上模拟结果可计算出：2022/7/2916理论计算2022/7/2917结果比较 虽然模拟结果与理论计算不完全一致，但它却能更加真实地表达实际战斗动态过程要使结果接近理论计算值，必须加大模拟次数，这就要求使用计算机模拟了。 用蒙特卡洛方法进行计算机模拟的步骤：1 设计一个逻辑框图，即建立模拟模型2 根据流程图编写程序，模拟随机现象可通过生成具有各种概率分布的随机数来模拟随机现象，进行模拟试验3 分析模拟结果，给出

10、所求问题的近似解（求解）.2022/7/2918注：rand(n)=rand(n,n)Matlab 中的随机数生成函数randperm(m) 生成一个由 1:m 组成的随机排列randn(m,n) 生成一个满足正态 m n 的随机矩阵rand(m,n) 生成一个满足均匀分布的 m n 随机矩阵，矩阵的每个元素都在 (0,1) 之间。perms(1:n) 生成由 1:n 组成的全排列，共 n! 个结果2022/7/2919 name 的取值可以是norm or Normalunif or Uniformpoiss or Poissonbeta or Betaexp or Exponentialg

11、am or Gammageo or Geometricunid or Discrete Uniform . .random(name,A1,A2,A3,M,N)Matlab 中的随机数生成函数2022/7/2920fix(x) : 截尾取整，直接将小数部分舍去floor(x) : 不超过 x 的最大整数ceil(x) : 不小于 x 的最小整数round(x) : 四舍五入取整Matlab中的取整函数2022/7/2921 模拟随机投掷均匀硬币，验证国徽朝上与朝下的概率是否都是 1/2 n=10000; % 给定试验次数m=0; % m 表示试验成功（国徽朝上）的次数for i=1:n x=r

12、andperm(2)-1; % randperm(2）生成1和2的随机排列 y=x(1); if y=0 % 0 表示国徽朝上，1 表示国徽朝下 m=m+1; endendfprintf(国徽朝上的频率为：%fn,m/n);小实例一：投掷硬币2022/7/2922 随机投掷骰子，验证各点出现的概率是否为 1/6 n=10000; m1=0; m2=0; m3=0; m4=0; m5=0;m6=0;for i=1:n x=randperm(6); y=x(1); switch y case 1, m1=m1+1; case 2, m2=m2+1; case 3, m3=m3+1; case 4,

13、 m4=m4+1; case 5, m5=m5+1; otherwise, m6=m6+1; endend. % 输出结果小实例二：投掷骰子2022/7/2923 用蒙特卡罗 ( Monte Carlo ) 投点法计算 的值 n=10000; a=2; m=0; % m表示落入圆内的次数for i=1:n x=rand(1)*a/2; y=rand(1)*a/2; if ( x2 + y2 = (a/2)2 ) m=m+1; endendfprintf(计算出来的 pi 为：%fn,4*m/n);小实例三：蒙特卡罗投点法2022/7/2924小实例三：蒙特卡罗投点法ezplot(x2+y2-1

14、,-1.1 1.1);hold onaxis equalplot(-1 -1 1 1 -1,-1 1 1 -1 -1);N=0;for k=1:100 N_point =10000; xy =(rand(2,N_point)-0.5)*2; d=sqrt(xy(1,:).2+xy(2,:).2); N(k) = length(find(d=1);endp1 = sum(N)/( N_point*k)p2 = pi*12/4pi0 =p1*4p_xy =(rand(500,2)-0.5).*2;plot(p_xy(:,1),p_xy(:,2),.);2022/7/2925 在画有许多间距为 d

15、的等距平行线的白纸上，随机投掷一根长为 l ( l d ) 的均匀直针，求针与平行线相交的概率，并计算的 值小实例四：蒲丰投针实验2022/7/2926n=100000; L=0.5; d=1; m=0; for i=1:n alpha=rand(1)*pi; y=rand(1)*d/2; if y=yy), 2) + 1; X(i) = x(d);endhist(X, x) % 思考：如何生产指定分布随机数?2022/7/2930随机游走模拟随机游走是一种基于运用0，1区间的均匀分布随机数序列来进行的计算。利用蒙特卡罗方法多次模拟醉汉行走，统计在 X 处占总次数的比值，即为概率 PN(x)。

16、2022/7/2931随机游走模拟设定朝右走的概率P, 总步数N, 模拟次数M , X , j=1 , S=0;i=1; x=0;i=Nj=M产生随机数确定行走方向计算所在位置xi=i+1j=j+1X=xS=S+1YNYNYN计算S/Mi 步数，j 模拟的次数，S为M次中成功次数x表示N步后实际位置, X表示指定的位置2022/7/2932程序P=0.5; N=10; M=100; X=0;S=0;for j =1:M x=0; for i=1:N if(rand(1)=P) x = x+1; else x = x-1; end end if (X = x) S=S+1; endendPP=S

17、/M2022/7/2933机械零件的可靠度计算可靠度计算主要研究机械零件在一定分布的随机应力 S（即零件的正常工作应力）作用下的可靠程度。应力和强度都是随机变量，影响应力和强度的各种因素也是随机变量。因此，机械零件的可靠度计算过程中，利用蒙特卡罗方法来处理对随机量的计算问题，具有特有的优势。问题： 设计一个拉杆，若受外力 P 均值为20000N，标准差 2000N，拉杆的半径 D 均值为20mm，标准差 0.5mm，材料强度 S 均值412 MPa，标准差15.6 MPa，计算其可靠度。2022/7/2934机械零件的可靠度计算问题分析：2022/7/2935机械零件的可靠度计算理论推导：20

18、22/7/2936机械零件的可靠度计算理论推导：2022/7/2937机械零件的可靠度计算蒙特卡罗模拟：通过多次试验，统计可靠零件个数，求出 其占总零件数的比值，即为可靠度。设定总零件数 N=1000，可靠零件数 St=0; 令 P_mean=20000; P_deta=2000; D_mean = 20; D_deta=0.5; Q_mean=412*106；Q_deta=15.6*106;i=1;i=N根据已知分布产生随机数p, d, q;计算应力s=(4*p)/(pi*d2)YNs=qSt=St+1,i=i+1;计算可靠度：St/N2022/7/2938机械零件的可靠度计算程序：N=10

19、0000; St=0;P_mean = 20000; P_deta=2000;D_mean = 20*10-3; D_deta=0.5*10-3;Q_mean=412*106; Q_deta=150.6*106;for i=1:N p=P_mean + P_deta*randn(1); % p为正态分布随机数 d=D_mean + D_deta*randn(1); q=Q_mean + Q_deta*randn(1); s=(4*p)/(pi*d2); if s=q St=St+1; endendSt_p =St/N;fprintf (可信度为：%f n, St_p );2022/7/2939

20、排队问题随机模拟排队论主要研究随机服务系统的工作过程。在排队系统中，服务对象的到达时间和服务时间都是随机的。排队论通过对随机服务现象的统计研究，找出反映这些随机现象平均特性的规律指标，如排队的长度、等待的时间及服务利用率。2022/7/2940系统的假设：（1） 顾客源是无穷的；（2） 排队的长度没有限制；（3）到达系统的顾客按先后顺序依次进入服务， “先到先服务”。 在某商店有一个售货员，顾客陆续来到，售货员逐个地接待顾客当到来的顾客较多时，一部分顾客便须排队等待，被接待后的顾客便离开商店设： 1顾客到来间隔时间服从参数为 0.1 的指数分布 2对顾客的服务时间服从，15上的均匀分布 3排队

21、按先到先服务规则，队长无限制 假定一个工作日为8小时，时间以分钟为单位。1模拟一个工作日内完成服务的个数及顾客平均等待时间t2模拟100个工作日，求出平均每日完成服务的个数及每日顾客的平均等待时间。单服务员的排队模型模拟2022/7/2941w：总等待时间；ci：第i个顾客的到达时刻；bi：第i个顾客开始服务时刻； ei：第i个顾客服务结束时刻；xi：第i-1个顾客与第i个顾客之间到达的间隔时间；yi ：对第i个顾客的服务时间。符号说明2022/7/2942c1b1c3c4c5c2e1b2e2b3e3b4e4b5ci = ci-1+ xi % 到达时刻、时间间隔ei = bi+yi % 结束服

22、务时刻、开始服务时刻、服务时间bi = max(ci, ei-1)t思路分析2022/7/2943初始化：令i=1，ei-1=0，w=0产生间隔时间随机数xi 参数为0.1的指数分布ci =xi , bi=xi 产生服务时间随机数yi4,15的均匀分布ei =bi + yi累计等待时间：w= w+bi-ci准备下一次服务：i=i+1产生间隔时间随机数xi 参数为0.1的指数分布ci =ci-1+ xi 确定开始服务时间：bi=max(ci, ei-1)bi480?YNi=i-1,t=w/i输出结果：完成服务个数：m=i 平均等待时间：t停止1模拟一日To Matlab(simu1)2模拟100

23、日To Matlab(simu2)b,服务时刻；c,到达时刻；e,结束时刻；2022/7/2944源代码simu1function i , t=simu1i=1; e=0; w=0;x(i) = random(exp,10); % lambda=1/10c(i) = x(i); b(i) = x(i);while ( b(i)MAXK或PMAXP时停止迭代2022/7/2948框 图初始化:给定MAXK,MAXP;k=0,p=0,Q:大整数xj = aj + R(bj-aj) j=1,2,nj=0j=j+1,p=p+1PMAXP?YNxj = aj + R(bj-aj)gi(X)0?i=1,2

24、nYNjMAXK?YN输出X,Q,停止YN2022/7/2949 在 Matlab中编程，共需三个文件: randlp.m ,mylp.m 和 lpconst.m . 主程序为randlp.m .% mylp.mfunction z = mylp(x) %目标函数z=2*x(1)2+x(2)2-x(1)*x(2)-8*x(1)-3*x(2); %转化为求最小值问题% lpconst.mfunction lpc =l pconst(x) %约束条件if 3*x(1)+x(2)2-10-0.5; %约束条件的误差为 lpc =1;else lpc =0;end 2022/7/2950% randl

25、p.m function sol,r1,r2=randlpdebug=1;a=0; %试验点下界b=10; %试验点上界n=1000; %试验点个数r1=unifrnd(a,b,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 endend2022/7/2951蒙特卡洛方法预测最佳搜救区域

26、问题介绍 根据搜救目标最后已知位置和时间,依据目标海上漂流风压特性,综合海区不同时段风速、流速等环境信息及其不确定性影响,预测目标漂移后的搜索区域,随时间推移计算目标漂流轨迹,同时求出某一具体时刻目标在某可能区域的概率。 准确地预测漂浮物的漂浮区域对于及早确定遇险对象搜救区域、成功实施援救具有决定意义。 准确的搜索区域划定包含两个要求: 搜救区域要能以最大概率包含搜救对象,不至于遗漏搜救对象; 搜救区域确定尽可能细致, 尽可能小,使搜救力量可以集中在最短时间内搜索可能性最高的区域。2022/7/2952蒙特卡洛方法预测最佳搜救区域影响漂流的因素 漂浮物漂流模型:漂浮物漂流位置受事发位置、总流压速度、风压、搜索目标自身风压特性等因素的影响。事发位置可由报告得知,包含一定半径误差,总流速可以通过海上流速测量设备获得,不同海区海流速度随时间而变化,包含一定的测量误差。波浪的影响,对于搜救对象外形大小远小于波浪波长的搜救对象可以忽略不计。而在一些特殊区域,如沿岸潮流、岛屿、潮汐、内河水流等因素则需要特殊考虑。2022/7/2