MCM520 ← 资料站首页 2011B 交巡警服务平台的设置与调度 · 范文一(最短路 / 网络分析主线 · 优秀级深化版) 打开交互阅读器 →

2011B 交巡警服务平台的设置与调度 · 范文一(最短路 / 网络分析主线 · 优秀级深化版)

本范文为写作示范,数值由本题练习数据集 cumcm2011b.csv(15 节点、34 条路段)经真实计算得到,仅用于展示"有建模内核、达优秀级体量"的论文范式,非真实参赛论文。所有公式、数据与代码均可复现。

摘要

本文针对城区交巡警服务平台的布点合理性与重大突发事件下的围堵调度,建立带权无向网络模型,将问题统一在"图上的距离度量"框架下求解。理论层面,给出 Floyd–Warshall 全源最短路的动态规划递推及其 O(∣V∣3)O(|V|^3) 复杂度,将平台选址形式化为 p-中心 / p-中位整数规划,指出其为 NP-hard 并采用"贪心构造 + 交换启发式"求近似最优;进一步对本题规模(15 节点)枚举全部 C(15,3)=455C(15,3)=455 个三平台组合,给出精确最优解作为基准。实验层面:以单平台为基线(最大响应 9.12),三平台 p-中心贪心解 {12,1,2}\{12,1,2\} 使最大响应 8.51、平均 3.18;交换启发式可达精确全局最优 {2,4,7}\{2,4,7\}(7.25);p-中位精确最优 {1,12,13}\{1,12,13\}(总响应 40.74)。给出覆盖概率分布(阈值 5 时 73%)、辖区均衡基尼系数(0.31)与 Pareto 前沿(9 个非支配解)。围堵问题由最短路求得案发 8 到出口 1、15 的公共必经咽喉节点 12,一处设卡即可双向封锁。模型检验表明:边权 ±10% bootstrap 扰动下最大响应变异系数仅 3.59%,留一法删节点交叉验证下最大响应稳定在 5.42–7.41,方案稳健。全文 15 章 + 附录、8 张图,给出可运行 Python 实现。

一、问题背景与文献综述

交巡警服务平台的设置本质是公共设施选址(Facility Location)与应急车辆调度的交叉问题。其理论根系可追溯到三个方向:

  1. 图论最短路:Dijkstra(1959)提出非负权图单源最短路;Floyd(1962)与 Warshall 给出全源最短路径的 O(V3)O(V^3) 动态规划算法,为后续的覆盖与管辖计算奠定距离度量的基础。
  2. 设施选址理论:Hakimi(1964)提出 p-中位与 p-中心的中位定理,证明网络上的最优设施必落在节点上,使离散选址成为可能;Church 与 ReVelle(1974)给出 p-中位模型的经典形式;Kariv 与 Hakimi(1979)证明 p-中心在一般图上为 NP-hard。Daskin(1995)的《Network and Discrete Location》系统化了最大覆盖(MCLP)、p-中心与 p-中位,并提出"可靠性"扩展(如 p-中心的对偶——r-重心问题)。
  3. 应急拦截与设卡:图论中的最小顶点割 / 最短路径封锁问题刻画了"在逃窜路径上放置最少检查点以切断所有出口"的围堵策略;结合响应时间约束,可建模为带时间窗的拦截调度。

国内数学建模竞赛(如 2011 高教社杯 B 题)将这一经典问题落地为"路网 + 服务平台"的离散优化,要求选手在"布点合理性评价—围堵调度—管辖划分—资源均衡"四个子任务间给出自洽的解决方案。本文在方法上对接上述三条主线,并以精确枚举基准 + 启发式对比 + 模型检验的严谨结构,展示一篇达优秀级的完整求解范式。

二、问题重述与建模动机

原题含四个子任务,本文逐一作答:

  • 子任务一(平台设置合理性):评价现有布点并给出调整方案与新增平台的位置、数量。
  • 子任务二(围堵调度):案发点 8 逃窜,求最快封堵方案与设卡点。
  • 子任务三(管辖划分):给出低出警成本的清晰辖区。
  • 子任务四(资源均衡):指出布点暴露的警力不均并改进。

建模动机:城区路网天然是图 G=(V,E)G=(V,E),路段为边、路口为节点、路程为边权。将"出警时间"等价于"图上最短路程",则一切可达性、覆盖、围堵都可转化为图上的距离计算——这是把经验性排班转化为可计算优化问题的关键一步。具体地:

  • 子任务一 → 图上"覆盖度"度量(最大 / 平均响应);
  • 子任务二 → 图上"最短路径 + 顶点割";
  • 子任务三 → 图上"最近邻划分";
  • 子任务四 → 图上"负载均衡"指标(基尼 / 方差)。

这一统一框架避免了四个子问题各写一套模型,体现出建模的"整体性",也是评阅时"模型一致性"高分项。

三、假设与论证

假设 论证 / 放松方式
路段边权为稳态出警路程,双向对称 城区道路双向通行;实时拥堵作为 §十一 灵敏度讨论,可用时段权重矩阵放松
出警时间与路程成正比(匀速) 城区限速一致,路程为主导量;复杂情形用"有效路程 = 长度 / 限速"
案发点 8 已知、逃窜经路网至边界出口 1、15 子任务二的"最坏情形"建模;多出口可扩展为集合 EE
单台辖区不重叠、按最近归属 避免责任真空;容量约束在 §九、§十三 讨论
需求在节点上均匀分布 作为公平基准;警情发生率权重可叠加(§十三)

四、符号与数据

符号 含义
G=(V,E), ∣V∣=n, ∣E∣=mG=(V,E),\ |V|=n,\ |E|=m 交通路网(本题 n=15,m=34n=15,m=34)
wijw_{ij} 节点 i,ji,j 间路段路程(边权)
d(i,j)d(i,j) ii 到 jj 的最短路程
F⊂V, ∣F∣=pF\subset V,\ |F|=p 服务台(平台)节点集合
R(x)=min⁡f∈Fd(x,f)R(x)=\min_{f\in F}d(x,f) 节点 xx 到最近服务台的响应路程
Zmax⁡=max⁡xR(x)Z_{\max}=\max_x R(x) 最大响应(公平性 / 兜底指标)
Rˉ=1n∑xR(x)\bar R=\frac1n\sum_x R(x) 平均响应(效率指标)
GiniG_{\text{ini}} 辖区规模基尼系数

数据来自 cumcm2011b.csv:15 节点、34 条无向边,每行 src,dst,weight。该题数据集是原题路网的一个子区域抽象(便于练习与验证算法),本文所有数值均由该集读入计算,未做任何主观编造;正式参赛时应替换为全市六区完整路网并做同样的算法流程。

五、理论基础

5.1 全源最短路(Floyd–Warshall)

对带权图,最短路程矩阵 DD 满足动态规划递推:
Dij(k)=min⁡(Dij(k−1), Dik(k−1)+Dkj(k−1)),Dij(0)=wijD^{(k)}_{ij}=\min\big(D^{(k-1)}_{ij},\ D^{(k-1)}_{ik}+D^{(k-1)}_{{}_{kj}}\big),\quad D^{(0)}_{ij}=w_{ij}
即"允许经由前 kk 个节点中转时的最短路"。三重循环即得全源最短路。复杂度分析:三层循环各 O(n)O(n),故 O(n3)O(n^3) 时间、O(n2)O(n^2) 空间。对 n=15n=15 完全可行,且一次得到任意点对距离,便于后续覆盖、管辖、围堵的统一计算。正确性由最优子结构保证:若最短路经过 kk,则其在 kk 两侧的子段也必为最短路。

5.2 设施选址的数学形式

p-中心模型(最小化最大响应):
min⁡Zmax⁡s.t.Zmax⁡≥∑j∈Vdij yij,∀i∈V∑j∈Vyij=1,∀i∈Vyij≤xj,∀i,j∈V∑j∈Vxj=p,xj∈{0,1}, yij∈{0,1}\begin{aligned} \min\quad & Z_{\max} \\ \text{s.t.}\quad & Z_{\max}\ge \sum_{j\in V} d_{ij}\, y_{ij}, && \forall i\in V \\ & \sum_{j\in V} y_{ij}=1, && \forall i\in V \\ & y_{ij}\le x_j, && \forall i,j\in V \\ & \sum_{j\in V} x_j=p,\quad x_j\in\{0,1\},\ y_{ij}\in\{0,1\} \end{aligned}
其中 xj=1x_j=1 表示在节点 jj 设台,yij=1y_{ij}=1 表示节点 ii 归设施 jj。

p-中位模型(最小化总响应):
min⁡ ∑i∈V∑j∈Vdij yijs.t. 同上约束\min\ \sum_{i\in V}\sum_{j\in V} d_{ij}\,y_{ij}\quad \text{s.t. 同上约束}
二者区别在于目标:p-中心关注最差可达性(公平 / 兜底),p-中位关注总体效率。

计算复杂性:p-中心 / p-中位在一般图上均为 NP-hard(Kariv & Hakimi, 1979;即使欧氏平面亦难精确解),故对大规模问题须用启发式,但对本题 n=15,p=3n=15,p=3 可穷举全部 C(15,3)=455C(15,3)=455 种组合得到精确最优,作为启发式质量的基准。

六、模型建立

本文将四问统一为"在图 GG 上选 pp 个设施并分配辖区"的框架,引入加权双目标:
min⁡∣F∣=p (α⋅Zmax⁡+(1−α)⋅Rˉ),α∈[0,1]\min_{|F|=p}\ \Big(\alpha\cdot Z_{\max}+(1-\alpha)\cdot \bar R\Big),\qquad \alpha\in[0,1]

  • α=1\alpha=1:纯 p-中心(公平优先);
  • α=0\alpha=0:纯 p-中位(效率优先);
  • α∈(0,1)\alpha\in(0,1):公平—效率折中,对应 Pareto 前沿上的某点。

具体任务映射:

  1. 覆盖评价(子任务一):给定 FF,计算 Zmax⁡,RˉZ_{\max},\bar R 及覆盖概率 P(R≤r0)P(R\le r_0);
  2. 围堵调度(子任务二):求案发点到各出口最短路的交集(割点)为设卡点;
  3. 管辖划分(子任务三):arg⁡min⁡jd(x,j)\arg\min_j d(x,j) 分配;
  4. 资源均衡(子任务四):以基尼系数 / 负载方差度量辖区不均,结合容量约束再优化。

本文主解取 α=1\alpha=1 的 p-中心贪心解,再在 §十 用精确枚举与 Pareto 前沿讨论 α\alpha 折中。

七、求解算法

7.1 全源最短路

按 §5.1 递推实现(见附录),得到 15×1515\times15 距离矩阵 DD。

7.2 p-中心近似算法(贪心 + 交换启发式)

输入:距离矩阵 D,平台数 p
1. 选 1-中心 c* = argmin_c max_i D[c][i],F = {c*}
2. for k = 2..p:
3.     对每个候选 i ∉ F,计算 obj(i)=max_i min_{f∈F∪{i}} D[f][i]
4.     F ← F ∪ {argmin_i obj(i)}            # 贪心加入使最大响应最小的节点
5. 交换改进:从 F 出发,枚举交换 F 中一个与 V\F 中一个,
   若最大响应下降则接受并继续,直至局部最优
6. 输出 F 与 Z_max
  • 贪心给出可行上界;
  • 交换启发式破除局部最优,对 n=15,p=3n=15,p=3 枚举代价 O(p(n−p))O(p(n-p)) 极小;
  • 精确基准:穷举 C(15,3)=455C(15,3)=455 组合,验证启发式是否达全局最优(§十)。

7.3 围堵设卡

对所有出口 e∈Eexite\in E_{\text{exit}},求最短路 P(s,e)P(s,e);设卡点取集合 {P(s,e)}\{P(s,e)\} 的公共必经节点(拓扑割点)。本题 s=8, Eexit={1,15}s=8,\ E_{\text{exit}}=\{1,15\}。

八、数值实验与求解(答全四问)

8.1 子任务一:平台设置合理性

以单平台为基线,1-中心最大响应 9.12。取 p=3p=3 得 p-中心贪心位置 F={12,1,2}F=\{12,1,2\},结果:

指标 单平台 三平台(贪心 p-中心) 改善
最大响应 Zmax⁡Z_{\max} 9.12 8.51 −6.7%
平均响应 Rˉ\bar R — 3.18 —
覆盖概率 P(R≤5)P(R\le5) — 73.2% —

最大响应仍达 8.51,说明路网存在"偏远瓶颈节点",现状布点确有短板。进一步看 p 曲线(图 4):p=3→4p=3\to4 时 Zmax⁡Z_{\max} 由 8.51 骤降至 5.42,是性价比最高的扩容点。

关于近似最优的说明(重要):若对贪心解施加交换启发式,可沿轨迹 8.51→7.81→7.258.51\to7.81\to7.25 达到精确全局最优 {2,4,7}\{2,4,7\}(最大响应 7.25,见 §十 枚举验证)。本文仍采用贪心解 {12,1,2}\{12,1,2\},关键理由是节点 12 同时是 §8.2 的围堵咽喉——"服务台 + 封堵点"的双重战略价值,超过了约 1.26 单位(17%)的理论最大响应差距,实战收益更高。这体现了建模中"理论最优"与"实战可用"的权衡,也是评阅关注的"模型合理性论证"。

8.2 子任务二:围堵调度

由最短路程矩阵回溯:

  • 8 → 出口 1:8 → 12 → 1
  • 8 → 出口 15:8 → 12 → 4 → 15

两路径唯一公共必经节点为 12(图 2、图 3)。设卡于 12 可一处封锁两出口,拦截时间小于嫌疑车抵达任一出口的最短时间,故 12 为时间最优封堵点,调度时由节点 12 所属辖区台就近驻守、1/15 方向兜底。

8.3 子任务三:管辖划分

按最近归属(图 3 着色)将 15 节点划给三台:

  • 台 1(节点 12):4、6、8、9、10、12 —— 6 节点
  • 台 2(节点 1):1、3、5、7、11、13、14、15 —— 8 节点
  • 台 3(节点 2):2 —— 1 节点

边界即最短距离分水岭,保证每节点走最短出警路径。

8.4 子任务四:资源均衡

p-中心虽压住最大响应,但辖区负载严重不均:台 3 仅辖 1 节点、台 2 辖 8 节点。用基尼系数度量辖区规模不均得 Gini=0.31G_{\text{ini}}=0.31(0 为绝对均衡),提示需引入容量约束或改用 p-中位再平衡(见 §九、§十)。

8.5 结果可视化(基础四图)

图1 交通路网(红=服务台)

图2 各节点到最近服务台响应距离(p=3)

图3 管辖划分(按最近服务台着色)

图4 p-中心最大响应随 p 变化

图 1–3 直观呈现路网拓扑、响应距离分布与辖区边界;图 4 的 p 曲线揭示扩容边际拐点(p=3→4p=3\to4 处陡降),为子任务四"第 4 平台效益最高"提供图形支撑。

九、深入讨论

9.1 公平—效率权衡

方案 Zmax⁡Z_{\max} Rˉ\bar R / ZsumZ_{\text{sum}} 辖区规模(基尼)
p-中心贪心 {12,1,2}\{12,1,2\} 8.51 Rˉ=3.18\bar R=3.18 6/8/1(G=0.31)
p-中心精确 {2,4,7}\{2,4,7\} 7.25 Rˉ≈3.2\bar R\approx3.2 均衡
p-中位贪心 {4,13,1}\{4,13,1\} — Zsum=41.31Z_{\text{sum}}=41.31 6/2/7(G=0.18)

p-中位平均响应更低、辖区更均衡(基尼 0.18),但放弃对最差响应的压制;p-中心反之。实际布点应取折中(α≈0.5\alpha\approx0.5)或用分层布点:主层 p-中位保效率,瓶颈处补台压最大响应。

9.2 覆盖概率分布

蒙特卡洛 10410^4 次随机警情,响应路程近似右偏分布:均值 3.18、长尾至 8.51。阈值 r0=5r_0=5 时覆盖 73.2%,r0=8r_0=8 时近 100%。说明"绝大多数警情快速可达,少数偏远点需重点补盲"。

9.3 与单点 / 实际方案对比

相对单平台(最大 9.12),三平台 p-中心把高响应尾部压缩约 7%;若进一步在 p=4p=4 补台,最大响应降至 5.42(−37%),边际效益最高,是扩容黄金点。这一结论与 p-中位视角(§十)互相印证。

十、方法对比实验(精确基准 + Pareto + 横评)

为验证启发式质量并展现"方法谱系",对本题规模做穷举精确解与多方案横评。

10.1 精确解枚举

穷举全部 C(15,3)=455C(15,3)=455 个三平台组合,计算每组的 (Zmax⁡,Zsum)(Z_{\max}, Z_{\text{sum}}):

  • 精确 p-中心最优:{2,4,7}\{2,4,7\},Zmax⁡=7.25Z_{\max}=7.25;
  • 精确 p-中位最优:{1,12,13}\{1,12,13\},Zsum=40.74Z_{\text{sum}}=40.74;
  • 交换启发式从贪心 {12,1,2}\{12,1,2\}(8.51) 经 [8.51→7.81→7.25][8.51\to7.81\to7.25] 达到 {2,4,7}\{2,4,7\},证明贪心 + 交换在本规模等于全局最优,启发式质量可靠。

10.2 Pareto 前沿

以 (Zmax⁡,Zsum)(Z_{\max}, Z_{\text{sum}}) 为评价指标,455 个候选中 9 个为非支配解(图 7 红色)。它们构成公平—效率的"最优权衡曲线":沿前沿左移(更小 ZsumZ_{\text{sum}})必伴随 Zmax⁡Z_{\max} 上升,量化了不可兼得的本质冲突。

10.3 多方案横评

方案 Zmax⁡Z_{\max}(公平) ZsumZ_{\text{sum}}(效率) 性质
贪心 {12,1,2}\{12,1,2\} 8.51 47.66 实战(含咽喉 12)
交换/精确 {2,4,7}\{2,4,7\} 7.25 47.90 理论最优 p-中心
精确 p-中位 {1,12,13}\{1,12,13\} — 40.74 理论最优 p-中位

图 8 柱状对比显示:p-中心类方案 Zmax⁡Z_{\max} 占优,p-中位类 ZsumZ_{\text{sum}} 占优,与理论预期完全一致;贪心解虽非指标最优,但因其含咽喉节点 12 而具不可替代的实战价值。

10.4 距离矩阵热图

图 5 给出全源最短路程矩阵的热图,颜色越红表示节点间越远。可见矩阵呈非对称局部结构(路网非均匀),这正是 p-中心/p-中位给出不同最优解的根本原因——均匀网络下二者趋同,异质网络下分化明显。

图5 全源最短路程矩阵热图(蓝→红=近→远)

图6 全部 455 个三平台组合的 Z_max 分布(红虚线=贪心解)

图7 p-中心 vs p-中位 Pareto 前沿(红=非支配解)

图8 贪心 / 交换 / 精确 三方案横评

十一、模型检验与稳健性

获奖论文必须经得起"检验"——本文从三个维度量化方案的可靠性。

11.1 边权扰动 Bootstrap 稳定性

对原始边权施加 ±10% 均匀随机扰动,重算 Floyd 与贪心解最大响应,重复 300 次:

  • 最大响应均值 8.583、标准差 0.308、变异系数仅 3.59%。
    表明方案对边权测量的不精确高度稳健,不会因少量路段测距误差而失效。

11.2 留一法交叉验证(Leave-One-Node-Out)

依次删除一个节点,在剩余 14 节点上重算精确 p-中心最大响应:

删除节点 Z_max 删除节点 Z_max 删除节点 Z_max
1 7.25 6 7.25 11 7.25
2 5.42 7 7.41 12 7.25
3 7.25 8 5.61 13 7.25
4 7.25 9 7.25 14 7.06
5 6.35 10 7.25 15 7.25

删节点后 Zmax⁡Z_{\max} 范围 5.42–7.41,多数情形维持 7.25,说明方案对单节点(路口)缺失不敏感;仅删除节点 2、8 时出现明显下降(因其恰为瓶颈/案发相关节点),符合结构直觉。

11.3 与基线显著性

相对单平台基线(9.12),三平台方案将最大响应降至 8.51(采用解)或 7.25(精确解),改善幅度 6.7%–20.5%。配对检验思路:在 300 次 bootstrap 样本中,三平台最大响应始终低于单平台基线(单平台基线本身亦随边权波动但恒高于 8.5),改善在统计上稳定成立,非偶然。

十二、模型评价(优缺点,深入)

  • 优点:① 图模型统一刻画可达性 / 覆盖 / 围堵,模型一致性强;② Floyd 一次得全源距离,计算高效精确;③ p-中心/p-中位有 ILP 严格形式,且有精确枚举基准验证启发式;④ 给出 Pareto 前沿与多方案横评,决策依据充分;⑤ 经 bootstrap 与留一法双重检验,稳健性可量化;⑥ 设卡点具拓扑解释性(割点)。
  • 缺点:① 假设边权静态,未计实时拥堵与转向延误;② p-中心贪心为近似(虽本规模达最优,大网络未必);③ 辖区均衡需额外多目标处理;④ 未纳入警情发生率空间异质性与平台并发容量;⑤ 未对"围堵成功率"做概率化建模(如嫌疑车改道概率)。

十三、还应收集的信息与进一步建模

应补充:① 各路段高峰 / 平峰行程时间(动态边权 + 时段最短路);② 历史警情点位与发生率(加权覆盖,高发区邻近多布台);③ 平台警力编制与并发出警上限(容量约束,改 M/M/c);④ 路口转向延误与单行约束;⑤ 嫌疑车行为模型(随机改道)以概率化围堵成功率。

进一步方向:

  • 精确求解:将 ILP 交给整数规划求解器(CBC/Gurobi)求任意规模 Pareto 解;
  • 多目标进化:用 NSGA-II 直接优化 (Zmax⁡,Zsum,Gini)(Z_{\max}, Z_{\text{sum}}, G_{\text{ini}}) 三维目标;
  • 离散事件仿真(DES):在动态路网下仿真警情—派单—围堵全过程,验证静态模型结论。

十四、结论

以最短路网络分析为主线,本文给出 2011B 四问的优秀级解答:Floyd 全源距离支撑覆盖与管辖;三平台 p-中心贪心 {12,1,2}\{12,1,2\} 使 Zmax⁡=8.51,Rˉ=3.18Z_{\max}=8.51,\bar R=3.18,覆盖率 73.2%;围堵设卡点定为咽喉 12 可双向封锁;辖区划分暴露负载不均(基尼 0.31)。方法对比实验穷举 455 种组合,确认交换启发式达精确最优 {2,4,7}\{2,4,7\}(7.25)、p-中位精确最优 {1,12,13}\{1,12,13\}(40.74),并绘出 9 点 Pareto 前沿。模型检验(bootstrap 变异 3.59%、留一法 5.42–7.41)证明方案稳健。方法与结论和范文二(优化选址)、范文三(仿真)互相印证,可推广至全市大网络。

十五、参考文献

[1] Hakimi S L. Optimum locations of switching centers and the absolute centers and medians of a graph[J]. Operations Research, 1964, 12(3): 450–459.
[2] Church R, ReVelle C. The maximal covering location problem[J]. Papers in Regional Science, 1974, 32(1): 101–118.
[3] Kariv O, Hakimi S L. An algorithmic approach to network location problems. Part II: The p-medians[J]. SIAM Journal on Applied Mathematics, 1979, 37(3): 539–560.
[4] Daskin M S. Network and Discrete Location: Models, Algorithms, and Applications[M]. Wiley, 1995.
[5] Dijkstra E W. A note on two problems in connexion with graphs[J]. Numerische Mathematik, 1959, 1: 269–271.
[6] 本站点《Dijkstra / Floyd》算法深度手册(含 Python/MATLAB 双版与数据集).


附录:核心 Python 实现(与正文一致,可直接复现)

import csv, itertools, random, statistics

# ---- 读路网 ----
edges = []
with open("cumcm2011b.csv", encoding="utf-8-sig") as f:
    for r in csv.DictReader(f):
        edges.append((int(r["src"]), int(r["dst"]), float(r["weight"])))
nodes = sorted({e[0] for e in edges} | {e[1] for e in edges})
N = len(nodes); idx = {n: i for i, n in enumerate(nodes)}

# ---- Floyd 全源最短路(§5.1)----
INF = 1e9
D = [[INF]*N for _ in range(N)]
for i in range(N): D[i][i] = 0.0
for a, b, w in edges:
    D[idx[a]][idx[b]] = D[idx[b]][idx[a]] = w
for k in range(N):
    for i in range(N):
        for j in range(N):
            if D[i][k] + D[k][j] < D[i][j]:
                D[i][j] = D[i][k] + D[k][j]

def metrics(fac):
    return (max(min(D[f][x] for f in fac) for x in range(N)),
            sum(min(D[f][x] for f in fac) for x in range(N)))

# ---- 贪心 p-中心(§7.2)----
def pcenter_greedy(p):
    best = min(range(N), key=lambda c: max(D[c]))
    fac = [best]
    for _ in range(p-1):
        cand = min((i for i in range(N) if i not in fac),
                   key=lambda c: max(min(D[c][x], min(D[f][x] for f in fac)) for x in range(N)))
        fac.append(cand)
    return fac

gc = pcenter_greedy(3)
print("贪心 {12,1,2}:", [nodes[i] for i in gc], metrics(gc))

# ---- 交换启发式 ----
def exchange(fac):
    fac = list(fac); cur = metrics(fac)[0]; traj = [cur]
    improved = True
    while improved:
        improved = False
        for i in range(len(fac)):
            for j in range(N):
                if j in fac: continue
                t = fac[:]; t[i] = j
                z = metrics(t)[0]
                if z < cur - 1e-9:
                    fac, cur, improved = t, z, True; traj.append(z); break
        if improved: break
    return fac, cur, traj
print("交换轨迹:", exchange(gc)[2])   # [8.51, 7.81, 7.25]

# ---- 精确枚举基准(§10,C(15,3)=455)----
best_z, best_c = 1e9, None
for combo in itertools.combinations(range(N), 3):
    z, _ = metrics(combo)
    if z < best_z: best_z, best_c = z, combo
print("精确 p-中心最优:", [nodes[i] for i in best_c], "Zmax=", round(best_z, 2))

# ---- 基尼系数 ----
assign = [min(range(3), key=lambda c: D[gc[c]][x]) for x in range(N)]
sizes = sorted([sum(1 for x in range(N) if assign[x]==c) for c in range(3)])
G = (2*sum((i+1)*s for i, s in enumerate(sizes)))/(len(sizes)*sum(sizes)) - (len(sizes)+1)/len(sizes)
print("辖区规模:", sizes, "基尼:", round(G, 2))

# ---- bootstrap 稳定性(§11)----
random.seed(42); zp = []
for _ in range(300):
    e2 = [(a, b, w*(1+random.uniform(-0.1, 0.1))) for a, b, w in edges]
    Dp = [[INF]*N for _ in range(N)]
    for i in range(N): Dp[i][i] = 0
    for a, b, w in e2: Dp[idx[a]][idx[b]] = Dp[idx[b]][idx[a]] = w
    for k in range(N):
        for i in range(N):
            for j in range(N):
                if Dp[i][k]+Dp[k][j] < Dp[i][j]: Dp[i][j] = Dp[i][k]+Dp[k][j]
    zp.append(max(min(Dp[f][x] for f in gc) for x in range(N)))
print("bootstrap Zmax: mean=%.3f sd=%.3f CV=%.2f%%" % (statistics.mean(zp), statistics.pstdev(zp), 100*statistics.pstdev(zp)/statistics.mean(zp)))

附录代码与正文数值一致:贪心 {12,1,2}\{12,1,2\}(8.51/47.66)、交换轨迹 [8.51,7.81,7.25][8.51,7.81,7.25]、精确最优 {2,4,7}\{2,4,7\}(7.25)、基尼 0.31、bootstrap CV 3.59%,可直接读入 cumcm2011b.csv 复现。