设施选址问题是运筹学中的经典问题之一，在生产、生活及物流等方面都有着非常广泛的应用，如工厂、仓库、急救中心、消防站及物流中心的选址等。而P中位问题（P-median problem，PMP）就是一种很常见的设施选址问题，该问题最初由Hakimi^[1-2]在1964年提出的，目的是要在给定没有容量限制的m个候选设施集中选择 $ p\left( {p < m} \right)$ 个作为开放设施，使得所有顾客到这些中位点的距离之和达到最小。

目前，用于求解此问题的方法多种多样，如动态规划方法^[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 = \left\{ {1,2, \cdots ,m} \right\}$ ，客户集 $N = $ $ \left\{ {1,2, \cdots ,n} \right\}$ ，对于任意给定的i∈M，j∈N，d_ij表示设施i与客户j之间的距离，所有的候选设施都没有容量限制，要求在m个候选位置点中选择p个位置作为开放设施（中位点），并且要求每个客户必须选择一个且只能选择一个开放设施来满足其需求，同时使总距离最小。PMP的数学模型可被描述为

$ \min \;z = \mathop \sum \limits_{i = 1}^m \mathop \sum \limits_{j = 1}^n {d_{ij}}{x_{ij}} $

(1)

$ \mathop \sum \limits_{i = 1}^m {x_{ij}} = 1,\quad \forall j \in N $

(2)

$ {x_{ij}} \leqslant {y_i},\quad\forall i \in M,j \in N $

(3)

$ \mathop \sum \limits_{i = 1}^m {y_i} = p,\quad\forall i \in M$

(4)

$ {x_{ij}} \in \left\{ {0,1} \right\},\forall i \in M,j \in N $

(5)

$ {y_i} \in \left\{ {0,1} \right\},\forall i \in M $

(6)

式中：z为目标函数，即总距离；p为允许开放的设施数目；x_ij表示客户j是否由设施i服务，如果客户j由设施i服务，则x_ij=1，否则，x_ij=0；y_i表示设施i是否开放，如果设施i开放，则y_i=1，否则，y_i=0。

2 基本蝙蝠算法

蝙蝠是一种神奇的动物，它靠自身的回声定位系统来探测猎物，避免障碍物，在黑暗中找到它们的栖息地，而BA是一种基于蝙蝠的回声定位系统而产生的智能启发式算法，将蝙蝠映射为问题空间中的可行解，将搜索和优化的过程模拟成蝙蝠寻找猎物的过程，利用适应度函数值的大小来衡量蝙蝠所处位置的优劣，将蝙蝠的优胜劣汰过程类比为搜索和优化过程中用好的可行解替代较差的可行解的迭代过程。

2.1 基本蝙蝠的搜索过程

BA的搜索过程实质上是通过频率的变化实现对蝙蝠速度的控制，从而更新蝙蝠的位置。假设在d维空间中，蝙蝠i在t时刻下的速度和位置分别为v_i^t和x_i^t，则该蝙蝠通过改变频率f_i来调整自身位置与当前所处位置最好的蝙蝠之间的距离，从而更新自身新位置x_i^t+1，向当前求解的最好解靠拢，则蝙蝠i在t+1时刻下的速度和位置更新公式定义为

$ {f_i} = {f_{\min}} + \left( {{f_{\max}} - {f_{\min}}} \right)\beta $

(7)

$ v_i^{t + 1} = v_i^t + \left( {x_i^t - {x^*}} \right){f_i} $

(8)

$ x_i^{t + 1} = x_i^t + v_i^{t + 1} $

(9)

式中： $\;\beta\in [0,\;1] $ 是服从均匀分布的随机数；f_min和f_max为频率的最小值和最大值；v_i^t和v_i^t+1为t和t+1时刻蝙蝠i的速度；x^*表示当前求得的最好解；x_i^t和x_i^t+1分别为t和t+1时刻蝙蝠i的位置。

2.2 基本蝙蝠的优化过程

蝙蝠的优化过程即局部搜索过程，在当前求解的最好解集中选择一个解x_old，那么，每只蝙蝠在该最好解附近按照随机游走法则产生局部新解x_new。

$ {x_{\rm new}} = {x_{\rm old}} + {\varepsilon} {A^t} $

(10)

式中： ${\varepsilon} $ 是属于[−1, 1]的p维随机向量；A^t表示t时刻当前代蝙蝠种群音量的平均值。

当蝙蝠找到猎物时，音量就会降低，同时脉冲发生率逐渐提高。具体按照式（11）和式（12）对音量和脉冲发生率进行更新。

$ A_i^{t + 1} = \alpha A_i^t $

(11)

$ r_i^{t + 1} = r_i^0\left[ {1 - {\rm exp}\left( { - \gamma t} \right)} \right] $

(12)

式中：α和γ是常量；A_i^t和A_i^t+1分别为t和t+1时刻蝙蝠i的音量；r_i⁰为初始脉冲发生率；r_i^t+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对应的编码为 $x_i=(5\;2\;3) $ ，即设施2，3，5为中位点。

3.2 相关定义

考虑到BA求解离散问题的局限性，为了将BA成功运用到PMP的求解中，并且有较好的运算效率，根据PMP模型的具体特征，设计了针对求解该问题的相关操作算子。

定义1 　速度的更新。假设t时刻下蝙蝠i的位置为 $ {x}_{i}^{t} = ( {x_{i1}^{t},{x}_{i2}^t, \cdots ,{x}_{ip}^t} )$ ，当前所处位置最好的蝙蝠的位置为 $ {{x}}^{*} = ( {{x}_{1}^{*},{x}_{2}^{*}, \cdots ,{x}_{p}^{*}} )$ ，则蝙蝠i的速度v_i^t+1=x^*−x_i^t。其中，v_i^t+1是由所有属于x^*且不属于x_i^t的元素组合而成的集合。当x^*和x_i^t相同时，则随机产生一个p维向量代替其中的x^*，进行相同的速度更新操作。

定义2 　位置的更新。蝙蝠i的位置为x_i^t+1= $Fea(x_i^t+x_i^{t+1}) $ ，可分为两部分。第一部分： $x_i^{t+1'}= $ $x_i^t+v_i^{t+1} $ ，其中， $x_i^{t+1'} $ 是由所有属于 $ x_i^{t+1'}$ 或 $v_i^{t+1} $ 的元素组合而成的集合；第二部分： $x_i^{t+1}=Fea(x_i^{t+1'}) $ ，其中，Fea（）函数是可行化函数。由于集合 $x_i^{t+1'} $ 的长度超过p，因此，需根据适应度值的大小逐次递减该向量中的一个元素，直至集合 $x_i^{t+1'} $ 的长度等于p，但同时属于x^*和x_i^t的元素应始终保持在 $x_i^{t+1'} $ 中，不得进行删减操作。

定义3 　局部搜索。基于遗传算法中交叉的思想，随机选择一个已开放的设施O₁和另一个未开放的设施C₁，交换O₁和C₁的位置，若产生的新解的目标函数值优于当前求得的最好解的目标函数值，则更新当前求得的最好解；否则，保持不变。

例如，假设存在一个m=5，p=3的PMP，距离矩阵

$ {{D }}=\left[ {\begin{array}{*{20}{c}} {{0}}&{{{424}}}&{{{569}}}&{{{254}}}&{{{275}}}\\ {{{424}}}&{{0}}&{{{391}}}&{{{452}}}&{{{298}}}\\ {{{569}}}&{{{391}}}&{{0}}&{{{167}}}&{{{528}}}\\ {{{254}}}&{{{452}}}&{{{167}}}&{{0}}&{{{456}}}\\ {{{275}}}&{{{298}}}&{{{528}}}&{{{456}}}&{{0}} \end{array}} \right] $

蝙蝠i在当前时刻t下的位置为 $ x_i^t=(3\;1\;4)$ ，当前所处位置最好的蝙蝠的位置为 $x^*=(1\;4\;5) $ ，则蝙蝠i的速度为 $v_i^{t+1}={(5)} $ 。第一部分所得到的 $x_i^{t+1'} = $ $(3\;1\;4\;5)$ ，再按照第二部分进行可行化，因为，元素1和4同时属于 $ x^*$ 和 $x_i^t $ ，则应始终保持在 $x_i^{t+1'} $ 中，不能进行删减，可得到2种不同的组合 $\left[ \begin{array}{*{35}{l}} 1 & 4 & 5 \\ 3 & 1 & 4 \\\end{array} \right]$ ，比较相应的适应度值 $ \left[ {\begin{array}{*{20}{c}} {{{0}} + {{298}} + {{167}} + {{0}} + {{0}}}\\ {{{0}} + {{391}} + {{0}} + {{0}} + {{275}}} \end{array}} \right] = $ $\left[ {\begin{array}{*{20}{c}} {{{465}}}\\ {{{666}}} \end{array}} \right] $ ，其中，适应度值465较小，则 $x_i^{t+1}= $ $(1\;4\;5) $ 。对当前所处位置最好的蝙蝠 $x^*=(1\;4\;5) $ 进行局部搜索时，随机选择一个已开放的设施O₁=1，一个未开放的设施C₁=3，则产生的局部新解 $x_{\rm new}=(3\;4\;5) $ ，计算该新解的适应度值为552，大于当前所处位置最好的蝙蝠的适应度值465，则不更新当前的最好蝙蝠。若该适应度值小于当前所处位置最好的蝙蝠的适应度值时，则进行更新。

3.3 混合蝙蝠算法的设计

混合蝙蝠算法的步骤：

步骤1 　各参数初始化。蝙蝠种群大小K；最大迭代次数 $ N_{\max}\left( {N\_iter = 1} \right)$ ；音量A；脉冲发生率r；距离矩阵D；适应度函数 $ f\left( x \right)$ 。

步骤2 　随机初始化种群，找出当前所处位置最好的蝙蝠 $x^*$ 。

步骤3 　按照定义1和定义2更新当前蝙蝠的速度v_i^t+1和位置x_i^t+1。

步骤4 　产生随机数rand 1，判定rand 1>r_i是否成立，如果成立，则根据定义3执行局部搜索。

步骤5 　产生随机数rand 2，判定rand 2>A_i且 $ f\left( {{x_i}} \right) < f\left( {{x^*}} \right)$ 是否同时成立，如果成立，则接受这个新解，并更新A，r和当前求得的最好解。

步骤6 　判断是否达到预设的最大迭代次数N_max或时间限制3 600 s，若达到，则输出当前求得的最好解；否则， $ N\_iter + + $ ，转至步骤3。

4 数值实验

为了验证混合蝙蝠算法的求解性能，在Intel（R）Core（TM）CPU@3.60 GHz处理器、7.87 GB内存、操作系统为64位Windows 10的实验环境下，应用Matlab R2016a编程实现了求解PMP的混合蝙蝠算法，对运筹学数据库OR-Library中的40个PMP算例进行仿真实验，并与文献[13]中的遗传算法（genetic algorithm，GA）和文献[14]中的粒子群算法（novel discrete particle swarm optimization，NDPSO）进行对比。求解这些算例的混合蝙蝠算法参数设置情况为：最大迭代次数N_max=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在迭代到265 s左右、图2到1 400 s左右时，2个算例的计算结果还在发生改变。但是，从运行总时间来看，HBA相对于NDPSO而言，需要耗用较多的时间，这是本文的不足之处。总的来说，HBA在求解这些规模不同的PMP算例时能够求得较好的实验结果，与已知最优解的误差较小，且绝大多数算例结果都可以达到已知最优解，故进一步验证了蝙蝠算法在求解PMP上的可行性与有效性。

5 结　论

针对P中位问题的特点以及蝙蝠算法的寻优机制，重新定义了原有蝙蝠算法的操作算子，并提出了一种可行化函数。同时在局部搜索部分引入了遗传算法的交叉思想，设计了一种针对求解该问题的局部搜索方式，提出了求解该问题的混合蝙蝠算法。

实验部分求解多个测试算例，并将混合蝙蝠算法的计算结果与文献中遗传算法和粒子群算法的计算结果进行对比，测试结果表明，混合蝙蝠算法在确保种群多样性的基础上，具有较好的全局优化性能，能够有效求解P中位问题，验证了该算法在求解质量上的优势。但是，与其他改进的智能算法相比，本文提出的混合蝙蝠算法在求解速度上还存在明显的不足，这将是进一步的研究方向。