楊兆軍，楊川貴，陳 菲，郝慶波，鄭志同，王 松

(1.吉林大學(xué) 機械科學(xué)與工程學(xué)院，長春130022;2.吉林大學(xué) 機械工業(yè)數(shù)控裝備可靠性技術(shù)重點實驗室，長春130022;3.空軍航空大學(xué) 飛行器與動力系，長春130022)

0 引 言

準(zhǔn)確地對加工中心可靠性模型進(jìn)行參數(shù)估計，有利于正確評估加工中心的可靠性水平，對加工中心制造企業(yè)實施可靠性增長具有重要的意義［1］?？煽啃怨烙嫷臏?zhǔn)確程度與所選擇的可靠性模型及其參數(shù)估計方法有關(guān)。已有研究結(jié)果表明，威布爾(Weibull)模型能較好地描述加工中心的可靠性分布，其準(zhǔn)確性依賴于對分布參數(shù)的估計精度。在參數(shù)估計方面，常用的方法有最小二乘法(LSM)、最大似然估計(MLM)［2］和貝葉斯方法［3］。其中，前兩種是基于經(jīng)驗風(fēng)險最小原則發(fā)展起來的經(jīng)典統(tǒng)計方法，但是在樣本數(shù)量較少的情況下，上述方法很難獲得可靠性模型的最優(yōu)參數(shù)。貝葉斯方法在小樣本實驗數(shù)據(jù)條件下體現(xiàn)出其優(yōu)勢，但其估計精度取決于所選擇的先驗分布。支持向量回歸(Support vector regression，SVR)模型［4］在小樣本情況下有很好的歸納能力，與性能優(yōu)越的優(yōu)化算法進(jìn)行耦合能得到比較準(zhǔn)確的參數(shù)估計值。

目前為止，鮮有文獻(xiàn)將支持向量回歸模型用于可靠性模型的參數(shù)估計。為此，本文嘗試將其用于估計加工中心可靠性模型的參數(shù)，并將其與本文提出的改進(jìn)的局部最優(yōu)粒子群優(yōu)化算法(Improved local best particle swarm optimization，Improved lbest PSO)耦合，實現(xiàn)支持向量回歸模型的參數(shù)優(yōu)化。Improved lbest PSO 算法在局部最優(yōu)粒子群優(yōu)化算法(lbest PSO)的基礎(chǔ)上引入了變異操作和自適應(yīng)調(diào)整慣性因子，以抑制算法的早熟并提高算法的全局最優(yōu)解搜索能力。

1 加工中心可靠性模型

加工中心是一種復(fù)雜的可修復(fù)系統(tǒng)。前期研究工作表明:威布爾模型能較好地描述加工中心的可靠性分布，同時加工中心的故障間隔時間是其可靠性水平的體現(xiàn)［5］。根據(jù)威布爾模型，加工中心故障間隔時間的密度函數(shù)(即故障密度函數(shù))可用式(1)表示:

式中:λ 為比例參數(shù);k 為形狀參數(shù)。

同時，加工中心的故障分布函數(shù)和可靠度函數(shù)(即可靠性模型)分別如式(2)和式(3)所示:

平均故障間隔時間(MTBF)是加工中心可靠性評價的重要指標(biāo)之一，平均故障間隔時間越長，則加工中心越可靠。其計算公式為:

通過式(3)可以看出，加工中心可靠性模型的準(zhǔn)確程度取決于其參數(shù)λ 和k 的估計精度。

2 可靠性模型參數(shù)估計方法

2.1 故障經(jīng)驗分布函數(shù)估計

為了估計可靠性模型的參數(shù)，需要對故障分布函數(shù)做出評估。假設(shè)加工中心提供了n 個故障間隔時間，依次為t1≤t2≤…≤tn，其中i 為故障間隔時間序號。當(dāng)處理不完全數(shù)據(jù)時，平均秩次法［6］能夠提供高精度的故障經(jīng)驗分布，因此本文采用該方法對故障經(jīng)驗分布函數(shù)進(jìn)行計算。該加工中心的故障經(jīng)驗分布函數(shù)計算過程如下:

式中:Ai為第i 個故障的平均秩次，并且A1=1。

2.2 用于參數(shù)λ 和k 估計的支持向量回歸模型

為了估計參數(shù)λ 和k，首先對式(2)進(jìn)行線性化處理，得出如下等式:

令y=lnln1/(1－F(t))，w=k，φ(x)=lnt，b=－klnλ，則式(7)可表示如下:

根據(jù)結(jié)構(gòu)風(fēng)險最小原則，SVR 模型通過對如下?lián)p失函數(shù)最小化以獲得w、b 的估計值。根據(jù)文獻(xiàn)［7-8］，用于估計參數(shù)λ、k 的支持向量回歸模型如下:

約束條件:

約束條件:

根據(jù)Lagrange 乘子法則，參數(shù)w 的計算公式如下:

利用Karush-Kuhn-Tucker(KKT)［7］條件，參數(shù)b 可通過如下公式計算:

綜上可知，式(3)中的參數(shù)可通過如下公式得出:

圖1 描述了參數(shù)C 和ε 對SVR 模型精度的影響。由圖1 可知，上述SVR 模型的估計精度完全受參數(shù)C、ε 的影響，因此SVR 模型的參數(shù)優(yōu)化是其面臨的首要任務(wù)。

2.3 支持向量回歸模型參數(shù)優(yōu)化

2.3.1 lbestPSO 算法概述

在本文中，lbestPSO［9-10］被用來優(yōu)化參數(shù)C、ε。粒子i 的位置向量和速度向量分別定義為xi=［xi1，xi2］和vi=［vi1，vi2］，其第一維和第二維分別對應(yīng)于參數(shù)C 和ε。同時，粒子的位置和速度更新公式如下:

圖1 參數(shù)C、ε 對SVR 模型精度的影響Fig.1 Influence of the parameters ε and C on the accuracy of SVR model

式中:γ 為慣性因子;c1、c2為學(xué)習(xí)因子;rand1、rand2為0 到1 之間服從均勻分布的隨機數(shù);為粒子i 第k 次迭代后的最優(yōu)位置;為粒子i 第k 次迭代后在其鄰域中的最優(yōu)位置。另外，粒子的全局最優(yōu)位置由Gbestk表示。

2.3.2 lbestPSO 算法改進(jìn)

為了抑制粒子早熟，本文引入變異操作:若粒子i 的位置向量xi與目前發(fā)現(xiàn)的全局最優(yōu)位置Gbestk的空間距離小于給定的最小距離Dmin時，對該粒子進(jìn)行如下變異操作［11］:

式中:rand3是0 到1 之間服從均勻分布的隨機數(shù)。

同時為了提高算法的全局搜索能力，引入一種非線性自適應(yīng)調(diào)節(jié)慣性因子［12］，如下式所示:

式中:f(·)為粒子的適應(yīng)度函數(shù);γmax和γmin分別為γ 的最大值和最小值;分別為第k 次迭代所有粒子的適應(yīng)度函數(shù)的最小值和平均值。

上述改進(jìn)后的算法在本文中稱為Improved lbest PSO 算法。

2.3.3 基于Improved lbest PSO 的SVR 模型參數(shù)選擇

利用參數(shù)(C，ε)，本文將Improved lbest PSO算法和SVR 模型耦合，耦合后的模型簡稱為Improved lbest PSO-SVR 模型。同時為了得到最優(yōu)的C、ε，本文將模型的均方根誤差作為Improved lbest PSO 算法的目標(biāo)函數(shù)(即粒子的適應(yīng)度函數(shù))，其定義如下:

式中:r 為樣本的數(shù)量;di為y 的實際值;fi為SVR模型對y 的估計值;RMSE 為均方根誤差。其中di和fi的表達(dá)式為:

支持向量回歸模型的參數(shù)可通過求解下面的優(yōu)化問題得出:

參數(shù)C、ε 的具體優(yōu)化過程如下所示。

Step1 參考表1，初始化Improved lbest PSO算法中的參數(shù);

Step2 隨機產(chǎn)生N1個粒子的位置和速度向量;

Step3 將每個粒子用于支持向量模型學(xué)習(xí)，確定每個粒子相應(yīng)的(w，b);

Step4 利用式(17)計算每個粒子的適應(yīng)度;

Step6 利用式(16)計算粒子群中每個粒子的自適應(yīng)調(diào)節(jié)慣性因子γ;

Step7 利用式(13)更新粒子群中每個粒子的速度向量，并將粒子的速度限制在規(guī)定的速度范圍內(nèi);

Step8 利用式(14)更新粒子群中每個粒子的位置向量，并將粒子的位置限制在規(guī)定的位置范圍內(nèi);

Step9 判斷各粒子是否滿足變異條件:如滿足則按照上述變異方式對該粒子進(jìn)行變異處理;

Step10 判斷是否滿足結(jié)束條件:若滿足則停止優(yōu)化過程，否則返回Step3。

3 數(shù)值仿真

假設(shè)加工中心的故障密度函數(shù)為:

首先，本文根據(jù)上述故障密度函數(shù)產(chǎn)生數(shù)量為20、30、50 的故障間隔時間樣本集。為了能夠?qū)χС窒蛄炕貧w模型的參數(shù)進(jìn)行優(yōu)化，本文將樣本集(含l 個元素)分為兩類，分別是訓(xùn)練組、校驗組［13］。校驗組樣本(含r 個元素)用于估計參數(shù)C、ε，訓(xùn)練組樣本(含l－r 個元素)用于估計可靠性模型的參數(shù)，其中校驗組的樣本數(shù)常為總樣本數(shù)的1/4。為了使Improved lbest PSO 算法能夠更高效地工作，其所需參數(shù)初始值設(shè)置如下:c1=c2=1.4962;wmax=0.9;wmin=0.4;Cmax=300;Cmin=1×10－5;εmax=0.5;εmin=1×10－5;Dmin=1;N1=40;迭代次數(shù)為300;誤差為1×10－5。

粒子群多樣性定義如下式所示，用以衡量改進(jìn)前、后算法的性能。

式中:N 為粒子數(shù)量。

為了比較Improved lbest PSO-SVR、LSM、MLE、局部最優(yōu)粒子群優(yōu)化算法優(yōu)化的支持向量回歸模型(lbest PSO-SVR)和遺傳算法優(yōu)化的支持向量回歸模型(GA-SVR)的性能，MTBF 估計值的相對誤差被作為各方法的綜合評價指標(biāo)，其計算公式如下:

式中:MTBF 通過式(4)計算。當(dāng)計算MTBFreal時，f(t)為實際故障密度函數(shù)(見式(19));當(dāng)計算MTBFest時，f(t)為估計得到的故障密度函數(shù)。

在上述基礎(chǔ)上，本文用LS、MLE、lbest PSOSVR、Improved lbest PSO-SVR 和GA-SVR 對上述加工中心可靠性進(jìn)行估計。不同樣本下，Improved lbest PSO、lbest PSO 和遺傳算法對參數(shù)C、ε 的優(yōu)化結(jié)果如表1 所示。從表1 可知，Improved lbest PSO 的均方根誤差小于lbest PSO和遺傳算法。

表1 參數(shù)C、ε 的優(yōu)化結(jié)果Table 1 Optimization results of the parameters C and ε

圖2 為適應(yīng)度函數(shù)和粒子群多樣性的變化情況。從圖2(a)可知，Improved lbest PSO 算法的收斂速度比lbest PSO 算法更快;從圖2(b)可知，Improved lbest PSO 算法保持著較高的多樣性，而lbest PSO 算法的多樣性則逐漸下降為0，這表明Improved lbest PSO 算法得到的參數(shù)更可信。上述現(xiàn)象證明了本文引入的變異操作、自適應(yīng)慣性因子使得Improved lbest PSO 算法的性能優(yōu)于lbest PSO 和遺傳算法。

表2 為LSM、MLE、lbest PSO-SVR、GA-SVR和Improved lbest PSO-SVR 對上述加工中心故障密度函數(shù)的參數(shù)估計結(jié)果。各方法估計誤差都在10%以下，但I(xiàn)mproved lbest PSO-SVR 模型的相對誤差小于5%，并且小于其他方法的估計誤差。

圖2 改進(jìn)前、后lbest PSO 算法的性能對比Fig.2 Performances comparison between the lbest PSO algorithm and improved lbest PSO algorithm

表2 LSM、MLE、SVR 的估計結(jié)果Table 2 Estimation results of LSM，MLE，SVR

圖3(a)、圖4(a)、圖5(a)分別表示l=20，l=30 和l=50 時，Improved lbest PSO-SVR 模型的估計結(jié)果，其中，最佳SVR 模型所需的支持向量機數(shù)量分別是13、4、12，并且得到的回歸直線與數(shù)據(jù)擬合良好。圖3(b)、圖4(b)、圖5(b)表示LSM、MLE、Improved lbest PSO-SVR 得到的故障時間間隔分布函數(shù)。從圖中可知，上述各方法都能較好地估計故障分布函數(shù)的參數(shù)，但是Improved lbest PSO-SVR 估計的故障分布函數(shù)與真實的故障分布函數(shù)最接近。因此，Improved lbest PSO-SVR 性能優(yōu)于LSM、MLE、lbest PSOSVR、GA-SVR。

圖3 樣本數(shù)量l=20 時，各方法的估計結(jié)果Fig.3 In the l=20 case，estimation results of every method

圖4 樣本數(shù)量l=30 時，各方法的估計結(jié)果Fig.4 In the l=30 case，the estimation results of every method

圖5 樣本數(shù)量l=50 時，各方法的估計結(jié)果Fig.5 In the l=50 case，the estimation results of every methods

4 應(yīng) 用

某系列加工中心的44 條故障間隔時間數(shù)據(jù)如圖6(a)所示。根據(jù)3.2 節(jié)的線性化處理方式，故障間隔時間數(shù)據(jù)的線性處理結(jié)果如圖6(b)所示。從圖6(b)可知，線性處理后的數(shù)據(jù)基本在一條直線上，因此可判斷該系列加工中心的故障間隔時間服從威布爾分布。

圖6 某加工中心的故障間隔時間Fig.6 Fault time intervals of a actual machining center

采用Improved lbest PSO-SVR 模型對該加工中心的可靠性模型進(jìn)行參數(shù)估計，其中Improved lbest PSO 算法的參數(shù)設(shè)置參見第3 節(jié)。此外，將故障間隔時間分為兩組:訓(xùn)練組(含33 個元素)和校驗組(含11 個元素)。Improved lbest PSOSVR 模型的估計結(jié)果如下:w =0.988;C =0.186;b=－6.521;ε=0.064;k=0.988;MTBF/h=735.693;λ/h=734.288。

圖7 為利用Improved lbest PSO-SVR 模型估計得到的回歸直線。從圖7 可知，最佳SVR 模型所需的支持向量機數(shù)量為24，同時SVR 模型能自動地剔除故障樣本中的異常樣本。由Improved lbest PSO-SVR 得到的該加工中心的故障分布函數(shù)和可靠度函數(shù)的估計結(jié)果如圖8 所示。同時利用Improved lbest PSO-SVR 模型估計出的參數(shù)值可以確定該加工中心的故障密度函數(shù)。利用式(5)可得，該機床的MTBF 估計值為735.693 h。這一估計結(jié)果將為該加工中心的可靠性增長提供基礎(chǔ)依據(jù)。

圖7 Improved lbest PSO-SVR 估計結(jié)果Fig.7 Estimation results of improved lbest PSO-SVR

圖8 某加工中心的可靠性估計結(jié)果Fig.8 Reliability assessment results of a actual machining center

5 結(jié) 論

(1)提出了改進(jìn)的粒子群優(yōu)化支持向量回歸模型，并將其應(yīng)用于可靠性領(lǐng)域，成功地對加工中心可靠性模型的參數(shù)進(jìn)行了估計，其估計精度優(yōu)于最小二乘法、最大似然估計法、原局部最優(yōu)粒子群優(yōu)化算法和遺傳算法優(yōu)化的支持向量回歸模型，由此建立的加工中心可靠性模型具有較高的可靠性評估精度。

(2)引入的變異操作、自適應(yīng)調(diào)整慣性因子提高了改進(jìn)的局部最優(yōu)粒子群優(yōu)化算法的全局尋優(yōu)能力和收斂速度，使其優(yōu)于原局部最優(yōu)粒子群優(yōu)化算法和遺傳算法。同時改進(jìn)的局部最優(yōu)粒子群優(yōu)化算法優(yōu)化后的支持向量回歸模型具有較高的參數(shù)估計精度。

［1］楊兆軍，陳傳海，陳菲，等.數(shù)控機床可靠性技術(shù)的研究進(jìn)展［J］.機械工程學(xué)報，2013，49(20):130-139.Yang Zhao-jun，Chen Chuan-hai，Chen Fei，et al.Progress in the research of reliability technology of machine tools［J］.Journal of Mechanical Engineering，2013，49(20):130-139.

［2］Balakrishnan N，Kateri M.On the maximum likelihood estimation of parameters of Weibull distribution based on complete and censored data［J］.Statistics＆Probability Letters，2008，78(17):2971-2975.

［3］Wu H C.Fuzzy reliability estimation using Bayesian approach［J］.Computers ＆ Industrial Engineering，2004，46(3):467-493.

［4］Vapnik V.Statistical Learning Theory［M］.New York:Wiley，1998.

［5］楊兆軍，李小兵，許彬彬，等.加工中心時間動態(tài)可靠性建模［J］.機械工程學(xué)報，2012，48(2):16-22.Yang Zhao-jun，Li Xiao-bing，Xu Bin-bin，et al.Time dynamic reliability modelling of machining center［J］.Chinese Journal of Mechanical Engineering，2012，48(2):16-22.

［6］Bryan D.The Weibull Analysis Handbook［M］.ASQ Quality Press，2006.

［7］Smola Alex J，Bernhard S.A tutorial on support vector regression［J］.Statistics and Computing，2004，14(3):199-222.

［8］吳堅，趙陽，何睿.基于支持向量機回歸算法的電子機械制動傳感器系統(tǒng)故障診斷［J］.吉林大學(xué)學(xué)報:工學(xué)版，2013，43(5):1178-1183.Wu Jian，Zhao Yang，He Rui.Fault detection and diagnosis of EMB sensor system based on SVR［J］.Journal of Jilin University(Engineering and Technology Edition)，2013，43(5):1178-1183.

［9］Kennedy J，Mendes R.Population structure and particle swarm performance［C］∥Proceedings of the 2002 Congress on Evolutionary Computation，Honolulu HI，USA，2002:1671-1676.

［10］寇曉麗，劉三陽.基于模擬退火的粒子群算法求解約束優(yōu)化問題［J］.吉林大學(xué)學(xué)報:工學(xué)版，2007，37(1):136-140.Kou Xiao-li，Liu San-yang.Particle swarm algorithm based on simulated annealing to solve constrained optimization［J］.Journal of Jilin University(Engineering and Technology Edition)，2007，37(1):136-140.

［11］王凌，劉波.微粒群優(yōu)化與調(diào)度算法［M］.北京:清華大學(xué)出版社，2008.

［12］Engelbrecht A P.Fundamentals of Computational Swarm Intelligence［M］.Hoboken:John Wiley ＆Sons，2006.

［13］Lins I D，Moura M C，Zio E，et al.A particle swarm‐optimized support vector machine for reliability prediction［J］.Quality and Reliability Engineering International，2012，28(2):141-158.