吳健飛,劉勤明,呂文元,葉春明
(上海理工大學(xué) 管理學(xué)院, 上海 200093)
隨著現(xiàn)代化工業(yè)技術(shù)的蓬勃發(fā)展,機械設(shè)備的可靠性和維護水平已成為保證系統(tǒng)正常運行的必要條件.在實際生產(chǎn)過程中,隨著設(shè)備使用時間的增加,設(shè)備的性能在不斷衰退,若不進行有效的維護,設(shè)備可能會失效或者發(fā)生故障,從而會影響企業(yè)的生產(chǎn)效益.所以在設(shè)備的健康狀態(tài)嚴(yán)重之前,對設(shè)備進行實時的處理,可以避免一些成本高昂的維修,降低設(shè)備的維護成本,提高設(shè)備的利用率.
設(shè)備的健康預(yù)測是有效實施維護策略的關(guān)鍵,由于數(shù)據(jù)監(jiān)測和計算機技術(shù)的進步,關(guān)于設(shè)備健康預(yù)測理論的研究引起了國內(nèi)外學(xué)者的關(guān)注[1-3].其中,基于馬爾可夫鏈的預(yù)測方法顯得尤為突出.例如:Du Ying等人將隱馬爾可夫模型(HMM)應(yīng)用到潤滑油的壽命預(yù)測中,在該方法中,假設(shè)狀態(tài)的駐留時間服從指數(shù)分布,但這與實際并不相符[4].針對HMM中的不足,許多學(xué)者又研究了基于隱半馬爾可夫模型(HSMM)的預(yù)測和診斷,HSMM模型沒有HMM模型中不符合實際的假設(shè),因此對實際問題具有更強的建模和分析能力.像Akram Khaleghei將HSMM模型用于部分觀測的故障系統(tǒng)的壽命預(yù)測中,并利用EM算法估計出HSMM中的參數(shù)[5].Dong等人使用隱半馬爾可夫模型(HSMM)進行機械故障診斷和故障預(yù)測,結(jié)果表明其故障識別率比HMM提高了,且基于HSMM的剩余壽命預(yù)測模型效果也較好[6,7].Wang等人針對設(shè)備運行狀態(tài)識別與故障預(yù)測問題,提出一種基于時變轉(zhuǎn)移概率的隱半Markov模型,并通過滾動軸承的實驗,證明了基于該模型的設(shè)備健康預(yù)測方法比傳統(tǒng)的隱半Markov模型方法更加有效[8].
然而目前的研究大都是基于樣本數(shù)據(jù)準(zhǔn)確的情況下對設(shè)備進行的健康預(yù)測[9-12],而關(guān)于在數(shù)據(jù)不確定情況下的設(shè)備健康預(yù)測文獻較少,或在數(shù)據(jù)不確定情況下對設(shè)備進行健康預(yù)測的研究很難得到令人滿意的結(jié)果[13-17].不確定數(shù)據(jù)是指在工程實際中由于噪聲、擾動、非理性測量儀器以及人為等因素導(dǎo)致數(shù)據(jù)在產(chǎn)生過程中出現(xiàn)異常的情況.
不確定的數(shù)據(jù)會導(dǎo)致樣本數(shù)據(jù)的不準(zhǔn)確,不僅增加了分析數(shù)據(jù)的難度,而且會導(dǎo)致數(shù)據(jù)挖掘結(jié)果產(chǎn)生偏差,不能準(zhǔn)確的對設(shè)備進行健康預(yù)測.因此,針對數(shù)據(jù)不確定的情況,本文基于馬爾可夫模型,首先利用DS證據(jù)理論構(gòu)建狀態(tài)識別框架.其次,用區(qū)間數(shù)來表示不確定的數(shù)據(jù),并利用區(qū)間數(shù)之間的距離和相似度作為產(chǎn)生BPA的證據(jù),為了使預(yù)測結(jié)果更加準(zhǔn)確,采用Pignistic概率轉(zhuǎn)換將BPA轉(zhuǎn)化為基礎(chǔ)狀態(tài)的概率分布.最后,通過一個案例研究驗證了提出方法的有效性.
設(shè)θ為一有限集,2θ為θ上的冪集,m是2θ到[0,1]上的函數(shù),若滿足:
m(φ)=0,∑A?2θm(A)=1
(1)
稱m為θ上的基本概率分配函數(shù),?A∈2θ,m(A)稱為A的基本概率賦值,反映了對A的信任程度.如果A的基本概率賦值m(A)>0,則把A稱為θ的一個焦元.θ通常也稱為識別框架,表示在條件E下所有可能結(jié)論的基礎(chǔ)命題的有限集,θ的一個子集A,即2θ中的元素,可以理解為一個命題.
隨機序列Xn在任一時刻n,它可以處在狀態(tài)O1,O2,…,ON,且它在l+k時刻所處的狀態(tài)為ql+k的概率,只與它在l時刻的狀態(tài)ql有關(guān),而與l時刻以前所處狀態(tài)無關(guān),即:
P(Xl+k=ql+k|Xl=ql,X1=q1)=P(Xl+k=ql+k|Xl=ql)
(2)
其中q1,q2,…,ql,ql+k∈(O1,O2,…,ON),則稱Xn為Markov鏈.
并且稱Pij(l,l+k)=P(ql+k=Oj|ql=Oi)為k步轉(zhuǎn)移概率,其中1≤i,j≤N,l,k為正整數(shù).當(dāng)Pij(l,l+k)與l無關(guān)時,稱這個Markov鏈為齊次Markov鏈,此時Pij(l,l+k)=Pij(k).當(dāng)k=1時,Pij(1)表示單步轉(zhuǎn)移概率,簡稱為轉(zhuǎn)移概率.
根據(jù)相同設(shè)備的歷史運行狀況可以建立狀態(tài)識別框架θ={θ1,θ2,…,θn},其中θ1,θ2,…,θn表示設(shè)備的n個基礎(chǔ)狀態(tài).識別框架θ的冪集為2θ,令2θ中的每個元素分別代表設(shè)備運行中的一個狀態(tài),則設(shè)備的狀態(tài)相應(yīng)的在原來的n個基礎(chǔ)狀態(tài)上增加了,然后根據(jù)實際情況可以確定每個狀態(tài)對應(yīng)的數(shù)據(jù)范圍.
另外,若A=[a1,a2],B=[b1,b2]分別表示兩個區(qū)間數(shù),則它們之間距離的平方為:
(3)
其中,D(A,B)表示區(qū)間數(shù)A和B之間的距離.雖然采集到的數(shù)據(jù)很多是確定的,但是仍可以用區(qū)間數(shù)來表示,例如10可以表示成[10,10].因為每個狀態(tài)都有對應(yīng)的數(shù)據(jù)范圍,基于公式(3),計算采集到的數(shù)據(jù)和狀態(tài)之間的距離.然后基于區(qū)間數(shù)之間的距離計算區(qū)間數(shù)之間的相似度.區(qū)間數(shù)之間的相似度定義為:若X=[x1,x2]和Y=[y1,y2]分別表示兩個區(qū)間數(shù),則這兩個區(qū)間數(shù)的相似度為:
(4)
最后,標(biāo)準(zhǔn)化計算出的區(qū)間數(shù)相似度就得到采集數(shù)據(jù)的基本概率賦值,具體計算過程見圖1.
圖1 基于區(qū)間數(shù)的BPA計算過程Fig.1 BPA calculation process based on interval number
通過一個例子來說明整個計算過程:若狀態(tài)E對應(yīng)的數(shù)據(jù)范圍為[0,5],狀態(tài)F對應(yīng)的數(shù)據(jù)范圍為(5,10],給定一個區(qū)間數(shù)C=[3,6],則區(qū)間C的基本概率賦值見表1.
獲得每個采集數(shù)據(jù)的基本概率賦值后,計算各個狀態(tài)之間的轉(zhuǎn)移概率,如下:
(5)
其中m(i)t表示狀態(tài)i在第t個時間點的基本概率賦值,n表示馬爾可夫鏈的長度.所以狀態(tài)轉(zhuǎn)移概率矩陣P=[Pij].在這里,狀態(tài)之間的轉(zhuǎn)移是指識別框架中的所有狀態(tài).
利用得到的狀態(tài)轉(zhuǎn)移概率矩陣P和采集到的最后一個時間點上數(shù)據(jù)的基本概率賦值m,可以計算出下一個時間點上數(shù)據(jù)的BPA,即:
m預(yù)測=m·P
(6)
其中m=m(i),i∈2θ表示采集到的最后一個時間點上數(shù)據(jù)的基本概率賦值.
雖然得到了下一個時間點上數(shù)據(jù)的基本概率賦值,但有時想通過基本概率賦值直接進行預(yù)測還是不容易的.為了更準(zhǔn)確的進行健康預(yù)測,本文采用Pignistic概率轉(zhuǎn)換將識別框架中所有狀態(tài)的基本概率賦值轉(zhuǎn)化為基礎(chǔ)狀態(tài)的概率分布.Pignistic概率轉(zhuǎn)換是為了重新分配系統(tǒng)已獲得的各命題的信度值,以得到更可靠的決策依據(jù).通常使用的分配方式是平均分配法,即認(rèn)為每個元素出現(xiàn)的概率相同,因此多元素命題的BPA值被平均分配到所包含的元素中:
(7)
m(A)為A的基本概率賦值,|A|表示A中元素的個數(shù).得到了基礎(chǔ)狀態(tài)的概率分布即可知道下一時間點設(shè)備狀態(tài).
通過機械設(shè)備的健康預(yù)測來驗證所提出方法的有效性.機械設(shè)備在運行使用過程中會經(jīng)歷不同狀態(tài)的演化,假設(shè)這個演化過程服從馬爾可夫過程,并且設(shè)備的狀態(tài)可以用下面三個狀態(tài)來表示:健康狀態(tài){0},退化狀態(tài){1},失效狀態(tài){2}.表2是在一段時間內(nèi)通過傳感器從設(shè)備上采集到的20個數(shù)據(jù),其中在第3、6、13、16個時間時由于傳感器的波動導(dǎo)致采集的數(shù)據(jù)是不確定的,我們用區(qū)間數(shù)來表示這些不確定的數(shù)據(jù).
表2 一段時間內(nèi)采集到的數(shù)據(jù)Table 2 Data collected over a period of time
根據(jù)相同設(shè)備的歷史信息可以得到三個基礎(chǔ)狀態(tài)對應(yīng)的數(shù)據(jù)范圍,見表3.
表3 3個狀態(tài)下對 應(yīng)的數(shù)據(jù)范圍Table 3 Corresponding data range in three states
考慮到在某時間點采集到的數(shù)據(jù)不能表示設(shè)備以概率1處于某個狀態(tài),即數(shù)據(jù)18可能表示設(shè)備以概率0.4處于健康狀態(tài){0},以概率0.6處于退化狀態(tài){1}.所以通過歷史信息得到的3個狀態(tài)的據(jù)范圍并不很合理.利用DS證據(jù)理論解決該問題.
因為采集到的所有數(shù)據(jù)都可以表示設(shè)備處于3個基礎(chǔ)狀態(tài):健康狀態(tài){0},退化狀態(tài){1},失效狀態(tài){2}.所以可以構(gòu)建識別框架θ={0,1,2},得到θ的子集為{0},{1},{2},{0,1},{0,2},{1,2},{0,1,2},φ,在一個證據(jù)識別框架中這些子集分別代表一個狀態(tài),這些狀態(tài)及對應(yīng)的數(shù)據(jù)范圍見表4.
表4 識別框架中的狀態(tài)及這些狀態(tài)對應(yīng)的數(shù)據(jù)范圍Table 4 Identify the states of the frame and the range of data corresponding to those states
考慮到在實際情況中,設(shè)備不可能既處于健康狀態(tài)0又處于失效狀態(tài)2,所以狀態(tài){0,2}沒有對應(yīng)的數(shù)據(jù)范圍,同理,狀態(tài){0,1,2}和φ也是.所以最后的有效狀態(tài)為{0},{0,1},{1},{1,2},{2}.在狀態(tài)和對應(yīng)的數(shù)據(jù)范圍確定之后,可以得到20個采樣時間點的狀態(tài)分布,見圖2.
圖2 20個采樣數(shù)據(jù)點的狀態(tài)分布Fig.2 State distribution of 20 sampling data points
根據(jù)圖1的計算過程,可以獲得20個數(shù)據(jù)的基本概率賦值.以第6個采樣時間點上的數(shù)據(jù)為例來表示計算過程.因為第6個時間點上的數(shù)據(jù)是不確定的,用區(qū)間數(shù)來表示.識別框架中有效狀態(tài)對應(yīng)的數(shù)據(jù)范圍也用區(qū)間數(shù)的形式表示,所以第6個時間點上的數(shù)據(jù)和有效狀態(tài)對應(yīng)的數(shù)據(jù)范圍之間的距離可以利用公式(3)計算出來;然后根據(jù)公式(4),可以計算區(qū)間數(shù)之間的相似度;最后標(biāo)準(zhǔn)化相似度就可以得到第6個時間點上的數(shù)據(jù)的BPA.結(jié)果見表5.
表5 第6個時間點上的數(shù)據(jù)的BPATable 5 BPA of data at sixth time points
對剩下的其余所有時間點上的數(shù)據(jù)重復(fù)上述過程,可以獲得每個時間點上的數(shù)據(jù)的基本概率賦值,見表6.
表6 20個時間點上數(shù)據(jù)的BPA Table 6 BPA of data at 20 time points
由表6可以看出每個數(shù)據(jù)的BPA都不為0.所以公式(5)可以推導(dǎo)為:
(8)
通過上述獲得的第20個時間點的BPA和狀態(tài)轉(zhuǎn)移概率矩陣,利用公式(6)可以計算第21個時間點的BPA:
m(21)=m(20)·P=(0.0922 0.1288 0.1487 0.3472 0.2831)
然后利用公式(7)將第21個時間點的BPA轉(zhuǎn)化為基礎(chǔ)狀態(tài)的概率分布,即:
BetP({0})=0.0922+0.1288/2=0.1566
BetP({1})=0.1288/2+0.1487+0.3472/2=0.3867
BetP({2})=0.3472/2+0.2831=0.4567
最后得到P({0},{1},{2})=(0.1566,0.3867,0.4567).所以在第21個時間點上設(shè)備最可能處于失效狀態(tài).在實際中,相同設(shè)備在第21個時間上的確已經(jīng)失效,這與預(yù)測結(jié)果相符合.
由于第3、6、13、16個時間點上的數(shù)據(jù)是不確定的,我們利用matlab軟件隨機生成4個數(shù)據(jù),然后基于本文的方法得到狀態(tài)轉(zhuǎn)移概率矩陣:
利用公式(6)可以計算第21個時間點的BPA:
m(21)=m(20)·P=(0.0803 0.1221 0.1938 0.3571 0.2427)
然后利用公式(7)得到P({0},{1},{2})=(0.1414,0.4334,0.4213).最后結(jié)果表明設(shè)備更可能處于健康狀態(tài)1,這與實際情況不符合,所以利用區(qū)間數(shù)來處理不確定數(shù)據(jù)能更精確的預(yù)測設(shè)備的健康狀態(tài).
另外,為了驗證本文所提方法的改進性,我們利用文獻[15]中的方法來進行比較.首先采用拉依達準(zhǔn)則法剔除樣本數(shù)據(jù)中不確定的數(shù)據(jù).拉依達準(zhǔn)則如下:有在線監(jiān)測數(shù)據(jù)序列X={x(1),x(2),…,x(n)},若采樣點x(i)滿足公式(9),則認(rèn)為x(i)為不確定數(shù)據(jù)應(yīng)剔除.
|x(i)-x|>3σ
(9)
本文提出了一種針對數(shù)據(jù)不確定下如何處理設(shè)備健康預(yù)測的新方法.提出的DS-MM框架能夠有效對設(shè)備的健康狀態(tài)進行預(yù)測.考慮到區(qū)間數(shù)的性質(zhì),利用區(qū)間數(shù)來表示不確定的數(shù)據(jù),并基于區(qū)間數(shù)之間的距離和相似度來作為產(chǎn)生BPA的證據(jù),為了使預(yù)測結(jié)果更加可靠,利用Pignistic概率轉(zhuǎn)換將BPA轉(zhuǎn)化為基礎(chǔ)狀態(tài)的概率分布.最后,通過設(shè)備健康預(yù)測的案例分析驗證了該方法的有效性.
本文的新方法能夠有效解決機械設(shè)備健康預(yù)測中存在的數(shù)據(jù)不確定情況.但是本文僅給出了基于DS-MM模型的設(shè)備健康預(yù)測,而設(shè)備健康預(yù)測的目的是進行設(shè)備維修決策,因此下一步是根據(jù)預(yù)測結(jié)果進行設(shè)備預(yù)測維修決策的研究.