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)是一种元启发式算法它的特点是使用禁忌表记录近期搜索历史避免重复搜索允许暂时接受较差的解以跳出局部最优通过邻域搜索机制探索解空间在选址问题中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 邻域搜索策略有效的邻域操作能显著提升搜索效率。我常用以下两种组合成对交换(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]单点变异(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 禁忌表管理我通常维护两个禁忌表全局禁忌表存储最近k代的最优解禁忌长度k50-100局部禁忌表存储当前邻域搜索中已尝试的解这种双重机制能有效平衡全局搜索和局部搜索。禁忌表实现示例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], indexrange(len(d_points)), columnsrange(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, replaceFalse ) 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_trials100): 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, replaceFalse) 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_iter500, tabu_tenure50): 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, history4.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], cblue, labelDemand Points ) plt.scatter( facility_points[:,0], facility_points[:,1], cred, markers, labelFacility Candidates ) plt.scatter( facility_points[selected,0], facility_points[selected,1], cgreen, marker*, s200, labelSelected 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, alpha0.3 ) plt.legend() plt.title(Facility Location Allocation) plt.show() # 运行算法 ts TabuSearch(dist_matrix, demand, capacity, p3) best_sol, best_val, history ts.solve() print(fBest 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个需求点可以分区域求解先聚类将需求点分组再分别求解并行计算使用多进程同时评估多个邻域解记忆化缓存已评估的解避免重复计算from multiprocessing import Pool def parallel_evaluate(chroms): with Pool() as p: return p.map(evaluate, chroms)5.3 常见问题排查遇到过的一些典型问题及解决方案算法陷入局部最优增加禁忌长度引入重启机制每N代重置当前解结合模拟退火的接受准则运行时间过长优化距离矩阵计算使用KDTree实现快速评估方法增量式计算设置早期终止条件无法满足容量约束惩罚函数法在目标函数中加入惩罚项修复算子专门处理不可行解的修复过程# 惩罚函数示例 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