考慮機器故障的作業(yè)車間調度方案魯棒測度方法

巴智勇，袁逸萍+，戴 毅，李曉娟，阿地蘭木·斯塔洪，劉金朵

(1.新疆大學 機械工程學院，新疆 烏魯木齊 830047；2.新疆生產力促進中心，新疆 烏魯木齊 830099)

0 引言

作業(yè)車間作為產品制造的直接執(zhí)行者，承載著大量的生產任務，生產任務間存在著復雜的關聯關系。隨著生產系統規(guī)模的擴大和復雜性的提高，生產系統運行中的不確定因素急劇增加，對于生產系統的影響難以估計，如果車間管理者無法提前或及時進行資源調度和計劃調整，極易導致生產進度的拖延和生產成本的上升，嚴重時引起生產過程混亂。

在實際生產過程中，機器故障是常見的不確定因素，它不僅會導致實際生產結果和計劃產生偏差，還可能引起物料配送等相關生產環(huán)節(jié)的混亂，對生產過程帶來重大影響。因此，在制定調度方案時考慮未來機器故障干擾，可以有效減少機器故障對調度性能的影響。如何有效評估故障環(huán)境下調度方案的性能成為魯棒性調度的一個難點。

近幾十年，大量學者對于生產過程中存在機器故障的魯棒調度問題進行了深入的研究。由于在生產之前無法得到實際生產數據，蒙特卡洛仿真(Monte Carlo simulation, MC)和替代測度(Surrogate Measure, SM)是常用的兩種魯棒性測度方法。

(1)基于蒙特卡羅仿真的測度方法

Zandieh等[1]通過免疫算法求解隨機機器故障下混合流水車間調度問題，通過故障模擬器評估可行調度方案的期望完工時間，將其作為選擇適應度值對種群進行選擇，引導算法生成期望時間較小的調度方案。張先超等[2]針對設備故障下多階段流水車間的調度問題，將生成的調度方案在故障環(huán)境下運行100次，將平均工期和工期標準差作為調度方案的性能和風險測度。Shen等[3]使用場景表示擾動的隨機性，進而提出兩種替代性測度方法，并通過仿真實驗證明其有效性。Amirl等[4]通過仿真軟件建立大量擾動場景，評估可行調度方案的性能，并通過對仿真數據進行訓練，構建調度決策神經網絡。Al-Hinai等[5]研究了隨機機器故障下柔性作業(yè)車間的魯棒調度問題，以工序完工期偏差為基礎提出3個魯棒性指標，通過模擬大量的故障場景來評估調度的魯棒性。顧澤平等[6]針對多種不確定因素環(huán)境下柔性作業(yè)車間的多目標調度優(yōu)化問題，提出使用離散仿真方法來評估調度方案的性能，引導混合遺傳算法的進化方向，并通過對比試驗證明了所得近似最優(yōu)解具有更高的魯棒性。

(2)替代測度方法

針對單機魯棒調度問題，Mehta等[7]詳細給出了5種替代測度，并進行了比較，提出了基于替代測度的魯棒調度算法。Goren等[8]根據機器加工時間和機器維修時間提出了兩種替代測度，并通過仿真驗證了替代測度對單機魯棒調度的有效性。陸志強等[9]針對考慮預防性維護的流水線魯棒調度問題，提出一種替代測度，并基于替代測度設計三階段啟發(fā)式算法對模型進行求解。Leon等[10]研究了作業(yè)車間魯棒性調度問題，以所有工序的平均總松弛(total slack)時間為替代測度，將該替代測度嵌入到遺傳算法中，引導算法生成魯棒性調度方案。Xiao等[11]結合調度方案和工序加工時間的分布信息，考慮關鍵路徑和非關鍵路徑對調度魯棒性劣化的影響，提出了兩種替代性的魯棒測度，并通過仿真實驗證明其有效性。Al-Fawzan等[12]用自由松弛(free slack)時間之和來測度資源受限項目調度的魯棒性。Xiong等[13]考慮機器負載與故障的相關性，用機器負載來表示工序總松弛時間的重要性，提出工序加權總松弛時間之和來近似調度方案的魯棒性。Wu等[14]結合工序松弛時間和故障信息提出一種調度魯棒替代測度方法，評估調度完工期的延遲風險。

以上成果都為本文的研究工作提供了借鑒和參考。蒙特卡洛仿真雖能有效評估調度魯棒性，但求解速度慢；同時多數替代性測度方法，由于未能充分利用故障信息，導致評估的準確性下降，若要有效利用故障信息需要解決以下兩個難點：①機器故障的隨機性和造成的影響在傳播過程中的不確定性；②多個機器故障在傳播過程中的疊加效應。

因此，本文結合調度方案結構、機器故障概率和維修時間提出一種基于期望影響效應的替代性測度方法，首先由機器故障分布函數計算工序加工過程中機器發(fā)生故障的概率，簡稱工序故障概率，同時得到工序加工過程中機器平均維修時間，簡稱工序期望維修時間，將機器故障映射到工序層面，分析單工序故障影響的傳播效應，進而評估多工序故障影響的綜合效應，獲得各工序的期望完工時間，最終求解出調度魯棒性。

1 機器故障下作業(yè)車間魯棒性調度問題描述與假設

為表述方便，表1列出了文中用到的符號及其說明。

表1 符號表示及其說明

續(xù)表1

1.1 機器故障下車間作業(yè)調度問題描述與假設

作業(yè)車間調度問題可描述為：n個工件在m臺機器上加工，每個工件有特定的加工工藝，使用機器的次序一定，加工時間為常數；工件在0時刻到達，所有機器在開工前均可用，每臺機器同一時間只能加工一個工件。在生產前調度計劃已經制定，且在加工過程中工序在機器上的加工順序不變。

本文的研究對象是可行的調度方案，各工件在機器上的加工時間和加工次序已確定。當發(fā)生故障時存在以下假設：

(1)機器故障只發(fā)生在機器加工期間；

(2)加工過程僅當機器發(fā)生故障時允許中斷，在維修后繼續(xù)加工，受到影響工序采用右移策略進行控制；

(3)維修操作不改變機器役齡。

1.2 隨機機器故障描述

隨機機器故障可以用故障發(fā)生概率與維修時間兩個參數進行描述，一般認為機器故障是服從一定的概率分布，機器的維修時間為定值。本文假設機器的可靠性服從Weibull分布。

機器在t時刻的可靠性為:

(1)

式中：λ(t)為故障率函數；β為形狀參數，當β>1時，隨著t的增加，機器故障概率逐漸增大；θ為尺寸參數。

1.3 調度方案魯棒性測度指標

調度的魯棒性指調度在不確定環(huán)境下保持原有狀態(tài)或性能的能力，通常分為性能魯棒性(Performance Robustness, PR)和調度穩(wěn)定魯棒性(Stability Robustness, SR)[15]。性能魯棒性是指初始調度的目標值與不確定因素擾動下實際調度目標值的接近程度。調度穩(wěn)定魯棒性指擾動因素影響下實際執(zhí)行的調度與初始調度的接近程度。

(1)性能魯棒性

用實際調度的σr與原調度σp期望最大完工時間偏差作為原調度σp的性能魯棒性測度，可表示為：

PR(σp)=E[|Cmax(σr)-Cmax(σp)|]。

(2)

由于原調度σp受到機器故障的影響，工序的完工時間不會減小，即Cmax(σr)-Cmax(σp)>0修改為：

PR(σp)=E[Cmax(σr)]-Cmax(σp)。

(3)

式中：原調度方案的Cmax(σp)已知，而E[Cmax(σr)]為實際調度σr中工序最大的期望完工時間。

(2)穩(wěn)定魯棒性

用實際調度σr與原調度σp各工序期望完工時間偏差之和表示穩(wěn)定魯棒性，可表示為：

(4)

同理，式(4)修改為：

(5)

式中，原調度方案中各工序的完工時間已知，E[Ci,j(σr)]表示實際調度σr中各工序的期望完工時間。

由以上分析可知，兩類調度魯棒性指標都轉化為工序期望完工時間的求解問題。但在生產完成之前，無法得到實際調度各工序完工時間，難以評價原調度方案的魯棒性。

2 基于故障影響傳播效應的魯棒性測度

2.1 作業(yè)車間工序關聯關系分析

工序在制造資源加工時序上存在次序性，在工藝上存在順序性，使得工序間存在復雜關聯關系。為描述工序間關聯關系，給出以下定義。

定義1關聯工序。由于機器共用和工藝的順序性，使得工序間存在時間約束關聯：先加工工序時間的減少或增加可能引起后加工工序時間的提前或延遲。

先加工的工序稱為前向關聯工序，后加工稱為后向關聯工序，若關聯工序加工的先后順序是緊鄰的，同一工件緊后加工的工序，記作工件維緊后工序，同一機器緊后加工的工序，記作機器維緊后工序。

以圖1中工序O4,2關聯工序分析為例，按工序加工次序可將O4,2關聯工序分為前向關聯工序集合{O3,1,O3,3,O4,3,O2,2}和后向關聯工序集合{O4,4,O4,1,O1,2}。前向關聯工序集合中任意工序的加工時間變化可能會導致O4,2的延遲，如前向關聯工序O3,1與O4,2既不是同一工件的工序也不在同一機器加工，但工序O3,1的加工時間變化會通過工序O3,3、O4,3間接影響工序O4,2的完工時間。

定義2影響傳播鏈。關聯工序及其之間工序組成的有序集合，集合中的工序至多存在一個緊前關聯工序和一個緊后關聯工序。如圖1中O3,1與O4,2之間的影響傳播鏈為{O3,1,O3,3,O4,3,O2,2}。

2.2 工序關聯下故障影響效應分析

機器故障導致工件的加工過程發(fā)生中斷，由于本文假設工件在機器修復后繼續(xù)加工，機器故障實質上增加了工件在機器上的停留時間[2]，在工序層面，可視作工序加工時間的增加。本文通過機器故障概率分布函數，將機器故障映射到工序層面，得到各工序的故障概率，如式(6)所示：

pri,j=1-e-[(ai,j/θ)β-(bi,j/θ)β]。

(6)

式中：ai,j為工序Oi,j開工時機器Mi的役齡；bi,j為工序Oi,j加工完成時機器Mi的役齡。

式(7)給出了工序Oi,j的期望維修時間：

(7)

2.3 魯棒性測度

作業(yè)車間調度問題中，各工序具有約束松弛度緊、約束內聯度高等特點[17]。當多個工序發(fā)生故障時，各工序故障將沿著約束鏈進行傳播并在該過程中出現疊加影響的效應。

如圖3所示，工序O1,1的前端關聯工序集合為{O2,1,O2,3,O1,2}；按工序間存在的時序約束，前向關聯工序可分為Path1和Path2兩條傳播路徑。

工序O1,1的完工時間期望延遲是兩條傳播路徑上工序造成延遲累計最大值為DT1,1=max(DT(1)1,1,DT(2)1,1)}。

根據以上分析，定義工序由前向關聯工序造成延遲時間

DTi,j=max(DT(1)i,j,…,DT(k)i,j)。

(8)

式中DT(k)i,j為第k個影響路徑上工序對Oi,j造成的累計期望延遲。

根據調度性能魯棒性的定義，將工序的期望完工時間帶入式(3),可得到調度的性能魯棒性為：

cmax(σp)。

(9)

同理可得調度穩(wěn)定魯棒性為:

(10)

3 算法設計

根據上述分析，本文設計一種基于故障影響效應的調度魯棒測度求解算法,符號說明如表2所示，具體算法如下：

輸入原調度σp(工件在機器上的加工順序，各工序加工時間、開工時間、完工時間),Weibull分布函數的形狀參數β，尺寸參數θ，故障維修時間tr。

步驟1確定未添加期望維修時間工序集合NDO和已添加的工序集FDO，并對以下集合、變量初始化FDO←?;NDO←TO;PAO←?;AO←?;STAO←?;ETAO←?;k=0;a=0。

步驟4若CO?AO,則AO←AO∪{CO},STAO←STAO∪{STN},ETAO←ETAO∪{ETN},a←a+1,否則，確定出當前工序CO在AO中編號a′,更新AO[a′]=CO;更新STAO[a′]=max{STN,STAO[a′]},ETAO[a′]=max{ETN,ETAO[a′]}。

步驟6從工序集ARO中刪除CO，轉步驟5。

步驟7遍歷受影響工序完工集合ETAO，計算受影響工序完工時間ETAO(Oij)與原調度方案工序完工時間Ci,j(δp)的差值并求和，獲得穩(wěn)定魯棒性指標；計算受影響工序完工時間ETAO(Oij)的最大值max(ETAO(Oi,j))與原調度方案Cmax(σp)的差值，獲得性能魯棒性。

表2 算法中用到的符號及說明

4 仿真驗證

4.1 仿真方案設計

為驗證本文測度方法的有效性，從典型的基準案例庫中選取25個作為實驗案例；Adams[18]設計的2個案例(abz5,abz9)、Fisher和Thompson[19]設計的2個案例(ft10、ft20)、Lawrence[20]設計的8個案例(la01、la06、la11、la16、la21、la26、la31、la35、la40)、Storer[21]設計的4個案例(swv01、swv06、swv11、swv16)、Vaccari[22]設計的8個案例(tai01、tai11、tai21、tai31、tai41、tai51、tai61、tai71)和Yamada[23]設計的案例yn1。

采用Windows 2008 Server，3.0 GHz CPU，4 G內存，Python3.5作為仿真語言。

使用Della Croce等[18]提出的遺傳算法生成25個案例的不考慮機器故障的調度方案，以最小化makespan為調度目標。設置種群規(guī)模為300，交叉率為0.7，變異概率為0.05，迭代次數300。終止條件滿足迭代次數。本文假設機器可靠性服從Weibull分布，其中形狀參數β=2，尺寸參數θ分為3類，分別為調度方案中機器最大負載的0.5倍、1倍和1.5倍;4類機器故障維修時間tr分別為10,20,30,60；共存在300個實例。求解各案例在隨機機器故障環(huán)境下調度方案σp的性能魯棒性PR(σp)、穩(wěn)定魯棒性SR(σp)和期望完工周期E[Cmax(δp)]。本文采用蒙特卡洛實驗模擬隨機機器故障場景，每個案例執(zhí)行5 000次，獲得期望完工期完工期偏差期望值(PRMC)和工序時間偏差期望總和(SRMC)和E[Cmax(σr)]作為真實值。

4.2 有效性分析

為檢驗測度方法的有效性，設計2個指標來檢驗測度方法得到結果的精度。PRD(PR(σp),PRDMC)表示PR(σp)與PRMC相對偏差，SRD(SR(σp),SRMC)表示SR(σp)與SRMC的相對偏差，相對偏差越小，說明所提方法的測度結果越精確。

(11)

(12)

表3是25個案例在θ,tr組合成的12種故障水平下的仿真結果，由于每種故障水平下存在25個案例，表3中數據均為25個案例在各故障水平下相應指標的平均值。由表3可以看出，PRD(PR(σp),PRMC)的值均小于0.21%，說明提出的方法對調度性能魯棒性的測度精度平均可達到99.79%以上，同樣的SRD(SR(σp),SRMC)的最大值僅為5.81%，說明提出的方法對調度穩(wěn)定魯棒性的測度精度平均可達到94.19%以上，可以有效地表征調度的穩(wěn)定魯棒性。由于本方法可以同時求得PR(σp)和SR(σp)，兩者求解時間相同，用T1表示；η在[0.58%～0.87%]之間，說明本文提出方法的運行時間遠遠小于蒙特卡羅仿真方法，證明了本文算法的高效性。

表3 有效性分析

圖5為不同故障水平下PRD(PR(σp),PRMC)和SRD(SR(σp),SRMC)分布的四分位圖。由圖5a可知，當維修時間tr相同時，隨著θ的增加，SRD(SR(σp),SRMC)略有下降，由式(6)可知，θ越大時，機器故障概率越低，說明測度方法在低機器故障率環(huán)境下的測度精度更高且波動范圍較小。當θ相同時，隨著維修時間tr的增大，SRD(SR(σp),SRMC)略有上升，說明維修時間的增加對穩(wěn)定魯棒性測度的精度有所下降。圖5b可得，PRD(PR(δp),PRMC)在各故障水平下保持極小的誤差，隨著維修時間tr增大，PRD(PR(δp),PRMC)的波動略有增大，但都保持在1%以內，說明本文方法對調度性能魯棒性的測度精度極高。

4.3 相關性分析

Leon等[10]認為當調度方案中存在更多總松弛時間時，可以更好地減少調度方案的延遲風險，從而提出了基于工序平均總松弛時間的魯棒性測度方法，如式(13)所示：

(13)

式中ti,j為工序間總松弛時間。

Al-Fawzan等[12]提出使用工序自由松弛時間總和的魯棒性測度方法，如式(14)所示：

(14)

式中tfi,j為工序間的自由松弛時間。

Xiong等[13]認為負載越大的機器發(fā)生故障的概率越大，該機器上的工序總松弛時間越重要，因此提出了工序總松弛時間加權和作為魯棒性的測度方法，如式(15)所示：

(15)

式中wi,j為Oi,j所在機器的負載，wtot為所有機器負載之和。

為了進一步說明提出測度方法的有效性，對3種替代方法和本文所提方法與實際調度魯棒性進行線性相關分析(如表4)，分析發(fā)現，判定系數R2越接近1，線性擬合效果越好。由于實際調度魯棒性不能提前獲得，使用蒙特卡洛仿真得到的PRMC和SRMC作為實際調度的性能魯棒性和穩(wěn)定魯棒性。

表4 相關性分析

從表4可知，SR(σp)的判定系數R2遠大于其他3種替代性方法，說明SR(σp)與實際調度穩(wěn)定魯棒性的線性相關性很強，同樣可得PR(σp)的判定系數R2均優(yōu)于其他三種替代性方法，且R2均大于99%,說明PR(σp)與實際調度性能魯棒性線性相關性極強。

圖6a表示不同的測度方法與SRMC的判定系數R2在不同故障水平下的變化情況,從圖6a中可以看出SR(σp)與SRMC判定系數R2都保持在90%以上，且隨著參數變化的幅度較小，保持較高的穩(wěn)定性。而另外3個替代性指標RM1、RM2、RM3與SRMC的相關性均在30%以下，表明三者與SRMC的相關性很差，不能替代表示調度方案的穩(wěn)定魯棒性。

圖6b表示對于不同測度方法與PRMC判定系數R2在不同故障水平下的變化情況，從圖6b中可得PR(σp)不隨著維修時間和機器故障概率變化，保持較高的穩(wěn)定性，同時R2均保持在99%以上，說明PR(σp)與PRMC幾乎是線性相關。而RM1和RM3與PRMC的判定系數R2保持在30%～60%之間，存在一定的線性相關，隨著參數變化存在一定的波動性。RM2與PRMC的判定系數R2在70%～90%之間，說明RM2與PRMC相關性較高，但隨著tr的增加相關性有所下降。

5 結束語

本文研究了設備故障環(huán)境下作業(yè)車間的調度魯棒性的測度問題，從工序關聯角度分析機器故障及其擴散效應對調度魯棒性的影響，提出了基于故障傳播效應的測度方法并設計求解算法，仿真結果顯示本文提出的測度方法對于調度性能魯棒性的測度精度可達到99.79%以上，對調度穩(wěn)定魯棒性的測度精度可達到94.19%以上，證明本文方法的測度有效性。通過相關性分析，得到本文方法在各種故障水平下與蒙特卡洛方法的結果均能保持較高的相關性，并與其他3種替代性測度方法比較，證明了本方法的優(yōu)越性。

本文提出的方法可為故障環(huán)境下魯棒性調度提供精確的測度方法,引導算法生成。在動態(tài)調度問題中，可通過評估機器故障對調度性能影響程度，準確選擇調度時機和重調度策略。未來將進一步研究集成預維護活動的調度方案的風險評估及預維護活動的插入策略，并對機器故障環(huán)境下作業(yè)車間重調度時機和策略進行研究。