李婷婷

(北京大學(xué) 城市規(guī)劃與設(shè)計(jì)學(xué)院，廣東 深圳518055)

軸輻網(wǎng)絡(luò)(hub-and-spoke network)廣泛應(yīng)用于通訊和交通領(lǐng)域[1]。其中，樞紐選址問題[2](hub location problem，HLP)經(jīng)常用于設(shè)計(jì)與分析軸輻網(wǎng)絡(luò)。研究表明[3]，軸輻網(wǎng)絡(luò)的樞紐很重要。然而，交通運(yùn)輸網(wǎng)絡(luò)的樞紐因自然災(zāi)害、惡劣天氣、設(shè)備故障、蓄意攻擊等可能被阻斷(interdiction，即禁止或停止運(yùn)行)，從而產(chǎn)生嚴(yán)重不良影響。因此，有必要識(shí)別軸輻網(wǎng)絡(luò)的關(guān)鍵樞紐以使阻斷發(fā)生后損失最小化。

軸輻網(wǎng)絡(luò)關(guān)鍵樞紐識(shí)別問題實(shí)質(zhì)為樞紐阻斷問題[4](hub interdiction problem，HIP)，目的在于識(shí)別軸輻網(wǎng)絡(luò)中哪些樞紐被阻斷后造成最大損失，即識(shí)別關(guān)鍵樞紐。與類似的網(wǎng)絡(luò)阻斷問題[5](network interdiction problem，NIP)、設(shè)施阻斷問題[6-7](facility interdiction problem，F(xiàn)IP)相比，HIP的研究較少。Lei[8]提出了樞紐阻斷中位問題(r-hub interdiction median problem，r-HIMP)，即從現(xiàn)有軸輻網(wǎng)絡(luò)中識(shí)別r個(gè)關(guān)鍵樞紐，使得r個(gè)樞紐被阻斷后軸輻網(wǎng)絡(luò)中的最小總需求加權(quán)路徑成本最大，并構(gòu)建了雙層規(guī)劃模型。Parvaresh等[9]研究了考慮惡意破壞的樞紐網(wǎng)絡(luò)設(shè)計(jì)問題，構(gòu)建了雙目標(biāo)的雙層規(guī)劃模型(下層模型是HIP)，提出了基于模擬退火與禁忌搜索的啟發(fā)式算法。Ghaffarinasab等[4]提出了3種類型的HIP，包括r-HIMP、樞紐阻斷最大覆蓋問題(r-hub interdiction maximal covering problem，r-HIMCP)、樞紐阻斷中心問題(r-hub interdiction center problem，r-HICP)，分別構(gòu)建雙層規(guī)劃模型并設(shè)計(jì)模擬退火算法求解。Ghaffarinasab等[10]對(duì)r-HIMP構(gòu)建了雙層規(guī)劃模型，提出了隱枚舉算法求解。Ramamoorthy等[11]對(duì)r-HIMP構(gòu)建了雙層規(guī)劃模型，比較了基于對(duì)偶與最近分配約束的將雙層規(guī)劃模型轉(zhuǎn)化為單層模型的方法，設(shè)計(jì)了Benders分解算法求解。

以上研究均采用雙層規(guī)劃模型解決HIP，然而均未考慮樞紐的能力限制。本文研究有能力限制的r-HIMP(capacitatedr-HIMP，Cr-HIMP)，考慮由于能力限制，樞紐被阻斷后部分需求可能得不到滿足的情況，構(gòu)建雙層規(guī)劃模型，通過算例驗(yàn)證方法的有效性，以期為有能力限制的軸輻網(wǎng)絡(luò)關(guān)鍵樞紐識(shí)別問題提供決策參考。

1模型

1.1 問題描述

Cr-HIMP考慮樞紐的能力限制，從現(xiàn)有軸輻網(wǎng)絡(luò)的p個(gè)樞紐中識(shí)別r個(gè)關(guān)鍵樞紐。當(dāng)r個(gè)關(guān)鍵樞紐被阻斷后，使得網(wǎng)絡(luò)的最小總需求加權(quán)路徑成本最大，部分需求由于能力限制不能得到滿足時(shí)則需要相應(yīng)懲罰。Cr-HIMP涉及襲擊者和防衛(wèi)者兩方，屬于Stackleberg博弈，襲擊者作為領(lǐng)導(dǎo)者，防衛(wèi)者作為追隨者，可以通過雙層規(guī)劃模型描述。其中，上層模型表示襲擊者從現(xiàn)有軸輻網(wǎng)絡(luò)的p個(gè)樞紐中挑選r個(gè)關(guān)鍵樞紐進(jìn)行襲擊，目標(biāo)為最大化阻斷發(fā)生后防衛(wèi)者的最小總需求加權(quán)路徑成本，下層表示防衛(wèi)者在阻斷發(fā)生后p-r個(gè)樞紐仍然正常運(yùn)行的情況下選擇最佳路線，目標(biāo)為最小化阻斷發(fā)生后的總需求加權(quán)路徑成本。

本文假設(shè)樞紐被襲擊后完全阻斷而不能運(yùn)行，同時(shí)遵循HLP的一般假設(shè)[12]。

1.2 模型構(gòu)建

N為所有節(jié)點(diǎn)的集合，H為現(xiàn)有樞紐的集合。參數(shù)說明如下。p為現(xiàn)有樞紐的數(shù)量；r為被阻斷的樞紐數(shù)量(即關(guān)鍵樞紐的數(shù)量)；wij為節(jié)點(diǎn)i、j之間的需求；dijkm為從節(jié)點(diǎn)i經(jīng)過樞紐k和m到達(dá)節(jié)點(diǎn)j的成本，其中dijkm=αcik+βckm+γcmj，α、β、γ分別為集合線路、樞紐線路、疏運(yùn)線路的折扣系數(shù)，cik、ckm、cmj分別為集合線路、樞紐線路、疏運(yùn)線路的運(yùn)輸成本；θ為單位需求沒有得到滿足的懲罰成本；Ck為樞紐k的能力。變量說明如下。zk=1表示樞紐k沒有被阻斷，否則為0；xijkm表示阻斷發(fā)生后從節(jié)點(diǎn)i經(jīng)過樞紐k和m到達(dá)節(jié)點(diǎn)j的流量占i、j之間需求的比例；T為下層模型的目標(biāo)函數(shù)值。

Cr-HIMP模型由上層模型(upper level problem,ULP)和下層模型(lower level problem,LLP)組成。

式(1)表示上層目標(biāo)為最大化阻斷發(fā)生后的最小總需求加權(quán)路徑成本；式(2)表示r個(gè)樞紐被阻斷；式(4)表示下層目標(biāo)為最小化阻斷發(fā)生后的總需求加權(quán)路徑成本；式(5)表示節(jié)點(diǎn)i、j間的需求可能不全部滿足；式(6)表示樞紐沒有被阻斷才能有流量通過；式(7)表示通過樞紐的流量不能超過其能力；式(3)、式(8)為變量約束。

1.3模型求解

式(4)等價(jià)于

式(5)～(7)分別引入對(duì)偶變量λij、ρijk、ηk，LLP的對(duì)偶問題如下。

LLP-DF(dual formulation of lower level problem)：

令φij=-λij、δijk=-ρijk、πk=-ηk，則LLP-DF等價(jià)于LLP-NDF(new dual formulation of lower level problem)。

將LLP-NDF與ULP合并，得到max-max問題從而把雙層規(guī)劃模型轉(zhuǎn)化為單層模型Cr-HIMP-S(singlelevel formulation of capacitatedr-HIMP)。

參考文獻(xiàn)[13-14]將Cr-HIMP-S線性化，用變量

綜上所述，Cr-HIMP轉(zhuǎn)化為單層的線性模型Cr-HIMP-SL(single-level and linearized formulation of capacitatedr-HIMP)。

Cr-HIMP-SL：

s.t.式(2)，(3)，(15)～(17)，(19)～(24)。

由于ρijk為約束(6)的影子價(jià)格，且 ρi jk≤0、δijk=-ρi jk，則δijk表示約束(6)右側(cè)減少1單位后目標(biāo)函數(shù)值變化的量。具體地，對(duì)于給定的i、j、k，當(dāng)zk為1，i、j可以經(jīng)過或者不經(jīng)過樞紐k，若i、j間的需求都分配到最長(zhǎng)路徑i-k1-m1-j，對(duì)應(yīng)的最大路徑成本為此情況下目標(biāo)函數(shù)值最大；當(dāng)zk變?yōu)?時(shí)不能經(jīng)過樞紐k，若i、j間的需求都分配到最短路徑i-k2-m2-j(k2≠k，m2≠k)，對(duì)應(yīng)的最小路徑成本為此情況下目標(biāo)函數(shù)值最小。因此，目標(biāo)函數(shù)值的最大變化量為故

由于ηk為約束(7)的影子價(jià)格，且ηk≤0、πk=-ηk，則πk表示約束(7)右側(cè)減少1單位后目標(biāo)函數(shù)值變化的量。具體地，對(duì)于給定的k，當(dāng)zk=1，由于增加1單位能力最多服務(wù)1單位需求：1)Ck減少1后，有1單位需求沒有滿足(被懲罰)，Ck減少1前，該1單位需求分配到k且經(jīng)過最短路徑i1-k-m3-j1(對(duì)應(yīng)的最

短數(shù)

路值

徑變

成化

本為

為

θ-

ddii11jj11kk mm33

；

=

mi2,

j,i

)mnC{

dk減i

jkm少})，1后

此，情有況1

下單

目位

標(biāo)需

函求不能分配到k且經(jīng)過最長(zhǎng)路徑i2-k4-m4-j2，對(duì)應(yīng)的最大路徑成本為di2j2k4m4=i,mj,l≠akx,m{di jlm}，Ck減少1前，該1單位需求分配到k且經(jīng)過最短路徑i1-k-m3-j1(對(duì)應(yīng)的最短路徑成本為di1j1km3=mi,j,imn{di jk m})，此情況下目標(biāo)函數(shù)值變化為di2j2k4m4-di1j1km3。因此，目標(biāo)函數(shù)值的最大變化為max{ θ-di1j1km3，di2j2k4m4-di1j1km3}，故πk≤max{θ-di1j1k m3,di2j2k4m4-di1j1km3}，可取Mk1=max{θ-di1j1k m3,di2j2k4m4-di1j1k m3}為πk的有效上界。

綜上所述，Cr-HIMP等價(jià)于Cr-HIMP-SL，后者為線性規(guī)劃模型，可以使用商業(yè)優(yōu)化軟件求解。

2算例

以CAB數(shù)據(jù)集(the civil aeronautics board(CAB)data set)[15]為例，該數(shù)據(jù)包含美國(guó)25個(gè)城市(如表1所示)之間的航空客流需求與成本。另外，p個(gè)樞紐及其能力采用文獻(xiàn)[16]的方法生成。為作比較，以下參數(shù)取與文獻(xiàn)[11]相同的值：p=7，α=γ=1，β ∈{0.9,0.5,0.1}，r∈{1,2,3,4,5}。θ參照文獻(xiàn)[17]分別取2 000、3 000、 4 000進(jìn)行計(jì)算。在主頻為3.20 GHz的Intel Core i5-6500處理器、內(nèi)存8 GB、操作系統(tǒng)Windows 10環(huán)境下，使用優(yōu)化軟件CPLEX 12.2[18]對(duì)模型求解，結(jié)果如表2所示。

表1 CAB數(shù)據(jù)集中25個(gè)城市的編號(hào)[19]Table 1 The ID number of 25 cities in the CAB data set[19]

表2模型r-HMIPDF[11]與C r-HMIP-SL的求解結(jié)果比較Table 2 Comparison of the results of models r-HMIPDF[11]with C r-HMIP-SL

對(duì)于25個(gè)城市、7個(gè)樞紐的軸輻網(wǎng)絡(luò)：文獻(xiàn)[11]不考慮樞紐能力限制，模型r-HMIPDF[11]的約束數(shù)量為 4 390，變量數(shù)量為 5 015；本文考慮樞紐能力限制，模型Cr-HMIP-SL的約束數(shù)量為30 668，變量數(shù)量為 5 029?？梢?，考慮樞紐能力限制使約束數(shù)量和變量數(shù)增多。如表2所示，Cr-HMIP-SL計(jì)算時(shí)間普遍比r-HMIPDF[11]長(zhǎng)，r-HMIPDF[11]計(jì)算時(shí)間均在15 s以下，Cr-HMIP-SL計(jì)算時(shí)間最長(zhǎng)超過1h(當(dāng)β=0.1、r=1、θ= 4 000時(shí))。求解模型r-HMIPDF[11]與Cr-HMIP-SL得到的關(guān)鍵樞紐編號(hào)在某些情況下不相同，說明考慮能力限制與不考慮能力限制的關(guān)鍵樞紐識(shí)別結(jié)果會(huì)有差異。例如，當(dāng)β=0.9、r=4，r-HMIPDF[11]對(duì)應(yīng)關(guān)鍵樞紐編號(hào)為“4”“7”“17”“20”，而Cr-HMIP-SL對(duì)應(yīng)“4”“14”“17”“20”。

對(duì)模型Cr-HMIP-SL，從表2看出，隨著β減小或r減小或θ增大，計(jì)算時(shí)間普遍增加；當(dāng)θ取不同值時(shí)，關(guān)鍵樞紐可能不相同，說明θ影響關(guān)鍵樞紐識(shí)別結(jié)果。例如，當(dāng)β=0.9、r=3、θ= 2 000，關(guān)鍵樞紐編號(hào)為“4”“17”“20”，而β=0.9、r=3、θ=3 000或4 000，關(guān)鍵樞紐編號(hào)為“4”“14”“17”。然而，某些關(guān)鍵樞紐比較固定，在管理中可對(duì)這些樞紐加強(qiáng)設(shè)防。例如，當(dāng)r=2，即使β、θ變化，“17”均為關(guān)鍵樞紐；當(dāng)r=3或r=4，即使β、θ變化，“4”“17”均為關(guān)鍵樞紐；當(dāng)r=5，即使β、θ變化，“4”“14”“17”均為關(guān)鍵樞紐。根據(jù)以上結(jié)果，判別出“14”“4”“17”關(guān)鍵程度依次增加，應(yīng)實(shí)行越來越高級(jí)別的設(shè)防。

3結(jié)論

1)軸輻網(wǎng)絡(luò)關(guān)鍵樞紐識(shí)別問題實(shí)質(zhì)為HIP。在現(xiàn)有模型r-HIMP的基礎(chǔ)上，增加樞紐能力約束，使用懲罰成本表示部分需求可能在樞紐阻斷發(fā)生后不能滿足，構(gòu)建了雙層規(guī)劃模型，通過下層模型求對(duì)偶問題將雙層規(guī)劃模型轉(zhuǎn)化為單層規(guī)劃并線性化，以CAB數(shù)據(jù)集為例，通過CPLEX優(yōu)化軟件求解，驗(yàn)證了模型的有效性。

2)在不考慮樞紐能力限制模型基礎(chǔ)上，增加樞紐能力約束使模型規(guī)模增大、求解時(shí)間變長(zhǎng)，且考慮能力限制與不考慮能力限制模型的關(guān)鍵樞紐識(shí)別結(jié)果有差異。有能力限制的模型計(jì)算時(shí)間普遍隨著β或r減小或θ增大而增加，且關(guān)鍵樞紐識(shí)別結(jié)果受θ影響。通過比較不同參數(shù)下的關(guān)鍵樞紐，可識(shí)別相對(duì)固定的關(guān)鍵樞紐并在管理中加強(qiáng)設(shè)防。例如，CAB數(shù)據(jù)集中，“14”“4”“17”為相對(duì)固定的關(guān)鍵樞紐，關(guān)鍵程度依次增加，應(yīng)實(shí)行越來越高級(jí)別的設(shè)防。