蒙特卡羅模擬應(yīng)用

1、3.5 隨機(jī)模擬與系統(tǒng)仿真 一. 隨機(jī)現(xiàn)象的模擬 1. 隨機(jī)變量及其分布 隨機(jī)事件：在一定條件下有可能發(fā)生的事件。 概率：隨機(jī)事件發(fā)生的可能性的度量 P(A)， 0 P(A) 1. 隨機(jī)變量：在一定的范圍內(nèi)隨機(jī)取值的變量, “X=ak” (k=1,2,n), 或 隨機(jī)發(fā)生. badttpbaXP)(),(,baX 隨機(jī)變量的分布： 若已知 則稱 ak a1 a2 an P(“X=ak”) p1 p2 pn 為隨機(jī)變量 X 的分布列, 簡(jiǎn)稱 X 的分布 若已知 則稱 p(x) 為隨機(jī)變量 X 的分布密度, 簡(jiǎn)稱X的分布1,) (1nkkkkppaXPbadttpbaXP)(),(1)(dttp

2、稱 為兩點(diǎn)分布 單次貝努里試驗(yàn)的結(jié)果 稱 為二項(xiàng)分布 n 重貝努里試驗(yàn)的結(jié)果 稱 為離散的均勻分布 以相同的概率取所有可能的數(shù)值 稱 為泊松分布 發(fā)生率較低次數(shù)無(wú)限增大時(shí)貝努里試驗(yàn)的極限）（222exp21)(xxP) 1( ,)1 ()(pppCaXPknkknknknaXPk, 2 , 1,/ 1)(, 2 , 1,!)(kekkXPk011)(xpxpxP 稱 為均勻分布密度 在a,b的任何相等的子區(qū)間上取值的概率相同 稱 為正態(tài)分布密度 許多偶然因素作用結(jié)果的總和。 稱 為指數(shù)分布密度 質(zhì)點(diǎn)于隨機(jī)時(shí)間陸續(xù)到達(dá)的時(shí)間間隔，其它,0,1)(baxabxP）（222exp21)(xxP0,0

3、0,)(xxexpx 2. 隨機(jī)數(shù)和隨機(jī)變量的模擬 10. 隨機(jī)數(shù)（RND）: 計(jì)算機(jī)隨機(jī)數(shù)發(fā)生器產(chǎn)生的數(shù)串, 它在（0, 1) 中的分布是均勻的。 一般稱之為偽隨機(jī)數(shù)。 20. 具特定分布隨機(jī)變量的模擬 變量 X 有分布列 令 則有 1,) (1nkkkkppaXP1,)()1()(nkkppp0,)0(1)(pppkiik1,)() 1()(nkkppp 以 p(k)為分點(diǎn)，將0，1分為 n 個(gè)小區(qū)間 取隨機(jī)數(shù) R , 則容易證明 P( “ p(k-1) R p(k) ” ) = pk = P ( “ X = ak” ) 隨機(jī)事件 “ p(k-1) R p(k) ” 與 “ X=ak” 有

4、相同的概率分布。 可以使用隨機(jī)數(shù)在各小區(qū)間出現(xiàn)的情況 來(lái)模擬隨機(jī)事件 “ X=ak” 發(fā)生的狀況。 例 隨機(jī)變量 x = 0，1，2表示每分鐘到達(dá)超市收款臺(tái)的人數(shù)，有分布列 xk 0 1 2 pk 0.4 0.3 0.3 模擬十分鐘內(nèi)顧客到達(dá)收款臺(tái)的狀況 用用MATLAB模擬隨機(jī)事件模擬隨機(jī)事件 rand(20, 1) %生成生成20個(gè)均勻隨機(jī)數(shù)向量。個(gè)均勻隨機(jī)數(shù)向量。 randn(20,1) %生成生成20個(gè)正態(tài)隨機(jī)數(shù)向量。個(gè)正態(tài)隨機(jī)數(shù)向量。 r=rand(1,10); for i=1:10; if r(i)0.4 n(i)=0; elseif 0.4=r(i)&r(i) r r =

5、0.5678 0.7942 0.0592 0.6029 0.0503 0.4565 0.0185 0.8214 0.4447 0.6154 n n = 1 2 0 1 0 1 0 2 1 1 r r =0.2311 0.6068 0.4860 0.8913 0.7621 0.7919 0.9218 0.7382 0.1763 0.4057 n n = 0 1 1 2 2 2 2 2 0 1 二. 系統(tǒng)仿真（Simulation） 1. 系統(tǒng)仿真：使用計(jì)算機(jī)對(duì)一個(gè)系統(tǒng)的結(jié)構(gòu)和行為進(jìn)行動(dòng)態(tài)模擬。 為決策提供必要的參考信息。 特點(diǎn)：對(duì)象真實(shí)、復(fù)雜，進(jìn)行模仿。 2. 仿真模型：由計(jì)算機(jī)程序控制運(yùn)行 從

6、數(shù)值上模仿實(shí)際系統(tǒng)的動(dòng)態(tài)行為。 3. 仿真過程 1. 現(xiàn)實(shí)系統(tǒng)的分析： 了解背景，明確目的，提出總體方案。 2. 組建模型： 確定變量, 明確關(guān)系, 設(shè)計(jì)流程,編制程序 3. 運(yùn)行檢驗(yàn)： 確定初始狀態(tài)，參量數(shù)值， 運(yùn)行程序，檢驗(yàn)結(jié)果，改進(jìn)模型。 4. 輸出結(jié)果 三. 動(dòng)態(tài)系統(tǒng)的仿真 1. 時(shí)間步長(zhǎng)法： 把整個(gè)仿真過程分為許多相等的時(shí)間間隔 每個(gè)間隔為一個(gè)時(shí)間單位時(shí)間步長(zhǎng)。 在每個(gè)時(shí)間步長(zhǎng)內(nèi)模擬系統(tǒng)的動(dòng)態(tài)。 仿真時(shí)鐘：用以控制時(shí)間步進(jìn)的過程 （每一次步進(jìn)一個(gè)步長(zhǎng)） 例 3.15 池水含鹽 池中有水 2000 m3，含鹽 2 kg， 以 6m3 / 分 的速率向池中注入濃度為 0.5 kg / m

7、3 的鹽水， 又以 4 m3 / 分的速率從池中流出混合后的鹽水 問欲使池中鹽水濃度達(dá)到 0.2 kg / m3，需要多長(zhǎng)時(shí)間？ 系統(tǒng)分析： 池中有鹽水， 勻速注入濃鹽水， 勻速流出混合后的鹽水， 池中鹽水的濃度變化。 目的：仿真池中鹽水濃度的變化，給出達(dá)到給定濃度的時(shí)間。 變量、參量 時(shí)間 t，體積 V(t), 鹽量 S(t), 濃度 p(t); 流入流速 rI, 流入濃度 pI, 流出流速 rO, 流出濃度 p(t), 給定濃度 p* 時(shí)間步長(zhǎng) t , 打印步長(zhǎng) T. 關(guān)系: 在 t, t+ t 內(nèi)有 trrtVttVOI)()()(ttprprtSttSOII)()()()(/ )()(

8、ttVttSttp 動(dòng)態(tài)系統(tǒng)模擬的偽代碼 運(yùn)算 池水含鹽動(dòng)態(tài)系統(tǒng)模擬 變量 V(n)=時(shí)刻 n 池中鹽水體積 p(n)=時(shí)刻 n 池中鹽水濃度 S(n)=時(shí)刻 n 池中鹽水含鹽量 t = 時(shí)間單位 N = 模擬時(shí)間長(zhǎng)度 輸入 t，V(0), p(0), S(0), N 過程過程 Begin for n=1 to N do Begin V(n+1)V(n)+(rI-r0) t S(n+1)S(n)+ripi-r0p(n)t p(n+1)S(n+1)/V(n+1) End End 輸出輸出 V(1), V(2), , V(n) S(1), S(2), , S(n) p(1), p(2), , p(

9、n) 系統(tǒng)仿真流程圖 N Y 初始化V(0),S(0)仿真時(shí)鐘 t=0打印時(shí)鐘T=0計(jì)算V(t+t),S(t+t),p(t+t)時(shí)鐘步進(jìn)t=t+t,T=T+1p(t)p*Tm 打印輸出 參數(shù) rI=6，rO=4，pI=0.5，p*=0.2 初始 V(0)=2000, S(0)=2 步長(zhǎng) t =1 , m=10 結(jié)果 表3.15，t=185分, 鹽水濃度為 0.2. 模型IIOprtStVrdtdS)()()0()(,VrttVrrrdtdVOIIIOprtStVrdtdS)()(clft=1; v=2000;s=2;p=1/1000;%初態(tài)初態(tài)V=v(end); S=s(end); P=p(e

10、nd); x=0;%終終態(tài)態(tài)while p(end)=0.2 %最終閾值最終閾值 T=0; while T0.2 T=20;%轉(zhuǎn)輸出 end; end; x=x t-1;%時(shí)間記錄 V=V v(end); S=S s(end); P=P p(end);end; V1=103.*V; a=x,V1,S,P %輸出調(diào)節(jié) MATLAB實(shí)現(xiàn) x=x/2+1， 0 t 10, x(0)=0 建立M文件 fun . M function y=fun(t,x) y=x/2+1; 運(yùn)行結(jié)果 t0=0;tf=10; x0=0; t,x=ode23(fun,t0,tf,x0); plot(t,x) MATLAB實(shí)

11、現(xiàn) y+(y2-1)y-y=0, y(0)=0.25, y(0)=0 化為一階方程組 y1=y1(1-y22)-y2, y2=y1, 建立M文件 rhf . M function fyy=rhf(t,x) fyy=y(1).*(1-y(2).2)-y(2);y(1); 運(yùn)行結(jié)果 t0=0;tfinal=15;Y0=0.25,0; T,Y=ode23(rhf,t0,tfinal,Y0); plot(T,Y)求表達(dá)式求表達(dá)式（符號(hào)運(yùn)算）S=dsolve(Dx=(3-6*x)/(2000+2*t);求數(shù)值解求數(shù)值解 建立M文件 fun . M function y=fun(t,x) y=(3-6*x

12、)/(2000+2*t);t0=0; tf=200; x0=2; t, x=ode23(fun,t0 tf,x0);Plot(t, x);drdrVtVpptptptrdttdptVOtIII)()()()0()()()()()(000 問題 1. 在池水含鹽的問題中令 rO=rI=6m3/分 10. 池中鹽水的濃度如何變化？ 20. 若當(dāng)p(t) = 0.3kg/m3 時(shí)令pI=0, 需要多少時(shí)間達(dá)到 p*= 0.2kg / m3? 2. 若池中鹽水的初始濃度為 p*， 對(duì)于不同的初始體積V0，當(dāng)pI=0時(shí)， 計(jì)算池中鹽水濃度降低一半所用的時(shí)間 例 3.16 市場(chǎng)服務(wù) 超市有兩個(gè)出口的收款臺(tái)

13、， 兩項(xiàng)服務(wù)：收款、裝袋。 兩名職工在出口處工作。 有兩種安排方案： 開一個(gè)出口，一人收款、一人裝袋； 開兩個(gè)出口，每個(gè)人既收款又裝袋。 問商店經(jīng)理應(yīng)選擇哪一種收款臺(tái)的服務(wù)方案。 假設(shè)： 1. 顧客的到達(dá)是隨機(jī)的。 2. 收款裝袋的時(shí)間是相同的。 3. 第一種方案中，收款與裝袋同時(shí)進(jìn)行。 參量、變量： n(t) 到達(dá)顧客人數(shù)，L(t) 隊(duì)列長(zhǎng), 服務(wù)時(shí)間 T1 排隊(duì)時(shí)間累加, T2 服務(wù)時(shí)間累加 模型 n(t)=0, L(t)=0 則 L(t+t)=0, T2 (t+t)=T2 (t) L(t)0 則 L(t+t)=L(t)+n-1,T2 (t+t)= n(t)=1, L(t)=0 則 L(t

14、+t)=0, T2 (t+t)=T2(t)+ L(t)0 則 L(t+t)=L(t)+n-1,T2 (t+t) = n(t)=2, L(t)=0 則 L(t+t)=n-1, T2 (t+t)= L(t)0 則 L(t+t)=L(t)+n-1,T2 (t+t) =clfL=zeros(1,31);% L 等待的顧客人數(shù)，等待的顧客人數(shù)，T1=zeros(1,31);%T1等待時(shí)間的累加，等待時(shí)間的累加，T2=zeros(1,31); %T2服務(wù)時(shí)間的累加服務(wù)時(shí)間的累加,L1=zeros(1,31); %L1到達(dá)顧客人數(shù)累加。到達(dá)顧客人數(shù)累加。t=1;tau=1;x=0:30;r=rand(1,3

15、0);%隨機(jī)數(shù)隨機(jī)數(shù)%隨機(jī)模擬隨機(jī)模擬for i=1:30; t=t+1; if 0=r(i) & r(i)0.4 n=0; elseif 0.4=r(i) & r(i)eL)L3=sum(L(L2)%未被服務(wù)的顧客等待時(shí)間總和未被服務(wù)的顧客等待時(shí)間總和g1=(T1(end)-L3)/eL %平均等待時(shí)間平均等待時(shí)間g2=g1+tau %平均逗留時(shí)間平均逗留時(shí)間g3=eL/30 %平均每分鐘服務(wù)的顧客人數(shù)平均每分鐘服務(wù)的顧客人數(shù) t RND n L T1 T2 1 .81 2 1 1 1 2 .02 0 0 0 1 3 .60 1 0 0 1 4 .04 0 0 0 0 5 .

16、46 1 0 0 1 6 .31 0 0 0 0 7 .67 1 0 0 1 8 .25 0 0 0 0 9 .24 0 0 0 0 10 .10 0 0 0 0 11 .40 1 0 0 1 12 .02 0 0 0 0 t RND n L T1 T213 .39 0 0 0 0 14 .68 1 0 0 115 .08 0 0 0 016 .59 1 0 0 117 .66 1 0 0 118 .90 2 1 1 119 .12 0 0 0 120 .64 1 0 0 121 .79 2 1 1 122 .31 0 0 0 123 .86 2 1 1 124 .68 1 1 1 1 2.

17、事件表法（面向事件法）：每處理一個(gè)事件就前進(jìn)一步（每步的時(shí)間可能不同），以事件為中心安排。 對(duì)系統(tǒng)中的一系列不同類型的事件按發(fā)生的前后順序逐個(gè)進(jìn)行分析、處理。 事件表：記錄今后將要發(fā)生的事件及其屬性（發(fā)生的時(shí)間，事件的類型等），調(diào)度事件執(zhí)行的順序。 它隨著程序的執(zhí)行過程不斷地更新和補(bǔ)充以保證仿真過程有序地進(jìn)行。系統(tǒng)仿真初始事件表找出最靠前事件事件屬性顧客到達(dá)服務(wù)結(jié)束仿真完畢？輸出結(jié)果下一顧客到來(lái)時(shí)刻，記事件表出納忙？為顧客服務(wù)產(chǎn)生服務(wù)結(jié)束時(shí)刻，記事件表統(tǒng)計(jì)數(shù)據(jù)顧客排隊(duì)有顧客排隊(duì)？為顧客服務(wù)收款臺(tái)閑置產(chǎn)生服務(wù)結(jié)束時(shí)刻，記事件表統(tǒng)計(jì)數(shù)據(jù)結(jié)束服務(wù)YNNY例 市場(chǎng)服務(wù)假設(shè)：1. 顧客的到達(dá)收款臺(tái)是隨機(jī)

18、的，平均時(shí)間間隔為0.5分鐘，即間隔時(shí)間服從=2的(負(fù))指數(shù)分布。2. 對(duì)不同的顧客收款和裝袋的時(shí)間服從正態(tài)分布 N(1，1/3)。參量，變量： t(i): 第i位顧客到達(dá)時(shí)間， t2(i):第i位顧客受到的服務(wù)時(shí)間(隨機(jī)變量)，T(i): 第i位顧客離去時(shí)間.將第i位顧客到達(dá)作為第i件事發(fā)生; t(i+1)- t(i)= t1(i) (隨機(jī)變量)平衡關(guān)系：當(dāng) t(i+1)T(i) 時(shí), T(i+1)=t(i+1)+t2(i);否則, T(i+1)=T(i)+t2(i)模擬20位顧客到收款臺(tái)前的排隊(duì)情況。clft=zeros(1,21);%每位顧客到達(dá)時(shí)間T=zeros(1,21);%每位顧客

19、離去時(shí)間w=zeros(1,21);%顧客等待時(shí)間累加ww=zeros(1,21);%收款臺(tái)空閑時(shí)間累加r=exprnd(2,1,21);%服從指數(shù)分布的隨機(jī)數(shù) t2=normrnd(1,1/3,1,21);%服從正態(tài)分布的隨機(jī)數(shù)for i=1:1:20 t(i+1)=t(i)+t1(i); if t(i+1)=T(i); w(i+1)=w(i); T(i+1)=t(i+1)+t2(i); ww(i+1)=t(i+1)-T(i)+ww(i); else w(i+1)=T(i)-t(i+1)+w(i); T(i+1)=T(i)+t2(i); ww(i+1)=ww(i); end;end;b=t,

20、T,w,ww;%求隊(duì)伍長(zhǎng)brow,bcol=size(b);b=b,zeros(brow,1);b(1,bcol+1)=0;for j=2:brow l=0; if j-l-10 while b(j,1)=b(j-l-1,2) l=l+1; end; b(j,bcol+1)=l;end;end; bg1=w(end)/20 %平均等待時(shí)間g2=sum(T-t)/20 %平均逗留時(shí)間g3=T(end)/20 %平均每分鐘服務(wù)的顧客人數(shù) 3. 系統(tǒng)仿真軟件 MATLAB-SIMULINK Simulink 是一個(gè)用來(lái)對(duì)動(dòng)態(tài)系統(tǒng)進(jìn)行建模、仿真和分析的軟件包。它支持線性和非線性系統(tǒng)，連續(xù)和離散時(shí)間模型等。