基于貝葉斯混合概率分布融合的系統(tǒng)可靠性分析與預(yù)測(cè)方法

楊樂(lè)昌, 郭艷玲

(1. 北京科技大學(xué)機(jī)械工程學(xué)院, 北京 100083;2. 北京航空航天大學(xué)自動(dòng)化科學(xué)與電氣工程學(xué)院, 北京 100191)

0 引 言

層次性結(jié)構(gòu)模型是一類具有普遍適用性的模型。基于“從整體到局部”的系統(tǒng)工程思想,大部分復(fù)雜機(jī)電設(shè)備(如衛(wèi)星),都可以視作多層次系統(tǒng),自上而下通常可劃分為系統(tǒng)層、子系統(tǒng)層,零部件層等,而其對(duì)應(yīng)的可靠性模型通常也具有顯著的層次性特征。在對(duì)這些可靠性模型進(jìn)行分析,做出可靠性預(yù)測(cè)的過(guò)程中,需要面對(duì)的問(wèn)題也與傳統(tǒng)可靠性模型有較大差別。

一方面,復(fù)雜機(jī)電設(shè)備通常包含眾多部件組件,作為底層的零部件多為標(biāo)準(zhǔn)件(如齒輪,軸等),通常已積累了較為充分的先驗(yàn)知識(shí)與實(shí)驗(yàn)數(shù)據(jù),且有歷史經(jīng)驗(yàn)數(shù)據(jù)可供參考,認(rèn)知較為完備;另一方面,對(duì)于系統(tǒng)整體而言,由于工況及實(shí)驗(yàn)條件所限,整體性實(shí)驗(yàn)一般代價(jià)高昂,開(kāi)展難度大,特別是在一些涉及航空航天、核工業(yè)的重大設(shè)備,由于客觀條件所限,幾乎不可能開(kāi)展整體性實(shí)驗(yàn)(如大型飛機(jī)整機(jī),空間站,核電站反應(yīng)堆等),因而系統(tǒng)整體極小樣本或無(wú)樣本,認(rèn)知極為匱乏。但在實(shí)際工程中,關(guān)注的重點(diǎn)恰恰是這些設(shè)備或系統(tǒng)整體的可靠性。

在本文中,將這種多層次信息分布不均衡(multi-level & information imbalanced, MLII)系統(tǒng)定義為MLII系統(tǒng)[1-2],圖1直觀地描述了MLII系統(tǒng)的特性。對(duì)于MLII系統(tǒng)的可靠性問(wèn)題,不能簡(jiǎn)單理解為傳統(tǒng)的“小樣本”(系統(tǒng)整體信息匱乏)或一般的“大數(shù)據(jù)”(底層零部件信息量大)問(wèn)題,而需針對(duì)其信息分布不均衡的特性,提出適用的系統(tǒng)可靠性分析方法。

圖1 MLII系統(tǒng)圖Fig.1 Graphical description of MLII system

1 問(wèn)題闡述

1.1 MLII系統(tǒng)特性分析

針對(duì)MLII系統(tǒng)的特性,一種自然而直觀的想法是充分利用系統(tǒng)底層完備的數(shù)據(jù)集,自下而上地“補(bǔ)償”頂層匱乏的信息,進(jìn)而對(duì)系統(tǒng)可靠性及其他物理性能參數(shù)做出準(zhǔn)確的分析與評(píng)估。貝葉斯方法作為基于條件概率推斷的一種統(tǒng)計(jì)方法,將主觀信息融入先驗(yàn)分布,以似然函數(shù)的形式利用客觀數(shù)據(jù),綜合了所有可用信息后作出概率推理結(jié)果,同時(shí)針對(duì)“小樣本”問(wèn)題同樣有效,可作為當(dāng)頂層信息匱乏時(shí)分析系統(tǒng)可靠性的一種手段[3-5],經(jīng)典貝葉斯推理方法的一般形式為

(1)

但傳統(tǒng)的貝葉斯方法在針對(duì)MLII系統(tǒng)時(shí)存在一些問(wèn)題。

(1) 在基于貝葉斯理論的統(tǒng)計(jì)推斷中,先驗(yàn)分布的選取對(duì)于后驗(yàn)分布的精度及最終可靠性分析結(jié)果的準(zhǔn)確性都起到至關(guān)重要的作用。在傳統(tǒng)的貝葉斯方法中,通常依據(jù)專家經(jīng)驗(yàn)給出適當(dāng)?shù)膮?shù)先驗(yàn)分布,但對(duì)于MLII系統(tǒng)而言,系統(tǒng)及子系統(tǒng)層通常信息匱乏,準(zhǔn)確的先驗(yàn)分布難以獲取,如果僅采用無(wú)信息先驗(yàn)分布,在小樣本的條件下,難以獲得足夠精度的可靠性分析結(jié)果[5]。對(duì)于具有MLII系統(tǒng)特性的復(fù)雜機(jī)電設(shè)備,需對(duì)經(jīng)典的貝葉斯推理算法做出改進(jìn),充分利用底層完備數(shù)據(jù)集的信息,給出具有足夠精度的先驗(yàn)分布,以提高分析結(jié)果精度。

(2) 復(fù)雜機(jī)電設(shè)備的MLII系統(tǒng),信息來(lái)源廣泛,數(shù)據(jù)多樣。即使對(duì)于同一事件(例如參數(shù)取值),也可能出現(xiàn)相異,甚至完全不同的認(rèn)知,即所謂的“信息沖突”問(wèn)題[6]。例如對(duì)于相同產(chǎn)品,不同專家根據(jù)各自經(jīng)驗(yàn)給出的故障率、平均故障間隔時(shí)間等參數(shù)的估測(cè)值可能有較大差異,可理解為可能存在的主觀信息沖突;又如,對(duì)速度、位移等物理量的動(dòng)力學(xué)分析結(jié)果可能與通過(guò)數(shù)理統(tǒng)計(jì)獲得的估測(cè)值存在一定偏差,可理解為可能存在的主客觀信息沖突,而經(jīng)典貝葉斯方法并無(wú)處理此類問(wèn)題的有效機(jī)制。

針對(duì)經(jīng)典貝葉斯方法在處理MLII系統(tǒng)可靠性問(wèn)題時(shí)的一些局限性[7-12],提出了基于貝葉斯與信息融合(Bayesian-based information extraction & aggregation, BIEA)的系統(tǒng)可靠性分析與預(yù)測(cè)方法。

1.2 BIEA方法理論框架

不失一般性,考慮不確定性傳遞模型Θ=M(θ)中的輸入?yún)?shù)向量θ和輸出參數(shù)向量Θ,分別具有獨(dú)立的直接先驗(yàn)分布(direct prior,DP)π(θ)和π(Θ)。此時(shí),由于參數(shù)向量θ具有獨(dú)立的先驗(yàn)分布,故參數(shù)向量Θ通過(guò)關(guān)系函數(shù)Θ=M(θ)會(huì)引入額外的概率分布π*(Θ)。由于信息來(lái)源不同,一般情況下,同一參數(shù)Θ的先驗(yàn)分布π(Θ)與π*(Θ)不同。

上述問(wèn)題的實(shí)質(zhì)在于如何綜合利用所有可用信息,融合概率分布π(Θ)與π*(Θ)。基于上述考量,定義概率分布融合運(yùn)算符“⊕”,提出BIEA方法,其基本形式為

(2)

相較于傳統(tǒng)的貝葉斯方法,BIEA方法通過(guò)構(gòu)建直接先驗(yàn)分布πD(Θ)與間接先驗(yàn)分布πI(Θ)的方式區(qū)分不同來(lái)源信息對(duì)參數(shù)先驗(yàn)分布的影響,定量計(jì)算不同信息源對(duì)參數(shù)先驗(yàn)概率分布的貢獻(xiàn)。通過(guò)融合先驗(yàn)分布πC(Θ),綜合所有可用信息,計(jì)算參數(shù)后驗(yàn)分布,進(jìn)而獲取所需的可靠性指標(biāo),即可靠度、平均故障間隔時(shí)間和故障率等。其中直接先驗(yàn)分布πD(Θ)一般基于主觀經(jīng)驗(yàn)獲得,而間接先驗(yàn)分布πI(Θ)通過(guò)分析模型傳遞后間接獲得。有關(guān)該方法與傳統(tǒng)貝葉斯方法的異同如圖2所示。

圖2 BIEA推理方法融合先驗(yàn)分布構(gòu)建示意圖Fig.2 Sketch map of BIEA reasoning method fusion priori distribution construction

2 基于自更新權(quán)重系數(shù)的BM概率分布融合方法

2.1 經(jīng)典BM方法

在BIEA方法中,構(gòu)建了融合先驗(yàn)分布以替代經(jīng)典貝葉斯方法中完全由主觀經(jīng)驗(yàn)給出的傳統(tǒng)先驗(yàn)分布,即

πC(Θ)=πD(Θ)⊕πI(Θ)

(3)

融合先驗(yàn)分布πC(Θ)由直接先驗(yàn)分布πD(Θ)和間接先驗(yàn)分布πI(Θ)共同構(gòu)成。其中直接先驗(yàn)分布基于研究對(duì)象的直接相關(guān)信息,如專家經(jīng)驗(yàn),相似產(chǎn)品的歷史數(shù)據(jù)等;間接先驗(yàn)分布則是對(duì)研究對(duì)象的間接相關(guān)信息進(jìn)行可靠性分析,并利用結(jié)構(gòu)函數(shù)間的不確定性傳遞關(guān)系后計(jì)算得到。由于信息來(lái)源與計(jì)算方式均不同,因此同一事件(或參數(shù))的直接先驗(yàn)分布和間接先驗(yàn)分布通常并不一致。一方面,對(duì)于研究對(duì)象本身是有直觀認(rèn)知的(直接先驗(yàn)分布);另一方面,系統(tǒng)或子系統(tǒng)層的參數(shù)與其父節(jié)點(diǎn)參數(shù)存在函數(shù)關(guān)系,而MLII系統(tǒng)的似然函數(shù)中也包含可靠性函數(shù)關(guān)系,這都相當(dāng)于添加了額外的約束,它表明通過(guò)結(jié)構(gòu)函數(shù)間的不確定性傳遞關(guān)系,上層單元引入了底層信息(額外約束)。而兩種不同的概率顯然無(wú)法直接相互疊加,故需引入一種具有數(shù)學(xué)一致性與完備性的概率融合方法。

貝葉斯混合(Bayesian melding, BM)是一種經(jīng)典的概率融合方法,該方法最早由文獻(xiàn)[13]提出,可有效地處理2種不同類型的概率分布,在各行各業(yè)均有應(yīng)用[14-15]。融合后的概率分布具有數(shù)學(xué)完備性,服從概率公理,且能夠繼承2種不同概率分布的統(tǒng)計(jì)特性。經(jīng)典BM方法通過(guò)Pooling的方法構(gòu)建,常見(jiàn)的線性融合方法為

πC(Θ)∝α·πD(Θ)+(1-α)·πI(Θ)

(4)

對(duì)數(shù)融合方法為

πC(Θ)∝πD(Θ)απI(Θ)(1-α)

(5)

式中,πC(Θ)、πD(Θ)和πI(Θ)分別為參數(shù)Θ的融合先驗(yàn)分布、直接先驗(yàn)分布和間接先驗(yàn)分布；權(quán)重系數(shù)α∈[0,1]用于平衡直接先驗(yàn)分布和間接先驗(yàn)分布對(duì)融合先驗(yàn)分布的貢獻(xiàn)量。

2.2 基于自更新權(quán)重系數(shù)的BM方法

在BM方法中,權(quán)重系數(shù)α的取值對(duì)最終融合先驗(yàn)分布的確定是至關(guān)重要的。文獻(xiàn)[16]認(rèn)為這種選取本質(zhì)上是靈活而非固定唯一的,并討論了關(guān)于權(quán)重系數(shù)α的取值。在沒(méi)有額外信息的情況下,同一模型(例如Θ=M(θ))中輸入?yún)?shù)(θ)和輸出信息(Θ)的可信程度是相同的,文獻(xiàn)[16]中例子將α設(shè)定為0.5,即直接先驗(yàn)分布和間接先驗(yàn)分布的置信水平是一致的。但事實(shí)上,均值0.5的選取方法并不是普遍適用的。文獻(xiàn)[17]嘗試不同的權(quán)重系數(shù)α的取值(0.2, 0.4, 0.5, 0.6, 0.8, 0.9),得到了不同的融合先驗(yàn)分布,但發(fā)現(xiàn)參數(shù)后驗(yàn)分布對(duì)α的取值并不敏感,究其原因在于該例子數(shù)據(jù)充足,樣本量較大,使得在推理計(jì)算過(guò)程中逐漸形成較強(qiáng)的似然函數(shù),足以修正初始先驗(yàn)分布中的偏差。事實(shí)上,該例子相較于先驗(yàn)分布,似然函數(shù)占據(jù)更為主導(dǎo)的位置。而直接先驗(yàn)分布選取無(wú)信息先驗(yàn)分布的做法對(duì)MLII系統(tǒng)并不適用,如前所述,對(duì)于MLII系統(tǒng)中的組件或部件而言,由于實(shí)際工況限制,往往數(shù)據(jù)匱乏,樣本量較小,似然函數(shù)在計(jì)算過(guò)程中無(wú)法占據(jù)主導(dǎo)位置;另一方面,對(duì)于底層部件而言,可能包含較多的專家經(jīng)驗(yàn)等主觀信息,因此先驗(yàn)分布的影響不可忽略。在這種情況下,合理選擇權(quán)重系數(shù)α就顯得尤為重要。

本論文針對(duì)MLII這一類存在信息不均衡特性的系統(tǒng),提出了一種基于自更新權(quán)重系數(shù)的BM方法。

對(duì)于線性情況,有

πC(Θ)∝πC(Θ|flinear(α,βD,βI))=

(6)

對(duì)于對(duì)數(shù)情況,有

πC(Θ)∝πC(Θ|flog(α,βD,βI))=

(7)

步驟1從參數(shù)θ的先驗(yàn)分布πD(Θ)中抽取I個(gè)樣本(θ1,θ2,…,θI)。從權(quán)重系數(shù)α的先驗(yàn)分布π(α)中抽取J個(gè)樣本 (α1,α2,…,αJ)。

步驟2對(duì)于抽取的每一個(gè)樣本θi,通過(guò)傳遞模型Θi=M(θi)計(jì)算對(duì)應(yīng)的樣本輸出值。

步驟3使用非參數(shù)估計(jì)方法(例如核密度估計(jì))計(jì)算間接先驗(yàn)分布πI(Θ)。

步驟4對(duì)每一組(θi,αj)的取值計(jì)算重要度抽樣概率為

步驟5以{wij:i=1,2,…,I;j=1,2,…,J}的概率從 (θi,αj)的離散分布中抽取L個(gè)樣本,即

π(θ,α)∝{wij:i=1,2,…,I;j=1,2,…,J}

步驟6由于輸入?yún)?shù)θ與權(quán)重系數(shù)α相互獨(dú)立,因此可以分別計(jì)算輸入?yún)?shù)θ的離散后驗(yàn)概率分布π(θ|D) 和權(quán)重系數(shù)α的離散后驗(yàn)概率分布π(α|D)，即

步驟7輸出參數(shù)Θ的離散后驗(yàn)概率分布π(Θ|D)可通過(guò)傳遞模型Θ(1,2,…,L)=M(θ1,2,…,L)估計(jì),即

本文所提出的基于自更新權(quán)重系數(shù)的BM方法不預(yù)先設(shè)定權(quán)重系數(shù)α的取值,而將其設(shè)置為未知超參數(shù),賦予其初始分布,再利用客觀數(shù)據(jù)對(duì)其取值修正,實(shí)時(shí)更新。權(quán)重系數(shù)α受到似然函數(shù)的影響,其取值會(huì)在推理過(guò)程隨著信息的積累而自動(dòng)變化調(diào)整,降低偏離實(shí)測(cè)數(shù)據(jù)的先驗(yàn)分布的權(quán)重,增加貼近實(shí)測(cè)數(shù)據(jù)的先驗(yàn)分布的權(quán)重,在反復(fù)迭代更新后,權(quán)重系數(shù)α達(dá)到穩(wěn)定值,此時(shí)的后驗(yàn)分布與實(shí)測(cè)數(shù)據(jù)吻合度最高。此方法可有效降低先驗(yàn)分布偏差對(duì)參數(shù)估計(jì)結(jié)果準(zhǔn)確性的影響,提高可靠性分析的精度。而提出的BM方法的另一個(gè)優(yōu)點(diǎn)則是,通過(guò)式(5)或式(6)得到融合先驗(yàn)分布是具有解析形式的。考慮到相同概率分布的直接疊加或乘積并不一定服從原始的概率分布形式，即式(3)和式(4)得到的融合先驗(yàn)分布通常并沒(méi)有解析形式,而這對(duì)于多層次的復(fù)雜貝葉斯模型,尤其是MLII系統(tǒng)的貝葉斯模型來(lái)說(shuō)很不方便。通過(guò)式(5)或式(6)得到的融合先驗(yàn)分布不僅具有明確的解析形式,還可以避免多峰概率分布的出現(xiàn),這也是在BM方法中經(jīng)常碰到的問(wèn)題。

現(xiàn)以簡(jiǎn)單的例子說(shuō)明上述問(wèn)題,分別選擇Beta(6,2)和Beta(2,6)作為某參數(shù)的直接先驗(yàn)分布和間接先驗(yàn)分布,通過(guò)經(jīng)典BM方法得到的融合先驗(yàn)分布和通過(guò)BIEA方法的BM在不同的權(quán)重系數(shù)α下的融合概率分布如圖3所示。

圖3 兩類BM方法對(duì)比Fig.3 Comparisons of the two kinds of BM approaches

可以發(fā)現(xiàn)，當(dāng)權(quán)重系數(shù)α取極限邊界值(0和1)時(shí),2種BM方法所得結(jié)果是相同的,這也意味著基于BIEA的BM方法與經(jīng)典貝葉斯方法具有相同融合邊界,即兼容性。但是,對(duì)于其他的權(quán)重系數(shù)α取值(0.2, 0.4,0.6,0.8),經(jīng)BM方法所得結(jié)果呈現(xiàn)出多峰值特征(雙峰特性),而B(niǎo)IEA的BM方法僅有一個(gè)峰值。在使用BIEA的BM后,參數(shù)后驗(yàn)分布為

p(Θ|β)∝L(Θ)πC(Θ|β,α)π(α)

(8)

通過(guò)馬爾可夫鏈蒙特卡羅(Markov chain Monto Carlo,MCMC)等抽樣方法進(jìn)行數(shù)值計(jì)算,此時(shí)可反解權(quán)重系數(shù)α的后驗(yàn)分布，即

p(α)=p(f-1(β))

(9)

通過(guò)應(yīng)用自更新權(quán)重系數(shù)的BM方法,可以更為精細(xì)地平衡基本先驗(yàn)的分布和間接先驗(yàn)分布對(duì)融合先驗(yàn)分布的貢獻(xiàn),而融合先驗(yàn)分布也會(huì)同時(shí)繼承2種先驗(yàn)分布的統(tǒng)計(jì)特征。該方法為權(quán)重系數(shù)α的定量量化提供了一種較為靈活的方法。該方法與經(jīng)典的BM方法有以下不同之處,如表1所示。

表1 經(jīng)典BM方法與自更新權(quán)重系數(shù)的BM方法的異同對(duì)比

3 MLII系統(tǒng)的BIEA方法

第2節(jié)給出了基于自更新權(quán)重系數(shù)的BM方法的數(shù)理算法,下面結(jié)合圖4所示的多層次結(jié)構(gòu)來(lái)闡述MLII系統(tǒng)的BIEA模型。

圖4 多層次結(jié)構(gòu)系統(tǒng)Fig.4 Multi-level structure system

3.1 MLII系統(tǒng)的BIEA模型

不失一般性,以l行第k個(gè)單元E(l,kl)為研究對(duì)象,其參數(shù)集θ(l,kl)。那么第l+1行的父節(jié)點(diǎn)E(l+1,kl+1)則有參數(shù)集θ(l+1,kl+1)。給定其參數(shù)直接先驗(yàn)分布πD(θ(l+1,kl+1)),其可靠度函數(shù)可一般地描述為R(l+1,kl+1)(t)=f(t|θ(l+1,kl+1))。其中，f()是由具體物理背景及失效機(jī)理確定的函數(shù)。那么,研究對(duì)象的可靠度函數(shù)及對(duì)應(yīng)的參數(shù)概率密度函數(shù)的一般形式為

R(l,kl)(t|θ(l,kl))=Ψ(l,kl)(R(l+1,kl+1)(t|θ(l+1,kl+1)):kl+1∈Q(l,kl))

f(l,kl)(t|θ(l,kl))=

(10)

式中,Ψ(l,kl)是由對(duì)象單元E(l,kl)及其父節(jié)點(diǎn)E(l+1,kl+1)確定的結(jié)構(gòu)函數(shù),Q(l,kl)是所有父節(jié)點(diǎn)的指標(biāo)集。對(duì)象單元參數(shù)的間接先驗(yàn)分布可通過(guò)式(14)的隨機(jī)變量轉(zhuǎn)換方法來(lái)計(jì)算，即

(11)

在給定參數(shù)直接先驗(yàn)分布πD(θ(l,kl))的基礎(chǔ)上,融合先驗(yàn)分布可由式(15)計(jì)算，即

πC(θ(l,kl))=πD(θ(l,kl))⊕πI(θ(l,kl))

(12)

融合先驗(yàn)分布中包含了底層數(shù)據(jù)與信息,并參與當(dāng)前層次的貝葉斯推理過(guò)程,重復(fù)這一過(guò)程,就可以將底層完備的信息逐漸傳遞至頂層,在系統(tǒng)級(jí)綜合所有可用信息,做出較為準(zhǔn)確的可靠性分析。

依據(jù)貝葉斯理論,在給定參數(shù)先驗(yàn)分布與似然函數(shù)的情況下,MLII系統(tǒng)的BIEA模型由(16)給出，即

π(Θ|D)∝πC(Θ)×L(D|Θ)=

(13)

式中,融合先驗(yàn)分布πC(Θ)考慮了直接先驗(yàn)分布πD(Θ)和間接先驗(yàn)分布πI(Θ)的2部分貢獻(xiàn);聯(lián)合似然函數(shù)由所有包含可用信息的父節(jié)點(diǎn)似然函數(shù)構(gòu)成,π(Θ|D)是考慮信息融合后模型參數(shù)集的后驗(yàn)分布。一般情況下,式(13)不存在解析解。通過(guò)使用各類MCMC方法,如Metropolis-Hastings方法(文獻(xiàn)[18])或Gibbs 抽樣方法(文獻(xiàn)[19])可獲得數(shù)值解。在本論文中,使用OpenBUGS(一種基于Gibbs抽樣做貝葉斯分析的開(kāi)源軟件包)進(jìn)行計(jì)算,詳細(xì)內(nèi)容可參考[20]。

3.2 可靠性指標(biāo)計(jì)算

通過(guò)式(16)可獲得所關(guān)心參數(shù)的聯(lián)合后驗(yàn)分布π(Θ|D),但在可靠性工程中,可能更關(guān)系某些特定的可靠性指標(biāo)，如可靠度、故障率等。

可靠度為

(14)

故障率為

(15)

平均故障時(shí)間(mean time to failure,MTTF)為

(16)

式中,D是MLII系統(tǒng)的可用數(shù)據(jù)集；R(tp|D)和λ(tp|D)分別是在tp時(shí)刻的可靠度與故障率。

現(xiàn)將BIEA方法的主要步驟羅列如下,其主要分析流程及與經(jīng)典貝葉斯方法的異同如圖5所示。

圖5 BIEA系統(tǒng)可靠性分析方法流程圖Fig.5 Flow chart of BIEA system reliability analysis method

步驟1依據(jù)各單元結(jié)構(gòu)特性,選取適當(dāng)?shù)拿枋瞿Ｐ?給出參數(shù)直接先驗(yàn)分布;

步驟2基于可用數(shù)據(jù)集建立似然函數(shù);

步驟3應(yīng)用貝葉斯更新計(jì)算參數(shù)后驗(yàn)分布;

步驟4基于系統(tǒng)結(jié)構(gòu)組成,計(jì)算系統(tǒng)可靠度函數(shù)表達(dá)式;

步驟5利用隨機(jī)變量轉(zhuǎn)換關(guān)系,獲得參數(shù)間接先驗(yàn)分布;

步驟6應(yīng)用BM算法計(jì)算參數(shù)融合先驗(yàn)分布;

步驟7基于融合先驗(yàn)分布計(jì)算更新后的參數(shù)后驗(yàn)分布;

步驟8模型輸出各類可靠性指標(biāo)。

4 方法驗(yàn)證與案例分析

4.1 模型描述

為了驗(yàn)證所提出的BIEA方法對(duì)MLII系統(tǒng)可靠性分析的有效性,選擇文獻(xiàn)[17]中的經(jīng)典雙層串聯(lián)系統(tǒng),其模型結(jié)構(gòu)與參數(shù)設(shè)置具有MLII系統(tǒng)的特性,如圖6所示。

圖6 雙層次結(jié)構(gòu)模型Fig.6 Two level structure model

具體地，C1、C2和C3分別為3類具有不同可靠性模型的組件。對(duì)于C1,其可靠性模型為logistic回歸模型,可用數(shù)據(jù)包含25個(gè)測(cè)試單元在11個(gè)時(shí)間點(diǎn)的二項(xiàng)型成敗記錄;對(duì)于C2,其可靠性模型服從雙參數(shù)的威布爾分布,可用數(shù)據(jù)包含25個(gè)測(cè)試單元的失效時(shí)間數(shù)據(jù);對(duì)于C3,選取的是性能退化模型,可用數(shù)據(jù)包含研究對(duì)象在10個(gè)時(shí)間點(diǎn)的性能退化數(shù)據(jù)。相關(guān)數(shù)據(jù)參如表2所示。為了保證對(duì)比的有效性,相關(guān)可靠性模型與模型參數(shù)的設(shè)置都與文獻(xiàn)[17]中例子一致。C1、C2和C3的可靠度函數(shù)分別由式(17)～式(19)給出。

R1(t|Θ1)=logit-1(θ1+η1t),Θ1=(θ1,η1)

(17)

(18)

(19)

基于表2中的有效數(shù)據(jù)集,可分別確定C1、C2和C3的似然函數(shù),而未知參數(shù)的后驗(yàn)分布可通過(guò)基于MCMC的抽樣方法計(jì)算得到。C1:成敗型數(shù)據(jù)(25樣本/時(shí)間點(diǎn)),C2:壽命數(shù)據(jù)(樣本量25),C3:退化數(shù)據(jù)(1樣本/單元/時(shí)間點(diǎn))。本文使用OpenBUGS軟件分別計(jì)算C1、C2和C3的可靠度隨時(shí)間變化曲線，結(jié)果如圖7所示,Mean為均值,Val 2.5%為置信度2.5%分位值,Val 97.5%為置信度97.5%分位值。

表2 C1、C2和C3的可用數(shù)據(jù)

圖7 C1、C2和C3的可靠度隨時(shí)間變化曲線Fig.7 Time variation curves of reliability for C1、C2 and C3

根據(jù)串聯(lián)系統(tǒng)結(jié)構(gòu)函數(shù),可得系統(tǒng)的可靠度函數(shù)為

R0(t|Θ0)=R1(t|Θ1)R2(t|Θ2)R3(t|Θ3)=

(20)

為了驗(yàn)證所提出的BIEA方法在信息分布不均衡情況下的有效性,分別對(duì)系統(tǒng)層參數(shù)的先驗(yàn)分布與可用數(shù)據(jù)考慮以下3類不同的情況:①先驗(yàn)分布為平坦的無(wú)信息先驗(yàn)分布且具有充足數(shù)據(jù)集(樣本量50);②相同的無(wú)信息先驗(yàn)分布,但僅為小樣本,稀疏數(shù)據(jù)集(樣本量10);③先驗(yàn)分布為非平坦的偏態(tài)分布,以此來(lái)模擬具有較強(qiáng)主觀意向的專家經(jīng)驗(yàn),同時(shí)樣本量保持為10的小樣本。為消除參數(shù)相關(guān)性的影響,假設(shè)系統(tǒng)的故障時(shí)間服從指數(shù)分布,故其可靠性模型僅包含唯一的參數(shù)λ。

對(duì)于每一種情況,都分別應(yīng)用簡(jiǎn)單貝葉斯方法(無(wú)BM),經(jīng)典BM方法和BIEA方法進(jìn)行參數(shù)估計(jì)與可靠性分析,并將結(jié)果做對(duì)比。需要注意的是,由于簡(jiǎn)單貝葉斯方法不區(qū)分直接先驗(yàn)分布和間接先驗(yàn)分布,故在進(jìn)行系統(tǒng)層的可靠性分析時(shí),是不包含底層信息的,而2種BM方法則包含所有信息。為了唯一確定系統(tǒng)可靠性預(yù)測(cè)曲線,假設(shè)系統(tǒng)的可靠度函數(shù)R0(t|λ)在某特定時(shí)間(t=20)服從Beta分布((beta(1,1) 或 beta(15,70)),那么根據(jù)λ和R之間的函數(shù)關(guān)系,可以得到關(guān)于參數(shù)λ的間接先驗(yàn)分布。這樣,融合先驗(yàn)分布就同時(shí)包含了直接先驗(yàn)分布和間接先驗(yàn)分布。系統(tǒng)層的數(shù)據(jù)是基于λ真值0.012抽取的仿真數(shù)據(jù)。3種情況中所涉及的先驗(yàn)分布與可用數(shù)據(jù)如表3所示。

表3 3種情況下系統(tǒng)層的可用數(shù)據(jù)

使用式(6)的形式構(gòu)建融合先驗(yàn)分布,由于直接先驗(yàn)分布πD(λ)和間接先驗(yàn)分布πI(λ)均服從beta分布,因此對(duì)于經(jīng)典BM方法和基于BIEA的BM方法,融合先驗(yàn)分布πC(λ)分別有

πC(λ)∝πD(λ)απI(λ)(1-α)

(21)

Beta((aD)α(aD)1-α,(bD)α(bD)1-α)

(22)

式中,a和b為beta分布參數(shù)。對(duì)于經(jīng)典BM方法,分別選擇不同的權(quán)重系數(shù)(α=0.2,0.4,0.6,0.8),混合后的融合先驗(yàn)分布如圖8所示。

圖8 不同權(quán)重因子(α=0.2,0.4,0.6,0.8)對(duì)應(yīng)的融合先驗(yàn)分布Fig.8 Combined prior of different weighing factor (α=0.2,0.4,0.6,0.8)

4.2 分析與討論

表4～表6分別給出了3種不同情況下,分別應(yīng)用簡(jiǎn)單貝葉斯方法、經(jīng)典BM,基于自更新權(quán)重系數(shù)的BIEA方法3種方法的參數(shù)估計(jì)結(jié)果;圖9給出了3種情況下,應(yīng)用3種不同方法所得的可靠性預(yù)測(cè)曲線對(duì)比結(jié)果,CB為簡(jiǎn)單貝葉斯,TM為經(jīng)典BM,AM為自更新BM。參數(shù)估計(jì)結(jié)果是利用OpenBUGS軟件中的統(tǒng)計(jì)工具箱計(jì)算完成后,在Matlab中繪圖得到的。

圖9 3種情況下不同方法的系統(tǒng)可靠性預(yù)測(cè)結(jié)果Fig.9 Prediction results of system reliability of different methodsunder three conditions

參數(shù)λ均值標(biāo)準(zhǔn)差誤差置信度2.5%分位值置信度50%分位值(中值)置信度97.5%分位值起始樣本點(diǎn)樣本量簡(jiǎn)單貝葉斯方法0.011 300.001 3751.417E-50.008 8160.011 240.014 181 00010 001經(jīng)典BM(α=0.2)0.011 870.001 5911.564E-50.008 9630.011 800.015 211 00010 001經(jīng)典BM(α=0.4)0.011 790.001 6031.613E-50.008 8800.011 730.015 121 00010 001經(jīng)典BM (α=0.6)0.011 760.001 6081.668E-50.008 8700.011 700.015 061 00010 001經(jīng)典BM(α=0.8)0.011 720.001 5951.646E-50.008 8190.011 660.014 991 00010 001BIEAλ0.011 730.001 5461.673E-50.008 8890.011 640.014 961 00010 001α0.322 400.258 7000.004 0720.011 3000.248 900.919 501 00010 001

表5 第2種情況下的參數(shù)統(tǒng)計(jì)結(jié)果

表6 第3種情況下的參數(shù)統(tǒng)計(jì)結(jié)果

由圖9可見(jiàn),在理想情況下(無(wú)偏先驗(yàn)分布&充足樣本)3種方法都可以得到令人滿意的結(jié)果(見(jiàn)圖9(a))。參數(shù)估計(jì)的結(jié)果對(duì)不同權(quán)重系數(shù)α的取值并不敏感(見(jiàn)表4)。此時(shí),3種方法得到的可靠性預(yù)測(cè)曲線幾乎是一致的,這一現(xiàn)象也與文獻(xiàn)[17]的描述相符。

當(dāng)樣本量逐漸減小后(從50減至10),應(yīng)用BM方法(經(jīng)典貝葉斯或BIEA)的分析與預(yù)測(cè)結(jié)果,其準(zhǔn)確性超過(guò)簡(jiǎn)單貝葉斯方法的結(jié)果(見(jiàn)圖9(b))。這說(shuō)明對(duì)于存在小樣本問(wèn)題的非理想情況,應(yīng)用BM可以改善分析精度,提高預(yù)測(cè)的準(zhǔn)確性。在第2種情況中,使用經(jīng)典BM方法時(shí),當(dāng)權(quán)重系數(shù)從0.2變化至0.8的過(guò)程中,參數(shù)估計(jì)的結(jié)果發(fā)生了變化(見(jiàn)表5),這是由于樣本量從50減至10后,似然函數(shù)的主導(dǎo)性顯著降低所致。

當(dāng)存在較強(qiáng)的專家經(jīng)驗(yàn)而可用數(shù)據(jù)又匱乏時(shí)，即信息非均衡的MLII系統(tǒng),使用經(jīng)典的BM方法也難以保證分析的準(zhǔn)確性(見(jiàn)圖9(c)的藍(lán)色曲線)。而對(duì)于基于BIEA的BM方法,即使先驗(yàn)分布中存在主觀偏差(beta(15,70)),其分析與預(yù)測(cè)結(jié)果仍然具有較高的精度(90%置信區(qū)間)。

事實(shí)上,在非理想情況下(偏差先驗(yàn)分布,小樣本有限數(shù)據(jù)),有限的數(shù)據(jù)集不足以構(gòu)建占據(jù)主導(dǎo)地位的似然函數(shù),故而先驗(yàn)分布對(duì)最終結(jié)果的影響較大。此時(shí),參數(shù)估計(jì)的結(jié)果對(duì)權(quán)重系數(shù)α的取值較為敏感,權(quán)重系數(shù)α須有較高精度。方法通過(guò)將權(quán)重系數(shù)α設(shè)為不確定性參數(shù),可以利用數(shù)據(jù)在推理過(guò)程中不斷更新α的取值,進(jìn)而自動(dòng)地平衡直接先驗(yàn)分布和間接先驗(yàn)分布對(duì)融合先驗(yàn)分布的影響。

經(jīng)典的BM方法將權(quán)重系數(shù)α設(shè)定為確定值,這對(duì)一般系統(tǒng)而言是可以接受的。但這對(duì)似然函數(shù)無(wú)法占據(jù)主導(dǎo)位置的MLII系統(tǒng)是不適用的。基于BIEA的BM方法通過(guò)將權(quán)重系數(shù)設(shè)為超參數(shù)的方法,在推理過(guò)程中使其與其他參數(shù)一起更新,通過(guò)這種方式,底層信息也參與到系統(tǒng)頂層對(duì)象的可靠性分析過(guò)程中。基于BIEA的BM方法在系統(tǒng)層面上重新求解了信息不均衡問(wèn)題,即使先驗(yàn)分布存在偏差,樣本量較小,參數(shù)估計(jì)結(jié)果的精度仍是可接受的。

5 結(jié) 論

本文就傳統(tǒng)可靠性分析與預(yù)測(cè)方法在處理含有MLII系統(tǒng)時(shí)的一些局限性,提出了基于BIEA的系統(tǒng)可靠性分析與預(yù)測(cè)方法。該方法主要針對(duì)在實(shí)際工程中存在的MLII系統(tǒng)的2類特點(diǎn)——多層次與信息分布不均衡，對(duì)經(jīng)典貝葉斯推理方法作出改進(jìn),使之可充分利用底層單元的完備數(shù)據(jù),自下而上地補(bǔ)償頂層匱乏的信息,獲得較為準(zhǔn)確的系統(tǒng)可靠性分析與預(yù)測(cè)結(jié)果。該方法可廣泛應(yīng)用于具有MLII系統(tǒng)特性的復(fù)雜機(jī)電設(shè)備的可靠性工程,降低分析結(jié)果的不確定性,提高預(yù)測(cè)精度。

參考文獻(xiàn)：

[1] LI M, LIU J, LI J, et al. Bayesian modeling of multi-state hierarchical systems with multi-level information aggregation[J]. Reliability Engineering & System Safety,2014,124(124):158-164.

[2] YANG L, ZHANG J, GUO Y, et al. Bayesian-based information extraction and aggregation approach for multilevel systems with multi-source data[J]. Journal of Systems Engineering and Electronics, 2017, 28(2): 385-400.

[3] RUI K, ZHANG Q, ZENG Z, et al. Measuring reliability under epistemic uncertainty: review on non-probabilistic reliability metrics[J]. Chinese Journal of Aeronautics,2016,29(3):571-579.

[4] PAN R, YONTAY P. Discussion of “Bayesian reliability: combining information”[J]. Quality Engineering,2017,29(1): 136-140.

[5] HAMADA M S, WILSON A G, REESE C S, et al. Bayesian Reliability[M]. New York: Springer, 2008.

[6] WILSON A G, FRONCZYK K M. Bayesian reliability: combining information[J]. Quality Engineering,2016,29(1):119-129.

[7] LI M, HU Q, LIU J. Proportional hazard modeling for hierarchical systems with multi-level information aggregation[J]. IIE Transactions, 2014, 46(2): 149-163.

[8] HAMADA M, MARTZ H F, REESE C S, et al. A fully Bayesian approach for combining multilevel failure information in fault tree quantification and optimal follow-on resource allocation[J]. Reliability Engineering & System Safety,2004,86(3):297-305.

[9] WILSON A G, REESE C S. Advances in data combination, analysis and collection for system reliability assessment[J]. Statistical Science, 2007, 21(4): 514-531.

[10] WOODS D C, OVERSTALL A M, ADAMOU M, et al. Bayesian design of experiments for generalized linear models and dimensional analysis with industrial and scientific application[J]. Quality Engineering, 2016, 29(1): 91-103.

[11] REESE C S, WILSON A G, GUO J, et al. A Bayesian model for integrating multiple sources of lifetime information in system-reliability assessments[J]. Journal of Quality Technology, 2011, 9(4): 489-495.

[12] YONTAY P, PAN R. A computational Bayesian approach to dependency assessment in system reliability[J]. Reliability Engineering & System Safety, 2016, 152: 104-114.

[13] POOLE D, RAFTERY A. Inference for deterministic simulation models: the Bayesian melding approach[J]. Publications of the American Statistical Association,2000,95(452):1244-1255.

[14] LIU Y, ZIDEK J V, TRITES A W, et al. Bayesian data fusion approaches to predicting spatial tracks: application to marine mammals[J].Annals of Applied Statistics,2016,10(3):1517-1546.

[15] LI Z S, GUO J, XIAO N C, et al. Multiple priors integration for reliability estimation using the Bayesian melding method[C]∥Proc.of the IEEE Reliability and Maintainability Symposium, 2017.

[16] FRENCH S. Group consensus probability distributions: a critical survey[J]. Bayesian Statistics, 1985, 2: 183-201.

[17] GUO J Q, WILSON A G. Bayesian methods for estimating system reliability using heterogeneous multilevel information[J]. Technometrics, 2013, 55(4): 461-472.

[18] GELFAND A E. Sampling-based approaches to calculating marginal densities[J]. Publications of the American Statistical Association, 1990, 85(410): 398-409.

[19] SMITH A F M, ROBERTS G O. Bayesian computation via the Gibbs sampler and related Markov chain Monte Carlo methods[J]. Journal of the Royal Statistical Society, 1993, 55(1): 3-23.

[20] KRUSCHKE J K. Doing Bayesian data analysis: a tutorial with R and BUGS[J]. Cognitivesciencesociety Org,2013,1(5):737-745.