多臺(tái)數(shù)控機(jī)床的時(shí)間截尾可靠性評(píng)估

王智明，楊建國(guó)，王國(guó)強(qiáng)，張根保

(1.上海交通大學(xué)機(jī)械與動(dòng)力工程學(xué)院，200240上海，wangzm301@sjtu.edu.cn; 2.重慶大學(xué)機(jī)械工程學(xué)院先進(jìn)制造技術(shù)實(shí)驗(yàn)室，400030重慶)

在數(shù)控機(jī)床的可靠性評(píng)估中，平均故障間隔時(shí)間(tMTBF)、給定時(shí)間的可靠度R(t)和一定可靠度下的工作時(shí)間t(R)是評(píng)估機(jī)床可靠性的重要指標(biāo)之一.可靠性試驗(yàn)是可靠性評(píng)估的基礎(chǔ)，一般常采用多臺(tái)定時(shí)截尾試驗(yàn)［1-2］.

分布模型的選擇是評(píng)估的難點(diǎn)和關(guān)鍵.常見的分布模型有指數(shù)分布和Weibull分布，其選擇的方法是先假設(shè)故障間隔時(shí)間(tTBF)服從某一分布，隨后用K-S或χ2檢驗(yàn)進(jìn)行檢驗(yàn)［3］，但有時(shí)會(huì)出現(xiàn)2種分布均被通過的情況［4］;對(duì)于可線性化的分布模型，可用最小二乘法求解模型，并用相關(guān)系數(shù)進(jìn)行檢驗(yàn)［5］，但該方法所得模型并非最優(yōu);另外，還有其他選擇模型的方法，如:綜合法［3］、誤差面積比［4］、灰關(guān)聯(lián)［6］、神經(jīng)網(wǎng)絡(luò)［7］等方法.但對(duì)于少樣本數(shù)據(jù)，上述方法所得模型精度一般較低.在機(jī)床可靠性評(píng)估中，參數(shù)估計(jì)多局限于參數(shù)的點(diǎn)估計(jì)，對(duì)于參數(shù)區(qū)間估計(jì)的報(bào)道則較少.

本文基于AIC、BIC信息準(zhǔn)則優(yōu)選tTBF模型和參數(shù)，給出了tTBF模型參數(shù)和可靠性指標(biāo)的點(diǎn)估計(jì)，用 Fisher信息矩陣法給出了模型參數(shù)及R(tMTBF)、t(0.9)等可靠性指標(biāo)的區(qū)間估計(jì).最后用蒙特卡洛法(MC)提高了區(qū)間估計(jì)的精度.

1 基于AIC、BIC信息準(zhǔn)則的模型和方法選擇

1.1 模型選擇

數(shù)控機(jī)床時(shí)間截尾數(shù)據(jù)通常為不完全數(shù)據(jù)，其中包含右刪失數(shù)據(jù).設(shè)某次試驗(yàn)投入k臺(tái)機(jī)床，所得數(shù)據(jù)為t=(ti，tcj)(i=1，2，…，n;j=1，2，…，k)，其中ti為失效數(shù)據(jù)，tcj為截尾數(shù)據(jù).一般每臺(tái)機(jī)床有一個(gè)截尾數(shù)據(jù)，且不相等.設(shè)截尾數(shù)據(jù)的概率密度和可靠度分別為f(t;θ)、R(t;θ)，則用極大似然法(MLE)優(yōu)選模型時(shí)，其對(duì)數(shù)似然函數(shù)為

令

則可得參數(shù)的估計(jì)值^θ.

不同模型對(duì)試驗(yàn)數(shù)據(jù)信息反映的準(zhǔn)確性各異，其中，只有對(duì)數(shù)據(jù)擬合度好、參數(shù)少、形式簡(jiǎn)單的模型才是最佳模型.AIC信息準(zhǔn)則將K-L距離和極大似然方法相結(jié)合，并利用了似然估計(jì)性質(zhì)［8-9］，其形式如下:

AIC=-2(maximum log likelihood)+2m. BIC信息準(zhǔn)則［10-11］是另一種模型選擇的有效方法之一，其形式如下:

BIC=-2(maximum log likelihood)+m×ln(n).其中m是參數(shù)數(shù)量，n是觀測(cè)數(shù)據(jù)數(shù).可見，BIC對(duì)增加的參數(shù)加上更大的懲罰項(xiàng)m×ln(n).因此，BIC更有可能選擇參數(shù)少的模型.AIC、BIC作為模型選擇的有效方法之一，其值越小者模型越好.

故障間隔時(shí)間分布模型一般為指數(shù)分布或Weibull分布，因此，以指數(shù)分布、Weibull分布或其他分布作為備選分布模型，用MLE法分別求取最大似然函數(shù)值，最后得到AIC、BIC值，該值最小者即為所求最佳模型.

1.2 方法選擇

常見的估計(jì)參數(shù)的方法有最小二乘法(LSM)和MLE法.其中，LSM根據(jù)因變量的不同可分X方向秩回歸(RRX)和Y方向秩回歸(RRY)［12］.

Weibull模型可經(jīng)對(duì)數(shù)變換線性化，然后用LSM求解.對(duì)于不完整數(shù)據(jù)可用Johnson法得到不完整數(shù)據(jù)中失效數(shù)據(jù)的失效順序號(hào)，再按中位秩估計(jì)得到相應(yīng)的可靠度［13］.Weibull模型用MLE法估計(jì)參數(shù)時(shí)，所得極大似然估計(jì)方程沒有封閉形式，可用迭代法解之.為了使其快速收斂，初值可取LSM所得估計(jì)值.最后用AIC、BIC準(zhǔn)則優(yōu)選最佳參數(shù)估計(jì)方法.

2 參數(shù)及可靠性指標(biāo)的區(qū)間估計(jì)

2.1 參數(shù)的區(qū)間估計(jì)

通常在大樣本情況下，極大似然估計(jì)量是漸進(jìn)正態(tài)分布的.在小樣本情況下，參數(shù)的對(duì)數(shù)更接近正態(tài)分布.所以，根據(jù)漸近正態(tài)原理有

進(jìn)而可得參數(shù)θ的估計(jì)區(qū)間值如下:

其中，var(^θ)為估計(jì)參數(shù)的方差，zα/2為置信度為1-α的正態(tài)分布的分位數(shù).因此，若求得參數(shù)的方差，即可得其估計(jì)區(qū)間.參數(shù)的方差可用Fisher信息矩陣法(FM)解決［14-15］，即參數(shù)方差矩陣為Fisher信息矩陣的逆.為了便于說明，下面以兩參數(shù)Weibull分布為例.

設(shè)時(shí)間截尾Weibull分布的失效數(shù)為n，截尾數(shù)為k，則由式(1)知其對(duì)數(shù)似然函數(shù)為

其中，β為形狀參數(shù)，η為尺度參數(shù)，均大于零，所以

由Fisher信息矩陣法知

2.2 可靠性指標(biāo)的區(qū)間估計(jì)

在獲得模型的參數(shù)估計(jì)值后，易得其點(diǎn)估計(jì)，平均故障間隔時(shí)間(tMTBF)、給定時(shí)間的可靠度R(t)和一定可靠度下的工作時(shí)間t(R)計(jì)算公式如下:

2.2.1 給定時(shí)間可靠度R(t)的區(qū)間估計(jì)

在Weibull分布中，設(shè)T=ln t，則T服從極小值分布，其可靠性函數(shù)為

μ、σ為該函數(shù)參數(shù).設(shè)μ=ln η，σ=1/β，u= β(ln t-ln η)，則極小值分布可靠性函數(shù)可變換為Weibull分布可靠性函數(shù)

所以，若u的置信區(qū)間已知，便可得出R(t)的置信區(qū)間.u的置信區(qū)間為

其中，

至此，由式(7)可得出給定時(shí)間可靠度的置信區(qū)間估計(jì)

2.2.2 一定可靠度下的工作時(shí)間t(R)的區(qū)間估計(jì)

在式(6)中，令

則

所以，v的置信區(qū)間為 CBv=［vL，vU］ =，最后得到t的置信區(qū)間為

3 案例分析

數(shù)據(jù)來源于文獻(xiàn)［1］，該文給出了10臺(tái)加工中心的不完整數(shù)據(jù)30個(gè)，其中失效數(shù)據(jù)20個(gè)，時(shí)間截尾數(shù)據(jù)10個(gè)，屬于小樣本數(shù)據(jù).具體為: 176.00，248.00，10.50，472.00，45.00，39.00，209.33， 261.25， 510.00，120.00， 224.00，267.50， 32.00， 50.00， 138.50， 398.00，478.00，353.00，348.00，137.06，332.50*，84.00*， 267.00*， 165.00*， 283.16*，700.00*，562.00*，387.00*，383.00*，1.50*(帶*號(hào)者為時(shí)間截尾數(shù)據(jù)).

3.1 模型選擇

以指數(shù)分布、Weibull分布作為備選分布模型，用MLE法求解，所得結(jié)果如表1所示.

表1 模型選擇結(jié)果

從表1可見，2參數(shù)Weibull分布的AIC、BIC值均較小，所以選取2參數(shù)Weibull分布作為時(shí)間截尾TBF可靠性評(píng)估模型.

3.2 方法選擇

2參數(shù)Weibull分布的可靠性函數(shù)如式(5)，線性化后令y=ln［-ln(R(t))］，x=ln(t)，則

R(t)可用如下中位秩公式估計(jì):

對(duì)于時(shí)間截尾數(shù)據(jù)，失效數(shù)據(jù)的失效順序ri可用Johnson法估計(jì)，公式如下:

得出ri后代替式(11)中的i即可得截尾數(shù)據(jù)的可靠度Ri.用LSM方法中的RRY估計(jì)參數(shù)時(shí)，由式(10)知β=1.024 4，η=407.942 6.為了選取最佳參數(shù)，還可用RRX和MLE法估計(jì)，結(jié)果如表2所示.由表2可見，當(dāng)2參數(shù)Weibull分布分別用RRX、RRY、MLE 3種不同方法估計(jì)參數(shù)時(shí)，MLE方法所得AIC、BIC最小，模型精度最高.

表2 方法選擇結(jié)果

至此，選擇MLE方法得模型參數(shù)、機(jī)床的tMTBF及該時(shí)刻的可靠度R(tMTBF)以及可靠度為0.9的工作時(shí)間t(0.9)的點(diǎn)估計(jì)分別如下:β= 1.200 9，η =376.989 3，tMTBF=354.545 9，R(tMTBF)=0.395 0，t(0.90)=57.88.

3.3 參數(shù)及可靠性指標(biāo)的區(qū)間估計(jì)

由式(3)、(4)知參數(shù)的協(xié)方差矩陣為

由式(2)得參數(shù)β、η的置信度為0.95的置信區(qū)間分別為［0.827 7，1.742 3］、［261.515 2，543.451 7］.另外，由式(8)得u的方差后知其估計(jì)區(qū)間為CBu=［-0.511 5，0.365 1］，其中方差var(u)=0.050 0.同樣，由式(9)得v的方差后可得其估計(jì)區(qū)間，最后得出R(tMTBF)和t(0.9)的估計(jì)區(qū)間分別為CBR(t)=［0.236 8，0.549 0］，CBt(0.9)=［26.91，124.47］.

由估計(jì)區(qū)間結(jié)果來看，該區(qū)間太大，估計(jì)精度較低，所以有必要提高區(qū)間估計(jì)的精度.事實(shí)上，上述區(qū)間估計(jì)采用漸進(jìn)正態(tài)分布原理，適合于大樣本.對(duì)于小樣本數(shù)據(jù)，其估計(jì)精度一般較低.通常將樣本單位數(shù)不少于30個(gè)的樣本稱為大樣本，不及30個(gè)的稱為小樣本.數(shù)控機(jī)床的故障數(shù)據(jù)一般為少樣本數(shù)據(jù)，因此，對(duì)于少樣本數(shù)據(jù)可用模擬抽樣的方法來擴(kuò)大樣本容量，從而提高區(qū)間估計(jì)精度.

3.4 MC仿真抽樣

2參數(shù)Weibull分布的tTBF的抽樣公式為

r為(0，1)區(qū)間上均勻分布的隨機(jī)變量.取樣本容量 N分別為 500、1 000、2 000、5 000、8 000、10 000，其參數(shù)和可靠性指標(biāo)抽樣區(qū)間估計(jì)分別如圖1、2所示.從圖可見，樣本容量為8 000和10 000時(shí)的估計(jì)結(jié)果已經(jīng)很接近，因此停止抽樣，并取容量為10 000時(shí)的抽樣結(jié)果為最終估計(jì)值.其結(jié)果如表3所示.從表3可以看出，估計(jì)區(qū)間均已大大減少.其中，t(0.9)、R(tMTBF)區(qū)間長(zhǎng)度分別減少93.21%、95.07%.

圖1 參數(shù)β、n的區(qū)間仿真

圖2 可靠性指標(biāo)R(tMTBF)、t(0.9)的區(qū)間仿真

表3 樣本容量為10 000時(shí)的區(qū)間模擬抽樣結(jié)果

4 討論

1)在Weibull分布模型中，當(dāng)形狀參數(shù)β=1時(shí)，Weibull分布就轉(zhuǎn)換為指數(shù)分布，模型處于浴盆曲線的偶然失效區(qū);當(dāng)β＜1或β＞1時(shí)，則模型對(duì)應(yīng)于浴盆曲線的早期失效或耗損失效.可見本案例中的加工中心處于耗損失效期.因此，Weibull分布模型因其靈活性作為故障間隔時(shí)間的分布模型而廣泛應(yīng)用于NC機(jī)床的可靠性分析中.在2-Weibull模型的參數(shù)估計(jì)中，MLE法較之LSE法有更高的估計(jì)精度.這是因?yàn)長(zhǎng)SE法在對(duì)數(shù)變形代換過程中改變了因變量的形態(tài)，變形后因變量殘差平方和最小并不能保證變形前因變量殘差平方和最小，同時(shí)，用Johnson法對(duì)不完整數(shù)據(jù)中的失效數(shù)據(jù)排序時(shí)，只考慮了其失效順序而沒有考慮截尾數(shù)據(jù)的相對(duì)大小.換言之，對(duì)于不同大小的截尾數(shù)據(jù)，當(dāng)用Johnson法分析時(shí)，失效數(shù)據(jù)可能有相同的失效順序.當(dāng)截尾數(shù)據(jù)較多時(shí)，只擬合重新排序后的失效數(shù)據(jù)(本案例只有20個(gè)失效數(shù)據(jù))的LSE法的評(píng)估精度自然會(huì)受到影響.MLE適用于完整數(shù)據(jù)和不完整數(shù)據(jù)，精度較之LSE法有較好的評(píng)估效果.但LSE法的評(píng)估結(jié)果可作為MLE法中迭代計(jì)算模型參數(shù)的初值，因此，使用MLE法時(shí)可結(jié)合LSE法而達(dá)到快速收斂的目的.本文所提方法雖然是針對(duì)時(shí)間截尾的，但仍適用于定數(shù)截尾實(shí)驗(yàn).當(dāng)采用定數(shù)截尾時(shí)，各有關(guān)公式中tcj=0.

5 結(jié)論

1)在NC機(jī)床的可靠性評(píng)估中，對(duì)不完整數(shù)據(jù)可結(jié)合LSE法的結(jié)果，優(yōu)先選用MLE法.

2)AIC、BIC信息準(zhǔn)則結(jié)合了MLE的計(jì)算結(jié)果和模型的簡(jiǎn)約原則，是機(jī)床可靠性評(píng)估中優(yōu)選模型和模型參數(shù)的有效方法.

3)評(píng)估模型參數(shù)和可靠性指標(biāo)的區(qū)間時(shí)，對(duì)于大樣本數(shù)據(jù)，可采用基于漸進(jìn)正態(tài)分布的FM法，但當(dāng)樣本數(shù)據(jù)較少時(shí)，可用MC仿真法擴(kuò)充樣本容量，提高區(qū)間估計(jì)精度.

［1］DAI Yi，ZHOU Y F，JIA Y Z.Distribution of time between failures of machining center based on type I censored data［J］.Reliability Engineering and System Safety，2003，79:375-377.

［2］張海波，賈亞洲，周廣文.數(shù)控系統(tǒng)故障間隔時(shí)間分布模型的研究［J］.哈爾濱工業(yè)大學(xué)學(xué)報(bào)，2005，37 (2):198-200.

［3］WANG Yiqiang，YAM R C M，ZUO Mingjian.A multicriterion evaluation approach to selection of the best statistical distribution［J］.Computers＆Industrial Engineering，2004，47:165-180.

［4］張英芝，申桂香，賈亞洲，等.數(shù)控車床故障分布規(guī)律及可靠性［J］.農(nóng)業(yè)機(jī)械學(xué)報(bào)，2006，37(1):156-159.

［5］ZHANG Haibo，JIA Yazhou，ZHOU Guangwen.Time between failures model and failure analysis of CNC system［J］.Journal of Harbin Institute of Technology (New Series)，2007，14(2):197-201.

［6］張英芝，賈亞洲，申桂香，等.數(shù)控機(jī)床故障分布的灰關(guān)聯(lián)分析［J］.農(nóng)業(yè)機(jī)械學(xué)報(bào)，2004，35(6):195-197.

［7］王金武，劉家福，許仲祥.履帶式拖拉機(jī)可靠性與維修性的分析［J］.農(nóng)業(yè)機(jī)械學(xué)報(bào)，2004，35(4):81-83.

［8］AKAIKE H.Information theory and an extension of the maximum likelihood principle［C］//Proc Second Int Symp on Information Theory.Budapest:Akademiai Kiado，1973:267-281.

［9］SAYAMA S，SEKINE M.Weibull distribution and K-distribution of sea clutter observed by X-band radar and analyzed by AIC［J］.IEICE Transactions on Communications，2000，E83B(9):1978-1982.

［10］SCHWARZ G.Estimating the dimension of a model［J］.Ann Statist，1978(6):461-464.

［11］ADRIAN E，RAFTERY.Bayes factors and BIC:comment on"a critique of the bayesian information criterion for model selection"［J］.Sociological Methods Research，1999，27(3):411-427.

［12］ZHANG L F，XIE M，TANG L C.A study of two estimation approaches for parameters of Weibull distribution based on WPP［J］.Reliability Engineering and System Safety，2007，92:360-368.

［13］JOHNSON L G.Theory and technique of variation research［R］.New York:Elsevier Publishing Co，1964.

［14］SANGUN P，BALAKRISHNAN N，GANG Zheng. Fisher information in hybrid censored data［J］.Statistics and Probability Letters，2008，78:2781-2786.

［15］ESCOBAR L A，MEEKER W Q.The asymptotic equivalence of the fisher information matrices for type I and type II censored data from location-scale families［J］. Communications in statistics theory and methods，2001，30(10):2211-2225.