基于蒙特卡羅模擬的航空發(fā)動(dòng)機(jī)故障風(fēng)險(xiǎn)預(yù)測(cè)

趙洪利，劉宇文

(1.中國(guó)民航大學(xué) 航空工程學(xué)院，天津300300;2.中國(guó)民航大學(xué) 中歐航空工程師學(xué)院，天津300300)

航空發(fā)動(dòng)機(jī)的風(fēng)險(xiǎn)水平指在規(guī)定的時(shí)間內(nèi)、正常飛行的條件下，航空發(fā)動(dòng)機(jī)出現(xiàn)故障的概率或次數(shù)［1］.故障風(fēng)險(xiǎn)預(yù)測(cè)方法能夠給出航空發(fā)動(dòng)機(jī)在未來(lái)一段時(shí)間周期內(nèi)可能發(fā)生故障的次數(shù).航空運(yùn)營(yíng)企業(yè)為了保證飛機(jī)運(yùn)行的可靠性，需要計(jì)劃好發(fā)動(dòng)機(jī)維修，并儲(chǔ)備相應(yīng)數(shù)量的需備件，這就需要對(duì)發(fā)動(dòng)機(jī)的故障風(fēng)險(xiǎn)水平做出預(yù)測(cè)，估計(jì)發(fā)動(dòng)機(jī)在未來(lái)一段時(shí)間內(nèi)各部件及發(fā)動(dòng)機(jī)系統(tǒng)出現(xiàn)故障的次數(shù)，從而為生產(chǎn)或運(yùn)營(yíng)計(jì)劃做出安排.

航空發(fā)動(dòng)機(jī)部件多，結(jié)構(gòu)復(fù)雜，不同部件都有各自的故障概率分布.經(jīng)過(guò)多年發(fā)展，航空發(fā)動(dòng)機(jī)部件故障概率的分析方法有很多［2-3］.其中，經(jīng)典的解析方法雖然能夠給出單個(gè)部件統(tǒng)計(jì)特性的精確表達(dá)式，但是用于確定多部件組成的復(fù)雜系統(tǒng)的故障風(fēng)險(xiǎn)等可靠性參數(shù)是十分困難的.而其他系統(tǒng)分析方法，如故障樹(shù)分析法、馬爾可夫狀態(tài)轉(zhuǎn)移法等，其問(wèn)題求解的規(guī)模往往會(huì)隨著部件數(shù)目的增多而呈指數(shù)量級(jí)增大［4-5］，且對(duì)所求解的系統(tǒng)往往有一定限制(如部件故障概率分布形式等).蒙特卡羅方法作為一種隨機(jī)模擬方法，通過(guò)對(duì)模型或系統(tǒng)的觀察或抽樣試驗(yàn)來(lái)計(jì)算所求參數(shù)的統(tǒng)計(jì)特征，能夠很好地避免以上問(wèn)題.

1 蒙特卡羅方法

蒙特卡羅方法的主要步驟［6］包括:針對(duì)具體問(wèn)題構(gòu)造相應(yīng)的概率過(guò)程;實(shí)現(xiàn)已知概率分布的隨機(jī)抽樣;利用隨機(jī)抽樣的結(jié)果計(jì)算具體的估計(jì)量.

1.1 構(gòu)造概率過(guò)程

構(gòu)造概率過(guò)程即找到具體問(wèn)題所服從的概率分布，如二項(xiàng)分布、指數(shù)分布、正態(tài)分布等.準(zhǔn)確描述問(wèn)題的概率過(guò)程是蒙特卡羅求解方法的關(guān)鍵，整個(gè)隨機(jī)模擬過(guò)程都是建立在已知的概率分布之上.如何確定問(wèn)題所服從的概率分布，這需要根據(jù)已有的數(shù)據(jù)并結(jié)合相應(yīng)的數(shù)據(jù)處理方法來(lái)確定.

1.2 隨機(jī)抽樣與統(tǒng)計(jì)

構(gòu)造完概率過(guò)程后，就知道了系統(tǒng)各個(gè)部分所服從的概率分布，如一輛汽車中儀表盤(pán)等電子元器件壽命可能服從指數(shù)分布，發(fā)動(dòng)機(jī)故障時(shí)間可能服從威布爾分布等.在已知概率分布的基礎(chǔ)上，產(chǎn)生滿足相應(yīng)分布的隨機(jī)變量即為隨機(jī)抽樣過(guò)程.隨機(jī)抽樣的核心在于隨機(jī)數(shù)的生成，常用的方法有工程上的物理方法和數(shù)學(xué)公式遞推法.工程上的物理方法無(wú)法重復(fù)便捷地生成隨機(jī)數(shù)，且價(jià)格昂貴.數(shù)學(xué)遞推公式法利用計(jì)算機(jī)就可重復(fù)產(chǎn)生大量的隨機(jī)數(shù)，雖然無(wú)法做到真正意義上的隨機(jī)性，但只要滿足相應(yīng)的隨機(jī)數(shù)檢驗(yàn)，即可滿足大部分要求［7］.

蒙特卡羅模擬方法能否得到好的結(jié)果關(guān)鍵在于整個(gè)隨機(jī)過(guò)程的構(gòu)造和計(jì)算過(guò)程中隨機(jī)數(shù)的性質(zhì).在隨機(jī)抽樣的基礎(chǔ)上，根據(jù)所求解的問(wèn)題，對(duì)隨機(jī)抽樣過(guò)程中產(chǎn)生的數(shù)據(jù)進(jìn)行記錄、統(tǒng)計(jì)，進(jìn)而確定問(wèn)題的估計(jì)量，得到問(wèn)題的解.

2 故障風(fēng)險(xiǎn)預(yù)測(cè)算法

2.1 故障概率模型

在航空發(fā)動(dòng)機(jī)故障數(shù)據(jù)處理中，威布爾分布是最常見(jiàn)的，且與數(shù)據(jù)符合程度相對(duì)較好的一種分布［8］.航空發(fā)動(dòng)機(jī)是一種高可靠性的產(chǎn)品，通常只有小樣本的故障數(shù)據(jù).威布爾分析在處理小樣本數(shù)據(jù)時(shí)，與其他分析方法相比通常都具有較好的效果［9-10］.因此，在航空發(fā)動(dòng)機(jī)故障風(fēng)險(xiǎn)預(yù)測(cè)中，本文采用了威布爾分布模型.

威布爾累積概率分布函數(shù)為

式中，F(xiàn)(t;β，η)為累積概率，t為故障時(shí)間，β 為形狀參數(shù)，η為尺度參數(shù);t0為起始點(diǎn).

對(duì)于小樣本航空發(fā)動(dòng)機(jī)故障數(shù)據(jù)，為了得到合適的尺度參數(shù)和形狀參數(shù)，采用中位秩回歸參數(shù)估計(jì)法［11］.

1)將故障時(shí)間數(shù)據(jù)T=［T1T2… Tn］T從小到大排列，并利用式(2)計(jì)算各個(gè)數(shù)據(jù)的中位秩，當(dāng)故障時(shí)間數(shù)據(jù)中存在刪失數(shù)據(jù)時(shí)，需利用式(3)調(diào)整數(shù)據(jù)排序值.

式中，i為數(shù)據(jù)序號(hào);Ri為第i個(gè)數(shù)據(jù)的中位秩;N為數(shù)據(jù)總數(shù).

式中，i′為調(diào)整后的排序值;j為前一個(gè)數(shù)據(jù)調(diào)整后的排序值;p為當(dāng)前刪失數(shù)據(jù)之后的數(shù)據(jù)個(gè)數(shù).

2)計(jì)算故障時(shí)間的自然對(duì)數(shù)，記為Y=ln T，并令

式中，I為單位矩陣;RY為Y的中位秩.

式中，xi，yi為向量 X，Y 的分量.

2.2 隨機(jī)變量的抽樣方法

滿足威布爾分布的隨機(jī)數(shù)生成算法有很多，反變換法是其中一種便捷、有效的方法［12-13］.設(shè)需產(chǎn)生分布函數(shù)F(x)的連續(xù)隨機(jī)數(shù)x，若已有［0，1］區(qū)間均勻隨機(jī)數(shù) r，則產(chǎn)生x的反變換公式［14］為

則威布爾的反變換公式為

式中，T為隨機(jī)故障時(shí)間.

反變換法的關(guān)鍵在于［0，1］區(qū)間上高品質(zhì)的均勻隨機(jī)數(shù).采用目前廣泛使用的乘同余組合發(fā)生器來(lái)產(chǎn)生［0，1］區(qū)間上的均勻隨機(jī)數(shù).遞推公式［15］為

式中，i=0，1，…;j=1，2，…，m;ri為［0，1］區(qū)間上隨機(jī)數(shù);mod為求余運(yùn)算符;m為組合數(shù).

選取組合數(shù)m=2，且令

2.3 蒙特卡羅模擬流程

設(shè)發(fā)動(dòng)機(jī)有n種故障模式，記為F1，F(xiàn)2，…，F(xiàn)n.這n種故障模式相互獨(dú)立，均服從威布爾分布，參數(shù)記為(η1，β1)，(η2，β2)，…，(ηn，βn)，故障時(shí)間抽樣為

式中，rk∈R，R為乘同余組合發(fā)生器生成的隨機(jī)數(shù)序列.

假定發(fā)動(dòng)機(jī)定期檢查時(shí)間為T，故障修復(fù)后，發(fā)動(dòng)機(jī)使用時(shí)間重新計(jì)，累積到時(shí)間T時(shí)，再做下次定檢，則蒙特卡羅模擬運(yùn)行流程如圖1所示.

圖1 蒙特卡羅模擬流程Fig.1 Monte Carlo simulation process

由圖1可知，模擬的基本步驟如下:

1)輸入機(jī)隊(duì)的原始數(shù)據(jù)，包括機(jī)隊(duì)總數(shù)，使用率，機(jī)隊(duì)年齡分布;

2)利用乘同余組合發(fā)生器構(gòu)建［0，1］區(qū)間均勻隨機(jī)數(shù)表;

3)針對(duì)每臺(tái)發(fā)動(dòng)機(jī)，從隨機(jī)數(shù)表中順序選取隨機(jī)數(shù)，利用式(13)計(jì)算各個(gè)故障模式的隨機(jī)故障時(shí)間 F1，F(xiàn)2，…，F(xiàn)n，首次故障時(shí)間需大于發(fā)動(dòng)機(jī)的已安全運(yùn)行時(shí)間，否則需重新產(chǎn)生一組隨機(jī)故障時(shí)間，直到均大于發(fā)動(dòng)機(jī)已安全運(yùn)行時(shí)間;

4)判斷 min{F1，F(xiàn)2，…，F(xiàn)n}是否小于定期檢查時(shí)間T，若是，則記錄該故障Fk，該故障模式發(fā)生次數(shù)累加1;

5)針對(duì)k故障模式，重新抽樣故障時(shí)間Fk;

6)重復(fù)步驟4)和步驟5)，直到min{F1，F(xiàn)2，…，F(xiàn)n}超過(guò)發(fā)動(dòng)機(jī)定期檢查時(shí)間T，則一次模擬結(jié)束;

7)遍歷機(jī)隊(duì)內(nèi)所有發(fā)動(dòng)機(jī)，每臺(tái)發(fā)動(dòng)機(jī)均進(jìn)行N次模擬，并記錄結(jié)果;

8)計(jì)算發(fā)動(dòng)機(jī)各個(gè)故障模式發(fā)生的概率;

9)結(jié)合機(jī)隊(duì)原始數(shù)據(jù)，計(jì)算機(jī)隊(duì)未來(lái)一段時(shí)間內(nèi)發(fā)動(dòng)機(jī)各個(gè)故障模式可能發(fā)生的次數(shù).

3 算例及結(jié)果分析

3.1 算例描述

為了驗(yàn)證該模型的準(zhǔn)確性與適用性，本文采用了某發(fā)動(dòng)機(jī)公司手冊(cè)中的數(shù)據(jù)作為輸入條件，利用上述故障風(fēng)險(xiǎn)預(yù)測(cè)算法進(jìn)行模擬仿真，并將結(jié)果與發(fā)動(dòng)機(jī)公司給出的軟件仿真結(jié)果進(jìn)行分析對(duì)比.手冊(cè)中提供了某噴氣式發(fā)動(dòng)機(jī)4種相互獨(dú)立的故障模式:F1表示過(guò)熱，F(xiàn)2表示葉片裂紋，F(xiàn)3表示油管裂紋，F(xiàn)4表示燃燒室裂紋.4種故障模式均服從不同參數(shù)的威布爾分布，具體參數(shù)見(jiàn)表1.同時(shí)，發(fā)動(dòng)機(jī)的定期檢查時(shí)間為1 000 h，整個(gè)機(jī)隊(duì)發(fā)動(dòng)機(jī)運(yùn)行時(shí)間分布如圖2所示.假定發(fā)動(dòng)機(jī)使用率為25 h/月.

表1 4種故障模式參數(shù)列表Table1 Parameters list of 4 failure modes

3.2 模擬流程

針對(duì)初始時(shí)間為0的發(fā)動(dòng)機(jī)，詳細(xì)闡述蒙特卡羅風(fēng)險(xiǎn)預(yù)測(cè)方法，模擬流程見(jiàn)圖3.

1)為4種模式生成隨機(jī)故障時(shí)間.利用隨機(jī)數(shù)表，順序選取前4個(gè)［0，1］區(qū)間上的隨機(jī)數(shù):0.007，0.028，0.517，0.603.根據(jù)式(13)，則有:

圖2 機(jī)隊(duì)發(fā)動(dòng)機(jī)運(yùn)行時(shí)間分布Fig.2 Operation time distribution of engine fleet

圖3 模擬流程Fig.3 Simulation process

2)4種故障模式中，最小的故障時(shí)間為951 h，并未達(dá)到定期檢查時(shí)間1000 h.

3)最先發(fā)生的為F1故障，記錄該故障發(fā)生時(shí)間為951 h，相當(dāng)于在未來(lái)第38個(gè)月發(fā)生故障.從隨機(jī)數(shù)表中選取下一個(gè)隨機(jī)數(shù)0.442，為F1故障重新生成隨機(jī)故障時(shí)間，F(xiàn)1=8827 h.此時(shí)，min{F1，F(xiàn)2，F(xiàn)3，F(xiàn)4}已超過(guò)定期檢查時(shí)間 1 000 h.至此，一次模擬結(jié)束.

針對(duì)每臺(tái)發(fā)動(dòng)機(jī)重復(fù)模擬N次，并遍歷所有發(fā)動(dòng)機(jī).根據(jù)大數(shù)定理，則發(fā)動(dòng)機(jī)故障發(fā)生概率近似為

式中，i=1，2，3，4 為故障模式序號(hào);j=1，2，…，12為發(fā)生故障的月份;Ni，j為第i種故障在j月份發(fā)生的次數(shù);Pi，j為第i種故障模式在j月份發(fā)生的概率;Na為額外抽樣次數(shù);Ne為發(fā)動(dòng)機(jī)數(shù)量.

則發(fā)動(dòng)機(jī)故障風(fēng)險(xiǎn)預(yù)測(cè)結(jié)果為

式中，F(xiàn)i，j為第 i種故障在 j月份發(fā)生的預(yù)測(cè)次數(shù).

發(fā)動(dòng)機(jī)一年內(nèi)的故障風(fēng)險(xiǎn)預(yù)測(cè)結(jié)果見(jiàn)表2，某發(fā)動(dòng)機(jī)公司提供的內(nèi)部軟件預(yù)測(cè)結(jié)果見(jiàn)表3，圖4給出了12月份發(fā)動(dòng)機(jī)故障次數(shù)的收斂結(jié)果.

表2 發(fā)動(dòng)機(jī)累積故障次數(shù)預(yù)測(cè)結(jié)果Table2 Forecasting results of cumulative failure number of engine

表3 某發(fā)動(dòng)機(jī)公司預(yù)測(cè)結(jié)果Table3 Forecasting results provided by an engine company

圖4 12月份發(fā)動(dòng)機(jī)故障次數(shù)收斂結(jié)果Fig.4 Convergence results of engine failure number on December

3.3 結(jié)果分析

在模擬過(guò)程中，為了得到穩(wěn)定的結(jié)果，需要實(shí)時(shí)觀察模擬輸出的數(shù)值，以便確定輸出結(jié)果是否穩(wěn)定.以12月份發(fā)動(dòng)機(jī)的故障次數(shù)預(yù)測(cè)結(jié)果為例，每10000次輸出一次模擬數(shù)值，并做100次抽樣，則總模擬次數(shù)為106.由圖4可知，最終結(jié)果已收斂.一般來(lái)說(shuō)，模擬次數(shù)越大，預(yù)測(cè)結(jié)果越逼近真實(shí)情況.在試驗(yàn)過(guò)程中，通過(guò)對(duì)比總模擬次數(shù)為106，107和108時(shí)的收斂圖形，發(fā)現(xiàn)結(jié)果差別不大.因此，總模擬次數(shù)設(shè)定為106已滿足收斂性要求.同時(shí)，由表2和表3可以看出，機(jī)隊(duì)內(nèi)發(fā)動(dòng)機(jī)一年內(nèi)各部件累積故障次數(shù)的預(yù)測(cè)結(jié)果與發(fā)動(dòng)機(jī)公司提供的數(shù)據(jù)相差不大，故障總次數(shù)僅相差0.97.預(yù)測(cè)結(jié)果準(zhǔn)確，驗(yàn)證了該算法的有效性與適用性.

4 結(jié)論

1)在確定了失效分布規(guī)律后，利用蒙特卡羅模擬算法進(jìn)行仿真，能夠比較準(zhǔn)確地估算出整個(gè)機(jī)隊(duì)發(fā)動(dòng)機(jī)在未來(lái)不同時(shí)間段內(nèi)的故障風(fēng)險(xiǎn)水平，從而為發(fā)動(dòng)機(jī)維修管理提供可靠性指標(biāo)參考;

2)蒙特卡羅模擬是多故障模式下風(fēng)險(xiǎn)預(yù)測(cè)的一個(gè)有效方法，并且在多故障模式下，該算法不但可以預(yù)測(cè)機(jī)隊(duì)整體風(fēng)險(xiǎn)水平，而且還可確定何種故障模式最易發(fā)生，從而可以有針對(duì)性地對(duì)該故障制定相應(yīng)的改進(jìn)措施，降低機(jī)隊(duì)的故障風(fēng)險(xiǎn);

3)當(dāng)某個(gè)易發(fā)故障得到改進(jìn)后，可重新進(jìn)行仿真模擬來(lái)找到下一個(gè)易發(fā)故障，不斷迭代，從而實(shí)現(xiàn)對(duì)機(jī)隊(duì)內(nèi)發(fā)動(dòng)機(jī)的動(dòng)態(tài)管理.

致謝 中國(guó)民航大學(xué)郭慶副教授提供了某發(fā)動(dòng)機(jī)廠商的技術(shù)資料，并為文章寫(xiě)作提出了建設(shè)性意見(jiàn)，在此表示感謝!

References)

［1］ 趙宇.可靠性數(shù)據(jù)分析［M］.北京:國(guó)防工業(yè)出版社，2011:14-20.Zhao Y.Data analysis of reliability［M］.Beijing:National Defense Industry Press，2011:14-20(in Chinese).

［2］ 楊宇航，馮允成.基于仿真的復(fù)雜系統(tǒng)可靠性、可用性和MTBF評(píng)估文獻(xiàn)綜述［J］.系統(tǒng)工程理論與實(shí)踐，2003，42(2):80-85.Yang Y H，F(xiàn)eng Y C.Survey of reliability and availability evaluation of complex system using Monte Carlo techniques［J］.Systems Engineering-Theory ＆ Practice，2003，42(2):80-85(in Chinese).

［3］ 劉曉平，唐益明，鄭利平.復(fù)雜系統(tǒng)與復(fù)雜系統(tǒng)仿真研究綜述［J］.系統(tǒng)仿真學(xué)報(bào)，2008，20(23):6303-6315.Liu X P，Tang Y M，Zheng L P.Survey of complex system and complex system simulation［J］.Journal of System Simulation，2008，20(23):6303-6315(in Chinese).

［4］ 邵偉.蒙特卡洛方法及在一些統(tǒng)計(jì)模型中的應(yīng)用［D］.濟(jì)南:山東大學(xué)，2012.Shao W.Monte Carlo methods theory and their applications in some statistical model［D］.Jinan:Shandong University，2012(in Chinese).

［5］ 袁明偉.復(fù)雜系統(tǒng)可靠性分析［D］.馬鞍山:安徽工業(yè)大學(xué)，2013.Yuan M W.Reliability analysis of complex system［D］.Maanshan:Anhui University of Technology，2013(in Chinese).

［6］ Dickman B H，Gilman M J.Technical note Monte Carlo optimization［J］.Journal of Optimization Theory and Applications，1989，60(1):149-157.

［7］ 周文彬.組合式偽隨機(jī)數(shù)發(fā)生器的研究與設(shè)計(jì)［D］.哈爾濱:哈爾濱工程大學(xué)，2013.Zhou W B.Design and research of the combined pseudo-random number generator［D］.Harbin:Harbin Engineering University，2013(in Chinese).

［8］ Baumshtein M V，Prokopenko A V，Ezhov V N.Probabilistic prediction of the fatigue life of gas-turbine engine compressor blades under two-level programmed loading［J］.Strength of Ma-terials，1985，17(5):587-592.

［9］ Nee A Y C，Song B，Ong S K.Re-engineering manufacturing for sustainability［M］.Singapore:Springer Singapore，2013:699-703.

［10］ Saralees N，Samuel K.Strength modeling using Weibull distributions［J］.Journal of Mechanical Science and Technology，2008，22(7):1247-1254.

［11］ 麻曉敏，張士杰，胡麗琴，等.可靠性數(shù)據(jù)威布爾分析中秩評(píng)定算法改進(jìn)研究［J］.核科學(xué)與工程，2007，27(2):152-155.Ma X M，Zhang S J，Hu L Q，et al.An improved rank assessment method for weibull analysis of reliability data［J］.Nuclear Science and Engineering，2007，27(2):152-155(in Chinese).

［12］ Yadolah D.The concise encyclopedia of statistics［M］.New York:Springer New York，2008:446-447.

［13］ Johnson G E.Constructions of particular random processes［J］.Proceedings of the IEEE，1994，82(2):270-285.

［14］ 金暢.蒙特卡羅方法中隨機(jī)數(shù)發(fā)生器和隨機(jī)抽樣方法的研究［D］.大連:大連理工大學(xué)，2005.Jin C.Study on random number generator and random sampling in Monte Carlo method［D］.Dalian:Dalian University of Technology，2005(in Chinese).

［15］ 楊自強(qiáng)，魏公毅.常見(jiàn)隨機(jī)數(shù)發(fā)生器的缺陷及組合隨機(jī)數(shù)發(fā)生器的理論與實(shí)踐［J］.數(shù)理統(tǒng)計(jì)與管理，2001，20(1):45-51.Yang Z Q，Wei G Y.Drawbacks in classical random number generators-theory and practice of combined generator［J］.Journal of Applied Statistics and Management，2001，20(1):45-51(in Chinese).