论文拾萃 | PISTS算法求解obnoxious p-median problem (附Python代码）

前 言

各位读者大家好，小编最近读到关于选址问题的一篇文章，读完感觉深有启发，特此来与大家分享~另外，该篇文章的作者也有将算法的代码进行公开，小编稍后也会分享给大家。

Chang J , Wang L , Hao J K , et al. Parallel iterative solution-based tabu search for the obnoxious p-median problem[J]. Computers & operations research, 2021(Mar.):127.

首先我们需要先搞清楚什么是obnoxious p-median problem，在此之前我们先了解p-median problem的概念。

P-中位问题（p-median problem）

P-中位问题（也叫P-中值问题）是研究在一个给定数量和位置的需求集合和一个给数量和候选位置的设施集合的前提下，分别为P个设施找到合适的位置，并指派每个需求点到一个特定的设施，使之达到在设施与需求点之间的运输费用最低。Hakimi提出该问题之后给出了 P-中位问题的 Hakimi 特性，他证明了 P-中位问题的服务站候选点限制在网络节点上时至少有一个最优解是与不对选址点限制时的最优解是一致的，所以将网络连续选址的 P-中位问题简化到离散选址问题不会影响到目标函数的最优值。

这篇文章虽然没有具体介绍该问题的数学模型，但是小编为了更加深入地理解，在其他资料中找到了p-median problems的数学模型：

其中，

N：在研究对象中的 n 个客户，N=（1,2,...,n）

M：在研究对象中的 m 个拟建设施的候选地点，M=（1,2,...,m）

：第个 i 客户的需求量

：从地点 i 到 j 的单位运输费用

p ：可以建立的设施总数

：如果在建立设施，则= 1；否则为0

：如果客户在，由设施来提供服务时，则其为1；否则为0

约束1：保证每个客户（需求点）只有一个设施来提供相应的服务

约束2：限制总的设施数目为 p 个

约束3：有效得保证没有设施的地点不会有客户对应

另外，除了p-median problems，选址研究中的典型问题，如Weber（韦伯） 问题、覆盖问题、中心问题（p-center problems）、多目标选址、竞争选址、选址-分配、选址-路线等，都是引起广泛关注和深入研究的热点课题，研究的也较为成熟，感兴趣的读者都可以去了解。

obnoxious p-median problem

不受欢迎的p中值问题与p-median problem相似，只是不再是为了最小化距离或者运输成本，而是为了确定一组设施的选址，使每个客户端与开放设施之间的距离之和最大化；即在这个问题中，需求点试图避开和远离设施。该问题的应用场景也很常见，典型的应用包括核反应堆、垃圾场等的最佳位置。在疫情形势下，需要建立大量的隔离场所供患者入院。考虑到潜在的感染风险，公共设施应远离居民区。

因此，obnoxious p-median problem的目标函数应该稍作修改。设 N 表示一组需求点集合，M 表示设施的选址集合，

为需求点i到选址j之间的距离。目标是找到 p 个设施的选址( p 已给定)，即一个子集

满足：

其中，

表示给定需求点 i 与所有设施之间的距离，定义为 i 与每个设施选址 j 的所有距离中的最小值。目标函数 f 最大化所有需求点到设施选址的距离之和。

Parallel iterative solution-based

tabu search（PISTS）

从名字我们就可以看出，该算法的特点有三个：“parallel”、“iterative”和具有创新性的“tabu search”，因此下面我们将从这三块来讲这个算法逐步拆分，便于大家理解。

首先是“并行”这一特点。并行迭代基于解的禁忌搜索(PISTS)的总体方案如算法1所示，其中ISTS函数表示的是Iterative solution-based tabu search，将在后面具体解说，因此这个伪代码主要体现的是并行过程。

该文章通过使用Python的多处理模块来利用多个处理器以实现“并行”。PISTS算法首先创建一个类Pool的对象ProcPool来存储

个进程。通过引入共享内存机制，使所有进程都可以访问禁忌列表，即先前某一进程得到并放入禁忌列表的解将不会被其他进程再次访问。在所有进程完成后，运行ProcPool来执行处理器上的每个进程并返回所有进程中的最优解

。

初步实验表明，在得到相同解的情况下，与典型的进程间不进行信息交互的并行实现方法相比，这种并行实现方法所需的计算时间更少，特别是在求解大型实例时。

然后是“迭代”。迭代基于解的禁忌搜索(ISTS)的伪代码如下图所示，对应于上图中的ISTS方法。下图中的STS方法表示的是solution-based tabu search，将在后面具体解说，因此这个伪代码主要体现的是迭代过程。

从伪代码中我们可以看出，ISTS在基于解的禁忌搜索和与VNS算法中的shaking相似的扰动过程之间交替使用。我们先介绍领域的定义，再对算法进行讲解。

邻域

：解

的邻域

由所有到解

距离为 k 的解组成，它被表述为

其中

表示解空间，那么

是什么意思呢？

定义

与

之间的距离为

其中|·|表示集合的基数。通俗点来说

与

之间的距离表示它们之间有多少元素不同，

需要从J\

中交换多少个元素才能得到

。

我们弄清楚领域的定义后，算法的过程也就易于理解了：首先使用基于解的禁忌搜索（STS）以达到局部最优解

。然后，算法2的外部“while”循环搜索邻域，直到达到给定的时间限制T。对于每个特定邻域

，内部的“while”循环迭代执行扰动过程以生成一个新的解

，并基于解的禁忌搜索过程通过改进

的质量来达到一个新的局部最优

。如果

得到了改进，就返回到邻域

；如果

在连续h次迭代中没有都得到改进，就搜索下一个领域

。

在扰动过程中，采用了一种随机策略，执行 k 次随机删除-添加移动。因此，k既是领域的下标，也是扰动强度。如果最优解连续h次迭代中没有得到改进，扰动过程动态调整扰动强度 k，系统地增加删除-添加移动的次数；当搜索改善时，则减小到最小次数。

当探索特定邻域

时，因为对于任何解，最多有

个元素可以删除和添加，所以当k超过上限时会被重置为2，以保证扰动解的可行性。

第三点是基于解的禁忌搜索算法。基于解的禁忌搜索算法的特点之一是将禁忌策略进行了创新，普通的禁忌策略通常会在一段时间内禁止之前执行过的移动，而基于解的禁忌策略则是直接禁忌之前得到过的解被重新访问。传统的禁忌策略需要设置禁忌长度(Tabu tenure)，这是一个参数，其调优通常会影响搜索性能。而基于解的禁忌搜索则消除了调优禁忌长度的需要。另外，在禁忌列表中记录搜索过程中得到的所有解需要太多的内存，因此不考虑保存每个解，而是保存它的哈希值。这样，要验证解是否是合格的候选解，只需检查该解的哈希值是否在禁忌列表中就足够了。

此外，为了解决obnoxious p-median problem，传统的交换邻域需要

次领域移动，因此在每次迭代中确定最佳交换移动是非常耗时的。为了克服这个问题，作者使用了删除-添加复合移动，首先从当前解

执行最佳删除移动，生成部分解

，然后对

执行最佳添加移动，以达到一个新的可行解决方案

。因此，在该算法中有两个禁忌列表TabuListDLT和TabuListADD。每次选择的删除和添加移动都应满足相应解

和

不在禁忌列表中。

算法3给出了基于解的禁忌搜索的伪代码。

实验证明，通过以上的改进措施，算法的速度得到了显著提升。

核心代码分享

在这里小编为大家搜集了文献中使用的实例用于测试代码。这些实例是根据OR-Library instances (Beasley, J. J Oper Res Soc (1990) 41: 1069. doi:10.1057/jors.1990.166)生成的。

算例OpM_LIB_2016\pmed17.txt.table.p25.A.txt的部分内容如图：

经小编注释的核心代码如下：

# =============================================================================
# neighbor search 
ISTS的内层while，对每个特定的领域
调用了tabu search 和 _perturb
# =============================================================================
def specific_neighbor_search(self, best_sol, best_obj, HTB, HTC):
    count = 0
    while True:
        if self._exit:
            break
        # perturbation
        res = set(range(self.dim)) - set(best_sol)
        min_l = min(self.level, self.dim - self.level)
        # 如果移动次数超过上限
        if self._step_num > min_l:
            self._move_flag = True
            break
        else:
            temp_sol = _perturb(self, res, best_sol)
            temp_obj = objective(self.data, temp_sol)
        # tabu search
        stag_sol, stag_obj = tabu_search(self, temp_sol, temp_obj, HTB, HTC)
        # condition: better solution found
        if stag_obj > best_obj:
            best_sol = stag_sol
            best_obj = stag_obj
            self._move_flag = True
            count = 0
            break
        else:
            count += 1
        # condition: reach range
        if count >= self.neighbor_depth:
            self._move_flag = False
            break
        pass
    return best_sol, best_obj

# =============================================================================
# tabu search（STS）
# _map_c2b_wt表示删除操作
# _map_b2c_wt表示添加操作
# =============================================================================
def tabu_search(self, curr_sol, curr_obj, HTB, HTC):
    best_sol, best_obj = curr_sol, curr_obj
    count = 0
    while True:
        # global exit condition
        if self.local_step >= 20 * self.level:
            self._exit = True
            break
        # update local step
        self.local_step += 1
        # delete-add local search
        back_sol, back_obj = _map_c2b_wt(self, curr_sol, HTB)
        if not back_obj:
            break
        curr_sol, curr_obj = _map_b2c_wt(self, back_sol, HTC)
        if not curr_obj:
            break
        # update tabu optimum
        if curr_obj > best_obj:
            best_sol = curr_sol
            best_obj = curr_obj
            count = 0
        else:
            count += 1
        # condition: search depth surpassed
        if count >= self.tabu_depth:
           break
        pass
    # update tabu_step
    self.tabu_step += 1
    return best_sol, best_obj

# =============================================================================
# with tabu 删除-添加操作的部分核心代码
# ============================================================================
def tabu_neighbor_current2back(curr_sol, HTB):
    back_sols = []
    for i in curr_sol:
        back_sol = list(curr_sol)
        back_sol.remove(i)
        back_sol = tuple(sorted(back_sol))
        if hash(back_sol) not in HTB: #如果删除操作后得到的解的哈希值不在禁忌列表中
            back_sols.append(back_sol)
        else:
            HTB[hash(back_sol)] += 1
        pass
    return back_sols


def tabu_neighbor_back2current(n_dim, back_sol, HTC):
    res = tuple(set(range(n_dim)) - set(back_sol))
    curr_sols = []
    for i in res:
        curr_sol = list(back_sol)
        curr_sol.append(i)
        curr_sol = tuple(sorted(curr_sol))
        if hash(curr_sol) not in HTC:#如果添加操作后得到的解的哈希值不在禁忌列表中
            curr_sols.append(curr_sol)
        else:
            HTC[hash(curr_sol)] += 1
        pass
    return curr_sols

欲获取源代码，请移步留言区！

-The End-

文案&代码注释：胡心瑶（华中科技大学管理学院本科四年级）

指导老师：秦虎老师 （华中科技大学管理学院）

排版：程欣悦

如对文中内容有疑问，欢迎交流。PS：部分资料来自网络。