MATLAB遗传算法求解VRP与VRPTW：从容量约束到时间窗的完整实现

简介资源为MATLAB环境下的车辆路径问题VRP求解代码集合覆盖经典VRP与带时间窗VRPTW场景适合物流优化、运筹学学习者及算法研究人员参考也适合作为本科或研究生相关课程的实验素材。压缩包共35个文件包含32个m函数/脚本、2个mat数据文件及1个asv备份文件整体仅25KB代码按遗传算法GA、模拟退火SA、禁忌搜索TS三类算法组织内附距离矩阵生成、初始解构造、2-opt/swap操作及目标函数计算等模块并提供数据文件便于直接运行验证。已有1521人学习下载。通过研读与运行这些代码可掌握三类全局优化算法在车辆调度中的编码实现与参数调优思路理解邻域搜索、禁忌表、退火降温等机制并可直接修改数据矩阵适配自身配送需求对于想快速上手组合优化和MATLAB算法编程的读者是一份兼具教学与实战价值的紧凑资源。1. 用MATLAB求解VRP和VRPTW一套代码从容量约束走到时间窗很多团队手里的“调度系统”本质只是一张订单表和几辆车的容量数离VRP只差一个距离矩阵和一个解码函数。用MATLAB来写求解器不是因为它能替代商用求解器或调包算法而是因为从模型到代码的距离最短坐标、需求量、时间窗是现成数组遗传算法的选择和交叉用矩阵索引几行就能表达调试时每个中间变量都在工作区里摊开。下面这套代码要给的是一套可直接复现的MATLAB求解VRP及VRPTW骨架先看容量约束下的VRP再扩展带时间窗的VRPTW最后落到怎么验证和调参。适合先要“把算法跑起来”的人更适合想知道一个可行解边界在哪的工程师。2. VRP与VRPTW的模型拆解容量约束和时间窗约束在MATLAB里怎么表达2.1 车辆路径问题的数据组织与解码思路VRP的标准描述是一个配送中心编号0和n个客户编号1到n多辆车从中心出发每辆车有载重上限Q每个客户恰好访问一次目标是总行驶距离最短。论文里常用三维决策变量表达但落到MATLAB里我更愿意围绕数组写代码坐标、需求量、时间窗各存一列距离矩阵用坐标一次算出来。这样换数据时只需要替换矩阵算法解码部分一行不用动。下面是准备数据的最小代码。% 客户数据第1行是配送中心后续每行一个客户 x [0; 36; 52; 18; 41; 28; 47; 32]; y [0; 28; 35; 60; 24; 52; 18; 38]; demand [0; 25; 15; 20; 30; 10; 12; 18]; % 仓库需求量置0 n length(demand) - 1; % 距离矩阵欧氏距离取整模拟公里数 distMat zeros(n1); for i 0:n for j 0:n distMat(i1, j1) round(sqrt((x(i1)-x(j1))^2 (y(i1)-y(j1))^2)); end enddistMat的维度是(n1)×(n1)索引对齐方式是仓库对应第1行第1列客户k对应第k1行第k1列这个约定要贯穿解码函数。欧氏距离取整只是原型阶段做法真实路网数据可以直接把distMat替换为OD矩阵。需求单位必须和车辆载重统一口径否则容量约束会失真距离单位与后文惩罚系数单位也要一致这决定惩罚项的量级。2.1.1 一条染色体怎么表示多辆车编码方式最影响实现难度。一种常见写法是“客户排列加分隔符0”比如[4 2 0 3 1]代表第一辆车服务客户4、2第二辆服务3、1。这种编码直观但交叉时0会跟着参与交换很容易产生两个0或没有0的子代。我更常用的做法是染色体只保存客户排列车辆切分留到解码阶段做从头累加需求量如果加入下一个客户会超过Q就关闭当前车、换一辆新车。这样交叉变异都只作用于一个排列0不进入基因位车辆数也是动态生成的不用提前写死。2.1.2 容量约束的两种策略硬分割与超载惩罚容量约束有两种落地方式。第一种是硬分割解码时严格保证每辆车不超过Q任何时候产生的路径都是容量可行的适应度里不用调容量惩罚系数第二种是软容量允许超载适应度里加上超载量与惩罚系数的乘积适合“车辆稍微超一点也能接受、但尽量不超”的业务。我一般默认硬分割因为多数物流场景超载不合法代码里仍然保留超载量输出口需要软容量时把解码里的容量判断改成load demand Q * 1.1再记录超出量就行。2.2 VRPTW的时间窗约束从“能不能装下”到“能不能按时到”VRPTW在VRP基础上给每个客户i增加三个参数最早可开始服务时间ready_i、最晚可开始服务时间due_i、服务时长service_i。车辆到达早于ready_i要原地等待到达晚于due_i就是迟到。判断一条路径是否可行要从仓库开始逐点递推离开前一个节点的时间加上行驶时间得到到达时间到达时间与ready_i取大得到实际开始服务时间再加上service_i就是离开该节点的时间。一路算下来会得到三个反映解质量的指标总行驶距离、总等待时间、总迟到量。对比维度VRPVRPTW决策内容路径顺序路径顺序 到达时间软硬约束载重 ≤ Q载重 ≤ Q 且 迟到量可控解码复杂度按容量切段按容量切段 时间递推适应度构成总距离总距离 迟到惩罚搜索空间排列 切分时间窗把解空间再切割从表格能看出VRPTW只是在VRP的解码循环里多了一个逐点时间累加所以第3章先实现不带时间窗的遗传算法第4章把decode替换成decodeTW即可。2.2.1 到达时间、等待时间与服务时间的递推关系设车辆从仓库出发时刻为0。从节点i到节点j到达时间 离开i的时刻 distance(i,j)。如果到达j的时刻早于ready_j开始服务时间 ready_j等待时间 ready_j - arrival_j否则开始服务时间就是到达时间。离开时刻 开始服务时间 service_j。注意标准VRPTW的时间窗通常指“最晚开始服务时间”如果业务给的是最晚离开时间录入数据前先换算成due - service否则整个惩罚体系都会偏移。2.2.2 硬时间窗与软时间窗怎么选硬时间窗要求任何迟到都不允许常见于冷链和即时配送软时间窗允许迟到但付出代价。遗传算法里不建议做纯硬时间窗只要一个体有一条路径迟到就淘汰种群会迅速退化到少数排列搜索停止。工程上更稳的做法是软时间窗配大惩罚让适应度里可行解一定优于不可行解算法一边压低迟到一边优化距离。这个思路和容量约束的软接口可以共用一套惩罚机制。3. 用遗传算法写一个MATLAB版VRP求解器从解码到主循环3.1 为什么选遗传算法而不是优化工具箱VRP的可行解是离散排列MATLAB优化工具箱里的ga更适合连续变量优化套到VRP上必须自己写编码、交叉和约束处理实际省不了多少事。自写遗传算法反而可控编码、解码、适应度、遗传算子都是独立函数后面从VRP扩展VRPTW只需要替换解码函数和适应度累加项。遗传算法的成本集中在一个解码函数和一个主循环对几十到几百个客户、几台车的规模能在几十秒内给出接近最优的解足够支撑原型验证和中小规模排线。3.2 解码函数与适应度把排列变成路径和距离解码函数接收一个客户排列perm返回按容量切段后的总行驶距离。切段规则从头累加需求量加入下一个客户不超Q就继续走超了就把当前客户放到下一辆新车。这段代码同时返回超载量overload硬分割下恒等于0但保留接口可以无缝切换软容量模式。function [dist, overload] decode(perm, distMat, demand, cap) n length(perm); dist distMat(1, perm(1)1); % 仓库到第一个客户 load demand(perm(1)1); overload 0; for i 2:n if load demand(perm(i)1) cap dist dist distMat(perm(i-1)1, perm(i)1); load load demand(perm(i)1); else % 当前车辆关闭回仓库新车辆从仓库出发 dist dist distMat(perm(i-1)1, 1); dist dist distMat(1, perm(i)1); load demand(perm(i)1); end end dist dist distMat(perm(n)1, 1); % 最后一辆车回仓库 endperm(i-1)1沿用仓库编号约定只要当前客户还能装入整车路径就是连通的一旦换车先把上一辆车拉回仓库再从仓库出发到新客户中间插入的仓库点就是路径分隔符。适应度取总距离加容量惩罚fitness dist penalty * overload。默认硬分割下overload为0惩罚项不干扰搜索若切换成软容量只需放宽容量判断条件惩罚项会自动起作用。提示overload接口是给软容量扩展用的硬分割下不要因为看到它为0就把代码删掉后面VRPTW的迟到量也要靠同样的返回机制传出来。3.3 遗传算子的两个关键顺序交叉和交换变异选择用锦标赛更稳每次从种群随机抽2个个体把适应度较小的那个放入交配池。交叉算子用顺序交叉Order Crossover而不是单点交叉因为排列编码要求每个客户只能出现一次普通单点交叉会产生重复和遗漏。顺序交叉第一步在父代1上选一个区间并复制给子代第二步从父代2里剔除区间内已选过的客户余下客户按原顺序填进子代空位。变异用两类方式交换两点位置或反转一段子序列。function child orderCrossover(p1, p2) n length(p1); a randi(n-1); b randi([a1, n]); child zeros(1, n); child(a:b) p1(a:b); % 复制父代1的区间 rest p2(~ismember(p2, p1(a:b))); % 父代2剔除重复客户 pos [1:a-1, b1:n]; child(pos) rest(1:length(pos)); enda和b分别是交叉区间的起止位置随机生成但保证abismember负责去掉父代2里已经出现在区间中的客户剩下的客户保持相对顺序填充空位子代仍是完整排列。变异函数类似pos randperm(n,2); c(pos) c(fliplr(pos))实现交换第二类用c(a:b) fliplr(c(a:b))实现反转。这里不需要处理分割点因为分割是解码阶段的事算子只需要维护客户顺序。3.4 主循环、早停条件与参数附表主循环每代算一遍所有个体的适应度记录全局最优然后锦标赛选父代、以概率pc执行顺序交叉、以概率pm执行变异最后把上一代最优个体原样放进下一代防止精英丢失。收敛判断用“最近40代最优适应度无变化”作为早停条件比固定迭代次数更适合现场调参。popSize 100; maxGen 400; pc 0.9; pm 0.1; pop zeros(popSize, n); for i 1:popSize, pop(i,:) randperm(n); end trace zeros(maxGen, 1); for gen 1:maxGen for i 1:popSize [d, o] decode(pop(i,:), distMat, demand, cap); fits(i) d penalty * o; end [bestVal, idx] min(fits); trace(gen) bestVal; if gen 40 trace(gen) trace(gen-40), break; end newPop zeros(popSize, n); newPop(1,:) pop(idx,:); for i 2:2:popSize p1 tournament(pop, fits); p2 tournament(pop, fits); if rand pc c1 orderCrossover(p1, p2); c2 orderCrossover(p2, p1); else c1 p1; c2 p2; end if rand pm, c1 mutateSwap(c1); end if rand pm, c2 mutateSwap(c2); end newPop(i,:) c1; newPop(i1,:) c2; end pop newPop; endpenalty取距离矩阵非0元素平均值乘30量级和总距离匹配popSize少于50容易早熟超过300迭代一次耗时明显增加pc太低会破坏优良片段组合太高则趋向随机搜索0.9是常用起点pm在0.05到0.2之间调整复杂数据从0.1起步。参数建议范围现场调法popSize80200早停过快就把种群放大maxGen300800配合40代无改善早停使用pc0.80.95收敛慢时适当降低交叉率pm0.050.2陷入局部最优时提高变异率penalty2030倍平均距离主要用于软容量扩展接口配套的群体算子也很短tournament用randi(size(pop,1),2,1)抽两个候选再比较适应度mutateSwap用randperm(numel(c),2)取两个位置并交换。function p tournament(pop, fits) k randi(size(pop,1), 2, 1); [~, j] min(fits(k)); p pop(k(j), :); end function c mutateSwap(c) pos randperm(numel(c), 2); c(pos) c(fliplr(pos)); end这两个算子都只操作客户排列不接触仓库分隔符因此不会破坏解码的切分逻辑。到这里不带时间窗的VRP求解器已完整可跑第4章在这个基础上只动解码和适应度。4. 从VRP到VRPTWMATLAB求解VRPTW的约束实现与惩罚调参4.1 解码时逐点推进到达时间VRPTW解码函数必须把时间累加进每条路径。新增输入是tw矩阵n行2列第1列ready第2列due和service列向量。递推变量有三个arrival是到达当前客户的时间start是实际开始服务时间depart是离开时刻。到达时间 上一节点离开时间 两节点距离开始服务时间 max(arrival, ready)离开时刻 start service。车辆切换时新一车从仓库出发depart重置为0这里隐含假设仓库在0时刻可用且无窗口关闭时间。function [dist, lateTotal] decodeTW(perm, distMat, demand, cap, tw, service) n length(perm); dist distMat(1, perm(1)1); load demand(perm(1)1); arr dist; lateTotal max(0, arr - tw(perm(1),2)); depart max(arr, tw(perm(1),1)) service(perm(1)); for i 2:n if load demand(perm(i)1) cap dist dist distMat(perm(i-1)1, perm(i)1); arr depart distMat(perm(i-1)1, perm(i)1); load load demand(perm(i)1); else dist dist distMat(perm(i-1)1, 1) distMat(1, perm(i)1); arr distMat(1, perm(i)1); load demand(perm(i)1); end start max(arr, tw(perm(i),1)); depart start service(perm(i)); lateTotal lateTotal max(0, arr - tw(perm(i),2)); end dist dist distMat(perm(n)1, 1); end这里有三个容易踩的点。一是tw第1列是“最早可开始时间”不是“最早到达时间”所以用max(arr, tw(...))计算开始服务二是due是“最晚开始时间”迟到判断用arr - due不能用start - due否则服务时间被重复算进迟到三是车辆切换时arr从仓库到新客户的距离开始算depart无需恢复因为下一段路程只依赖arr。运行前检查tw列顺序很多现场数据把最早到达时间当成了最早开始服务时间差一个service会让整段路径时间计算失真。4.2 软时间窗惩罚与可行解的引导把decodeTW接进适应度循环容量约束由硬分割保证适应度公式变成总距离 迟到惩罚。其中时间惩罚 omega1 × lateTotal。omega1的取值直接决定搜索方向取太小算法会为省路程而允许大迟到取太大又会让距离优化停滞。我常用的标定方法先用第3章不带时间窗的代码跑一遍得到参考最优距离D0再根据业务能容忍的总迟到量T比如30分钟把omega1设为D0/T。这样“晚到T分钟”的代价约等于“多跑一个完整最优路线的距离”量级匹配算法会自然把迟到压到可接受范围。omega1 D0 / 30; % 容忍总迟到30分钟 for i 1:popSize [d, late] decodeTW(pop(i,:), distMat, demand, cap, tw, service); fits(i) d omega1 * late; enddecodeTW内部已经是硬分割容量不需要额外容量惩罚如果业务改成软容量把overload从decodeTW里一起返回并乘omega2即可。运行后看一个关键现象如果适应度里late项在前50代没有明显下降检查omega1是否被距离量级压住如果late归零但总距离比参考D0大很多说明omega1偏大算法过早牺牲了路径优化。4.3 时间窗引起的两个跑偏问题早熟与长等待第一个问题是早熟。当某个排列恰好把时间窗排通lateTotal变成0以后惩罚项失去梯度种群很快转向距离优化不再探索其它时间窗可行区域。对策是让omega1随迭代衰减每50代乘以0.98拿前250代做“先压时间、后压距离”的搜索。第二个问题是等待时间不计入适应度算法会排出“车辆早早到场、停着等好几个小时”的路径这在标准VRPTW论文里合理但现场很难接受。处理方式是在适应度里加一个小等待系数取0.1~0.3倍的平均边距离把等待时间折成距离损失。权重不能大否则算法会为了躲等待而绕远路让总距离恶化。这两个问题在数据里表现为前者适应度曲线在70代左右突然走平后者总行驶距离看着合理但路径甘特图有大段空白。加一个逐代的[bestVal, meanLate, meanWait]三列输出观察三个数随迭代的变化能快速判断是哪种情况再决定衰减系数和等待系数往哪个方向调。5. 从能跑到跑快MATLAB求解VRPTW的验证与调参技巧5.1 先用小规模穷举证明代码没写错跑真实数据前先构造一个7个客户、2辆车的实例节点数小到可以用全排列枚举验证。枚举时每个排列还要考虑所有可能的车辆切分点因为不同切分点对应不同距离。7! × 2^6约4万次计算MATLAB几秒出结果得到全局最优值。然后把同一实例交给第3章的GA跑10次如果10次都能在50代内复现穷举数值说明解码和距离矩阵索引没有错如果稳定差一个固定值优先查distMat里仓库来回两段是否都算了。这一步只花十分钟能避免后续调参时把代码错误误判成算法问题。5.2 三个诊断指标看遗传算法有没有在干活调参不能只盯总适应度。每10代记录三个数全局最优适应度、平均迟到总量、平均等待时间。正常形态是总适应度阶梯下降lateTotal在50代内快速压到很低或归零如果lateTotal长期不降说明omega1小到被距离淹没按5.1的最优距离重新标定。种群规模先从100起步早停要是出现在前100代把pm提到0.2如果100代内还没出现早停、单代耗时又高把popSize降到80。所有对比都在同一个随机种子下重跑三次取中位数避免单次随机波动误导判断。5.3 收尾加一个2-optGA结果还能再降5%GA的排列搜索擅长找好顺序但不够精细。把每辆车的路径单独拿出来做一轮2-opt遍历路径上所有边对尝试把(i,i1)和(j,j1)重连成(i,j)和(i1,j1)如果重连后总代价变短就接受。对VRPTW接受条件改写成newD omega1*newLate oldD omega1*oldLate否则2-opt可能省几公里却造成大迟到。实现如下。function route twoOptTW(route, D, tw, service) % route: 含仓库端点的单条路径如 [0,3,1,0] improved true; while improved improved false; for i 2:length(route)-2 for j i1:length(route)-1 oldD D(route(i-1)1, route(i)1) D(route(j)1, route(j1)1); newD D(route(i-1)1, route(j)1) D(route(i)1, route(j1)1); oldLate lateOfRoute(route, D, tw, service); newRoute route; newRoute(i:j) fliplr(route(i:j)); newLate lateOfRoute(newRoute, D, tw, service); if newD newLate oldD oldLate route newRoute; improved true; end end end end end function late lateOfRoute(route, D, tw, service) late 0; depart 0; for k 2:length(route)-1 arr depart D(route(k-1)1, route(k)1); start max(arr, tw(route(k),1)); late late max(0, arr - tw(route(k),2)); depart start service(route(k)); end endlateOfRoute独立于GA解码函数只服务2-opt因此循环里使用arr和start分离的写法保证迟到判断口径与decodeTW一致。2-opt结束后把每条车路径拼回去再算一遍总适应度确认局部搜索没有把某辆车的时间窗漏洞放大。改进幅度与数据分布有关客户越密集、距离矩阵方差越大的数据收益越明显通常能挤出5%~10%的距离。本文还有配套的精品资源点击获取