蘭澤康，何世偉，黎浩東

(北京交通大學(xué)城市交通復(fù)雜系統(tǒng)理論與技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室，北京100044)

0 引言

當(dāng)列車出發(fā)晚點(diǎn)、運(yùn)行晚點(diǎn)、臨時(shí)加開列車或列車早點(diǎn)到達(dá)等情況發(fā)生時(shí)，都會(huì)導(dǎo)致列車運(yùn)行偏離圖定時(shí)刻，為了恢復(fù)列車安全有序運(yùn)行，需要采取列車運(yùn)行調(diào)整措施.目前列車運(yùn)行調(diào)整研究主要集中于單條線路，限制了調(diào)整的空間，若可以擴(kuò)展至一定區(qū)域的線路網(wǎng)絡(luò)，則可以充分利用其他線路上空閑的資源，達(dá)到整體的優(yōu)化效果.但鐵路網(wǎng)絡(luò)求解規(guī)模較大，不易快速求解，而列車運(yùn)行調(diào)整實(shí)時(shí)性要求較高，故需要設(shè)計(jì)高效的算法.

在鐵路網(wǎng)絡(luò)列車運(yùn)行圖方面，彭其淵等[1]在列車選擇線路固定條件下，提出采用加邊法求解；Meng等[2]建立鐵路網(wǎng)絡(luò)列車運(yùn)行調(diào)整的MILP模型，利用拉格朗日松弛算法化解成單列車的最短路徑問(wèn)題；陳雍君等[3]考慮微觀路網(wǎng)的復(fù)雜性采用了序優(yōu)化算法；Mu等[4]提出了多種啟發(fā)式方法，包括初始化候選路徑集，以及從路網(wǎng)或列車數(shù)量角度分解成小規(guī)模問(wèn)題滾動(dòng)求解的方法.可以看出研究方法大多采取不同的分解策略，達(dá)到簡(jiǎn)化求解的目的.本文首先在文獻(xiàn)[2]研究成果的基礎(chǔ)上，建立基于列車到發(fā)時(shí)刻的MILP模型，再構(gòu)建基于列車時(shí)空路徑的ILP模型，設(shè)計(jì)分支定價(jià)方法求解.

1 問(wèn)題描述與模型假設(shè)

1.1 問(wèn)題描述

相對(duì)于單條線路，鐵路網(wǎng)絡(luò)的列車運(yùn)行調(diào)整問(wèn)題，每一列列車可能有多條行駛線路可以選擇.如圖1所示，是由2條單線鐵路和3個(gè)車站組成的調(diào)度區(qū)段，標(biāo)有數(shù)字的節(jié)點(diǎn)指軌道電路的分界點(diǎn)，連接相鄰兩個(gè)節(jié)點(diǎn)的弧段指區(qū)間內(nèi)閉塞分區(qū)或車站內(nèi)股道線路.圖中矩形方框表示相應(yīng)弧段在某一時(shí)間段內(nèi)進(jìn)行維修作業(yè)，期間不可被列車占用.正常情況下，列車可按照既定的計(jì)劃有序地在該網(wǎng)絡(luò)上運(yùn)行.當(dāng)列車晚點(diǎn)進(jìn)入調(diào)度區(qū)段時(shí)，就可能與其他列車和維修作業(yè)發(fā)生沖突.假設(shè)可以變更列車運(yùn)行線路的條件下，需要重新調(diào)整列車運(yùn)行計(jì)劃，使得列車盡可能按圖定時(shí)刻到達(dá)終點(diǎn).

現(xiàn)有研究大多將車站看作是一個(gè)點(diǎn)，到發(fā)線數(shù)量約束作為模型的約束條件，如文獻(xiàn)[5].本文考慮車站內(nèi)進(jìn)路，將到發(fā)線運(yùn)用和列車時(shí)刻表一同優(yōu)化，因此可以不必單獨(dú)考慮車站能力約束，但同時(shí)也造成求解規(guī)模增大.

1.2 模型假設(shè)

假設(shè)1將列車看作是一個(gè)質(zhì)點(diǎn)，列車可選擇任意1條從始發(fā)節(jié)點(diǎn)到終到節(jié)點(diǎn)的可達(dá)路徑，包括更改運(yùn)行線路和變更站內(nèi)到發(fā)線.

圖1 鐵路網(wǎng)絡(luò)示意圖Fig.1 Sketch map of railway network

假設(shè)2不考慮停站變動(dòng)的影響.正線不允許停車，只準(zhǔn)在側(cè)線上停車，例如圖1中的2-4、15-16和9-10弧段為側(cè)線.

假設(shè)3已知維修作業(yè)的開始時(shí)間、結(jié)束時(shí)間及作業(yè)地點(diǎn).每一個(gè)維修作業(yè)所占用的弧段是緊鄰的.

2 基于列車到發(fā)時(shí)刻的模型建立

2.1 基本符號(hào)定義

給定鐵路網(wǎng)絡(luò)G=(V,A)，其中V為節(jié)點(diǎn)集合，A為弧段集合，v∈V，a∈A；G內(nèi)所有運(yùn)行的列車集合為F，f∈F，列車f有唯一索引If=0,1,…,|F|-1；列車f的始發(fā)節(jié)點(diǎn)為of，終到節(jié)點(diǎn)為df；Af為列車f可能通過(guò)的弧段集合，As為允許停車的弧段集合；A o(v)和A d(v)分別為以v起點(diǎn)和終點(diǎn)的弧段集合；Rf(a)為列車f在弧段a上的最小運(yùn)行時(shí)間分別為列車f在弧段a(a∈As)上的最小和最大停留時(shí)間；h為列車車頭安全間隔時(shí)間，即統(tǒng)一所有安全間隔時(shí)間為h；Sf為列車f最早可能始發(fā)時(shí)間；Ef為列車f圖定終到時(shí)間；K為維修作業(yè)集合，k∈K；M為足夠大的正整數(shù).

決策變量sf(a)和ef(a)分別為列車f到達(dá)和離開弧段a的時(shí)間；yf(a)表示列車f是否通過(guò)弧段a的二元變量，是則取值為1，否為0；λ(f,f′,a)為二元變量，當(dāng)列車f′在f之后通過(guò)弧段a時(shí)，取值為1，否為0.

2.2 約束條件和目標(biāo)函數(shù)

(1)流平衡約束.

式(1)和式(2)分別保證列車在始發(fā)節(jié)點(diǎn)和終到節(jié)點(diǎn)的流平衡，式(3)保證列車在中間節(jié)點(diǎn)的流平衡.

(2)最早發(fā)車時(shí)間及時(shí)間變量相關(guān)約束.

式(4)為列車最早允許出發(fā)時(shí)間約束；式(5)表示列車離開前一弧段與到達(dá)下一相鄰弧段的時(shí)間相等；式(6)和式(7)表示當(dāng)列車f不占用弧段a時(shí)，則sf(a)和ef(a)均取值為 0.

(3)最小運(yùn)行時(shí)間和停站時(shí)間約束.

式(8)為最小運(yùn)行時(shí)間約束；式(9)和式(10)分別將最小停站時(shí)間和最大停站時(shí)間附加到弧段運(yùn)行時(shí)間上，以滿足停站時(shí)間約束.

(4)安全間距約束.

式(11)表示當(dāng)列車f和f′都占用弧段a，且f′在f之后占用時(shí)，f′到達(dá)弧段a時(shí)間與f離開弧段a的時(shí)間應(yīng)間隔h，(f,f′)?F表示F中滿足If′＞If的兩列車；式(12)與式(11)同理；式(13)表示當(dāng)列車f和f′都占用弧段a時(shí)，有λ(f,f′,a)+λ(f′,f,a)≥ 1，結(jié)合

式(14)，得到λ(f,f′,a)+λ(f′,f,a)=1，即f在f′之前或之后占用a；式(15)和(16)表示列車f或f′不占用a時(shí)，λ(f,f′,a)和λ(f′,f,a)均取值為 0.

(5)維修作業(yè)約束.

式中：Ak為維修作業(yè)k所占用弧段集合分別為維修作業(yè)k的開始時(shí)間和結(jié)束時(shí)間；zf(k)是二元變量，當(dāng)列車f在維修作業(yè)k結(jié)束之后通過(guò)弧段a(a∈Ak)時(shí)，取值為1，f在維修作業(yè)k開始之前通過(guò)a或不通過(guò)a時(shí)均取值為0.

式(17)表示當(dāng)列車f在維修作業(yè)k結(jié)束之后通過(guò)弧段a時(shí)，應(yīng)滿足sf(a)≥moteknd；式(18)表示當(dāng)列車f在維修作業(yè)k還未開始時(shí)通過(guò)弧段a，或者不通過(guò)a，應(yīng)滿足ef(a)≤motstartk.

(6)目標(biāo)函數(shù)及相關(guān)約束.

3 基于列車時(shí)空路徑的模型建立

為刻畫列車在鐵路網(wǎng)絡(luò)上運(yùn)行的時(shí)空路徑，如圖2所示，以鐵路網(wǎng)絡(luò)為空間平面，并加上時(shí)間軸，形成三維的立體空間.以固定時(shí)間長(zhǎng)度離散化時(shí)間軸，為平衡求解精度和時(shí)間，本文取固定時(shí)長(zhǎng)為1 min.用時(shí)間拓展的時(shí)空聯(lián)弧(v1,v2,t,s)連接節(jié)點(diǎn)v1和v2，t、s分別表示到達(dá)節(jié)點(diǎn)v1和v2的時(shí)間，順次連接列車經(jīng)過(guò)的所有節(jié)點(diǎn)，即形成1條時(shí)空路徑p，p在空間平面上的投影即為列車在鐵路網(wǎng)絡(luò)中實(shí)際通過(guò)的物理路徑.

圖2 列車在鐵路網(wǎng)絡(luò)中運(yùn)行的時(shí)空路徑Fig.2 The time-space path of train on a railway network

假設(shè)開始調(diào)整的時(shí)間為0，調(diào)整時(shí)段長(zhǎng)為T.令Pf為列車f所有可能的時(shí)空路徑集.cf(p)為p(p∈Pf)的費(fèi)用，若p的終到時(shí)間為則cf(p)可按式(23)計(jì)算.

xf(p)為列車f是否選擇p(p∈Pf)的二元變量，若選擇則取1，否取0.為滿足列車安全間距要求，需要保證任意弧段在時(shí)間段[τ,τ+h-1](τ=1,2,…,T-h)內(nèi)至多被1列列車占用，定義Pf中在該時(shí)間段內(nèi)占用弧段a的時(shí)空路徑集為Pf(a,τ,τ+h-1).建立基于列車時(shí)空路徑的模型M2為

式(24)為目標(biāo)函數(shù)，表示最小化所需費(fèi)用；式(25)表示列車f選擇Pf中唯一的1條時(shí)空路徑；式(26)確保被選擇的任意2條時(shí)空路徑滿足安全間距要求.

模型M2與文獻(xiàn)[5]模型主要有以下幾點(diǎn)不同：①本文定義了每1列列車的可選路徑集，而文獻(xiàn)[5]中只定義了1個(gè)可選路徑集；②由于本文考慮了車站內(nèi)軌道電路設(shè)置，不必單獨(dú)考慮到發(fā)線數(shù)量約束，且不同于文獻(xiàn)[5]以車站到發(fā)間距描述安全間隔，本文以任一弧段不被2列列車同時(shí)占用進(jìn)行描述.

4 求解模型M2的分支定價(jià)算法

分支定價(jià)算法是一種將列生成算法和分支定界算法相結(jié)合的組合算法.由于列生成算法所得到的解一般不是整數(shù)解，需要借助分支定界算法對(duì)非整數(shù)變量分支以獲取整數(shù)解，每個(gè)子結(jié)點(diǎn)仍采用列生成算法求解.以下是基于目標(biāo)函數(shù)為最小化，0-1整數(shù)規(guī)劃問(wèn)題為前提進(jìn)行分析的.

4.1 求解線性松弛模型的列生成算法

上述模型M2結(jié)構(gòu)簡(jiǎn)單，但缺點(diǎn)在于列車可選時(shí)空路徑數(shù)量是非常巨大的，無(wú)法一一枚舉，因此采用列生成算法定價(jià)生成路徑的方法.首先將M2分解成限制主問(wèn)題RMP和|F|個(gè)子問(wèn)題，列車f對(duì)應(yīng)的子問(wèn)題為SPf.令P′f表示Pf的子集，替換M2中的Pf，且將式(27)線性松弛，即構(gòu)成RMP.在RMP求解之后，可以得到對(duì)偶最優(yōu)解(μ,λ)，其中μ=(μf)，λ=(λ(τ,a))分別為式(25)和式(26)的對(duì)偶變量，則SPf求解模型為

式中：AP為所包含的弧段集合為上弧段a所對(duì)應(yīng)的時(shí)空聯(lián)弧在時(shí)間軸上投影所包含的時(shí)間點(diǎn)集合.

將λ(τ,a)看作是在τ時(shí)刻占用弧段a的價(jià)格，結(jié)合時(shí)空路徑的形成，子問(wèn)題可轉(zhuǎn)化成求每1列列車的最短路問(wèn)題，本文最短路算法采用的是文獻(xiàn)[2]中離散時(shí)空標(biāo)號(hào)更新算法.若SPf的目標(biāo)值，則在P′f中添加當(dāng)前最短路，否則不需要添加.RMP和子問(wèn)題不斷迭代求解，直至所有子問(wèn)題都沒(méi)有添加新的路徑時(shí)終止，可獲得模型M2線性松弛后的最優(yōu)解，它是模型M2的下界值.

4.2 獲取整數(shù)解的分支定界算法

分支求解過(guò)程較為耗時(shí)，需采用加快算法收斂的分支策略和結(jié)點(diǎn)搜索策略.本文采用的分支策略為偽費(fèi)用分支[6-7](Pseudocost Branching).其主要思想是給每個(gè)變量定義相應(yīng)的評(píng)估值，每次選擇出分支后最有可能使得下界值更快上升變量進(jìn)行分支.對(duì)于一個(gè)非整數(shù)變量xi，若對(duì)其分支可得目標(biāo)值對(duì)xi變化的敏感度為

式中：c(q)為父結(jié)點(diǎn)的目標(biāo)值分別為對(duì)變量xi上、下取整之后得到子結(jié)點(diǎn)的目標(biāo)值.

由于xi可能多次被分支，記的累計(jì)值分別為累計(jì)次數(shù)分別為若子結(jié)點(diǎn)的RMP問(wèn)題不可行則不予累計(jì)，最終xi的評(píng)估值按式(30)計(jì)算，每次選擇評(píng)估值最大的變量進(jìn)行分支.由于在分支樹的開始階段，大部分變量沒(méi)有被分支過(guò)，無(wú)法計(jì)算評(píng)估值，可采用強(qiáng)分支(Strong Branching)方法對(duì)評(píng)估值進(jìn)行初始化，這里不在贅述，可以參見文獻(xiàn)[6-7].

式 中 ：score(s+,s-)=(1-μ)?min(s+,s-)+μ?max(s+,s-)，μ=1 6.

分支中上取整對(duì)目標(biāo)值的影響較大，也可能是無(wú)解的，而下取整對(duì)目標(biāo)值的影響較小，故本文僅對(duì)上取整得到的子結(jié)點(diǎn)用列生成算法求解，而下取整的子結(jié)點(diǎn)添加取整約束后用對(duì)偶單純形算法求解以獲取目標(biāo)值，進(jìn)一步加快算法收斂.為便于后續(xù)論述，記上取整得到的子結(jié)點(diǎn)為右結(jié)點(diǎn)，而下取整得到的子結(jié)點(diǎn)為左結(jié)點(diǎn).結(jié)點(diǎn)搜索采用最佳優(yōu)先(Best-node First)搜索策略，即每次選擇目標(biāo)值最小的結(jié)點(diǎn)，使得下界可以更快地上升.

4.3 處理眾多約束條件的行生成方法

4.4 算法流程

算法流程中涉及到的符號(hào)定義如表1所示.

Step 1獲得初始可行解和初始化.

按照列車最早可能出發(fā)時(shí)間，依次對(duì)列車按照最短路算法進(jìn)行規(guī)劃，獲得可行解的目標(biāo)值作為UB的初始值，令LB=0.創(chuàng)建1個(gè)結(jié)點(diǎn)，以初始可行解建立該結(jié)點(diǎn)的RMP問(wèn)題，將其加入ANL.

Step 2終止條件.

表1 算法流程中的符號(hào)定義Table 1 Notations for the algorithm flow

Step 3結(jié)點(diǎn)選擇.

以最佳優(yōu)先搜索策略從ANL中選擇1個(gè)結(jié)點(diǎn)作為當(dāng)前結(jié)點(diǎn)n，將n從ANL中移除.若n是左結(jié)點(diǎn)，用對(duì)偶單純形算法求解RMP(n)，否則用列生成算法求解RMP(n).

Step 4結(jié)點(diǎn)操作判斷.

若RMP(n)不可行或Z(n)≥UB，轉(zhuǎn)Step2；若B(n)=1，更新UB，轉(zhuǎn)Step 2；若B(n)=0，轉(zhuǎn)Step5.

Step 5分支.

按照偽費(fèi)用分支策略分支，生成右結(jié)點(diǎn)和左結(jié)點(diǎn)，并在RMP(n)的基礎(chǔ)上添加分支變量的取整約束，構(gòu)建相應(yīng)子結(jié)點(diǎn)的RMP問(wèn)題，將2個(gè)子結(jié)點(diǎn)加入ANL，轉(zhuǎn)Step2.

5 算例分析

5.1 算例數(shù)據(jù)

本文算例設(shè)計(jì)鐵路網(wǎng)絡(luò)共有76個(gè)節(jié)點(diǎn)和85個(gè)弧段，總長(zhǎng)度137.6 km，如圖3所示.維修作業(yè)和部分列車相關(guān)數(shù)據(jù)分別如表2和表3所示，列車晚點(diǎn)進(jìn)入?yún)^(qū)段的時(shí)間分布在[0,20]min范圍內(nèi).計(jì)劃運(yùn)行圖中所有列車均未通過(guò)交叉渡線(圖3中54-58，56-57，44-47，45-46)，而列車運(yùn)行調(diào)整后，列車可能通過(guò)交叉渡線從一條線路換軌至另一條線路.

表2 維修作業(yè)相關(guān)數(shù)據(jù)Table 2 The data for maintenances

圖3 算例的鐵路網(wǎng)絡(luò)Fig.3 The railway network for numerical experiment

表3 部分列車相關(guān)數(shù)據(jù)Table 3 The data for some trains

5.2 算法性能分析

為分析本文分支定價(jià)算法的性能，用2種方法分別求解算例，測(cè)試計(jì)算機(jī)參數(shù)為Intel Pentium CPUB9602.20GHZ，6GB內(nèi)存.求解結(jié)果如表4所示.

(1)IP_GUROBI.運(yùn)用求解器GUROBI6.5.0求解基于列車到發(fā)時(shí)刻的模型.

(2)BAP_C#.運(yùn)用C#語(yǔ)言實(shí)現(xiàn)分支定價(jià)方法求解模型M2，以GUROBI作為L(zhǎng)P求解器.

求解結(jié)果表明，隨著列車數(shù)量的增多，即求解規(guī)模的增大，IP_GUROBI方法求解越來(lái)越困難，列車數(shù)為4列時(shí)，5.3 s即可得到最優(yōu)解，當(dāng)列車數(shù)增加到20列時(shí)，求解時(shí)間達(dá)到3 468 s.在規(guī)模較小時(shí)，IP_GUROBI比BAP_C#更具有優(yōu)勢(shì)，隨著規(guī)模增大，BAP_C#可以在較短的時(shí)間內(nèi)獲得滿意解，如列車數(shù)量為20列時(shí)，求解時(shí)間減少(3 463.7-290)/3 463.7=91.6%，得到的最終可行解距離最優(yōu)解的間隔為(508-463)/463=9.72%.

在分支策略中，一種較為常用的分支策略是最為分?jǐn)?shù)分支(Most Fractional Branching)策略，即每次選擇最接近0.5的變量進(jìn)行分支[6-7].現(xiàn)以列車數(shù)為20列的情況為例，比較本文分支策略和最為分?jǐn)?shù)分支策略，兩者的上下界收斂曲線如圖4所示.本文分支策略在290 s處收斂，最終收斂目標(biāo)值為508 min，而最為分?jǐn)?shù)分支策略在389 s處收斂，最終收斂目標(biāo)值為515 min，說(shuō)明了本文分支策略具有一定優(yōu)勢(shì).需要指出，最終收斂時(shí)上下界相等并不是得到了最優(yōu)解，因?yàn)榉种渲械哪承┙Y(jié)點(diǎn)沒(méi)有用列生成算法求解.

表4 IP_GUROBI和BAP_C#的求解時(shí)間和目標(biāo)值Table 4 The computational time and objective value for IP_GUROBI and BAP_C#

5.3 可變更線路與按原始線路行駛的比較

本文按原始線路行駛指的是限定列車不能使用交叉渡線，同樣用GUROBI求解得到該種調(diào)整方式的目標(biāo)值，結(jié)果如圖5所示.可以看出，可變線路相比于按原始行駛線路得到優(yōu)化，目標(biāo)值降低31.4%～41.9%，平均降低37.4%.

圖4 不同分支策略上下界收斂曲線Fig.4 The upper and lower bound for different branch strategies

圖5 2種列車運(yùn)行調(diào)整方式的目標(biāo)值Fig.5 The objective value for two kinds of train rescheduling

6 結(jié)論

本文分別構(gòu)建了基于列車到發(fā)時(shí)刻和基于列車時(shí)空路徑的鐵路網(wǎng)絡(luò)列車運(yùn)行調(diào)整模型，前者適用于商業(yè)軟件求解，可獲取最優(yōu)解.結(jié)合后者模型的特點(diǎn)，設(shè)計(jì)了分支定價(jià)算法求解，對(duì)比得出隨著求解規(guī)模的增大，本文算法求解時(shí)間變化較為平緩，且能夠得到較優(yōu)解，表現(xiàn)出更好地適應(yīng)性.同時(shí)得出本文分支策略比最為分?jǐn)?shù)分支策略可以更快地使算法收斂，可變更線路相對(duì)于按照原始線路行駛的調(diào)整方式可降低目標(biāo)值31.4%～41.9%.

[1]彭其淵，朱松年，王培.網(wǎng)絡(luò)列車運(yùn)行圖的數(shù)學(xué)模型及算法研究[J].鐵道學(xué)報(bào)，2001，23(1):1-8.[PENG Q Y，ZHU S N，WANG P.Study on a general optimization model and its solution for railway network train diagram[J].Journal of the China Railway Society，2001，23(1):1-8.]

[2]MENG L Y，ZHOU X S.Simultaneous train rerouting and rescheduling on an N-track network:a model reformulation with network-based cumulative flow variables[J].Transportation Research Part B Methodological，2014，67(3):208-234.

[3]陳雍君，周磊山.基于序優(yōu)化方法的列車運(yùn)行調(diào)整算法研究[J].鐵道學(xué)報(bào)，2010，32(3):1-8.[CHEN Y J，ZHOU L S.Study on train operation adjustment algorithm based on ordinal optimization[J].Journal of the China Railway Society，2010，32(3):1-8.]

[4]MU S，DESSOUKY M.Scheduling freight trains traveling on complex networks[J].Transportation Research Part B Methodological，2011，45(7):1103-1123.

[5]CACCHIANI V，CAPRARA A，TOTH P.A column generation approach to train timetabling on a corridor[J].4or Quarterly Journal of the Belgian French&Italian Operations Research Societies，2008，6(2):125-142.

[6]ACHTERBERG T，KOCH T，MARTIN A.Branching rules revisited[J].Operations Research Letters，2005，33(1):42-54.

[7]FRIBERG D.An implementation of the branch-and price algorithm applied to opportunistic maintenance planning[D].Gothenburg:University of Gothenburg，2015.