基于排隊網(wǎng)絡(luò)的共享單車壞車運維決策優(yōu)化

郭忠玉 羅 俊 夏 俊

(1.上海交通大學(xué) 安泰經(jīng)濟與管理學(xué)院, 上海 200030; 2.上海交通大學(xué) 中美物流研究院, 上海 200030)

0 引言

共享單車因其便捷性在日常通勤中扮演著重要的角色。根據(jù)交通運輸部2019年公布的數(shù)據(jù),截至2019年8月底,我國互聯(lián)網(wǎng)租賃自行車共有1950萬輛,覆蓋全國360個城市,注冊顧客數(shù)超過3億人次,日均訂單數(shù)達(dá)到4700萬單[1]。然而自從共享單車借助移動互聯(lián)網(wǎng)從傳統(tǒng)的有樁模式演進(jìn)成無樁模式開始,故障單車(后文簡稱“壞車”)的運維管理始終是運營者面臨的管理難題。企鵝智庫2018年發(fā)布的《共享單車數(shù)據(jù)報告:解讀摩拜ofo們的顧客與未來》中顯示,55.2%的ofo顧客認(rèn)為共享單車未能滿足其騎行需求,其中上報車輛故障的比例高達(dá)39.3%[2]。

相比傳統(tǒng)的固定車樁模式,無樁共享單車模式(后文簡稱“無樁模式”)極大地提高了顧客的使用便利,但同時也給城市治理帶來了挑戰(zhàn)。共享單車在投入使用的過程中會產(chǎn)生各種故障或損壞。單車一旦發(fā)生故障即無法有效服務(wù)顧客,需要經(jīng)過運維才能再次投入使用。當(dāng)前壞車的運維主要按照區(qū)域進(jìn)行批量回收、集中維修和批量重新投放的方式來進(jìn)行。具體地,運維人員駕駛運載工具收集所轄范圍內(nèi)的壞車,將其運送至特定的維修中心,然后裝載維修完畢的好車,重新投放至所轄區(qū)域。隨著城市治理政策的出臺和行業(yè)競爭格局日趨穩(wěn)定,當(dāng)前共享單車的運維越來越受到其運營商的重視。為了提升品牌競爭力和增加收入,運營方開始投入大量的運維資源用于壞車的回收、維修和重新投放。這些資源包括用于回收和投放的運載車輛資源,以及用于維修的服務(wù)資源等。從運維效率的角度看,維修資源的配置決策會影響壞車修復(fù)的速率,而運載資源的配置決策會影響壞車回收速率和好車投放速率。因此,兩種服務(wù)資源的配置決策會共同影響單車系統(tǒng)的運維效率。本文旨在針對共享單車系統(tǒng)的運維實踐構(gòu)建封閉排隊網(wǎng)絡(luò)和系統(tǒng)仿真模型,并基于該仿真模型優(yōu)化運維資源的配置,從而最小化顧客損失比例。

當(dāng)前共享單車的文獻(xiàn)中,考慮壞車運維的研究工作較少。Kaspi等人[3]提出一種基于貝葉斯估計的數(shù)學(xué)模型來估算有樁模式下單車的故障概率。Kaspi等人[4]指出壞車會嚴(yán)重降低顧客對共享單車系統(tǒng)的滿意度,需要提高壞車運維力度。徐國勛等人[5]以運營成本最小化為目標(biāo),針對壞車回收問題建立混合整數(shù)線性規(guī)劃模型,并使用禁忌搜索算法對較大規(guī)模的問題進(jìn)行求解。針對無樁模式的共享單車系統(tǒng),Wang和Szeto[6]在考慮壞車的情形下研究共享單車的調(diào)度策略。王涵霄等[7]提出一個比例風(fēng)險模型來實時監(jiān)控單車的使用情況,進(jìn)而制定單車的維修和調(diào)度策略。

針對大型復(fù)雜的共享單車系統(tǒng),部分文獻(xiàn)通過建立排隊網(wǎng)絡(luò)模型來研究。Li等[8-10]最先使用封閉排隊網(wǎng)絡(luò)來描述和分析有樁模式的共享單車系統(tǒng)。該排隊網(wǎng)絡(luò)中,單車看作“虛擬顧客”,停車區(qū)和道路看作“虛擬服務(wù)臺”。該研究假設(shè)道路網(wǎng)絡(luò)結(jié)構(gòu)是不可約圖,顧客到達(dá)是馬爾可夫到達(dá)過程,且單車騎行時間是獨立同分布的,并基于以上假設(shè)應(yīng)用流體和擴散極限來分析共享單車系統(tǒng)的性能表現(xiàn)。排隊網(wǎng)絡(luò)也經(jīng)常被用來研究其他服務(wù)系統(tǒng)。Adelman等人[11]使用封閉排隊網(wǎng)絡(luò)來對不同節(jié)點之間的集裝箱運輸過程進(jìn)行建模。George等人[12]針對汽車租賃系統(tǒng)構(gòu)建了封閉排隊網(wǎng)絡(luò)模型,用以確定租賃公司的最優(yōu)車隊規(guī)模,并評估公司各個租賃站汽車的可用性。

在共享單車相關(guān)研究中,部分學(xué)者使用系統(tǒng)仿真或仿真優(yōu)化方法來研究不同區(qū)域間單車的調(diào)度問題。Wang和Wang[13]運用仿真方法比較了在是否擁有實時或歷史運營數(shù)據(jù)的情形下單車的調(diào)度策略。韓琦等[14]運用系統(tǒng)仿真方法分析共享單車系統(tǒng)中車在不同區(qū)域之間的轉(zhuǎn)移狀態(tài)。Caggiani和Ottomanelli[15]基于仿真方法通過優(yōu)化路線選擇來確定最佳調(diào)度流程,從而降低共享單車運營商的車輛重置成本,提高顧客滿意度。孫怡璇等[16]結(jié)合實際調(diào)研得到的數(shù)據(jù),運用仿真方法求解共享單車的調(diào)度問題。Jian等[17]使用一種基于梯度的啟發(fā)式仿真優(yōu)化算法來優(yōu)化紐約城區(qū)的自行車和車樁分布。排序擇優(yōu)算法是仿真優(yōu)化中一類通過比較仿真實驗結(jié)果,從有限解空間中選擇出最優(yōu)或近似最優(yōu)解的方法。其中考慮無差別區(qū)域(indifference zone,當(dāng)兩個解之間的差異小于該值時,認(rèn)為二者無差異)[18]的典型的方法包括兩階段[19]、兩階段帶篩選[20]和完全序貫(fully sequential)的方法[21],以及在并行計算環(huán)境下解決大規(guī)模問題的新序貫算法[22]。排序擇優(yōu)算法可以以較少的計算資源,以給定概率選擇出最優(yōu)或接近最優(yōu)的系統(tǒng)。

近年來共享單車的相關(guān)工作主要研究系統(tǒng)的運行狀態(tài),例如單車需求量預(yù)測、顧客騎行規(guī)律,和單車調(diào)度等等。這些研究對壞車(處于損壞狀態(tài)的共享單車,后文同)的回收、維修和好車投放等運維環(huán)節(jié)涉及較少,尚有一些關(guān)鍵問題未得到解決。例如運維環(huán)節(jié)會如何影響系統(tǒng)整體表現(xiàn),如何進(jìn)行運維資源的優(yōu)化配置等問題都是共享單車運營商亟待解決的問題。本文通過對共享單車運維系統(tǒng)建立封閉排隊網(wǎng)絡(luò)模型,結(jié)合連續(xù)時間馬爾可夫鏈和離散事件系統(tǒng)仿真方法,探究了單車運營各環(huán)節(jié)對系統(tǒng)整體服務(wù)表現(xiàn)的影響。進(jìn)一步,研究總運維資源有限約束下的運維能力分配決策問題,使用排序擇優(yōu)這類仿真優(yōu)化算法對該問題進(jìn)行求解。最后在數(shù)值實驗部分,首先,我們通過小規(guī)模問題的實驗,將封閉排隊網(wǎng)絡(luò)模型的理論計算結(jié)果和離散事件系統(tǒng)仿真結(jié)果對比,同時,將理論計算結(jié)果與排序擇優(yōu)算法結(jié)果進(jìn)行對比,分別進(jìn)行交叉驗證。最后在較大規(guī)模問題上進(jìn)行了試驗,驗證了排序擇優(yōu)算法的有效性。

后文章節(jié)內(nèi)容安排如下:第一章構(gòu)建了共享單車整體系統(tǒng)的封閉排隊網(wǎng)絡(luò)模型,并基于該模型定義本文要求解的決策問題;第二章介紹基于連續(xù)時間馬爾可夫鏈(continuoustime markov chain, CTMC)和基于離散事件系統(tǒng)仿真的方法來求解系統(tǒng)性能指標(biāo),并描述了用于求解決策問題的排序擇優(yōu)算法;第三章通過數(shù)值實驗驗證前述建模及求解方法的正確性,在設(shè)定的算例上進(jìn)行系統(tǒng)性能指標(biāo)求解和決策問題的求解,并分析實驗結(jié)果;第四章對本文的研究工作進(jìn)行總結(jié)。

1 數(shù)學(xué)模型

本節(jié)結(jié)合當(dāng)前共享單車的實際運行情況和特點,建立一個封閉排隊網(wǎng)絡(luò)模型,并基于該模型提出決策問題。其中,1.1節(jié)介紹了模型所需的假設(shè),以及相關(guān)變量和參數(shù)。1.2節(jié)對模型進(jìn)行詳細(xì)描述。1.3節(jié)描述了針對維修和運載資源優(yōu)化配置的決策問題。

1.1 模型假設(shè)和符號標(biāo)記

記A為區(qū)域數(shù)量,i為區(qū)域的腳標(biāo);M為單車總數(shù)。

1) 共享單車系統(tǒng)劃分成A個不同的區(qū)域;

2) 區(qū)域i的顧客到達(dá)服從參數(shù)λi的泊松過程,各區(qū)域的顧客到達(dá)互相獨立;

3) 記P為顧客在各個區(qū)域間騎行的概率轉(zhuǎn)移矩陣,顧客到達(dá)區(qū)域i,如果有好車,則以轉(zhuǎn)移概率Pij騎往區(qū)域j,否則該顧客立即離開,i記為起始區(qū)域,j記為目標(biāo)區(qū)域。若起始與目標(biāo)區(qū)域相同,則表示顧客在區(qū)域內(nèi)騎行;

4) 顧客從區(qū)域i至區(qū)域j的單車騎行時間服從參數(shù)ρij的指數(shù)分布,各騎行過程相互獨立;

5) 顧客騎行完畢后,以1-β的概率進(jìn)入目標(biāo)區(qū)域,否則進(jìn)入損壞狀態(tài);

6) 進(jìn)入損壞狀態(tài)的單車通過回收、維修、再投放這三個過程,重新服務(wù)顧客;

7) 運載車間隔參數(shù)為τ的指數(shù)分布時間,對運載車的最大載荷C以內(nèi)盡可能多的壞車進(jìn)行回收作業(yè);然后間隔相同參數(shù)τ的指數(shù)分布時間,對載荷C以內(nèi)盡可能多的維修完畢的好車進(jìn)行投放作業(yè);

8) 同一輛運載車持續(xù)不斷地交替進(jìn)行回收和投放作業(yè);

9) 該特定維修中心有N個維修工,維修工修復(fù)一輛單車的服務(wù)時間服從獨立同分布且參數(shù)μ的指數(shù)分布。

表1為系統(tǒng)狀態(tài)變量的各分量描述。記t為系統(tǒng)時間變量。

表1 狀態(tài)隨機變量分量表Table 1 Random variables of system states

1.2 封閉排隊網(wǎng)絡(luò)模型

共享單車系統(tǒng)的運行可以解構(gòu)為兩大部分:(1)正常服務(wù)狀態(tài)的部分;(2)因故障進(jìn)入運維狀態(tài)的部分。本文中,我們將損壞單車的回收、維修、投放統(tǒng)稱為共享單車的運維。

共享單車日常會根據(jù)供需情況對區(qū)域間的單車進(jìn)行平衡調(diào)度,比如在需求低谷期將單車從過剩區(qū)域轉(zhuǎn)移至不足區(qū)域。共享單車調(diào)度策略是一個復(fù)雜的研究課題,目前已有相關(guān)研究工作[5,13,23]。本文聚焦于共享單車系統(tǒng)的運維。由于單車的重新投放過程一定程度上發(fā)揮了再平衡的作用,因此我們不具體考慮平衡調(diào)度對系統(tǒng)的影響。

圖1展示了由兩個區(qū)域組成的共享單車系統(tǒng),系統(tǒng)各過程中等待或正在被服務(wù)的共享單車構(gòu)成隊列。因為同一隊列中的單車沒有區(qū)別,所以不需要考慮排隊的順序。某區(qū)域i的顧客到達(dá)服從參數(shù)λi的泊松過程,并且各區(qū)域間的顧客到達(dá)相互獨立。顧客到達(dá)某區(qū)域i時,選擇一輛單車,并根據(jù)概率轉(zhuǎn)移矩陣P中第i行,確定目標(biāo)區(qū)域j,進(jìn)入騎行隊列。各個騎行過程相互獨立,因此相當(dāng)于進(jìn)入騎行隊列的單車都能夠馬上得到服務(wù)。騎行的時間服從參數(shù)為ρij的指數(shù)分布。當(dāng)顧客到達(dá)而該區(qū)域沒有可騎單車時,該顧客立即流失。騎行完畢后,單車以概率1-β排至目標(biāo)區(qū)域j的單車隊列,否則以概率β進(jìn)入一個回收池BrokenPool(簡稱BP),即單車進(jìn)入待回收隊列。損壞的單車需要經(jīng)過運載車回收,運維中心的維修工維修和運載車的投放來重新進(jìn)入各個區(qū)域的單車隊列,并服務(wù)顧客。運載車有回收和投放兩種工作狀態(tài)。處于搬運狀態(tài)時,其經(jīng)過服從參數(shù)為τ的指數(shù)分布時間后,將不多于載荷C的盡可能多的壞車運至維修中心。運送至維修中心的壞車進(jìn)入等待維修隊列。運載車隨即進(jìn)入投放狀態(tài)。運載車再經(jīng)過參數(shù)為τ的指數(shù)分布時間,將不多于C的盡可能多的維修后的好車重新投放至各區(qū)域的單車隊列。運載車隨即進(jìn)入回收狀態(tài)?；厥蘸屯斗哦呓惶孢M(jìn)行。這里以DistributePool(簡記為DP),表示已維修待投放的單車隊列。維修中心有N個獨立、并行的維修工,每個維修工每次維修一輛車,每個維修工的維修時間服從相同參數(shù)μ的指數(shù)分布。維修完畢的單車進(jìn)入待投放隊列。系統(tǒng)中單車的總量不變,并在各個隊列中轉(zhuǎn)移,于是構(gòu)成一個封閉排隊網(wǎng)絡(luò)。

圖1 區(qū)域數(shù)為2時單車封閉排隊系統(tǒng)示意圖Figure 1 Diagram of closed queueing network consisting of two areas

本文采用一種固定的投放策略,具體步驟為:(1)對某區(qū)域,以其顧客到達(dá)率除以各個區(qū)域到達(dá)率之和,再乘以單車總數(shù)并向下取整,得到初始各區(qū)域的目標(biāo)數(shù)量;由于取整,可能會出現(xiàn)初始目標(biāo)數(shù)量未達(dá)到單車總量,進(jìn)一步以單車總量減去各區(qū)域初始目標(biāo)數(shù)量之和得到剩余數(shù)量;按到達(dá)率由大至小,給前剩余數(shù)量個區(qū)域的目標(biāo)數(shù)量加一,以此作為每個區(qū)域的目標(biāo)投放數(shù)量。(2)投放時,如果某區(qū)域的好車數(shù)量已達(dá)目標(biāo)量,則不向該區(qū)域投放;否則以當(dāng)前運載車上的好車數(shù)為上限,盡可能多地投放至該目標(biāo)量。

由于該排隊網(wǎng)絡(luò)中各過程的指數(shù)時間性質(zhì),系統(tǒng)中相鄰事件發(fā)生的間隔時間服從指數(shù)分布,因此該系統(tǒng)的狀態(tài)變化構(gòu)成一個CTMC。由于事件發(fā)生時間的隨機性,以及單車狀態(tài)的隨機性,系統(tǒng)中任意狀態(tài)之間可以以一定概率和一定時間進(jìn)行轉(zhuǎn)移。因此,這是一個遍歷(ergodic)的CTMC,系統(tǒng)存在長程穩(wěn)定狀態(tài)。由此,我們在給定系統(tǒng)參數(shù)設(shè)置下,可以通過系統(tǒng)狀態(tài)概率轉(zhuǎn)移方程,求解系統(tǒng)處于各個狀態(tài)的概率,從而得到系統(tǒng)性能指標(biāo)。

1.3 決策模型

壞車運維需要配置合適數(shù)量的運維資源,如運載車、維修工等。通常運載車和維修工的單位成本是已知的,運營方有固定的預(yù)算來執(zhí)行運維作業(yè)。在總預(yù)算限制下,運營方需要優(yōu)化運載車和維修工的配置數(shù)量,從而最小化系統(tǒng)中期望顧客損失比例。模型所需的參數(shù)和變量如表2所示。

表2 決策問題相關(guān)參數(shù)及變量表Table 2 Parameters and variables of the decision problem

數(shù)學(xué)模型如下所示:

模型目標(biāo)函數(shù)中的顧客損失比例為各個區(qū)域因沒有可供騎行單車而流失的顧客數(shù)量總和,占所有到達(dá)系統(tǒng)的顧客數(shù)量的比例。

2 模型求解

本節(jié)描述封閉排隊網(wǎng)絡(luò)模型性能指標(biāo)的求解方法和決策模型的求解算法。2.1節(jié)介紹兩種性能指標(biāo)的求解方法。其中2.1.1節(jié)介紹基于CTMC長程穩(wěn)態(tài)概率的系統(tǒng)性能求解方法。2.1.2節(jié)介紹基于離散事件系統(tǒng)仿真的求解方法。2.2節(jié)介紹用于求解決策模型的排序擇優(yōu)算法。

2.1 系統(tǒng)性能指標(biāo)求解方法

2.1.1 CTMC狀態(tài)方程求解方法

對于遍歷的CTMC,當(dāng)系統(tǒng)時間趨向于無窮大時,系統(tǒng)處于各個狀態(tài)的時間比例收斂于其極限概率[24]。因此,我們可以由系統(tǒng)狀態(tài)轉(zhuǎn)移的偏微分方程構(gòu)建對任意系統(tǒng)而言狀態(tài)轉(zhuǎn)出率與轉(zhuǎn)入率相等的方程組。求解該狀態(tài)轉(zhuǎn)移方程組可得系統(tǒng)處于各狀態(tài)的概率,即系統(tǒng)處于每個狀態(tài)的時間比例?；谠摃r間比例可計算得到系統(tǒng)性能指標(biāo)。以上基于CTMC的狀態(tài)微分方程組求解結(jié)果理論上是精確解。

記S為系統(tǒng)所有可能的狀態(tài)的集合,s為集合中元素。記s為系統(tǒng)狀態(tài)向量如下:

因為系統(tǒng)中的車的總數(shù)固定,所以狀態(tài)變量符合以下等式:

記πs為系統(tǒng)處于狀態(tài)s的概率。記INi如下:

對任意系統(tǒng)狀態(tài),其轉(zhuǎn)出速率為該狀態(tài)可以到達(dá)的任意下一狀態(tài)的長程比例乘以相應(yīng)轉(zhuǎn)移速率之和。即:

由于系統(tǒng)狀態(tài)的繁雜,為方便實際計算,我們使用以下示性函數(shù)表示一個狀態(tài)是否可能出現(xiàn)在當(dāng)前系統(tǒng)狀態(tài):

穩(wěn)態(tài)情況下,對任意狀態(tài)其轉(zhuǎn)入速率為任意可以到達(dá)該狀態(tài)的前繼狀態(tài)的長程比例乘以相應(yīng)轉(zhuǎn)移速率之和。具體為,任意系統(tǒng)狀態(tài)可由另一個狀態(tài)經(jīng)由(一)顧客騎走一輛好車;(二)騎行到達(dá);(三)單車損壞;(四)壞車回收至維修中心;(五)一輛單車維修完畢;(六)投放至各區(qū)域;這六種過程的某一種而進(jìn)入。

其中,過程(五)的計算需要考慮到當(dāng)前狀態(tài)中C,BP和RC的值。從前一狀態(tài)進(jìn)入當(dāng)前狀態(tài)所搬運的數(shù)量不超過C和RC之間的較小值。同時,當(dāng)該搬運數(shù)量小于C時,BP必為0,否則由搬運的規(guī)則可知,搬運的數(shù)量應(yīng)當(dāng)更多。于是,我們使用以下示性函數(shù)判斷一個狀態(tài)是否可能為通過搬運進(jìn)入當(dāng)前狀態(tài)的系統(tǒng)狀態(tài):

此外,過程(六)需要根據(jù)1.3節(jié)中描述的投放策略得到所有可能的投放前系統(tǒng)狀態(tài)集合Q。前述六種過程分別對應(yīng)表達(dá)式(7)中的1-6行,從而確定系統(tǒng)任意狀態(tài)的轉(zhuǎn)入速率之和為:

因為系統(tǒng)必須處在某個狀態(tài),因此必須滿足:

上述(2)～(8)式構(gòu)成了系統(tǒng)狀態(tài)轉(zhuǎn)移概率方程組。由于遍歷的CTMC極限穩(wěn)態(tài)概率的唯一性,該方程組具有唯一解[25]。通過求解上述方程組,我們可求解得到系統(tǒng)處于每個狀態(tài)的長程穩(wěn)態(tài)概率,也即長程比例。根據(jù)系統(tǒng)的每個狀態(tài)的長程穩(wěn)態(tài)概率,可以計算得到一些系統(tǒng)性能指標(biāo)。本文所選取的三個系統(tǒng)性能指標(biāo)主要是為了反映共享單車系統(tǒng)的運行性能狀況,與實際中單車企業(yè)的收入、成本和運行效率緊密相關(guān)。其中,期望顧客損失比例表示到達(dá)系統(tǒng)而未被服務(wù)的顧客數(shù)占到達(dá)系統(tǒng)顧客總數(shù)的比例的期望,與用戶使用體驗直接相關(guān)。高的顧客損失比例對應(yīng)著差的用戶體驗,同時也意味著用戶的減少,進(jìn)而帶來單車企業(yè)收入的損失,是一個同時關(guān)系到共享單車供需雙方的指標(biāo);期望好車比例為系統(tǒng)中好車數(shù)占總車數(shù)的比例的期望,與共享單車企業(yè)的單車成本投入相關(guān)聯(lián),好車比例越高所需要的維修和替換投入越少;期望維修工空閑比例則為空閑維修工數(shù)占總維修工數(shù)比例的期望,反映了單車企業(yè)進(jìn)行維修的人力成本投入的利用效率。根據(jù)解決問題的需要,我們還可以考慮其他類似的性能指標(biāo),例如各區(qū)域的平均可使用的單車數(shù)。為了簡便起見,本文中僅舉出三個代表性的指標(biāo),并以期望顧客損失比例作為決策目標(biāo)進(jìn)行詳細(xì)論證(類似地,我們都可以用同樣的方法得出其他性能指標(biāo)的結(jié)論)。

(1)期望顧客損失比例

當(dāng)系統(tǒng)中某個區(qū)域沒有好車時,新到達(dá)的顧客會因騎行需求不能立即得到滿足而損失。因此,我們用沒有好車的區(qū)域的顧客到達(dá)速率除以總的顧客到達(dá)速率,來計算任意狀態(tài)的顧客損失比例。然后以所有狀態(tài)比例為權(quán)重,計算顧客損失率的期望值:

(2)期望好車比例

系統(tǒng)中好車和壞車的比例可以由任意狀態(tài)中的好車數(shù)量占車輛總數(shù)的比例,以系統(tǒng)狀態(tài)比例為權(quán)重加權(quán)平均得到:

(3)期望維修工空閑比例

維修中心處可能會有多個維修工。當(dāng)維修中心處待維修的單車數(shù)少于維修工數(shù)量時,即出現(xiàn)維修工空閑。因此,維修中心的車數(shù)除以維修工數(shù)量得到某一狀態(tài)的空閑比例,再以所有狀態(tài)比例為權(quán)重加權(quán)平均,即可計算期望維修工空閑比例。

2.1.2 離散事件系統(tǒng)仿真

對本文所建立的封閉排隊網(wǎng)絡(luò)模型所對應(yīng)的CTMC而言,其狀態(tài)數(shù)量可通過系統(tǒng)中的區(qū)域數(shù)量、單車總數(shù)以及運載車的數(shù)量計算得到,總狀態(tài)數(shù)量為(此處!為階乘符號)。

由此可知,系統(tǒng)中的狀態(tài)數(shù)量與區(qū)域數(shù)量、單車總數(shù)和運載車數(shù)量呈指數(shù)關(guān)系,也就是說系統(tǒng)對應(yīng)的系統(tǒng)狀態(tài)轉(zhuǎn)移概率方程組將隨共享單車系統(tǒng)的規(guī)模增加而迅速增多。由于時間和計算機內(nèi)存的限制,2.1.1節(jié)介紹的方程求解方法難以解決系統(tǒng)狀態(tài)數(shù)較多的情況。相比之下,基于仿真的方法可以無需遍歷所有狀態(tài)來獲得系統(tǒng)性能指標(biāo)期望值的估計值。因此,我們提出基于離散事件的系統(tǒng)仿真方法來計算系統(tǒng)性能指標(biāo)。

我們基于前述構(gòu)建的封閉排隊網(wǎng)絡(luò)模型,使用計算機程序編寫離散事件仿真系統(tǒng),再生成相應(yīng)分布的隨機數(shù)來運行離散事件的仿真。在仿真過程中,我們對系統(tǒng)狀態(tài)進(jìn)行記錄或計數(shù)來輔助計算系統(tǒng)性能指標(biāo)。通過記錄系統(tǒng)中發(fā)生的事件,可以計算系統(tǒng)性能指標(biāo)的均值。隨著系統(tǒng)仿真的進(jìn)行,系統(tǒng)逐漸進(jìn)入穩(wěn)態(tài),由大數(shù)定律,仿真得到的性能指標(biāo)均值也逐漸收斂到系統(tǒng)性能指標(biāo)的理論期望值。具體地,顧客損失比例可通過到達(dá)顧客總數(shù)和未被服務(wù)而丟失的顧客總數(shù)相除得到;好車比例和維修工空閑比例可通過記錄系統(tǒng)狀態(tài)和狀態(tài)持續(xù)時間來計算。

2.2 決策問題求解-排序擇優(yōu)算法

決策模型(1)的解空間有限、離散。由于系統(tǒng)狀態(tài)數(shù)量較多,系統(tǒng)性能指標(biāo)難以通過方程直接求解得到;此外,每一組特定參數(shù)設(shè)置下的仿真結(jié)果具有隨機性,因而無法進(jìn)行不同參數(shù)設(shè)置之間的比較,除非運行大量仿真實驗。Kim和Nelson在2001年提出的KN算法(KN procedure)[21],即可在正確率確保選出最優(yōu)參數(shù)設(shè)置的同時,通過序列式篩選,節(jié)省計算開支,適用于求解本文的決策問題。Hong等人[26]的研究表明,該算法的計算時間復(fù)雜度與待比較的方案個數(shù)K的關(guān)系為O(KlogK)。

具體而言,本文所采用的KN排序擇優(yōu)算法可保證以設(shè)定的1-α概率水平從候選系統(tǒng)中挑選出最優(yōu)的系統(tǒng)。參數(shù)δ表示無差別區(qū)間長度,當(dāng)某個待選系統(tǒng)在性能表現(xiàn)與最優(yōu)的系統(tǒng)差距在δ以內(nèi)時,認(rèn)為兩個系統(tǒng)的性能表現(xiàn)無差別。表現(xiàn)在結(jié)果上是當(dāng)某些系統(tǒng)的均值與最優(yōu)系統(tǒng)對應(yīng)的均值之間的距離小于該值時,多次執(zhí)行算法,最優(yōu)解無差別范圍內(nèi)的系統(tǒng)也將被當(dāng)做最優(yōu)系統(tǒng)而以一定頻率被選出。這種解記作δ-最優(yōu)解。該方法的優(yōu)點在于:(一)可以設(shè)置無差別區(qū)間,從而可以認(rèn)為決定什么樣程度的性能差別是沒有影響的;(二)其篩選過程中,可以逐步剔除明顯差的方案,從而節(jié)省比較所需要的總樣本數(shù);(三)該方法雖然基于樣本分布為正態(tài)的假設(shè),但是對樣本分布非正態(tài)同樣有著良好的適用性,可以應(yīng)對在實際使用過程中,樣本分布未知的情形[21]。算法的輸入?yún)?shù)如表3所示。

表3 KN算法參數(shù)表Table 3 Parameters of KN procedure

具體到求解本文的運維決策問題中,K表示待比較的運維方案的個數(shù);1-α為人為設(shè)定的,實際上的最優(yōu)運維方案被算法選出的概率大小;δ為人為設(shè)定的判斷閾值,當(dāng)兩個運維方案的性能指標(biāo)差異小于該值時,認(rèn)為兩個方案的性能指標(biāo)沒有差異;n0為算法初始比較時,每個運維方案的仿真樣本個數(shù);c為算法參數(shù),與實際運維問題無關(guān);k為運維方案的標(biāo)號。

記Xkp為系統(tǒng)k的第p個觀測值。

算法流程如算法1所示。

算法1 KN算法輸入:K,α,δ,n0,c輸出:最優(yōu)系統(tǒng)編號1 初始化:2 I={1,2,…,K}3 η=1 2[( )-2/(n0-1)-1 2α K-1]4 h2=2(n0-1)cη 5 S2kl=1 n0-1 n0 l(n0)])2,?k∈Iandk≠l 6 Nkl=∑p=1(Xkp-Xlp-[X k(n0)-Ximages/BZ_223_1480_1843_1503_1942.pngh2S2k l δ2images/BZ_223_1586_1843_1608_1942.png,?k∈I and k≠l 7 Nk=max l≠k Nkl,?k∈I 8 Xk(r)=1 r ∑r p=1Xkp 9 初始判斷:10 if n0＞maxkNk then 11 return argmaxkXk(n0)12 else 13 r=n0 14 進(jìn)入篩選過程15 end if 16 篩選過程:17 while |I|＞1: 18 ?系統(tǒng)k∈I中產(chǎn)生一個觀測值Xk,r+1 19 r=r+1 20 if r==maxkNk+1 then 21 return argmaxkX k(r)22 else 23 Iold=I 24 Xk(r)=1 r ∑r p=1Xkp,?k∈I 25 Wkl(r)=max 0,δ 2cr h2S2 kl { },?k∈Iandk≠l ( )δ2-r

26 I={k:k∈Iold and X k(r)≥X l(r)-Wkl(r),?l∈Iold,l≠k}27 end if 28 end while 30 return I

系統(tǒng)k與l觀測值之差的樣本方差,Nk為系統(tǒng)k的最大樣本量。在將該排序擇優(yōu)算法應(yīng)用到本文的決策問題時,注意到算法原適用于最大化問題,因此在本研究實際應(yīng)用該算法時,將目標(biāo)函數(shù)做取反處理。

3 數(shù)值實驗

本節(jié)我們驗證第2章中封閉排隊網(wǎng)絡(luò)建模和離散事件系統(tǒng)仿真方法的正確性并設(shè)置算例測試排序擇優(yōu)算法的效果。3.1節(jié)介紹系統(tǒng)性能指標(biāo)求解的過程,在設(shè)定的算例上驗證其正確性,并分析結(jié)果。3.2節(jié)基于小規(guī)模問題通過比較方程計算得到的結(jié)果、仿真計算得到的期望值和排序擇優(yōu)算法的結(jié)果,驗證本文所采用算法的正確性。然后在不同算例上運行排序擇優(yōu)算法,驗證其優(yōu)化效果。

3.1 運維模型數(shù)值實驗

3.1.1 求解方法設(shè)置

在如表4所示的默認(rèn)系統(tǒng)參數(shù)設(shè)置下,使用方程求解和系統(tǒng)仿真兩種方法來計算系統(tǒng)性能指標(biāo)。

表4 系統(tǒng)仿真默認(rèn)參數(shù)表Table 4 Default parameters of simulation

使用系統(tǒng)狀態(tài)轉(zhuǎn)移方程求解主要經(jīng)過三大步驟:(一)根據(jù)系統(tǒng)參數(shù)設(shè)置,生成所有可能的系統(tǒng)狀態(tài);(二)使用稀疏矩陣存儲方程組;(三)求解方程組。本研究中,我們基于Python調(diào)用Scipy包的稀疏矩陣方程求解算法來對方程進(jìn)行求解。

通過對仿真模型的簡單擴展,可以仿真系統(tǒng)中包含多輛運載車同時運行的情況。根據(jù)仿真結(jié)果,我們對結(jié)果進(jìn)行排序,從而挑選出最優(yōu)解或近似最優(yōu)解,也即δ-最優(yōu)解,繼而求解目標(biāo)決策問題。由于初始系統(tǒng)狀態(tài)對長程穩(wěn)定狀態(tài)下的系統(tǒng)性能指標(biāo)沒有影響,我們統(tǒng)一設(shè)置初始狀態(tài)下系統(tǒng)中所有車均為好車,且均勻分布于各個區(qū)域中。仿真運行過程中,觀察系統(tǒng)指標(biāo)變化可判斷系統(tǒng)是否進(jìn)入穩(wěn)態(tài)。取系統(tǒng)進(jìn)入穩(wěn)態(tài)后的一段時間,計算這段時間某性能指標(biāo)均值,作為該指標(biāo)的期望值的估計[25]。為了快速得到CTMC模型方程的解,我們首先考慮一個小規(guī)模問題。特別要指出的是,在現(xiàn)實中不同場景下會有不同的顧客騎行轉(zhuǎn)移概率和到達(dá)速率。本文中不同的Pij、λi相當(dāng)于現(xiàn)實中不同的地區(qū)或時間等不同的單車系統(tǒng)的運行場景。在實際中可以通過實際運行情況進(jìn)行統(tǒng)計、估計。在此,為了研究便利,我們采用隨機數(shù)。3.1節(jié)中取

并且,在后文的數(shù)值實驗中生成不同的隨機數(shù),以比較不同的場景,反映計算方法對不同情景的魯棒性。為了文章簡潔,后文不再具體給出取值。仿真結(jié)果如圖2所示,縱軸為比例值,橫軸為離散系統(tǒng)仿真的系統(tǒng)時間,三條曲線分別代表了性能指標(biāo)仿真進(jìn)行至當(dāng)前時間的平均值。如圖所示,2.1.1節(jié)中介紹的三種性能指標(biāo)隨仿真進(jìn)程推進(jìn),其均值逐漸趨于穩(wěn)定。如2.1.2節(jié)所述,后續(xù)將根據(jù)穩(wěn)定值估計系統(tǒng)性能指標(biāo)。

圖2 A=2,M=6時系統(tǒng)性能指標(biāo)-時間30000圖Figure 2 Performance indicators-simulation time graph with A=2, M=6

3.1.2 系統(tǒng)狀態(tài)參數(shù)數(shù)值算例

由于輸入?yún)?shù)的隨機性,離散事件系統(tǒng)仿真結(jié)果存在波動。為準(zhǔn)確得到結(jié)果,我們以20次仿真結(jié)果的均值作為一個仿真得到的系統(tǒng)性能指標(biāo)結(jié)果。基于表4中的默認(rèn)參數(shù)設(shè)置,我們在不同的系統(tǒng)運維速率下比較方程求解和系統(tǒng)仿真所得的結(jié)果。如下圖3所示。

圖3 A=2,M=6時方程與仿真結(jié)果對比圖-變動運載速率Figure 3 Comparison of ODE and simulation result varying speed of carrier with A=2, M=6

圖4 A=2,M=6時方程與仿真結(jié)果對比圖-變動維修速率Figure 4 Comparison of ODE and simulation result varying repairing speed with A=2, M=6

通過對比方程求解和仿真結(jié)果可以看到,仿真方法可以逼近系統(tǒng)性能指標(biāo)的CTMC方程計算值。此外,單獨提升運維系統(tǒng)中的運載車或維修工的速率,可以降低顧客的損失比例,提升系統(tǒng)中好車的比例,并提高維修工空閑的比例。但從圖中也可以看到隨著速率的進(jìn)一步提升,系統(tǒng)性能指標(biāo)的曲線變得平緩,也即單獨提升運維系統(tǒng)中的運載車或維修工的速率,提升改善效果存在邊際遞減效應(yīng)。

圖5展示了不同運載車數(shù)量和維修工數(shù)量下顧客損失比例的變化。從圖5中可知,增加維修工或運載車的數(shù)量均可降低顧客的損失率,改善系統(tǒng)性能表現(xiàn)。但當(dāng)數(shù)量持續(xù)增加時,曲線的下降趨勢漸趨平緩,說明改善效果邊際遞減。當(dāng)二者增加到一定程度后對系統(tǒng)的性能改善效果幾乎消失。由此可知,對一個共享單車系統(tǒng)而言,單純通過增加維修或投放的能力對改善系統(tǒng)的性能表現(xiàn)是存在上限的。因此在相同資源限制條件下,可能存在性能表現(xiàn)相近,而資源投入存在差異的方案。實驗結(jié)果符合常識推論。

圖5 A=2,M=6時變動維修工及運載車數(shù)量仿真結(jié)果Figure 5 Simulation result with varying number of repairing workers and carriers server with A=2, M=6

3.2 決策模型數(shù)值實驗

3.2.1 仿真優(yōu)化算法驗證

我們首先在A=2,M=6的情況下使用CTMC方程方法求解CostR=1,CostC=1,CostTotal=5,nC=1的小規(guī)模決策問題中的每一個解的期望顧客損失比例。

通過方程計算得到如表5所示結(jié)果,顯示(N,nc)=(4,1)為最優(yōu)解。在δ=0.01時,(N,nc)=(2,1)和(N,nc)=(3,1)在最優(yōu)解的無差別范圍之內(nèi),為δ-最優(yōu)解。

為了排序擇優(yōu)算法的適用性,仿真結(jié)果應(yīng)當(dāng)盡可能符合正態(tài)分布。圖6為不同決策問題解的200次仿真結(jié)果的頻數(shù)分布分布統(tǒng)計圖。4個解對應(yīng)仿真結(jié)果的分布均能通過正態(tài)性的KS檢驗(Kolmogorov-Smirnov test)。按照表3中的參數(shù)設(shè)置,運行1000次排序擇優(yōu)算法,結(jié)果(N,nc)=(4,1) 被選中974次,(N,nc)=(3,1)選中26次。所有結(jié)果中,在最優(yōu)解無差別區(qū)域內(nèi)的解被選中的比例超過了設(shè)定的95%。該實驗驗證了算法的正確性。

圖6 小規(guī)模決策問題解對應(yīng)仿真結(jié)果頻數(shù)分布統(tǒng)計圖Figure 6 Frequency distribution graph for the simulation result of the solutions of the small-scale decision problem

3.2.2 決策問題數(shù)值算例

分別在A=2,M=6,A=10,M=100和A=10,M=1000的情況下,取CostR=1,CostC=1,CostTotal=10;不考慮N=0或nC=0的情況;相應(yīng)的KN算法的參數(shù)K=45。以不同的轉(zhuǎn)移概率和各區(qū)域的顧客到達(dá)率作為不同場景,在各場景中求解其對應(yīng)的決策問題。在如表4所給出的參數(shù)設(shè)置下,對前述場景分別運行重復(fù)試驗。結(jié)合3.1.2的結(jié)果,和從表6和表7中看到的50次仿真結(jié)果的置信區(qū)間,我們可以認(rèn)為50次仿真結(jié)果的均值可以較準(zhǔn)確地估計實際值。以每個場景中的每個解的50次仿真結(jié)果的均值和標(biāo)準(zhǔn)差作為比較標(biāo)準(zhǔn),統(tǒng)計算法結(jié)果。其中正確選擇比例(PCS,probability of correct selection)為KN算法選出最優(yōu)解或δ-最優(yōu)解的次數(shù)占總實驗次數(shù)的比例。展示KN算法選出的結(jié)果如表6、表7和表8所示。

表6 決策問題數(shù)值算例,A=2,M=6,500次重復(fù)Table 6 Numerical experiments of decision problem, A=2, M=6, 500 times repetition

表7 決策問題數(shù)值算例,A=10,M100,500次重復(fù)Table 7 Numerical experiments of decision problem, A=10, M=100, 500 times repetition

表8 決策問題數(shù)值算例,A=10,M=1000,10次重復(fù)Table 8 Numerical experiments of decision problem, A=10, M=1000, 10 times repetition

三張表中的結(jié)果顯示,排序擇優(yōu)算法在不同場景下,可以給定概率選出性能指標(biāo)與最優(yōu)解在給定的無差別區(qū)域內(nèi)的方案,試驗顯示較少次數(shù)的排序擇優(yōu)算法重復(fù)試驗,也已可以選出較好的方案。從而保證了問題求解的最優(yōu)性,并具有一定的魯棒性。該方法適用于在難取樣的復(fù)雜系統(tǒng)中挑選可行方案。此外,排序擇優(yōu)算法的運行是通過逐步剔除明顯差解。本研究中,運行KN算法需要900次左右仿真得到一個結(jié)果。而通常情況下,基于仿真樣本均值估計的求解方式往往難以準(zhǔn)確估計所需的仿真次數(shù)。因此,在算法效率方面,KN算法對于解決本文的管理決策問題具有一定的優(yōu)勢。

4 結(jié)論

本文結(jié)合共享單車的實際運維特點,將無樁共享單車的整體運行過程建模成封閉排隊網(wǎng)絡(luò),并構(gòu)建對應(yīng)的CTMC模型。該模型可通過狀態(tài)轉(zhuǎn)移方程求得系統(tǒng)長程穩(wěn)態(tài)概率和系統(tǒng)性能指標(biāo),具有理論上的精確性。在大規(guī)模共享單車系統(tǒng)運行環(huán)境下,系統(tǒng)中的狀態(tài)數(shù)量呈指數(shù)速度增長,導(dǎo)致狀態(tài)轉(zhuǎn)移方程難以求解。并且該計算方法,依賴于指數(shù)分布假設(shè),難以適用于實際問題。本文提出基于離散事件系統(tǒng)仿真的方法來估計系統(tǒng)性能指標(biāo)。該方法適用于不同的系統(tǒng)假設(shè),例如不同的顧客到達(dá)分布假設(shè)等,極大增加了研究的可擴展性,為研究共享單車系統(tǒng)運行狀態(tài)提供參考。通過觀察單車系統(tǒng)參數(shù)變化對性能指標(biāo)的影響,我們發(fā)現(xiàn)增加系統(tǒng)中運載車的速率和數(shù)量,或提升維修環(huán)節(jié)速度、增加維修工數(shù)量可提高系統(tǒng)性能表現(xiàn),但對系統(tǒng)的改善作用都存在邊際遞減現(xiàn)象。因此,需要解決將總數(shù)有限的運維資源在運載和維修環(huán)節(jié)之間進(jìn)行合理分配的決策問題。排序擇優(yōu)算法適用于求解該決策問題,并且可以從理論上保證,在給定的正確概率水平下,從不同性能指標(biāo)接近的方案中挑選得到可接受的最優(yōu)方案或接近最優(yōu)的δ-最優(yōu)方案。本文通過小規(guī)模數(shù)值算例對比理論計算和仿真優(yōu)化算法結(jié)果,交叉驗證了算法的正確性。并進(jìn)一步在設(shè)置的算例上展示算法運行效果。為共享單車實際運營中進(jìn)行運維能力的分配決策提供解決方案參考。