作業車間調度的多工序精確聯動鄰域結構混合進化算法

巴智勇,袁逸萍,裴國慶,王 波

(1.新疆大學 機械工程學院,新疆 烏魯木齊 830017;2.卓郎新疆智能機械有限公司,新疆 烏魯木齊 830057)

0 引言

作業車間調度問題(Job Shop Scheduling Problem, JSP)不僅是典型的NP-hard問題,還是管理領域中許多優化問題的通用模型,如項目計劃調度、岸橋調度優化[3]、軌道列車時刻表調度優化[4]等。另外,JSP模型能以一種明確的方式對問題添加各種側約束和目標,使其應用范圍不限于調度優化問題[5]。

學者們針對JSP提出了大量高效的求解方法,如數學規劃[6]、調度規則[7]、啟發式和元啟發式算法。由于JSP的NP-hard特性,精確求解方法只適合求解小規模問題,難以應用于實際生產,當前主流研究方向是采用近似方法求解JSP。近似求解方法包括大量經典和新型的元啟發式方法,如粒子群優化[8]、遺傳算法[9]、差分進化[10]、蟻群算法[11]、細菌覓食算法[12]、教與學優化算法[13]和灰狼優化算法[14]等。然而,受限于基本理論和算法框架,單一的元啟發式算法在JSP上的求解效果并不理想,越來越多學者利用兩個或多個算法的互補特征開發更為先進的混合算法。REN等[15]提出一種混合局部搜索策略的遺傳算法求解JSP,通過引入多樣化的種群更新機制防止算法過早收斂;REHAB等[16]設計了基于鄰域多樣化強化策略的粒子群優化算法,實現了對JSP解空間的廣泛搜索;PENG等[17]提出基于路徑重連的禁忌算法,采用路徑重連方法產生更多樣的候選解,同時在禁忌搜索中應用更復雜的N7鄰域,改進了多個JSP基準算例的上界;CHENG等[5]提出一種基于機器最長公共子序列的新型交叉算子,并設計基于質量與相似性的種群替代策略來保持種群多樣性,更新了多個JSP基準算例的已知最優解;LIANG等[18]提出一種基于新型鄰域結構的自適應遺傳算法求解JSP,該算法的交叉概率和變異概率可以根據種群適應度的離散程度進行非線性和自適應調整,實驗結果表明算法在求解精度和收斂效率方面均有顯著提升;ALKHATEEB等[19]設計了一種集成模擬退火搜索的混合離散布谷鳥算法求解JSP,該算法首先采用基于對抗學習的方法生成多樣性的初始解,通過集成變鄰域搜索與萊維飛行方法實現廣泛的局部搜索,并設計了基于精英對抗的學習方法跳出局部最優解;SHAHED等[20]提出性能驅動的元啟發式切換方法求解JSP,該方法將多算子差分進化和粒子群優化集成在單一算法框架中,并引入性能驅動的切換機制將搜索性能不佳的算法切換為其他可能的算法。

大多數先進的JSP混合算法通過局部搜索組件來強化對有潛力解空間的深度搜索,而鄰域結構作為局部搜索的核心,極大地影響著算法的效率和質量。JSP局部搜索中較為常用的是基于關鍵路徑的工序塊鄰域結構,通過交換關鍵工序塊中任意兩工序可以形成新解。當前研究者大多從減少鄰域結構的規模和保證可行鄰域移動兩方面構建有效的鄰域結構。VAN LAARHOVEN等[25]交換同一臺機器上相鄰關鍵工序的鄰域結構,記作N1,然而N1鄰域結構不但規模大,而且包含大量非改進的移動操作;GRABOWSKI等[26]首先提出塊的概念,并將鄰域結構定義為把一個內部關鍵工序插入其塊首之前或塊尾之后的一種機制。此后,NOWICKI等[27]、BALAS等[28]和ZHANG等[29]提出更具競爭力的N4,N5,N6,N7鄰域結構;趙詩奎等[30-32]通過分析機器上空閑時間的利用機理,擴大了關鍵工序的移動范圍,并在此基礎上進一步提出多工序聯動鄰域結構和近似評價方法。

本文結合JSP的特性,設計了一種3對工序精確聯動的鄰域結構。針對進化算法中存在的種群早熟現象,設計了基于鄰域懲罰的交叉父本配對選擇算子與基于動態懲罰閾值的種群更新策略,并將基于新型鄰域結構的禁忌搜索嵌入進化算法框架,設計了兼具全局與局部搜索能力的JSP求解算法。最后,將本文所提算法與其他先進算法在基準算例上進行對比分析,驗證了所提算法的有效性。

1 問題描述

JSP問題可描述為:n個工件J={J1,J2,…,Jn}在m臺機器M={M1,M2,…,Mn}上加工,每個工件Ji由m個連續工序{Oi,1,Oi,2,…,Oi,m}組成,相應的加工機器集合為{σi,1,σi,2,…,σi,m},工序Oi,j在機器σi,j上的加工時間為pti,j,以最小化最大完工時間(Cmax)為優化目標。

對JSP問題作以下假設:

(1)所有工件零時刻達到,所有機器零時刻可用。

(2)同一工件同一時刻最多只能在一臺機器上加工,一臺機器同一時刻最多只能加工一個工件。

(3)工序在機器上一旦開始加工就不能中斷。

(4)工序之間的運輸時間和生產準備時間都包含在加工時間中。

以最小化最大完工時間為優化目標建立JSP混合整數規劃模型,模型構建過程中涉及的參數符號及說明如下:

sti,j為工序Oi,j的開工時間;

eti,j為工序Oi,j的完工時間;

JSP混合整數規劃模型為:

min(Cmax)=min(max(eti,j)),?Oi,j;

(1)

sti,j≥sti,j-1+pti,j-1,?{Oi,j,Oi,j-1};

(2)

eti,j=sti,j+pti,j,?Oi,j;

(3)

(sti,j-etr,s)yr,s,i,j(1-yi,j,r,s)+(str,s-eti,j)·

yi,j,r,s(1-yr,s,i,j)≥0,?{Oi,j,Or,s};

(4)

yi,j,r,s+yr,s,i,j=0;

(5)

sti,j≥0。

(6)

其中:式(1)為目標函數;式(2)為同一工件的緊鄰工序的順序約束;式(3)為計算工序的完工時間;式(4)與式(5)定義了同一機器上緊鄰兩道工序的加工時間不能重疊,即同一臺機器上后加工工序的開始時間不能早于前一道工序的完工時間;式(6)定義了工序開工時間的非負性。

2 多工序精確聯動的鄰域結構

鄰域結構的本質在于如何引導關鍵工序更充分地利用機器空閑時間,作為最有效的鄰域結構之一,N7鄰域結構雖然通過在較大范圍內移動單個關鍵工序對原方案產生較強擾動來生成鄰域解[33],但因缺乏足夠的理論分析,可能會產生大量無效的移動操作。為此,本文對N7鄰域結構進行理論分析,提出關鍵工序無效移動的條件并進行證明,進而對N7鄰域結構進行擴展優化。相關符號說明如表1所示。

表1 符號表示及其說明

2.1 N7鄰域無效移動的理論分析

N7鄰域結構存在4類關鍵工序移動操作,如圖1所示,當工序移動后,采用近似求解[34]方法重新計算受影響工序新的頭尾長度,則鄰域解的近似值C′max為受影響工序的頭尾長度之和的最大值,若C′max大于原方案的最大完工時間Cmax,則可證明該工序移動是無效的。

定理1關鍵塊中存在兩個工序u和v,分別為塊首與塊中工序,對于u與MS[v]之間的任意工序wi,若rJP[wi]+dJP[wi]+du≥rwi,則工序u移動到工序v之后不會減少最大完工時間。

證明如圖2所示,將工序u移動至v后,得到受影響工序{w1,…,wi,…,wn,v,u},采用近似評估方法重新計算各受影響工序的頭尾長度和C′max,計算如下:

(7)

C′max=max{r′w1+q′w1,…,r′wi+q′wi,…,
r′wn+q′wn,r′v+q′v,r′u+q′u}。

(8)

移動前,從工序wi角度,

Cmax=rwi+dwi+…+dwn+dv+qMS[v]。

(9)

移動后,工序MP[u]與v之間的工序將直接或間接利用塊前原有機器空閑時間以及u移動后產生的新的機器空閑時間,從工序wi的角度,近似最大完工時間C′max(wi)=r′wi+q′wi,由式(7)遞推可得:

r′wi≥rJP[wi]+dJP[wi];

(10)

q′wi≥qMS[v]+du+dv+dwn+…+dwi。

(11)

整理式(10)和式(11)可得

C′max(wi)≥rJP[wi]+dJP[wi]+du+dv+dwi+…
+dwn+qMS[v]。

(12)

若rJP[wi]+dJP[wi]+du≥rwi成立,則必有

C′max(wi)≥

rwi+dwi+…+dwn+dv+qMS[v]。

(13)

與式(9)相比,有C′max(wi)≥Cmax,則C′max≥Cmax必成立。證畢。

定理2關鍵塊中存在兩個工序u和v,分別為塊中與塊尾工序,若qJS[u]≥qJS[v],則工序u移動到v之后不會減少最大完工時間。

證明如圖3所示,移動后,從工序u的角度,近似最大完工時間C′max(u)=r′u+q′u,由式(7)遞推可得:

r′u≥rMP[u]+dMP[u]+dw1+…+dwn+dv;

(14)

q′u≥qJS[u]+du。

(15)

整理式(14)和式(15)可得C′max(u)≥rMP[u]+dMP[u]+du+dw1+…+dwn+dv+qJS[u],若qJS[u]≥qJS[v]成立,則

C′max(u)≥rMP[u]+dMP[u]+du+dw1+…
+dwn+dv+qJS[v]=Cmax,

(16)

C′max≥Cmax必成立。證畢。

定理3關鍵塊中存在兩個工序u和v,為塊首與塊尾工序,若rJP[wi]+dJP[wi]+du≥rwi且qJS[v]≤qJS[u],則工序u移動到v之后不會減少最大完工時間。

證明如圖4所示,移動后,從工序u的角度,近似最大完工時間C′max(u)=r′u+q′u,由式(7)遞推可得

r′u≥rJP[wi]+dJP[wi]+dwi+…+dwn+dv。

(17)

當rJP[wi]+dJP[wi]+du≥rwi成立時,由式(17)可得

r′u≥rwi+(dwi+…+dwn+dv)-du

=rv+dv-du。

(18)

由式(7)推得q′u≥qJS[u]+du必成立,結合式(18)可得C′max(u)≥rv+dv+qJS[u],若qJS[v]≤qJS[u],則C′max(u)≥rv+dv+qJS[v]>Cmax,有C′max≥Cmax必成立。證畢。

定理4關鍵塊中存在兩個工序u和v,分別為塊中與塊尾工序,對于u與MS[v]之間的任意工序wi,若qJS[wi]+dv≥qwi-dwi,則工序v移動到u之前不會減少最大完工時間。

證明如圖5所示,將工序v移動到u之前,得到受影響工序{v,u,w1,…,wi,…,wn},各自近似頭尾長度計算如下:

(19)

從工序wi的角度來看,移動后的最大完工時間C′max(wi)=r′wi+q′wi。由式(19)推導得到工序wi的近似頭尾長度:

r′wi≥qMP[u]+dMP[u]+dv+du+dw1+…
+dwi-1=rwi+dv;

(20)

q′wi≥qJS[wi]+dwi。

(21)

結合式(20)和式(21)可得C′max(wi)≥rwi+dv+qJS[wi]+dwi,若qJS[wi]+dv≥qwi-dwi成立,則C′max(wi)≥rwi+(qwi-dwi)+dwi=rwi+qwi=Cmax必成立,進一步可得C′max≥Cmax。證畢。

定理5關鍵塊中存在兩個工序u和v,分別為塊首與塊中工序,若rJP[v]+dJP[v]≥ru,則將工序v移動到工序u之前不會減少最大完工時間。

證明如圖6所示,從工序v的角度,移動后的近似最大完工時間C′max(v)=r′v+q′v,由式(19)推導可得:

r′v≥rJP[v]+dJP[v];

(22)

q′v≥qMS[v]+dwn+…+dw1+du+dv。

(23)

當rJP[v]+dJP[v]≥ru成立時,結合式(22)和式(23)可得C′max(v)≥ru+du+dw1+…+dwn+dv+qMS[v]=Cmax,進一步可得C′max≥Cmax。證畢。

定理6關鍵塊中存在兩個工序u和v,分別為塊首與塊尾工序,對于MP[u]與v之間的任意工序wi,如果qJS[wi]+dv≥qwi-dwi且rJP[v]+dJP[v]≥ru,則將工序v移動至工序u之前不會減少Cmax。

證明如圖7所示,移動后,從工序v的角度,近似最大完工時間C′max(v)=r′v+q′v,由式(19)推導得

q′v≥qJS[wi]+(dwi+…+dw1+du)+dv。

(24)

若qJS[wi]+dv≥qwi-dwi成立,則q′v≥qwi+(dwi-1+…+dw1+du)=qu必成立。由式(19)推導得r′v≥rJP[v]+dJP[v],若rJP[v]+dJP[v]≥ru成立,則r′v≥ru必成立。由以上結果可得C′max(v)≥ru+qu=Cmax必成立,進一步可得C′max≥Cmax。證畢。

可見,N7鄰域結構中參與移動的工序數量多且移動范圍廣,能從多個角度利用機器空閑時間,但也存在不能減少Cmax的無效移動操作。以上基礎理論分析為識別N7鄰域結構中的無效移動操作提供了理論支撐,同時也為更加復雜的多工序聯動鄰域結構設計指明了方向。

2.2 多工序精確聯動鄰域

本文將文獻[32]提出的多工序聯動的觸發條件與N7鄰域無效移動判定條件結合,構建更為強化的鄰域結構PN7+2MT,具體實現如下:

情形1塊首工序u移動至塊尾工序v后,聯動工序的移動方式如圖8所示,具體步驟如下:

(1)令msu=MS[u],如果JP[msu]≠?,且滿足rJP[msu]+dJP[msu]>rMP[u]+dMP[u],則按照工序的最早完工時間查找工序msu的某一道工件前序工序msu′,如果存在工序msu′滿足條件MP[msu′]≠?,且rMP[msu′]+dMP[msu′]=r′msu,則交換工序MP[msu′]和msu′。

(2)如果JS[u]≠?,且滿足qJS[u]>qMS[v],則按照工序的最晚開工時間查找工序u的工件維某一道后序工序u′,如果存在工序u′滿足條件MS[u′]≠?,且q′u-d′u=qMS[u′],則交換工序u′和MS[u′]。

情形2塊首工序u移動到塊中工序v之后,聯動工序的移動方式如圖9所示,具體步驟如下:

(1)令msu=MS[u],如果JP[msu]≠?,且滿足rJP[msu]+dJP[msu]>rMP[u]+dMP[u],則按照工序的最早完工時間查找工序msu的某一道工件前序工序msu′,如果存在工序msu′滿足條件MP[msu′]≠?,且rMP[msu′]+dMP[msu′]=r′msu,則交換工序MP[msu′]和msu′。

(2)若MS[u]與MS[v]之間的任意工序w均滿足rJP[w]+dJP[w]+duqMS[v]均成立,則按照工序的最晚開工時間查找工序u的某一道工件后序工序u′,如果存在工序u′滿足條件MS[u′]≠?,且qu′-du′=qMS[u′],則交換工序u′和MS[u′]。

情形3塊中工序u后移到塊尾工序v之后,聯動工序的移動方式如圖10所示,具體步驟如下:

(1)令msu=MS[u],如果JP[msu]≠?,且滿足rJP[msu]+dJP[msu]>ru,則按照工序的最早完工時間查找工序msu的某一道工件前序工序msu′,如果工序msu′滿足條件MP[msu′]≠?,且rMP[msu′]+dMP[msu′]=r′msu,則交換工序MP[msu′]和msu′。

(2)如果JS[u]≠?,且滿足qJS[u]>qMS[v],則按照工序的最晚開工時間查找工序u的某一道工件的后序工序u′,如果存在工序u′滿足條件MS[u′]≠?,且qu′-du′=qMS[u′],則交換工序u′和MS[u′]。

情形4塊尾工序v移動到塊首工序u之前,聯動工序的移動方式如圖11所示,具體步驟如下:

(1)令mpv=MP[v],如果JS[u]≠?,且滿足qJS[mpv]>qMS[v],則按照工序的最晚開工時間查找工序mpv的某一道工件后序工序mpv′,如果存在工序mpv′滿足條件MS[mpv′]≠?,且qmpv′-dmpv=qMS[mpv′],則交換工序mpv′和MS[mpv′]。

(2)如果JP[v]存在,且滿足rJP[v]+dJP[v]>rMP[u]+dMP[u],則按照工序的最早完工時間查找工序v的某一道工件前序工序v′,如果存在工序v′滿足條件MP[v′]≠?,且rMP[v′]+dMP[v′]=r′v,則交換工序MP[v′]和v′。

情形5塊尾工序v移動到塊中工序u之前,聯動工序的移動方式如圖12所示,具體步驟如下:

(1)令mpv=MP[v],如果JS[u]≠?,且滿足qJS[mpv]>qMS[v],則按照工序的最晚開工時間查找工序mpv的某一道工件的后序工序mpv′,如果存在工序mpv′滿足條件MS[mpv′]≠?,且qmpv′-dmpv=qMS[mpv′],則交換工序mpv′和MS[mpv′]。

(2)對于MP[u]與v之間的任意工序w滿足qJS[w]+dvrMP[u]+dMP[u]同時成立,則按照工序的最早完工時間查找工序v的某一道工件的前序工序v′,如果存在工序v′滿足條件MP[v′]≠?,且rMP[v′]+dMP[v′]=r′v,則交換工序MP[v′]和v′。

情景6塊中工序v移動到塊首u之前,聯動工序的移動方式如圖13所示,具體步驟如下:

(1)令mpv=MP[v],如果JS[u]≠?,且滿足qJS[mpv]>qMS[v],則按照工序的最晚開工時間查找工序mpv的某一道工件的后序工序mpv′,如果存在工序mpv′滿足條件MS[mpv′]≠?,且qmpv′-dmpv=qMS[mpv′],則交換工序mpv′和MS[mpv′]。

(2)如果JP[v]存在,且滿足rJP[v]+dJP[v]>rMP[u]+dMP[u],則按照工序的最早完工時間查找工序v的某一道工件的前序工序v′,如果存在工序v′滿足條件MP[v′]≠?,且rMP[v′]+dMP[v′]=r′v,則交換工序MP[v′]和v′。

2.3 鄰域結構復雜度分析

對于n個工件m臺機器的JSP,根據N7鄰域結構的4類移動方式,當工序塊有n個工序時,總的交換次數最多,為(n-1)×4×m,其計算復雜度為O(nm)。對于PN7+2MT,在N7鄰域結構的基礎上查找另外兩對交換工序,兩次交換的最大查找次數為2×(m-1),同時因為塊首工序后移與塊尾工序前移產生的聯動工序交換只需計算一次,所以PN7+2MT的最多查找次數為(n-1)×2×m×2×(m-1)+(n-1)×2×m+2×m×2×(m-1),計算復雜度為O(nm2)。

3 算法設計

3.1 算法框架

進化算法不但具有良好的理論研究基礎,而且在工程領域應用廣泛。然而進化算法偏重全局搜索,需要與基于問題特征的局部搜索集成才能獲得高質量的解。禁忌搜索算法側重于局部搜索,可以很好地結合JSP特征信息的鄰域結構,引導算法進行有效搜索。為了將全局搜索和局部搜索相結合,本文提出一種基于多工序精確聯動鄰域結構的混合進化算法(Hybrid Evolutionary Algorithm with Multi-operation precise joint movement Neighborhood Structure,HEA-MNS)。在進化算法部分引入基于新型個體距離的配對選擇算子,防止近親繁殖;然后,采用基于動態懲罰閾值的種群更新機制構建下一代種群,通過動態調整懲罰閾值實現集中性與多樣性搜索的動態平衡。禁忌搜索采用多工序聯動鄰域結構PN7+2MT對大鄰域空間進行精確搜索,提高了HEA-MNS的局部搜索能力。

算法1為HEA-MNS的基本框架,算法首先從一個初始隨機種群開始,采用禁忌搜索優化每個染色體(4行～6行);然后,采用配對選擇算子選擇交叉父本,通過交叉算子生成新的子代染色體(8行～9行);再次,采用禁忌搜索對新生成的子代染色體進行優化(10行～11行);最后,采用種群更新算子從原種群與子代染色體的集合中選擇下一代種群(12行)。下面詳細介紹HEA-MNS的主要組件。

算法1HEA-MNS偽代碼。

1 輸入:J,M,D。

2 輸出:算法終止時的最優調度方案s*與完工期Cmax。

3 初始化種群P={s1,…,sp}

4 For i={1,…,p} do

5 si←Tabu_Search(si)

6 end for

7 While終止條件未滿足do

8 (sj,sk)= Select_Parents(P)

9 (sp+1,sp+2)←Crossover(sj,sk)

10 sp+1← Tabu_Search(sp+1)

11 sp+2← Tabu_Search(sp+2)

12 {s1,…,sp}←Population_Updating(s1,…,sp,sp+1,sp+2)

13 end while

14 s*=argmim{f(si)|i=1,…,p}

15 Return s*

3.2 編碼與解碼

本文采用文獻[35]的編碼方法,即一條染色體由m個基因串組成,每個基因串表示各工件在機器上的加工序列,基因串的工序基因用相應的工件序號表示,因此染色體的長度為n×m。

對于JSP,以最大完工時間為優化目標的最優調度必在主動調度集合中[36],因此本文選擇GONG等[37]提出的插入式貪婪主動解碼算法。

3.3 禁忌搜索

作為車間調度問題最有效的求解方法之一,禁忌搜索能夠有效引導算法深入搜索有潛力的區域,為了加強局部搜索的擾動強度和有效性,本文采用第2章提出的精確多工序聯動鄰域結構PN7+2MT生成鄰域解。

(1)移動評估 當一組可行鄰域解被構造出來后,每個鄰域解都需要被重新計算或評估其最大完工時間,本文采用趙詩奎等[32]提出的快速評估方法對鄰域解進行近似計算。

(2)禁忌表 禁忌表作為禁忌搜索的短期記憶,可以有效避免搜索過程出現死循環。算法采用與TSSA(tabu search and simulatecl annealing)算法[38]相同的禁忌表結構。具體來說,如果一次移動包含交換工序u和v,則u和v之間的所有工序及其位置都存儲在禁忌表中。禁忌表長度根據決策變量規模動態調整[38],即在兩個給定的極值[Lmin,Lmax]之間隨機選取。

(3)鄰域解選擇 在每次迭代中,禁忌搜索過程首先考慮不在禁忌表中的移動,并從中選擇產生最優目標值的鄰域解。若兩種或兩種以上的移動都具有相同的最優目標值,則從中隨機選擇一個移動;如果一個移動能夠帶來比當前最優解更好的解,且該移動在禁忌表中,則對該移動進行特赦。

(4)終止準則 當找到最優解或禁忌搜索的連續未改善迭代次數超過給定極限值時,禁忌搜索過程停止。

3.4 基于鄰域懲罰的配對選擇與交叉操作

在進化算法中,配對選擇操作的主要任務是從父代種群中選擇用于生成子代個體的配對父本,其選擇標準在很大程度上影響算法的收斂性與多樣性。為了增強算法產生多樣化后代的能力,本文提出基于鄰域懲罰的配對選擇操作,其基本思想為:算法隨機從父代種群中選取一個染色體,另一個染色體從其懲罰鄰域外隨機選擇,如果懲罰鄰域外無其他染色體,則從懲罰鄰域內選擇距離最遠的一個染色體。為方便后續論述,本文給出以下定義:

定義1兩個解之間的距離。對于給定的兩個解x=(x1,…,xi,…,xm)和y=(y1,…,yi…,ym),解x與y的距離為各機器上工序排列的偏差距離的總和,記作Dpos(x,y),

Dpos(x,y)=

(25)

式中:n和m分別為工件數量與機器數量;xi和yi分別為解x與y的第i機器上的工序序列;pos(a,b[])為a在b序列中的位置。

定義2解與種群的距離。對于給定的種群P={s1,…,sp}與一個解I,I與種群P的距離為I與種群P中所有個體距離的最小值,記為DCI(I,P),

DCI(I,P)=min{Dpos(I,sj)|sj∈P}。

(26)

定義3種群的密度。對于給定的種群P={s1,…,sp},種群P的密度為種群的各個解si相對于P{si}距離的平均值,記作DCIav(P),

(27)

算法2中給出了交叉算子的偽代碼,對于當前種群P={s1,…,sp},采用以下規則選取交叉父本{sj,sk}。首先,根據當前種群的密度設定父本選擇的懲罰距離D,并初始化候選集合Candidate與懲罰集Penalized(3行),從當前種群Candidate中隨機選擇一個染色體作為第1個父本sj,將其從Candidate刪除(4行～5行);然后,通過不斷迭代,將與父本sj距離小于D的染色體從Candidate刪除,并添加至Penalized中(6行～11行);最后,若Candidate不為空,則從Candidate中隨機選擇一個染色體作為第2個父本,否則從Penalized中選擇距離sj最遠的染色體作為第2個父本(12行～16行)。

算法2基于動態閾值的父本選擇算子。

1 輸入:種群P={s1,…,sp}。

2 輸出:{sj,sk}。

4 從Candidate中隨機選擇一個染色體sj

5 Candidate=Candidate{sj}

6 foreach I∈P{sj} do

7 if Dpos(I,sj)

8 Candidate=Candidate{I},

9 Penalized=Penalized∪{I}

10 end if

11 end foreach

12 if Candidate≠? then

13 從Candidate中隨機選擇一個個體sk

14 else

15 sk= argmax{Dpos(si,sj)|si∈Penalized}

16 end if

17 return{sj,sk}

作為混合進化算法的重要組件,交叉算子不僅要生成多樣化的子代,還要將父輩的精英成分遺傳給子代,因此本文不再采用單一的交叉方式,而是將基于最長相同子序列的交叉算子[5]與基于機器工序串的交叉算子[15]混合使用。每次交叉均隨機選擇一種交叉方法進行染色體重組。

3.5 種群更新算子

為進一步增強種群的多樣性,設計了一種基于動態懲罰閾值的種群更新策略。該策略的主要設計原則是在算法的初始階段通過設定較大的懲罰閾值誘導更多樣化的搜索,而在算法后期通過設定較小的懲罰閾值來強化集中搜索。

算法3給出了種群更新算法的偽代碼,對于禁忌算法優化的兩個子代{sp+1,sp+2}與當前種群P={s1,…,sp},本文采用以下規則來確定sp+1,sp+2是否應該替換原有種群的個體:首先,計算種群更新的懲罰閾值(3行),將{sp+1,sp+2}與P合并構建初始候選集合Candidate,選擇目標值最小的個體添加到新種群中,并將該個體從Candidate中刪除(4行～7行),初始化懲罰解集合Penalized(8行);然后,不斷迭代,直到新種群的規模與原種群相同(9行～26行)。新種群個體選擇過程如下:首先,計算候選集合中個體與新種群的距離(11行),然后將距離小于選擇閾值D的個體從Candidate中刪除并添加至Penalized(12行～16行);最后,若Candidate不為空,則從Candidate中選取目標值最優的個體進入新種群,否則從Penalized選擇與新種群距離最遠的個體進入新種群(17行～26行)。

算法3種群更新算子。

1 輸入:種群P={s1,…,sp}和{sp+1,sp+2}。

2 輸出:新種群NewPop。

3 計算當前種群的懲罰閾值D

4 初始化候選集合Candidate=P∪{sp+1,sp+2}

5 選擇最優解Best=argmin{f(si)|si∈Candidate}

6 構建新種群NewPop={best}

7 Candidate=Candidate{best}

8 初始化懲罰解集合Penalized=?

9 while|NewPop|

10 foreach I∈Candidate do

11 I.dci=DCI(I,NewPop)

12 if I.dci

13 Candidate=Candidate{I}

14 Penalized=Penalized∪{I}

15 end if

16 end foreach

17 if Candidate≠? then

18 NewSelected =argmin{f(si)|si∈Candidate}

19 else

20 foreach I∈Penalized do

21 I.dci=DCI(I,NewPop)

22 end foreach

23 NewSelected=argmax{s.dci|s∈Penalized}

24 end if

25 NewPop=NewPop∪{NewSelected}

26 end while

27 Return NewPop

當懲罰閾值D較大時有利于增強搜索的廣度,D值較小時有利于在有潛力的區域進行集中搜索。綜合考慮以上特征,將D設定為隨算法運行時間線性減少的動態值,使算法逐漸從多樣化搜索轉化為集中化搜索,

(28)

式中:Tcur為當前運行時間;Ttotal為算法設定的最大運行時間;β為縮放系數,β=0.5;P0為初始種群。

3.6 算法復雜度分析

假設種群個體數為N、工件數量為n、機器數量為m。算法步驟包括種群的初始化、交叉父本匹配選擇操作、交叉操作、禁忌搜索操作、解碼操作、種群更新操作。算法采用隨機生成方法構建初始種群,該過程的時間復雜度為O(n×m×N);由于交叉父本匹配選擇操作與種群更新操作均需計算種群密度,根據個體之間的距離計算公式可得其在最壞情況下的復雜度為O(n2×m),個體與種群其他個體的距離最小值計算的時間復雜度為O(N×n2×m),進而可得種群密度的計算復雜度為O(N2×n2×m);交叉匹配選擇操作將種群隨機分為兩組的時間復雜度為O(N),最壞情況下,查找另一組中距離第1個父本最遠個體的時間復雜度為O(N);算法隨機執行兩種交叉操作,每種交叉操作均需訪問父染色體兩次,同時每次迭代只執行1次交叉操作,因此交叉操作的時間復雜度為O(n×m);禁忌搜索操作中鄰域操作的時間復雜度為O(n×m2),近似評價的時間復雜度為O(n×m),因此禁忌搜索操作的時間復雜度為O(n2×m3);在解碼過程中,每道工序需要從前后獲取空閑時間段,因此最壞情況下,該操作的時間復雜度為O(n2×m);在種群更新過程中,由于每次迭代只產生兩個子代個體,需要重新計算子代個體與原種群個體之間的距離,該過程的時間復雜度為O(N×n2×m),依次添加個體進入新種群過程的時間復雜度為O(N2)。由于初始種群構建與初始種群密度計算只執行一次,對算法整體性能影響不大,從以上分析可知,算法迭代過程中在禁忌搜索操作與種群更新操作上消耗的時間最多。

4 算法測試

HEA-MNS采用C++編程實現,并在Linux操作系統的計算機(CPU主頻2.3 GHz,內存16 GB)上使用g++編譯運行。算法參數設置如下:算法執行的總時間為10 min,種群規模為50,禁忌搜索的終止條件為未改進迭代次數達到1 000次;禁忌長度在兩個極值之間隨機選取,極值的設置分別為Lmin=0.5(10+n/m)和Lmax=1.5(10+n/m),其中n和m分別為工件數和機器數?？紤]到HEA-MNS的隨機特性,本文獨立求解每個問題算例10次。測試實驗包括兩組JSP基準算例,均來源于JSP基準算例庫(http://jobshop.jjvh.nl/)。第1組為LAWRENCE[39]設計的40個算例(LA01～LA40),第2組為ADAMS等[40]設計的5個基準算例(ABZ05～ABZ09)。

4.1 鄰域結構測試

ImpNS=

(29)

采用N5,sCET+2MT,N6,N7,PN7+2MT 5種鄰域結構對45個JSP算例的測試結果進行統計,如表2所示。表中算例順序按照PN7+2MT鄰域結構的相對改進百分比Imp降序排列得到,并以該順序對算例進行編號。按照該順序,繪制sCET+2MT,N6,N7,PN7+2MT 4種鄰域結構的相對改進百分比曲線,如圖14所示。對于所有45個JSP算例,PN7+2MT的相對改進百分比曲線均位于N6與N7鄰域結構之上,說明多工序精確聯動可有效改善鄰域結構性能。與sCET+2MT鄰域結構相比,PN7+2MT的相對改進百分比曲線具有明顯優勢,說明更大范圍的工序參與聯動操作可以大幅提升鄰域結構的搜索能力。

表2 4種鄰域結構的相對改進百分比對比

同一規模下,對各算例應用不同鄰域結構產生鄰域解的數量(SNS)av進行統計分析,結果如表3所示,其中:CR表示PN7+2MT相對于N7鄰域結構的鄰域解數量的消減率,計算公式如式(30)。因為sCET+2MT鄰域結構是在N5的基礎上聯動另外兩對工序,其鄰域解的數量相同,所以只列出N5的鄰域解數量。從表3可見,PN7+2MT的鄰域解數量大于N5和N6,表明PN7+2MT對解空間具有廣泛搜索的能力。相對于N7,PN7+2MT的鄰域規模減小約14%,說明所提無效移動判定條件可有效減少鄰域規模。

表3 4種鄰域結構的鄰域解數量對比

(30)

4.2 算法性能測試

從LA算例集中選取10個求解比較困難的算例測試算法性能,結果如表4所示,表中參數的詳細描述如下:Best為10次求解結果的最優值;Mav為10次求解結果的平均值;Tav(s)為10次獲得求解結果的平均用時(單位:s);MRE為求解結果UBsolve與已知下界LB相對偏差的平均值,其中相對偏差RE的計算如式(31)所示。表4的4個對比算法分別為文獻[41]具有局部搜索策略和多交叉算子的改進遺傳算法(Multi-Crossover Local Search Genetic Algorithm,MXLSGA)、文獻[42]基于多算子通信的序列禁忌搜索差分進化(Multi-operator Communication based Differential Evolution with sequential Tabu Search approach,MCDE/TS)算法、文獻[21]基于無延遲調度路徑重連與回溯禁忌搜索算法(Path Relinking based on Non-Delay scheduling and Backtracking Tabu Search algorithm,PRND/BTS)、文獻[22]基于切換策略的混合進化算法(Switching strategy-based Hybrid Evolutionary Algorithms,SHEA)。

表4 10個困難LA算例測試結果與其他算法結果比較

RE=(UBsolve-LB)/LB×100

(31)

從表4可知,HEA-MNS求得的最優結果Best與LB相對偏差的平均值b-MRE為0.01,且均低于其余4種算法的b-MRE。同時,HEA-MNS對其中9個算例的最優結果均優于或等于其他4種算法,僅對LA29算例的求解結果劣于SHEA算法。HEA-MNS求得的最優結果平均值Mav與LB相對偏差的平均值av-MRE均小于其他4種算法的b-MRE。

表5總結了HEA-MNS與其他4種算法求解全部LA算例的比較結果。實驗結果表明,HEA-MNS對每組算例求得的b-MRE優于或等于MXLSGA,MCDE/TS,PRND/BTS 3種算法,由于HEA-MNS未能求得LA29算例的已知最優解,使其在LA26-30算例組獲得的b-MRE值劣于SHEA算法,同時HEA-MNS對全部LA算例的av-MRE均小于其他4種算法的b-MRE,驗證了本文算法的穩定性。

表5 LA算例測試結果與其他算法結果比較

表6所示為對5個ABZ基準算例的測試結果,4種對比算法包括文獻[41]的MXLSGA算法、文獻[43]基于增強Lévy飛行和差分進化的混合鯨魚優化算法(hybrid Whale Optimization Algorithm enhanced with Lévy Flight and Differential Evolution,WOA-LFDE)、文獻[22]的SHEA算法、文獻[44]基于頻率分析算子的遺傳算法(Genetic Algorithm based on Frequency Analysis operator,FAGA)。從MRE角度進行統計對比分析,由于SHEA 算法未給出ABZ7算例的測試結果,在計算b-MRE時,設置該基準算例的相對偏差RE=0。HEA-MNS的b-MRE和av-MRE分別為0.62和0.73,其余3種算法的b-MRE分別為4.46,0.92,1.07,2.10,均大于本文算法。HEA-MNS求得了3個基準算例(ABZ5,ABZ6,ABZ9)的下界值,算例ABZ8的求解最優結果為667,其余4種算法分別為713,669,669,713,算例ABZ9的求解最優結果為678,其余4種算法分別為721,683,685,680,驗證了本文算法的有效性。圖15和圖16所示分別為HEA-MNS獲得ABZ8和ABZ9基準算例最優調度結果的甘特圖。

表6 ABZ算例測試結果與其他算法結果比較

4.3 參數敏感度分析

為了維持種群多樣性,通常需要在進化算法設計中添加額外的參數。在HEA-MNS中需要設置懲罰閾值D,該值越高,對決策空間的多樣化探索能力越強。由式(29)可知,懲罰閾值D由算法執行時間、初始種群密度與縮放系數β共同決定,上述測試中的β均設置為0.5。為了探討β對HEA-MNS性能的影響,將其設置6個級別{0,0.2,0.4,0.6,0.8,1.0},選用比較困難的ABZ8,LA29,LA40作為測試算例,算法其他參數設置與4.2節相同,圖17所示為不同β值下HEA-MNS求解3個測試算例的MRE曲線圖。

從圖17可知,3個測試算例的MRE值在β=0時最差,說明當HEA-MNS缺乏有效的種群多樣性控制策略時,算法性能會受到明顯影響;當β≥0.4時,算法均能保持較優的性能,說明在保證算法性能的前提下,β的取值范圍很大,HEA-MNS具有較強的穩定性。

5 結束語

本文提出一種HEA-MNS算法,通過將基于多工序精確聯動鄰域結構的禁忌搜索算法與多樣化控制的進化算法相結合,有效地解決了傳統進化算法在求解JSP中出現的局部搜索能力差、易早熟等問題。

本文從理論上給出了N7鄰域結構中工序無效移動的判定條件,據此設計了3對工序精確聯動的鄰域結構,有效擴展了鄰域搜索的范圍。為了解決傳統進化算法的種群多樣性隨進化過程快速下降的問題,本文提出基于鄰域懲罰的父本選擇算子與基于動態懲罰閾值的種群更新策略,有效提高了子代與種群的多樣性。

通過對比HEA-MNS與多種先進算法對不同規模JSP基準算例的測試結果,證明了算法的穩定性與有效性。然而,當前工作僅對算法進行了改進,并在靜態基準算例上進行了測試,初步驗證了HEA-MNS的可行性與有效性,下一步擬將HEA-MNS應用于實際生產調度場景,并考慮不確定因素的影響。