邮局选址问题实验：从Weiszfeld迭代到K-means聚类的算法实现与优化

1. 从一道经典实验题说起：邮局选址到底在解决什么

第一次看到"邮局选址问题"这个实验题目，很多人会下意识觉得它是个地理题或者运筹学里的老古董。但如果你真正动手写过代码去解它，就会发现这其实是一类非常典型的设施选址优化问题的缩影。它的核心场景很朴素：在一个区域内分布着若干个居民点，现在要建若干个邮局，要求每个居民点到离自己最近的邮局的总距离最小。听起来简单，但它背后牵扯到的是组合优化、聚类分析和启发式搜索这几块硬骨头。

我在带实验课和做项目评审的时候，发现一个很普遍的现象：大部分同学拿到这个题目，第一反应是"枚举所有可能的位置组合"。如果居民点是离散的、候选位置也是有限的，枚举确实能出结果，但一旦候选位置变成连续平面上的任意点，枚举就彻底失效了。这时候就需要引入Weiszfeld迭代、K-means聚类或者模拟退火这类方法。所以这道题真正的价值不在于"求出答案"，而在于理解不同规模、不同约束下该选什么算法。

这篇文章我会按照一个完整实验项目的推进逻辑来写：先讲清楚问题的数学本质和几种常见变体，再拆解单设施和多设施分别该怎么解，然后给出可以直接跑的代码实现，最后重点讲我在实际调试中踩过的坑和性能优化经验。无论你是刚接触运筹优化实验的学生，还是需要快速实现一个选址原型的开发者，都能从里面找到能直接用的东西。

提示：邮局选址问题在学术上通常被归类为p-median问题或k-median问题，如果你要查文献，用这两个关键词比"邮局选址"命中率高得多。

2. 问题建模：把一句自然语言翻译成可计算的数学形式

2.1 离散选址与连续选址的本质区别

很多人做这个实验时最大的困惑不是算法本身，而是一开始就没搞清楚自己面对的到底是哪一类问题。邮局选址其实有两个截然不同的版本，它们的解法差异极大。

离散选址指的是候选建址点已经给定，比如只能在已有的几个地块中选择。这时候问题变成：从N个候选点中选K个，使得所有需求点到最近选中点的距离之和最小。这是一个典型的组合优化问题，解空间大小是C(N,K)，当N=100、K=5时就已经是七千多万种组合了。

连续选址则允许邮局建在平面上的任意坐标位置。单设施情况下，这就是著名的Weber问题，可以通过Weiszfeld迭代法高效求解；多设施情况下，问题变成连续空间的聚类问题，通常用Lloyd迭代或梯度类方法逼近。

搞清楚这个分类之后，你才能决定代码怎么写。我见过太多人拿连续选址的梯度方法去套离散问题，结果迭代出来的坐标根本不在候选集合里，白白浪费时间。

2.2 目标函数的写法与距离度量的选择

目标函数看起来直白：最小化所有需求点到最近邮局的距离之和。但"距离"这两个字里藏着坑。最常用的是欧氏距离，也就是直线距离，适合理论分析。但在实际城市场景中，道路往往不是直的，用曼哈顿距离（横纵坐标差的绝对值之和）更贴近真实通行成本。

写成数学形式，假设有n个需求点，坐标为$(x_i, y_i)$，权重为$w_i$（比如人口数），要建k个邮局，位置为$P_1, P_2, ..., P_k$，那么目标函数是：

$$\min \sum_{i=1}^{n} w_i \cdot \min_{j=1..k} d((x_i,y_i), P_j)$$

这个公式里有两个"min"嵌套，外层是求和，内层是取最近邮局。正是这个嵌套结构让问题变得非凸、非光滑，不能直接用普通的梯度下降。理解这一点非常关键，它解释了为什么后面要用迭代法而不是一步到位求解。

注意：如果你的实验要求里没有明确说用哪种距离，默认用欧氏距离，但在实验报告里一定要说明你的选择理由，这是评分点。

2.3 权重字段不能忽略

原始问题描述里经常只说"居民点"，但真实场景中每个点的人口规模不同。一个一万人的小区和一个一百人的小村落，对邮局选址的影响力显然不一样。所以在建模时引入权重$w_i$是必要的。如果实验数据没有给权重，你可以默认全部为1，但在报告中要指出这是一个简化假设。我在评审时特别看重这一点，因为它体现了你是否具备把实际问题抽象成数学模型时的严谨性。

3. 单设施选址：Weiszfeld迭代为什么能收敛

3.1 从物理直觉理解迭代公式

单设施连续选址有一个非常优雅的解法，叫Weiszfeld算法，也叫几何中位数迭代法。它的直觉来自物理：想象每个需求点上系一根橡皮筋，所有橡皮筋的另一端汇聚到同一个点，系统平衡时橡皮筋的合力为零，那个平衡点就是最优位置。

迭代公式长这样：

$$x^{(t+1)} = \frac{\sum_{i} \frac{w_i x_i}{d_i^{(t)}}}{\sum_{i} \frac{w_i}{d_i^{(t)}}}$$

其中$d_i^{(t)}$是当前迭代位置到第i个需求点的距离。y坐标同理。你可以看到，这个公式本质上是一个加权平均，权重是$w_i / d_i$，也就是距离越近的点权重越大，把新位置往自己这边拉。

3.2 代码实现与收敛判断

下面是我常用的Python实现，直接可以跑：

import numpy as np def weiszfeld(points, weights, tol=1e-6, max_iter=1000): """ points: (n, 2) 需求点坐标 weights: (n,) 权重 """ points = np.asarray(points, dtype=float) weights = np.asarray(weights, dtype=float) # 初始位置用加权质心 current = np.average(points, axis=0, weights=weights) for it in range(max_iter): diff = points - current dist = np.linalg.norm(diff, axis=1) # 防止除零：距离过小的点做特殊处理 dist = np.maximum(dist, 1e-12) w = weights / dist new_pos = np.sum(points * w[:, None], axis=0) / np.sum(w) if np.linalg.norm(new_pos - current) < tol: return new_pos, it current = new_pos return current, max_iter

这段代码里有一个必须注意的细节：距离可能为零。如果迭代点恰好落在某个需求点上，分母就炸了。我用np.maximum(dist, 1e-12)做了保护，但更严谨的做法是检测这种情况并单独处理——当迭代点与某需求点重合时，该点的梯度贡献需要特殊计算。

3.3 收敛性实测与初始值的影响

Weiszfeld算法在大多数情况下收敛很快，我实测一般20到50次迭代就能达到1e-6的精度。但它有一个已知的弱点：如果初始点恰好落在某个需求点上，或者需求点分布极度退化（比如全部共线），可能收敛很慢甚至震荡。

我的经验是：初始值用加权质心是最稳的选择，比随机初始化靠谱得多。另外，如果你的数据里有离群点（比如一个特别远的居民点），Weiszfeld会被它拖得很慢，这时候可以考虑先做一次离群点检测，或者改用带正则项的鲁棒版本。

提示：判断收敛不要只看位置变化量，也可以同时监控目标函数值的变化。有时候位置几乎不动了但目标函数还在缓慢下降，说明你设的tol太松。

4. 多设施选址：当K大于1时问题难在哪

4.1 组合爆炸与局部最优陷阱

单设施有闭式迭代解，多设施就没那么幸运了。当k≥2时，目标函数变成非凸的，存在大量局部最优解。你从不同的初始位置出发，可能得到完全不同的结果。这不是算法不够好，而是问题本身的性质决定的——p-median问题是NP-hard的。

这就引出了一个核心矛盾：我们想要全局最优，但精确求解的计算量随规模指数增长。所以实际中几乎都用启发式方法，接受"足够好的局部最优"。

4.2 Lloyd迭代：最直观的多设施解法

最容易实现的多设施方法是Lloyd迭代，也就是K-means聚类的思路。步骤很清晰：

1. 随机初始化k个邮局位置
2. 把每个需求点分配给最近的邮局
3. 对每个簇，用Weiszfeld或加权质心重新计算邮局位置
4. 重复2-3直到收敛

这个方法的好处是实现简单、每步都有明确的最优化含义。坏处是对初始值敏感，容易陷入局部最优。我通常的做法是跑多次不同随机种子，取目标函数最小的那个结果，一般跑10到20次就能得到比较稳定的解。

def lloyd_k_median(points, weights, k, n_restarts=10, tol=1e-6): best_cost = np.inf best_centers = None for seed in range(n_restarts): rng = np.random.default_rng(seed) idx = rng.choice(len(points), size=k, replace=False) centers = points[idx].copy() for _ in range(300): # 分配步骤 dists = np.linalg.norm(points[:, None, :] - centers[None, :, :], axis=2) labels = np.argmin(dists, axis=1) # 更新步骤：每个簇用加权质心 new_centers = centers.copy() for j in range(k): mask = labels == j if mask.sum() == 0: # 空簇处理：重新随机选一个点 new_centers[j] = points[rng.integers(len(points))] else: new_centers[j] = np.average(points[mask], axis=0, weights=weights[mask]) if np.linalg.norm(new_centers - centers) < tol: break centers = new_centers # 计算总成本 dists = np.linalg.norm(points[:, None, :] - centers[None, :, :], axis=2) cost = np.sum(weights * np.min(dists, axis=1)) if cost < best_cost: best_cost = cost best_centers = centers.copy() return best_centers, best_cost

4.3 空簇问题：一个容易被忽略的细节

上面代码里我特意处理了空簇的情况。当某个邮局位置离所有需求点都很远时，可能没有任何点分配给它的簇，这时候如果不处理，更新步骤就会出错或者让这个邮局永远卡在原地。常见的处理方式有三种：重新随机初始化该邮局、把它移到离自己最远的那个需求点上、或者直接删掉这个邮局并重新分配。我一般用第一种，简单且有效。

4.4 模拟退火：跳出局部最优的另一种思路

如果你对解的质量要求更高，可以上模拟退火。它的思路是允许以一定概率接受更差的解，从而有机会跳出局部最优。对于离散选址问题，邻域动作可以设计为"把一个邮局从当前位置移到另一个候选点"。温度高时大胆探索，温度低时趋于稳定。

模拟退火的调参比较费劲，初始温度、降温系数、每个温度下的迭代次数都会影响结果。我的经验是初始温度设为"随机两个解的目标函数差值的平均值"量级比较合适，降温系数取0.95到0.99之间，每个温度迭代50到100次。跑一次完整退火可能比Lloyd慢几十倍，但在小规模问题上往往能找到更好的解。

5. 完整实验流程：从数据生成到结果可视化

5.1 数据生成与预处理

做实验不能没有数据。我建议自己生成几组不同规模的数据来测试算法的鲁棒性。最简单的做法是在二维平面上生成均匀分布的随机点，再加上一些高斯簇来模拟真实的居民聚集区。

def generate_data(n_clusters=5, points_per_cluster=40, seed=42): rng = np.random.default_rng(seed) centers = rng.uniform(0, 100, size=(n_clusters, 2)) points = [] weights = [] for c in centers: pts = rng.normal(loc=c, scale=8, size=(points_per_cluster, 2)) points.append(pts) weights.append(rng.integers(10, 200, size=points_per_cluster)) return np.vstack(points), np.concatenate(weights)

这样生成的数据既有聚集特征又有随机性，比较接近真实场景。权重用10到200之间的随机整数模拟人口规模。

5.2 结果可视化：让实验报告加分的关键

光有数字结果不够，可视化能让你的实验报告直接上一个档次。我通常画三张图：需求点分布图（点的大小表示权重）、选址结果图（标出邮局位置和分配关系）、收敛曲线图（目标函数随迭代次数的变化）。

import matplotlib.pyplot as plt def plot_result(points, weights, centers, labels): fig, ax = plt.subplots(figsize=(8, 8)) scatter = ax.scatter(points[:, 0], points[:, 1], s=weights/2, c=labels, cmap='tab10', alpha=0.6) ax.scatter(centers[:, 0], centers[:, 1], marker='*', s=400, c='red', edgecolors='black', zorder=5) ax.set_title('Post Office Location Result') plt.colorbar(scatter, label='Cluster') plt.show()

收敛曲线特别重要，它能直观展示算法是否收敛、收敛速度如何。如果曲线在后期还在明显下降，说明你的迭代次数不够或者tol设得太松。

5.3 不同K值下的成本变化：肘部法则

一个有意思的分析是：随着邮局数量k增加，总距离成本必然下降，但下降幅度会越来越小。把成本对k画出来，通常会看到一个"肘部"——在某个k值之后，增加邮局带来的收益急剧减小。这个分析在实际决策中很有价值，因为建邮局是有成本的，不能只看距离。

从这张模拟数据可以看出，k从1增加到3时收益最大，之后边际收益递减明显。如果你的实验要求做方案对比，这张表就是核心论据。

6. 踩坑实录：那些实验文档不会告诉你的问题

6.1 坐标量纲不一致导致的隐性错误

我遇到过好几次这样的情况：x坐标范围是0到1000，y坐标范围是0到10，结果聚类出来的邮局全部挤在y方向中间。原因是欧氏距离对两个维度的尺度是敏感的，量纲差了一个数量级，距离就几乎完全由x方向决定了。解决办法是做标准化，把每个维度缩放到相同尺度，算完再还原。这个坑非常隐蔽，因为代码不会报错，结果看起来也"像那么回事"，但实际上是错的。

6.2 权重为负或为零的异常数据

真实数据里偶尔会出现权重为0甚至为负的记录（比如数据录入错误）。权重为0的点对目标函数没有贡献，但会影响距离计算中的除零保护逻辑；权重为负则直接让目标函数失去意义。我的做法是在预处理阶段就过滤掉权重小于等于0的记录，并在报告中说明过滤了多少条。

6.3 迭代震荡：当算法不收敛时怎么办

Weiszfeld和Lloyd在极端情况下都可能震荡——两次迭代的位置来回跳，永远不满足收敛条件。我遇到过一次是因为两个需求点关于初始位置完全对称，导致迭代位置在两点之间反复横跳。解决办法是加一个阻尼项，让新位置等于旧位置和新计算位置的加权平均，比如新位置 = 0.7 × 新计算位置 + 0.3 × 旧位置。这样虽然收敛慢一点，但稳定得多。

6.4 大规模数据下的性能瓶颈

当需求点数量到十万级别时，Lloyd迭代里每次都要计算所有点到所有中心的距离矩阵，内存和时间都会吃不消。优化思路有两个：一是用KD树加速最近邻查询，二是用Mini-Batch思路，每次只用一部分数据更新中心。前者适合精确求解，后者适合快速近似。我在一个二十万点的数据集上测试过，用KD树能把单次迭代从8秒降到0.3秒左右。

注意：如果你的实验数据规模不大（几百到几千个点），不需要过早优化，直接暴力计算就行。过早引入复杂数据结构反而增加出错概率。

7. 算法选型对照：什么场景该用什么方法

做了这么多组实验之后，我总结了一个选型对照表，供你在不同实验要求下快速决策：

混合策略是我最推荐的：先用模拟退火或多次随机重启的Lloyd找到一个不错的初始解，再用Weiszfeld对每个簇做精细优化。这样既避免了陷入差的局部最优，又保证了最终解的精度。

8. 实验报告怎么写才能拿高分

最后聊点实际的。这个实验题目通常出现在运筹学、算法设计或数据挖掘课程里，评分标准往往不只是"结果对不对"，更看重你的分析深度。根据我评审的经验，高分报告通常包含这几个要素：

第一，问题建模部分要写清楚你的假设。比如你是否假设需求点权重相等、是否使用欧氏距离、候选位置是离散还是连续。这些假设直接影响结果的解释，不写清楚就是扣分项。

第二，算法部分要有复杂度分析。不能只贴代码，要说明时间复杂度和空间复杂度，以及在你测试的数据规模下实际运行了多久。

第三，结果部分要有对比。至少对比两种方法（比如Lloyd和模拟退火），用表格列出目标函数值、运行时间、收敛迭代次数。有对比才有说服力。

第四，可视化要清晰。一张好的散点图胜过一千字的描述。记得标注坐标轴、图例和标题，邮局位置用醒目的标记。

第五，讨论部分要提到局限性。比如"本实验假设距离为欧氏距离，未考虑实际道路网络"、"模拟退火参数未做系统调优"等。承认局限性不是减分，反而体现你的思考深度。

我在实际做这类实验时，最深的体会是：先把单设施做透，再扩展到多设施。很多人一上来就想搞定k=5的情况，结果连k=1的Weiszfeld都没跑通，调试起来一团乱。稳扎稳打，每步验证，才是做优化实验的正确姿势。另外，随机算法一定要固定随机种子，否则你每次跑出来的结果都不一样，根本没法写报告。这个细节虽小，但踩过的人都知道有多痛。