资讯动态

超市建模程序解析:从OD矩阵到整数规划选址

发布时间:2026/10/22 8:08:21 来源:尧图企业网站定制
简介一份面向数学建模学习者与参赛者的rar压缩包围绕2004年全国大学生数学建模竞赛A题“奥运会临时超市网点布局优化”提供完整MATLAB解题程序。程序涵盖数据处理、选址模型构建、优化求解与结果输出可用于理解超市选址中的需求预测、网络优化及约束条件设置等关键环节。压缩包仅含1个m文件文件大小1KB内容精简但核心逻辑完整适合已掌握MATLAB基础、希望研读建模实战代码的读者参考。目前已有154人学习下载。通过阅读和运行这份代码可以直观掌握如何将地理信息、交通流等因素转化为可求解的数学规划模型并借助MATLAB实现从建模到仿真的全过程对准备数学建模竞赛或零售网点布局问题研究均有启发。1. 从大型赛事临时超市到日常网点这份超市建模程序能解决什么问题如果你也在头疼超市建模问题那这包资源算是踩中了要害。某年大学生数学建模竞赛A题把“一座大型赛事场馆周边该放几个临时超市、放在哪、配多大面积”变成了一整套可计算的数学题而这个压缩包里装的就是那份赛题的MATLAB完整解题程序。别指望打开就能跑出最优解真正拉开差距的是前两步怎么把抽样问卷扩成全量出行矩阵怎么把步行路径折算成每个候选点的潜在需求。这套源码适合两类人——备赛的学生以及做零售网点选址但不想一直靠经验拍板的从业者。超市建模问题看着是选址题本质是流量分配加整数规划。2. 数据预处理与OD矩阵构建把问卷表变成稀疏矩阵2.1 从题面到建模假设先回答三个问题拿到题面第一件事不是找数据文件而是把文字描述变成三个能写进代码的量。第一个是需求点的划分场馆观众席、出入口、商业街连廊被切成若干交通小区每个小区的质心就是一个需求点需求量用“全天从该小区出发去往其他小区的人次”表示。第二个是候选点的集合临时超市不能随便放只能放在允许搭建的商业通道节点上这些节点就是候选点。第三个是流动路径题目给的问卷数据是抽样调查不是全量统计需要按抽样权重放大否则后面所有数字都偏小一个量级。这套源码在开头部分就体现了这个顺序先是几个注释块明确变量作用域再做数据导入。如果自己重写建议照抄这个结构——先定义n_zone、n_cand、Rmax三个全局量再往下展开。很多人拿到代码就直接跑结果后面想改一个参数要全局替换非常痛苦。问卷数据通常长这样每一行是一次被访者的出行记录列包括起点分区编号、终点分区编号、交通方式、消费档位、周出行频次、抽样权重。抽样权重是本阶段最关键的列比如调查了总人流的5%那么每一条记录要放大20倍才代表全量。把这个放大做在数据导入阶段而不是做到模型里后面算覆盖、算收益才不会乱。2.2 用accumarray快速构建OD矩阵OD矩阵是所有后续计算的基石。它的行是出发小区列是到达小区单元格的值是全天从i到j的人次。用MATLAB的accumarray可以一次性把几千行问卷聚合出这个矩阵不需要写三重循环。% 读取原始抽样表分区编号、交通方式、出行频次、抽样权重 raw xlsread(survey_data.xlsx); from_zone raw(:,1); % 出行起点小区编号 to_zone raw(:,2); % 出行终点小区编号 samp_wt raw(:,6); % 抽样比倒数用于扩样成全天总量 n_zone max([max(from_zone), max(to_zone)]); OD accumarray([from_zone, to_zone], samp_wt, [n_zone, n_zone]); % 剔除对角线区内出行不产生经过商业街的跨区流 OD(1:n_zone1:end) 0; figure; imagesc(OD); colorbar; title(全天跨区出行OD);这段代码的核心是accumarray([from_zone, to_zone], samp_wt, [n_zone, n_zone])。第一个参数是下标矩阵每一行的两个值决定结果矩阵里的行列位置第二个参数是对应值如果有多条记录落在同一个(i,j)位置accumarray会自动累加。第三个参数指定输出尺寸防止出现编号正好缺一位导致矩阵维度错位。这里最容易踩的坑是编号必须是从1开始的连续正整数。原始问卷里如果起点编号是0、100、200这种跳号accumarray会直接报错或生成超大矩阵。比赛数据里经常出现编号带着楼栋号、看台号一起编的情况所以要先做一次编号映射把原始编号压缩成1:n的连续序列% 将非连续分区编号映射为1:n all_zone [raw(:,1); raw(:,2)]; [~, ~, ic] unique(all_zone); n_zone max(ic); from_zone ic(1:size(raw,1)); to_zone ic(size(raw,1)1:end);unique返回的ic是每个元素在去重列表里的位置天然就是1开头的连续整数。这套映射做完之后OD矩阵的行列含义才真正对齐。对角线置零这一步也别省区内出行不会从看台绕到商业街不剔除的话后面需求归集会虚高。2.3 把距离矩阵算出来再谈选址OD矩阵只回答了“有多少人从哪到哪”还没回答“这些人是否会经过某个候选超市”。要把两者连接起来需要构造一张路网邻接矩阵算每个需求点到每个候选点的最短步行时间得到距离矩阵dist。dist(i,j)表示第i个需求点到第j个候选点要走的分钟数这个矩阵就是下一章的主角。很多人在这一步直接把欧氏距离当步行时间这种近似在图纸上看着合理实际上会漏掉“商业街连廊必须从看台出口绕行”这种关键约束。后面所有模型结论都会被一个错误的距离矩阵带偏而且很难察觉因为最优解的图看起来还挺漂亮。我一般建议至少把道路拓扑表读进来跑一遍最短路径再决定要不要用简化距离。3. 人流模拟与需求预测用最短路径把出行量折算成潜在需求3.1 为什么要先模拟人流而不是直接套统计回归数模评审里经常能看到一种做法把各小区常住人口密度和消费能力直接回归出需求量然后用这个需求量做选址。这个思路在常态社区商业里说得通但在大型赛事临时超市场景下是错位的——观众不是居民他们不会专程去某家超市而是从看台散场时顺路经过。所以需求量的正确来源是“路径流量”不是“人口密度”。这就要把人流模拟放在选址之前先建立路网计算所有OD对的最短路径再把每条路径上的流量按照“是否会经过某个候选点”归集到各候选点身上。这样得到的demand向量才是超市建模问题里真正该放进优化模型的需求量它天然带着动线信息后面做“如果封闭某个出口”的情景分析时也能直接改路网重算。3.2 路网邻接矩阵与最短路计算路网数据一般给的是拓扑表每条边三个字段端点A、端点B、步行时间。构造邻接矩阵时有一个最常见的错误——直接用W(a,b)t赋值。如果拓扑表里同一条边出现两次比如双向道路各占一行后面的值会覆盖前面的而且两条平行路径也没法取最小值。正确做法是初始化全Inf赋值时取min。% 路网拓扑表三列分别为端点A、端点B、步行时间分钟 edge xlsread(road_network.xlsx); n_node max(max(edge(:,1:2))); % 构造无向图邻接矩阵不连通的节点之间保持Inf W inf(n_node, n_node); for k 1:size(edge,1) a edge(k,1); b edge(k,2); t edge(k,3); W(a,b) min(W(a,b), t); W(b,a) W(a,b); end W(1:n_node1:end) 0; % 对角线为0表示原地点到点不需要时间 % 需求点到候选点的最短步行时间 dist zeros(n_zone, n_cand); for i 1:n_zone for j 1:n_cand [~, d] graphshortestpath(sparse(W), zone_node(i), cand_node(j)); dist(i,j) d; end endgraphshortestpath是旧版MATLAB的图论函数新版建议换成shortestpath(graph(W), s, t)用法类似返回值里的距离含义不变。把W转成稀疏矩阵是因为路网节点多但边少稀疏格式能把计算时间砍掉一个量级。这段代码里的zone_node和cand_node是两个映射数组需求点质心在路网里的节点编号、候选超市所在的路网节点编号。当需求点正好落在节点上时最理想不在的话要找最近节点做投影这个投影逻辑通常由距离判断完成。3.3 需求量归集就近分配与距离衰减两种口径拿到dist矩阵之后需求归集就变成一个分配问题。最简单的口径是“就近分配”每一条OD出行从其起点出发在所有步行可达的候选点里选最近的那个把这段出行人次加到该候选点的需求量上。Rmax 15; % 最大步行半径单位分钟 demand zeros(1, n_cand); for i 1:n_zone for k 1:n_zone if OD(i,k) 0 continue; end % 从起点小区i出发找所有步行可达候选点 reach_cand find(dist(i,:) Rmax); if isempty(reach_cand) % 所有候选点都太远视为流失需求 continue; end [~, idx] min(dist(i, reach_cand)); j reach_cand(idx); demand(j) demand(j) OD(i,k); end end这个逻辑把每条i→k的出行按出发点i归属到离它最近的候选点隐含的假设是“观众从看台出来第一眼看到哪家超市就去哪家”。Rmax15是步行极限的常见经验值比赛时可以扫一组值做敏感性分析不要拍脑袋定死。更细的做法是用距离衰减代替硬截断把OD(i,k)按1/(1dist(i,j)/d0)的比例分摊给所有可达候选点其中d0是衰减尺度一般取5到10分钟。距离衰减口径更接近真实消费行为但解释起来比就近分配费劲。我自己的习惯是两种口径都跑一遍如果最优选址结果差异不大说明模型稳健如果换了口径点位就完全变了那要警惕是不是距离矩阵本身有硬伤。3.4 把中间结果落盘方便后续回代OD、dist、demand、Rmax这四个量是整套流水线的中间产物后面做敏感性分析、覆盖回代都要反复用到。建议在源码里加一行保存save(mid_result.mat, OD, dist, demand, Rmax, n_zone, n_cand);这样做的好处是后面改选址模型参数时不用重新读问卷、重新算最短路直接load回来就能接着跑。比赛里时间紧这步能省下大量重复计算。4. p-中位选址模型与整数规划用intlinprog锁定超市点位4.1 为什么是p-中位而不是p-中心或覆盖模型网点选址的经典模型有三个p-中心、集合覆盖、p-中位。p-中心的优化目标是让最远的需求点到最近设施的距离最小适合消防站、急救中心这种强调公平性的场景集合覆盖的目标是最少设施数覆盖全部需求点适合应急资源布局而p-中位让所有需求点到最近设施的加权总距离最小本质是“贴近大多数人的动线”。临时超市选址属于典型的p-中位场景。超市不需要保证每个观众都在5分钟内到达但需要尽量多地拦截主流人流让加权总步行时间最小。目标函数长这样每个需求点i的需求量h_i去往它被分配的候选点j的步行时间d_ij把所有组合的人次乘分钟累加起来作为总代价。这个代价越小说明超市越贴近大多数观众的行动路径。理解了为什么选p-中位你就能明白为什么题目给的压缩包里没有用模拟退火或遗传算法——因为决策变量全是0-1整数约束是线性的本身就是一个标准的整数规划问题直接用分支定界求解器能拿到全局最优解启发式算法在这个规模下既不必要也不可靠。4.2 intlinprog求解p-中位变量展平与约束构造用MATLAB的intlinprog解这个模型最麻烦的不是建模本身而是把二维的分配变量x_ij展平成求解器需要的一维列向量。这里约定决策变量排列顺序是先放所有x_ij按列优先再放所有y_j所以变量总数为n*m m。function [x_opt, y_opt, fval] solve_pmedian(demand, dist, p) % demand: n_zone×1 需求向量列向量 % dist: n_zone×n_cand 距离矩阵单位步行分钟 % p: 计划建设的超市数量 n size(dist,1); m size(dist,2); h demand(:); % 统一为列向量 % 决策变量排列先n*m个x按列展平再m个y f [reshape(dist .* h, [], 1); zeros(m, 1)]; intcon 1:(n*m m); % 约束1每个需求点只能分配给一个候选点 Aeq [kron(eye(n), ones(1,m)), sparse(n, m)]; beq ones(n, 1); % 约束2x_ij y_j只有被选中的候选点才能承接需求 A2 sparse(n*m, n*m m); row 0; for j 1:m for i 1:n row row 1; A2(row, (j-1)*n i) 1; A2(row, n*m j) -1; end end b2 zeros(n*m, 1); % 约束3恰好选中p个候选点 Aeq2 [sparse(1, n*m), ones(1, m)]; Aeq [Aeq; Aeq2]; beq [beq; p]; % 0-1变量边界 lb zeros(n*m m, 1); ub ones(n*m m, 1); opts optimoptions(intlinprog, Display, iter); [x, fval] intlinprog(f, intcon, A2, b2, Aeq, beq, lb, ub, opts); x_opt reshape(x(1:n*m), n, m); y_opt x(n*m1:end); end目标函数里dist .* h利用隐式扩展把需求向量广播到每一列再按列展平顺序和x的排列严格对应。这里一旦展平顺序搞反结果就是错位的最优解表面上exitflag1实际上完全不能用。约束1用kron(eye(n), ones(1,m))生成块对角矩阵每一行对应一个需求点保证每个需求点必须且只能分配给一个候选超市。约束2的循环构造是给每一个(i,j)对写一条x_ij - y_j 0只有y_j1的候选点才能承接需求。约束3把选中数量锁死在p。循环构造约束在n200、m50时会生成一万条约束虽然能用但效率不高。建议先写成循环方便核对逻辑确认无误后再用稀疏矩阵向量化拼接。比赛里时间有限正确性高于优雅性。Displayiter会把分支定界的迭代过程打到命令行如果看到节点数持续暴涨说明变量规模太大需要回头削减候选点数量。4.3 求解结果怎么用求解结束后y_opt里值为1的下标就是最终选出的超市点位x_opt(:,j)给出哪些需求点被分配给了第j个超市。要把结果变成报告里能直接用的表格需要把每个选中超市承接的需求量算出来sel find(y_opt 0.5); h demand(:); for s 1:length(sel) j sel(s); volume sum(x_opt(:, j) .* h); % 分配给该超市的总需求量 fprintf(候选点%d 承接需求 %.0f 人次/天\n, j, volume); endx_opt(:,j) .* h是把分配到第j个超市的所有需求点需求量加起来得到这个超市的“服务量”。这个值可以直接换算成超市面积常见做法是按人次和客单价估算营业额再按每平米日销售额反推面积。源码里如果只给出了点位没给面积这一步就是你自己补全的部分。5. 复现Sams-副本.m的避坑指南四个最容易翻车的细节5.1 拿到源码先别跑先确认这四个状态这份压缩包里的Sams-副本.m是完整流程的串联但直接双击运行大概率不会一次通过。我建议动代码之前先检查四件事。第一数据单位是否统一。源码里距离是分钟OD是“人次/全天”。如果有其他脚本把某个表换成了小时或千米最优解会彻底变样。这个错误最隐蔽因为图照样能画出来曲线照样长得像模像样但所有数字都差了一个倍率。第二编号是否从1开始且连续。accumarray和kron都要求下标是连续正整数如果原始分区号是0、100、200这种必须先映射。跳号导致的不是报错就是矩阵维度对不上排查起来很费劲。第三求解器的exitflag是否等于1。intlinprog返回exitflag0或负值意味着没拿到最优解甚至模型不可行。这时候不要急着看x_opt先把约束矩阵打出来核对。第四可视化叠加的坐标系。很多人在最后画图时把需求点质心、候选点、路网画在三张图里看起来叠上了实际缩放比例不一致导致点位错判。用同一套坐标轴画或者直接用hold on叠加。5.2 常见问题排查记录下面几条是复现这类源码时最高频的问题按“现象→原因→解决”记录。现象一intlinprog返回exitflag0或者直接提示约束不可行。原因约束2里x_ij y_j的行列对应关系写错导致某个需求点在可达半径内没有任何被选中的候选点另一个常见原因是dist矩阵中存在整行全是Inf也就是某个需求点根本到不了任何候选点。解决先检查dist的每一行是否有小于等于Rmax的值用find(min(dist,[],2)Inf)快速定位。再用m3, n3的小样例手算A2的某一行确认(j-1)*ni这个索引与reshape的展平顺序一致。现象二最优解里某个超市承接的需求量是0。原因需求归集用了Rmax硬截断把该候选点附近的所有OD都判成了不可达。或者是demand归集时dist(i,:)用的下标范围不对候选点编号和dist列索引错位。解决打印该候选点到所有需求点的距离分布看阈值是否设置得太小。换距离衰减口径再跑一次如果结果差异超过20%说明模型对阈值敏感需要在报告里做敏感性分析。现象三最短路径矩阵出现Inf或者对角线不是0。原因邻接矩阵初始化时用了0而不是Inf导致不连通的节点之间“距离为0”或者n_node算小了一个号最大节点编号的边被截断。解决初始化全部为Inf赋值时取min(W(a,b), t)。n_node用max(max(edge(:,1:2)))不要用length(unique())后者在节点编号不连续时会直接算错。现象四脚本跑了十几分钟还没出结果。原因决策变量数量为n*mm当需求点200个、候选点50个时整数变量超过一万个分支定界节点数爆炸。这不是代码写错了是规模超出了求解器的舒适区。解决先把候选点数量压到30以内方法是对候选点按需求热度聚类每组只保留质心。或者先用LP松弛跑一遍看下界再决定值不值得用启发式。这是比赛里最常见的血泪经验——不要一上来就上全量先用小规模验证算法正确再逐步扩大。6. 把结果从图做成结论敏感性分析与方案回代6.1 用P曲线确定“到底放几个超市”intlinprog的输入p是超市数量但题目通常不会直接告诉你放几个最合理。这个值要靠敏感性分析扫出来。把p从1取到10逐个跑一遍solve_pmedian记录每次的总加权步行时间fval画一条随p变化的曲线p_list 1:10; fval_list zeros(length(p_list),1); for k 1:length(p_list) [~,~, fval_list(k)] solve_pmedian(demand, dist, p_list(k)); end plot(p_list, fval_list, -o); xlabel(超市数量 p); ylabel(加权总步行时间人次·分钟);这条曲线通常先是快速下降到某个值之后边际收益明显变缓这个拐点就是推荐数量。比单点结果有说服力得多。写报告时把这张图放进去比只贴一张选址地图更有讨论价值。6.2 回代验证选中点位之后重新核算覆盖比例选完点位后还有一个习惯必须养成把选址结果放回OD矩阵做覆盖回代。做法是固定sel重新遍历所有需求点看每个需求点距离最近的选中超市是否在Rmax以内统计覆盖了多少总出行量sel find(y_opt 0.5); cover 0; for i 1:n_zone if min(dist(i, sel)) Rmax cover cover sum(OD(i,:)); end end fprintf(覆盖比例 %.1f%%\n, cover / sum(OD(:)) * 100);覆盖率是选址方案的最终验收指标。如果覆盖率很高但fval很大说明超市数量够但位置不够贴动人流如果覆盖率低就要回到需求归集那段检查Rmax。从那以后我每次复现这类选址代码都强制把“单位→矩阵维度→求解器状态→覆盖回代”四步走一遍。省掉任何一步后面所有结论都可能要推翻重来。这份压缩包里的源码值得你打开MATLAB一行行跑通再换一组自己的数据试试希望帮到你。本文还有配套的精品资源点击获取

读完文章,也想定制专属网站?

尧图设计师 24 小时内与您沟通定制方案

免费获取报价 →
↑