數(shù)學建模-計算機模擬技術(shù)

1模擬基礎(chǔ)

2MonteCaHo模擬

3模擬模型案例

1模擬基礎(chǔ)

1.0模擬的背景、思路

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

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

子擴散”的科研項目取名為"Monte-Carlo"

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

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

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

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

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

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

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

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

模擬思路：

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

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

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

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

1.1模擬的基本知識

1.1.1模擬的概念及作用

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

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

中的一些隨機性因素很難用準確的數(shù)學公

式表述，從而也無法對整個系統(tǒng)采用解析

法求解。模擬是處理這類實際問題的有力

工具。

?模擬通常借助于計算機進行。

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

型的基礎(chǔ)之上，通過計算機試驗，對一個

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

一個狀態(tài)變換為另一個狀態(tài)的行為進行描

述和分析。

模擬的作用：

?對于很難用解析方法加以處理的問題，模

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

?對建模過程中的假設(shè)進行鑒定，對理論研

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

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

定的目標函數(shù)對不同方案進行比較，從中

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

1.1.2模擬的分類

通常，模擬時間是模擬的主要自變量。

設(shè)計正確的模擬時間推進機理：模擬過程中

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

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

推進模擬時間的基本方法：

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

的時間點推進到緊接著的下一次事件發(fā)生

的時間點。

?固定時間步長法：模擬時間每次均以相等

的固定步長向前推進，每到達一個新的模

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

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

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

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

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

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

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

時間推進方法：一般采用下次事件法

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

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

因

變

量

A

事件時間時間

圖1離散型模擬

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

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

F。，

T客達二

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

有.

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

。上

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

次，

程描述為圖幼扃而疑就講財比

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

化。在大多數(shù)計算機模擬過程中，按固定的步長

推進模擬時間。

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

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

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

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

1.1.3模擬的方式

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

達到T時，模擬終止。其性能指標明顯取決于系

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

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

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

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

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

定）。

1.1.4模擬的一般步驟

?明確問題，建立模型。正確描述待研究問題，明

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

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

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

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

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

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

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

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

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

?編制程序，模擬運行。

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

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

最優(yōu)方案。

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

征未知的隨機變量，每次運行的結(jié)果僅僅是對該

隨機變量所有觀察值總體的一次抽樣，對總體的

代表性很差，雖然可以增加模擬運行的時間從而

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

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

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

1.2隨機模擬案例：趕上火車的概率

【問題】如圖，一列火車從A站開往B站，某人每天

趕往B站上這趟火車.

工/某人

火車運行方向/

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

思考：請研究他能否趕上這趟火車。

他已了解到：

1）火車從A站到B站的運行時間是均值為30分

鐘，標準差為2分鐘的隨機變量；

2）火車在下午大約1點離開A站，離開時刻

的頻率分布如下：

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

頻率0.70.20.1

他到達B站的時刻的頻率分布為

時刻午后1:28午后1:30午后1:32午后1:34

頻率0.30.40.20.1

他能否及時趕上火車?

明確問題：他能及時趕上火車的概率是多少?

建模方向（思路）：

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

析解。（思考）

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

求近似解。

即先建立模擬模型，然后通過計算機模擬得到

問題的近似解。在同樣條件下多次試驗，計算

他能及時趕上火車的頻率。

問題分析：能及時趕上火車的充要條件是:

TH

即心―(<T2

其中

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

T2一火車的運行時間；

走什么變

丁3一他到達B站的時刻。重？如何

模擬?

基本假設(shè):

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

JLJ24T

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

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

建立模型：為了簡化計算，將下午L點、記為初始時

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

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

律分別為：

T1（分）0510

P(t)0.70.20.1

丁(分)

T3（分）2288303033223434

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

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

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

機變量%和T3,可以通過如下方法：

0,0<r<0.7,

4=<5,0.7<r<0.95

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

用來模擬隨機變

〃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è)計

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

輸出：趕上火車的頻率

兩種不同風格的算法描述

主要變量說明：

n模擬次數(shù)

k臨時變量，存儲當前累計模擬次數(shù)

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

算法1（分步驟描述）

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

第2步k=l,count=0

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

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

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

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

第7步由r及公式確定T3模擬人達到時刻

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

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

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

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

算法2（偽代碼描述）

i）初始化：

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

count=0;

道）模擬11次

fori=lton,

模擬隨機變量￡,丁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;當模擬隨機變量tl（火車從A站出發(fā)的時刻）

程序ifrtl<0.7

Tl=0;

elseifrtl>=0.7&rtl<0.9

Tl=5;

else?exid

T2=30+randn*2;%模取隨機菱量t2（火車的運行時間）

%模擬隨機變量t3（他到達B站的時亥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次

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

110000.6280

210000.6920

310000.6530

450000.6490

550000.6260

650000.6288

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

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

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

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

析。

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

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

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

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

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

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

的選擇。

1.3隨機變量的建模

?利用理論分布，基于對問題的實際的、合理的假

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

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

析和處理。但此方法僅限于十分簡單的情況，當問題越復

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

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

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

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

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

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

1.4均勻分布隨機變量模擬

1.4.1平方取中法

。1.從一個4位數(shù)對開始，稱為種子；

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

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

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

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

數(shù)。

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

采用這種方法

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

■functionsimrandl

%平方取中法:會退化為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個正整數(shù)包仇gio

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

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

?特點（問題）：

。1.序列不斷重復（循環(huán)周期至多「個數(shù)）：

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

。3,序列是個確定性

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

七7計1=（5。/）mod-31）

工片+531

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

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

。產(chǎn)生n個隨機整數(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)

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

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

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

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

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

。思考解答

(設(shè)randint為線性同余法隨機數(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隨機模擬函數(shù):

1.5.1常見分布隨機變量的模擬

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

陣：rand(m,n)

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

隨機數(shù)？

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

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

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

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

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

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

/⑺=$expT(號丹/wR

V27T(7/°

?randn(m,n)產(chǎn)生標準正態(tài)分布的隨機數(shù)

。當研究對象視為大量相互獨立的隨機變量之和，且其

中每一種變量對總和的影響都很小時，可以認為該對

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

。機械加工得到的零件尺寸的偏差、射擊命中點與目標

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

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

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

陣：exprnd（出m,n）

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

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

布。

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

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

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

布。

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

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

令為exprnd(l/A)

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

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

平均間隔時間是10個單位時間.即平均10個單位時間到

達1個顧客.顧客到達的間隔時間可用exprnd(lO)模

擬。

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

陣：poissrnd(AFmFn)

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

，且取各個值的概率為

Xke~x

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

A-!

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

布。

9泊松分布的期望值為人

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

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

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

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

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

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

為：

.Xke~x

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

卜！

反之亦然。

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

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

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

。說明：

。(1)指兩個顧客到達商店的平均間隔時間是io個單位時

間.即平均10個單位時間到達1個顧客.

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

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

松分布，平均每分鐘到達4輛.

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

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

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

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

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

擬。入=4

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

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

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

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

1.5.2其它隨機變量的模擬

X01

概率0.50.5

。模擬方法:

()</?<().5.

r?U(0,l).

0.5<r<1.

o實現(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)方法，也稱為隨機模擬

(randomsimulation)

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

面的問題

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

于問題的解；

?2.然后通過對模型或過程的觀察或抽樣試驗來計算所

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

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

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

。對于非線性規(guī)劃問題:

min/(/).￡eE"

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

s.t.

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

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

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

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

誓騁若f實驗點，然后從試驗點中找出可行

再從可行點中選擇最小點.

9基本假設(shè)

。試驗點的第/個分量叼服從對.%］內(nèi)的均勻分布.

。符號假設(shè)

。P：試驗點總數(shù)：MAXP：最大試驗點總數(shù)；

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

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

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

表示的最大數(shù).

。算法思路

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

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

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

隨機數(shù)分量.當K>MAXK或P>MAXP時停止迭代.

例子1

maxz———+&ii+

(3J：I+12=10

S.t<N1>0

X2>0

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

件：

。主程序；

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

。約束條件測試程序.

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

序）

2.2.2估算圓周率

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

。?可題:請建立確定性模型估算國周率.

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

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

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

圓的面積6=7F12=7T

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

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

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

°即工仁如今兀仁坦

'4nn

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

9n=input('請輸入產(chǎn)生點的個數(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更簡潔實現(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

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

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

100003.1028

100003.1395

1000003.1418

1000003.1403

。小結(jié)

。運用蒙特卡羅法求解問題要先建立概率模型或者隨機

過程。

。思考

?請建立更精確的估計圓周率的數(shù)學模型。

2.2.3估算定積分

例子3

求定積分J32dl=宗

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

1)頻率法

。分析：，

?條件:

?？紤]〃上)20的情形

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

(1>max/(z)

a<i<b

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

的概率。

。具體思路：

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

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

?、艘?竺

1—57

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

。則所求定積分

rb

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

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

a=0;%積分下限

b=l;%積分上限

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

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

m=0;

fori=l:n,

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

y=d*rand;%隨機產(chǎn)生y坐標

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

m=m+l;%計算增力口

end

end

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

。運行程序:

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

Jo

。思考題

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

嗎？

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

呢？

2)平均值法

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

來計算定積分1

I=f3dH(3.1)

Ja

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

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

一系列相互獨立同分布的隨機變量，有

1rbI

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

I)—aLb—a

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

n

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

72T8,=]u?

9因此當〃足夠大時，可得近似公式

an

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

Jai=l

。由此得到計算/的平均值法；

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

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

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

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

學期望的估計工

9%program:simJni.mean

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

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

ifisempty(n)|(n<10000),

n=10000;

end

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

a=0;%積分下限

b=l;%積分上限

fori=l:n,

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

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

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

量s

end

I=(b-a)*s/n

。運行結(jié)果：

sim-int_inean

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

I=

0.33795894941753

練習1

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

如下:一開始D、工作:若Di損壞則D2立即開始工

作;若D2損壞則D3立即開始工作,其余類推。

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

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

令丁為30個器件的總使用時間，試用M—C（蒙特卡

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

3模擬模型案例

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

在僅有一個可供裝卸的泊位，船舶先到則先進行

裝卸，如果船舶得不到及時裝卸而造成的滯期費

為每小時100元?，F(xiàn)要弄清該系統(tǒng)的性能，重點

考察船舶進入該港后等待裝卸的滯留時間以及卸

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

效益分析。

首先，對進入該港口的100艘船進行實際調(diào)查，記錄其活動情況，得到

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

率分布。

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

到達間隔56789101112131415161718

（小時）

到達船舶136791011111197654

數(shù)（艘）

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

裝卸時間（小時）910111213141516

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

表3船舶到達港口的時間間隔累積分布表

到達間隔56789101112131415161718

累積頻率0.010.040.100.170.260.360.470.580.690.780.850.910.961.00

表4船舶裝卸時間累積分布表

裝卸時間910111213141516

累積頻率0.200.420.610.770.870.950.981.00

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

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

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

置為8640小時?？紤]到船舶到達港口的事件是影

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

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

所用變量說明：

t:一艘船到達港口的時間；

tl:前一艘船駛離港口的時間；

td:兩艘船到達港口的時間間隔；

ts：當前船舶裝卸所需時間；

tw：t時刻所有已到達港口的船舶等待裝卸時

間^總、和；

tf:t時刻裝卸位已空閑時間總和；

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

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

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

1736748010.1674800035396%

272243696.0543690051594%

3737910912.3690190042695%

473657587.8257580038196%

5741750110.1275010034196%

673968409.2668400042295%

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

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

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

滯留6843小時，因此港口每年約需支付68萬元的

滯留費，這顯然是一筆不小的開支。為此，可修

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

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

否需要再新建一個裝卸泊位。

練習：某公共汽車站每隔30分鐘到達一輛汽車，但

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

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

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

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

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

客按每30分鐘到達12個人的泊松分布到達汽車站,

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