游達章，李秋實，肖 哲

(湖北工業(yè)大學(xué)機械工程學(xué)院，湖北武漢430068)

在數(shù)控系統(tǒng)［1］和機床的可靠性研究上，許多學(xué)者做了大量的研究［2］，而故障時間分布、參數(shù)估計、壽命評估等是數(shù)控機床或系統(tǒng)可靠性的一個重要內(nèi)容。數(shù)控系統(tǒng)故障間隔時間服從威布爾分布［3］，極大似然估計方法下的三參數(shù)威布爾分布，估算結(jié)果高效，準(zhǔn)確，但計算量大，過程復(fù)雜［4］。威布爾分布參數(shù)估計有如下方法:常用的圖估計法直觀簡單，通俗易懂，使用方便，但依賴于個人對數(shù)據(jù)曲線的視覺審查，所得結(jié)果往往因人而異，精確度較差;最小二乘法是最早的參數(shù)估計法，因方法簡單而使用較廣泛，但精確度不及極大似然估計法;極大似然估計是一種十分有效且精確度很高的算法，但其求解效率卻一直受困于龐雜的數(shù)據(jù)和復(fù)雜的超越方程求解［5］。國內(nèi)外很多學(xué)者一直在探索一種高效的極大似然解法，文獻［7］即通過給定形狀參數(shù)或尺寸參數(shù)等方法來達到似然方程降階的目的。這些方法雖然一定程度上簡化了計算，但計算過程依舊繁雜。

Newton-Raphson迭代法較梯度法、二分法等簡單迭代法具有更快的二階收斂速度，尤其是在方程根附近具有較高的收斂速度，因此是將近似根精確化的一種相當(dāng)有效的方法，但其缺點在于過分依賴于初始值，一旦選取的迭代初值偏離解析解太遠，將得不到收斂結(jié)果。

極大似然估計是一種基于出現(xiàn)概率最大原理下的參數(shù)估計方法，其算法精度高，適應(yīng)性廣，在參數(shù)估計問題中占有重要地位，尤其在處理不完全壽命的情況下，極大似然估計具有明顯優(yōu)勢。針對極大似然估計中超越方程求解困難的問題，采用穩(wěn)定快速的Newton-Raphson迭代法建立極大似然參數(shù)估計求解模型。同時，對于Newton-Raphson迭代法高度依賴于其初始值選取的問題，則利用Matlab軟件首先繪制出似然函數(shù)的曲線，并根據(jù)曲線在零點附近的位置，大致確定初始值的選取區(qū)間。然后依據(jù)Newton-Raphson迭代收斂的充分條件，對選取區(qū)間內(nèi)的初值，判定其迭代結(jié)果是否收斂。最后由Matlab繪制出不同初始值迭代過程的三維趨勢圖，以此驗證Newton-Raphson迭代法的解題結(jié)果。

1 基于極大似然估計算法的威布爾分布模型

根據(jù)課題組從2007年4月開始收集的某數(shù)控系統(tǒng)故障時間間隔實驗的數(shù)據(jù):

值得說明的是，該故障數(shù)據(jù)為工廠直接收集到的現(xiàn)場數(shù)據(jù)，其數(shù)控系統(tǒng)2002年左右投入使用。

采用利用極大似然估計算法進行威布爾分布參數(shù)的點估計。把樣本 { ti}(i=1，2，…n) 代入

對密度函數(shù)取對數(shù)求和得似然函數(shù)。

對似然函數(shù)求關(guān)于參數(shù)β，α偏導(dǎo)數(shù)得似然方程:簡化上述似然方程，可得

式(3)不含參數(shù)α，只含有參數(shù)β，因此，先由似然方程(3)求出參數(shù)β∧，再把β∧

代入似然方程()，即可求得另一個參數(shù)α∧，最后求得兩個參數(shù)的極大似然估計值。

2 Matlab和 Newton-Raphson迭代法的優(yōu)化算法

威布爾分布參數(shù)的極大似然估計在求解未知參數(shù)解中，似然方程式(3)是超越方程，用初等方法無法求解，此處借助計算機用Newton-Raphson迭代法，在迭代過程中初始值可以用圖估法或最小二乘估計給出。

Newton-Raphson迭代法雖然具有二階收斂速度，但其前提是選定的初值接近似然方程的解析解，否則有可能得不到收斂結(jié)果，因此迭代初值的確定直接影響Newton-Raphson迭代法求解的結(jié)果，甚至由于發(fā)散而得不到解。

首先利用Matlab繪出似然函數(shù)F(m)的圖像，初步確定迭代初值的區(qū)域，然后利用Matlab繪出F(m)的一階導(dǎo)函數(shù)和二階導(dǎo)函數(shù)圖像，結(jié)合Newton-Raphson迭代法收斂的充分條件最終確定其迭代初值的區(qū)間。

2.1 采用圖估法確定迭代初值的取值區(qū)間

首先根據(jù)似然方程(3)，建立似然函數(shù)F(m)。為了方便說明和計算，將參數(shù)β由字母m代替。據(jù)文獻［7］，一般情況下，對于機械系統(tǒng)可靠性工程威布爾分布的形狀參數(shù)的范圍在1～10之間，考慮到數(shù)控系統(tǒng)會發(fā)生初期故障，取m的范圍(0，10］。

似然函數(shù):

通過Matlab，以m為橫坐標(biāo)，F(xiàn)(m)為縱坐標(biāo)，繪制出m在0～10區(qū)間對應(yīng)F(m)的函數(shù)曲線(圖1)。由圖1可以確定，當(dāng)F(m)=0時m的值位于1.5附近，則大致取Newton-Raphson迭代過程的初始值區(qū)間為［1，2］。

圖1 F(m)的函數(shù)圖像

2.2 依賴于迭代初值的關(guān)于收斂的判定條件

Newton-Raphson迭代法局部收斂性定理:

設(shè)x*是方程f(x)=0的根，在包含x*的某個開區(qū)間內(nèi)f″(x)連續(xù)，且f'(x)≠0，則存在x*的鄰域Bδ(x*)= ［x*－ δ，x*+ δ］，使得任取初值 x0∈Bδ(x*)，由牛頓迭代法產(chǎn)生的迭代序列{xn}，以不低于二階的收斂速度收斂于x*。

故選取的初值如果在收斂領(lǐng)域內(nèi)則牛頓迭代可快速收斂，若初始值偏離收斂域較遠則收斂速度很慢，甚至不發(fā)生收斂。

Newton-Raphson迭代法收斂的充分條件:

設(shè) F(x) ∈［a，b］，且在區(qū)間［a，b］上滿足:

1)F(a)F(b)＜0

2)F'(x)≠0

3)F ″(x) 在區(qū)間［a，b］上恒正或恒負

4)x0∈［a，b］，滿足條件 F(x0)F ″(x0) ＞ 0

由上述4個條件可歸結(jié)為4種情形:

若在區(qū)間［a，b］上F(x)及其導(dǎo)函數(shù)的特性滿足以上任意一種情形，則Newton-Raphson迭代序列{xn}，單調(diào)的收斂于方程F(x)=0的唯一解。

圖2 F'(m)和F″(m)的函數(shù)圖像

為求證在區(qū)間［1，2］上似然函數(shù)F(m)的收斂性，利用Matlab在區(qū)間［1，2］上繪制出F'(m)和F″(m)的函數(shù)圖像，結(jié)合圖1和圖2得:F(1) ＞0，F(xiàn)(2) ＜0，F(xiàn)'(m) ＜0，F(xiàn) ″(m) ＞0，與上述情形3)相符，故判定Newton-Raphson迭代初值在區(qū)間［1，2］內(nèi)是收斂的，在該區(qū)域內(nèi)選取初值可行。

2.3 Matlab環(huán)境下的Newton-Raphson迭代循環(huán)

Newton-Raphson迭代法的基本思想:

設(shè)xk是f(x)=0的一個近似根，把f(x)在xk處作泰勒展開

取前兩項來近似代替f(x)，為f(x)的線性化，則得到近似的線性方程

設(shè)f'(xk)≠0，令其解為xk+1，得到Newton-Raphson迭代公式:

迭代初值m0為區(qū)間［1，2］內(nèi)以0.01為步長的100個數(shù)，循環(huán)終止條件通過 Matlab程序下按圖3算法流程完成Newton-Raphson迭代，得到結(jié)果為1.3119。

圖3 循環(huán)流程圖

圖4 不同初始值的迭代圖

為了驗證Newton-Raphson迭代法的收斂性、有效性，采用Matlab繪出了迭代次數(shù)、迭代初值和迭代值為橫、縱坐標(biāo)的three-dimensions圖。圖中每塊區(qū)域由多條曲線組成，迭代初值范圍為［1，2］，每條曲線對應(yīng)一個迭代初值，曲線走向代表收斂趨勢，曲線上的節(jié)點數(shù)代表迭代次數(shù)，同一顏色的曲線表示經(jīng)過相同的迭代次數(shù)后收斂，可以發(fā)現(xiàn)迭代次數(shù)多則4次，少則2次即收斂，與Newton-Raphson迭代具有二階收斂速度的特性相符合。同時發(fā)現(xiàn)越接近解析解的初始值，其迭代次數(shù)越少，也就表示收斂速度越快。只要能夠準(zhǔn)確定位初始值的選取區(qū)間，不論在區(qū)間內(nèi)選取何初始值，最終經(jīng)過若干次迭代都收斂于同一結(jié)果1.3119。

3 該模型與最小二乘法算法模型的比較

對實驗數(shù)據(jù)采用文獻［3］中的最小二乘估計算法得到形狀參數(shù) β為 1.0226，尺度參數(shù) α 為429.956，則由最小二乘算法得到的故障間隔時間的概率密度函數(shù)f(t)、分布函數(shù)F(t)分別為:

采用文獻［3］中故障間隔時間經(jīng)驗分布函數(shù)圖的繪制方法，在Matlab中繪制經(jīng)驗分布函數(shù)曲線、極大似然算法分布函數(shù)曲線、最小二乘算法分布函數(shù)曲線。

圖5 擬合曲線圖

對比兩種算法得到的分布函數(shù)曲線，基于極大似然算法擬合的分布函數(shù)曲線更接近于經(jīng)驗分布函數(shù)圖，擬合度更高，體現(xiàn)出了極大似然估計的有效性和高精度性，基于極大似然參數(shù)估計的模型和Newton-Raphson迭代求解法，更加精確有效。

4 結(jié)束語

采用了雙參數(shù)威布爾分布描述數(shù)控系統(tǒng)的故障時間分布，并介紹了基于極大似然法的參數(shù)估計。針對雙參數(shù)威布爾分布極大似然估計的數(shù)值求解算法困難的問題，應(yīng)用了穩(wěn)定快速的Newton-Raphson迭代法，同時利用Matlab圖形，很好解決了Newton-Raphson迭代初值的選取及數(shù)值解的收斂判斷問題。極大似然估計法具有有效性、高精度性，而Newton-Raphson迭代法穩(wěn)定并二階收斂。該模型結(jié)合二者特點，理論上科學(xué)可行，解題方法快速、穩(wěn)定、精確，適用于可靠性工程和壽命分析的威布爾分布極大似然參數(shù)估計。

［1］盛伯浩.我國數(shù)控機床現(xiàn)狀與技術(shù)發(fā)展策略［J］.現(xiàn)代金屬加工，2005(06):38-24.

［2］賈亞洲.數(shù)控系統(tǒng)可靠性國內(nèi)外現(xiàn)狀及對策［J］.中國制造業(yè)信息化，2006(06):51-53.

［3］ 游達章.數(shù)控系統(tǒng)加速壽命試驗方法及可靠性評估技術(shù)研究［D］.武漢:華中科技大學(xué)，2011.

［4］ 范 英，王順坤，晉民杰.多種數(shù)據(jù)狀態(tài)下三參數(shù)Weibull分布的極大似然估計［J］.機械強度，2012，34(01):53-57.

［5］趙呈建.威布爾分布參數(shù)估計的研究［D］.天津:南開大學(xué)，2007.

［6］ 徐 薇，胡偉明，孫 鵬，等.基于兩參數(shù)威布爾分布的設(shè)備可靠性預(yù)測研究［J］.中國工程機械學(xué)報，2013，11(02):112-116.

［7］ 楊存謀，聶 宏.三參數(shù)Weibull分布參數(shù)的極大似然估計數(shù)值解法［J］.南京航空航天大學(xué)學(xué)報，2007，29(01):23-25.

［8］ 肖光強，方 壯，余顯志.對牛頓迭代法條件的一個改進［J］.湖北民族學(xué)院學(xué)報，2008，26(04):395-397.

［9］ 李麗容.對牛頓迭代法的改進［J］.中國水運，2006，4(05):204-206.

［10］李慶楊，王能超，易大義.數(shù)值分析［M］.武漢:華中科技大學(xué)出版社，2001.