求解P中位問題的混合蝙蝠算法

（上海理工大學 管理學院，上海 200093）

設施選址問題是運籌學中的經典問題之一，在生產、生活及物流等方面都有著非常廣泛的應用，如工廠、倉庫、急救中心、消防站及物流中心的選址等。而P中位問題（P-median problem，PMP）就是一種很常見的設施選址問題，該問題最初由Hakimi[1-2]在1964年提出的，目的是要在給定沒有容量限制的m個候選設施集中選擇p(p＜m)個作為開放設施，使得所有顧客到這些中位點的距離之和達到最小。

目前，用于求解此問題的方法多種多樣，如動態規劃方法[3]、拉格朗日松弛算法[4-5]等經典求解方法。但是，這些經典方法僅對小規模的問題有效,而當求解問題的規模較大時，其計算復雜度也將極大地增加，乃至無法求解。由于Garey和Johnson[6]在1979年通過計算復雜性理論已經證明了該問題是NP-hard問題，因此，除非P=NP，否則不存在多項式時間算法求解該問題。所以，絕大多數學者都致力于啟發式算法的研究，以期得到問題的滿意解或近似最優解。許多傳統的優化算法已被用來求解PMP，如局部搜索算法[7-8]、隨機搜索算法[9]等。此外，人們也提出了多種現代啟發式算法，如模擬退火算法[10]、遺傳算法[11-13]及粒子群算法[14]等。這些啟發式算法能夠生成較好的最優解或近似最優解，特別是以遺傳算法為代表的仿生算法，以其自組織和自適應等良好特征，被廣大學者研究并應用到組合優化問題的求解過程中。而本文正是將蝙蝠算法這種具有良好優化性能和發展前景的仿生算法用于求解PMP，以期能夠克服已有算法的局限性，獲得較好的優化性能。

蝙蝠算法（bat algorithm，BA）是 Yang[15]在2010年提出的一種新型啟發式算法，這是一種搜索全局最優解的有效方法。自提出以來，該算法已經被證明可以應用于多種實際問題的求解。Nakamura等[16]提出了一種二進制蝙蝠算法來解決特征選擇問題；李枝勇等[17]在元胞自動機原理和基本蝙蝠算法的基礎上提出了一種元胞蝙蝠算法用于0-1規劃問題的研究；Rizk-Allah等[18]給出了一種二元蝙蝠算法實現0-1背包問題的求解。但是，至目前尚未有文獻將BA應用于PMP的求解。因此，本文基于BA求解離散問題的局限性，結合PMP的具體特征，重新定義了BA的操作算子，提出了求解PMP的混合蝙蝠算法（hybrid bat algorithm，HBA），并與其他智能算法如遺傳算法[13]和粒子群算法[14]進行比較，不僅驗證了該混合蝙蝠算法用來求解PMP的有效性，而且拓寬了BA的應用領域。

1 P中位問題

假設給定設施集M={1,2,···,m}，客戶集N={1,2,···,n}，對于任意給定的i∈M，j∈N，dij表示設施i與客戶j之間的距離，所有的候選設施都沒有容量限制，要求在m個候選位置點中選擇p個位置作為開放設施（中位點），并且要求每個客戶必須選擇一個且只能選擇一個開放設施來滿足其需求，同時使總距離最小。PMP的數學模型可被描述為

式中：z為目標函數，即總距離；p為允許開放的設施數目；xij表示客戶j是否由設施i服務，如果客戶j由設施i服務，則xij=1，否則，xij=0；yi表示設施i是否開放，如果設施i開放，則yi=1，否則，yi=0。

2 基本蝙蝠算法

蝙蝠是一種神奇的動物，它靠自身的回聲定位系統來探測獵物，避免障礙物，在黑暗中找到它們的棲息地，而BA是一種基于蝙蝠的回聲定位系統而產生的智能啟發式算法，將蝙蝠映射為問題空間中的可行解，將搜索和優化的過程模擬成蝙蝠尋找獵物的過程，利用適應度函數值的大小來衡量蝙蝠所處位置的優劣，將蝙蝠的優勝劣汰過程類比為搜索和優化過程中用好的可行解替代較差的可行解的迭代過程。

2.1 基本蝙蝠的搜索過程

BA的搜索過程實質上是通過頻率的變化實現對蝙蝠速度的控制，從而更新蝙蝠的位置。假設在d維空間中，蝙蝠i在t時刻下的速度和位置分別為vit和xit，則該蝙蝠通過改變頻率fi來調整自身位置與當前所處位置最好的蝙蝠之間的距離，從而更新自身新位置xit+1，向當前求解的最好解靠攏，則蝙蝠i在t+1時刻下的速度和位置更新公式定義為

式中：β∈[0,1]是服從均勻分布的隨機數；fmin和fmax為頻率的最小值和最大值；vit和vit+1為t和t+1時刻蝙蝠i的速度；x*表示當前求得的最好解；和xit+1分別為t和t+1時刻蝙蝠i的位置。

2.2 基本蝙蝠的優化過程

蝙蝠的優化過程即局部搜索過程，在當前求解的最好解集中選擇一個解xold，那么，每只蝙蝠在該最好解附近按照隨機游走法則產生局部新解xnew。

式中： ε是屬于[-1,1]的p維隨機向量；At表示t時刻當前代蝙蝠種群音量的平均值。

當蝙蝠找到獵物時，音量就會降低，同時脈沖發生率逐漸提高。具體按照式（11）和式（12）對音量和脈沖發生率進行更新。

式中：α和γ是常量；Ait和Ait+1分別為t和t+1時刻蝙蝠i的音量；ri0為初始脈沖發生率；rit+1為t+1時刻蝙蝠i的脈沖發生率。

3 混合蝙蝠算法

BA與蟻群算法、粒子群算法等其他群體智能算法類似，是利用蝙蝠的回聲定位特性來模擬蝙蝠的捕食獵物行為而建立的仿生算法。它的速度和位置更新步驟有些類似于標準粒子群優化，在一定程度上，BA可以看作是標準粒子群優化和強化的局部搜索的平衡結合，其平衡受音量和脈沖發生率的控制。BA的提出是用來求解連續域的函數優化問題，不能直接用來求解組合優化問題，因此，本文根據PMP的具體特征和BA的基本思想，提出了一種適用于求解PMP的混合蝙蝠算法。

3.1 編碼方式

PMP的主要任務是在m個候選位置點中選擇p個位置作為中位點，以這些中位點來服務客戶，而僅當每個客戶都被分配到離其最近的中位點時才能保證總距離最小，則可認為PMP的求解是通過搜索不同組合的p個中位點來尋找問題最優解。因而本文采用基于中位數的方式對個體進行編碼，每只蝙蝠對應一個p維向量，每個維度上的值為[1，m]內的隨機整數，該p維向量內每一維度上的值都是唯一的，且沒有大小順序。例如，對于一個m=5，p=3的PMP，與蝙蝠i對應的編碼為xi=(5 2 3)，即設施2，3，5為中位點。

3.2 相關定義

考慮到BA求解離散問題的局限性，為了將BA成功運用到PMP的求解中，并且有較好的運算效率，根據PMP模型的具體特征，設計了針對求解該問題的相關操作算子。

3.3 混合蝙蝠算法的設計

混合蝙蝠算法的步驟：

步驟1各參(數初始化)。蝙蝠種群大小K；最大迭代次數NmaxN_iter=1；音量A；脈沖發生率r；距離矩陣D；適應度函數f(x)。

步驟2隨機初始化種群，找出當前所處位置最好的蝙蝠x?。

步驟3按照定義1和定義2更新當前蝙蝠的速度vit+1和位置。

步驟 4產生隨機數rand1，判定rand1＞ri是否成立，如果成立，則根據定義3執行局部搜索。

步驟 5產生隨機數rand2，判定rand2＞Ai且f(xi)＜f(x?)是否同時成立，如果成立，則接受這個新解，并更新A，r和當前求得的最好解。

步驟6判斷是否達到預設的最大迭代次數Nmax或時間限制3600s，若達到，則輸出當前求得的最好解；否則，N_iter++，轉至步驟3。

4 數值實驗

為了驗證混合蝙蝠算法的求解性能，在Intel（R）Core（TM）CPU@3.60GHz 處理器、7.87GB內存、操作系統為64位Windows10的實驗環境下，應用MatlabR2016a編程實現了求解PMP的混合蝙蝠算法，對運籌學數據庫OR-Library中的40個PMP算例進行仿真實驗，并與文獻[13]中的遺傳算法（genetic algorithm，GA）和文獻[14]中的粒子群算法（novel discrete particle swarm optimization，NDPSO）進行對比。求解這些算例的混合蝙蝠算法參數設置情況為：最大迭代次數Nmax=100，蝙蝠音量A=0.9，脈沖發生率r=0.25，音量衰減系數a=0.95，脈沖發生率增加系數g=0.05。同時考慮到測試算例大小的不同和運行時間的限制，根據開放中位點的大小來設置蝙蝠種群的數目，當中位點p＜50時，種群規模K為15；否則，種群規模為7。文獻[13]通過C++編程實現了求解PMP的遺傳算法，其中，最大迭代次數為150 000，交叉概率為0.8，變異概率為0.4，種群規模為100；文獻[14]通過C編程實現了粒子群算法，其中，最大迭代次數為1 000，學習因子為0.5，種群規模設置為算例中設施數目的2倍。

表1不僅給出了40個PMP測試算例相應的設施數和中位數，而且給出了這些算例的已知最優解和應用混合蝙蝠算法10次運行測試求得的最好解，并給出最好解與已知最優解之間的Gap值，Gap=（本研究算法所求得的最好解-已知最優解）/已知最優解×100%。其中，GA和NDPSO求得的最好解和Gap值分別來自文獻[13-14]。另外，表1中還給出了HBA和NDPSO的運行時間，由于GA的運行時間在原文獻中未提及，這里就不加以考慮。除此之外，表1的最后還給出了所有算例已知最優解和3種算法求得解的平均值、Gap的平均值和運行時間的平均值。

由表1可知，HBA在求解40個PMP算例時，其中30個算例可以達到已知最優解，另外10個算例求得的最好解盡管稍遜于已知最優解，但其最大Gap值僅為0.654%；而GA僅可以求得7個算例的最優解，NDPSO可以求得16個算例的最優解。另外，從3種算法的平均Gap值來看，GA的平均Gap=3.231%最大，而HBA的平均Gap=0.06%最小，非常接近已知最優解，因此，HBA在與GA和NDPSO相比時，該算法在求解精度上具有明顯的優勢。

為了進一步分析HBA在求解PMP時的有效性和計算復雜度，給出了一個解得最優解的算例pmed14和另一個解得近優解的算例pmed15的求解迭代過程圖，如圖1和圖2所示。由于算例的規模不同而導致運行迭代次數不同，因此，這里以運行時間來作相關分析。圖1在迭代到265s左右、圖2到1400s左右時，2個算例的計算結果還在發生改變。但是，從運行總時間來看，HBA相對于NDPSO而言，需要耗用較多的時間，這是本文的不足之處。總的來說，HBA在求解這些規模不同的PMP算例時能夠求得較好的實驗結果，與已知最優解的誤差較小，且絕大多數算例結果都可以達到已知最優解，故進一步驗證了蝙蝠算法在求解PMP上的可行性與有效性。

5 結 論

針對P中位問題的特點以及蝙蝠算法的尋優機制，重新定義了原有蝙蝠算法的操作算子，并提出了一種可行化函數。同時在局部搜索部分引入了遺傳算法的交叉思想，設計了一種針對求解該問題的局部搜索方式，提出了求解該問題的混合蝙蝠算法。

表1 OR-Library 算例運行結果Tab.1 Operating results of the examples from OR-Library

圖1 HBA 求解 pmed14 的最優迭代圖Fig.1 Optimal iterative graph of pmed14 solved by HBA

圖2 HBA 求解 pmed15 的最優迭代圖Fig.2 Optimal iterative graph of pmed15 solved by HBA

實驗部分求解多個測試算例，并將混合蝙蝠算法的計算結果與文獻中遺傳算法和粒子群算法的計算結果進行對比，測試結果表明，混合蝙蝠算法在確保種群多樣性的基礎上，具有較好的全局優化性能，能夠有效求解P中位問題，驗證了該算法在求解質量上的優勢。但是，與其他改進的智能算法相比，本文提出的混合蝙蝠算法在求解速度上還存在明顯的不足，這將是進一步的研究方向。