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

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

0 引 言

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

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

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

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

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

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

1 系統(tǒng)描述與假設

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

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

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

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

威布爾分布通常被用于描述系統(tǒng)的故障過程,其概率密度函數(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時,即為指數分布。為了更好地描述預測模型,系統(tǒng)中信息源和單機的數量分別用n和m表示。在各個信息源和單機層面,用i(1≤i≤n)表示信息源的序號,用j(1≤j≤m)表示系統(tǒng)中單機的序號。在仿真中,用k表示根據單機壽命分布獲得的樣本量,l(1≤l≤k)表示單機壽命樣本xi, j(l)的序號。

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

2 融合建模與預測

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

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

圖2 系統(tǒng)級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之間生成連續(xù)均勻分布的k個隨機數α,使Fi, j(t)=α,得到k個單機壽命樣本。

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

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

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

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

(4)

(5)

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

(6)

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

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

算法 2抽取冷備系統(tǒng)的壽命樣本

步驟 1對于第j(j=1,2,…,m)個單機的第i(i=1,2,…,n)個信息源,從冷備系統(tǒng)中的單機壽命樣本{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)次,直到系統(tǒng)失效部件超過r個,得到系統(tǒng)壽命為Yi,G(l)=min CLi, j(l);

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

2.2 多源驗前分布推導

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

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

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

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

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

(14)

(15)

式中,

(16)

(17)

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

(18)

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

(19)

(20)

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

(21)

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

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

(22)

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

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

(23)

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

2.4 系統(tǒng)RUL預測

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

(24)

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

(25)

(26)

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

算法 3求解RUL預測結果

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

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

3 實例驗證

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

3.1 算例1

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

表1 單機的失效分布信息

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

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

表2 系統(tǒng)RUL預測結果

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

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

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

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

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

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

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

3.2 算例2

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

圖6 對比系統(tǒng)的可靠性框圖Fig.6 Reliability block diagram of comparison system

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

表3 對比系統(tǒng)的失效信息

4 結 論

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

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