蒙特卡羅方法教材課件

1、數(shù)學(xué)建模專題之蒙特卡羅方法主講：李培巒時(shí)間：2014-08-032022/7/291內(nèi)容提綱1.引言2.Monte Carlo模擬基本思想3.隨機(jī)數(shù)生成函數(shù)4.應(yīng)用實(shí)例舉例5.隨機(jī)游走模擬6.機(jī)械零件的可靠度計(jì)算7.排隊(duì)論模擬8.Monte Carlo方法求解規(guī)劃問題9.Monte Carlo方法預(yù)測(cè)搜救區(qū)域2022/7/292引言(Introduction)Monte Carlo方法：蒙特卡羅方法，又稱隨機(jī)模擬方法，屬于計(jì)算數(shù)學(xué)的一個(gè)分支，它是在上世紀(jì)四十年代中期為了適應(yīng)當(dāng)時(shí)原子能事業(yè)的發(fā)展而發(fā)展起來的。亦稱統(tǒng)計(jì)模擬方法，statistical simulation method 利用隨機(jī)數(shù)

2、進(jìn)行數(shù)值模擬的方法Monte Carlo名字的由來：Monte Carlo是摩納哥（monaco)最大的城市，該城以賭博聞名。Monte-Carlo, Monaco2022/7/293引言(Introduction)Monte Carlo方法的起源：二十世紀(jì)四十年代中期，由于科學(xué)技術(shù)的發(fā)展和電子計(jì)算機(jī)的發(fā)明，蒙特卡羅方法作為一種獨(dú)立的方法被提出來，并首先在核武器的試驗(yàn)與研制中得到了應(yīng)用(中子的鏈鎖反應(yīng))。但其基本思想并非新穎，可追溯到18世紀(jì)后半葉的蒲豐（Buffon）隨機(jī)投針實(shí)驗(yàn)（1777年）。Nicholas Metropolis (1915-1999)John Von Neumann (

3、1903-1957) 2022/7/294引言(Introduction)Monte Carlo方法的應(yīng)用： 物理：核物理，熱力學(xué)與統(tǒng)計(jì)物理，粒子輸運(yùn)問題等 數(shù)學(xué)：多重積分、解微分方程、非線性方程組求解等 工程領(lǐng)域：真空技術(shù)，水力學(xué)，激光技術(shù)等 經(jīng)濟(jì)學(xué)領(lǐng)域：期權(quán)定價(jià)、項(xiàng)目管理、投資風(fēng)險(xiǎn)決策等 其他領(lǐng)域：化學(xué)、醫(yī)學(xué)，生物，生產(chǎn)管理、系統(tǒng)科學(xué)、公用事業(yè)等方面。隨著科學(xué)技術(shù)的發(fā)展，其應(yīng)用范圍將更加廣泛。2022/7/295Monte Carlo方法的基本思想Buffon投針實(shí)驗(yàn)1768年，法國數(shù)學(xué)家Comte de Buffon利用投針實(shí)驗(yàn)估計(jì)的值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. 蒲豐投針問題利用關(guān)系式：求出值 其中 為投計(jì)次數(shù)，n 為針與平行線相交次數(shù)。這就是古典概率論中著名的蒲豐氏問題。2022/7/298一些人進(jìn)行了實(shí)驗(yàn)，其結(jié)果列于下表 ：實(shí)驗(yàn)者年份投計(jì)次數(shù)的實(shí)驗(yàn)值沃爾弗(Wolf)185050003.1596斯密思(Smith)185532043.1553福克斯(Fox)189411203.1419拉查里尼(Lazzarini)1901340837/299基本思想 MS基本思想：為了求解數(shù)學(xué)、物理、工程技術(shù)或隨機(jī)服務(wù)系統(tǒng)等方面的問題，首先構(gòu)造一個(gè)模型（概率統(tǒng)計(jì)模型），使所求問題的解正好

6、是該模型的參數(shù)或特征量或有關(guān)量，然后通過模擬（統(tǒng)計(jì)試驗(yàn)），給出模型參數(shù)或特征量的估計(jì)值，最后得出所求問題的近似解。 MS特點(diǎn)： 1.方法新穎、應(yīng)用面廣、實(shí)用性強(qiáng)； 2.隨機(jī)模擬方面的算法簡單，但計(jì)算量大； 3.模擬結(jié)果具有隨機(jī)性，精度較低； 4.模擬結(jié)果的收斂過程服從概率規(guī)律。2022/7/2910例1在我方某前沿防守地域，敵人以一個(gè)炮排（含兩門火炮）為單位對(duì)我方進(jìn)行干擾和破壞為躲避我方打擊，敵方對(duì)其陣地進(jìn)行了偽裝并經(jīng)常變換射擊地點(diǎn) 經(jīng)過長期觀察發(fā)現(xiàn)，我方指揮所對(duì)敵方目標(biāo)的指示有50是準(zhǔn)確的，而我方火力單位，在指示正確時(shí)，有1/3的射擊效果能毀傷敵人一門火炮，有1/6的射擊效果能全部毀傷敵人火

7、炮 現(xiàn)在希望能用某種方式把我方將要對(duì)敵人實(shí)施的20次打擊結(jié)果顯現(xiàn)出來，確定有效射擊的比率及毀傷敵方火炮的平均值。分析：這是一個(gè)概率問題，可以通過理論計(jì)算得到相應(yīng)的概率和期望值.但這樣只能給出作戰(zhàn)行動(dòng)的最終靜態(tài)結(jié)果，而顯示不出作戰(zhàn)行動(dòng)的動(dòng)態(tài)過程. 為了能顯示我方20次射擊的過程，現(xiàn)采用模擬的方式。舉例2022/7/2911 需要模擬出以下兩件事： 2 當(dāng)指示正確時(shí)，我方火力單位的射擊結(jié)果情況1 觀察所對(duì)目標(biāo)的指示正確與否模擬試驗(yàn)有兩種結(jié)果，每一種結(jié)果出現(xiàn)的概率都是1/2 因此，可用投擲一枚硬幣的方式予以確定，當(dāng)硬幣出現(xiàn)正面時(shí)為指示正確，反之為不正確 模擬試驗(yàn)有三種結(jié)果：毀傷一門火炮的可能性為1/

8、3(即2/6)，毀傷兩門的可能性為1/6，沒能毀傷敵火炮的可能性為1/2(即3/6) 這時(shí)可用投擲骰子的方法來確定：如果出現(xiàn)的是、三個(gè)點(diǎn)：則認(rèn)為沒能擊中敵人；如果出現(xiàn)的是、點(diǎn)：則認(rèn)為毀傷敵人一門火炮；若出現(xiàn)的是點(diǎn)：則認(rèn)為毀傷敵人兩門火炮問題分析2022/7/2912i：要模擬的打擊次數(shù)； k1：沒擊中敵人火炮的射擊總數(shù)； k2：擊中敵人一門火炮的射擊總數(shù)；k3：擊中敵人兩門火炮的射擊總數(shù)；E：有效射擊比率； E1：20次射擊平均每次毀傷敵人的火炮數(shù)符號(hào)說明2022/7/2913模擬框圖初始化:i=0,k1=0,k2=0,k3=0i=i+1骰子點(diǎn)數(shù)?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模擬結(jié)果2022/7/2915從以上模擬結(jié)果可計(jì)算出：2022/7/2916理論計(jì)算2022/7/2917結(jié)果比較 雖然模擬結(jié)果與理論計(jì)算不完全一致，但它卻能更加真實(shí)地表達(dá)實(shí)際戰(zhàn)斗動(dòng)態(tài)過程要使結(jié)果接近理論計(jì)算值，必須加大模擬次數(shù)，這就要求使用計(jì)算機(jī)模擬了。 用蒙特卡洛方法進(jìn)行計(jì)算機(jī)模擬的步驟：1 設(shè)計(jì)一個(gè)邏輯框圖，即建立模擬模型2 根據(jù)流程圖編寫程序，模擬隨機(jī)現(xiàn)象可通過生成具有各種概率分布的隨機(jī)數(shù)來模擬隨機(jī)現(xiàn)象，進(jìn)行模擬試驗(yàn)3 分析模擬結(jié)果，給出

10、所求問題的近似解（求解）.2022/7/2918注：rand(n)=rand(n,n)Matlab 中的隨機(jī)數(shù)生成函數(shù)randperm(m) 生成一個(gè)由 1:m 組成的隨機(jī)排列randn(m,n) 生成一個(gè)滿足正態(tài) m n 的隨機(jī)矩陣rand(m,n) 生成一個(gè)滿足均勻分布的 m n 隨機(jī)矩陣，矩陣的每個(gè)元素都在 (0,1) 之間。perms(1:n) 生成由 1:n 組成的全排列，共 n! 個(gè)結(jié)果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 中的隨機(jī)數(shù)生成函數(shù)2022/7/2920fix(x) : 截尾取整，直接將小數(shù)部分舍去floor(x) : 不超過 x 的最大整數(shù)ceil(x) : 不小于 x 的最小整數(shù)round(x) : 四舍五入取整Matlab中的取整函數(shù)2022/7/2921 模擬隨機(jī)投擲均勻硬幣，驗(yàn)證國徽朝上與朝下的概率是否都是 1/2 n=10000; % 給定試驗(yàn)次數(shù)m=0; % m 表示試驗(yàn)成功（國徽朝上）的次數(shù)for i=1:n x=r

12、andperm(2)-1; % randperm(2）生成1和2的隨機(jī)排列 y=x(1); if y=0 % 0 表示國徽朝上，1 表示國徽朝下 m=m+1; endendfprintf(國徽朝上的頻率為：%fn,m/n);小實(shí)例一：投擲硬幣2022/7/2922 隨機(jī)投擲骰子，驗(yàn)證各點(diǎn)出現(xiàn)的概率是否為 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. % 輸出結(jié)果小實(shí)例二：投擲骰子2022/7/2923 用蒙特卡羅 ( Monte Carlo ) 投點(diǎn)法計(jì)算 的值 n=10000; a=2; m=0; % m表示落入圓內(nèi)的次數(shù)for i=1:n x=rand(1)*a/2; y=rand(1)*a/2; if ( x2 + y2 = (a/2)2 ) m=m+1; endendfprintf(計(jì)算出來的 pi 為：%fn,4*m/n);小實(shí)例三：蒙特卡羅投點(diǎn)法2022/7/2924小實(shí)例三：蒙特卡羅投點(diǎn)法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、的等距平行線的白紙上，隨機(jī)投擲一根長為 l ( l d ) 的均勻直針，求針與平行線相交的概率，并計(jì)算的 值小實(shí)例四：蒲豐投針實(shí)驗(yàn)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) % 思考：如何生產(chǎn)指定分布隨機(jī)數(shù)?2022/7/2930隨機(jī)游走模擬隨機(jī)游走是一種基于運(yùn)用0，1區(qū)間的均勻分布隨機(jī)數(shù)序列來進(jìn)行的計(jì)算。利用蒙特卡羅方法多次模擬醉漢行走，統(tǒng)計(jì)在 X 處占總次數(shù)的比值，即為概率 PN(x)。

16、2022/7/2931隨機(jī)游走模擬設(shè)定朝右走的概率P, 總步數(shù)N, 模擬次數(shù)M , X , j=1 , S=0;i=1; x=0;i=Nj=M產(chǎn)生隨機(jī)數(shù)確定行走方向計(jì)算所在位置xi=i+1j=j+1X=xS=S+1YNYNYN計(jì)算S/Mi 步數(shù)，j 模擬的次數(shù)，S為M次中成功次數(shù)x表示N步后實(shí)際位置, 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機(jī)械零件的可靠度計(jì)算可靠度計(jì)算主要研究機(jī)械零件在一定分布的隨機(jī)應(yīng)力 S（即零件的正常工作應(yīng)力）作用下的可靠程度。應(yīng)力和強(qiáng)度都是隨機(jī)變量，影響應(yīng)力和強(qiáng)度的各種因素也是隨機(jī)變量。因此，機(jī)械零件的可靠度計(jì)算過程中，利用蒙特卡羅方法來處理對(duì)隨機(jī)量的計(jì)算問題，具有特有的優(yōu)勢(shì)。問題： 設(shè)計(jì)一個(gè)拉桿，若受外力 P 均值為20000N，標(biāo)準(zhǔn)差 2000N，拉桿的半徑 D 均值為20mm，標(biāo)準(zhǔn)差 0.5mm，材料強(qiáng)度 S 均值412 MPa，標(biāo)準(zhǔn)差15.6 MPa，計(jì)算其可靠度。2022/7/2934機(jī)械零件的可靠度計(jì)算問題分析：2022/7/2935機(jī)械零件的可靠度計(jì)算理論推導(dǎo)：20

18、22/7/2936機(jī)械零件的可靠度計(jì)算理論推導(dǎo)：2022/7/2937機(jī)械零件的可靠度計(jì)算蒙特卡羅模擬：通過多次試驗(yàn)，統(tǒng)計(jì)可靠零件個(gè)數(shù)，求出 其占總零件數(shù)的比值，即為可靠度。設(shè)定總零件數(shù) N=1000，可靠零件數(shù) 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根據(jù)已知分布產(chǎn)生隨機(jī)數(shù)p, d, q;計(jì)算應(yīng)力s=(4*p)/(pi*d2)YNs=qSt=St+1,i=i+1;計(jì)算可靠度：St/N2022/7/2938機(jī)械零件的可靠度計(jì)算程序：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為正態(tài)分布隨機(jī)數(shù) 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、排隊(duì)問題隨機(jī)模擬排隊(duì)論主要研究隨機(jī)服務(wù)系統(tǒng)的工作過程。在排隊(duì)系統(tǒng)中，服務(wù)對(duì)象的到達(dá)時(shí)間和服務(wù)時(shí)間都是隨機(jī)的。排隊(duì)論通過對(duì)隨機(jī)服務(wù)現(xiàn)象的統(tǒng)計(jì)研究，找出反映這些隨機(jī)現(xiàn)象平均特性的規(guī)律指標(biāo)，如排隊(duì)的長度、等待的時(shí)間及服務(wù)利用率。2022/7/2940系統(tǒng)的假設(shè)：（1） 顧客源是無窮的；（2） 排隊(duì)的長度沒有限制；（3）到達(dá)系統(tǒng)的顧客按先后順序依次進(jìn)入服務(wù)， “先到先服務(wù)”。 在某商店有一個(gè)售貨員，顧客陸續(xù)來到，售貨員逐個(gè)地接待顧客當(dāng)?shù)絹淼念櫩洼^多時(shí)，一部分顧客便須排隊(duì)等待，被接待后的顧客便離開商店設(shè)： 1顧客到來間隔時(shí)間服從參數(shù)為 0.1 的指數(shù)分布 2對(duì)顧客的服務(wù)時(shí)間服從，15上的均勻分布 3排隊(duì)

21、按先到先服務(wù)規(guī)則，隊(duì)長無限制 假定一個(gè)工作日為8小時(shí)，時(shí)間以分鐘為單位。1模擬一個(gè)工作日內(nèi)完成服務(wù)的個(gè)數(shù)及顧客平均等待時(shí)間t2模擬100個(gè)工作日，求出平均每日完成服務(wù)的個(gè)數(shù)及每日顧客的平均等待時(shí)間。單服務(wù)員的排隊(duì)模型模擬2022/7/2941w：總等待時(shí)間；ci：第i個(gè)顧客的到達(dá)時(shí)刻；bi：第i個(gè)顧客開始服務(wù)時(shí)刻； ei：第i個(gè)顧客服務(wù)結(jié)束時(shí)刻；xi：第i-1個(gè)顧客與第i個(gè)顧客之間到達(dá)的間隔時(shí)間；yi ：對(duì)第i個(gè)顧客的服務(wù)時(shí)間。符號(hào)說明2022/7/2942c1b1c3c4c5c2e1b2e2b3e3b4e4b5ci = ci-1+ xi % 到達(dá)時(shí)刻、時(shí)間間隔ei = bi+yi % 結(jié)束服

22、務(wù)時(shí)刻、開始服務(wù)時(shí)刻、服務(wù)時(shí)間bi = max(ci, ei-1)t思路分析2022/7/2943初始化：令i=1，ei-1=0，w=0產(chǎn)生間隔時(shí)間隨機(jī)數(shù)xi 參數(shù)為0.1的指數(shù)分布ci =xi , bi=xi 產(chǎn)生服務(wù)時(shí)間隨機(jī)數(shù)yi4,15的均勻分布ei =bi + yi累計(jì)等待時(shí)間：w= w+bi-ci準(zhǔn)備下一次服務(wù)：i=i+1產(chǎn)生間隔時(shí)間隨機(jī)數(shù)xi 參數(shù)為0.1的指數(shù)分布ci =ci-1+ xi 確定開始服務(wù)時(shí)間：bi=max(ci, ei-1)bi480?YNi=i-1,t=w/i輸出結(jié)果：完成服務(wù)個(gè)數(shù)：m=i 平均等待時(shí)間：t停止1模擬一日To Matlab(simu1)2模擬100

23、日To Matlab(simu2)b,服務(wù)時(shí)刻；c,到達(dá)時(shí)刻；e,結(jié)束時(shí)刻；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時(shí)停止迭代2022/7/2948框 圖初始化:給定MAXK,MAXP;k=0,p=0,Q:大整數(shù)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中編程，共需三個(gè)文件: randlp.m ,mylp.m 和 lpconst.m . 主程序?yàn)閞andlp.m .% mylp.mfunction z = mylp(x) %目標(biāo)函數(shù)z=2*x(1)2+x(2)2-x(1)*x(2)-8*x(1)-3*x(2); %轉(zhuǎn)化為求最小值問題% 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; %試驗(yàn)點(diǎn)下界b=10; %試驗(yàn)點(diǎn)上界n=1000; %試驗(yàn)點(diǎn)個(gè)數(shù)r1=unifrnd(a,b,n,1); %a, b均勻分布隨機(jī)數(shù)矩陣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蒙特卡洛方法預(yù)測(cè)最佳搜救區(qū)域

26、問題介紹 根據(jù)搜救目標(biāo)最后已知位置和時(shí)間,依據(jù)目標(biāo)海上漂流風(fēng)壓特性,綜合海區(qū)不同時(shí)段風(fēng)速、流速等環(huán)境信息及其不確定性影響,預(yù)測(cè)目標(biāo)漂移后的搜索區(qū)域,隨時(shí)間推移計(jì)算目標(biāo)漂流軌跡,同時(shí)求出某一具體時(shí)刻目標(biāo)在某可能區(qū)域的概率。 準(zhǔn)確地預(yù)測(cè)漂浮物的漂浮區(qū)域?qū)τ诩霸绱_定遇險(xiǎn)對(duì)象搜救區(qū)域、成功實(shí)施援救具有決定意義。 準(zhǔn)確的搜索區(qū)域劃定包含兩個(gè)要求: 搜救區(qū)域要能以最大概率包含搜救對(duì)象,不至于遺漏搜救對(duì)象; 搜救區(qū)域確定盡可能細(xì)致, 盡可能小,使搜救力量可以集中在最短時(shí)間內(nèi)搜索可能性最高的區(qū)域。2022/7/2952蒙特卡洛方法預(yù)測(cè)最佳搜救區(qū)域影響漂流的因素 漂浮物漂流模型:漂浮物漂流位置受事發(fā)位置、總流壓速度、風(fēng)壓、搜索目標(biāo)自身風(fēng)壓特性等因素的影響。事發(fā)位置可由報(bào)告得知,包含一定半徑誤差,總流速可以通過海上流速測(cè)量設(shè)備獲得,不同海區(qū)海流速度隨時(shí)間而變化,包含一定的測(cè)量誤差。波浪的影響,對(duì)于搜救對(duì)象外形大小遠(yuǎn)小于波浪波長的搜救對(duì)象可以忽略不計(jì)。而在一些特殊區(qū)域,如沿岸潮流、島嶼、潮汐、內(nèi)河水流等因素則需要特殊考慮。2022/7/2