以運(yùn)載火箭、風(fēng)電機(jī)組、新型裝備等為代表的產(chǎn)品由于可靠性試驗(yàn)周期短,實(shí)驗(yàn)成本高,實(shí)驗(yàn)數(shù)據(jù)短缺,在對其進(jìn)行可靠性評估等量化工作時面臨小子樣問題[1]。最大似然估計法等傳統(tǒng)可靠性評估方法在面對小子樣可靠性評估問題時容易產(chǎn)生較大偏差[2]。而貝葉斯推斷能融合并利用各種先驗(yàn)信息(包括歷史數(shù)據(jù)、相似產(chǎn)品數(shù)據(jù)以及專家經(jīng)驗(yàn))與試驗(yàn)信息,并用試驗(yàn)信息不斷更新先驗(yàn)信息,最終得到可靠性評估模型中待估參數(shù)和可靠性指標(biāo)的估計值與置信區(qū)間,能有效控制小樣本數(shù)據(jù)帶來的不確定性[3]。因此,貝葉斯推斷在小子樣可靠性評估中應(yīng)用廣泛。
在已有研究多集中于指數(shù)分布型產(chǎn)品或簡化威布爾型產(chǎn)品的情況下[4-6],對更一般雙參數(shù)的威布爾型產(chǎn)品建立基于貝葉斯推斷的可靠性評估模型。威布爾分布函數(shù)無共軛先驗(yàn)分布,假設(shè)待估參數(shù)服從Gamma分布[7],保證其正定性,利用先驗(yàn)信息與專家經(jīng)驗(yàn)確定Gamma分布中超參數(shù),得到待估參數(shù)先驗(yàn)概率密度函數(shù),結(jié)合似然函數(shù)生成后驗(yàn)概率密度函數(shù),應(yīng)用馬爾可夫鏈蒙特卡洛法具體為NUTS抽樣模擬待估參數(shù)的后驗(yàn)概率密度函數(shù),求得待估參數(shù)估計值與置信區(qū)間。以某型號風(fēng)力發(fā)電機(jī)齒輪箱的8次故障數(shù)據(jù)為算例,證明該方法的正確性與有效性,并完成對其可靠性評估。
貝葉斯推斷能將先驗(yàn)信息與觀察到的故障數(shù)據(jù)相結(jié)合,以推斷待估參數(shù)ξ(ξ1,ξ2,…,ξn)的數(shù)學(xué)特征。貝葉斯推斷的核心是貝葉斯定理,貝葉斯定理可以實(shí)現(xiàn)通過收集到的數(shù)據(jù)y(y1,y2,…,yn)對已有信息更新[8,9]。如下所示:
其中y(y1,y2,…,yn),表示觀測到隨機(jī)變量y的數(shù)據(jù)向量;p(ξ|y)為后驗(yàn)密度函數(shù);p(ξ)為先驗(yàn)密度函數(shù);f(y|ξ)為數(shù)據(jù)的抽樣密度函數(shù)即似然函數(shù);m(y)為數(shù)據(jù)的邊緣密度函數(shù)。根據(jù)貝葉斯定理,后驗(yàn)分布可成比例表示為似然函數(shù)與先驗(yàn)分布的乘積,即
其中m(y)與參數(shù)ξ無關(guān),被視為比例常數(shù)的一部分。獲得觀測數(shù)據(jù)前,對參數(shù)ξ(ξ1,ξ2,…,ξn)的信息集中于p(ξ),在獲得觀測數(shù)據(jù)y(y1,y2,…,yn)后生成似然函數(shù),并更新生成后驗(yàn)分布。
Weibull分布作為常用可靠性分析統(tǒng)計分布,對機(jī)械產(chǎn)品的故障時間數(shù)據(jù)擬合能力強(qiáng)[10,11],故障時間t的概率密度函數(shù)為
其中,β是形狀參數(shù);θ是尺度參數(shù),也稱為特征壽命。待估參數(shù)為β與θ。
可靠度函數(shù)為
累積分布函數(shù)為
平均故障間隔時間為
對于故障時間T(t1,t2,…,tN),T~Weibull(θ,β),似然函數(shù)為各故障時間概率密度函數(shù)乘積,即
當(dāng)Weibull分布的形狀參數(shù)和尺度參數(shù)均未知時,不存在適用的共軛先驗(yàn)分布。對于θ與β,假設(shè)其先驗(yàn)分布為Gamma分布,且兩個參數(shù)相互獨(dú)立,即
θ~Gamma(α1,λ1),β~Gamma(α2,λ2)
α1,α2>0為形狀參數(shù);λ1,λ2>0為尺寸參數(shù)。對于超參數(shù)α1,α2>0與λ1,λ2>0可根據(jù)產(chǎn)品可靠性先驗(yàn)信息,以矩估計與專家經(jīng)驗(yàn)相結(jié)合的方法求解。
此時θ概率密度函數(shù)(條件先驗(yàn))為f(θ|α1,λ1)=,均值對θ概率密度函數(shù)(條件先驗(yàn))化簡可得p(θ|α1,λ1)∝θα1-1exp(-λ1θ)。θ條件后驗(yàn)密度函數(shù)為:
對于以上條件后驗(yàn)密度函數(shù)由于無法用解析方法求得,需要應(yīng)用數(shù)值模擬方法—馬爾可夫鏈蒙特卡洛方法[12,13]。
馬爾可夫鏈蒙特卡洛法(MCMC)也稱為馬爾可夫鏈蒙特卡洛抽樣器,其將馬爾可夫鏈與蒙特卡洛法相結(jié)合。其 中 馬爾可夫鏈?zhǔn)侵敢幌盗须S 機(jī)樣本ξ(ξ1,ξ2,…,ξn),P(ξn+1|ξ1,ξ2,ξ3,…,ξn)=P(ξn+1|ξn)。即ξn+1取值只與當(dāng)前取值ξn有關(guān),而與之前取值無關(guān)。在馬爾可夫鏈中相鄰樣本(ξi,ξi+1)之間可以發(fā)生轉(zhuǎn)換,其轉(zhuǎn)換概率P(ξi→ξi+1)由轉(zhuǎn)換概率分布函數(shù)確定。在齊次馬爾可夫鏈中,生成的樣本分布將逐漸收斂,并趨于平穩(wěn)[14,15]。在貝葉斯推斷背景下,該平穩(wěn)分布即為待估參數(shù)后驗(yàn)分布[16,17]。但馬爾可夫鏈的初始樣本并不來源于平穩(wěn)分布,因此也不能代表該平穩(wěn)分布,該過程稱為老煉階段(burn-in),應(yīng)從馬爾可夫鏈中剔除。
設(shè)老練期生成樣本數(shù)Nburn-in,馬爾可夫鏈生成總樣本數(shù)N。當(dāng)算法進(jìn)行了Nburn-in次后,產(chǎn)生的仿真序列可以認(rèn)為是來自后驗(yàn)分布的樣本,而Nburn-in次之前的樣本屬于算法的老煉階段,不能認(rèn)為來自后驗(yàn)分布。因此,MCMC創(chuàng)造了一個馬爾可夫鏈,該馬爾可夫鏈的不變分布即為所求的后驗(yàn)概率分布。常用MCMC算法有Metropolis-Hastings算 法、Gibbs抽 樣 與Hamiltonian蒙特卡洛方法[18]。其中Hamiltonian蒙特卡羅(Hamiltonian Monte Carlo,HMC)方法也稱為混合蒙特卡羅(Hybrid Monte Carlo)方法,其采用系統(tǒng)動力學(xué)而不是概率分布來建議馬爾可夫鏈中的未來狀態(tài),使得馬爾可夫鏈能更有效地探索目標(biāo)分布從而實(shí)現(xiàn)更快的收斂。無掉頭抽樣器(No-U-Turn Sampler,NUTS)作為HMC的拓展,使用成對平均算法[19],可以在完全不進(jìn)行任何手動調(diào)整的情況下運(yùn)行,生成的樣本至少與精細(xì)的手動調(diào)整HMC一樣。
對威布爾型產(chǎn)品,參數(shù)向量中包含2個元素(θ,β),由1.3節(jié)知每個元素的條件后驗(yàn)分布為p(θ|T,β,α1,λ1),p(β|T,θ,α2,λ2)。
以某型號風(fēng)力發(fā)電機(jī)齒輪箱的8次故障數(shù)據(jù)為算例[20],驗(yàn)證上述可靠性評估模型。風(fēng)電機(jī)組齒輪箱故障狀態(tài) 數(shù) 據(jù)(4707,6167,9887,9991,13712,16543,18007,20793),均為停機(jī)故障,單位:時。使用最小二乘法計算得到,θ=19000,β=2.3。
對于先驗(yàn)分布,應(yīng)用矩估計法并借助專家經(jīng)驗(yàn),可確定以下方程式:
解得,λ1=3.6;α1=64800,λ2=1,α2=1。將結(jié)果代入θ條件后驗(yàn)密度函數(shù)與β條件后驗(yàn)密度函數(shù)。
采用基于Python平臺的PyMC3進(jìn)行NUTS抽樣,仿真環(huán)境為PyCharm Edu,抽樣次數(shù)設(shè)置為5000次,Nburn-in=1000,抽樣結(jié)果如表1所示。
表1 NUTS抽樣結(jié)果
即θ后驗(yàn)均值17998.02,標(biāo)準(zhǔn)差72.72;β后驗(yàn)均值2.11,標(biāo)準(zhǔn)差0.63。
生成后驗(yàn)密度直方圖和最高密度區(qū)域圖如圖1。
圖1 θ、β后驗(yàn)密度直方圖與最高密度區(qū)域
生成密度估計圖和痕跡圖如圖2、圖3。
圖2 θ密度估計圖和痕跡圖
圖3 β密度估計圖和痕跡圖
文獻(xiàn)[20]分別使用最小二乘法根據(jù)8次故障停機(jī)數(shù)據(jù)及文獻(xiàn)改進(jìn)方法根據(jù)8次停機(jī)故障數(shù)據(jù)加12次非停機(jī)故障數(shù)據(jù)對雙參數(shù)進(jìn)行估計,結(jié)果如表2所示。
表2 三種方法參數(shù)估計對比
通過對比可知,提出基于貝葉斯推斷的威布爾型產(chǎn)品可靠性評估方法在利用8次停機(jī)故障數(shù)據(jù)情況下,估計所得θ均值介于最小二乘法與文獻(xiàn)[20]改進(jìn)方法。在可利用故障數(shù)據(jù)一致情況下,本文方法得到的尺度參數(shù)θ估計值比最小二乘法更精確,相較文獻(xiàn)[20]方法更加合理。對于形狀參數(shù)β,估計值大于1,說明失效率隨時間增加,根據(jù)本文方法得出的故障規(guī)律與最小二乘法、文獻(xiàn)[20]改進(jìn)方法一致。
通過生成抽樣能量圖,可以看出能量(Marginal Energy Distribution)及能量轉(zhuǎn)換(Energy Transition Distribution)兩個分布范圍接近,寬窄相近,說明NUTS抽樣中失去的信息不重要,后驗(yàn)估計準(zhǔn)確。
圖4 抽樣能量圖
將NUTS抽樣產(chǎn)生的形狀參數(shù)與尺度參數(shù)帶入可靠性評估模型(4)(5)(6)(7),可得故障時間t的概率密度函數(shù)為:
可靠度函數(shù)為:
R(t)=exp[-(t/17998)2.11]
累積分布函數(shù)為:
F(t)=1-exp[-(t/17998)2.11]
平均故障間隔時間為:
由上述可靠性模型計算得出某型號風(fēng)力發(fā)電機(jī)齒輪箱的平均故障間隔時間為15943h,完成了可靠性評估。
本文針對更一般的雙參數(shù)威布爾型產(chǎn)品建立基于貝葉斯推斷的可靠性評估模型,在無共軛先驗(yàn)分布情況下,均假設(shè)其待估參數(shù)服Gamma從先驗(yàn)分布,推導(dǎo)出待估參數(shù)后驗(yàn)概率密度,而非簡化為指數(shù)分布型可靠性模型。充分利用已有數(shù)據(jù)與專家經(jīng)驗(yàn),應(yīng)用MCMC法中NUTS抽樣求解待估參數(shù)點(diǎn)估計值與置信區(qū)間,實(shí)現(xiàn)復(fù)雜密度函數(shù)的仿真模擬。經(jīng)過算例分析可知該模型的輸出結(jié)果正確可靠,既擴(kuò)大貝葉斯推斷的應(yīng)用范圍又能保證求解精度。模型通過調(diào)整先驗(yàn)分布可以求解不同威布爾型產(chǎn)品的待估參數(shù)與可靠度函數(shù),具有可重用性。對于不同威布爾型產(chǎn)品,借助先驗(yàn)信息與專家經(jīng)驗(yàn)確定待估參數(shù)先驗(yàn)分布后,通過修改PyMC3中參數(shù)信息,可用于求解待估參數(shù),有助于對小子樣產(chǎn)品的可靠性評估與可靠性管理。
圖5 故障時間概率密度函數(shù)
圖6 可靠度函數(shù)
在本文基礎(chǔ)上,可以進(jìn)一步研究融合多源數(shù)據(jù)的基于貝葉斯推斷可靠性評估,量化專家經(jīng)驗(yàn)控制主觀性,充分利用各種先驗(yàn)信息,進(jìn)行更準(zhǔn)確的小子樣產(chǎn)品可靠性評估與優(yōu)化。