一種新的威布爾比例失效模型參數(shù)估計(jì)方法

張路朋，王 娟，王 攀，陳 亮

(1.國(guó)網(wǎng)河北省電力公司邯鄲供電分公司，河北 邯鄲 056000；2.國(guó)網(wǎng)河北省電力公司石家莊供電分公司，石家莊 050051)

2017-06-05

張路朋(1988-)，男，工程師，主要從事電力系統(tǒng)繼電保護(hù)工作。

一種新的威布爾比例失效模型參數(shù)估計(jì)方法

張路朋1，王 娟1，王 攀1，陳 亮2

(1.國(guó)網(wǎng)河北省電力公司邯鄲供電分公司，河北 邯鄲 056000；2.國(guó)網(wǎng)河北省電力公司石家莊供電分公司，石家莊 050051)

針對(duì)比例失效模型在參數(shù)估計(jì)中存在計(jì)算速度慢、收斂效果差以及初始數(shù)據(jù)要求高等缺陷，提出了基于極大似然估計(jì)方法的分步估計(jì)算法，首先對(duì)形狀參數(shù)和尺度參數(shù)進(jìn)行極大似然估計(jì)，然后將計(jì)算結(jié)果代入比例失效模型，再對(duì)狀態(tài)指標(biāo)回歸系數(shù)進(jìn)行極大似然估計(jì)，經(jīng)算例分析表明，該算法在計(jì)算速度、計(jì)算時(shí)間、數(shù)據(jù)樣本選擇方面具有較大的優(yōu)勢(shì)，證實(shí)了所提算法的有效性。

比例失效模型；參數(shù)估計(jì)；極大似然估計(jì)；分步估計(jì)算法

比例失效模型(Proportional Hazards Model，PHM)由Cox首先提出，最初用于醫(yī)學(xué)研究，隨著設(shè)備可靠性技術(shù)的不斷發(fā)展，逐步成為設(shè)備狀態(tài)評(píng)估與狀態(tài)維修的典型模型，PHM能有效地將設(shè)備狀態(tài)信息用于設(shè)備的狀態(tài)維修，從而在考慮經(jīng)濟(jì)性和可靠性的前提下獲得設(shè)備的最大可用壽命，并進(jìn)行預(yù)知性維修或者更換決策。

目前應(yīng)用最為廣泛的是威布爾比例失效模型(Weibull Proportional Hazards Model, WPHM)，文獻(xiàn)[1-2]將WPHM應(yīng)用于風(fēng)電機(jī)組的狀態(tài)檢修中，通過采集風(fēng)機(jī)的溫度和振動(dòng)等數(shù)據(jù)，制定維修閾值，通過最小成本法計(jì)算出最佳維修時(shí)機(jī)，取得了較好的效果。

WPHM在應(yīng)用過程中，參數(shù)的估算始終是一大難題，目前應(yīng)用最多的方法是極大似然估計(jì)算法，文獻(xiàn)[3-5]利用極大似然估計(jì)算法，將WPHM各參數(shù)進(jìn)行整體估計(jì)，在參數(shù)估計(jì)過程中，采取的做法是利用歷史故障數(shù)據(jù)，將形狀參數(shù)、尺度參數(shù)、狀態(tài)指標(biāo)回歸系數(shù)同時(shí)計(jì)算，這一做法大大增加了參數(shù)估計(jì)的計(jì)算量。但這種算法對(duì)初始數(shù)據(jù)的要求較高，同時(shí)可能會(huì)出現(xiàn)結(jié)果不收斂的情況。為了避免這種情況的發(fā)生，對(duì)WPHM的參數(shù)估計(jì)算法進(jìn)行了改進(jìn)，提出了分步估計(jì)算法。所謂分步估計(jì)算法，即將WPHM形狀參數(shù)、尺度參數(shù)、狀態(tài)指標(biāo)回歸系數(shù)分開估計(jì)，這種方法收斂效果更好，計(jì)算速度更快，且對(duì)初始數(shù)據(jù)的要求較低。

1 威布爾比例失效模型

威布爾比例失效模型兼顧設(shè)備自然老化和狀態(tài)協(xié)變量沖擊兩方面的影響，在設(shè)備狀態(tài)檢修中得到廣泛應(yīng)用，如式(1)。

(1)

式中：h0(t)為基本失效函數(shù)，由運(yùn)行時(shí)間t唯一決定，在威布爾比例失效模型中，用威布爾故障率函數(shù)來表示，反應(yīng)變壓器隨時(shí)間的自然老化[6]；η是壽命參數(shù)；β>1是形狀參數(shù)；φ(Z)為狀態(tài)協(xié)變量函數(shù)，由狀態(tài)協(xié)變量集合Z唯一決定，反應(yīng)設(shè)備受狀態(tài)協(xié)變量集合Z的沖擊，設(shè)狀態(tài)協(xié)變量集合Z包含p個(gè)狀態(tài)協(xié)變量，分別用Z1～Zp表示，其回歸系數(shù)分別為γ1～γp。

提出的分步估計(jì)算法就是通過極大似然估計(jì)算法，對(duì)WPHM的參數(shù)η，β與[γ1，γ2，…，γp]分開估計(jì)。

首先對(duì)式(1)所示的威布爾比例失效模型進(jìn)行分析發(fā)現(xiàn)：

a. WPHM的基本失效函數(shù)h0(t)為威布爾分布，描述的是設(shè)備在理想條件下的自然劣化過程，其分布只取決與運(yùn)行時(shí)間t。

b. 狀態(tài)指標(biāo)函數(shù)φ(Z)描述的是設(shè)備在實(shí)際工況下受運(yùn)行環(huán)境影響的程度，其分布只取決于選取的狀態(tài)指標(biāo)集Z=[Z1，Z2，…，Zp]，即 φ(Z)受實(shí)際工況的影響。

c.h0(t)和φ(Z)在數(shù)學(xué)上沒有聯(lián)系。

以上結(jié)論進(jìn)一步表明，可將基本失效函數(shù)h0(t)中涉及到的形狀參數(shù)β和尺度參數(shù)η與狀態(tài)指標(biāo)函數(shù)φ(Z)中涉及到的γ=[γ1，γ2，…，γp]分開估計(jì)。采取先求取β和η，從而得到基本失效函數(shù)h0(t)，然后將其作為已知量代入式(1)，對(duì)γ=[γ1，γ2，…，γp]進(jìn)行計(jì)算，最終得到威布爾比例失效模型。

在可靠性模型中，進(jìn)行參數(shù)估計(jì)所用到的數(shù)據(jù)主要有故障數(shù)據(jù)和結(jié)尾數(shù)據(jù)，在參數(shù)估計(jì)過程中，暫不考慮結(jié)尾數(shù)據(jù)，只采用故障數(shù)據(jù)。參數(shù)求取分為β，η的求取和γ=[γ1，γ2，…，γp]求取，采用的算法是極大似然估計(jì)算法。

2 參數(shù)估計(jì)方法

2.1 形狀參數(shù)和尺度參數(shù)的求取

基本失效函數(shù)用威布爾失效率函數(shù)表示，表達(dá)式為

(2)

β，η的求取采取的是連續(xù)函數(shù)的極大似然估計(jì)方法，數(shù)理統(tǒng)計(jì)原理中強(qiáng)調(diào)，連續(xù)函數(shù)的極大似然估計(jì)采用概率密度函數(shù)，所以首先要求解威布爾分布的概率密度函數(shù)。根據(jù)可靠性知識(shí)[7-8]，可推出威布爾分布的概率密度函數(shù)為

(3)

基本失效函數(shù)h0(t)只取決于運(yùn)行時(shí)間，所以對(duì)其估計(jì)時(shí)，只需要統(tǒng)計(jì)設(shè)備的歷史故障時(shí)間，然后將其從小到大排列，得到歷史故障時(shí)間集合T={t1，t2，…，tp}， 的極大似然估計(jì)算法如下。

為了便于計(jì)算，令m=ηβ，得到似然函數(shù)如式(4)所示：

(4)

對(duì)式(4)兩端同時(shí)取對(duì)數(shù)，并分別求β,m的偏導(dǎo)，代入m=η-β，令其等于零可得：

(5)

采用牛頓法對(duì)上式進(jìn)行收斂計(jì)算，可以求解出β，η的具體值。

2.2 狀態(tài)協(xié)變量參數(shù)的估計(jì)

求取γ=[γ1，γ2，…，γp]，不但需要考慮故障時(shí)的運(yùn)行工況，還需要考慮正常時(shí)的運(yùn)行工況。有別于β，η的求取，γ=[γ1，γ2，…，γp]的求取采用的是散點(diǎn)極大似然估計(jì)，因此，不需要先求取概率密度函數(shù)。采用固定間隔時(shí)間采集的狀態(tài)指標(biāo)檢測(cè)值來對(duì)γ=[γ1，γ2，…，γp]進(jìn)行散點(diǎn)估計(jì)。設(shè)在ti時(shí)刻設(shè)備的狀態(tài)監(jiān)測(cè)值為Z(ti)，從開始投入運(yùn)行到故障截止，一共采集了n個(gè)樣本點(diǎn)數(shù)據(jù)，則似然函數(shù)可以寫成如式(6)的形式：

(6)

(7)

對(duì)上式兩邊同時(shí)取對(duì)數(shù)可得

(8)

上式中含有積分項(xiàng)，對(duì)γ求偏導(dǎo)較難，考慮到在連續(xù)監(jiān)測(cè)中，部件的狀態(tài)協(xié)變量Z(ti)在較小的一段時(shí)間內(nèi)變化幅度不大，可以把{Z(t)|t≥0}看作是右連續(xù)的階躍過程，于是采用把積分項(xiàng)寫成分段積分和的形式。假設(shè)0

(9)

分部積分可得

(10)

定義ev=exp[γ·Z(tv)](v=0,2,…,j)，則式(10)的積分項(xiàng)運(yùn)用分部積分求解，可進(jìn)一步推出

(11)

帶入式(8)可得

(12)

式(12)對(duì)γ=[γ1，γ2，…，γp]中的γk求偏導(dǎo)并令其等于零，由于β，η此時(shí)都為常數(shù)，于是可得

(13)

令上式等于零，采用牛頓法對(duì)上式進(jìn)行收斂計(jì)算，可以求解γ=[γ1，γ2，…，γp]的具體值。

各參數(shù)求出之后，代入式(1)可以得到威布爾比例失效模型。

3 算例分析

變壓器長(zhǎng)期運(yùn)行于過負(fù)荷狀態(tài)，內(nèi)部會(huì)發(fā)生油氣變化，此時(shí)需要時(shí)刻監(jiān)測(cè)變壓器的運(yùn)行狀態(tài)，以便制定合理的檢修策略。利用油色譜分析技術(shù)可以得到變壓器油中溶解的氫氣(H2)、一氧化碳(CO)、二氧化碳(CO2)、甲烷(CH4)、乙烯(C2H4)、乙烷(C2H6)、乙炔(C2H2)等氣體體積分?jǐn)?shù)[9-10]。在進(jìn)行算例分析時(shí)，利用變壓器的油色譜分析數(shù)據(jù)來計(jì)算變壓器的威布爾比例失效模型，可為制定變壓器檢修策略提供依據(jù)。

對(duì)某型號(hào)主變壓器的歷史故障數(shù)據(jù)和油色譜監(jiān)測(cè)數(shù)據(jù)進(jìn)行統(tǒng)計(jì)，表1為該型號(hào)主變壓器的8次歷史故障時(shí)間，表2為該型號(hào)第8次維修周期中的油色譜分析數(shù)據(jù)。運(yùn)用MATLAB軟件進(jìn)行計(jì)算。利用表1的歷史故障時(shí)間并結(jié)合2.1節(jié)方法，可對(duì)形狀參數(shù)和尺度參數(shù)進(jìn)行計(jì)算，計(jì)算結(jié)果如表3所示，利用表2的數(shù)據(jù)結(jié)合2.2節(jié)方法，可對(duì)該型號(hào)變壓器的以及H2、CO、CO2、CH4、C2H4、C2H6、C2H2的回歸系數(shù)進(jìn)行計(jì)算，計(jì)算結(jié)果如表3所示。利用分步計(jì)算法進(jìn)行計(jì)算時(shí)，計(jì)算參數(shù)γ=[γ1，γ2，…，γp]用時(shí)2 s，計(jì)算 用時(shí)14 s，整個(gè)計(jì)算過程用時(shí)16 s。

表1 某型號(hào)變壓器的歷史故障數(shù)據(jù)部分樣本

主變壓器序號(hào)12345678故障天數(shù)/d856106895514281222113113971646

表2 某型號(hào)變壓器第8次維修周期內(nèi)油色譜監(jiān)測(cè)數(shù)據(jù)部分樣本 μL/L

采樣時(shí)間φ(H2)φ(CO)φ(CO2)φ(CH4)φ(C2H4)φ(C2H6)φ(C2H2)2016-01-1510881563.00.30.90.42016-07-16161002893.10.41.03.72016-10-15291192223.50.61.04.32016-11-18311283974.20.81.14.82016-12-16391372114.51.11.06.42017-01-15441383285.31.01.77.2

表3 參數(shù)計(jì)算值(分步計(jì)算法)

參數(shù)名稱βηγ1γ2γ3γ4γ5γ6γ7參數(shù)值9.0372853.30.13420.09260.07430.16320.21530.17200.1682計(jì)算用時(shí)/s 2 14

利用整體估計(jì)方法計(jì)算結(jié)果如表4所示，計(jì)算用時(shí)為37 s。

表4 參數(shù)計(jì)算值(整體計(jì)算法)

參數(shù)名稱βηγ1γ2γ3γ4γ5γ6γ7參數(shù)值9.0423853.50.13380.09270.07490.16520.21430.17260.1692計(jì)算用時(shí)/s37

表4與表3對(duì)比發(fā)現(xiàn)，利用提出的分步估計(jì)算法求解結(jié)果與整體估計(jì)算法結(jié)果相差較小，可滿足工程誤差要求，但是分步計(jì)算法總體用時(shí)16 s，而整體計(jì)算方法用時(shí)37 s，是分步計(jì)算方法的兩倍多，可見分步估計(jì)算法在計(jì)算速度、計(jì)算時(shí)間方面具有較大的優(yōu)勢(shì)，而這種優(yōu)勢(shì)隨著數(shù)據(jù)量的增加會(huì)逐步擴(kuò)大。同時(shí)在數(shù)據(jù)樣本選擇方面，分步計(jì)算法對(duì)數(shù)據(jù)要求不高，計(jì)算更容易收斂，而整體計(jì)算方法對(duì)數(shù)據(jù)要求較高，計(jì)算結(jié)果容易發(fā)散。

將表3計(jì)算結(jié)果代入式(1)中，可以得到變壓器的威布爾比例失效函數(shù)，可根據(jù)該失效函數(shù)制定相應(yīng)的狀態(tài)檢修策略，為主變運(yùn)行維護(hù)人員提供檢修建議。

4 結(jié)束語

對(duì)傳統(tǒng)比例失效模型參數(shù)估計(jì)方法進(jìn)行了改進(jìn)，提出一種基于極大似然估計(jì)算法的新的參數(shù)估計(jì)方法，即分步估計(jì)算法，該算法認(rèn)為威布爾比例失效模型的參數(shù)與在數(shù)學(xué)上互不相關(guān)，可以利用極大似然估計(jì)算法進(jìn)行分開計(jì)算。算例分析表明，分步估計(jì)算法在計(jì)算速度、計(jì)算時(shí)間、數(shù)據(jù)樣本選擇方面具有較大的優(yōu)勢(shì)，證實(shí)了所提方法的有效性。

[1] Ardine，Banjevicd．Optimization a mine haul truck wheel motors' condition monitoring program use for proportional hazards modeling [J]．Journa l of Quality in Maintenance Engineering，2001，7(4)：286 - 301．

[2] Xiang Wu，Sarah M. Ryan．Optimal Replacement in the Proportional Hazards Model With Semi-Markovian Covariate Process and Continuous Monitoring[J]． IEEE Tanscctions on Reliability，2011，60(3)：580-589．

[3] 張繼全，趙洪山．基于比利失效模型的風(fēng)機(jī)齒輪箱軸承檢修[J]．電源技術(shù)，2011，35(5)： 570-573．

[4] 陳寒雨．基于設(shè)備故障率評(píng)估的視情維修研究[D]．成都：電子科技大學(xué)，2011．

[5] 馬維青，于瑤章．基于比例失效模型的配電網(wǎng)設(shè)備狀態(tài)檢測(cè)[J]．信息技術(shù)，2015，(15)：206-209．

[6] 趙洪山，張路朋．基于可靠度的風(fēng)電機(jī)組預(yù)防性機(jī)會(huì)維修策略[J]．中國(guó)電機(jī)工程學(xué)報(bào)，2014，34(22)：3777-3783．

[7] 張路朋．風(fēng)電機(jī)組的狀態(tài)機(jī)會(huì)維修策略研究[D]．保定：華北電力大學(xué)，2015．

[8] 張路朋，趙洪山．基于時(shí)間延遲的風(fēng)機(jī)齒輪箱狀態(tài)優(yōu)化維修[J]．中國(guó)電力，2014，47(11)：108-111．

[9] DL/T 722-2014．變壓器油中溶解氣體分析和判斷導(dǎo)則[S]．

[10] 李守學(xué)，司昌健，張春豐，等．某500 kV變壓器油色譜異常監(jiān)測(cè)及缺陷診斷[J]．吉林電力，2017，45(1)：44-46．

A New Method of Parameter Estimation for Proportional Hazards Model

Zhang Lupeng1，Wang Juan1,Wang Pan1,Chen Liang2

(1.State Grid Hebei Electric Power Corporation,Handan Power Supply Branch,Handan 056000,China；2.State Grid Hebei Electric Power Company Shijiazhuang Power Supply Branch,Shijiazhuang 050051,China)

Aiming at the slow computation speed,poor convergence effect and high initial data requirements in parameter estimation of proportional hazards model,this paper proposes step estimation algorithm based on maximum likelihood estimation.Firstly,this paper has shape parameter and scale parameter estimated by maximum likelihood estimation,then it substitutes the result into proportional hazards model,and then it has regression coefficient of condition index estimated by maximum likelihood estimation.Through the calculation example analysis,step estimation algorithm is more effective in computation speed,computation time and data sample selection.

proportional hazards model;parameter estimation;maximum likelihood estimation;step estimation algorithm

TM41

B

1001-9898(2017)05-0018-03

本文責(zé)任編輯：靳書海