基于性能退化預測的數控機床預防維修方法

代愽超 張英杰 李陽帆 陳 波

西安交通大學機械工程學院，西安，710049

0 引言

隨著制造業整體實力的提高，通過壽命試驗或加速壽命試驗獲取機床故障數據日趨困難且成本昂貴，為了保證數控機床加工質量的持續穩定性，采用合理的預防性維修方法已成為提高設備利用率和實現資產效率最大化的有效途徑[1]。

生產控制方面，統計過程控制(statistical process control, SPC)作為一種集測定、記錄、分析、控制等于一體的實用技術[2-4]，廣泛應用于機床等設備的可靠性監測。其中普通控制圖[2,5-6]、預控圖[7]和選控圖[8-9]直接關注評測目標的狀態參數變化，而在狀態參數檢測設備不完善的生產線上，則可利用不合格品率控制圖進行生產過程的控制[10-11]。

維護維修方面，大多數研究者將 “維修如新”作為假設，其原理簡單、操作方便[12-14]。考慮到維修程度不同所帶來的效果差異，“不完全維修”方面的研究也取得了一定的進展[15-16]。但“不完全維修”的維修程度難以掌控，因此一些學者認為將數控機床的維修視為“維修如舊”更為實用[17-18]。

目前在性能退化分析和預防維修決策方面的研究大多為單一方法研究，針對從狀態預測到維修決策的全過程的方法依然較少，難以滿足企業的實際需求。為此，本文提出了一種基于系統性能退化預測的數控機床預防維修方法。利用不合格品率控制圖和Markov狀態轉移矩陣預測發生異常的時間、位置和原因，建立了Wiener退化模型來預測設備剩余壽命，計算并比較了維修和換新的后期收益效率，以實現對維修策略的指導。

1 基于SPC-Markov模型的機床狀態分析

1.1 不合格品率控制圖相關參數的計算

不合格品率控制圖是建立在二項分布基礎上的一種統計過程控制技術[10]。設樣本的不合格品率為p，采用統計方法對其進行估計，表達式如下：

(1)

考慮到各樣本中的產品數量可能不同，故對式(1)進行修正，得到

(2)

根據六西格瑪(6σ)質量管理方法，結合二項分布的性質可計算得到中心線、上控制限和下控制限，其表達式分別如下：

(3)

1.2 基于Markov過程的異常概率預測

考慮到當前維修部位和維修水平會對下一次異常的發生造成影響，采用Markov狀態轉移矩陣進行各種異常模式的概率預測。

Markov狀態轉移矩陣如下：

(4)

其中，z為同一時刻最多可能的異常種類數；PIJ為t時刻處于I狀態，到t+1時刻變為J狀態的概率，且滿足：

(5)

t時刻的狀態空間向量

Et=(e1,e2,…,ez)

(6)

其中，e1,e2，…，ez的取值分別表示每一種異常是否發生，當取值為1時，表示該種異常發生；當取值為0時，該種異常不會發生。

則可求得t+1時刻各種異常情況的發生概率向量：

Gt+1=EtM

(7)

2 性能退化模型的建立及維修策略分析

2.1 基于Wiener過程的性能退化模型

Wiener退化過程模型[19-21]定義為

y(t)=a+μt+σBB(t)

(8)

式中，y(t)為t時刻的退化量測量值；a為初始退化量；μ為漂移系數；σB為擴散系數；B(t)為標準布朗運動值。

設i(i=1,2,…,N)為退化試驗樣本序號，N為退化試驗樣本總數；j(j=1,2,…,mi)為測量次數序號，mi為各樣本的測量次數。進而可得到如下試驗數據模型：

yij=ai+μitij+σBiBij

(9)

式中，tij、Bij、yij分別為第i個樣本第j次的測量時間、布朗運動值和退化量測量值；ai、μi、σBi分別為第i個樣本的初始退化量、漂移系數和擴散系數。

采用各樣本分離的方法計算參數值，基本流程見圖1。

圖1 Wiener退化計算流程Fig.1 Evaluation process of Wiener degradation model

利用Wiener過程的高斯獨立增量性可簡化計算。令

(10)

j=2,3,…,mi

則

(11)

式中，Δtij、Δyij分別為樣本i第j-1次和第j次測量之間的時間間隔和退化量增量。

建立似然函數：

(12)

先求得對數似然函數L(ai,μi,σBi)，再根據如下邊界條件：

(13)

最終解得各樣本參數的估計值分別為

(14)

(15)

(16)

將失效閾值記為Df，將y(t)首次達到Df的時間定義為機床壽命T，即

T=inf(t|y(t)=Df)t≥0

(17)

式中,inf(·)表示求集合中最大下界的操作。

退化初值a通常小于閾值Df，可靠度可描述為機床壽命T大于規定時間t的概率，即

(18)

實際求解時，選取參數a、μ、σB中變異系數較大的參數服從正態分布，以避免出現多重積分[20]。變異系數計算表達式如下：

(19)

式中,X為具體的變量(即本例中的a、μ和σB)；var(X)為X的方差；E(X)為X的期望。

2.2 基于剩余壽命和收益的維修策略分析

以“維修如舊”為前提，選取可靠度退化閾值Dr，認為系統平均壽命為退化曲線上的可靠度降低至Dr所需的總時間tA。則剩余壽命s(t)可表示為

s(t)=tA-t-t-

(20)

式中，t-為當前研究時段之前的總運行時間。

假設在機床正常工作期間，其日常運行的基礎費用cm(單位時間)不變，則進行零部件維修和換新的后期收益效率分別為

(21)

(22)

式中，v為產品的出廠價格；w為單位時間產量；cr為每件產品所需原材料成本；Cf為維修費用；Cd為換新費用；tR為平均維修時間(MTTR)；tD為零部件更換平均用時(MTTD)。

通過比較式(21)、式(22)，將出現如下4種情況：

(1)ηf≥ηd且tR≤tD。ηf≥ηd表明在后續壽命內，相較于零部件換新，零部件維修產生的單位時間內的收益更高；又滿足tR≤tD，表明維修造成的誤工影響小。因此，宜進行關鍵零部件維修。

(2)ηf≥ηd且tR>tD。若當前生產任務緊迫則可先進行零部件更換，并在后續生產壓力小的時候，對換下的零部件進行維修以備用。

(3)ηf<ηd且tR≤tD。若生產任務輕松，則直接更換零部件；若急需趕工且預期修一次能完成趕工任務，則建議維修，否則建議盡早更換零部件。

(4)ηf<ηd且tR>tD。此時更適合進行關鍵零部件的換新。

當tR和tD分別比s(t)和tA小很多時，利用式(21)、式(22)計算收益效率時可略去tR和tD以簡化計算。此時若要滿足ηf≥ηd，只需滿足：

CftA≤Cds(t)

(23)

3 實例應用

3.1 生產線機床加工過程控制及狀態預測

國內某活塞生產線的主要加工流程見圖2。

圖2 某活塞生產線加工流程Fig.2 Machining process of piston production line

本研究收集了該生產線30天(實際工作24天，第24天因檢測出了大量返修品和報廢品而被迫進行停機維修)的不合格品數據，選取前22天的數據，并計算得到不合格品率控制圖參數，見表1，其中中心線值為0.049 6。

將各參數的計算結果繪制成控制圖見圖3，可以看出，由于每日產量的波動導致控制圖的上下限不斷變化，因此選擇如下3個通用準則進行判異：①樣本點距中心線大于3個標準差；②連續9點在中心線同一側；③連續14點上下交錯。

由圖3可以看出，出現了4個異常點，其中點4、點21和點22均為數據超界點(情況①)，而點14則是因為連續超過9個點同側(情況②)。由分析可知，點4為超下界點，而點6～14同樣處于中心線的下側，它們在統計規律上異常，但不合格品率偏低以及長期處于偏低態勢都是對實際情況有利的現象，因此應查明引起這種情況的原因，并力圖維持該狀態以保證生產質量，但也要檢測和分析是否是由檢測偏差引起的。此外，點21為超上界點，表明系統可能出現了某種故障，點22的超界又加劇了這一可能性，因此預測系統將在第23個工作日發生異常。

表1 不合格品率控制圖參數計算結果Tab.1 The calculation result of parameters for unqualified product rate control chart

圖3 不合格品率控制圖的預測結果Fig.3 Prediction results of unqualified product rate control chart

企業并未及時采取措施，且已將后續樣本點添加到已有控制圖中，結果見圖4。

圖4 不合格品率控制圖的效果驗證Fig.4 Effect validation of unqualified product rate control chart

圖4是在圖3的基礎上繪制得到的，且計算中心線保持不變。由圖4可以看出，點23回到了上下控制限范圍內，但點22到點23是從超上界回落到中心線LC以下，跨度很大，生產過程極不穩定；點24則再次超出了上控制限。整個生產過程的不合格品率達到25%，報廢量和返修量過大，企業無法接受。

利用Markov狀態轉移矩陣進行異常問題預測性診斷時，考慮到每臺機床都有很多種潛在的異常類型，因數據量十分龐大，故無法詳盡羅列完整的矩陣，本文僅分析下一次的異常最可能發生在哪一臺機床。

收集得到該生產線的372條歷史異常數據(實際狀態轉移371次)，見表2，其中，第Mi行Mj列的元素表示前一次異常發生在機床Mi，隨后一次異常發生在機床Mj的情況出現的頻數。

表2 各機床發生異常的轉移關系Tab.2 Abnormity transmission relationship between machine tools

計算得到Markov狀態轉移矩陣：

M=

(24)

因收集不合格品數據前最近一次的異常發生在機床M7，故令該時間節點為t，則狀態空間向量可表示為

Et=(0,0,0,0,0,0,1,0)

(25)

則t+1時刻各機床發生異常的概率向量可表示為

Gt+1=EtM=
(0.043,0.130,0.239,0.022,0.043,0.022,0.130,0.370)

(26)

依據式(26)，機床M8的異常概率最高，達到37%；機床M3的異常概率為23.9%以及機床M2、M7的異常概率均為13%，上述異常概率相對其他機床的異常概率也較高。由此可知，通過圖3預測得到第23個工作日可能發生異常時，應優先對機床M8進行排查，若機床M8正常，則再檢查機床M3、M2、M7，最后才檢查其他機床。

在本生產實例中，正是M8機床的定位問題導致點22～24之間不合格品率的連續大幅度波動，這與模型預測的大概率事件相符合，表明該方法具有一定的準確性，可以引導異常排查，提高檢修效率。

3.2 零部件剩余壽命預測及維修策略分析

根據預測第23個工作日發生異常，實際已工作22天(每天三班輪換，每班8 h)，即528 h。且預測最可能出現問題的是機床M8，若及時采取了相應檢修措施，就會發現機床M8的定位元件出現了問題(但實際上是被迫停機后才檢測出了定位問題)。

采用M8同類機床的定位精度數據進行分析。現有5臺定位失效閾值Df=300 μm的樣本機床，分別記作Si(i=1,2,3,4,5)，其定位精度檢測結果見圖5。

圖5 定位精度的性能檢測數據Fig.5 Performance testing data for positioning accuracy

利用式(14)～式(16)計算得到各樣本初始退化量、漂移系數、擴散系數的估計值，見表3。表3中，樣本5的初始退化量出現了負值，但考慮到有的機床需先使用一段時間后，定位精度才會提高到最佳狀態，因此這里不做改動。

表3 各樣本參數估計Tab.3 Parameter estimation of each sample

表4 參數估計值的變異性分析Tab.4 Variability analysis of the estimation values for parameters

(27)

(28)

式中，μa、σa分別為a的均值和標準差。

依據式(18)計算可靠度函數：

(29)

(30)

式中，φ(a)為a的概率密度函數。

取式(29)的積分區間為(-3σa,+3σa)，為補償被忽略的區間，增加修正系數：

(31)

在MATLAB軟件中繪制得到機床M8的可靠度退化曲線，見圖6。

圖6 機床M8的可靠度退化曲線Fig.6 Reliability degradation curve of M8

根據企業要求選取可靠度閾值Dr=0.6，由圖6可知tA=1 300 h，又查得該定位元件在數據收集前實際已運行約21天，取t-=504 h，則有

s(t)=tA-t-t-=268 h

(32)

維修和換新所需時間均在3 h左右，可依據式(23)進行比較，得到優先維修的條件為

Cd/Cf≥4.85

(33)

對于引起本次故障的關鍵零部件，更換一次的平均費用約為2 360元，維修費用約為400元，滿足式(33)。由2.2節中情況(1)可知：從經濟性和時間效率兩方面考慮，此時對關鍵定位元件進行維修比換新具有更好的綜合效益。

求得tA后，根據式(20)，在當前成本條件下若要滿足式(23)，只需滿足：

t+t-≤1 079.66 h

(34)

這表明引起本次定位問題的關鍵零部件在總運行時間不超過1 079.66 h的情況下，優先選擇維修是更加經濟的選擇。

4 結論

(1)提出了一種基于系統性能退化預測的數控機床預防維修方法，利用不合格品率控制圖實現了生產線上機床設備的統計過程控制，降低了對專用檢測工具的依賴性。

(2)基于Wiener退化模型，利用其增量的高斯獨立性，采用各樣本分離的方法求解模型參數，可較為準確地評估關鍵零部件的剩余壽命。

(3)基于“維修如舊”的假設，結合生產效益指導預防維修策略，且符合企業的實際利益需求。

(4)對國內某汽車活塞生產線數控機床的分析結果表明，所提方法的預測結果與實際情況基本一致，所構建的決策模型可為設備的預防性維修提供一定的參考。