第4章蒙特卡羅仿真方法

1、航空航天中的計(jì)算方法航空航天中的計(jì)算方法授課教師：陳琪鋒授課教師：陳琪鋒中南大學(xué)航空航天學(xué)院中南大學(xué)航空航天學(xué)院第第4章章 蒙特卡羅蒙特卡羅仿真方法仿真方法航空航天中的計(jì)算方法Page 32021-12-31內(nèi)容提要內(nèi)容提要4.1蒙特卡羅法蒙特卡羅法的的基本概念基本概念4.2蒙特卡羅法蒙特卡羅法的的理論依據(jù)理論依據(jù)4.3隨機(jī)變量的隨機(jī)變量的抽樣抽樣4.4隨機(jī)模擬隨機(jī)模擬結(jié)果的統(tǒng)計(jì)分析結(jié)果的統(tǒng)計(jì)分析參考資料：參考資料：1 閆曉東閆曉東,許志許志. 飛行器系統(tǒng)仿真實(shí)訓(xùn)教程飛行器系統(tǒng)仿真實(shí)訓(xùn)教程, 西北工業(yè)大學(xué)出西北工業(yè)大學(xué)出版社，版社，2013. Chp. 6.2 David L. Darmofa

2、l, Computational Methods in Aerospace Engineering (Lecture Notes), MIT, 2005. Chp. 15, 16.3 Averill M. Law, 仿真仿真建模與建模與分析分析 (Fourth Edition), 清華清華大學(xué)出版社大學(xué)出版社, 2009. Chp. 1.5.航空航天中的計(jì)算方法Page 42021-12-314.1 蒙特卡羅法蒙特卡羅法的的基本基本概念概念問題問題 干擾和偏差客觀存在，無法消除：干擾和偏差客觀存在，無法消除：風(fēng)、推力誤差、結(jié)構(gòu)和裝配誤差風(fēng)、推力誤差、結(jié)構(gòu)和裝配誤差如何分析干擾和偏差對飛行過程的

3、影響？如何分析干擾和偏差對飛行過程的影響？隨機(jī)因素對系統(tǒng)運(yùn)動(dòng)影響的統(tǒng)計(jì)規(guī)律？隨機(jī)因素對系統(tǒng)運(yùn)動(dòng)影響的統(tǒng)計(jì)規(guī)律？蒙特卡羅（蒙特卡羅（Monte Carlo）法）法 通過隨機(jī)變量的統(tǒng)計(jì)試驗(yàn)（或隨機(jī)模擬），求解問題近通過隨機(jī)變量的統(tǒng)計(jì)試驗(yàn)（或隨機(jī)模擬），求解問題近似解的數(shù)值方法，又似解的數(shù)值方法，又稱統(tǒng)計(jì)試驗(yàn)法或隨機(jī)稱統(tǒng)計(jì)試驗(yàn)法或隨機(jī)模擬法。模擬法。4.1 蒙特卡羅法蒙特卡羅法的的基本概念基本概念航空航天中的計(jì)算方法Page 52021-12-31基本思想：基本思想：對隨機(jī)現(xiàn)象總體進(jìn)行抽樣，通過樣本對隨機(jī)現(xiàn)象總體進(jìn)行抽樣，通過樣本的的仿真和分析，仿真和分析，對結(jié)果進(jìn)行統(tǒng)計(jì)。對結(jié)果進(jìn)行統(tǒng)計(jì)?；具^程：

4、基本過程：對影響系統(tǒng)的隨機(jī)現(xiàn)象根據(jù)其概率分布產(chǎn)生隨機(jī)樣本對影響系統(tǒng)的隨機(jī)現(xiàn)象根據(jù)其概率分布產(chǎn)生隨機(jī)樣本（抽樣）；（抽樣）；用確定性方法對樣本進(jìn)行分析；用確定性方法對樣本進(jìn)行分析；對隨機(jī)樣本分析結(jié)果進(jìn)行統(tǒng)計(jì)，得到統(tǒng)計(jì)規(guī)律。對隨機(jī)樣本分析結(jié)果進(jìn)行統(tǒng)計(jì)，得到統(tǒng)計(jì)規(guī)律。4.1 蒙特卡羅法蒙特卡羅法的的基本概念基本概念航空航天中的計(jì)算方法Page 62021-12-31蒙特卡羅法的基本步驟蒙特卡羅法的基本步驟（1）建立系統(tǒng)數(shù)學(xué)模型；）建立系統(tǒng)數(shù)學(xué)模型；（2）確定系統(tǒng)中的各種隨機(jī)擾動(dòng)因素及概率分布；）確定系統(tǒng)中的各種隨機(jī)擾動(dòng)因素及概率分布；（3）產(chǎn)生各隨機(jī)擾動(dòng)變量的抽樣值；）產(chǎn)生各隨機(jī)擾動(dòng)變量的抽樣值；（

5、4）將隨機(jī)擾動(dòng)變量的抽樣值輸入系統(tǒng)數(shù)學(xué)模型，進(jìn)行）將隨機(jī)擾動(dòng)變量的抽樣值輸入系統(tǒng)數(shù)學(xué)模型，進(jìn)行確定性仿真計(jì)算，得到擾動(dòng)下系統(tǒng)運(yùn)動(dòng)過程；確定性仿真計(jì)算，得到擾動(dòng)下系統(tǒng)運(yùn)動(dòng)過程；（5）重復(fù)（）重復(fù)（4），進(jìn)行多次抽樣仿真，得到系統(tǒng)運(yùn)動(dòng)過），進(jìn)行多次抽樣仿真，得到系統(tǒng)運(yùn)動(dòng)過程子樣；程子樣；（6）對子樣仿真結(jié)果進(jìn)行處理，得到系統(tǒng)運(yùn)動(dòng)的統(tǒng)計(jì)特）對子樣仿真結(jié)果進(jìn)行處理，得到系統(tǒng)運(yùn)動(dòng)的統(tǒng)計(jì)特征量。征量。4.1 蒙特卡羅法蒙特卡羅法的的基本概念基本概念航空航天中的計(jì)算方法Page 72021-12-314.2 蒙特卡羅法蒙特卡羅法的理論依據(jù)的理論依據(jù)大數(shù)定理大數(shù)定理 如果隨機(jī)變量如果隨機(jī)變量X1，X2，XN獨(dú)

6、立同分布，且具有有限獨(dú)立同分布，且具有有限期望值期望值E(X)，則：，則：即隨機(jī)變量即隨機(jī)變量X的簡單子樣的算術(shù)平均值，當(dāng)子樣數(shù)的簡單子樣的算術(shù)平均值，當(dāng)子樣數(shù)N充分充分大時(shí)，以概論大時(shí)，以概論1收斂于它的期望值。收斂于它的期望值。子樣均值與數(shù)學(xué)期望的偏差？子樣均值與數(shù)學(xué)期望的偏差？4.2 蒙特卡羅法蒙特卡羅法的的理論依據(jù)理論依據(jù) lim()1NNPXE X11NNiiXXN 航空航天中的計(jì)算方法Page 82021-12-31方差：方差：得：得：4.2 蒙特卡羅法蒙特卡羅法的的理論依據(jù)理論依據(jù) 2221212212NEXEXXXXN 22NEXN NE X 航空航天中的計(jì)算方法Page 92

7、021-12-31中心極限定理中心極限定理 如果隨機(jī)變量如果隨機(jī)變量X1，X2，XN獨(dú)立同分布，且具有有限獨(dú)立同分布，且具有有限非零的方差非零的方差2，則當(dāng)，則當(dāng)N充分大時(shí)，有：充分大時(shí)，有：其中其中稱為置信度，稱為置信度，1- 稱為置信水平。稱為置信水平。蒙特卡羅方法（均值）的誤差：蒙特卡羅方法（均值）的誤差： 0.50.050.003 t0.67451.9634.2 蒙特卡羅法蒙特卡羅法的的理論依據(jù)理論依據(jù)2202()12ttNtPXE XedtN N 查標(biāo)準(zhǔn)正態(tài)查標(biāo)準(zhǔn)正態(tài)分布表確定分布表確定精度提高一個(gè)數(shù)量級，試驗(yàn)次數(shù)精度提高一個(gè)數(shù)量級，試驗(yàn)次數(shù)N需增加需增加2個(gè)數(shù)量級個(gè)數(shù)量級,NXNN

8、航空航天中的計(jì)算方法Page 102021-12-314.3 隨機(jī)變量隨機(jī)變量的的抽樣抽樣4.3.1 隨機(jī)數(shù)隨機(jī)數(shù)的的產(chǎn)生產(chǎn)生隨機(jī)數(shù)序列隨機(jī)數(shù)序列0,1上服從均勻分布的隨機(jī)變量的抽樣序上服從均勻分布的隨機(jī)變量的抽樣序列。列。 用專門的算法和程序可在計(jì)算機(jī)上生成用專門的算法和程序可在計(jì)算機(jī)上生成偽隨機(jī)數(shù)偽隨機(jī)數(shù)。 平方取中法、線性同余法、加同余法、平方同余法平方取中法、線性同余法、加同余法、平方同余法確定的算法、給定初值（種子），確定的序列！確定的算法、給定初值（種子），確定的序列！需對偽隨機(jī)數(shù)進(jìn)行需對偽隨機(jī)數(shù)進(jìn)行獨(dú)立性獨(dú)立性、均勻性等統(tǒng)計(jì)性質(zhì)的檢驗(yàn)、均勻性等統(tǒng)計(jì)性質(zhì)的檢驗(yàn)。周期長、統(tǒng)計(jì)性質(zhì)好的

9、偽隨機(jī)數(shù)，用于模擬試驗(yàn)造成的周期長、統(tǒng)計(jì)性質(zhì)好的偽隨機(jī)數(shù)，用于模擬試驗(yàn)造成的誤差不大。誤差不大。 4.3 隨機(jī)變量隨機(jī)變量的的抽樣抽樣航空航天中的計(jì)算方法Page 112021-12-314.3.2 離散隨機(jī)變量離散隨機(jī)變量的抽樣的抽樣 設(shè)設(shè)隨機(jī)變量隨機(jī)變量X，其概率分布為：，其概率分布為： 可按輪盤賭的方式可按輪盤賭的方式，由產(chǎn)生的隨機(jī)數(shù)序列進(jìn)行，由產(chǎn)生的隨機(jī)數(shù)序列進(jìn)行抽樣抽樣。 令：令： 依次產(chǎn)生隨機(jī)數(shù)依次產(chǎn)生隨機(jī)數(shù)ri，若，若 則第則第i次抽樣時(shí)隨機(jī)變量次抽樣時(shí)隨機(jī)變量X的抽樣值為的抽樣值為xk。 4.3 隨機(jī)變量隨機(jī)變量的的抽樣抽樣 1,1,2,1kknkkP Xxpknp (0)(

10、)10,1,2,kkjjPPpkn (1)( ),1,2,kkiPrPkn 航空航天中的計(jì)算方法Page 122021-12-314.3.3 連續(xù)隨機(jī)變量連續(xù)隨機(jī)變量的抽樣的抽樣直接法（反變換法）直接法（反變換法） 0,1上服從均勻分布的隨機(jī)變量上服從均勻分布的隨機(jī)變量R的分布函數(shù)：的分布函數(shù)： 連續(xù)隨機(jī)變量連續(xù)隨機(jī)變量X的分布函數(shù)的分布函數(shù)為為F(x)是是值域在值域在0,1上單調(diào)增上單調(diào)增的連續(xù)函數(shù)，利用的連續(xù)函數(shù)，利用ri得到得到X的抽樣值：的抽樣值： 函數(shù)函數(shù)F-1(R)可顯示表示時(shí)，才能應(yīng)用?？娠@示表示時(shí)，才能應(yīng)用。 4.3 隨機(jī)變量隨機(jī)變量的的抽樣抽樣0,0( ),011,1rG r

11、rrr 1( )iixFr 航空航天中的計(jì)算方法Page 132021-12-31例：指數(shù)分布隨機(jī)變量的生成例：指數(shù)分布隨機(jī)變量的生成 參數(shù)為參數(shù)為a的指數(shù)分布的概率密度函數(shù)：的指數(shù)分布的概率密度函數(shù)： 累積分布函數(shù)累積分布函數(shù)： 反函數(shù)：反函數(shù)： 用隨機(jī)數(shù)產(chǎn)生服從指數(shù)分布的隨機(jī)變量用隨機(jī)數(shù)產(chǎn)生服從指數(shù)分布的隨機(jī)變量X的抽樣：的抽樣：4.3 隨機(jī)變量隨機(jī)變量的的抽樣抽樣( ),0axf xaex 0( )( )1,0 xaxF xf t dtex 11( )ln(1),0,1Frrra 1ln ,0,1xrRUa 航空航天中的計(jì)算方法Page 142021-12-31函數(shù)變換函數(shù)變換法法 已知

12、已知隨機(jī)變量隨機(jī)變量 的聯(lián)合概率密度函數(shù)為的聯(lián)合概率密度函數(shù)為 。 隨機(jī)變量的變換隨機(jī)變量的變換 若逆變換存在且具有一階連續(xù)若逆變換存在且具有一階連續(xù)偏導(dǎo)數(shù)，即偏導(dǎo)數(shù)，即4.3 隨機(jī)變量隨機(jī)變量的的抽樣抽樣 1112221212,nnnnnYgXXXYgXXXYgXXX 12,nXXX 12,nfxxx 1111212212112,nnnnnXgY YYXgY YYXgY YY 航空航天中的計(jì)算方法Page 152021-12-31 則隨機(jī)變量則隨機(jī)變量 的聯(lián)合概率密度函數(shù)的聯(lián)合概率密度函數(shù)為：為：4.3 隨機(jī)變量隨機(jī)變量的的抽樣抽樣1111222212120nnnnnnXXXYYYXXXYY

13、YJXXXYYY 12,nY YY 11111221212,nnnnfgyyygyyygyyy航空航天中的計(jì)算方法Page 162021-12-31例：正態(tài)分布隨機(jī)變量的生成例：正態(tài)分布隨機(jī)變量的生成 設(shè)設(shè)X，Y獨(dú)立服從標(biāo)準(zhǔn)正態(tài)分布獨(dú)立服從標(biāo)準(zhǔn)正態(tài)分布： 變換：變換： 隨機(jī)變量隨機(jī)變量?和和的分布函數(shù)：的分布函數(shù)： 用反變換法得用反變換法得?和和的抽樣：的抽樣： 得到隨機(jī)變量得到隨機(jī)變量X，Y的抽樣：的抽樣：4.3 隨機(jī)變量隨機(jī)變量的的抽樣抽樣 2221( , )2xyf x ye cossinxy 22( )1,( )2FeF 122ln,2rr 22221( , )( ) ( )21( )

14、,( )2fefffef 12122lncos(2)2lnsin(2)xrryrr 航空航天中的計(jì)算方法Page 172021-12-314.4 隨機(jī)隨機(jī)模擬結(jié)果的統(tǒng)計(jì)分析模擬結(jié)果的統(tǒng)計(jì)分析隨機(jī)模擬試驗(yàn)隨機(jī)模擬試驗(yàn)的目的是獲得所研究的隨機(jī)現(xiàn)象的的目的是獲得所研究的隨機(jī)現(xiàn)象的統(tǒng)計(jì)特統(tǒng)計(jì)特性。性。 把所要研究的隨機(jī)現(xiàn)象用隨機(jī)變量把所要研究的隨機(jī)現(xiàn)象用隨機(jī)變量X表示，表示，N次隨機(jī)模擬次隨機(jī)模擬試驗(yàn)試驗(yàn)得到的結(jié)果（抽樣）序列為得到的結(jié)果（抽樣）序列為 。 X的數(shù)學(xué)期望的數(shù)學(xué)期望和方差分別用樣本均值和樣本方差和方差分別用樣本均值和樣本方差來統(tǒng)計(jì)：來統(tǒng)計(jì)：4.4 隨機(jī)隨機(jī)模擬結(jié)果的統(tǒng)計(jì)分析模擬結(jié)果的統(tǒng)計(jì)

15、分析12,Nxxx 1221111NiiNiixxNSxxN 航空航天中的計(jì)算方法Page 182021-12-31方差方差的的估值為估值為： 數(shù)學(xué)期望的置信區(qū)間（置信水平數(shù)學(xué)期望的置信區(qū)間（置信水平為為1- ）為：為： 事件概率的估計(jì)為事件概率的估計(jì)為： 其中其中m為為N次試驗(yàn)中事件出現(xiàn)的次數(shù)。次試驗(yàn)中事件出現(xiàn)的次數(shù)。 4.4 隨機(jī)隨機(jī)模擬結(jié)果的統(tǒng)計(jì)分析模擬結(jié)果的統(tǒng)計(jì)分析 22211(1)NxiiSxxNN N ,SSxtxtNNxmpN 航空航天中的計(jì)算方法Page 192021-12-31作業(yè)：作業(yè)：考慮炮彈在二維鉛錘平面內(nèi)的運(yùn)動(dòng)，假設(shè)地面考慮炮彈在二維鉛錘平面內(nèi)的運(yùn)動(dòng)，假設(shè)地面為一為一水平水平面，重力加速度為面，重力加速度為9.8m/s2，大氣密度為大氣密度為1.29kg/m3。 炮彈炮彈質(zhì)量質(zhì)量30kg，阻力系數(shù)為，阻力系數(shù)為0.1，迎風(fēng)截面積為，迎風(fēng)截面積為0.05m2。初速。初速度為度為95,105m/s區(qū)間的均勻分布，發(fā)射角（初速度方向與區(qū)間的均勻分布，發(fā)射角（初速度方向與水平面的夾角）服從均值為水平面的夾角）服從均值為40deg、均方差為、均方差為0