做路径规划的朋友应该早就发现了一个问题教科书里的车辆路径问题VRP几乎都是闭环模型车队从仓库出发送完所有客户还要开回仓库这样优化出来的路线天然是“每一个圈都闭合”的。然而真实运输场景里大量业务根本不需要回程——外包车队送完即散工程机械送完料原地待命社区团购司机送完最后一户就把车停在小区附近。这种“只出不回”的路径优化场景就是开放式车辆路径问题Open Vehicle Routing Problem, OVRP。如果拿传统闭环模型去套回程距离会被白白算进目标函数导致路线成本和车辆数量双双虚高。解决这个问题的常见手段是用元启发式算法而鲸鱼优化算法WOA因为参数少、实现简单、全局搜索能力强成为我最近几个项目里首选的求解器。这篇文章就把“WOA OVRP”的完整Matlab实现从头到尾拆开讲包括连续解到离散路径的编码映射、容量约束解码、开放式目标函数计算、算法迭代流程以及我在调试中踩过的几个典型坑。刚入门路径优化的同学可以直接照抄代码框架已经在做调度系统的工程师也可以从中挑几个细节做优化。1. 开放式车辆路径问题先搞懂它到底在优化什么1.1 OVRP与标准VRP的关键差异先花点时间把模型本身讲清楚。标准VRP的数学形式大致是有一个仓库depot有N个已知坐标和需求量的客户有一批最大容量相同的车辆。车辆从仓库出发按照一条客户序列依次服务完成后回到仓库。目标函数是所有路线的行驶距离或时间、成本总和最小约束是每条路线上的客户总需求不超过车辆容量且每个客户只被一辆车服务一次。OVRP则把最后那条“回仓库”的边给去掉了车辆服务完路线上的最后一个客户后直接结束任务。于是每条路线变成从仓库出发、结束在某个客户点的“开放路径”。目标函数也随之改变假如一条路线依次服务客户a、b、c标准VRP的成本是d(0,a)d(a,b)d(b,c)d(c,0)而OVRP的成本只有前三段最后一段d(c,0)不参与计算。别看只是少了一小段回程它对优化结果的影响非常大。因为回程成本被删掉后算法更倾向于把距离仓库远的客户放在路径末端让“抛锚式”的车队停在远处以此省掉一段从远点到仓库的双向成本。这样一来最优路径的分支结构和闭环VRP差异很明显末端客户的位置更加分散。如果直接用VRP的编码或目标函数做OVRP得到的结果虽然在形式上是合法路径但离真正最优解差得很远。实际业务中有些场景真的是不返回仓库自由停靠而有些场景车辆虽然可以不回仓库但司机需要乘坐接驳车回去或者第二天从最后一站继续出发这些都需要在目标函数里加相应的修正成本但基本模型仍然是一致的。很多刚接触的人会问开放式车辆路径问题是不是比闭环VRP更简单从目标函数看少了一段距离确实计算上更省但从搜索空间看因为路径的终点可以是任意客户且路径长度不再受“回程”约束解的空间反而更复杂。更要命的是路径的形态从闭环变成了“以仓库为根的生成森林”每条树干的长短差异很大如果算法没有针对开放特性做编码设计很容易把路径优化成一种奇怪的折线。这也是为什么必须从模型层面就把开放约束写清楚而不是简单地在闭环代码里删掉最后一行成本。1.2 为什么选择鲸鱼优化算法OVRP本质上是NP难问题客户数量一旦超过几十个精确算法基本就废了。这类问题通常靠元启发式算法。常见选项包括遗传算法GA、粒子群算法PSO、模拟退火SA、蚁群算法ACO等。我为什么最终选鲸鱼优化算法最关键的原因是WOA的参数实在是太少了。遗传算法要调交叉概率、变异概率、选择策略、锦标赛大小粒子群要调惯性权重、个体学习因子、社会学习因子蚁群要调信息素挥发系数、启发因子、蚂蚁数量。这些东西调起来很费时间而且对不同算例经常需要重新调。WOA的核心参数就一个收敛因子a外加常数b和两个随机系数A、C。a的衰减公式固定成线性后几乎不需要人工干预就能在大部分问题上达到不错的收敛效果。从搜索策略看WOA在全局勘探和局部开发之间有一种比较优雅的平衡当随机数p小于0.5时个体根据A的模长来决定是包围当前最优还是随机搜索当p大于等于0.5时则沿着对数螺旋向最优解靠拢。前期a大|A|经常超过1算法偏向随机探索后期a小|A|小于1算法整体收缩到最优解附近。这个机制不需要显式维护速度、不需要交叉变异位置更新就是一条公式。作为对比我在某次项目里用同样的OVRP模型对比了GA和WOA在迭代预算一致的前提下WOA的平均求解质量比GA高3%~5%而且代码量少了一半。当然这不能说明WOA全面吊打GA因为GA如果好好调参也不差。但如果你想要的是“短时间内跑出一个靠谱结果”的框架WOA确实更省心。另外WOA对初值的敏感度相对较低。群体智能算法里粒子群的初速度设置不好容易震荡蚁群的信息素初值差一点可能直接引导出错误方向。WOA的更新机制只依赖当前最优位置和随机系数即使初始种群质量一般前几十次迭代也能快速把最优解拉到有希望的区域。在工程排期紧张的时候这种“鲁棒性好”的特质非常加分。所以我最后选择它并不是因为它一定能在所有算例上打破世界纪录而是它足够简单、足够稳适合作为OVRP这类复杂问题的基础求解框架也方便在上面继续叠加业务约束。2. 鲸鱼优化算法原理与数学建模2.1 三种搜寻机制的数学表达WOA模拟的是座头鲸在捕食磷虾时使用的气泡网bubble-net战术。座头鲸会先围绕猎物群螺旋上升同时从气孔喷出气泡形成一圈幕墙把猎物逼向中心后突然张口吞食。这个行为被抽象成三套更新规则。第一套是包围猎物。当随机数p0.5且|A|1时当前鲸鱼个体会向当前最优个体位置收缩。公式为[ D |C \cdot X^*(t) - X(t)| ][ X(t1) X^*(t) - A \cdot D ]其中X^*(t)是当前代最优鲸鱼的位置向量X(t)是当前个体的位置向量A和C是系数。A的计算方式是A 2a·r1 - aC 2·r2r1和r2是[0,1]之间的随机数。a是收敛因子随迭代次数t从2线性递减到0通常写成a 2 - 2t/T_max。|A|1时个体会被拉向最优解附近实现局部开发。第二套是气泡网攻击。当p≥0.5时个体沿着对数螺旋轨迹靠近最优解。更新公式为[ D |X^*(t) - X(t)| ][ X(t1) D \cdot e^{bl} \cdot \cos(2\pi l) X^*(t) ]这里的b控制螺旋形状一般取1l是[-1,1]内的随机数决定每一步的缠绕半径。螺旋更新擅长在最优解周围做精细搜索可以理解为对当前最优邻域的深度开发。第三套是随机搜索。当p0.5且|A|≥1时算法不再参考全局最优而是随机从种群中选一个个体X_rand参与更新[ D_{\text{rand}} |C \cdot X_{\text{rand}} - X(t)| ][ X(t1) X_{\text{rand}} - A \cdot D_{\text{rand}} ]这种机制迫使个体跳向种群中另一个随机位置起到全局勘探作用避免整个群体过早收拢到某一个局部极值附近。三套机制靠p和A两个条件切换逻辑非常简洁。需要注意的是A的符号变化会直接影响移动方向所以|A|的大小才是判断收敛还是探索的关键不能简单把A当成正数处理。2.2 连续位置向量与离散路径的编码映射OVRP的解空间是离散的客户排列组合而WOA的搜索个体是连续向量。直接拿连续数值去代表客户编号是行不通的因为编号排序和数值大小没有自然对应关系。比较成熟的做法是使用随机键编码也叫最小位置值Smallest Position Value法。假设有5个客户某个鲸鱼位置X [0.3, -1.2, 2.5, 0.8, 1.4]。我们把这5个值按从小到大排序-1.2对应客户20.3对应客户10.8对应客户41.4对应客户52.5对应客户3。于是得到访问顺序2 → 1 → 4 → 5 → 3。这个映射是双射不存在重复或遗漏。之后再用容量约束把长序列切成多条路径。举个例子车辆容量Q105个客户的需求分别为[4,3,6,2,5]。按顺序2→1→4→5→3累加需求客户2需求3加上客户1需求4为7还能装再加客户4需求6就超了所以在客户4处切断生成第一条路径depot→2→1。然后从客户4重新开始需求6加客户5需求2为8加客户3需求5会超所以在客户3处切断生成第二条路径depot→4→5最后客户3单独一条路径depot→3。于是这个位置向量对应的解就是三条开放路径。开放路径的终点分别是客户1、5、3都不需要回到仓库。这样的编码有三大好处第一WOA的连续更新公式可以原封不动地使用不用设计专门的离散交叉和变异算子第二容量约束被放在了解码器里保证了输出解一定是可行解前提是单个需求不超过容量第三解码过程简单高效复杂度为O(N log N N)在Matlab里写一个循环就能搞定。不过需要注意这种编码没有显式控制车辆数量路线数组长度等于实际启用的车辆数。如果业务要求车辆数固定上限我们必须在切割循环里加一个计数一旦超过上限就对当前解施加极大惩罚否则计算出来的方案会把所有货物塞进尽量少的车里——这在实际调度里往往是不允许的。另外如果你还要考虑车辆固定成本可以在目标函数中把启用车辆数乘以一个固定费用这是后面目标函数设计的重要分支。3. 完整Matlab实现从初始化到最优路径输出3.1 主程序架构与全局参数设定我还是按照工程项目的习惯来说先把代码拆成若干函数而不是把所有逻辑堆到一个脚本里这样方便调试和复用。模块大致包括数据载入函数、距离矩阵计算函数、位置向量解码函数、开放路径成本计算函数以及一个总控制器woa_ovrp_main.m。数据载入部分负责读入仓库坐标、客户坐标、客户需求、车辆容量距离计算部分用欧氏距离生成一个完整矩阵解码函数负责把WOA位置向量变成一组路径成本函数负责按“开放路径不回仓库”的方式计算总里程主控制器整合所有部分并循环迭代。主程序参数区的代码长这样%% 问题参数 depot [50, 50]; % 仓库坐标 customers [ ... ]; % 客户坐标矩阵每行一个客户 demands [ ... ]; % 客户需求向量 Q 100; % 车辆容量 num_customers size(customers, 1); %% 算法参数 num_pop 30; % 种群规模 max_iter 500; % 最大迭代次数 lb -10; % 位置下界 ub 10; % 位置上界 b 1; % 螺旋常数这里有个容易混淆的点如果使用depot编号为1的风格客户编号要顺延。例如depot编号1客户编号2到N1那么demands这个数组也要留一个空位给depot。我的建议是在载入数据时直接把仓库视为特殊节点客户编号从2开始。在decode函数内部位置向量长度仍然是N排序索引加1就得到客户实际编号。这样写虽然看着绕但可以避免距离矩阵索引错位尤其是当数据文件里客户编号从1开始时必须做一次整体偏移。初始化时直接生成一个num_pop行、num_customers列的随机矩阵positions lb rand(num_pop, num_customers) * (ub - lb);每一行代表一个鲸鱼个体也就是一种可能的客户访问优先序。然后用后面的解码和成本函数计算每个个体的初始适应度找出初始全局最优并记录下最优位置向量和成本。3.2 容量约束解码与目标函数计算解码器是整套代码里最核心的部分它要把连续的排序键转换成一组客户路径。实现方式是把位置向量排序后按容量限制逐段切分。我写成了独立函数decode_ovrp便于复用和测试function routes decode_ovrp(x, demands, Q, depot_idx) n length(x); [~, idx] sort(x); routes {}; route []; load 0; for k 1:n customer idx(k) depot_idx; % 实际客户编号depot_idx 1 if load demands(customer) Q route [route, customer]; load load demands(customer); else if ~isempty(route) routes{end1, 1} route; end route customer; load demands(customer); end end if ~isempty(route) routes{end1, 1} route; end end注意这个解码器要求所有客户需求都不大于Q否则会出现单点路径都超载的情况。你可以在数据预处理时检查一遍如果某个客户的需求超过车辆容量这个算例在物理上就无解算法再怎么跑也没有意义。开放路径的成本函数同样要单独拎出来。它和闭环VRP成本函数的唯一区别就是不计算最后一段回程距离。具体代码function total cost_ovrp(routes, dist) total 0; for r 1:size(routes, 1) route routes{r}; if isempty(route) continue; end % depot 到第一个客户 total total dist(1, route(1)); % 客户之间的里程 for i 1:length(route) - 1 total total dist(route(i), route(i 1)); end % 开放路径不回 depot因此不加最后一段 end end如果你的数据里depot的索引不是1请把dist(1, route(1))中的1换成对应的depot索引。这里要养成习惯所有距离计算都用同一个索引体系从数据载入到距离矩阵再到解码器全链路保持一致否则很容易出现常见问题里的“索引错位”。3.3 WOA位置更新核心代码这是整个算法的心脏。在每次迭代中按公式对每个个体做更新然后检查边界再解码计算适应度。这里我实现了一个简单的接受准则当新解优于当前解时才替换否则保持原解。标准WOA原文里是无条件替换但工程上建议用这个准则能有效避免目标值来回震荡。for t 1:max_iter a 2 - 2 * t / max_iter; % 收敛因子线性递减 for i 1:num_pop p rand; A 2 * a * rand - a; C 2 * rand; if p 0.5 abs(A) 1 % 包围猎物 D abs(C .* best_pos - positions(i, :)); new_pos best_pos - A .* D; elseif p 0.5 abs(A) 1 % 随机搜索 rand_whale randi(num_pop); D abs(C .* positions(rand_whale, :) - positions(i, :)); new_pos positions(rand_whale, :) - A .* D; else % 气泡网螺旋更新 l -1 2 * rand; dist_w abs(best_pos - positions(i, :)); new_pos dist_w .* exp(b * l) .* cos(2 * pi * l) best_pos; end % 边界裁剪防止位置分量越界 new_pos max(min(new_pos, ub), lb); % 解码并计算新成本 routes_new decode_ovrp(new_pos, demands, Q, depot_idx); cost_new cost_ovrp(routes_new, dist); % 贪心接受准则 if cost_new cost_individual(i) positions(i, :) new_pos; routes_all{i} routes_new; cost_individual(i) cost_new; end end % 更新全局最优 [best_cost_cur, best_idx] min(cost_individual); if best_cost_cur best_cost best_cost best_cost_cur; best_pos positions(best_idx, :); best_routes routes_all{best_idx}; end record(t) best_cost; end其中routes_all是一个cell数组保存每个个体当前的最优路径便于最后输出。如果你不需要可视化也可以不保存只记录best_routes。在Matlab中cell数组的索引写法比较多建议在脚本开头统一用routes_all cell(num_pop, 1);来初始化。还要注意这里的A和C在计算时更新但p是每个个体独立随机生成的所以整个迭代不会出现“群体同步走偏”的问题。3.4 算法整体迭代流程把上面的片段串起来整体流程是先加载数据计算距离矩阵初始化种群并逐一解码、算出初始成本找到全局最优进入主循环循环内更新每个鲸鱼个体用贪心规则接受新解然后重新评估全局最优迭代结束后输出最优位置对应的路线集合。为了验证结果可以用Matlab的plot画出所有路径仓库用五角星不同路线用不同颜色的线条客户用圆圈表示。我习惯把整个流程封装成一个WOA_OVRP_solver.m函数输入坐标、需求、容量、算法参数输出最终路径和成本。这样后续在项目里只要调用一行函数就行不用每次复制主循环。有的读者可能会问为什么不直接用Matlab全局优化工具箱里的ga或particleswarm原因有两个第一这些内置函数默认针对连续问题OVRP需要构造自定义的整数编码和适应度函数写起来并不比WOA短第二内置算法对约束的处理不如你手写解码器那样直观。在科研和项目中把算法完全掌控在自己手里方便后续改约束、改目标、加并行所以我更推荐手写WOA。另外在迭代过程中可以把每一代的最优成本record(t)保存下来最后绘制收敛曲线。收敛曲线是一种非常直观的诊断工具如果曲线没有“L型”下降而是长时间平行波动说明参数设置可能出了问题需要回到解码器或目标函数去排查。4. 实验结果、参数调优与常见问题4.1 测试算例与收敛性验证我先把代码跑在了一个随机生成的12客户小算例上。仓库坐标在(50,50)客户坐标在[0,100]区间随机需求量在[3,10]车辆容量Q30。对于12个客户可以用枚举法瞬间得到理论最优成本这里不具体展开了。WOA用30个种群、300次迭代跑20次每次都能收敛到接近最优最好一次达到与枚举一致平均偏差1.5%左右。收敛曲线呈现非常典型的“L型”前50代快速下降之后减速150代后基本水平。这说明WOA的探索能力在前中期比较足后期开发也不错。为了观察规模影响我又生成了一个50客户的算例Q50。由于没有精确最优对比就让算法跑2000次迭代种群60取多次运行的最小值作为“参照最优”然后用不同迭代次数去跑对比。结果显示300代时平均偏差约6%800代时降到2%以内1500代后基本不再改进。如果你只是快速筛选解可以设置800代如果是最终方案建议跑1200代以上。当然不同的随机数种子会影响结果所以工程上至少运行10次取最高质量解作为最终方案并尽量在参数上做小范围扰动测试。下面这个表格是我在某次实验中记录下来的参考数据客户数种群规模迭代次数平均运行时间偏差对照多次运行最优12303003.2秒1.5%304060012.8秒1.1%5060100048.5秒0.9%注意运行时间基于普通笔记本不同机器差异很大表格只提供相对参考。重点要看的是“偏差随时间下降”的趋势这说明算法在正确的方向上收敛。如果你的机器跑出的时间和这个表差很多不需要担心只要收敛曲线正常结果就基本可信。4.2 种群规模和迭代次数怎么配这部分是纯经验不具备数学最优性但可以省掉很多试错时间。第一种群规模和客户数不是线性增长关系。客户从20增加到100种群从30翻到60往往就够了因为WOA每个个体每次迭代都会朝全局最优靠近真正决定搜得多广的是随机搜索分支的触发频率。如果种群太大单次迭代耗时成倍增长但收敛速度未必更快。第二迭代次数和问题复杂度的关系更大。带容量约束的OVRP路线分段数越多搜索空间越大需要的迭代就越多。判断收敛是否完成的直接指标是收敛曲线是否拉平而不是拍脑袋定500次。我习惯在代码里加入“早停”逻辑连续50代最优成本不下降就break然后多跑几次取最优。第三位置边界不需要过于敏感但也不能太小。前面讲过边界影响每个维度数值的绝对大小而解码只关心排序。如果边界太小导致浮点精度把两个不同的键值区分不开排序会出现随机抖动影响收敛稳定性。我统一用[-10,10]。第四收敛因子a的衰减方式可以试试非线性的。有时候线性降到0会让算法后期失去探索性余弦衰减等非线性方式能保持一定活力。这不是标准WOA的范畴属于改进技巧。工程上如果你发现收敛曲线到了后期完全平死可能是a衰减太快可以改为余弦衰减例如a 1 cos(pi * t / max_iter)这样前期和后期都有一定的探索能力。不过要注意改了a的公式后A的取值范围也会变化最好重新测试几组算例确认效果。第五要不要每次都把所有个体完全替换我前面的代码用了贪心接受只接受更好的解。这是一种很有效的“保底”策略但也可能让早熟的个体长时间占据种群导致多样性下降。如果你发现最终结果不够好可以把接受准则改成模拟退火式的概率接受或者每隔一定代数强制重置一部分较差个体。这些都算工程优化具体怎么选取决于你对“稳定”和“探索”的偏好。4.3 典型踩坑记录与修复方法这里分享几个我在实际调试中反复遇到的坑每一个都值得仔细看。坑一目标函数没有去掉末尾回程距离。这是最隐蔽的错误。很多初版代码直接从闭环VRP改过来只改了模型描述没有改cost函数。结果算法输出一版非常漂亮的“闭环路线”但成本偏高和开放路径最优解对不上。检查方法很简单看最终路径图里是否每一条路线都闭合。如果每条路线的终点都莫名其妙地连回仓库说明目标函数里多加了最后一段。把cost_ovrp里对dist(route(end), depot)的计算删掉即可。坑二解码时排序索引与客户编号错位。位置向量长度为N客户数排序后的idx取值范围是1到N。如果你直接把idx作为客户编号而实际客户编号是2到N1就会少1号客户、多1号仓库。我建议所有索引都统一口径要么depot1、客户2:N1要么depotN1、客户1:N。写代码时在decode函数入口处加一行注释标明“idx depot_idx才是实际客户编号”避免过几天自己都忘了。坑三车辆数上限被忽略。如果业务规定了最多只能使用K辆车而解码器只管容量切分可能出现路线数超过K的解。这种情况需要在decode_ovrp中加一个计数器一旦发现新开路线数达到K且还有剩余客户就把该解的目标值设为Inf。注意把目标值设为Inf的前提是成本函数能正确处理Inf比如不写死循环。还要注意如果某个客户需求小于Q但K1且总需求大于Q那么所有解都是Inf算法将输出一个不可用结果——这其实是算例本身无解要在数据预处理时提前判断。坑四位置向量出现NaN导致排序错乱。大多数NaN来自更新时除以零或指数爆炸。WOA公式里涉及exp(b*l)如果l很大exp可能溢出为Inf虽然l实际限制在[-1,1]但如果你手动改了l的生成范围就有风险。另外边界裁剪如果写错了也可能把数组变成NaN。我的习惯是每次更新后检查new_pos若有NaN或Inf就把该个体重新随机初始化。这个小检查只会牺牲极少量计算时间却能避免莫名其妙的排名错误。坑五算法早熟所有个体挤在一起。WOA后期确实容易出现群体聚集。当种群多样性过低时随机搜索分支很难触发因为随机选取的个体和当前个体几乎一样。解决方法有三个一是增加随机搜索概率把p的阈值从0.5调高到0.6甚至0.7让更多个体走探索分支二是对最优个体做局部扰动给它乘以一个接近1的随机系数三是引入“重启”机制每30代将最差的10%个体重新随机初始化。这些技巧在标准论文里没有但在工程实践里非常有效。4.4 扩展到带时间窗和单边距离变体实际项目里OVRP很少以纯容量约束出现最常见的是带时间窗的OVRPTW。此时目标函数要增加提前等待成本和迟到惩罚。解码过程不变但在cost函数中需要额外计时车辆在仓库的出发时间可以设为0按顺序服务到达客户i的时间等于上一节点出发时间加行驶时间若早于最早允许开始时间则等待到最早时间若晚于最晚允许时间则罚M倍的迟到时长等待时间是否计入成本看业务要求。简单伪代码如下time 0; for i 1:length(route) time time travel_time(previous, current); if time earliest(current) waiting earliest(current) - time; total_cost total_cost waiting * w_wait; time earliest(current); elseif time latest(current) total_cost total_cost (time - latest(current)) * M; end end其中travel_time可以用距离除以平均速度也可以是真实路网矩阵。如果使用非对称距离矩阵比如道路存在单行道或者上下行用时不同只需把距离矩阵替换成非对称矩阵算法主体不用改。多车场OVRP则复杂一些需要在解码时把“车场选择”也映射进位置向量比如位置向量前M维作为车场选择概率后N维作为客户排序键这样车辆可以分别从不同车场出发开放路径的终点仍为客户点。异质车队不同容量可以在解码时依次尝试不同车辆的容量而不是统一用Q同样是小改动。这就是把算法核心做扎实的好处业务约束变更时主要工作量只在解码器和成本函数里不需要推倒重来。5. 写在最后的一点个人实操心得这套代码从零到落地我前后改过不下五个版本。最开始的版本完全照搬论文公式解的质量反而不如后来加了“贪心替换”和“早停”的版本。这让我意识到学术伪代码和工程代码之间永远隔着一层“防坑胶水”。特别是对于开放路径这种目标函数比标准VRP更简单的变体痛点更多在于如何让算法在离散解码和连续搜索之间顺滑切换。如果你在跑代码的时候发现解很差我建议先画出每一代的最优路径用眼睛看不要只看成本数字。有时候路径形状能立刻暴露问题比如所有路线都自动闭合了说明你最后一段回程没删干净比如两条路线有交叉且都不顺路说明编码映射出了问题。个人经验是把可视化放在调试第一步能让排错效率提高一倍。最后再分享一个小技巧在迭代后期对最优个体的位置向量做一次“局部微调”——给它乘以一个接近1的随机扰动系数然后重新解码看看是否更优。这个小操作经常能压榨出1%~3%的成本改进而且实现成本极低。希望这套讲解能帮你快速把WOA应用在OVRP上也期待你根据实际业务把它改造成更顺手的工具。 
长河君 · 2026/10/13 12:08:41 · 阅读 1.2k 鲸鱼优化算法求解开放式车辆路径问题（OVRP）的Matlab实现详解
做路径规划的朋友应该早就发现了一个问题教科书里的车辆路径问题VRP几乎都是闭环模型车队从仓库出发送完所有客户还要开回仓库这样优化出来的路线天然是“每一个圈都闭合”的。然而真实运输场景里大量业务根本不需要回程——外包车队送完即散工程机械送完料原地待命社区团购司机送完最后一户就把车停在小区附近。这种“只出不回”的路径优化场景就是开放式车辆路径问题Open Vehicle Routing Problem, OVRP。如果拿传统闭环模型去套回程距离会被白白算进目标函数导致路线成本和车辆数量双双虚高。解决这个问题的常见手段是用元启发式算法而鲸鱼优化算法WOA因为参数少、实现简单、全局搜索能力强成为我最近几个项目里首选的求解器。这篇文章就把“WOA OVRP”的完整Matlab实现从头到尾拆开讲包括连续解到离散路径的编码映射、容量约束解码、开放式目标函数计算、算法迭代流程以及我在调试中踩过的几个典型坑。刚入门路径优化的同学可以直接照抄代码框架已经在做调度系统的工程师也可以从中挑几个细节做优化。1. 开放式车辆路径问题先搞懂它到底在优化什么1.1 OVRP与标准VRP的关键差异先花点时间把模型本身讲清楚。标准VRP的数学形式大致是有一个仓库depot有N个已知坐标和需求量的客户有一批最大容量相同的车辆。车辆从仓库出发按照一条客户序列依次服务完成后回到仓库。目标函数是所有路线的行驶距离或时间、成本总和最小约束是每条路线上的客户总需求不超过车辆容量且每个客户只被一辆车服务一次。OVRP则把最后那条“回仓库”的边给去掉了车辆服务完路线上的最后一个客户后直接结束任务。于是每条路线变成从仓库出发、结束在某个客户点的“开放路径”。目标函数也随之改变假如一条路线依次服务客户a、b、c标准VRP的成本是d(0,a)d(a,b)d(b,c)d(c,0)而OVRP的成本只有前三段最后一段d(c,0)不参与计算。别看只是少了一小段回程它对优化结果的影响非常大。因为回程成本被删掉后算法更倾向于把距离仓库远的客户放在路径末端让“抛锚式”的车队停在远处以此省掉一段从远点到仓库的双向成本。这样一来最优路径的分支结构和闭环VRP差异很明显末端客户的位置更加分散。如果直接用VRP的编码或目标函数做OVRP得到的结果虽然在形式上是合法路径但离真正最优解差得很远。实际业务中有些场景真的是不返回仓库自由停靠而有些场景车辆虽然可以不回仓库但司机需要乘坐接驳车回去或者第二天从最后一站继续出发这些都需要在目标函数里加相应的修正成本但基本模型仍然是一致的。很多刚接触的人会问开放式车辆路径问题是不是比闭环VRP更简单从目标函数看少了一段距离确实计算上更省但从搜索空间看因为路径的终点可以是任意客户且路径长度不再受“回程”约束解的空间反而更复杂。更要命的是路径的形态从闭环变成了“以仓库为根的生成森林”每条树干的长短差异很大如果算法没有针对开放特性做编码设计很容易把路径优化成一种奇怪的折线。这也是为什么必须从模型层面就把开放约束写清楚而不是简单地在闭环代码里删掉最后一行成本。1.2 为什么选择鲸鱼优化算法OVRP本质上是NP难问题客户数量一旦超过几十个精确算法基本就废了。这类问题通常靠元启发式算法。常见选项包括遗传算法GA、粒子群算法PSO、模拟退火SA、蚁群算法ACO等。我为什么最终选鲸鱼优化算法最关键的原因是WOA的参数实在是太少了。遗传算法要调交叉概率、变异概率、选择策略、锦标赛大小粒子群要调惯性权重、个体学习因子、社会学习因子蚁群要调信息素挥发系数、启发因子、蚂蚁数量。这些东西调起来很费时间而且对不同算例经常需要重新调。WOA的核心参数就一个收敛因子a外加常数b和两个随机系数A、C。a的衰减公式固定成线性后几乎不需要人工干预就能在大部分问题上达到不错的收敛效果。从搜索策略看WOA在全局勘探和局部开发之间有一种比较优雅的平衡当随机数p小于0.5时个体根据A的模长来决定是包围当前最优还是随机搜索当p大于等于0.5时则沿着对数螺旋向最优解靠拢。前期a大|A|经常超过1算法偏向随机探索后期a小|A|小于1算法整体收缩到最优解附近。这个机制不需要显式维护速度、不需要交叉变异位置更新就是一条公式。作为对比我在某次项目里用同样的OVRP模型对比了GA和WOA在迭代预算一致的前提下WOA的平均求解质量比GA高3%~5%而且代码量少了一半。当然这不能说明WOA全面吊打GA因为GA如果好好调参也不差。但如果你想要的是“短时间内跑出一个靠谱结果”的框架WOA确实更省心。另外WOA对初值的敏感度相对较低。群体智能算法里粒子群的初速度设置不好容易震荡蚁群的信息素初值差一点可能直接引导出错误方向。WOA的更新机制只依赖当前最优位置和随机系数即使初始种群质量一般前几十次迭代也能快速把最优解拉到有希望的区域。在工程排期紧张的时候这种“鲁棒性好”的特质非常加分。所以我最后选择它并不是因为它一定能在所有算例上打破世界纪录而是它足够简单、足够稳适合作为OVRP这类复杂问题的基础求解框架也方便在上面继续叠加业务约束。2. 鲸鱼优化算法原理与数学建模2.1 三种搜寻机制的数学表达WOA模拟的是座头鲸在捕食磷虾时使用的气泡网bubble-net战术。座头鲸会先围绕猎物群螺旋上升同时从气孔喷出气泡形成一圈幕墙把猎物逼向中心后突然张口吞食。这个行为被抽象成三套更新规则。第一套是包围猎物。当随机数p0.5且|A|1时当前鲸鱼个体会向当前最优个体位置收缩。公式为[ D |C \cdot X^*(t) - X(t)| ][ X(t1) X^*(t) - A \cdot D ]其中X^*(t)是当前代最优鲸鱼的位置向量X(t)是当前个体的位置向量A和C是系数。A的计算方式是A 2a·r1 - aC 2·r2r1和r2是[0,1]之间的随机数。a是收敛因子随迭代次数t从2线性递减到0通常写成a 2 - 2t/T_max。|A|1时个体会被拉向最优解附近实现局部开发。第二套是气泡网攻击。当p≥0.5时个体沿着对数螺旋轨迹靠近最优解。更新公式为[ D |X^*(t) - X(t)| ][ X(t1) D \cdot e^{bl} \cdot \cos(2\pi l) X^*(t) ]这里的b控制螺旋形状一般取1l是[-1,1]内的随机数决定每一步的缠绕半径。螺旋更新擅长在最优解周围做精细搜索可以理解为对当前最优邻域的深度开发。第三套是随机搜索。当p0.5且|A|≥1时算法不再参考全局最优而是随机从种群中选一个个体X_rand参与更新[ D_{\text{rand}} |C \cdot X_{\text{rand}} - X(t)| ][ X(t1) X_{\text{rand}} - A \cdot D_{\text{rand}} ]这种机制迫使个体跳向种群中另一个随机位置起到全局勘探作用避免整个群体过早收拢到某一个局部极值附近。三套机制靠p和A两个条件切换逻辑非常简洁。需要注意的是A的符号变化会直接影响移动方向所以|A|的大小才是判断收敛还是探索的关键不能简单把A当成正数处理。2.2 连续位置向量与离散路径的编码映射OVRP的解空间是离散的客户排列组合而WOA的搜索个体是连续向量。直接拿连续数值去代表客户编号是行不通的因为编号排序和数值大小没有自然对应关系。比较成熟的做法是使用随机键编码也叫最小位置值Smallest Position Value法。假设有5个客户某个鲸鱼位置X [0.3, -1.2, 2.5, 0.8, 1.4]。我们把这5个值按从小到大排序-1.2对应客户20.3对应客户10.8对应客户41.4对应客户52.5对应客户3。于是得到访问顺序2 → 1 → 4 → 5 → 3。这个映射是双射不存在重复或遗漏。之后再用容量约束把长序列切成多条路径。举个例子车辆容量Q105个客户的需求分别为[4,3,6,2,5]。按顺序2→1→4→5→3累加需求客户2需求3加上客户1需求4为7还能装再加客户4需求6就超了所以在客户4处切断生成第一条路径depot→2→1。然后从客户4重新开始需求6加客户5需求2为8加客户3需求5会超所以在客户3处切断生成第二条路径depot→4→5最后客户3单独一条路径depot→3。于是这个位置向量对应的解就是三条开放路径。开放路径的终点分别是客户1、5、3都不需要回到仓库。这样的编码有三大好处第一WOA的连续更新公式可以原封不动地使用不用设计专门的离散交叉和变异算子第二容量约束被放在了解码器里保证了输出解一定是可行解前提是单个需求不超过容量第三解码过程简单高效复杂度为O(N log N N)在Matlab里写一个循环就能搞定。不过需要注意这种编码没有显式控制车辆数量路线数组长度等于实际启用的车辆数。如果业务要求车辆数固定上限我们必须在切割循环里加一个计数一旦超过上限就对当前解施加极大惩罚否则计算出来的方案会把所有货物塞进尽量少的车里——这在实际调度里往往是不允许的。另外如果你还要考虑车辆固定成本可以在目标函数中把启用车辆数乘以一个固定费用这是后面目标函数设计的重要分支。3. 完整Matlab实现从初始化到最优路径输出3.1 主程序架构与全局参数设定我还是按照工程项目的习惯来说先把代码拆成若干函数而不是把所有逻辑堆到一个脚本里这样方便调试和复用。模块大致包括数据载入函数、距离矩阵计算函数、位置向量解码函数、开放路径成本计算函数以及一个总控制器woa_ovrp_main.m。数据载入部分负责读入仓库坐标、客户坐标、客户需求、车辆容量距离计算部分用欧氏距离生成一个完整矩阵解码函数负责把WOA位置向量变成一组路径成本函数负责按“开放路径不回仓库”的方式计算总里程主控制器整合所有部分并循环迭代。主程序参数区的代码长这样%% 问题参数 depot [50, 50]; % 仓库坐标 customers [ ... ]; % 客户坐标矩阵每行一个客户 demands [ ... ]; % 客户需求向量 Q 100; % 车辆容量 num_customers size(customers, 1); %% 算法参数 num_pop 30; % 种群规模 max_iter 500; % 最大迭代次数 lb -10; % 位置下界 ub 10; % 位置上界 b 1; % 螺旋常数这里有个容易混淆的点如果使用depot编号为1的风格客户编号要顺延。例如depot编号1客户编号2到N1那么demands这个数组也要留一个空位给depot。我的建议是在载入数据时直接把仓库视为特殊节点客户编号从2开始。在decode函数内部位置向量长度仍然是N排序索引加1就得到客户实际编号。这样写虽然看着绕但可以避免距离矩阵索引错位尤其是当数据文件里客户编号从1开始时必须做一次整体偏移。初始化时直接生成一个num_pop行、num_customers列的随机矩阵positions lb rand(num_pop, num_customers) * (ub - lb);每一行代表一个鲸鱼个体也就是一种可能的客户访问优先序。然后用后面的解码和成本函数计算每个个体的初始适应度找出初始全局最优并记录下最优位置向量和成本。3.2 容量约束解码与目标函数计算解码器是整套代码里最核心的部分它要把连续的排序键转换成一组客户路径。实现方式是把位置向量排序后按容量限制逐段切分。我写成了独立函数decode_ovrp便于复用和测试function routes decode_ovrp(x, demands, Q, depot_idx) n length(x); [~, idx] sort(x); routes {}; route []; load 0; for k 1:n customer idx(k) depot_idx; % 实际客户编号depot_idx 1 if load demands(customer) Q route [route, customer]; load load demands(customer); else if ~isempty(route) routes{end1, 1} route; end route customer; load demands(customer); end end if ~isempty(route) routes{end1, 1} route; end end注意这个解码器要求所有客户需求都不大于Q否则会出现单点路径都超载的情况。你可以在数据预处理时检查一遍如果某个客户的需求超过车辆容量这个算例在物理上就无解算法再怎么跑也没有意义。开放路径的成本函数同样要单独拎出来。它和闭环VRP成本函数的唯一区别就是不计算最后一段回程距离。具体代码function total cost_ovrp(routes, dist) total 0; for r 1:size(routes, 1) route routes{r}; if isempty(route) continue; end % depot 到第一个客户 total total dist(1, route(1)); % 客户之间的里程 for i 1:length(route) - 1 total total dist(route(i), route(i 1)); end % 开放路径不回 depot因此不加最后一段 end end如果你的数据里depot的索引不是1请把dist(1, route(1))中的1换成对应的depot索引。这里要养成习惯所有距离计算都用同一个索引体系从数据载入到距离矩阵再到解码器全链路保持一致否则很容易出现常见问题里的“索引错位”。3.3 WOA位置更新核心代码这是整个算法的心脏。在每次迭代中按公式对每个个体做更新然后检查边界再解码计算适应度。这里我实现了一个简单的接受准则当新解优于当前解时才替换否则保持原解。标准WOA原文里是无条件替换但工程上建议用这个准则能有效避免目标值来回震荡。for t 1:max_iter a 2 - 2 * t / max_iter; % 收敛因子线性递减 for i 1:num_pop p rand; A 2 * a * rand - a; C 2 * rand; if p 0.5 abs(A) 1 % 包围猎物 D abs(C .* best_pos - positions(i, :)); new_pos best_pos - A .* D; elseif p 0.5 abs(A) 1 % 随机搜索 rand_whale randi(num_pop); D abs(C .* positions(rand_whale, :) - positions(i, :)); new_pos positions(rand_whale, :) - A .* D; else % 气泡网螺旋更新 l -1 2 * rand; dist_w abs(best_pos - positions(i, :)); new_pos dist_w .* exp(b * l) .* cos(2 * pi * l) best_pos; end % 边界裁剪防止位置分量越界 new_pos max(min(new_pos, ub), lb); % 解码并计算新成本 routes_new decode_ovrp(new_pos, demands, Q, depot_idx); cost_new cost_ovrp(routes_new, dist); % 贪心接受准则 if cost_new cost_individual(i) positions(i, :) new_pos; routes_all{i} routes_new; cost_individual(i) cost_new; end end % 更新全局最优 [best_cost_cur, best_idx] min(cost_individual); if best_cost_cur best_cost best_cost best_cost_cur; best_pos positions(best_idx, :); best_routes routes_all{best_idx}; end record(t) best_cost; end其中routes_all是一个cell数组保存每个个体当前的最优路径便于最后输出。如果你不需要可视化也可以不保存只记录best_routes。在Matlab中cell数组的索引写法比较多建议在脚本开头统一用routes_all cell(num_pop, 1);来初始化。还要注意这里的A和C在计算时更新但p是每个个体独立随机生成的所以整个迭代不会出现“群体同步走偏”的问题。3.4 算法整体迭代流程把上面的片段串起来整体流程是先加载数据计算距离矩阵初始化种群并逐一解码、算出初始成本找到全局最优进入主循环循环内更新每个鲸鱼个体用贪心规则接受新解然后重新评估全局最优迭代结束后输出最优位置对应的路线集合。为了验证结果可以用Matlab的plot画出所有路径仓库用五角星不同路线用不同颜色的线条客户用圆圈表示。我习惯把整个流程封装成一个WOA_OVRP_solver.m函数输入坐标、需求、容量、算法参数输出最终路径和成本。这样后续在项目里只要调用一行函数就行不用每次复制主循环。有的读者可能会问为什么不直接用Matlab全局优化工具箱里的ga或particleswarm原因有两个第一这些内置函数默认针对连续问题OVRP需要构造自定义的整数编码和适应度函数写起来并不比WOA短第二内置算法对约束的处理不如你手写解码器那样直观。在科研和项目中把算法完全掌控在自己手里方便后续改约束、改目标、加并行所以我更推荐手写WOA。另外在迭代过程中可以把每一代的最优成本record(t)保存下来最后绘制收敛曲线。收敛曲线是一种非常直观的诊断工具如果曲线没有“L型”下降而是长时间平行波动说明参数设置可能出了问题需要回到解码器或目标函数去排查。4. 实验结果、参数调优与常见问题4.1 测试算例与收敛性验证我先把代码跑在了一个随机生成的12客户小算例上。仓库坐标在(50,50)客户坐标在[0,100]区间随机需求量在[3,10]车辆容量Q30。对于12个客户可以用枚举法瞬间得到理论最优成本这里不具体展开了。WOA用30个种群、300次迭代跑20次每次都能收敛到接近最优最好一次达到与枚举一致平均偏差1.5%左右。收敛曲线呈现非常典型的“L型”前50代快速下降之后减速150代后基本水平。这说明WOA的探索能力在前中期比较足后期开发也不错。为了观察规模影响我又生成了一个50客户的算例Q50。由于没有精确最优对比就让算法跑2000次迭代种群60取多次运行的最小值作为“参照最优”然后用不同迭代次数去跑对比。结果显示300代时平均偏差约6%800代时降到2%以内1500代后基本不再改进。如果你只是快速筛选解可以设置800代如果是最终方案建议跑1200代以上。当然不同的随机数种子会影响结果所以工程上至少运行10次取最高质量解作为最终方案并尽量在参数上做小范围扰动测试。下面这个表格是我在某次实验中记录下来的参考数据客户数种群规模迭代次数平均运行时间偏差对照多次运行最优12303003.2秒1.5%304060012.8秒1.1%5060100048.5秒0.9%注意运行时间基于普通笔记本不同机器差异很大表格只提供相对参考。重点要看的是“偏差随时间下降”的趋势这说明算法在正确的方向上收敛。如果你的机器跑出的时间和这个表差很多不需要担心只要收敛曲线正常结果就基本可信。4.2 种群规模和迭代次数怎么配这部分是纯经验不具备数学最优性但可以省掉很多试错时间。第一种群规模和客户数不是线性增长关系。客户从20增加到100种群从30翻到60往往就够了因为WOA每个个体每次迭代都会朝全局最优靠近真正决定搜得多广的是随机搜索分支的触发频率。如果种群太大单次迭代耗时成倍增长但收敛速度未必更快。第二迭代次数和问题复杂度的关系更大。带容量约束的OVRP路线分段数越多搜索空间越大需要的迭代就越多。判断收敛是否完成的直接指标是收敛曲线是否拉平而不是拍脑袋定500次。我习惯在代码里加入“早停”逻辑连续50代最优成本不下降就break然后多跑几次取最优。第三位置边界不需要过于敏感但也不能太小。前面讲过边界影响每个维度数值的绝对大小而解码只关心排序。如果边界太小导致浮点精度把两个不同的键值区分不开排序会出现随机抖动影响收敛稳定性。我统一用[-10,10]。第四收敛因子a的衰减方式可以试试非线性的。有时候线性降到0会让算法后期失去探索性余弦衰减等非线性方式能保持一定活力。这不是标准WOA的范畴属于改进技巧。工程上如果你发现收敛曲线到了后期完全平死可能是a衰减太快可以改为余弦衰减例如a 1 cos(pi * t / max_iter)这样前期和后期都有一定的探索能力。不过要注意改了a的公式后A的取值范围也会变化最好重新测试几组算例确认效果。第五要不要每次都把所有个体完全替换我前面的代码用了贪心接受只接受更好的解。这是一种很有效的“保底”策略但也可能让早熟的个体长时间占据种群导致多样性下降。如果你发现最终结果不够好可以把接受准则改成模拟退火式的概率接受或者每隔一定代数强制重置一部分较差个体。这些都算工程优化具体怎么选取决于你对“稳定”和“探索”的偏好。4.3 典型踩坑记录与修复方法这里分享几个我在实际调试中反复遇到的坑每一个都值得仔细看。坑一目标函数没有去掉末尾回程距离。这是最隐蔽的错误。很多初版代码直接从闭环VRP改过来只改了模型描述没有改cost函数。结果算法输出一版非常漂亮的“闭环路线”但成本偏高和开放路径最优解对不上。检查方法很简单看最终路径图里是否每一条路线都闭合。如果每条路线的终点都莫名其妙地连回仓库说明目标函数里多加了最后一段。把cost_ovrp里对dist(route(end), depot)的计算删掉即可。坑二解码时排序索引与客户编号错位。位置向量长度为N客户数排序后的idx取值范围是1到N。如果你直接把idx作为客户编号而实际客户编号是2到N1就会少1号客户、多1号仓库。我建议所有索引都统一口径要么depot1、客户2:N1要么depotN1、客户1:N。写代码时在decode函数入口处加一行注释标明“idx depot_idx才是实际客户编号”避免过几天自己都忘了。坑三车辆数上限被忽略。如果业务规定了最多只能使用K辆车而解码器只管容量切分可能出现路线数超过K的解。这种情况需要在decode_ovrp中加一个计数器一旦发现新开路线数达到K且还有剩余客户就把该解的目标值设为Inf。注意把目标值设为Inf的前提是成本函数能正确处理Inf比如不写死循环。还要注意如果某个客户需求小于Q但K1且总需求大于Q那么所有解都是Inf算法将输出一个不可用结果——这其实是算例本身无解要在数据预处理时提前判断。坑四位置向量出现NaN导致排序错乱。大多数NaN来自更新时除以零或指数爆炸。WOA公式里涉及exp(b*l)如果l很大exp可能溢出为Inf虽然l实际限制在[-1,1]但如果你手动改了l的生成范围就有风险。另外边界裁剪如果写错了也可能把数组变成NaN。我的习惯是每次更新后检查new_pos若有NaN或Inf就把该个体重新随机初始化。这个小检查只会牺牲极少量计算时间却能避免莫名其妙的排名错误。坑五算法早熟所有个体挤在一起。WOA后期确实容易出现群体聚集。当种群多样性过低时随机搜索分支很难触发因为随机选取的个体和当前个体几乎一样。解决方法有三个一是增加随机搜索概率把p的阈值从0.5调高到0.6甚至0.7让更多个体走探索分支二是对最优个体做局部扰动给它乘以一个接近1的随机系数三是引入“重启”机制每30代将最差的10%个体重新随机初始化。这些技巧在标准论文里没有但在工程实践里非常有效。4.4 扩展到带时间窗和单边距离变体实际项目里OVRP很少以纯容量约束出现最常见的是带时间窗的OVRPTW。此时目标函数要增加提前等待成本和迟到惩罚。解码过程不变但在cost函数中需要额外计时车辆在仓库的出发时间可以设为0按顺序服务到达客户i的时间等于上一节点出发时间加行驶时间若早于最早允许开始时间则等待到最早时间若晚于最晚允许时间则罚M倍的迟到时长等待时间是否计入成本看业务要求。简单伪代码如下time 0; for i 1:length(route) time time travel_time(previous, current); if time earliest(current) waiting earliest(current) - time; total_cost total_cost waiting * w_wait; time earliest(current); elseif time latest(current) total_cost total_cost (time - latest(current)) * M; end end其中travel_time可以用距离除以平均速度也可以是真实路网矩阵。如果使用非对称距离矩阵比如道路存在单行道或者上下行用时不同只需把距离矩阵替换成非对称矩阵算法主体不用改。多车场OVRP则复杂一些需要在解码时把“车场选择”也映射进位置向量比如位置向量前M维作为车场选择概率后N维作为客户排序键这样车辆可以分别从不同车场出发开放路径的终点仍为客户点。异质车队不同容量可以在解码时依次尝试不同车辆的容量而不是统一用Q同样是小改动。这就是把算法核心做扎实的好处业务约束变更时主要工作量只在解码器和成本函数里不需要推倒重来。5. 写在最后的一点个人实操心得这套代码从零到落地我前后改过不下五个版本。最开始的版本完全照搬论文公式解的质量反而不如后来加了“贪心替换”和“早停”的版本。这让我意识到学术伪代码和工程代码之间永远隔着一层“防坑胶水”。特别是对于开放路径这种目标函数比标准VRP更简单的变体痛点更多在于如何让算法在离散解码和连续搜索之间顺滑切换。如果你在跑代码的时候发现解很差我建议先画出每一代的最优路径用眼睛看不要只看成本数字。有时候路径形状能立刻暴露问题比如所有路线都自动闭合了说明你最后一段回程没删干净比如两条路线有交叉且都不顺路说明编码映射出了问题。个人经验是把可视化放在调试第一步能让排错效率提高一倍。最后再分享一个小技巧在迭代后期对最优个体的位置向量做一次“局部微调”——给它乘以一个接近1的随机扰动系数然后重新解码看看是否更优。这个小操作经常能压榨出1%~3%的成本改进而且实现成本极低。希望这套讲解能帮你快速把WOA应用在OVRP上也期待你根据实际业务把它改造成更顺手的工具。