基于高維多目標候鳥優(yōu)化算法的柔性作業(yè)車間調度

王秋蓮 段星皓

南昌大學經濟管理學院，南昌，330031

0 引言

柔性作業(yè)車間調度問題(flexible job shop scheduling problem，F(xiàn)JSP)是指在一定的約束條件下，同時確定工件加工機器和加工順序，以達到理想的調度目標[1]。FJSP是NP-hard問題，大規(guī)模的FJSP很難通過數(shù)學規(guī)劃方法、拉格朗日松弛法等精確方法求得理論解，因此更多的研究采用啟發(fā)式算法對FJSP進行求解。如何改進啟發(fā)式算法使得其收斂性、穩(wěn)定性達到更優(yōu)是當前研究的熱點。

FJSP根據(jù)調度目標的數(shù)量分為單目標FJSP和多目標FJSP(multi-objective FJSP，MOFJSP)。其中，MOFJSP可進一步細分為低維(調度目標不超過3個)MOFJSP和高維MOFJSP。目前，大部分學者主要對低維MOFJSP進行研究，如孫麗珍等[2]以最小化最大完工時間、機器負載的均方差為目標，設計了一種改進解碼方法和初始化規(guī)則的遺傳算法；孟冠軍等[3]提出一種混合人工蜂群算法求解以最大完工時間、最大機器負荷和機器總負荷為目標的低維MOFJSP。實際生產中，制造企業(yè)不僅需要考慮上述目標，還需要考慮能源消耗。李進宇等[4]研究了基于遞歸分析和圖像處理技術的混流車間能耗。李明等[5]提出一種新型帝國競爭算法來求解以最小化完工時間、最大拖期、最大機器負荷和總能耗為目標的高維MOFJSP，但該方法對高維多目標解集個體的選擇壓力不足，可能導致算法的過早收斂。黃夏寶等[6]提出一種改進的NSGAⅢ算法，引入基于參考點的小生境選擇機制來提高算法的多樣性，求解以最小化最大完工時間、機器總負荷、最大機器負荷和機器能耗為目標的高維MOFJSP，但該算法的局部搜索能力還有待提高。由此可見，考慮能耗目標的高維MOFJSP研究是當前車間調度問題研究的熱點及難點，但處于起步階段[7]。

候鳥優(yōu)化(migrating birds optimization，MBO)算法[8]可用于求解二次分配問題，且所得的最終結果優(yōu)于遺傳算法、模擬退火算法等常用算法。QUAN等[9]將MBO算法用于求解混合流水車間的調度問題；任彩樂等[10]提出兩階段解碼方法及4種鄰域結構的改進MBO算法，求解混合流水車間調度問題。此外，有學者通過MBO算法來求解FJSP，如姚妮[11]以最小化最大完工時間為目標，結合MBO算法和變鄰域搜索策略來提高算法局部搜索能力；劉雪紅等[12]設計一種改進的MBO算法，引入精英分批策略和可行鄰域結構策略來解決以最小化最大完工時間和最大批次數(shù)目為目標的低維MOFJSP。

盡管這些學者證明MBO算法對FJSP能獲得高質量解集，但用MBO算法進行高維MOFJSP的研究較少，因此本文建立一種以最小化完工時間、機器總負荷、延期時間、生產總能耗為優(yōu)化目標的MOFJSP模型，在改進MBO算法的基礎上，引入Pareto思想和基于參考點的選擇算子，設計了一種高維多目標候鳥優(yōu)化(multi-objective migrating birds optimization，MOMBO)算法，通過2種選擇算子給予鳥群選擇壓力，以提高鳥群的多樣性和收斂性，并基于屬性層次模型和灰色關聯(lián)度分析法的組合權重決策選擇最優(yōu)解。最后通過算例和實驗驗證算法的有效性。

1 高維多目標柔性作業(yè)車間調度模型

FJSP可以描述為：工件集J={J1，J2，…，Jn}中的n個工件在機器集M={M1，M2，…，Mm}中的m臺機器上加工，其中，第i個工件Ji有ei個工序。MOFJSP主要需要解決兩類子問題：機器選擇和工序排序，其中，機器選擇是指從可選機器集Mij(i∈[1，n]，j∈[1，ei])中選擇一個機器完成工序Oij，工序排序是確定每個機器的加工順序。

工件加工過程中，有以下假設：①所有機器在0時刻均可開始加工工件，所有工件在0時刻均被釋放；②每臺機器一次只能完成一道工序，每道工序在同一時刻只能由一臺機器完成；③工件的加工需嚴格按照工序依次完成；④加工一旦開始就不能中止，直到完成；⑤機器完成最后一道加工就停止運行；⑥加工過程中，工件的轉移時間忽略不計。

本文提出的多目標柔性作業(yè)車間調度模型的調度目標的具體描述如下：

(1)最大完工時間。最大完工時間是指所有工件加工完成的時間，即最后一個工件加工的完工時間，其計算公式為

f1=max(Ci|i=1，2，…，n)

(1)

式中，Ci為工件i最后一道工序的完工時間。

(2)總拖期。客戶和企業(yè)會約定產品的交貨期，交貨時間超過交貨期，企業(yè)會損失信譽并被罰款?？偼掀诳梢苑从称髽I(yè)的信譽情況，其計算公式為

(2)

式中，Di為工件i的交貨期。

(3)機器總負荷。機器總負荷指所有加工機器的負荷總和，可以反映企業(yè)資源的利用效率，其計算公式為

(3)

式中，tijk為工序Oij在機器k上的加工時間；xijk為0-1變量，xijk=1表示工序Oij在機器k上加工，xijk=0表示Oij不在機器k上加工。

(4)總能耗?？偰芎闹饕ㄜ囬g固定能耗、機器空載能耗、機器加工能耗，其計算公式為

(4)

(5)

(6)

Eg=Pgmax(Ci|i=1，2，…，n)

(7)

2 高維多目標候鳥優(yōu)化算法

候鳥在飛行過程中會自然排成V字形隊列，這是因為領飛鳥通過翅膀扇動形成一股氣流，氣流向后方流動從而加大后方壓力，跟飛鳥利用壓力差飛行得更加輕松。MBO算法模仿候鳥的遷徙過程，將表現(xiàn)最優(yōu)的可行解作為領飛鳥，將其他可行解作為跟飛鳥。MBO算法主要分為鳥群初始化、領飛鳥進化、跟飛鳥進化和替換領飛鳥四個階段。本文在此基礎上，引入基于Pareto支配和基于參考點的選擇算子，對鳥群初始化、飛鳥進化階段進行改進，提出一種高維多目標候鳥優(yōu)化算法，該算法流程如圖1所示。

圖1 高維多目標候鳥優(yōu)化算法流程圖

2.1 基于Pareto支配和參考點的選擇規(guī)則

本文引入基于Pareto支配的選擇算子和基于參考點的選擇算子來增強對Pareto前沿面的選擇壓力，進而對解進行篩選和替換。

2.1.1基于Pareto支配的選擇算子

Pareto支配的定義為：對于解p和q，當且僅當解p的所有目標值都不大于解q的目標值且至少存在一個目標值小于解q的目標值，則稱p支配q(pq)。

定義si為鳥群中被個體i支配的個體集合，zi為鳥群中支配個體i的個體數(shù)量。

基于上述定義，提出一種快速非支配排序算法，具體步驟如下：

(1)將種群中不被其他個體支配的個體標記為第1級非支配前端F1；

(2)F1中的個體p支配的個體集合為sp，對于sp中的個體q，將支配q的個體數(shù)zq減1即zq←zq-1后，若zq=0即不存在可以支配q的其他個體，則將q保存至集合Q中。對F1的每個支配個體集合中的全部個體執(zhí)行上述操作后，將最終得到的集合Q作為鳥群的第2級非支配前端F2；

(3)遍歷F2中的個體，參照步驟(2)中的操作，找到第3級非支配前端F3。以此類推，直到種群中的每一個個體都分配到對應的非支配前端等級。

通過基于Pareto思想的快速非支配排序算法將鳥群分為不同等級，等級越低說明該個體越優(yōu)。

2.1.2基于參考點的選擇算子

基于參考點的選擇算子首先在目標空間內采用單層參考點生成法[13]得到均勻分布的多個參考點，即將A個目標均勻劃分成H份，然后根據(jù)參考點依附的個體數(shù)量選擇較優(yōu)個體，具體流程如下：

(8)

然后通過下式歸一化鳥群個體的目標值：

(9)

歸一化鳥群個體的目標值后，此時理想點變?yōu)樵c，然后將理想點與各參考點的連線作為參考線，計算歸一化后的目標點與所有參考線的距離，選擇最近的參考線作為個體依附對象；

(10)

最后依次遍歷鳥群個體，統(tǒng)計每個參考點的依附個體數(shù)量即該參考點的小生境數(shù)。在MOMBO算法鳥群進化操作時，對于同一級非支配前端中的個體集合，優(yōu)先選擇小生境數(shù)更少的個體，保證算法的多樣性。

2.2 多目標候鳥優(yōu)化算法的主要步驟

2.2.1編碼

運用MOMBO算法求解MOFJSP時首先要解決編碼問題，即將FJSP的潛在解轉換至可被MOMBO算法搜索的空間中。本算法采用二段式編碼規(guī)則[14]，將解分為機器選擇(machine selection，MS)和工序排序(operation sequence，OS)兩部分。MS按工件的工序順序依次排列，每個位置的數(shù)值為該道工序選擇的加工機器號；OS中，每個位置的數(shù)值為工件號，從左往右遍歷OS，同一工件號出現(xiàn)第j次數(shù)表示對應工件的第j道工序。MS和OS如圖2所示。

圖2 二段式編碼示意圖

2.2.2解碼

解碼則是將解空間轉換為問題空間。本算法采用貪婪插入式解碼法[15]，以保證每道工序在所選加工機器上盡可能早地開始，對于前插的工序，它在OS上的位置也要相應調整。

2.2.3鳥群初始化

初始化是車間調度問題的一個重要步驟，它生成解集的好壞對算法的收斂性和多樣性有很大的影響。為保證初始解的質量和多樣性，本文分別對機器選擇和工序排序采用多種初始化規(guī)則。

2.2.3.1 MS初始化規(guī)則

(1)改進的全局選擇規(guī)則[16]。設定一個初始值為0的機器負荷數(shù)組，從工件集合中隨機選擇一個工件，從該工件的工序中隨機選擇一個工序，將該工序的可選機器加工時間與機器負荷數(shù)組的已有負荷相加，得到臨時機器負荷數(shù)組。從臨時機器負荷數(shù)組中選擇臨時負荷最小的機器作為該工序的加工機器；若臨時負荷最小的機器不只一個，則選擇加工時間最短的機器作為該工序的加工機器；更新機器負荷數(shù)組。當前工件的所有工序選擇完加工機器后，隨機選擇下一個工件，依次類推，直到所有工件的工序都選擇完機器。

(2)改進的局部選擇規(guī)則[16]。設定一個初始值為0的機器負荷數(shù)組，從工件集合中依次選擇一個工件，從該工件的工序集合中隨機選擇一個工序，將該工序的可選機器加工時間與機器負荷數(shù)組的已有負荷進行相加，得到臨時機器負荷數(shù)組。從臨時機器負荷數(shù)組中選擇臨時負荷最小的機器作為該工序的加工機器；若臨時負荷最小的機器不只一個，則選擇加工時間最短的機器作為該工序的加工機器；更新機器負荷數(shù)組。當前工件的所有工序選擇完加工機器后，將機器負荷數(shù)組重新置0，選擇下一個工件，依次類推，直到所有工件的工序都選擇完機器。

(3)隨機選擇。在加工工序的可選機器中隨機選擇一臺機器進行加工。

2.2.3.2 OS初始化規(guī)則

(1)剩余負荷最大規(guī)則[17]。將工件按剩余工序加工時間總和進行降序排列，優(yōu)先加工剩余加工時間長的工件。

(2)最短加工時間規(guī)則[18]。將工序按照加工時間升序排列，優(yōu)先完成加工時間短的工序。

(3)隨機選擇。隨機選擇工序進行加工。

參考文獻[16，19]，首先按照鳥總數(shù)(即個體總數(shù))4∶4∶2的比例，分別通過改進的全局選擇規(guī)則、改進的局部選擇規(guī)則和隨機選擇規(guī)則生成初始鳥群的MS，再按照上述比例生成初始鳥群的OS。

初始化之后，按如下步驟將鳥群分為領飛鳥和左右跟飛鳥群：首先根據(jù)基于Pareto支配的選擇算子進行非支配排序，找到初始種群中的非支配解集；然后通過基于參考點的選擇算子，選擇小生境數(shù)最小的非支配解作為領飛鳥，若小生境數(shù)最小的非支配解有多個，則隨機選擇一個作為領飛鳥，并將剩下個體按非支配等級均勻分為左右跟飛鳥群。

2.2.4鳥群進化

首先定義三種鄰域結構：

N1：隨機產生i(i=2，3，4)個變異工序。將i個工序順序重新排列組合，然后計算不同排列組合的目標值，通過上文中的兩種選擇算子選擇其中最好的新解作為鄰域解。

N2：基于工序變異，交換兩個不同工件的任意工序位置。

N3：基于機器變異，隨機選擇兩個變異工件，并在這兩個變異工件上分別隨機選擇一道工序，在不考慮原加工機器的情況下，重新選擇這兩道工序加工時間最短的機器。

基于上述3種鄰域結構，增加鳥群個體的局部尋優(yōu)能力，鳥群進化具體步驟如下。

2.2.4.1 領飛鳥進化

按上述3種鄰域結構產生領飛鳥的3個鄰域解后，將領飛鳥與3個鄰域解合并作為進化解集，接著從該進化解集中選擇1個個體作為新的領飛鳥，具體的選擇規(guī)則如下：①若鄰域解支配領飛鳥，則替換領飛鳥；②若領飛鳥支配鄰域解，則不替換；③若領飛鳥與鄰域解互不支配，則進行領飛鳥選擇，具體的領飛鳥選擇操作如下。

首先將原鳥群中的進化個體剔除，計算剔除進化個體后的鳥群的小生境數(shù)；然后找出具有最小小生境數(shù)的參考點集J，遍歷J中的每個參考點。假設Ja是J中的一個參考點，判斷進化解集(包括進化個體、進化個體的鄰域解)的非支配解集中是否有個體依附于參考點Ja，若進化解集的非支配解集中沒有個體依附于參考點Ja，則將參考點Ja的小生境數(shù)ρa設為無窮大。若進化解集的非支配解集中有個體依附于參考點Ja，則進一步判斷ρa是否為0，若ρa為0，則從進化解集的非支配解集中依附Ja的個體中將與參考點Ja對應參考線距離最短的個體作為新的領飛鳥加入鳥群；若ρa不為0，則在進化解集的非支配解集中隨機選擇1個依附該參考點的個體作為新的領飛鳥加入鳥群。進化過程中選擇領飛鳥的流程如圖3所示。

圖3 進化過程中基于參考點的選擇算子流程圖

更新完領飛鳥之后，選擇一個次優(yōu)解作為第一個跟飛鳥的共享解，以增強跟飛鳥群的尋優(yōu)能力。次優(yōu)解的選擇過程如下：判斷進化解集中剩余個體的非支配等級，選擇非支配等級最低的個體作為次優(yōu)解；若非支配等級最低的個體數(shù)大于1，則比較其依附參考點在鳥群中的小生境數(shù)，優(yōu)先選擇小生境數(shù)更小的個體，若小生境數(shù)相同，則隨機選擇一個體作為次優(yōu)解。

2.2.4.2 跟飛鳥進化

跟飛鳥進化包括對左右跟飛鳥群的進化，即首先將領飛鳥進化解集中的次優(yōu)解作為左右跟飛鳥群中第一只跟飛鳥的共享解，將共享解與生成的3個鄰域解組成該跟飛鳥的進化解集。然后，選擇進化解集中的最優(yōu)解作為新的跟飛鳥，同時將進化解集中的次優(yōu)解作為下一只跟飛鳥的共享解，替換跟飛鳥和選擇次優(yōu)解規(guī)則與上文中領飛鳥的進化過程相同。按上述步驟依次遍歷所有跟飛鳥，完成跟飛鳥的進化。

2.2.5跟飛鳥交叉

每次對左右跟飛鳥群執(zhí)行完進化操作后，從左右跟飛鳥群中隨機選擇Cr(Cr為交叉?zhèn)€數(shù))對個體，并對每對個體依次進行交叉操作。

MS、OS的交叉操作分別采用均勻交叉策略和POX(precedence operation crossover)交叉策略[15]。MS的均勻交叉策略是：隨機產生一個長度為工序總數(shù)的0、1數(shù)組，交換父代數(shù)字1對應位的機器號，產生2個新的子代。OS的POX交叉策略是：首先將工件隨機地分為子集合J1和J2；然后將父代1中J1部分的工件號復制到子代1對應位置，將父代2中J2部分的工件號復制到子代2中對應位置；最后將父代2中J2部分的工件號按順序依次填入子代1中的剩余位，將父代1中J1部分的工件號按順序依次填入子代2中的剩余位。MS交叉如圖4所示，OS交叉如圖5所示。

圖4 MS均勻交叉圖

圖5 OS POX交叉圖

2.2.6篩選個體產生新鳥群

個體交叉生成子代鳥群Q后，將子代加入至原鳥群Y(數(shù)量為n)中，將Q和Y合并為R，對R重新進行非支配排序，在非支配解集中隨機選擇小生境數(shù)最小的個體作為新鳥群S的領飛鳥，同時更新S的小生境數(shù)；然后按照非支配層等級由低到高將每層個體加入至S中，直至S的數(shù)量等于n；若S的飛鳥數(shù)量大于n，則需要重新在最后一層非支配層中通過基于參考點的選擇算子選擇一定數(shù)量個體保留在S中，未被選擇個體則剔除出S，使S的數(shù)量等于n。新鳥群生成如圖6所示。

圖6 新鳥群生成圖

生成新鳥群之后，若未達到最大迭代次數(shù)，則返回至鳥群進化階段進行下一次迭代；若達到最大迭代次數(shù)，則將新鳥群的非支配解集作為最終優(yōu)化結果輸出。

2.2.7基于AHM和灰色關聯(lián)度分析法的組合權重決策

實際生產需要選擇一個確定的方案。MOMBO最終的非支配解集可能存在多個解，本節(jié)采用組合賦權法進行方案選擇。組合賦權是主觀賦權和客觀賦權的線性加權，其中，主觀賦權采用屬性層次模型，客觀賦權采用灰色關聯(lián)度分析法。

(1)屬性層次模型(attribute hierarchical model，AHM)先將各個指標劃分為相互關聯(lián)的有序結構，再結合專家的主觀意見和分析者的判斷來確定指標權重?；贏HM 的權重計算方法如下：

①建立判斷矩陣A=[aij]，比較同一等級指標之間的重要性，其中，aij=1/aji(i≠j)，aii=1。

②設μ=[μij]為測度矩陣，μij是由aij換算而成的測度指標，令

(11)

k=1，2，…

這里設置β=1。

③計算指標權重。指標主觀權重向量Wa=(w1，w2，…，wm)按下式計算：

(12)

其中，w1+w2+…+wm=1且0≤wi≤1。

(2)灰色關聯(lián)度分析(gray relation analysis，GRA)法的思想是：先根據(jù)某個問題的實際情況確定理想的最優(yōu)序列；然后，通過方案的序列曲線和幾何形狀與理想最優(yōu)序列的曲線和幾何形狀的相似程度來判斷方案的序列與理想最優(yōu)序列之間的關聯(lián)程度，曲線和幾何形狀越接近說明關聯(lián)度越大，方案越接近理想最優(yōu)，則該方案的權重越大。GRA的主要步驟如下：

①數(shù)據(jù)規(guī)范化處理。假設有n個方案，每個方案有m個目標值，首先對各個指標進行規(guī)范化處理，規(guī)范化處理計算式為

(13)

②由下式計算灰色關聯(lián)系數(shù)：

(14)

灰色關聯(lián)系數(shù)反映yij與理想值之間的關聯(lián)程度。

③權重計算。根據(jù)灰色關聯(lián)系數(shù)矩陣計算得到客觀權重向量Wb，其中，第j個目標的權重為

(15)

(16)

(3)綜合權重確定。將AHM得出的權重Wa與灰色關聯(lián)度分析法的權重Wb結合，求得綜合權重：

w*=θWa+(1-θ)Wb

(17)

其中，θ為主觀偏重系數(shù)，θ∈[0，1]，由決策者依據(jù)偏好給出。組合賦權法既能充分利用專家長期積累的生產加工經驗，使結果更貼近實際生產，又能在一定程度上從數(shù)據(jù)本身特性來確定其重要性。

3 性能測試與案例研究

為驗證本算法的有效性，首先將MOMBO分別與三個變體、NSGAⅡ及NSGAⅢ的最終解集進行對比，然后通過仿真對MOMBO求得的非支配解集進行分析。本節(jié)所有算法均由MATLAB R2014實現(xiàn)，MATLAB在Intel i5-8300H、內存8.00 GB的計算機上運行。

3.1 算法性能測試

本小節(jié)通過10個BRdata測試問題[20]來驗證MOMBO的性能，將測試問題的加工時間作為本調度問題中各個工序的基本加工時間，算例的其他數(shù)據(jù)如交貨期、機器加工功率、空載功率等數(shù)據(jù)如表1所示。

表1 算例數(shù)據(jù)

多目標優(yōu)化算法的解集有多種評價指標，本文采用的評價指標是反向世代距離(inverted generational distance，IGD)[21]和超體積(hypervolume，HV)[22]。

IGD能綜合評價解集的收斂性和多樣性，其計算公式為

(18)

其中，P*為一組理想解集，本文設定為所有算法多次運算后得到的全部解中的非支配解集；|P*|為解集P*中解的個數(shù)；x為解集P*中的某個解；解集P是某算法求出的該前沿面的一個近似解集，y是解集P中的某個解；d(x，y)為解x和y在目標空間中的歐氏距離。IGD值越小，算法性能越好。

超體積能綜合反映解集的收斂性和多樣性，其計算公式為

(19)

其中，解集P是Pareto前沿面的近似解集；fm為解集P中某個非支配解f=(f1，f2，…，fm)T的第m個目標值；r=(r1，r2，…，rm)T為目標空間中的一個參考向量，它被解集P中所有目標向量支配；volume(*)表示近似Pareto前沿面與參考向量所圍成的超立方體的超體積。那么關于參考向量r的超體積指的是被解集P所支配且以參考向量r為邊界的目標空間的體積。

將超體積定義為被解集P所支配且以參考點r為邊界的目標空間的體積，HV(P，r)越大，算法性能越好。對于4個目標的柔性作業(yè)車間調度問題，本文采用BADER等[23]提出的蒙特卡羅模擬方法計算HV(P，r)。

3.1.1MOMBO與MOMBO變體的比較

為驗證所做改進對MOMBO算法性能的影響，構建以下幾種變體類型：

MOMBO1：采用半主動解碼。

MBMBO2：種群初始化階段采用隨機初始化方式。

MOMBO3：飛鳥進化過程僅采用N1鄰域結構產生鄰域解。

MOMBO算法參數(shù)設置如表2所示。多次運行MOMBO及其3種變體，求得多個近似解集，然后計算IGD值和HV值，最終取平均值進行對比，結果如表3所示，其中，最優(yōu)的結果用加粗數(shù)字表示。由表3可知，除了算例MK04和MK05，MOMBO算法的IGD值、HV值都優(yōu)于3種變體算法，由此可知基于MBO算法提出的貪婪插入式解碼法、改進的多種初始化規(guī)則，以及改進的鄰域搜索策略在大部分情況下更適用于高維多目標柔性作業(yè)車間調度模型。

表2 MOMBO參數(shù)設置

表3 MOMBO及其變體的IGD、HV

3.1.2MOMBO與NSGAⅡ、NSGAⅢ的比較

NSGAⅡ、NSGAⅢ的參數(shù)設置如表4所示。多次運行MOMBO、NSGAⅡ和NSGAⅢ，求得多個近似解集，然后計算IGD值和HV值，最終取平均值進行對比，結果如表5所示，其中，最優(yōu)的結果用加粗數(shù)字表示。由表5可知，MOMBO的IGD值和HV值都顯著優(yōu)于NSGAⅡ和NSGAⅢ，驗證了MOMBO求解該問題的有效性。

表4 NSGAⅡ、NSGAⅢ的參數(shù)設置

表5 MOMBO、NSGAⅡ和NSGAⅢ的IGD、HV

3.1.3算法收斂性分析

為驗證算法運行過程中的收斂性，繪制MOMBO、3種變體、NSGAⅡ和NSGAⅢ的HV值隨迭代次數(shù)變化的進化軌跡圖。由圖7可以看出，6種算法的HV值都隨著迭代次數(shù)的增大而增大，這意味著它們的收斂性能穩(wěn)定；MOMBO每代的HV值在MK08、MK09算例上遠大于其他算法的HV值，說明MOMBO的進化效果最好；MOMBO2、NSGAⅡ、NSGAⅢ的HV值明顯小于MOMBO、MOMBO1、MOMBO的HV值，這說明本文提出的多種初始化方式是有效的。接下來考慮非支配解集占比情況，MOMBO求解MK01的初始解集中，非支配等級共有5層，其中，非支配解集占所有解集的比例為28.6%；MOMBO求解MK01的最終解集的非支配等級為4層，其中，非支配解集占所有解集的68.6%。這說明所提出的兩種選擇算子能給予算法有效的進化壓力，能有效促進算法的收斂。

(a)MK04 (b)MK08 (c)MK09

3.2 實例研究

為驗證本文算法在實際生產加工中的效果，在某工廠進行加工實驗。實驗內容是，5臺加工機器上加工4個工件，通過Precision Power Analyzer LMG671收集機床加工工件過程中的功率。對得到的數(shù)據(jù)進行預處理，獲得每臺機器的加工功率和空載功率曲線圖。為方便計算，通過取均值求得各機器的加工功率和空載功率。車間固有功率為2kW，每個工件的交貨期和加工時間如表6所示，其中，“-”表示該工序不能在此機器上完成。

表6 實驗數(shù)據(jù)

3.2.1MOMBO與單目標MBO的比較

通過MOMBO算法求得Pareto前沿解集，如表7所示。通過MBO算法分別求解出以最大完工時間f1、總拖期f2、機器最大負荷f3、總能耗f4為目標的單目標最優(yōu)解，如表8所示。觀察表8可以發(fā)現(xiàn)，雖然單目標最優(yōu)解在某個目標值上最優(yōu)，但其他3個目標值較差，單目標最優(yōu)解與通過MOMBO算法求出的最優(yōu)解集中的某個解的目標值之差即Δf1～Δf4如表9所示。由表9可見，在MOMBO求出的最優(yōu)解集中，總能找出優(yōu)于單目標最優(yōu)解的解，因此通過MOMBO求出的最優(yōu)解集相對于單目標最優(yōu)解具有更好的質量。

表7 實例Pareto前沿解集

表8 實例單目標最優(yōu)解

表9 單目標最優(yōu)解與多目標最優(yōu)解的目標值之差

3.2.2MOMBO與NSGAⅡ、NSGAⅢ的比較

NSGAⅡ、NSGAⅢ求得的最優(yōu)解集分別如表10、表11所示。對比表7可以看出，MOMBO求出的最優(yōu)解大部分都優(yōu)于NSGAⅡ、NSGAⅢ求得的解。為進一步對比MOMBO與NSGAⅡ、NSGAⅢ的收斂性和多樣性，通過式(18)、式(19)求得3種算法的IGD值和HV值，如表12所示。在求解過程中，本文設定P*為三種算法多次運算后得到的非支配解集。如表12所示，對于本實例，MOMBO取得了優(yōu)于另外2個算法的結果。

表10 NSGAⅡ的實例Pareto前沿解集

表11 NSGAⅢ的實例Pareto前沿解集

表12 MOMBO、NSGAⅡ和NSGAⅢ求解實例的結果

圖8、圖9所示為3種算法求解實例的HV值和IGD值隨迭代次數(shù)的進化曲線。3種算法的HV值都隨著迭代次數(shù)的增大而增大，IGD值均隨著迭代次數(shù)的增大而減小，這說明它們的收斂性能穩(wěn)定。本文提出的MOMBO的HV值大于另外兩種算法的HV值，IGD值小于另外兩種算法的IGD值，因此MOMBO的進化效果最好。

圖8 3種算法求解實例所得HV值隨迭代次數(shù)的進化軌跡

圖9 3種算法求解實例所得IGD值隨迭代次數(shù)的進化軌跡

3.2.3基于組合權重法的最優(yōu)方案決策

首先根據(jù)專家對不同目標的打分情況建立AHM需要的判斷矩陣，然后通過計算得到各個目標的主觀權重向量Wa=(0.264，0.361，0.264，0.111)；接著根據(jù)表7中的Pareto前沿解集計算出灰色關聯(lián)度分析法的各個數(shù)據(jù)：

Wb=(0.261，0.189，0.281，0.269)

進一步取θ=0.5，求得綜合權重w*=(0.263，0.275，0.272，0.190)。由于本文是逆指標，故需要對y采用減法一致法處理得到y(tǒng)′，即利用每一列的最大值Mj依次減去該列的原始數(shù)據(jù)，然后將y′×w*，得到最終得分矩陣：

由R可知第二組方案的得分最高，其加工時間為53.42 s，拖期為0 s，機器負載為178.50 s，總能耗為364 681 J。圖10為按第二個調度方案加工的甘特圖。

圖10 解2的甘特圖

4 結論

本文提出以最小化最大完工時間、總拖期、機器總負荷、總能耗為目標的高維多目標柔性作業(yè)車間調度模型，并運用一種MOMBO進行求解。算法性能測試表明MOMBO相對于其3種變體以及NSGAⅡ、NSGAⅢ來說，其最優(yōu)解集在IGD和HV兩種性能指標上表現(xiàn)更優(yōu)。實例研究證明本文提出的MOMBO算法相對于單目標MBO算法、NSGAⅡ和NSGAⅢ有著更好的解集質量，在實際生產加工中能夠給予決策者更好的選擇。