基于ITLS和DE的加工中心三參數(shù)威布爾分布*

王曉峰 張英芝 申桂香 龍哲 張立敏

(1.吉林大學(xué) 建設(shè)工程學(xué)院，吉林 長春130022;2.吉林大學(xué) 機械科學(xué)與工程學(xué)院，吉林 長春130022)

近些年來，威布爾分布模型已經(jīng)廣泛應(yīng)用于加工中心可靠性研究. 在進行加工中心的可靠性建模時，常假設(shè)數(shù)據(jù)服從位置參數(shù)γ=0 的二參數(shù)威布爾分布［1-2］.但加工中心并非總在t=0 時刻發(fā)生故障，并且二參數(shù)威布爾分布的處理過程中人為地丟掉了一個參數(shù)，增大了其參數(shù)估計結(jié)果的誤差. 因此，選用位置參數(shù)γ≠0 的三參數(shù)威布爾分布更有利于反映加工中心可靠性的真實情況. 三參數(shù)威布爾分布的參數(shù)估計比較復(fù)雜，有關(guān)學(xué)者進行了許多相關(guān)研究［3-5］.在不斷探討其在加工中心可靠性建模應(yīng)用的過程中發(fā)現(xiàn)以下問題［6］:①對加工中心三參數(shù)威布爾分布模型的參數(shù)優(yōu)化多采用迭代法，即通過一定步長不斷增大優(yōu)化參數(shù)值，計算每次迭代得到的目標(biāo)函數(shù)值，最后取最優(yōu)目標(biāo)函數(shù)所對應(yīng)的參數(shù)值為估計值.這種迭代法如果步長過大可能會丟失最優(yōu)解，步長過小則以犧牲時間為代價，有時迭代法可能會耗時很長.②在對三參數(shù)威布爾分布進行參數(shù)估計時多采用最小二乘法確定回歸直線方程. 由于選擇x 方向或者y 方向擬合所得的回歸直線是不同的［7］，所以當(dāng)在y 方向上擬合直線時，只考慮到y(tǒng) 坐標(biāo)上的誤差而沒有考慮到x 坐標(biāo)上的誤差. 也就是說只考慮了故障發(fā)生概率F 的誤差，沒有考慮間隔時間t 的誤差. 但是故障間隔時間t 的誤差必然存在，即x 坐標(biāo)誤差必然存在，所以x 坐標(biāo)的誤差也應(yīng)該予以充分考慮.③在三參數(shù)威布爾線性變換后，并非是尺度參數(shù)α 直接使得目標(biāo)函數(shù)達到最佳，而且目標(biāo)函數(shù)所考慮的誤差是y 坐標(biāo)上的誤差，并沒有直接考慮到故障發(fā)生概率F 的誤差，但實際上兩組誤差并不相同.

針對上述問題，文中采用一種新興的優(yōu)化算法(即差分進化算法)進行參數(shù)優(yōu)化，它可以準(zhǔn)確快速地尋得目標(biāo)函數(shù)全局最優(yōu)解，克服了迭代法的缺點.其次，結(jié)合三參數(shù)威布爾分布特點，文中提出一種改進的整體最小二乘法. 該方法一方面具有整體最小二乘法的固有特點，即同時兼顧x 坐標(biāo)和y 坐標(biāo)誤差，彌補了最小二乘法的不足，使所得結(jié)果更加真實合理;另一方面以故障發(fā)生概率F 的偏差平方和作為目標(biāo)函數(shù)進行優(yōu)化，確保最終所得的三參數(shù)威布爾曲線與原始數(shù)據(jù)點最大程度地接近.

1 試驗數(shù)據(jù)

對23 臺加工中心進行了現(xiàn)場試驗，根據(jù)所得數(shù)據(jù)整理出加工中心整機的105個故障間隔時間數(shù)據(jù)，具體見表1.

表1 加工中心故障間隔時間1)Table 1 Time between failures of machining center h

2 基于改進的整體最小二乘法的參數(shù)估計

三參數(shù)威布爾分布函數(shù)為

式中:α 為尺度參數(shù)，α ＞0;β 為形狀參數(shù)，β ＞0;γ 為位置參數(shù)，γ≥0.

對式(1)兩次取對數(shù)后線性變換為

設(shè)一元線性回歸模型為

最小二乘(LS)法的約束準(zhǔn)則為

整體最小二乘(TLS)法的約束準(zhǔn)則為

由此可見:①最小二乘法的實質(zhì)是求一條直線，使所有已知點到這條直線的縱坐標(biāo)的距離平方和最小，因此，LS 法只顧及了y 坐標(biāo)的誤差沒有涉及x坐標(biāo)的誤差.②整體最小二乘法的實質(zhì)是求一條直線，使所有的已知點到此直線的距離和的平方最小，所以TLS 法同時顧及了x 坐標(biāo)和y 坐標(biāo)的誤差. 因此將整體最小二乘法代替最小二乘法引入加工中心可靠性建模后，既可以考慮到故障發(fā)生概率F 的誤差，又可以考慮到故障間隔時間t 的誤差，使所得三參數(shù)威布爾模型更加真實合理.

2.1 對β 的參數(shù)估計

整體最小二乘法最早由Golub 等［8］在1980年提出.Byungsoo［9］總結(jié)了TLS 法的重要理論結(jié)果和計算方法，發(fā)現(xiàn)在一些典型應(yīng)用中用TLS 法相較于用LS 法可以使參數(shù)估計精度顯著提高.本節(jié)主要通過TLS 法對威布爾分布中的形狀參數(shù)β 進行估計.

根據(jù)式(5)取QTLS關(guān)于a 的偏導(dǎo)并令其為0:

由上式得到a 的表達式如下:

把式(6)代入式(5)中，得到

根據(jù)式(7)求QTLS關(guān)于b 的偏導(dǎo)并令其為0:

根據(jù)上式可得

因為經(jīng)三參數(shù)威布爾分布線性變換后yi=ln ln所以

又因為b=β，所以此時β 的估計值可以表示成僅含有未知參數(shù)γ 的表達式β(γ).

2.2 對α 和γ 的參數(shù)估計

根據(jù)三參數(shù)威布爾線性變換過程中的特點，可以發(fā)現(xiàn):首先，式(2)中a= -βlnα，b=β，其中β 是可以直接求的，但是求出的α 并非直接使QTLS達到最小，而是a 使QTLS達到最小;其次，以上方法考慮的是x和y 坐標(biāo)上的誤差，沒有直接考慮到故障間隔時間t和故障發(fā)生概率F 的誤差，但兩組誤差并不相同.

針對以上問題，在上節(jié)得到β 的估計值β(γ)的基礎(chǔ)上，對α 和γ 的估計過程如下.

將β(γ)代入原三參數(shù)威布爾分布函數(shù)(式(1))中，此時分布函數(shù)F^中只含有未知參數(shù)α 與γ. 然后根據(jù)以下約束準(zhǔn)則進行優(yōu)化:

QITLS實際上是威布爾分布函數(shù)F 在節(jié)點ti處的偏差平方和.代入β 表達式β(γ)后QITLS為自變量α與γ 的函數(shù)，結(jié)合表1 中數(shù)據(jù)，應(yīng)用差分進化算法對α 與γ 進行優(yōu)化使QITLS達到最小值.

根據(jù)以上分析過程可知，改進的整體最小二乘法繼承了整體最小二乘法的特點，即兼顧自變量和因變量的誤差.根據(jù)三參數(shù)威布爾分布函數(shù)線性變換的特點，為給出參數(shù)α 的更佳估計值，改進的整體最小二乘法以F 的偏差平方和最小作為約束準(zhǔn)則進行優(yōu)化，考慮到了初始數(shù)據(jù)t 和F 的誤差，使所得曲線在擬合效果上達到最佳.

3 差分進化算法

差分進化(DE)算法最早由Storn 等［10］在1995年提出.DE 可以動態(tài)跟蹤當(dāng)前的搜索情況以調(diào)整搜索策略，具有較好的全局收斂能力和魯棒性，適用于一些常規(guī)的數(shù)學(xué)規(guī)劃方法所無法求解的復(fù)雜環(huán)境中的優(yōu)化問題.

3.1 DE 算法的基本步驟

在使用DE 算法之前需要設(shè)定以下幾個參數(shù):種群大小N，即群體中個體的數(shù)量;最大迭代代數(shù)G;變異因子M，通常取值范圍為M∈［0，2］;交叉概率C，通常取值范圍為C∈［0，1］. 實施DE 算法時包含以下4個基本步驟［11］:

(1)初始化.個體初始化公式為

式中，r1，r2，r3∈{1，2，…，N}，為互不相同的整數(shù)且r1，r2，r3與i 不同，因此種群規(guī)模N≥4.

rand(j)∈［0，1］，為均勻分布隨機數(shù);rand(i)∈［1，2，…，D］，是隨機選擇的維數(shù)變量索引.

(4)選擇操作.以最小化優(yōu)化為例，選擇操作為

3.2 DE 對三參數(shù)威布爾模型的優(yōu)化

根據(jù)2.2 節(jié)分析結(jié)果，結(jié)合表1 中數(shù)據(jù)，利用DE對參數(shù)α 與γ 進行優(yōu)化，使式(9)達到最小值.文中的DE 優(yōu)化過程是通過在Matlab 中編程實現(xiàn)的.其優(yōu)化α 與γ 的基本流程圖如圖1 所示，具體步驟如下:

圖1 DE 優(yōu)化流程圖Fig.1 Process of DE optimization

(1)初始化種群. 根據(jù)文獻［12-15］研究結(jié)論，文中設(shè)置種群規(guī)模N =10;最大迭代代數(shù)G =150;變異因子M=0.5;交叉概率C=0.9;count=0;因為有α 與γ 兩個參數(shù)，所以每個個體有兩個分量，即種群維數(shù)D=2;γ 的取值應(yīng)該為0≤γ≤tmin，根據(jù)該加工中心故障數(shù)據(jù)(見表1)，γ 的界限為γL=0，γU=18.29;αL=0，αU=500. 根據(jù)式(10)對α 與γ 進行個體初始化.

(2)變異、交叉.根據(jù)式(11)、(12)對種群中的每個個體進行變異和交叉操作.

(3)選擇.選取式(9)為目標(biāo)函數(shù)，根據(jù)式(13)進行下一子代的選擇操作.

(4)判斷. 如果達到規(guī)定的進化代數(shù)則停止程序輸出結(jié)果，否則返回變異操作步驟繼續(xù)運行.

根據(jù)以上流程和表1 中數(shù)據(jù)，計算得到該系列加工中心可靠性三參數(shù)威布爾分布模型在γ=13.7816、α=406.6531時QITLS達到最小，最小值為0.0402，根據(jù)式(8)有β=0.9720.故得到該加工中心可靠性三參數(shù)威布爾分布函數(shù)為

4 參數(shù)估計結(jié)果對比

根據(jù)普通最小二乘法和文中改進的整體最小二乘法分別得到兩個不同的三參數(shù)威布爾分布模型，對這兩個模型的相關(guān)參數(shù)進行對比分析，結(jié)果如表2所示.

表2 所得模型相關(guān)參數(shù)對比Table 2 Comparison of parameters of models

從表2 可以看出:對于尺度參數(shù)α，改進的整體最小二乘法得到的α 比最小二乘法得到的α 有所增大;對于形狀參數(shù)β，由改進的整體最小二乘法得到的β 比最小二乘法有所增大，β 值更加接近1，說明該型號加工中心處于偶然故障期;對于位置參數(shù)γ，由改進的整體最小二乘法得到的γ 小于由最小二乘法得到的γ 值.

對于QITLS，改進的整體最小二乘法所得QITLS小于最小二乘法所得QITLS，這說明由改進的整體最小二乘法所得三參數(shù)威布爾分布函數(shù)F 的偏差平方和最小，所得分布函數(shù)F 的曲線擬合優(yōu)度最佳.

5 結(jié)語

針對以往加工中心可靠性建模中存在的一些問題，提出了一種對三參數(shù)威布爾分布模型進行建模的新方法，并且結(jié)合某型號加工中心故障數(shù)據(jù)進行計算.最后通過對比分析可以發(fā)現(xiàn)文中所提出的方法具有以下特點:所采用的DE 算法具有快速收斂、全局收斂、指導(dǎo)記憶能力、穩(wěn)定性、通用性等特點，較普通的迭代算法更適合解決這類復(fù)雜目標(biāo)函數(shù)的尋優(yōu)問題;改進的整體最小二乘法兼顧到了故障間隔時間t 和故障發(fā)生概率F 的誤差，使分析結(jié)果更加合理.改進的整體最小二乘法以三參數(shù)威布爾分布函數(shù)F 在節(jié)點ti處的偏差平方和作為約束準(zhǔn)則進行優(yōu)化，所得曲線在整體上最大程度地接近原始數(shù)據(jù)點，對三參數(shù)威布爾分布曲線擬合效果最好. 另外，應(yīng)用文中所提出的方法對該型號加工中心的三參數(shù)威布爾分布進行了參數(shù)估計，根據(jù)所得結(jié)果，可以發(fā)現(xiàn)形狀參數(shù)β 近似等于1，說明該型號加工中心處于偶然故障期.

［1］申桂香，孟書，張英芝，等. 平均秩次法在子系統(tǒng)可靠性建模中的應(yīng)用［J］.吉林大學(xué)學(xué)報:工學(xué)版，2014，44(1):101-105.Shen Gui-xiang，Meng Shu，Zhang Ying-zhi，et al. Application of average rank time method in reliability modeling for subsystems［J］. Journal of Jilin University:Engineering and Technology Edition，2014，44(1):101-105.

［2］王曉峰，申桂香，張英芝，等. 可靠性模型參數(shù)估計方法的對比［J］.華南理工大學(xué)學(xué)報:自然科學(xué)版，2011，39(6):47-52.Wang Xiao-feng，Shen Gui-xiang，Zhang Ying-zhi，et al.Comparison of parameter estimation methods for reliability model［J］.Journal of South China University of Technology:Natural Science Edition，2011，39(6):47-52.

［3］Hideki Nagatsuka，Toshinari Kamakura.A consistent method of estimation for the three-parameter Weibull distribution［J］. Computational Statistics and Data Analysis，2013，58:210-226.

［4］Moeini Asghar，Jenab Kouroush.Fitting the three parameter Weibull distribution with ross entropy ［J］. Applied Mathematical Modelling，2013，37:6354-6363.

［5］孫麗玢，湯銀才. 在定時截尾樣本下三參數(shù)威布爾分布的矩估計［J］. 數(shù)學(xué)的實踐與認識，2010，40(22):156-161.Sun Li-bin，Tang Yin-cai.Moment estimation of the threeparameter Weibull distribution under type I censored samples［J］. Mathematics in Practice and Theory，2010，40(22):156-161.

［6］王曉峰.加工中心可信性影響度分析及增長技術(shù)研究［D］.長春:吉林大學(xué)機械科學(xué)與工程學(xué)院，2012.

［7］李雄軍. 對X 和Y 方向最小二乘線性回歸的討論［J］.誤差與數(shù)據(jù)處理，2005(1):50-52.Li Xiong-jun.The least squares discussion of X and Y direction ［J］. Modern Manufacturing Engineering，2005(1):50-52.

［8］Golub G H，Van Loan C F.An analysis of the total least squares problem［J］.SIAM Journal on Numerical Analysis，1980，17(6):883-893.

［9］Byungsoo Kim.Total least square method applied to rating curves［J］.Hydrol Process，2014，28(13):4057-4066.

［10］Storn R，Price K. Differential evolution—a simple and efficient adaptive scheme for global optimization over continuous spaces ［R］. Berkley:Technical Report International Computer Science Institute，1995.

［11］Neri Ferrante，Tirronen Ville. Recent advances in differential evolution:a survey and experimental analysis［J］.Artif Intell Rev，2010，33(1):61-106.

［12］吳亮紅.差分進化算法及應(yīng)用研究［D］. 長沙:湖南大學(xué)電氣與信息工程學(xué)院，2007.

［13］Qin A K，Huang V L，Suganthan P N.Differential evolution algorithm with strategy adaptation for global numerical optimization［J］.IEEE Trans Evol Comput，2009，13(2):398-417.

［14］張寶珍，張堯，林凌雪. 基于改進差分進化算法的估計等值法［J］. 華南理工大學(xué)學(xué)報:自然科學(xué)版，2014，42(4):7-12.Zhang Bao-zhen，Zhang Yao，Lin Ling-xue. Estimation equivalence method based on modified differential evolution algorithm［J］. Journal of South China University of Technology:Natural Science Edition，2014，42(4):7-12.

［15］魏德敏，陳貴濤.自適應(yīng)子群體米母算法及其在混凝土框架位移性能優(yōu)化中的應(yīng)用［J］.華南理工大學(xué)學(xué)報:自然科學(xué)版，2013，41(10):108-116.Wei De-min，Chen Gui-tao. Self-daptive subpopulationsbased memetic algorithm and its application to displacement performance optimization of concrete frames［J］.Journal of South China University of Technology:Natural Science Edition，2013，41(10):108-116.