求解零空閑流水車間調(diào)度問題的離散正弦優(yōu)化算法

趙 芮， 顧幸生

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

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

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

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

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

1 零空閑流水車間調(diào)度問題

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 基本正弦優(yōu)化算法

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

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

(4)

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

(5)

(6)

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

基于以上分析，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 離散正弦優(yōu)化算法

3.1 編碼及初始化評價

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

3.2 位置更新策略

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

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

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

(7)

(8)

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

因此，在DSOA中，個體的位置更新公式轉(zhuǎn)變?yōu)?/p>

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 交叉操作

為了避免個體過早地陷入局部收斂狀態(tài)，需要對個體Xi實施交叉操作[21].此交叉操作的基本思想是：利用種群中的當前最優(yōu)解Xbest作為參考，產(chǎn)生包含有最優(yōu)解部分信息的新解，并用新解替換當前解.具體地，首先隨機地選擇兩個位置P1和P2，并將Xi中P1和P2之間的工件序列作為子序列sub1取出，再另外規(guī)定sub2是Xbest去除序列sub1中工件后所得的子序列.最后生成一個[0,1]之間的隨機數(shù)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 選擇策略

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

(10)

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

3.5 迭代局部搜索

當找到的最優(yōu)解不夠理想時，需要對最優(yōu)解進行局部搜索，其目的是在相對最優(yōu)解附近增強密集搜索，以期望找到更好的解.具體使用迭代局部搜索(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所示，終止條件由迭代次數(shù)進行設置.算法經(jīng)過初始化、位置更新策略、交叉操作、選擇策略以及ILS局部搜索等步驟之后，最終返回最優(yōu)解.

圖3 DSOA框架

4 仿真實驗

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

(11)

(12)

4.1 參數(shù)設置

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

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

圖4 α的測試結果

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

4.2 結果分析

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

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

(13)

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

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

表2 SOA與DSOAi的比較

圖5 算例Ta051的收斂曲線

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

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

表3 3種算法的結果對比

表3顯示了3種算法求解不同規(guī)模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求解不同規(guī)模問題的尋優(yōu)精度和穩(wěn)定性均優(yōu)于兩種對比算法.此外，觀察可發(fā)現(xiàn)DSOA在n=20,m=5、n=50,m=5以及n=100,m=5這3個問題上的SD幾近于0，隨著問題規(guī)模的增大，其SD值相對于IWO和DFWA-LS的更小，這表明提出的DSOA在解決相對大規(guī)模的問題時仍然具有出色的穩(wěn)定性.這些性能受益于DSOA的搜索策略以及與IG算法的結合.因此，DSOA可較好地應用于NIFSP，能在保證算法穩(wěn)定性的情況下找到比IWO和DFWA-LS更優(yōu)的解.

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

圖6 算例Ta058的收斂曲線

5 結語

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

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