基于極大似然估計的數控系統故障分布模型

游達章，李秋實，肖 哲

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

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

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

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

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

根據課題組從2007年4月開始收集的某數控系統故障時間間隔實驗的數據:

值得說明的是，該故障數據為工廠直接收集到的現場數據，其數控系統2002年左右投入使用。

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

對密度函數取對數求和得似然函數。

對似然函數求關于參數β，α偏導數得似然方程:簡化上述似然方程，可得

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

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

2 Matlab和 Newton-Raphson迭代法的優化算法

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

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

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

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

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

似然函數:

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

圖1 F(m)的函數圖像

2.2 依賴于迭代初值的關于收斂的判定條件

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

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

故選取的初值如果在收斂領域內則牛頓迭代可快速收斂，若初始值偏離收斂域較遠則收斂速度很慢，甚至不發生收斂。

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

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

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

2)F'(x)≠0

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

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

由上述4個條件可歸結為4種情形:

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

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

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

2.3 Matlab環境下的Newton-Raphson迭代循環

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

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

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

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

迭代初值m0為區間［1，2］內以0.01為步長的100個數，循環終止條件通過 Matlab程序下按圖3算法流程完成Newton-Raphson迭代，得到結果為1.3119。

圖3 循環流程圖

圖4 不同初始值的迭代圖

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

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

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

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

圖5 擬合曲線圖

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

4 結束語

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

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

［2］賈亞洲.數控系統可靠性國內外現狀及對策［J］.中國制造業信息化，2006(06):51-53.

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

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

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

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

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

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

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

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