隨機(jī)模擬(仿真)-simulation

1、隨機(jī)模擬（仿真）-simulation仿真仿真（也稱為模擬模擬）：就是用計(jì)算機(jī)程序在計(jì)算機(jī)上模仿各種實(shí)際系統(tǒng)的運(yùn)行過程，并通過計(jì)算了解系統(tǒng)隨時(shí)間變化的行為或特性。計(jì)算機(jī)仿真計(jì)算機(jī)仿真：是在已經(jīng)建立的數(shù)學(xué)、邏輯模型之上，通過計(jì)算機(jī)實(shí)驗(yàn)，對一個(gè)系統(tǒng)按照一定的決策原則或作業(yè)規(guī)則，由一個(gè)狀態(tài)變換為另外一個(gè)狀態(tài)的行為進(jìn)行描述和分析。實(shí)際問題數(shù)學(xué)、邏輯模型計(jì)算機(jī)模型數(shù)學(xué)、計(jì)算機(jī)解實(shí)際解為什么要進(jìn)行仿真 實(shí)際系統(tǒng)建立之前，要對系統(tǒng)的行為或結(jié)果進(jìn)行分析研究； 有些真實(shí)系統(tǒng)做實(shí)驗(yàn)會(huì)影響系統(tǒng)運(yùn)行，例如，在生產(chǎn)中任意改變工藝系數(shù)可能導(dǎo)致廢品，在經(jīng)濟(jì)活動(dòng)中隨意將一個(gè)決策付諸行動(dòng)會(huì)導(dǎo)致經(jīng)濟(jì)混亂； 在系統(tǒng)上做多次試驗(yàn)，很

2、難保證每次試驗(yàn)的操作條件相同，因而對實(shí)驗(yàn)結(jié)果好壞很難作出正確的判斷； 當(dāng)人是系統(tǒng)的一部分時(shí)，他的行為往往實(shí)驗(yàn)結(jié)果有所影響，這時(shí)，最好進(jìn)行模擬研究； 實(shí)驗(yàn)時(shí)間太長，費(fèi)用太大，或者有危險(xiǎn)，使得試驗(yàn)不容易進(jìn)行； 有些系統(tǒng)一旦建立起來后無法復(fù)原，例如，建立大型企業(yè)，要分析社會(huì)和經(jīng)濟(jì)效益，不能用建立起來試試看的辦法。那些問題適合計(jì)算機(jī)仿真解決 難以用數(shù)學(xué)公式表示的系統(tǒng)，或者沒有求解數(shù)學(xué)模型的有效方法； 雖然可以用解析的方法解決問題，但是問題的分析與計(jì)算過于復(fù)雜，這時(shí)計(jì)算機(jī)仿真可能提供簡單可行的求解方法； 希望在較短時(shí)間內(nèi)觀察到系統(tǒng)的發(fā)展全過程，以估計(jì)某參數(shù)對系統(tǒng)行為的影響； 難以在實(shí)際環(huán)境中進(jìn)行試驗(yàn)和觀

3、察，計(jì)算機(jī)仿真是唯一的方法； 需要對系統(tǒng)或過程進(jìn)行長期運(yùn)行的比較，從大量方案中尋找最優(yōu)。仿真的分類 模擬是系統(tǒng)狀態(tài)隨時(shí)間而變化的動(dòng)態(tài)寫照，因此，通常時(shí)間是模擬的主要自變量，其它的變量為因變量。 （1）按照模擬過程中因變量的變化情況，可以將模擬分為離散、連續(xù)、混合3種類型； （2）如果采用模擬計(jì)算機(jī)、采用數(shù)字計(jì)算機(jī)以及聯(lián)合使用則分為模擬仿真、數(shù)字仿真以及混合仿真； （3）根據(jù)仿真變量的特征分為隨機(jī)模擬隨機(jī)模擬仿真和模糊模擬仿真。模擬的方法模擬的方法 設(shè)計(jì)正確的模擬時(shí)間推進(jìn)機(jī)理是進(jìn)行模擬的一個(gè)非常重要的問題，模擬過程應(yīng)該根據(jù)系統(tǒng)的特征正確推進(jìn)模擬時(shí)間，使系統(tǒng)中各項(xiàng)要素與發(fā)生的事件保持同步，推進(jìn)時(shí)間

4、模擬的基本方法有： （1）下次事件法：是將模擬時(shí)間由一個(gè)事件發(fā)生時(shí)間點(diǎn)推進(jìn)到緊接著下一個(gè)事件發(fā)生的時(shí)間點(diǎn)。既時(shí)間變化幅度由事件變化確定； （2）固定時(shí)間步長法：此種方法模擬時(shí)間每次均以相等的固定步長向前推進(jìn)，每到達(dá)一個(gè)新的模擬時(shí)間點(diǎn)需要檢查相應(yīng)的時(shí)間段內(nèi)是否發(fā)生了事件；（有時(shí)候需要?jiǎng)討B(tài)調(diào)整步長）模擬的一般步驟模擬的一般步驟 明確問題，建立模型：正確描述研究的問題，明確規(guī)定模擬的目標(biāo)和任務(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í)現(xiàn)往往離不開大量數(shù)據(jù)的輸入，且需要確定隨機(jī)因

5、素的概率分布特性，并以此為抽樣的根據(jù)； 編制程序：模擬運(yùn)行，選擇適當(dāng)?shù)挠?jì)算機(jī)語言，按照系統(tǒng)數(shù)學(xué)、邏輯模型編寫計(jì)算機(jī)程序。 分析模擬輸出結(jié)果：一般包括如下幾個(gè)方面 （1）模擬結(jié)果的統(tǒng)計(jì)特性：均值、方差以及置信區(qū)間； （2）靈敏度分析； （3）根據(jù)確定的目標(biāo)函數(shù)，在眾多的實(shí)現(xiàn)方案中選取最優(yōu)方案。仿真的基礎(chǔ)仿真的基礎(chǔ) 概率論的大數(shù)定理； 微積分基礎(chǔ)； 各種分布：在matlab模擬中，常用的隨機(jī)分布為unif(均勻分布),exp(指數(shù)分布),norm(正態(tài)分布),chi2(2分布),t(t分布),f(F分布),bino(二項(xiàng)分布),poiss(泊松分布),unid(整數(shù)均勻分布)； 有關(guān)分布的計(jì)算功能

6、：pdf(概率密度),cdf(分布函數(shù) ),inv(逆概率分布),stat(均值與方差),rnd(隨機(jī)數(shù)生成)； 將分布與計(jì)算功能結(jié)合在一起，就是實(shí)現(xiàn)某種及算功能，例如，expinv就是計(jì)算指數(shù)分布的逆分布。用隨機(jī)模擬計(jì)算積分用隨機(jī)模擬計(jì)算積分例1 如下圖所示，在正方形內(nèi)有1/4單位圓。向正方形內(nèi)投小石頭，假設(shè)每次都能夠投進(jìn)正方形內(nèi)且可以落在正方形內(nèi)任何一點(diǎn)。問，小石頭落在1/4單位圓內(nèi)(包含邊界)的概率多大？01x1y分析：假設(shè)頭入正方形內(nèi)的石頭有n塊，有k塊落入了1/4單位圓內(nèi)。P為小石頭落入1/4單位圓內(nèi)的概率。那么根據(jù)Bernoull(伯努利)大數(shù)定理，有1|pnk|plimn即，當(dāng)實(shí)

7、驗(yàn)次數(shù)n充分大時(shí)，頻率和概率之差小于任意數(shù)的概率趨于1。而另外一方面由幾何概型有4Ap樣本空間的測度的測度事件在實(shí)際操作中，實(shí)驗(yàn)次數(shù)n不可能趨于無窮大，所以有nk4p44p(數(shù)學(xué)模型)對于估計(jì)nk4只有不斷重復(fù)做實(shí)驗(yàn)，這種試驗(yàn)可以具體去操作，(均勻投石塊，然后數(shù)數(shù)，這樣需要較高成本)。也可以讓計(jì)算機(jī)去重復(fù)試驗(yàn)，但是需要將數(shù)學(xué)模型轉(zhuǎn)化為計(jì)算機(jī)模擬模型(讓計(jì)算機(jī)完成均勻投石塊，自動(dòng)計(jì)數(shù)，也需要成本)。用計(jì)算機(jī)模擬投石塊過程和步驟如下：1、自動(dòng)生成隨機(jī)點(diǎn)0,1x0,1，模擬石塊在正方形內(nèi)的任意位置，用(xi,yi)表示，共n個(gè)點(diǎn)；2、判別(xi,yi)是否滿足xi2+yi21，即判別石塊是否落在1/

8、4單位圓內(nèi)，共k個(gè)點(diǎn)滿足；3、整理、統(tǒng)計(jì)模擬結(jié)果，用4k/n估計(jì)。利用matlab編程計(jì)算過程1、編寫M文件：function pai=fangzhenpai(n)%生成均勻分布;X=unifrnd(0,1,n,2);%判別是否落在1/4圓內(nèi);k=0;for k1=1:n y=X(k1,1)2+X(k1,2)2; if y for k=1:length(n)pai(k)=fangzhenpai(n(k);end paipai =2.0000 3.2000 3.3400 3.1740 3.1464 3.1460 3.1393 3.1423值得注意注意的是，每次模擬即或程序相同，不同人不同計(jì)算時(shí)間

9、，結(jié)果一般不同(因?yàn)殡S機(jī))。例2計(jì)算定積分dxx1322分析：若隨機(jī)變量X的概率分布密度是P(x)(axb)，則隨機(jī)變量Y=f(X)的數(shù)學(xué)期望為badx)x(p)x(f)X(f E)Y(E當(dāng)X是a,b上的均勻分布時(shí)，有bxa,ab1)x(p將p(x)代入上面的期望算子，有)X( f E)ab(dx)x( fdxab)x( f)X( f Ebaba 而另外一方面，根據(jù)大數(shù)定理，設(shè)y1,y2,yn為來自總體Y的一組樣本，有1| )y(Ey|limn其中)X(f E)Y(E)y(E,yn1yn1ii所以有niin1iiba)x(fnabny)ab(y)ab()X(f E)ab()Y(E)ab(dx)

10、x(f其中,x1,x2,xn為來自服從a,b上均勻分布總體X的一組樣本。且yi=f(xi)，i=1,2,n。dxx1322那么計(jì)算的步驟如下：1、在2，3上抽樣x1,x2,xn；2、計(jì)算n,.,2, 1i ,x1y2ii3、整理估計(jì)定積分；n1i2i322x1n23dxx1將上述過程編寫M文件如下function I=jifen1(a,b,n)X=unifrnd(2,3,n,1);L=(b-a)/n;for k=1:n Y(k)=(1+X(k)2)(1/2);endI=L*sum(Y); a=2;b=3;n=30000; I=jifen1(a,b,n)I = 2.6912計(jì)算結(jié)果如下：而該積分

11、的準(zhǔn)確值為：syms x int(1+x2)(1/2),2,3) ans = 3/2*10(1/2)-1/2*log(-3+10(1/2)-5(1/2)-1/2*log(2+5(1/2) double(ans)ans = 2.6948例3 估計(jì)事件發(fā)生的概率9XX, 3XXPr243221 其中X1服從均勻分布U2,5，X2服從指數(shù)分布exp(3)，X3和X4分別服從正態(tài)分布 N(3,2)和N(1,1)。1、在總體X=(X1,X2,X3,X4)中抽樣，得到樣本(x1k,x2k,x3k,x4k)，k=1,2,N；2、判斷隨機(jī)樣本是否滿足不等式組，假設(shè)滿足不等式組的有n個(gè)；3、用n/N估計(jì)。步驟步

12、驟：編寫M文件：function p=fangzhenguji(N)x1=unifrnd(2,5,N,1);x2=exprnd(3,N,1);x3=normrnd(3,2,N,1);x4=normrnd(1,1,N,1);n=0;for k=1:N y1=x1(k)+x2(k)2; y2=x3(k)+x4(k)2; if y1=3&y2b-0.00001); rr=bb(1); f1=f1;f(rr); f(rr)=;end bar=f1(n); fbar=jieji(10000,0.8)fbar = 4.9698計(jì)算結(jié)果：例5一個(gè)排隊(duì)系統(tǒng)的仿真 某店只有一個(gè)收款臺，顧客到達(dá)收款臺的時(shí)間間隔服

13、從均值為4.5的負(fù)指數(shù)分布，每個(gè)顧客的服務(wù)時(shí)間服從均值為3.2，標(biāo)準(zhǔn)差為0.6的正態(tài)分布，對100位顧客去收款臺繳款的排隊(duì)過程進(jìn)行仿真。分析分析： 排隊(duì)系統(tǒng)涉及的參考變量為：顧客的等待時(shí)間Wq，顧客的逗留時(shí)間W,排隊(duì)長Lq,隊(duì)長L，服務(wù)臺的空閑時(shí)間I，服務(wù)臺的繁忙時(shí)間B，排隊(duì)系統(tǒng)的狀態(tài)概率P。最后需要統(tǒng)計(jì)出有關(guān)指標(biāo)的平均值均方差。合理假設(shè)合理假設(shè)1、進(jìn)入系統(tǒng)的顧客沒有人因?yàn)椴辉敢獾榷x開系統(tǒng)；2、先到先服務(wù)；仿真步驟：仿真步驟：1、產(chǎn)生初始事件表，并從初始事件表中找出最靠前的 事件；2、啟動(dòng)仿真時(shí)鐘；3、判別當(dāng)前事件是哪一個(gè)事件，如果是顧客到達(dá)事 件，則啟動(dòng)顧客到達(dá)子程序；如果是服務(wù)事件，則

14、啟動(dòng)服務(wù)結(jié)束子程序；4、重復(fù)2-3，直到仿真完畢；5、整理、分析并輸出結(jié)果。顧客到達(dá)的子程序顧客到達(dá)的子程序1、產(chǎn)生下一個(gè)顧客到來的時(shí)刻，記入事件表；2、判別服務(wù)臺是否忙s=1?，若是，則顧客則排隊(duì)長L=L+1;否則，服務(wù)臺空閑,L=0;3、產(chǎn)生一個(gè)服務(wù)結(jié)束時(shí)刻，記入事件表；4、統(tǒng)計(jì)所需的數(shù)據(jù)；服務(wù)結(jié)束子程序服務(wù)結(jié)束子程序1、結(jié)束服務(wù)，產(chǎn)生結(jié)束時(shí)間，記入事件表；2、判別是否有顧客排隊(duì)，L=0?，若L=0，則服務(wù)臺空閑；否則，排隊(duì)長L=L-1;3、產(chǎn)生服務(wù)結(jié)束時(shí)刻，記入事件表；4、統(tǒng)計(jì)所需的數(shù)據(jù)；產(chǎn)生初始事件表從事件表中找出最靠前的事件仿真時(shí)鐘步進(jìn)是哪一類事件顧客到達(dá)服務(wù)結(jié)束產(chǎn)生下一個(gè)顧客到達(dá)紀(jì)

15、錄，記入事件表結(jié)束服務(wù)，產(chǎn)生結(jié)束時(shí)間，記入事件表出納是否忙？出納空閑，s=0顧客排隊(duì)，Lq=Lq+1否是產(chǎn)生服務(wù)結(jié)束時(shí)刻，記入事件表是否顧客排隊(duì)隊(duì)長減1，Lq-1出納空閑，s=0 否產(chǎn)生服務(wù)結(jié)束時(shí)刻，記入事件表是統(tǒng)計(jì)所需數(shù)據(jù)仿真是否完畢是結(jié)束程序否M文件程序:paiduisimu.mfunction ndaoda,nlikai,t,LLq,LL,BI,tdaoda,tfuwu=paiduisimu(n,a,b,c)%L時(shí)鐘進(jìn)步步長；n需要處理的人數(shù)；a表示負(fù)指數(shù)分布參數(shù)；b,c表示正態(tài)分布參數(shù);t=0;%起步時(shí)鐘;LL=;%記錄隊(duì)長狀態(tài);LLq=;%記錄排隊(duì)長狀態(tài);BI=;%記錄服務(wù)窗口空閑或

16、繁忙狀態(tài);tdaoda=;%記錄顧客到達(dá)時(shí)刻;tfuwu=;%記錄顧客接受服務(wù)離開的時(shí)刻;ndaoda=0;%記錄顧客到達(dá)累計(jì)總數(shù);nlikai=0;%記錄顧客離開累計(jì)總數(shù);L=0; %隊(duì)長，系統(tǒng)中的顧客數(shù);Lq=0;%排隊(duì)長，系統(tǒng)中排隊(duì)等待的顧客數(shù);t1=exprnd(a);t2=+inf;%t2=normrnd(b,c);%t1,t2初始事件表;s=0; %記錄窗口繁忙狀態(tài),s=0表示空閑;s=1表示窗口繁忙;while nlikain; if t1t2 %啟動(dòng)服務(wù)結(jié)束子程序; t=t2;%仿真時(shí)鐘進(jìn)步; tfuwu=tfuwu,t; if Lq0 s=1; L=L-1; Lq=Lq-1;

17、 r2=normrnd(b,c); t2=t+r2; else s=0; L=0;Ls=max(L-1,0); t2=+inf; end LL=LL,L; LLq=LLq,Lq; BI=BI,s; nlikai=nlikai+1; ndaoda=ndaoda; elseif t1=t2 %兩個(gè)子程序同時(shí)啟動(dòng); t=t1; %仿真時(shí)鐘進(jìn)步; tdaoda=tdaoda,t; tfuwu=tfuwu,t; ndaoda=ndaoda+1; nlikai=nlikai+1; L=L; Lq=Lq; s=s; LL=LL,L; LLq=LLq,Lq; BI=BI,s; t1=t+exprnd(a);

18、t2=t+normrnd(b,c); endend ndaoda,nlikai,t,LLq,LL,BI,tdaoda,tfuwu=paiduisimu(n,a,b,c);ndaoda = 100nlikai = 100t = 489.3746 mean(LLq)ans = 1.0950 mean(LL)ans = 1.9300 mean(BI)ans = 0.8350在例5中，如果窗口數(shù)為2，那么排隊(duì)情況會(huì)怎樣？ clear all n=100;a=4.5;b=3.2;c=0.6; lambda=0.1 ndaoda,nlikai,t,LLq,LL,BI1,BI2,tdaoda,tfuwu1=

19、paiduiL2(n,a,b,c,lambda);ndaoda = 100t = 497.5000 mean(LLq)ans = 0.1200 mean(LL)ans = 1.2300 mean(BI1)ans = 0.5600 mean(BI2)ans = 0.5500從結(jié)果看，情況有所改善。例6公交車調(diào)度問題公交車調(diào)度問題 在一條公共汽車上，有n個(gè)車站（兩個(gè)方向各n/2個(gè)），m輛汽車按照發(fā)車時(shí)刻表循環(huán)往復(fù)地由發(fā)車、中途運(yùn)行、到站、終點(diǎn)調(diào)頭、反向運(yùn)行等環(huán)節(jié)運(yùn)送等待上車的旅客，希望模擬此路線的交通運(yùn)行情況。為了簡化問題，對于交通路口、并行路線、和其它車輛的影響不予考慮。系統(tǒng)的實(shí)體有汽車、乘客、車站、發(fā)車場等；