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

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

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

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

1 試驗(yàn)數(shù)據(jù)

對(duì)23 臺(tái)加工中心進(jìn)行了現(xiàn)場(chǎng)試驗(yàn)，根據(jù)所得數(shù)據(jù)整理出加工中心整機(jī)的105個(gè)故障間隔時(shí)間數(shù)據(jù)，具體見表1.

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

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

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

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

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

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

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

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

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

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

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

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

由上式得到a 的表達(dá)式如下:

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

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

根據(jù)上式可得

因?yàn)榻?jīng)三參數(shù)威布爾分布線性變換后yi=ln ln所以

又因?yàn)閎=β，所以此時(shí)β 的估計(jì)值可以表示成僅含有未知參數(shù)γ 的表達(dá)式β(γ).

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

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

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

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

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

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

3 差分進(jìn)化算法

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

3.1 DE 算法的基本步驟

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

5 結(jié)語

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

［1］申桂香，孟書，張英芝，等. 平均秩次法在子系統(tǒng)可靠性建模中的應(yīng)用［J］.吉林大學(xué)學(xué)報(bào):工學(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ì)方法的對(duì)比［J］.華南理工大學(xué)學(xué)報(bào):自然科學(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í)截尾樣本下三參數(shù)威布爾分布的矩估計(jì)［J］. 數(shù)學(xué)的實(shí)踐與認(rèn)識(shí)，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é)機(jī)械科學(xué)與工程學(xué)院，2012.

［7］李雄軍. 對(duì)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］吳亮紅.差分進(jìn)化算法及應(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ìn)差分進(jìn)化算法的估計(jì)等值法［J］. 華南理工大學(xué)學(xué)報(bào):自然科學(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é)報(bào):自然科學(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.