計算機仿真原理及應用第三講 離散事件系統(tǒng)仿真方法ppt課件

單位 物理電子學院主講人 王剛Email buncan wang 計算機仿真原理與應用 3離散事件系統(tǒng)仿真方法 離散事件系統(tǒng)示例例 單人理發(fā)店系統(tǒng) 設上午9 00開門 下午5 00關(guān)門 顧客到達時間一般是隨機的 為每個顧客服務的時間長度也是隨機的 系統(tǒng)的狀態(tài) 服務臺的狀態(tài) 忙或閑 顧客排隊等待的隊列長度 狀態(tài)量的變化也只能在離散的隨機時間點上發(fā)生 系統(tǒng) 一些具有特定功能 相互之間以一定的規(guī)律聯(lián)系著的物體所組成的總體 系統(tǒng)邊界 為了限制所研究問題涉及的范圍 用系統(tǒng)邊界把所研究的系統(tǒng)與影響系統(tǒng)的環(huán)境區(qū)分開來 實體 系統(tǒng)的對象 系統(tǒng)的組成元素都可以稱為實體 是仿真系統(tǒng)中可單獨識別和刻劃的構(gòu)成要素 屬性 屬性反映實體的某些性質(zhì) 是實體特征的描述 用特征參數(shù)或變量表示 它可以是文字型 數(shù)字型或邏輯型 活動 實體在一段時間內(nèi)持續(xù)進行的操作或過程 通常用于表示兩個可以區(qū)分的事件之間的過程 標志著系統(tǒng)狀態(tài)的轉(zhuǎn)移 狀態(tài) 系統(tǒng)的狀態(tài)是指在某一時刻實體及其屬性值的集合 事件 事件是引起離散事件系統(tǒng)狀態(tài)發(fā)生變化的行為 進程 由若干個有序事件及若干有序活動組成 描述了它所包括的事件及活動間的相互邏輯關(guān)系及時序關(guān)系 專用術(shù)語 仿真策略 要將系統(tǒng)模型轉(zhuǎn)換為計算機模型 首先要從總體上確定仿真模型的控制邏輯和仿真時鐘推進機制 即確定仿真策略 仿真策略是仿真模型的核心 反映了仿真模型的本質(zhì) 從根本上決定了仿真模型的結(jié)構(gòu) 迄今為止 離散事件系統(tǒng)形成了三種基本仿真策略 事件調(diào)度法 eventschedule ES 活動掃描法 activityscanning AS 進程交互法 ProcessInteractive PI 事件調(diào)度法由蘭德公司在1963年提出 在美國廣泛采用 歐洲不很流行 基本思想 將事件例程作為仿真模型的基本模型單元 按照事件發(fā)生的先后順序不斷執(zhí)行相應的事件例程 每一個有確定發(fā)生時間的事件 都有一個事件例程 用事件例程來處理事件發(fā)生后對實體狀態(tài)所產(chǎn)生的影響 并安排后續(xù)事件 事件調(diào)度法用事件的觀點分析真實系統(tǒng) 通過定義事件及每個事件的發(fā)生 引起系統(tǒng)狀態(tài)的變化 按時間順序 在每個事件發(fā)生時 確定并執(zhí)行有關(guān)的邏輯關(guān)系 按這種策略建立模型時 所有的事件均放在事件表中 模型中設有一個時間控制器 從事件表中選擇最早發(fā)生的事件并將仿真時鐘改到該事件發(fā)生的時間 然后調(diào)用與該事件相應的事件處理模塊 這樣 事件的選擇與處理不斷地進行 直到仿真終止或程序事件產(chǎn)生為止 事件調(diào)度法的仿真策略 1 初始化置仿真的開始時間t0和結(jié)束時間tf 置實體的初始狀態(tài) 置初始事件及其發(fā)生時間ts 2 仿真時鐘TIME ts 3 確定在當前時鐘TIME下發(fā)生的事件類型Ei i 1 2 n 并按規(guī)則排序 4 如果TIME tf 執(zhí)行 caseEiofE1 執(zhí)行E1的事件例程 產(chǎn)生后續(xù)事件類型及發(fā)生時間 En 執(zhí)行En的事件例程 產(chǎn)生后續(xù)事件類型及發(fā)生時間 endcase 否則 轉(zhuǎn) 6 事件調(diào)度法的仿真策略 5 將仿真時鐘TIME推進到下一個最早事件發(fā)生時間 這一步體現(xiàn)了仿真時鐘的推進機制 是將仿真時鐘推進到下一個最早事件的發(fā)生時刻 它與連續(xù)系統(tǒng)仿真中的時間推進方法 固定時間增量法不同 反映了離散事件系統(tǒng)狀態(tài)僅在離散時刻點上發(fā)生變化的特點 這種時間推進方法在離散事件系統(tǒng)仿真中普遍采用 稱為下一事件增量法 簡稱事件增量法 6 仿真結(jié)束 1 系統(tǒng)建模由觀測數(shù)據(jù)確定隨機變量的分布和參數(shù) 一般可用流程圖或網(wǎng)格圖的方式描述 反映臨時實體在系統(tǒng)內(nèi)部歷經(jīng)的過程 永久實體對臨時實體的作用以及它們相互之間的邏輯關(guān)系 2 確定仿真算法兩個方面內(nèi)容 如何產(chǎn)生所需求的隨機變量 采用什么方法對離散事件系統(tǒng)仿真 確定仿真策略 3 建立仿真模型仿真時鐘在各種算法中的推進方法 根據(jù)仿真算法建立仿真系統(tǒng)的計算機模型 變量定義及程序流程 4 設計仿真程序?qū)崿F(xiàn)仿真模型 5 仿真結(jié)果分析每次仿真結(jié)果只是隨機變量的一次取樣 仿真結(jié)果的可信度 離散事件系統(tǒng)仿真步驟 已知的基本信息 等待區(qū)足夠大 排隊規(guī)則先進先出FIFO 到達間隔服從負指數(shù)分布 1 1 10 臺 天 修理時間服從負指數(shù)分布 2 1 15 臺 天 仿真時間長度為365天 編程序求解 機器的平均等待時間 機器的平均逗留時間 修理臺利用率 實例 機器修理車間的仿真 這是一個典型的單服務員單隊列的排隊系統(tǒng)仿真模型 這類排隊系統(tǒng)主要包括兩個要素 顧客 即服務對象 和服務員 即服務設備 該系統(tǒng)由到達模式 服務模式 并行服務員數(shù)目 系統(tǒng)容量 排隊規(guī)則來表示 由命題可知 被修理的機器為 顧客 而修理臺為 服務員 該排隊系統(tǒng)的到達模式用機器到達間隔時間的負指數(shù)分布表示 服務模式由修理時間的負指數(shù)分布表示 系統(tǒng)中并行服務員數(shù)目為1 系統(tǒng)容量足夠大 排隊規(guī)則采用先進先出FIFO方式 系統(tǒng)描述 系統(tǒng)建模 采用事件調(diào)度法 具體的仿真步驟如下 1 初始化 給出當前仿真時鐘 系統(tǒng)狀態(tài)量及統(tǒng)計量的初始值 2 掃描事件表 將當前仿真時鐘增加到下一個最早發(fā)生事件的時間上 3 處理該事件 相應地改變系統(tǒng)狀態(tài) 4 收集統(tǒng)計數(shù)據(jù) 5 若仿真時間未結(jié)束 則返回2 否則 執(zhí)行下一步 6 分析收集的統(tǒng)計數(shù)據(jù) 產(chǎn)生報告 通過分析可知 該仿真模型只存在兩類事件 第一類事件為 到達事件 第二類事件為 離開事件 那么下一事件的類型由變量EVTFLAG給出 仿真模型的總體結(jié)構(gòu)圖如下頁所示 其中INIT為系統(tǒng)初始化子程序 TIMEDV為時間推進子程序 ARRIVE為到達事件處理子程序 DEPART為離開事件處理子程序 REPORT為報告生成子程序 仿真模型的總體結(jié)構(gòu)圖 建模變量說明 到達事件的處理流程 離開事件的處理流程 計算機仿真結(jié)果 由已知條件可知 到達間隔時間服從 1 1 10 臺 天 的負指數(shù)分布 修理時間服從 2 1 15 臺 天 的負指數(shù)分布 仿真時間長度為365天 故到達間隔時間均值EATI 1 1 10 天 修理時間均值ERT 1 2 15 天 仿真結(jié)束時間TIME 365 天 給定隨機數(shù)發(fā)生器種子SEED 113 則仿真結(jié)果機器在系統(tǒng)中平均逗留時間為33天 機器在隊列中平均等待時間為40天 修理臺利用率為78 9 基本思想 針對待求問題 根據(jù)物理現(xiàn)象本身的統(tǒng)計規(guī)律 或人為構(gòu)造一合適的依賴隨機變量的概率模型 使某些隨機變量的統(tǒng)計量為待求問題的解 進行大統(tǒng)計量 N 的統(tǒng)計實驗方法或計算機隨機模擬方法 理論依據(jù) 大數(shù)定理 均勻分布的算術(shù)平均收斂于真值中心極限定理 置信水平下的統(tǒng)計誤差 兩個例子 Buffen投針實驗求 射擊問題 打靶游戲 蒙特卡羅方法概述 基本思想 Buffon投針實驗 1777年 求 N 大數(shù)定理 蒙特卡羅方法概述 基本思想 針與平行線垂直方向夾角為 則相交概率為 一些人進行了實驗 其結(jié)果列于下表 1 設r表示射擊運動員彈著點到靶心的距離 r 表示擊中r處相應的得分數(shù) 環(huán)數(shù) f r 為該運動員彈著點的分布密度函數(shù) 它們反映運動員的射擊水平 該運動員的射擊成績?yōu)?用概率語言來說 是隨機變量 r 的數(shù)學期望 即 假設該運動員進行了 次射擊 每次射擊的彈著點依次為r1 r2 rN 則 次得分g r1 g r2 g rN 的算術(shù)平均值 代表了該運動員的成績 射擊問題 打靶游戲 用 次試驗所得成績的算術(shù)平均值作為數(shù)學期望的估計值 例如 設射擊運動員的彈著點分布為 用計算機作隨機試驗 射擊 的方法為 選取一個隨機數(shù) 按右邊所列方法判斷得到成績 這樣 就進行了一次隨機試驗 射擊 得到了一次成績 r 作 次試驗后 得到該運動員射擊成績的近似值 收斂性 大數(shù)定理 作為所求解的近似值 由大數(shù)定律可知 如果X1 X2 XN獨立同分布 且具有有限期望值 E X 則 即隨機變量X的簡單子樣的算術(shù)平均值 當子樣數(shù) 充分大時 以概率1收斂于它的期望值E X 由前面介紹可知 蒙特卡羅方法是由隨機變量X的簡單子樣X1 X2 XN的算術(shù)平均值 f x 是X的分布密度函數(shù) 則當N充分大時 有如下的近似式 蒙特卡羅方法的近似值與真值的誤差問題 概率論的中心極限定理給出了答案 該定理指出 如果隨機變量序列X1 X2 XN獨立同分布 且具有有限非零的方差 2 即 其中 稱為置信度 1 稱為置信水平 這表明 不等式 近似地以概率1 成立 且誤差收斂速度的階為 O N 1 2 蒙特卡洛法中的誤差 中心極限定理 上式中 與置信度 是一一對應的 根據(jù)問題的要求確定出置信水平后 查標準正態(tài)分布表 就可以確定出 通常 蒙特卡羅方法的誤差 定義為 給出幾個常用的 與 的數(shù)值 兩點說明 1 MC方法的誤差為概率誤差 這與其他數(shù)值計算方法是有區(qū)別的 2 誤差中的均方差 是未知的 必須使用其估計值來代替 在計算所求量的同時 可計算出估計值 2 減小估計的均方差 比如降低一半 那誤差就減小一半 這相當于N增大四倍的效果 減小方差的各種技巧 1 增大試驗次數(shù)N 在 固定的情況下 要把精度提高一個數(shù)量級 試驗次數(shù)N需增加兩個數(shù)量級 因此 單純增大N不是一個有效的辦法 顯然當給定置信度 后 誤差 由 和N決定 要減小 效率 一般來說 降低方差的技巧 往往會使觀察一個子樣的時間增加 在固定時間內(nèi) 使觀察的樣本數(shù)減少 所以 一種方法的優(yōu)劣 需要由方差和觀察一個子樣的費用 使用計算機的時間 兩者來衡量 這就是MC方法中效率的概念 它定義為 2c 其中c是觀察一個子樣的平均費用 顯然 2c越小 方法越有效 1 能夠比較逼真地描述具有隨機性質(zhì)的事物的特點及物理實驗過程 蒙特卡羅方法的優(yōu)點 從這個意義上講 蒙特卡羅方法可以部分代替物理實驗 甚至可以得到物理實驗難以得到的結(jié)果 用蒙特卡羅方法解決實際問題 可以直接從實際問題本身出發(fā) 而不從方程或數(shù)學表達式出發(fā) 它具有直觀 形象的特點 2 受幾何條件限制小 計算s維空間中的任一區(qū)域Ds上的積分 無論區(qū)域Ds的形狀多么特殊 只要能給出描述Ds的幾何特征的條件 就可以從Ds中均勻產(chǎn)生N個點 因此 在具有隨機性質(zhì)的問題中 如考慮的系統(tǒng)形狀很復雜 難以用一般數(shù)值方法求解 而使用蒙特卡羅方法 不會有原則上的困難 其中Ds為區(qū)域Ds的體積 這是數(shù)值方法難以作到的 得到積分的近似值 3 收斂速度與問題的維數(shù)無關(guān) 由誤差定義可知 在給定置信水平情況下 MC方法的誤差為O N 1 2 與問題本身的維數(shù)無關(guān) 維數(shù)的變化 只引起抽樣時間及估計量計算時間的變化 不影響誤差 這一特點 決定了蒙特卡羅方法對多維問題的適應性 而一般數(shù)值方法 比如計算定積分時 計算時間隨維數(shù)的冪次方而增加 而且由于分點數(shù)與維數(shù)的冪次方成正比 需占用相當數(shù)量的計算機內(nèi)存 這些都是一般數(shù)值方法計算高維積分時難以克服的問題 4 具有同時計算多個方案與多個未知量的能力 2 使用蒙特卡羅方法還可以同時得到若干個所求量 1 對于那些需要計算多個方案的問題 使用蒙特卡羅方法有時不需要像常規(guī)方法那樣逐個計算 而可以同時計算所有的方案 其全部計算量幾乎與計算一個方案的計算量相當 例如對于屏蔽層為均勻介質(zhì)的幾何平板 要計算若干種厚度的穿透概率時 只需計算最厚的一種情況 其他厚度的穿透概率在計算最厚一種情況時稍加處理便可同時得到 例如在模擬粒子過程中 可以同時得到不同區(qū)域的通量 能譜 角分布等 而不像常規(guī)方法那樣 需要逐一計算所求量 5 誤差容易確定 根據(jù)蒙特卡羅方法的誤差公式 可以在計算所求量的同時計算出誤差 6 程序結(jié)構(gòu)簡單 易于實現(xiàn) 在計算機上進行蒙特卡羅方法計算時 程序結(jié)構(gòu)簡單 分塊性強 易于實現(xiàn) 1 收斂速度慢 2 誤差具有概率性 蒙特卡羅方法的收斂速度為O N 1 2 一般不容易得到精確度較高的近似結(jié)果 對于維數(shù)少 三維以下 的問題 不如其他方法好 由于蒙特卡羅方法的誤差是在一定置信水平下估計的 所以它的誤差具有概率性 而不是一般意義下的誤差 1蒙特卡羅方法概述 MC缺點 3 計算結(jié)果與系統(tǒng)大小有關(guān) 例如在粒子輸運問題中 經(jīng)驗表明 只有當系統(tǒng)的大小與粒子的平均自由程可以相比較時 一般在十個平均自由程左右 蒙特卡羅方法計算的結(jié)果較為滿意 但對于大系統(tǒng)或小概率事件的計算問題 計算結(jié)果往往比真值偏低 在使用蒙特卡羅方法時 可以考慮把蒙特卡羅方法與解析 或數(shù)值 方法相結(jié)合 取長補短 既能解決解析 或數(shù)值 方法難以解決的問題 也可以解決單純使用蒙特卡羅方法難以解決的問題 真隨機數(shù) 不可重復 物理方法產(chǎn)生 隨機數(shù)是實現(xiàn)由已知分布抽樣的基本量 在由已知分布的抽樣過程中 將隨機數(shù)作為已知量 用適當?shù)臄?shù)學方法可以由它產(chǎn)生具有任意已知分布的簡單子樣 由具有已知分布的總體中抽取簡單子樣 在蒙特卡羅方法中占有非常重要的地位 總體和子樣的關(guān)系 屬于一般和個別的關(guān)系 或者說屬于共性和個性的關(guān)系 隨機數(shù)的產(chǎn)生和檢驗 蒙特卡羅模擬就是邊產(chǎn)生隨機數(shù)邊進行隨機模擬的方法 準隨機數(shù) 不具隨機性質(zhì) 只要處理問題能得到正確結(jié)果 如放射性衰變 電子設備的熱噪音 宇宙射線的觸發(fā)時間等 偽隨機數(shù) 可重復 數(shù)學方法產(chǎn)生 必須通過統(tǒng)計檢驗 區(qū)分 數(shù)列的隨機性和隨機數(shù)的分布 一個完美的隨機數(shù)序列可能具有某種分布 如均勻分布 高斯分布等 但具有某種分布的數(shù)列卻可能完全不是隨機的 分布函數(shù)為 最簡單 最基本 也是最重要的隨機數(shù)是在單一區(qū)間 0 1 上的均勻分布的隨機數(shù) 其分布密度函數(shù)為 由于隨機數(shù)在蒙特卡羅方法中占有極其重要的位置 我們用專門的符號 表示 用 1 2 代表相互獨立且具有相同單位均勻分布的隨機數(shù)序列 獨立性 均勻性是隨機數(shù)必備的兩個特點 如 擲篩子游戲 投擲硬幣 用來作為隨機數(shù)發(fā)生器的物理源主要有兩種 一種是根據(jù)放射性物質(zhì)的放射性 另一種是利用計算機的固有噪聲 用物理方法產(chǎn)生的隨機數(shù)序列無法重復實現(xiàn) 不能進行程序復算 給驗證結(jié)果帶來很大困難 而且 需要增加隨機數(shù)發(fā)生器和電路聯(lián)系等附加設備 費用昂貴 因此 該方法也不適合在計算機上使用 在計算機上用物理方法產(chǎn)生隨機數(shù)的基本原理是 利用某些物理現(xiàn)象 在計算機上增加某些特殊設備 可以在計算機上直接產(chǎn)生隨機數(shù) 這些特殊設備稱為隨機數(shù)發(fā)生器 1 馮 諾伊曼平方取中法 遞推公式 以十進制數(shù)為例 平方取中法是把一個2S位的十進制自平方后 去頭截尾只保留中間2S個數(shù)字 然后用102S來除 這樣就可以得到在 0 1 上均勻分布的偽隨機數(shù)序列 例如 設十進制數(shù)的2S 4 并取x1 6406 則有 相應的偽隨機數(shù)序列是0 6406 0 3680 0 1354 0 8333 0 4388等 具有周期性 有些數(shù)甚至會緊接著重復出現(xiàn) 很少使用 由Lehmer在1951年提出來的 它的一般形式是 對于任一初始值x0 偽隨機數(shù)序列由下面遞推公式確定 2 Lehmer線性同余法 例如 乘同余法 x0稱為種子 改變它的值就得到基本序列的不同區(qū)段隨機數(shù) a 乘子 c 增量 m 模 乘同余法具有在計算機上容易實現(xiàn) 快速等特點 所以乘同余法已被廣泛采用 偽隨機數(shù)的均勻性 偽隨機數(shù)的獨立性 均勻性是指在 0 1 區(qū)間內(nèi)等長度子區(qū)間中隨機數(shù)的數(shù)量是一樣的 按先后順序出現(xiàn)的隨機數(shù)中 每個隨機數(shù)的取值與其相距一定間隔的隨機數(shù)取值之間無關(guān) 判斷偽隨機數(shù)序列是否滿足均勻和相互獨立的要求 要靠統(tǒng)計檢驗的方法實現(xiàn) 對于偽隨機數(shù)的統(tǒng)計檢驗 一般包括兩大類 均勻性檢驗和獨立性檢驗 偽隨機數(shù)的均勻性 將區(qū)間 0 1 分為K個子區(qū)間 統(tǒng)計隨機數(shù)落在第k個子區(qū)間的實際頻數(shù)nk 它應當趨近于理論頻數(shù)mk 計算統(tǒng)計量 如果 2值很大 表示遠遠偏離理想值 因此要求 2值盡可能小 但如果它趨于0則有可能N已進入循環(huán) 通常求和中的每一項的大小約為1 因此 2的值約為K K的取值不能太大也不能太小 太大反映不出 小區(qū)間 的均勻性 太小反映不出 大區(qū)間 的均勻性 限制條件 1 概率論中的Pearson定理說明 上式的極限概率分布是 2分布 給出了 2 x的概率 整數(shù) 是系統(tǒng)自由度表示獨立測量的次數(shù) 由于存在一個限制條件 故 K 1 給出了 2 x的概率 余函數(shù) 因此 當給定顯著水平 或置信度1 后 由方程Q 2 或P 2 1 解出 值 或從已有的 2表中查得 值 如果由 1 式計算出來的 小于 則認為在此置信度下 偽隨機數(shù)序列在 0 1 中是均勻分布的 1 順序相關(guān)法 它用相鄰兩個隨機數(shù)的自相關(guān)函數(shù) 或相關(guān)系數(shù) 來標識偽隨機數(shù)序列的獨立性情況 間距為l的自相關(guān)函數(shù)是 相關(guān)系數(shù)越小 獨立性越好 2 多維頻率檢驗 1 將偽隨機數(shù)序列用任意一種辦法進行組合 每S個隨機數(shù)作為S維空間中的一個點的坐標值 于是可以構(gòu)成一個點序列 2 把S維空間中的單位方體分成為K個子方體 方體邊長 3 統(tǒng)計落在第k個子方體中的實際頻數(shù)nk 它應當趨近于理論頻數(shù) 隨機數(shù)的獨立性統(tǒng)計檢驗 隨機變量的抽樣 離散型 隨機變量X x1 x2 x3 xN 例如 x可取3個值x1 x2 x3 它們出現(xiàn)的幾率分別為2 8 5 8 1 8 則隨機數(shù)小于2 8時實現(xiàn)x1 在區(qū)間2 8 7 8中時實現(xiàn)x2 大于7 8時實現(xiàn)x3 概率密度f p1 p2 p3 pN 如果從 0 1 區(qū)間中均勻抽樣得到的隨機數(shù) 滿足下式時 則物理量x取值為xn 實際需要的大多數(shù)隨機變量并不是 0 1 區(qū)間均勻分布的 而是有各種不同形式分布密度函數(shù)的隨機變量 因此 對不均勻的隨機變量抽樣的關(guān)鍵問題是如何從均勻分布的偽隨機變量樣本中 抽取符合所要求的分布密度函數(shù)的簡單子樣 設連續(xù)型變量x在區(qū)間 a b 中取值 可視為將上述的離散情形取連續(xù)極限 要求變量x 可由上式解析反解出x 的函數(shù)表達式 即求反函數(shù) 這對一些簡單的幾率密度函數(shù)解析表達式是很容易做到的 則求反函數(shù)后得 隨機變量的抽樣 連續(xù)型 累積函數(shù) 注意 1 和 同樣服從 0 1 的均勻分布 一維 變換抽樣法的基本思想是將一個比較復雜的分布p x 的抽樣 變換為已知的簡單分布g y 的抽樣 我們希望找到x y之間的對應關(guān)系 使得幾率密度守恒 變換抽樣法 例如 黑體輻射的譜密度按頻率 表示時為 當希望將譜密度用波長 2 c 表示時 有 顯然 當g y 取 0 1 均勻分布時 問題即化為 尋找y x 使其導數(shù)為p x 然后在 0 1 區(qū)間中對變量y抽樣得到均勻分布的隨機數(shù) 再由x y 關(guān)系得到對應幾率密度函數(shù)p x 的隨機抽樣x 二維 有兩個變量x和y的聯(lián)合分布密度函數(shù)為p x y 欲變換至變量u和v 它們的聯(lián)合分布密度函數(shù)為g u v