基于中智數(shù)及ε-支配的加工時間不確定IPPS問題多目標(biāo)優(yōu)化研究

金亮亮 張超勇 GEORGE GERSHOM CHRISTOPHER 文笑雨

(1.紹興文理學(xué)院 機電學(xué)院，浙江 紹興 312000；2.華中科技大學(xué) 機械科學(xué)與工程學(xué)院，湖北 武漢 430074；3.紹興文理學(xué)院 土木工程學(xué)院，浙江 紹興 312000；4.鄭州輕工業(yè)大學(xué) 河南省機械裝備智能制造重點實驗室，河南 鄭州 450002)

0 引言

集成式工藝規(guī)劃與車間調(diào)度(IPPS)問題是在柔性作業(yè)車間調(diào)度問題[1,2]的基礎(chǔ)上引入序列柔性、加工柔性得到的；通過整合工藝規(guī)劃及車間調(diào)度消除資源沖突、提高設(shè)備利用率[3-7]，因此，更適合于多周期、小批量、個性化的生產(chǎn)制造模式.

IPPS問題的本質(zhì)在于工藝規(guī)劃和車間調(diào)度模塊的有機融合.目前，IPPS問題中兩大模塊常見的融合方法可分為兩類：(1)順序融合模式：先生成若干備選工藝規(guī)劃，在車間調(diào)度環(huán)節(jié)中根據(jù)優(yōu)化目標(biāo)選取最合理的工藝規(guī)劃[8],但由于此模式常常無法生成全部可行的工藝規(guī)劃，導(dǎo)致僅能得到局部最優(yōu)解，同時該模式忽略了工藝規(guī)劃和車間調(diào)度內(nèi)在的聯(lián)系；(2)集成融合模式：將工藝規(guī)劃和車間調(diào)度混為一體進行建模及優(yōu)化，避免了順序融合模式中的缺陷，真正實現(xiàn)了集成優(yōu)化.

IPPS問題的優(yōu)化方法大致可分為三類：(1)元啟發(fā)式方法：由于IPPS屬于NP-Hard問題，因此這類方法應(yīng)用較為廣泛，從單目標(biāo)到多目標(biāo)，從靜態(tài)到動態(tài)，從確定到不確定問題等都有涉及.韓國全南國立大學(xué)Kim等人使用共生進化算法求解了基于網(wǎng)絡(luò)圖的IPPS實例，Kim等人提出的測試實例也成為經(jīng)典的實例[9].華中科技大學(xué)李新宇和高亮等針對Kim實例進行了相關(guān)的研究并取得了一系列成果[7-8]；Lian等人基于Kim實例的特點，提出了基于網(wǎng)絡(luò)圖中OR節(jié)點的編碼方式，使用新穎的殖民競爭算法對Kim實例進行了優(yōu)化[10].巴黎等人研究了基于工件批量柔性劃分的IPPS問題，提出了新的編碼方法并使用粒子群算法進行了求解[6]；Lv和Qiao提出了新的個體初始化方法及相應(yīng)的遺傳操作，使用混合遺傳算法優(yōu)化了IPPS實例并獲得了較好的結(jié)果[11].近年來Zhang和Wong針對IPPS實例提出了對象編碼方法，在求解Kim實例時取得了較好的效果，改進了許多已有的結(jié)果[12].文獻[13]提出了一種基于變鄰域搜索的文化基因算法，對于Kim的實例進行測試，得到了許多實例的最優(yōu)解.此外，也有學(xué)者使用元啟發(fā)式算法求解了多目標(biāo)IPPS問題.例如Li等人受博弈論中的納什平衡啟發(fā)，用納什平衡處理各個優(yōu)化目標(biāo)間的關(guān)系，并提出了一種混合算法求解多目標(biāo)IPPS問題[14].文獻[15]針對最大完工時間、機床最大負(fù)荷、機床總負(fù)荷三個指標(biāo)使用多目標(biāo)文化基因算法進行優(yōu)化.(2)基于Agent的優(yōu)化方法.這類方法優(yōu)化后的結(jié)果與使用元啟發(fā)式算法得到的結(jié)果有一定的差距；Wong等使用基于Agent的談判方法求解IPPS問題并提出了工件Agent和機器Agent間的談判協(xié)議[16].(3)其他方法：例如基于啟發(fā)式規(guī)則的方法[17]及數(shù)學(xué)規(guī)劃方法，前者常常無法得到較優(yōu)的結(jié)果，后者無法在合理的時間內(nèi)得到最優(yōu)解.文獻[3]建立了一個IPPS問題的集成融合模型，并進行了求解.

目前，IPPS問題的研究已取得了相應(yīng)的進展和研究成果，但實際生產(chǎn)制造過程中某一工序的加工時間常常不固定，前續(xù)工序加工時間波動必將影響后續(xù)工序的加工，從而影響整個調(diào)度方案的實施，這是現(xiàn)有IPPS問題研究的一大挑戰(zhàn).因此，研究加工時間不確定的IPPS問題有現(xiàn)實的意義.相對于確定的IPPS問題，加工時間不確定的IPPS問題的研究尚未成熟.常見的優(yōu)化方法諸如機會約束規(guī)劃方法[18]無法應(yīng)用于大規(guī)模優(yōu)化問題.基于模糊理論的優(yōu)化方法是一種較好的優(yōu)化方法.Wang等人通過用模糊數(shù)表征不確定加工時間，提出了一種基于模糊數(shù)的人工蜂群算法求解柔性作業(yè)車間調(diào)度問題[19]；雷德明將不確定加工時間映射為三角模糊數(shù)，在遺傳算法框架下優(yōu)化了最大完工時間[20].最近，Jamrus等人將加工時間用三角模糊數(shù)模擬，使用粒子群算法對柔性車間調(diào)度問題進行優(yōu)化[21].筆者前期在這方面做了相關(guān)的研究[22]，特別是在中智數(shù)建模的基礎(chǔ)上使用線性加權(quán)的方法優(yōu)化了IPPS實例，但這類方法無法給出最優(yōu)Pareto前沿,且現(xiàn)有方法在表達決策者偏好方面存在一定不足.

綜上，基于模糊數(shù)的方法大多將加工時間假設(shè)為三角模糊數(shù)，但現(xiàn)實中常常無法準(zhǔn)確估計加工時間服從何種分布或某種規(guī)律，通常僅知曉大致的波動范圍，這為傳統(tǒng)基于模糊數(shù)優(yōu)化方法的使用帶來了困難.相反，使用中智數(shù)描述一個不確定量只需要知道極大和極小值即可，而不必求得具體的分布或規(guī)律，這為刻畫不確定加工時間帶來了方便.基于中智數(shù)的優(yōu)點，本文使用中智數(shù)對不確定加工時間進行模擬并研究相應(yīng)的優(yōu)化方法.此外，決策者常常僅對目標(biāo)空間某部分的非支配解感興趣而非全部非支配解，然而現(xiàn)有文獻中提出的方法大多無法利用偏好信息引導(dǎo)Pareto前沿趨向于決策者感興趣的區(qū)域.因此，本文利用中智數(shù)刻畫不確定加工時間，提出基于偏好的ε-支配優(yōu)化算法求解多目標(biāo)不確定IPPS問題，并對名義最大完工時間及魯棒性兩個指標(biāo)進行優(yōu)化.

1 基本概念與問題描述

1.1 中智數(shù)及其運算

中智數(shù)最早由Smarandache教授提出[23]，可表示為N=a+bI，其中a和b為兩個實數(shù)，I為不確定區(qū)間，滿足I2=I,0·I=0.例如：N=2.2+0.3I其中I=[-0.1,0.5]等價于N=[2.17,2.35]，也可表示為N=2.17+0.18I,I=[0,1].通過調(diào)整參數(shù)a和b可將兩中智數(shù)的不確定部分調(diào)整為同一個I.中智數(shù)的加減法運算可表示為N1±N2=a1±a2+(b1±b2)I.對兩個中智數(shù)Ni=ai+biI,Nj=aj+bjI，I=[β-,β+]，則Ni≥Nj的概率可表示為[24]：

(1)

對于n個中智數(shù)，兩兩之間通過式(1)計算后可得到概率矩陣P=(Pij)n×n，根據(jù)(2)式計算每個中智數(shù)排序順序值ri，則最大ri值對應(yīng)的中智數(shù)即為n個中智數(shù)中最大者.

(2)

中智數(shù)的“取大”算子可定義為：

1.2 問題描述與建模

IPPS問題可定義為：給定n個工件在m臺機器上加工，每個工件均含有若干道工序，通過選取合理的工序、機器、加工順序，得到最優(yōu)調(diào)度方案.如圖1所示，IPPS問題中可選的工序、加工資源及加工時間可用網(wǎng)絡(luò)圖進行描述[9].圖中開始節(jié)點和終止節(jié)點為虛擬節(jié)點，表示一個工件加工的開始與結(jié)束.每個普通節(jié)點代表一個工序，提供了該工序的ID、備選加工設(shè)備及相應(yīng)的加工時間.若某個分支旁標(biāo)有“OR”，則兩個分支中只需選擇一個分支上的工序即可.各工序只要滿足箭頭所指的先后關(guān)系(直接或間接)即可.相對于傳統(tǒng)柔性作業(yè)車間問題，IPPS問題中存在著三種柔性：工序柔性、序列柔性、加工柔性.工序柔性是指一個工序可在任意一臺備選的設(shè)備上加工；序列柔性是指同樣一組工序可以有若干種排列方式，只要符合網(wǎng)絡(luò)圖中的先后順序；加工柔性是指可以選擇不同的“OR”分支上的工序來完成工件加工.參考現(xiàn)有的IPPS問題的數(shù)學(xué)模型[15]，對中智數(shù)的確定部分和不確定部分分離，可以建立基于中智數(shù)的IPPS數(shù)學(xué)模型.

首先對模型中所需的相關(guān)符號、變量等進行定義.

n：工件集合

M：機器集合

Ti：工件i的所有工藝規(guī)劃集合

NOil：工件i的第j個工藝規(guī)劃所有工序集合

Rijl：工序Oijl的可用機器集合

A：一個很大的正整數(shù)

Xil：=1，若工件i的第l個工藝規(guī)劃選中；=0，否則

Zijkl：=1，若工序Oijl在機器k上加工；=0，否則

Yijli′j′l′：=1，若工序Oijl直接或間接地在工序Oi′j′l′之后加工；=0，否則

優(yōu)化目標(biāo)：

(3)

(4)

約束方程：

(5)

(6)

(7)

(8)

(9)

(10)

(11)

(12)

(13)

(14)

(15)

(16)

1.3 ε-支配與多目標(biāo)優(yōu)化

多目標(biāo)優(yōu)化算法中，常得到一組折中解(非支配解)，而不是全局最優(yōu)解，以便決策者從中選取較為合理的方案.具有代表性的算法為Deb等人提出的NSGA-II多目標(biāo)遺傳算法[25].然而，現(xiàn)有多目標(biāo)進化算法常將算法得到的非支配解全部呈現(xiàn)于決策者，而無法在進化中融入決策者的偏好信息，因此可以在算法中融入偏好信息，引導(dǎo)算法的搜索過程集中于感興趣的區(qū)域.本文通過目標(biāo)空間中的一控制點引入偏好信息，將目標(biāo)空間劃分為偏好區(qū)域和非偏好區(qū)域，提出基于偏好及中智數(shù)的ε-支配多目標(biāo)進化算法.

定義1：設(shè)p和q為群體中任意兩個個體，則當(dāng)fi(p)≤fi(q),?i且fj(p)

定義2：若一個個體p沒有被當(dāng)前所有個體支配，稱p為非支配解；群體中所有非支配解構(gòu)成Pareto前沿.

在ε-支配中先將目標(biāo)空間劃分為若干尺度為εj的格子，通過格子間的支配關(guān)系決定個體間的支配關(guān)系.通過控制格子的尺度進而控制非支配解的分布.

2 算法設(shè)計

IPPS問題為NP-Hard問題，整數(shù)線性規(guī)劃難以對其求解.本研究由于引入了中智數(shù)，導(dǎo)致模型更加復(fù)雜，求解更加困難.故本文提出一種ε-支配多目標(biāo)遺傳算法.以下介紹算法的編碼、選擇交叉等操作并給出該算法的流程圖.

2.1 編碼方案

如圖2所示，一個個體中包含調(diào)度串、工藝規(guī)劃串、工序串.調(diào)度串采用基于工序的編碼方法，數(shù)字i第j次出現(xiàn)對應(yīng)工序為工件i工序串中第j個位置的工序，調(diào)度串中含有的格子數(shù)恰好為實際一個調(diào)度問題中最大可能的工序數(shù).工件i的序號在調(diào)度串中正好出現(xiàn)|NOil|次，當(dāng)各工件實際工序數(shù)之和小于各工件最大可能的工序數(shù)之和，則在調(diào)度串多余位置用數(shù)字0補足.圖2a)中各工件分別有3,2,4個工序，因此，無須在調(diào)度串中補零.工序串反映了工藝規(guī)劃的信息，每一個格子代表一個工序，括號內(nèi)分別是工序ID及對應(yīng)的機床.例如，圖2a)中第一個工件三個工序加工順序為3,1,2且分別在機床2,1,4上加工.工藝規(guī)劃串僅含一位數(shù)字，代表該工件所選擇的工藝規(guī)劃(工序組合)代號.對于圖1中的工件，根據(jù)“OR”節(jié)點，分為兩個工藝規(guī)劃(工序組合)：工序1-11、14和工序1-8、12-14.

圖2 編碼方案與遺傳操作

2.2 解碼方法

傳統(tǒng)調(diào)度問題采用主動調(diào)度方法進行解碼并生成甘特圖，相對于半主動調(diào)度，基于主動調(diào)度的解碼方法能夠生成更優(yōu)的調(diào)度方案.但本研究中工序加工時長和最大完工時間均為中智數(shù)，無法直接使用現(xiàn)有的解碼方法.改進的解碼方法主要步驟如下：

(N3為機床k上空余時間段結(jié)束時間)時可將當(dāng)前工序安排到機器上當(dāng)前空余時間段.

(3)若無法在空余的時間段內(nèi)插入工序，只能將工序安排到相應(yīng)機床的當(dāng)前末尾的位置，工序開工時間為：

(4)按照個體調(diào)度串中的序號，重復(fù)上述步驟，逐個處理完所有工序.

2.3 遺傳操作

考慮到非支配解的數(shù)量多少以及個體進化的快慢，本文提出的算法部分參考了NSGA-II算法及現(xiàn)有的ε-支配多目標(biāo)算法，算法包括的遺傳操作主要有交叉和變異操作.圖2b)為交叉操作的過程，選取兩個父代個體P1、P2中被選中的工序串與工藝規(guī)劃串并復(fù)制一份分別放到兩個子代個體O2、O1對應(yīng)位置；將P1(P2)中未選中工件的工藝規(guī)劃串和工序串復(fù)制一份到O1(O2)相應(yīng)位置.由于同一工件不同工藝規(guī)劃包含不同數(shù)目的工序，需要調(diào)整兩個個體的調(diào)度串.將個體P1(P2)調(diào)度串中未被選中個體的ID按照原位置填入O1(O2)調(diào)度串對應(yīng)位置中，將個體P2(P1)調(diào)度串中被選中個體的ID依次填入O1(O2)調(diào)度串對應(yīng)空位中.所有工件的工序均填入調(diào)度串后空缺位置用數(shù)字0填充.圖2b)中由于P1、P2中第三個工件分別有3、4個工序，交換后個體O2的調(diào)度串中需要補充一個數(shù)字0.此外，為充分利用IPPS網(wǎng)絡(luò)圖中序列柔性，如圖2c)所示，設(shè)計了工藝規(guī)劃間的單點交叉操作.

變異操作包含三種：隨機交換調(diào)度串中兩個任意位置的數(shù)字；為某工序隨機選取不同的機床；為某工件重新生成一新的工藝規(guī)劃.算法運行變異操作時在三種變異中隨機選取一種.

2.4 選擇操作與算法流程

對于決策者給定的點P(坐標(biāo)為(Pf1,Pf2))，最小化問題中偏好區(qū)域即為坐標(biāo)原點與點P圍成的矩形，如圖3所示.除偏好區(qū)域外的其余目標(biāo)空間為非偏好區(qū).將兩部分區(qū)域劃分不同尺度的網(wǎng)格，ε-支配即為網(wǎng)格間的支配.顯然，偏好區(qū)網(wǎng)格較非偏好區(qū)更密，這表明偏好區(qū)域能夠容納更多的非支配解，而在非偏好區(qū)由于網(wǎng)格尺度較粗，僅能容納少量的非支配解，這樣引導(dǎo)算法的非支配解更多地分布于偏好區(qū).

圖3 ε-支配示意

隨機產(chǎn)生初始個體后，構(gòu)造普通Pareto支配的非支配解集合：將初始個體利用NSGA-II類似算法得到非支配個體(此時可能有不止一個個體處于同一格子中)，按照第一個指標(biāo)f1的數(shù)值由小到大排列各個非支配解，依次將排列后的非支配解按上述定義3的方法放入指定的格子中，從第二個非支配解開始，判斷若當(dāng)前處理的個體與前一個個體均處于偏好區(qū)或非偏好區(qū)則進一步按以下處理：

1)若當(dāng)前個體C(如圖3)在已有個體A下方，則刪除個體C上方所有個體；

2)若當(dāng)前個體B(如圖3)在已有個體B右方，則刪除當(dāng)前個體B；

3)若當(dāng)前個體E(如圖3)和已有個體D在同一格，則只保留距離該格子左下角頂點歐氏距離最短的點(圖中為E點).

否則接受當(dāng)前個體.

這樣，處理完所有非支配解后得到了初始?xì)w檔集.初始?xì)w檔集中的個體相互間均不構(gòu)成ε-支配關(guān)系.之后算法按照選擇、交叉、變異進行優(yōu)化，但傳統(tǒng)ε-支配多目標(biāo)優(yōu)化算法中每次僅進行單個個體間的交叉操作，效率與優(yōu)化效果較差.此外，為提高偏好區(qū)域中的個體被選入下一代的概率，對選擇操作進行了改進，具體如下：

1)任取兩個個體若均在非偏好區(qū)，則去掉非支配的個體，保留另一個個體；

2)若兩個個體均在偏好區(qū)，則同時接受這兩個個體；

3)若兩個個體所屬區(qū)域不同，確保處于偏好區(qū)的個體被選中，若另一個體支配處于偏好區(qū)的個體，則也被選中進入下一代；

4)重復(fù)以上步驟直到所選個體數(shù)目達到規(guī)定要求.

為改善個體多樣性，選擇操作中允許一小部分新生成的個體代替老的個體.算法每一輪進化后允許15%的新生成個體代替老的個體.

歸檔集更新主要考慮用當(dāng)前群體中的非支配解更新原歸檔集中的ε-非支配解.因此，先將當(dāng)前群體中非支配解和歸檔集中的個體合并，并按照第一個指標(biāo)f1的數(shù)值由小到大排列各個體，參考初始?xì)w檔集生成方法，得到當(dāng)前更新后的歸檔集.需要注意的是：若前后兩個個體均處于(非)偏好區(qū)，按ε-支配關(guān)系決定是否保留個體，若前后兩個個體處于不同區(qū)域的直接按普通的Pareto支配關(guān)系決定是否保留個體.

本文使用中智數(shù)來模擬現(xiàn)實中不確定的加工時長及完工時間，定義了兩種不同的工序加工時長波動范圍：1)小范圍波動：工序加工時長最大波動量為其名義值的±5%；2)較大范圍波動：工序加工時長最大波動量為其名義值的±15%.算法主要流程如圖4所示.

圖4 算法流程

本文主要考慮兩個優(yōu)化指標(biāo)：名義最大完工時間指標(biāo)及魯棒性指標(biāo).由于實際的最大完工時間是一個中智數(shù)，完工時間存在圍繞名義值上下波動的情況，故本文中魯棒性指標(biāo)定義為最大完工時間上下波動的最大偏差.顯然，完工時間波動的偏差越小則魯棒性越好.

3 計算實例及分析

3.1 參數(shù)設(shè)置

計算實例來自Kim等人提出的一組實例[9]，共有24個實例.工序加工時長精確到小數(shù)點后三位.名義最大完工時間和魯棒性指標(biāo)分別對應(yīng)目標(biāo)空間中的橫、縱坐標(biāo).在偏好區(qū)，橫、縱坐標(biāo)的偏差εj分別為0.005和0.004；在非偏好區(qū)，橫、縱坐標(biāo)的偏差分別為0.05和0.04.個體總數(shù)設(shè)置為800個，迭代次數(shù)800次，交叉和變異概率分別為0.7和0.08.所有計算在一臺裝有I5-9600(3.7GHz)CPU及16GB內(nèi)存的計算機上進行.

3.2 計算結(jié)果與分析

首先求解Kim實例中最復(fù)雜的第24個實例，該實例中含有18個工件，在15臺機床上加工.比較在工序加工時長處于較小波動范圍時得到的非支配解；圖5為相應(yīng)的Pareto前沿，本例中偏好區(qū)域為0≤f1≤650及0≤f2≤20圍成的矩形區(qū)域.由圖可知，初始非支配解經(jīng)過不斷進化形成了最終的Pareto前沿，并且大多數(shù)ε-非支配解集中于偏好區(qū)域(圖中方框內(nèi)區(qū)域為偏好區(qū)域).此外，非偏好區(qū)域內(nèi)的非支配解只有一個，這說明提出的ε-非支配多目標(biāo)算法確實能夠根據(jù)決策者偏好獲得較多的非支配解.圖6為相同條件下，偏好區(qū)域為0≤f1≤550,0≤f2≤35時的非支配解分布.顯然，對于不同的偏好區(qū)域，ε-非支配解分布的區(qū)域不同且隨偏好區(qū)域改變而改變.本例中ε-非支配解仍然集中于偏好區(qū)且只有兩個非支配解處于偏好區(qū)之外，因此，這表明提出的算法確實能夠使非支配解盡可能集中于偏好區(qū).

圖7為當(dāng)加工時長在較大范圍波動下第24個實例的Pareto前沿.對比前幾個圖可以發(fā)現(xiàn)當(dāng)加工時間的波動范圍增加時，得到的最大工期波動范圍也隨著增加，但對名義最大工期影響不是很大.此外，大多數(shù)非支配解均分布于給定的偏好區(qū)域(0≤f1≤600,0≤f2≤45).

圖5第24個實例小波動范圍非支配解分布

圖6 第24個實例小波動范圍非支配解分布(不同偏好區(qū)域)

圖7第24個實例較大波動范圍非支配解分布

圖8為圖7中偏好區(qū)域內(nèi)某非支配解對應(yīng)的甘特圖(名義加工時間).圖9為第12個實例優(yōu)化前后的非支配解分布(加工時間較大范圍波動下)，優(yōu)化后非支配解數(shù)量增加且主要分布于偏好區(qū)內(nèi).

選取Kim實例中第1、12、24三個實例(此三個實例分別代表典型的小、中、大規(guī)模實例)，得到每個實例5次計算時長的平均值并將該平均值與Kim實例相應(yīng)實例的計算時長進行比較，相關(guān)結(jié)果列于表1中.

圖8 第24個實例的甘特圖

圖9第12個實例的非支配解分布

表1 計算時長對比

由表1可知，求解大規(guī)模實例時，本文提出的算法所需計算時間較短.說明本文算法在求解這類實例時具有較高的效率，能快速找到合理的非支配解.同時，對于中小規(guī)模的實例，文獻[9]的算法計算較快，這是由于本文使用了中智數(shù)模擬不確定加工時間，在解碼等過程中中智數(shù)的運算增加了計算量.此外，為得到偏好區(qū)與非偏好區(qū)內(nèi)的非支配解，特別是保留并增加偏好區(qū)中的非支配解，本文提出了較為復(fù)雜的非支配排序方法，這也使得計算過程變得相對復(fù)雜，因此增加了計算時間.

4 結(jié)語

工序加工時間的不確定性給調(diào)度的魯棒性帶來了挑戰(zhàn)，鑒于中智數(shù)在描述不確定信息方面具有優(yōu)勢，本文將其引入集成式工藝規(guī)劃與調(diào)度問題的優(yōu)化中，建立了基于中智數(shù)的數(shù)學(xué)模型，提出了基于ε-支配的多目標(biāo)優(yōu)化算法.通過控制網(wǎng)格尺度及參考點的設(shè)定控制非支配解的分布.算法對名義最大完工時間及魯棒性兩個指標(biāo)進行了優(yōu)化，計算結(jié)果證明了算法的有效性，大多數(shù)非支配解集中于偏好區(qū)域.