1、第四章 系統仿真w 4.1 基本概念w 4.2 蒙特卡羅仿真w 4.3 均勻隨機數的產生方法w 4.4 隨機變量的產生方法w 4.5 離散事件隨機系統的仿真方法w 4.6 4.6 隨機分布的確定隨機分布的確定w 4.7 系統動力學仿真w 4.8 4.8 仿真語言仿真語言w 4.9 4.9 可視化仿真技術可視化仿真技術w 4.10 4.10 虛擬現實技術虛擬現實技術w 4.11 4.11 仿真技術在礦業中的應用仿真技術在礦業中的應用w 作業4.1 4.1 基本概念基本概念 在下圖中，已知采區來煤量q1是一個隨機變量，根據統計結果，它服從正態分布，已知其均值=1.45t/ min，其標準差=0.0

2、95t/min。大巷是膠帶輸送機，其運輸量是常數：q2=1.33t/min。求8h后，采區最大煤倉容量。 q2q1煤倉投擲硬幣投擲硬幣隨機模擬實驗隨機模擬實驗 建立一個計算機能夠運算的模擬模型,問題的關鍵是用隨機數去模擬投擲硬幣的實驗。 所謂隨機數就是一組隨機出現的數列，通常用的是0，1區間的隨機數。這時要求0，1區間的隨機數必須具備下列性質：一是它的均勻性，即這些隨機數落在0，1區間內任一位置的概率是相同的；另一個是它的隨機性，即這些數在 0，1區間內任一位置出現與否是隨機的、獨立的。如以R表示隨機數，則出現0R0.5和 0.5R1的概率相等，都為0.5.這樣如以隨機數R0.5模擬出現正面，

3、以R0.5模擬出現反面，則模擬模型就建立起來了。 設設 I I 是投擲次數是投擲次數 N N 為正面出現次數為正面出現次數 R R 為隨機數為隨機數 P P是正面出現頻率是正面出現頻率 I0.50.5實驗開始 I=I+1產生隨機數R記錄出現正面次數N=N+1計算并輸出結果N，P模擬結束初始狀態I=0，N=0YNYN圖3-2 投擲硬幣實驗模擬框圖返回4.2 4.2 蒙特卡羅仿真蒙特卡羅仿真 設f(x)是區間（a,b）內的連續函數，且已知在該區間內f(x)=0,要求計算在（a，b）和（0，c）區間中分別產生兩個均勻分布的隨機數xi和yi，由此構成的坐標點 Pi(xi，yi)必定在矩形abde內。若

4、yif(xi)，則接受該點并作計數統計；否則舍棄該點不作統計. badxxf)(拒絕接受dxxbxjxia0cef(xi)f(xj) 設共產生N個隨機仿真點，如果共接受個m點，（即恰在曲線上或在曲線以下的點），當N充分大時有 曲線下面積 badxxfcabNm)()(返回4.3 4.3 均勻隨機數的產生方法均勻隨機數的產生方法 4.3.1 4.3.1 均勻隨機數均勻隨機數 0，1區間上的均勻隨機數，它的密度函數f(x)和分布F(x)分別為 其他0101)(xxf111000)(xxxxxF 4.3.2 4.3.2 均勻隨機數的產生方法均勻隨機數的產生方法 用數學公式產生的均勻隨機數，稱為偽隨機

5、數,常用的公式方法是線性同余法 . 若有a、b兩個整數被某個整數M除后，它們的余數相同，則稱a、b兩數對模M同余，記作 a=b(ModM) 求余數的公式為例如, MMyyc1510505151010515515c產生隨機數的線性同余法是一個遞推公式： 即 通過下式得到落在0，1區間內的隨機數： 兩個較好的計算公式如下 niMcaZZii,2, 1)(mod(1MMcaZcaZZiii11Mxrn1)2)(mod15(35115iiZZ)2)(mod453806246314159269(311iiZZ4.3.3 4.3.3 均勻隨機數檢驗方法均勻隨機數檢驗方法 設N個隨機數從小到大排列起來并分成

6、互不重疊的等區間的k段.從理論上講，這些隨機數應該均勻分布在每一段上，在每一段上理論頻數為如果統計出隨機數實際落在每段的頻數為ni，則令 x2值越大表明隨機數列越不均勻。 kikNmi,2,1,kiiiimmnx122)( 給定判定標準稱為顯著度,根據自由度r=k-1，從x2表中查得 值，如果計算所得的x250），若取零假設H0：=0，并取統計量 (近似服從N（0，1）標準正態分布)，同時選定顯著度=0.05，則根據概率論，當|u|1.96時為差異顯著，拒絕H0假設；當|u|1.96時接受假設，即檢驗的隨機數符合獨立性要求。4.3.4 4.3.4 計算機函數計算機函數 計算機高級語言有產生隨機

7、數的標準函數，可直接調用。 202/)(1SZZZNjNijiijjNuj返回4.4 4.4 隨機變量的產生方法隨機變量的產生方法 4.4.1 4.4.1 逆變換法逆變換法 已知某隨機變量的分布函數F（x） ，并知其定義域為0，1區間，這時可用0，1區間的均勻分布隨機數R來表示F（x），即F（x）=R，由此可得到隨機變量的樣本值x=F-1(R)。F-1為逆變換的表達式。xx0R1F(x)例如，已知負指數分布的密度函數為：例如，已知負指數分布的密度函數為：其分布函數為 令上式兩邊取對數后有 000)(xxexfx01)(0 xedtexFxxtuReexFRxx11)(uxuxlnln4.4.2

8、 4.4.2 近似法近似法 正態分布的密度函數為 分布函數為 當=0，=1時稱標準正態分布: xeaxfx, 021)(222)(xtdtexF222)(21)(2221)(xexfxtdtexF2221)( 任何非標準的正態分布隨機變量Y，通過下式都可以化為標準正態分布的隨機變量李雅普諾夫的中心極限定理 設 是相互獨立的隨機變量，它們的數學期望和方差為令 當滿足 yxn21,nkDEkkkk2 , 1,)(,)(2nknkkknkk11221,0)(lim212nkkkEE則隨機變量 的分布近似服從標準正態分布 如為 為0，1區間均勻隨機數，其期望 ， 方差 xnkknknkkkx1211n

9、21,2/ 12/ )()(abEk)1,0(12112)()(2baabD有 若已知一非標準正態，其期望值為 ，均方差為 ，按前式，則非標準的正態分布隨機變量為將標準正態分布的隨機變量代入后為 當n=5時 )2(121221212111111nkknkknknknkknnnnxxy)2(121nkknny) 3(261kky返回4.5 4.5 離散事件隨機系統的仿真方法離散事件隨機系統的仿真方法4.5.1 4.5.1 仿真時鐘仿真時鐘 仿真時鐘是隨著仿真的進程而不斷更新的時間機構。通常，在仿真開始時將仿真時鐘置為零。隨后，仿真時鐘不斷地給出仿真時間的當前值。仿真時間是仿真模型中的時間指示，代

10、表模型運行的真實時間。4.5.2 4.5.2 時間步長仿真法時間步長仿真法 一個系統的各種變量和參數往往隨著時間變化而變化，因此可把時間分成段（如1天、1小時、1分鐘或1秒鐘），仿真時鐘按時間段等距推進，每次推進需要掃描所有的活動，模擬各段時間內系統各種變量和參數的變化，這種仿真方法稱為時間步長法。 YN初始條件和數據增加一個時間步長在當前的時間步長內考察、分析、計算和記錄系統的活動。時鐘是否到達模 擬總時間?對模擬數據進行分析，輸出結果。圖3-5 時間步長仿真法 以井底煤倉系統為例構造仿真模型，分析煤倉容量變化情況. 煤煤倉箕斗 主井煤倉系統示意圖 根據實測數據統計分析，單位時間內（根據實測

11、數據統計分析，單位時間內（10min10min） 煤倉來煤量為符合正態分布的隨機變量，其均值為1=125t/10min，均方差為1=12.5t/10min；同時單位時間內箕斗提升次數也是隨機的，統計結果見表3-1，符合正態分布，每箕斗的提升量為10t。根據表3-1中的數據可求出其均值為2=46.5t/10min，均方差為2=9.63t/10min。現在要求模擬300天，求出： 每天最大煤倉容量和300天的平均值； 300天的平均日產量。 設10min為仿真時鐘的時間推進步長，每隔10min考察一次來煤量、煤倉容量提升量，構造仿真模型，其結構見圖3-7。 產生來煤量隨機數R1，產生箕斗提升隨機數

12、R2， 模擬天數D=0煤倉容量V=0，最大煤倉容量VD=0，日產量QD=0，仿真時鐘T=0。時鐘推進：T=T+1令 V=V+R1-10R2，Rd= Rd+ R1。判斷VD是否小于V，若小于V則令VD=V D=D+1T300d ?輸出煤倉容量每天最大值和300d平均值以及日產量平均值。NNY圖3-7 主井煤倉仿真結構4.5.3 4.5.3 事件步長仿真法事件步長仿真法 事件步長法是以事件的發生來推進模擬時鐘的增長。 初始狀態 選取最短時間事件 考察此事件所處的狀態 改變事件的狀態，并計算有關的數據。 累加事件的當前時間是否模擬完畢 ？ 輸出結果NY 圖3-8 事件步長仿真法汽車運輸系統汽車運輸系

13、統 已知一汽車隊共有N輛汽車，每輛車的載重量為Qt，每天工作7h。圖3-9為運輸系統示意圖。車庫和倉庫在一起，工作開始時汽車在倉庫前排隊裝車，每次一輛；裝完后汽車重載運行到火車站卸貨；卸車后空載運行回到倉庫再裝車。經統計知道，汽車運輸的裝車、重載運行、卸車、空載運行時間都是隨機的，符合正態分布，其均值和均方差見表3-2。因汽車運行一個周期要花一個多小時，因此汽車到達倉庫時間等于或超過6點半時就不再發車。 現在需要通過建立仿真模型，了解以下情況： 在汽車數N一定的情況下，每班可運輸多少噸貨物（QA）？ 汽車排隊總等待時間（F）； 汽車的利用率（FS）； 倉庫的空閑時間（A）及倉庫裝車的利用率（A

14、S）； 每輛車平均運行時間（Ti） 倉庫 3 N 2 1 汽車空載運行重載運行火車站卸車圖3-9 汽車運輸系統 設置初始狀態計算第一輪裝車及有關參數求最短時間事件（即最先回到倉庫的汽車當前時間及號碼）求裝車、卸車和空、重載運行的正態分布運行時間。汽車到達時間倉庫當前時間？計算汽車等待時間及有關參數計算卸車噸位數是否到收車時間？計算收車數目及運行時間汽車是否全部回到車庫？計算有關參數及輸出每天數據模擬天數是否完畢？輸出分析結果計算倉庫空閑時間及有關參數YNNNYYNY4.5.4 4.5.4 仿真結果的分析仿真結果的分析 仿真結果的平均值 4.5.5 誤差分析誤差分析 設用 作為未知參數的估計量，

15、兩者誤差不超過某一正數的概率為1-，即 區間 稱為置信區間 , 1-稱為置信概率，稱為顯著性水平（或置信度） nxxxxn211)(p),(由中心極限定理可知，當試驗次數足夠大時服從正態分布，由于2未知，需用樣本方差S2來表示，它服從t分布，即：對給定置信度，有nxU/) 1()(ntSnx1)1()(2/ntSnxp即x在置信度下的置信區間是例如，5次仿真結果：12500，12650，12600，12750。試求置信概率為0.99的置信區間。 ) 1(,) 1(2/2/nSntxnSntx125905151iixx45700)(11122niixxnS試驗次數n=5，n-1=41-=0.99

16、, /2=0.005 因此有：置信區間為（12344.19，12835.81） 6041. 4) 4() 1(005. 02/tnt81.128355457006041. 412590) 1(19.123445457006041. 412590) 1(2/2/nSntxnSntx 定置信區間（即允許波動范圍）及置信度，欲求重復仿真次數n值，即尋找最小的n，滿足 可利用試湊法，即令n=0，n=1，n=2，直到滿足上式為止。 例如，已知=0.1，=1，S=1.78，分別將n=8、9、10、11、12代入上式，可發現n=11為符合要求的最少仿真次數。 nSnt) 1(2/返回4.6 隨機分布的確定隨

17、機分布的確定4.6.1 分布函數的假設分布函數的假設樣本直方圖 有一組樣本值x1，x2，xn。樣本值中的最大值xmax和最小值xmin。選b略小于xmax，選a略大于xmin，將區間a，b等分為m個互不相交的子區間，即有下列分割點ti： a = t0 t1 t2 ti tm = b其中 mabttii1 樣本值落在區間ti-1，ti）中的個數，記作Vi ,則樣本落入本區間的頻率 以組距（ti-1，ti）為底邊，以 為高度繪制直方矩形,構成直方圖.畫一條曲線大致經過各個直方矩形的上邊，便得出近似的概率密度的函數。 mVfii1iiiittfyx0yti-1 ti4.6.2 參數估計參數估計 設總

18、體X的概率密度的形式f（x，）為已知，它只含一個未知參數。若x1，x2，xn為X的一組觀察值，則它們的聯合概率密度等于記作若存在的一個值 ，使得似然函數在時= 時則稱 是的極大似然估計。 ),(),(121niinxfxxxf),(),(121niinxfxxxLLmax),(21nxxxL 求極大: 令或例如，設X服從指數分布 0Ldd0lnLdd000)(xxexfxniniinxnxexxxLi1121)exp(),()(ln)(lnln1xnxnLniiniixnx11令 得0lnLdxdniixnx114.6.3 X2檢驗檢驗 用n個觀測值來檢驗總體分布的假設H0。 H0：總體X的分

19、布函數為F（x）。 將x軸分為k個不相交的區間（ai，ai+1），i =1，2，k,其中ai、ak+1可分別取-及+ ,這時有 pi = F（ai+1）- F（ai）= p（ai50）時，則不論總體屬什么分布，皮爾遜統計量總是近似服從自由度為k-1的分布，其中k是所分區間個數，是估計參數的個數。 于是，若在假設下算得的 ，則在水平下拒絕H0，否則就接受H0。) 1(22rkxx返回4.10 系統動力學系統動力學4.10.1 基本概念及特點 系統動力學（System Dynamics，簡稱SD）是一種研究復雜系統的方法1-4。該方法基于系統論，并融會了控制論、信息論、組織理論等多種理論，通過對系

20、統結構因果關系的分析，設計出反映系統行為的反饋回路，建立相應的方程式，最終建立系統動態模型，以此為基礎對系統進行仿真實驗。 特點:w 系統動力學問題是動態的問題，這些問題通常是用隨時間連續變化的量來表示的，如人口變化、資源變化等。w 應用系統動力學研究復雜系統，能夠容納大量變量，一般可達數千個以上。w 系統動力學模型，既有描述系統各要素之間因果關系的結構模型，以此來認識和把握系統結構，又有專門形式表現的數學模型，據此進行仿真實驗和計算，以掌握系統的未來動態行為。因此，系統動力學是一種定性分析和定量分析相結合的仿真技術。 w 在系統動力學模型中，能夠設定各種控制因素，當改變輸入的控制因素（例如不

21、同的組織狀態、經濟參數或不同的政策因素）時，可觀察系統的行為和變化的趨勢，對系統進行動態模擬。w 系統動力學通過模型進行仿真計算的結果，都用預測未來一定時期內各種變量隨時間而變化的曲線來表示，也就是說，系統動力學能處理高階次、非線性、多重反饋的復雜時變的復雜系統的有關問題。4.10.2 回路w 復雜系統總是由許多因果反饋回路偶合而成的。系統動力學就是使用反饋來揭示原因和尋找解決問題的辦法，因果關系是構成系統動力學模型的基礎。所謂反饋回路是指由兩個以上具有因果關系的變量,以因果關系彼此連接而成的閉合結構.為了建立系統動力學方程，通常用反饋決策回路表示系統各元素間的關系. 回路分正回路和負回路:

22、生產增加商品減少收入增加+-訂貨速度庫存差額庫存量-+-a-正反饋回路b-負反饋回路圖 2-15 正、負反饋回路+-+人口密度平均壽命死亡率人口生育率出生率+-+- 圖 2-16 復雜因果關系回路決策回路： w因果反饋回路表達了系統發生變化的原因，但這種定性描述還不能確定使回路中的變量發生變化的機制。為了建立系統動力學方程，可進一步用反饋決策回路表示系統各元素間的關系，其基本結構可用下圖表示。 信息水平(系統的表現水平)水平(系統狀態)決策源決策反饋回路基本結構4.10.3 系統動態行為 w系統動態行為是指系統特征變量隨時間變化的情況，典型情況：w曲線A表示的是簡單的反饋系統，變量隨時間增加而

23、增加，但增加速率隨變量值增加而逐漸減小，最后趨于某數值，屬負反饋。曲線B則為復雜系統的動態行為，隨時間而振蕩，然后趨于某一目標值，屬負反饋。曲線C表示系統變量隨時間作指數增長，屬正反饋。曲線D表示的動態行為開始一段時間如同曲線C，而后同曲線B呈振蕩狀，屬組合型系統動態行為。 4012312345BADC 動態行為曲線4.10.4 模型方程模型方程 w 系統是由許多反饋回路組合而成的，而每一個反饋回路模型則是通過一組微分方程來表示的，這就是模型方程。w 模型方程能反映系統行為隨時間改變的動態行為，其中最重要的兩個概念是：水平變量（也稱狀態變量）和速率變量（也稱決策變量）。所謂水平變量是指能表征系

24、統某種屬性的變量，一般它應該是一個積累量，如礦產資源地質儲量、庫存量等；而速率變量是指狀態變量變化的速度，在系統中描述的是物質的實際流動，如礦產資源的增長與消耗、庫存的入庫與出庫等。下下圖反映了系統動態行為過程: 水平方程。水平方程標識符為L，如 L LL.K=LL.J+（DT）（RA.JK-RS.JK） 速率方程。速率表示單位時間內水平的變化，標號為R。以訂貨速率為例： R OR.KL=（1/AT）（DI-IK） w 輔助方程。 輔助方程的標號為A。上式的速率方程中，DI要求的庫存量假如為一變量（等于要求庫存量平均銷售的周數WID與平均銷售率AS.RK的乘積），通過輔助方程獨立列出：w R

25、OR.KL=（1/AT）（DI.K-IK）w A DI.K=（WID）（AS.RK） J DT K DT LL1.LL2.LL1.KL2.KL1.jL2.jR1.JKR2.JK時間 圖2-19 系統動態行為過程描述返回系統動力學建模步驟 w 首先，確定系統的邊界，畫出因果圖；w 選擇模型的基本變量水平變量；w 以水平變量為中心構造各自的子系統；w 然后，根據因果圖，連接各子系統；w 根據以上的描述，寫出方程式；w 最后，進行仿真運算，進行仿真結果及政策分析。 礦區規劃系統動力學模型 從礦區煤炭生產系統的整體性出發，研究礦區發展的規模、礦區產量、投資、效益等彼此相互影響和制約的關系。通過對礦區生

26、產和自然資源內部結構和變量之間的因果關系、社會行為和發展過程的仿真實驗，對礦區未來的發展趨勢進行多種可行方案的預測。 礦區規劃的總模型 技術進步職工培訓固定資產生產效率投資工業產值職工總人數礦區生產能力施工工期自投資礦區人員礦井衰老物資消耗利潤留成礦區環境礦區總產量銷售利潤銷售總收入+-共有5個正反饋回路和2個負反饋回路。 正回路： 投資-礦區生產能力-產量-收入利潤-利潤留成-自投資-投資。 投資-固定資產-產值在企業收入中的貢獻-收入-利潤-利潤留成- 自投資-投資。 投資-職工培訓-生產效率-職工數減少-物質消耗-利潤-利潤留成- 自投資-投資。 投資-技術進步-生產效率-工期縮短-物質

27、消耗-利潤-利潤留成- 自投資-投資。 投資-技術進步-礦區環境-收入-利潤-利潤留成-自投資-投資。負回路： 礦區生產能力-產量-礦井衰老與報廢-礦區生產能力。 礦區生產能力-產量-礦區環境-收入-利潤-利潤留成-自投資-投資 -礦區生產能力。礦區總模型劃分為礦區總模型劃分為5個子模型個子模型:w礦區生產能力與產量子模型：從已有的礦井生產能力增加或減少、新建礦井的增加兩個方面考慮。 w礦區經濟子模型:礦區經濟子模型以固定資產和洗精煤產量的變化作為變量，主要反映礦區收入、利潤和積累的狀況，綜合體現礦區不同規劃方案的經濟效果。w礦區職工人數和效率子模型:礦區職工人數子模型通過工人數和效率反映礦區

28、生產發展、工資增長和需要進行生活福利設施建設等方面的情況。w礦區人口和住宅子模型:礦區人口和住宅子模型主要反映礦區人口發展狀況及對住房的要求。w礦區環境子模型:礦區環境子模型主要考慮礦區規劃中社會效益方面的內容，按照煤礦開采的特點，從矸石、廢水、廢氣、地表塌陷破壞等方面對礦區的環境狀況進行評價。礦區生產能力與產量子模型： DELAY3DELAY3量增加老井儲開采量老井老井可采儲量新井可采儲量量增加新井儲開采量新井計劃產量礦區總產量老井報廢需求差額改擴建投資新井投資礦區總生產能力老井能力增加老井能力減少改擴建開工率新井開工率新井能力老井能力圖2-16 礦區生產能力與產量模型流程圖系統模擬：w方案

29、：對8個生產礦井進行改擴建，并陸續新建3對礦井。w方案：對5個生產礦井進行改擴建，并陸續新建3對礦井。w方案：只對8個生產礦井進行改擴建，不建新礦井。w方案：維持現有狀況，既不改擴建也不新建礦井。w利用系統動力學模型對礦區30年內的生產能力、產量、投資、經濟效益、生產效率、職工人數、人口、生態環境等各項參數的動態發展和變化進行模擬計算，得到大量數據結果供決策層分析決策。 煤產量模擬結果比較煤產量模擬結果比較 方案礦區產量，萬噸各方案產量比較1990199520002005201020151990199520002005201020151882.42069.52134.72345.72476.4

30、2256.90000001822.41936.01925.22119.42257.22048.30-133.5-209.5-226.3-219.2-208.61822.42069.52016.42018.11868.81479.300-118.3-327.6-327.6-777.61822.41936.01806.91791.81649.61270.60-133.5-328.1-553.9-826.8-986.34.8 4.8 仿真語言仿真語言 開發仿真系統的專用語言,能夠提供離散事件仿真系統需要的各種功能：產生隨機數和具有特定分布性質的隨機變量；推進仿真時間；根據事件表確定下次事件，并把控制

31、轉移到相應的程序模塊；往事件表中增加記錄或從表中刪除記錄；收集和分析數據；生成仿真結果報表；檢查各種錯誤等。 具有較長發展歷史的仿真語言有GPSS（GPSSII、III、360、V）、GASP、和NGPSS以及SIMSCRIPTII，II-PL-US和SIMULA67等。 發展趨向：連續、離散混合仿真; 多種程序結構; 綜合化、一體化。 返回4.9 4.9 可視化仿真技術可視化仿真技術 可視化仿真的實質是采用圖形或圖象處理方式對仿真計算過程的跟蹤、控制和結果的后處理，同時實現仿真軟件界面的可視化。 OpenGL是包括Microsoft、SGI、Intergragh、IBM以及Intel等在內的ARB機構推出的3D圖形應用程序接口的工業標準。Ope