基于約束規劃的煤炭出港作業調度算法

李 偉，鄭瀾波

（武漢理工大學 物流工程學院，湖北 武漢 430063）

1 引言

我國是世界上最大的煤炭生產與消費國，龐大的消費量以及產銷在地理上的差異，導致煤炭運輸周期長、供應鏈結構復雜。港口作為煤炭多式聯運中水陸模式的轉運節點，其協調上、下游供需的作用直接影響著整條供應鏈的作業能力。據國家統計局數據，2014-2018年我國主要港口煤炭及制品吞吐量增長了10.26%，年均增長達3 811.75萬t。在我國煤炭需求旺季，港口經常出現煤炭場存供不應求的現象，因此，對港口的堆場資源進行合理調度，最大程度地滿足船舶需求是煤炭港口面臨的重要問題[1]。目前，港口調度的研究集中在集裝箱碼頭，而散貨碼頭，更確切地說對于煤炭碼頭的研究較少[2]。輸出型煤炭港口的運作模式可以分為拉動式（Pull-based）和推動式（Push-based），拉動式以客戶需求為導向進行原煤運輸和混配[2]，采用開放式的堆場，港口需要為每個到港煤堆分配堆存空間；而推動式則采用固定垛位模式，通過預測、安全庫存等方式，事先在堆場堆存一定數量的原煤來滿足客戶的需求。

拉動式港口的研究可分為精確算法和啟發式算法兩個方向，由于港口調度問題的動態性和復雜性，精確算法想要達到較高的效率往往依賴于模型簡化。文獻[3]以巴西散貨碼頭為研究對象，針對港口的主要作業環節建立了整體的混合整數規劃（MIP，Mixed Integer Programming）模型，并提出一種結合列生成和分支定價算法的精確算法。考慮到煤堆矩形在堆場的時空布局與二維裝箱問題的相似性，有學者將煤炭堆場調度問題轉化成帶特殊約束的二維裝箱問題，并用約束規劃（CP，Constraint Programming）范式對其進行建模。文獻[4]分別建立了煤炭堆場空間調度問題的MIP模型和CP模型，并通過計算實驗證明CP的求解效果更好。文獻[5]結合兩種規劃范式的優勢，將整體問題進行分解，分別用MIP和CP對港口調度主問題和堆場調度子問題進行建模，并提出一種基于benders分解的精確算法來提高模型求解效率。文獻[3]和[5]的決策時間步長均為1h，相比而言，啟發式算法能在更細的時間粒度下實現對堆場的精準調度。以典型的拉動式煤炭供應鏈HVCC（Hunter Valley Coal Chain）為例，文獻[6]提出一種基于貪婪思想的遺傳算法，該算法從全局和系統角度考慮了更全面的港口調度問題，并實現了煤炭出港作業在時空維度的連續調度。結合精確算法和啟發式算法的優點，文獻[7]以堆場為核心節點分別建立了煤炭供應的MIP和CP模型，通過滾動求解獲得初始解，并設計了大規模鄰域搜索算法對解進行改進，計算實驗表明，CP的自定義搜索策略機制讓其在求解中表現出更優異的性能。

我國大部分的煤炭港口采用推動式進行管理，其調度問題的文獻研究主要基于規則和仿真方法，而基于數學模型和優化技術的研究較少。文獻[8]將連續的堆場空間進行網格化處理，并提出了兩種物料堆存方案，使用WITNESS軟件對堆場離散時間進行仿真建模。文獻[9]結合貨物特性和堆場情況，提出三種堆存策略，并使用Plant-Simulation軟件構建系統仿真模型對堆存方案進行實驗和評價。文獻[10]應用WITNESS軟件對黃驊港的堆場進港作業系統進行仿真實驗，并建立多目標MIP模型與之對比，發現MIP的求解效果更好。

在煤炭吞吐量逐年上漲的背景下，現有的基于規則和仿真的調度方法逐漸無法滿足我國港口的日常運作要求。本文借鑒CP在拉動式港口應用的成功經驗，以數學模型和優化技術為基礎，為港口日常運作管理提供更合理、高效的調度算法。本文的研究對象為國內某大型出口配煤港口，該港的堆場分為四個部分，其中一期工程的建立時間最早，堆場結構也最復雜。因此，選取一期工程作為建模場景建立煤炭出港作業調度的CP模型，該模型充分考慮了固定垛位模式下的堆場作業特點，對于同類港口的調度問題具有普適性。為了解決問題動態性和組合性帶來的求解困難，采用分而治之的思想將問題進行分解，并設計基于變鄰域搜索的求解算法優化求解結果，從而實現求解效率和精度的平衡。

2 出港作業調度的約束規劃模型

2.1 問題描述

一期工程共有6個條形堆場和8臺取料機，堆場的空間布局如圖1所示。每個條形堆場中坐落著多個垛位，垛位在堆場中占據一定的長度并存有一定數量的原煤。每個垛位的位置、煤種和存量信息已知，相鄰垛位為防止污染會留有一段距離。取料機位于條形堆場兩側的軌道上，取料時，取料機沿著軌道行走到指定垛位旁，通過傳送帶將原煤轉運至裝船機，最后由裝船機裝載到靠泊船的對應艙位。

圖1 堆場的空間布局示意圖

客戶到港前，會先向港口發出提名（Nomination）以告知船的預計到達時間和需求，一艘船一般需要一至三種合同煤。合同煤是指發熱量、灰度等煤礦特性達到一定質量要求的煤種，由于推動式煤炭港口只存有數量有限的原煤，因此需要在裝船過程中通過原煤混裝形成客戶需求的合同煤。每種合同煤一般由一至兩種原煤按比例混合而成，稱為配煤方案。此外，為了保持船體在水中的物理平衡，裝船前還需要預先設定不同艙位的作業順序。配煤方案選擇和裝船順序可通過將問題轉化為網絡流模型求得[1]，本文模型不涉及配煤方案的選擇且假定裝船順序已知。一艘船的裝船任務實際上就是按順序分輪、定量地裝載原煤，每個輪次指定了所需原煤的種類和數量，由于單個垛位的存量有限，一種原煤時常需要在多個垛位取料。

圖2是一個裝船作業的過程示例，船舶v的裝船作業分為3個輪次，輪次U1需要裝載10 000t原煤C1，存有 C1的垛位為N105、N408和N604，實際作業在8：00到9：40從N105取料10 000t；輪次 U2需要20 000t原煤C2，由于單個垛位存量有限，因此從多個垛位取料完成該輪次作業；最后，輪次U3從唯一存有C3的N201取料完成作業。每次作業都是一個決策單元，稱為裝船單元，裝船作業調度就是在船舶泊位和需求已知的情況下，為每個裝船單元分配合理的取料時間、取料量和取料設備。由于堆場中存放的原煤種類較多且港口設備之間存在復雜的可達性約束，多條船舶同時在港作業可能會搶奪堆場的場存和設備資源，造成裝船作業的延遲。裝船作業的目的就是科學合理地調度堆場資源，讓所有船能盡早完成作業離港。

圖2 船舶v取料裝船作業示意圖

2.2 模型假設

（1）所有到港船只的需求和停靠泊位已知，且所有決策不會超出預定的規劃期。

（2）堆場中垛位的位置、煤種和存量信息確定且已知，場存充足。

（3）單次取料過程由取料機、傳送帶和裝船機組成一條整體的作業線，其作業瓶頸為取料機，故以取料機能力為標準計算裝船單元的取料時長。

（4）由于CP模型只能處理整數，故對數據進行取整處理，煤炭數量都取整為t，時間取整為min，距離取整為m，數據步長均為1。

2.3 符號說明

（1）集合

V為船舶集合，按預計到達時間排序；

U為輪次集合，按所屬船舶的預計到達時間及裝船順序排序；

S為裝船單元集合，按所屬輪次排序。

（2）參數

etav為船舶v的預計到達時間，min；

Du為輪次u的取料量（時長），min；

vs為裝船單元s所屬的船舶；

us為裝船單元s所屬的輪次；

ns為裝船單元s的取料垛位；

pn∈{0,...,1 300}為垛位n在條形堆場的位置，測算起點為距離泊位最遠的堆場邊界，m；

Rn為垛位n可用的取料機集合；

Qn為垛位n的場存量，t；

speedR=30為取料機的行走速度，m/s；

rateR=100為取料機的作業效率，t/min，所有設備的行走速度和取料效率一致。

（3）決策變量

rs∈Rns為裝船單元s使用的取料機；

ts∈為裝船單元s的取料開始時間，其中T為規劃期的時間終點，min；

ds∈為裝船單元s的取料時長，由于設備的取料效率一致，為方便建模，用取料時長代表取料量，min。

（4）輔助變量

依托鄉土體育資源開展的陽光體育活動，改變了以往學生疲于應付，被動活動的局面，變“要我運動”為“我要運動”，使學生在文明和諧的氛圍中陶冶情操，鍛煉身體，讓學生轉瞬即逝的青春留下照亮終生的火花，讓學生的人生充滿了陽光，使他們無論是處于順境還是逆境，都能樂觀而不消極。

departv為船舶v完成作業任務離港的時間，即該船最后一個裝船單元完成作業的時間，min。

2.4 約束規劃模型

裝船時間約束：裝船單元的取料開始時間不能早于所屬船的預計到港時間。

裝船順序約束：同一艘船分屬不同輪次的裝船單元需要按順序裝船。

船舶作業不重疊約束：同一艘船裝船單元的作業時間不能重疊。

垛位作業不重疊約束：在同一個垛位取料的裝船單元的作業時間不能重疊。

取料量約束：同屬一個輪次的裝船單元的取料量（時長）之和等于該輪次的總取料量（時長）。

垛位場存約束：任一垛位的取料量不能超過其場存量。

取料機防撞約束：條形堆場之間一般只有一條軌道，且按假設每條軌道只有一臺取料機作業，而堆場1和堆場2之間有兩條軌道，R1位于靠近堆場2的軌道上，R2靠近堆場1，當設備在各自靠近的堆場取料時，不會對彼此造成影響，否則將對另一臺設備的行走產生限制。如圖3所示，R2正在垛位N202取料，則位于其右側的R1無法越過R2到左側的N101和N201進行作業。

圖3 取料機R2在N202取料時對R1行走產生限制

先考慮R2限制R1行走的情況，假設裝船單元s使用設備R2在堆場2作業，則可能會對R1的作業產生影響，對任意一對先后由R1作業的裝船單元t1、t2需要至少滿足以下條件之一：t1、t2的設備移動發生在s作業之前或之后，即R1在時間上不會和R2產生沖突；t1、t2位于堆場的同一側，即R1在空間上不會和R2產生沖突；在裝船單元s作業開始之前或作業完成之后，R1有足夠的時間從t1的垛位移動到t2的垛位，同理可得R1限制R2行走的情況。

為了計算目標函數值，給出派生變量departv的計算公式，裝船作業調度的目標是讓所有船盡早完成作業離港，因此，目標函數為最小化船舶離港時間departv的累和。

3 基于變鄰域搜索的求解算法

在上述CP模型中，決策變量ts和ds共同構建了模型中大量的不重疊約束，導致一次性求解完整模型非常困難。因此，考慮將求解過程進行分解，先確定決策變量ds的取值，即裝船單元的取料方案，再求解完整模型，并設計一種改進的變鄰域搜索（VNS，Variable Neighborhood Search）算法優化取料方案。VNS是一種基于多鄰域結構的元啟發式算法，已被用于求解多個經典的組合或全局優化問題[12]，其主要思想是在搜索過程中系統地改變鄰域結構，從而實現局部收斂和全局搜索的平衡。

3.1 算法主要組成部分

VNS只提供了求解問題的一般框架，本文根據問題結構，將概率和延遲思想引入到當前解更新的過程中，設計的VNS算法主要由五個部分組成：鄰域結構的設計，這是基礎和關鍵；初始解與鄰域解的生成，實現解空間內的隨機搜索；評價值計算與接受準則，用來更新當前解；延遲策略，實現解的局部收斂；終止條件。

3.1.1 鄰域結構的設計。當有多個裝船單元同時作業時，不同裝船單元可能會搶奪堆場資源，造成作業延遲，鄰域結構的設計思路是改變部分裝船單元的取料量，通過調整不同輪次的取料方案消除作業沖突。應用VNS時，為了避免消耗過多計算時間，鄰域結構的數量不宜過多，建議為2[13]。本文設計了N0和N1兩種鄰域結構，N1的結構由N0的搜索結果決定。

（1）鄰域結構N0。假設裝船單元的集合為S，裝船單元i對應的取料量為di，則當前解可定義為，從 X1中 隨 機 選 擇 K1（如個變量進行擾動（改變變量的取值），固定其他變量取值，隨機生成鄰域解

（2）鄰域結構N1。給定集合U′?U，假設當前解中 隨 機 選 擇 K2（如個變量進行擾動，要求被選變量所屬的輪次u∈U′，固定其他變量的取值，隨機生成鄰域解

（3）N0和N1的關系。算法首先在鄰域結構N0內搜索，如果當前解有所改進，說明改變部分裝船單元的取料方案可以改進解的質量，將這部分裝船單元所屬的輪次U′記錄下來，進入鄰域結構N1，將其他輪次的取料方案固定，僅在更小的輪次U′范圍內搜索鄰域解。

3.1.2 初始解與鄰域解的生成。初始解與鄰域解的生成有兩個要求，一是快速求解，二是保證解生成的隨機性。考慮將CP模型中含有大量邏輯運算的約束（2）、（3）、（4）、（7）、（8）移除，并將最小化目標改為求可行解，形成簡化的CP模型來求解取料方案。由于移除的約束不涉及決策變量ds，所以可以保證解的可行性。如前所述，鄰域解的生成依賴于擾動變量的選擇，假設當前解為X1，如果擾動變量選擇不當，則鄰域解X2可直接由約束（5）推出，此時X1=X2。為了保證鄰域解的隨機性，要求每次生成的鄰域解滿足，否則重新生成，直至滿足條件。

簡化的CP模型使用chuffed求解器的自由搜索（FS，Free Search）策略求解模型，FS的自適應搜索算法能進一步加快求解速度，并保證求解結果的隨機性[11]。

3.1.3 評價值計算與接受準則。本文使用式（10）表示的目標函數值評價解的質量。用解表示的取料方案實例化CP模型中的取料時長變量ds，求出目標函數值作為解的評價值。在VNS的基本框架中，算法只接受改進的解來更新當前解。為了提高算法的全局搜索能力，本文借鑒模擬退火（SA，Simulated Annealing）的Metropolis準則，以一定概率接受較差的解。給定比例系數α，假設評價函數為 f()x，當前解為i，則接受鄰域解 j的概率為：

3.1.4 延遲策略。由于鄰域解的產生具有隨機性，因此每次當前解更新后，選擇相同的變量繼續擾動，尋找局部范圍內的更優解。循環迭代，直至解在最近m次局部搜索沒有得到改進，說明已經達到局部收斂，跳出循環進入下一步。

3.1.5 終止條件。使用兩種條件判斷算法是否終止運行：一是達到最大的迭代次數I，二是解足夠接近全局最優。引入文獻[2]中的理想離開時間作為問題的理論下界，理想離開時間（EPDT，Earliest Possible Departure Time）是指船v在港口資源無限的情況下最早完成作業離港的時間，即船的預計到達時間加上其裝船單元的取料時長：EPDTv=etav+

3.2 算法框架

給定VNS的最大迭代次數I，VNS算法步驟如下：

Step 1 將場存數據Q和需求數據D載入到主程序。

Step 2 設置在鄰域結構N0和N1內執行鄰域搜索的迭代次數為n0和n1。設置k=0，i=0，n=0，其中k、i和n分別表示第k鄰域結構、第i次迭代和第n次鄰域內迭代。

step 3 生成初始解X0，設置當前解X=X0，最優解Xb=X0。

step 4 設置當前鄰域結構為Nk，X執行迭代次數為nk的鄰域搜索，在每次迭代中，使用接受準則判斷是否接受鄰域解，若接受，根據延遲策略對X執行局部搜索，轉step5；否則，轉step7。

step 5 如果X優于Xb，則更新Xb，轉step6；否則，轉step7。

step 6 設置n=-1、k=k+1，如果k＞1，則設置k=0。

step 7 更新參數i、n，判斷是否達到終止條件，若是，則輸出Xb，結束程序；否則，轉step4。

4 數值實驗

為了驗證模型和算法的性能，本文根據港口歷史數據生成了10組測試用例，每組用例的船舶數量為2～7，每艘船舶需求量為10 000～35 000t，輪次數量為6～25。圖4是一個4艘船舶、11個輪次的測試用例求解的甘特圖，每個矩形代表一個裝船單元，矩形長度表示取料時長，標注分別表示取料垛位和設備。

本次實驗在64位Windows10操作系統，intel?雙核2.30GHz處理器，6G運行內存的電腦上完成。VNS算法使用Python3.6語言編寫，VNS的參數設置為：最大迭代次數I=55，鄰域N0內迭代次數n0=5，鄰域N1內迭代次數n1=3，局部搜索參數m=3。CP模型使用專用的約束規劃語言Minizinc2.3.2建模[11]，使用chuffed求解器求解。為了提高模型的求解效率，模型變量采用分組求解的策略，先從最重要的變量取料開始時間ts開始分支求解，之后依次確定取料時長ds以及取料設備rs。其中取料開始時間優先嘗試可行域內的最小值進行搜索，取料設備選擇可行域最小的變量先求解。由表1可以看出，制定求解策略后的模型求解效率得到大幅提升，CP的定制求解策略機制有利于充分發揮建模者的先驗知識，幫助求解器更快找到好的解[7]。

圖4 測試用例的甘特圖示例

表1 制定求解策略對CP模型的影響

為了驗證求解算法的有效性，將純CP模型的求解結果和加入VNS的求解算法的求解結果進行對比，并使用前文所述的理論下界評價解的質量，求解時間設定為600s。從表2可看出，CP在數據規模較小的用例中求解效果好，用例1、2、3、4均能求得最優解，且僅有用例4的耗時比VNS長。在其他數據規模較大用例中，CP無法在限定的時間內求出結果，VNS則可以求得滿意解，其在用例5、6的求解結果與理論下界的差距不到1%，其他用例的差距不到10%。

綜合實驗結果可知，對于小規模算例，純CP的求解效果更好，但當數據規模較大時，CP已無法在有效時間內求解，而VNS則能得到較優的滿意解。以上分析表明，基于VNS的求解算法是一種求解裝船作業調度的有效算法。

表2 CP與VNS的對比結果

5 結語

本文針對固定垛位模式下的煤炭出港作業調度問題，設計了結合CP和VNS的求解算法，算法求解分為兩步，先確定裝船的取料方案，再求解完整模型。將VNS引入到取料方案的優化中，其中當前解的更新采用Metropolis準則以一定概率接受劣解，并設計了一種延遲策略在每次解更新后加速算法的局部收斂。該方法綜合考慮了出港作業調度問題的組合性和動態性，為提高煤炭港口的運作管理效率和水平提供了一種新的思路。