基于量子計算和威布爾分布的混合CHIO算法求解JSP 問題*

亓祥波 趙品威 王 潤

（沈陽大學機械工程學院，遼寧 沈陽 110044）

在制造業中，作業車間調度問題（job-shop scheduling problem，JSP）被證明是一種極具挑戰性的約束組合優化問題，也是一個典型的NP-hard 問題。這類問題的特點是，無法找到一種有效的算法在多項式時間內準確地求解其最優解[1]。它涉及如何合理地安排工作流程和資源分配，以最大化生產效率和利潤。隨著全球制造業的發展和競爭的加劇，對調度的優化變得越來越重要，設計科學合理的調度方案可以提高生產效率和經濟效益，因此，對調度優化算法的研究成為調度領域研究熱點之一。

近年來，受自然界各種生物群體行為、物理現象和自然現象等啟發設計出來各種智能優化算法，因其具有較強的全局搜索能力、自適應性、魯棒性、并行計算能力以及易實現等優點而被廣泛應用于離散調度的問題求解中。劉麗娜等提出了一種改進麻雀算法（ISSA）[2]，通過引入高斯變異擾動、正余弦搜索等策略提高了麻雀算法對離散調度問題的求解精度。為了解決狀態轉移算法在過分依賴初值、局部搜索能力差、收斂速度慢等方面存在的不足，吳貝貝等提出了一種量子狀態轉移算法（QSTA）[3]，該算法引入了量子旋轉門的思想，通過將最小化最大完工時間作為目標函數來解決離散調度問題，最后通過試驗驗證了該算法能夠有效地求解離散調度問題。陳斌等提出了一種多策略引導的電磁場優化算法（MGEFO）[4]，該算法通過引入了自適應電磁移動系數，并采用不同的引導策略來動態調整粒子間的正電距離和負電距離，通過實驗結果驗證了該算法在解決離散調度問題方面的可行性。亓祥波等針對離散調度問題的特性提出了一種基于交叉選擇的變鄰域蜂群算法（VABC）[5]，增強了人工蜂群算法在離散調度問題上的搜索能力。

以上針對調度問題設計的改進算法在一定程度上取得了成果，為進一步研究提供了參考。然而，這些改進算法大多是對較為傳統的優化算法的改良和設計，很少嘗試使用新型智能優化算法進行對調度問題進行求解。冠狀病毒群免疫優化算法[6]是2020 年Al-Betar M A 等提出的一種新的元啟發式優化算法，靈感來源于造成世界性傳播的新冠病毒，該算法模仿了群體免疫策略和社會距離概念，其基于3 種類型的個體病例用于群體免疫：易感、感染和免疫。每個易感人群與感染個體接觸而感染病毒，再通過個體是否免疫而更新引導群體免疫以達到算法收斂。

針對CHIO 后期搜索能力差、收斂精度低等缺點，已有不少學者對其進行了改進并應用在不同實際問題上，例如路徑規劃[7]、醫療檢測[8]和設計優化[9]等，但對于離散調度問題應用相對較少。本文以CHIO 為基礎并對其進行了改進將其應用到車間調度問題上，與其他幾種改進的優化算法進行了仿真對比實驗，證明了算法在求解JSP 方面有較強的能力。最后，將提出算法應用到發動機生產調度實例上，進一步驗證了所提出算法的可行性和優越性。

1 作業車間調度模型

1.1 JSP 問題的描述

JSP 問題可以理解為在具有特定工序的n個工件在m臺機器上進行加工。每個工序的加工機器以及加工每道工序所需的時間都是已知的。調度的目標是安排每臺機器上工件的加工順序，根據排序確定每個工序的開始和結束時間，以實現最小的完工時間。JSP 問題滿足以下假設：

①機器是獨立的。

② 加工準備時間已被計入加工時間。

③所有機器從0 時刻開始可用。

④ 各工件加工的優先級相同，同一工件不同工序之間存在先后約束。

⑤ 每道工序在同一時刻只能由一臺機器進行加工。

⑥ 加工過程中沒有缺料或機器故障的情況。

1.2 作業車間調度的數學模型

本文使用整數規劃模型[10]對JSP 進行描述，其模型如下：

其中：式（1）為最小化最大完工時間的目標函數；式（2）為工藝約束；式（3）為機器約束；cik為工件i在機器k上的完成時間；pik為工件i在機器k上的加工時間；M是一個充分大的正數；xihk是指示系數，xijk是指示變量，其含義為

1.3 基于工序的編碼

選擇恰當的編碼規則對于獲得優秀的調度方案至關重要，尤其是在表示JSP 的解時。本文采用基于工序的編碼方式來實現這一目的。工序編碼由工件號組成，編碼的總長度等于任務的總加工工序數。在解碼過程中，按照從左到右的順序對工序編碼進行安排加工。當一個工件號第一次出現時，表示該工件正在進行第一道工序的加工；當同一個工件號第二次出現時，表示該工件正在進行第二道工序的加工；以此類推，直到所有工序都被安排完成。以一個3×2 規模問題為例，具體見表1。

表1 編碼對應規則

對離散變量除以機器總數向上取整，即得到基于工序的編碼[1,3,3,2,1,2]。

2 量子混合CHIO

2.1 標準CHIO 算法

在CHIO 算法中，有3 種類型的個體病例：易感、感染和免疫。傳播冠狀病毒的速度取決于感染者如何與其他社會成員直接接觸。群體免疫是當大多數群體具有免疫力時，群體達到的一種狀態，這種狀態可以防止疾病的傳播。這些概念是根據優化概念建模的，既模仿群體免疫策略，又模仿社會距離概念。CHIO 對群體免疫策略進行了建模：在初始化種群HIP 后，隨機產生與HIS一樣多的一組病例（即個體）。生成的病例以二維矩陣n×HIS形式存儲，且每個個體均對應一個狀態向量Si和年齡向量Aj，HIP 為

式中：HIS為種群個數；n為問題維度。

在算法迭代過程中，xi病例的j基因要么保持不變，要么受到社會距離的影響，其影響根據基本繁殖率BRr，取值為0.05，見式（8）。

式中：r表示(0,1)的隨機數；(t+1) 表示(t)更新后的基因值。

式中：r表示(0,1)的隨機數；(t)表示從易感個體中選擇的個體，此時m={i|Si=0}。

式中：r表示(0,1)的隨機數；(t)表示從免疫個體中選擇的最好的個體，此時v={i|Si=2}。

隨后更新群體免疫群體，計算新個體xi(t+1)的免疫率（即適應度值），如果f(xi(t+1))>f(xi(t))，則xi(t+1) 取代xi(t)。且如果該個體為感染個體，則年齡向量Aj增加1。其中個體對應狀態向量Si也隨著適應度值改變，見式（12）。

式中：?f(x) 表示當前種群的平均免疫率；is_Corona表示當前個體是否繼承過感染個體基因。如果繼承過感染個體的基因、當前個體為易感人群且免疫率小于平均免疫率，則該個體感染病毒，此時該個體的 Si=1；如果當前個體免疫率大于平均免疫率且當前個體為感染個體，則該個體為免疫個體，此時該個體的 Si=2。

最后，當個體年齡大于最大年齡后，該個體將重置：

式中：ub、lb分別表示問題的上下界；rand表示(0,1)的隨機數。

2.2 QCHIO 算法

2.2.1 量子計算

由于CHIO 在求解JSP 時存在收斂精度低、尋優能力差等缺點，因此將量子計算[11]的思想引入到CHIO 算法。其中單量子比特的疊加狀態為

在量子計算中，實現邏輯變換的基礎是量子邏輯門。這些門操作主要用于改變量子態的性質。常見的幾種量子門包括Hadamard 門、相位門和量子旋轉門等。量子旋轉門表達式為

量子旋轉門的更新表達式如下：

在更新過程中對工序序列進行量子優化，則工序序列經過單鏈編碼量子旋轉門變換后表達式為

式中：Qxi表示xi通過量子計算產生的新個體。

2.2.2 威布爾分布變異

威布爾分布又稱為韋伯分布、韋氏分布，是一種連續概率分布，一般用于工業制造、預測天氣、可靠性和失效分析、量化壽險模型的重復索賠等。其概率密度函數（probability density function，PDF）為

式中：A為尺度參數，表示函數的走勢；B為形狀參數，表示函數的形狀。圖1 表示尺度為A、形狀為B的威布爾分布的PDF。Sahoo S K 等[12]將威布爾分布算子應用到飛蛾撲火算法上證實了其可以當前位置搜索到有用的信息，因此將改進威布爾分布算子引入CHIO 以增加算法的多樣性，通過改進威布爾分布算子的大步長和小步長來平衡算法的局部搜索和全局開發能力，以增強算法的收斂精度。具體實現見式（21）～式（23）。

圖1 威布爾分布示意圖

式中：S TEP表 示大步長，step表 示小步長，sign(rand-0.5) 表示隨機取得的-1 或1 的向量，wbl(1,1)表示A、B均為1 時生成的威布爾分布的隨機數，norm為范數函數，j表示個體的第j維，xi表示第i個個體，xbest表示當前種群中最優個體，xnew表示xi生成的新個體。

2.2.3 β-登山局部搜索

登山（hill climbing，HC）算法是一種簡單的貪心鄰域搜索算法。該算法每次從當前解的臨近解空間中選擇一個最優解作為當前解，直到達到一個局部最優解。但其主要缺點是易陷入局部最優解。為克服該缺點，Al-Betar M A 等提出了一種 β-登山算法[13]。在CHIO 中引入了 β-登山算子以平衡局部搜索和全局開發，從而提高CHIO 的收斂精度。具體實現見式（24）和式（25）。

式中：xj表 示當前個體的第j維，j是從當前個體 中隨機選擇的一個維度；r1表示(0,1)的隨機數；表示一個常數，取值為0.001。

式 中：xi表示當前個體的第i維，lbj、ubj分別表示問題求解范圍的下界和上界在第j維的位置，r2和R均表示(0,1)的隨機數。最后采用貪婪選擇，如新個體的適應度值大于原個體，則新個體取代原個體。

2.2.4 多鄰域搜索

在優化算法中，鄰域搜索是一種通過使用鄰域函數來創建解的鄰域，并在鄰域中找到比當前解更優的解來進行替換的方法。在本文中，經過迭代過程后，對種群中的最優個體依次進行了交換、倒置和倒序操作。每次操作都會保留適應度值較優的個體，以增強算法的挖掘能力。

具體的多鄰域搜索方法如下。

（1）交換操作：假設加工任務總加工工序數為N，在1～N中隨機選擇兩個不同的位置，然后交換這兩個位置的加工工序。

（2）倒置操作：假設加工任務總加工工序數為N，在1～N中隨機選擇一個位置Oa，然后將該位置與其后面的工序位置Oa+1進行交換。如果選擇的位置是最后一個位置，則將該位置的工序與工序列表中的第1 個位置的工序進行交換。

（3）倒序操作：假設加工任務總加工工序數為N，在1～N中隨機選擇兩個不同的整數為a和a+p，假設未進行倒序操作前的工序編碼為O＝[O1,O2,O3,···,Oa,Oa＋p,···,ON]，則執行倒序操作后的序編碼為Onew＝[O1,O2,···,Oa＋p,Oa+p-1,···,Oa,···,ON]。

2.3 QCHIO 算法的流程

QCHIO 的流程圖如圖2 所示。

圖2 QCHIO 流程圖

3 實驗驗證與結果分析

為了驗證所提出算法QCHIO 的性能，采用了車間調度上典型的FT 算例與LA 算例進行仿真測試。實驗PC 設備參數為：CPU 為AMD Ryzen7 5700U、主頻為1.8 GHz、運行內存為16 GB、操作系統為Windows 11。編程環境為Matlab 2018b。

3.1 對比算法及參數設置

為了便于對比分析，選取了幾種較為先進的算法作為對比算法，其中有改進麻雀算法（ISSA）[2]、量子狀態轉移算法（QSTA）[3]、多策略引導的電磁優化算法（MGEFO）[4]、標準CHIO[6]算法。為了驗證改進策略的有效性，對算法進行了消融實驗，見表2，其中ICHIO1、ICHIO2、ICHIO3、ICIO4 分別表示在QCHIO 的基礎上無量子計算算子、無威布爾分布算子、無 β-登山算子、無多鄰域搜索算子。為了便于對比，將實驗參數均設置種群個數為50，最大迭代次數為500，算法獨立運行30 次，記錄其平均值（avg）、最優值（best）作為實驗結果，結果見表3，QCHIO 相對于對比算法的改進提升水平見表4，同時選取了幾個具有代表性的算例給出了平均迭代曲線和甘特圖，如圖3 和圖4 所示。

圖3 在部分算例上的平均迭代曲線

圖4 在部分算例上的甘特圖

表2 消融實驗結果

表3 QCHIO 與部分改進算法在算例上的對比（%）

表4 QCHIO 相對于CHIO、ISSA、QSTA、MGEFO 的提升水平（%）

3.2 仿真結果

從表2 可以看出，QCHIO 相對于CHIO 在所有算例上無論是最優值還是平均值都得到了較大程度的提升，相較于ICHIO1、ICHIO2、ICHIO3、ICHIO4 也得到了不同程度的提升，從而證明了改進策略的有效性。其中相較于ICHIO1，QCHIO 在最優值上的提升大于平均值的提升，由此可以推斷加入量子計算算子有助于提升算法的全局開發能力。相較于CHIO2，QCHIO 在大規模問題的LA21、LA26、LA31、LA36 上均得到了較大的提升，由此推斷威布爾分布算子在求解大規模問題上具備一定優勢。文獻[14]指出全局最優解是在所有鄰域結構中的局部最優解中找到的，并且某種鄰域結構的局部最優解通常接近于另一種鄰域結構的局部最優解。相對于ICHIO3、ICHIO4，QCHIO 分別在當前最優解和全局最優解附近進行了鄰域搜索，基于實驗結果也同樣可得出，算法在最優解附近確實可以搜索到有用的信息，從而提升算法的局部搜索能力。

在表3 中，CHIO(50,500) 意為種群大小為50，總迭代次數為500，例如ISSA(100,500)表示在原文中種群大小為100，總迭代次數為500。在表4 中，?best表示QCHIO 相較于其他算法最優值的提升率，?avg表示相較于其他算法平均值的提升率，“－”表示無提升，以QCHIO 相對于CHIO 的提升水平為 例，?best=(CHIObest-QCHIObest)/ CHIObest,?avg=(CHIOavg-QCHIOavg)/ CHIOavg。

由表3 和表4 可知，QCHIO 相較于CHIO，最優值和平均值均得到了較大提升，尤其是在大規模問題上均取得了10% 左右的提升。相較于ISSA，雖然其設置了更大的種群規模，但QCHIO 在平均值上有7 個算例上均具有不同程度的提升，說明了QCHIO 在一定程度較ISSA 更穩定。相較于QSTA，QCHIO 與其實驗參數設置均保持了一致，從平均值角度看，QCHIO 在9 個算例上得到了提升，在最優值上有4 個算例得到提升，3 個算例保持一致，由此得出引入了量子計算和威布爾分布算子的QCHIO 在解決調度問題上相較于QSTA 具有更強的全局尋優能力和穩定性。與MGEFO 相比，該算法總迭代次數為QCHIO 的2 倍，但在10 個算例的最優值上QCHIO 有5 個算例得到提升，3 個算例保持一致，平均值上有9 個算例得到提升，由此可以推斷出QCHIO 較于MGEFO 在求解JSP 問題上更具優越性。

3.3 統計學分析

為了進一步驗證本文算法與其他算法在仿真結果在統計學上是否存在顯著性差異，通過使用95%置信度的配對t檢驗對各算法在10 個算例上的平均值進行試驗并記錄其實驗結果，見表5。在四組配對數據中顯著性水平p均小于0.05，另外，Cohen's d 值是用來衡量效應量大小的，即算法之間差異的幅度大小。當Cohen's d 值越大時，說明算法之間的差異越大。根據常用的分級標準，效應量可分為小、中、大三個等級，對應的臨界點分別是0.2、0.5 和0.8。根據表中的數據可知，4 組數據的Cohen's d 值均大于0.8，表示算法之間存在較大的差異性。因此可以得出結論，QCHIO 與所有對比算法在統計學上均存在顯著性差異。

表5 各算法配對t 檢驗結果

4 實例驗證

為了進一步驗證QCHIO 求解車間調度問題JSP 的可行性和有效性，本文在某公司汽車發動機生產的作業排產問題上進行了例證，根據文獻[15]整理了某型號的發動機的生產工件工藝和數據，各工件的加工順序見表6，加工時長見表7。

表6 各工件加工工藝

表7 各工件加工時長 min

根據表6 可以看出需要排產的發動機的工件生產過程需要10 個工件，每個工件需要6 道工序，為10×6 規模的問題，解集規模為(10!)6。目前該公司所采用的人工經驗調度方式的原則為：優先加工總時長最短的工件，同時在整個調度中盡可能將總加工時長相對較短的工件先加工，同時盡可能保證在調度中排在前面的加工機器能夠連續工作，不出現空閑狀況，即優先保證粗加工車M1 優先加工，所得到的該型號發動機的主要10 個工件總工時為199 min。

運用本文所提算法QCHIO 對其進行求解，設置種群規模為100，最大迭代次數為100，算法獨立運行10 次，所得結果見表8。最優加工方案甘特圖如圖5 所示。

圖5 最優加工方案甘特圖

表8 實例求解結果

所得最優加工方案為[3,3,9,10,3,5,9,3,5,4,10,9,7,9,5,6,4,9,10,1,5,1,4,9,7,2,3,10,8,4,6,1,1,5,10,4,7,7,3,8,1 0,5,6,2,8,2,1,4,6,1,8,7,2,6,8,7,2,6,8,2]，加工總時長為122 min。較原始人工調度方法提升了38.7%，較原文提出算法DFACO 提升了10.3%。盡管只是對工件加工工時進行初步改善，但這是在該公司所生產的發動機的其中一種型號，對10 個工件進行的一次10×6 規模的調度優化?？紤]到整個公司的產品種類和數量較多，每年實際的生產成品量巨大，并且包括單獨外售的零部件，所以從整體上看，這種改善效果和產生的效益是較為可觀的。

5 結語

隨著全球制造業的發展和競爭的加劇，車間調度作為提高生產車間效率的關鍵，不斷對其進行優化變得越來越重要。本文針對最小化最大完工時間為求解JSP 目標，提出QCHIO 算法，并通過在10個算例進行仿真實驗和消融實驗，與其他幾種算法上進行了對比分析，并對結果進行了顯著性分析，最后將其應到一個實際生產案例上，得出了以下結論：

通過對10 個算例的消融結果進行分析，驗證了所提算法QCHIO 在改進策略的有效性。通過對算例仿真結果的分析得出QCHIO 具有較強的全局搜索能力，收斂精度較高等特點。引入量子計算算子，通過量子相關性實現全局搜索和快速收斂的目標，增強了算法的全局開發能力，威布爾分布算子進一步增加了算法的多樣性和搜索能力，在搜索空間的尾部有更好的擬合性能，能夠更好地探索整個搜索空間。β-登山算子和多鄰域搜索算子是針對局部最優解問題的改進措施。這兩種策略的引入有效地避免了算法陷入局部最優解而無法跳出的問題，提高了算法的局部搜索能力。通過對實例結果分析進一步驗證了所提算法的可行性和優越性。

本文的研究對象主要為調度問題中JSP 問題，作為實際生產的簡化模型，除了本文應用求解的發動機生產，金屬沖壓件和彈簧[16]等生產都符合JSP的生產特點，由此看出本文的研究具有一定的代表性，所提算法也可為求解調度問題提供一種思路。但隨著綠色制造的提出，越來越多的企業重視企業生產能耗，以求由此來降低成本，因此在同時考慮最大完工時間和能耗的多目標問題上有待進一步研究。