基于多目標離散正弦優化算法的混合零空閑置換流水車間調度

趙 芮， 郎 峻， 顧幸生

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

零空閑置換流水車間調度問題（No-idle Permutation Flow Shop Scheduling Problem, NPFSP）是廣泛存在于集成電路制造、玻璃纖維加工和鑄造等現代工業的組合優化問題[1]。NPFSP 通常要求機器都處在零空閑狀態，即每臺機器要不間斷地加工工件[2]。然而在實際的生產過程中，往往只有部分特定的機器因設備昂貴或技術限制等受到零空閑條件的約束。混合零空閑置換流水車間調度問題（MNPFSP）是NPFSP 的擴展，其要求部分機器受零空閑條件的約束，其余機器為常規機器。2014 年，Pan 等[3]首先提出了混合零空閑置換流水車間調度問題，并證明其為NP 難問題。然而隨著規模的增大，MNPFSP 的復雜性呈指數增長，精確算法不再適用，因而，通過高效的智能優化算法解決MNPFSP 是學術界和工業界的重要課題。

目前，將智能優化算法應用于MNPFSP 的研究并不多。Pan 等[3]提出了一種改進的迭代貪婪（Iterated Greedy, IG）算法以最小化最大完工時間（makespan）。張曉霞等[4]提出了一種改進的分布估計算法（Distribution Estimation Algorithm, EDA），將禁忌搜索（Tabu Search, TS）算法引入了EDA，通過禁忌列表控制EDA 的進化方向，并以最小化makespan 為優化目標，對比了在7 種不同零空閑情況的MNPFSP 上的求解效果。Cheng 等[5]基于云理論設計了一種新的IG 算法對MNPFSP進行求解。Rossi 等[6]開發了一種構造性啟發式算法，求解以最小化總流水時間（Total Flow Time, TFT）為目標的具有序列相關準備時間的MNPFSP。

本文針對車間生產效率和客戶滿意度這兩個不同的生產需求，將makespan 和最大拖期（Maximum Tardiness）作為優化目標。兩者在多目標優化問題的研究中已有相關成果，如Santosa 等[7]針對最小化makespan、總拖期時間和總機器空閑時間為目標的Flow shop 調度問題，提出了一種離散粒子群算法（Discrete Particle Swarm Optimization, DPSO）。Pan 等[8]為解決零等待Flow shop 問題，提出了一種離散差分進化（Discrete Differential Evolution, DDE）算法以實現最小化makespan 和最大拖期的目標，同時設計了幾種加速方法提高算法效率。Lei 等[9]針對以最小化總 拖 期 時 間（Total Tardiness, TTd）、最 大 拖 期 和makespan 為目標的多目標混合Flow shop 問題，提出了一種具有全局交換功能的鄰域搜索策略（Neighborhood Search with Global Exchange, NSG），通過全局交換和鄰域搜索相結合的方式提高解的質量。Murilo等[10]提出了一種多目標進化算法（Multi-Objective Evolutionary Algorithm, MOEA），在分布估計算法中引入對數據排名的Mallows 概率模型，解決了以makespan、總流程時間、總拖期時間為目標的多目標PFSP。

正 弦 優 化 算 法（Sine Optimization Algorithm,SOA）[11]是通過改進正余弦算法（Sine Cosine Algorithm,SCA）[12]的一種利用正弦函數對個體進行位置更新的全局優化方法。本文在SOA 的基礎上設計了一種多目標離散正弦優化算法（Multi-objective Discrete Sine Optimization Algorithm, MDSOA）。引入IG 算法的破壞重構機制對SOA 的位置更新進行了改進，并設計適用于離散調度問題的位置更新策略。構建外部檔案集（AS）儲存Pareto 解集，通過AS 選取的當前最優解與個體之間的交叉操作保留最優解的部分信息。利用快速非支配排序（FNDS）和擁擠距離篩選更新種群，以此提高種群的分布性。

1 多目標混合零空閑置換流水車間調度問題

MNPFSP 可以被描述如下：n個工件J={1,2,···,n}按照相同的順序在m臺機器（M={1,2,···,m} )上進行加工，所有機器上工件的加工次序相同[13]。約束條件有：（1）在任意時刻，每個機器僅加工一個工件，每個工件僅在一臺機器上加工；（2）在加工過程中部分機器處于零空閑狀態，其余機器為置換流水車間中的常規機器，零空閑狀態的機器在加工兩相鄰工件時不間斷，即第j個工件在機器i上的開始時間必須等于第j?1 個工件在機器i上的完工時間。

在生產過程中，最小化最大完工時間（Cmax）可以提高生產效率和機器利用率，最小化最大拖期（Tmax）可以使工廠最大程度地避免延期交貨，滿足客戶對產品工期的要求。因此，本文將最小化最大完工時間和最大拖期作為MNPFSP 的優化目標，兩者之間的關系已經被證明是互相矛盾的[14]。用 π={π1,π2,···,πl,···,πn} 表示工件加工序列，用 [l] 表示排序中處于l的工件，即 πl；用pi,[l]表示工件 πl在機器i上的加工處理時間；用Si,[l]和Ci,[l]分別表示工件 πl在機器i上加工的開始時間和完成時間；T[l]和d[l]分別為工件 πl的拖期時間和交期；ai表示受零空閑條件的約束，工件 πl在機器i上延遲加工的總時間；將具有零空閑狀態的機器集合記為M′?M。加工序列 π的最大完工時間和最大拖期的計算公式如下[3]：

式（1）計算 π1在首臺機器上的開始和完成時間；式（2）計算 π1在（除首臺機器外）其他各臺機器上的開始和完成時間；式（3）計算每個工件在首臺機器上的開始和完成時間；式（4）表示若機器2 受零空閑條件的約束，則要保證各個操作之間是零空閑的，πl之 前 的 工 件 開 始 時 間 都 需 要 右 移 或 延 遲a2=max{C1,[l]?C2,[l?1],0} ，若機器2 為常規機器，則不需要延遲加工，a2=0 ；式（5）類似，工件因機器i為零空閑機器而延遲加工的時間為 m ax{Ci?1,[l]?(Ci,[l?1]+ai?1),0} ，ai?1為上游零空閑機器導致的延遲時間，因此ai為總延遲時間。式（6）和式（7）規定了工件序列 π 的最大完工時間和最大拖期。則優化的目標函數表示為

2 正弦優化算法（SOA）

SOA 是基于SCA 的一種新的群智能優化算法，對SCA 的位置更新策略進行了改進，僅使用正弦函數對初始候選解進行優化。相較于SCA，SOA 的參數更少，算法更簡便，同時收斂速度更快，精度更高[11]。

SOA 使用一組隨機解作為初始種群，假設在n維搜索空間中有Psize個個體，第i個個體的位置表示為

3 求解多目標MNPFSP 的多目標離散正弦優化算法

3.1 編碼方案與初始化

對工件排序采用自然數編碼的方法，排序的序列由一串不重復的整數表示，并且每個整數都對應一道工序，這串序列代表一個解的編碼，例如解X={4,1,3,6,2,5} 表示工件加工順序為 π={4,1,3,6,2,5} 。初始種群均隨機產生，隨機的初始種群更有可能在迭代之初遍布可行域的搜索空間，提高算法的多樣性。

3.2 外部檔案集的構建、更新與篩選

在求解多目標調度問題時，往往無法獲得一個最優解，而是在相互競爭的目標中進行取舍從而獲得一組Pareto 解集。在MDSOA 中構建一個外部檔案集（AS）存儲Pareto 解，避免算法中的優良個體在種群迭代更新中被遺失。在種群初始化后，根據個體子目標的適應度值進行非支配排序，得到Pareto 解集作為初始AS。同時在每代新種群產生后要對AS 進行篩選和更新，具體方法是將當前種群X中不受AS 支配的解放入AS 中，并把AS 中被種群X支配的解剔除出去，使AS 中始終存放著算法所求的Pareto 解集。AS 創建及更新的偽代碼如下：

由于在求解多目標問題時可以得到大量滿足約束條件的Pareto 解，為了控制計算的復雜度，將外部檔案集的規模限制為K。若某次更新時，AS 的長度超過了K，則需要刪除AS 中擁擠距離更小的解，以此控制外部檔案集的規模。

3.3 基于破壞重構的位置更新策略

MDSOA 中個體的位置更新公式如下：

其中：DestructConstruct（X，d）表示對個體X進行破壞重構操作，破壞工件個數為d。當r3<0.5 時，利用隨機個體進行位置更新，否則使用最優個體進行位置更新。位置更新策略的偽代碼如下：

3.4 基于FNDS 和擁擠距離的選擇策略

在各類處理多目標問題的優化方法中，基于Pareto 解的方法具有很大的優勢，Deb 等[15]提出的非支配排序遺傳算法II（Non-dominated Sorting Genetic Algorithm II，NSGA-II）具有很好的求解效果。NSGAII 根據對種群的非支配排序和擁擠距離計算來選擇下一代種群，可以在保留種群精英解的同時，保持解具有良好的分布性和收斂性。

快速非支配排序（FNDS）的思想是根據支配關系將種群劃分為若干個不同等級的Pareto 前沿。在多目標搜索空間中，兩個個體之間的支配關系需要針對每個目標進行比較。在這種情況下，當且僅當個體x的所有目標值都比個體y更好或相等，且至少有一個目標值優于y時，才有x支配y[16]。

Deb 等[17]給出了擁擠距離的定義：在同等級的Pareto 前沿上，與個體i相鄰的兩個個體在各子目標上的距離差之和，計算公式如下：

圖1 基于FNDS 和擁擠距離的選擇策略Fig. 1 Selection strategy based on fast non-dominated sorting and crowding distance

3.5 多目標離散正弦優化算法

MDSOA 的流程如圖2 所示，終止條件為滿足最大迭代次數。算法首先初始化參數并隨機生成初始種群，再對初始解集進行非支配排序，構建初始AS。當r3<0.5 時，選擇最優個體根據位置更新策略對種群進行更新，否則選擇隨機個體對種群進行位置更新，得到種群X。交叉策略采用的是將種群的最優解Xbest作為參考，與種群中的其他解進行兩點交叉，使產生的新解包含最優解的部分信息，并用新解代替當前解。將交叉操作得到的種群Y與位置更新策略得到的種群X進行非支配排序和擁擠距離計算，通過選擇策略保留最優解。同時更新AS 和種群的最優解，完成一次迭代。最終返回最優的Pareto解集。

圖2 MDSOA 算法流程圖Fig. 2 Flowchart of the MDSOA

4 仿真及分析

算法基于Matlab 2017a 實現，在處理器主頻為2.40 GHz，內存為64.0 GB 的PC 機上運行。為了檢驗MDSOA 求解MNPFSP 的有效性，本文選取Taillard Benchmark[18]中不同問題規模的第一個算例作為測試對象，其中，規模用工件數(n)×機器數(m)表示，并將仿真結果與NSGA-II 算法[15]和NSGA-III[19]算法進行比較。

由MNPFSP 的實際生產過程可知，一般第一臺機器處于零空閑狀態，最后一臺機器不處于零空閑狀態，并且處于零空閑狀態的機器占比較小，因此，本文設置兩臺機器處于零空閑狀態：除了第一臺機器外，再隨機選擇一臺（除最后一臺機器外）機器，并用Mset表示。按算例序號的升序隨機生成的Mset數據依次為 [ 2 5 11 3 8 15 4 7 10 3 6] 。對于最大拖期Tmax這個性能指標，需要規定每個工件的交期，交期的構造公式如下[20]：

4.1 參數設置

MDSOA 有4 個主要參數：最大迭代次數（ M axGen ）、種群規模（Psize）、控制參數（ β ）和外部檔案集長度(K)。為確定各個參數的取值，為每個參數設置3 個不同水平，選擇L9(34) 正交試驗表對參數的各種組合進行討論，如表1 所示。對于每種參數組合，算法均獨立運行20 次，并將仿真結果合并為該參數組合下的Pareto 解集，再合并每種參數組合下的解集，最終得到該算例下的Pareto 解集。將每組參數組合下Pareto 解的數量占算例Pareto 解集的百分比作為響應變量（Response Variable, RV），RV 值越大，參數越好。正交表及算例Ta051 的正交結果見表2（括號中的數值為選擇的等級(Level)）。

表1 正交設計因子水平表Table 1 Factors and levels for orthogonal design

表2 算例Ta051 的正交試驗結果Table 2 Orthogonal test results for the instance of Ta051

由表2 可以看出，對算法影響較大的是 MaxGen和Psize。較多的迭代次數有利于種群的充分進化，合適的種群規模有利于平衡算法每次迭代的搜索廣度。其次是K，K太小可能造成收斂緩慢，K太大可能造成早熟收斂，故K的取值要適中。 β 對MDSOA性能的影響較小。基于以上分析，MDSOA 的參數如下： M axGen=300,Psize=50, β=0.5,K=40 。

為驗證MDSOA 求解多目標MNPFSP 的有效性，將其與兩種典型的多目標算法NSGA-II 和NSGAIII 進行比較，其中NSGA-II 和NSGA-III 的交叉操作和變異操作均采用的是最常見的兩點交叉和單點變異。NSGA-II 和NSGA-III 的參數設置見表3。

表3 NSGA-II 和NSGA-III 的參數設置Table 3 Parameter settings of NSGA-II and NSGA-III algorithms

4.2 算法性能測試

本文采用均勻性指標（SM）[21]、反世代距離（IGD）[21]以及非支配解的數量(NNDS[22])和距離指標（ D IR）[22]來評價3 種算法得到的Pareto 解集的質量，并分別統計其平均值（Avg）和標準偏差（SD）。

（1）SM 用于度量算法求得的近似Pareto 解集中每個解到其他解的最小距離的標準差，計算公式如下：

（2）IGD 是世代距離（GD）的一種變體，它不僅可以度量算法求得的近似Pareto 解集與最優Pareto解集的距離，還可以衡量解集的多樣性、分布性。IGD 的計算公式如下：

（4） D IR用于度量Sh相對于參考解集S?的距離：

表4 列出了3 種算法支配解個數的比較結果，最優結果以黑體表示。表5 列出了3 種算法均勻性指標的比較結果。

表5 SM 的計算結果Table 5 Computational results of SM

從表4 可以看出，本文提出的MDSOA算法在11 個算例中都獲得了非支配解個數的最優值。雖然在算例Ta021、Ta041、Ta051、Ta081 和Ta101上沒有獲得更小的標準偏差（SD），但從平均結果（Avg)來看，MDSOA 的SD 值要優于其他兩個算法。因此可以說明，相較于另兩種多目標優化算法，MDSOA能夠找到更多的非支配解。

表4 NNDS 的計算結果Table 4 Computational results of NNDS

從表5 可以看出，算例Ta031、Ta061、Ta091 和Ta101 中較優的SM 值是由NSGA-II 和NGSA-III 獲得的，其余算例均是MDSOA 獲得了更優的SM 值。因此，對于大部分算例，MDSOA 求得的解在解空間中分布得更均勻。

表6 列出了3 種算法距離指標的對比結果。可以看出，除了在算例Ta041 中MDSOA 求解的SD 值稍遜于NSGA-II 算法外，剩余算例中MDSOA 都獲得了更優的結果，因此MDSOA 可以獲得更接近真實Pareto 前沿的解集。

表6 D IR 的計算結果Table 6 Computational results of DIR

表7 列出了3 種算法反世代距離的比較結果。對于測試的11 個算例，MDSOA 求解到的 IGD 值均明顯小于NSGA-II 和NGSA-III 算法，并且隨著算例規模的增大，IGD 值之間的差距也逐漸增大。同時，除了算例Ta001 和Ta031 的SD 最小值由另外兩個算法獲得，其余算例的SD 最小值均由MDSOA 算法獲得。由此可見，MDSOA 算法具有良好的綜合性能，且算例的規模越大，其性能越優越。

表7 IGD 的計算結果Table 7 Computational results of IGD

為直觀地對比3 種算法性能，選擇部分獲得最優IGD 值的非支配解集繪制成Pareto 前沿，如圖3所示。顯然在解的質量和覆蓋均勻性方面，MDSOA要優于其他兩個算法，并且隨著問題規模的逐漸增大，算法性能差距越大，MDSOA 越能得到更優的Pareto 解集。

圖3 MDSOA、NSGA-II 和NSGA-III 獲得的Pareto 前沿Fig. 3 Pareto front obtained by MDSOA, NSGA-II and NSGA-III algorithms

基于以上分析，可以得出，在解的數量、分布性以及收斂性上，MDSOA 都顯示了比NSGA-II 和NGSAIII 更優的算法性能，進一步表明了MDSOA 算法求解MNPFSP 的有效性。

5 結 論

本文基于正弦優化算法提出了一種多目標離散正弦優化算法（MDSOA）以提高求解MNPFSP 的效率和解的質量。通過在SOA 位置更新策略中引入IG 算法的破壞重構機制，平衡算法的全局搜索和局部搜索能力。設計交叉操作增加算法多樣性，幫助算法跳出局部最優，同時，構建外部檔案集，通過基于FNDS 和個體擁擠距離的選擇策略保留精英解。使用Taillard Benchmark 算例，對比MDSOA、NSGAII 和NGSA-III 在4 個評價指標上的性能表現，相較于NSGA-II 和NGSA-III 算法，MDSOA 獲得的解數量更多，更接近真實Pareto 前沿，分布更均勻，算法的收斂性也更好。生產調度作為一類多約束問題，求解十分困難，本文算法所針對的模型更接近實際生產過程，同時，滿足生產過程中的多個生產目標的求解需要。