MATLAB求解LRP：节约算法与禁忌搜索的启发式实现

简介：基于MATLAB实现的本科毕业设计源码包，将节约算法与禁忌搜索算法用于LRP（定位-路径）问题的求解。面向物流工程、计算机、自动化等专业的在校学生与毕设开发者，既适合作为算法课程设计参考，也可用于科研前期的算法对比实验与算法改进验证。资源共16个文件，以12个.m格式的MATLAB源程序为核心，覆盖初始解构造、距离矩阵计算、节约值排序、禁忌搜索迭代等关键环节，另有3个.asv自动备份文件和1个README说明文档，便于梳理代码结构和运行流程。目前已有148人学习下载。源码全部测试通过，答辩评审平均分达96分，压缩包内包含可直接运行的测试脚本与演示程序，从数据载入、算法参数调整、优化迭代到结果输出构成完整闭环，便于在此基础上替换测试数据、修改变邻域策略或扩展其他启发式算法。

1. LRP 难在哪：为什么纯节约算法不够，还要叠禁忌搜索

把设施选址和车辆路径放在同一个模型里解，是物流优化里典型的“硬骨头”——设施开在哪、开几个、车从哪个设施出、沿着什么顺序送，四个决策互相耦合，单独拉出来任何一个都是 NP-hard。而这个毕设项目把两条启发式路线都做全了：先让节约算法（Clarke-Wright Savings）快速构造一个不差的解，再用禁忌搜索（Tabu Search）在这个解附近持续翻找更优结构。这个“先构造、再改进”的组合在工程上非常常见，因为你总不能拿小规模 MIP 硬解上百个客户。这套 MATLAB 代码适合两类人：一类是物流工程、工业工程相关专业准备毕设答辩的学生，另一类是想快速在 MATLAB 里落地启发式算法、拿真实规模数据跑实验的研究者。读完这篇，你会清楚每个.m文件在流水线里的位置，也知道调参时该动哪几个旋钮。

2. 节约算法的数学动机与 Gmcws.m 的 MATLAB 实现

2.1 从 depot 出发的合并逻辑

节约算法的出发点很朴素：如果一辆车单独给客户 i 跑一趟、再单独给客户 j 跑一趟，总里程是两条往返之和；如果让同一辆车从 depot 出发依次服务 i 和 j，就能把这两趟的公共路段省掉。省下来的那段就是节约值：

s(i, j) = d(depot, i) + d(depot, j) - d(i, j)

s(i, j) 越大，说明把 i 和 j 拼在同一条线路上越划算。算法把所有客户对的节约值从大到小排序，依次尝试合并两条线路，只要合并后不超车辆容量、不超线路长度约束，就执行合并。

在 LRP 场景里，这个公式要改一个地方：线路的起点不一定只有一个 depot，而是有多个候选设施。每个客户在初始化时会先被分配到某个设施下，后续合并也只是在同一个设施覆盖的客户集合内部进行。

% Gmcws.m 中节约值矩阵建立的常见写法 n = size(coord, 1); savings = zeros(n, n); for i = 2:n for j = i+1:n % dist(1, i) 表示设施(索引1)到客户i的距离 savings(i, j) = dist(1, i) + dist(1, j) - dist(i, j); savings(j, i) = savings(i, j); end end % 按节约值降序排列，得到一个待合并的候选序列 [~, idx] = sort(savings(:), 'descend');

这段代码把节约值写成一个对称矩阵，再拍平排序。注意dist(1, i)里的索引 1 是设施编号，实际项目中设施可能不止一个，那就要把目标客户先归类到某个设施下，再以该设施作为公式里的“0 点”。排序后的idx是后续合并循环的遍历顺序，这一步决定了算法是从最划算的合并开始试，而不是按客户编号顺序。

类 Weighted Savings 的做法更适用于 LRP：给第二个距离项乘上一个权重 λ，s(i,j) = d(0,i) + d(0,j) - λ·d(i,j)，λ 大于 1 时倾向合并地理位置近的客户，λ 小于 1 时倾向保住大的辐射路径。在项目里调整 λ 是改变解形态最直接的手段。

2.2 initial.m 与 Gmcws 的配合

初始解生成通常分两段。initial.m负责把一个粗糙的可行解搭出来：每个客户分配到最近且有剩余容量的设施，然后每个设施下挂一条或多条单车线路；随后Gmcws.m在这个基础上执行节约合并。注意不要跳层——直接拿纯随机解喂给节约算法，合并效果会被随机性冲淡。

这个工程里Gmcws的关键参数一般在文件头部集中定义，常见设计是：

pic.m里的图形会直接反映这一类参数的效果：lambda 偏大时，散点图上会看到明显簇状分布，因为算法更偏好局部抱团线路。

核心合并循环如下：

while ~isempty(cand) % 取出当前最大的节约值及其客户对 s_val = val_top; i = ci; j = cj; if isSameRoute(i, j) % 已在同一线路，跳过 cand = cand(2:end); continue; end if feasMerge(i, j, capacity) % 检查容量 mergeRoute(i, j); % 修改线路链表 end cand = cand(2:end); end

循环里每次只取当前最高节约值。isSameRoute判断两个客户是否已经在同一条物理线路上，feasMerge做容量复核。容易栽的坑是忘记在mergeRoute后更新设施剩余容量——因为一个客户不再占用原设施的资源，而是要由新线路对应设施服务，LRP 里这个数据不一致会让后面的countloc.m统计出错误结果。

2.3 节约算法给禁忌搜索留的起点

跑完 Gmcws，得到的解通常已经比纯最近邻好出一截，下一步是把它喂给tabusearch.m。在对接时要注意：保存解的数据结构必须统一。项目里常见是用结构体数组，每个元素包含facility、route、cost三段，tabusearch才能直接读取和修改。如果 saved solution 里只存了客户序列而丢了设施归属，禁忌搜索的邻域交换就无从谈起。

另外，节约算法是贪婪的，它给出的线路集合里常常存在单个客户被孤零零挂在线路末尾的情况，这正是禁忌搜索最值得改进的地方——通过跨线路搬迁把它塞进别的车次里。这段先记住结论：初始解的质量决定禁忌搜索的起点，但不决定最终质量；真正决定质量的是禁忌搜索能不能跳出局部最优。

3. 禁忌搜索的三件套：邻域、禁忌表、特赦规则

3.1 解编码与邻域算子

tabusearch.m负责在两个层面同时搜索：路径层面和选址层面。路径层面用的是经典算子——2-opt 反转一段客户序列、单点搬迁把一个客户从 A 线路移到 B 线路、交换算子互换两个客户的位置。选址层面则是翻转某个候选设施的开关状态：一个设施从“关闭”变成“开放”，原来挂在别的设施下的客户可以搬过来，以节省长途运输。

编码用两段式：设施向量facility（每个客户属于哪个设施）+ 客户序列向量route（每条线路的访问顺序）。这样设计的好处是邻域算子在两个段上操作不会相互污染，算成本和校验容量也容易拆开。

3.2 tabusearch.m 的迭代骨架

% tabusearch.m 的核心循环 curr = initSol; best = initSol; bestCost = totaldis(best); tabuList = zeros(n, n); % 禁忌表：记录边(i,j)的剩余禁忌代数 iter = 0; while iter < MaxIter nbr = generateNeighborhood(curr); % 生成候选邻域 [nbr, move] = selectBest(nbr, tabuList); % 挑未禁忌中最好的 if totaldis(nbr) < bestCost - 1e-6 best = nbr; bestCost = totaldis(best); end % 禁忌当前做过的边交换，防止回跳 tabuList(move.i, move.j) = TabuLen; tabuList = max(tabuList - 1, 0); % 每轮衰减 curr = nbr; iter = iter + 1; end

totldis在这里是总行驶距离或总成本计算函数（工程里对应totaldis.m）。generateNeighborhood一般不会枚举全邻域——n 到 200 时全枚举会直接卡死，常见做法是随机抽 30～60 个候选移动，挑其中目标函数最小的那个。tabuList存的是边对 (i, j)，表示最近若干代禁止把客户 i 和 j 重新接到一起；这里用矩阵存储、每轮全体减 1 的实现简单直接，缺点是内存开销随 n 平方增长，n 超过 500 时就要改用稀疏哈希表。

3.3 禁忌长度与迭代次数怎么定

禁忌长度是这套算法里最敏感的参数，没有之一。它太小，搜索会绕回刚走过的地方；它太大，候选解全被锁死，算法退化成局部爬山。经验规律是：禁忌长度设成客户数量的 1/5～1/3，并且在迭代过程中动态波动（比如每 20 代在 base 和 2×base 之间随机变化），效果比固定长度稳定得多。

MaxIter的设定要和禁忌长度成比例。如果发现执行到最后几十代bestCost还在下降，说明迭代次数不够；如果前 1/3 就到瓶颈、后面全是无效搜索，优先缩减迭代次数并扩大邻域抽样数，而不是继续加迭代。

3.4 特赦规则与解改进记录

特赦规则用于处理一种尴尬情况：某个候选移动虽然被禁忌，但它带来的目标值比当前全局最优还好。这时的标准做法是忽略禁忌直接采纳，因为走到这一步说明原来的路线已经被充分探索，这次移动可能通向真正的新区域。实现上只需要在selectBest里加一条判断：候选解的 cost 小于 bestCost 时绕过 tabuList。

一个容易被忽略的细节是记录历史最优解的序列，而不仅仅是 cost 值。项目里的best必须保存完整的route和facility字段，否则收敛后发现结果成本下降，却拿不出对应路径来绘图验证，问题排查会变得很困难。

4. 从 initial 到 cplex_function：数据流与精确解对照

4.1 文件清单与调用关系

这套代码里每个.m文件的职责边界比较清晰，跑通整个流程需要厘清它们之间的依赖关系。

Gmcws.asv、demo.asv、tabusearch.asv是 MATLAB 自动保存的备份文件，内容比同名.m旧一个版本。调试时如果.m被改坏了，可以用这些.asv找回上一个能跑的版本，这也是一个免费的“后悔药”。

4.2 demo.m 的完整调用序列

% demo.m 主流程（结构说明） coord = load('coords.txt'); % 每行: id x y demand = load('demands.txt'); % 每行: id qty facCap = load('facility_cap.txt'); % 候选设施的容量 dist = disr(coord); % 用坐标生成距离矩阵 initSol = initial(coord, demand, dist, facCap); savingSol = Gmcws(initSol, demand, dist, facCap, 1.1); finalSol = tabusearch(savingSol, demand, dist, facCap, ... 'MaxIter', 400, 'TabuLen', 25); pic(finalSol, coord); % 画图 fprintf('Total cost: %.2f\n', totaldis(finalSol));

第 4 行的disr只是坐标转矩阵的纯计算——注意两点间距离要取欧氏距离，坐标单位不一致的结果会直接偏掉。第 6 行的Gmcws多传了一个 1.1，那就是加权节约参数 λ。整条链路是串行的：初始解 → 节约优化 → 禁忌搜索 → 可视化。在tabusearch内部，会用findis在候选设施里刷新当前最优设施集合，再在集合上做路径邻域搜索。

4.3 cplex_function.m 的 MIP 建模骨架

精确解对照是这份代码里很有价值的一块。cplex_function.m做的是把 LRP 写成混合整数规划，交给 CPLEX 求解。目标函数包含三块：设施开设固定成本、从设施到客户的运输成本、客户点之间的路径连接成本。

% cplex_function.m 的建模骨架 prob = struct; prob.Name = 'LRP'; % x(i,j,k): 车辆k从设施i到客户j的弧; y(i): 设施i是否开放 prob.vtype = [repmat('B', 1, nFacility), repmat('B', 1, nArcs)]; prob.obj = [facOpenCost, arcCost]; % 目标函数系数 prob.Aineq = [...]; % 容量约束、度约束 prob.bineq = [...]; prob.sense = 'minimize'; sol = cplexmilp(prob);

这组约束里最关键的是“每个客户恰好被服务一次”和“线路起点必须对应一个开放设施”这两条。后者是把选址和路径绑定起来的核心约束：如果没有这条，CPLEX 会给出一个“从关闭设施出发送所有客户”的奇葩解，这在数据层面合法、在实际业务里完全没意义。

demo.m可以加一个对比实验：取 10 个候选设施、20 个客户的小样本，分别跑cplex_function和tabusearch，输出两者的成本差。常见结果是 CPLEX 给出精确最优，禁忌搜索给出差距在 3%～8% 的解，但耗时差两个数量级。这个数字写进毕设论文里，说服力远比“我跑了两个算法”强。

现在像 Codex 这类编码助手已经能直接操作 MATLAB 任务，把这段 MIP 骨架转成实际的 CPLEX 接口调用并不难，前提是你自己先清楚数据流和约束语义——结构体字段名写错，它只能报语法错误，纠不了模型错误。

5. 参数整定、pic 画图和三类典型运行错误

5.1 tabuLen 与 MaxIter 的搭配整定

先把禁忌长度和迭代次数的关系想清楚：MaxIter是搜索的时间上限，TabuLen是记住历史移动的时间窗口。两者的合理配比是 TabuLen 占 MaxIter 的 5%～10%。比如 MaxIter=500，TabuLen 取 25～50；设到 80 以上，搜索很容易守着一个区域原地打转。判定方法是打开totaldis的收敛曲线：如果曲线像锯齿一样抖个不停，优先加大 TabuLen；如果曲线是平滑下降但后半段完全不降，说明邻域算子太弱，优先增加抽样数。

countloc.m输出的设施数量和totaldis有很强的负相关性——设施开得越多，运输距离越短，但设施固定成本越高。调参时盯着总成本看，别只看运输距离。

5.2 运行前先在内存里验约束

在跑完整流程之前，用下面这段检查脚本确认初始解和搜索结果都没有破约束：

% 验证解是否满足容量约束 for s = finalSol loadSum = sum(demand(s.route)); assert(loadSum <= facCap(s.facility), ... '容量越界: 设施%d 超载 %.1f', s.facility, loadSum); end

这类断言函数放在ceshi.m里很合适。如果断言被触发，错误信息里会直接给出是哪个设施超载多少，省去逐条线路排查的功夫。

5.3 三个容易踩的运行时错误

第一类，距离矩阵不对称。disr.m返回的矩阵偶尔因为浮点舍入或数据读入顺序出现微小不对称，导致最终解成本偏差。跑完第一时间验证max(abs(dist - dist'))，超过 1e-9 就该原地对称化。第二类，initial.m生成的解本身不可行——未读入设施容量就开始分配客户。这类问题通常在countloc输出异常时露出马脚。第三类，禁忌表长度不符合规模，小样本上用大禁忌表导致搜索过早停滞，表现为bestCost在前 50 代就锁定。遇到这种情况别急着删代码，先打印每 20 代的 cost 变化率，通常能明显看出是哪一类故障。

最后留一个操作清单：跑通demo.m之后，把pic.m画出的图导出成矢量格式，再在图上标注 open facility 的编号和每条线路的总需求。这样既完成了结果验证，也直接有了论文里的成品插图。

本文还有配套的精品资源，点击获取