代愽超 張英杰 李陽帆 陳 波

西安交通大學(xué)機械工程學(xué)院，西安，710049

0 引言

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

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

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

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

1 基于SPC-Markov模型的機床狀態(tài)分析

1.1 不合格品率控制圖相關(guān)參數(shù)的計算

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

(1)

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

(2)

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

(3)

1.2 基于Markov過程的異常概率預(yù)測

考慮到當(dāng)前維修部位和維修水平會對下一次異常的發(fā)生造成影響，采用Markov狀態(tài)轉(zhuǎn)移矩陣進行各種異常模式的概率預(yù)測。

Markov狀態(tài)轉(zhuǎn)移矩陣如下：

(4)

其中，z為同一時刻最多可能的異常種類數(shù)；PIJ為t時刻處于I狀態(tài)，到t+1時刻變?yōu)镴狀態(tài)的概率，且滿足：

(5)

t時刻的狀態(tài)空間向量

Et=(e1,e2,…,ez)

(6)

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

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

Gt+1=EtM

(7)

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

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

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

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

(8)

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

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

yij=ai+μitij+σBiBij

(9)

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

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

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

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

(10)

j=2,3,…,mi

則

(11)

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

建立似然函數(shù)：

(12)

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

(13)

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

(14)

(15)

(16)

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

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

(17)

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

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

(18)

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

(19)

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

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

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

s(t)=tA-t-t-

(20)

式中，t-為當(dāng)前研究時段之前的總運行時間。

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

(21)

(22)

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

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

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

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

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

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

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

CftA≤Cds(t)

(23)

3 實例應(yīng)用

3.1 生產(chǎn)線機床加工過程控制及狀態(tài)預(yù)測

國內(nèi)某活塞生產(chǎn)線的主要加工流程見圖2。

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

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

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

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

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

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

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

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

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

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

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

表2 各機床發(fā)生異常的轉(zhuǎn)移關(guān)系Tab.2 Abnormity transmission relationship between machine tools

計算得到Markov狀態(tài)轉(zhuǎn)移矩陣：

M=

(24)

因收集不合格品數(shù)據(jù)前最近一次的異常發(fā)生在機床M7，故令該時間節(jié)點為t，則狀態(tài)空間向量可表示為

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

(25)

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

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

(26)

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

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

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

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

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

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

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

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

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

(27)

(28)

式中，μa、σa分別為a的均值和標(biāo)準差。

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

(29)

(30)

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

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

(31)

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

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

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

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

(32)

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

Cd/Cf≥4.85

(33)

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

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

t+t-≤1 079.66 h

(34)

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

4 結(jié)論

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

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

(3)基于“維修如舊”的假設(shè)，結(jié)合生產(chǎn)效益指導(dǎo)預(yù)防維修策略，且符合企業(yè)的實際利益需求。

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