劉天虎，許維勝，吳啟迪

（同濟大學(xué) 電子與信息工程學(xué)院，上海201804）

車輛路徑問題（vehicle routing problem，VRP）是在滿足一定約束條件下使適宜的行車路徑達到既定的目標（如費用最小、耗時最短、行程最短等）.而帶時間窗的車輛路徑問題（vehicle routing problem with time windows，VRPTW）則是在VRP問題上加上了訪問需求的時間窗口.Desaulniers等[1]研究了具有等待成本的多資源站的VRPTW問題，得到了最優(yōu)解.Qureshi等[2]提出了在軟時間窗約束下的VRP搜索的精確最優(yōu)解決方案，從而識別最短的行駛路徑，控制最小的成本消耗.Ghoseiri等[3]對多目標的VRPTW搜索問題進行了系統(tǒng)研究，利用目標規(guī)劃及遺傳算法得到了最優(yōu)路徑解.為了求解VRPTW問題，過去的研究中最常用的算法有禁忌搜索算法[4]、插入啟發(fā)式算法[5]、粒子群算法[6]、人工智能算法[7]、拉格朗日松弛算法、近似算法等，本文在以上研究的基礎(chǔ)上，利用混合遺傳算法（hybrid genetic algorithm，HGA）[8]研究當醫(yī)療救援站僅有1個時，如何調(diào)度多輛醫(yī)療救援車輛的最優(yōu)VRPTW問題.

1 問題描述

突發(fā)災(zāi)害事件的發(fā)生往往對醫(yī)療救援活動提出急迫的需求，如何在有限的資源條件下實現(xiàn)最有效的救援問題越來越引起人們的重視，特別是如發(fā)生地震這樣的災(zāi)害事件，要求醫(yī)療急救資源迅速做出響應(yīng)，這就需要各個醫(yī)療救援站對救援車輛進行合理調(diào)度，然而一個醫(yī)療救援站不可能為所有的災(zāi)害點提供救援服務(wù)，而且每個災(zāi)害點在不同階段的醫(yī)療救助需求也可能不同.本文需要解決以下2個問題：①當發(fā)生突發(fā)災(zāi)害事件時如何最小化救援的總成本.總成本的控制對救援站資源的使用效率將帶來重要影響.②如何評價HGA算法在解決帶軟時間窗VRPTW搜索問題方面的有效性.需要深入分析尋優(yōu)模型，并通過與精確算法數(shù)例仿真比較分析來說明其有效性.

2 模型的建立

實際救援中需要依據(jù)災(zāi)害點的特征及時調(diào)整救援路徑，利用軟時間窗作為約束條件比硬時間窗更加具有彈性，可以得到更加可行的解決方案.當災(zāi)害點不能在軟時間窗的時限內(nèi)得到及時救助時，引入懲罰成本，懲罰成本參數(shù)可以依據(jù)實際情況加以調(diào)整，通過權(quán)衡懲罰成本和使用額外車輛的固定成本，可以對醫(yī)療救援車輛進行合理調(diào)度.假定條件：①假定僅有1個醫(yī)療救援站為所屬災(zāi)害點提供救援服務(wù)，救援車輛完成服務(wù)后需返回醫(yī)療救援站.②假定醫(yī)療救援站的車輛數(shù)及所有車輛的運載能力已知.③假定救援車輛與醫(yī)療救援站之間可以實現(xiàn)無障礙溝通，醫(yī)療救援站可監(jiān)控所有救援車輛的所在位置，可以評估出救援車輛的行駛時間信息.

這里可以通過圖1來說明軟時間窗下的多車路徑搜索問題.圖中的第1個括號內(nèi)的數(shù)值表示災(zāi)害點的x，y坐標方位，而第2個括號內(nèi)的數(shù)值表示災(zāi)害點對醫(yī)療救援的需求量，而單車的最大承載量不能超過13個單位量；A0表示救援站；Di表示災(zāi)害點，i＝1，…，13.

從圖1a中可以看出，最初的災(zāi)害點有10個，有2部醫(yī)療救援車輛提供救援服務(wù)，車輛1的行駛路徑為：A0→D5→D1→D9→D8→D6→A0.而車輛2的行駛路徑為：A0→D4→D10→D3→D7→D2→A0.

當T1時刻車輛1行駛到災(zāi)害點D5，而車輛2行駛到D10時，此時收到醫(yī)療救援站的最新信息，又有4個災(zāi)害點D11，D12，D13，D14需要救助，從圖1b可以看出，在規(guī)劃的路徑范圍內(nèi)，災(zāi)害點D11和D13可由車輛1提供救助，而災(zāi)害點D12和D14可由車輛2提供救助，在不超過車載容量的前提下，對2部醫(yī)療救援車輛的后續(xù)行駛路徑進行了優(yōu)化，車輛1后續(xù)行駛路徑優(yōu)化為：D5→D13→D1→A0，A0→D11→D9→D8→D6→A0，車輛2后續(xù)行駛路徑優(yōu)化為：D10→D3→D14→A0，A0→D7→D12→D2→A0.

從以上的分析可以看出，軟時間窗下的路徑搜索問題實際是一個混合整數(shù)線性規(guī)劃，于是可以建立相應(yīng)的數(shù)學(xué)模型，模型的目標是為了最小化救援的總成本（包括車輛成本、路徑成本和懲罰成本）.而該模型的限制條件由車輛、路徑、災(zāi)害點時間窗所約束，于是可以建立模型的最優(yōu)方程z如下：

限制條件如下：

式中：Cf為車輛使用的固定成本；i，j為災(zāi)害點；M為所有車輛的出發(fā)節(jié)點；Ψ為所有災(zāi)害及救援節(jié)點；τ，t0，tn為時刻點；θij，τ為從τ 時刻開始救援車輛從災(zāi)害 點i到 達 災(zāi) 害 點j， θij，τ＝；G為所有車輛的到達節(jié)點；Tij，τ為從τ時刻開始從災(zāi)害點i到達災(zāi)害點j的時間；We為過早到達災(zāi)害點的懲罰成本；Ei為在災(zāi)害點i的等待時間；Wl為延遲到達災(zāi)害點的懲罰成本；Li為延遲到達災(zāi)害點i的時間；M0為車輛的出發(fā)節(jié)點0；M1為車輛的出發(fā)節(jié)點1；G0為車輛的到達節(jié)點0；G1為車輛的到達節(jié)點1；Q為使用車輛的集合；αi，k為 災(zāi) 害 點i由 車 輛k實 施 救 援，αi，k＝，k為車輛；d為災(zāi)j害點j的需求量；ωi為車輛離開災(zāi)害點i后的載重量；K為所有車輛總和；Sk為車輛k的承載能力；R為車輛路徑總成本；Pi為第個災(zāi)害點被選上的概率；λi為期望到達災(zāi)害點i的時間。

式（1）是目標函數(shù)，用于最小化總成本，包括醫(yī)療救援車輛成本、路徑成本和懲罰成本.之后的約束條件分為3類，第1類是關(guān)于醫(yī)療救援車輛的，由式（2）～（9）決定；第2類是關(guān)于災(zāi)害點特征的，由式（10）～（17）決定；第3類是關(guān)于路徑特征的，由式（18）決定.式（2）和（3）表示醫(yī)療救援車輛在軟時間窗的約束下從救援站或者其中一個災(zāi)害點開出；式（4）和（5）表示所有車輛都會回到救援站；式（6）～（8）表示車輛的承載能力；式（9）表示在救援計劃完成前所有車輛都會返回到救援站；式（10）表示所有災(zāi)害點都會得到救助；式（11）和（12）表示每個災(zāi)害點僅由1部醫(yī)療救援車輛實施救助；式（13）決定了醫(yī)療救援車輛對任務(wù)分配的執(zhí)行條件；式（14）和（15）表示2個相鄰的災(zāi)害點由同一部救援車輛實施救助；式（16）和（17）用于計算由于違背軟時間窗而帶來的懲罰成本.方程式（13）決定路徑的約束條件.

3 模型的求解

3.1 初始化

醫(yī)療救援路徑通過一定的編碼來分析，以便災(zāi)害點能夠以一定的軟時間窗順序得到醫(yī)療救助，例如：有6個災(zāi)害點需要醫(yī)療救助，分別為D1，D2，D3，D4，D5，D6，如果令A(yù)0為救援站，而路徑編碼為（A0，D3，D5，D2，A0，D6，D1，D4，A0），可以看出，在軟時間窗的約束下有2條路徑可選，安排2部車輛分別實施救援，車輛1救援路徑為：A0→D3→D5→D2→A0，車輛2救援路徑為：A0→D6→D1→D4→A0，如果存在n部救援車輛，則每條染色體將有n條鏈接.在初始化階段，通過以下3個步驟產(chǎn)生可行的初始化結(jié)果.①分配災(zāi)害點到每個鏈接上，實質(zhì)是分組問題，災(zāi)害點將按與救援站的物理區(qū)域進行分組，目標是使總體的車輛行駛時間最小化.例如有2部車輛v1，v2用于醫(yī)療救援，救援站為A0，將每個災(zāi)害點i（xi，yi）按如下的原則進行分組：以救援站A0為圓心設(shè)定x，y軸，將x軸正值區(qū)的災(zāi)害點i（xi，yi）分配給v1，而x軸負值區(qū)的災(zāi)害點i（xi，yi）分配給v2，此時的災(zāi)害點i（xi，yi）與救援站的絕對距離為D（i，A0）＝[（xi－xA0）2＋（yi－yA0）2]－2.②對救援車輛的路徑進行分配，利用Clarke等[9]的節(jié)約矩陣E（i，j）可以構(gòu)建災(zāi)害點i，j與救援車輛的關(guān)聯(lián).首先假設(shè)在T0時刻2個災(zāi)害點分配在同一個路徑鏈接上，再比較成本節(jié)約值E（i，j）＝D（A0，i）＋D（A0，j）－D（i，j），在保證車輛承載能力的前提下，成本節(jié)約值大的災(zāi)害點將得到保留.③利用近鄰啟發(fā)算法解決車輛路徑問題，其原則是首先隨機選擇一個災(zāi)害點i，然后再選擇另一個最近的災(zāi)害點j形成路徑，直至所有的災(zāi)害點都覆蓋到為止.

通過以上步驟可以實現(xiàn)問題的初始化，參照圖1中的數(shù)據(jù)，可以得到具體初始化結(jié)果如圖2所示.

圖2 混合遺傳算法（HGA）初始化過程Fig.2 Initialization process of hybrid genetic algorithm

3.2 遺傳迭代交換

此時引入最優(yōu)搜索運算，對每2個災(zāi)害點進行遺傳交換比較，如果產(chǎn)生的結(jié)果比初始狀態(tài)好，則初始狀態(tài)將被取代從而形成父代，直至所有的災(zāi)害點的遺傳交換比較完成，最終產(chǎn)生新的父代。運算過程將消耗一定的時間，而迭代交換運算可以改進父代與子代之間的鏈接，需要通過以下5個步驟實現(xiàn)：①從父代中隨機選擇2個遺傳基因；②交換2個遺傳基因的位置從而形成子代；③互換2個遺傳基因的相鄰基因從而形成更多的子代；④評估所有子代并找到最好的一個；⑤如果最好的子代比父代好，則替換父代并重復(fù)步驟①，否則停止.圖3展示這種迭代交換過程.

3.3 救援時間評價

這里需要最小化車輛的救援時間，由于每輛車的救援完成時間各不相同，所以必需以最晚完成救援的車輛作為救援結(jié)束的時刻。令Tva為車輛v所量，而Tvm（Xω） 為 車 輛vi所 耗 費 的最長時間，Xω代表第ω個染色體，則有Tva＝，n為救援車輛的數(shù)max（Tv1a，Tv2a，…Tvna），式中：nv為所有救援車輛數(shù)目；t（i，j）為 從 災(zāi) 害 點i到j(luò)所 耗 費 時 間，為車輛行駛速度；nτ為在一條救援路徑上災(zāi)害點的數(shù)量.

圖3 混合遺傳算法（HGA）迭代交換過程Fig.3 Iterative exchange process of hybrid genetic algorithm

3.4 染色體選擇

這里需要選擇一些染色體用于遺傳運算，采用輪盤賭法選擇健康的染色體，越健康的染色體被選上的機會越大，也就是說，染色體被選上的機率與其健康程度成正比。染色體選擇通過以下步驟實現(xiàn)：①計算所有染色體的健康程度psize為被選染色體的規(guī)模；②計算每個染色體Xω被選擇上的概率pω＝[Ufit－Tvm（Xω）]／[Ufit（psize－1）]，ω＝1，2，…psize；③計算每個染色體Xω的累積概率Pω＝∑ωi＝1pi；④在范圍（0，1]中產(chǎn)生任意的隨機數(shù)λ，如果Pω≥λ＞Pω－1，則染色體Xω將被選擇.

3.5 遺傳運算

遺傳交叉算子可以找到更好的解決方案，而遺傳變異算子可以擴展更大的搜索空間。下文用交叉算子、啟發(fā)式變異和反轉(zhuǎn)突變來產(chǎn)生改進的子代染色體。

3.5.1 順序交叉

通過遺傳交叉每次可產(chǎn)生2個子代染色體，可以通過以下步驟實現(xiàn)：①從第1組父代中隨機選擇1組基因子串；②復(fù)制該基因子串到原子串對應(yīng)的位置，產(chǎn)生1組新的子串；③刪除第2組父代中相應(yīng)位置的基因，該位置的基因形成序列；④從左到右將原子串中空出位置按前步的序列順序填滿基因，從而形成一個新的子代；⑤通過交換2組父代的位置，重復(fù)步驟①～④的過程可產(chǎn)生另外一個子代.具體的遺傳交叉過程如圖4所示.

3.5.2 啟發(fā)式變異

可以通過以下步驟實現(xiàn)：①從父代中隨機選出3個遺傳基因；②產(chǎn)生被選擇基因所有可能的鄰里排列，并生成所有的子代.具體的啟發(fā)式變異過程如圖5所示.

3.5.3 反轉(zhuǎn)突變

逆算子作為一種變異運算僅僅用來增加樣本的多樣性，而不是加強樣本的質(zhì)量.具體的反轉(zhuǎn)突變過程如圖6.

4 數(shù)據(jù)仿真分析

2008年發(fā)生在我國四川汶川地區(qū)的里氏8.0級特大地震給災(zāi)區(qū)人民帶來嚴重的生命和財產(chǎn)損失，此時需要各個醫(yī)療救援站積極響應(yīng)，使醫(yī)療急救物資迅速運輸?shù)轿?而各個救援站的車輛資源是有限的，需要合理調(diào)度，故此以非常規(guī)突發(fā)災(zāi)害事件為案例，依據(jù)前面的數(shù)學(xué)模型通過數(shù)據(jù)仿真分析HGA算法的有效性并與 Azi等[10]的精確算法（exact algorithm，EA）進行比較，說明其在處理帶軟時間窗的路徑搜索問題方面的優(yōu)勢.EA算法可以精確到20個災(zāi)害點的情形，而災(zāi)害點大于20個時EA算法的計算時間會很長，在這里取某救援站所屬區(qū)域有20個災(zāi)害點的情形進行比較分析，HGA算法采用MATLAB 7.1編程、Intel Core P8400／2.26GHz處理器 運 算，而 EA 算 法 利 用 AMPL Plus 1.5／CPLEX求解.在整個仿真測試過程中，假設(shè)有2部救援車輛用于醫(yī)療救助，各車輛的承載能力與災(zāi)害點的坐標及軟時間窗都采用計算機隨機生成，而HGA的初始數(shù)據(jù)也通過計算機隨機生成，HGA迭代次數(shù)為200次，車輛行駛的時間間隔為2個單位，超出時間窗的懲罰成本W(wǎng)e和Wl各取20個單位.

首先比較HGA和EA算法的總成本缺口率Cgap＝（CHGA－CEA）／CEA×100%，其中CHGA，CEA分別為HGA和EA算法的總體成本消耗；運算時間的比率REH＝TEA／THGA，其中TEA，THGA分別為EA 和HGA算法的運算時間，于是可以計算比較2種算法的總成本消耗和運算時間，由于災(zāi)害點小于4的比較意義不大，故隨機采樣數(shù)據(jù)計算從4個災(zāi)害點開始，如表1所示.

從表1可以看出，當災(zāi)害點的數(shù)量小于7個時，HGA和EA算法的總成本是相當?shù)?，而隨著災(zāi)害點數(shù)量的增大，HGA和EA算法之間的總成本缺口發(fā)生變化，在20個災(zāi)害點時，總成本缺口達到了4.2%，在10個災(zāi)害點以內(nèi)時，HGA和EA算法的運算時間差別不太大，而隨著災(zāi)害點數(shù)量的增大，HGA算法體現(xiàn)出明顯的時間優(yōu)勢，EA算法的運算時間將大量增加，而在20個災(zāi)害點時，EA算法運算所耗費的時間是HGA算法的133.3倍.

表1 HGA和EA算法的總成本及時間消耗比較Tab.1 Comparison of the total cost and time by HGA and EA algorithm

隨著HGA算法迭代次數(shù)的增大醫(yī)療救援車輛的行駛時間進一步降低，從而優(yōu)化總成本消耗和救援運輸時間.下面分析HGA算法在20個災(zāi)害點的情況下，迭代運算200次，交叉率取0.4，突變率取0.2，隨機取3對染色體進行遺傳交叉，其中3個染色體作遺傳變異和反轉(zhuǎn)突變，于是每次迭代將產(chǎn)生18個子代，這樣重復(fù)迭代200次可以得到如圖7的結(jié)果.

圖7 HGA算法200次迭代收斂Fig.7 200 iterations converge of HGA algorithm

通過圖7可以看出，隨著HGA算法迭代次數(shù)的增加，救援車輛的行駛時間的優(yōu)化趨于穩(wěn)定，經(jīng)過200次的迭代運算后，救援車輛的實際行駛時間從最初的85個單位縮減到37個單位.

5 結(jié)語

面對突發(fā)性非常規(guī)災(zāi)害事件，醫(yī)療急救物資的及時運輸顯得尤為重要，將有益于減少突發(fā)事件下的傷亡.針對軟時間窗約束下的多輛醫(yī)療急救車輛的VRPTW問題建立了數(shù)學(xué)模型，得到最優(yōu)總成本消耗下對各個災(zāi)害點救援的車輛行駛路徑，同時有效縮短救援車輛的實際行駛時間.通過與EA算法進行數(shù)據(jù)仿真計算比較，結(jié)果表明HGA算法在總體成本消耗上與EA算法相當，但運算時間方面具有明顯的優(yōu)勢.

[1] Desaulniers G，Lavigne J，Soumis F.Multi－depot vehicle scheduling problems with time windows and waiting costs[J].European Journal of Operational Research， 1998， 111（3）：479.

[2] Qureshi A G，Taniguchi E，Yamada T.An exact solution approach for vehicle routing and scheduling problems with soft time windows[J].Transportation Research Part E：Logistics and Transportation Review，2009，45（6）：960.

[3] Ghoseiri K，Ghannadpour S F.Multi－objective vehicle routing problem with time windows using goal programming and genetic algorithm[J].Applied Soft Computing，2010，10（4）：1096.

[4] Taillard E D，Laporte G，Gendreau M.Vehicle routing with multiple use of vehicles[J].Journal of the Operational Research Society，1996，47：1065.

[5] Campbell A M，Savelsbergh M.Efficient insertion heuristics for vehicle routing and scheduling problems[J].Transportation Science，2004，38：369.

[6] Marinakis Y，Marinaki M，Dounias G.A hybrid particle swarm optimization algorithm for the vehicle routing problem[J].Engineering Applications of Artificial Intelligence，2010，23（4）：463.

[7] Tan K C，Lee L H，Ou K.Artificial intelligence heuristics in solving vehicle routing problems with time window constraints[J].Engineering Applications of Artificial Intelligence，2001，14（6）：825.

[8] Park Y B.A hybrid genetic algorithm for the vehicle scheduling problem with due times and time deadlines[J].International Journal of Production Economics，2001，73（2）：175.

[9] Clarke G，Wright J W.Scheduling of vehicles from a central depot to a number of delivery points[J].Operations Research，1964，12：568.

[10] Azi N，Gendreau M，Potvin J Y.An exact algorithm for a vehicle routing problem with time windows and multiple use of vehicles[J].European Journal of Operational Research，2010，202（3）：756.