1. 项目背景与核心问题拆解2020年数维杯数学建模C题题目全称是“垃圾转运优化模型设计”这是一个典型的运筹学与物流优化问题。当时拿到这个题目很多队伍的第一反应是这不就是个车辆路径问题VRP吗套个模型调个算法跑个结果就完事了。但真正深入进去你会发现从“知道这是个VRP”到“建出一个能求解、能解释、能应对各种现实约束的模型”中间隔着十万八千里。这道题之所以经典就是因为它把一个看似标准的VRP问题嵌套在了城市垃圾收运这个非常具体的、充满现实掣肘的场景里。它考察的绝不仅仅是你会不会用粒子群算法或者遗传算法而是你如何将一堆模糊的、矛盾的现实条件转化为清晰的、可计算的数学模型并且最终给出一个在数学上严谨、在现实中可行的调度方案。简单来说题目给了一个虚构的城市区域里面有若干个垃圾产生点居民区、商业区等和若干个垃圾转运站。每天垃圾收运车从停车场出发去各个产生点收集垃圾然后运往转运站进行处理最后空车返回停车场。这里面的核心矛盾是垃圾产生点的垃圾量每天不同转运站有处理能力上限收运车有载重和行驶距离限制而我们的目标是以最低的成本通常包括固定车辆成本、行驶距离成本等完成每天的收运任务。听上去规则很清晰对吧但魔鬼藏在细节里车辆到了转运站卸货时间算不算如果某个产生点的垃圾量超过一辆车的容量是必须一次运完还是可以分多次车辆在转运站卸货后是必须返回停车场还是可以继续去下一个点收集这些看似“理所当然”或者“可以假设”的细节恰恰是模型能否贴近实际、求解是否可行的关键。我当时的思路是不能一上来就埋头建模型而是要先当一回“垃圾收运调度员”在脑子里把整个流程跑一遍把每一个可能卡住我们的现实环节都揪出来。2. 模型构建从现实约束到数学方程建模的第一步是定义清晰。我们把所有垃圾产生点称为“需求点”每个点有确定的垃圾量把转运站称为“设施点”每个点有确定的最大处理能力把收运车称为“车辆”每辆车有最大载重和最大行驶距离。停车场是车辆的起点和终点。目标函数很明确最小化总成本一般包括动用车辆的固定成本和所有车辆行驶路径的总距离成本假设距离成本与油耗、时间成正比。接下来就是重头戏约束条件的数学化。这是最考验建模者功力的地方也是论文拉开差距的地方。### 2.1 核心约束的数学表达流量平衡约束这是路径问题的基石。对于任何一个非停车场、非转运站的需求点进去的车辆数等于离开的车辆数并且只能有一辆车服务一次假设一辆车一次装完该点垃圾。用数学语言说就是对于每个需求点i其相关的决策变量表示车辆k是否从点i行驶到点j的流入和流出相等且总和为1。这个约束保证了路径的连续性不会出现车“瞬移”或者“断头路”。载重约束任何时刻车辆k上的垃圾装载量不能超过其最大载重。这需要引入一个额外的变量来跟踪车辆在访问每个点后的累计装载量。当车辆访问一个需求点时装载量增加当访问转运站时装载量清零。这个约束是模型非线性的一个重要来源因为装载量与路径选择是耦合的。转运站处理能力约束所有车辆运到转运站s的垃圾总量不能超过该转运站日处理能力。这需要统计所有以转运站s为终点的路径段即车辆前来卸货所携带的垃圾量。这个约束将各个独立的车辆路径耦合在了一起因为它们共享有限的转运站资源。车辆使用与路径唯一性约束一辆车如果被使用它必须从停车场出发最后返回停车场。同时要防止出现“子回路”——即车辆没有服务所有指派给它的点而是自己绕了一个小圈圈就回去了。这是VRP建模的经典难题通常通过引入“MTZ约束”或“流约束”来消除。### 2.2 关键细节的模型化处理题目没有明确说明的细节需要我们自己做出合理且自洽的假设并在模型中体现卸货时间如果考虑卸货时间相当于增加了车辆在转运站的“服务时间”。这会影响车辆的最大行驶时间约束如果有时限要求或者可以将其折算为额外的固定成本。在我们的模型中为了简化可以将转运站视为一个特殊的“需求点”其服务时间卸货时间为固定值。大垃圾点拆分如果一个需求点的垃圾量超过单车容量必须允许被多辆车服务或一辆车多次访问。这会将问题从标准的VRP每个点访问一次变为更复杂的“可拆分需求的VRP”。实现方式有两种一是在预处理阶段将该大点虚拟拆分为多个位置相同、垃圾量之和为原总量的子点二是在模型层面允许决策变量表示“车辆k从点i运走多少量的垃圾”这大大增加了模型复杂度。车辆是否必须返回停车场再出发这决定了模型的网络结构。如果必须返回那么每次从转运站出来都相当于开始一次新的行程停车场-需求点-转运站-停车场构成一个闭合回路。如果允许直接前往下一个需求点那么路径就是停车场-需求点-转运站-需求点-...-转运站-停车场。后者通常能得到成本更低的解但模型更复杂。我们选择后者因为它更符合实际运营中提高车辆利用率的逻辑。把这些约束和假设用数学方程主要是0-1整数规划或混合整数规划写下来一个完整的垃圾转运优化模型就初具雏形了。但模型建得好不代表能解得出来。这类问题属于NP-hard问题随着问题规模点、车数量增大精确算法如分支定界法会在短时间内失去可行性。3. 算法选择与求解策略为什么是启发式算法面对一个大规模的混合整数规划模型直接丢给MATLAB的intlinprog或者商业求解器Gurobi、CPLEX很可能几个小时都得不到一个可行解。这时候就必须借助启发式算法。2020年参赛时粒子群算法PSO是很多队伍的选择因为它概念相对直观实现起来不算太复杂而且对于连续和离散优化问题都有不错的适配性。但这里有一个关键的思维转换如何用粒子群算法来求解一个路径规划问题粒子群算法中的“粒子”其位置通常代表优化问题的一个解。在连续优化中位置是实数向量而在VRP中一个解是一组车辆的路径序列。因此编码与解码是首要难题。### 3.1 粒子编码设计一种有效的路径表示法我们采用了一种称为“基于客户点排列的编码”方式。具体来说首先生成一个包含所有需求点的随机排列序列。例如有10个需求点生成序列 [3, 7, 1, 9, 4, 10, 2, 8, 5, 6]。然后按照这个序列的顺序结合车辆的载重约束和转运站能力约束对这个序列进行“分割”分配车辆和确定转运站访问点。从第一辆车开始按序列顺序将需求点加入当前车辆的路径。每加入一个点检查车辆当前累计载重是否超限。如果加入当前点后超限则该点不能加入当前车辆路径结束。车辆需要前往一个转运站卸货选择哪个转运站可以基于距离最近或处理能力余量等规则然后从转运站开始用同一辆车继续从序列中装载下一个点如果允许车辆连续作业或者换下一辆车从停车场出发继续装载。同时需要检查所选转运站的处理能力余量是否足够接收这批垃圾如果不够则需要选择其他转运站。这个过程一直持续到所有需求点都被服务完毕。最终这个序列以及分割规则就唯一确定了一组车辆路径方案解。在粒子群算法中每个粒子的“位置”就是这个需求点的排列序列。但粒子群更新公式速度位置是针对实数向量的如何对排列序列进行“加减”和“移动”呢这就需要定义特殊的算子。### 3.2 离散粒子群更新交换子与速度的含义我们将粒子的“速度”定义为一系列“交换操作”的概率或强度。例如速度可以是一个实数向量其长度与序列长度相同值的大小表示对应位置元素需要发生交换的意愿强度。 更新时根据个体历史最优解和全局历史最优解计算出每个位置的建议交换概率。然后通过一个随机过程比如将序列中某些位置的点以一定概率进行交换如两点交换、片段逆序等来模拟“位置”的更新。这样粒子就能在解空间所有可能的排列中进行搜索。### 3.3 目标函数与约束处理解码得到路径方案后就可以计算总成本固定成本行驶距离。但这里还有一个关键解码过程中生成的方案一定能满足载重和转运站能力约束吗我们上面描述的分割过程是“贪婪”的它主动保证了载重约束并通过规则选择转运站来尝试满足处理能力约束。但如果转运站能力非常紧张这种贪婪构造法可能失败产生不可行解。处理约束的常用方法有罚函数法将约束违反程度乘以一个很大的惩罚系数加到目标函数上。这样算法在优化时会倾向于远离不可行区域。但惩罚系数设置需要技巧太大则搜索僵化太小则约束形同虚设。可行解保持法设计特殊的解码规则和更新算子确保产生的解始终是可行的。这种方法更优雅但设计难度更大。在实际编程中我们采用了罚函数法与修复策略相结合的方式。首先在解码时尽量遵守约束。如果仍然违反了转运站能力约束比如某个站接收的垃圾超了则在目标函数中加上一个巨大的惩罚项。同时我们会记录这个不可行解但它的适应度会很差在粒子群的选择压力下会被淘汰。为了帮助算法跳出局部最优我们允许迭代过程中存在少量不可行解带着惩罚它们可能包含着通往更优可行解的结构信息。4. MATLAB实现核心步骤与踩坑实录理论说得再漂亮代码跑不起来都是白搭。用MATLAB实现上述算法有几个核心模块和一堆坑等着你。### 4.1 数据准备与距离矩阵计算第一步永远是数据。题目会给出各点的坐标或地址。我们需要计算所有点两两之间的行驶距离。这里第一个坑就来了是计算直线欧氏距离还是考虑实际道路的曼哈顿距离街区距离对于城市垃圾收运道路通常是网格状的曼哈顿距离更真实。我们假设了道路网格采用曼哈顿距离公式d |x1-x2| |y1-y2|。将这个距离矩阵预先算好并存储后续所有路径距离计算都通过查表完成这是巨大的性能优化点。% 假设 points 是一个 Nx2 的矩阵第一列是x坐标第二列是y坐标 num_points size(points, 1); dist_matrix zeros(num_points, num_points); for i 1:num_points for j 1:num_points if i ~ j dist_matrix(i, j) abs(points(i,1)-points(j,1)) abs(points(i,2)-points(j,2)); end end end### 4.2 粒子群算法主框架搭建主循环结构是标准的PSO框架但关键在updatePosition和calculateFitness这两个自定义函数。% 参数初始化 pop_size 50; % 粒子数量 max_iter 200; % 最大迭代次数 w 0.729; % 惯性权重 c1 1.49445; % 个体学习因子 c2 1.49445; % 社会学习因子 % 初始化粒子位置随机排列和速度 particle_position cell(pop_size, 1); % 每个粒子位置是一个序列 particle_velocity rand(pop_size, num_demand_nodes); % 速度矩阵用于决定交换概率 % ... 初始化每个particle_position为随机排列 ... pbest_position particle_position; % 个体历史最优位置 pbest_value inf(pop_size, 1); % 个体历史最优值 gbest_position []; % 全局历史最优位置 gbest_value inf; % 全局历史最优值 % 主循环 for iter 1:max_iter for i 1:pop_size % 1. 更新粒子速度针对离散问题这里是更新交换概率向量 r1 rand(1, num_demand_nodes); r2 rand(1, num_demand_nodes); % 这里需要自定义一个函数将pbest和gbest与当前位置的差异转化为“速度”增量 velocity_increment c1*r1.*(position_diff_to_probability(particle_position{i}, pbest_position{i})) ... c2*r2.*(position_diff_to_probability(particle_position{i}, gbest_position)); particle_velocity(i, :) w * particle_velocity(i, :) velocity_increment; % 2. 根据速度更新粒子位置序列 particle_position{i} update_sequence_by_velocity(particle_position{i}, particle_velocity(i, :)); % 3. 解码位置计算适应度总成本 [total_cost, routes, constraint_violation] decode_sequence(particle_position{i}, dist_matrix, vehicle_capacity, station_capacity); fitness total_cost penalty_factor * constraint_violation; % 4. 更新个体最优和全局最优 if fitness pbest_value(i) pbest_value(i) fitness; pbest_position{i} particle_position{i}; end if fitness gbest_value gbest_value fitness; gbest_position particle_position{i}; end end % 可以在这里动态调整惯性权重w实现前期探索后期收敛 end### 4.3 解码函数decode_sequence算法的心脏这是整个程序最复杂、最容易出bug的部分。它输入一个序列输出总成本、详细的路径列表以及约束违反程度。function [total_cost, routes, violation] decode_sequence(sequence, dist_matrix, cap_vehicle, cap_station) num_demands length(sequence); demands ...; % 各需求点垃圾量数组 station_ids ...; % 转运站编号列表 depot_id ...; % 停车场编号 total_cost 0; routes {}; % 用于存储每辆车的路径如 {[depot, 3, 7, station1, depot], ...} violation 0; current_route [depot_id]; current_load 0; remaining_station_cap cap_station; % 各转运站剩余能力数组 idx 1; % 指向序列中待处理的需求点 vehicle_count 0; while idx num_demands next_node sequence(idx); next_demand demands(next_node); % 情况1当前车辆能装下这个点 if current_load next_demand cap_vehicle % 尝试将点加入当前路径 current_route [current_route, next_node]; current_load current_load next_demand; idx idx 1; % 处理下一个点 else % 情况2当前车辆装不下了需要先去卸货 % 选择一个转运站例如选择当前路径末端距离最近的、且有容量的站 [selected_station, dist_to_station] select_station(current_route(end), station_ids, remaining_station_cap, dist_matrix); if isempty(selected_station) % 没有找到可用的转运站严重违反约束 violation violation 1e9; % 施加巨大惩罚 % 应急处理强行结束当前路径返回停车场这是一个不可行解的部分 total_cost total_cost dist_matrix(current_route(end), depot_id); current_route [current_route, depot_id]; routes{end1} current_route; vehicle_count vehicle_count 1; total_cost total_cost vehicle_fixed_cost; % 加上车辆固定成本 % 开始新的路径 current_route [depot_id]; current_load 0; % 注意当前点(next_node)还没有被服务idx没有增加下一轮循环会重新尝试 else % 找到可用转运站前往卸货 current_route [current_route, selected_station]; total_cost total_cost dist_to_station; remaining_station_cap(selected_station) remaining_station_cap(selected_station) - current_load; current_load 0; % 卸货后不一定要回停车场可以直接继续服务下一个点 % 这里我们选择继续服务所以路径继续延伸 end end % 额外判断如果当前路径已经很长或者所有点都服务完了则结束当前车辆路径返回停车场 if idx num_demands || some_other_condition if current_route(end) ~ depot_id total_cost total_cost dist_matrix(current_route(end), depot_id); current_route [current_route, depot_id]; end if length(current_route) 2 % 不只是[depot, depot] routes{end1} current_route; vehicle_count vehicle_count 1; total_cost total_cost vehicle_fixed_cost; end break; end end total_cost total_cost sum(dist_matrix_for_routes(routes)); % 加上所有路径内部的行驶距离 end### 4.4 实际踩坑与调试经验初始解质量至关重要完全随机的初始序列解码出来的路径可能极其糟糕甚至大部分都是不可行解违反转运站能力。这会导致算法初期在极差的解空间里徘徊。我们的改进方法是采用一些简单的启发式规则生成初始序列比如按照需求点距离停车场的远近排序或者按照“垃圾量/距离”的密度排序。这能提供一个不错的起点。罚因子的动态调整固定的大罚因子如1e9虽然能保证最终解可行但可能会像一堵墙把搜索限制在可行域边界很小的区域内。我们尝试了动态罚因子初期设置较小的罚因子允许算法探索一些不可行区域随着迭代进行逐渐增大罚因子迫使粒子收敛到可行域。这有点像模拟退火的思想效果提升明显。局部搜索的引入单纯的粒子群算法在路径优化上容易早熟。我们在每次迭代后对全局最优解gbest_position对应的序列进行局部搜索。例如随机交换序列中的两个点或者逆序一段子序列看看能否得到更好的解。这个步骤虽然增加了计算量但能显著提升解的质量。MATLAB性能瓶颈解码函数decode_sequence会被调用成千上万次粒子数×迭代次数。一定要做好预分配、向量化操作避免在循环中动态增长数组。dist_matrix查表比实时计算距离快得多。对于大规模问题可以考虑用更快的语言如C编写核心解码函数通过MEX接口在MATLAB中调用。结果的可视化与验证一定要把求出的最优路径画出来用plot函数把点、转运站、停车场标出来然后用线条把每辆车的路径画出来。肉眼一看就能发现很多问题路径交叉严重通常不是最优、某辆车绕了远路、某个转运站过于拥挤。可视化是调试和验证模型逻辑最直观的工具。5. 模型检验、灵敏度分析与论文写作要点模型和算法跑通了得到了一个“最优”解工作只完成了一半。如何让人相信你的解是合理的、稳健的这就需要模型检验和灵敏度分析。### 5.1 模型检验与简单规则对比一个强有力的说服方法是将你的优化模型得到的结果与一两种简单的、符合直觉的调度规则进行对比。例如最近邻规则车辆总是从当前位置前往最近的未服务需求点装不下就去最近的转运站。先到先得规则车辆按照需求点编号顺序服务。用同样的数据和成本计算方法运行这些规则算出总成本。通常你的优化模型结果应该显著优于这些简单规则成本降低10%-30%或更多。这个对比能直观地体现你模型的价值。### 5.2 灵敏度分析关键参数的影响数学建模论文非常看重灵敏度分析即研究当某些关键参数发生变化时你的模型结果如总成本、所需车辆数如何变化。这展示了模型的鲁棒性和你对问题深度的理解。主要分析点包括垃圾产生量的波动将所有需求点的垃圾量统一增加或减少10%重新求解观察总成本和路径方案的变化。结果可能会发现成本增加并非线性的可能存在规模效应。转运站处理能力逐个收紧某个转运站的处理能力观察系统是否变得脆弱成本是否急剧上升。这可以帮助识别系统中的瓶颈设施。车辆载重改变车辆的载重上限分析其对车辆使用数量和总成本的影响。通常会有一个经济载重区间。单位距离成本与车辆固定成本的比值这个比值会影响模型的决策倾向。比值高时模型会更倾向于减少总行驶距离哪怕多用几辆车比值低时模型会更倾向于减少用车数量哪怕有些车跑得远一点。分析这个比值的影响能为车队配置是大车少辆还是小车多辆提供决策依据。做灵敏度分析时不需要对每个参数变化都重新运行完整的优化算法太耗时。可以基于得到的最优解进行“如果-那么”的局部调整和估算或者对关键参数进行小范围的重新优化。### 5.3 论文写作的核心讲好一个解决问题的故事数模论文的本质是技术报告但好的论文是在讲故事。你的故事线应该是问题重述与分析不是简单抄题目而是用你的理解提炼出核心矛盾、决策变量和目标。模型假设清晰列出所有假设并说明其合理性。这是你简化现实世界的依据。模型建立这是核心章节。建议分小节1符号说明用表格2目标函数3约束条件每条约束用文字描述其实际意义再给出数学公式4模型总结指出这是一个什么类型的数学模型。算法设计详细说明为什么选择PSO或其他算法如何将你的模型与算法适配编码、解码、约束处理并给出清晰的算法流程图或伪代码。求解与结果分析给出具体的数据、参数设置、最终求出的路径方案可以用表格列出每辆车的行驶顺序、总成本。一定要有结果的可视化图如路径网络图。模型检验与灵敏度分析如上所述展示模型的优越性和鲁棒性。优缺点与推广客观评价自己模型的优点如贴近实际、求解高效和缺点如未考虑交通拥堵、假设较理想并提出可能的改进方向或模型在其他领域的应用如快递配送、校车路线规划。写作时公式要清晰编号图表要规范美观MATLAB出图后可以用Visio或PPT稍作美化文字表述要专业且流畅。避免口语化但也要避免过于晦涩。最终你的论文应该让一个没参与建模的读者也能清晰地理解你面对的问题、你的解决思路、以及你的方案为什么好。回过头看这道题之所以让人印象深刻就是因为它完美地诠释了数学建模的全过程从现实抽象出模型用算法求解模型再回到现实去解释和检验结果。每一个环节都有无数细节可以打磨也都有无数的坑等着你去踩。而真正让你成长的恰恰是处理这些细节和爬出这些坑的过程。