求解異構并行機調度問題的混合煙花算法

石 慶 民

(新鄉學院人事處 河南 新鄉 453000)

0 引 言

異構并行機調度(Unrelated parallel machine scheduling，UPMS)是實際生產制造過程中常見的一類調度問題，具有廣泛的工程應用背景，例如：云計算、半導體加工、紡織生產等[1-2]。傳統的UPMS問題多以時間指標為優化目標，如makespan、提前/拖期總和等。現如今，低碳節能調度已成為實現低碳制造的重要途徑，并已引起社會各界的廣泛關注[3]。與此同時，快速響應客戶需求也構成了企業競爭力的重要組成部分。因此，研究以能耗和延遲成本為目標的UPMS問題具有較高的理論意義和應用價值。

UPMS已經被證明是具有NP-hard性質的組合優化問題，因而設計有效的求解算法成為這一領域的研究重點。目前，求解UPMS問題的算法分為三類：精確求解算法、啟發式算法和現代優化算法[4]。精確求解算法包括分枝定界、動態規劃等，這類算法能夠獲得該問題的精確解，但實用性較弱，無法在合理的時間內求解中大規模算例[5]。啟發式算法是利用調度規則對問題進行快速求解，具有運算時間短、求解精度較高等優點，但該類算法對問題的設置具有嚴格限制，通用性較差[6]。近年來，隨著現代優化算法的快速發展，學者們逐步開始使用智能算法來求解這一問題，相關算法包括遺傳算法、禁忌搜索算法、模擬退火算法等。研究表明，這類算法能夠在較短的時間內獲得中大規模算例的滿意解，且具有較強的通用性[7]。因此，本文基于基本煙花算法(Firework algorithm,FWA)設計了改進算法用于求解考慮節能的異構并行機調度問題。

FWA是Tan等[8]提出的新型現代優化算法，屬于群智能優化算法，該算法受煙花爆炸行為的啟發通過煙花爆炸生成火花這一過程實現迭代進化。FWA具有收斂速度快、求解精度高等優點，目前已經在數據挖掘、神經網絡訓練、車輛路徑規劃等多個方面獲得了成功應用[9]。即便如此，FWA在中大規模測試問題的求解過程中仍然存在易陷入局部最優、收斂速度慢等缺點。因此，提升FWA的尋優性能仍然是該鄰域的一個重難點。

基于以上研究內容，本文考慮一類加工時間可控的并行機調度問題，并以能耗和拖期懲罰成本為優化目標構建了混合整數線性規劃(Mixed Integer Linear Programming,MILP)模型。算法設計中創建了特定的編解碼方法以表示問題的解，同時融入了反向學習初始化方法以提升初始解的質量，并構建了基于變鄰域搜索算法的局部優化流程用以強化基本算法的尋優性能。仿真測試中，利用正交實驗對算法參數進行校驗，并將HFWA的優化結果與多種經典的優化算法進行比較，測試結果驗證了HFWA的可行性和有效性。

1 問題模型

1.1 問題描述

在UPMS-ETC問題中，假設有n個任務需要在m臺機器上加工，各個機器存在多種不同的加工速度，各機器在不同加工速度下的能耗存在差異性。此外，已知各個訂單的交貨期，晚于規定的時間交貨將會產生延遲懲罰成本。

本文以最小化所有任務的加工能耗成本和延遲懲罰成本總和為目標，決策內容包含以下三點：1) 各個機器所加工的任務集合；2) 各個機器上加工任務的排序；3) 各任務的加工速度。

1.2 數學建模

為了構建UPMS-ETC問題的MILP模型，定義以下數學符號。

問題參數：

m：機器總數，i∈{1,2,…,m}。

n：任務總數，v∈{1,2,…,n}。

l：加工位置編號，l∈{1,2,…,n}，各臺機器均設有n個加工位置。

Qi：機器i可選的速度類別總數，s∈{1,2,…,Qi}。

lv：任務v的標準加工時間。

ris：機器i以第s種加工速度的取值。

tvis：任務v在機器i上以速度ris加工時相應的加工時間，tvis=lv/ris。

eis：機器i在單位時間內以速度ris運行所產生的能耗。

dv：任務v的交貨期。

Cv：任務v的完成時間。

τv：任務v的延誤時間，τv=max{0,Cv-dv}。

μv：任務v單位時間的延誤懲罰成本。

決策變量：

xvils：0-1變量，任務v在機器i以速度ris進行加工，xvils取值為1；否則，xvils取值為0。

Bil：機器i上位置l任務開始加工的時間。

基于以上問題描述和符號定義，構建了UPMS-ETC問題的MILP模型。

目標函數：

(1)

約束條件：

(2)

(3)

(4)

τv≥Bi,l+tvis·xvils-(1-xvils)·M-dj
i∈{1,2,…,m},s∈{1,2,…,Qi},v∈{1,2,…,n},
l∈{1,2,…,n}

(5)

τv≥0v∈{1,2,…,n}

(6)

xvils∈{0,1},i∈{1,2,…,m},s∈{1,2,…,Qi},
v∈{1,2,…,n},l∈{1,2,…,n}

(7)

式(1)表示目標函數，即所有任務的加工能耗成本和延誤懲罰成本之和；式(2)用于確保各個加工任務所選的加工機器、加工位置和加工速度的唯一性；式(3)表示各臺機器上每個加工位置處最多只能安排一個加工任務；式(4)用于計算機器i上位置l處加工任務的開始時間；式(5)-式(6)用于計算各個任務的延誤時間；式(7)定義了決策變量的取值范圍。

2 基本煙花算法

FWA以煙花表示待優化問題的解，并通過煙花爆炸生成火花這一過程迭代進化。首先依據各煙花個體的評價函數值計算相應的爆炸火花數目和爆炸半徑數值；其次，采用均分分布和高斯分布生成新個體；最后，利用選擇規則構建新種群以進入下一次迭代。具體實現步驟歸納如下：

Step1初始化算法參數，利用隨機方法構建初始種群并對每個煙花個體進行評價。

Step2依據評價函數計算當前種群中各個煙花的爆炸火花數目和爆炸半徑。

Step3利用均分分布和高斯分布生成爆炸火花和高斯變異火花。

Step4利用煙花、爆炸火花和高斯變異火花構建候選種群，并依據篩選規則生成新種群。

Step5判斷是否滿足終止條件。若滿足則停止迭代搜索并輸出當前最優解；否則，轉Step 2。

FWA的核心內容包含以下四點：(1) 計算爆炸火花數目和爆炸半徑；(2) 生成爆炸火花；(3) 生成高斯變異火花；(4) 基于選擇規則構建新種群。

計算爆炸火花數目和爆炸半徑：以N表示煙花數目，以xi表示當前種群中的第i個煙花，N個煙花數目生成的爆炸火花總數記為SN。以最小化問題為例，對于煙花xi，以符號fi表示其評價函數值，以Si和Ri表示由xi生成的爆炸火花數目和相應的爆炸半徑。Si和Ri的計算公式歸納如下：

(8)

(9)

式中：fmax和fmin分別表示當前煙花種群中最大和最小的評價函數值；A表示基本爆炸半徑，需要預先設定；ε為機器最小值，用以避免除零操作。同時，為了對上述Si計算公式做修正以確保優勢煙花個體生成的火花數不過多、劣勢煙花個體生成的火花數不過少，修正公式為：

(10)

式中：常數a和b滿足0

生成爆炸火花：依據計算所得Si和Ri值，生成當前煙花個體xi對應的爆炸火花。以D表示待優化問題的維度，隨機選取維度集合{1,2,…,D}的一個子集Dz，利用式(11)更新xi的相應維度的編碼數值取值生成新個體。

xiz←xiz·rand(-1,1) ?Z∈Dz

(11)

式中：rand(-1,1)表示區間[-1,1]的均分分布。

生成高斯變異火花：首先，利用隨機的方法從當前種群中選擇GN個煙花；其次，對所選的各個煙花，隨機選擇一些維度并利用高斯分布進行更新，以Dz表示所選煙花xi的一個維度子集。相應的更新公式如下：

xiz←xiz·N(1,1) ?Z∈Dz

(12)

式中：N(1,1)表示均值為1、方差為1的高斯分布。

基于選擇規則構建新種群：首先，將當前煙花種群、爆炸火花種群和高斯變異火花種群合并，將其中最優個體直接選入下一代種群；其次，依據輪盤賭規則選出其余N-1個解構建新種群。以K表示候選種群規模，則其中第i個煙花解xi被選擇的概率取值為：

(13)

式中：d(xi,xj)表示個體xi和xj的距離。

此外在算法進化中，變異操作往往使得新生成的煙花個體中某些維度的取值在規定的范圍之外，為此采用如下公式修復不可行解：

(14)

3 混合煙花算法

為了有效地求解UPMS-ETC問題的數學模型，本文基于標準FWA構建了HFWA。算法設計了雙層編碼方式和解碼方法，用于表示UPMS-ETC問題的解；融入了反向學習初始化方法，致力于提升初始解的質量；構建了基于變鄰域搜索算法的局部優化流程，用以強化基本FWA的尋優性能。

3.1 編碼與解碼

標準FWA采用實數編碼方式，無法直接用于求解離散組合優化問題。為此，本文結合問題的性質，構建了UPMS-ETC問題的編碼與解碼機制。

編碼過程：采用雙層實數編碼，兩層編碼長度值均為n；編碼第一層用于各個機器的加工任務序列，各維度編碼的取值范圍為[1,m+1)；編碼第二層用于確定各任務在機器上的加工速度，各維度編碼的取值范圍為[0,1]。

UPMS-ETC問題需要決策以下三點：(1) 各個機器所加工的任務集合；(2) 各個機器上加工任務的排序；(3) 各訂單的加工速度。相應的解碼過程歸納如下：

? 取第一層編碼的整數部分，用于確定各任務的加工機器編號，進而得到各個機器的加工任務集合。

? 基于各個機器的加工任務集合，取第一層編碼的小數部分用于確定各個機器上的任務加工順序，小數部分取值小所對應的加工任務優先加工。

? 針對各個機器，構建選擇加工速度的輪盤賭模型，并依據輪盤賭規則確定各個機器所分得的加工任務的加工速度。如圖1所示，假設某一機器共有三種不同的加工速度，則當該任務對應第二編碼的數值在區間[0,1/3)上，則采用加工速度1；若其數值在區間[1/3,2/3)上，則采用加工速度2；若其數值在區間[2/3,1]上，則采用加工速度3。

圖1 基于輪盤賭規則的加工速度選擇

為了清晰地說明該編解碼方法，給出如下算例：任務總數為8，機器總數為3，機器1和機器3具有3種加工速度，機器2具有兩種加工速度。已知如下編碼：第一層(3.19,1.80,2.97,2,90,1.47,1.35,2.45,3.78)，第二層(0.13,0.05,0.84,0.93,0.20,0.48,0.15,0.83)，圖2給出了完整的解碼過程。

(a)

(b)

(c)

3.2 反向學習初始化方法

基本FWA采用隨機的方法構建初始種群，一定程度上限制了算法的優化性能。為此，本文利用反向學習初始化方法構建初始解，力求提升初始種群的質量[10]。該方法的基本思路為：借助隨機方法生成多個初始解和相應的反向解，結合評價函數進行篩選，選取較優的解構造初始種群。結合以上描述，反向學習初始化方法的實現步驟歸納如下：

Step1令i←1，d←1，轉Step 2。

Step2若i≤N，轉Step 3；否則，轉Step 7。

Step3若d≤D，轉Step 4；否則，轉Step 5。

Step5令d←d+1，若d>D，轉Step 6,否則，轉Step 3。

Step6令i←i+1，轉Step 2。

3.3 基于變域搜索算法的局部優化流程

在算法迭代后期，基本FWA存在尋優性能不佳、易陷入局部最優的缺點。為此，本文基于變域搜索算法構建了局部優化流程，力求增強基本FWA的性能[11]。

變域搜索算法通過系統地改變當前解以拓展算法的搜索范圍，進而獲得待優化問題的局部最優解；同時，基于此局部最優解，再次系統地進行解空間的拓展，力求獲得另一局部最優解。圖3為變鄰域搜索算法的示意圖。

圖3 變鄰域搜索算法示意圖

鄰域結構設計構成了變鄰域搜索算法的一項核心內容，對于當前研究的PMS-CPT問題，本文采用交換和翻轉兩種鄰域結構實現從當前解到新解的變換，具體描述如下：

? 交換變異：隨機選擇當前解編碼的兩處位置，交換這兩處位置上的編碼數值生成新解。

? 翻轉變異：隨機選擇當前解編碼的兩處位置，翻轉這兩處位置間的編碼數值生成新解。

圖4所示的鄰域結構示意圖清晰地說明了以上兩種變異算子。當前編碼方法中，編碼第一層用于各個機器的加工任務序列，編碼第二層用于確定各任務在機器上的加工速度。因此，以上鄰域變換能夠對改善的解產生擾動，探索其周圍的解。

圖4 鄰域結構示意圖

由編解碼方法可知，編碼第一層用于各個機器的加工任務序列，編碼第二層用于確定各任務在機器上的加工速度。因此，以上鄰域變換能夠對當前解產生擾動，從而探索當前解附近的其他解。

3.4 混合煙花算法流程

Step1初始化算法各參數，利用反學習初始化方法生成初始解，并對每個煙花個體進行評價。

Step2依據評價函數計算當前種群中各個煙花的爆炸火花數目和爆炸半徑。

Step3利用均分分布和高斯分布生成爆炸火花和高斯變異火花。

Step4利用煙花、爆炸火花和高斯變異火花構建候選種群，并依據篩選規則生成新種群。

Step5對于各個煙花個體，利用基于變域搜索算法的局部優化流程進一步改善其性能。

Step6判斷是否滿足終止條件。若滿足則停止迭代搜索并輸出當前最優解；否則，轉Step 2。

4 仿真測試

4.1 實驗準備

在MATLAB 2015a平臺上進行仿真測試，設備參數為1.6 GHz、內存8 GB、Intel Core i5-8250U CPU，并利用CPLEX 12.4和CPLEX-MATLAB接口求解UPMS-ETC問題的MILP模型，CPLEX最大求解時間設置為3 600 s。考慮到難以從現有文獻中獲取與當前UPMS-ETC問題相似的測試算例，本文依據相關機器調度問題生成方法采用隨機方法構建測試算例集，具體包括小規模和中大規模兩類算例[12]。

在小規模算例中，機器數目m∈{2,3}，任務數目n∈{5,8,10,12,15}，共計10個算例；在中大規模算例中，機器數目m∈{2,3,5}，任務數目n∈{20,40,60,80,100}，共計15個算例。對于算例(m,n)，其他參數設置如下：

? 任務的標準加工時間lv在均分分布U(50,200)上隨機產生；

? 任務v單位時間的拖期成本μv的數值在均分分布U(0.1,1)上隨機產生；

? 加工速度類型總數Qi均設置為5，加工速度ris在均分分布U(1,20)上隨機產生；

? 各機器在不同速度下的能耗參數eis在均分分布U(0.1,2)上隨機產生，且對于同一機器而言，加工速度越快，其能耗越高；

? 任務v的交貨期dv在均分分布U(D,3·D)上隨機產生，參數D的計算公式為：

(15)

為了驗證算法的有效性，將本文提出的HFWA與多種優化算法進行比較，包括：FWA、粒子群算法(Particle swarm optimization, PSO)和蜂群算法(Artificial bee colony, ABC)。其中，與基本FWA進行對比用于加入改進方法后對標準算法性能的改善；與此同時，PSO與ABC算法作為目前較優秀的優化算法，其進化機理和HFWA存在些許差異，因而同這兩種算法進行對比能夠凸顯HFWA的競爭力。對于各個測試算例，以上四種算法均獨立運行15次，對測試結果進行統計分析，采用最優相對偏差(RPDb)和平均相對偏差(RPDm)作為評價指標[13]，公式定義如下：

(16)

(17)

4.2 參數校驗

表1 正交試驗參數設置

表2 正交試驗仿真結果

表3 正交試驗極差分析

由以上測試結果可知，局部搜索參數LS和高斯變異參數GN對HFWA的性能影響最大，其次是參數爆炸火花參數SN和基本爆炸半徑A，最后是種群規模N。據此可知這表明局部優化算子和高斯變異對算法的優化性能具有重要影響，因而需要合理設置LS和GN，若取值過大將使算法過分注重局部搜索，若取值過小則弱化了算法的挖掘能力。與此同時，其余三個參數也需要進行合理設置，才能使HFWA達到最優的搜索性能。各參數取值歸納如下：N=30，SN=10，GN=30，A=0.6，LS=10。

類似地，借助正交實驗進行小規模情況下的HFWA參數校驗，算法取值歸納如下：N=10，SN=5，GN=20，A=0.6，LS=5。為了確保對比的公平性，FWA各參數取值與FWFA相同，PSO和ABC的種群規模與HFWA相同，其余參數采用相似方法進行校驗。小規模情形下PSO和ABC的其他參數設置如下：PSO權重系數為1.0，全局學習參數取值為1.5，自身歷史最優學習參數2.0，ABC偵察蜂步數為10；中大規模情形下PSO和ABC的其他參數設置如下：PSO權重系數為0.95，全局學習參數取值為1.8，自身歷史最優學習參數2.0，ABC偵察蜂步數為20。

4.3 算法對比

結合以上參數校驗的結果，分別采用HFWA、FWA、PSO和ABC四種優化算法對各個小規模算例進行求解，測試結果見表4和表5。其中Minsol表示CPLEX所得目標函數的最優解，計算時長則為CPLEX求解各個算例的耗時。表格最后一行優勝率表示某一測試算的RPDm指標在所有測試算例中獲得最優性能的比率，例如“2/10”表示對應的算法在共計10組測試算例中兩次獲得最優RPDm值。此外，對于各個算例，以符號表示某一算法15次運行結果相對偏差百分比的標準差，圖5給出了各算法在不同算例中的δ取值。

表4 小規模算例測試結果1

表5 小規模算例測試結果2

圖5 小規模算例方差對比

由測試結果可知，在小規模算例中，各個測試算法的性能都較為良好，RPDb指標均能達到0，即能夠獲得小規模算例的精確解。就RPDm指標而言，HFWA的性能最優，其優勝率為10/10。由圖5可知，HFWA在多數情形下具有較小的δ值，表明其具有較強的魯棒性，算法優化結果更為穩定。此外，由CPLEX計算時長參數可知，隨著問題規模的擴大，求解MILP模型所需的時間急劇增加，算例(3,15)的計算時長達到2 670.20 s，因此有必要設計高效的智能算法來求解UPMS-ETC問題。

進一步地，利用HFWA、FWA、PSO和ABC四種測試算法求解各種大規模算例，測試結果見表6和表7。其中Minsol表示表示四種算法在共計60次仿真中所得最優解。對于各個算例，圖6給出了各算法在不同算例中的δ取值。

表6 中大規模算例測試結果1

續表6

表7 中大規模算例測試結果2

圖6 中大規模算例方差對比

由測試結果可知，四種算法中，HFWA表現值最優秀，其RPDb指標均為0，表明HFWA能夠獲得各個算例的已知最優解。就RPDm指標而言，HFWA的性能最優，其優勝率為12/15，表明本文算法在15組測試算例中12次獲勝；在剩余三組算例中，HFWA的表現同樣較為優秀，與最優RPDm值較為接近。此外，由圖6可知，HFWA在多數情形下具有較小的δ值，表明算法優化結果更為穩定，即算法具有較強的魯棒性。

5 結 語

本文以加工時間可控的UPMS問題為研究對象，構建了該問題的MILP模型，并提出了HFWA求解算法以最小化作業期間的能耗和延遲懲罰成本總和。算法設計中，創建了特定的編解碼方法以表示問題的解，同時融入了反向學習初始化方法以提升初始解的質量，并構建了基于變鄰域搜索算法的局部優化流程用以強化基本算法的尋優性能。在實驗部分，借助正交實驗對算法參數進行校正，并將HFWA與多種智能優化算法進行對比，結果表明本文構建的HFWA能夠快速、有效地求解UPMS-ETC問題，并具有顯著的優越性。