劉喜梅,潘立軍
湖南工程學(xué)院 管理學(xué)院,湖南 湘潭411104
公共自行車(chē)再平衡問(wèn)題(Bike Sharing Rebalancing Problem,BRP)是一類(lèi)NP-難問(wèn)題[1]。其可描述為利用有相同載重的車(chē)輛從單車(chē)維護(hù)補(bǔ)給中心出發(fā)完成一定區(qū)域內(nèi)各單車(chē)停放點(diǎn)自行車(chē)的再分配,使各停放點(diǎn)達(dá)到預(yù)先設(shè)置容量后回到中心,如何使所有的車(chē)輛行駛的距離或花費(fèi)的時(shí)間成本最少。20 世紀(jì)60 年代,公共自行車(chē)系統(tǒng)最早在荷蘭阿姆斯特丹出現(xiàn),隨著社會(huì)經(jīng)濟(jì)的不斷發(fā)展,特別是交通擁堵、環(huán)境污染等城市問(wèn)題不斷突顯,公共自行車(chē)作為健康、環(huán)保、便捷的出行方式逐步在世界各地普及。目前,全世界約有49個(gè)國(guó)家共計(jì)超過(guò)500個(gè)公共自行車(chē)系統(tǒng),如法國(guó)巴黎(20 600 輛)、中國(guó)杭州(78 000 輛)、中國(guó)武漢(90 000 輛)、中國(guó)株洲(20 000輛)[2]。公共自行車(chē)系統(tǒng)的普及使各地公共自行車(chē)運(yùn)營(yíng)部門(mén)單車(chē)再平衡任務(wù)加劇,通常公共自行車(chē)運(yùn)營(yíng)公司組織自有運(yùn)輸車(chē)輛從維護(hù)補(bǔ)給中心出發(fā),按照事先各停靠點(diǎn)分配單車(chē)取送量依次訪(fǎng)問(wèn)各停靠點(diǎn),再回到維護(hù)補(bǔ)給中心,在該過(guò)程中,運(yùn)輸車(chē)輛既可以在各停放點(diǎn)間調(diào)配自行車(chē),也可在出發(fā)前從維護(hù)補(bǔ)給中心載入一部分自行車(chē)調(diào)配到各停放點(diǎn),或?qū)⒏魍7劈c(diǎn)多余的自行車(chē)運(yùn)回維護(hù)補(bǔ)給中心。
國(guó)外有關(guān)公共自行車(chē)再平衡的研究報(bào)道目前主要關(guān)注二類(lèi)BRP,一是靜態(tài)BRP[3-9],即在車(chē)輛開(kāi)始服務(wù)各停放點(diǎn)后,各點(diǎn)平衡數(shù)量不發(fā)生變化,如在夜間(22點(diǎn)至第二天凌晨6點(diǎn)間)用戶(hù)幾乎不使用自行車(chē)的時(shí)段進(jìn)行再平衡作業(yè);二是動(dòng)態(tài)BRP[10-12],即在車(chē)輛開(kāi)始平衡作業(yè)服務(wù)后,各停靠點(diǎn)單車(chē)平衡數(shù)量依舊發(fā)生變化,如在公共自行車(chē)使用頻率較高的區(qū)域,再平衡期間仍有用戶(hù)用車(chē)。從BRP的問(wèn)題特點(diǎn)來(lái)看,求解目標(biāo)有最小化車(chē)輛行駛費(fèi)用[11]、最小化再平衡總費(fèi)用(該費(fèi)用包括車(chē)輛運(yùn)輸費(fèi)用、裝卸單車(chē)費(fèi)用)[10,12]、最小化作業(yè)時(shí)間[3-8]、最大化客戶(hù)滿(mǎn)意度(滿(mǎn)意度表示為站點(diǎn)初始單車(chē)保有量與用戶(hù)取還車(chē)頻率的函數(shù),用戶(hù)在站點(diǎn)取還車(chē)失敗的可能性越高,滿(mǎn)意度越低,反之則越大)[9-10,12]。從BRP求解方法來(lái)看,目前報(bào)道較多的主要有三類(lèi),一類(lèi)是精確算法,如分枝定界法[3-5];二是各類(lèi)經(jīng)典啟發(fā)式算法,如節(jié)約算法[6]、貪婪算法[7],可變鄰域下降搜索[8]等;三是智能啟發(fā)式算法,如禁忌搜索算法[9]。表1 對(duì)國(guó)外主要的研究報(bào)道進(jìn)行了梳理總結(jié)。
我國(guó)雖然是公共自行車(chē)保有量最大的國(guó)家,但目前國(guó)內(nèi)公共自行車(chē)再平衡的研究報(bào)道比較少,白雪等[13]研究了考慮維修車(chē)輛的公共自行車(chē)系統(tǒng)再平衡問(wèn)題,提出了基于動(dòng)態(tài)規(guī)劃的精確算法。前期也對(duì)求解BRP的遺傳算法[14]進(jìn)行了研究,運(yùn)用標(biāo)準(zhǔn)算例測(cè)試表明其具有較好的求解效率。此外,葉麗霞[15]、喬曉[16]對(duì)我國(guó)南京市江寧區(qū)、西安市雁塔區(qū)的公共自行車(chē)再平衡作業(yè)進(jìn)行了案例研究。
綜合分析來(lái)看,國(guó)外對(duì)BRP 的研究起步早,成果較多,對(duì)BRP問(wèn)題特點(diǎn)及精確算法、經(jīng)典啟發(fā)式算法作了較深入的研究,但應(yīng)用智能啟發(fā)式算法求解BRP的研究比較少,尚沒(méi)有求解BRP的人工免疫克隆選擇算法的報(bào)道,而國(guó)內(nèi)研究多為案例研究,研究成果難以推廣應(yīng)用,因此均不適合開(kāi)發(fā)具有一定通用性的計(jì)算機(jī)調(diào)度軟件來(lái)輔助公共自行車(chē)系統(tǒng)再平衡。在實(shí)踐中,公共自行車(chē)的再平衡是一項(xiàng)周期性的調(diào)度工作,用戶(hù)的用車(chē)行為具有慣性(即一段時(shí)間內(nèi)各??奎c(diǎn)取車(chē)、用車(chē)數(shù)量不會(huì)出現(xiàn)劇烈波動(dòng)),前一天的再平衡調(diào)度方案對(duì)第二天的調(diào)度方案優(yōu)化具有重要的參考價(jià)值。
人工免疫算法是模仿生物免疫系統(tǒng)的一種智能優(yōu)化方法。生物免疫系統(tǒng)通過(guò)構(gòu)建具有動(dòng)態(tài)性、自適應(yīng)性和自組織性的信息防御系統(tǒng),以此來(lái)抵御外部無(wú)用、有害和干擾信息的侵入,從而保證系統(tǒng)接受信息的有效性與無(wú)害性。代表性的人工免疫算法是巴西人工免疫學(xué)專(zhuān)家Castro 等人[17]借鑒生物免疫系統(tǒng)的克隆選擇原理提出的克隆選擇算法,該算法通過(guò)模擬生物免疫系統(tǒng)的對(duì)抗原初次免疫應(yīng)答和二次免疫應(yīng)答機(jī)制,能較好地處理優(yōu)化求解中基于過(guò)往方案優(yōu)化現(xiàn)有方案的場(chǎng)景[18],作為一種全局性鄰域搜索智能啟發(fā)式算法,在處理NP-難問(wèn)題時(shí)有較好效果,Kim[19]與Coello[20]分別運(yùn)用該方法處理了多目標(biāo)的數(shù)值優(yōu)化問(wèn)題,文獻(xiàn)[21]則運(yùn)用該方法求解帶工作時(shí)間與時(shí)間窗的開(kāi)放式車(chē)輛路徑問(wèn)題,均顯示出較好的求解效果。
表1 國(guó)外BRP研究代表性文獻(xiàn)
借鑒文獻(xiàn)[5],BRP 的數(shù)學(xué)模型表示如下:設(shè)有一完備圖G=(V,E),V={1,2,…,n}表示頂點(diǎn)集,其中1表示自行車(chē)維護(hù)與補(bǔ)給中心,V′ {2,3,…,n}表示單車(chē)停放點(diǎn)集,E={(i,j)|i,j ∈V,i ≠j}表示弧集或邊集。其符號(hào)定義為:M 表示實(shí)施再平衡運(yùn)輸?shù)膶?zhuān)用車(chē)輛數(shù);Q表示專(zhuān)用運(yùn)輸車(chē)輛的最大載重量,di表示車(chē)輛在第i 個(gè)單車(chē)停放點(diǎn)取送量,di的取值范圍為[-Q,Q],di小于0,表示該停放點(diǎn)需補(bǔ)充單車(chē),di大于0,表示該停放點(diǎn)要回收單車(chē);cij表示專(zhuān)用車(chē)輛訪(fǎng)問(wèn)弧(i,j)產(chǎn)生的運(yùn)輸成本或時(shí)間;θj表示專(zhuān)用車(chē)輛訪(fǎng)問(wèn)第j 個(gè)停放點(diǎn)后的載重量,即車(chē)輛線(xiàn)路載重量。決策變量xij=1,表示車(chē)輛訪(fǎng)問(wèn)了弧(i,j),否則,xij=0。
式(1)為目標(biāo)函數(shù),為最小化運(yùn)輸總成本或時(shí)間,式(2)、式(3)確保除維護(hù)補(bǔ)給中心點(diǎn)外其余單車(chē)停放點(diǎn)均被訪(fǎng)問(wèn)且只訪(fǎng)問(wèn)一次,式(4)、式(5)確保所有的專(zhuān)用運(yùn)輸車(chē)輛在完成任務(wù)后均回到維護(hù)與補(bǔ)給中心,式(6)為消除子回路約束,式(7)、式(8)、式(9)為BRP車(chē)輛載重約束,確保各專(zhuān)用運(yùn)輸車(chē)輛在實(shí)施單車(chē)回收與投放過(guò)程中,車(chē)輛線(xiàn)路載重量不超過(guò)額定載重量,其假設(shè)在某一??奎c(diǎn)回收的自行車(chē)可被投放到任意其他需要的點(diǎn),且車(chē)輛出發(fā)或者回到補(bǔ)給中心時(shí),載重可不為零,因?yàn)榭梢允孪妊b入部分單車(chē),或者帶回部分單車(chē)。
求解BRP 的克隆選擇算法基本步驟如圖1 所示。其中抗體的編碼方式采用多維整數(shù)編碼方法[14],初始解生成借鑒文獻(xiàn)[22]的方法,采用隨機(jī)選擇??奎c(diǎn)、插入位置插入。算法的終止條件設(shè)定為迭代次數(shù)Max_iter或當(dāng)前最優(yōu)解連續(xù)不改進(jìn)的次數(shù)Max_Not_Improve 達(dá)到事先設(shè)定的常數(shù)。
圖1 克隆選擇算法基本流程圖
算法運(yùn)行中設(shè)定了一定容量的記憶細(xì)胞群體Memory,算法每迭代100 次,將抗體群中的最優(yōu)抗體存入Memory,算法終止后,從當(dāng)前Memory 輸出最優(yōu)抗體。當(dāng)有新的問(wèn)題求解時(shí),也可直接從Memory 中選取一定數(shù)量的記憶細(xì)胞,直接計(jì)算是否滿(mǎn)足新的問(wèn)題約束條件,若滿(mǎn)足,可直接計(jì)算抗體與抗原的親和力,進(jìn)入算法后續(xù)步驟迭代,其過(guò)程如圖1 中虛線(xiàn)部分所示,這部分反映了人工免疫系統(tǒng)中特有的二次免疫應(yīng)答機(jī)制,即當(dāng)有新的抗原侵入免疫系統(tǒng)時(shí),可直接借助成熟的記憶細(xì)胞展開(kāi)免疫應(yīng)答,以達(dá)到快速形成免疫能力的目標(biāo)。
在實(shí)踐中,由于部分公共自行車(chē)用戶(hù)需求以通勤出行為主,因此公共自行車(chē)各??奎c(diǎn)的再平衡量作業(yè)量在相鄰的幾天內(nèi)變動(dòng)較小,前一天的再平衡調(diào)度方案對(duì)形成新的再平衡調(diào)度方案有重要的參考價(jià)值。同時(shí)分析BRP模型來(lái)看,約束條件式(9)為兩不等式,表明可行線(xiàn)路中車(chē)輛通過(guò)每點(diǎn)的裝載量有一定的變動(dòng)空間。因此,在克隆選擇算法中引入二次免疫應(yīng)答機(jī)制對(duì)提升算法求解效率有幫助。
在本文算法中,求解目標(biāo)為使目標(biāo)函數(shù)式(1)最小。因此將抗原與抗體的親和力定義為目標(biāo)函數(shù)式(1)的倒數(shù),即目標(biāo)函數(shù)值越小,則該抗體與抗原的親和力越大。
克隆選擇算子設(shè)計(jì)參考文獻(xiàn)[21]的方法,主要有以下兩步:
步驟1 計(jì)算當(dāng)前抗體群中每個(gè)抗體與抗原的親和力,并按親和力降序排列。
步驟2 按照下式對(duì)種群中親和力高的抗體進(jìn)行克隆,得到新的抗體群Nc:
式中,Scale 代表抗體群的規(guī)模,Pos 表示抗體群在降序排列后該抗體在群中的序位,Pos 為整數(shù)且滿(mǎn)足:0 <Pos <;Nc代表克隆后新的抗體群,Nc >Scale;β 為克隆系數(shù),其取值區(qū)間為:0.2 <β <0.3。
抗體經(jīng)過(guò)克隆后則經(jīng)歷高頻變異,高頻變異是克隆選擇算法進(jìn)行鄰域搜索的主要手段。本文采用混合鄰域結(jié)構(gòu),即隨機(jī)選擇點(diǎn)對(duì)操作對(duì)當(dāng)前解的鄰域進(jìn)行搜索,設(shè)計(jì)了以下五類(lèi)變異算子。為了便于說(shuō)明,設(shè)S1=(1356 147928),編號(hào)1 代表維修保養(yǎng)中心,選取進(jìn)行變換的點(diǎn)對(duì)為i=5,j=9,具體變換方法如下:
(1)頂點(diǎn)前向重新指派
將所挑選的頂點(diǎn)i 從線(xiàn)路上當(dāng)前的位置中取出,并將其插入到所挑選頂點(diǎn)j 的位置之前。例如:S1=(1356 147928)→S2=(136 1475928)。
(2)頂點(diǎn)后向重新指派
將所挑選的頂點(diǎn)i 從線(xiàn)路上當(dāng)前的位置中取出,并將其插入到所挑選頂點(diǎn)j 的位置之后。例如:S1=(1356 147928)→S2=(136 1479528)。
(3)頂點(diǎn)交換
將所挑選的頂點(diǎn)i、j 交換位置。例如:S1=(1356 147928)→S2=(1396 147528)。
(4)尾巴交換
若所挑選的兩個(gè)頂點(diǎn)i、j 位于不同的線(xiàn)路上,則將所挑選的兩個(gè)頂點(diǎn)后面的“尾巴”(從被選頂點(diǎn)至線(xiàn)路末尾)互換。例如:S1=(1356 147928)→S2=(13928 14756)。
(5)2-opt
若所挑選的兩點(diǎn)i、j 位于同一條線(xiàn)路上,則將i、j兩點(diǎn)間所有頂點(diǎn)的排列順序逆轉(zhuǎn)。若所挑選的兩點(diǎn)i、j 位于不同線(xiàn)路上,則執(zhí)行以下四種變換:變換1,S1=(1356 147928)→S2=(13829 14765);變換2,S1=(1356 147928)→S2=(139741 6528);變換3,S1=(1356 147928)→S2=(97416 53128);變換4,S1=(1356 147928)→S2=(8296 147531)。
以上幾種變異方法隨機(jī)采用,且為兼顧算法搜索的質(zhì)量與速度,算法前期(當(dāng)前迭代次數(shù)<Max_iter/2)時(shí)采取隨機(jī)變異,即變異后的抗體可行,且與原抗體目標(biāo)值不同即接受變異,后期(當(dāng)前迭代次數(shù)≥Max_iter/2)時(shí)采取尋優(yōu)變異,即每種變異取不同的點(diǎn)對(duì)執(zhí)行5 次,若變異后的抗體可行,且親和力優(yōu)于原抗體,則接收變異。
為進(jìn)行抗體抑制,定義抗體相似性L(fǎng)xy的度量方法如下:
設(shè)x 與y 為同一抗體種群中兩個(gè)不同抗體,p 為抗體邊矩陣(n 維方陣,n 為完備圖G 中頂點(diǎn)數(shù)量),若抗體x 中有邊E(i,j),則=1,否則為0(i,j ∈V,i ≠j)。
式(11)中,函數(shù)sum(p)的功能為計(jì)算抗體邊矩陣中各元素之和,若抗體x 與y 完全相同,則其邊矩陣px=py,則Lxy=0;若抗體x 與y 完全不相同,則Lxy=1。
抗體群經(jīng)歷克隆與高頻變異后,產(chǎn)生了一些親和力低或與其他抗體結(jié)構(gòu)相近的抗體,通??贵w抑制過(guò)程即將這些抗體從種群中刪除并補(bǔ)充新的抗體到原抗體群中,該過(guò)程有利于增加算法的搜索區(qū)域,避免算法陷入局部收斂。本文抗體抑制的具體步驟如下:
步驟1 定義抗體相似度系數(shù)λ ,其取值范圍為[0.3,0.6]。
步驟2 計(jì)算抗體與抗原的親和力,按親和力降序排列,得抗體群Pop(a1,a2,…,aNc),Nc 為抗體群規(guī)模。
步驟3 將a1加入新的抗體群Popnew。
步驟4 從a2開(kāi)始,依次將Pop 中的抗體與Popnew中當(dāng)前抗體比較,若相似度小于等于λ,則將其加入到Popnew中,直到Popnew的規(guī)模達(dá)到Scale,若Pop 中的抗體數(shù)量不足,則從記憶細(xì)胞庫(kù)中選擇記憶細(xì)胞加入,若還不足,則生成新抗體加入到Popnew中。
算法采用matlab2014 編程實(shí)現(xiàn),運(yùn)行于CPU(core i3,3.1 GHz)、ROM(4 GB)的PC機(jī)上。
首先運(yùn)用BRP標(biāo)準(zhǔn)測(cè)試算例,測(cè)試人工免疫克隆選擇算法初次應(yīng)答求解效率,主要參數(shù)設(shè)置為:變異概率取0.9,種群規(guī)模Scale=50,β=0.3,λ=0.4 ,最大循環(huán)次數(shù)Search_Max 為4 000,Not_Improve 為800。每個(gè)算例運(yùn)行10 次,取最好解與平均運(yùn)行時(shí)間。將本文算法與分支定界法(B&C)[5]、遺傳算法[14]進(jìn)行比較(三種算法的測(cè)試硬件環(huán)境基本相同)。
表2、表3 中,Si代表各算法的求解結(jié)果,CPUi表示各算法求解CPU消耗,N/C 分別代表算例的??奎c(diǎn)數(shù)量與運(yùn)送卡車(chē)額定載重量,g_B&C=(S2-S1)/S2、g_GA=(S3-S1)/S3,表示本文算法相較于B&C、GA 在求解質(zhì)量上的差距,g_t_B&C=(CPU2-CPU1)/CPU2、g_t_GA=(CPU3-CPU1)/CPU3,為本文算法相較于B&C、GA在求解時(shí)間上的差距。(1)在規(guī)模小于50個(gè)點(diǎn)算例上,本文算法能找到所有算例的當(dāng)前最好解,平均CPU 消耗為9.6 s,比B&C 的快96.80%,但比GA 慢,如表2所示;(2)在規(guī)模為50至100個(gè)點(diǎn)的算例,本文算法找到6 個(gè)算例的最優(yōu)解(表3 中加粗表示),平均求解質(zhì)量比B&C 低7.43%,但比GA 高1.02%,平均CPU 消耗相較于B&C快96.8%。
為了測(cè)試人工免疫克隆選擇算法二次應(yīng)答的求解效率,本文對(duì)BRP 標(biāo)準(zhǔn)測(cè)試問(wèn)題Roid(N=55/C=20)的停靠點(diǎn)取送車(chē)量(d)進(jìn)行修改,生成3組新的測(cè)試問(wèn)題,用以模擬不同場(chǎng)景下公共自行車(chē)運(yùn)營(yíng)系統(tǒng)連續(xù)兩天的再平衡需求。修改方法為:分別取p=1、2、3,隨機(jī)選取Roid 中80%的??奎c(diǎn),對(duì)原??奎c(diǎn)取送單車(chē)數(shù)量按d±p 進(jìn)行修改(p 的正負(fù)號(hào)隨機(jī)選取,但需確保新的??奎c(diǎn)取送車(chē)數(shù)量不等于0),每組算例生成4個(gè)新的問(wèn)題。算法測(cè)試的參數(shù)設(shè)置除Not_Improve=400 外,其他參數(shù)設(shè)置與第一部分測(cè)試相同。二次應(yīng)答測(cè)試方法為算法先運(yùn)行Roid問(wèn)題,并將求解過(guò)程中的記憶細(xì)胞保存,再分別運(yùn)行改進(jìn)后的測(cè)試問(wèn)題,記錄運(yùn)行結(jié)果及CPU消耗。
測(cè)試算例Roid中各點(diǎn)再平衡車(chē)輛數(shù)的取值區(qū)間為[-7,7],表4~6分別模擬測(cè)試連續(xù)兩天內(nèi),再平衡量變化幅度為14.3%(p=1)、28.6%(p=2)、42.9%(p=3)的情況下,算法二次應(yīng)答效率。測(cè)試結(jié)果表明:二次應(yīng)答的求解質(zhì)量分別比一次應(yīng)答的高0.65% (p=1) 、0.24%(p=2) 、0.16% (p=3) ,CPU 消耗比一次應(yīng)答分別快39.78%(p=1)、46.18%(p=2)、40.66%(p=3)。表4~6中:g_32=(S1-S2)/S1,t_32=(CPU1-CPU3)/CPU3。
表2 小于50個(gè)點(diǎn)的算例測(cè)試結(jié)果
表3 大于50個(gè)點(diǎn)小于100個(gè)點(diǎn)的算例測(cè)試結(jié)果
公共自行車(chē)系統(tǒng)在我國(guó)普及度高,系統(tǒng)的再平衡調(diào)度方法研究有著非常廣泛的應(yīng)用價(jià)值。本文在已有研究基礎(chǔ)上,設(shè)計(jì)了求解BRP的人工免疫克隆選擇算法,算法采用多維整數(shù)編碼方法,結(jié)合問(wèn)題特點(diǎn)設(shè)計(jì)了新的抗體相似性定義方法及抗體抑制策略,并結(jié)合BRP周期性調(diào)度特點(diǎn)引入二次應(yīng)答求解機(jī)制。運(yùn)用標(biāo)準(zhǔn)算例測(cè)試表明:(1)一次應(yīng)答測(cè)試中,本文算法在求解規(guī)模小于50個(gè)點(diǎn)的問(wèn)題中,均能找到最優(yōu)解,求解質(zhì)量與分枝定界法相當(dāng),但速度更快;(2)在求解規(guī)模為50個(gè)點(diǎn)到100個(gè)點(diǎn)間的算例上,算法求解質(zhì)量相較于GA提高1.02%;(3)二次應(yīng)答測(cè)試中,二次應(yīng)答的求解質(zhì)量比一次應(yīng)答的略高,但CPU消耗比一次應(yīng)答快39%以上。
表4 二次應(yīng)答測(cè)試(p=1)
表5 二次應(yīng)答測(cè)試(p=2)
表6 二次應(yīng)答測(cè)試(p=3)
結(jié)合本文算法一次應(yīng)答、二次應(yīng)答的求解測(cè)試結(jié)果,本文認(rèn)為,求解BRP的人工免疫克隆選擇算法的求解速度較精確算法更有優(yōu)勢(shì),求解質(zhì)量與遺傳算法相差不大,其二次應(yīng)答機(jī)制能加快周期性問(wèn)題的求解速度,在實(shí)現(xiàn)公共自行車(chē)再平衡實(shí)時(shí)調(diào)度方面能體現(xiàn)更好的實(shí)用性。公共自行車(chē)再平衡是一項(xiàng)系統(tǒng)工程,即與??奎c(diǎn)的規(guī)劃布置有關(guān),也與用戶(hù)取用車(chē)習(xí)慣有關(guān),還與運(yùn)營(yíng)平臺(tái)再平衡調(diào)度管理策略有關(guān),將這些因素融入BRP模型,設(shè)計(jì)新的問(wèn)題類(lèi)型與對(duì)應(yīng)求解方法將是后續(xù)BRP問(wèn)題的重要研究方向之一。