孫 亮潘全科鄒溫強(qiáng)王亞敏

(1.上海大學(xué)機(jī)電工程與自動(dòng)化學(xué)院,上海 200072;2.山東理工大學(xué)交通與車輛工程學(xué)院,山東淄博 255049;3.南京審計(jì)大學(xué)信息工程學(xué)院,江蘇南京 211815)

1 引言

車輛路徑問題(vehicle routing problem,VRP)是指在滿足一定的約束條件(如需求量、車輛容量等)下,配送中心構(gòu)造有限條車輛運(yùn)行路徑為客戶提供配送服務(wù),最終實(shí)現(xiàn)客戶服務(wù)水平和運(yùn)輸成本等決策目標(biāo)的最優(yōu)化.配送完成時(shí)間是衡量客戶服務(wù)水平的一個(gè)重要因素.但配送完成時(shí)間受交通擁堵等多種因素的影響,很難通過精確的概率分布進(jìn)行描述,使得隨機(jī)規(guī)劃方法[1]處理不確定旅行時(shí)間條件下的車輛路徑問題(vehicle routing problem with uncertain travel times,VRP–UT)的效果不夠理想.與隨機(jī)規(guī)劃方法不同,魯棒優(yōu)化方法[2]利用不確定數(shù)據(jù)集的邊界關(guān)系,可以將不確定的優(yōu)化模型轉(zhuǎn)化為確定的優(yōu)化模型進(jìn)行求解.目前常用的最壞場景魯棒優(yōu)化方法[3],要求所有車輛均需要在客戶指定的時(shí)間內(nèi)將貨物送達(dá)客戶,但從實(shí)際營運(yùn)角度看,稍微偏離客戶允許時(shí)間范圍的配送服務(wù)并不對服務(wù)水平構(gòu)成影響.輕魯棒優(yōu)化方法[4–7]通過對魯棒優(yōu)化模型中的某些約束進(jìn)行松弛,以約束違背最小化為優(yōu)化目標(biāo),從而改善最壞實(shí)現(xiàn)下目標(biāo)函數(shù)值與預(yù)期值之間的偏差.所以采用輕魯棒優(yōu)化方法研究VRP–UT,更具有理論意義和實(shí)際價(jià)值.

從求解算法看,求解此類問題的算法主要有鄰域搜索算法[8–10]、分支定界算法[11–13]和超啟發(fā)式算法[14]等.現(xiàn)有算法主要關(guān)注總變動(dòng)成本的降低,但對如何在總變動(dòng)成本降低的同時(shí),使盡可能多的客戶在其允許的時(shí)間范圍完成配送服務(wù)關(guān)注較少.基于上述分析,建立了一種描述VRP–UT特征的輕魯棒優(yōu)化模型(light-robust-optimization model,LRO),提出了一種超啟發(fā)式粒子群算法(hyper particle swarm optimization,HPSO).大量的仿真實(shí)驗(yàn)表明了所提模型和算法的有效性.

2 問題描述

VRP–UT描述如下:給定完全賦權(quán)圖G(V,E)表示配送中心與客戶點(diǎn)之間的相對位置關(guān)系,其中V{0,1,2,……,n}表示點(diǎn)集,

表示邊集,0 表示配送中心;CV{0}表示顧客集.配送中心擁有相同容量的m輛車,每輛車的最大負(fù)載均為Q,車輛集合用K{1,2,……,m}表示.配送車輛的單位懲罰成本ω1和單位等待成本ω2,第i個(gè)顧客的需求量di,客戶允許的服務(wù)時(shí)間窗(客戶時(shí)間窗)[ETi,LTi].關(guān)于旅行時(shí)間的不確定數(shù)據(jù)集UT定義如下:

車輛從客戶i到客戶j的旅行時(shí)間在之間任意取值.

根據(jù)已知決策信息,獲取一個(gè)調(diào)度方案,該調(diào)度方案具有兩個(gè)特征:①對UT內(nèi)的任意實(shí)現(xiàn)均保持可行;②該調(diào)度方案對應(yīng)的TVC最小.

條件假設(shè)如下:車輛從配送中心出發(fā),服務(wù)完最后一個(gè)客戶后,必須回到配送中心;車輛進(jìn)入某個(gè)客戶點(diǎn),也必須從該點(diǎn)離開;一輛車可以為多個(gè)客戶服務(wù),但一個(gè)客戶只能被一輛車服務(wù);車輛的實(shí)際負(fù)載不得超過Q.

3 VRP–UT輕魯棒優(yōu)化模型

在文獻(xiàn)[15]給出的魯棒優(yōu)化模型的基礎(chǔ)上,利用新增加的松弛變量表示變動(dòng)成本,構(gòu)建描述該問題特征的輕魯棒優(yōu)化模型.模型中M表示一個(gè)充分大的正數(shù).LRO模型如下:

目標(biāo)函數(shù)(3)表示TVC;約束(4)–(5)表示每輛車必須從配送中心出發(fā),服務(wù)完最后一個(gè)客戶后,回到配送中心;約束(6)–(8)表示車輛服務(wù)完某個(gè)客戶后,必須從該客戶離開;約束(9)表示配送中心允許使用的最大車輛數(shù)限制;約束(10)表示每個(gè)客戶只能由一輛車進(jìn)行服務(wù);約束(11)表示每一輛車的實(shí)際負(fù)載不能超過Q;約束(12)為避免出現(xiàn)子循環(huán)(客戶點(diǎn)之間形成的環(huán)路)約束;約束(13)–(14)表示決策變量之間的關(guān)系;約束(15)表示可行解對于UT的任意實(shí)現(xiàn)均保持可行;約束(16)–(17)表示車輛到達(dá)客戶的時(shí)間幾乎處處在[ETi,LTi]內(nèi);約束(18)–(24)表示決策變量的取值范圍.

4 LRO的求解方法

從LRO可行解滿足的條件可以看出,LRO可行解的每一個(gè)子回路都是一個(gè)環(huán),而深度優(yōu)先搜索(depth first search,DFS)[16]是搜索圖中是否存在回路的有效的遍歷方式.因此,以是否存在零環(huán)(可行解中某條子回路的TVC為零,則稱該子回路為零環(huán))為搜索目的的DFS可以發(fā)現(xiàn)零環(huán).與此同時(shí),使用DFS生成零環(huán)時(shí),當(dāng)遍歷到某個(gè)結(jié)點(diǎn)的前驅(qū)結(jié)點(diǎn)為已經(jīng)訪問過的結(jié)點(diǎn)時(shí),將不再訪問這個(gè)結(jié)點(diǎn),因此可以完全避免子循環(huán)的產(chǎn)生.

圖1 HPSO算法流程圖Fig.1 Flow chart of HPSO

因?yàn)镻SO可以兼顧尋優(yōu)過程的全局探測和局部開采能力并且沒有許多的參數(shù)調(diào)節(jié),所以采用PSO作為外層搜索算法.HPSO的算法流程如圖1所示.

4.1 DFS的實(shí)現(xiàn)方法

1) 可行解的編碼方式(encoding of feasible solutions).

可行解采用基于客戶順序的編碼方法,如圖2所示.0表示配送中心.其他數(shù)字表示客戶點(diǎn).例如,路徑(0 2 6 7 9 4 0)表示車輛服務(wù)客戶的順序是0→2→6→7→9→4→0.

圖2 可行解編碼結(jié)構(gòu)示意圖Fig.2 Feasible solution representation

2) DFS的算法步驟(pseudo-code of DFS).

DFS由兩層循環(huán)構(gòu)成.內(nèi)層循環(huán)利用DFS,產(chǎn)生一個(gè)零環(huán).在內(nèi)層循環(huán)中,從0點(diǎn)開始,如果發(fā)現(xiàn)一個(gè)未訪問的及時(shí)點(diǎn),該點(diǎn)不是0點(diǎn)并且該點(diǎn)不違背約束(9),該點(diǎn)加入第k輛車對應(yīng)的服務(wù)序列中,并從該點(diǎn)繼續(xù)搜索,直到形成第k個(gè)零環(huán)zk.將發(fā)現(xiàn)的及時(shí)點(diǎn)從客戶服務(wù)序列MSi1,……,ik,……,in中移除.在外層循環(huán)中,將zk加入所有子回路的集合R中.重復(fù)此過程,直到無法在MS 中發(fā)現(xiàn)及時(shí)點(diǎn).剩余的點(diǎn)根據(jù)約束式(11)插入到R中.sumdemand(zk)表示zk的總需求.Γ(·)表示某個(gè)點(diǎn)的鄰接點(diǎn).DFS(MS)算法的偽代碼如下:

4.2 低層啟發(fā)式規(guī)則變換方式設(shè)計(jì)

根據(jù)幾何不等式,如果每條子路徑上的TVC盡可能相同,那么TVC將得到進(jìn)一步改善.從這個(gè)角度出發(fā),對文獻(xiàn)[14]中提出的啟發(fā)式規(guī)則變換方式進(jìn)行重新設(shè)計(jì).新的啟發(fā)式規(guī)則變換方式的具體實(shí)施步驟如下:

2) LRO–insert.分別從TVC 最大和最小的兩個(gè)環(huán)中隨機(jī)選擇位置qm和qh,將客戶服務(wù)序列中對應(yīng)的元素從總變動(dòng)成本最大的路徑中移除,然后將該位置的元素插入到客戶服務(wù)序列中qm對應(yīng)的元素前面;

3) LRO–inverse.首先,隨機(jī)選擇一個(gè)環(huán).其次,從該環(huán)中隨機(jī)選擇兩個(gè)位置kr和kq,將該序列中kr和kq之間的所有位置對應(yīng)的元素kr,kr+1,……,kq?1,kq變換為kq,kq?1,……,kr+1,kr.

4.3 粒子的編碼與解碼

4.3.1 粒子的編碼與解碼

粒子是一個(gè)1×3q維向量.采用1,2,3分別表示LRO–swap,LRO–inverse,LRO–insert這3種變換方式.第1,4,……,3q ?2維表示變換方式,其他維度表示變換位置.q表示啟發(fā)式規(guī)則的數(shù)量.例如,(1,2,8)表示將初始可行解對應(yīng)的客戶服務(wù)序列中的第2位和第8位采用LRO–swap的方式進(jìn)行交換.

4.3.2 粒子的解碼方法

記P(r1,……,rq)表示位置向量.依次將啟發(fā)式規(guī)則ri作用于該粒子對應(yīng)的可行解,得到最終的客戶服務(wù)序列后,再使用DFS將該序列轉(zhuǎn)化為一個(gè)可行解計(jì)算的TVC.

4.4 更新規(guī)則

因?yàn)槲恢孟蛄勘硎镜牡蛯訕?gòu)造算法是根據(jù)給定的啟發(fā)式規(guī)則的變換方式生成,而不是對給定啟發(fā)式規(guī)則的進(jìn)行對換,所以現(xiàn)有用于VRP求解的粒子群優(yōu)化算法(particle swarm optimization,PSO)的更新規(guī)則不能直接應(yīng)用于LRO的求解.因此需要對更新公式進(jìn)行重新設(shè)計(jì).重新設(shè)計(jì)的更新式如下:

從結(jié)構(gòu)變動(dòng)度來看，至2016年，真正體現(xiàn)醫(yī)務(wù)人員服務(wù)價(jià)值的費(fèi)用項(xiàng)目，如手術(shù)費(fèi)、護(hù)理費(fèi)，這兩項(xiàng)明細(xì)費(fèi)用項(xiàng)目的構(gòu)成比較往年有所增加并達(dá)到五年來的最高值，反映了該院在提升醫(yī)者服務(wù)價(jià)值、提高醫(yī)務(wù)人員積極性方面已有一定進(jìn)步；但灰色關(guān)聯(lián)數(shù)據(jù)顯示，手術(shù)費(fèi)及護(hù)理費(fèi)兩項(xiàng)指標(biāo)的關(guān)聯(lián)系數(shù)逐年下降，這說明二者與住院次均費(fèi)用之間關(guān)系的緊密程度有所減弱。建議衛(wèi)生管理者在制定完善醫(yī)療服務(wù)價(jià)格調(diào)整政策時(shí)，參考國內(nèi)相關(guān)的標(biāo)化價(jià)值方法學(xué)體系[11]，明確技術(shù)勞務(wù)、物資耗材等項(xiàng)目的價(jià)值構(gòu)成，逐步理順醫(yī)療服務(wù)價(jià)格，充分調(diào)動(dòng)醫(yī)務(wù)人員工作積極性；而醫(yī)療機(jī)構(gòu)應(yīng)繼續(xù)提高醫(yī)生業(yè)務(wù)收入中技術(shù)勞務(wù)性收入比重，進(jìn)一步優(yōu)化費(fèi)用結(jié)構(gòu)。

其中:t表示當(dāng)前迭代次數(shù);ceil(·)表示向上取整函數(shù);mod(·)為取余函數(shù);Xi(t){xi1(t),……,xi3m(t)}為粒子i在第t代的位置向量;

為粒子i在第t代的速度向量;為當(dāng)前迭代次數(shù)下,種群中適應(yīng)度函數(shù)值中最好的個(gè)體位置向量;XG為截至到當(dāng)前迭代次數(shù),種群中所有中最好的個(gè)體位置向量;ω為慣性權(quán)重;T為迭代的最大次數(shù);M為種群規(guī)模;a1,a2為加速因子;D為隨機(jī)生成[?(T ?t),(T ?t)]之間的3m個(gè)正整數(shù).

4.5 HPSO的算法步驟

首先,使用DFS方法生成一組初始可行解.其次,隨機(jī)生成初始粒子群.將可行解分配到每一個(gè)粒子中.然后,進(jìn)入迭代階段,在該階段,利用重新設(shè)計(jì)的更新公式,進(jìn)行迭代,根據(jù)TVC最小的原則,更新當(dāng)前最優(yōu)解.當(dāng)滿足迭代終止條件時(shí),輸出最終的最優(yōu)解.HPSO(LRO)算法的偽代碼如下所示:

利用DFS方法得到一組初始可行解

5 仿真分析

5.1 算例選取與參數(shù)設(shè)置

從文獻(xiàn)[16]中選取若干標(biāo)準(zhǔn)算例,引入U(xiǎn)T的參數(shù)信息,構(gòu)造新算例用于分析LRO和HPSO的相關(guān)性能.算例的基本信息詳見表1.

表1 標(biāo)準(zhǔn)算例的基本信息Table 1 Configure information of benchmarks

tij表示標(biāo)準(zhǔn)算例中到的旅行時(shí)間.令ω1,a11,a21,M30,T30,q10,ω11,ω21,Vmin?30,Vmax30,ρij1.參數(shù)取上述值的原因是便于觀察合理的計(jì)算時(shí)間內(nèi)算法優(yōu)化解的相關(guān)性能.為了避免正態(tài)分布和方差齊性假設(shè)對分析結(jié)果的影響,本文選擇Kruskal-Wallis ANOVA檢驗(yàn)和趨勢卡方檢驗(yàn)對算法優(yōu)化解的相關(guān)性能進(jìn)行分析.

5.2 仿真環(huán)境

采用MATLAB R2016a(64bit)實(shí)現(xiàn)了求解LRO 的HPSO及其他對比算法,在Intel(R)Core(TM)i5–6200U CPU@2.3 GHz,8 GB內(nèi)存的計(jì)算機(jī)上進(jìn)行了仿真實(shí)驗(yàn).相關(guān)統(tǒng)計(jì)學(xué)檢驗(yàn)使用SPSS(version 26)實(shí)現(xiàn)(SPSS的全稱為statistical product and service solutions,即:統(tǒng)計(jì)產(chǎn)品與服務(wù)解決方案).

5.3 LRO與RO比較分析

文獻(xiàn)[15]給出的魯棒優(yōu)化模型(robust optimization model),記為RO.兩者在模型特征上的差異,可能造成兩者最優(yōu)解(TVC為0的解)對UT參數(shù)的靈敏度不同.由于LRO與RO的目標(biāo)函數(shù)不同,對于RO,只利用其約束條件發(fā)現(xiàn)最優(yōu)解,對兩者的最優(yōu)解進(jìn)行靈敏度分析.

從文獻(xiàn)[16]中選取35個(gè)標(biāo)準(zhǔn)算例進(jìn)行最優(yōu)解的靈敏度分析.選取算例及其特征如表2所示.參數(shù)信息如表3所示.采用HPSO對表2中的每個(gè)算例運(yùn)算30次,記錄各自的TVC.表2中加粗部分表示該算例兩種魯棒優(yōu)化方法均發(fā)現(xiàn)理想最優(yōu)解.分別將LRO和RO獲得的最優(yōu)解代入不同的Λ和θij,觀察其TVC的變化.Λ從0.25n增加到n,增加幅度為0.125n.ρij1,maxt表示算例中tij的最大值.θij從maxt增加到4 maxt,增加幅度為0.5 maxt.

表3–4說明,當(dāng)Λ和θij發(fā)生變化時(shí),RO在最壞實(shí)現(xiàn)下的最優(yōu)解,其TVC隨著Λ和θij的增加而改變,并且這種改變具有統(tǒng)計(jì)學(xué)意義(P <0.05).LRO在最壞情況下的最優(yōu)解,雖然也發(fā)生一定改變,但這種改變并沒有統(tǒng)計(jì)學(xué)意義(P >0.05).與RO不同,LRO是通過不斷調(diào)整偏差值產(chǎn)生的最優(yōu)解,這可能是LRO對參數(shù)變化不敏感的一個(gè)重要原因.

表2 算例類型及其算例對應(yīng)表Table 2 Feature of benchmark

表3 Λ的變化對TVC的影響(θijρij 1.5 max t,M30,T 30)Table 3 Λ impact on TVC (θijρij 1.5 max t,M 30,T 30)

表3 Λ的變化對TVC的影響(θijρij 1.5 max t,M30,T 30)Table 3 Λ impact on TVC (θijρij 1.5 max t,M 30,T 30)

表4 θij的變化對TVC的影響(Λ0.75?n,M30,T 30)Table 4 θij impact on TVC(Λ0.75?n,M30,T30)

表4 θij的變化對TVC的影響(Λ0.75?n,M30,T 30)Table 4 θij impact on TVC(Λ0.75?n,M30,T30)

5.4 HPSO的求解性能分析

因?yàn)槲墨I(xiàn)[9]中的求解算法(簡記為AMP)是目前普遍認(rèn)為求解此類魯棒優(yōu)化模型性能非常理想的啟發(fā)式算法,所以本文采用AMP與HPSO進(jìn)行對比試驗(yàn).為了避免初始解生成策略對分析結(jié)果的影響,兩者均采用DFS方法生成初始解.兩個(gè)算法遍歷的可行解總數(shù)相同.HPSO的參數(shù)設(shè)置詳見第4.1節(jié).針對表5中的測試數(shù)據(jù),使用HPSO和AMP分別獨(dú)立求解30次LRO.將30次運(yùn)算結(jié)果的運(yùn)行時(shí)間(t),TVC的95%置信區(qū)間記錄于表5中.針對表5中描述的數(shù)據(jù)進(jìn)行Kruskal-Wallis ANOVA檢驗(yàn),相關(guān)結(jié)論見表6.

表6中關(guān)于兩種算法的TVC的平均秩次說明,HPSO獲取的TVC顯著小于AMP獲取的TVC(p<0.05).HPSO獲取的TVC能夠顯著改善的原因是:AMP根據(jù)給定的啟發(fā)式規(guī)則產(chǎn)生可行解;HPSO在搜索過程中不斷生成新規(guī)則產(chǎn)生可行解.因此,HPSO搜索范圍更加廣泛,從而提升了算法優(yōu)化解的性能.

表5 TVC與運(yùn)算時(shí)間分析Table 5 Detailed results of TVC and Time

表6中關(guān)于兩種算法運(yùn)行時(shí)間的平均秩次說明,兩種算法在求解LRO時(shí)運(yùn)行時(shí)間不存在統(tǒng)計(jì)學(xué)差異(p>0.05),造成這種現(xiàn)象的原因是兩者的算子本質(zhì)上都屬于對換,在時(shí)間復(fù)雜度上是相同的.

表6 針對TVC和時(shí)間的平均秩次比較分析Table 6 Results of mean ranks with regard to TVC and time

6 結(jié)論

本文的核心內(nèi)容如下:①模型方面,以TVC最小為優(yōu)化目標(biāo),構(gòu)建描述該問題特征的LRO;②算法方面,設(shè)計(jì)HPSO求解LRO,該算法利用DFS策略生成初始解,然后再利用重新設(shè)計(jì)的更新公式,生成構(gòu)造算法,進(jìn)一步降低算法優(yōu)化解對應(yīng)的TVC;③仿真實(shí)驗(yàn)表明,HPSO是一種求解LRO的有效算法.