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

數(shù)學(xué)建模算法之蒙特卡羅方法 原理 編程及應(yīng)用 一 前言一 前言 1946 年 美國(guó)拉斯阿莫斯國(guó)家實(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ù) 或更常見(jiàn)的偽隨機(jī)數(shù) 來(lái)解決很多計(jì)算問(wèn)題的方法 由于傳統(tǒng)的經(jīng)驗(yàn)方法由于不能逼近真實(shí)的物理過(guò)程 很難得到滿意的結(jié)果 而蒙特卡羅方法由于能夠真實(shí)地模擬實(shí)際物理過(guò)程 故解決問(wèn)題與實(shí)際非常符合 可以得到很圓滿的結(jié)果 二 二 蒙特卡羅方法的蒙特卡羅方法的基本原理以及思想基本原理以及思想 1 蒲豐投針實(shí)驗(yàn) 其基本思想源于法國(guó)數(shù)學(xué)家蒲豐提出著名的蒲豐投針實(shí)驗(yàn) 并以該方法求圓 周率 為了求得圓周率 值 在十九世紀(jì)后期 有很多人作了這樣的試驗(yàn) 將長(zhǎng)為 2l 的一根針任意投到地面上 用針與一組相間距離為 2a l a 的平行線相交 的頻率代替概率 P 再利用準(zhǔn)確的關(guān)系式 a l P 2 22 n N a l aP l 求出 值 其中 為投針次數(shù) n 為針與平行線相交次數(shù) 這就是古典概率論中 著名的蒲豐氏問(wèn)題 2 射擊問(wèn)題 設(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)員的射擊成績(jī)?yōu)?0 gg r f r dr 用概率語(yǔ)言來(lái)說(shuō) 是隨機(jī)變量 r 的數(shù)學(xué)期望 即 rgEg 當(dāng)所求解問(wèn)題是某種隨機(jī)事件出現(xiàn)的概率 或者是某個(gè)隨機(jī)變量的期望值時(shí) 通過(guò)某種 實(shí)驗(yàn) 的方法 以這種事件出現(xiàn)的頻率估計(jì)這一隨機(jī)事件的概率 或 者得到這個(gè)隨機(jī)變量的某些數(shù)字特征 并將其作為問(wèn)題的解 有一個(gè)例子可以使 你比較直觀地了解蒙特卡洛方法 假設(shè)我們要計(jì)算一個(gè)不規(guī)則圖形的面積 那么圖形的不規(guī)則程度和分析性計(jì) 算 比如 積分 的復(fù)雜程度是成正比的 蒙特卡洛方法是怎么計(jì)算的呢 假想 你有一袋豆子 把豆子均勻地朝這個(gè)圖形上撒 然后數(shù)這個(gè)圖形之中有多少顆豆 子 這個(gè)豆子的數(shù)目就是圖形的面積 當(dāng)你的豆子越小 撒的越多的時(shí)候 結(jié)果 就越精確 在這里我們要假定豆子都在一個(gè)平面上 相互之間沒(méi)有重疊 蒙特卡羅方法通過(guò)抓住事物運(yùn)動(dòng)的幾何數(shù)量和幾何特征 利用數(shù)學(xué)方法來(lái)加 以模擬 即進(jìn)行一種數(shù)字模擬實(shí)驗(yàn) 它是以一個(gè)概率模型為基礎(chǔ) 按照這個(gè)模型 所描繪的過(guò)程 通過(guò)模擬實(shí)驗(yàn)的結(jié)果 作為問(wèn)題的近似解 蒙特卡羅方法與一般計(jì)算方法有很大區(qū)別 一般計(jì)算方法對(duì)于解決多維或因 素復(fù)雜的問(wèn)題非常困難 而蒙特卡羅方法對(duì)于解決這方面的問(wèn)題卻比較簡(jiǎn)單 其 特點(diǎn)如下 直接追蹤粒子 物理思路清晰 易于理解 采用隨機(jī)抽樣的方法 較真切的模擬粒子輸運(yùn)的過(guò)程 反映了統(tǒng)計(jì)漲落的規(guī)律 不受系統(tǒng)多維 多因素等復(fù)雜性的限制 是解決復(fù)雜系統(tǒng)粒子輸運(yùn)問(wèn)題的好方 法 三 蒙特卡羅方法的編程三 蒙特卡羅方法的編程 計(jì)算陰影部分面積 一個(gè)古人要求一個(gè)圖形的面積 他把圖形畫(huà)在一塊方形布上 然后找來(lái)一袋 豆子 然后將所有豆子灑在布上 落在圖形內(nèi)豆子的重量比上那塊布上所有豆子 的重量再乘以布的面積就是他所要求的圖形的面積 兩種編程思路來(lái)計(jì)算這個(gè)面積 方法一 將整個(gè)坐標(biāo)軸看成一個(gè)邊長(zhǎng)為 12 的正方形 然后均勻的這個(gè)正方 形分成 N N 的大小取決于劃分的步長(zhǎng) 個(gè)點(diǎn) 然后找出 N 個(gè)點(diǎn)中有多少個(gè)點(diǎn)是 屬于陰影部分中 假設(shè)這個(gè)值為 k 則陰影部分的面積為 k N 12 2 方法二 將整個(gè)坐標(biāo)軸看成一個(gè)邊長(zhǎng)為 12 的正方形 然后在 6 6 中隨 機(jī)出 N N 越大越好 至少超過(guò) 1000 個(gè)點(diǎn) 然后找出這 N 個(gè)點(diǎn)中有多少個(gè)點(diǎn)在 陰影區(qū)域內(nèi) 假設(shè)這個(gè)值為 k 則陰影部分的面積為 k N 12 2 然后重復(fù)這個(gè) 過(guò)程 100 次 求出 100 次面積計(jì)算結(jié)果的均值 這個(gè)均值為陰影部分面積 對(duì)比分析 以上兩個(gè)方法都是利用蒙特卡羅方法計(jì)算陰影部分面積 只是在 處理的細(xì)節(jié)有一點(diǎn)區(qū)別 前者是把豆子均勻分布在布上 后者則是隨機(jī)把豆子仍 在布上 就計(jì)算結(jié)果的精度而言 前者取決點(diǎn)的分割是否夠密 即 N 是否夠大 后者不僅僅通過(guò) N 來(lái)控制精度 因?yàn)殡S機(jī)的因素會(huì)造成單次計(jì)算結(jié)果偏高和偏小 所以進(jìn)行反復(fù)多次計(jì)算最后以均值來(lái)衡量陰影部分面積 附上 MATLAB 程序 方法一 clear x 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 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 end end mj 12 2 k zs 運(yùn)行結(jié)果 mj 7 2150 方法二 clear N 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 a3 x i 2 2 y i 1 2 if a1 1 if a2 1 if a3 9 k k 1 end end end end m j 12 2 k N end mj mean m 運(yùn)行結(jié)果 mj 7 2500 四 蒙特卡羅方法在數(shù)學(xué)建模中的應(yīng)用舉例四 蒙特卡羅方法在數(shù)學(xué)建模中的應(yīng)用舉例 1 問(wèn)題的提出 某食品加工廠主要生產(chǎn)即食產(chǎn)品 一般當(dāng)天生產(chǎn)的產(chǎn)品必須當(dāng)天售出 否則 就會(huì)出現(xiàn)不能保質(zhì) 或變質(zhì) 造成一定的經(jīng)濟(jì)損失 如果市場(chǎng)需求量大而生產(chǎn)量 不足 則也會(huì)影響工廠的銷售收入 該產(chǎn)品的單位成本為 1 5 元 單位產(chǎn)品售價(jià) 為 4 元 工廠為了避免產(chǎn)品滯銷存貨過(guò)多而造成的經(jīng)濟(jì)損失 提出了如何制定合 理的生產(chǎn)與庫(kù)存數(shù)量的方案問(wèn)題 能夠使得工廠能有盡可能多的收益 經(jīng)初步考 慮擬從以下兩種生產(chǎn)與庫(kù)存方案中選出一個(gè)較好的方案 方案 1 按前一天的銷售量作為當(dāng)天的生產(chǎn)庫(kù)存量 方案 2 按前兩天的平均銷售量作為當(dāng)天的生產(chǎn)庫(kù)存量 2 問(wèn)題的分析與假設(shè)及參數(shù)定義 實(shí)際中的問(wèn)題 生產(chǎn)與庫(kù)存多了 銷售不出去會(huì)造成經(jīng)濟(jì)損失 生產(chǎn)與庫(kù)存 少了不能滿足需求也會(huì)造成一定的損失 工廠需要依據(jù)實(shí)際不確定的需求量來(lái)制 定合理的生產(chǎn)與庫(kù)存方案 使得能有盡量大的經(jīng)濟(jì)收益 解決問(wèn)題的基本思路 利用蒙特卡羅方法隨機(jī)模擬市場(chǎng)對(duì)該產(chǎn)品需求量 統(tǒng)計(jì)計(jì)算出按照兩種不同方 案 T 天后工廠的經(jīng)濟(jì)值 比較不同方案經(jīng)濟(jì)效益的大小 選出一個(gè)較好的方案 假設(shè)市場(chǎng)對(duì)該產(chǎn)品的每天需求數(shù)量是一個(gè)隨機(jī)變量 從統(tǒng)計(jì)學(xué)的角度分析得 知 該隨機(jī)變量服從正態(tài)分布 為了編程實(shí)現(xiàn)問(wèn)題的目標(biāo) 引入如下的變量 本文中參數(shù)符號(hào)定義 T 表示 模擬天數(shù) C 表示每天的需求量 KC1 表示方案 1 當(dāng)天的生產(chǎn)與庫(kù)存量 KC2 表 示方案 2 當(dāng)天的生產(chǎn)與庫(kù)存量 S1 表示方案 1 前一天的銷售量 S21 表示方案 2 前一天的銷售量 S22 表示方案 2 前二天的銷售量 ST1 表示方案 1 當(dāng)天的 實(shí)際銷售量 ST2 表示方案 2 當(dāng)天的實(shí)際銷售量 L1 表示方案 1 當(dāng)天的實(shí)際利 潤(rùn) L2 表示方案 2 當(dāng)天的實(shí)際利潤(rùn) LS1 表示方案 1 實(shí)際累計(jì)總利潤(rùn) LS2 表 示方案 2 實(shí)際累計(jì)總利潤(rùn) 3 模型的建立與求解 根據(jù)上面的分析 利用蒙特卡羅方法編程實(shí)現(xiàn) 主要隨機(jī)模擬前一天和前兩 天的各種不同的銷售量 來(lái)確定當(dāng)天的生產(chǎn)與庫(kù)存量 依據(jù)可能的實(shí)際銷售量 計(jì)算出當(dāng)天的銷售利潤(rùn) 選擇使連續(xù)幾天利潤(rùn)盡可能大的方案 下面給 MATLAB 程序 1 建立蒙特卡羅方法的 M 文件 函數(shù)名 mcun m Function LS1 LS2 mcun T S1 S21 S22 LS1 0 LS2 0 k 1 While k T KC1 S1 KC2 S21 S22 2 C normrnd 1500 30 30 ifC KC1 ST1 KC1 Else ST1 C end if C KC2 ST2 KC2 else ST2 C end L1 4 ST1 1 5 KC1 L2 4 ST2 1 5 KC2 LS1 LS1 L1 LS2 LS2 L2 k k 1 end S1 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ì)算 并對(duì)結(jié)果進(jìn)行分析 由若干次模擬實(shí)演的 結(jié)果可以看出 方案 1 的利潤(rùn)總是小于方案 2 的利潤(rùn) 所以該工廠實(shí)際按方案 2 進(jìn)行組織生產(chǎn)與庫(kù)存 工廠會(huì)有更好的經(jīng)濟(jì)效益 五 參考文獻(xiàn)五 參考文獻(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 蒙特卡羅方法與問(wèn)題的維數(shù)