Matlab遗传算法求解共享单车调度问题:从建模到完整源代码 简介自行车调度问题在数学建模竞赛中属于常见的优化应用题核心是结合站点分布、用户需求预测、调度策略与运输成本在供需平衡和时间窗限制下找到成本最小的自行车调配方案。这套MATLAB源代码压缩包正是一套可直接运行的完整程序集适合数学建模参赛者、运筹优化方向的学生以及需要快速搭建调度模型的开发者使用。包内共46个文件主要包含24个m脚本负责数据预处理、模型定义、算法求解、7个mat数据文件用于存取站点坐标与中间结果、6张jpg图像展示调度路径和迭代收敛曲线另有7个asv自动备份和txt说明文档整体打包大小仅为237KB结构清晰便于按模块学习与二次开发。代码实现覆盖了从数据导入、初始解生成到遗传算法中选择、交叉、变异以及免疫优化等启发式算子的完整流程并带有时间窗检查与结果可视化等辅助函数。已有852人查看学习对于希望掌握MATLAB建模实践、算法编写或进行课程设计的人来说具有较好的参考价值和可复用性。1. 自行车调度问题为什么是数学建模竞赛的常客自行车调度问题本质上是让一批调度车辆在有限时间内把共享单车从“淤积站”搬到“空桩站”同时让总成本最低。这类题在数学建模竞赛里出镜率极高因为它同时包含资源分配、路径规划和多目标权衡非常适合考察参赛者从问题定义到算法落地的完整能力。标题里的“matlab程序源代码”意味着读者不只需要一个漂亮的目标函数更需要一套能直接改、能跑出结果的完整代码框架。现实中做这件事比看上去难得多早高峰地铁口半小时内涌入几百辆车调度车还没到站点已经溢出晚高峰反过来小区门口一排空锁桩。调度车辆数量有限、载车量有限、站点距离和道路拥堵又给每一条路径叠加了时间成本。问题一旦超过五六个站点手算和穷举就已经失效必须用数学模型加启发式算法求解。这篇文章按“建模 → 写代码 → 调参数 → 排坑 → 进阶”的顺序展开目标是让你拿到一套能在 Matlab 里跑通的最小可用方案。新手照着抄能出图出数熟手可以跳过基础部分直接看参数边界和进阶方向。下面列出的代码都不是伪代码是可以在 R2018a 之后的版本里直接运行的完整片段你只需要准备全局优化工具箱。2. 建模先行把站点盈缺写成一个可计算的目标函数2.1 问题描述与符号定义站点盈缺、调度车与时间窗在做任何代码之前先把问题“数学化”。一个典型的自行车调度场景是这样的某区域里有 N 个站点每个站点在某个时间段的初始车辆数为 S_i预测需求为 D_i站点的桩位容量为 C_i。调度任务就是用 M 辆调度车从车场出发按一定路径访问站点补充或运走自行车使得最终各站点的车辆数趋向于一个理想区间 [L_i, U_i]同时总花费最小。先定义符号体系后面写代码和改代码都依赖这套符号N站点总数编号 1 到 N车场编号为 0M可用的调度车数量Q单辆调度车的最大装载量S_i站点 i 的初始自行车数量C_i站点 i 的桩位容量d_ij站点 i 到站点 j 的距离单位为百米或公里需要统一t_ij站点 i 到站点 j 的行驶时间st_i在站点 i 进行一次装/卸操作所需的服务时间tw_i [a_i, b_i]站点 i 的可访问时间窗车必须在这个窗口内到达才能作业x_ijk调度车 k 是否从站点 i 行驶到站点 j取 0 或 1y_ik调度车 k 是否访问站点 i取 0 或 1p_i站点 i 的调度需求正值表示需要运走负值表示需要补充。有了这套符号问题的目标就很清楚了设计 M 条从车场出发再回到车场的路径在满足时间窗和容量约束的前提下最小化总成本。成本不但包含行驶距离和车辆使用数还要包含对“调度不足”的惩罚否则算法一定会选择一辆车都不派成本为零但那不是我们想要的解。2.2 三种建模思路对比整数规划、遗传算法与仿真优化同一个调度问题可以选择不同的数学框架。第一种是混合整数线性规划MILP把路径决策、装载量决策都写成线性约束直接用 intlinprog 求解。对于站点数小于 15、车辆数小于 5 的小规模问题MILP 能给出全局最优解适合作为对照基准。但站点一旦超过 30 个MILP 的求解时间会指数级上升竞赛时间不允许你等。第二种是遗传算法GA把每条路径编码成一个“个体”通过选择、交叉、变异迭代搜索。GA 不能保证全局最优但能在合理时间内找到工程上可用的近似解而且实现灵活容易加入时间窗和容量约束。第三种是仿真优化用离散事件模拟评估每个调度策略在动态需求下的表现适合研究需求波动很大的场景但计算量最大一般放到进阶阶段。我的建议很直接竞赛题目如果没有明确要求精确解就用遗传算法。原因有三点——Matlab 自带的 ga 函数封装程度高不需要从零实现选择算子罚函数法可以非常自然地处理复杂的业务约束求解规模可以撑到上百个站点足以覆盖绝大多数赛题。后面所有代码都基于遗传算法展开。2.3 目标函数的构成行驶成本、固定成本与惩罚项目标函数是建模的核心我把它拆成三个部分相加。第一部分是行驶成本即所有调度车的总行驶距离乘以单位距离成本。第二部分是固定成本即每启用一辆调度车就产生一个固定费用这能促使算法尽量少用车避免出现“一个站点派一辆车”的奢侈方案。第三部分是惩罚项对每个站点最终车辆数偏离理想区间的量进行惩罚。用数学式子表达就是min Z α · ΣΣΣ x_ijk · d_ij β · Σ sign(y_ik) λ · Σ penalty_i其中 α 是单位行驶成本β 是单车的启用车成本λ 是惩罚系数penalty_i 需要根据最终车辆数落在哪个区间来分段计算。如果调度后站点 i 的车辆数在 [L_i, U_i] 内惩罚为 0如果低于 L_i按缺车辆数乘一个短缺惩罚系数如果高于 U_i按溢出量乘一个溢出惩罚系数。两个方向的惩罚系数通常不一样因为“空桩导致用户无法还车”和“车辆堆积导致用户无法取车”在业务上严重程度不同这点在后续代码里会体现。设计适应度函数时直接把目标函数取倒数或取负值。注意一个常见误区不要把惩罚系数设得过大否则遗传算法会把“少派一辆车”的优化权重压到几乎为零所有个体都在满足约束的窄带里搜索过早收敛。后面避坑章节会专门讨论这个问题。3. Matlab 源码实现从数据生成到 ga 求解的最小完整方案3.1 生成模拟站点数据用随机种子固定一套可复现的测试集写算法之前需要一套可靠的测试数据。这里我用固定随机种子生成站点坐标、容量和需求这样每次运行代码得到的结果可复现方便你检查算法改动带来的影响而不是把随机波动误当成了优化效果。rng(42); % 固定随机种子保证结果可复现 numStations 20; % 站点数量 stationPos rand(numStations, 2) * 10; % 站点坐标范围 0~10 公里 capacity randi([20, 60], numStations, 1); % 每个站点桩位容量 20~60 initialBikes randi([5, capacity], numStations, 1); % 初始车辆数不超过容量 demandBikes randi([-10, 10], numStations, 1); % 正为多余负为短缺 % 计算站点间距离矩阵车场坐标为 (0,0) depotPos [0, 0]; allPos [depotPos; stationPos]; distMat sqrt((allPos(:,1) - allPos(:,1)).^2 (allPos(:,2) - allPos(:,2)).^2);这段代码的目的是构造一个“站点网络”坐标矩阵生成 20 个随机站点容量在 20 到 60 之间波动初始车辆数在容量范围内随机需求同时存在正负值——正值代表该站车辆过多需要运走负值代表需要补充。距离矩阵用欧氏距离计算如果你有真实地图数据可以替换成路网距离矩阵代码其余部分不需要改动。参数说明站点数量、容量范围、需求幅度都影响求解难度建议先用这套参数跑通再逐步加大站点数到 50、100 测试算法极限。3.2 编码方式与约束处理报童思路下的路径串与罚函数设计遗传算法的核心是“怎么用一串数字表示一个调度方案”。这里采用最常用的整数编码每个个体是一个向量其长度等于调度车数量乘站点数每一段表示一辆车的访问站点序列。站点编号从 1 到 20编号 0 表示车场序列中连续的非零数字代表车辆按顺序访问的站点。例如个体片段 [3 5 0 2 7] 表示一辆车先去站点 3再去站点 5回到车场另一辆车去站点 2 和 7。这种编码方式有一个天然问题遗传算子的交叉和变异可能把序列改得支离破碎比如同一个站点被访问两次或某些站点从未被访问。处理方法是罚函数而不是硬性修复。每出现一次重复访问就在适应度中扣掉一个较大的惩罚值引导算法自行淘汰这类非法个体。function score fitnessFunc(chromosome, model) % 输入染色体路径编码、模型参数结构体 % 输出得分需最小化 numVeh model.numVeh; numSta model.numSta; totalDist 0; vehUsed 0; penalty 0; stationLoad model.initialBikes; % 各站点当前的车辆数 % 按车辆数切分染色体 segLen length(chromosome) / numVeh; visitedCount zeros(numSta, 1); for k 1:numVeh seg chromosome((k-1)*segLen1 : k*segLen); seg seg(seg 0); % 去掉 00 表示空闲 if isempty(seg) continue; end vehUsed vehUsed 1; curPos 1; % 从车场出发车场索引为 1 load 20; % 调度车初始装载 20 辆车 for s seg if s numSta continue; end visitedCount(s) visitedCount(s) 1; % 累加行驶距离 totalDist totalDist distMat(curPos, s1); % 模拟车辆装卸需求为正就把车装走为负就放下车 bikeDiff model.demandBikes(s); if bikeDiff 0 takeAway min(bikeDiff, load); stationLoad(s) stationLoad(s) - takeAway; load load - takeAway; else putBack min(-bikeDiff, model.capacity(s) - stationLoad(s)); stationLoad(s) stationLoad(s) putBack; load load putBack; end curPos s 1; end % 返回车场的距离 totalDist totalDist distMat(curPos, 1); end % 检查站点访问次数 for i 1:numSta if visitedCount(i) 1 penalty penalty 1000 * (visitedCount(i) - 1); end end % 检查站点车辆数是否在理想区间 [0.2*C, 0.8*C] for i 1:numSta idealLow 0.2 * model.capacity(i); idealHigh 0.8 * model.capacity(i); if stationLoad(i) idealLow penalty penalty 500 * (idealLow - stationLoad(i)); elseif stationLoad(i) idealHigh penalty penalty 300 * (stationLoad(i) - idealHigh); end end score totalDist 200 * vehUsed penalty; end这个函数是整套代码最关键的模块每一行都在定义一个“好方案”的偏好。逻辑说明先按车辆数把染色体切段每段代表一条路径路径中的站点逐个访问根据需求正负模拟装车或卸车同时维护站点实时车辆数和调度车实时载重最后计算三类惩罚——重复访问站点惩罚、站点车辆数低于理想下限惩罚、高于理想上限惩罚。参数说明重复访问惩罚 1000 是最大级别的约束违背代价应当显著高于正常行驶成本短缺惩罚 500 和溢出惩罚 300 的差异体现了业务上“缺车比堆车更严重”的优先级你可以按赛题要求调整这个比例。3.3 主程序调用 ga 函数的最小可跑通代码主程序负责把数据组装成模型结构体、设置遗传算法选项、调用 ga 并输出结果。这里没有引入复杂的自定义算子ga 自带的交叉和变异就能处理只要编码格式合理。% 组装模型结构体 model.numSta numStations; model.numVeh 4; % 调度车数量 model.distMat distMat; model.capacity capacity; model.initialBikes initialBikes; model.demandBikes demandBikes; model.vehCapacity 20; % 单辆调度车最大载车数 % 染色体长度 车辆数 * 8每辆车最多访问 8 个站点 segMax 8; chromLen model.numVeh * segMax; % 设置 ga 选项 options optimoptions(ga, ... PopulationSize, 200, ... MaxGenerations, 500, ... Display, iter, ... PlotFcn, gaplotbestf, ... UseParallel, false); % 定义变量边界允许 0~numSta 的整数0 代表该位置空闲 lb zeros(1, chromLen); ub numStations * ones(1, chromLen); IntCon 1:chromLen; % 所有变量都是整数 % 调用 ga 求解 [x, fval] ga((x) fitnessFunc(x, model), chromLen, [], [], [], [], ... lb, ub, [], IntCon, options); % 输出结果 fprintf(最优成本: %.2f\n, fval); disp(最优染色体:); disp(x);逻辑说明染色体长度设为车辆数乘 8即每辆车最多访问 8 个站点超出部分用 0 填充表示空闲。这种定长编码的好处是 ga 的初始化、交叉、变异都不需要特殊处理坏处是当车辆路径过长时信息被截断需要根据站点规模调整 segMax。参数说明种群 200、代数 500 是一组折中值适合 20 个站点的规模站点数加大到 50 时建议种群 500、代数 800UseParallel 设为 true 可以启用并行计算加速但需要额外配置并行池。IntCon 声明所有变量为整数这是调度问题正确性的关键漏掉它 ga 会把站点编号变成小数。3.4 结果重构与可行性检查从染色体还原调度方案ga 输出的只是数字序列必须还原成可读的调度路径才能判断结果是否合理。这一步放在主程序之后建议独立成一个函数便于反复调用。function printSchedule(x, model) numVeh model.numVeh; segLen length(x) / numVeh; for k 1:numVeh seg x((k-1)*segLen1 : k*segLen); seg seg(seg 0); if isempty(seg) fprintf(车辆 %d: 未使用\n, k); else routeStr 0; for s seg routeStr [routeStr, - , num2str(s)]; end routeStr [routeStr, - 0]; fprintf(车辆 %d: %s\n, k, routeStr); end end % 检查是否存在重复访问 allVisited x(x 0); if length(allVisited) ~ length(unique(allVisited)) fprintf(警告: 存在重复访问的站点结果不合法!\n); else fprintf(结果合法性检查通过: 所有站点最多被访问一次\n); end end逻辑说明输出的格式是“0 - 3 - 5 - 0”直观展示每辆车的访问序列。合法性检查用 unique 判断是否有重复站点这是最基础的合理性验证。参数说明如果结果频繁出现重复访问说明罚函数强度不够需要把重复惩罚从 1000 提高到 3000而不是去调整遗传算子。这一步看似简单但在实际项目中很多参赛队伍拿到解后没有还原就直接写进论文结果路径出现环路评审扣分惨重。4. 参数调优的三个硬指标种群规模、进化代数与罚函数系数4.1 种群规模与进化代数先看收敛曲线再定预算遗传算法的参数没有统一标准但有一个可靠的调试顺序先固定一个中等规模比如 200 种群、300 代运行后查看 gaplotbestf 给出的收敛曲线。如果曲线在最后 50 代还在明显下降说明代数不够加大到 800如果曲线在 100 代左右就平了说明种群多样性不足应加大种群而不是继续加代数。一个经验值站点数为 N 时种群规模取 10N 到 15N进化代数取 30N 到 50N。比如 20 个站点种群 200 到 300代数 600 到 1000。这个范围比大多数竞赛队伍的默认设置要激进但换来的是更稳定的解。如果运行时间过长优先降代数而不是降种群因为种群深度决定搜索广度代数只决定搜索终点。4.2 交叉与变异概率保持多样性和收敛性的平衡ga 默认的交叉概率是 0.8变异概率是 0.01对于调度问题通常要调高变异概率。原因是调度问题的解空间高度离散变异概率过低会导致种群快速同质化所有个体都收敛到同一个次优解。我一般设为交叉 0.85、变异 0.05然后观察早熟现象是否缓解。判断早熟有个简单方法连续 50 代最优适应度完全没有变化且当前值离你预期最优值还有相当差距说明算法陷入了局部最优。处理手段有两个一是把变异概率提高到 0.1二是改用 adaptive feasible 变异函数让 ga 自动调整变异步长。注意变异概率也不是越大越好超过 0.2 会把好解频繁破坏收敛曲线变成锯齿状。4.3 惩罚系数的灵敏性怎么判断约束是否过紧或过松罚函数系数是调度问题里最讲究“玄学”的部分。系数太小算法会产出大量重复访问或站点供需严重失衡的解系数太大算法把所有精力都放在满足约束上行驶距离反而失控。我建议做一个简单的灵敏度实验把惩罚系数从 100 按 10 倍递增到 10000各跑三次记录平均最优值和约束违背次数。表惩罚系数灵敏度实验记录模板惩罚系数平均最优成本约束违背次数平均运行时间100略低很多较短1000中等少量中等10000偏高几乎为零较长我的做法是取“约束违背次数降到 0 的最小惩罚系数”再留 20% 的余量。这样既保证解的可行性又不至于因为惩罚太大压缩搜索空间。如果你发现结果中频繁出现“某站点被访问 3 次但其他站点完全没被访问”的情况优先怀疑的不是算子问题而是重复访问惩罚设低了。4.4 用多次运行固定随机种子让结果更稳定遗传算法是随机算法同一份代码跑两次结果可能不同。做灵敏度实验时每一组参数至少要跑三次三次的最优值取中位数来比较而不是取最好的一次否则容易把运气当成参数优势。每次运行前调用 rng(seed)用不同的种子记录结果数据才可信。竞赛写论文时有一个容易被质疑的点只给一组最优解而不说明计算过程。用三组不同种子运行并列出最优值区间能有效回应“这个解是不是碰运气跑出来的”这类问题。这也是 matlab 程序源代码类题目的一个隐藏加分项评审老师看到稳定的多组数据比看到一个孤零零的最优解更有说服力。5. 避坑指南遗传算法求解自行车调度最常见的五个翻车现场5.1 运行结果每次都不一样论文没法复现现象同一份代码前后两次运行得到的最优成本相差 10% 以上把结果写进论文后自己再跑都对不上。原因ga 函数默认使用随机数生成器初始化种群只要没有固定随机种子每次运行都是不同的搜索起点最终收敛到不同的局部最优。解决在主程序开头强制 rng(固定数字) 或 rng(default)并在论文中注明使用的种子值。当需要展示算法稳定性时固定三到五个种子分别运行把结果中位数写入正文完整的种子列表附在附录里。这个方法成本为零却能避免“结果不可复现”这个竞赛论文最致命的硬伤。5.2 站点编号出现小数路径完全无法解读现象ga 返回的染色体包含如 7.32、3.81 这样的小数printSchedule 输出的路径中出现无效站点编号。原因调用 ga 时没有设置 IntCon 参数。ga 默认把决策变量当作连续变量处理只有显式声明整数约束才会在交叉和变异后做整数化处理。解决在 ga 参数中传入 IntCon 1:chromLen声明全部变量为整数。这是一个非常容易遗忘的细节很多同学在 fitnessFunc 里写了 round(x) 来取整这不是不行但会破坏遗传算子的连续性导致搜索效率下降。正确做法是让 ga 自己维护整数性。5.3 所有个体都收敛到“不派车”的最优解现象最优成本非常低但检查发现 vehicleUsed 为 0没有任何车辆被派出站点供需完全未处理。原因目标函数里缺少对调度不足的惩罚或者惩罚系数远低于启用一辆车的固定成本算法发现在数学上“什么都不做”反而是最优解。解决在 fitnessFunc 中对站点车辆数偏离理想区间施加高额惩罚确保“不作为”的得分高于“派车但不够完美”的得分。一个稳妥的做法是先手动设置一个最简单的可行调度方案计算其成本然后把惩罚系数调到该成本的 1.5 倍以上迫使算法至少尝试派车。5.4 时间窗约束形同虚设路径顺序完全忽视调度时效现象输出的路径虽然访问了所有站点但在时间上根本不可行——第 3 站的服务时间还没结束路径里已经安排第 7 站因为第 7 站的时间窗早已过期。原因适应度函数里根本没有加入时间累积逻辑只计算了距离并检查了容量和访问次数时间窗口完全没有被编码。解决在 fitnessFunc 中维护一个当前时间变量每次访问站点后累加行驶时间和该站的服务时间如果到达时间超出该站点时间窗累加一个时间惩罚。时间约束是调度问题中最容易被忽略、但在现实中最影响可行性的约束新增一个 100 到 200 行的子函数并不复杂建议在基础版本跑通后立刻补上。5.5 站点数量增加到 50 以后求解时间从 2 分钟暴涨到 2 小时现象20 个站点时算法运行流畅把站点数改为 50 后每次迭代耗时猛增完全无法在竞赛时间内完成。原因染色体长度随站点数线性增长fitnessFunc 中的循环次数也随之增加但更大的问题是搜索空间变成了原来的指数倍500 代的进化次数不足以支撑更大规模的搜索。解决把站点数分阶段增加每次增加 10 个站点观察运行时间和解质量同时调整种群规模和代数。如果 50 站点仍然太慢考虑把单次运行拆成两步先用聚类算法把站点分成若干区域每个区域单独用 ga 求解再合并路径。这是工程上常见的“分而治之”策略通用性很强。6. 从静态到动态把单时段调度扩展成滚动时域方案上面所有的求解都建立在“一次规划、一次执行”的静态假设上也就是每辆车出发后全程不再接收新订单。真实业务里需求是随时间变化的早高峰结束前的最后十分钟某个地铁站突然涌出几百辆车静态方案无法响应。解决这个问题的标准做法是滚动时域调度每 30 分钟重新获取各站点的实时车辆数把“当前状态 短时预测需求”作为输入重新跑一次 ga生成未来 30 分钟的新路径。实现时把基质模型结构体改成动态更新每轮滚动开始时更新 initialBikes、demandBikes把上一轮已经执行完的路径站点从待访问集合中移除。ga 的初始种群可以用上一轮的最优染色体做根基加入少量扰动来加速收敛。这个“热启动”技巧值得单独提一下它通常能把第二轮及之后的求解时间压缩到第一轮的三分之一以下因为新方案与旧方案的大部分路径是重合的。代码层面的改动并不大核心是包一层外部循环horizon 8; % 总共滚动 8 个时段每时段 30 分钟 for t 1:horizon % 更新站点状态 model.initialBikes getCurrentBikes(t); model.demandBikes predictDemand(t); % 如果 t 1用上一轮结果构造初始种群 if t 1 exist(prevX, var) options.InitialPopulationMatrix createWarmStart(prevX, model); end % 调用 ga 求解本时段调度方案 [x, fval] ga((x) fitnessFunc(x, model), chromLen, ... [], [], [], [], lb, ub, [], IntCon, options); prevX x; % 执行本时段的调度路径 executeSchedule(x, model); end这个方案把单次求解变成了一个可以持续运行的闭环系统代价是求解频率提高后计算压力增加但对 30 分钟的调度周期来说20 个站点规模的 ga 求解时间在几秒到几十秒之间完全够用。如果你做的是竞赛论文在最后加入这一章“从静态到动态”的扩展讨论通常能展现你对实际问题复杂度的理解深度。我个人的习惯是动态调度方案的验证不要只看总成本下降了多少而是看“站点溢出的总时长”和“空桩总时长”这两个业务指标。成本曲线下降很漂亮如果溢出时长没有明显减少说明算法优化的方向和真实需求是错位的。这是我带项目时踩过最大的一个坑——把太多精力放在调参让目标函数数值变好却忘了目标函数本身可能没有反映业务痛点。希望这篇笔记能帮你在建模和写代码的每个环节都先想清楚“这个数算出来到底是为了什么”少走这些弯路。希望帮到你。本文还有配套的精品资源点击获取