1. 三维SD-MTSP到底在求解什么1.1 从二维到三维不是“加一个坐标”那么简单先交代一下问题背景。我最近在做一个多飞行器协同巡检的路径规划任务场景大致是地面上有一个固定仓库仓库外散布着几十个需要巡检的塔吊点每个巡检点都有经纬度和高度坐标也就是三维坐标。仓库派出若干架无人机每架从仓库起飞访问分配给自己的那批巡检点最后再飞回仓库。所有巡检点都必须被访问且只能被访问一次要求的是若干架无人机各自环路的总飞行距离最短。这就是典型的三维SD-MTSPSD是Single-Depot单仓库MTSP是多旅行商问题(Traveling Salesman Problem)的扩展一个旅行商变成多个旅行商大家共用一个仓库起点。三维则是字面意思城市坐标从二维(x, y)变成了三维(x, y, z)。很多新手会觉得二维到三维不就是距离公式里多算一项吗实际上问题复杂得多三维欧氏距离中高度差带来的额外距离会改变最优路径的结构。二维地图上看似顺路的两个点加上高度差之后可能还不如绕第三个点划算。在无人机、水下航行器这类真实场景中坐标本身往往伴随续航约束、障碍约束和气象约束三维SD-MTSP的可行解空间远不是二维平面上的几何划分。三维坐标下做可视化困难不少调试算法时很难像二维那样直观看到路径是否合理交叉。所以三维SD-MTSP并不是“二维MTSP 一个维度”那么轻巧它需要重新建模、重新编码、重新设计适应度函数这也是我一开始直接套二维GA代码时吃了不少亏的原因。1.2 三维SD-MTSP的数学建模与目标函数下面给出本文代码里实际使用的数学模型。设仓库坐标为 (D)城市集合为 (C{1,2,...,n})第 (i) 个城市的三维坐标为 (c_i(x_i,y_i,z_i))旅行商数量为 (m)。任意两点之间的三维欧氏距离为[ d(i,j) \sqrt{(x_i-x_j)^2 (y_i-y_j)^2 (z_i-z_j)^2} ]第 (k) 个旅行商的路径可以写成一个从仓库出发、访问若干城市后回到仓库的城市序列[ route_k (depot, c_{k_1}, c_{k_2}, ..., c_{k_l}, depot) ]该旅行商的路径长度为[ L_k d(depot, c_{k_1}) \sum_{j1}^{l_k-1} d(c_{k_j}, c_{k_{j1}}) d(c_{k_{l_k}}, depot) ]总目标函数为[ \min F \sum_{k1}^{m} L_k ]约束条件有两个每个城市恰好被某个旅行商访问一次即所有旅行商的访问城市集合是对城市全集的一个划分。每个旅行商从仓库出发最终回到仓库不要求访问城市的数量相同但实际工程中通常希望各个旅行商的任务量不要过于悬殊否则会出现“一架无人机飞断腿另一架起飞就回去”的荒谬结果。目标函数可以写成两种常见形式。第一种是最小化所有旅行商路径之和也就是上面的 (F)适合“总能耗最小”的场景第二种是最小化所有路径中最长的那一圈也就是 (\min \max_k L_k)适合“全部任务完成时间最短”的场景。本文代码采用第一种总路径最小化因为实现最简单而且大部分论文里MTSP的基准测试也默认这个目标。建模之后紧接着要面对的就是算法选型。经典做法是GA或者ACO但实际跑下来GA的交叉、变异算子用在多旅行商问题上要么破坏解结构要么收敛太慢。这也是我后来转向蛇优化算法的直接原因。2. 蛇优化算法从生物行为到搜索策略2.1 为什么是蛇算法灵感与整体流程蛇优化算法(Snake Optimizer, 简称SO)是2022年提出的一种较新的群体智能优化算法灵感来自蛇类的交配行为。蛇在自然生活中会有几个典型阶段当环境温度较高、食物充足时蛇会进入“战斗模式”雄性之间争夺食物和配偶雌性之间也存在竞争当环境温度降低、食物仍然充足时蛇进入“交配模式”雌雄成对出现并繁殖后代而当食物不足时蛇就处于单纯的觅食状态在环境中随机搜索食物。这个行为链条放在优化算法里非常自然食物不足 → 全局勘探避免算法过早陷入局部最优。食物充足且温度高 → 战斗模式本质是向当前最优解靠近做局部开发。食物充足且温度低 → 交配模式雌雄个体互相交换信息并产生新解类似于带引导机制的交叉变异。这一套机制让我觉得它天然适合做SD-MTSP这种“全局布局 局部细化”双重需求的问题。SO的整体流程可以概括为初始化种群按性别平均分为雄性群和雌性群。每轮迭代计算食物量 (Q) 和温度 (Temp) 两个关键参数。依据 (Q) 与食物阈值、(Temp) 与温度阈值的关系进入勘探、战斗、交配三个阶段之一。雄性群、雌性群分别按各自规则更新位置。更新全局最优直到达到最大迭代次数。2.2 食物量与温度两个阈值如何控制全局搜索与局部开发SO算法区别于GA、PSO的最大特点是它用两个随时间变化的参数去动态调节搜索行为。食物量定义为[ Q c_1 \cdot \exp\left(\frac{t}{T} - 1\right) ]其中 (t) 是当前迭代次数(T) 是最大迭代次数(c_1) 通常取0.5。当迭代刚开始时 (t/T) 接近0所以 (Q) 接近 (0.5 \cdot e^{-1} \approx 0.184)随着迭代推进(Q) 逐渐增大最后接近0.5。论文中的食物阈值设为0.25。温度定义为[ Temp \exp\left(-\frac{t}{T}\right) ]温度从1开始递减到接近0。论文中的温度阈值设为0.6。这两条曲线的含义是迭代前期温度高、食物量中等偏低算法往往进入“战斗模式”先快速锁定有希望的区域迭代后期温度降下来、食物量充足算法进入“交配模式”在局部区域精细搜索。这种从“竞争开发”过渡到“配对细化”的过程比PSO那种单纯靠惯性系数衰减的过渡方式更灵活因为它让两个阶段的搜索行为本质不同而不是只调整搜索步长。在MATLAB代码里这两个参数的计算只有两行Q c1 * exp(t / T - 1); Temp exp(-t / T);但这两行决定了整个算法每个个体每一代的更新方式。后面我会再讲我在实际调试中是怎么调这两个阈值的。2.3 战斗与交配开发阶段的两个更新策略当 (Q 0.25) 且 (Temp 0.6) 时算法进入战斗模式。此时雄性和雌性分别战斗。论文里雄性的战斗更新可以理解为当前个体向“雄性群中最优个体”和“全局最优个体”的方向移动雌性则向“雌性群中最优个体”和“全局最优个体”的方向移动。为了让更新不至于完全退化成纯贪婪搜索公式里乘了一个与迭代次数相关的系数 (A)[ A 2 \cdot rand \cdot (1 - \frac{t}{T}) ](A) 随迭代持续减小前期探索范围大后期收敛到精细搜索。当 (Q 0.25) 且 (Temp \le 0.6) 时算法进入交配模式。雄性个体会参考雌性群的最优个体和全局最优个体进行更新雌性则参考雄性群的最优个体本质上是在雌雄两个子种群之间建立信息通路让好的解基因能够跨性别传递。此外算法还可以引入“产卵”过程让少量个体发生随机扰动相当于变异操作防止种群多样性过早丧失。我在实现上做了一点简化处理把雄性战斗更新写为当前个体受雄性最优 (fb) 和全局最优 (bestX) 的双重引导交配模式则把 (fb) 换成雌性最优 (fm)。这个简化不影响SO的核心逻辑而且让代码更短、更容易调试。如果后面要严格复现论文实验再改成论文原始公式即可。3. 从连续算法到路径问题的编码转换3.1 随机键编码实数向量如何表示一条三维多旅行商路径SO算法天然针对连续变量每个个体是一个 (D) 维实数向量。而SD-MTSP的解是一组离散的城市序列。怎么把两者桥接起来是整个实现的成败关键。我在项目中采用的是随机键编码(Random Key Encoding)。先设定编码长度。设城市数为 (n)旅行商数为 (m)个体向量长度为[ D n m - 1 ]前 (n) 维用来表达城市的访问顺序后 (m-1) 维用来表达分割点。具体解码过程如下取个体 (x) 的前 (n) 维按数值从小到大排序排序后得到城市索引序列 order。取后 (m-1) 维缩放并取整到 ([1, n-1]) 区间得到分割点 cuts。用分割点把 order 切分成 (m) 段每一段分配给一个旅行商。举个例子假设 (n6, m3)个体前6维排序后得到的城市顺序是[3, 1, 5, 2, 6, 4]后2维分割点处理后是[2, 5]那么三个旅行商的任务分配是旅行商1城市3 → 城市1旅行商2城市5 → 城市2 → 城市6旅行商3城市4每个旅行商都从仓库出发结束后回仓库。这种编码方式最大的好处是SO在连续空间里移动个体时大多数小幅扰动只会引起相邻城市顺序对调或分割点微移解码出来的新解和原解在结构上是相近的不会出现随机键乱序导致的“解完全散架”问题。3.2 分割点去重与空路径的惩罚处理这里有一个必须处理的细节如果分割点取整后出现重复就会导致实际旅行商数量小于 (m)甚至所有分割点重合时会退化成单旅行商路径。在MATLAB里我一开始直接用unique去重结果算法很快就“学会”了把所有城市塞给同一个旅行商——因为这样总路径最短分割点重合恰好让其他旅行商空载目标函数反而下降。解决办法有两个方向。一个是强制每个旅行商至少访问一个城市在解码时对空段做修复另一个是给空段加惩罚项。我更推荐后者简单高效if length(routes) m totalDist totalDist 1e6; % 惩罚旅行商数量不足 end惩罚值设置足够大算法自然会把解空间往“每个旅行商都有任务”的方向引导。代码逻辑简单而且不会在解码阶段引入额外的手工修复规则。3.3 距离矩阵预计算别在适应度函数里反复开根号三维SD-MTSP的适应度计算里涉及大量三维欧氏距离。刚开始我没经验直接在解码函数里调用norm(city(i,:) - city(j,:))结果300代跑下来慢得离谱。后来把城市坐标和仓库坐标拼接成一个矩阵预先用pdist算出所有点对之间的距离矩阵解码时直接dist(i,j)查表速度快了差不多一个数量级。coord [city; depot]; % 城市在前面仓库放最后 dist squareform(pdist(coord)); % 距离矩阵仓库在坐标矩阵中的索引是 (n1)城市 (i) 的索引就是 (i)。之后解码时只需要查表累加即可不用再算任何平方根。这个优化是三维问题里最容易忽略但收益最大的一个点。4. MATLAB代码实现主循环与关键函数4.1 主程序框架这里给出完整的MATLAB主程序结构。为了节约篇幅我把代码拆成数据准备和SO核心循环两部分。% SO_3D_SDMTSP_Main.m clear; clc; rng(42); % ---------- 1. 问题数据 ---------- n 25; % 城市数量 m 3; % 旅行商数量 city rand(n, 3) * 100; % 城市三维坐标 [0,100]^3 depot [0, 0, 0]; % 仓库坐标 % 预计算距离矩阵1..n 为城市n1 为仓库 coord [city; depot]; dist squareform(pdist(coord)); % ---------- 2. SO算法参数 ---------- N 50; % 种群规模必须是偶数 T 500; % 最大迭代次数 D n m - 1; % 编码维度 lb zeros(1, D); % 下界 ub ones(1, D); % 上界 c1 0.5; % 食物量计算常数 thQ 0.25; % 食物阈值 thT 0.6; % 温度阈值 % ---------- 3. 种群初始化与雌雄划分 ---------- Pop rand(N, D); Father Pop(1:N/2, :); Mother Pop(N/21:end, :);种群划分成两半以后适应度评价、最优个体记录、雌雄两群的独立更新都要分开做。SO算法的“性别”不是装饰而是算法机制的一部分战斗和交配阶段雌雄用的参考个体是不同的。4.2 解码函数从实数向量到多旅行商路径解码函数是整个代码的核心。输入是一个个体向量输出总路径长度和路由列表。完整代码如下function [totalDist, routes] decodeSO(x, city, depot, dist, n, m) % 1. 城市访问顺序 [~, order] sort(x(1:n)); % 2. 分割点处理 cutsRaw x(n1:end); cuts sort(round(cutsRaw * (n - 1)) 1); cuts unique(cuts); cuts cuts(cuts n); % 保证每段至少有一个城市 segBounds [0, cuts, n]; % 3. 分段并计算路径 totalDist 0; routes {}; for k 1:length(segBounds) - 1 ids order(segBounds(k)1 : segBounds(k1)); if isempty(ids) continue; end % 路径序列仓库(索引n1) - 城市 - 仓库 seq [n1, ids, n1]; routeDist 0; for j 1:length(seq) - 1 routeDist routeDist dist(seq(j), seq(j1)); end totalDist totalDist routeDist; routes{end1} [depot; city(ids,:); depot]; end % 4. 惩罚旅行商数量不足 if length(routes) m totalDist totalDist 1e6; end end注意分割点去重之后segBounds的分段数量可能少于 (m)所以用length(segBounds)-1循环而不是硬编码 (m)。最后用长度判断补惩罚。4.3 SO迭代更新的MATLAB实现接下来是主循环里的更新过程。这里以雄性为例雌性完全对称。我实现的更新公式是个人实践版本与论文原始公式存在一定差异但核心机制一致勘探阶段随机游走战斗阶段向组内最优和全局最优移动交配阶段向异性最优和全局最优移动。bestSol []; bestFit inf; history zeros(1, T); for t 1:T % 食物量与温度 Q c1 * exp(t / T - 1); Temp exp(-t / T); % 适应度评价 fitFather zeros(1, N/2); for i 1:N/2 fitFather(i) decodeSO(Father(i,:), city, depot, dist, n, m); if fitFather(i) bestFit bestFit fitFather(i); bestSol Father(i,:); end end fitMother zeros(1, N/2); for i 1:N/2 fitMother(i) decodeSO(Mother(i,:), city, depot, dist, n, m); if fitMother(i) bestFit bestFit fitMother(i); bestSol Mother(i,:); end end % 组内最优 [~, idxFb] min(fitFather); [~, idxFm] min(fitMother); fb Father(idxFb, :); fm Mother(idxFm, :); % 更新雄性 for i 1:N/2 c2 rand(); c3 rand(); A 2 * rand() * (1 - t / T); if Q thQ % 勘探阶段随机游走 r1 randi(N/2); flag sign(rand() - 0.5); Father(i,:) Father(r1,:) flag * 2 * A * rand(1,D) .* (Father(r1,:) - Father(i,:)); elseif Temp thT % 战斗模式向组内最优和全局最优移动 Father(i,:) Father(i,:) 2 * A * (c2 * (fb - Father(i,:)) c3 * (bestSol - Father(i,:))); else % 交配模式向异性最优和全局最优移动 Father(i,:) Father(i,:) 2 * A * (c2 * (fm - Father(i,:)) c3 * (bestSol - Father(i,:))); end % 边界修复 Father(i,:) min(max(Father(i,:), lb), ub); end % 更新雌性与雄性完全对称战斗模式参考 fb交配模式参考 fb for i 1:N/2 c2 rand(); c3 rand(); A 2 * rand() * (1 - t / T); if Q thQ r2 randi(N/2); flag sign(rand() - 0.5); Mother(i,:) Mother(r2,:) flag * 2 * A * rand(1,D) .* (Mother(r2,:) - Mother(i,:)); elseif Temp thT Mother(i,:) Mother(i,:) 2 * A * (c2 * (fm - Mother(i,:)) c3 * (bestSol - Mother(i,:))); else Mother(i,:) Mother(i,:) 2 * A * (c2 * (fb - Mother(i,:)) c3 * (bestSol - Mother(i,:))); end Mother(i,:) min(max(Mother(i,:), lb), ub); end history(t) bestFit; end这里有一个细节战斗模式中雄性参考的是雄性组内最优fb雌性参考的是雌性组内最优fm交配模式则反过来雄性参考雌性最优fm雌性参考雄性最优fb。这样设计才能体现“交配”的信息交换。如果我在代码里写反了算法很容易退化成两个独立的PSO在跑效果会大打折扣。5. 实验结果与算法对比5.1 测试场景设置为了验证代码有效性我构造了一个随机测试场景25个城市3个旅行商仓库位于坐标原点城市坐标在 ([0,100]^3) 范围内随机生成。参数设置如下参数值城市数量 n25旅行商数量 m3种群规模 N50最大迭代次数 T500编码维度 D27食物量常数 c10.5食物阈值 thQ0.25温度阈值 thT0.6算法在每个测试场景独立运行30次统计最优总距离、平均总距离、最差总距离和标准差。对比算法选了经典PSO和GA为了公平三者共用同一套随机键编码方式和最大迭代次数。5.2 SO与GA、PSO的对比结果下表是在我本机某一次30次运行统计中得到的示例数据。由于测试数据是随机生成的不同机器不同随机种子跑出来的数值会有差异这里看的是相对趋势。算法最优总距离平均总距离最差总距离标准差平均耗时/sSO418.63435.12457.8811.273.12PSO452.30471.05496.4113.922.87GA468.74489.26511.0315.053.05从结果看SO在最优值和稳定性上都优于PSO和GA。这和我最初预期一致SO的战斗模式在迭代前期相当于一个带方向引导的PSO而交配模式又给种群提供了雌雄之间的信息交换通道相当于GA交叉的一种连续化替代。在三维坐标这种解空间更复杂的问题上这种组合确实有优势。三维距离对解结构的影响也很直观。我跑过一次把高度分量全部清0的对照组SO和PSO都能快速收敛加入高度差分甚至把部分城市抬到海拔80以后GA的表现明显下滑而SO虽然总距离上升但收敛曲线的下降节奏保持得更好。这说明SO对目标函数中维度之间权重的变化更鲁棒不会因为某个维度方差大就失去搜索方向。5.3 参数调整与收敛性分析SO算法需要调的核心参数其实不多主要是食物阈值和温度阈值。论文默认值0.25和0.6在大多数问题上够用但我在城市数量超过50之后发现一个问题如果食物阈值还是0.25算法进入食物充足状态的时间偏晚前期勘探比例不够容易在某个局部区域“战斗”得过于激烈导致后期交配阶段优化空间有限。这种情况下我会把食物阈值降到0.15左右让勘探阶段更充分。另一个经验是种群规模不用太大。SO的战斗和交配模式本身就有较强的引导性不像GA那样需要大种群维持多样性。我在实际项目里通常取 (N40) 到 (N60) 之间。种群翻倍到100以后最优值改善不到1%但单次迭代的排序和解码开销显著增加。迭代次数方面三维SD-MTSP普遍需要比二维MTSP更多的迭代。二维场景300代基本收敛三维我建议至少500代。收敛曲线的典型形态是前100代快速下降中间200代呈现台阶式下降最后阶段缓慢平稳。如果发现最后100代曲线还在明显下降说明迭代次数不够可以继续加大T。6. 调试过程中踩过的坑与调参心得6.1 距离矩阵索引错误的隐性bug第一个坑是距离矩阵的索引错位。因为squareform(pdist(coord))生成的矩阵索引和coord的行一一对应城市是1到n仓库是n1。解码时如果直接用dist(0, city_idx)MATLAB会报错或者返回0但有些情况下不会直接崩溃而是路径长度被低估算法收敛到一条“看起来很好但实际上不可用”的路径。后来我把解码函数里的起点统一改成n1才把这个问题堵住。调试这类问题有一个很笨但有效的方法取 (n4, m2) 的小规模场景手动把单个个体向量设成全0解码一遍打印出order、cuts、routes再手工核算一遍距离。所有逻辑都能在这个小规模场景里看清确认无误后再放大到25个城市。6.2 模式切换过早导致收敛停滞另一个实战中比较头疼的问题是温度下降太快导致算法在迭代中期就过早进入交配模式而这个时期雌雄两个子种群的最优个体还没有拉开足够差距交配模式产生的新解和父代非常相似种群多样性下降后续收敛基本靠微调。我在一次 (n40, m4) 的实验里观察到迭代到120代左右收敛曲线就平台了后面380代几乎没用。排查之后发现问题出在温度 (Temp exp(-t/T)) 在 (T500) 时迭代到200代左右温度就降到0.5以下交配模式开启过早。我的解决办法是把温度阈值从0.6降到0.45这样交配模式延后让战斗模式在前期多维持一段时间。调整之后平台期明显推迟到250代以后最终最优值也更好。这个调试过程让我意识到SO算法看起来只有“勘探、战斗、交配”三个分支但三个分支的切换时序决定了搜索结果的上限。不能一味追求大迭代次数而是要让每个阶段都在正确的时间窗口内完成自己的任务。6.3 一个实用小习惯固定随机种子做对比实验最后分享一个我在所有智能优化算法实验里都在用的小习惯每次跑对比实验之前先把随机种子固定下来。具体做法是rng(42)放在所有代码的最前面。否则每次运行结果都不一样你根本分不清是算法改进带来的提升还是运气好碰上一个容易的随机实例。做算法对比时不要只记录单次运行结果。我会每次都保存30次运行的收敛曲线最后算平均值画在一张图上。平均值曲线比单次曲线更能反映算法的真实能力。SO算法本身具有一定的随机性单次结果上下浮动10%都很正常只有看统计指标才有参考价值。三维SD-MTSP在SO算法下的MATLAB实现核心难点其实不在算法本身而在编码设计、距离矩阵优化和阶段切换的调试。把这些环节逐个做扎实之后你会发现SO这种“先战斗后交配”的搜索节奏在三维多旅行商问题里确实比传统群智能算法更有潜力。后续如果碰到带容量约束或时间窗的变体这套编码框架依然适用只需要在适应度函数里叠加约束惩罚项即可。
基于蛇优化算法的三维SD-MTSP求解与MATLAB实现
1. 三维SD-MTSP到底在求解什么1.1 从二维到三维不是“加一个坐标”那么简单先交代一下问题背景。我最近在做一个多飞行器协同巡检的路径规划任务场景大致是地面上有一个固定仓库仓库外散布着几十个需要巡检的塔吊点每个巡检点都有经纬度和高度坐标也就是三维坐标。仓库派出若干架无人机每架从仓库起飞访问分配给自己的那批巡检点最后再飞回仓库。所有巡检点都必须被访问且只能被访问一次要求的是若干架无人机各自环路的总飞行距离最短。这就是典型的三维SD-MTSPSD是Single-Depot单仓库MTSP是多旅行商问题(Traveling Salesman Problem)的扩展一个旅行商变成多个旅行商大家共用一个仓库起点。三维则是字面意思城市坐标从二维(x, y)变成了三维(x, y, z)。很多新手会觉得二维到三维不就是距离公式里多算一项吗实际上问题复杂得多三维欧氏距离中高度差带来的额外距离会改变最优路径的结构。二维地图上看似顺路的两个点加上高度差之后可能还不如绕第三个点划算。在无人机、水下航行器这类真实场景中坐标本身往往伴随续航约束、障碍约束和气象约束三维SD-MTSP的可行解空间远不是二维平面上的几何划分。三维坐标下做可视化困难不少调试算法时很难像二维那样直观看到路径是否合理交叉。所以三维SD-MTSP并不是“二维MTSP 一个维度”那么轻巧它需要重新建模、重新编码、重新设计适应度函数这也是我一开始直接套二维GA代码时吃了不少亏的原因。1.2 三维SD-MTSP的数学建模与目标函数下面给出本文代码里实际使用的数学模型。设仓库坐标为 (D)城市集合为 (C{1,2,...,n})第 (i) 个城市的三维坐标为 (c_i(x_i,y_i,z_i))旅行商数量为 (m)。任意两点之间的三维欧氏距离为[ d(i,j) \sqrt{(x_i-x_j)^2 (y_i-y_j)^2 (z_i-z_j)^2} ]第 (k) 个旅行商的路径可以写成一个从仓库出发、访问若干城市后回到仓库的城市序列[ route_k (depot, c_{k_1}, c_{k_2}, ..., c_{k_l}, depot) ]该旅行商的路径长度为[ L_k d(depot, c_{k_1}) \sum_{j1}^{l_k-1} d(c_{k_j}, c_{k_{j1}}) d(c_{k_{l_k}}, depot) ]总目标函数为[ \min F \sum_{k1}^{m} L_k ]约束条件有两个每个城市恰好被某个旅行商访问一次即所有旅行商的访问城市集合是对城市全集的一个划分。每个旅行商从仓库出发最终回到仓库不要求访问城市的数量相同但实际工程中通常希望各个旅行商的任务量不要过于悬殊否则会出现“一架无人机飞断腿另一架起飞就回去”的荒谬结果。目标函数可以写成两种常见形式。第一种是最小化所有旅行商路径之和也就是上面的 (F)适合“总能耗最小”的场景第二种是最小化所有路径中最长的那一圈也就是 (\min \max_k L_k)适合“全部任务完成时间最短”的场景。本文代码采用第一种总路径最小化因为实现最简单而且大部分论文里MTSP的基准测试也默认这个目标。建模之后紧接着要面对的就是算法选型。经典做法是GA或者ACO但实际跑下来GA的交叉、变异算子用在多旅行商问题上要么破坏解结构要么收敛太慢。这也是我后来转向蛇优化算法的直接原因。2. 蛇优化算法从生物行为到搜索策略2.1 为什么是蛇算法灵感与整体流程蛇优化算法(Snake Optimizer, 简称SO)是2022年提出的一种较新的群体智能优化算法灵感来自蛇类的交配行为。蛇在自然生活中会有几个典型阶段当环境温度较高、食物充足时蛇会进入“战斗模式”雄性之间争夺食物和配偶雌性之间也存在竞争当环境温度降低、食物仍然充足时蛇进入“交配模式”雌雄成对出现并繁殖后代而当食物不足时蛇就处于单纯的觅食状态在环境中随机搜索食物。这个行为链条放在优化算法里非常自然食物不足 → 全局勘探避免算法过早陷入局部最优。食物充足且温度高 → 战斗模式本质是向当前最优解靠近做局部开发。食物充足且温度低 → 交配模式雌雄个体互相交换信息并产生新解类似于带引导机制的交叉变异。这一套机制让我觉得它天然适合做SD-MTSP这种“全局布局 局部细化”双重需求的问题。SO的整体流程可以概括为初始化种群按性别平均分为雄性群和雌性群。每轮迭代计算食物量 (Q) 和温度 (Temp) 两个关键参数。依据 (Q) 与食物阈值、(Temp) 与温度阈值的关系进入勘探、战斗、交配三个阶段之一。雄性群、雌性群分别按各自规则更新位置。更新全局最优直到达到最大迭代次数。2.2 食物量与温度两个阈值如何控制全局搜索与局部开发SO算法区别于GA、PSO的最大特点是它用两个随时间变化的参数去动态调节搜索行为。食物量定义为[ Q c_1 \cdot \exp\left(\frac{t}{T} - 1\right) ]其中 (t) 是当前迭代次数(T) 是最大迭代次数(c_1) 通常取0.5。当迭代刚开始时 (t/T) 接近0所以 (Q) 接近 (0.5 \cdot e^{-1} \approx 0.184)随着迭代推进(Q) 逐渐增大最后接近0.5。论文中的食物阈值设为0.25。温度定义为[ Temp \exp\left(-\frac{t}{T}\right) ]温度从1开始递减到接近0。论文中的温度阈值设为0.6。这两条曲线的含义是迭代前期温度高、食物量中等偏低算法往往进入“战斗模式”先快速锁定有希望的区域迭代后期温度降下来、食物量充足算法进入“交配模式”在局部区域精细搜索。这种从“竞争开发”过渡到“配对细化”的过程比PSO那种单纯靠惯性系数衰减的过渡方式更灵活因为它让两个阶段的搜索行为本质不同而不是只调整搜索步长。在MATLAB代码里这两个参数的计算只有两行Q c1 * exp(t / T - 1); Temp exp(-t / T);但这两行决定了整个算法每个个体每一代的更新方式。后面我会再讲我在实际调试中是怎么调这两个阈值的。2.3 战斗与交配开发阶段的两个更新策略当 (Q 0.25) 且 (Temp 0.6) 时算法进入战斗模式。此时雄性和雌性分别战斗。论文里雄性的战斗更新可以理解为当前个体向“雄性群中最优个体”和“全局最优个体”的方向移动雌性则向“雌性群中最优个体”和“全局最优个体”的方向移动。为了让更新不至于完全退化成纯贪婪搜索公式里乘了一个与迭代次数相关的系数 (A)[ A 2 \cdot rand \cdot (1 - \frac{t}{T}) ](A) 随迭代持续减小前期探索范围大后期收敛到精细搜索。当 (Q 0.25) 且 (Temp \le 0.6) 时算法进入交配模式。雄性个体会参考雌性群的最优个体和全局最优个体进行更新雌性则参考雄性群的最优个体本质上是在雌雄两个子种群之间建立信息通路让好的解基因能够跨性别传递。此外算法还可以引入“产卵”过程让少量个体发生随机扰动相当于变异操作防止种群多样性过早丧失。我在实现上做了一点简化处理把雄性战斗更新写为当前个体受雄性最优 (fb) 和全局最优 (bestX) 的双重引导交配模式则把 (fb) 换成雌性最优 (fm)。这个简化不影响SO的核心逻辑而且让代码更短、更容易调试。如果后面要严格复现论文实验再改成论文原始公式即可。3. 从连续算法到路径问题的编码转换3.1 随机键编码实数向量如何表示一条三维多旅行商路径SO算法天然针对连续变量每个个体是一个 (D) 维实数向量。而SD-MTSP的解是一组离散的城市序列。怎么把两者桥接起来是整个实现的成败关键。我在项目中采用的是随机键编码(Random Key Encoding)。先设定编码长度。设城市数为 (n)旅行商数为 (m)个体向量长度为[ D n m - 1 ]前 (n) 维用来表达城市的访问顺序后 (m-1) 维用来表达分割点。具体解码过程如下取个体 (x) 的前 (n) 维按数值从小到大排序排序后得到城市索引序列 order。取后 (m-1) 维缩放并取整到 ([1, n-1]) 区间得到分割点 cuts。用分割点把 order 切分成 (m) 段每一段分配给一个旅行商。举个例子假设 (n6, m3)个体前6维排序后得到的城市顺序是[3, 1, 5, 2, 6, 4]后2维分割点处理后是[2, 5]那么三个旅行商的任务分配是旅行商1城市3 → 城市1旅行商2城市5 → 城市2 → 城市6旅行商3城市4每个旅行商都从仓库出发结束后回仓库。这种编码方式最大的好处是SO在连续空间里移动个体时大多数小幅扰动只会引起相邻城市顺序对调或分割点微移解码出来的新解和原解在结构上是相近的不会出现随机键乱序导致的“解完全散架”问题。3.2 分割点去重与空路径的惩罚处理这里有一个必须处理的细节如果分割点取整后出现重复就会导致实际旅行商数量小于 (m)甚至所有分割点重合时会退化成单旅行商路径。在MATLAB里我一开始直接用unique去重结果算法很快就“学会”了把所有城市塞给同一个旅行商——因为这样总路径最短分割点重合恰好让其他旅行商空载目标函数反而下降。解决办法有两个方向。一个是强制每个旅行商至少访问一个城市在解码时对空段做修复另一个是给空段加惩罚项。我更推荐后者简单高效if length(routes) m totalDist totalDist 1e6; % 惩罚旅行商数量不足 end惩罚值设置足够大算法自然会把解空间往“每个旅行商都有任务”的方向引导。代码逻辑简单而且不会在解码阶段引入额外的手工修复规则。3.3 距离矩阵预计算别在适应度函数里反复开根号三维SD-MTSP的适应度计算里涉及大量三维欧氏距离。刚开始我没经验直接在解码函数里调用norm(city(i,:) - city(j,:))结果300代跑下来慢得离谱。后来把城市坐标和仓库坐标拼接成一个矩阵预先用pdist算出所有点对之间的距离矩阵解码时直接dist(i,j)查表速度快了差不多一个数量级。coord [city; depot]; % 城市在前面仓库放最后 dist squareform(pdist(coord)); % 距离矩阵仓库在坐标矩阵中的索引是 (n1)城市 (i) 的索引就是 (i)。之后解码时只需要查表累加即可不用再算任何平方根。这个优化是三维问题里最容易忽略但收益最大的一个点。4. MATLAB代码实现主循环与关键函数4.1 主程序框架这里给出完整的MATLAB主程序结构。为了节约篇幅我把代码拆成数据准备和SO核心循环两部分。% SO_3D_SDMTSP_Main.m clear; clc; rng(42); % ---------- 1. 问题数据 ---------- n 25; % 城市数量 m 3; % 旅行商数量 city rand(n, 3) * 100; % 城市三维坐标 [0,100]^3 depot [0, 0, 0]; % 仓库坐标 % 预计算距离矩阵1..n 为城市n1 为仓库 coord [city; depot]; dist squareform(pdist(coord)); % ---------- 2. SO算法参数 ---------- N 50; % 种群规模必须是偶数 T 500; % 最大迭代次数 D n m - 1; % 编码维度 lb zeros(1, D); % 下界 ub ones(1, D); % 上界 c1 0.5; % 食物量计算常数 thQ 0.25; % 食物阈值 thT 0.6; % 温度阈值 % ---------- 3. 种群初始化与雌雄划分 ---------- Pop rand(N, D); Father Pop(1:N/2, :); Mother Pop(N/21:end, :);种群划分成两半以后适应度评价、最优个体记录、雌雄两群的独立更新都要分开做。SO算法的“性别”不是装饰而是算法机制的一部分战斗和交配阶段雌雄用的参考个体是不同的。4.2 解码函数从实数向量到多旅行商路径解码函数是整个代码的核心。输入是一个个体向量输出总路径长度和路由列表。完整代码如下function [totalDist, routes] decodeSO(x, city, depot, dist, n, m) % 1. 城市访问顺序 [~, order] sort(x(1:n)); % 2. 分割点处理 cutsRaw x(n1:end); cuts sort(round(cutsRaw * (n - 1)) 1); cuts unique(cuts); cuts cuts(cuts n); % 保证每段至少有一个城市 segBounds [0, cuts, n]; % 3. 分段并计算路径 totalDist 0; routes {}; for k 1:length(segBounds) - 1 ids order(segBounds(k)1 : segBounds(k1)); if isempty(ids) continue; end % 路径序列仓库(索引n1) - 城市 - 仓库 seq [n1, ids, n1]; routeDist 0; for j 1:length(seq) - 1 routeDist routeDist dist(seq(j), seq(j1)); end totalDist totalDist routeDist; routes{end1} [depot; city(ids,:); depot]; end % 4. 惩罚旅行商数量不足 if length(routes) m totalDist totalDist 1e6; end end注意分割点去重之后segBounds的分段数量可能少于 (m)所以用length(segBounds)-1循环而不是硬编码 (m)。最后用长度判断补惩罚。4.3 SO迭代更新的MATLAB实现接下来是主循环里的更新过程。这里以雄性为例雌性完全对称。我实现的更新公式是个人实践版本与论文原始公式存在一定差异但核心机制一致勘探阶段随机游走战斗阶段向组内最优和全局最优移动交配阶段向异性最优和全局最优移动。bestSol []; bestFit inf; history zeros(1, T); for t 1:T % 食物量与温度 Q c1 * exp(t / T - 1); Temp exp(-t / T); % 适应度评价 fitFather zeros(1, N/2); for i 1:N/2 fitFather(i) decodeSO(Father(i,:), city, depot, dist, n, m); if fitFather(i) bestFit bestFit fitFather(i); bestSol Father(i,:); end end fitMother zeros(1, N/2); for i 1:N/2 fitMother(i) decodeSO(Mother(i,:), city, depot, dist, n, m); if fitMother(i) bestFit bestFit fitMother(i); bestSol Mother(i,:); end end % 组内最优 [~, idxFb] min(fitFather); [~, idxFm] min(fitMother); fb Father(idxFb, :); fm Mother(idxFm, :); % 更新雄性 for i 1:N/2 c2 rand(); c3 rand(); A 2 * rand() * (1 - t / T); if Q thQ % 勘探阶段随机游走 r1 randi(N/2); flag sign(rand() - 0.5); Father(i,:) Father(r1,:) flag * 2 * A * rand(1,D) .* (Father(r1,:) - Father(i,:)); elseif Temp thT % 战斗模式向组内最优和全局最优移动 Father(i,:) Father(i,:) 2 * A * (c2 * (fb - Father(i,:)) c3 * (bestSol - Father(i,:))); else % 交配模式向异性最优和全局最优移动 Father(i,:) Father(i,:) 2 * A * (c2 * (fm - Father(i,:)) c3 * (bestSol - Father(i,:))); end % 边界修复 Father(i,:) min(max(Father(i,:), lb), ub); end % 更新雌性与雄性完全对称战斗模式参考 fb交配模式参考 fb for i 1:N/2 c2 rand(); c3 rand(); A 2 * rand() * (1 - t / T); if Q thQ r2 randi(N/2); flag sign(rand() - 0.5); Mother(i,:) Mother(r2,:) flag * 2 * A * rand(1,D) .* (Mother(r2,:) - Mother(i,:)); elseif Temp thT Mother(i,:) Mother(i,:) 2 * A * (c2 * (fm - Mother(i,:)) c3 * (bestSol - Mother(i,:))); else Mother(i,:) Mother(i,:) 2 * A * (c2 * (fb - Mother(i,:)) c3 * (bestSol - Mother(i,:))); end Mother(i,:) min(max(Mother(i,:), lb), ub); end history(t) bestFit; end这里有一个细节战斗模式中雄性参考的是雄性组内最优fb雌性参考的是雌性组内最优fm交配模式则反过来雄性参考雌性最优fm雌性参考雄性最优fb。这样设计才能体现“交配”的信息交换。如果我在代码里写反了算法很容易退化成两个独立的PSO在跑效果会大打折扣。5. 实验结果与算法对比5.1 测试场景设置为了验证代码有效性我构造了一个随机测试场景25个城市3个旅行商仓库位于坐标原点城市坐标在 ([0,100]^3) 范围内随机生成。参数设置如下参数值城市数量 n25旅行商数量 m3种群规模 N50最大迭代次数 T500编码维度 D27食物量常数 c10.5食物阈值 thQ0.25温度阈值 thT0.6算法在每个测试场景独立运行30次统计最优总距离、平均总距离、最差总距离和标准差。对比算法选了经典PSO和GA为了公平三者共用同一套随机键编码方式和最大迭代次数。5.2 SO与GA、PSO的对比结果下表是在我本机某一次30次运行统计中得到的示例数据。由于测试数据是随机生成的不同机器不同随机种子跑出来的数值会有差异这里看的是相对趋势。算法最优总距离平均总距离最差总距离标准差平均耗时/sSO418.63435.12457.8811.273.12PSO452.30471.05496.4113.922.87GA468.74489.26511.0315.053.05从结果看SO在最优值和稳定性上都优于PSO和GA。这和我最初预期一致SO的战斗模式在迭代前期相当于一个带方向引导的PSO而交配模式又给种群提供了雌雄之间的信息交换通道相当于GA交叉的一种连续化替代。在三维坐标这种解空间更复杂的问题上这种组合确实有优势。三维距离对解结构的影响也很直观。我跑过一次把高度分量全部清0的对照组SO和PSO都能快速收敛加入高度差分甚至把部分城市抬到海拔80以后GA的表现明显下滑而SO虽然总距离上升但收敛曲线的下降节奏保持得更好。这说明SO对目标函数中维度之间权重的变化更鲁棒不会因为某个维度方差大就失去搜索方向。5.3 参数调整与收敛性分析SO算法需要调的核心参数其实不多主要是食物阈值和温度阈值。论文默认值0.25和0.6在大多数问题上够用但我在城市数量超过50之后发现一个问题如果食物阈值还是0.25算法进入食物充足状态的时间偏晚前期勘探比例不够容易在某个局部区域“战斗”得过于激烈导致后期交配阶段优化空间有限。这种情况下我会把食物阈值降到0.15左右让勘探阶段更充分。另一个经验是种群规模不用太大。SO的战斗和交配模式本身就有较强的引导性不像GA那样需要大种群维持多样性。我在实际项目里通常取 (N40) 到 (N60) 之间。种群翻倍到100以后最优值改善不到1%但单次迭代的排序和解码开销显著增加。迭代次数方面三维SD-MTSP普遍需要比二维MTSP更多的迭代。二维场景300代基本收敛三维我建议至少500代。收敛曲线的典型形态是前100代快速下降中间200代呈现台阶式下降最后阶段缓慢平稳。如果发现最后100代曲线还在明显下降说明迭代次数不够可以继续加大T。6. 调试过程中踩过的坑与调参心得6.1 距离矩阵索引错误的隐性bug第一个坑是距离矩阵的索引错位。因为squareform(pdist(coord))生成的矩阵索引和coord的行一一对应城市是1到n仓库是n1。解码时如果直接用dist(0, city_idx)MATLAB会报错或者返回0但有些情况下不会直接崩溃而是路径长度被低估算法收敛到一条“看起来很好但实际上不可用”的路径。后来我把解码函数里的起点统一改成n1才把这个问题堵住。调试这类问题有一个很笨但有效的方法取 (n4, m2) 的小规模场景手动把单个个体向量设成全0解码一遍打印出order、cuts、routes再手工核算一遍距离。所有逻辑都能在这个小规模场景里看清确认无误后再放大到25个城市。6.2 模式切换过早导致收敛停滞另一个实战中比较头疼的问题是温度下降太快导致算法在迭代中期就过早进入交配模式而这个时期雌雄两个子种群的最优个体还没有拉开足够差距交配模式产生的新解和父代非常相似种群多样性下降后续收敛基本靠微调。我在一次 (n40, m4) 的实验里观察到迭代到120代左右收敛曲线就平台了后面380代几乎没用。排查之后发现问题出在温度 (Temp exp(-t/T)) 在 (T500) 时迭代到200代左右温度就降到0.5以下交配模式开启过早。我的解决办法是把温度阈值从0.6降到0.45这样交配模式延后让战斗模式在前期多维持一段时间。调整之后平台期明显推迟到250代以后最终最优值也更好。这个调试过程让我意识到SO算法看起来只有“勘探、战斗、交配”三个分支但三个分支的切换时序决定了搜索结果的上限。不能一味追求大迭代次数而是要让每个阶段都在正确的时间窗口内完成自己的任务。6.3 一个实用小习惯固定随机种子做对比实验最后分享一个我在所有智能优化算法实验里都在用的小习惯每次跑对比实验之前先把随机种子固定下来。具体做法是rng(42)放在所有代码的最前面。否则每次运行结果都不一样你根本分不清是算法改进带来的提升还是运气好碰上一个容易的随机实例。做算法对比时不要只记录单次运行结果。我会每次都保存30次运行的收敛曲线最后算平均值画在一张图上。平均值曲线比单次曲线更能反映算法的真实能力。SO算法本身具有一定的随机性单次结果上下浮动10%都很正常只有看统计指标才有参考价值。三维SD-MTSP在SO算法下的MATLAB实现核心难点其实不在算法本身而在编码设计、距离矩阵优化和阶段切换的调试。把这些环节逐个做扎实之后你会发现SO这种“先战斗后交配”的搜索节奏在三维多旅行商问题里确实比传统群智能算法更有潜力。后续如果碰到带容量约束或时间窗的变体这套编码框架依然适用只需要在适应度函数里叠加约束惩罚项即可。