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

(數(shù)控裝備可靠性教育部重點(diǎn)實(shí)驗(yàn)室∥吉林大學(xué) 機(jī)械與航空航天工程學(xué)院，吉林 長(zhǎng)春 130022)

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

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

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

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

1 威布爾模型

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

(1)

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

對(duì)式(1)兩邊連取兩次自然對(duì)數(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ù)估計(jì)

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

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

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

(5)

具體計(jì)算步驟如下：

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

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

均方根誤差計(jì)算公式為

(6)

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

2.2 粒子群算法求解流程

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

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

(7)

(8)

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

粒子群算法的步驟如下：

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

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

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

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

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

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

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

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

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

(9)

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

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

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

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

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

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

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

(11)

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

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

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

(12)

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

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

(13)

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

(14)

3 演算驗(yàn)證

3.1 完全樣本

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

表1 軸承壽命時(shí)間與經(jīng)驗(yàn)分布值

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

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

表3 均方根誤差和K-S檢驗(yàn)值

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

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

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

圖1 完全樣本可靠度曲線

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

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

表4 元件失效時(shí)間與經(jīng)驗(yàn)分布值

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

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

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

表6 均方根誤差與K-S檢驗(yàn)值

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

圖2 定時(shí)截尾樣本可靠度曲線

表7 樣品失效時(shí)間與經(jīng)驗(yàn)分布值

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

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

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

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

表9 均方根誤差與K-S檢驗(yàn)值

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

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

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

4 結(jié)論

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

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

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