简介这份资源是一篇完整的数学建模竞赛论文聚焦公园内道路设计的优化问题面向备战数学建模竞赛的高校学生与指导教师也适合作为图论算法落地实际场景的自学范例。论文以矩形公园的8个入口为对象要求在任意两入口最短路径不超过直线距离1.4倍的前提下最小化新建道路总长并按给定交叉点、自由布点、含矩形湖约束三种情形层层递进。包内共1个doc文件约495KB正文结构完整包含问题重述、模型假设、符号说明、数据分析与模型建立求解等章节并附Prim最小生成树、Floyd最短路、Fermat点与斯坦纳树的求解过程及结果数据。目前已有330人学习读者可借此掌握从抽象建模、算法组合选型到局部优化与结果验证的完整思路也可对照复现代码与数值结论自查适合需要拆解赛题套路、积累可复用模型框架的读者参考。1. 公园里修路这件事卡住人的从来不是最短路算法一个 200 米 × 100 米的矩形公园8 个入口全部钉在四条边上内部只允许在指定交叉点问题一给了 A、B、C、D 四个或自行选定的位置上修路要求任意两个入口之间的实际走行距离不超过两点直线距离的 1.4 倍同时内部新修道路的总长度最小。多数人第一反应是套 Dijkstra 或者 Floyd把最短路跑出来就交差但真正卡人的从来不是最短路算法本身最小生成树优化的是边权和压根不管直径而 1.4 倍恰恰是一个直径约束公园四周已经建好的路不计入成本等于系统白送了一张零权子图用不用、怎么接直接决定最终数字第三问再塞进一个 25 米见方的矩形人工湖当禁区路不能穿湖只能贴着走。这份数学建模论文把三问串成一条完整链路Prim 先建候选路网Floyd 只做合规校验超标的边用局部换边或者 Fermat 点重排遇到障碍再做可见性判断并用斯坦纳树收口最终落到 394.5 米、358.4 米、359.6 米三个结果。下文按“怎么把公园抽象成图、校验怎么跑、点怎么插、坑在哪”的顺序拆开讲。2. 把公园压成加权图12 个节点、两张距离表和一条 1.4 倍红线动手写代码之前先把物理场景翻译成图论语言。这一步做不干净后面所有算法的输出都是垃圾。公园道路设计问题里能作为图节点的只有两类东西8 个入口以及允许修路时打结的交叉点。矩形四条边上的其它点一律不能成为节点——题目写得很死新修道路与四周的连接只能与 8 个入口相通不能连到边上的其它位置。这条约束在很多人的模型里是缺失的结果程序跑出来一条从 (200, 30) 这种位置拐进公园的路看上去总长更短实际上题意都不满足。2.1 节点编号与坐标表问题一把四个交叉点直接给定了编号接在 8 个入口后面凑成 12 个点。这个编号顺序后面会反复用到尤其是距离矩阵的行列下标最好一次定死别中途改。编号名称x 坐标y 坐标位置说明1P1200下边2P2500下边3P31600下边4P420050右边5P5120100上边6P635100上边7P710100上边8P8025左边9A5075内部交叉点10B4040内部交叉点11C12040内部交叉点12D11570内部交叉点问题二把 9 到 12 号节点全部放开交叉点位置变成决策变量问题三在此基础上再加一个矩形禁区。所以整篇文章的骨架其实只有一套变的只是“哪些点固定、哪些点自由、哪些区域不能碰”。2.2 欧氏直连矩阵12 个点两两之间的直线距离第一张表是 12 个点两两之间的欧氏距离。它描述的是“如果这两个点之间直接拉一条路需要多长”。注意这里的“直连”是建模意义上的允许直连不是物理上一定可修——比如入口 1 和入口 8 之间直连要穿过公园内部那是要花钱的而沿边界走是免费的这两件事必须分开算。% 12 个节点的坐标行序即节点编号 1..12 XY [ 20 0; % P1 50 0; % P2 160 0; % P3 200 50; % P4 120 100; % P5 35 100; % P6 10 100; % P7 0 25; % P8 50 75; % A 40 40; % B 120 40; % C 115 70]; % D n size(XY, 1); D zeros(n); for i 1:n for j 1:n D(i,j) hypot(XY(i,1)-XY(j,1), XY(i,2)-XY(j,2)); end endhypot比sqrt(dx^2dy^2)在极端量级下更稳虽然这道题里坐标都是几十上百差异体现不出来但养成习惯没坏处。得到的D是一个对称矩阵对角线上是 0。这个矩阵有两个用途一是作为 Prim 的输入权重二是作为 1.4 倍判定的分母——判据里的“两点直线距离”指的就是D(i,j)不是沿边距离这一点特别容易搞混。2.3 沿边免费距离把 8 个入口投影到周长参数上第二张表是 8 个入口之间“只走公园四周已建道路”的最短距离。这张表决定了哪些入口对根本不需要在公园内部修路直接从候选集里排除能大幅缩小后面的搜索空间。计算方法不复杂但如果手工一个个数格子十个里有八个会出错。稳妥的做法是把矩形周长拉成一条数轴选定一个起点比如 P1沿顺时针方向累计弧长每个入口得到一个参数 t沿边最短距离就是min(|ti-tj|, L-|ti-tj|)其中 L 是周长 600 米。入口顺时针累计弧长 t米P120P250P3160P4250P5380P6465P7490P8575% 8 个入口按顺时针排列的周长参数 t [20 50 160 250 380 465 490 575]; L 600; % 2*(200100) E zeros(8); for i 1:8 for j 1:8 d abs(t(i) - t(j)); E(i,j) min(d, L - d); % 两个方向取短 end end这段代码只有四行但它把“沿边走多远”这件事彻底形式化了。只要t的顺序和方向固定结果就是唯一的不会出现两个人算出两个答案的情况。2.4 1.4 倍筛选哪 10 对入口必须走内部道路有了E和D筛选就是一次逐对比较如果E(i,j) 1.4 * D(i,j)说明走边界已经满足要求这一对入口之间不需要在公园内修任何东西否则这一对必须通过内部道路连通并且内部路径长度还得再受一次 1.4 倍约束。把 28 对入口全部过一遍落在“边界不够用”这一侧的有 10 对入口对直线距离1.4 倍上限沿边距离判定1→5141.4198.0240不可走边界1→6101.1141.5155不可走边界1→832.044.845不可走边界2→5122.1170.9270不可走边界2→6101.1141.5185不可走边界2→7107.7150.8160不可走边界3→464.089.690不可走边界3→5107.7150.8220不可走边界3→6160.1224.1295不可走边界3→7180.3252.4270不可走边界1→8 和 3→4 这两对是典型的“擦线”案例沿边 45 米对上限 44.8 米沿边 90 米对上限 89.6 米差距都不到半米。如果坐标或者弧长参数有半点马虎这两对就会被误判成合规最后结果差出一整条边。这也是后文避坑章节要重点说的一条。提示判据里的比较方向是“内部路径 ≤ 1.4 × 直线距离”不是“≤ 1.4 × 沿边距离”。见过不少建模稿把分母写成沿边距离那样约束会松得离谱1→5 直接变成合规问题一的最优解凭空少掉一大截。3. Prim 加 Floyd 的双层校验最小生成树为什么天生过不了直径约束把 12 个点、所有允许的直连边凑成一个完全图之后第一反应通常是求最小生成树因为“总长度最小”和“边权和最小”听起来是一回事。但这两件事只在没有直径约束的时候等价一旦加上 1.4 倍最小生成树就只是候选方案不是可行方案。这一章的核心就是把这个区别用代码跑出来。3.1 在完全图上跑 Prim拿到的是候选边集Prim 的逻辑很朴素从一个点出发每次挑一条连接“已选集合”和“未选集合”的最短边扩进去直到所有点都被覆盖。对 12 个点来说输出是 11 条边构成的树。function [T, W, total] prim_mst(D) % D 为对称距离矩阵T 为 (n-1)x2 的树边表total 为树的总权重 n size(D, 1); INF 1e9; visited false(1, n); visited(1) true; T zeros(n-1, 2); total 0; for k 1:n-1 best INF; bi 0; bj 0; for i 1:n if ~visited(i), continue; end for j 1:n if visited(j), continue; end if D(i,j) best best D(i,j); bi i; bj j; end end end if bi 0 error(图不连通检查距离矩阵中的 INF 项); end visited(bj) true; T(k,:) [bi bj]; total total best; end W total; endbi 0那个分支不是摆设。问题三加入湖之后如果把穿湖的边权重置成INF图有可能在某些参数下变得不连通这时候程序必须明确报错而不是悄悄返回一棵缺边的树让你在后面 Floyd 阶段看到一堆莫名其妙的 1e9。Prim 输出的树边表T才是真正要修的路。剩下的边在这个方案里不修但它们不是“不存在”而是“不划算”。这个区分在后面做局部优化时很关键换边本质上就是把某条不修的边加进来、把某条已修的边踢出去。3.2 Floyd 在树上只是复述不是优化很多人对 Floyd 有一个误解以为把它套在最小生成树上能“再优化一次”。不会的。树里任意两点之间的路径唯一Floyd 跑完得到的距离矩阵就是这条唯一路径的长度它没有任何别的选择可以比较。所以在这套流程里Floyd 的角色是校验器不是优化器。function [dist, nxt] floyd_all(W) % W 为邻接权重矩阵不存在的边填 infdist 为全源最短距离nxt 用于回溯路径 n size(W, 1); dist W; nxt repmat(1:n, n, 1); for i 1:n nxt(i,i) i; end for k 1:n for i 1:n for j 1:n if dist(i,k) dist(k,j) dist(i,j) - 1e-12 dist(i,j) dist(i,k) dist(k,j); nxt(i,j) nxt(i,k); end end end end end function p rebuild_path(nxt, i, j) % 根据 nxt 表回溯 i 到 j 的完整节点序列 p i; while p(end) ~ j p(end1) nxt(p(end), j); end end那个1e-12的容差是必须的。浮点加法不满足严格结合律dist(i,k)dist(k,j)和dist(i,j)在数学上相等的时候机器上可能差一个最低位不加容差会导致某些本不该更新的路径被更新nxt表被写乱回溯出来的路径出现环。这种 bug 极难查因为距离数值看起来完全正常。把树边填进邻接矩阵的时候记得对称填充树是无向的W inf(n); for k 1:size(T,1) i T(k,1); j T(k,2); W(i,j) D(i,j); W(j,i) D(i,j); % 无向图必须对称 end [dist, nxt] floyd_all(W);只填一半是另一个高频错误。填完之后矩阵不对称Floyd 跑出来的dist(i,j)和dist(j,i)不相等后面按入口对逐个校验时会得到自相矛盾的结果。3.3 逐对校验问题一里有两对超标校验对象只有那 10 对“边界不够用”的入口其余 18 对走边界天然合规不用管。跑完 Floyd 之后逐对取值和 1.4 倍上限比入口对树上路径路径长度1.4 倍上限是否合规1→51-2-10-9-12-5203.2198.0否1→61-2-10-9-6136.8141.6是1→81-832.044.8是2→52-10-9-12-5173.2170.9否2→62-10-9-6106.8141.6是2→72-10-9-6-7131.8150.8是3→43-464.089.6是3→53-11-12-5140.0150.8是3→63-11-12-9-6203.9222.1是3→73-11-12-9-6-7228.9252.3是只有 1→5 和 2→5 超标而且超得不多203.2 对 198.0、173.2 对 170.9。这种“擦边超标”是最有希望用最小代价修好的。观察这两条路径1→5 是 1-2-10-9-12-52→5 是 2-10-9-12-5后半段完全重合都经过 9-12 这条边。于是优化思路出来了把 9-12 这条边换成 9-5。9 号是 A(50,75)12 号是 D(115,70)9-12 长 65.19 米9 号到 5 号 P5(120,100) 长 74.33 米。换边之后总长多出 9.14 米但两条超标的路径都被截短1→5 变成 1-2-10-9-5长度 30 41.2 36.4 74.33 181.9 米小于 198.02→5 变成 2-10-9-5长度 41.2 36.4 74.33 151.9 米小于 170.9。顺便还得确认另外 8 对没有因为删掉 9-12 而变长。这一步不能省换边是全局操作动一条边会影响经过它的所有路径。检查完发现其余 8 对的最短路径都不经过 9-12或者即便经过也有替代走法且不超标换边才成立。最终内部道路总长 394.5 米。这段推导里有一个容易被忽略的点局部优化的前提是“找出超标路径的公共边”。如果两条超标路径没有公共边换边就得一次动两处代价和复杂度都不一样。实际做题时可以先按公共边排序优先处理重合度最高的那几条。4. 避坑与排查这套流程里翻车最多的五个地方前面两章是顺利路径这一章写不顺利的部分。下面五条都是我在复现类似模型时真实撞过的每条按“看到什么现象、根子在哪、怎么改”来说能省掉不少通宵。4.1 沿边距离表和程序算的对不上差了十几米现象很典型手算出来 P3 到 P4 沿边走 60 米程序跑出来 90 米两边对不上后面的筛选结论跟着全乱。原因几乎总是周长参数化的起点或者方向不一致——有人从 P1 出发顺时针累计有人从 P8 出发逆时针累计还有人把上边的 x 方向搞反了从上边右端往左走x 是递减的。解决办法是把参数化过程写成不可辩驳的代码而不是靠人肉数格子。具体做法是先确定一个遍历顺序我一般用顺时针下边从左到右、右边从下到上、上边从右到左、左边从上到下然后写一个坐标到弧长的映射函数把 8 个入口的 t 值打印出来人工核对一遍对上了再往下走。这一步花五分钟能省掉后面两小时的重跑。还有一种更隐蔽的情况入口恰好落在矩形角点上。这道题里没有角点入口但如果换一道题遇到t会出现两个候选值角点在两条边上都算必须显式规定取哪一个否则沿边距离会算成 0判定全错。4.2 直接把最小生成树当最终方案交上去现象Prim 跑完图画出来挺好看总长也确实不大但一校验就有好几对入口超标而且超得很离谱。原因在于优化目标和约束目标错位——最小生成树最小化的是 11 条边的权重和它对“树上最远两点有多远”这件事完全不敏感。一棵很瘦长的树边权和可以很小但直径可以很大。这个问题不是靠换算法解决的Prim 换成 Kruskal 结果一模一样。正确的做法是把最小生成树当成起点而不是终点先用它拿到一个下界和一个初始拓扑再用 Floyd 算出所有入口对的实际走行距离把超标的挑出来针对性地加边、删边、挪点。整个过程是一个“构造—校验—修补”的循环而不是一次求解。问题二里这个现象更明显8 个入口的最小生成树跑完7 对入口超标局部修补之后总长从原始树长涨到 427.5 米。如果一开始就把树长当成答案交上去数字会好看很多但方案根本不成立。4.3 Fermat 点用 fmincon 跑飞落在三角形外面现象明明三角形三个内角都小于 120°理论上 Fermat 点应该在三角形内部但fmincon给出的解跑到几十米开外目标函数值反而更大。原因通常有两个一是初值给得太随意直接从原点或者某个顶点起步sqp算法在非凸的目标函数上滑到了局部极小二是三角形本身有一个内角大于等于 120°这时候 Fermat 点的定义直接退化成那个钝角顶点根本不该用数值优化去求。改成两步走就稳了。第一步先算三个内角只要有一个 ≥ 120°直接返回对应顶点不走优化器第二步用重心做初值同时加一个多点重启取目标函数最小的那个解。function [q, fval] fermat_point(A, B, C) V [A; B; C]; ang zeros(3,1); for k 1:3 p V(k,:); u V(mod(k,3)1,:) - p; v V(mod(k1,3)1,:) - p; ang(k) acosd(dot(u,v) / (norm(u) * norm(v))); end [mx, idx] max(ang); if mx 120 q V(idx,:); fval sum(vecnorm(V - q, 2, 2)); return; end obj (x) sum(vecnorm(V - x, 2, 2)); opts optimoptions(fmincon, Display, off, Algorithm, sqp, ... OptimalityTolerance, 1e-10, StepTolerance, 1e-10); cand [mean(V,1); (AB)/2; (BC)/2; (CA)/2]; % 四个初值兜底 best inf; q mean(V,1); for k 1:size(cand,1) [xk, fk] fmincon(obj, cand(k,:), [],[],[],[], [0 0], [200 100], [], opts); if fk best best fk; q xk; end end fval best; endacosd返回的是角度制和 120 直接比省掉一次pi/180换算。上下界[0 0]和[200 100]是公园的矩形范围把搜索空间限死在场地内也能挡掉一部分跑飞。四个初值分别取重心和三条边的中点代价是四次优化对这道题的规模完全无所谓换取的是稳定性。4.4 把湖的四个角当成普通节点现象问题三加入矩形湖之后把 R1 到 R4 四个顶点也编进节点表一起送进 Prim结果生成的路径直接横穿湖面。原因是对“障碍”的建模错了——湖的四个角是障碍物的边界点不是可以随便打结的交叉点。在这类问题里允许新修道路与湖的边界相接题目说了“可以到达湖四周的边”但不允许穿过湖的内部。正确的处理是在距离矩阵层面把“穿越湖体”的边权重置成INF让 Prim 和 Floyd 在选边时自动避开。判断一条边是否穿湖用的是线段与轴对齐矩形的相交检测只要线段和矩形四条边中的任意一条真正相交端点接触不算就认为是穿越。function hit seg_hit_rect(p, q, rect) % rect [xmin ymin xmax ymax]判断线段 pq 是否穿过矩形内部 hit false; corners [rect(1) rect(2); rect(3) rect(2); rect(3) rect(4); rect(1) rect(4)]; for k 1:4 a corners(k,:); b corners(mod(k,4)1,:); if seg_cross(p, q, a, b) % 严格相交端点接触不算 hit true; return; end end end这个函数本身很短真正的坑在seg_cross的实现上用叉积判方向的时候必须把共线和端点重合的情况分开处理否则“道路贴着湖边走”会被误判成穿湖白白丢掉一批合法解。4.5 浮点容差吃掉一次合规判定现象某一对入口的路径长度算出来是 198.0000000001上限是 198.0程序判定不合格于是多做了一次优化总长无谓地增加了。原因就是纯浮点比较a b在两者数学上相等时返回 false。这道题里最危险的是 1→8 和 3→4 这两对沿边距离和 1.4 倍上限只差 0.2 到 0.4 米虽然不是浮点级别但如果坐标或者参数有一点偏差就会翻盘。稳妥的做法是给所有合规判定加一个绝对容差比如len 1.4*d 1e-6并且在打印结果时保留 3 位小数让人眼能复核。另一个相关问题是单位。整道题都是米但如果有人在中间某一步用了“格子”做单位或者把厘米混进来所有数量级全乱。我一般在脚本开头就把单位写死成注释所有函数只接受米。5. Fermat 点枚举从 427.5 米压到 358.4 米的完整推演问题一里交叉点是给定的问题二把约束放开交叉点位置变成决策变量搜索空间一下大了一个维度。直接上全局优化不现实可行的路子是“先构造一个合规的初始拓扑再用 Fermat 点逐个替换三角形结构每替换一次比较一次总长留下更短的那个”。5.1 八个入口的初版拓扑与局部修补先只拿 8 个入口建最小生成树不加任何内部节点。跑完 Floyd 校验7 对入口超标入口对树上路径长度1.4 倍上限超标幅度1→5217.7198.019.72→5247.7170.976.82→6162.7141.621.13→5158.3150.87.53→6243.3222.121.23→7268.3252.415.9超标的都以 5、6、7 这三个上边入口为端点问题集中在上半部分。局部修补的思路和问题一一样找公共边、加替代边、删多余边。具体做了四步先试连 1-6发现 2→5 会变成 226.1 米、超过 170.9 的上限放弃改连 2-62→6 直接合规但 2→5 还是超标于是补一条 2-58-7 那条直连删掉因为 8 到 7 走边界只有 50 多米比直连还短3-4-5 这条链超标删掉 4-5改连 3-5这样 3→6 走 3-5-6、3→7 走 3-5-6-7 都能压到上限以内。修补完总长 427.5 米。这个数字是问题二的第一版可行解也是后面插入 Fermat 点的比较基准。5.2 Fermat 点的插入判据省下来的必须比加进去的多一个三角形三条边如果能在内部找一个点 Q使得 QA QB QC 小于原来三条边之和那把这个三角形换成“Q 连三个顶点”的星形结构就划算。Fermat 点就是这个使距离和最小的 Q。但“划算”这件事必须每次现算不能拍脑袋。判据写成代码就一行function gain fermat_gain(A, B, C) % 正数表示插入 Fermat 点后总长会缩短值得插入 [q, fval] fermat_point(A, B, C); base norm(A-B) norm(B-C) norm(C-A); gain base - fval; end三角形越“胖”省得越多三角形越扁省得越少扁到某个程度就变成负收益。前面 4.3 节处理的内角 ≥ 120° 退化成顶点的情况就是扁到极致的表现。问题二里2、3、5 三个入口构成一个比较胖的三角形插入一个 Fermat 点后总长从 427.5 米降到 392.9 米省了 34.6 米。这个降幅相当可观说明初版拓扑里那三条边确实绕了远路。5.3 一点、两点、三点枚举为什么三个反而更差既然一个点有效自然会想到试两个点。两个 Fermat 点的做法是在初版拓扑里找两个不相交的三角形各自插入一个点。这里选了 △P2P5P6 和 △P3P4P5分别求 Fermat 点得到 Q2 和 Q3。跑完总长 358.4 米比单点又降了 34.5 米。再往下试三个点把 Q1、Q2、Q3 同时作为交叉点插入结果总长 416.9 米比两点方案还长。这个结果乍看反直觉其实很合理第三个三角形在当时的拓扑下已经不够“胖”插入 Q1 之后把原本一条较短的边拆成了三条辐条反倒增加了长度。用 5.2 节的fermat_gain单独算一下那个三角形的收益大概率是负的。把三种方案放一起对比方案交叉点数量总长米初版局部优化0427.5插入 Q11392.9插入 Q2、Q32358.4插入 Q1、Q2、Q33416.9所以最终答案是两个交叉点总长 358.4 米。这里的方法论意义比具体数字更重要Fermat 点的数量不是越多越好每插一个都要独立算收益收益为负的直接拒绝。贪心地“能插就插”在这道题上会得到一个比两点方案差 58 米的结果。还有一点值得说两个 Fermat 点的解不是从三个点的解里删出来的而是从初版拓扑里独立构造的。如果先构造三点方案再删点很可能删不到 358.4 这个局部最优因为三点方案的整体拓扑已经和两点方案不一样了。6. 湖面障碍与收尾359.6 米怎么守住并验证第三问加了一个 25 米 × 25 米的矩形湖四个角点 R1(140,70)、R2(140,45)、R3(165,45)、R4(165,70)。问题二那条 358.4 米的最优路径恰好从这个位置穿过去所以必须先重排。重排的枚举方式是绕行方向路径要么从湖的左侧过要么从右侧过。左侧方案算下来 362.8 米右侧方案 397.95 米差距 35 米实际上就是右侧要多绕一个湖的宽度。取左侧方案作为新的基准。在这个基础上再对 △P3P4R4 求一次 Fermat 点总长压到 360.6 米。最后用斯坦纳树的方法对整体拓扑做一次收口得到三个交叉点 h1(68.825, 64.7)、h2(112.216, 75.34)、h3(161.8, 30.7)总长 359.6 米。这个数字比无障碍时的 358.4 米只多了 1.2 米说明湖的位置恰好落在原路径的边缘附近绕行代价很小。验证这一步千万别省。把最终路径的每一条边都过一遍seg_hit_rect确认没有一条穿湖再把 8 个入口的 28 对距离全部跑一遍 Floyd确认没有一对超过 1.4 倍最后检查每条边的端点是不是都落在 8 个入口或者选定的交叉点上别出现“路修到边界中间某个位置”的违规情况。验证项检查方法通过标准不穿湖seg_hit_rect逐边检测全部返回 false1.4 倍约束Floyd 全源距离对比28 对全部满足边界落点端点集合比对只落在 P1–P8节点连通邻接矩阵连通性检查8 个入口互相可达这套检查我现在的习惯是写成一个函数本文还有配套的精品资源点击获取