江 峰，倪少權(quán)*,b,c，呂紅霞,b,c

(西南交通大學(xué)a.交通運輸與物流學(xué)院;b.全國鐵路列車運行圖編制研發(fā)培訓(xùn)中心;c.綜合交通運輸智能化國家地方聯(lián)合工程實驗室,成都610031)

0 引言

列車運行圖調(diào)整是根據(jù)運輸需求完善運行圖的過程，分局部調(diào)整與全局調(diào)整.局部調(diào)整在既有框架下根據(jù)增開列車需求增鋪運行線，涉及新增運行線鋪畫和既有運行線調(diào)整；全局調(diào)整則是全圖運行線的重新鋪畫.既有高速鐵路運行圖框架是列車開行方案長期優(yōu)化的結(jié)果，通常不進(jìn)行全局調(diào)整.隨著路網(wǎng)完善及運輸需求變化，高速鐵路運行圖局部調(diào)整增多，現(xiàn)有的人工推線求解方法難以保證質(zhì)量，存在很大局限性，因此有必要研究列車運行圖局部調(diào)整模型與算法.

列車運行圖調(diào)整屬列車運行圖編制(Train Timetabling Problem,TTP)范疇，是 NP-hard問題[1-3]，大量國內(nèi)外學(xué)者對相關(guān)問題進(jìn)行了研究.國內(nèi)方面，文獻(xiàn)[2]以單線區(qū)段為背景，優(yōu)化各運行線旅行速度.文獻(xiàn)[4]構(gòu)造了適用于單雙線區(qū)段的運行圖時空窗口，滾動求解各運行線.文獻(xiàn)[5]考慮上下行列車到發(fā)站屬性差異，求解單線非追蹤平行運行圖最小周期.文獻(xiàn)[6]在周期性運行圖基礎(chǔ)上，以深度搜索法添加非周期運行線.文獻(xiàn)[7]總結(jié)了基于周期事件規(guī)劃問題(PESP)的周期性運行圖模型及其應(yīng)用情況.其中文獻(xiàn)[4,6]屬啟發(fā)式方法，缺乏對所得解的優(yōu)劣評價；文獻(xiàn)[2,5]研究單線區(qū)段，不完全適用于高速鐵路雙線區(qū)段；文獻(xiàn)[6-7]研究對象為周期性運行圖.國外方面，文獻(xiàn)[3]通過給定各列車初始利潤，以總收益最大化為目標(biāo)建立整數(shù)規(guī)劃模型，對原模型進(jìn)行拉格朗日松弛，以動態(tài)規(guī)劃方法求解各運行線.文獻(xiàn)[1]根據(jù)給定列車開行方案，基于時空網(wǎng)絡(luò)鋪畫周期性運行圖.文獻(xiàn)[8]以路網(wǎng)中列車總旅行時間最小為目標(biāo)，研究貨物列車最優(yōu)徑路問題.上述研究基于拉格朗日松弛技術(shù)，適用于周期[1]及非周期列車運行圖[7-8]，并能有效評價所得解的優(yōu)化程度[1,7-8].文獻(xiàn)[9]歸納總結(jié)了不同運行圖編制模型的特點及適用情況.

上述研究為本文奠定了理論基礎(chǔ).我國高速鐵路運行圖呈非周期性，局部調(diào)整以新增列車開行方案為基本輸入，需調(diào)整的既有運行線可視為一類特殊的新增列車，從而將局部調(diào)整轉(zhuǎn)化為給定框架下的新增運行線求解問題.現(xiàn)有基于拉格朗日松弛的運行圖調(diào)整模型簡化了列車在區(qū)間運行狀態(tài)，未考慮停站引起的起停車附加時分[7-9]，無法實際應(yīng)用.本文在文獻(xiàn)[1,7-8]的基礎(chǔ)上，考慮列車停站需求，細(xì)化停站對列車區(qū)間運行時分的影響，給定各新增列車一個由理想始發(fā)時刻及始發(fā)時刻可調(diào)整度、全程停時總延長決定的始發(fā)時間域，將運行線鋪畫轉(zhuǎn)化為最短路求解問題，建立整數(shù)規(guī)劃模型并進(jìn)行拉格朗日松弛求解，根據(jù)對偶信息設(shè)計啟發(fā)式算法迭代求解運行圖最優(yōu)可行解.

1 問題描述

運行圖可由有向時空網(wǎng)絡(luò)描述[1,7-9]，其中縱軸表示車站、橫軸表示時間.將1天以分鐘為單位離散為1～q，每個車站包含1 440個節(jié)點，每條列車運行線所經(jīng)節(jié)點即可由車站及所經(jīng)節(jié)點對應(yīng)時刻確定.以G=(N,A)表示該有向時空網(wǎng)絡(luò)，其中點集N中元素代表時空域上車站所在位置，包括到達(dá)節(jié)點U及出發(fā)節(jié)點W；弧集A中元素代表相鄰節(jié)點間有向弧，為列車運行線組成部分.對每條運行線設(shè)置一個虛擬出發(fā)點σj及虛擬終到點εj，則有N=

設(shè)S為全體車站集合，給定列車集以表示列車j∈T所經(jīng)車站集合(共s站).運行圖局部調(diào)整即為在G=(N,A)中，在原有框架下，根據(jù)給定新增列車始發(fā)時間域鋪畫新增運行線的過程，如圖1所示，其中T1為新增列車.

上述問題為一定約束下G中自fj起，途徑fj+1,…,lj-1至lj止的最短路徑，可表示為，其中Nj,j∈T表示列車j所經(jīng)車站節(jié)點構(gòu)成的點子集，由列車j所經(jīng)車站及在該站的到達(dá)、出發(fā)時刻決定；Aj,j∈T表示與列車j有關(guān)的弧子集，由列車j的區(qū)間運行時分或在站停時決定，其中的求解范圍由列車始發(fā)時間域及全程停時總延長范圍決定.列車j的運行線即為與其有關(guān)的弧子集Aj中的元素，包含有：從虛擬出發(fā)點σj至列車j始發(fā)站出發(fā)節(jié)點v∈Wfj的始發(fā)??；在中間站，由到達(dá)節(jié)點u∈Ui至出發(fā)節(jié)點v∈Wi的停留?。挥绍囌緄至i+1區(qū)間內(nèi)經(jīng)由車站i出發(fā)節(jié)點v∈Wi至車站i+1到達(dá)節(jié)點u∈Ui+1的運行弧(v,u)；從列車j終點站到達(dá)節(jié)點u∈Ulj至虛擬終到點εj的終止弧(u,ε).其中，在站停留弧(u,v)的占用時間為列車j在車站i的停站時間，當(dāng)列車在站通過不停車時，其占用時間為0；運行弧(v,u)的占用時間為列車j在區(qū)間(i,i+1)的總運行時分，由列車區(qū)間運行時分決定.

2 列車運行圖局部調(diào)整模型

2.1 變量說明

2.2 模型約束

新增運行線j需滿足一定約束.對兩相鄰列車j、k而言，在任意站，其構(gòu)成弧需滿足車站間隔約束、運行線唯一性約束、列車停站約束及相關(guān)連接變量約束等.

2.2.1 車站間隔約束

(1)列車在車站i的出發(fā)時刻應(yīng)不小于該站發(fā)車間隔di，即Δ(v1,v2)≥di,v1＜v2，如圖2所示.

圖2 發(fā)車間隔示意圖Fig.2 Departure interval

(2)列車在車站i+1的到達(dá)時刻應(yīng)不小于該站到達(dá)間隔ai+1，即Δ(u1,u2)≥ai+1,u1＜u2，如圖3所示.

圖3 到達(dá)間隔示意圖Fig.3 Arrival interval

約束可表示為對任意i∈S/{1},u∈Ui+1，在Δ(u1,u2)＜ai+1內(nèi)所有節(jié)點u至多有1個被某運行線經(jīng)過，即

(3)兩運行線在區(qū)間內(nèi)不相交，設(shè)v1＜v2,u2＜u1，如圖4所示.

圖4 區(qū)間不相交約束示意圖Fig.4 No overpass in sector constraint

設(shè)區(qū)間內(nèi)列車k的速度高于列車j，約束可表示為對任意相鄰車站i和i+1，運行線j、k，圖4所示的4個節(jié)點中，至多3個被這兩條運行線經(jīng)過，設(shè)v1＜v＜v2，u2＜u＜u1，有

2.2.2 運行線唯一性約束

(4)始發(fā)弧數(shù)量約束.

(5)節(jié)點v進(jìn)入、離開弧數(shù)量約束.

(6)節(jié)點v選擇次數(shù)約束.

2.2.3 列車停站約束

列車停站約束為列車j在經(jīng)停站停時不少于該站計劃停時；在非經(jīng)停站出發(fā)時刻與到達(dá)時刻相等.設(shè)θ(v)=θ(u)+tij,i∈Scj?Sj,j∈T，Scj為列車j的經(jīng)停站，tij為列車j在車站i的計劃停時，包括：

(7)經(jīng)停站出發(fā)時刻選擇限制.

(8)非經(jīng)停站出發(fā)時刻選擇限制.

2.2.4 連接變量約束

連接變量約束為各0-1變量間相互約束.

(9)弧構(gòu)成約束.

(10)節(jié)點選擇約束.

2.3 目標(biāo)函數(shù)

各運行線始發(fā)時間域及總停時限制決定了其所有可選節(jié)點及弧的子圖Gj=(Nj,Aj)，其可行解即為在Gj=(Nj,Aj)中滿足上述約束條件的1條路徑.為評價解的質(zhì)量，對每一條運行線j初始點賦予初始利潤πj，對其相對計劃始發(fā)時刻的調(diào)整及總停時延長賦予罰數(shù)，該運行線的最終利潤即為各構(gòu)成弧的利潤之和[1].

設(shè)αj、βj為相應(yīng)系數(shù)，對每個初始弧(σ,v),v∈Wfj，設(shè)υ(v)為始發(fā)點v對應(yīng)時刻與理想始發(fā)時刻θ(fj)的差值，有其利潤為p(σ,v)=πj-αjυ(v)；對每一個停留弧(u,v),u∈Ui,v∈Wi,i∈Sj/{fj,lj}，設(shè)μ(u,v)為實際停時與計劃停時θ(u,v)的差值，其利潤為p(u,v)=-βjμ(u,v)，以pa表示運行線j中所有弧a∈Aj利潤之和，有pa=p(σ,v)+p(u,v)，以全圖所有運行線利潤最大構(gòu)建目標(biāo)函數(shù)[1].

3 基于拉格朗日松弛的啟發(fā)式算法

3.1 約束松弛

上述模型中行車組織約束屬于簇約束，數(shù)量隨問題規(guī)?？焖僭鲩L，會引起模型規(guī)模的急劇擴(kuò)大，對其進(jìn)行拉格朗日松弛可加快求解速度[1,8-9].基于文獻(xiàn)[1,8-9]所述方法，對原模型中式(1)～式(3)分別增加拉格朗日乘子后添加至目標(biāo)函數(shù)，將其轉(zhuǎn)化為

式中：C——由行車組織約束決定的與該節(jié)點不相容弧集[10]；

λC1、λC2——對應(yīng)的拉格朗日乘子.

對式(11)進(jìn)行恒等變換，即

λh——第h個松弛約束對應(yīng)的拉格朗日乘子，其含義為，①當(dāng)λh對應(yīng)約束為式(1)或式(2)，即yw時，有bh=1，任何經(jīng)過點w的運行線將被賦予1個罰數(shù)λh.②當(dāng)λh對應(yīng)約束為式(3)，即zjw時，有bh=3，任何經(jīng)過點w的列車j將被賦予1個罰數(shù)λh.

原模型松弛后，各運行線間約束關(guān)系被弱化，可由線性規(guī)劃方法快速求解[1,7,10]，所得結(jié)果即為各運行線的拉格朗日松弛解.

3.2 可行解啟發(fā)式求解

由于松弛了原模型中部分約束，LD中各運行線拉格朗日松弛解可能是實際非可行解.為求得實際可行解FD，設(shè)計下述啟發(fā)式算法：

(2)新增運行線按拉格朗日松弛解利潤大小排序確定鋪畫順序.

(3)在新增運行線始發(fā)時間域內(nèi)基于所有可行出發(fā)節(jié)點并考慮p(σ,v)、p(u,v)，選擇λh(N)=0的節(jié)點，以min(p(σ,v)+p(u,v))為目標(biāo)，基于深度搜索算法[5]求解實際利潤最大的運行線，并記錄每個車站的最大停留列車數(shù)，若大于該站到發(fā)線數(shù)量，則該運行線無實際可行解，所得解為各運行線的實際可行解.

(4)求解全圖實際利潤大于零的運行線實際利潤之和，結(jié)果記錄為FD的解.

所求得的FD為D的1個實際可行解[5,10].由于考慮了λh(N)=0的限制，因此有FD≤D，即FD為D的1個下界.

3.3 拉格朗日乘子更新

由于FD≤D≤LD，故D最優(yōu)解處于區(qū)間(FD,LD)內(nèi)，通過更新拉格朗日乘子，可優(yōu)化FD的目標(biāo)函數(shù)值，使其趨近于D最優(yōu)解[7-10].拉格朗日乘子的更新由LD中式(1)～式(3)的違反情況確定[10,11]，其步驟為：

(1)初始化各約束次梯度要素viol(λh)=-1.

(2)對各運行線構(gòu)成弧檢查是否滿足第h個約束條件，若違反該約束則對第h個約束的次梯度要素更新為viol(λh)=viol(λh)+1.

(4)根據(jù)Fisher提出的式(15)更新第h個約束對應(yīng)的拉格朗日乘子λh[10-11].

式中：ρ——給定的非負(fù)步長[10-11]，初始值設(shè)為ρ=1，當(dāng)?shù)?0次min(LD)未發(fā)生改變則減半ρ取值.

通過上述過程可優(yōu)化LD及FD的目標(biāo)函數(shù)值，從而對FD進(jìn)行迭代優(yōu)化，當(dāng)LD與FD滿足一定條件時可認(rèn)為已求得最優(yōu)可行解[1,7-10].

3.4 算法流程

Step 0初始化，記錄起始時間，設(shè)迭代次數(shù)為0，當(dāng)前最優(yōu)可行解為0，當(dāng)前最優(yōu)上界為全圖各運行線初始利潤之和.

Step 1求解LD，記錄當(dāng)前最優(yōu)上界為min(LD).

Step 2根據(jù)LD所得各節(jié)點罰數(shù)λh(N)，以3.2節(jié)所述方法求解運行線實際可行解，計算全圖各運行線實際可行解利潤和FD，記錄當(dāng)前最優(yōu)下界為max(F D).

Step 3設(shè)gap=(min(LD)-max(FD))/max(FD)?100%，判斷max(FD)與min(LD)間差值是否小于設(shè)定gap值，若是，算法終止，最優(yōu)下界所對應(yīng)的解即為最優(yōu)可行解；若否，迭代次數(shù)+1，轉(zhuǎn)Step4.

Step 4以3.3節(jié)所述方法更新拉格朗日乘子λh.

Step 5檢查算法運行時間、迭代次數(shù)是否達(dá)到設(shè)定最大值，若是，算法終止，最優(yōu)下界對應(yīng)解即為最優(yōu)可行解；若否，轉(zhuǎn)Step1.

4 算例驗證

以2015年5月京滬高鐵運行圖為例，在原圖302條運行線框架下增鋪24條G類高速運行線，新增列車信息如表1所示.

表1 新增列車信息Table 1 New trains’information

本文主要驗證模型算法的有效性，求解過程中暫不要求列車上下行數(shù)量相等.

首先以現(xiàn)有調(diào)整方法進(jìn)行實驗：以初始始發(fā)時刻推線求解，結(jié)果表明僅能增加AD1015/AD1016、AD1019、AD1023/AD1024次等5條運行線；擴(kuò)大始發(fā)時刻可調(diào)整度至60 min，可增加AD1009、 AD1015/AD1016、 AD1018/AD1019、AD1021/AD1022、AD1023/AD1024次等 9條運行線.

然后驗證本文所提出模型：給定各運行線初始利潤10 000，設(shè)αj=10、βj=20、新增列車總停時延長上限為10 min、算法迭代次數(shù)上限為200次、運行時間上限為60 min、gap為1%.列車運行參數(shù)根據(jù)實際情況確定，列車始發(fā)時間域以給定的理想始發(fā)時刻及可調(diào)整度確定.

算例對應(yīng)規(guī)模下，可逐一列出所有簇約束所包含的約束條件.為驗證拉格朗日松弛算法效率，調(diào)用ILOG CPLEX 12.5求解原模型D進(jìn)行對比驗證.相關(guān)算法以C語言實現(xiàn)，實驗在1臺Ubuntu14.04 Xeon 2.66GHz 4G內(nèi)存工作站上進(jìn)行，結(jié)果如表2所示，其中C代表CPLEX求解結(jié)果，L代表拉格朗日算法求解結(jié)果.

可見可增鋪運行線數(shù)量隨列車始發(fā)時刻可調(diào)整度增加而增加，說明列車始發(fā)時間域是影響列車運行圖局部調(diào)整的重要因素.簇約束式(1)～式(3)使得模型規(guī)模急劇增加，CPLEX雖然求解精度較高，但求解速度下降明顯；而拉格朗日算法求解時間增加較為平緩，目標(biāo)函數(shù)值、新增列車運行線平均旅速與CPLEX求解所得結(jié)果十分接近(表2)，在保證求解質(zhì)量的同時維持了較高的計算效率.統(tǒng)計不同始發(fā)時刻可調(diào)整度下拉格朗日算法所得結(jié)果如表3所示，其中始發(fā)時刻調(diào)整總時長、平均時長、最大調(diào)整時長分別為相應(yīng)始發(fā)時刻可調(diào)整度下各新增運行線較給定理想始發(fā)時刻調(diào)整值之和、各新增運行線始發(fā)時刻調(diào)整值均值、各新增運行線最大始發(fā)時刻調(diào)整度；總停時延長值為相應(yīng)始發(fā)時刻可調(diào)整度下各新增運行線較給定最短停時延長值之和.

表2 實驗計算結(jié)果Table 2 Experiment results

表3 拉格朗日算法計算結(jié)果Table 3 Lagrangian algorithm results

實驗中隨著始發(fā)時刻可調(diào)整域的增加，相較于根據(jù)給定的理想始發(fā)時刻推線，由于考慮了列車在站停時可在一定程度上延長，使得運行線求解空間增大，例如在始發(fā)時刻可調(diào)整度為60 min時，新增運行線在站停時可以在給定最短停時的基礎(chǔ)上適當(dāng)延長(總停時延長不超過10 min)，某站停時的適當(dāng)延長可能會避免后續(xù)鋪畫時與其他運行線的潛在沖突，在此基礎(chǔ)上遍歷所有的運行線可行解，增加了新增運行線的鋪畫概率，因此，相較按理想始發(fā)時刻及最短停時推線，額外增鋪了6條運行線；但同時由于全程停時的增加，引起了旅速在一定程度上的下降.

進(jìn)一步驗證京滬高鐵現(xiàn)有運行圖框架下的最大能力，將各運行線始發(fā)時刻可調(diào)整度放寬至4 h，計算滿足天窗限制(0:00-6:00)的可新增運行線數(shù).由于簇約束影響，CPLEX無法在給定時間內(nèi)得到最優(yōu)解；拉格朗日算法經(jīng)過200次迭代、509 s后由于迭代次數(shù)達(dá)到設(shè)定上限終止，共鋪畫了有效運行線18條，說明在算例框架下京滬高鐵通過能力已接近飽和，進(jìn)一步增開列車需對原有框架做出一定調(diào)整.

5 結(jié)論

針對給定框架下增鋪運行線問題，基于時空網(wǎng)絡(luò)構(gòu)建整數(shù)規(guī)劃模型，采用基于拉格朗日松弛的啟發(fā)式算法求解，實例驗證得到結(jié)論如下：

(1)采用始發(fā)時刻推線方法，算例所述情況下能額外鋪畫9條運行線，而拉格朗日算法能夠鋪畫15條運行線.

(2)除始發(fā)時刻可調(diào)整度±10 min情景外，拉格朗日算法的求解速度均優(yōu)于CPLEX.

(3)隨著始發(fā)時刻可調(diào)整度的增加，模型規(guī)模急劇增長，CPLEX所需的求解時間快速增加，而拉格朗日算法計算時間增加較平緩，計算效率優(yōu)勢明顯.

(4)始發(fā)時刻可調(diào)整度由10 min增加至60 min時，由于各運行線可行解范圍增加，可增鋪的運行線由6條增加至15條.

(5)始發(fā)時刻可調(diào)整度4 h條件下，CPLEX無法在給定時間內(nèi)求得最優(yōu)解，拉格朗日算法在現(xiàn)有框架下最多增鋪了18條運行線.

在原列車運行圖框架不變條件下，新增運行線鋪畫受既有運行圖框架限制較大，下一步將研究既有列車運行圖框架可部分調(diào)整條件下所需調(diào)整的運行線并給出調(diào)整方案.