基于禁忌搜索的带容量约束p-中值选址问题求解与Python实现

1. 什么是带容量约束的p-中值选址问题

我第一次接触选址问题是在一个物流优化项目中。客户需要在某个区域新建3个配送中心，要求既满足周边超市的配送需求，又要控制运输成本。这其实就是典型的p-中值问题。

p-中值问题属于设施选址问题的经典类型，它的核心目标是：从多个备选位置中选择p个设施点，使得所有需求点到其最近设施点的加权距离之和最小。这里的"加权"很重要，因为不同需求点的需求量（比如超市的订单量）可能差异很大。

带容量约束的意思是每个设施都有服务上限。比如一个仓库最多只能处理5000件商品/天的配送量，超过这个限制就会导致运营压力过大。在实际项目中，忽略容量约束往往会得到无法落地的"纸上谈兵"方案。

举个例子：某连锁便利店要在城市开设5个配送中心，已知：

• 20个候选仓库位置（租金、面积不同）
• 100家门店的分布位置和每日需求量
• 每个配送中心的最大处理能力 目标是选择5个位置建仓，使得总运输成本最低，同时不超出任何仓库的处理能力。

2. 问题建模与数学表达

2.1 基本符号定义

先明确问题中的关键要素：

• 需求点集合：用I表示，比如I={1,2,...,100}代表100家门店
• 备选设施集合：用J表示，比如J={1,2,...,20}代表20个候选仓库
• 需求点的需求量：d_i表示第i个需求点的需求量
• 设施的服务能力：c_j表示第j个设施的最大服务能力
• 距离或运输成本：w_ij表示从设施j到需求点i的单位运输成本

2.2 数学模型构建

建立数学模型需要定义决策变量：

• y_j：是否选择第j个位置建仓（1选择，0不选）
• x_ij：需求点i是否由设施j服务（1是，0否）

目标函数和约束条件如下：

minimize ∑(i∈I)∑(j∈J) w_ij * d_i * x_ij subject to: ∑(j∈J) y_j = p # 选择p个设施 ∑(j∈J) x_ij = 1, ∀i∈I # 每个需求点只由一个设施服务 x_ij ≤ y_j, ∀i∈I, ∀j∈J # 只有被选中的设施才能服务 ∑(i∈I) d_i * x_ij ≤ c_j, ∀j∈J # 容量约束 x_ij, y_j ∈ {0,1} # 0-1变量

这个模型看起来简单，但当I和J规模较大时（比如都有上百个点），求解会变得非常困难。我在第一次实现时就遇到了性能问题。

3. 禁忌搜索算法设计

3.1 算法核心思想

禁忌搜索(TS)是一种元启发式算法，它的特点是：

1. 使用禁忌表记录近期搜索历史，避免重复搜索
2. 允许暂时接受较差的解，以跳出局部最优
3. 通过邻域搜索机制探索解空间

在选址问题中，TS的表现通常优于遗传算法和模拟退火，特别是在处理容量约束时。根据我的项目经验，TS的求解质量比传统启发式算法高15-20%。

3.2 编码方案设计

编码直接影响搜索效率。我推荐使用混合编码方案：

• 前|J|位：0/1表示是否选择该设施（y_j）
• 后|I|位：自然数表示服务该需求点的设施编号

例如有5个候选设施和8个需求点，选3个设施：

[1,0,1,0,1, 1,3,2,1,3,2,1,3,1] ↑设施选择部分 ↑需求点分配部分

这种编码既保留了设施选择信息，也明确了需求分配。

3.3 邻域搜索策略

有效的邻域操作能显著提升搜索效率。我常用以下两种组合：

1. 成对交换(Swap)：随机选择两个设施交换其状态（选变不选/不选变选）

# 示例：交换第2和第4个设施 原解：[1,0,1,0,1, 1,3,2,1,3,2,1,3,1] 新解：[1,1,1,1,1, 1,3,2,1,3,2,1,3,1]

2. 单点变异(Mutation)：随机改变某个需求点的服务设施

# 示例：将第6个需求点的服务设施从3改为2 原解：[1,0,1,0,1, 1,3,2,1,3,2,1,3,1] 新解：[1,0,1,0,1, 1,2,2,1,3,2,1,3,1]

实际项目中，我会设置变异概率（通常0.1-0.3），避免过度扰动。

3.4 禁忌表管理

我通常维护两个禁忌表：

1. 全局禁忌表：存储最近k代的最优解（禁忌长度k=50-100）
2. 局部禁忌表：存储当前邻域搜索中已尝试的解

这种双重机制能有效平衡全局搜索和局部搜索。禁忌表实现示例：

tabu_list = { '[1,0,1,0,1,1,3,2,1,3,2,1,3,1]': 5, # 禁忌剩余代数 '[1,1,1,0,0,1,2,2,1,3,2,1,3,1]': 3 }

4. Python完整实现

4.1 数据准备与初始化

首先定义问题数据。这里用我之前项目的一个简化案例：

import numpy as np import pandas as pd from matplotlib import pyplot as plt # 需求点坐标和需求量 demand_points = np.array([ [10, 20], [30, 40], [50, 60], [70, 80], [90, 10], [20, 30], [40, 50], [60, 70] ]) demand = [5, 8, 6, 9, 7, 4, 6, 5] # 备选设施坐标和能力 facility_points = np.array([ [15, 25], [35, 45], [55, 65], [75, 85], [25, 35] ]) capacity = [20, 25, 30, 15, 18] # 计算距离矩阵 def calc_distance_matrix(d_points, f_points): return pd.DataFrame( [[np.linalg.norm(d-f) for f in f_points] for d in d_points], index=range(len(d_points)), columns=range(len(f_points)) ) dist_matrix = calc_distance_matrix(demand_points, facility_points)

4.2 禁忌搜索核心算法

完整实现禁忌搜索算法：

class TabuSearch: def __init__(self, dist_matrix, demand, capacity, p): self.dist_matrix = dist_matrix self.demand = demand self.capacity = capacity self.p = p self.n_demand = len(demand) self.n_facility = len(capacity) def initialize(self): # 随机选择p个设施 selected = np.random.choice( self.n_facility, self.p, replace=False ) chrom = [0]*self.n_facility for s in selected: chrom[s] = 1 # 随机分配需求点 for _ in range(self.n_demand): chrom.append(np.random.choice(selected)+1) return chrom def evaluate(self, chrom): selected = [i for i in range(self.n_facility) if chrom[i]==1] demand_assigned = [0]*self.n_facility total_cost = 0 for i in range(self.n_demand): fac = chrom[self.n_facility + i] - 1 if fac not in selected: return float('inf') # 无效解 # 检查容量 demand_assigned[fac] += self.demand[i] if demand_assigned[fac] > self.capacity[fac]: return float('inf') # 容量超限 total_cost += self.dist_matrix.iloc[i, fac] * self.demand[i] return total_cost def neighborhood_search(self, chrom, tabu_list, max_trials=100): best_neighbor = None best_value = float('inf') trials = 0 while trials < max_trials: neighbor = chrom.copy() # 随机选择邻域操作 if np.random.random() < 0.7: # 70%概率选择交换 f1, f2 = np.random.choice(self.n_facility, 2, replace=False) neighbor[f1], neighbor[f2] = neighbor[f2], neighbor[f1] else: # 30%概率选择变异 d = np.random.randint(self.n_demand) new_fac = np.random.choice( [i for i in range(self.n_facility) if neighbor[i]==1] ) neighbor[self.n_facility + d] = new_fac + 1 # 检查禁忌状态 str_rep = str(neighbor) if str_rep not in tabu_list and self.evaluate(neighbor) < best_value: best_neighbor = neighbor best_value = self.evaluate(neighbor) trials += 1 return best_neighbor, best_value def solve(self, max_iter=500, tabu_tenure=50): current = self.initialize() best = current.copy() best_value = self.evaluate(best) tabu_list = {} # 禁忌表：解 -> 剩余禁忌代数 history = [] for _ in range(max_iter): # 邻域搜索 neighbor, neighbor_value = self.neighborhood_search( current, tabu_list ) if neighbor is None: # 未找到可行邻域解 continue # 更新当前解 current = neighbor current_value = neighbor_value # 更新最优解 if current_value < best_value: best = current.copy() best_value = current_value # 更新禁忌表 tabu_list[str(current)] = tabu_tenure tabu_list = {k: v-1 for k, v in tabu_list.items() if v > 1} history.append(best_value) return best, best_value, history

4.3 结果可视化

绘制优化过程和最终选址方案：

def visualize_result(ts, best_solution): # 提取选中的设施 selected = [i for i in range(ts.n_facility) if best_solution[i]==1] # 绘制需求点和设施 plt.figure(figsize=(10,6)) plt.scatter( demand_points[:,0], demand_points[:,1], c='blue', label='Demand Points' ) plt.scatter( facility_points[:,0], facility_points[:,1], c='red', marker='s', label='Facility Candidates' ) plt.scatter( facility_points[selected,0], facility_points[selected,1], c='green', marker='*', s=200, label='Selected Facilities' ) # 绘制分配关系 for i in range(ts.n_demand): fac = best_solution[ts.n_facility + i] - 1 plt.plot( [demand_points[i,0], facility_points[fac,0]], [demand_points[i,1], facility_points[fac,1]], 'gray', alpha=0.3 ) plt.legend() plt.title('Facility Location Allocation') plt.show() # 运行算法 ts = TabuSearch(dist_matrix, demand, capacity, p=3) best_sol, best_val, history = ts.solve() print(f"Best solution value: {best_val}") visualize_result(ts, best_sol) # 绘制收敛曲线 plt.plot(history) plt.title('Convergence Curve') plt.xlabel('Iteration') plt.ylabel('Total Cost') plt.show()

5. 实战技巧与优化建议

5.1 参数调优经验

经过多个项目实践，我总结出以下参数设置经验：

• 禁忌长度：通常设为问题规模的10-20%。比如100次迭代，禁忌长度10-20
• 邻域大小：建议每次搜索生成50-100个邻域解
• 变异概率：0.1-0.3效果较好，过高会导致搜索不稳定
• 迭代次数：至少500次，复杂问题需要1000-5000次

可以通过网格搜索寻找最优参数组合：

param_grid = { 'tabu_tenure': [10, 20, 50], 'max_iter': [300, 500, 1000], 'p': [2, 3, 4] }

5.2 处理大规模问题

当问题规模较大时（如超过100个需求点），可以：

1. 分区域求解：先聚类将需求点分组，再分别求解
2. 并行计算：使用多进程同时评估多个邻域解
3. 记忆化：缓存已评估的解，避免重复计算

from multiprocessing import Pool def parallel_evaluate(chroms): with Pool() as p: return p.map(evaluate, chroms)

5.3 常见问题排查

遇到过的一些典型问题及解决方案：

1. 算法陷入局部最优：

   • 增加禁忌长度
   • 引入重启机制（每N代重置当前解）
   • 结合模拟退火的接受准则

2. 运行时间过长：

   • 优化距离矩阵计算（使用KDTree）
   • 实现快速评估方法（增量式计算）
   • 设置早期终止条件

3. 无法满足容量约束：

   • 惩罚函数法：在目标函数中加入惩罚项
   • 修复算子：专门处理不可行解的修复过程

# 惩罚函数示例 def evaluate_with_penalty(chrom): base_cost = evaluate(chrom) if base_cost == float('inf'): # 计算约束违反程度 violation = calc_capacity_violation(chrom) return base_cost + 1e6 * violation # 大惩罚系数 return base_cost