求解無容量設施選址問題的混合蝙蝠算法

劉春苗，張惠珍

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

1 引言

無容量設施選址問題可以應用于許多不同的領域。無容量設施選址[1]問題雖然相對來說易于理解描述，但是卻很難對它進行求解，它是一個NP難題。為了更好地選擇相應的地址來建造一些設施，如醫院、超市、學校、地鐵站等，從而使得客戶的需求得到最大化的滿足，消耗成本最小。

目前為止，可用在求解無容量設施選址問題的算法可概括性地分成：精確型算法（分支定界法[2-3]、列生成算法、割平面算法等），近似算法（1.488-近似算法[4]），智能優化算法（蟻群算法[5]、禁忌搜索算法、遺傳算法、粒子群算法[6]等）。用精確算法來求解無容量設施選址問題雖然能求得最優解，但這只能用于求解規模較小的無容量設施選址問題，并且計算速度較慢。近似算法是指在多項式時間內能夠求得問題的一個解，并且其目標函數值與最優解的目標函數值之比不超過一個常數的算法[7]。智能優化算法雖適用于求解規模大的無容量設施選址問題[8]，但往往很容易陷入局部最優，一般最后求得的是滿意解而不是問題的最優解。綜合精確型算法的優缺點和智能優化算法的優缺點來看：一方面，改進精確算法的復雜度和應用范圍，提高其運行效率。另一方面，迫切需要設計一種可以融合其他智能算法的優點，并彌補自身不足的混合智能優化算法，這不僅是求解無容量設施選址問題的研究趨勢，也是求解其他任何一個NP難題的趨勢。

蝙蝠算法[9]在求解全局優化問題時，具有求解速度快、效率高、通用性強、可以和其他算法很好地結合等特點，但也存在一些缺點。盛孟龍等人針對蝙蝠容易陷入局部最優和尋優精度不高的問題，通過引入一種交叉變換的方式重新更新蝙蝠種群的位置，減少了算法陷入局部極值的可能，增強了算法的尋優精度。尹進田等人使用混沌映射的方法來對算法進行初始化，增強了種群的多樣性。李枝勇等人根據蝙蝠算法容易與其他算法結合的特點，將量子進化算法和蝙蝠算法相結合，提高蝙蝠的全局搜索能力。本文通過分析無容量設施選址的具體特征，改變原有的游走規則，將三種局部搜索策略、和聲搜索機制與蝙蝠算法進行合理的結合，用混合蝙蝠算法求解無容量設施選址問題。本文通過求解算例，并與其他算法作比較來驗證該算法的可行性和有效性。

2 無容量設施選址問題的數學模型

無容量設施選址問題一般描述為：假設待選設施的容量都沒有限制，客戶用I={1,2,…,n}來表示，待選設施的位置用J={1,2,…,m}來表示，cij表示設施 j到客戶i之間的運輸成本，fj表示設施 j的建造費用。該模型的目的是確保在每個客戶的需求最大化滿足的情況下確保最低的成本。

其中，I表示客戶集；i表示第幾個客戶，?i∈I；n表示客戶總數；J表示設施集；j表示第幾個設施，?j∈J；m表示設施總數；xij=1表示設施位置 j為客戶i提供服務，否則xij=0；yj=1表示在 j處建設服務設施，否則yj=0。

上述模型中，式（1）為目標函數，表示最小化運輸成本與設施建造費之和；約束（2）表示每個客戶必須被服務且只能被服務一次；約束（3）表示只有開放的設施才能為客戶服務；約束（4）和（5）表示變量的取值范圍。

記 X*是UFL問題（1）～（5）的最優解集，通過具體分析UFL問題的具體特征，得到如下定理。

定理1假設(x*,y*)∈X*為最優解集中的一個解，記對于那么存在：

證明（1）假設存在，使下式成立：

那么定義給定的無容量設施選址問題的解(x?,y?)為：

計算使得下式成立的 j″：

這與(x*,y*)∈X*是矛盾的。

（2）對于 j0∈J,假設存在i0∈Ij0,使計算使得下式成立的 j″：

并定義給定的無容量設施選址問題的解(x?,y?)為：

(x*,y*)是求解無容量設施選址問題的一個最優解，Ij表示設施 j為哪些客戶提供服務。通過定理1可以看出，無容量設施選址問題的最優解具備兩個特征：

（2）在已經開放的設施里，客戶都會選擇與之相對應消耗費用最少的設施為他服務。

3 基本蝙蝠算法

蝙蝠算法[10-11]是美國學者楊新社在2010年受到蝙蝠捕捉獵物的啟發而提出的一種新型智能啟發式算法，蝙蝠通過自己特有的回聲定位方法來檢測距離，采用獨特的方式來區分障礙物與獵物，蝙蝠通過發射出響亮的聲音脈沖，根據周圍物體反射回來的回聲響度和自身的雙耳時間差額去建立一個三維環境的場景。蝙蝠算法[12]中尋找最優目標函數值的過程，就是將虛擬蝙蝠類比成當前的可行域內所分布的搜索點，蝙蝠通過不斷飛行來尋找獵物。

3.1 虛擬蝙蝠的運動

通常需要定義虛擬蝙蝠i在d維搜索獵物的空間中它的位置xi和運行速度vi的更新方式，t時刻蝙蝠i的速度和位置的更新公式為：

其中，β∈[ ]0,1是一個服從均勻分布的隨機向量；fi為蝙蝠i的頻率；x*表示當前全局最優位置（解）。在實際的問題求解過程中，可以根據需要搜索范圍的大小來確定 fi的值，比如令 fmin=0，fmax=100。初始時，每只蝙蝠按間的均勻分布隨機賦給一個頻率。

與其他仿真算法相類似的是，運行過程中出現了最優解集，從這個最優解集中選擇了一個最優解，此時每只蝙蝠又隨機產生了局部解[13-15]。局部搜索時，每只蝙蝠的新位置可通過式（9）產生：

3.2 音量和脈沖發生率

通常情況下，在接近獵物的過程中，響度會逐漸降低，脈沖發射的速率會逐漸提高。音量Ai和脈沖發射速率ri的更新公式如下：其中，α和γ分別表示音量衰減系數和脈沖頻度增加系數，隨著蝙蝠接近獵物的過程，音量是不斷降低，脈沖頻率是不斷增加的[16]，二者均為常量。實際上，α類似于模擬退火算法中冷卻進程中的冷卻因素。不難發現，對于任意的0＜α＜1，γ＞0，都有：

實際應用中，通常使用α=γ=0.9。

4 求解無容量設施選址問題混合蝙蝠算法

4.1 開放設施的改進方法

局部搜索方式1（交換法）采用實數編碼的方式，在已經開放的設施中，通過隨機交換一部分不同客戶所對應的設施來改進當前解。

例如：如圖1所示，假設客戶i1由設施1服務，客戶i2由設施2服務，采用交換法對圖1中的設施進行變換，以(i1,1)、(i2,2)替代(i1,2)、(i2,1)。

圖1 交換法示意圖

需要注意的是：當利用交換法變化前后的成本滿足ci1,1+ci2,2＜ci1,2+ci2,1，則說明交換設施服務得到了改進。

局部搜索方式2、3以定理1為基礎，本文除了運用上述局部搜索方法，還將定理1中的兩大有效的局部搜索策略與蝙蝠算法結合，求解無容量設施選址問題。通過下面的例子來對定理1的這兩種策略進行簡要介紹：

給定一個m=5和n=5的無容量設施選址問題，矩陣C表示客戶到設施之間的運輸費用，F表示設施的建造費用。

上述例子的最優解是101。假設蝙蝠算法中L只蝙蝠搜索的滿意解為x14=1,x23=1,x32=1,x42=1,x55=1,y1=0,y2=1,y3=1,y4=1,y5=1,即設施開放集合為{2，3，4，5}；客戶1選擇了設施4，客戶2選擇了設施3，客戶3和客戶4選擇了設施2，客戶5選擇了設施5；消耗的總費用之和為19+14+12+10+122+174+122+178+114=765。

下面使用兩種局部搜索方法對蝙蝠算法的求解結果進行改進：

搜索策略1本文以設施2服務客戶3和客戶4為例，客戶3和客戶4分別同時被其他各個設施服務所消耗的成本為：客戶3和客戶4被設施1服務的費用為284，客戶3和客戶4被設施2服務的費用為319，客戶3和客戶4被設施3服務的費用為321，客戶3和客戶4被設施4服務的費用為158，客戶3和客戶4被設施5服務的費用為159。根據計算結果以及定理1得出，如果客戶3和客戶4選擇了為其服務成本最低的設施4，則可以使得蝙蝠算法的求解結果得到改進。那么令x32=0，x42=0,x34=1,x44=1,y2=0,y4=1。同理其他的也都可以改進，最后得出結果為x14=1,x23=1,x34=1,x44=1,x55=1,y1=0,y2=0,y3=1,y4=1,y5=1,總的成本為592。

搜索策略2 min{c13,c14,c15}=19，設施5應為客戶1服務，min{c23,c24,c25}=13,設施4應為客戶2服務，min{c33,c34,c35}=12，設施4應為客戶3服務，min{c43,c44,c45}=11，設施5應為客戶4服務，min{c53,c54,c55}=10，設施3應為客戶5服務，v=101。

經過兩個策略計算，蝙蝠算法的求解結果明顯得到了改進。

4.2 游走法則的改進

通過分析蝙蝠算法的主要特性，發現通過式（9）的游走法則來尋找目標的搜索能力較弱，因此本文給出了一個新的游走法則，增強蝙蝠的搜索能力，公式如下：

其中，μ是慣性因子，μ=1-ω,ω=ωmax-(ωmax-ωmin)×t/MaxIter；ψ∈[-1,1]；xold是最佳蝙蝠的位置；xk是從種群中隨機選取的任一只蝙蝠的位置。

4.3 和聲算法

為了進一步改進蝙蝠算法的局部搜索能力，平衡蝙蝠算法的開發能力和探索能力，在蝙蝠算法中融入并改進了和聲搜索機制。

2001年，Geem等人首次提出了和聲搜索（Harmony Search，HS）算法，這種算法原理來自于樂隊演奏，每個樂器的音符都代表著目標函數中的變量，演奏目的是使得音樂更加動聽，同樣算法的優化目的就是尋找到全局最優值。在HS算法中，各樂器聲調的每個和聲相當于一個解，并通過一個n維實向量表示。在和聲搜索算法中存在一個和聲記憶庫，用來存放解集，和聲記憶庫的大小表示種群的大小，存儲在和聲記憶庫中的和聲向量和目標函數值用矩陣來表示。在和聲算法中初始解是隨機產生的，初始化之后，算法通過三種方式產生新的和聲向量xnew，包括記憶考慮、音調調整和隨機選擇，公式如下：

其中，r1,r2∈[0,1]是服從均勻分布的隨機向量；r3是區間[-1,1]之間的隨機數；HMCR∈[0,1]表示和聲記憶庫考慮概率，用它來決定是否從和聲記憶庫中選擇xnew；PAR∈[0,1]表示微調概率，它用來決定變量是否進行調整；BW表示的是音調微調帶寬。

一個新的和聲向量xnew產生以后，和聲記憶將通過比較xnew和最差和聲向量xw的值來更新和聲記憶庫，如果xnew的值優于xw的值，那么xw將被xnew取代。

在本文算法中，采用和聲搜索機制產生候選解xnew，通過與和聲記憶庫中最差的和聲向量xw比較，來決定是否接受產生的新解。為了更好地利用當前最優解，重新更改了和聲搜索機制，公式如下：

其中，τ=1-ω；r1,r2∈[0,1]是一個服從均勻分布的隨機向量；r3、r4是區間[-1,1]之間的隨機數；xk1,j、xk2,j、xk3,j是隨機選取種群中的解；表示當前最優解中的第 j個值，j=1,2,…,n。從上式中可以看出，如果r1≤HMCR,xnew,j是由當前最優解產生的，這個過程類似于和聲算法中的記憶考慮；如果r2≤PAR,xnew,j是通過擾動進行微調，這個過程類似于和聲算法中的音調調整。擔當的角色類似于和聲算法中的BW，如果r1＞HMCR,xnew,j則通過隨機選擇產生，這個過程也類似于和聲算法中的隨機更新規則。在式（16）中，HMCR和RAR是不斷變化的，更新公式如下：

其中，t是當前迭代次數；max cycle是最大迭代次數；HMCR(t)是隨著t變化的和聲記憶庫取值概率；HMCRmax和HMCRmin表示最大保留概率和最小保留概率。

其中，PAR(t)是隨著t變化的和聲記憶庫的微調概率；PARmax和PARmin表示最大擾動概率和最小擾動概率。從式（16）可以看出，當前最優解是從新的候選解中產生的，HMCRmax和HMCRmin是不斷線性增加的。因此，蝙蝠在早期進行隨機搜索，隨著迭代次數的增加，主要集中在當前最優解附近搜索。

4.4 混合蝙蝠算法求解無容量設施選址問題的實現過程

步驟1設置蝙蝠的種群大小為N，最大迭代次數Max Iter，蝙蝠i的初始響度，脈沖發射率，最小脈沖頻率 fmin和最大脈沖頻率 fmax，脈沖音強衰減系數α和脈沖頻度增加系數γ。

步驟2隨機產生每只蝙蝠所對應的初始解，產生的初始解對應的是每個客戶都得到一個開放的設施為他們服務，并且每個客戶僅僅被一個設施服務，如一組解[2 1 3 5 3]則表示5個客戶分別被設施1、2、3、5服務。

步驟3計算與每只蝙蝠所對應的目標函數值fitness(i)(i=1,2,…,N)，并記當前最優解為X*及其對應的目標函數值為gfvalue。

步驟4判斷該算法是否取得了最優解或者達到了最大迭代次數MaxIter，如果取得了最優解或者達到了最大迭代次數，算法結束，此時輸出它的最優解；否則，轉步驟5。

步驟5根據式（6）～（8）計算得出蝙蝠i的一個待定解Xnew。通過頻率和速度更新位置Xnew，不一定每個客戶都能得到正確的分配，此時需要對產生的解進行處理。首先找出未能得到正確分配的客戶，對他們采用運輸成本最低的方式進行重新分配。用交換法對Xnew進行局部優化，直到交換出最優解，則停止使用交換法。

步驟6產生一個隨機數rand1，如果rand1＞R0(i),則用式（14）得到一個新的待定解Xnew，此時運用定理1的兩大局部搜索策略、和聲搜索進一步提高蝙蝠的搜索能力和尋優精度。

步驟7產生一個隨機數rand2，如果rand2＜A(i)且 fitness(Xnew)＜gfvalue，則轉入步驟8。

步驟8接受這個蝙蝠飛行產生的新解，并根據式（10）和（11）更新響度A(i)和脈沖發射率R0(i)。

步驟9排列蝙蝠的位置并找到最佳位置X*，并返回到步驟4。

混合蝙蝠算法流程如圖2所示。

5 算例分析

為了驗證混合蝙蝠算法在求解無容量設施選址問題上的可行性、有效性和優越性，本文采用了UFL基準問題庫中的12個算例數據對目標函數進行求解，并對基本蝙蝠算法、基本蟻群算法、混合蟻群算法和混合蝙蝠算法求解無容量設施選址的計算結果進行了對比分析。利用MATLAB（R2008b）編寫程序，算法的實驗環境為Intel?CoreTMi7-6500U CPU@2.50 GHz，4 GB內存的Windows 7操作系統。

參數的選取可以根據具體的實際情況進行調整，本文測試算法時的參數取值包括：蝙蝠的個數N為20；音量衰減系數用α表示，取值為0.9；脈沖頻度增加系數用γ表示，取值為0.9；最大慣性權重用ωmax表示，取值為0.9；最小慣性權重用ωmin表示，取值為0.9；最大迭代次數 MaxIter為200次；和聲算法中的最大保留概率HMCRmax為0.99，最小保留概率 HMCRmin為0.9，最大擾動概率PARmax為0.5，最小擾動概率PARmin為0.1。如表1所示。

表2將基本蝙蝠算法、混合蝙蝠算法與文獻[5]中的基本蟻群算法、混合蟻群算法的求解結果進行了對比，包括迭代200次的總時間和目標函數值。圖3、圖4和圖5分別是選擇三種規模算例例1（gs00250a1）、例2（gs00500a1）、例3（gs00750a1）用基本蝙蝠算法和混合蝙蝠算法求解的結果對比圖。

表1 參數取值

表2 計算結果

圖2 混合蝙蝠算法流程圖

圖3 BA與Hybrid BA用例1時目標函數的收斂曲線

圖4 BA與Hybrid BA用例2時目標函數的收斂曲線

圖5 BA與Hybrid BA用例3時目標函數的收斂曲線

通過表2和收斂曲線圖得出如下結果：基本蟻群算法、混合蟻群算法、基本蝙蝠算法和混合蝙蝠算法求解12個算例的平均運行總時間分別是651.08、671.63、5.28和81.13。混合蝙蝠算法求解12個算例的平均目標函數值為530 897，無論從時間上或者求解結果上來說都比基本蝙蝠算法、基本蟻群算法和混合蟻群算法的結果要好。加入局部搜索策略、更改種群游走規則的混合蝙蝠算法運行總時間仍然明顯小于基本蟻群算法和混合蟻群算法，其結果也明顯優于基本蝙蝠算法、基本蟻群算法和混合蟻群算法。在求解大規模選址算例時，蝙蝠算法計算出的結果與蟻群算法求出的結果隨著規模的逐步變大，差值也越來越大，從而表明蝙蝠算法的研究價值不容小覷。總之，相對于其他三種算法而言，用混合蝙蝠算法求解大規模問題的可行性更高，收斂性能更好。

6 結束語

本文將蝙蝠算法用于求解無容量設施選址問題，根據UFL問題的具體特性，加入了三種局部搜索方法、和聲搜索機制改進算法，并更新了種群隨機游走規則，設計出了求解UFL問題的混合蝙蝠算法。改進后的算法極大程度地改善了蝙蝠易陷入局部最優，尋優精度不高，后期收斂速度慢的缺點。本文的求解算例結果顯示出了混合蝙蝠算法在求解UFL問題上的可行性、有效性和優越性。主要創新如下：

（1）交換法的加入使得蝙蝠每次搜索解的時候，能夠遍歷整個鄰域空間，每次得出的新解都能得到最大化的改進。

（2）重新定義的游走規則提高了蝙蝠尋找獵物的能力，在求解無容量設施選址問題時，最優解的精度得到明顯的改進。

（3）根據無容量設施選址模型設計的兩大局部搜索策略，不僅符合求解模型的特點，還能很好地與蝙蝠算法結合，有效避免了蝙蝠過早陷入局部最優，使得設施的分配更加合理。

（4）和聲算法增強了蝙蝠的全局搜索能力，平衡了算法的開發能力和探索能力，提高了蝙蝠的收斂速度。

通過分析本次測試算例的整個過程，可以知道蝙蝠算法參數少，比其他算法易于實現，其運行時間大大小于其他算法的運行時間，同時也適應了現在快節奏生活的要求。蝙蝠算法求解的選址問題規模越大，其優點越發突顯，在未來解決復雜性、大規模選址問題上用途十分廣泛。為了擴展蝙蝠算法的應用領域，并進一步研究蝙蝠算法求解相關問題的可行性和優越性，利用混合蝙蝠算法求各類選址問題將是下一步的研究工作。