張英芝 翟粉莉 鄭玉彬 朱繼微 侯勝冬

(數(shù)控裝備可靠性教育部重點實驗室∥吉林大學(xué) 機械與航空航天工程學(xué)院，吉林 長春 130022)

自威布爾模型在產(chǎn)品的壽命研究中應(yīng)用，關(guān)于兩參數(shù)威布爾模型的參數(shù)估計問題就得到廣泛關(guān)注。常用的威布爾模型參數(shù)估計方法有最小二乘估計、矩估計、極大似然估計、貝葉斯估計等。其中最小二乘法因?qū)颖玖繘]有要求且原理及計算簡單，故在工程實際中得到較廣泛的應(yīng)用[1]。

應(yīng)用中發(fā)現(xiàn)，經(jīng)典最小二乘法估計參數(shù)雖滿足無偏性，但因異方差的存在，使得回歸模型不再顯著[2]。故相關(guān)學(xué)者提出加權(quán)最小二乘估計，既滿足無偏性又為最佳線性無偏估計，但余曉燕[3]通過數(shù)學(xué)證明該方法雖能在權(quán)重為各樣本點擾動項方差的倒數(shù)時誤差方差達到最小，但仍不滿足有效性。

王永山[4]分析不同先驗信息融合方法對貝葉斯參數(shù)估計精度的影響。劉正賢等[5]將威布爾分布轉(zhuǎn)化為極值分布，再利用貝葉斯公式進行計算，雖簡化了計算過程，但工作量仍然很大。Sawadogo等[6]針對運用牛頓迭代求解似然函數(shù)時的缺點，提出用遺傳算法求解似然函數(shù)的方法，但因不同樣本下似然函數(shù)表達式不同導(dǎo)致計算復(fù)雜度增加。Ehab等[7]研究了最大似然估計和貝葉斯估計在威布爾模型參數(shù)估計中的應(yīng)用，發(fā)現(xiàn)不同樣本下的參數(shù)估計模型推導(dǎo)難度較大，且不具有通用性。董力等[8]將遺傳算法與最小二乘法結(jié)合進行威布爾模型參數(shù)估計，此方法雖減少數(shù)學(xué)計算過程，但不能彌補最小二乘法的不足。許偉等[9]提出基于混沌模擬退火的PSO算法求威布爾模型參數(shù)，該方法避免了精確解直接求解時的大量計算，但是其目標(biāo)沒有考慮經(jīng)驗分布函數(shù)，降低了估計精度。

綜上可知，當(dāng)前威布爾模型參數(shù)估計研究多考慮計算工作量或求解精確性，忽略模型總偏差。因此，借鑒智能算法與數(shù)學(xué)方程結(jié)合進行求解的研究方法，文中提出一種基于累積誤差平方和最小的兩參數(shù)威布爾模型參數(shù)估計方法，應(yīng)用粒子群算法實現(xiàn)參數(shù)快速求解，采用柯爾漠哥洛夫-斯米爾洛夫檢驗(K-S檢驗)進行模型檢驗，用均方根誤差(RMSE)來衡量參數(shù)估算精度。結(jié)合相關(guān)文獻數(shù)據(jù)進行方法應(yīng)用有效性證明。

1 威布爾模型

兩參數(shù)威布爾模型的故障分布函數(shù)為

(1)

式中，β>0，α>0，β為形狀參數(shù)，α為尺度參數(shù)。

對式(1)兩邊連取兩次自然對數(shù)得：

ln[-ln[1-F(t)]]=β[ln(t)-ln(α)]

(2)

令

y=ln[-ln(1-F(t))]，x=ln(t)

B=β，A=-βln(α)

(3)

可得回歸方程：

y=Bx+A

(4)

2 累積誤差平方和最小參數(shù)估計

2.1 目標(biāo)函數(shù)及參數(shù)計算流程

文中以累積誤差平方和最小作為目標(biāo)構(gòu)建參數(shù)估計模型，并結(jié)合完全樣本、定時截尾樣本與定數(shù)截尾樣本情況進行具體參數(shù)估計。

目標(biāo)函數(shù)為：

(5)

具體計算步驟如下：

步驟2應(yīng)用K-S法進行擬合性檢驗。

步驟3根據(jù)式(6)計算均方根誤差，判斷所提方法的優(yōu)劣。判斷標(biāo)準(zhǔn)為均方根誤差值越大則表明參數(shù)估計精度越差[10]。

均方根誤差計算公式為

(6)

文中對于同一類型可靠性數(shù)據(jù)，分別采用經(jīng)典最小二乘法、加權(quán)最小二乘法以及基于累積誤差平方和最小法3種方法進行威布爾模型參數(shù)估計，并計算均方根誤差，以此為指標(biāo)驗證文中所提方法的有效性。

2.2 粒子群算法求解流程

粒子群算法(PSO)從鳥類捕食行為特征中得到啟發(fā)并用于求解優(yōu)化問題，算法中每個粒子代表問題的一個潛在解，每個粒子對應(yīng)一個由適應(yīng)度函數(shù)決定的適應(yīng)度值。在算法中每個粒子的速度都隨自身及其它粒子移動經(jīng)驗進行動態(tài)調(diào)整，從而實現(xiàn)個體在可解空間中的最優(yōu)。

算法中每個粒子代表一個潛在解，粒子數(shù)目一般取20-60，文中取40個粒子。根據(jù)參考文獻[11- 12]，設(shè)定粒子最大速度為0.8，設(shè)定粒子初始速度為0.8rand(1，2)，rand(1，2)表示每個粒子隨機產(chǎn)生一組2維速度初值。依據(jù)式(7)、(8)更新其速度與位置。在算法中設(shè)置算法的停止條件為迭代次數(shù)大于某一個數(shù)(MaxNum)。

(7)

(8)

式中，V表示速度，X表示位置，i代表第i個粒子，i=1，2，…，40，c1、c2為學(xué)習(xí)因子，分別代表粒子自身學(xué)習(xí)能力與全局學(xué)習(xí)能力，一般為非負(fù)常數(shù)，常取0-2，文中根據(jù)文獻[12]，取c1=c2=1.494。k為當(dāng)前迭代次數(shù)，1≤k≤MaxNum；r1、r2為[0，1]內(nèi)的隨機數(shù)，d代表每個粒子的第d個估計參數(shù)，估計參數(shù)個數(shù)為2，故d=1，2，分別代表尺度參數(shù)與形狀參數(shù)，Pi、Pg分別為個體極值與群體極值。

粒子群算法的步驟如下：

步驟1初始化粒子群的數(shù)目，設(shè)定學(xué)習(xí)因子值、粒子位置與速度。

步驟2根據(jù)式(5)計算個體適應(yīng)度值，將個體最優(yōu)值(p_ibest)設(shè)為群體最優(yōu)值(p_gbest)。

步驟3依據(jù)式(7)、(8)更新粒子速度與位置，并計算新位置上個體的適應(yīng)度值。

步驟4取粒子當(dāng)前適應(yīng)度與個體極值p_ibest適應(yīng)度的最小值，將其賦值給p_ibest，并記錄當(dāng)前參數(shù)值。

步驟5對比個體極值p_ibest與群體最優(yōu)值p_gbest值，若個體適應(yīng)度值更小，則將其賦值給p_gbest，反之保持p_gbest不變。

步驟6若迭代次數(shù)小于設(shè)定值，則返回步驟3，否則，算法結(jié)束，同時輸出最優(yōu)參數(shù)值。

2.3 經(jīng)驗分布函數(shù)

(1)完全樣本下經(jīng)驗分布函數(shù)

將n個完全樣本的故障數(shù)據(jù)按從小到大排序，則t1

(9)

(2)不完全樣本下經(jīng)驗分布函數(shù)

1)定時截尾樣本經(jīng)驗分布函數(shù)

設(shè)N為參加試驗樣品個數(shù)，定時截尾試驗下產(chǎn)生n個故障數(shù)據(jù)和N-n個截尾數(shù)據(jù)。

為修正截尾數(shù)據(jù)對故障數(shù)據(jù)秩次影響，采用Johnson法，引入秩增量修正秩次變動。

設(shè)At(i-1)為第i-1個故障數(shù)據(jù)t(i-1)秩次，定義At0=0；j為所有N個數(shù)據(jù)從小至大的排列秩次，即j=1，2，…，N。

第i個故障數(shù)據(jù)ti的秩次Ati為

秩次修正后的經(jīng)驗分布函數(shù)計算公式為

(11)

2)定數(shù)截尾試驗下經(jīng)驗分布函數(shù)

設(shè)n為參加試驗樣品個數(shù)，定數(shù)截尾試驗下產(chǎn)生r個故障數(shù)據(jù)。文中采用“殘存比率法”來估計定數(shù)截尾樣本下產(chǎn)品在一定時間區(qū)間內(nèi)殘存概率，各故障時刻產(chǎn)品可靠度值，進而計算經(jīng)驗分布函數(shù)。

S(ti)為產(chǎn)品在時間區(qū)間(ti-1，ti)內(nèi)殘存概率，計算式為

(12)

因此，產(chǎn)品在ti時刻的可靠度R(ti)可表示為

R(ti)=R(ti-1)S(ti)

(13)

式中，R(ti-1)為產(chǎn)品在ti-1時刻的可靠度，且R(t0)=1。

(14)

3 演算驗證

3.1 完全樣本

文中采用文獻[13]中滾動軸承208在載荷為 606 N，轉(zhuǎn)速為4 000 r/min的條件下，21個樣本壽命時間與經(jīng)驗分布值如表1所示。

表1 軸承壽命時間與經(jīng)驗分布值

采用經(jīng)典最小二乘法、加權(quán)最小二乘法及基于累積誤差平方和最小法得到參數(shù)估計值，如表2所示.3種參數(shù)估計方法的K-S檢驗值(Dmax)和均方根誤差如表3所示。

表2 完全樣本威布爾模型參數(shù)估計值

表3 均方根誤差和K-S檢驗值

顯著性水平α=0.1，查表得K-S檢驗臨界值為0.258 6，即3種估計方法均通過檢驗，表明滾動軸承壽命服從兩參數(shù)威布爾分布。

由表3可知，基于累積誤差平方和最小的參數(shù)估計法明顯優(yōu)于其他兩種方法，且較最小二乘法、加權(quán)最小二乘法估計精度分別提高24.94%、8.23%。

將3種參數(shù)估計法得到的可靠度函數(shù)與近似中位秩公式所求經(jīng)驗可靠度函數(shù)繪圖，如圖1所示。

圖1 完全樣本可靠度曲線

3.2 定時截尾樣本數(shù)據(jù)

文中參考文獻[13]，16個機電元件試驗持續(xù) 1 510 個工作周期，其中有10個元件失效。根據(jù)平均秩次法計算得各元件的經(jīng)驗分布值，如表4所示。

表4 元件失效時間與經(jīng)驗分布值

同理計算得到參數(shù)估計值如表5所示，RMSE與Dmax如表6所示。查表得3種估計方法均能通過假設(shè)檢驗，表明機電元件的壽命服從兩參數(shù)威布爾分布。

表5 定時截尾樣本威布爾模型參數(shù)估計值

由表6可知，基于累計誤差平方和最小法的參數(shù)估計精度最高，較其他兩種算法精度分別提高97.44%和37.18%。

表6 均方根誤差與K-S檢驗值

將3種方法的可靠度函數(shù)與經(jīng)驗可靠度函數(shù)繪圖，如圖2所示。由圖2可知，3種估計方法得到的可靠度隨著時間的推移偏差逐漸增大，但應(yīng)用累積誤差平方和最小的參數(shù)估計法得到的可靠度與經(jīng)驗可靠度是最接近的。

圖2 定時截尾樣本可靠度曲線

表7 樣品失效時間與經(jīng)驗分布值

3.3 定數(shù)截尾樣本數(shù)據(jù)

文中參考文獻[13]，13件樣品做壽命試驗，當(dāng)失效數(shù)r=10時停止試驗。失效時間以及根據(jù)殘存比率法計算得到的經(jīng)驗分布值如表7所示。同理計算得到參數(shù)估計值如表8所示，RMSE與Dmax如表9所示。查表得3種估算方法均能通過假設(shè)檢驗，表明該樣本的壽命服從兩參數(shù)威布爾分布。

表8 定數(shù)截尾樣本威布爾模型參數(shù)估計值

從表9可知，基于累積誤差平方和最小參數(shù)估計法的精度最高，較經(jīng)典最小二乘法和加權(quán)最小二乘法精度分別提高1.72%和1.15%。

表9 均方根誤差與K-S檢驗值

根據(jù)定數(shù)截尾樣本下3種參數(shù)估計法得到的可靠度函數(shù)與近似中位秩公式所求經(jīng)驗可靠度函數(shù)繪圖，如圖3所示。

圖3 定數(shù)截尾樣本可靠度曲線

圖3表明，定數(shù)截尾樣本下3種參數(shù)估計方法得到的可靠度函數(shù)隨著時間的推移差別逐漸增大，基于累積誤差平方和最小法的參數(shù)估計精度最高。

4 結(jié)論

(1)文中提出的以累積誤差平方和最小為目標(biāo)的威布爾分布參數(shù)估計法考慮樣本擬合分布函數(shù)與經(jīng)驗分布值的誤差，可彌補經(jīng)典最小二乘法與加權(quán)最小二乘法的不足。

(2)引入粒子群算法進行基于累積誤差平方和最小的威布爾分布參數(shù)估計，避免大量數(shù)學(xué)推導(dǎo)過程，原理清晰易于理解，便于推廣應(yīng)用。

(3)文中將研究方法結(jié)合參考文獻的3種樣本數(shù)據(jù)進行應(yīng)用，K-S檢驗值與均方根誤差值表明，基于累積誤差平方和最小的威布爾分布參數(shù)估計方法是合理與有效的。