简介面向全国大学生数学建模竞赛参赛者及运筹优化学习者这份资源围绕2004年高教社杯CUMCM A题“奥运会临时超市网点设计”给出完整参赛解决方案。压缩包内共6个文件包含4个Matlab脚本.m、1个Access数据库.mdb和1个赛题文档.doc整体仅1.5MB轻量但内容紧凑同时兼顾建模代码、原始数据与题目说明。已有292人浏览学习适合备赛参考或复现经典赛题。资源不仅提供可运行的求解程序还覆盖线性规划、遗传算法、模拟退火等关键优化思想并结合需求预测、选址建模、结果可视化等环节一步步还原从原始数据到最终布局方案的建模过程帮助读者深入理解运筹优化模型的构建与Matlab实现切实提升数学建模与编程实战能力。1. 从CUMCM2004A题看Matlab求解临时超市选址的完整路径拿到高教社杯全国大学生数学建模竞赛(CUMCM2004A题)时最常见的误区是立刻在网上搜“奥运会临时超市网点设计Matlab源码”然后直接改数跑结果。真正让这道题从“解出来”变成“能获奖”的是把题目给出的观众小区分布、场馆坐标、购物概率和超市容量约束转成一组可以被Matlab优化工具箱处理的稀疏矩阵和整数规划模型。这道题表面是商业选址本质上是带容量约束的集合覆盖与需求分配问题。下面按“建模—数据预处理—求解—验证”的顺序给出一套用Matlab复现时可落地的代码路径包含每个参数为什么这么设以及替换数据后需要改哪些变量。2. 先把题面翻译成数学模型变量、目标函数与约束条件2.1 从“网点设计”到集合覆盖定义下标集合常见做法是把问题拆成两类实体需求点和候选网点。题目通常给出若干个观众居住小区或观众来源区每个小区有固定人数并给出每个小区到各个场馆的人流概率临时超市则安排在候选位置比如场馆周边的空地或道路节点。为了建立可求解模型需要先定义以下集合(S)观众需求点集合大小为 (N_s)每个点 (s) 有需求量 (d_s)。(J)候选超市网点集合大小为 (N_j)每个候选点 (j) 有最大服务容量 (c_j)。(K)场馆或热点区域集合用于把需求点与场馆之间的距离纳入目标。题目给的“人流数据”往往不是现成的需求点而是“小区—场馆”的OD量。所以第一步是把OD矩阵按概率折算成每个小时的总购物人数再聚合到需求点坐标上。这个过程不只是读Excel还要做空间重采样。在Matlab里建议用结构体或table保存原始数据而不是散落的工作区变量% 读取题目附带的观众小区数据 raw readtable(demand_points.csv); % 假设有两列坐标x,y和人数pop s.x raw.x; s.y raw.y; s.demand raw.pop; % 若人数是全天总量需要按题目给的购物时间窗口折算 N_s height(raw);这里把需求点直接存成结构体 (s)后面的距离计算和矩阵构建都基于这个结构体。注意需求量不是人口普查人数而是“有购物意愿且在该时段到场馆的人数”题目往往给出比例常见做法是乘以一个购物系数 (\alpha)这个系数在后面会作为敏感性分析的变量。2.2 目标函数覆盖需求最大化与总步行距离最小化CUMCM2004A题常见的评阅标准会让选手权衡两个指标一个是满足的购物需求总量另一个是观众到最近超市的平均距离。在数学上可以写成双目标但为了用Matlab的整数线性规划求解通常把其中一个转成约束或者线性加权。如果以“最小化总步行距离”为主目标可以写成[ \min \sum_{s \in S}\sum_{j \in J} t_{s,j} \cdot y_{s,j} \cdot d_s ]其中 (y_{s,j}) 是0-1变量表示需求点 (s) 是否被分配给候选网点 (j)(t_{s,j}) 是两者步行距离。同时需要引入另一个0-1变量 (x_j) 表示是否在候选点 (j) 建超市。如果以“最大化覆盖需求”为目标则目标函数为[ \max \sum_{s \in S}\sum_{j \in J} y_{s,j} \cdot d_s ]实际建模中我一般取前者因为它对“临时超市”场景更直观临时网点本身是服务设施观众不会跑太远。下面的推导都按最小化总距离展开但切换目标时只要改Matlab的f向量即可。2.3 三类核心约束容量、唯一分配与覆盖半径约束条件决定了模型能不能用线性规划解。第一类是每个候选超市的容量上限[ \sum_{s \in S} d_s \cdot y_{s,j} \le c_j \cdot x_j, \quad \forall j \in J ]这里 (c_j) 是最大服务人数和超市面积、平均购物时间、营业窗口时长相关。如果题目给出的容量单位不是“人/小时”需要先换算。第二类是每个需求点只能被一个已开放网点服务[ \sum_{j \in J} y_{s,j} 1, \quad \forall s \in S ]第三类是覆盖半径约束只有距离小于阈值的需求点才能被分配[ y_{s,j} \le a_{s,j} ]其中 (a_{s,j}) 是0-1可达矩阵由距离阈值 (\tau) 生成。虽然第三类约束可以由距离下限隐含但显式写出能减少变量组合让整数规划求解更快。这三类约束组合起来就是一个典型的带容量约束的选址-分配问题。因为决策变量 (x_j) 和 (y_{s,j}) 都是整数所以需要调用整数规划求解器。下面的表总结了所有符号和对应Matlab变量名符号含义Matlab变量备选值(S)需求点集合s.x,s.y题目给出(J)候选网点集合j.x,j.y网格剖分或题目给(d_s)需求点需求量s.demand人口×购物概率(c_j)候选点容量cap超市面积/人均面积(t_{s,j})步行距离dist欧氏距离×路网系数(x_j)是否建点x(j)0/1(y_{s,j})是否分配y(s,j)0/12.4 用Matlab构建模型矩阵的骨架代码把上面的模型转成intlinprog能接受的矩阵形式最关键的是把双下标变量 (y_{s,j}) 展成一维向量。常见做法是按列展开即idx(s,j) (j-1)*N_s s。下面这段代码构建目标函数和容量约束% 假设已有 dist 矩阵、 cap 向量、 demand 向量 N_s length(s.demand); N_j length(cap); % 决策变量顺序 x(1)..x(N_j) 为建点变量 y(s,j) 展开在后面 n_x N_j; n_y N_s * N_j; f zeros(n_x n_y, 1); % 目标距离 * 需求量只作用于 y 变量 for j 1:N_j for k 1:N_s idx_y n_x (j-1)*N_s k; f(idx_y) dist(k,j) * s.demand(k); end end % 容量约束 sum_s demand(s) * y(s,j) cap(j) * x(j) A zeros(N_j, n_x n_y); b zeros(N_j, 1); for j 1:N_j % 对应 x(j) 的系数是 -cap(j)因为 cap(j)*x(j) 移到左边为 -cap(j) A(j, j) -cap(j); for k 1:N_s idx_y n_x (j-1)*N_s k; A(j, idx_y) s.demand(k); end end这段代码的意义是用矩阵A同时表达容量约束和0-1变量的关系。注意A(j, j) -cap(j)是让不等式变为 (\sum demand \cdot y - cap \cdot x \le 0)。这样写比直接在循环里拼intlinprog的输入更容易检查错误。实际调试时我会先打印A(1,:)前10个非零元素确认idx_y映射没有错位。3. 数据预处理从原始人流数据到可分块的稀疏矩阵3.1 题目数据读入Matlab先做清洗再谈建模CUMCM2004A题通常附带Excel或CSV文件里面往往混合了文本表头、单位行和空行。如果直接readtable会出现NaN或字符串变量混入数值列。推荐的做法是先读成table再用ismissing找到缺测点最后按列类型转换。下面代码处理常见格式T readtable(CUMCM2004A_data.csv, VariableNamingRule, preserve); % 去掉全空的行 T rmmissing(T); % 把小区编号列转成字符串如果原来是文本 if iscell(T{:,1}) id string(T{:,1}); end % 取出坐标和人数列列名按实际文件调整 coords [T{:, x_coord}, T{:, y_coord}]; pop T{:, population}; % 有些题目会给出“购物概率”这里直接乘上去 shopping_prob 0.65; % 这个值来自题目附表或假设 demand pop .* shopping_prob;这里有一个容易踩的坑题目给的人流量可能是“全天总客流”但临时超市只在比赛前后营业3到4小时所以必须把全天量除以营业时段数否则容量约束会严重失配。我一般会先画出demand的直方图如果最大值超过候选点容量的几十倍就要检查是不是忘记换算时间窗口。3.2 坐标转距离矩阵pdist2、矩阵广播和循环的取舍需求点到候选点的距离矩阵是后续建模的核心。如果候选点数量不大比如几百个直接用pdist2最方便% 假设 demand_coords 是 N_s x 2 candidate_coords 是 N_j x 2 dist_matrix pdist2(demand_coords, candidate_coords, euclidean);pdist2计算的是直线距离实际步行距离需要乘以一个绕行系数常见取值在1.2到1.5之间。如果数据量达到几万pdist2可能内存不够这时用分块计算或者直接考虑用路网矩阵。三种方式的对比方法代码复杂度内存占用适用场景pdist2低(O(N_s N_j))数据量在 1e4 以内广播dist sqrt((x_s - x_j.).^2 ...)中同上但更底层需要额外处理维度时双重循环高只需两行数据内存极小的教学环境速度慢实际参赛时不必追求极致的距离精度因为题目给的坐标已经是平面投影坐标线性距离和真实步行距离差别不会影响网点数量级。重点是把dist_matrix转成0-1可达矩阵reachablecover_radius 800; % 单位米来自题目对“步行5分钟”的约束 reachable dist_matrix cover_radius;这里阈值是模型里最重要的一个超参数。后来做敏感性分析时可以用这个阈值作为横轴观察结果变化。3.3 生成候选网点网格剖分与坐标聚类题目不一定直接给候选点列表有时候需要自己生成。最常见做法是把场馆周边区域划分成正方形网格取每个网格中心作为候选点。这样做的好处是网点规模可控且能覆盖整个区域。下面是生成候选点的代码% 由场馆坐标给出区域边界 x_min min(stadium(:,1)) - 500; x_max max(stadium(:,1)) 500; y_min min(stadium(:,2)) - 500; y_max max(stadium(:,2)) 500; grid_step 300; % 每300米一个网格 [gx, gy] meshgrid(x_min:grid_step:x_max, y_min:grid_step:y_max); candidate_coords [gx(:), gy(:)];网格步长直接决定候选点数量和求解难度。如果步长从300改为200候选点数量会从约 (N_s) 级变成 (N_s \times 1.5^2)整数规划耗时通常指数上升。所以建议先跑粗网格快速验证模型正确性再细化网格。还有一种做法是用kmeans对需求点聚类把聚类中心作为候选点。它能保证候选点贴近人流密集处但缺点是会漏掉一些交通枢纽位置。如果题目中明确提到“可以在场馆出入口附近建点”那么网格法中再加入出入口坐标即可。3.4 预处理阶段最容易错的三个细节第一是坐标单位。有的题给经纬度有的给直角坐标如果直接混用距离会差几个数量级。处理方法是用deg2km或utm2deg转换但更简单的是先看题目配图里的比例尺把单位统一成“百米”或“米”。第二是重复点。需求点里可能有多个小区在同一个坐标如果不去重会在距离矩阵里产生大量全同行导致intlinprog出现退化。用unique(coords,rows)合并需求量累加。第三是demand为0的点。如果某小区没有购物需求保留它会白增加变量数量。我通常把demand 1e-6的行直接过滤掉这样在后期求解数据规模更小。4. 求解策略从intlinprog到遗传算法的切换4.1 线性整数规划用Matlab优化工具箱的intlinprog落地当模型只有线性目标和线性约束时优先用intlinprog。它来自Matlab优化工具箱内置分支定界算法对小规模问题变量数几千个能在几十秒内返回全局最优。调用时需要指定决策变量中哪些是整数。在CUMCM2004A题里所有决策变量都是0-1变量所以intcon 1:num_vars。intcon 1:(n_x n_y); lb zeros(n_x n_y, 1); ub ones(n_x n_y, 1); % 唯一分配约束每个需求点只被一个网点服务 Aeq zeros(N_s, n_x n_y); beq ones(N_s, 1); for k 1:N_s for j 1:N_j idx_y n_x (j-1)*N_s k; Aeq(k, idx_y) 1; end end % 调用求解器 options optimoptions(intlinprog, Display, iter, MaxTime, 120); [x_opt, fval, exitflag] intlinprog(f, intcon, A, b, Aeq, beq, lb, ub, options);这段代码里A,b沿用上一章的容量约束Aeq,beq实现唯一分配。运行后会得到x_opt向量前N_j个是建点标志后面是分配关系。注意exitflag不等于1时需要看x_opt中的NaN值通常是没有可行解而不是求解失败。4.2 模型太大时换成遗传算法ga找近似解如果候选点超过500个intlinprog可能会在分支定界时卡住。此时改用全局优化工具箱的ga来求次优解更实际。遗传算法的核心是写一个适应度函数输入建点方案输出总距离并对违反容量约束的方案施加惩罚。function total_cost fitness_func(selected) % selected 是 0/1 向量长度 N_j if sum(selected) 0 total_cost 1e8; return; end % 每个需求点找最近的开放网点 open_idx find(selected); total_cost 0; for s_idx 1:N_s dist_to_open dist_matrix(s_idx, open_idx); [mind, ~] min(dist_to_open); total_cost total_cost demand(s_idx) * mind; end % 容量惩罚开放网点总容量不够时加大惩罚 total_capacity sum(cap(open_idx)); if total_capacity sum(demand) total_cost total_cost 1e6 * (sum(demand) - total_capacity); end end调用ga时需要设置种群大小和代数。适应度函数里的1e6是一个很大的惩罚系数目的是让算法避开不可行解。注意ga默认求解最小化问题所以适应度值越小越好。下面表格列出ga的关键参数及其推荐值参数推荐值作用PopulationSize200种群越大越不容易早熟MaxGenerations300迭代上限看收敛曲线CrossoverFraction0.8交叉比例太高会破坏好解MigrationFraction0.1子种群间迁移率Displayiter打印每代最优值4.3 整数规划与遗传算法的分工建议从工程角度我一般分两步走先用intlinprog在放松网格密度下求一个精确解作为上界参考再把这个解转成initial_population喂给ga让遗传算法在更细的网格上继续搜索。这样比直接让ga从随机解开始快很多。需要注意ga的初始种群必须是矩阵一行一个个体。可以把intlinprog的x_opt(1:N_j)作为第一个个体其余个体用随机01向量。这样遗传算法在早期就能集中在高质量解附近。4.4 求解失败的排查清单遇到intlinprog提示“No feasible solution”先检查三处第一所有需求点的总需求是否大于所有候选点容量之和如果是问题本身无解需要增加候选点或放宽容量第二覆盖半径是否太小导致某些需求点不在任何候选点覆盖范围内打印sum(reachable,2)可以找到那些全0行第三容量约束矩阵A是否写反了不等号方向。如果ga收敛很快但结果明显不合理往往是惩罚系数太大或太小。太大导致算法只关心可行性太小则容易产生容量超限的解。常见的做法是让惩罚系数为容量超限量的1000倍然后看多次运行是否稳定。5. 结果验证与网点调整的实用技巧5.1 用热力图检查覆盖盲区求解完成后把需求点和网点画在一起颜色表示需求量圈表示服务半径可以立刻发现哪些区域没有被覆盖。下面的代码画出覆盖热力图figure; scatter(demand_coords(:,1), demand_coords(:,2), 30, demand, filled); hold on; open_j find(x_opt(1:N_j) 0.5); viscircles(candidate_coords(open_j,:), cover_radius * ones(size(open_j)), ... Color, red, LineWidth, 1.2); colorbar; xlabel(x 坐标); ylabel(y 坐标); title(需求点与开放超市覆盖范围);如果发现某个需求点被多个圆圈覆盖说明容量分配可能不均匀。可以通过查看Aeq对应的y变量来检查每个开放网点的实际负荷并把手动调整候选点坐标作为下一步的输入。5.2 用敏感性分析选定最优覆盖半径CUMCM2004A题评分时很看重对参数的讨论。最直接的做法是把覆盖半径从600米逐步改到1200米步长100米每次调用同一个求解脚本记录总步行距离和开放网点数。把这个循环写成脚本最后用plot画折线图。如果曲线在某一半径后变得平缓那这个点就是合理的阈值。5.3 把求解结果整理成可直接放进论文的表格最后一步是把x_opt转成表格输出选中网点的坐标和服务人数占比。可以用writetable生成Excel方便在论文里直接截图selected_j find(x_opt(1:N_j) 0.5); out_table table(candidate_coords(selected_j,1), candidate_coords(selected_j,2), ... cap(selected_j), VariableNames, {x, y, capacity}); out_table.served zeros(height(out_table), 1); for jj 1:length(selected_j) out_table.served(jj) sum(demand(find_allocated_to_j(selected_j(jj)))); end writetable(out_table, resolved_networks.xlsx);这里find_allocated_to_j是从y变量中还原每个网点服务了哪些需求点的小函数。整个流程跑通后换一套数据只需要改第2章的readtable文件名和第3章的cover_radius。本文还有配套的精品资源点击获取 Matlab求解CUMCM2004A临时超市选址：整数规划与优化实践
简介面向全国大学生数学建模竞赛参赛者及运筹优化学习者这份资源围绕2004年高教社杯CUMCM A题“奥运会临时超市网点设计”给出完整参赛解决方案。压缩包内共6个文件包含4个Matlab脚本.m、1个Access数据库.mdb和1个赛题文档.doc整体仅1.5MB轻量但内容紧凑同时兼顾建模代码、原始数据与题目说明。已有292人浏览学习适合备赛参考或复现经典赛题。资源不仅提供可运行的求解程序还覆盖线性规划、遗传算法、模拟退火等关键优化思想并结合需求预测、选址建模、结果可视化等环节一步步还原从原始数据到最终布局方案的建模过程帮助读者深入理解运筹优化模型的构建与Matlab实现切实提升数学建模与编程实战能力。1. 从CUMCM2004A题看Matlab求解临时超市选址的完整路径拿到高教社杯全国大学生数学建模竞赛(CUMCM2004A题)时最常见的误区是立刻在网上搜“奥运会临时超市网点设计Matlab源码”然后直接改数跑结果。真正让这道题从“解出来”变成“能获奖”的是把题目给出的观众小区分布、场馆坐标、购物概率和超市容量约束转成一组可以被Matlab优化工具箱处理的稀疏矩阵和整数规划模型。这道题表面是商业选址本质上是带容量约束的集合覆盖与需求分配问题。下面按“建模—数据预处理—求解—验证”的顺序给出一套用Matlab复现时可落地的代码路径包含每个参数为什么这么设以及替换数据后需要改哪些变量。2. 先把题面翻译成数学模型变量、目标函数与约束条件2.1 从“网点设计”到集合覆盖定义下标集合常见做法是把问题拆成两类实体需求点和候选网点。题目通常给出若干个观众居住小区或观众来源区每个小区有固定人数并给出每个小区到各个场馆的人流概率临时超市则安排在候选位置比如场馆周边的空地或道路节点。为了建立可求解模型需要先定义以下集合(S)观众需求点集合大小为 (N_s)每个点 (s) 有需求量 (d_s)。(J)候选超市网点集合大小为 (N_j)每个候选点 (j) 有最大服务容量 (c_j)。(K)场馆或热点区域集合用于把需求点与场馆之间的距离纳入目标。题目给的“人流数据”往往不是现成的需求点而是“小区—场馆”的OD量。所以第一步是把OD矩阵按概率折算成每个小时的总购物人数再聚合到需求点坐标上。这个过程不只是读Excel还要做空间重采样。在Matlab里建议用结构体或table保存原始数据而不是散落的工作区变量% 读取题目附带的观众小区数据 raw readtable(demand_points.csv); % 假设有两列坐标x,y和人数pop s.x raw.x; s.y raw.y; s.demand raw.pop; % 若人数是全天总量需要按题目给的购物时间窗口折算 N_s height(raw);这里把需求点直接存成结构体 (s)后面的距离计算和矩阵构建都基于这个结构体。注意需求量不是人口普查人数而是“有购物意愿且在该时段到场馆的人数”题目往往给出比例常见做法是乘以一个购物系数 (\alpha)这个系数在后面会作为敏感性分析的变量。2.2 目标函数覆盖需求最大化与总步行距离最小化CUMCM2004A题常见的评阅标准会让选手权衡两个指标一个是满足的购物需求总量另一个是观众到最近超市的平均距离。在数学上可以写成双目标但为了用Matlab的整数线性规划求解通常把其中一个转成约束或者线性加权。如果以“最小化总步行距离”为主目标可以写成[ \min \sum_{s \in S}\sum_{j \in J} t_{s,j} \cdot y_{s,j} \cdot d_s ]其中 (y_{s,j}) 是0-1变量表示需求点 (s) 是否被分配给候选网点 (j)(t_{s,j}) 是两者步行距离。同时需要引入另一个0-1变量 (x_j) 表示是否在候选点 (j) 建超市。如果以“最大化覆盖需求”为目标则目标函数为[ \max \sum_{s \in S}\sum_{j \in J} y_{s,j} \cdot d_s ]实际建模中我一般取前者因为它对“临时超市”场景更直观临时网点本身是服务设施观众不会跑太远。下面的推导都按最小化总距离展开但切换目标时只要改Matlab的f向量即可。2.3 三类核心约束容量、唯一分配与覆盖半径约束条件决定了模型能不能用线性规划解。第一类是每个候选超市的容量上限[ \sum_{s \in S} d_s \cdot y_{s,j} \le c_j \cdot x_j, \quad \forall j \in J ]这里 (c_j) 是最大服务人数和超市面积、平均购物时间、营业窗口时长相关。如果题目给出的容量单位不是“人/小时”需要先换算。第二类是每个需求点只能被一个已开放网点服务[ \sum_{j \in J} y_{s,j} 1, \quad \forall s \in S ]第三类是覆盖半径约束只有距离小于阈值的需求点才能被分配[ y_{s,j} \le a_{s,j} ]其中 (a_{s,j}) 是0-1可达矩阵由距离阈值 (\tau) 生成。虽然第三类约束可以由距离下限隐含但显式写出能减少变量组合让整数规划求解更快。这三类约束组合起来就是一个典型的带容量约束的选址-分配问题。因为决策变量 (x_j) 和 (y_{s,j}) 都是整数所以需要调用整数规划求解器。下面的表总结了所有符号和对应Matlab变量名符号含义Matlab变量备选值(S)需求点集合s.x,s.y题目给出(J)候选网点集合j.x,j.y网格剖分或题目给(d_s)需求点需求量s.demand人口×购物概率(c_j)候选点容量cap超市面积/人均面积(t_{s,j})步行距离dist欧氏距离×路网系数(x_j)是否建点x(j)0/1(y_{s,j})是否分配y(s,j)0/12.4 用Matlab构建模型矩阵的骨架代码把上面的模型转成intlinprog能接受的矩阵形式最关键的是把双下标变量 (y_{s,j}) 展成一维向量。常见做法是按列展开即idx(s,j) (j-1)*N_s s。下面这段代码构建目标函数和容量约束% 假设已有 dist 矩阵、 cap 向量、 demand 向量 N_s length(s.demand); N_j length(cap); % 决策变量顺序 x(1)..x(N_j) 为建点变量 y(s,j) 展开在后面 n_x N_j; n_y N_s * N_j; f zeros(n_x n_y, 1); % 目标距离 * 需求量只作用于 y 变量 for j 1:N_j for k 1:N_s idx_y n_x (j-1)*N_s k; f(idx_y) dist(k,j) * s.demand(k); end end % 容量约束 sum_s demand(s) * y(s,j) cap(j) * x(j) A zeros(N_j, n_x n_y); b zeros(N_j, 1); for j 1:N_j % 对应 x(j) 的系数是 -cap(j)因为 cap(j)*x(j) 移到左边为 -cap(j) A(j, j) -cap(j); for k 1:N_s idx_y n_x (j-1)*N_s k; A(j, idx_y) s.demand(k); end end这段代码的意义是用矩阵A同时表达容量约束和0-1变量的关系。注意A(j, j) -cap(j)是让不等式变为 (\sum demand \cdot y - cap \cdot x \le 0)。这样写比直接在循环里拼intlinprog的输入更容易检查错误。实际调试时我会先打印A(1,:)前10个非零元素确认idx_y映射没有错位。3. 数据预处理从原始人流数据到可分块的稀疏矩阵3.1 题目数据读入Matlab先做清洗再谈建模CUMCM2004A题通常附带Excel或CSV文件里面往往混合了文本表头、单位行和空行。如果直接readtable会出现NaN或字符串变量混入数值列。推荐的做法是先读成table再用ismissing找到缺测点最后按列类型转换。下面代码处理常见格式T readtable(CUMCM2004A_data.csv, VariableNamingRule, preserve); % 去掉全空的行 T rmmissing(T); % 把小区编号列转成字符串如果原来是文本 if iscell(T{:,1}) id string(T{:,1}); end % 取出坐标和人数列列名按实际文件调整 coords [T{:, x_coord}, T{:, y_coord}]; pop T{:, population}; % 有些题目会给出“购物概率”这里直接乘上去 shopping_prob 0.65; % 这个值来自题目附表或假设 demand pop .* shopping_prob;这里有一个容易踩的坑题目给的人流量可能是“全天总客流”但临时超市只在比赛前后营业3到4小时所以必须把全天量除以营业时段数否则容量约束会严重失配。我一般会先画出demand的直方图如果最大值超过候选点容量的几十倍就要检查是不是忘记换算时间窗口。3.2 坐标转距离矩阵pdist2、矩阵广播和循环的取舍需求点到候选点的距离矩阵是后续建模的核心。如果候选点数量不大比如几百个直接用pdist2最方便% 假设 demand_coords 是 N_s x 2 candidate_coords 是 N_j x 2 dist_matrix pdist2(demand_coords, candidate_coords, euclidean);pdist2计算的是直线距离实际步行距离需要乘以一个绕行系数常见取值在1.2到1.5之间。如果数据量达到几万pdist2可能内存不够这时用分块计算或者直接考虑用路网矩阵。三种方式的对比方法代码复杂度内存占用适用场景pdist2低(O(N_s N_j))数据量在 1e4 以内广播dist sqrt((x_s - x_j.).^2 ...)中同上但更底层需要额外处理维度时双重循环高只需两行数据内存极小的教学环境速度慢实际参赛时不必追求极致的距离精度因为题目给的坐标已经是平面投影坐标线性距离和真实步行距离差别不会影响网点数量级。重点是把dist_matrix转成0-1可达矩阵reachablecover_radius 800; % 单位米来自题目对“步行5分钟”的约束 reachable dist_matrix cover_radius;这里阈值是模型里最重要的一个超参数。后来做敏感性分析时可以用这个阈值作为横轴观察结果变化。3.3 生成候选网点网格剖分与坐标聚类题目不一定直接给候选点列表有时候需要自己生成。最常见做法是把场馆周边区域划分成正方形网格取每个网格中心作为候选点。这样做的好处是网点规模可控且能覆盖整个区域。下面是生成候选点的代码% 由场馆坐标给出区域边界 x_min min(stadium(:,1)) - 500; x_max max(stadium(:,1)) 500; y_min min(stadium(:,2)) - 500; y_max max(stadium(:,2)) 500; grid_step 300; % 每300米一个网格 [gx, gy] meshgrid(x_min:grid_step:x_max, y_min:grid_step:y_max); candidate_coords [gx(:), gy(:)];网格步长直接决定候选点数量和求解难度。如果步长从300改为200候选点数量会从约 (N_s) 级变成 (N_s \times 1.5^2)整数规划耗时通常指数上升。所以建议先跑粗网格快速验证模型正确性再细化网格。还有一种做法是用kmeans对需求点聚类把聚类中心作为候选点。它能保证候选点贴近人流密集处但缺点是会漏掉一些交通枢纽位置。如果题目中明确提到“可以在场馆出入口附近建点”那么网格法中再加入出入口坐标即可。3.4 预处理阶段最容易错的三个细节第一是坐标单位。有的题给经纬度有的给直角坐标如果直接混用距离会差几个数量级。处理方法是用deg2km或utm2deg转换但更简单的是先看题目配图里的比例尺把单位统一成“百米”或“米”。第二是重复点。需求点里可能有多个小区在同一个坐标如果不去重会在距离矩阵里产生大量全同行导致intlinprog出现退化。用unique(coords,rows)合并需求量累加。第三是demand为0的点。如果某小区没有购物需求保留它会白增加变量数量。我通常把demand 1e-6的行直接过滤掉这样在后期求解数据规模更小。4. 求解策略从intlinprog到遗传算法的切换4.1 线性整数规划用Matlab优化工具箱的intlinprog落地当模型只有线性目标和线性约束时优先用intlinprog。它来自Matlab优化工具箱内置分支定界算法对小规模问题变量数几千个能在几十秒内返回全局最优。调用时需要指定决策变量中哪些是整数。在CUMCM2004A题里所有决策变量都是0-1变量所以intcon 1:num_vars。intcon 1:(n_x n_y); lb zeros(n_x n_y, 1); ub ones(n_x n_y, 1); % 唯一分配约束每个需求点只被一个网点服务 Aeq zeros(N_s, n_x n_y); beq ones(N_s, 1); for k 1:N_s for j 1:N_j idx_y n_x (j-1)*N_s k; Aeq(k, idx_y) 1; end end % 调用求解器 options optimoptions(intlinprog, Display, iter, MaxTime, 120); [x_opt, fval, exitflag] intlinprog(f, intcon, A, b, Aeq, beq, lb, ub, options);这段代码里A,b沿用上一章的容量约束Aeq,beq实现唯一分配。运行后会得到x_opt向量前N_j个是建点标志后面是分配关系。注意exitflag不等于1时需要看x_opt中的NaN值通常是没有可行解而不是求解失败。4.2 模型太大时换成遗传算法ga找近似解如果候选点超过500个intlinprog可能会在分支定界时卡住。此时改用全局优化工具箱的ga来求次优解更实际。遗传算法的核心是写一个适应度函数输入建点方案输出总距离并对违反容量约束的方案施加惩罚。function total_cost fitness_func(selected) % selected 是 0/1 向量长度 N_j if sum(selected) 0 total_cost 1e8; return; end % 每个需求点找最近的开放网点 open_idx find(selected); total_cost 0; for s_idx 1:N_s dist_to_open dist_matrix(s_idx, open_idx); [mind, ~] min(dist_to_open); total_cost total_cost demand(s_idx) * mind; end % 容量惩罚开放网点总容量不够时加大惩罚 total_capacity sum(cap(open_idx)); if total_capacity sum(demand) total_cost total_cost 1e6 * (sum(demand) - total_capacity); end end调用ga时需要设置种群大小和代数。适应度函数里的1e6是一个很大的惩罚系数目的是让算法避开不可行解。注意ga默认求解最小化问题所以适应度值越小越好。下面表格列出ga的关键参数及其推荐值参数推荐值作用PopulationSize200种群越大越不容易早熟MaxGenerations300迭代上限看收敛曲线CrossoverFraction0.8交叉比例太高会破坏好解MigrationFraction0.1子种群间迁移率Displayiter打印每代最优值4.3 整数规划与遗传算法的分工建议从工程角度我一般分两步走先用intlinprog在放松网格密度下求一个精确解作为上界参考再把这个解转成initial_population喂给ga让遗传算法在更细的网格上继续搜索。这样比直接让ga从随机解开始快很多。需要注意ga的初始种群必须是矩阵一行一个个体。可以把intlinprog的x_opt(1:N_j)作为第一个个体其余个体用随机01向量。这样遗传算法在早期就能集中在高质量解附近。4.4 求解失败的排查清单遇到intlinprog提示“No feasible solution”先检查三处第一所有需求点的总需求是否大于所有候选点容量之和如果是问题本身无解需要增加候选点或放宽容量第二覆盖半径是否太小导致某些需求点不在任何候选点覆盖范围内打印sum(reachable,2)可以找到那些全0行第三容量约束矩阵A是否写反了不等号方向。如果ga收敛很快但结果明显不合理往往是惩罚系数太大或太小。太大导致算法只关心可行性太小则容易产生容量超限的解。常见的做法是让惩罚系数为容量超限量的1000倍然后看多次运行是否稳定。5. 结果验证与网点调整的实用技巧5.1 用热力图检查覆盖盲区求解完成后把需求点和网点画在一起颜色表示需求量圈表示服务半径可以立刻发现哪些区域没有被覆盖。下面的代码画出覆盖热力图figure; scatter(demand_coords(:,1), demand_coords(:,2), 30, demand, filled); hold on; open_j find(x_opt(1:N_j) 0.5); viscircles(candidate_coords(open_j,:), cover_radius * ones(size(open_j)), ... Color, red, LineWidth, 1.2); colorbar; xlabel(x 坐标); ylabel(y 坐标); title(需求点与开放超市覆盖范围);如果发现某个需求点被多个圆圈覆盖说明容量分配可能不均匀。可以通过查看Aeq对应的y变量来检查每个开放网点的实际负荷并把手动调整候选点坐标作为下一步的输入。5.2 用敏感性分析选定最优覆盖半径CUMCM2004A题评分时很看重对参数的讨论。最直接的做法是把覆盖半径从600米逐步改到1200米步长100米每次调用同一个求解脚本记录总步行距离和开放网点数。把这个循环写成脚本最后用plot画折线图。如果曲线在某一半径后变得平缓那这个点就是合理的阈值。5.3 把求解结果整理成可直接放进论文的表格最后一步是把x_opt转成表格输出选中网点的坐标和服务人数占比。可以用writetable生成Excel方便在论文里直接截图selected_j find(x_opt(1:N_j) 0.5); out_table table(candidate_coords(selected_j,1), candidate_coords(selected_j,2), ... cap(selected_j), VariableNames, {x, y, capacity}); out_table.served zeros(height(out_table), 1); for jj 1:length(selected_j) out_table.served(jj) sum(demand(find_allocated_to_j(selected_j(jj)))); end writetable(out_table, resolved_networks.xlsx);这里find_allocated_to_j是从y变量中还原每个网点服务了哪些需求点的小函数。整个流程跑通后换一套数据只需要改第2章的readtable文件名和第3章的cover_radius。本文还有配套的精品资源点击获取