基于貝葉斯融合與仿真的系統剩余壽命預測

宋兆理, 賈 祥, 郭 波, 程志君

(國防科技大學系統工程學院, 湖南 長沙 410073)

0 引 言

在工程實踐中,復雜系統的可靠性和剩余壽命(remaining useful life,RUL)預測是重要的研究問題。尤其是對于大型工程,例如衛星研制與發射等,可靠的系統級RUL預測可以為上層決策提供強有力的支持,從而科學化地降低整個系統工程的運營成本,提高系統的可用性。

單個設備的RUL稱為單機級RUL,而由多個單機或更多層級組成的復雜系統的RUL稱為系統級RUL。對于處于工作狀態的設備,無論是在單機層面還是在系統層面,都需要合適的RUL預測方法來監管設備的使用方式,制定將來的維護和更換計劃。與單機級RUL相比,系統級的RUL預測往往更加重要,也更加困難。單機級RUL預測的關鍵在壽命分布、數據特征、融合策略等方面[1-2],而系統級的RUL預測還需要考慮系統結構、信息折算等[3-5]。因此,探索有效的系統級RUL預測方法是可靠性研究領域的重要問題之一。近年來已經涌現出一些研究成果,大致分為3大類:① 基于貝葉斯理論的方法;② 基于隨機過程的方法;③ 統計學方法。

對于具有高可靠性和少量數據的設備,貝葉斯理論可以充分利用驗前信息和試驗數據。文獻[6]提出一種貝葉斯非參數估計方法,以獲取具有單機和系統驗前的系統驗后可靠性函數。基于貝葉斯框架,最小方差無偏估計器(minimum-variance unbiased estimator,MVUE)[7]或數據驅動的預測器[8]可以準確地進行系統狀態估計。在典型的系統可靠性結構中,串聯和并聯結構的工作原理相對簡單,因此基于貝葉斯理論的相關研究成果和公式推論較多[9-12];而對于表決系統和冷備系統,由于其結構的復雜性,簡便易操作的方法較少。

對于設備使用過程中產生的性能監測數據,很多方法通過建立基于隨機過程的可靠性模型進行系統級RUL預測。一定程度上來說,這些方法能夠還原系統在工作狀態下的真實退化過程,從而保證預測精度。目前,利用隨機過程預測系統級RUL的模型很多,例如Markov過程[13]、Gamma過程[14]、Wiener過程[15]、Petri網[16]。文獻[17]將系統視為一個整體考慮其性能退化過程,提出了一種基于解析模型的RUL預測方法。這些方法運用場景廣泛,能夠在無法獲得失效數據的情況下預測系統的RUL,但數據來源相對單一,在預測對象具備多源信息的條件下,難以基于不同類型的信息進行融合預測。

統計學方法內涵豐富,包括非參數分析、回歸分析、蒙特卡羅(Monte-Carlo,MC)仿真等一系列具體方法。文獻[18-19]使用極大似然估計方法進行了衛星和衛星子系統的非參數可靠性分析和威布爾擬合。統計學方法對不確定性問題具有良好的魯棒性,能夠比較準確地預測系統RUL,所以被越來越多地用于隨機疲勞分析和可靠性分析中[20-22]。然而,由于隨機仿真預測的精確度是建立在足夠的時間成本之上,所以難以滿足工程中存在的實時預測的需求。

通常,系統級RUL預測往往趨于在單機級融合所有信息,然后通過系統結構將其轉換到系統級,導致數據損耗。如果將融合過程推遲到系統級[23],則原始的單機級信息可以直接折合到系統,減少信息折損。因此,本文提出一種將貝葉斯融合和隨機仿真相結合的方法,在系統級融合現場信息和單機提供的多源驗前信息,進而預測系統的RUL。該方法既能借助貝葉斯理論將多源信息融合運用,又能發揮仿真方法便捷準確的優勢,為系統級RUL預測提供一套完整可行的流程。

1 系統描述與假設

系統通常由具有特定結構的若干個單機組成,例如圖1所示的衛星平臺上的某功能系統。該系統可以分為3個子系統,分別由S1、S2和S3表示。具體地,S1子系統是包含5個C1單機的5中取3冷備結構;S2子系統是由3個C2單機,2個C3單機,1個C4單機和6個C6單機組成的混聯結構;S3子系統是包含4個C5單機的2/4(G)表決結構。

圖1 示例系統的可靠性框圖Fig.1 Reliability block diagram of the prototype system

根據可靠性框圖,可以通過逐層劃分,將復雜的系統分解為具有特定關聯形式的單機,從而呈現出系統的邏輯結構。

本問題聚焦于系統級的RUL預測,所以不考慮單機級的計算過程,而是基于所有單機的多源信息已知的前提構建模型。當任何單機的信息源存在少量缺失時,將利用易于獲得的專家信息(如可靠度估計值、壽命估計值等)來填補這些丟失的數據。這樣可以保證每個信息源都能獲取到系統的無信息驗前分布,并與系統級的試驗信息相結合來預測RUL,從而避免該方法不適用于不完整信息的情況,增強其適用性。

威布爾分布通常被用于描述系統的故障過程,其概率密度函數(probability density function,PDF)、可靠性函數和累積概率分布函數(cumulative probability function,CDF)分別為

f(t)=λβtβ-1exp(-λtβ)

(1)

R(t)=exp(-λtβ)

(2)

F(t)=1-exp(-λtβ)

(3)

式中,λ和β為威布爾分布的參數。值得注意的是,當β=1時,即為指數分布。為了更好地描述預測模型,系統中信息源和單機的數量分別用n和m表示。在各個信息源和單機層面,用i(1≤i≤n)表示信息源的序號,用j(1≤j≤m)表示系統中單機的序號。在仿真中,用k表示根據單機壽命分布獲得的樣本量,l(1≤l≤k)表示單機壽命樣本xi, j(l)的序號。

考慮到系統RUL預測的復雜性,本文建立的模型適當簡化了實際問題,提出以下3個假設:① 系統的結構清晰,可以用可靠性框圖來描述;② 系統中所有單機都是獨立的;③ 除工作單機外,連接部件均完全可靠,如冷備系統中的開關等。

2 融合建模與預測

根據多源信息,可以通過單機級的壽命預測方法[1-2,24]獲得與每個單機的不同信息源相關的壽命分布。本模型不考慮單機級的計算過程,將單機壽命分布作為預測模型的原始信息。此外,系統級還提供了少量的現場試驗數據,其中包含被監測系統的壽命數據和工作狀態。對于具有極少失效特征的高可靠性系統,極高的測試成本導致測試數量很少,所以現場試驗信息是系統RUL預測的重要參考。

本文提出的RUL預測流程如圖2所示,每種信息都可以獨立計算,并基于貝葉斯理論在系統級融合來自每個信息源的驗后分布,提高最終預測結果的準確性。應當注意的是,當可以確定單機的多源信息來自同一設備時,在貝葉斯估計中不需要進行一致性檢驗。

圖2 系統級RUL預測流程Fig.2 Prediction procedure of system RUL

2.1 壽命樣本隨機仿真

對于每個信息源,已知單機壽命的CDF和PDF分別為Fi, j(t)和fi, j(t),其中i=1,2,…,n,j=1,2,…,m。借助反函數法和從0到1的隨機數,可以通過仿真獲得單機的壽命樣本xi, j(l)(l=1,2,…,k)。具體的采樣步驟如算法1所示。

算法 1抽取單機的壽命樣本

步驟 1對于第j(j=1,2,…,m)個單機的第i(i=1,2,…,n)個信息源,在0～1之間生成連續均勻分布的k個隨機數α,使Fi, j(t)=α,得到k個單機壽命樣本。

步驟 2對于系統中的所有m個單機,重復步驟1,得到第i個信息源對應的所有單機的壽命樣本。

步驟 3對于已知的所有n個信息源,重復步驟1和步驟2,得到每種信息源對應的所有單機的壽命樣本{xi, j(l)}(i=1,2,…,n;j=1,2,…,m;l=1,2,…,k)。

接下來,需要借助這些壽命樣本將單機級的信息傳遞到系統級。以4種常見的可靠性結構為例進行說明,用Yi,S、Yi,P、Yi,G和Yi,C分別表示串聯結構、并聯結構、表決結構和冷備結構的壽命樣本集。實際上,這些傳遞原則不僅限于單機級,同樣也適用于子系統級和系統級。

從仿真的角度來看,根據各個單元的壽命樣本集合,串聯結構、并聯結構和r/m: G表決結構的系統壽命樣本分別為

(4)

(5)

Yi,G(l)=xi,m-r+1(l)

(6)

式(6)中的單機壽命樣本xi,m-r+1(l)來自將表決系統中所有樣本升序排列后得到的樣本集合{xi,1(l)≤xi,2(l)≤…≤xi,m(l)}。

在m中取r冷備結構,假設所有備件都相同,則每個位置的故障過程都是獨立的[25]。對于任一部件位置,每當工作單機發生故障時,都會有一個冷備單機被激活,并接替其繼續工作,從而延長了該部件位置的壽命。以此類推,當所有冷備單機都被激活后,一旦有單機發生故障,系統將停止工作。因此,使用壽命最短的部件位置決定了冷備系統的壽命。具體抽樣過程如算法2所示。

算法 2抽取冷備系統的壽命樣本

步驟 1對于第j(j=1,2,…,m)個單機的第i(i=1,2,…,n)個信息源,從冷備系統中的單機壽命樣本{xi, j(l)}(l=1,2,…,k)中找到m個樣本作為最開始工作的單機,記為CLi, j(l)=xi, j(l),已失效的單機數p=0。

步驟 2找到min CLi, j(l),令CLi, j(l)=CLi, j(l)+xi, j+1+p(l),更新p=p+1。

步驟 3重復步驟2(m-r)次,直到系統失效部件超過r個,得到系統壽命為Yi,G(l)=min CLi, j(l);

步驟 4重復步驟1～步驟3k次,得到冷備系統的壽命樣本集Yi,C。

2.2 多源驗前分布推導

根據單機級信息獲得系統級不同信息源對應的壽命樣本Si(1≤i≤n)后,采用貝葉斯方法將這些信息視為驗前信息,與系統級信息融合進行RUL預測。在貝葉斯融合中,主要步驟包括將多源信息分類為驗前信息和現場信息,由驗前信息確定驗前分布以及融合現場試驗數據得出驗后分布。

各個信息源的系統壽命樣本集合可以近似描述壽命的分布特征,被視為系統的驗前信息。同時,系統的已記錄壽命數據被視為現場試驗信息。以下是從一個信息源的角度進行分析計算的方法,不同信息源的過程是相同的。

通過對系統壽命樣本進行排序,假設顯著性水平為α,易得壽命的100(1-α)%置信區間為[QL,QH]。置信區間反映了壽命樣本的分布,可以用來計算系統的驗前分布,從而達到將單機級信息折合到系統級的效果。

為了根據不完全信息獲得參數的驗前分布,通常采用最大熵方法(maximum entropy method,MEM),其中信息熵可以作為目標函數[26]。顯然,RUL的置信區間[μL,μH]滿足Fτ(μL)=α/2和Fτ(μH)=1-α/2,其中Fτ(t)表示在τ時刻系統RUL的CDF。在系統壽命服從威布爾分布的情況下,置信區間的下限μL和上限μH分別為

(7)

根據MEM原理[27],約束問題可以表示為

(8)

(9)

由工程經驗假設λ服從Gamma分布,β服從均勻分布,λ和β互相獨立,則其融合驗前分布為

(10)

式中,a和b是需要計算的超參數,λ∈[0,∞),β∈[βL,βH]。

(11)

根據式(7)～式(11),可以通過

(12)

計算得到a和b。式中,M1和M2是懲罰因子,需要通過數值方法計算[27],且

(13)

2.3 融合驗后分布推導

隨著多源信息驗前分布的獲得,必須將單機傳遞的驗前信息與系統的現場數據結合運用以減少不確定性。通常,多源信息的融合驗前分布表示為

(14)

(15)

式中,

(16)

(17)

根據MLE-II原理,融合權重πi隨L(D|πi)增加而增加,表示為

(18)

用P(D|λ,β)表示系統現場數據的似然函數,則系統的驗后分布為

(19)

(20)

則融合驗后分布可以表示為

(21)

式中,πi(λ,β|D)(i=1,2,…,n)是每個信息源的驗后分布。不難發現,融合驗后分布即為由多源驗前分布獲得的驗后分布的加權總和。

根據貝葉斯理論,第i個驗后分布可以表示為

(22)

對于一個壽命服從威布爾分布的系統,類似于式(15),現場試驗數據的似然函數為

P(D|λ,β)=λgβgMβ-1e-λN

(23)

因此,第i個信息源的融合權重wi及其驗后分布可以分別由式(20)、式(22)和式(23)得到,從而通過式(21)加權融合獲得最終的驗后分布。

2.4 系統RUL預測

根據式(1)～式(3),在獲得參數驗后分布的條件下,系統RUL的點估計[5]為

(24)

同時,由置信區間的滿足條件可以獲得系統RUL的置信區間下限μL和上限μH分別為

(25)

(26)

得到融合驗后分布之后,由于其形式復雜無法直接求解,故采用MC仿真[22]近似計算,過程如算法3所示。

算法 3求解RUL預測結果

步驟 1根據融合驗后分布π(λ(ε),β(ε)|D)抽樣,得到?組參數樣本{λ(ε),β(ε)}(ε=1,2,…,?)。

步驟 3對樣本求均值作為系統RUL的無偏估計,得到預測結果。

3 實例驗證

為進行方法驗證,列舉兩個算例具體分析。算例1通過與仿真方法對比,驗證了本文方法的正確性,算例2通過與現有方法對比,驗證了本文方法的可行性。

3.1 算例1

本算例的系統結構如圖1所示。根據工程實際問題提供的數據,在給定單機的相似產品壽命信息、歷史壽命信息和性能監測信息3個信息源(n=3)的條件下,系統中6種單機的壽命服從不同的指數分布和威布爾分布。為了簡化計算,本例認為來自相同單機的不同信息源提供相同的壽命分布。根據模型假設,單機的壽命分布相互獨立。具體地,單機C1、C2、C3、C4、C5服從威布爾分布,單機C6服從指數分布,參數如表1所示。根據圖2中的系統級RUL預測流程,給定來自多個信息源的單機壽命分布信息,每個分布生成k個樣本。為分析k的不同取值對結果的影響,將k分別設置為1 000、2 000、5 000、10 000、20 000和50 000。然后,根據圖1中的系統可靠性框圖,可以將單機壽命信息折合為系統的驗前信息,并獲得系統樣本集。

表1 單機的失效分布信息

根據工程經驗,可以認為威布爾型系統形狀參數的取值范圍是2～10。單機級數據被視為系統的驗前信息,該信息將通過參數模型結合系統級現場數據確定驗后分布。給定的系統現場試驗數據為不同時間開始工作的設備的定時截尾數據,分別是7 503 h、4 854 h、5 642 h、2 638 h、4 097 h,將其視為無失效數據。最終,當系統工作到τ時刻時,將獲得參數的驗后PDF,從而估計系統RUL。在本實例中,對于同一單機,由于不同的信息源來自相同的壽命分布,因此預測結果中每個信息源所占的權重幾乎相等。

為了證明本文方法的有效性和可行性,將用仿真方法與其進行比較。在仿真中,不考慮多源信息融合,而是由表1中給出的信息直接提取樣本,并將其按照第2.1節中的折算原則折合到系統級,從而作為驗前信息計算RUL。仿真結果可以被認為是實際值,用作驗證本文方法的參考。雖然本方法基于大量的仿真抽樣,但算法的耗時都較小,以秒為單位,所以此例沒有對算法的運行時間展開分析。表2匯總了兩種方法在τ=2年時的RUL預測結果,包括點估計和80%置信區間,以及對應點估計的均方誤差(mean square error,MSE)。圖3展示了不同k取值情況下通過兩種方法獲得的PDF曲線。

表2 系統RUL預測結果

以仿真方法所得的系統RUL點估計預測結果為對比標準,本文方法所得點估計預測結果的相對誤差如圖4所示。

由表2可以看出,本文方法獲得的系統RUL點估計和置信區間結果與仿真方法所得結果近似吻合,體現了本文方法的可行性。對于樣本量k的不同取值,本文方法的MSE均小于仿真方法的MSE,說明本文提出的預測模型所得結果分布更集中,具有更好的收斂性。由圖3可以看出,隨著樣本量k的增大,本文方法所得PDF曲線與仿真方法所得PDF曲線的重合度增加,反映了仿真實驗“樣本量越大,結果越精確”的特性。然而,當k不斷增大時,仿真抽樣的時間成本也不斷提高,雖然本例的時間成本可忽略不計,但對于實際中的大型復雜系統來說,樣本量的增大可能會使算法的運行時間成倍增長。因此,需要在保證結果精確性和較低時間成本的前提下,設置合適的樣本量大小。通過分析圖3和圖4可以判斷,當k=10 000時,預測結果已經基本穩定,且點估計預測的相對誤差較低(小于5%)。

圖3 不同k取值下本文方法與仿真方法所得的PDFFig.3 PDFs from the proposed method and simulation method with different k

圖4 系統RUL點估計的相對預測誤差Fig.4 Relative prediction error of system RUL point estimation

為了探究本文提出的RUL預測方法對于不同評估時期的適用性,在實例給定的τ=2年的前提下,將預測時間推遲,得到RUL點估計和置信區間的預測結果(設k=10 000)如圖5所示。

圖5 不同τ取值下系統RUL預測結果Fig.5 Prediction results of system RUL with different τ

由圖5可以看出,系統RUL預測結果隨著τ增大而減小,并且τ越大,τ與RUL預測值之和(即期望的系統壽命)越小。對于正常工作的設備,壽命期望值應維持不變,而此例中壽命隨τ增大而減小,表明本方法對于工作到中后期的設備RUL預測趨于保守。這一現象與系統壽命服從威布爾分布的特性有關,其故障率會隨著工作時間的延長而增大,導致系統加速退化,屬于“浴盆曲線”的耗損故障階段[29]。

3.2 算例2

本算例將本文中的預測方法應用于文獻[30]中的實際問題(已省略時間單位),以驗證本方法的可行性。該對比系統是一個簡化的衛星系統,由7個單機組成,其結構如圖6所示。

圖6 對比系統的可靠性框圖Fig.6 Reliability block diagram of comparison system

通過比較發現,利用文獻[30]提供的系統和失效信息(見表3),及其壽命數據5.4, 8.59, 6.34, 6.08, 9.26, 5.11, 10.55, 6.37, 10.56, 4.73，可以應用本方法預測該衛星系統的RUL,并且預測結果與原方法所得的結果近似。另外,本方法還能預測RUL的置信區間,從而提供更豐富的評估結果。

表3 對比系統的失效信息

4 結 論

本文的問題來自工程實踐中具有高可靠性和極少試驗數據的衛星平臺,在單機級RUL研究相對成熟的基礎上,結合貝葉斯理論和隨機仿真,提出融合多源信息的系統級RUL預測方法。當然,除了衛星平臺系統,本文方法也適用于具有類似特征的其他系統。本文提出的方法具有較強的適應性,介紹的實例僅為驗證提出方法的一種可能情況。在工程實踐中,如果單機的壽命分布未知,可以通過性能退化監測來獲得本文方法所需的壽命樣本,同樣能利用本方法進行系統RUL預測。本方法的優勢還在于利用仿真將單機的失效信息傳遞到系統級,減少分析推導和計算的復雜操作;在系統級融合多源信息,減小信息折損,保證已知數據的充分運用;幫助決策者了解不同信息源在RUL預測中所占的權重,為后期的大型設備運行策略和任務規劃提供支撐。

基于現有研究,后續將對更復雜的系統展開RUL分析,不局限于由獨立單機組成的系統,深入考慮單機的相關性失效,以更貼近工程實際的應用需求。