陳揆能，袁小芳

（1.復(fù)雜環(huán)境特種機(jī)器人控制技術(shù)與裝備湖南省工程研究中心（湖南理工職業(yè)技術(shù)學(xué)院），湖南湘潭 411104；2.湖南大學(xué)電氣與信息工程學(xué)院，長沙 410082）

0 引言

開放車間調(diào)度問題（Open-shop Scheduling Problem，OSP）作為智能制造領(lǐng)域的重點(diǎn)研究問題，廣泛存在于檢測、維修、加工等生產(chǎn)制造行業(yè)［1］，對該問題的深入研究有重要理論和應(yīng)用價值。目前，大多數(shù)OSP 相關(guān)文獻(xiàn)只考慮了固定的工序加工時間。然而，隨著加工制造技術(shù)的發(fā)展，以數(shù)控加工（Computer Numerical Control，CNC）機(jī)床為代表的先進(jìn)加工設(shè)備可以通過調(diào)節(jié)加工功率控制加工時間［2-3］，并且縮短加工時間通常需要更高的加工能耗。因此，同時優(yōu)化完工時間和總加工能耗對高效、節(jié)能的開放車間生產(chǎn)具有重要意義。為此，本文研究了可控加工時間的多目標(biāo)開放車間調(diào)度問題（Multi-Objective Open-shop Scheduling Problem with Controllable Processing Time，MOOSPCPT）問題。

OSP 屬于NP-hard 問題［4］，精確算法很難在合理的時間內(nèi)找到滿足精度要求的近似解。因此，近年來研究人員多采用群體進(jìn)化算法對相關(guān)問題進(jìn)行求解。Rahmani Hosseinabadi 等［5］和Shareh 等［6］分別提出了一種擴(kuò)展遺傳算法和改進(jìn)蝙蝠優(yōu)化算法求解OSP。Abreu 等［7］針對連續(xù)加工工序之間存在設(shè)置時間的OSP 進(jìn)行了研究并提出一種混合遺傳算法對問題進(jìn)行求解。Mejía 等［8］和Zhuang 等［9］對考慮工序設(shè)置時間和機(jī)器運(yùn)輸時間的OSP 進(jìn)行了研究。Sheikhalishahi 等［10］針對考慮人為失誤與機(jī)器維修生產(chǎn)因素在內(nèi)的多目標(biāo)OSP 進(jìn)行了研究。然而，目前對OSP 的研究多集中在理論層面，提出的基本算法框架通常無法完全適用于實(shí)際問題的求解。因此，依據(jù)問題特性設(shè)計(jì)改進(jìn)的混合進(jìn)化算法成為行業(yè)內(nèi)研究的熱點(diǎn)問題之一，開發(fā)高效的全局搜索算子和局部搜索算子是研究的核心內(nèi)容。

生物地理學(xué)優(yōu)化（Biogeography-Based Optimization，BBO）算法是由Simon［11］提出的基于種群的新型元啟發(fā)式方法，常用于全局搜索。BBO 算法通過遷移策略實(shí)現(xiàn)不同解之間的信息交流，通過變異策略實(shí)現(xiàn)解的內(nèi)部進(jìn)化，具有參數(shù)少、收斂速度快的優(yōu)點(diǎn)［12］，廣泛應(yīng)用于車間調(diào)度等NP-hard問題求解。張國輝等［13］將BBO 算法應(yīng)用于柔性作業(yè)車間調(diào)度問題的求解。Wang 等［14］以最小化總延遲時間為目標(biāo)，提出了一種混合BBO 算法用于解決流水車間調(diào)度問題。吳定會等［15］將BBO 算法應(yīng)用于求解多目標(biāo)柔性車間調(diào)度問題并取得了良好的效果。然而，由于BBO 中解的遷移概率與解之間的優(yōu)劣關(guān)系直接相關(guān)，容易使種群多樣性降低，導(dǎo)致算法陷入局部最優(yōu)。為此，本文改進(jìn)了傳統(tǒng)的遷移策略和變異策略用于平衡算法的開發(fā)能力與探索能力。

變鄰域搜索（Variable Neighborhood Search，VNS）通過系統(tǒng)化地移動當(dāng)前解到其鄰域解空間中的“鄰居”位置以實(shí)現(xiàn)局部搜索。變鄰域搜索顯著的特點(diǎn)之一是可以很好地結(jié)合問題知識進(jìn)行鄰域結(jié)構(gòu)的設(shè)計(jì)，從而實(shí)現(xiàn)對當(dāng)前解的規(guī)則化的移動，降低搜索的隨機(jī)性，極大地提高搜索的效率和精度。關(guān)鍵路徑是生產(chǎn)調(diào)度問題重要的知識，因此，基于關(guān)鍵路徑的變鄰域搜索得到了廣泛的研究。Xia 等［16］針對動態(tài)集成工藝規(guī)劃與調(diào)度問題設(shè)計(jì)了一種混合遺傳算法，并使用基于關(guān)鍵路徑的變鄰域搜索算法增強(qiáng)了算法的搜索能力。Yuan等［17］在差分進(jìn)化算法中加入了基于關(guān)鍵路徑的鄰域搜索策略，并對柔性作業(yè)車間調(diào)度問題進(jìn)行求解。趙詩奎［18］開發(fā)了一種基于關(guān)鍵路徑的兩級鄰域搜索算法，該算法用于求解時間最小化的柔性作業(yè)車間調(diào)度問題。目前，應(yīng)用變鄰域搜索算法求解開放車間調(diào)度問題的研究相對匱乏。為此，本文設(shè)計(jì)了一種基于關(guān)鍵路徑的自調(diào)整變鄰域搜索用于提高算法的局部搜索性能。

本文的主要工作包括：1）在開放車間環(huán)境下考慮可控的加工時間，同時優(yōu)化了完工時間和總額外能耗，既有效地提高了生產(chǎn)效率又滿足當(dāng)下綠色制造的重大需求；2）設(shè)計(jì)了高效的多目標(biāo)混合進(jìn)化算法（Multi-Objective Hybrid Evolutionary Algorithm，MOHEA）求解MOOSPCPT，在MOHEA 中提出了一系列改進(jìn)策略，豐富了生產(chǎn)調(diào)度優(yōu)化問題的求解方法。

1 問題描述與模型

本文提出的MOOSPCPT 可以描述為：有n個待加工工件，第i個工件為Ji（i=1，2，…，n）；共在m臺機(jī)器上加工，第j臺機(jī)器為Mj（j=1，2，…，m）；每個工件包含與機(jī)器數(shù)量相等的工序，且每個工序唯一對應(yīng)一個機(jī)器進(jìn)行加工；Oij表示Ji由Mj負(fù)責(zé)加工的工序。同工件的工序之間不存在優(yōu)先關(guān)系。工序的加工時間可以通過調(diào)節(jié)機(jī)器的輸出功率控制，對于同一臺機(jī)器加工同一個工序而言，較短的加工時間對應(yīng)較高的加工能耗，以最長加工時間對應(yīng)的能耗為基礎(chǔ)能耗，則實(shí)際能耗與基礎(chǔ)能耗的差值為額外能耗。

因此，MOOSPCPT 包含2 個子問題：安排最優(yōu)的工序加工順序，確定各工序最優(yōu)的加工時間。本文通過解決上述兩個子問題同時優(yōu)化完工時間和總額外能耗。模型中所用到的符號及其含義如表1 所示。

MOOSPCPT 的目標(biāo)函數(shù)為：

其中：式（1）為目標(biāo)函數(shù)——最小化完工時間；式（2）為目標(biāo)函數(shù)——最小化總額外能耗；式（3）表示Oij的實(shí)際加工時間在取值范圍內(nèi)；式（4）表示加工一旦開始就不允許中斷；式（5）表示同一個工件不能同時在兩臺機(jī)器上加工；式（6）表示同一臺機(jī)器在同一時刻只能加工一個工件。

2 多目標(biāo)優(yōu)化問題及評價指標(biāo)

與單目標(biāo)決策可以獲得單一最優(yōu)解不同，多目標(biāo)優(yōu)化問題中各個目標(biāo)函數(shù)之間通常是相互沖突的，即一個目標(biāo)的改進(jìn)會導(dǎo)致另一個目標(biāo)的退化。一般多目標(biāo)問題可以表示為：

如果解向量X不被其他任何解向量所支配，則稱X為Pareto 最優(yōu)解，一組Pareto 最優(yōu)解稱為Pareto 最優(yōu)解集，問題的所有Pareto 最優(yōu)解稱為Pareto 前沿。多目標(biāo)優(yōu)化問題的目標(biāo)在于盡可能多地找到接近真實(shí)Pareto 前沿的解。

本文使用三個通用的性能指標(biāo)定量衡量算法獲得的Pareto 前沿的優(yōu)劣：

1）GD（Generational Distance）：用于表征算法的收斂性。GD 值越小，算法的收斂性越好，其定義為：

其中：N為Pareto 前沿中解的個數(shù)；di為Pareto 前沿中第i個解與真實(shí)的Pareto 前沿中距離其最近解之間的歐氏距離。

2）IGD（Inverse Generational Distance）：用于表征算法包括收斂性和分散性在內(nèi)的綜合性能。IGD 越小，算法的綜合性能越好，其定義為：

其中：N為真實(shí)Pareto 前沿中解的個數(shù)；Dj為真實(shí)Pareto 前沿中第j個解與算法得到的Pareto 前沿中與其距離最近的解之間的歐氏距離。

3）Δ（Spread）：用于表征Pareto 前沿的分布性與多樣性，Δ越小，代表Pareto 前沿的分布性與多樣性越好。Δ定義為：

其中：N是Pareto 前沿中解的個數(shù)；M為目標(biāo)函數(shù)個數(shù)；di、代表Pareto 前沿中第i個解與最近解之間的歐幾里得距離，dˉ代表所有di的平均值；dl為Pareto 前沿中第l個目標(biāo)上的極值解與Pareto 前沿邊界解的歐氏距離。

由于MOOSPCPT 的真實(shí)Pareto 前沿未知，本文將所有單獨(dú)求得的Pareto 前沿合并，得到的新Pareto 前沿視為真實(shí)Pareto 前沿。

3 MOHEA求解MOOSPCPT

3.1 MOHEA的基本流程

本文所提出的MOHEA 的基本流程如下：

輸入?yún)?shù)：種群數(shù)量n，最大遷入率I，最大遷出率E，最大變異率mmax，最大迭代數(shù)gmax，變鄰域搜索初始概率pv，精英庫大小nA。

步驟1 初始化：隨機(jī)產(chǎn)生初始種群。

步驟2 全局搜索：對種群進(jìn)行快速非支配排序與擁擠度排序，計(jì)算得到棲息地的遷入遷出率。對種群進(jìn)行遷移變異操作產(chǎn)生新解，如果新解優(yōu)于原來的解則用新解替代舊解。

步驟3 局部搜索：對種群執(zhí)行基于關(guān)鍵路徑的自調(diào)整變鄰域搜索操作，并更新變鄰域搜索概率pv。

步驟4 時間重置并更新精英庫：對當(dāng)前種群Pareto 前沿解進(jìn)行時間重置。將得到的解與當(dāng)前精英庫中的解合并，并保留新的Pareto 前沿；

步驟5 終止判斷：判斷是否達(dá)到最大迭代次數(shù)，是則結(jié)束算法，輸出精英庫；否則轉(zhuǎn)步驟2。

3.2 編碼與解碼

與普通的開放車間調(diào)度問題相比，本文需要額外考慮工序的加工時間。因此，本文采用的編碼由加工順序（Processing Sequence，PS）和工序時間（Processing Time，PT）兩部分組成。PS 為1 到n×m的一個隨機(jī)排列，PS 中的元素代表元素中的數(shù)值與固定編號對應(yīng)的工序；PT 中的元素表示與固定編號對應(yīng)的工序的加工時間。下面通過一個示例展示編碼過程。表2 為一個3×3（n×m）的算例，針對該算例的一個個體編碼如圖1 所示。PS 中包含一個1～9 的隨機(jī)排列。PS 第1 個元素值為6，代表編號6 對應(yīng)的工序O23第1 個加工；同理，如PS 第5 個元素值為9，代表編號9 對應(yīng)的工序O33第5 個加工。PT 第1 個元素為30，代表編號1 對應(yīng)的工序O11的加工時間為30 個單位時間；同理，如第9 個元素代表編號9 對應(yīng)的工序O33的加工時間為28 個單位時間。

表2 示例的加工信息Tab.2 Processing information of example

圖1 個體編碼示意圖Fig.1 Schematic diagram of individual encoding

本文采用常用的主動調(diào)度［19］解碼方式，基本過程如下：按照從左至右的順序遍歷PS 對應(yīng)的工件；在不推遲其他當(dāng)前已經(jīng)安排的工序的前提下，將工序安排在其加工機(jī)器上盡可能早的時間開始加工；該工件的加工時間為PT 序列中對應(yīng)的時間。

為了更加清晰地描述解碼過程，對圖1 所示編碼進(jìn)行解碼。對于同一工序在同一機(jī)器上加工，其加工時間越短額外加工能耗越高。因此，本文使用反比例函數(shù)刻畫機(jī)器額外加工能耗隨工序加工時間變化的特性，加工Oij的額外加工能耗由式（12）得到。圖2 為解碼獲得的甘特圖?？梢缘玫狡渫旯r間為154 個單位時間，總額外能耗為39 個單位能耗。

圖2 個體解碼甘特圖Fig.2 Gantt chart of individual decoding

3.3 基于BBO的改進(jìn)遷移策略和變異策略

BBO 通過模擬物種在不同棲息地間的遷移行為對種群進(jìn)行優(yōu)化，以棲息地作為個體進(jìn)行操作。BBO 算法主要包括遷移策略和變異策略。

遷移策略模擬了物種在棲息地之間的遷移行為，物種生存的適宜程度決定了棲息地間的物種遷移概率：適宜程度高的棲息地物種數(shù)量趨于飽和，物種遷出率高而遷入率低；適宜程度低的棲息地物種數(shù)量較少，物種遷入率高而遷出率低。如圖3 所示為單個棲息地的物種線性遷移模型，其中：λ為棲息地的遷入率，最大值為I；μ為棲息地的遷出率，最大值為E，S為棲息地物種適宜度；Smax為棲息地物種適宜度最大值。當(dāng)棲息地遷入率與遷出率相等時，棲息地種群數(shù)量達(dá)到動態(tài)平衡點(diǎn)S0。

變異策略是BBO 算法中另一個重要的搜索策略。變異過程模擬了自然界中由于棲息地環(huán)境發(fā)生改變引起的棲息地內(nèi)部的物種進(jìn)化，用于對個體產(chǎn)生小范圍的擾動。積極的變異（如：新物種的出現(xiàn)、附近物種的遷入等）能夠有效提升種群的多樣性，而消極事件（如：疾病、貪婪捕食者入侵等）引起的突變將降低種群多樣性。因此，在設(shè)計(jì)變異策略時應(yīng)該更多地引入積極的變異并盡量降低消極事件對種群的影響。圖3 為線性遷移模型。當(dāng)棲息地遷入率等于遷出率時，該棲息地物種種群數(shù)量達(dá)到動態(tài)平衡點(diǎn)S0，遷入該棲息地的種群數(shù)量等于遷出該棲息地的種群數(shù)量；當(dāng)棲息地的環(huán)境變化平衡點(diǎn)發(fā)生偏移（S0增大為正偏移，減小為負(fù)偏移）：S0增大是由于其他種群突然遷入，例如大量生物從附近棲息地遷入或由于基因突變導(dǎo)致大量的新物種出現(xiàn)；S0減小是由于疾病、貪婪捕食者入侵或一些自然大災(zāi)害導(dǎo)致的棲息地種群數(shù)急劇減少。

圖3 線性遷移模型Fig.3 Linear migration model

本文采用上述BBO 算法的遷移策略對PS 編碼進(jìn)行遷移操作，采用變異策略對PS 編碼和PT 編碼進(jìn)行變異操作。其中，每個個體代表一個棲息地，個體優(yōu)劣的排名代表?xiàng)⒌匚锓N適應(yīng)度，遷移概率用于選擇個體進(jìn)行遷移，變異概率用于控制個體是否進(jìn)行變異。同時，針對BBO 算法易陷入局部最優(yōu)的缺點(diǎn)，本文對遷移概率和變異概率進(jìn)行了改進(jìn)。

為適應(yīng)多目標(biāo)問題的求解，本文在計(jì)算棲息地的遷入遷出率之前，先對棲息地進(jìn)行非支配關(guān)系和擁擠度排序［20］，并根據(jù)計(jì)算結(jié)果將種群按照優(yōu)劣關(guān)系進(jìn)行排序。其中，最劣的個體排序?yàn)?，最優(yōu)的個體排序?yàn)榉N群數(shù)n。以排序?yàn)閭€體i的適應(yīng)度Si，則Smax=n。

3.3.1 改進(jìn)遷移策略

本文對傳統(tǒng)遷移策略進(jìn)行改進(jìn)：在種群進(jìn)化早期將遷移率設(shè)為常數(shù)保證種群的多樣性，從而獲得較好的全局搜索性能；在進(jìn)化后期遷移率與個體適應(yīng)度的優(yōu)劣程度相關(guān)，使算法能夠進(jìn)行充分的局部搜索，從而獲得更加精確的解。遷入率和遷出率隨迭代次數(shù)的變化如式（13）～（14）所示：

其中：λi和μi分別為Si個體的遷入率和遷出率；c1和c2分別為種群進(jìn)化早期遷入率和遷出率常數(shù)（c1，c2∈[0，1]）；iter為當(dāng)前迭代次數(shù)；T為常數(shù)。T值較小時會獲得較快的收斂速度但全局搜索能力較差，且T與c1、c2之間相互影響，為平衡遷移策略的全局搜索能力和局部搜索能力，實(shí)際應(yīng)用時通過正交實(shí)驗(yàn)確定c1、c2和T的具體取值。

遷移策略的偽碼如下所示：

第10）行所述的遷移操作過程如圖4 所示：1）隨機(jī)產(chǎn)生兩個遷移位點(diǎn)將遷入個體Hi和遷出個體Hj的PS 編碼分成3段；2）將遷入個體Hi的中間一段遷入新個體Hk中；3）遷出個體Hj中第1 段和第3 段的元素按照遷出個體Hj中對應(yīng)元素的順序填入Hk中的剩余位置。

圖4 遷移過程示意圖Fig.4 Schematic diagram of migration process

3.3.2 改進(jìn)變異策略

變異策略算法對防止過早收斂具有重要作用。本文對固定變異概率的方式進(jìn)行了改進(jìn)：按照個體的優(yōu)劣采用不同的變異概率，對于特別優(yōu)或者特別劣的個體變異概率大，從而實(shí)現(xiàn)較大范圍的局部搜索。Si個體的變異概率定義為：

其中：mmax為最大變異率。

本文通過比較隨機(jī)數(shù)r（0＜r＜1）與Si個體的變異概率mi控制Si個體是否進(jìn)行變異：若r＜mi，則對Si個體執(zhí)行變異操作。變異操作包括PS 編碼變異和PT 編碼變異。其中，PS 編碼變異：隨機(jī)選擇PS 編碼中的兩個元素進(jìn)行互換；PT 編碼變異：在PT 編碼中隨機(jī)選擇一個元素，在該元素對應(yīng)的工序加工時間取值范圍內(nèi)隨機(jī)選擇一個與當(dāng)前加工時間不相等的值。

3.4 基于關(guān)鍵路徑的自調(diào)整變鄰域搜索

變鄰域搜索能快速改善解的質(zhì)量，加快算法收斂速度。然而，在算法迭代過程中全局搜索的進(jìn)度往往難以匹配鄰域搜索的進(jìn)度，導(dǎo)致變鄰域搜索很難繼續(xù)對解進(jìn)行有效的改進(jìn)，反而由于無效的鄰域搜索耗費(fèi)了大量的時間，影響了算法最終的性能。為了解決這一問題，本文提出了一個自調(diào)整選擇策略：定義一個變鄰域搜索概率pv，當(dāng)產(chǎn)生的隨機(jī)數(shù)小于pv時，對當(dāng)前種群執(zhí)行變鄰域搜索；同時，變鄰域搜索算法可以根據(jù)對鄰域搜索結(jié)果的分析，自動地調(diào)整pv的值，即在鄰域搜索對解改進(jìn)大時增大pv的值，在鄰域搜索對解改進(jìn)小時減小pv的值。

自調(diào)整pv的過程為：設(shè)定變鄰域搜索初始概率pv（本文設(shè)為0.4）；當(dāng)前種群的Pareto 前沿PX包含x個非支配解，經(jīng)過變鄰域搜索操作得到新的Pareto 前沿PY包含y個非支配解，定義計(jì)數(shù)器c=0；對于PY中的每一個解，如果在PX中找到被其支配的解，則c=c+1。定義變量Ra=c/（xy），則下一次迭代時pv由式（16）計(jì)算得到：

關(guān)鍵路徑是在變鄰域搜索中常用的問題知識。析取圖中從起點(diǎn)到終點(diǎn)的最長路徑稱為關(guān)鍵路徑［21］，關(guān)鍵路徑上的工序稱為關(guān)鍵工序，同一機(jī)器上相鄰工序組成工序塊，工序塊中第一道工序稱為塊首工序，最后一道工序稱為塊尾工序，其余為塊內(nèi)工序。研究表明：對關(guān)鍵路徑上的工序塊進(jìn)行鄰域搜索是改善解質(zhì)量的最有效方式之一［22］。因此，本文采用了基于關(guān)鍵路徑的三種鄰域結(jié)構(gòu)［23］進(jìn)行鄰域搜索。對于關(guān)鍵工序塊L=｛l1，l2，…，ln-1，ln｝［23］，這三種鄰域結(jié)構(gòu)包括：

1）對塊首工序l1與塊尾工序ln，將其插入工序塊L的內(nèi)部；

2）對塊內(nèi)工序｛l2…ln-1｝，將其插入塊首工序l1之前或塊尾工序ln之后；

3）交換工序塊前兩道工序｛l1，l2｝或后兩道工序｛ln-1，ln｝。

綜上所述，本文所提出的基于關(guān)鍵路徑的自調(diào)整變鄰域搜索步驟如下：

步驟1 生成隨機(jī)數(shù)r，如果r＜pv，則進(jìn)行變鄰域搜索，否則結(jié)束變鄰域搜索；

步驟2 確定選定解的關(guān)鍵路徑及關(guān)鍵工序塊；

步驟3 隨機(jī)選擇一個未被選中過的關(guān)鍵工序，選取與該關(guān)鍵工序?qū)?yīng)的鄰域結(jié)構(gòu)對其進(jìn)行鄰域搜索操作，產(chǎn)生新解；

步驟4 若新解支配當(dāng)前解，則用新解替換當(dāng)前解，并回到步驟2；否則進(jìn)入步驟3 直到所有關(guān)鍵工序都被選中過；

步驟5 結(jié)束當(dāng)前解的鄰域搜索并選定種群中的下一個解，循環(huán)步驟2～5 直到種群中的Pareto 前沿的解都被搜索。

步驟6 根據(jù)變鄰域搜索結(jié)果按照式（16）更新pv的值。

3.5 加工時間重置算子

當(dāng)一個解中所有工序的起始加工時間確定后，只需調(diào)整每個工序的加工時間，就可以在不增大完工時間的前提下得到當(dāng)前解最小的總額外能耗［1］。為此，本文提出一種工序加工時間重置算子，用于確定解在當(dāng)前工序的起始加工時間下最優(yōu)的總額外能耗。

由于本文采取主動解碼方式對編碼方案進(jìn)行解碼，因此各工序的起始加工時間不可后移。加工時間重置算子在保證所有工序的起始加工時間不變的前提下將所有工序的加工時間盡可能增大，具體擴(kuò)展規(guī)則為：

步驟1 選擇選定解的一個未被選中過的工序O；

步驟2 判斷：1）編碼中工序O同工件上緊鄰的后一道工序的開始時間是否等于工序O的結(jié)束時間；2）編碼中工序O同機(jī)器上緊鄰的后一道工序的開始時間是否等于工序O的結(jié)束時間。若上述兩個判斷都為否，則工序O存在可擴(kuò)展的加工時間間隙，轉(zhuǎn)至步驟3；否則轉(zhuǎn)至步驟4。

步驟3 在工序O的加工時間范圍內(nèi)增大其加工時間，直到步驟2 中的判斷至少有一個滿足時，將此時的加工時間設(shè)為工序O的加工時間。

步驟4 若Pareto 前沿中的解都被選中，則終止；否則，轉(zhuǎn)步驟1。

圖2 所示個體通過加工時間重置算子優(yōu)化后得到的結(jié)果如圖5 所示。將工序O11、O12、O13、O32的加工時間向后擴(kuò)展（灰色塊所示為擴(kuò)展時間），在總加工時間不變的前提下，加工總額外能耗從39 下降到17，有效地提高了解的質(zhì)量。

圖5 加工時間重置算子示意圖Fig.5 Schematic diagram of processing time resetting operator

4 實(shí)驗(yàn)與結(jié)果分析

為了驗(yàn)證比較算法性能，本文仿真實(shí)驗(yàn)所用到的測試算例基于常用的Taillard 系列OSP Benchmarks［24］測試算例集進(jìn)行拓展。其中，工序的加工時間上界tr為OSP Benchmarks 中的工序加工時間，加工時間下界tl由式（17）產(chǎn)生。

從而產(chǎn)生了與OSP Benchmarks 對應(yīng)的6 種規(guī)模新測試算例，分別為（n×m）：4×4、5×5、7×7、10×10、15×15 和20×20，其中每個規(guī)模各包含4 個算例，共計(jì)24 個算例。為降低隨機(jī)性的影響，將算法在每個算例運(yùn)行20 次。所有算法使用Matlab R2018b編程，在 CPU 為 Intel Core i7-8750H（2.20 GHz），內(nèi)存為4 GB 的環(huán)境下運(yùn)行。

c1、c2和T為MOHEA 的重要參數(shù)，本文通過正交實(shí)驗(yàn)確定其具體取值。在參數(shù)的取值范圍內(nèi)均勻選擇了5 個水平，如表3 所示。其中，參數(shù)T用于表征進(jìn)化初期的迭代次數(shù)，因此在［0.1×Maxiter，0.5×Maxiter］內(nèi)取值較為合適。選擇了正交實(shí)驗(yàn)表L25 進(jìn)行實(shí)驗(yàn)，實(shí)驗(yàn)在算例5×5、10×10、15×15 上分別運(yùn)行20 次，并用獲得非支配解級的IGD 值作為評價指標(biāo)，如表4 所示。對各組IGD 結(jié)果進(jìn)行處理，得到各參數(shù)對結(jié)果的影響，如表5 所示?？梢钥闯鰞?yōu)勢水平為c1的第4 水平、c2的第4 水平和T的第3 水平。因此，最終確定的參數(shù)值為c1=0.7、c2=0.7，T=0.3×MaxIter。

表3 參數(shù)各水平取值Tab.3 Values of parameters with different levels

表4 正交實(shí)驗(yàn)表及對應(yīng)平均IGD值Tab.4 Orthogonal experiment table and corresponding average IGD values

表5 參數(shù)水平對應(yīng)的平均IGD值Tab.5 Average IGD values corresponding to parameter levels

其他參數(shù)設(shè)置為：種群數(shù)n=100，最大遷入率I=1，最大遷出率E=1，最大變異率mmax=0.7。

4.1 新測試算例的收斂時間

算法求解時間應(yīng)設(shè)置在收斂時間附近以得到最優(yōu)的計(jì)算結(jié)果。由于新算例的收斂時間未知，首先進(jìn)行預(yù)實(shí)驗(yàn)以確定不同測試算例使用MOHEA 求解的收斂時間。預(yù)實(shí)驗(yàn)基本過程為：1）設(shè)置足夠長的算法運(yùn)行終止時間；2）每隔10 s記錄一次算法得到的Pareto 前沿；3）計(jì)算終止后，計(jì)算每個時間間隔采集到的Pareto 前沿解的GD 和IGD 值（算法最終得到的Pareto 前沿看作是算例的真實(shí)Pareto 前沿），當(dāng)GD 和IGD 值不再隨計(jì)算時間的增加而明顯變化則說明計(jì)算收斂。

如圖6 所示為算例5×5_2 在不同運(yùn)行時間下Pareto 前沿的GD 和IGD 值?？梢钥闯觯?dāng)時間到80 s 時，兩個指標(biāo)值基本接近0 并且不再明顯變化，表明80 s 時算法得到的結(jié)果已經(jīng)收斂。同理，得到了其他算例的收斂時間如表6 所示。

圖6 算例5×5_2的收斂情況Fig.6 Convergence of example 5×5_2

表6 其他算例的收斂時間Tab.6 Convergence times of other examples

4.2 改進(jìn)策略的有效性驗(yàn)證

為了驗(yàn)證本文所提出的三種改進(jìn)策略對MOHEA 性能提升的有效性，設(shè)計(jì)了三個MOHEA 的變體進(jìn)行對比：1）MOHEA-1，使用傳統(tǒng)的遷移變異策略，同時包含自調(diào)整變鄰域搜索算子與加工時間重置算子；2）MOHEA-2，不包含自調(diào)整變鄰域搜索算子，僅包含改進(jìn)遷移與變異策略及加工時間重置算子；3）MOHEA-3，不包含加工時間重置算子，僅包含改進(jìn)遷移與變異策略及自調(diào)整變鄰域搜索算子。

將MOHEA 及MOHEA-1、MOHEA-2 和MOHEA-3 在24 個算例上獨(dú)立運(yùn)行20 次，GD、IGD 和Δ的統(tǒng)計(jì)結(jié)果分別如圖7～9 所示。圖中：A0、A1、A-Ⅱ和A3 分別表示MOHEA、MOHEA-1、MOHEA-2 和MOHEA-3；菱形表示20 次運(yùn)行結(jié)果得到的評價指標(biāo)平均值；子圖名為算例工件數(shù)-算例序號。

從圖7～9 可以看出：

圖7 GD結(jié)果的箱線圖Fig.7 Boxplots of GD results

1）MOHEA 在除4-4、10-3 和15-2 外的所有算例上都獲得了最優(yōu)的GD 均值，因此三種改進(jìn)策略有效地提升了算法的收斂性；

2）MOHEA 在所有算例上都獲得了最優(yōu)的IGD 均值，因此三種改進(jìn)策略對算法綜合性能有顯著提升；

3）可以看出四種算法在Δ方面的差異不大，MOHEA 在除5-2 和20-2 外的所有算例上都獲得了最優(yōu)的Δ均值，因此三種改進(jìn)策略在對算法求解質(zhì)量和綜合性能提升的同時沒有明顯影響分布性與多樣性；

4）在IGD 方面（圖8 所示），MOHEA-1 在小規(guī)模算例上（4×4、5×5、7×7、10×10）具有明顯的劣勢，說明本文提出的改進(jìn)遷移變異算子對算法綜合性能的提升相較于其他兩個策略更為顯著，歸功于改進(jìn)遷移策略對算法全局搜索能力的提升以及改進(jìn)變異策略對避免算法過早收斂起到的重要作用；

圖8 IGD結(jié)果的箱線圖Fig.8 Boxplots of IGD results

5）MOHEA-3 在大規(guī)模算例上（15×15 和20×20）具有非常明顯的劣勢，表明當(dāng)問題規(guī)模增大或工序加工時間范圍增大時，加工時間重置算子能在不增大加工時間的前提下顯著降低總額外能耗，從而高效地改進(jìn)解的質(zhì)量。

圖9 Δ結(jié)果的箱線圖Fig.9 Boxplots of Δ results

4.3 MOHEA性能驗(yàn)證

為了驗(yàn)證本文提出的MOHEA 的性能，將其與目前廣泛應(yīng)用的多目標(biāo)求解優(yōu)化算法NSGA-Ⅱ（Non-dominated Sorting Genetic Algorithm Ⅱ）［25］、NSGA-Ⅲ（Non-dominated Sorting Genetic Algorithm Ⅲ）［26］和 SPEA2（Strength Pareto Evolutionary Algorithm 2）［27］進(jìn)行比較。對比算法的各參數(shù)在文獻(xiàn)常用的范圍內(nèi)選取多個取值，然后通過在算例上的仿真進(jìn)行選優(yōu)，最終確定的參數(shù)值如表7 所示。

表7 不同算法的參數(shù)設(shè)置Tab.7 Parameter setting of different algorithms

從表8 所示的評價指標(biāo)統(tǒng)計(jì)結(jié)果可以看出：在所有算例上，MOHEA 在GD 和IGD 兩個方面均獲得了全部的最優(yōu)均值，因此MOHEA 在求解MOOSPCPT 時要優(yōu)于另外3 種多目標(biāo)優(yōu)化算法；從標(biāo)準(zhǔn)差可以看出，MOHEA 的穩(wěn)定性也優(yōu)于另外3 種算法；同時，隨著問題規(guī)模的增大，MOHEA 的優(yōu)勢更加顯著，說明MOHEA 特別適合求解大規(guī)模的MOOSPCPT。各算法得到的Pareto 前沿解如圖10 所示。

表8 不同算法的計(jì)算結(jié)果Tab.8 Computational results of different algorithms

從圖10 中可以直觀地看出：隨著問題規(guī)模增大，MOHEA 的優(yōu)勢更加明顯，且在總額外資源消耗方面，MOHEA 能夠求得能耗最小的解。綜上所述，相較于NSGA-Ⅱ、NSGA-Ⅲ和SPEA2，本文提出的MOHEA 可以更有效地解決MOOSPCPT。

圖10 不同算法得到的Pareto前沿Fig.10 Pareto frontiers obtained by different algorithms

4.4 算法時間復(fù)雜度分析

對照MOHEA 的基本流程，下面給出MOHEA 各操作策略的時間復(fù)雜度：

1）快速非支配排序O（2N2）［25］；

2）遷入遷出率計(jì)算O（2N）［28］；

3）改進(jìn)遷移變異機(jī)制O（N）［28］；

4）自適應(yīng)變鄰域搜索O（N4）［29］；

5）精英解時間重置策略O(shè)（N）。

因此，MOHEA 的時間復(fù)雜度為O（N4）。同時，理論上耗時最長的策略為變鄰域搜索，與仿真實(shí)驗(yàn)的結(jié)果一致。結(jié)合圖7～8 可以得到：在大規(guī)模問題中，由于變鄰域搜索策略的貢獻(xiàn)相對較少，可以通過適當(dāng)減小變鄰域搜索的強(qiáng)度加速算法的收斂。

5 結(jié)語

本文針對現(xiàn)實(shí)生產(chǎn)場景中機(jī)床的加工時間可控的現(xiàn)象，研究了加工時間可控的開放車間調(diào)度問題。以最小化完工時間與總額外能耗為目標(biāo)建立了調(diào)度模型，并設(shè)計(jì)了MOHEA 進(jìn)行求解。通過仿真實(shí)驗(yàn)，首先確定了測試算例的合理收斂時間；驗(yàn)證了本文所提出的三種搜索策略對MOHEA 的搜索性能具有明顯的提升作用，尤其對算法的綜合性能提升最為顯著，同時，改進(jìn)遷移變異算子對算法綜合性能的提升在小規(guī)模算例上優(yōu)于其他算子，加工時間重置算子能高效地降低總額外能耗；最后，將MOHEA 與NSGA-Ⅱ、NSGA-Ⅲ和SPEA2 進(jìn)行對比，表明本文提出的MOHEA 可以更有效地解決加工時間可控的開放車間調(diào)度問題，并且隨著問題規(guī)模的增大，MOHEA 的優(yōu)勢更加顯著。