劉喜梅，潘立軍

湖南工程學(xué)院 管理學(xué)院，湖南 湘潭411104

1 引言

公共自行車再平衡問(wèn)題（Bike Sharing Rebalancing Problem，BRP）是一類NP-難問(wèn)題[1]。其可描述為利用有相同載重的車輛從單車維護(hù)補(bǔ)給中心出發(fā)完成一定區(qū)域內(nèi)各單車停放點(diǎn)自行車的再分配，使各停放點(diǎn)達(dá)到預(yù)先設(shè)置容量后回到中心，如何使所有的車輛行駛的距離或花費(fèi)的時(shí)間成本最少。20 世紀(jì)60 年代，公共自行車系統(tǒng)最早在荷蘭阿姆斯特丹出現(xiàn)，隨著社會(huì)經(jīng)濟(jì)的不斷發(fā)展，特別是交通擁堵、環(huán)境污染等城市問(wèn)題不斷突顯，公共自行車作為健康、環(huán)保、便捷的出行方式逐步在世界各地普及。目前，全世界約有49個(gè)國(guó)家共計(jì)超過(guò)500個(gè)公共自行車系統(tǒng)，如法國(guó)巴黎（20 600 輛）、中國(guó)杭州（78 000 輛）、中國(guó)武漢（90 000 輛）、中國(guó)株洲（20 000輛）[2]。公共自行車系統(tǒng)的普及使各地公共自行車運(yùn)營(yíng)部門單車再平衡任務(wù)加劇，通常公共自行車運(yùn)營(yíng)公司組織自有運(yùn)輸車輛從維護(hù)補(bǔ)給中心出發(fā)，按照事先各?？奎c(diǎn)分配單車取送量依次訪問(wèn)各停靠點(diǎn)，再回到維護(hù)補(bǔ)給中心，在該過(guò)程中，運(yùn)輸車輛既可以在各停放點(diǎn)間調(diào)配自行車，也可在出發(fā)前從維護(hù)補(bǔ)給中心載入一部分自行車調(diào)配到各停放點(diǎn)，或?qū)⒏魍７劈c(diǎn)多余的自行車運(yùn)回維護(hù)補(bǔ)給中心。

國(guó)外有關(guān)公共自行車再平衡的研究報(bào)道目前主要關(guān)注二類BRP，一是靜態(tài)BRP[3-9]，即在車輛開始服務(wù)各停放點(diǎn)后，各點(diǎn)平衡數(shù)量不發(fā)生變化，如在夜間（22點(diǎn)至第二天凌晨6點(diǎn)間）用戶幾乎不使用自行車的時(shí)段進(jìn)行再平衡作業(yè)；二是動(dòng)態(tài)BRP[10-12]，即在車輛開始平衡作業(yè)服務(wù)后，各停靠點(diǎn)單車平衡數(shù)量依舊發(fā)生變化，如在公共自行車使用頻率較高的區(qū)域，再平衡期間仍有用戶用車。從BRP的問(wèn)題特點(diǎn)來(lái)看，求解目標(biāo)有最小化車輛行駛費(fèi)用[11]、最小化再平衡總費(fèi)用（該費(fèi)用包括車輛運(yùn)輸費(fèi)用、裝卸單車費(fèi)用）[10，12]、最小化作業(yè)時(shí)間[3-8]、最大化客戶滿意度（滿意度表示為站點(diǎn)初始單車保有量與用戶取還車頻率的函數(shù)，用戶在站點(diǎn)取還車失敗的可能性越高，滿意度越低，反之則越大）[9-10，12]。從BRP求解方法來(lái)看，目前報(bào)道較多的主要有三類，一類是精確算法，如分枝定界法[3-5]；二是各類經(jīng)典啟發(fā)式算法，如節(jié)約算法[6]、貪婪算法[7]，可變鄰域下降搜索[8]等；三是智能啟發(fā)式算法，如禁忌搜索算法[9]。表1 對(duì)國(guó)外主要的研究報(bào)道進(jìn)行了梳理總結(jié)。

我國(guó)雖然是公共自行車保有量最大的國(guó)家，但目前國(guó)內(nèi)公共自行車再平衡的研究報(bào)道比較少，白雪等[13]研究了考慮維修車輛的公共自行車系統(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ū)的公共自行車再平衡作業(yè)進(jìn)行了案例研究。

綜合分析來(lái)看，國(guó)外對(duì)BRP 的研究起步早，成果較多，對(duì)BRP問(wèn)題特點(diǎn)及精確算法、經(jīng)典啟發(fā)式算法作了較深入的研究，但應(yīng)用智能啟發(fā)式算法求解BRP的研究比較少，尚沒有求解BRP的人工免疫克隆選擇算法的報(bào)道，而國(guó)內(nèi)研究多為案例研究，研究成果難以推廣應(yīng)用，因此均不適合開發(fā)具有一定通用性的計(jì)算機(jī)調(diào)度軟件來(lái)輔助公共自行車系統(tǒng)再平衡。在實(shí)踐中，公共自行車的再平衡是一項(xiàng)周期性的調(diào)度工作，用戶的用車行為具有慣性（即一段時(shí)間內(nèi)各停靠點(diǎn)取車、用車數(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é)專家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í)間窗的開放式車輛路徑問(wèn)題，均顯示出較好的求解效果。

表1 國(guó)外BRP研究代表性文獻(xiàn)

2 BRP數(shù)學(xué)模型

借鑒文獻(xiàn)[5]，BRP 的數(shù)學(xué)模型表示如下：設(shè)有一完備圖G=(V,E)，V={1,2,…,n}表示頂點(diǎn)集，其中1表示自行車維護(hù)與補(bǔ)給中心，V′ {2,3,…,n}表示單車停放點(diǎn)集，E={(i,j)|i,j ∈V,i ≠j}表示弧集或邊集。其符號(hào)定義為：M 表示實(shí)施再平衡運(yùn)輸?shù)膶Ｓ密囕v數(shù)；Q表示專用運(yùn)輸車輛的最大載重量，di表示車輛在第i 個(gè)單車停放點(diǎn)取送量，di的取值范圍為[-Q,Q]，di小于0，表示該停放點(diǎn)需補(bǔ)充單車，di大于0，表示該停放點(diǎn)要回收單車；cij表示專用車輛訪問(wèn)弧(i,j)產(chǎn)生的運(yùn)輸成本或時(shí)間；θj表示專用車輛訪問(wèn)第j 個(gè)停放點(diǎn)后的載重量，即車輛線路載重量。決策變量xij=1，表示車輛訪問(wèn)了弧(i,j)，否則，xij=0。

式（1）為目標(biāo)函數(shù)，為最小化運(yùn)輸總成本或時(shí)間，式（2）、式（3）確保除維護(hù)補(bǔ)給中心點(diǎn)外其余單車停放點(diǎn)均被訪問(wèn)且只訪問(wèn)一次，式（4）、式（5）確保所有的專用運(yùn)輸車輛在完成任務(wù)后均回到維護(hù)與補(bǔ)給中心，式（6）為消除子回路約束，式（7）、式（8）、式（9）為BRP車輛載重約束，確保各專用運(yùn)輸車輛在實(shí)施單車回收與投放過(guò)程中，車輛線路載重量不超過(guò)額定載重量，其假設(shè)在某一?？奎c(diǎn)回收的自行車可被投放到任意其他需要的點(diǎn)，且車輛出發(fā)或者回到補(bǔ)給中心時(shí)，載重可不為零，因?yàn)榭梢允孪妊b入部分單車，或者帶回部分單車。

3 求解BRP的克隆選擇算法

3.1 算法主要框架

求解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ì)算是否滿足新的問(wèn)題約束條件，若滿足，可直接計(jì)算抗體與抗原的親和力，進(jìn)入算法后續(xù)步驟迭代，其過(guò)程如圖1 中虛線部分所示，這部分反映了人工免疫系統(tǒng)中特有的二次免疫應(yīng)答機(jī)制，即當(dāng)有新的抗原侵入免疫系統(tǒng)時(shí)，可直接借助成熟的記憶細(xì)胞展開免疫應(yīng)答，以達(dá)到快速形成免疫能力的目標(biāo)。

在實(shí)踐中，由于部分公共自行車用戶需求以通勤出行為主，因此公共自行車各?？奎c(diǎn)的再平衡量作業(yè)量在相鄰的幾天內(nèi)變動(dòng)較小，前一天的再平衡調(diào)度方案對(duì)形成新的再平衡調(diào)度方案有重要的參考價(jià)值。同時(shí)分析BRP模型來(lái)看，約束條件式（9）為兩不等式，表明可行線路中車輛通過(guò)每點(diǎn)的裝載量有一定的變動(dòng)空間。因此，在克隆選擇算法中引入二次免疫應(yīng)答機(jī)制對(duì)提升算法求解效率有幫助。

3.2 親和力定義與克隆選擇算子

在本文算法中，求解目標(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ù)且滿足：0 ＜Pos ＜；Nc代表克隆后新的抗體群，Nc ＞Scale；β 為克隆系數(shù)，其取值區(qū)間為：0.2 ＜β ＜0.3。

3.3 變異算子

抗體經(jīng)過(guò)克隆后則經(jīng)歷高頻變異，高頻變異是克隆選擇算法進(jìn)行鄰域搜索的主要手段。本文采用混合鄰域結(jié)構(gòu)，即隨機(jī)選擇點(diǎn)對(duì)操作對(duì)當(dāng)前解的鄰域進(jìn)行搜索，設(shè)計(jì)了以下五類變異算子。為了便于說(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 從線路上當(dāng)前的位置中取出，并將其插入到所挑選頂點(diǎn)j 的位置之前。例如：S1=（1356 147928）→S2=（136 1475928）。

（2）頂點(diǎn)后向重新指派

將所挑選的頂點(diǎn)i 從線路上當(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 位于不同的線路上，則將所挑選的兩個(gè)頂點(diǎn)后面的“尾巴”（從被選頂點(diǎn)至線路末尾）互換。例如：S1=（1356 147928）→S2=（13928 14756）。

（5）2-opt

若所挑選的兩點(diǎn)i、j 位于同一條線路上，則將i、j兩點(diǎn)間所有頂點(diǎn)的排列順序逆轉(zhuǎn)。若所挑選的兩點(diǎn)i、j 位于不同線路上，則執(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)于原抗體，則接收變異。

3.4 抗體相似性定義與抗體抑制

為進(jìn)行抗體抑制，定義抗體相似性Lxy的度量方法如下：

設(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開始，依次將Pop 中的抗體與Popnew中當(dāng)前抗體比較，若相似度小于等于λ，則將其加入到Popnew中，直到Popnew的規(guī)模達(dá)到Scale，若Pop 中的抗體數(shù)量不足，則從記憶細(xì)胞庫(kù)中選擇記憶細(xì)胞加入，若還不足，則生成新抗體加入到Popnew中。

4 算法測(cè)試

算法采用matlab2014 編程實(shí)現(xiàn)，運(yùn)行于CPU（core i3，3.1 GHz）、ROM（4 GB）的PC機(jī)上。

4.1 初次應(yīng)答求解效率測(cè)試

首先運(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)送卡車額定載重量，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%。

4.2 算法二次應(yīng)答的求解效率測(cè)試

為了測(cè)試人工免疫克隆選擇算法二次應(yīng)答的求解效率，本文對(duì)BRP 標(biāo)準(zhǔn)測(cè)試問(wèn)題Roid(N=55/C=20)的?？奎c(diǎn)取送車量（d）進(jìn)行修改，生成3組新的測(cè)試問(wèn)題，用以模擬不同場(chǎng)景下公共自行車運(yùn)營(yíng)系統(tǒng)連續(xù)兩天的再平衡需求。修改方法為：分別取p=1、2、3，隨機(jī)選取Roid 中80%的?？奎c(diǎn)，對(duì)原停靠點(diǎn)取送單車數(shù)量按d±p 進(jìn)行修改（p 的正負(fù)號(hào)隨機(jī)選取，但需確保新的停靠點(diǎn)取送車數(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)再平衡車輛數(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é)果

5 結(jié)束語(yǔ)

公共自行車系統(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)公共自行車再平衡實(shí)時(shí)調(diào)度方面能體現(xiàn)更好的實(shí)用性。公共自行車再平衡是一項(xiàng)系統(tǒng)工程，即與?？奎c(diǎn)的規(guī)劃布置有關(guān)，也與用戶取用車習(xí)慣有關(guān)，還與運(yùn)營(yíng)平臺(tái)再平衡調(diào)度管理策略有關(guān)，將這些因素融入BRP模型，設(shè)計(jì)新的問(wèn)題類型與對(duì)應(yīng)求解方法將是后續(xù)BRP問(wèn)題的重要研究方向之一。