數(shù)學(xué)建模-計(jì)算機(jī)模擬技術(shù)

1模擬基礎(chǔ)

2MonteCaHo模擬

3模擬模型案例

1模擬基礎(chǔ)

1.0模擬的背景、思路

應(yīng)用領(lǐng)域：

■第二次世界大戰(zhàn)期間，J.V.Neumann等人將進(jìn)行的“中

子擴(kuò)散”的科研項(xiàng)目取名為"Monte-Carlo"

?運(yùn)輸系統(tǒng)模擬

■摩天大樓安全疏散系統(tǒng)模擬

?國(guó)民經(jīng)濟(jì)發(fā)展模擬

?人口增長(zhǎng)系統(tǒng)模擬

?供水系統(tǒng)模擬

?管理系統(tǒng)模擬

?雷達(dá)系統(tǒng)模擬

?戰(zhàn)爭(zhēng)系統(tǒng)模擬

模擬思路：

?“模擬”一對(duì)系統(tǒng)抽象建模

?“試驗(yàn)”一根據(jù)模型設(shè)計(jì)算法，編程進(jìn)行反復(fù)試驗(yàn)

?“估計(jì)”一根據(jù)試驗(yàn)數(shù)據(jù)

?“收集”一根據(jù)試驗(yàn)結(jié)果作出判斷

1.1模擬的基本知識(shí)

1.1.1模擬的概念及作用

?現(xiàn)實(shí)系統(tǒng)的數(shù)學(xué)或邏輯模型可能十分復(fù)雜,

例如大多數(shù)具有隨機(jī)因素的復(fù)雜系統(tǒng)，其

中的一些隨機(jī)性因素很難用準(zhǔn)確的數(shù)學(xué)公

式表述，從而也無(wú)法對(duì)整個(gè)系統(tǒng)采用解析

法求解。模擬是處理這類實(shí)際問(wèn)題的有力

工具。

?模擬通常借助于計(jì)算機(jī)進(jìn)行。

計(jì)算機(jī)模擬：在已經(jīng)建立的數(shù)學(xué)、邏輯模

型的基礎(chǔ)之上，通過(guò)計(jì)算機(jī)試驗(yàn)，對(duì)一個(gè)

系統(tǒng)按照一定的決策原則或作業(yè)規(guī)則，由

一個(gè)狀態(tài)變換為另一個(gè)狀態(tài)的行為進(jìn)行描

述和分析。

模擬的作用：

?對(duì)于很難用解析方法加以處理的問(wèn)題，模

擬是一種有效的技術(shù)；

?對(duì)建模過(guò)程中的假設(shè)進(jìn)行鑒定，對(duì)理論研

究的結(jié)論加以檢驗(yàn)；

?對(duì)不同的實(shí)現(xiàn)方案進(jìn)行多次模擬，按照既

定的目標(biāo)函數(shù)對(duì)不同方案進(jìn)行比較，從中

選擇最優(yōu)方案。

1.1.2模擬的分類

通常，模擬時(shí)間是模擬的主要自變量。

設(shè)計(jì)正確的模擬時(shí)間推進(jìn)機(jī)理：模擬過(guò)程中

應(yīng)根據(jù)系統(tǒng)的特性正確推進(jìn)模擬時(shí)間，使

系統(tǒng)中各要素與發(fā)生的事件保持同步。

推進(jìn)模擬時(shí)間的基本方法：

?下次事件法：將模擬時(shí)間由一個(gè)事件發(fā)生

的時(shí)間點(diǎn)推進(jìn)到緊接著的下一次事件發(fā)生

的時(shí)間點(diǎn)。

?固定時(shí)間步長(zhǎng)法：模擬時(shí)間每次均以相等

的固定步長(zhǎng)向前推進(jìn)，每到達(dá)一個(gè)新的模

擬時(shí)間點(diǎn)需檢查相應(yīng)時(shí)間段內(nèi)是否發(fā)生了

事件。需根據(jù)實(shí)際問(wèn)題合理設(shè)置模擬時(shí)間

發(fā)生改變的步長(zhǎng)。

根據(jù)模擬過(guò)程中因變量的變化情況進(jìn)行分類:

1）離散型模擬：因變量在與事件時(shí)間有關(guān)的具體模

擬時(shí)間點(diǎn)呈離散性變化。大多數(shù)系統(tǒng)（如排隊(duì)服

務(wù)系統(tǒng)）可采用離散型模擬。

時(shí)間推進(jìn)方法：一般采用下次事件法

應(yīng)當(dāng)重點(diǎn)對(duì)系統(tǒng)狀態(tài)可能發(fā)生改變的事件進(jìn)行描述,

并確定這些事件之間的邏輯關(guān)系。

因

變

量

A

事件時(shí)間時(shí)間

圖1離散型模擬

采用離散型模擬模型其中發(fā)生系

邦r隊(duì)系統(tǒng)通常

F。，

T客達(dá)二

統(tǒng)狀態(tài)變化的事件兩個(gè)一是有顧到；

有.

是服務(wù)員完成服務(wù)將最近?發(fā)生述兩種事件之

。上

一的時(shí)刻設(shè)置為下事件發(fā)生點(diǎn)就可將服務(wù)過(guò)

次，

程描述為圖幼扃而疑就講財(cái)比

2）連續(xù)型模擬：因變量隨時(shí)間的改變呈連續(xù)性變

化。在大多數(shù)計(jì)算機(jī)模擬過(guò)程中，按固定的步長(zhǎng)

推進(jìn)模擬時(shí)間。

通常需建立一系列的由系統(tǒng)狀態(tài)變量組成的

狀態(tài)方程組，以描述狀態(tài)變量與模擬時(shí)間的關(guān)系。

3）混合型模擬：因變量隨時(shí)間的推移而作連續(xù)性

的變化并具有離散性的突變，如庫(kù)存控制系統(tǒng)。

1.1.3模擬的方式

終態(tài)模擬：在規(guī)定的時(shí)間T內(nèi)進(jìn)行模擬運(yùn)行，時(shí)間

達(dá)到T時(shí)，模擬終止。其性能指標(biāo)明顯取決于系

統(tǒng)的初始狀態(tài)。

穩(wěn)態(tài)模擬：隨著模擬時(shí)間的推移，系統(tǒng)的性能逐漸

趨于平穩(wěn)。其目的是研究非終態(tài)系統(tǒng)長(zhǎng)期運(yùn)行條

件下的穩(wěn)態(tài)性能，模擬時(shí)間的長(zhǎng)短取決于能否獲

得系統(tǒng)性能的優(yōu)良估計(jì)（可由模擬輸出的精度確

定）。

1.1.4模擬的一般步驟

?明確問(wèn)題，建立模型。正確描述待研究問(wèn)題，明

確規(guī)定模擬的目的和任務(wù)，確定衡量系統(tǒng)性能或

模擬輸出結(jié)果的目標(biāo)函數(shù)，然后根據(jù)系統(tǒng)的結(jié)構(gòu)

及作業(yè)規(guī)則，分析系統(tǒng)各狀態(tài)變量之間的關(guān)系，

以此為基礎(chǔ)建立所研究的系統(tǒng)模型。

?收集和整理數(shù)據(jù)資料。模擬技術(shù)的正確運(yùn)用，往

往由大量的輸入數(shù)據(jù)作依靠。在隨機(jī)模擬中，應(yīng)

認(rèn)真分析具體收集到的隨機(jī)性數(shù)據(jù)資料，確定系

統(tǒng)中隨機(jī)性因素的概率分布特性，以此為依據(jù)產(chǎn)

生模擬過(guò)程所必需的抽樣數(shù)據(jù)。

?編制程序，模擬運(yùn)行。

?分析模擬輸出結(jié)果：模擬結(jié)果的統(tǒng)計(jì)特性（樣本

均值、方差、置信區(qū)間等），靈敏性分析，選擇

最優(yōu)方案。

注：模擬結(jié)果的統(tǒng)計(jì)分析模擬的輸出結(jié)果是分布特

征未知的隨機(jī)變量，每次運(yùn)行的結(jié)果僅僅是對(duì)該

隨機(jī)變量所有觀察值總體的一次抽樣，對(duì)總體的

代表性很差，雖然可以增加模擬運(yùn)行的時(shí)間從而

增加抽樣次數(shù)，但這些數(shù)據(jù)總是由一個(gè)“種子”

經(jīng)過(guò)一定的算法而獲得的偽隨機(jī)序列，它們是自

相關(guān)的，并不能構(gòu)成統(tǒng)計(jì)上獨(dú)立的隨機(jī)樣本。

1.2隨機(jī)模擬案例：趕上火車的概率

【問(wèn)題】如圖，一列火車從A站開(kāi)往B站，某人每天

趕往B站上這趟火車.

工/某人

火車運(yùn)行方向/

A°--------------------°B

思考：請(qǐng)研究他能否趕上這趟火車。

他已了解到：

1）火車從A站到B站的運(yùn)行時(shí)間是均值為30分

鐘，標(biāo)準(zhǔn)差為2分鐘的隨機(jī)變量；

2）火車在下午大約1點(diǎn)離開(kāi)A站，離開(kāi)時(shí)刻

的頻率分布如下：

出發(fā)時(shí)刻午后1:00午后1:05午后1:10

頻率0.70.20.1

他到達(dá)B站的時(shí)刻的頻率分布為

時(shí)刻午后1:28午后1:30午后1:32午后1:34

頻率0.30.40.20.1

他能否及時(shí)趕上火車?

明確問(wèn)題：他能及時(shí)趕上火車的概率是多少?

建模方向（思路）：

i）分析法：用概率統(tǒng)計(jì)知識(shí)建立分析模型，求解

析解。（思考）

ii）模擬法：用概率統(tǒng)計(jì)知識(shí)建立模型，通過(guò)模擬

求近似解。

即先建立模擬模型，然后通過(guò)計(jì)算機(jī)模擬得到

問(wèn)題的近似解。在同樣條件下多次試驗(yàn)，計(jì)算

他能及時(shí)趕上火車的頻率。

問(wèn)題分析：能及時(shí)趕上火車的充要條件是:

TH

即心―(<T2

其中

彳；一火車從A站出發(fā)的時(shí)刻;

T2一火車的運(yùn)行時(shí)間；

走什么變

丁3一他到達(dá)B站的時(shí)刻。重？如何

模擬?

基本假設(shè):

i)假設(shè)T1，T,T?都是相互獨(dú)立的隨機(jī)變量;

JLJ24T

ii)將午后1時(shí)記為t=0,設(shè)火車運(yùn)行時(shí)間T?

服從正態(tài)分布：T2?N(30,22)o

建立模型：為了簡(jiǎn)化計(jì)算，將下午L點(diǎn)、記為初始時(shí)

刻。得到隨機(jī)變量「和人的分布律如下：

火車出發(fā)時(shí)刻！\和人到達(dá)B站時(shí)刻T3的分布

律分別為：

T1（分）0510

P(t)0.70.20.1

丁(分)

T3（分）2288303033223434

~~P⑴(tj00^.33004.400^.22O0H.1

能及時(shí)趕上火車的概率〃=P{T3VTi+T2}

如果r為在(0J)均勻分布的隨機(jī)數(shù)，為了模擬隨

機(jī)變量%和T3,可以通過(guò)如下方法：

0,0<r<0.7,

4=<5,0.7<r<0.95

10,0.9<r<1.貝和，3可分另U

用來(lái)模擬隨機(jī)變

〃28,0<r<0,3

9量％和T3。

30,0.3<r<0,7,

/3=<

32，0.7<r<0,9?

34,0.9<r<1.0.

模擬算法設(shè)計(jì)

輸入：趕火車次數(shù)（天數(shù)）

輸出：趕上火車的頻率

兩種不同風(fēng)格的算法描述

主要變量說(shuō)明：

n模擬次數(shù)

k臨時(shí)變量，存儲(chǔ)當(dāng)前累計(jì)模擬次數(shù)

count存儲(chǔ)趕上火車的次數(shù)

算法1（分步驟描述）

第1步輸入模擬次數(shù)n

第2步k=l,count=0

第3步當(dāng)k〈=n,執(zhí)行Step4-9,否則執(zhí)行第11步

第4步生成均勻分布隨機(jī)數(shù)賦給r

第5步由r及公式確定T1模擬火車出發(fā)時(shí)刻

爭(zhēng)6步生成均勻分布隨機(jī)數(shù)賦給r;.

第7步由r及公式確定T3模擬人達(dá)到時(shí)刻

第8步生成正態(tài)分布隨機(jī)數(shù)T2模擬火車運(yùn)行時(shí)間

第9步IFT1+T2>T3,count=count+l,END

第10步執(zhí)行第3步

第H步輸出趕上火車頻率p=count/n

算法2（偽代碼描述）

i）初始化：

輸入模擬次數(shù)n;

count=0;

道）模擬11次

fori=lton,

模擬隨機(jī)變量￡,丁2，丁3,

分別賦給tl,t2,t3;

iftl+t2>t3,

count=count+l

endif

endfor

app_prb=count/n;

n=input「輸入模擬次數(shù)：1）;

count=0;

fori=l:n,

模擬rtl=rand;當(dāng)模擬隨機(jī)變量tl（火車從A站出發(fā)的時(shí)刻）

程序ifrtl<0.7

Tl=0;

elseifrtl>=0.7&rtl<0.9

Tl=5;

else?exid

T2=30+randn*2;%模取隨機(jī)菱量t2（火車的運(yùn)行時(shí)間）

%模擬隨機(jī)變量t3（他到達(dá)B站的時(shí)亥lj）

rt3=rand;

ifrt3<0.3

T3=28;

elseifrt3>=0.3&rt3<0.7

T3=30;

elseifrt3>=0.7&rt3<0.9

T3=32;

else

T3=34;

end

ifT3<T1+T2zcount=count+l;end

end%for

prob=count/n

模擬結(jié)果：每次模擬1000次或5000次

序號(hào)模擬次數(shù)近似概率p

110000.6280

210000.6920

310000.6530

450000.6490

550000.6260

650000.6288

系統(tǒng)模擬注意事項(xiàng)：

?一次模擬結(jié)果毫無(wú)意義！

?模擬是試驗(yàn)性的,是思維結(jié)果的驗(yàn)證。

?必須進(jìn)行足夠多次的模擬，并對(duì)結(jié)果進(jìn)行統(tǒng)計(jì)分

析。

系統(tǒng)模擬特點(diǎn)：系統(tǒng)模擬是研究系統(tǒng)，特別是動(dòng)態(tài)

系統(tǒng)的重要方法，對(duì)于：

?結(jié)構(gòu)復(fù)雜的系統(tǒng)；

?很難用解析方法求出變量關(guān)系的系統(tǒng)；

?內(nèi)部機(jī)理不明的“黑箱”系統(tǒng)；

?為驗(yàn)證用其他方法建立的模型及結(jié)果,應(yīng)是較好

的選擇。

1.3隨機(jī)變量的建模

?利用理論分布，基于對(duì)問(wèn)題的實(shí)際的、合理的假

設(shè)，選擇適當(dāng)?shù)睦碚摲植寄M隨機(jī)變量

優(yōu)點(diǎn)是給出了各種理論結(jié)果出現(xiàn)的概率，便于進(jìn)行數(shù)學(xué)分

析和處理。但此方法僅限于十分簡(jiǎn)單的情況，當(dāng)問(wèn)題越復(fù)

雜，數(shù)學(xué)處理變得越困難，并且丟失了試驗(yàn)數(shù)據(jù)的信息。

?基于實(shí)際數(shù)據(jù)的頻率做近似模擬

優(yōu)點(diǎn)在于完全與觀察數(shù)據(jù)相符，并且隨實(shí)際問(wèn)題的復(fù)雜程

度增大不會(huì)產(chǎn)生更大的困難，僅增大工作量而已。缺點(diǎn)是

不便于進(jìn)行數(shù)學(xué)分析，不得不依賴于模擬得到的統(tǒng)計(jì)結(jié)果。

?應(yīng)用常將兩種模擬方法結(jié)合起來(lái)使用

1.4均勻分布隨機(jī)變量模擬

1.4.1平方取中法

。1.從一個(gè)4位數(shù)對(duì)開(kāi)始，稱為種子；

。2.將它平方得到一個(gè)8位數(shù)(必要時(shí)前面加0)；

。3.取中間的4位數(shù)作為下一個(gè)隨機(jī)數(shù)。

。按照如上方式，就能得到0—9999隨機(jī)出現(xiàn)的整

數(shù)。因而可以用來(lái)模擬(0,1)上均勻分布的隨機(jī)

數(shù)。

。缺點(diǎn)：會(huì)退化為0或幾個(gè)數(shù)的循環(huán)，所以現(xiàn)在一般不

采用這種方法

。思考：如何皆斷隨機(jī)數(shù)生成算法的好壞？

■functionsimrandl

%平方取中法:會(huì)退化為0

x0=6532161%2040,2041不好

%653217出現(xiàn)循環(huán)

n=6;%n位整數(shù)

count=0;disp(sprintf(n'))

whilexO?=0%fori=l:100,

t=x02;t=mod(t,10(n+n/2));

t=fix(t/10(n/2));disp(sprintf('%10d\t))

x0=t;

count=count+1;

end

1.4.2線性同余法

。給定4個(gè)正整數(shù)包仇gio

。產(chǎn)生規(guī)則：

一葉1=（。/九+6）modc,c可以取很大，如c=231。

?特點(diǎn)（問(wèn)題）：

。1.序列不斷重復(fù)（循環(huán)周期至多「?jìng)€(gè)數(shù)）：

?2.序列各數(shù)之間缺少統(tǒng)計(jì)獨(dú)立性

。3,序列是個(gè)確定性

O兩組在實(shí)踐中效果較好的參數(shù)設(shè)置：

七7計(jì)1=（5。/）mod-31）

工片+531

1=(8<rn)mod(2—1)

線性同余法程序：產(chǎn)生o—c整數(shù)

。產(chǎn)生n個(gè)隨機(jī)整數(shù)

?functionr=randint(x,aFb,c,n)

%輸入x為種子

%線性同余法

fori=l:n,

x=mod((a*x+b),c);

?=x;

end

。調(diào)用例子

。x=7;a=l,b=7,c=10;n=20;

r—randint(x,a,b.c,n)

。x=7;a=l,b=7,c=231;n=20;

r=randint(x,a,b,c,n)

。思考題：利用線性同余法實(shí)現(xiàn)如下功能

。產(chǎn)生(0,1)上的隨機(jī)數(shù)；

。產(chǎn)生U(a,b)均勻分布的隨機(jī)數(shù)；

。產(chǎn)生0,1，一n的隨機(jī)整數(shù)；

。產(chǎn)生從j到k(jWk)上的隨機(jī)整數(shù)。

。思考解答

(設(shè)randint為線性同余法隨機(jī)數(shù)生成函數(shù))

?(1)randint/c;

■(2)a+(b-a)*randint/c;

。(3)INT((n+l)*randint/c)

?(4)j+INT((k-j+l)*randint/c)

1.5Matlab隨機(jī)模擬函數(shù):

1.5.1常見(jiàn)分布隨機(jī)變量的模擬

?1.產(chǎn)生mxn階［0,1］均勻分布的隨機(jī)數(shù)矩

陣：rand(m,n)

。思考題：請(qǐng)問(wèn)如何用rand產(chǎn)生均勻分布U(a,b)的

隨機(jī)數(shù)？

?2.產(chǎn)生mxn階［a,b］均勻分布U(a,b)的隨

機(jī)數(shù)矩陣：unifrnd(a,b,m,n)

。思考題：在Matlab中用rand和unifrnd分別產(chǎn)生10個(gè)均

勻分布U(0,5)的隨機(jī)數(shù)。

。3.產(chǎn)生mxn階均值為〃，標(biāo)準(zhǔn)差為。的正態(tài)分布的

隨機(jī)數(shù)矩陣：normrndn)

/⑺=$expT(號(hào)丹/wR

V27T(7/°

?randn(m,n)產(chǎn)生標(biāo)準(zhǔn)正態(tài)分布的隨機(jī)數(shù)

。當(dāng)研究對(duì)象視為大量相互獨(dú)立的隨機(jī)變量之和，且其

中每一種變量對(duì)總和的影響都很小時(shí)，可以認(rèn)為該對(duì)

象服從正態(tài)分布。

。機(jī)械加工得到的零件尺寸的偏差、射擊命中點(diǎn)與目標(biāo)

的偏差、各種測(cè)量誤差、人的身高、體重等，都可近

似看成服從正態(tài)分布。

?4.產(chǎn)生mxn階期望值為〃的指數(shù)分布的隨機(jī)數(shù)矩

陣：exprnd（出m,n）

。若連續(xù)型隨機(jī)變量X的概率密度函數(shù)為

。其中X>0為常數(shù)，則稱X服從參數(shù)為人的指數(shù)分

布。

0指數(shù)分布的期望值為1/八

。排隊(duì)服務(wù)系統(tǒng)中顧客到達(dá)率為常數(shù)時(shí)的到達(dá)間

隔、故障率為常數(shù)時(shí)零件的壽命都服從小指數(shù)分

布。

。指數(shù)分布在排隊(duì)論、可靠性分析中有廣泛應(yīng)用。

?注意：Matlab中，產(chǎn)生參數(shù)為八的指數(shù)分布的命

令為exprnd(l/A)

顧客到達(dá)某商店的間隔時(shí)間服從參數(shù)為0」的指數(shù)分布

指數(shù)分布的均值為1/0.1=10。指兩個(gè)顧客到達(dá)商店的

平均間隔時(shí)間是10個(gè)單位時(shí)間.即平均10個(gè)單位時(shí)間到

達(dá)1個(gè)顧客.顧客到達(dá)的間隔時(shí)間可用exprnd(lO)模

擬。

。5.產(chǎn)生mxn階參數(shù)為人的泊松分布的隨機(jī)數(shù)矩

陣：poissrnd(AFmFn)

e設(shè)離散型隨機(jī)變量X的所有可能取值為0J,2,...

，且取各個(gè)值的概率為

Xke~x

P(X=A-)="=0J2???，

A-!

?其中》＞0為常數(shù)，則稱X服從參數(shù)為人的泊松分

布。

9泊松分布的期望值為人

?泊松分布在排隊(duì)系統(tǒng)、產(chǎn)品檢驗(yàn)、天文、物理等

領(lǐng)域有廣泛應(yīng)用。

。指數(shù)分布與泊松分布的關(guān)系：

如相繼兩個(gè)事件出現(xiàn)的間隔時(shí)間服從參數(shù)為A的指數(shù)

分布，則在單位時(shí)間間隔內(nèi)事件出現(xiàn)的次數(shù)服從參數(shù)

為人的泊松分布.即單位時(shí)間內(nèi)該事件出現(xiàn)k次的概率

為：

.Xke~x

p.X=卜）=———.卜=0.12…，

卜！

反之亦然。

(1)顧客到達(dá)某商店的間隔時(shí)間服從參數(shù)為0.1的指數(shù)

分布㈡(2)該商店在單位時(shí)間內(nèi)到達(dá)的顧客數(shù)服從參

數(shù)為0」的泊松分布

。說(shuō)明：

。(1)指兩個(gè)顧客到達(dá)商店的平均間隔時(shí)間是io個(gè)單位時(shí)

間.即平均10個(gè)單位時(shí)間到達(dá)1個(gè)顧客.

。(2)指一個(gè)單位時(shí)間內(nèi)平均到達(dá)0.1個(gè)顧客

敵坦克分隊(duì)對(duì)我方陣地實(shí)施突襲,其到達(dá)規(guī)律服從泊

松分布，平均每分鐘到達(dá)4輛.

(1)模擬敵坦克在3分鐘內(nèi)到達(dá)目標(biāo)區(qū)的數(shù)量,以

及在第1、2、3分鐘內(nèi)各到達(dá)幾輛坦克.

(2)模擬在3分鐘內(nèi)每輛敵坦克的到達(dá)時(shí)刻。

?求解思路：先明確隨機(jī)變量及其分布

。(1)每分鐘到達(dá)車輛數(shù)：4用poissmd(4)進(jìn)行模

擬。入=4

。(2)兩輛車到達(dá)的平均間隔時(shí)間0.25分鐘.

一坦克到達(dá)的間隔時(shí)間應(yīng)服從參數(shù)為4的負(fù)指數(shù)分

布，用exprnd(1/4)模擬。

。練習(xí)：如何用Matlab編程實(shí)現(xiàn)上述問(wèn)題的模擬？

1.5.2其它隨機(jī)變量的模擬

X01

概率0.50.5

。模擬方法:

()</?<().5.

r?U(0,l).

0.5<r<1.

o實(shí)現(xiàn)：

。u=rand

ifu<0.5,

x=o；

else

X0123

概率p0.30.20.40.1

9模擬方法：

'0.0<7-<0.3,

1.0.3<r<0,5,

?，'=j2.0,5<r<0,9,r?U(0,l).

3,0.9<r<1,

?u=rand

ifu<0.3,

X=0;

elseifu<0.5%0.3+0.2

X=l;

elseifu<0.9%0.3+0.2+0.4

X=2;

else

X=3;

end

2Monte-Carlo模擬

2.1Monte-Carlo原理

。蒙特卡羅(MonteCarlo)方法，也稱為隨機(jī)模擬

(randomsimulation)

?基本思想：為了解決數(shù)學(xué)、物理、工程技術(shù)等方

面的問(wèn)題

。1.首先建立一個(gè)概率模型或隨機(jī)過(guò)程，使它的參數(shù)等

于問(wèn)題的解；

?2.然后通過(guò)對(duì)模型或過(guò)程的觀察或抽樣試驗(yàn)來(lái)計(jì)算所

求參數(shù)的統(tǒng)計(jì)特征，最后給出所求解的近似值。

2.2蒙特卡羅法應(yīng)用

2.2.1求解非線性規(guī)劃

。對(duì)于非線性規(guī)劃問(wèn)題:

min/(/).￡eE"

((//(X)>0J=1,2,-??,m.

s.t.

\cij<Xj<%,j=1,2,???

?用蒙特卡羅法求解的基本思想是：在估計(jì)的區(qū)域

{(“1"：2,,??6

.，。幾升勺[ajj}j]j1,2,...,n

誓騁若f實(shí)驗(yàn)點(diǎn)，然后從試驗(yàn)點(diǎn)中找出可行

再?gòu)目尚悬c(diǎn)中選擇最小點(diǎn).

9基本假設(shè)

。試驗(yàn)點(diǎn)的第/個(gè)分量叼服從對(duì).%］內(nèi)的均勻分布.

。符號(hào)假設(shè)

。P：試驗(yàn)點(diǎn)總數(shù)：MAXP：最大試驗(yàn)點(diǎn)總數(shù)；

。K:可行點(diǎn)總數(shù)：MAXK：最大可行點(diǎn)數(shù)；

。X*：迭代產(chǎn)生的最優(yōu)占?

。Q：迭代產(chǎn)生的最小宿i('x*),其初始值為計(jì)算機(jī)所能

表示的最大數(shù).

。算法思路

。先產(chǎn)生一個(gè)隨機(jī)數(shù)作為初始試驗(yàn)點(diǎn)，以后則將上一個(gè)試

驗(yàn)點(diǎn)的第j個(gè)分量隨機(jī)產(chǎn)生，其它分量不變而產(chǎn)生一新

的試驗(yàn)點(diǎn).這樣，每產(chǎn)生一個(gè)新試驗(yàn)點(diǎn)只需一個(gè)新的

隨機(jī)數(shù)分量.當(dāng)K>MAXK或P>MAXP時(shí)停止迭代.

例子1

maxz———+&ii+

(3J：I+12=10

S.t<N1>0

X2>0

9在Matlab軟件包中編程，一般需要設(shè)計(jì)三個(gè)M-文

件：

。主程序；

。目標(biāo)函數(shù)程序；

。約束條件測(cè)試程序.

例子2（飛行管理模型的蒙特卡羅法求解程

序）

2.2.2估算圓周率

?(1)采用確定性模型估算圓周率

。?可題:請(qǐng)建立確定性模型估算國(guó)周率.

。利用定積分的幾何意義：估算曲邊梯形面積

。(2)采用蒙特卡羅法估算圓周率

平行四邊形ABCD的面積為22=4

圓的面積6=7F12=7T

?，F(xiàn)在模擬產(chǎn)生在正方形ABCD中均勻分布的點(diǎn)n個(gè)，如

果這n個(gè)點(diǎn)中有m個(gè)點(diǎn)在該圓內(nèi)，則圓的面積與正方

形ABCD的面積之比可近似為m/n；

°即工仁如今兀仁坦

'4nn

。模擬程序：估算圓周率：

9n=input('請(qǐng)輸入產(chǎn)生點(diǎn)的個(gè)數(shù)：)m=0;

fori=l:n,

x=-l+2*rand;y=-l+2*rand;

ifx~2+y-2<=Lm=m+l;end

end

mypi=4*m/n

a更簡(jiǎn)潔實(shí)現(xiàn)：

9n=10000;xy=-l+2*rand(2,n);

x=xy(l,:);y=xy(2,:);

mypi=4*sum(x「2+y「2j=l)/n

。下面是一組運(yùn)行結(jié)果：

模擬次數(shù)”7T的近似

100003.1028

100003.1395

1000003.1418

1000003.1403

。小結(jié)

。運(yùn)用蒙特卡羅法求解問(wèn)題要先建立概率模型或者隨機(jī)

過(guò)程。

。思考

?請(qǐng)建立更精確的估計(jì)圓周率的數(shù)學(xué)模型。

2.2.3估算定積分

例子3

求定積分J32dl=宗

圖3.1:直線〃=2,直線1=1與曲線〃=x2

1)頻率法

。分析：，

?條件:

?？紤]〃上)20的情形

?構(gòu)造一個(gè)矩形包含了曲邊梯形，矩形區(qū)域高"滿足：

(1>max/(z)

a<i<b

。利用大數(shù)定理，事件發(fā)生頻率依概率收斂與事件發(fā)生

的概率。

。具體思路：

。在矩形區(qū)域內(nèi)產(chǎn)生〃(足夠大)個(gè)點(diǎn)

。落在曲邊梯形內(nèi)的點(diǎn)為m個(gè)，則

?、艘?竺

1—57

其中矩形區(qū)域面積S=(6-a)cl

。則所求定積分

rb

If(x)dx?—(fc—d)d

?n=106;%產(chǎn)生隨機(jī)點(diǎn)的個(gè)數(shù)

a=0;%積分下限

b=l;%積分上限

maxfx=l;%函數(shù)最大值(估計(jì))

d=maxfx+l;%確定一個(gè)d(作為大矩形的高)

m=0;

fori=l:n,

x=a+rand*(b-a);%隨機(jī)產(chǎn)生x坐標(biāo)

y=d*rand;%隨機(jī)產(chǎn)生y坐標(biāo)

if丫〈=*八2,%在曲線下方

m=m+l;%計(jì)算增力口

end

end

s=m/n*d*(b-a)

。運(yùn)行程序:

/fd.r70.3332160()0()()0()0

Jo

。思考題

。(1)模型中的參數(shù)(如矩形區(qū)域的高d)對(duì)結(jié)果有影響

嗎？

。(2)這種方法對(duì)于函數(shù)值可正可負(fù)的情形該怎么應(yīng)用

呢？

2)平均值法

o平均值法是利用隨機(jī)變量的平均值(數(shù)學(xué)期望)

來(lái)計(jì)算定積分1

I=f3dH(3.1)

Ja

?？紤]一系列在(0,1)上均勻分布的獨(dú)立隨機(jī)變

量：△,&「??&，則｛/(&)｝，”1,2,…，7?是

一系列相互獨(dú)立同分布的隨機(jī)變量，有

1rbI

同/(&)=----/，3心=-----

I)—aLb—a

9由強(qiáng)大數(shù)定律知以概率為1成立

n

山】】！￡/(&)=￡(3.2)

72T8,=]u?

9因此當(dāng)〃足夠大時(shí)，可得近似公式

an

1=//(1)"七他一°葉￡/出)(3.3)

Jai=l

。由此得到計(jì)算/的平均值法；

。(1)產(chǎn)生在(0.1)上均勻分布隨機(jī)數(shù)4.孫…，九；

9(2)令Uq—a+(6—1限2：=1,2，7，；

9(3)計(jì)算寧作為/的估計(jì)值。

。分析式(3.1)?缶1(3.3),可發(fā)現(xiàn)本質(zhì)上平均值法是用樣

學(xué)期望的估計(jì)工

9%program:simJni.mean

%蒙特卡羅模擬法求定積分(平均值法)

「input('輸入產(chǎn)生隨機(jī)數(shù)個(gè)數(shù):')；

ifisempty(n)|(n<10000),

n=10000;

end

s=0;%存儲(chǔ)函數(shù)值之和

a=0;%積分下限

b=l;%積分上限

fori=l:n,

x=a+(b-a)*rand;%rand(l)

funvalue=x2;%計(jì)奠函數(shù)宿

S=s+funvalue;%結(jié)合函數(shù)計(jì)算函數(shù)值，并累計(jì)到變

量s

end

I=(b-a)*s/n

。運(yùn)行結(jié)果：

sim-int_inean

。輸入產(chǎn)生隨機(jī)數(shù)個(gè)數(shù)：iuo

I=

0.33795894941753

練習(xí)1

有30個(gè)電子部件D「D?,…,。30，它們的使用情況

如下:一開(kāi)始D、工作:若Di損壞則D2立即開(kāi)始工

作;若D2損壞則D3立即開(kāi)始工作,其余類推。

沒(méi)D的壽命X*i=l,2,…,30）服從參數(shù)

為入=0.1（小時(shí)）的指數(shù)分布且相互獨(dú)立。

令丁為30個(gè)器件的總使用時(shí)間，試用M—C（蒙特卡

羅）模擬法求解P{T230}和石（T尸。

3模擬模型案例

例1：某港口提供有足夠的泊位供船舶?？?，但現(xiàn)

在僅有一個(gè)可供裝卸的泊位，船舶先到則先進(jìn)行

裝卸，如果船舶得不到及時(shí)裝卸而造成的滯期費(fèi)

為每小時(shí)100元。現(xiàn)要弄清該系統(tǒng)的性能，重點(diǎn)

考察船舶進(jìn)入該港后等待裝卸的滯留時(shí)間以及卸

位（即裝卸用的泊位）的利用率，從而進(jìn)行經(jīng)濟(jì)

效益分析。

首先，對(duì)進(jìn)入該港口的100艘船進(jìn)行實(shí)際調(diào)查，記錄其活動(dòng)情況，得到

這100艘船到達(dá)港口的時(shí)間間隔和裝卸時(shí)間的分布情況的頻數(shù)和累積頻

率分布。

表1100艘船到達(dá)港口的時(shí)間間隔頻數(shù)表

到達(dá)間隔56789101112131415161718

（小時(shí)）

到達(dá)船舶136791011111197654

數(shù)（艘）

表2100艘船裝卸時(shí)間頻數(shù)表

裝卸時(shí)間（小時(shí)）910111213141516

裝卸船舶數(shù)（艘）2022191610832

表3船舶到達(dá)港口的時(shí)間間隔累積分布表

到達(dá)間隔56789101112131415161718

累積頻率0.010.040.100.170.260.360.470.580.690.780.850.910.961.00

表4船舶裝卸時(shí)間累積分布表

裝卸時(shí)間910111213141516

累積頻率0.200.420.610.770.870.950.981.00

為了比較準(zhǔn)確地反映系統(tǒng)的性能，我們采用穩(wěn)態(tài)

模擬方式，可以考察系統(tǒng)運(yùn)行360天的情況。假

定港口每天24小時(shí)連續(xù)工作，因此模擬時(shí)間T設(shè)

置為8640小時(shí)?？紤]到船舶到達(dá)港口的事件是影

響整個(gè)系統(tǒng)狀態(tài)的主要因素，可以將模擬事件設(shè)

置為該事件的發(fā)生時(shí)刻。

所用變量說(shuō)明：

t:一艘船到達(dá)港口的時(shí)間；

tl:前一艘船駛離港口的時(shí)間；

td:兩艘船到達(dá)港口的時(shí)間間隔；

ts：當(dāng)前船舶裝卸所需時(shí)間；

tw：t時(shí)刻所有已到達(dá)港口的船舶等待裝卸時(shí)

間^總、和；

tf:t時(shí)刻裝卸位已空閑時(shí)間總和；

表5某港口對(duì)船舶服務(wù)的模擬結(jié)果

模擬到達(dá)船到港船舶滯留船舶平均滯港口支付滯裝卸泊位空裝卸泊位

次數(shù)舶數(shù)總時(shí)間（小時(shí)）留時(shí)間（小時(shí)）留費(fèi)（元）閑時(shí)間（小時(shí)）利用率

1736748010.1674800035396%

272243696.0543690051594%

3737910912.3690190042695%

473657587.8257580038196%

5741750110.1275010034196%

673968409.2668400042295%

經(jīng)過(guò)6次模擬計(jì)算，裝卸泊位的平均利用率為

95.3%,但到達(dá)的船舶的卻因得不到及時(shí)的服務(wù)

而造成平均每艘滯留9.295小時(shí)，每年到港船舶總

滯留6843小時(shí)，因此港口每年約需支付68萬(wàn)元的

滯留費(fèi)，這顯然是一筆不小的開(kāi)支。為此，可修

改模擬模型，考慮增設(shè)一個(gè)裝卸泊位，并重新進(jìn)

行模擬運(yùn)行，將這些數(shù)據(jù)提供給決策者一確定是

否需要再新建一個(gè)裝卸泊位。

練習(xí)：某公共汽車站每隔30分鐘到達(dá)一輛汽車，但

可能有［0,3］分鐘的誤差，此誤差大小與前一輛汽

車的運(yùn)行無(wú)關(guān)。汽車最多容納50名旅客，到達(dá)該

汽車站時(shí)車內(nèi)旅客人數(shù)服從［20,50］的均勻分布。

到站下車的旅客人數(shù)服從［3,7］的均勻分布，每

名旅客下車的時(shí)間服從［1,7］秒的均勻分布。旅

客按每30分鐘到達(dá)12個(gè)人的泊松分布到達(dá)汽車站,

單列排隊(duì)等車，先到先上，如果某位旅客未能上