蔣華偉 郭 陶 楊 震 趙麗科
(糧食信息處理與控制教育部重點實驗室(河南工業大學) 鄭州 450001)
(河南工業大學信息科學與工程學院 鄭州 450001)
自然災害及重大公共衛生事件的發生會嚴重威脅人們的生命財產安全,并對社會生產造成不利影響[1]。為最大限度降低災害帶來的損失,就需要研究和構建科學合理的物資(如糧食)應急調度模型,來獲取最優路徑,減少物資配送時間,以保證災區救援物資的充足供應,從而為開展及時高效的救援工作提供物資保障。
災后物資應急調度的本質是多種約束條件下的車輛路徑問題(Vehicle Routing Problem, VRP)[2],即具有多配送中心帶時間窗的VRP。作為經典的組合優化問題,VRP已被證明為NP-hard問題[3],其求解方法主要有精確算法[4,5]和啟發式算法[6,7]兩類。精確算法針對具體問題建立相應的數學模型,利用數學方法求出問題的最優解,但其計算時間隨問題規模的增大呈爆炸式增長,因此只能解決規模相對較小的問題。啟發式算法是相對于精確算法提出的,在可接受范圍內給出問題的解,相比于精確算法,在處理大規模VRP時,具有更高的魯棒性、可行性。因此國內外學者主要采用啟發式算法對VRP及其變體問題進行研究,如在求解帶時間窗VRP時,Marinakis等人[8]采用3種不同的自適應策略優化粒子群算法,分別用于初始解的生成、解的移動以及算法參數的自適應調整,以提高算法的求解性能;此外,Ramachandranpillai等人[9]將改進螢火蟲算法與脈沖神經系統結合,使其可以快速地搜索解空間,從而提高算法求解的收斂速度。為求解具有多配送中心VRP,Lahyani等人[10]提出了基于混合自適應大鄰域搜索算法,算法結合3種插入、5種移除啟發式算法和4個后優化局部搜索過程,來靈活解決多配送中心VRP;另外,胡蓉等人[11]提出結合聚類分解策略的增強蟻群算法,將多配送中心問題分解為多個單配送中心問題,以控制問題的求解規模。求解多目標VRP時,Zhang等人[12]設計了兩個多目標局部搜索算法,結合兩個算法并采用附加技術增強局部搜索算法,提高算法的求解性能;此外,Long等人[13]結合遺傳算法與局部搜索策略,提出一種基于帕累托的進化算法,用于求解多目標規劃問題。在求解其他變體問題時,駱劍平等人[14]將改進冪律極值動力學優化引入混合蛙跳算法,以提高算法求解容量約束VRP的性能;李國明等人[15]將禁忌搜索與最近鄰算法相結合,并對禁忌長度等進行自適應調整,引入自適應懲罰系數,以提高禁忌搜索算法在求解隨機VRP時的尋優能力和魯棒性;Xiang等人[16]提出基于成對鄰近學習的蟻群算法,用于求解動態VRP,采用成對鄰近學習方法研究變化前的最優路徑,并預測變化后最優路徑中客戶的局部訪問順序。上述算法在求解VRP及其變體問題時,針對不同問題的特點設計相應的算法進行求解。在解決算法易陷入局部最優這一問題時,主要借助變異操作跳出局部極值,具有很強的隨機性,無法保證變異后個體的優劣,并且未充分考慮種群多樣性與算法陷入局部極值間的密切關系,僅使用適應度函數選擇個體,因而無法在算法迭代周期內維持高種群多樣性以最大限度發揮啟發式算法的優勢。
針對啟發式算法存在的缺點,以及實際物資應急調度問題龐大的解空間和整體編解碼的復雜性,本文提出一種改進離散鯨魚群算法(Improved Discrete Whale Swarm Algorithm, IDWSA),對具有多配送中心帶時間窗的物資應急調度問題(Material Emergency Scheduling Problem with Mul- tiple Distribution Centers and Time Windows, MESPMDCTW)進行求解。首先分析問題中各種復雜的約束條件,為其建立嚴謹的數學模型;其次,為獲得高質量的初始種群,提出混合初始化策略生成初始解,即動態模糊聚類(Dynamic Fuzzy Clustering, DFC)與隨機生成相結合;然后設計分別以相似配送順序和相同配送中心為比較項的個體移動規則,采用自適應柯西變異算子對個體進行變異,并提出路徑選擇策略;此外,構造全局評價函數,用于評價子代個體對種群的貢獻度,以避免僅使用適應度函數對個體進行評價時的局限性,使得算法在迭代后期仍可保持高種群多樣性;最后在Solomon標準測試集[17]上進行仿真實驗,并與其他算法進行對比分析,驗證IDWSA的有效性。
實際物資應急調度中存在多種復雜并相互制約的因素,為了簡化模型本文給出以下約束條件:(1)每個受災點僅能由一個配送中心的一輛車完成配送任務;(2)每輛車所訪問的各受災點的總需求量不可超過該車輛的最大載重;(3)每輛車的行駛距離不可超出車輛最大行駛里程;(4)所有車輛在完成配送任務后需返回其始發配送中心;(5)車輛對其配送受災點的訪問時間不可超過該受災點的最晚送達時間。
MESPMDCTW問題可由有向圖描述,在有向圖中頂點表示配送中心及受災點,有向邊表示車輛的可行駛路徑,如圖1所示。

圖1 MESPMDCTW問題的有向圖表示
假設該問題的配送中心集合D= {d1,d2,···,dm};每個配送中心擁有vi輛車,其中i表示第i個配送中心,每輛車的最大載重為Q;受災點集合R={r1,r2,··,rn},每個受災點的物資需求量為qj(j=1,2,···,n)。
根據所述,本文構建MESPMDCTW問題的數學模型如下:
(1)目標函數
MESPMDCTW問題的目標是使總運輸成本最低和使用的車輛數目最少,其中總運輸成本最低可以表示為總運輸路徑最短、總服務時間最小等,車輛數目最少可使用車輛費用最小化表示。該目標函數用于衡量所生成車輛路徑的優劣程度,目標函數值越大表明采用該路徑對各受災點進行物資配送時所耗費的代價越大。

由于各受災點對物資配送時間具有一定限制,當車輛提前或者延時到達時,要對該車輛進行懲罰,即增加其成本。其中ti表示車輛到達受災點i的實際時間,ETi,LTi分別表示受災點i的最早和最晚訪問時間;a,b分別表示車輛提前到達和延時到達的懲罰系數,由于不允許延時到達,所以b為一個很大的正數。

式(4)表示每個受災點只能由一個配送中心的一輛車訪問,式(5)表示第t個配送中心服務的受災點的物資總需求量,式(6)表示訪問受災點j的車輛k流入流出的平衡條件。
鯨魚群算法(Whale Swarm Algorithm, WSA)是2017年由Zeng等人[18]基于群體智能提出的一種新元啟發式算法,經大量實驗證明,與遺傳算法(Genetic Algorithm, GA)、差分進化算法(Differential Evolution Algorithm, DEA)等綜合對比,WSA的求解性能更優。其基本思想如下:首先在定義域內隨機生成一組解作為初始種群,種群中每條鯨魚代表解空間中的一個候選解;然后根據種群中鯨魚的適應度值和位置,依次為每條鯨魚搜索其“更優且最近”目標鯨魚,即周圍適應度更優的鯨魚中距離其最近的鯨魚;最后每條鯨魚均以其目標鯨魚為導向以某種方式進行移動,從而產生新一代種群,其向目標鯨魚移動的方式如式(7)。

由式(7)可知,當鯨魚X與其“更優且最近”鯨魚Y之間的距離很近時,鯨魚X會積極地向鯨魚Y隨機移動;反之,鯨魚X會消極地向鯨魚Y隨機移動。
WSA主要適用于求解連續性問題,而MESPMDCTW屬于離散問題,無法直接使用它對問題進行求解,因此本文對基本WSA進行優化,主要改進包括:(1)設計新的鯨魚個體編解碼方式,并提出相應的個體間距離計算方式;(2)為提高初始種群的質量及其多樣性提出混合初始化策略;(3)提出自適應柯西變異算子,用于更新種群中無引導個體的個體;(4)設計路徑選擇策略用于個體中各受災點的重新規劃;(5)提出以相似配送順序和相同配送中心為比較項的兩種個體移動規則;(6)構造全局評價函數衡量個體對種群的貢獻度用于選擇個體組成子代種群。
IDWSA算法的完整步驟如下:
(1)采用混合初始化策略對種群進行初始化得到大小為n的初始種群pop;
(2)采用式(14)計算種群中個體的適應度值,并從中選擇適應度值最大的個體作為當前種群的最優個體,記為best。
(3)尋找種群中每個個體的“更優且最近”個體Y,對于個體X,若Y存在,則分別根據以相似配送順序和相同配送中心為比較項的移動規則對X進行移動,獲得兩個子代個體;若Y不存在,則對X分別進行兩次自適應柯西變異,獲得兩個子代個體;對種群中所有個體操作后,獲得大小為2n的子代種群new_pop;
(4)分別采用式(14)和式(16)計算子代種群new_pop中每個個體的適應度值和貢獻度,從中選擇適應度值最大的個體記為new_best并與種群pop中的best作比較,保留適應度值大的個體記為best,然后根據個體的貢獻度,從new_pop中選擇貢獻度最大的前n–1個個體與個體best組成子代種群pop。
(5)判斷是否滿足終止條件,若不滿足則返回步驟(2),反之則進行步驟(6)。
(6)獲得最終種群pop,從中選擇適應度值最大的個體作為問題的最終解。
MESPMDCTW問題中每個受災點可以由任意配送中心的任意車輛配送,且其在車輛中的配送順序是隨機的,因此該問題具有龐大的解空間,為了減小搜索空間,提高算法的搜索效率,本文提出3層編碼方式,包括配送中心編碼(Distribution Coding, DC)、車輛編碼(Vehicle Coding, VC)和順序編碼(Sequential Coding, SC),如圖2所示。其中,DC值為對應受災點所屬配送中心,VC值為與DC對應的受災點的車輛號,SC為與VC段對應的該受災點在車輛中的配送順序。
解碼時,首先由DC確定受災點所屬配送中心,然后由VC確定同一配送中心下各受災點的配送車輛,最后根據SC確定車輛訪問各受災點的順序。圖2中該示例表示受災點1,2,4,6,7由配送中心1的兩輛車配送,第1輛車訪問順序為1,6,7,第2輛車訪問順序為4,2;受災點3,5,8,9由配送中心2的兩輛車配送,第1輛車訪問順序為3,8,第2輛車訪問順序為9,5。

圖2 3層編碼示例
由于基本WSA中的距離公式僅適用于求解連續性問題,而不能直接用于求解離散問題,因此根據MESPMDCTW問題的特點及本文設計的3層編碼方式,提出一種計算鯨魚個體間相對距離的方法,如式(8)。

高質量的初始種群不僅能加快算法的收斂速度,而且可產生質量更優的最終解。目前,對種群進行初始化的方法有全局最小化處理時間規則[19]、MinEnd啟發式規則[20]、全局搜索和局部搜索方法[21]等。但這些種群初始化方法均是基于車間調度問題提出的,并不適用于本文問題。因此為提高初始種群的質量并保持其多樣性,防止種群在迭代中過早喪失多樣性而陷入局部極值,根據MESPMDCTW問題的特點,本文提出混合初始化策略,即結合DFC和隨機生成方法,按照一定比例生成初始種群,如圖3所示。
圖3中,初始種群中35%的個體由DFC算法生成,65%的個體隨機生成。在隨機生成個體時,首先生成一個兩倍大的臨時種群,然后根據適應度函數計算臨時種群中個體的適應度值,選擇適應度值最大的前12.5%的個體組成隨機生成部分25%的個體,最后從剩余87.5%的個體中隨機選擇個體作為隨機生成部分40%的個體。

圖3 混合初始化策略
DFC主要根據各受災點位置、訪問時間窗及服務時間,進行相似性分析,使同一類別的受災點距離最近且各受災點訪問時間重疊最小,以此將所有受災點劃分為與配送中心數相同的幾類。此外,為防止各配送中心所分配受災點數極度不均勻導致整體配送效率降低,本文在對各受災點進行聚類后,對形成的受災點集合進行均衡化處理,使得每類受災點數大致相同,即當某一類受災點數較多時,根據式(9)選擇部分受災點將其移向數目相對較少的受災點集合,循環移動直至達到各受災點集合大小均衡。

WSA中個體的優化是以其“更優且最近”個體為引導的,但種群中存在某些不具有引導個體的個體,為了對它們進行更新,根據式(12)的柯西變異算子,本文改進提出了式(13)所示的自適應柯西變異算子。


其中,xij為個體i的第j維位置,n為種群大小,C是由t=1的柯西分布函數產生的隨機數,[xmin,xmax]是各受災點的編號區間。
WSA中個體的運動以其“更優且最近”個體為引導,隨著個體的不斷移動,種群中個體不斷收斂。在算法迭代前期,個體位置分散,種群平均位置較大,隨著個體不斷地向最優解方向移動,種群平均位置逐漸變小,算法逐漸收斂,因此種群平均位置變化與算法收斂特性是一致的,采用種群平均位置作為控制變異步長的參數有利于提高算法在迭代前期的搜索能力,并能在后期加快算法的收斂速度。
鯨魚個體向其“更優且最近”個體移動時,為了在其引導個體周圍進行細致的搜索,以最大概率尋找區域內最優個體,本文提出路徑選擇策略,使個體根據該策略向其引導個體移動,其基本思想如表1。首先根據引導個體即“更優且最近”個體確定配送中心及各受災點之間的路徑權值矩陣W,若節點相鄰則將兩節點之間的路徑權值置為1,反之則置0;根據當前已生成路徑的最后一個節點,從未被訪問的受災點集合S中選擇與最后一個節點相連權值最大的受災點,作為下一個要訪問的受災點,以此循環直至未被訪問的受災點集合S為空。

表1 路徑選擇策略
IDWSA中,鯨魚個體向其“更優且最近”個體移動,采用適應度函數式(14)衡量其質量。該適應度函數是基于MESPMD CTW的目標函數式(1)構建的,主要有3個影響因素,分別為個體X所表示的路徑總距離、X所使用的總車輛數以及X中車輛訪問所有受災點的總延遲時間,用于從每代種群中選擇質量最優即適應度值最大的個體。

由于MESPMDCTW屬于離散問題,基本WSA中的個體移動方式無法使用,因此本文提出以相似配送順序和相同配送中心為比較項的兩種個體移動規則,具體的移動規則如圖4所示。

圖4 個體移動規則
為了說明個體移動規則,圖5給出了基于相似配送順序的個體移動示例。圖5中,個體Y為個體X的引導個體,在個體X和Y中,受災點3,9都位于對應車輛的第1個訪問順序,受災點5,6都位于對應車輛的第2個訪問順序,它們具有相同的配送順序,因此將受災點3,5,6,9根據其在引導個體Y中的位置復制到新個體Z中;然后隨機生成兩個整數P1=2和P2=5,將引導個體Y中索引P1和P2之間(不包括P2)的受災點復制到Z中相應位置;最后將X中剩余受災點依次復制到Z中。

圖5 基于相似配送順序的個體移動示例
對于種群中的個體X,在對其進行1次移動后會產生2個子代個體,則完成1輪搜索后,種群規模將變成原來的兩倍。由于個體是向其“更優且最近”個體移動的,當采用適應度函數選擇個體時,隨著迭代次數增加,個體之間距離逐漸縮小,使個體逐漸聚集在某一區域,種群多樣性降低,從而使算法陷入局部極值,難以求出問題最優解。因此為了維持種群多樣性,使算法在進行迭代求解時可以在問題的較大解空間上進行搜索,本文基于個體的適應度值構造了如式(16)所示的全局評價函數,用于衡量個體對整個種群的貢獻程度,從而在算法迭代期間對個體進行選擇以組成新的子代種群。

全局評價函數涉及子代個體、父代個體及整個種群的適應度值,在計算個體x對種群的貢獻度時,綜合考慮其父代個體在父代種群中的影響程度、個體x在子代種群中的優劣程度以及父代個體向其引導個體移動生成個體x時的優化程度。根據貢獻度選擇個體時,由于不僅考慮當前個體x的適應度值,同時考慮其父代個體所產生的影響,使得適應度值大的個體其貢獻度不一定大,因此對于適應度值較小但可能位于最優解區域的個體x仍有機會被選擇,從而擴大算法的搜索空間。采用全局評價函數選擇個體,根據個體的貢獻度在每次迭代時對個體進行選擇,使生成的子代種群中個體間具有較大的差異性,即種群個體分布在問題解空間的較大區域,從而不會使算法過早陷入局部最優,且每次迭代中都會保留當前種群中適應度值最大的個體,因此可以在最大限度維持種群多樣性的情況下求得問題的最好解。
為驗證IDWSA求解MESPMDCTW的有效性,本文在Solomon標準數據集上進行仿真實驗,該數據集主要分為6類:C1,C2,R1,R2,RC1和RC2,共56個測試集。其中C類是聚類數據,R類數據是隨機分布的,RC類數據則是C類和R類的混合數據。由于本文問題涉及多配送中心這一約束條件,因此本文選取隨機分布的R類測試集中的R101進行仿真實驗,并在上述數據集中添加3個配送中心,其余數據保持不變,添加的配送中心位置信息如表2。

表2 數據集配送中心信息
本文所有實驗均使用python 3.7語言編寫,在pyCharm上編譯,并在Intel(R) Core(TM) i5-8500T CPU @ 2.10GHz 2.11 GHz Windows 10操作系統上運行。為了衡量算法的性能,文中使用了多個評價指標,其含義如表3所示。

表3 評價指標及含義
采用IDWSA對問題求解時,參數的選取對算法性能具有一定影響,本文所提出的IDWSA的參數主要有:種群大小和算法迭代次數。為最大限度發揮IDWSA的優勢,需要對算法參數進行分析,以選取最優參數。為此本文分別在種群大小為20,30,40和迭代次數為10,15,20,25,30,35時進行仿真實驗,每組實驗分別進行20次,其結果如表4,圖6,7所示。
理論上,算法所求問題最好解的質量應該與種群規模和迭代次數呈嚴格正相關,但由于WSA具有不穩定性,在同一實驗條件下對問題進行多次求解時,其所求解的質量具有一定差別,這一問題從表4也可看出。此外由圖6可以看出,隨著種群規模增加,IDWSA所求最好解的距離逐漸減小,但其降低的幅度較?。煌瑫r隨著迭代次數增加,算法所求最好解的質量有所增加,但兩者之間不具有嚴格的正相關性。

表4 不同參數下算法所求解的質量
此外,從圖7可以看出,隨著迭代次數增加,不同種群大小下算法的運算時間都呈明顯遞增趨勢。結合圖6和圖7可知,不同種群大小下所求解的質量差異較小,但在相同迭代次數下,其運算時間相差較大,且隨著迭代次數增加,算法的運算時間大幅度增長。在種群大小為20和迭代次數為30時,算法平均運算時間處于46.79 s左右,在可接受范圍內可求出問題最好解,因而在進行后續實驗時,本文選擇種群大小20、迭代次數30作為最優參數。

圖6 不同種群規模和不同迭代次數下所求最好解

圖7 不同迭代次數和不同種群規模下的求解時間
初始解的質量對后續問題的求解至關重要,為驗證本文所提混合初始化策略的有效性,本文構建了式(17)用于衡量算法迭代過程中種群的多樣性。

其中,i,j表示種群中第i和第j個個體,dij表示個體i與j之間的距離,D表示種群多樣性,其值越大表示種群多樣性越好。
分別采用混合初始化策略、DFC與隨機生成方式生成初始解,結果如表5所示。
由表5可知,采用混合初始化策略生成初始種群,其種群多樣性位于僅采用動態模糊聚類和隨機生成的初始種群之間,且其初始種群的平均最好解、平均最差解及平均解均優于動態模糊聚類和隨機生成;此外,采用混合初始化策略所求問題最終解的平均最好解、平均最差解及平均解也均優于動態模糊聚類和隨機生成,由此可表明混合初始化策略有助于算法求得質量更優的可行解。

表5 3種種群初始化方式結果對比
為驗證全局評價函數在IDWSA中的有效性,分別使用適應度函數和全局評價函數選擇個體并對其種群多樣性進行計算分析,結果如圖8,表6所示。
從圖8可以看出,采用適應度函數選擇個體,在迭代5次時,種群多樣性就已大幅度降低,表明種群中個體在迭代初期便快速聚集在某一較小區域,且在迭代過程中種群多樣性一直降低,由此可知采用適應度函數選擇個體時算法易過早陷入局部最優,從而難以求出全局最優解。采用全局評價函數選擇個體,與適應度函數相比,在迭代5次時,其種群多樣性降低幅度較小,且在迭代后期種群多樣性逐漸趨于穩定即算法逐漸收斂時,其種群多樣性仍遠高于適應度函數,由此可知采用全局評價函數選擇個體時,算法在問題解空間的較大區域內尋找可行解,使算法求得全局最優解的可能性大幅度增加。

圖8 不同迭代次數下種群多樣性
同時,由表6可知,在運算時間大致相同的情況下,采用全局評價函數所求問題最好解優于采用適應度函數所求的最好解,且其最大偏差和平均偏差均大于適應度函數,表明全局評價函數在維持高種群多樣性的同時可搜索到高質量的可行解,由此驗證了本文所提全局評價函數的有效性。

表6 兩種函數求得可行解
為驗證IDWSA求解MESPMDCTW的有效性,本文采用IDWSA進行仿真計算,同時使用文獻[8]中的多自適應粒子群優化(Multi Adaptive Particle Swarm Optimi- zation, MAPSO)算法、文獻[22]中的遺傳算法(Genetic Algorithm, GA)、文獻[23]中的混合蟻群(Hybrid Ant Colony Optimi- zation, HACO)算法以及文獻[24]中的人工蜂群(Artificial Bee Colony, ABC)算法對本文問題進行求解,并將5種算法的計算結果進行對比分析,結果如圖9所示。
由圖9可知,在20次實驗中,IDWSA所求最好解的質量與MAPSO, GA, HACO和ABC相比,具有明顯提升,其最好解距離與MAPSO, GA,HACO, ABC相比分別減少了2.25%, 13.4%, 6%,1.46%,且所求平均最好解、平均最差解以及平均解均優于MAPSO, GA, HACO和ABC,由此可以看出IDWSA在求解物資應急調度問題時可以縮短車輛行駛的總距離,從而減少運輸成本和配送時間。此外,從圖9還可以看出,IDWSA所求解的平均偏差與最大偏差均大于MAPSO, GA, HACO和ABC,這表明IDWSA所求最好解與平均解以及最差解之間的差距較大,種群在迭代后期仍保持較高的個體間差異性,從而提高算法求得全局最優解的概率。

圖9 5種算法實驗結果對比
在相同實驗條件下,IDWSA, MAPSO, GA,HACO和ABC的平均運算時間分別為46.16 s,22.6 s, 3.8 s, 5.43 s和7.7 s。IDWSA的平均運算時間遠大于MAPSO, GA, HACO和ABC,其中種群移動更新所用時間為43.52 s,占總運算時間的94.3%,由此可知,父代個體在向其引導個體移動時,為獲得高質量的子代個體,在引導個體周圍進行細致的搜索,因此耗費了大量時間,這也是后續研究中需要解決的問題。
本文針對物資應急調度背景下的MESPMDCTW問題,提出了一種改進離散鯨魚群算法(IDWSA)。綜合考慮了影響算法性能的可能因素,分別從種群初始化方式、鯨魚子代個體生成方式及個體選擇方式3個方面改進鯨魚群算法,結果表明,本文研究的IDWSA可以對車輛行駛路徑進行合理規劃,縮短車輛行駛總距離,減少物資配送時間,從而有效解決物資應急調度問題。
通過在Solomon標準測試集上進行仿真計算,結論如下:(1)根據不同種群規模和迭代次數下所求解的質量與運算時間,選取種群大小20、迭代次數30作為算法的最優參數。(2)混合初始化策略所生成初始種群的多樣性介于動態模糊聚類和隨機生成之間,且初始種群和其所求問題最終解的平均最好解、平均最差解及平均解均優于二者,驗證了混合初始化策略的有效性。(3)采用全局評價函數選擇個體,所生成子代種群的多樣性降低幅度較小,且在迭代后期其種群多樣性仍高于適應度函數,表明全局評價函數可有效維持種群多樣性。(4)與MAPSO, GA, HACO和ABC相比,IDWSA所求問題最好解的距離分別減少了2.25%, 13.4%, 6%,1.46%,并且可在一定程度上降低所求解為局部最優解的概率。