基于PSO 算法和SVR 模型的加工中心可靠性模型參數估計

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

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

0 引 言

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

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

1 加工中心可靠性模型

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

式中:λ 為比例參數;k 為形狀參數。

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

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

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

2 可靠性模型參數估計方法

2.1 故障經驗分布函數估計

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

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

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

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

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

根據結構風險最小原則，SVR 模型通過對如下損失函數最小化以獲得w、b 的估計值。根據文獻［7-8］，用于估計參數λ、k 的支持向量回歸模型如下:

約束條件:

約束條件:

根據Lagrange 乘子法則，參數w 的計算公式如下:

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

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

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

2.3 支持向量回歸模型參數優化

2.3.1 lbestPSO 算法概述

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

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

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

2.3.2 lbestPSO 算法改進

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

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

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

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

上述改進后的算法在本文中稱為Improved lbest PSO 算法。

2.3.3 基于Improved lbest PSO 的SVR 模型參數選擇

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

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

支持向量回歸模型的參數可通過求解下面的優化問題得出:

參數C、ε 的具體優化過程如下所示。

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

Step2 隨機產生N1個粒子的位置和速度向量;

Step3 將每個粒子用于支持向量模型學習，確定每個粒子相應的(w，b);

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

Step6 利用式(16)計算粒子群中每個粒子的自適應調節慣性因子γ;

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

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

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

Step10 判斷是否滿足結束條件:若滿足則停止優化過程，否則返回Step3。

3 數值仿真

假設加工中心的故障密度函數為:

首先，本文根據上述故障密度函數產生數量為20、30、50 的故障間隔時間樣本集。為了能夠對支持向量回歸模型的參數進行優化，本文將樣本集(含l 個元素)分為兩類，分別是訓練組、校驗組［13］。校驗組樣本(含r 個元素)用于估計參數C、ε，訓練組樣本(含l－r 個元素)用于估計可靠性模型的參數，其中校驗組的樣本數常為總樣本數的1/4。為了使Improved lbest PSO 算法能夠更高效地工作，其所需參數初始值設置如下: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;迭代次數為300;誤差為1×10－5。

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

式中:N 為粒子數量。

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

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

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

表1 參數C、ε 的優化結果Table 1 Optimization results of the parameters C and ε

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

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

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

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

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

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

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

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

4 應 用

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

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

采用Improved lbest PSO-SVR 模型對該加工中心的可靠性模型進行參數估計，其中Improved lbest PSO 算法的參數設置參見第3 節。此外，將故障間隔時間分為兩組:訓練組(含33 個元素)和校驗組(含11 個元素)。Improved lbest PSOSVR 模型的估計結果如下: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 模型所需的支持向量機數量為24，同時SVR 模型能自動地剔除故障樣本中的異常樣本。由Improved lbest PSO-SVR 得到的該加工中心的故障分布函數和可靠度函數的估計結果如圖8 所示。同時利用Improved lbest PSO-SVR 模型估計出的參數值可以確定該加工中心的故障密度函數。利用式(5)可得，該機床的MTBF 估計值為735.693 h。這一估計結果將為該加工中心的可靠性增長提供基礎依據。

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

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

5 結 論

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

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

［1］楊兆軍，陳傳海，陳菲，等.數控機床可靠性技術的研究進展［J］.機械工程學報，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］楊兆軍，李小兵，許彬彬，等.加工中心時間動態可靠性建模［J］.機械工程學報，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］吳堅，趙陽，何睿.基于支持向量機回歸算法的電子機械制動傳感器系統故障診斷［J］.吉林大學學報:工學版，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］寇曉麗，劉三陽.基于模擬退火的粒子群算法求解約束優化問題［J］.吉林大學學報:工學版，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］王凌，劉波.微粒群優化與調度算法［M］.北京:清華大學出版社，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.