基于極大似然估計(jì)的數(shù)控系統(tǒng)故障分布模型

游達(dá)章，李秋實(shí)，肖 哲

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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)。為了方便說(shuō)明和計(jì)算，將參數(shù)β由字母m代替。據(jù)文獻(xiàn)［7］，一般情況下，對(duì)于機(jī)械系統(tǒng)可靠性工程威布爾分布的形狀參數(shù)的范圍在1～10之間，考慮到數(shù)控系統(tǒng)會(huì)發(fā)生初期故障，取m的范圍(0，10］。

似然函數(shù):

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

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

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

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

設(shè)x*是方程f(x)=0的根，在包含x*的某個(gè)開(kāi)區(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)則牛頓迭代可快速收斂，若初始值偏離收斂域較遠(yuǎn)則收斂速度很慢，甚至不發(fā)生收斂。

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

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

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

2)F'(x)≠0

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

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

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

若在區(qū)間［a，b］上F(x)及其導(dǎo)函數(shù)的特性滿(mǎn)足以上任意一種情形，則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的一個(gè)近似根，把f(x)在xk處作泰勒展開(kāi)

取前兩項(xiàng)來(lái)近似代替f(x)，為f(x)的線(xiàn)性化，則得到近似的線(xiàn)性方程

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

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

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

圖4 不同初始值的迭代圖

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

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

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

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

圖5 擬合曲線(xiàn)圖

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

4 結(jié)束語(yǔ)

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

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

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

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

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

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

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

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

［8］ 肖光強(qiáng)，方 壯，余顯志.對(duì)牛頓迭代法條件的一個(gè)改進(jìn)［J］.湖北民族學(xué)院學(xué)報(bào)，2008，26(04):395-397.

［9］ 李麗容.對(duì)牛頓迭代法的改進(jìn)［J］.中國(guó)水運(yùn)，2006，4(05):204-206.

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