數(shù)學建模算法之蒙特卡羅方法——原理、編程及應(yīng)用

1、.數(shù)學建模算法之蒙特卡羅方法原理、編程及應(yīng)用一、前言1946年，美國拉斯阿莫斯國家實驗室的三位科學家John von Neumann,Stan Ulam和Nick Metropolis共同發(fā)明了蒙特卡羅方法。此算法被評為20世紀最偉大的十大算法之一。蒙特卡羅方法（Monte Carlo method），又稱隨機抽樣或統(tǒng)計模擬方法，是一種以概率統(tǒng)計理論為指導的一類非常重要的數(shù)值計算方法。此方法使用隨機數(shù)（或更常見的偽隨機數(shù)）來解決很多計算問題的方法。由于傳統(tǒng)的經(jīng)驗方法由于不能逼近真實的物理過程，很難得到滿意的結(jié)果，而蒙特卡羅方法由于能夠真實地模擬實際物理過程，故解決問題與實際非常符合，可以得到很

2、圓滿的結(jié)果。二、蒙特卡羅方法的基本原理以及思想1、蒲豐投針實驗其基本思想源于法國數(shù)學家蒲豐提出著名的蒲豐投針實驗，并以該方法求圓周率。為了求得圓周率值，在十九世紀后期，有很多人作了這樣的試驗：將長為2l的一根針任意投到地面上，用針與一組相間距離為2a（la）的平行線相交的頻率代替概率P，再利用準確的關(guān)系式：求出值。其中為投針次數(shù)，n為針與平行線相交次數(shù)。這就是古典概率論中著名的蒲豐氏問題。2、 射擊問題設(shè)r表示射擊運動員的彈著點到靶心的距離，(r)表示擊中r處相應(yīng)的得分數(shù)（環(huán)數(shù)），f(r)為該運動員的彈著點的分布密度函數(shù)，它反映運動員的射擊水平。則該運動員的射擊成績?yōu)?用概率語言來說，<

3、g>是隨機變量(r)的數(shù)學期望，即 當所求解問題是某種隨機事件出現(xiàn)的概率，或者是某個隨機變量的期望值時，通過某種“實驗”的方法，以這種事件出現(xiàn)的頻率估計這一隨機事件的概率，或者得到這個隨機變量的某些數(shù)字特征，并將其作為問題的解。有一個例子可以使你比較直觀地了解蒙特卡洛方法：假設(shè)我們要計算一個不規(guī)則圖形的面積，那么圖形的不規(guī)則程度和分析性計算（比如，積分）的復雜程度是成正比的。蒙特卡洛方法是怎么計算的呢？假想你有一袋豆子，把豆子均勻地朝這個圖形上撒，然后數(shù)這個圖形之中有多少顆豆子，這個豆子的數(shù)目就是圖形的面積。當你的豆子越小，撒的越多的時候，結(jié)果就越精確。在這里我們要假定豆子都在一個平面上

4、，相互之間沒有重疊。蒙特卡羅方法通過抓住事物運動的幾何數(shù)量和幾何特征，利用數(shù)學方法來加以模擬，即進行一種數(shù)字模擬實驗。它是以一個概率模型為基礎(chǔ)，按照這個模型所描繪的過程，通過模擬實驗的結(jié)果，作為問題的近似解。蒙特卡羅方法與一般計算方法有很大區(qū)別，一般計算方法對于解決多維或因素復雜的問題非常困難，而蒙特卡羅方法對于解決這方面的問題卻比較簡單。其特點如下： 直接追蹤粒子，物理思路清晰，易于理解。 采用隨機抽樣的方法，較真切的模擬粒子輸運的過程，反映了統(tǒng)計漲落的規(guī)律。不受系統(tǒng)多維、多因素等復雜性的限制，是解決復雜系統(tǒng)粒子輸運問題的好方法。三、蒙特卡羅方法的編程計算陰影部分面積。一

5、個古人要求一個圖形的面積，他把圖形畫在一塊方形布上，然后找來一袋豆子，然后將所有豆子灑在布上，落在圖形內(nèi)豆子的重量比上那塊布上所有豆子的重量再乘以布的面積就是他所要求的圖形的面積。兩種編程思路來計算這個面積：方法一：將整個坐標軸看成一個邊長為12的正方形，然后均勻的這個正方形分成N（N的大小取決于劃分的步長）個點，然后找出N個點中有多少個點是屬于陰影部分中，假設(shè)這個值為k，則陰影部分的面積為：k/N*122方法二：將整個坐標軸看成一個邊長為12的正方形，然后在（-6，6）中隨機出N（N越大越好，至少超過1000）個點，然后找出這N個點中有多少個點在陰影區(qū)域內(nèi)，假設(shè)這個值為k，則陰影部分的面積為

6、：k/N*122。然后重復這個過程100次，求出100次面積計算結(jié)果的均值，這個均值為陰影部分面積。對比分析：以上兩個方法都是利用蒙特卡羅方法計算陰影部分面積，只是在處理的細節(jié)有一點區(qū)別。前者是把豆子均勻分布在布上；后者則是隨機把豆子仍在布上。就計算結(jié)果的精度而言，前者取決點的分割是否夠密，即N是否夠大；后者不僅僅通過N來控制精度，因為隨機的因素會造成單次計算結(jié)果偏高和偏小，所以進行反復多次計算最后以均值來衡量陰影部分面積。附上MATLAB程序：方法一：clearx=-6:0.01:6;y=x;s=size(x);zs=s(1,2)2;k=0;for i=1:s(1,2) for j=1:s(

7、1,2) a1=(x(i)2)/9+(y(j)2)/36; a2=(x(i)2)/36+y(j)2; a3=(x(i)-2)2+(y(j)+1)2; if a1<1 if a2<1 if a3<9 k=k+1; end end end endendmj=(122)*k/zs;運行結(jié)果：mj =7.2150方法二：clearN=10000;n=100;for j=1:n k=0;for i=1:N a=12*rand(1,2)-6; x(i)=a(1,1); y(i)=a(1,2); a1=(x(i)2)/9+(y(i)2)/36; a2=(x(i)2)/36+y(i)2; a

8、3=(x(i)-2)2+(y(i)+1)2; if a1<1 if a2<1 if a3<9 k=k+1; end end endendm(j)=(122)*k/N;endmj=mean(m);運行結(jié)果：mj =7.2500四、蒙特卡羅方法在數(shù)學建模中的應(yīng)用舉例1.問題的提出某食品加工廠主要生產(chǎn)即食產(chǎn)品，一般當天生產(chǎn)的產(chǎn)品必須當天售出，否則就會出現(xiàn)不能保質(zhì)、或變質(zhì)、造成一定的經(jīng)濟損失，如果市場需求量大而生產(chǎn)量不足，則也會影響工廠的銷售收入，該產(chǎn)品的單位成本為1.5元，單位產(chǎn)品售價為4元。工廠為了避免產(chǎn)品滯銷存貨過多而造成的經(jīng)濟損失，提出了如何制定合理的生產(chǎn)與庫存數(shù)量的方案問題

9、，能夠使得工廠能有盡可能多的收益，經(jīng)初步考慮擬從以下兩種生產(chǎn)與庫存方案中選出一個較好的方案：方案(1)：按前一天的銷售量作為當天的生產(chǎn)庫存量。方案(2)：按前兩天的平均銷售量作為當天的生產(chǎn)庫存量。2.問題的分析與假設(shè)及參數(shù)定義實際中的問題，生產(chǎn)與庫存多了，銷售不出去會造成經(jīng)濟損失，生產(chǎn)與庫存少了不能滿足需求也會造成一定的損失，工廠需要依據(jù)實際不確定的需求量來制定合理的生產(chǎn)與庫存方案，使得能有盡量大的經(jīng)濟收益。解決問題的基本思路：利用蒙特卡羅方法隨機模擬市場對該產(chǎn)品需求量，統(tǒng)計計算出按照兩種不同方案T天后工廠的經(jīng)濟值，比較不同方案經(jīng)濟效益的大小，選出一個較好的方案。假設(shè)市場對該產(chǎn)品的每天需求數(shù)量

10、是一個隨機變量，從統(tǒng)計學的角度分析得知，該隨機變量服從正態(tài)分布。為了編程實現(xiàn)問題的目標，引入如下的變量：本文中參數(shù)符號定義：T表示模擬天數(shù)；C表示每天的需求量；KC1表示方案(1)當天的生產(chǎn)與庫存量；KC2表示方案(2)當天的生產(chǎn)與庫存量；S1表示方案(1)前一天的銷售量；S21表示方案(2)前一天的銷售量；S22表示方案(2)前二天的銷售量；ST1表示方案(1)當天的實際銷售量；ST2表示方案(2)當天的實際銷售量；L1表示方案(1)當天的實際利潤；L2表示方案(2)當天的實際利潤；LS1表示方案(1)實際累計總利潤；LS2表示方案(2)實際累計總利潤。3.模型的建立與求解根據(jù)上面的分析，利

11、用蒙特卡羅方法編程實現(xiàn)，主要隨機模擬前一天和前兩天的各種不同的銷售量，來確定當天的生產(chǎn)與庫存量，依據(jù)可能的實際銷售量，計算出當天的銷售利潤，選擇使連續(xù)幾天利潤盡可能大的方案，下面給MATLAB程序。(1)建立蒙特卡羅方法的M文件，函數(shù)名：mcun.mFunction LS1,LS2=mcun(T,S1,S21,S22)LS1=0;LS2=0;k=1；While k<TKC1=S1;KC2=(S21+S22)/2；C=normrnd(1500,30*30)if C<KC1ST1=KC1；ElseST1=C；endif C<KC2ST2=KC2；elseST2=C；endL1=4

12、*ST11.5*KC1；L2=4*ST21.5*KC2；LS1=LS1+L1；LS2=LS2+L2；k=k+1；endS1=ST1:S22=S21;S21=ST2；(2)調(diào)用函數(shù)mcun(T,S1,S21,S22)4.模型的結(jié)果分析與推廣在MATLAB命令窗口多次輸入?yún)?shù)T,S1,S21和S22數(shù)值，分別調(diào)用函數(shù)meun(T,S1,S21,S22)，進行求解計算，并對結(jié)果進行分析，由若干次模擬實演的結(jié)果可以看出，方案(1)的利潤總是小于方案(2)的利潤，所以該工廠實際按方案(2)進行組織生產(chǎn)與庫存，工廠會有更好的經(jīng)濟效益。五、參考文獻1蒙特卡羅方法及應(yīng)用J. 尹增謙,管景峰,張曉宏,曹春梅. 物理與工程. 2002(03).45-49.2蒙特卡羅方法應(yīng)用研究J. 王巖,尹海麗,竇在祥. 青島理工大學學報. 2006(02).111-113.3高職數(shù)學中蒙特卡羅方法的應(yīng)用J. 陳楊林,羅婷. 九江職業(yè)技術(shù)學院學報. 2013(02).22-23.4蒙特卡羅方法在教學中的應(yīng)用J. 康件麗,黃俊杰. 電腦學習. 2008(02).46-47.5蒙特卡羅方法與問題的維數(shù)J. 裴鹿成. &#