突發(fā)災(zāi)害下帶軟時(shí)間窗多車路徑搜索建模

劉天虎，許維勝，吳啟迪

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

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

1 問(wèn)題描述

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

2 模型的建立

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

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

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

當(dāng)T1時(shí)刻車輛1行駛到災(zāi)害點(diǎn)D5，而車輛2行駛到D10時(shí)，此時(shí)收到醫(yī)療救援站的最新信息，又有4個(gè)災(zāi)害點(diǎn)D11，D12，D13，D14需要救助，從圖1b可以看出，在規(guī)劃的路徑范圍內(nèi)，災(zāi)害點(diǎn)D11和D13可由車輛1提供救助，而災(zāi)害點(diǎn)D12和D14可由車輛2提供救助，在不超過(guò)車載容量的前提下，對(duì)2部醫(yī)療救援車輛的后續(xù)行駛路徑進(jìn)行了優(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í)間窗下的路徑搜索問(wèn)題實(shí)際是一個(gè)混合整數(shù)線性規(guī)劃，于是可以建立相應(yīng)的數(shù)學(xué)模型，模型的目標(biāo)是為了最小化救援的總成本（包括車輛成本、路徑成本和懲罰成本）.而該模型的限制條件由車輛、路徑、災(zāi)害點(diǎn)時(shí)間窗所約束，于是可以建立模型的最優(yōu)方程z如下：

限制條件如下：

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

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

3 模型的求解

3.1 初始化

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

通過(guò)以上步驟可以實(shí)現(xiàn)問(wèn)題的初始化，參照?qǐng)D1中的數(shù)據(jù)，可以得到具體初始化結(jié)果如圖2所示.

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

3.2 遺傳迭代交換

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

3.3 救援時(shí)間評(píng)價(jià)

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

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

3.4 染色體選擇

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

3.5 遺傳運(yùn)算

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

3.5.1 順序交叉

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

3.5.2 啟發(fā)式變異

可以通過(guò)以下步驟實(shí)現(xiàn)：①?gòu)母复须S機(jī)選出3個(gè)遺傳基因；②產(chǎn)生被選擇基因所有可能的鄰里排列，并生成所有的子代.具體的啟發(fā)式變異過(guò)程如圖5所示.

3.5.3 反轉(zhuǎn)突變

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

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

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

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

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

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

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

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

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

5 結(jié)語(yǔ)

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

[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.