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

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

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

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

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

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

6、：k/N*122。然后重復(fù)這個(gè)過程100次，求出100次面積計(jì)算結(jié)果的均值，這個(gè)均值為陰影部分面積。對比分析：以上兩個(gè)方法都是利用蒙特卡羅方法計(jì)算陰影部分面積，只是在處理的細(xì)節(jié)有一點(diǎn)區(qū)別。前者是把豆子均勻分布在布上；后者則是隨機(jī)把豆子仍在布上。就計(jì)算結(jié)果的精度而言，前者取決點(diǎn)的分割是否夠密，即N是否夠大；后者不僅僅通過N來控制精度，因?yàn)殡S機(jī)的因素會(huì)造成單次計(jì)算結(jié)果偏高和偏小，所以進(jìn)行反復(fù)多次計(jì)算最后以均值來衡量陰影部分面積。附上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;運(yùn)行結(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);運(yùn)行結(jié)果：mj =7.2500四、蒙特卡羅方法在數(shù)學(xué)建模中的應(yīng)用舉例1.問題的提出某食品加工廠主要生產(chǎn)即食產(chǎn)品，一般當(dāng)天生產(chǎn)的產(chǎn)品必須當(dāng)天售出，否則就會(huì)出現(xiàn)不能保質(zhì)、或變質(zhì)、造成一定的經(jīng)濟(jì)損失，如果市場需求量大而生產(chǎn)量不足，則也會(huì)影響工廠的銷售收入，該產(chǎn)品的單位成本為1.5元，單位產(chǎn)品售價(jià)為4元。工廠為了避免產(chǎn)品滯銷存貨過多而造成的經(jīng)濟(jì)損失，提出了如何制定合理的生產(chǎn)與庫存數(shù)量的方案問題

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

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

11、用蒙特卡羅方法編程實(shí)現(xiàn)，主要隨機(jī)模擬前一天和前兩天的各種不同的銷售量，來確定當(dāng)天的生產(chǎn)與庫存量，依據(jù)可能的實(shí)際銷售量，計(jì)算出當(dāng)天的銷售利潤，選擇使連續(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)，進(jìn)行求解計(jì)算，并對結(jié)果進(jìn)行分析，由若干次模擬實(shí)演的結(jié)果可以看出，方案(1)的利潤總是小于方案(2)的利潤，所以該工廠實(shí)際按方案(2)進(jìn)行組織生產(chǎn)與庫存，工廠會(huì)有更好的經(jīng)濟(jì)效益。五、參考文獻(xiàn)1蒙特卡羅方法及應(yīng)用J. 尹增謙,管景峰,張曉宏,曹春梅. 物理與工程. 2002(03).45-49.2蒙特卡羅方法應(yīng)用研究J. 王巖,尹海麗,竇在祥. 青島理工大學(xué)學(xué)報(bào). 2006(02).111-113.3高職數(shù)學(xué)中蒙特卡羅方法的應(yīng)用J. 陳楊林,羅婷. 九江職業(yè)技術(shù)學(xué)院學(xué)報(bào). 2013(02).22-23.4蒙特卡羅方法在教學(xué)中的應(yīng)用J. 康件麗,黃俊杰. 電腦學(xué)習(xí). 2008(02).46-47.5蒙特卡羅方法與問題的維數(shù)J. 裴鹿成. &#