求解零空閑流水車間調度問題的離散正弦優化算法

趙 芮， 顧幸生

(華東理工大學 化工過程先進控制和優化技術教育部重點實驗室， 上海 200237)

流水車間調度問題(FSP)是被廣泛研究的組合優化問題，在制造系統和工業過程中起著非常重要的作用[1].在Johnson[2]的開創性工作之后，學者們提出了許多解決FSP的方法.作為FSP的重要分支，零空閑流水車間調度問題(NIFSP)假定機器不間歇地進行工作，而不必在相鄰的工件之間等待時間，廣泛存在于實際生產過程[3-4].在傳統行業如陶瓷熔塊生產、玻璃纖維加工、鑄造及集成電路生產中，由于某些制造環境中的技術要求或是使用了昂貴的機械設備，使得正在加工工件的機器不允許被停止.

早在1997年，Baptiste等[5]就證明了當機器數大于2時，使makespan準則最小的NIFSP是一個NP(非確定性多項式)難問題.因此如果只用精確算法和構造性啟發式算法進行求解，往往難以獲得理想的求解效果.在之前的研究中，國內外的學者提出了一些智能優化算法來求解NIFSP.為最小化最大完工時間，文獻[6]提出了一種混合離散粒子群(HDPSO)算法來求解NIFSP.文獻[7]提出混合離散差分進化(HDDE)算法.二者都采用了評估插入鄰域的加速方法，以降低時間復雜度.針對最大完工時間和總流水時間(TFT)兩個目標準則，文獻[8]利用差分進化算法來優化可變迭代貪婪算法的參數，提出了一種被稱為vIG_DE(variable Iterated Greedy Algorithm with Differential Evolution)的算法來求解NIFSP，并在Ruiz基準測試集上進行了測試.結果表明，vIG_DE算法與HDDE、HDPSO算法相比性能有所提高.文獻[9]利用離散人工蜂群DABC)算法來求解以總拖期時間(TTd)為目標的NIFSP，并首次在Taillard基準測試上為TTd準則提供了已知最優解.文獻[10]提出入侵雜草優化(IWO)算法，實驗結果表明，在Taillard測試集上，IWO算法的性能優于迭代貪婪(IG)算法和HDPSO算法，但其沒有使用任何局部搜索方法.文獻[11]提出雙種群分布估計算法(Bi-population EDA, BEDA)，以最小化TTd，BEDA由兩個子種群協同工作，通過對不同的概率模型進行采樣來生成兩個子種群.實驗表明，BEDA優于DABC算法.最近，針對以最小化最大完工時間為目標的NIFSP，文獻[12]提出了一種帶有局部搜索的離散煙花算法(DFWA-LS)，實驗結果表明DFWA-LS算法在尋優精度和穩定性上均不劣于HDPSO和IWO算法.文獻[13]提出具有混合節點和邊緣直方圖矩陣的模因算法，以最小化最大完工時間；文獻[14]提出基于元啟發式的混合離散學習算法，以最小化TTd.

受正弦余弦波形的啟發，Mirjalili[15]提出了一種基于智能群體的正弦余弦算法(SCA)，用以求解優化問題.此后，為了提高SCA的性能，Meshkat等[16]在SCA的基礎上提出了一種新的位置更新策略，該方法僅用正弦函數對個體位置進行更新，因此稱為正弦優化算法.研究結果表明，相比SCA，正弦優化算法(SOA)具有更快的收斂速度以及更高的尋優精度.到目前為止，對SOA的研究還很少，尚未應用于實際優化問題.但SCA在各種優化問題上的成功應用，例如參數優化[17]、風速預測[18]、電力系統的最優潮流問題[19]，可以揭示出SOA也有解決優化問題的潛力.

據知，目前尚無關于SOA求解NIFSP的研究成果.本文提出一種有效的基于IG算法的離散正弦優化算法來求解NIFSP，其目標是最小化最大完工時間.本文特別之處在于，將IG與SOA的位置更新策略相結合，并提出了新的位置更新策略來提高算法的搜索能力.加入了個體與當前最優解之間的交叉操作以增加算法多樣性，并根據種群數的增加定義了選擇策略，保留了種群中優秀的個體信息.此外，還采用了一種基于插入的局部搜索算法以增強局部開發能力.實驗表明，采用本文提出的離散正弦優化算法(DSOA)解決NIFSP是可行且有效的.

1 零空閑流水車間調度問題

NIFSP描述如下：n個工件集合J={1,2,…,n}按照相同的順序在m臺機器集合M={1,2,…,m}上進行加工.在任意時刻，每臺機器最多加工一個工件，每個工件最多只能被一臺機器加工[20].為了遵循零空閑的約束，機器i加工第j個工件的開始時間必須等于第j-1個工件的加工完成時間，即每臺機器加工任意兩相鄰工件時沒有空閑時間.

若用π={π(1),π(2),…,π(n)}來表示一個工件排序，用πE(j)={π(1),π(2),…,π(j)}來表示工件排序π的部分序列，其中1

F(πE(1),i,i+1)=p(π(1),i+1)

(1)

i=1,2,…,m-1

F(πE(j),i,i+1)=p(π(j),i+1)+

max{F(πE(j-1),i,i+1)-p(π(j),i),0}

(2)

i=1,2,…,m-1;j=2,3,…,n

(3)

式中：最大完工時間Cmax=p(π(1),1)+p(π(2),1)+F(πE(2),1,2)+F(πE(2),2,3).為了便于理解，將2個工件、3臺機器的Cmax計算過程在圖1中給出，圖中t為機器加工工件所花費的時間.

2 基本正弦優化算法

Meshkat等[16]在2017年提出的SOA是一種新穎的基于智能群體的隨機優化算法，單目標實參數數值優化的結果表明，與SCA相比，SOA的收斂速度更快，精度更高.SOA構造多個隨機的初始候選解，并使用基于正弦函數的數學模型對它們進行優化.該算法包括一些隨機變量和自適應變量，以平衡優化過程中的全局探索和局部開發.由于隨機優化算法將優化問題視為黑箱，使得SOA具有更大的靈活性，這意味著SOA可輕松應用于不同領域的問題.

SOA使用一組隨機解開始優化過程，假設在n維搜索空間中有NP個個體，則第i個個體的位置表示為

(4)

在SOA中，使用正弦函數來更新個體的位置：

(5)

(6)

式中:Tmax為最大迭代次數；a為常數，一般取值為2.優化過程開始時，r1的值較大，導致r1sinr2有較強的波動，使位置的隨機變化量更大，可以更好地進行全局探索.隨著時間的推移，r1的值減小，相應個體位置的隨機變化量也會減小，可以更好地進行局部開發；r2是[0,2π]之間的隨機數，它規定了個體下一個位置移動的方向；r3是[0,1]之間的隨機數，決定了當前個體下一位置的移動起點.當r3∈[0,0.5)時，使用式(5)中第一個公式更新位置，這時個體的下一位置將等于隨機個體的修改位置.換言之，在位置更新時首先隨機地選擇一個個體作為移動起點，然后用r1sinr2對其進行修改，得到新的位置.這樣一來，就能夠很好地覆蓋所有搜索區域，從而增強了算法的全局探索能力.當r3∈[0.5,1]時，使用式(5)中第二個公式更新位置，這時個體的下一位置將等于最優個體的修改位置.如此，能夠在最優個體周圍的區域進行搜索，從而增強了算法的局部開發能力.

基于以上分析，SOA的偽代碼為

Initialize a set of search agents

While (T

Evaluate each of the search agents

Update the best solution

Updater1,r2andr3

Update the position of search by Formula (5)

End while

Return the best solution as the global optimum

3 離散正弦優化算法

3.1 編碼及初始化評價

DSOA采用基于工件排序的自然數編碼，例如，個體Xi={3,5,1,4,2}表示工件在機器上的加工順序為π={3,5,1,4,2}.初始化種群均隨機產生.

3.2 位置更新策略

在SOA中，每一次的迭代過程，都需要對所有個體的位置進行更新.而個體位置的更新可以看作是對隨機個體(或最優個體)位置的一次修改，位置修改的變化量取決于r1sinr2波動的幅度.r1sinr2波動的幅度越大，其鄰域搜索空間越大.

為了解決組合優化問題，IG算法作為一種簡單、有效且易于適應的算法是非常可取的.基于NEH(Nawaz-Enscore-Ham)的基本原理，IG算法可以通過不斷地進行破壞操作和重構操作來獲得更好的解.在破壞(Destruction)階段，首先從π中移除隨機選擇的d個工件，產生兩個排序：πR和πD，其中πR包含d個移除工件，πD包含剩余的n-d個工件.在重構(Construction)階段，將πR的第一個工件插入πD中，產生n-d+1個部分解，選擇其中最優解用于下一次迭代；接下來將第二個工件插入πD中，……，直至πD中包含n個工件，形成完整的調度解.

在DSOA中，將IG算法的思想用于個體的位置更新，把IG的迭代過程視為個體的位置修改操作.每次迭代后，參數d更新如下：

(7)

(8)

式(7)中，d為對r1sinr2四舍五入取整、再取絕對值后的值.式(8)中，n為工件數，α為控制位置更新變化量大小的參數，由于參數d≤n，因此α∈(0,1).

因此，在DSOA中，個體的位置更新公式轉變為

Xi(T+1)=

(9)

完整的位置更新策略的偽代碼為

Input: removing size:d

1r3=a random number in[0,1]

2 if (r3<0.5)

3Xrand=select one solution fromXrandomly

4π=Xrand

5 else

6π=Xbest

7 end if

8πR=removedjob fromπrandomly

9πD=π-πR

10 fori=1:d

11 insert jobπR(i) into the optimal location which gives the minimum makespan in all possible positions ofπD

12 end for

13Xi=πD

Output: solutionXi

3.3 交叉操作

為了避免個體過早地陷入局部收斂狀態，需要對個體Xi實施交叉操作[21].此交叉操作的基本思想是：利用種群中的當前最優解Xbest作為參考，產生包含有最優解部分信息的新解，并用新解替換當前解.具體地，首先隨機地選擇兩個位置P1和P2，并將Xi中P1和P2之間的工件序列作為子序列sub1取出，再另外規定sub2是Xbest去除序列sub1中工件后所得的子序列.最后生成一個[0,1]之間的隨機數r，若r<0.5，則新解Yi由[sub1,sub2]組成；否則，Yi由[sub2,sub1]組成.交叉操作如圖2所示.

圖2 交叉操作

交叉操作偽代碼為

1 generate two positionsP1andP2randomly, 1≤P1≤P2≤n

3 sub2=Xbest-sub

4 if (r<0.5)

5Yi={sub1,sub2}

6 else

7Yi={sub2,sub1}

8 end if

Output: solutionYi

3.4 選擇策略

為了將種群中優秀的信息傳遞到下一代種群中，在交叉操作產生新的NP個個體Yi后，由原種群X與新種群Y共同組成候選解集合S，從S中選擇NP個個體組成下一代的種群X.選擇策略包括兩部分：① 采用精英保留策略保留最優個體，即候選解集合中適應度值最小的個體會被確定性地保留到下一代種群中；② 對剩下的NP-1個個體采用輪盤賭方式進行選擇，對于候選解中個體Si，其被選擇的概率滿足：

(10)

式中：f(Si)為第i個個體的適應度值.在候選解集合中，若個體適應度值較高，該個體被選擇的概率會提高，反之，被選擇的概率會降低.

3.5 迭代局部搜索

當找到的最優解不夠理想時，需要對最優解進行局部搜索，其目的是在相對最優解附近增強密集搜索，以期望找到更好的解.具體使用迭代局部搜索(ILS)算法[9]：

Input: sequenceπ;

1π′=π;s=1; counter=1

2 while (counter

3 insertjob=π′(s)

4 remove insertjob fromπ′

5π″=best sequence obtained by inserting insertjob in all possible positions ofπ′

6 ifCmax(π″)

7π=π″

8 counter=1

9 else

10 counter=counter+1

11 end if

12s=mod(s+1,n)

13 end while

Output: sequenceπ

3.6 DSOA的框架

DSOA的總體框架如圖3所示，終止條件由迭代次數進行設置.算法經過初始化、位置更新策略、交叉操作、選擇策略以及ILS局部搜索等步驟之后，最終返回最優解.

圖3 DSOA框架

4 仿真實驗

算法基于MATLAB 2017a實現，在處理器為Intel Xeon E5-2640，主頻為2.40 GHz，內存為64.0 GB的PC機下運行.為了檢驗DSOA算法求解NIFSP的性能，本文采用Taillard Benchmark問題[22]中的11種不同規模的共110個典型算例進行測試.為了證明所提出的算法的有效性并與其他方法進行性能比較，采用平均相對百分比偏差(ARPD)和標準偏差(SD)來衡量算法性能：

(11)

(12)

4.1 參數設置

算法的性能很大程度上取決于參數的設置，因此通過實驗來找到最佳的參數至關重要.為確定DSOA中參數α的取值，在算例Ta051、Ta061、Ta071及Ta081上進行參數測試，α取0.1、0.3及0.5，各進行10次仿真，每次仿真的最大迭代次數設為300，測試結果如表1所示，表中avg和min分別為算法運行結果的平均值和最優值.此外，為了更直觀地表示3種參數設置下DSOA的運行結果，用迭代收斂圖進行對比，如圖4所示.

由參數測試結果可見，α=0.1時，DSOA的運行結果較差，特別是當問題規模的機器數較大時，算法性能的差距更為明顯；α=0.3時，算法雖然在4個算例上都可以得到最優解，但其運行結果的平均值并不是最優的；而α=0.5時，算法在算例Ta061、Ta071上能獲得最優解，在算例Ta051、Ta081上雖然沒有獲得最優解，但運行結果的最小值與最優解的值相差不大，且算法運行結果的平均值最小，平均相對百分比偏差也最小，這表明α=0.5時算法的平均尋優精度最高，能使算法性能達到最優.基于以上分析，在后續實驗中參數α設置為0.5.確定了α的取值后，d和r1的取值可由式(7)、(8)得到.

圖4 α的測試結果

DSOA的其余參數設置如下：NP為30，Tmax設置為300，r2為[0,2π]之間的隨機數，r3為[0,1]之間的隨機數.

4.2 結果分析

為了驗證DSOA中各個策略的有效性，采用SOA、DSOA1、DSOA2、DSOA3及DSOA分別對11個規模大小不同的算例進行求解，每個算例獨立求解10次，最大迭代次數為300.DSOA1表示只改變了位置更新策略的離散正弦優化算法；而DSOA2不僅改變了位置更新策略，還加入了交叉操作和選擇策略；DSOA3則表示改變了位置更新策略，并使用了ILS局部搜素，但沒有加入交叉操作和選擇策略的離散正弦優化算法.需要說明的是，SOA求解NIFSP的編碼方法采用的是基于最大值排序(LRV)規則.

表2列出了這11個算例的運行結果，選取了10次運行結果中的平均值進行對比，其中性能提高百分比(PIP)定義為

(13)

顯然，當i的取值為1時，表示的是DSOA1相比于SOA提高的性能百分比.

為了更直觀地觀察各算法的求解效果，以Ta051為例，畫出各個算法求解的收斂曲線，如圖5所示.結合表2和圖5，不難看出，相比于SOA，DSOA1、DSOA2、DSOA3以及DSOA的性能都有大幅度提升.對比SOA和DSOA1的實驗結果可以看出，重新定義的位置更新策略更適合于求解調度問題，算法的探索能力明顯增強.相比于SOA，DSOA2和DSOA1的性能提升相差不大，但交叉操作和選擇策略卻并非無效的，因為從DSOA3來看，DSOA3沒有經過交叉操作和選擇策略，在種群更新策略之后直接進行局部搜索，求解效果卻不如DSOA，導致這一結果的原因可能是種群多樣性被破壞，算法難以跳出局部最優.而DSOA2和DSOA的實驗結果表明，加入ILS局部搜索后能夠加快收斂速度，提高尋優精度.

表2 SOA與DSOAi的比較

圖5 算例Ta051的收斂曲線

從表2中還可以看出，對每一個算例，DSOA都能找到比SOA更令人滿意的最優解，算例Ta051的PIP達到最大，為15.415%，平均性能提高了9.260%.這表明對SOA的改進總體來說是有效的，對于NIFSP的求解，DSOA相比于SOA，具有更好的尋優精度.

為了進一步測試所提出的DSOA求解NIFSP的有效性，將DSOA與IWO[10]和DFWA-LS[12]算法進行了比較.由于在另一臺計算機上實現同一算法的結果可能不同，所以，為了確保操作環境的一致性，算法比較的公平性，重新實現了對比算法，并在相同的PC機上運行.對比算法的參數設置與文獻[10,12]中一致.此外，為了避免測試的隨機性帶來的誤差，每次不同的仿真都獨立運行10次.測試結果見表3，其中計算ARPD的式(11)中，C*采用這3種算法所能找到的最優解進行計算.

表3 3種算法的結果對比

表3顯示了3種算法求解不同規模Taillard算例的性能.如表3所示，所提出的DSOA的平均ARPD值為0.15，小于DFWA-LS和IWO的對應值0.31、2.21；DSOA的平均SD值為0.050，也小于由DFWA-LS和IWO得到的0.090、0.351.以上結果說明，在最小化最大完工時間的目標上，DSOA求解不同規模問題的尋優精度和穩定性均優于兩種對比算法.此外，觀察可發現DSOA在n=20,m=5、n=50,m=5以及n=100,m=5這3個問題上的SD幾近于0，隨著問題規模的增大，其SD值相對于IWO和DFWA-LS的更小，這表明提出的DSOA在解決相對大規模的問題時仍然具有出色的穩定性.這些性能受益于DSOA的搜索策略以及與IG算法的結合.因此，DSOA可較好地應用于NIFSP，能在保證算法穩定性的情況下找到比IWO和DFWA-LS更優的解.

圖6顯示了3種算法在求解算例Ta058時的收斂性.可以看出，與其他算法相比，DSOA具有相對較快的收斂速度和更高的精度，并且能夠快速跳出局部最優.

圖6 算例Ta058的收斂曲線

5 結語

本文提出了一種新穎的DSOA來求解以最小化最大完工時間為目標的零空閑置換流水車間調度問題.所提出算法的主要思想是將IG算法嵌入到SOA的位置更新策略中.在DSOA中，隨機選擇的個體和不斷變化的去除工件數提高了算法的全局搜索能力，避免了算法陷入局部最優狀態.同時，IG算法增強了算法的局部搜索能力.此后，交叉操作可以增加算法多樣性，幫助算法跳出局部最優，選擇策略能夠選擇相對優解構成新的種群.最后，ILS局部搜索的重點在于圍繞最優解進行搜索，從而提高了算法的收斂速度.這4個基本策略的集成為算法有效性提供了保證.

在Taillard Benchmark算例上測試了提出的DSOA，并通過與IWO和DFWA-LS的對比，驗證了DSOA在求解NIFSP上的有效性和優越性.下一步的研究將包括：① ILS是局部搜索的一種基本策略，考慮用更高效的局部搜索，算法的性能將更為優越；② 考慮用DSOA求解其他復雜的調度問題.