基于魯棒接續時間的鐵路空車調運優化

羅 霞 ，胡劍鵬 ，甘易玄

（1.西南交通大學交通運輸與物流學院, 四川 成都 611756；2.西南交通大學綜合交通運輸智能化國家地方聯合工程實驗室, 四川 成都 611756）

空車調運是實現鐵路資源優化配置的重要內容，空車調運的優劣直接影響鐵路車輛的利用效率、市場需求的滿足程度以及運輸組織的現代化水平.

近年來空車調運優化問題的研究多從供需及環境的隨機性出發建立動態隨機規劃模型.文獻[1]基于廣義旅速建立時空服務網絡，提出了空車調配多階段動態優化方法；在此基礎上，文獻[2]著重考慮了實際運輸生產中的能力約束，建立了考慮車種替代動態規劃模型；文獻[3]則將空車調配分為基于固定需求的優化調整和基于客戶實時需求的策略再優化兩個階段，解決了空車需求的動態變化問題；文獻[4]基于供需的不確定性，為約束條件和目標函數設置了一定的置信水平，建立了多車種調配優化模型；文獻[5]考慮了路網中轉能力和通行能力不確定性的影響，建立多目標隨機期望約束模型；文獻[6]則針對列車走行時間的不確定性，提出了魯棒優化模型，并利用對稱和對偶變化降低了模型求解難度，提高了求解效率.

在上述研究的基礎上，考慮到實際鐵路生產中空車的產生以及空車的需求是一個時空問題，除了物理層面的連接外，還應當結合空車與列車接續時間關系進行列車之間的空車配流.同時，網絡中列車的旅行時間以及車站的技術作業時間具有不確定性，時間的波動會直接對供需網絡的接續關系造成影響.因此，本文基于魯棒優化理論對多車種、可替代條件下的空車調運問題進行了研究.

1 問題描述與建模

1.1 問題描述

供應站與需求站的車站技術作業時間以及站間旅行時間不確定性時，如何確定需求站的空車來源、類型與空車數量是魯棒接續時間要求下空車調運的核心問題.模型中包含下述參數.

定義集合及索引：I、i分別為空車供應站集合和索引，i∈I；J、j分別為空車需求站集合及索引，j∈J；E、e分別為供應站內產生空車的到達列車集合及索引，e∈E；F、f分別為供應站發出的可掛運空車的列車集合及索引，f∈F；G、g分別為需求站內需使用到達空車裝車的發出列車集合及索引，g∈G；η為實際的空車類型， χ 為考慮車種替代后按用途歸類的空車類型，空車類型集合H包括平車（N）、棚車（P）和敞車（C），H={N,P,C} ，且 η , χ ∈H.

定義供應站相關變量：ci、ti分別為供應站i的等待車小時費用及車站技術作業時間；Sie,η、Sif分別為供應站i內第e列到達列車產生 η 型空車的數目及第f列發出列車的最大空車掛運數目；Aie、Lif分別為供應站i第e列到達列車的到達時刻及第f列發出列車的最晚編組時刻.

定義需求站相關變量：pj、cj和tj分別為需求站j的空車裝車后產生收益、空車等待車小時費用及車站技術作業時間；Sjg,χ為空車需求站j第g列發出列車需要的 χ 型空車數量；Ljg為空車需求站j第g列發出列車的最晚編組時刻.

定義0-1變量：xi,e,f表示供應站i的第e列到達列車與供應站i的第f列發出列車間是否存在接續關系，xi,e,f= 1 為是，xi,e,f= 0 為否；yif,jg= 0 表示供應站i的第f列發出列車與需求站j的第g列發出列車間是否存在接續關系，yif,jg= 1 為是，yif,jg= 0 為否.

定義決策變量：di,e,f,η為供應站i的第e列到達列車為供應站i發出的第f列列車供應的η型空車的數量；dif,jg,χ為供應站i第f列發出列車為需求站j發出的第g列列車供應的χ型空車的數量；dif,jg,η為供應站i第f列發出列車為需求站j第g列發出列車供應的η型空車的數量.

另外定義了空車供應站i和空車需求站j間的空車輸送費用cij及站間旅行時間tij.

1.2 名義模型構建

名義模型是將供應站與需求站車站技術作業時間以及站間旅行時間看作確定值時的空車調運模型.

1）目標函數

以空車調運效益最大化為目標函數，如式（1）所示，其中：等號右端第1項表示空車到達需求站后由需求站裝車發出所產生收益；第2項表示空車調運時產生的路徑消耗費用[7]；第3項表示空車在需求站的等待費用；第4項表示空車在供應站的等待費用.

2）約束條件

式（2）表示供應站i的第e列到達列車給該站各發出列車f提供的η型空車總數不得超過第e列到達列車產生η型空車總數.

式（3）表示供應站i的各次到達列車e為該站第f列發出列車供應的空車總數不得超過第f列發出列車的最大空車掛運數目.

式（4）表示dif,jg,η中 η 型箱的數量等于其實際替換的各車型dif,jg,χ的數量總和.設Rη為 η 型箱可替代的車型集合，參考文獻[8]，平車代替敞車有90%概率，棚車代替平車和敞車有75%的可能，敞車代替棚車有30%的可能，故有RN= {N,C}，RP= {N,P,C}，RC= {P,C}.

式（5）表示供應站i發出的第f列列車中掛運的各型空車均用于滿足需求站的空車需求[9-10].

式（6）表示各供應站發出的列車為需求站j發出的第g列列車提供的χ型空車數量等于該列車對χ型空車的需求量.

式（7）表示列車間空車輸送數量非負.

式（8）表示如果供應站i的第e列到達列車在供應站i完成相關技術作業后，到達調車場的時刻不晚于i站第f列發出列車的最晚編組時刻，則該組供應站到達列車與發出列車之間滿足最小接續時間要求，即xi,e,f= 1，否則，xi,e,f= 0.供應站到達列車產生的空車來自重車卸車時，ti= (ti,dj+ti,zx) + (ti,jt+ti,zx) +(ti,xc+ti,zx)；到達列車本身掛有空車時，ti= (ti,dj+ti,zx) +(ti,jt+ti,zx)；其中：ti,dj、ti,zx、ti,jt和ti,xc分別為i站的到達技術作業時間標準、轉線作業時間標準、解體作業時間標準（含推峰時間）及卸車作業時間標準.

式（9）表示若供應站i發出的第f列列車在需求站j完成相關技術作業后，到達調車場的時刻不晚于j站發出第g列列車的最晚編組時刻，即yie,jg= 1，否則，yie,jg= 0.列車在需求站j的技術作業時間tj=(tj,dj+tj,zx) + (tj,jt+tj,zx) + (tj,zc+tj,zx)，其中：tj,dj、tj,zx、tj,jt、tj,zc分別為j站到達技術作業時間標準、轉線作業時間標準、解體作業時間標準及裝車作業的時間標準.

式（10）表示只有當供應站i第f列發出列車存在滿足接續時間要求的空車來源時，即xi,e,f= 1時，才考慮該發出列車與需求站發出列車的接續關系.

為便于說明，本文定義了屬于關系、連接關系和車流關系，如圖1所示，圖中：三角框為供應站，正方形框為需求站，圓圈為列車，后同.

圖1 空車調運網絡Fig.1 Empty wagon allocation network

定義1屬于關系指供應站到達列車e、發出列車f與供應站i之間的關系以及需求站發出列車g與需求站j之間的關系.

定義2連接關系為供應站發出的列車f與該列車到達的需求站j之間的關系，供應站發出的任一可掛運空車的列車有且只與一個需求站建立連接關系.

定義3車流關系是指當某一列車中的車輛可以成為另一列車中的車流來源時，兩列列車之間存在的關系.車流關系可以定義為兩部分：

1）供應站到達列車e與該站發出列車f之間的車流關系定義為第1階段車流關系；

2）供應站發出的列車f與需求站發出的列車e之間的車流關系定義為第2階段車流關系.

車流關系由約束條件式（8）～（10）確定.

1.3 魯棒優化模型構建

實際運輸中由于設備因素、環境因素和人為因素的影響，列車站間旅行時間以及車站技術作業時間是不確定的，若以固定的旅行時間和技術作業時間作為空車接續時間的判斷依據，當接續條件Lif-(Aie+ti)≥0或Ljg-(Lif+tij+tj)≥0 ，但取值較小時，技術作業時間和旅行時間的波動可能使得既有的接續關系失效，導致求得方案不可行.因此，基于Dimitris等[11]提出的魯棒優化理論建立了空車調運魯棒優化模型.

1）目標函數

魯棒優化模型的目標函數如式（12），約束關系如式（13）所示，結合波動程度與波動下限約束，式（12）中的C2如（14）所示，矩陣表達如式（15）所示.

式（13）、（15）中： δij、 δi、 δj、 θij、 θi、 θj為松弛系數；n為供應站數目 ；m為需求站數目；En為n×n的單 位 矩 陣 ；Zn為 元 素 全 為 1的 1 ×n矩 陣 ； δn=(δ1,δ2,···,δn)T； θn=(θ1,θ2,···,θn)T；Unm=(α11,α12,···,αnm)T；Vn=(β1,β2,···,βn)T；Wm=(γ1,γ2,···,γm)T.

證明一個矩陣為完全幺模矩陣的充分條件為：① 矩陣中元素僅包含 0、1、-1；② 矩陣每列最多兩個非0元素；③ 矩陣的行可分劃成兩個子集，使得同列中兩個非零元素符號相同時，對應的兩行在不同的行子集中，當符號不同時，對應的兩行在同一行子集中.易得：約束矩陣式（15）滿足條件 ①、②，并且當矩陣行劃分為(1,4)和(2,3,5,6)兩個行子集時滿足條件 ③，故該約束矩陣為完全幺模矩陣.結合完全幺模矩陣的性質，當 Γ1、 Γ2和 Γ3均取整數時，式（13）、（14）的最優解 α*、 β*、 λ*也為整數，因αij,βi,γj≥0，且 αij,βi, γj≤1 ，故其最優值取0或1.

2）約束條件

原約束條件中式（2）～（7）、（10）不變，式（8）、（9）分別變更為式（16）、（17）.

2 求解算法

2.1 名義模型求解算法

名義模型為非線性整數規劃模型，求解困難，因此，求解前結合約束關系式(8) ～ (10)對模型進行簡化.

以算例中的供應站3為例，為便于表示，此處僅保留三列供應站發出列車，變化前后網絡結構如圖2所示.與原網絡相比，新網絡中僅保留了xi,e,f=1 和yie,jg=1時所對應的決策變量di,e,f,η和dif,jg,η，進而去除模型中的約束條件式（8）～（10），轉化為線性規劃模型，可調用CPLEX進行求解.

圖2 簡化前后網絡圖Fig.2 Network diagrams of pre- and post-simplification

2.2 魯棒優化模型求解算法

魯棒優化模型中引入了變量 αij、 βi、 γj，這些變量會導致供需網絡中既有車流關系的變化，但對于每組 αij、 βi、 γj存在唯一的網絡結構，求解時可以看作名義模型.

遍歷滿足約束的每組 αij、 βi、 γj是最準確的求解方法，但當 Γ1、 Γ2、 Γ3取值較高時模型的求解復雜度會劇增，導致求解困難.研究發現，只有當網絡結構中的車流關系減少時，最優調運方案才可能改變，故設計下述步驟求解.求解時 Γ1、 Γ2、 Γ3和為已知值.

步驟1結合式(8)～(10)，得出固定出行時間下供需網絡中存在的車流關系，令供應站i的第1階段車流關系集合為 Ω1,i，供應站i的第2階段車流關系集合為 Ω2,i，形式如下：

式中：Ω1,i中(e,f)表示i站第e列到達列車可以為i站第f列發出列車供應空車；Ω2,i中(f,j,g)表示i站發出的第f列列車可以為j站發出的第g列列車供應空車.

步驟 2令集合U、V、W中的各項元素均為 1，即所有的tij、ti、tj都存在波動，結合式（10）、（16）、（17），得出絕對魯棒下的網絡中存在的車流關系 Π1,i、Π2,i，形式同 Ω 1,i、 Ω 2,i.

步驟 3比較 Π 1,i、 Π 2,i和 Ω 1,i、 Ω 2,i中的車流關系數量，確定此時減少的車流關系.若 Π1,iΩ1,i，說明第1階段車流關系發生了變化，若 Π2,iΩ2,i，說明第2階段車流關系發生變化.假設此時減少的車流關系為 Ω1,1的 (1,1)， Ω1,2的 (1,2)， Ω2,2的(1,1,1)，Ω2,3的(2,2,1).

步驟4找出減少的車流關系中包含的供應站集合I1、需求站集合J1以及連接關系集合O，以步驟3中假設為例，減少的車流關系中包含的供應站為I1= {1,2,3}，包含的需求站為J1= {1,2}，包含的連接關系為O={(2,1),(3,2)}.

步驟5令集合U、V、W中的各項元素均為0.同時，根據供應站n1、需求站n2、連接關系n3及波動下限 Γ1、Γ2、Γ3，按下述方法確定 βi、 γj、 αij的取值：

1）若 Γ1≥|n1| （ |n1| 表示n1中元素個數，其余同），則令I1內各供應站的系數 βi為 1，以I1={1,2,3} 為例

即 β1=1 、 β2=1 、 β3=1 ；若 Γ1＜|I1| 則得到從I1中取Γ1個元素的全排列集合N1，令每組排列的系數為1；

2）針對 Γ2的調整方法同T1，當 Γ2＜|J1| 時，得到從J1中取 Γ2個元素的全排列集合N2；

3）若 Γ3≥|O| ，則令O中各元素的系數 αij=1 ，以步驟3中假設為例， α21=1 ， α32=1 ；否則，得到從O中取 Γ3個元素的全排列組合N3,令每組排列的系數為1.

步驟 6輸入由N1、N2、N3中元素組成的 αij、 βi、γj排列，調用CPLEX求解各網絡結構下的最佳效益值，取其中的最小值作為魯棒解.

該算法在給定波動量t?ij、t?i、t?j的情況下，找出最劣網絡結構較理想網絡結構中減少的車流關系，以此為入手點找出波動下限 Γ1、 Γ2、 Γ3下最可能引起網絡結構變化的 αij、 βi、 γj排列，避免了無效遍歷，減少了求解時間.

3 模擬計算及結果分析

3.1 參數輸入

以5個需求站和4個供應站構成的網絡進行分析.供應站列車的到達時刻Aie、產生的各型空車數量Sie,η見表1.供應站發出列車f的最晚編組時刻Lif、發往需求站j和最大空車掛運數量Sif見表2.需求站發出列車g的最晚編組時刻Ljg及空車需求量Sjg,χ見表3.站間空車輸送費用cij、名義旅行時間tij見表4.各車站等待車小時費用ci、cj，固定作業時間ti、tj，空車裝車后的發送效益pj見表5.規劃階段起始時刻為 0，1 天為 1 440 min，表中負號對應前一天的時間.前一階段殘存空車也按照其實際列車到達情況加入供應站到達列車數據.

表1 供應站列車到達時刻及空車產生數量Tab.1 Time of train arrival and number of empty wagons produced in supply stations

表2 供應站列車最晚編組時刻及最大空車掛運數量Tab.2 Train marshalling time limit and maximum number of empty wagons in supply stations

表3 需求站列車最晚編組時刻及所需空車數量Tab.3 Train marshalling time limit and number of empty trains required in demand stations

表4 站間空車輸送費用及名義旅行時間Tab.4 Empty wagon transportation cost and nominal travel time between stations

表5 車站參數數據Tab.5 Station parameter data

3.2 名義模型求解結果

通過MATLAB 2018a調用數學優化軟件CPLEX編程求解，得到該供需網絡下的最優空車匹配方案，供應站各發出列車掛運的空車來源見表6.各需求站發出列車的空車來源見表7.該網絡下存在474個第1階段車流關系，204個第2階段車流關系，空車輸送效益值為53 708元.

表6 供應站發出列車的空車來源Tab.6 Empty wagon source of trains leaving from supply stations

表7 需求站發出列車空車來源Tab.7 Empty wagon source of trains leaving from demand station

空車替代情況：需求站2發出的列車2中由列車(4,1,1)供應的1輛敞車為平車用作敞車.需求站4發出的列車2中由列車(1,5,4)供應的棚車中，有1輛為敞車用作棚車.需求站5發出的列車1中由列車(1,3,7)供應的棚車中，2輛為敞車用作棚車.

3.3 魯棒優化模型分析

波動量和波動下限造成的網絡結構變化可能導致無可行解，為保證可行解的存在對網絡進行下述變化.

1）為各供應站增加一輛虛擬到達列車，該列車與供應站所有發出列車f均存在接續關系，且掛有充足的各類型空車；

2）任一供應站與所有需求站間均增加一列虛擬發出列車，該列車與需求站所有發出列車g均存在接續關系且發出列車的空車掛運能力充足；

3）當供應站發出列車空車來源為供應站虛擬到達列車時，供應站庫存費用取較大常數，當需求站發出列車車流來源為供應站虛擬發出列車時，需求站發車效益為0.

4）求得最優解后，除去供應站虛擬到達列車和虛擬發出列車對應車流關系，結合式（1）求解實際效益，作為當前方案下的效益值.

3.3.1 單因素波動影響分析

模型中的不確定因素包括站間旅行時間、供應站及需求站作業時間3項.各因素的不確定程度由波動量和波動下限 Γ1、 Γ2、 Γ3共同決定，令不同波動下限下效益值隨的變化即可看作效益值隨波動率q1、q2、q3的變化.單因素波動導致效益值和網絡結構的變化如圖3所示.

圖3 單因素波動率影響分析Fig.3 Fluctuation influence analysis from single factor

由圖3可知：

1）單因素的不確定性導致效益值出現較大變化時，必然伴隨車流關系的變化.只有當減少的車流關系屬于當前最優調運方案時效益值才會有較大變化.

2）相同波動率下效益值差距較大時，波動下限越大效益值越小.此時，相同波動率下的最優調運方案不同，波動下限越大，網絡中存在車流關系越少，網絡結構越劣，效益越小.

3）相同波動率下效益值很接近時，波動下限越大效益值越大.此時最優調運方案一致，波動下限越大，受到影響的作業或旅行時間越多，車輛的等待時間越小，效益值越大.

4）從單因素的絕對魯棒性來看，波動率相同時，效益值越小則受波動率影響越大，故各因素對效益的影響程度為tj＞tij＞ti.

3.3.2 多因素波動率影響分析

給定 Γ1=15 ， Γ2=2 ， Γ3=3 ，q1 、q2 、q3 之間互動關系如圖4.q2的大小決定第1階段連接關系的數量，也決定了優化結果的上限.隨著q2的增大，效益值整體波動幅度逐漸減小，q2不變時效益值隨q1、q3的增大而減小.

圖4 多因素波動率影響分析Fig.4 Fluctuation influence analysis from multiple factors

3.3.3 波動下限影響分析

當q1=5% ，q2=4% ，q3=5% 時，不同Γ1、Γ2、Γ3下方案效益值如圖5所示.

圖5 波動下限影響分析Fig.5 Influence analysis of lower limit of fluctuation

1）站間旅行時間波動下限影響分析

Γ2、 Γ3取值相同， Γ1∈[0,2]時， Γ1對效益有顯著影響，并且 Γ1越大效益越小.Γ1∈[3,20]時， Γ1對效益無顯著影響， Γ1越大效益值越大.

2）供應站作業時間波動下限影響分析

Γ1、 Γ3取值相同， Γ2的取值對效益無較大影響，并且 Γ2越大效益值越大.

3）需求站作業時間波動下限影響分析

Γ1、Γ2取值相同，Γ3∈[0,1]時的效益明顯大于 Γ3∈[2,5]時的效益，且效益與 Γ3的呈負相關.但當Γ3∈[2,5]時，雖然效益仍與 Γ3呈負相關，但無顯著差異.

綜上，供應站作業時間波動下限對效益值的影響最小.旅行時間波動下限和需求站作業時間波動下限取值較小時，它們對效益值有較大影響，并且二者間存在相互影響.

4 結 論

1）提出供應站到達列車、發出列車以及需求站發出列車間的空車接續時間關系判斷方法，并據此分別建立了確定環境與不確定環境下的空車調運優化模型.

2）考慮構建模型的特征，將確定環境下的模型轉化為了線性模型，并結合車流關系的變化是效益值變化的主要原因這一特點，設計了魯棒優化模型的求解算法，求得了考慮車種替代的空車配流方案.

3）需求站技術作業時間波動率對效益值影響最大.供應站技術作業時間波動下限對效益值影響最小.方案效益值隨方案魯棒性的提高而下降.