薛 鋒 ，劉泳博 ，戶佐安 ，陳逸飛

（1. 西南交通大學(xué)交通運輸與物流學(xué)院，四川 成都 611756；2. 西南交通大學(xué)唐山研究生院，河北 唐山 063000；3. 西南交通大學(xué)綜合交通大數(shù)據(jù)應(yīng)用技術(shù)國家工程實驗室， 四川 成都 611756）

路網(wǎng)的車流分配及徑路安排是鐵路運輸組織的重要內(nèi)容，其優(yōu)化問題主要是指在現(xiàn)有的路網(wǎng)環(huán)境、設(shè)施設(shè)備及貨流條件下，合理安排車流徑路，使車流在路網(wǎng)中的運營指標到達最優(yōu). 鐵路車流分配優(yōu)化能夠在滿足運輸需求的同時，確保路網(wǎng)中各設(shè)施設(shè)備協(xié)調(diào)配合，合理使用車流徑路，進而提高整個路網(wǎng)的運輸效率.

目前，不少學(xué)者對鐵路車流分配及徑路優(yōu)化問題做過研究，并取得了豐碩的成果. 紀麗君等[1-3]基于多商品流理論建立了0-1 整數(shù)規(guī)劃車流分配與徑路優(yōu)化模型；嚴余松等[4]通過機會約束理論與車流量波動，構(gòu)建了列車編組計劃與車流徑路綜合優(yōu)化模型；王文憲等[5]利用0-1 非線性規(guī)劃模型來解決重載直達車流的分配問題；Upadhyay 等[6]基于動態(tài)規(guī)劃理論優(yōu)化了列車徑路模型；Bornd?rfer 等[7]從戰(zhàn)略層面出發(fā)，在鐵路網(wǎng)建立了以運行時間和延誤最小為目標的非線性優(yōu)化模型；Yaghini 等[8]使用混合整數(shù)規(guī)劃模型求解列車徑路問題；Khaled 等[9]以鐵路網(wǎng)為基礎(chǔ)，構(gòu)建了以系統(tǒng)總費用最小為目標函數(shù)的列車徑路優(yōu)化模型. 鐵路車流的分配屬于組合優(yōu)化中的NP-hard 問題，目前的研究大多采用CPLEX、LINGO 等進行求解[10]，雖有學(xué)者采用現(xiàn)代啟發(fā)式算法，如遺傳算法、模擬退火算法等，但這些算法在求解的質(zhì)量和算法參數(shù)控制方面各有差異，當遇到難約束問題時，算法初始解難以產(chǎn)生，影響算法的求解效率. 2014 年Mirjalili 等[11]提出灰狼優(yōu)化算法（grey wolf algorithm，GWO）. 該算法通過模擬狼群捕食獵物過程求解優(yōu)化問題，具有原理簡單、可調(diào)參數(shù)少、容易實現(xiàn)且全局搜索能力強的優(yōu)點，在函數(shù)尋優(yōu)、機器學(xué)習(xí)等方面得到了有效應(yīng)用[12]，但目前很少有學(xué)者用該算法求解鐵路車流分配及徑路優(yōu)化問題.

本文在文獻[2]的基礎(chǔ)上，引入虛擬弧和動態(tài)多段映射懲罰函數(shù)對模型進行改進，優(yōu)化其約束條件，解決模型中的難約束問題，使模型在應(yīng)對不可行流時更加靈活. 同時，運用佳點集理論對傳統(tǒng)灰狼算法初始種群的生成進行改進，并設(shè)計收斂因子的非線性變化策略，形成一種改進的灰狼算法（improved gray wolf algorithm，IGWO），用于求解基于動態(tài)罰函數(shù)的鐵路車流分配及徑路優(yōu)化模型.

1 鐵路車流分配的基本模型

1.1 基本模型的描述

鐵路車流分配問題與多商品流問題類似，即鐵路的每股車流視為一種商品，那么鐵路網(wǎng)中所有車流的分配問題可以將其視為多種商品的網(wǎng)絡(luò)配流問題[1]. 鐵路車流分配問題的本質(zhì)是在滿足路網(wǎng)約束條件的前提下，以最小的成本將車流分配到路網(wǎng)中. 因此，基于多商品流的鐵路車流分配及徑路優(yōu)化基本模型通常以貨物走行公里為目標函數(shù)，結(jié)合鐵路路網(wǎng)的線路能力要求、車流不可拆散原則等構(gòu)建約束條件. 圖1 是由5 個車站構(gòu)成的簡化網(wǎng)絡(luò)，以此為基礎(chǔ)描述鐵路車流分配及徑路優(yōu)化基本模型的特點. 假定技術(shù)站1 到技術(shù)站5 的車流量OD （origindestination）為8 （單位：車），各弧段參數(shù)為(B，M)，B為弧段容量，M為弧段里程，車流的繞行距離不得超過最短距離的兩倍.

圖1 簡化路網(wǎng)Fig. 1 Simplified railway network

如圖1 所示，從技術(shù)站1 到技術(shù)站5 的路徑有4 條（1→2→5，1→3→5，1→5，1→4→5），其中，路徑1→5 盡管為最短路徑，但不滿足弧段容量的約束，路徑1→2→5 滿足弧段容量約束，但不滿足車流的繞行距離的限制. 因此，可行路徑為1→3→5，1→4→5. 考慮以貨物走行公里最小為目標，最優(yōu)的車流徑路為1→3→5. 需要注意的是，在車流分配的過程中，必須滿足同一支車流不拆散的原則.

1.2 基本模型的假設(shè)和相關(guān)參數(shù)

考慮鐵路車流分配及徑路優(yōu)化問題的完整性和簡潔性[2]，基本模型包含如下假設(shè)：

1） 模型目標函數(shù)僅考慮貨物在弧段中能力的約束，不考慮通過支點站的走行、改編和能力限制；

2） 每支車流在通過支點站時流量不發(fā)生變化.

鐵路網(wǎng)簡化為無向圖G=(V,E) ，V為站點的集合，E為弧段的集合. 站點i，j，s，t∈V，弧段 (s,t)∈E.Ni,j為始發(fā)站i到終點站j的車流量；Bs,t為弧段(s,t)的通過能力（單位：車）；M(s,t)為弧段(s,t)的里程（單位：km）. 引入0-1 型決策變量yij,st表示車流Ni,j是否經(jīng)過弧段(s,t)，即

1.3 基本模型的構(gòu)建

根據(jù)鐵路車流的不可拆分原則，考慮線路能力、運輸服務(wù)水平等約束，以車流的總走行里程最小為目標，基于多商品流模型構(gòu)建鐵路車流分配及徑路優(yōu)化基本模型P 如式（1） ～ （5）所示.

式中：Z為路網(wǎng)中所有車輛總的走行里程；Dij為車流Ni,j在沒有能力限制下通行的最短距離；λ為最大繞行率，且 λ ≥1 .

式（1）為目標函數(shù)，表示路網(wǎng)中所有車輛總的走行里程最小. 式（2） ～ （5）為約束條件，式（2）體現(xiàn)了路網(wǎng)中的同一支車流不可拆分原則，式（3）為經(jīng)過路網(wǎng)中的每個支點站的流量守恒約束，式（4）為路網(wǎng)中的弧段能力限制的約束，式（5）為運輸服務(wù)水平約束，此處主要考慮走行距離，設(shè)置該約束可以避免少數(shù)車輛過度繞行.

2 基于動態(tài)罰函數(shù)的改進模型

鐵路車流分配按照運輸組織原則將多支車流分配到給定的路網(wǎng)中，包含具有NP-hard 特性的多商品流結(jié)構(gòu)[13]，屬于帶容量約束的徑路問題（capacitated vehicle routing problem，CVRP）[14]. 為方便對CVRP 問題的求解，需要對模型進行針對性的改進，優(yōu)化約束條件.

2.1 懲罰函數(shù)對模型的改進

懲罰函數(shù)具有結(jié)構(gòu)簡單、參數(shù)少的優(yōu)點，是解決約束優(yōu)化問題的一種有效方法. 借助懲罰函數(shù)可以將難約束吸收到目標函數(shù)中，便于問題求解[15]. 其懲罰策略為：對在無約束的求解過程中企圖違反約束的迭代點給予很大的目標函數(shù)值，迫使無約束的迭代點向可行域靠近，或者一直保持在可行域內(nèi)移動，直到收斂到原來約束最優(yōu)化的極小值點.

分析模型P 可知，路網(wǎng)的弧段能力約束式（4）是問題的難約束，參考文獻[15]，并結(jié)合懲罰函數(shù)的思想，將弧段(s,t)分為真實弧和虛擬弧，并定義新的決策變量wij,st，若車流由于弧段能力限制，只能選擇弧段(s,t)的虛擬弧運行，則wij,st取1，否則取0. 把新定義的變量加到目標函數(shù)中，構(gòu)造增廣目標函數(shù)：

式中：Z1為車流在真實弧中總的走行里程；σH(wij,st)為懲罰項，其中， σ 為懲罰力度，H(wij,st) 為懲罰因子.

當分配的車流滿足弧段的容量約束時，虛擬弧上沒有分配流量，H(wij,st) =0 ，目標函數(shù)不受額外懲罰；當分配的車流不滿足弧段的容量約束，需要在虛擬弧上分配流量，此時H(wij,st) ＞0 ，目標函數(shù)受額外懲罰； σH(wij,st) 越大，懲罰越重. 要使Z極小，H(wij,st)應(yīng)該充分小，當其趨近于0 時，Z的極小值逼近了Z1的極小值. 于是，有容量約束的鐵路車流分配及徑路優(yōu)化問題即轉(zhuǎn)化成了沒有容量約束的優(yōu)化問題. 基于懲罰函數(shù)的鐵路車流分配及徑路優(yōu)化模型Q 如式（7）所示.

2.2 懲罰項的動態(tài)更新策略

2.2.1 懲罰力度的更新

在懲罰函數(shù)的應(yīng)用過程中，需要選擇合適的懲罰力度. 若懲罰力度過小，則會大大降低懲罰項對目標函數(shù)的影響，導(dǎo)致算法搜索時間作用于不可行域，嚴重影響算法求解效率；如果懲罰力度太大，則會造成增廣目標函數(shù)Z過大，趨于病態(tài)，也會給模型的求解造成很大困難，甚至無法求解. 因此，選取合適懲罰力度的更新策略對算法的求解非常關(guān)鍵.

懲罰力度在優(yōu)化過程中的變化依賴于當前群體中可行解所占的比例[16]，本文在文獻[16]基礎(chǔ)上提出一種動態(tài)調(diào)整懲罰因子的策略. 第k次迭代過程中，懲罰力度的更新規(guī)則為

式中： ρk為在第k次迭代過程中，滿足弧段容量約束的個體所占種群數(shù)量的比值.

懲罰力度的更新規(guī)則具有以下性質(zhì)：

1） 函數(shù)結(jié)構(gòu)簡單，其參數(shù)不需要人為調(diào)整.

2） σk隨著種群中可行解的比例而變化，種群中的可行解比例高，則 σk較?。环粗?，種群中的可行解比例低，則 σk較大. 該函數(shù)的初值并不是人為給出，是由當前種群的狀態(tài)決定.

3） 在算法優(yōu)化的過程中，如果 σk較大，則會有更多的可行解進入種群，此時的 σk將會隨之減?。环粗?，如果 σk較小，則會有更多的不可行解進入種群， σk也將會隨之增大. 由此可以看出，本文設(shè)計σk的更新函數(shù)是動態(tài)的，自適應(yīng)的.

4） 在可行解的比值 ρk從0 變化到1 的過程中，前期 σk急劇減小，后期 σk減小緩慢，如此可以使算法從搜索可行解的過程中快速轉(zhuǎn)移到最優(yōu)目標函數(shù)的搜索.

5） 式（8）中的 ( 2k-1)/k是值域為[1, 2)的增函數(shù)，這樣設(shè)計可以在算法迭代的初期降低懲罰函數(shù)的影響，從而保留較多的不可行解來增加種群的多樣性，在后期增加懲罰函數(shù)的影響，從而保留較多的可行解來確保目標函數(shù)的尋優(yōu).

2.2.2 懲罰因子的更新

通過非固定階段映射懲罰力度的方法可以解決約束優(yōu)化問題，避免固定罰函數(shù)中懲罰因子難以確定的缺陷[17]. 基于該思路，結(jié)合文獻[17]的多階段映射罰函數(shù)的表達式，對優(yōu)化模型Q 中的懲罰因子修正為

式中：q為對容量約束的違背程度，其值等于虛擬弧中總車流量與真實弧中總車流量的比值，如式（10）所示；b為懲罰函數(shù)的強度大小，如式（11）所示；θ為分段映射函數(shù)，通過和真實弧的流量對比確定多階段映射罰函數(shù)各個區(qū)間的函數(shù)值，如式（12）所示.

3 灰狼算法的改進

3.1 GWO

3.1.1 狼群的社會等級

灰狼算法在迭代過程中，根據(jù)適應(yīng)度值的大小，將狼群劃分 α 、 β 、 δ 、ω4 個等級，其中， α 、 β 、 δ 為適應(yīng)度值排名前3 的頭狼，剩余的狼群為ω. 計算過程中，狼群ω實現(xiàn)整個算法的尋優(yōu)過程，3 只高等級的狼 α 、 β 、 δ 被假設(shè)擁有獲取獵物位置的潛在能力，共同指揮狼群ω的移動，然后狼群ω將信息反饋給上層的3 只狼，由他們決定是否需要更新信息. 當達到算法的迭代次數(shù)時， α 、 β 、 δ 分別對應(yīng)所求解問題的最優(yōu)解、次優(yōu)解、次次優(yōu)解.

3.1.2 狼群的捕食過程

狼群捕食過程分為包圍獵物和攻擊獵物，狼群先通過包圍獵物尋找到進行狩獵的最佳路線，然后對獵物進行攻擊，最后達到捕食獵物的目的. 狼群的捕食過程用式（13） ～ （15）進行描述.

式中：X1,k、X2,k、X3,k分別為狼群ω相對3 只頭狼 α 、 β 、δ在第k次迭代時的位置向量；Xm,k為頭狼m在第k次迭代時的位置向量，m∈{α，β，δ}；Dm為頭狼m到狼群ω間的距離；Xk為第k次迭代時狼群ω的位置向量；Am、Cm為系數(shù)向量，分別為

其中：r1、r2為隨機向量，且0≤ |r1| ≤1，0≤ |r2| ≤1；ak= 2(1-k/K)，為收斂因子，K為最大迭代次數(shù).

3.2 IGWO

3.2.1 初始種群的優(yōu)化

對群智能算法而言，初始種群的好壞影響算法的求解效率和性能，多樣性較好的初始種群有助于提高群智能算法的尋優(yōu)能力. 根據(jù)文獻[18]可知：均勻取點能夠較好地保證初始種群的多樣性，并且佳點集序列的取點方式比其他方式的取點更為均勻. 但是，GWO 在形成初始種群的過程中是通過隨機數(shù)產(chǎn)生，無法保證初始種群的個體在搜索空間中均勻分布. 因此，為克服這一不足，本文采取佳點集序列的方式對GWO 初始種群的形成進行優(yōu)化. 圖2 給出了隨機序列和佳點集序列在二維搜索空間(x,y) （ 1 00×100 ）中生成的100 個個體的初始種群分布. 從圖2 可以看出：佳點集方式形成的初始種群個體分布更加均勻，并且具有較好的多樣性.

圖2 兩種方式下初始種群的分布Fig. 2 Distribution of initial population in two modes

3.2.2 收斂因子的改進策略

在GWO 算法中，隨著算法的不斷迭代，收斂因子ak從2 線性遞減至0. 但是，線性變化的收斂因子不能很好地平衡算法的全局搜索能力和局部搜索能力[19]. 對此，本文根據(jù)余弦函數(shù)前期下降慢后期下降快的特點，設(shè)計一種基于余弦函數(shù)的收斂因子，如式（18）所示.

ak隨著迭代次數(shù)呈非線性遞減. 在迭代初期，ak下降較慢，算法保持較大的搜索步長，可以更好地尋找全局最優(yōu)解；而到了后期，ak下降較快，搜索步長快速減小，從而使狼群更加集中，能夠更加精確地尋找到局部最優(yōu)解.

4 求解改進模型的IGWO 算法設(shè)計

4.1 編碼策略

本文提出的IGWO 采用先路由后分組的策略生成多支車流的可行路徑集，為了滿足式（2）、（3），采用整數(shù)編碼，編碼規(guī)則為：將狼群中的第e個灰狼個體的位置向量Xe與各支車流序號形成映射關(guān)系，Xe= (xe,1,xe,2,···,xe,p,···,x1,n)，xe,p為第e個灰狼個體第p個OD 車流在其所在的可行路徑集中選擇的徑路編號，p= 1，2， · ··,n，n為OD 車流個數(shù)，然后按照車流的順序依次訪問所對應(yīng)的可行路徑集，可行路徑集的大小由式（5）決定，每支車流在其所在的路徑集中有且只能選擇一條徑路. 具體的編碼操作如圖3 所示，圖中：mp為第p個OD 車流的可行路徑集中的徑路個數(shù)；x1,x2,···,xn為灰狼位置分量.

圖3 編碼操作Fig. 3 Encoding operation

4.2 求解步驟

IGWO 在求解模型Q 的步驟如圖4 所示.

圖4 IGWO 在求解模型Q 的流程Fig. 4 Flowchart of IGWO solving model Q

5 算例分析

5.1 基礎(chǔ)數(shù)據(jù)

本文采用模擬的車流OD 數(shù)據(jù)，在區(qū)域鐵路網(wǎng)[2]上對模型和算法進行驗證. 模型涉及的站點和弧段的鄰接關(guān)系如圖5 所示，該路網(wǎng)由14 個支點站和20 條弧組成.

圖5 某區(qū)域簡化路網(wǎng)Fig. 5 Simplified railway network of a certain area

路網(wǎng)其他屬性參數(shù)、各區(qū)段之間的里程以及線路的容量見表1. 線路的容量為該弧段上、下行車流年通過能力的總和，20 支模擬車流OD 量及發(fā)到站見表2.

表1 路網(wǎng)相關(guān)參數(shù)Tab. 1 Railway network parameters

表2 年車流OD 量Tab. 2 Annual OD volume of cargo flow

5.2 模型求解

基于上述路網(wǎng)結(jié)構(gòu)和參數(shù)，在處理器為i7-7 700 HQ 四核2.8 GHz 的個人計算機上，使用MATLAB 2018 進行編程求解. 算法的參數(shù)設(shè)置如下：種群規(guī)模為80，最大迭代次數(shù)為100 次，最大繞行率為2.經(jīng)過37 s 的計算，各支車流優(yōu)化后的走行徑路如表3所示. 在20 支OD 車流中，按最短徑路運輸?shù)挠? 支車流，分別為N7,14、N10,3、N2,12、N14,5、N10,1、N14,5，其余車流因為線路容量的限制均發(fā)生了不同程度的繞行.

5.3 結(jié)果分析

5.3.1 弧段容量限制檢驗

為驗證表3 中的最佳配流方案是否滿足容量的約束，需要對各個弧段的流量進行統(tǒng)計. 經(jīng)計算，各弧段的通過車流量及其通過能力占用情況如表4 所示. 在表4 中，部分弧段的能力利用率高達90%以上，如弧段(2,6)、(5,8)、(13,14)、(8,10)，而弧段(4,5)、(5,6)、(9,12)的能力利用率低于50%，并且所有弧段車流的通過量均小于弧段能力的上限，滿足模型中弧段容量限制約束.

表3 車流徑路的優(yōu)化方案Tab. 3 Optimization scheme of cargo flow route

表4 區(qū)間通過流量統(tǒng)計Tab. 4 Interval traffic statistics

5.3.2 改進前后模型的結(jié)果對比分析

為驗證動態(tài)懲罰函數(shù)對模型求解的影響，應(yīng)用GWO 和IGWO 分別對改進前、后的模型進行求解，對其求解10 次，并取車流總走行公里的均值，結(jié)果如表5 所示.

表5 改進前后模型的車流總走行公里Tab. 5 Cargo flow kilometers before and after improving model

從表5 可以看出：由于模型P中包含的弧段容量約束的限制，直接用GWO 和IGWO 難以在限定范圍內(nèi)找到可行解，這說明，群智能算法在直接求解難約束問題時有著自身的局限性.

對模型進行約束優(yōu)化后，通過GWO 和IGWO可以求得滿足約束條件的可行解，這說明本文基于懲罰函數(shù)提出的約束優(yōu)化策略，克服了群智能算法在求解難約束問題時的局限性，并且可以應(yīng)用于CVRP 問題的求解.

5.3.3 算法的性能分析

1） 解的質(zhì)量

本文選取20 個OD 車流徑路的平均繞行率、選擇最短路徑的OD 個數(shù)以及目標函數(shù)作為衡量兩種算法解的質(zhì)量指標，結(jié)果如表6 所示.

表6 兩種算法的質(zhì)量指標Tab. 6 Quality metrics for two algorithms

從表6 可以看出：相比于GWO，IGWO 路徑平均繞行率下降了2.6%，車流總走行公里下降了5.2%，且IGWO 選擇最短路徑的OD 數(shù)也多于GWO.

2） 收斂性能

如圖6 所示，IGWO 在算法迭代的初期就開始收斂，在第13 代的時候達到最優(yōu)，其收斂速度與求解的結(jié)果均優(yōu)于GWO.

圖6 GWO 和IGWO 的求解示意Fig. 6 Solution illustration of GWO and IGWO

綜上，本文提出的初始種群和收斂因子的改進策略在一定程度上避免了GWO 容易陷入局部最優(yōu)的不足，使其在解空間的搜索能力上有了較大提高.

6 結(jié) 論

本文引入動態(tài)更新的懲罰函數(shù)對鐵路車流分配及徑路優(yōu)化模型中的弧段容量約束進行優(yōu)化，并設(shè)計改進灰狼算法對其進行求解，通過算例驗證，得到以下結(jié)論：

1） 通過弧段容量限制的檢驗分析發(fā)現(xiàn)，本文基于懲罰函數(shù)提出的改進模型求得的配流方案滿足弧段容量約束.

2） 通過對改進前后模型的結(jié)果分析發(fā)現(xiàn)，本文提出的對容量約束進行優(yōu)化的策略可以有效解決GWO 無法處理難約束問題的不足，使其可以在限定的條件下找到滿足約束條件的可行解.

3） 通過GWO 和IGWO 分別對改進后的模型進行求解可知，與GWO 的求解結(jié)果相比，IGWO 求得的配流方案使OD 車流的平均繞行率下降了2.6%，貨物走行公里數(shù)下降了5.2%，并且收斂速度也優(yōu)于GWO，從而體現(xiàn)了本文所提出的對GWO 改進策略的有效性.