郭 森，王大為，張紹偉，姚永超

(上海機(jī)電工程研究所，上海 201109)

0 引言

隨著裝備復(fù)雜化程度不斷提升，綜合保障成為影響其作戰(zhàn)使用的重要因素，而通過裝備全壽命周期數(shù)據(jù)的應(yīng)用分析挖掘出故障及使用維修規(guī)律，確定預(yù)防性維修間隔，是提升其綜合保障能力的關(guān)鍵環(huán)節(jié)[1-2]。其中，根據(jù)裝備失效維護(hù)數(shù)據(jù)確定其壽命分布是重要途徑。威布爾分布是一種連續(xù)型概率分布模型，它能夠描述各種類型機(jī)械零部件失效數(shù)據(jù)的分布規(guī)律，在壽命數(shù)據(jù)分析、可靠評(píng)估和預(yù)測(cè)等方面得到了一定程度的應(yīng)用[3]。對(duì)于復(fù)雜裝備而言，失效數(shù)據(jù)往往不會(huì)對(duì)各失效模式進(jìn)行區(qū)分，用簡(jiǎn)單威布爾分布并不能有效地對(duì)設(shè)備的壽命進(jìn)行描述[4]。多重威布爾混合模型具有良好的擬合性能，更適合于復(fù)雜裝備失效規(guī)律的描述，但也會(huì)引入更多的未知參數(shù)，增加了模型參數(shù)估計(jì)算法的復(fù)雜度[5]。

常見的威布爾分布參數(shù)估計(jì)方法包括圖解法[6]、期望極大化(EM)法、極大似然估計(jì)(MLE)法[7]和貝葉斯法[8]等。其中：圖解法是粗略的估計(jì)，精度較低；EM算法和MLE算法需求解超越方程，在大樣本下估計(jì)精度較高，但是求解過程比較復(fù)雜，有時(shí)不易收斂，且面對(duì)截尾小樣本時(shí)偏差較大；貝葉斯法可以充分利用專家經(jīng)驗(yàn)、設(shè)計(jì)參數(shù)等信息提高參數(shù)估計(jì)的準(zhǔn)確性，并且隨著使用過程中新信息的出現(xiàn)不斷更新結(jié)果，是解決小樣本參數(shù)估計(jì)的有效方法，但對(duì)設(shè)備先驗(yàn)信息要求較高。

為了取得更精確的壽命分布模型，并減少對(duì)先驗(yàn)信息的依賴性，許多學(xué)者引入了非線性擬合優(yōu)化方法來對(duì)威布爾模型參數(shù)進(jìn)行估計(jì)。文獻(xiàn)[9]將兩重二參數(shù)威布爾混合分布參數(shù)估計(jì)轉(zhuǎn)化為非線性最小二乘優(yōu)化問題，并采用Quasi-Newton法進(jìn)行迭代求解，能夠得到更加精確的結(jié)果，但理論要求比較高；文獻(xiàn)[10]采用灰色模型與支持向量機(jī)相結(jié)合的方式，解決了小樣本情況下三參數(shù)威布爾分布參數(shù)估計(jì)精度較低的問題；文獻(xiàn)[11]采用混合粒子群(PSO)算法對(duì)模型似然函數(shù)進(jìn)行求解，避免了傳統(tǒng)極大似然估計(jì)方法的短板，一定程度上解決了威布爾混合分布模型的參數(shù)估計(jì)問題；文獻(xiàn)[12]引入模擬退火(SA)算法，直接對(duì)威布爾分布模型進(jìn)行參數(shù)估計(jì)，取得了較好的結(jié)果。

針對(duì)威布爾混合分布模型復(fù)雜、參數(shù)估計(jì)難度大的問題，本文提出采用基于自適應(yīng)改進(jìn)模擬退火粒子群(SA-PSO)算法加以解決。將威布爾混合分布模型參數(shù)估計(jì)視為一個(gè)非線性優(yōu)化問題，采用粒子群算法求解，并對(duì)粒子速度與位置更新進(jìn)行自適應(yīng)改進(jìn)，在全局最優(yōu)值的選擇上引入模擬退火機(jī)制，進(jìn)一步改善粒子群算法的全局搜索能力。將改進(jìn)的算法應(yīng)用到某型柴油機(jī)噴油器的失效數(shù)據(jù)分析中，利用圖解法對(duì)參數(shù)初值進(jìn)行粗估，用于自適應(yīng)SA-PSO優(yōu)化算法的種群初始化，對(duì)威布爾混合分布進(jìn)行參數(shù)估計(jì)，結(jié)果表明該算法能夠有效提高參數(shù)估計(jì)的精度和效率。

1 模擬退火粒子群算法及其改進(jìn)

1.1 PSO算法基本原理

PSO算法是一種用于非線性函數(shù)優(yōu)化的集群智能優(yōu)化算法[13]，它在解決大規(guī)模非線性問題中具有更高的優(yōu)化效率和更好的優(yōu)化結(jié)果。該算法將各個(gè)備選解稱為“粒子”，全部的粒子集合稱作“種群”。PSO算法首先生成初始種群，即在可行解空間中隨機(jī)初始化生成一群粒子，并根據(jù)要優(yōu)化的目標(biāo)函數(shù)確定與之相對(duì)應(yīng)的適應(yīng)度值。優(yōu)化過程就是初始種群在可行解空間內(nèi)部飛行，并尋找最大或最小適應(yīng)度值的過程，粒子飛行的方向與距離由速度決定。

其數(shù)學(xué)過程描述如下所述。

vid(t+1)=vid(t)+c1r1(pid(t)-xid(t))+

c2r2(gd(t)-xid(t))

(1)

xid(t+1)=xid(t)+vid(t+1)

(2)

1≤i≤N,1≤d≤D

c1、c2為正數(shù)，稱作加速因子,c1調(diào)節(jié)粒子飛向自身最好位置方向上的步長，c2調(diào)節(jié)粒子飛向全局最好位置方向上的步長;r1、r2為[0,1]的隨機(jī)數(shù)。

通過比較更新前后各粒子的適應(yīng)度值，對(duì)單個(gè)粒子最優(yōu)位置和全局最優(yōu)位置進(jìn)行更新，直到滿足迭代終止條件。

1.2 PSO算法的自適應(yīng)改進(jìn)

與其他群體智能算法相比，PSO算法概念簡(jiǎn)單，易于編程實(shí)現(xiàn)求解，但在優(yōu)化過程中同樣會(huì)面臨易于陷入局部極值點(diǎn)、搜索精度不高、收斂速度慢的問題。為此，引入慣性權(quán)重因子ω，對(duì)粒子飛行速度進(jìn)行調(diào)節(jié)，提高PSO算法的收斂速度及全局搜索能力，使其能夠快速得到全局最優(yōu)結(jié)果。慣性權(quán)重因子ω對(duì)速度的調(diào)節(jié)原理為

vid(t+1)=ωvid(t)+c1r1(pid(t)-xid(t))+

c2r2(gd(t)-xid(t))

(3)

由式(3)可以看出，引入的ω代表的含義主要是上一代粒子狀態(tài)對(duì)當(dāng)前粒子狀態(tài)的影響程度，ω取值越大，則上一代粒子狀態(tài)對(duì)與當(dāng)前粒子的狀態(tài)影響越大；反之，影響則越小。同時(shí)，算法的全局搜索能力也會(huì)隨著ω的增大而提高。因此，通過控制慣性權(quán)重因子ω的大小，可以改進(jìn)粒子群算法的全局搜索能力，當(dāng)ω值取為常數(shù)1時(shí)，式(3)即退化為基本粒子群算法的速度更新公式。

為了進(jìn)一步提升優(yōu)化性能，并加快粒子群算法的收斂速度，引入距離控制因子C(t)對(duì)慣性權(quán)重因子及加速系數(shù)進(jìn)行自適應(yīng)調(diào)整[14]，原理為

(4)

D為各粒子到全局最優(yōu)值的平均距離；max(D)為至今平均距離的最大值；C(t)為距離控制因子，即各粒子到全局最優(yōu)點(diǎn)的平均距離與平均距離最大值的比值。

1.3 基于模擬退火的自適應(yīng)粒子群優(yōu)化策略

模擬退火的思想來自常用于金屬材料提升性能的熱處理工藝退火，通過模擬高溫下金屬逐漸降溫直至熱平衡的過程，引入概率接收準(zhǔn)則，并適當(dāng)?shù)乜刂茰囟葏?shù)的下降，反復(fù)迭代尋找能量最小值點(diǎn)，最后得到問題近似的全局最優(yōu)解[12]。

將模擬退火思想引入PSO算法的優(yōu)化過程，改善算法在訓(xùn)練過程精度降低容易發(fā)散的問題。在最優(yōu)粒子的取舍中加入Metropolis準(zhǔn)則，在接收更優(yōu)粒子的同時(shí)以一定的概率接收較差粒子，能夠從局部極值中跳出，隨著溫度的降低逐漸收斂到全局最優(yōu)位置。

基于模擬退火的自適應(yīng)粒子群優(yōu)化過程如下：

a.初始化種群中各粒子的位置和速度。

b.計(jì)算各粒子的適應(yīng)度值，將當(dāng)前各粒子位置和適應(yīng)度值存儲(chǔ)在各粒子的pi中，將所有Pb中最優(yōu)個(gè)體位置和適應(yīng)度值存儲(chǔ)在pb中。

c.根據(jù)式(5)確定初始溫度t0，即

t0=f(pb)/ln5

(5)

d.根據(jù)式(6)確定當(dāng)前溫度下各pi的適配值，即

(6)

f.計(jì)算各粒子新的適應(yīng)度值，更新各粒子的pi值及群體的pb值。

g.采用衰減系數(shù)法進(jìn)行退溫操作，即

tk+1=λtk

(7)

h.若滿足終止條件搜索停止，輸出結(jié)果，否則返回步驟d繼續(xù)搜索。

2 威布爾分布及其參數(shù)估計(jì)方法

2.1 二參數(shù)威布爾分布

由于威布爾分布既能描述失效率上升的產(chǎn)品壽命，又能描述失效率下降的產(chǎn)品壽命，因此威布爾分布是可靠性工程中適于描述壽命分布規(guī)律的統(tǒng)計(jì)模型之一。

若設(shè)備壽命t服從二參數(shù)威布爾分布，則其累積失效概率密度函數(shù)、可靠度和失效率為：

(8)

(9)

(10)

β、θ分別為威布爾模型的形狀參數(shù)、尺度參數(shù)。

2.2 威布爾混合分布模型

復(fù)雜機(jī)械系統(tǒng)包含多個(gè)零部件，每個(gè)零部件的故障可能是在多種失效機(jī)理共同作用下發(fā)生的，而且在不同的壽命階段，不同的失效機(jī)理對(duì)系統(tǒng)的失效起主導(dǎo)作用。上述情況下，采集到的可靠性試驗(yàn)數(shù)據(jù)在威布爾概率圖(WWP)上表現(xiàn)為曲線，如果用標(biāo)準(zhǔn)的單威布爾分布模型來描述時(shí)，會(huì)出現(xiàn)較大的誤差，威布爾混合分布提供了更符合實(shí)際的壽命數(shù)據(jù)擬合模型[3]。

如果1個(gè)總體由n個(gè)子體組成，設(shè)各子體的概率密度函數(shù)分別為f1(t)，f2(t),…,fn(t)，各子體的混合權(quán)重分別為p1(t),p2(t),…,pn(t)，權(quán)重之和為1，則總體概率密度函數(shù)f(t)為

(11)

若各子體都服從威布爾分布，即

(12)

βi為第i個(gè)威布爾分布的形狀參數(shù)；θi為第i個(gè)威布爾分布的尺度參數(shù)；pi為第i個(gè)威布爾分布的權(quán)重，它表示第i種失效原因或失效模式的失效數(shù)據(jù)在產(chǎn)品總的失效數(shù)據(jù)中所占的比例。

2.3 基于自適應(yīng)SA-PSO算法的參數(shù)估計(jì)優(yōu)化

對(duì)于單一威布爾分布而言，可以經(jīng)過簡(jiǎn)單的數(shù)學(xué)變換將樣本數(shù)據(jù)的非線性關(guān)系轉(zhuǎn)換線性關(guān)系，并采取線性回歸方法對(duì)數(shù)據(jù)進(jìn)行擬合。但威布爾混合分布模型較為復(fù)雜，一般的數(shù)學(xué)變換無法將模型向線性問題進(jìn)行轉(zhuǎn)換，因此采用非線性優(yōu)化方式對(duì)其參數(shù)估計(jì)問題進(jìn)行求解[15]。針對(duì)兩重威布爾混合分布模型，基于誤差平方和最小的思想，建立多參數(shù)威布爾混合分布的非線性最小二乘(NLS)優(yōu)化模型。

將樣本容量為m的失效數(shù)據(jù)按壽命大小進(jìn)行順序排列(t1,t2,…,tm)，采用中位秩公式計(jì)算累計(jì)失效概率的觀測(cè)值[3]，如式(13)：

(13)

(14)

根據(jù)誤差平方和最小思想構(gòu)建威布爾混合模型參數(shù)非線性最小二乘優(yōu)化模型，如式(15)：

(15)

采用自適應(yīng)改進(jìn)的模擬退火粒子群算法對(duì)威布爾混合分布參數(shù)優(yōu)化模型進(jìn)行求解，流程如圖1所示。在威布爾分布的參數(shù)估計(jì)方法中，圖解法參數(shù)估計(jì)值的精度較低，往往不能滿足實(shí)際使用需求，但是可以作為其他方法的迭代初值。因此，先對(duì)設(shè)備失效數(shù)據(jù)用圖解法進(jìn)行求解，得到迭代初始值，之后按照基于模擬退火的自適應(yīng)粒子群優(yōu)化策略進(jìn)行迭代優(yōu)化，直到滿足終止條件。

圖1 基于自適應(yīng)SA-PSO的威布爾混合分布參數(shù)估計(jì)流程

3 實(shí)例分析

為了驗(yàn)證本文所提方法的準(zhǔn)確性，采用文獻(xiàn)[3]中某型柴油機(jī)噴油器服從兩重威布爾混合分布的故障里程進(jìn)行擬合優(yōu)化，其數(shù)據(jù)如表1所示。根據(jù)式(13)構(gòu)造樣本集，其可靠度函數(shù)為

(16)

不可靠度為

F(t)=1-R(t)

(17)

優(yōu)化目標(biāo)函數(shù)為

(18)

圖解法粗估值為p=0.620 5，θ1=9 996.6，β1=1.408 5，θ2=2 250.3，β2=11.007 9。在此結(jié)果基礎(chǔ)上構(gòu)造自適應(yīng)SA-PSO優(yōu)化算法的初始種群，以優(yōu)化目標(biāo)函數(shù)式(18)作為適應(yīng)度函數(shù)按圖1所示流程進(jìn)行優(yōu)化。同時(shí)，為了驗(yàn)證自適應(yīng)SA-PSO方法的優(yōu)勢(shì)，將其結(jié)果與標(biāo)準(zhǔn)PSO算法、自適應(yīng)PSO算法結(jié)果進(jìn)行對(duì)比。

表1 某型柴油機(jī)噴油器的故障里程

3.1 自適應(yīng)SA-PSO參數(shù)優(yōu)化性能分析

不同優(yōu)化算法對(duì)模型解決參數(shù)估計(jì)問題的優(yōu)化過程如圖2所示。從圖2可以看出，標(biāo)準(zhǔn)PSO算法在迭代704步之后適應(yīng)度值收斂至0.031 70，自適應(yīng)PSO算法在迭代至250步時(shí)適應(yīng)度值收斂至0.030 60，說明二者在收斂精度上相差不大，但后者收斂速度明顯提高；與前2種算法相比，自適應(yīng)SA-PSO算法在迭代過程中，適應(yīng)度值于第113步便快速達(dá)到局部極值點(diǎn)0.033 25，并隨后在第527步、第780步產(chǎn)生突變，最終收斂到0.021 00，這表明，引入模擬退火機(jī)制并經(jīng)權(quán)重自適應(yīng)更新的粒子群算法能夠有效跳出局部極值點(diǎn)收斂到更低的適應(yīng)度值，具有更優(yōu)的全局搜索能力以及更快的收斂速度。綜合分析，3種算法中，自適應(yīng)SA-PSO算法優(yōu)化性能最優(yōu)，自適應(yīng)PSO算法次之，標(biāo)準(zhǔn)PSO算法由于未經(jīng)過任何改進(jìn)，收斂速度較慢且易于陷入局部最優(yōu)，優(yōu)化性能最弱。

圖2 不同算法參數(shù)優(yōu)化過程

3.2 威布爾分布參數(shù)估計(jì)結(jié)果分析

采用圖解法、基于Levenberg-Marquardt (L-M)的非線性最小二乘法(NLS)、標(biāo)準(zhǔn)PSO算法、自適應(yīng)PSO算法和自適應(yīng)SA-PSO算法，得到的模型參數(shù)估計(jì)結(jié)果及相對(duì)均方根誤差如表2所示。采用基于L-M的NLS方法得到的模型相對(duì)均方根誤差為0.089 1，通過智能優(yōu)化算法擬合的模型誤差都要更小，而且隨著算法的改進(jìn)，誤差也變得越來越小，其中經(jīng)過自適應(yīng)改進(jìn)的SA-PSO算法結(jié)果最好，擬合模型的相對(duì)均方根誤差僅為0.024 8，與優(yōu)化性能分析結(jié)果相符。

表2 不同方法參數(shù)估計(jì)結(jié)果

根據(jù)各方法參數(shù)估計(jì)結(jié)果畫出相應(yīng)的模型曲線，如圖3所示。圖中能直觀看出噴油器失效數(shù)據(jù)中位秩點(diǎn)比較均勻的分布于相應(yīng)擬合曲線兩邊。

結(jié)果表明，與傳統(tǒng)方法相比，基于群體智能優(yōu)化算法擬合的模型具有更好的精度，其中自適應(yīng)SA-PSO算法擬合精度最好。根據(jù)該算法結(jié)果得到某型柴油機(jī)噴油器的可靠度函數(shù)為

(19)

圖3 不同方法擬合曲線結(jié)果對(duì)比

4 結(jié)束語

本文對(duì)威布爾混合分布模型傳統(tǒng)參數(shù)估計(jì)方法進(jìn)行了分析。針對(duì)模型復(fù)雜、求解難度較大的問題，在殘差最小化思想基礎(chǔ)上建立了多重兩參數(shù)威布爾混合分布的非線性最小二乘優(yōu)化模型，并提出了一種基于自適應(yīng)SA-PSO算法加以求解，最后采用某型柴油機(jī)噴油器失效數(shù)據(jù)，對(duì)本文所構(gòu)造的參數(shù)估計(jì)模型及求解方法進(jìn)行了驗(yàn)證，可以得到如下結(jié)論：

a.與標(biāo)準(zhǔn)PSO算法及經(jīng)過自適應(yīng)改進(jìn)的PSO算法相比，通過引入模擬退火機(jī)制對(duì)優(yōu)化過程全局最優(yōu)值選取進(jìn)行控制，可以有效避免陷入局部最優(yōu)，同時(shí)可以提高算法的收斂速度。

b.與傳統(tǒng)用于解決參數(shù)估計(jì)的圖解法、非線性最小二乘法相比，采用基于自適應(yīng)SA-PSO算法求解威布爾混合分布模型參數(shù)的最小二乘估計(jì)值，能夠最大程度上減少先驗(yàn)信息缺乏的影響，有效簡(jiǎn)化了計(jì)算的復(fù)雜度，并且具有更高的求解精度。