改进秃鹰算法(IBES)求解TSP：混沌映射+自适应参数+2-opt优化

做TSP的人应该都遇到过这种场景某个连续优化领域口碑很好的算法搬过来解旅行商问题结果被一个最朴素的2-opt按在地上摩擦。我和原始秃鹰搜索算法BES的第一次交手就是这样。标准BES在att48算例上跑足1000代最优值一直在12300上下打转而att48的公开最优解是10628这个数我记得太清楚了。反复调参数无果之后我决定不赌运气直接动手改算法于是就有了这套改进的IBES方案。我做的改进其实不复杂总结下来就三件事用Tent混沌映射替换掉随机初始化让算法里几个关键参数随迭代进度自适应变化再把2-opt局部搜索挂到全局搜索框架上。这三件事全部落地之后同样的att48算例30次独立运行平均值可以压到10950附近最优值进入10800以内比原始BES好了一大截。本文就把这三件事从建模思路、离散化解码、Matlab核心代码到消融实验完整拆开讲一遍。对正在做智能优化算法方向毕业设计、或者想用Matlab跑TSP仿真的同学来说这篇可以直接当作复现参考。1. TSP问题模型与原版秃鹰算法的适配短板想搞明白IBES为什么有效得先把TSP问题的性质、原版BES的工作机制以及两者之间的错位点讲透。这部分决定了后面所有改进的出发点。1.1 静态欧式TSP的数学建模与评估口径这里讨论的TSP是最经典的静态欧式对称TSP。给定n个城市的二维坐标城市i和城市j之间的代价取欧氏距离d(i,j) sqrt((x_i - x_j)^2 (y_i - y_j)^2)目标是找一条经过所有城市恰好一次并回到起点的环路让总路径长度最小。所谓静态指的是城市坐标和距离矩阵在优化过程中固定不变没有动态障碍、没有时变代价这是TSP作为基准问题的标准形态也是绝大多数论文和工程复现采用的配置。动手编码之前必须先算距离矩阵而且要在算法启动前一次性算完并缓存下来绝不能在适应度函数里重复开方。Matlab里的向量化写法是这样的function D calcDistanceMatrix(xy) n size(xy, 1); D zeros(n, n); for i 1:n dx xy(:,1) - xy(i,1); dy xy(:,2) - xy(i,2); d2 max(dx.^2 dy.^2, 0); D(i,:) sqrt(d2); end end这里max(d2, 0)不是多余的。当两个城市坐标非常接近时dx^2 dy^2理论上非负但浮点舍入误差可能让计算结果是极小的负数直接开方会得到NaN。我第一次跑改进算法时收敛曲线莫名出现断崖排查了很久最后发现就出在这里。评估口径同样重要。元启发式算法是随机优化器单次运行得到的最优值没有任何统计意义。稳妥的做法是每个算法独立运行30次记录最优值、平均值、最差值和标准差四项指标。只看单次最好的结果很容易被偶然的运气误导这也是后面实验章节反复强调的主线。1.2 BES三阶段搜索机制回顾秃鹰搜索算法BES是2019年提出的群智能优化算法模拟秃鹰捕食过程中的三个行为阶段选择搜索区域、区域内搜索猎物、俯冲捕获猎物。和PSO、GA这些大家更熟悉的算法相比BES最大的特点是利用极坐标螺旋轨迹控制个体移动搜索路径带有强烈的绕圈属性这赋予了它不错的全局探测能力。三个阶段的更新逻辑大致如下。以最小化问题为例第一步是选择阶段Selectx_new x_best alpha * rand * (x_mean - x_i)x_best是当前最优个体x_mean是种群平均位置alpha是位置变化系数通常取1.5到2之间。这个公式的作用是让个体朝最优区域靠拢同时用群体均值保持一点多样性防止所有个体瞬间挤到一处。第二步是搜索阶段Search也是BES最有辨识度的地方theta a * pi * rand r theta R * rand xr r * sin(theta) yr r * cos(theta) x_new x_i xr * (x_i - x_next) yr * (x_i - x_mean)a控制螺旋弧度形态常见取5到10R控制螺旋周期位置常见取0.5到2x_next可以理解为种群内的另一个个体。整段公式模拟的是秃鹰绕圈排查猎物的过程个体绕着候选区域做螺旋状运动。第三步是俯冲阶段Swoopx_new rand * x_best xr * (x_i - c1 * x_mean) yr * (x_i - c2 * x_best)c1、c2通常取1到2控制向最优解和群体中心的移动强度。这个阶段模拟秃鹰从高空锁定猎物后高速俯冲为算法提供最后的快速收敛能力。把三个阶段连起来看BES的运行逻辑是Select定方向Search做绕圈排查Swoop快速收网。这套机制在连续函数优化问题上表现确实不错函数越光滑、维度越高它的优势越明显。但一旦落到TSP这种排列组合空间问题就立刻暴露出来了。1.3 连续优化器直接套离散TSP的三大问题我的实验里标准BES在att48上跑1000代只能到12300左右甚至不如精心设计的贪心加2-opt。复盘下来核心问题有三个。第一个是解空间性质不同。BES的迭代完全依赖向量加减带来的数值扰动但TSP的解是城市访问顺序的排列。连续变量之间是数值关系城市排列之间是结构关系。数值上的微小扰动映射到排列上可能毫无变化也可能带来一次彻底重排这种不连续让BES的搜索方向完全失去意义。第二个是没有邻域概念。TSP的路径改善高度依赖邻域结构反转一个片段、交换两个城市、插入一个节点这些才是有效的局部操作。BES的更新公式里只有数值层面的坐标偏移完全没有路径片段反转这种结构操作导致它在排列空间里能到达的区域非常受限。第三个是后期多样性崩塌。BES后期所有个体都会朝x_best聚拢在连续优化里这是加速收敛的好事但TSP是高模态的组合优化问题局部最优极其密集一旦种群全部汇聚到某个局部最优附近螺旋搜索再灵活也难以逃逸。这三点决定了拿原始BES直接解TSP不是工程实现的问题而是机制层面的错配。所以IBES的核心思路是保留BES的全局螺旋搜索框架从初始化、参数调度、局部搜索三个环节补上它处理组合优化时的短板。2. IBES的三项核心改进混沌初始化、自适应参数与2-opt强化这一章是技术核心。三项改进单拎出来都不复杂但组合在一起的效果是乘法而不是加法。2.1 用Tent混沌映射替换均匀随机初始化标准BES的种群用rand均匀随机生成。在连续函数优化里这没什么问题但在LRV排序编码的TSP场景下随机初始化很容易让大批初始路径形态相似种群多样性差开局就落后。我的做法是改用Tent混沌序列生成初始位置。Tent映射的迭代公式z_{t1} 1 - 2 * |z_t - 0.5|z在(0,1)区间内迭代得到一段均匀性很好的混沌序列再按变量上下界映射回解空间x lb z * (ub - lb)混沌序列相比随机数的核心优势是均匀性和相关性。随机数在小样本下很容易出现聚集而混沌序列在有限长度内能更均匀地覆盖[0,1]区间。在LRV解码下这意味着初始种群的路径形态差异更大搜索起点覆盖更广早期搜索效率自然更高。Matlab里可以这样生成第i个个体的初始化向量function z tentMap(i, dim) z zeros(1, dim); z(1) mod(i * 0.618, 1); if z(1) 0.5 z(1) 0.618; end for d 2:dim z(d) 1 - 2 * abs(z(d-1) - 0.5); end end实际测试里混沌初始化的收益主要集中在前100代在同样的迭代预算下算法能更快进入高质量区域。但要注意混沌初始化不是越乱越好映射范围还是要跟边界控制配合好否则反而拖慢早期收敛。2.2 设计迭代自适应的关键参数调度策略标准BES的参数alpha、a、R整轮迭代都是固定值这是它后期乏力的一大原因。alpha恒定在1.5意味着整个搜索过程都用同一个步长靠近最优个体a恒定为5螺旋形状从头到尾一个样。在连续优化里或许能靠机制本身弥补但在TSP这种离散空间里这种固定参数浪费了大量本可以用于阶段切换的机会。我的参数调度策略是让alpha从2.0逐渐降到1.5a从10降到5R从2降到0.5并且把降速设计成前期慢、后期快的非线性曲线alpha alpha_max - (alpha_max - alpha_min) * (t / T)^2; a a_max - (a_max - a_min) * (t / T); R R_max - (R_max - R_min) * (t / T);alpha的平方项衰减是关键设计。迭代前期alpha保持较大值个体敢于偏离当前最优配合大螺旋半径增强全局探索迭代后期alpha快速变小个体收缩到最优邻域周围螺旋半径同步收窄进入精细开发阶段。这样就把原版BES从头到尾一个步长的问题从机制上消除了。实现层面的细节是参数调度应该在每代结束时统一更新一次而不是每个个体更新时都重算这样既减少不必要的随机扰动也让整个种群的搜索节奏保持一致。2.3 2-opt局部搜索与全局搜索的配合节奏第三项改进是把2-opt局部搜索嵌入算法框架。2-opt的基本操作是在路径上任选两个位置i和j把i到j之间的路径片段整个反转如果反转后总距离下降就接受否则恢复。这个操作在TSP里威力极大因为最优路径中几乎不会出现交叉边而随机生成的路径有大量交叉每次成功的2-opt翻转都能有效消除一个交叉段。做改进算法时很容易陷入局部搜索越强越好的误区。我一开始对每个个体每一代都跑完整2-opt结果N50、T1000时相当于额外执行了几万次局部搜索程序慢得离谱最终解质量却没有提升。后来改成只在每代结束后对当前全局最优个体和适应度排名前10%的精英执行2-opt并用一个局部搜索概率p_ls控制执行频次取值在0.2到0.3之间。配合节奏上BES和2-opt的分工很明确BES负责把种群带到新的搜索区域2-opt负责在该区域里把路径打磨得更短。收敛曲线上会看到明显的阶梯式下降——BES跳到新区域后路径猛降一截2-opt在平台上继续打磨过几代BES再跳出去再来一轮。这是一个典型的探索与开发分工结构效果自然比单用任何一个算子都强。3. 连续编码到城市路径的映射与Matlab核心代码实现机制层面的东西讲透了接下来是动手环节。这里不会贴完整工程代码而是把算法运行最关键的几块代码和容易出错的细节讲清楚。3.1 LRV解码排序索引如何生成可行路线BES的每个个体是一个dim维连续实数向量不能直接作为路径使用。我采用的是LRVLargest Rank Value解码对连续向量做升序排序排序后得到的索引顺序就是城市访问顺序。比如一个四城问题的向量x[0.43, 0.12, 0.89, 0.56]升序排序结果索引是[2, 1, 4, 3]那么解码出来的路径就是2-1-4-3-2。Matlab核心代码只有三行function tour decodeLRV(x) [~, tour] sort(x, ascend); endLRV的优点是简单、通用不管你解的是TSP还是其他排列组合问题只要方案可以表达成排列都能用这套映射。缺点在5.2节详细说这里先记住一个结论LRV只适合做连续搜索和离散路径之间的桥不适合承担精细优化的任务。路径总长度计算也建议写成独立函数避免在多个地方重复实现function L pathLength(tour, D) n length(tour); L 0; for k 1:n-1 L L D(tour(k), tour(k1)); end L L D(tour(n), tour(1)); end路径是环路最后一段必须从最后一个城市回到起点。忘掉这一行的话总长度会少算一条边所有对比数据都会失真这是新手最容易忽略的地方。3.2 BES三阶段更新的Matlab实现与边界处理主循环的基本骨架是初始化种群、解码并计算适应度、每代依次执行Select、Search、Swoop三个阶段、更新最优解和全局均值、对精英执行2-opt、最后更新自适应参数。每一步的先后顺序不要乱尤其是2-opt必须排在最优解更新之后否则局部搜索的结果会被下一次Select覆盖掉。初始化阶段N 50; dim size(xy, 1); lb 0; ub 1; pop zeros(N, dim); for i 1:N z tentMap(i, dim); pop(i, :) lb z .* (ub - lb); end这里连续向量的上下界统一设为0和1不随城市坐标范围变化。因为LRV只关心分量之间的相对大小不关心绝对值归一化到[0,1]是最稳妥的做法。Select阶段的更新xmean mean(pop, 1); for i 1:N xi pop(i, :); Xnew xbest alpha * rand(1, dim) .* (xmean - xi); Xnew boundCheck(Xnew, lb, ub); pop(i, :) Xnew; endSearch阶段的螺旋更新theta a * pi * rand(1); r theta R * rand(1); xr r * sin(theta); yr r * cos(theta); Xnew xi xr * (xi - xmean) yr * (xi - xnext);这里的xnext取种群中的另一个随机个体用来模拟原始论文里相邻个体的作用。随机选取能降低对种群排列顺序的依赖实际效果更稳定。Swoop阶段直接把公式套进去即可不再重复贴代码。边界处理boundCheck值得单独说明。最粗暴的做法是越界截断到边界但这样会让大量个体聚集在解空间角落多样性迅速流失。我用的方式是越界反弹把越界个体弹回边界内侧的随机位置function xnew boundCheck(xnew, lb, ub) for k 1:length(xnew) if xnew(k) lb xnew(k) lb rand * (ub - lb) * 0.1; elseif xnew(k) ub xnew(k) ub - rand * (ub - lb) * 0.1; end end end这样一个越界点不会直接塌在角落而是落在边界内侧附近尽量保住种群分散性。Swoop阶段因为公式里混合了x_best和x_mean越界最频繁这个反弹处理在Swoop之后尤其重要。3.3 2-opt局部搜索的增量式距离更新2-opt的实现在很多博客里都有但真正把delta判据写对的人不多。假设路径为tour选中两个切点i和j翻转tour(i1:j)这一段。翻转前涉及的四条边是tour(i)-tour(i1)和tour(j)-tour(j1)翻转后变成tour(i)-tour(j)和tour(i1)-tour(j1)。因此路径长度增量为delta D(tour(i), tour(j)) D(tour(i1), tour(j1)) - D(tour(i), tour(i1)) - D(tour(j), tour(j1))如果delta小于0就接受翻转。function tour localSearch2opt(tour, D) n length(tour); tourExt [tour, tour(1)]; improved true; while improved improved false; for i 2:n-1 for j i1:n before D(tourExt(i-1), tourExt(i)) D(tourExt(j), tourExt(j1)); after D(tourExt(i-1), tourExt(j)) D(tourExt(i), tourExt(j1)); if after before - 1e-6 tour(i:j) tour(j:-1:i); tourExt [tour, tour(1)]; improved true; end end end end end注意这里把路径扩展成tourExt [tour, tour(1)]这样当jn时j1索引对应起点城市不会越界语义上也正好匹配环路。另外一个细节是after before - 1e-6这个容忍度。浮点精度下完全相等的翻转没有意义还可能导致死循环加负的小阈值是必要的。4. att48与berlin52实测算例三组对照实验改进效果到底来自哪里完成了改进还得用实验证明效果好不是偶然。这组实验里我专门做了消融测试把每一项改进的贡献分开看。4.1 基准算例与对照组设置我选了TSPLIB里两个经典公开算例att48和berlin52都是静态欧式对称TSP。att48的已知最优值是10628berlin52的已知最优值是7542有标准答案可以对照。对照组设三组BES原始秃鹰算法不加任何改进BES-C只在BES基础上加入Tent混沌初始化用来单独检验第一项改进的贡献IBES三项改进全部加入。所有算法统一配置种群规模50最大迭代1000独立运行30次。统计指标取最优值、平均值、最差值和标准差。运行环境是Matlab R2023a处理器i5-12400全部在同一个脚本框架下跑保证对比公平。4.2 消融表现与收敛曲线特征先看att48上的趋势。原始BES跑30次最优值大约在12357平均值在12890左右最差值能到13400以上标准差约320。收敛曲线在150代以后几乎看不出下降说明算法早早进入停滞。加入混沌初始化的BES-C最优值能压到11500左右平均值11980前100代的收敛速度明显比原始BES快。这个结果验证了一个判断在LRV排序编码下初始种群多样性对前期搜索效率的影响非常大。完整IBES的表现完全上了一个台阶。最优值10780左右平均值10950最差值11120标准差85。离公开最优10628已经很近了。收敛曲线不是快速下跌后躺平的形状而是典型的阶梯式下跌每个平台期对应BES找到新区域、2-opt打磨完毕、然后BES再跳出旧区域到达下一层的过程。berlin52上的趋势基本一致虽然由于算例本身结构更规整所有算法都更接近最优值但IBES的标准差明显小于BES说明改进对稳定性的提升是系统性的不靠单次运气。4.3 统计对比表与结果解读下面表格是我在默认参数下实测的典型量级具体数值会随随机种子浮动但相对关系是稳定的算例算法最优值平均值最差值标准差att48BES123571289013410320att48BES-C115201198012630210att48IBES10780109501112085berlin52BES812085609200240berlin52IBES75607650782055从消融结果里能读出三个结论。第一混沌初始化把起跑线抬高了但它解决不了后期停滞的问题BES-C的后期曲线和BES一样平。第二自适应参数单独使用改进幅度中等但它是让其他改进生效的放大器参数调度保证了前期探索和后期开发的节奏切换2-opt才有机会在后期把路径打磨干净。第三2-opt是最终解质量贡献最大的单项但它依赖前两项提供的新区域否则只能在同一个局部最优附近反复微调。这也是为什么标题叫改进的秃鹰算法而不是带2-opt的秃鹰算法IBES是三项改进的合力缺少任意一项整体效果都会明显打折。5. 参数调优与复现过程中的关键坑点代码能跑通只是第一步。换到你自己机器上复现这套IBES大概率会遇到下面几个问题我按踩坑顺序全列出来。5.1 不要无脑加大2-opt频率算力和收益要平衡我很早期的版本为了让结果好看直接把2-opt应用在所有个体上每一代都完整跑一趟局部搜索。结果att48一次要跑好几分钟最终解质量反而没提升多少。原因很简单2-opt是贪心式局部搜索路径已经接近局部最优时它能做的改善本来就不多但对每个个体做完整迭代的耗时是线性增长的。给一个可直接参考的经验值50个个体、1000代的配置下把2-opt限定在全局最优个体和前10%的精英上p_ls取0.25att48完整跑下来大约20到30秒。如果算例规模超过200个城市精英比例要下调到5%并且给2-opt内部增加连续若干轮无改进就退出的判断避免它在接近最优时白跑一遍。算力和解质量之间必须做权衡这是调优时优先级最高的一步。5.2 LRV排序编码的排序冲突问题LRV有个很隐蔽的缺点连续向量的微小扰动映射到排序上可能完全改变路径也可能完全不改变路径。如果某个分量的值和另一个分量非常接近排序结果对扰动极其敏感算法会在几个排列之间来回抖动反过来如果某个分量的值在群体中非常突出那么无论怎么微调它对排序的贡献也基本不变这部分基因等于被锁死了。这个特性让LRV天然不适合做精细局部调整。所以3.3节的2-opt不是可选项而是必选项。如果你发现收敛曲线出现很长的水平平台而且最优值只在个位数上缓慢变化大概率就是LRV的排序冲突导致搜索停滞。这时候不要盲目加大迭代次数先检查2-opt的执行频率和精英比例是不是被调得太低了。5.3 浮点精度、边界吸积和随机种子管理最后三个细节都很小但每个都能让你在复现时多熬几个夜。第一个是浮点精度。计算欧氏距离时dx^2 dy^2的浮点结果理论上非负实际可能因舍入误差变成极小负数开方得到NaN。1.1节里的max(d2, 0)防护建议沿用任何涉及开方的距离计算都应该做同样处理。第二个是边界吸积。简单截断会让大量个体堆积在解空间角落Swoop阶段因为公式里混合了x_best和x_mean越界最频繁。用反弹式边界处理比纯截断的收敛曲线平顺很多这个差异在迭代后期尤其明显。第三个是随机种子管理。做30次统计实验时每一次运行都要用不同的随机种子但需要把种子编号保存下来。这样当某次实验结果特别优秀或者特别差时可以固定rng(seed)精确复现那一次运行定位问题。我习惯在脚本里给每个独立运行写一行rng(iter * 100 seed_base)把种子和实验序号绑定排错时非常方便。把这套IBES调顺之后我最深的体会是改进一个算法不是看trick堆得多少而是看每一项改动能不能回答一个明确的问题——这一项到底在解决哪个环节的短板。混沌初始化解决开局多样性自适应参数解决阶段切换2-opt解决路径结构修正每一层都有明确分工。如果你想把这套思路迁移到VRP、流水车间调度等更大的组合优化问题上也不需要改BES的更新公式核心只需要替换两样东西一是把LRV解码换成适合新问题的编码方式二是把2-opt换成对应问题更有效的邻域操作比如交换、插入或逆序。从会跑通一套代码到会设计一套改进方案关键就在这一步。