王婷婷 張惠珍 趙玉蘋

摘要無容量設(shè)施選址問題（Uncapacitated Facility Location Problem， UFLP）是一類經(jīng)典的組合優(yōu)化問題，被證明是一種NP-hard問題，易于描述卻難于求解.首先根據(jù)UFLP的數(shù)學(xué)模型及其具體特征，重新設(shè)計了蝙蝠算法的操作算子，給出了求解UFLP的蝙蝠算法.其次構(gòu)建出三種可行化方法，并將其與求解UFLP的蝙蝠算法和拉格朗日松弛算法相結(jié)合，設(shè)計了求解該問題的拉格朗日蝙蝠算法.最后通過仿真實例和與其他算法進(jìn)行比較的方式，驗證了該混合算法用來求解UFLP的可行性，是解決離散型問題的一種有效方式.

關(guān)鍵詞管理科學(xué)與工程；無容量設(shè)施選址問題；拉格朗日蝙蝠算法；拉格朗日松弛算法；蝙蝠算法

中圖分類號TP301.6文獻(xiàn)標(biāo)識碼A

Lagrangian Bat Algorithm for Solving

The Uncapacitated Facility Location Problem

Tingting Wang1， Huizhen Zhang1， Yuping Zhao2

（1.School of Management， University of Shanghai for Science & Technology， Shanghai200093， China；

2. State Grid Shanghai Procurement Company， Shanghai200093， China）

AbstractThe uncapacitated facility location problem is a classical combinatorial optimization problem and has been proved an NPhard problem， easy to describe but difficult to solve. Based on the mathematical model and specific features of UFLP， the operators of bat algorithm and proposed a bat algorithm for solving it has been redesigned. Secondly， three feasible methods were designed and lagrangian bat algorithm was presented by combining these three methods， the bat algorithm for UFLP and lagrangian relaxation algorithm. Lastly， the experimental results showed that the hybrid algorithm can feasibly solve UFLP through simulation examples and comparisons with other algorithms and it is an effective algorithm in solving the discrete problem.

Key wordsmanagement science and engineering； uncapacitated facility location problem； lagrangian bat algorithm； lagrangian relaxation algorithm； bat algorithm

1引言

設(shè)施選址問題是在各需求點已給定的條件下，選擇設(shè)施的最佳位置，使設(shè)施的總成本降到最低[1].無容量設(shè)施選址問題是從給定沒有容量限制的設(shè)施集中選擇要開放的設(shè)施，使其以最小的總費用來服務(wù)所有客戶.到目前為止，求解該問題的方法主要有近似算法、精確算法和智能優(yōu)化算法.近似算法可以在多項式時間內(nèi)求得問題的一個解，并且其目標(biāo)函數(shù)值與最優(yōu)解的目標(biāo)函數(shù)值之比不超過一個常數(shù)[2]；精確算法是一種可以求得最優(yōu)解的算法，但僅適用于小規(guī)模問題，且計算速度較慢.智能優(yōu)化算法更期望在可接受時間內(nèi)獲得問題的最優(yōu)解或近似最優(yōu)解.近年來，許多學(xué)者運用各種智能優(yōu)化算法來解決UFLP：Kole（2014）[3]等首次將蟻群算法用于UFLP的求解，并與粒子群算法進(jìn)行對比，拓寬了蟻群算法的應(yīng)用范圍；Damgacioglu（2015）[4]等為求解該問題提出了一種遺傳算法，并驗證該算法的可行性，同時將未來的研究工作致力于多設(shè)施選址和有容量選址問題；李倩（2016）[5]等針對無容量設(shè)施選址問題的特點，設(shè)計了兩種局部搜索方法，并與基本蟻群算法相結(jié)合，提出了求解該問題的混合蟻群算法，結(jié)果表明該混合算法具有明顯的有效性.

蝙蝠算法（Bat Algorithm， BA）是Yang（2010）[6]提出的一種新型啟發(fā)式算法.BA的提出是用來求解連續(xù)域的函數(shù)優(yōu)化問題，不能直接用來求解組合優(yōu)化問.通過人們的不斷改進(jìn)，先后對BA進(jìn)行了改進(jìn).Adarsh（2016）[7]等在蝙蝠算法的基礎(chǔ)上加入了混沌序列，提出了求解經(jīng)濟(jì)調(diào)度問題的混沌蝙蝠算法，提高了算法的求解性能.劉春苗（2016）[8]等將蝙蝠算法的操作算子重新進(jìn)行了定義，設(shè)計了求解車輛路徑問題的離散蝙蝠算法，這是首次將蝙蝠算法用于車輛路徑問題的研究.

拉格朗日松弛的思想起源于20世紀(jì)60年代.20世紀(jì)70年代初才對拉格朗日松弛算法（Lagrangian Relaxation Algorithm， LR）提出定義.后來，F(xiàn)isher[9]使用該算法來解決調(diào)度問題，同時提出了更新拉格朗日乘子的方法.

深入分析UFLP模型的具體特點，結(jié)合BA和LR各自的優(yōu)點，提出了一種拉格朗日蝙蝠算法，并通過求解系列經(jīng)典UFLP算例驗證了該算法的尋優(yōu)性能.

2無容量設(shè)施選址問題

假設(shè)給定m個設(shè)施、n個客戶，且所有的候選設(shè)施都沒有容量限制，cij表示設(shè)施i與客戶j之間的服務(wù)費用，hi表示設(shè)施i的開放費用，要求每個客戶都應(yīng)由一個且只能一個設(shè)施提供服務(wù)，同時使總費用最小.UFLP的數(shù)學(xué)模型可被描述如下：

min z=∑mi=1∑nj=1cijxij+∑mi=1hiyi，（1）

∑mi=1xij=1，j∈N，（2）

xij

SymbolcB@ yi，i∈M，j∈N，（3）

xij∈0，1，i∈M，j∈N，（4）

yi∈0，1，i∈M.（5）

模型中，z是目標(biāo)函數(shù)，即總費用；xij表示客戶j是否由設(shè)施i服務(wù)，如果客戶j由設(shè)施i服務(wù)，則xij=1，否則xij=0；yi表示設(shè)施i是否開放，如果設(shè)施i開放，則yi=1，否則yi=0.

3 蝙蝠算法

3.1基本蝙蝠算法

蝙蝠算法是一種基于蝙蝠的回聲定位系統(tǒng)而產(chǎn)生的智能啟發(fā)式算法，將蝙蝠類比為當(dāng)前可行域內(nèi)分布的搜索點，用目標(biāo)函數(shù)值的大小來判斷當(dāng)前蝙蝠所處位置的優(yōu)劣，通過蝙蝠的飛行移動來搜索目標(biāo)函數(shù)最優(yōu)解[10].蝙蝠飛行移動過程中的關(guān)鍵操作包括：速度、位置、音量和脈沖發(fā)生率的更新.

3.1.1蝙蝠速度和位置的更新

在d維空間中，蝙蝠i在t+1時刻下的速度和位置更新公式定義如下：

fi=fmin +（fmax -fmin ）β，（6）

vt+1i=vti+（xti-x*）fi，（7）

xt+1i=xti+vt+1i，（8）

其中，β∈[0，1]是服從均勻分布的隨機(jī)數(shù)；fmin和fmax是頻率的最小值和最大值；vit和vit+1是t和t+1時刻蝙蝠i的速度；x*表示當(dāng)前全局最優(yōu)解；xit和xit+1分別是t和t+1時刻蝙蝠i的位置.

在當(dāng)前最優(yōu)解集中選中一個解xold，那么每只蝙蝠在最優(yōu)解附近按照隨機(jī)游走法則可產(chǎn)生局部新解xnew.

xnew=xold+εAt， （9）

其中，ε是屬于[-1，1]的隨機(jī)數(shù)；At表示t時刻當(dāng)前代蝙蝠種群音量的平均值.

3.1.2蝙蝠音量和脈沖發(fā)生率的更新

當(dāng)蝙蝠接近獵物時，音量降低，脈沖發(fā)生率逐漸提高.音量和脈沖發(fā)生率按照以下公式更新：

At+1i=αAti，（10）

rt+1i=r0i1-exp -γt.（11）

其中，α和γ是常量；Ait和Ait+1分別是t和t+1時刻蝙蝠i的音量；ri0是初始脈沖發(fā)生率；rit+1是t+1時刻蝙蝠i的脈沖發(fā)生率.

3.2求解UFLP的蝙蝠算法

3.2.1編碼方式

對m個設(shè)施、n個客戶的UFLP，每只蝙蝠對應(yīng)一個n維向量，每個維度上的值為[1，m]內(nèi)的一個隨機(jī)整數(shù).例如：給定m=3、n=5的UFLP，若某只蝙蝠的位置向量為[2， 1， 2， 1， 1]，則表示設(shè)施1和2開放，3關(guān)閉；設(shè)施1服務(wù)客戶2、4、5；設(shè)施2服務(wù)客戶1和3.用這種方式構(gòu)建UFLP的解，會使每個客戶均被考慮到并限制每個客戶僅能由一個設(shè)施提供服務(wù)，大大降低了計算復(fù)雜度.

3.2.2種群初始化

考慮到UFLP的特點，當(dāng)設(shè)施開放費用遠(yuǎn)高于服務(wù)費用時，可以開放較少的設(shè)施，否則可以通過減少服務(wù)費用來降低總費用.在初始化種群時就不是完全隨機(jī)的，可以帶有一定的方向性.將設(shè)施的平均開放費用與平均服務(wù)費用的比值記為分類指標(biāo)t，設(shè)定t1為t的一個參考值，p1、p2∈[0，1]且p1t1時，如果r

3.2.3速度和位置更新規(guī)則

為簡化算法復(fù)雜度，沒有考慮BA的頻率，而用當(dāng)前蝙蝠和最優(yōu)蝙蝠位置間的漢明距離s來代替速度，再用精英保留策略更新蝙蝠位置.為了保持種群多樣性，用最優(yōu)蝙蝠位置向量中長度為floor（s/2）的基因段來替代當(dāng)前蝙蝠位置中的同一個基因段.其中，第一個基因段是從第1維至第floor（s/2）維，第二個基因段從第2維至第floor（s/2）+1維，按這種長度直至能替換到位置中的第n維，即分別替換n+1-floor（s/2）次，找出多種替換后適應(yīng)度值最小的一種方式，則用該基因片段替換得到的位置即為更新后的蝙蝠位置.

例如：給定m=n=5的UFLP，其服務(wù)費用矩陣C=[1696 1309 1488 1235 1621； 1666 1980 1227 1783 1135； 1354 1249 1240 1831 108； 1224 1320 121 1340 1004； 192 1439 1375 112 1009]、開放費用矩陣H=[130 199 139 127 103]，蝙蝠i的位置為[3， 2， 2， 5， 3]，適應(yīng)度值為5222，x*為[5， 3， 3， 4， 4]，s為5，則每次替換的基因片段長度為floor（s/2）=2，用精英保留策略分別替換4次得到的位置為[5， 3， 2， 5， 3]、[3， 3， 3， 5， 3]、[3， 2， 3， 4， 3]、[3， 2， 2， 4， 4]，對應(yīng)的適應(yīng)度值分別為3329、4305、6487、7370，其中最小的適應(yīng)度值為3329，小于5222，則該蝙蝠的新位置為[5， 3， 2， 5， 3].

3.2.4局部搜索

為進(jìn)一步提高算法的搜索精度，根據(jù)UFLP的最優(yōu)解所具有的特點，在最優(yōu)蝙蝠附近進(jìn)行局部搜索，在當(dāng)前開放的所有設(shè)施中隨機(jī)挑出一個設(shè)施p∈[1，m]，將其服務(wù)的所有客戶由其它m-1個設(shè)施分別來服務(wù)，其中設(shè)施q∈[1，m]且q≠p所需的服務(wù)總費用最小，則最優(yōu)蝙蝠中由設(shè)施p服務(wù)的客戶都改為由設(shè)施q來服務(wù).

4拉格朗日蝙蝠算法

拉格朗日蝙蝠算法（Lagrangian Bat Algorithm， LR_BA）是由LR和求解UFLP的蝙蝠算法（Discrete Bat Algorithm， DBA）結(jié)合而成，嘗試將這兩種算法相結(jié)合，通過計算上下界值來同時逼近問題的最優(yōu)解.在求解UFLP時，DBA得到原問題的一個上界值，而LR通過松弛其中的等式約束得到原問題的下界值，再利用可行化方法將該下界值可行化為原問題的另一個上界值.通過比較兩個上界值，從而得到UFLP的滿意值.

4.1UFLP的拉格朗日松弛

對約束式（2）進(jìn)行松弛，松弛后的目標(biāo)函數(shù)變?yōu)槭剑?2）.

min z=∑mi=1∑nj=1（cij-λj）xij+∑mi=1hiyi+∑nj=1λj.（12）

通過標(biāo)準(zhǔn)次梯度優(yōu)化方法，按式（13）計算步長T，再根據(jù)式（14）更新拉格朗日乘子λ.

Tt=πZUB-ZLB∑nj=1∑mi=1xij-12，（13）

λt+1j=max 0，λtj+Tt1-∑mi=1xij.（14）

其中，ZUB為原問題的上界值，ZLB為松弛問題的下界值；xij為松弛問題的解；π是常數(shù).

根據(jù)UFLP的具體特點，對于任意客戶，其應(yīng)當(dāng)由所需服務(wù)費用最少的設(shè)施服務(wù).對于任意設(shè)施，將其服務(wù)的所有客戶分別由其他設(shè)施來服務(wù)，而該設(shè)施所需的總費用應(yīng)當(dāng)最小，提出了松弛問題的3種可行化方法.

下面通過一個UFLP算例對此進(jìn)行簡要介紹.

給定m=n=5，C和H的值與3.2.3中的值相同.假設(shè)松弛問題的解為：xm×n=[0 0 1 0 1； 1 1 1 1 0； 1 1 0 1 0； 1 0 0 0 0； 0 0 0 0 1].

可行化方法一：每個客戶應(yīng)選擇所需服務(wù)費用最小的設(shè)施為其服務(wù).以客戶1為例，其分別由設(shè)施2、3、4服務(wù)，比較三種服務(wù)費用，設(shè)施4所需費用最少，則客戶1由設(shè)施4服務(wù)，即x41=1，x21=x31=0.可行化后為x23=x24=x32=x41=x55=1，目標(biāo)函數(shù)值為7060.

可行化方法二：對于已開放的設(shè)施，將其服務(wù)的客戶分別由所有設(shè)施來服務(wù)，最后確認(rèn)由所需總費用最小的設(shè)施來服務(wù).以服務(wù)客戶3和4的設(shè)施2為例，計算客戶3和4由所有設(shè)施服務(wù)的總費用，min[1488+1235+130，1227+1783+199，1240+1831+139，121+1340+127，1375+112+103]=1588，設(shè)施4所需費用最少，則客戶3和4由設(shè)施4服務(wù).最后x32=x51=x53=x54=x55=1，目標(biāo)函數(shù)值為4179.

可行化方法三：每個客戶應(yīng)從已開放的設(shè)施集合中選擇所需服務(wù)費用最小的設(shè)施為其服務(wù).以客戶1為例，當(dāng)前開放設(shè)施3和5，比較客戶1分別由設(shè)施3和5服務(wù)的費用，設(shè)施5所需費用較少，則客戶1由設(shè)施5服務(wù).最后的結(jié)果為x32=x33=x35=x51=x54=1，目標(biāo)函數(shù)值為3143.

4.2拉格朗日蝙蝠算法設(shè)計

Step1：各參數(shù)初始化.蝙蝠種群大小K；最大迭代次數(shù)Nmax（N_iter=1）；音量A；脈沖發(fā)生率r；服務(wù)費用C；開放費用H，適應(yīng)度函數(shù)f（x）.初始化種群找出最優(yōu)蝙蝠x*.

Step2：根據(jù)漢明距離和精英保留策略更新蝙蝠位置.

Step3：產(chǎn)生隨機(jī)數(shù)rand1，判定rand1>ri是否成立，如果成立，則執(zhí)行局部搜索.

Step4： 產(chǎn)生隨機(jī)數(shù)rand2，判定rand2>Ai且f（xi）

Step5：求解松弛問題，得到原問題的下界值.

Step6：根據(jù)式（14-15）更新拉格朗日乘子.

Step7：下界可行化，得到原問題的可行解UB2，并按以下兩個標(biāo)準(zhǔn)更新π：如果連續(xù)5次沒有提高解，則π將減半；如果連續(xù)循環(huán)50次后，該值被重置為初始值.

Step8：將UB2與所有蝙蝠的適應(yīng)度值進(jìn)行比較，替換其中最大的適應(yīng)度值.

Step9：比較UB1和UB2，將兩者較小值保存為UB=min（UB1，UB2）.

Step10：判斷是否達(dá)到預(yù)設(shè)的最大迭代次數(shù)Nmax或時間限制3600s，若達(dá)到，則輸出最優(yōu)解UB，否則N_iter++，轉(zhuǎn)至Step2.

其中，去除Step5Step9為求解UFLP的蝙蝠算法.

5實驗與分析

為驗證拉格朗日蝙蝠算法的可行性及其求解性能，在操作系統(tǒng)為64位Windows7的實驗環(huán)境下，通過MATLAB編程對兩組標(biāo)準(zhǔn)算例求解測試，結(jié)果見圖1.第一組算例來自O(shè)R-Library，從中選取12個算例，m=16/25/50，n=50；第二組算例來自UFL基準(zhǔn)問題庫，同樣選取12個算例，m=n=250/500/750.同時，將其與改進(jìn)蟻群算法、基本蟻群算法、混合蟻群算法、DBA進(jìn)行比較.DBA和LR_BA的所有目標(biāo)函數(shù)值都是運行測試10次中的最好值.DBA和LR_BA的具體參數(shù)設(shè)置如下：K=20，Nmax=200，A=0.9，r=0.25，α=0.95，γ=0.05，t1=0.9，π=2，λj1=max cij，同時借鑒文獻(xiàn)[11]中p1和p2的參數(shù)設(shè)置，當(dāng)求解OR-Library的12個算例時，p1=0.3，p2=0.5；當(dāng)m=n=250時，p1和p2仍設(shè)定相同值；當(dāng)m=n=500時，p1=0.15，p2=0.3；當(dāng)m=n=750時，p1=0.1，p2=0.2.表1和2分別給出了OR-Library和UFL基準(zhǔn)問題庫中算例的計算結(jié)果.

在表1中，Na和Nb分別表示DBA和LR_BA的十次運行結(jié)果中優(yōu)于改進(jìn)蟻群算法的次數(shù).表中包括每個算例的已知最優(yōu)值、各算法的所得值、Gap值（Gap是各算法所得值與已知最優(yōu)值的相對偏移百分比）和平均運行時間.DBA和LR_BA的計算結(jié)果都遠(yuǎn)優(yōu)于改進(jìn)蟻群算法，其中DBA的計算結(jié)果非常接近最優(yōu)值，而LR_BA更可以得到全部算例的最優(yōu)值.但在運行時間方面，DBA和LR_BA都耗時較長，這是文章的不足之處.總的來說，DBA和LR_BA都可以用來解決UFLP，但LR_BA在處理時具有更明顯的可行性和有效性.

在表2中，Na和Nb分別表示DBA和LR_BA的十次運行結(jié)果中優(yōu)于混合蟻群算法的次數(shù).同時，為進(jìn)一步分析DBA和LR_BA在求解UFLP時的有效性和計算復(fù)雜度，選擇算例gs00250a給出求解的迭代過程圖，如圖1所示，雖然LR_BA中加入了拉格朗日松弛算法和三種可行化方法，致使其運行總時間遠(yuǎn)高于DBA，但是LR_BA的計算結(jié)果明顯優(yōu)于DBA，而且LR_BA的收斂性能明顯得到改善.就UFL基準(zhǔn)問題庫的12個算例而言，DBA和LR_BA在超過一半的算例中得到的計算結(jié)果都優(yōu)于混合蟻群算法.對比Na和Nb的值，發(fā)現(xiàn)LR_BA比DBA求解更為有效.LR_BA的最小Gap值和最大Gap值分別為0.790%和6.066%，其明顯小于由DBA計算得到的最小Gap值1.256%和最大Gap值11.573%，然而混合蟻群算法的最大Gap值為35.976%，基本蟻群算法的最大Gap值竟達(dá)123.572%.四種算法的平均Gap值分別為83.532%、16.707%、4.909% 和3.092%，從中可以看出DBA和LR_BA的計算結(jié)果非常接近已知最優(yōu)值，并且LR_BA的計算結(jié)果更優(yōu).但就運行時間而言，無論是求解小規(guī)模算例，還是大規(guī)模算例，DBA和LR_BA總會消耗較長時間.

6結(jié)論

為解決無容量設(shè)施選址問題，給出了一種結(jié)合BA和LR的混合優(yōu)化算法.針對BA求解離散問題的局限性，重新定義了BA的操作算子，引入了漢明距離和精英保留策略，提出了DBA.在此基礎(chǔ)上，通過對UFLP的具體特點進(jìn)行深入分析，得到其最優(yōu)解所具有的基本特征，結(jié)合LR，設(shè)計三種可行化方法，使兩種算法取長補(bǔ)短，進(jìn)一步引導(dǎo)LR_BA在有希望的解區(qū)域中強(qiáng)化搜索，明顯改進(jìn)了DBA的計算性能和求解結(jié)果.

通過求解24個UFLP算例，可以發(fā)現(xiàn)LR_BA無論對小規(guī)模算例還是大規(guī)模算例都具有較好的全局優(yōu)化性能，能夠有效求解UFLP，驗證了該算法在求解質(zhì)量上的優(yōu)勢，但也存在一定的缺陷，所需消耗的運行時間相對于其他算法來說較長，很難在合理的計算時間內(nèi)求得問題的滿意解，需要對此作進(jìn)一步的改進(jìn)，這是作者今后的研究方向.

參考文獻(xiàn)

[1]吳燁. 物流配送網(wǎng)絡(luò)選址的模糊數(shù)學(xué)模型及其算法[J]. 經(jīng)濟(jì)數(shù)學(xué)， 2008， 25（3）：277-282.

[2]堵丁柱. 近似算法的設(shè)計與分析[M]. 高等教育出版社， 2011.

[3]KOLE A， CHAKRABARTI P， BHATTACHARYYA S. An Ant Colony Optimization Algorithm for Uncapacitated Facility Location Problem[J]. Artificial Intelligence & Applications， 2014， 1（1）：55-61.

[4]DAMGACIOGLU H， DINLER D， EVIN O， et al. A genetic algorithm for the uncapacitated single allocation planar hub location problem [J]. Computers & Operations Research， 2015， 62（C）：224-236.

[5]李倩， 張惠珍， Cesar Beltran-Royo. 求解無容量設(shè)施選址問題的混合蟻群算法[J]. 上海理工大學(xué)學(xué)報， 2016， 38（4）：367-372.

[6]YANG X S. A new metaheuristic bat-inspired algorithm[C]//Nature Inspired Cooperative Strategies for Optimization （NISCO 2010）（Eds.J R Gonzalez et al.）， SCI 284， 2010： 65-74.

[7]ADARSH B R， RAGHUNATHAN T， JAYABARATHI T， et al. Economic dispatch using chaotic bat algorithm [J]. Energy， 2016（96）：666-675.

[8]劉春苗， 張惠珍， 馬祥麗. 求解車輛路徑問題的離散蝙蝠算法[J]. 經(jīng)濟(jì)數(shù)學(xué)， 2016， 33（4）：91-95.

[9]FISHER M L. Optimal Solution of Scheduling Problems Using Lagrange Multipliers： Part II[M]// Symposium on the Theory of Scheduling and Its Applications. Springer Berlin Heidelberg， 1973：1114-1127.

[10]高珊， 馬良， 張惠珍. 函數(shù)優(yōu)化的小生境蝙蝠算法[J]. 數(shù)學(xué)的實踐與認(rèn)識， 2014， 44（15）：253-260.

[11]李岳佳. 基于遺傳算法的設(shè)施選址問題算法研究[D]. 濟(jì)南：山東科技大學(xué)信息科學(xué)與工程學(xué)院， 2012.