基于抽樣和變鄰域搜索的隨機(jī)共享單車重平衡問題

賈永基， 惲博文， 許媛媛

(東華大學(xué) 旭日工商管理學(xué)院， 上海 200051)

隨著我國經(jīng)濟(jì)的迅猛發(fā)展，交通運輸需求快速增加，未來仍將長期呈現(xiàn)增長態(tài)勢，這導(dǎo)致交通運輸系統(tǒng)的碳排放所占份額越來越大。在2021年全國兩會上，“碳達(dá)峰”與“碳中和”首次被寫入政府工作報告。為了如期實現(xiàn)“雙碳”目標(biāo)，降低交通系統(tǒng)碳排放，須大力推廣零排放的共享單車系統(tǒng)。近20年來，共享單車作為一種低碳、環(huán)保的短途出行方式在全球范圍內(nèi)得到廣泛應(yīng)用，在緩解交通擁堵、構(gòu)建綠色出行體系方面做出了極大貢獻(xiàn)[1]。然而，隨著共享單車規(guī)模的不斷擴(kuò)大，供需不匹配現(xiàn)象愈加嚴(yán)重，在加大企業(yè)運營成本的同時，降低用戶滿意度，這將極大阻礙共享單車行業(yè)的健康發(fā)展。因而實施有效的共享單車供需重新匹配操作，即共享單車重平衡問題(bike-sharing rebalance problem，BRP)[2]，顯得尤為重要。

在目前的研究成果中，所有參數(shù)都已知的確定性BRP研究已較為成熟[3]，而在實際運營場景中，客戶需求、行駛時間等參數(shù)都不可避免地會出現(xiàn)波動，因而隨機(jī)BRP逐漸引起研究人員的重視。目前隨機(jī)BRP的主要求解方法大致可分為4類：需求預(yù)測、馬爾可夫決策過程、滾動時域和隨機(jī)規(guī)劃。Zhang等[4]在預(yù)測單車庫存和用戶需求的基礎(chǔ)上建立BRP模型。靳文舟等[5]基于用戶預(yù)約數(shù)據(jù)進(jìn)行需求預(yù)測，構(gòu)建了考慮用戶滿意度的最小化總成本的BRP模型。陳菁等[6]基于小波神經(jīng)網(wǎng)絡(luò)預(yù)測區(qū)域需求量構(gòu)建BRP模型。Federico等[7]利用歷史數(shù)據(jù)來預(yù)測共享單車網(wǎng)絡(luò)狀況并在必要時采取重平衡操作。Pan等[8]建立了基于馬爾可夫決策過程的BRP模型，并提出一種基于確定性策略梯度算法的深度學(xué)習(xí)算法。Legros[9]則使用馬爾可夫決策過程開發(fā)一種可實施的決策支持工具，幫助運營商在任何時間決定站點優(yōu)先級及單車庫存。Zhai等[10]同樣利用馬爾可夫決策過程和線性規(guī)劃方法來解決共享單車車隊規(guī)模和BRP問題。董紅召等[11]加入模糊時間窗約束，以站點滿意度最大為目標(biāo)建立公共自行車系統(tǒng)調(diào)度模型，采用滾動時域算法來求解該模型。Shui等[12]也采用滾動時域獲得多個階段的靜態(tài)子問題，使用改進(jìn)的人工蜂群算法和路線截斷算法來優(yōu)化每個階段的重平衡車輛路徑。Maggioni等[13]考慮共享單車隨時間變化的隨機(jī)需求，提出兩階段和多階段隨機(jī)規(guī)劃模型，以確定在次日服務(wù)開始時每個站點的最優(yōu)單車數(shù)量。

現(xiàn)實BRP系統(tǒng)中，客戶隨機(jī)取車、還車，使得各站點的單車數(shù)量時刻處于變化中，因而無法提前獲得較準(zhǔn)確的預(yù)測數(shù)據(jù)，導(dǎo)致基于需求預(yù)測的重平衡模型的計算結(jié)果偏離實際值較大。馬爾可夫決策過程計算復(fù)雜，因而其不適合大規(guī)模的實際問題求解。滾動時域采用“wait and see”策略等待隨機(jī)需求的確定值出現(xiàn)，導(dǎo)致重平衡的時效性較差。隨機(jī)規(guī)劃是一種事前策略，同時不要求事先知道所有參數(shù)的確定值，因而具有更大的使用范圍。據(jù)筆者所知，目前還沒有文獻(xiàn)使用隨機(jī)規(guī)劃方法來研究需求隨機(jī)的共享單車重平衡問題(bike-sharing rebalancing problem with stochastic demands，BRPSD)。

使用共享單車企業(yè)的歷史數(shù)據(jù)，可以估計客戶隨機(jī)需求的分布概率。為了簡化計算，本文假設(shè)客戶需求的分布概率是已知的，因而可以使用拉丁超立方抽樣(Latin hypercube sampling，LHS)來生成一系列具有給定概率的離散場景來表示客戶需求的隨機(jī)性。在每個場景下，客戶需求是固定不變的，因而可以使用常規(guī)啟發(fā)式算法(例如變鄰域搜索算法)進(jìn)行求解。LHS是一種基于蒙特卡洛模擬的方法。L?hndorf[15]研究流行的隨機(jī)場景生成方法，指出蒙特卡洛方法是目前最流行的方法，但其需要大量樣本以確保隨機(jī)變量概率分布的準(zhǔn)確近似，該缺陷可以使用LHS加以緩解。Nikzad等[16]已證實基于LHS獲得的結(jié)果更接近基于大樣本獲得的結(jié)果，且在相同數(shù)量的場景中，LHS比蒙特卡洛方法可獲得更準(zhǔn)確的近似值。因此，本文采用LHS來生成一系列隨機(jī)場景。此外，本文重新定義共享單車站點的平衡狀態(tài)，提出了表示站點平衡狀態(tài)的平衡區(qū)間概念，給出了衡量站點重平衡程度的重平衡效應(yīng)度函數(shù)。在此基礎(chǔ)上，建立BRPSD的兩階段隨機(jī)規(guī)劃模型，并提出了求解該模型的基于LHS的變鄰域搜索算法。

1 問題描述與數(shù)學(xué)模型

1.1 問題描述

在共享單車系統(tǒng)中，由于單車的單程性所導(dǎo)致的“潮汐現(xiàn)象”會使單車供需不匹配，從而導(dǎo)致用戶滿意度和單車?yán)寐实南陆担焕诠蚕韱诬囆袠I(yè)的可持續(xù)發(fā)展。共享單車重平衡是指利用運輸車輛將單車從過剩站點運往不足站點，即重新在站點之間分配單車以達(dá)到供需平衡。如何經(jīng)濟(jì)、高效地實現(xiàn)共享單車重平衡也就成為共享單車行業(yè)亟待解決的關(guān)鍵問題。

在共享單車企業(yè)實際運營中，為明確工作職責(zé)，往往將單車服務(wù)區(qū)域劃分成多個不關(guān)聯(lián)的子區(qū)域，每個子區(qū)域只由一輛運輸車單獨進(jìn)行重平衡操作?？紤]到隨機(jī)重平衡問題的復(fù)雜性，不失一般性，本文只考慮一輛運輸車輛的情形，且假設(shè)每個站點只能被服務(wù)一次。BRPSD問題可以描述為：車輛從車庫出發(fā)，依次訪問處于不平衡狀態(tài)的共享單車站點，通過取車或卸車操作使得所有站點趨于平衡狀態(tài)，最后返回車庫，同時考慮車輛行駛距離的最小化和重平衡效用度期望的最大化。

本文所使用的決策變量和符號如表1所示。

表1 決策變量和符號

1.2 平衡區(qū)間和重平衡效應(yīng)度

現(xiàn)有研究多將共享單車站點的平衡狀態(tài)設(shè)為一個平衡點，但由于客戶取車、還車活動的隨機(jī)性，站點的單車數(shù)量是波動的，因而平衡點的設(shè)置使得重平衡操作不具有穩(wěn)定性，提高了重平衡成本。因此，本文將站點的平衡狀態(tài)由平衡點改為平衡區(qū)間[14]，如圖1所示，其中aj是站點j的當(dāng)前單車數(shù)量，bj和[(1-θ)bj,(1+θ)bj]分別是其平衡點和平衡區(qū)間。當(dāng)站點的單車數(shù)量在其平衡區(qū)間內(nèi)時，該站點就處于平衡狀態(tài)。即使站點單車數(shù)量由于客戶活動出現(xiàn)了隨機(jī)波動，只要波動范圍不超過平衡區(qū)間閾值，就無需進(jìn)行重平衡操作，這可以大大減少共享單車系統(tǒng)的重平衡需求，從而提高共享單車企業(yè)的重平衡效率。

圖1 共享單車系統(tǒng)的平衡區(qū)間

由于站點單車數(shù)量的隨機(jī)性，本模型不強(qiáng)制要求所有站點在被服務(wù)后都處于平衡狀態(tài)，而采用重平衡效用度[14]來衡量重平衡操作的效益，如圖2所示。在重平衡操作后，若站點j的單車數(shù)量處于平衡區(qū)間[(1-θ)bj,(1+θ)bj]內(nèi)，那么其重平衡效用度為1。當(dāng)單車數(shù)量小于平衡區(qū)間下界(1-θ)bj或大于平衡區(qū)間上界(1+θ)bj時，其重平衡效用度隨aj與平衡區(qū)間差值的逐漸增大而逐漸減小。

圖2 站點的重平衡效用度

1.3 兩階段隨機(jī)規(guī)劃模型

基于LHS的兩階段隨機(jī)規(guī)劃模型的第一階段為運輸車輛路徑優(yōu)化模型，最大化共享單車系統(tǒng)的總效益；第二階段基于第一階段所獲得的車輛路徑信息，最大化重平衡效用度期望。第二階段基于第一階段的路徑信息，又影響著第一階段的路徑?jīng)Q策。

第一階段模型(P1)：

(1)

s.t.

(2)

(3)

?S?V,2≤|S|≤n-1

(4)

xij∈{0,1} ?i,j∈V,i≠j

(5)

其中：式(1)表示共享單車系統(tǒng)的總效益，重平衡操作產(chǎn)生正效益，系數(shù)為ε1，而車輛行駛成本產(chǎn)生負(fù)效益，系數(shù)為ε2；式(2)表示每個站點只能被訪問一次；式(3)表示流量平衡；式(4)用于消除子回路；式(5)表示決策變量的取值范圍。

第二階段模型(P2)：

(6)

(7)

?j∈N,ω∈Ω

(8)

lj,ω2=lj,ω1-sj,ω++sj,ω-?j∈N,ω∈Ω

(9)

sj,ω+×sj,ω-=0 ?j∈N,ω∈Ω

(10)

(11)

(12)

(13)

sj,ω+，sj,ω-，fij,ω≥0 ?i,j∈V,i≠j,ω∈Ω

(14)

其中：式(6)表示最大化隨機(jī)場景下的重平衡效用度期望；式(7)表示場景ω下運輸車輛的裝載量不大于車輛容量；式(8)表示場景ω下車輛服務(wù)站點前后的裝載量之差等于站點被滿足的單車數(shù)量；式(9)表示場景ω下站點被服務(wù)前后的單車數(shù)量之差等于站點被滿足的單車數(shù)量；式(10)表示場景ω下站點不能同時具有取車和卸車需求；式(11)表示場景ω下的站點重平衡效用度；式(12)表示決策變量sj,ω+不能超過車輛訪問站點時的剩余容量；式(13)表示決策變量sj,ω-不能超過車輛訪問站點時的裝載量；式(14)表示決策變量的取值范圍。

2 基于LHS的變鄰域搜索算法

2.1 算法流程

BRPSD模型的解由兩部分決定：第一階段的車輛路徑和第二階段一系列隨機(jī)場景下的重平衡需求。其中，車輛路徑可表示為(0,i1,i2,…,in,0)，其中0是車庫，而i1,i2,…,in是車輛所訪問站點的一個排列。通過改變站點的順序，可以獲得車輛的不同行駛路徑。顯然，可以通過求解一個旅行商問題來獲得第一階段的最優(yōu)車輛路徑。本文設(shè)計了基于拉丁超立方抽樣的變鄰域搜索(Latin hypercube sampling based variable neighbourhood search，LHS-VNS)算法來求解BRPSD問題，該算法流程如圖3所示。

圖3 LHS-VNS算法流程圖

首先，采用LHS生成一系列隨機(jī)場景，每個場景都具有相同的概率。接著，使用商業(yè)優(yōu)化軟件CPLEX求解第一階段模型P1，獲得車輛最優(yōu)路徑。然后，將車輛最優(yōu)路徑代入第二階段模型P2，求解該路徑在所有場景下完成重平衡操作后達(dá)到的重平衡效用度期望的最大值。如果沒有達(dá)到停止條件，則利用變鄰域搜索(variable neighbourhood search，VNS)算法來獲得新的車輛路徑，并將其代入模型P2進(jìn)行計算，否則結(jié)束迭代，輸出車輛最優(yōu)路徑及目標(biāo)函數(shù)F的值。

2.2 隨機(jī)場景生成

選擇適當(dāng)?shù)膱鼍吧煞椒ú⒋_定最佳場景數(shù)是求解BRPSD的首要問題。為了獲得好的場景分布，本文使用LHS來生成隨機(jī)場景。LHS是一種分層隨機(jī)抽樣，具有均勻分層的特性，可以在較少抽樣的情況下保證每個變量范圍的全覆蓋。假設(shè)有k個變量x1,x2,…,xk,要從規(guī)定的區(qū)間中取出N個樣本，首先計算每個變量的累積分布，并將其分成N個相同的小區(qū)間，隨后從每個區(qū)間中隨機(jī)選擇一個值，如圖4所示，這樣每一個變量都會獲得N個值，最后將不同變量的值進(jìn)行隨機(jī)組合，從而獲得所需樣本。在進(jìn)行抽樣時，為了簡化計算，本文假設(shè)每個抽樣的概率都是相同的。

圖4 拉丁超立方抽樣

2.3 站點重平衡需求計算

在第一階段結(jié)果的基礎(chǔ)上，第二階段求解一系列隨機(jī)場景下共享單車站點的重平衡需求。設(shè)站點j的最佳單車數(shù)量L=0.5Hj，其在場景ω下的重平衡需求計算過程如圖5所示。如果lj,ω1L，站點j需要取走單車，若此時車輛剩余容量足夠，則將多余單車移走使得站點j的單車降低到最佳數(shù)量L，否則將車輛裝滿即可。

圖5 場景ω下重平衡需求計算過程

2.4 VNS算法

Mladenovic等[17]在1997年提出VNS算法，該算法容易實現(xiàn)，通用性強(qiáng)，目前已經(jīng)被廣泛用于求解各類組合優(yōu)化問題[18]。本文采用包括3種鄰域結(jié)構(gòu)(見圖6)的VNS算法來生成新的運輸車輛路徑，每次隨機(jī)選擇其中一種鄰域結(jié)構(gòu)來執(zhí)行。(1)點插入：隨機(jī)選擇一個站點，將其從原始位置上移除，并插入到另外一個隨機(jī)位置，如圖6(a)所示。(2)節(jié)點交換：隨機(jī)選擇兩個站點，交換其位置，如圖6(b)所示。(3)2-opt：隨機(jī)選擇兩個站點，將它們之間的站點逆序，如圖6(c)所示。

圖6 鄰域結(jié)構(gòu)

3 仿真測試

使用MATLAB R2014b軟件編程實現(xiàn)LHS-VNS算法，所有試驗都是在Intel Core i5-8265U和8 GB RAM的筆記本電腦上運行。

3.1 測試算例與參數(shù)設(shè)置

為測試LHS-VNS算法的有效性，隨機(jī)生成站點規(guī)模分別為15、20、25和30的4組測試實例，每組3個算例。站點的坐標(biāo)是在緯度31.17°～31.23°和經(jīng)度121.37°～121.45°內(nèi)隨機(jī)生成的，車庫位于區(qū)域中心，其坐標(biāo)為(31.20°，121.41°)。站點容量均為100。假設(shè)客戶需求服從正態(tài)分布，所有算例中客戶需求的均值和標(biāo)準(zhǔn)差如表2所示。其他參數(shù)如下：平衡區(qū)間閾值θ=0.2，站點平衡區(qū)間均為[40,60]；運輸車輛容量為50，出庫時的裝載量為10；系數(shù)ε1和ε2分別取0.5和0.8；抽取的樣本數(shù)分別為75、100、125和150，迭代次數(shù)為10 000。

表2 客戶需求的均值與標(biāo)準(zhǔn)差

3.2 場景生成效果分析

本節(jié)計算LHS所獲得樣本的均值和標(biāo)準(zhǔn)差，以及與所設(shè)定的均值及標(biāo)準(zhǔn)差的絕對差，并與蒙特卡洛抽樣的結(jié)果進(jìn)行了對比。算例1的對比結(jié)果如表3所示，樣本量均為75。由表3可以看出，LHS所獲樣本的均值絕對差和標(biāo)準(zhǔn)差絕對差的絕大部分都優(yōu)于蒙特卡洛抽樣，且LHS樣本的均值與設(shè)定均值的絕對差值都小于0.10，標(biāo)準(zhǔn)差的絕對差值都小于0.40，而蒙特卡洛抽樣的波動較大，其均值絕對差值和標(biāo)準(zhǔn)差絕對差值最高都超過了1.00。由此說明LHS的抽樣效果優(yōu)于蒙特卡洛抽樣。

表3 LHS與蒙特卡洛抽樣對比結(jié)果

將算例6的樣本量分別設(shè)為50、75、100、125、150和175，LHS的抽樣結(jié)果如表4所示。由表4可以看出，隨著樣本量的增大，F(xiàn)的最大值趨于穩(wěn)定，所有樣本的車輛路徑均相同，但抽樣運行時間卻隨著樣本量的增加而增加。這也驗證了LHS的優(yōu)勢，可以使用較小的樣本量來生成隨機(jī)場景，在縮短計算時間的同時獲得大致相同的計算結(jié)果。

表4 樣本量分析

3.3 求解結(jié)果分析

基于樣本量為100的一系列隨機(jī)場景，每個算例均運行10次，所求得目標(biāo)F的最大值、最小值、均值和標(biāo)準(zhǔn)差如表5所示。由表5可以看出，同一規(guī)模算例的F均值和運行時間均相差不大，而規(guī)模越大的算例，其F值和運行時間越大。此外，所有算例的標(biāo)準(zhǔn)差都小于1.00，說明算法具有良好的穩(wěn)健性，其中小規(guī)模算例1和算例3的標(biāo)準(zhǔn)差為0，說明這兩個算例10次的運算結(jié)果均相同。

表5 求解結(jié)果

為了進(jìn)一步驗證本文算法求解大規(guī)模算例的性能，采用同樣的方法生成站點規(guī)模分別為50(算例13～15)和60(算例16～18)的兩組算例。設(shè)樣本量為100，每組算例均運行10次的求解結(jié)果如表6所示。由表6可以看出，所求得結(jié)果的穩(wěn)定性可以接受，但運行時間差異較大。

表6 大規(guī)模算例求解結(jié)果

4 結(jié) 語

針對需求隨機(jī)的共享單車重平衡問題，引入平衡區(qū)間代替平衡點來表示站點的平衡狀態(tài)，提出了衡量站點平衡程度的重平衡效應(yīng)度，通過LHS得到一系列具有相同概率的隨機(jī)場景來表示客戶的隨機(jī)需求，提出了BRPSD兩階段隨機(jī)規(guī)劃模型的數(shù)學(xué)表達(dá)。由于該問題是NP-hard問題，因此提出了求解該模型的LHS-VNS啟發(fā)式算法。為了驗證該算法的有效性，設(shè)計了多種規(guī)模的測試算例，測試結(jié)果表明：(1)與蒙特卡洛抽樣相比，LHS的抽樣效果更優(yōu)也更穩(wěn)定，其運行時間更短；(2)LHS-VNS可在較短時間內(nèi)獲得BRPSD的穩(wěn)定結(jié)果，具有良好的穩(wěn)健性。

未來的研究方向：一方面可以關(guān)注提高LHS-VNS的運行效率，以縮短較大規(guī)模算例的運行時間；另一方面可以研究使用多輛運輸車同時進(jìn)行重平衡操作，從而擴(kuò)展本文提出的BRPSD模型。