賈永基， 惲博文， 許媛媛

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

隨著我國經(jīng)濟(jì)的迅猛發(fā)展，交通運(yùn)輸需求快速增加，未來仍將長(zhǎng)期呈現(xiàn)增長(zhǎng)態(tài)勢(shì)，這導(dǎo)致交通運(yùn)輸系統(tǒng)的碳排放所占份額越來越大。在2021年全國兩會(huì)上，“碳達(dá)峰”與“碳中和”首次被寫入政府工作報(bào)告。為了如期實(shí)現(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ùn)營成本的同時(shí)，降低用戶滿意度，這將極大阻礙共享單車行業(yè)的健康發(fā)展。因而實(shí)施有效的共享單車供需重新匹配操作，即共享單車重平衡問題(bike-sharing rebalance problem，BRP)[2]，顯得尤為重要。

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

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

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

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

1.1 問題描述

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

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

本文所使用的決策變量和符號(hào)如表1所示。

表1 決策變量和符號(hào)

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

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

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

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

圖2 站點(diǎn)的重平衡效用度

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

基于LHS的兩階段隨機(jī)規(guī)劃模型的第一階段為運(yùn)輸車輛路徑優(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)表示每個(gè)站點(diǎn)只能被訪問一次；式(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ī)場(chǎng)景下的重平衡效用度期望；式(7)表示場(chǎng)景ω下運(yùn)輸車輛的裝載量不大于車輛容量；式(8)表示場(chǎng)景ω下車輛服務(wù)站點(diǎn)前后的裝載量之差等于站點(diǎn)被滿足的單車數(shù)量；式(9)表示場(chǎng)景ω下站點(diǎn)被服務(wù)前后的單車數(shù)量之差等于站點(diǎn)被滿足的單車數(shù)量；式(10)表示場(chǎng)景ω下站點(diǎn)不能同時(shí)具有取車和卸車需求；式(11)表示場(chǎng)景ω下的站點(diǎn)重平衡效用度；式(12)表示決策變量sj,ω+不能超過車輛訪問站點(diǎn)時(shí)的剩余容量；式(13)表示決策變量sj,ω-不能超過車輛訪問站點(diǎn)時(shí)的裝載量；式(14)表示決策變量的取值范圍。

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

2.1 算法流程

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

圖3 LHS-VNS算法流程圖

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

2.2 隨機(jī)場(chǎng)景生成

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

圖4 拉丁超立方抽樣

2.3 站點(diǎn)重平衡需求計(jì)算

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

圖5 場(chǎng)景ω下重平衡需求計(jì)算過程

2.4 VNS算法

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

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

3 仿真測(cè)試

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

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

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

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

3.2 場(chǎng)景生成效果分析

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

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

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

表4 樣本量分析

3.3 求解結(jié)果分析

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

表5 求解結(jié)果

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

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

4 結(jié) 語

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

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