简介Matlab混合粒子群算法HPSO求解TSP问题的完整代码实例面向需要掌握启发式优化与组合优化实践的开发者、研究生及相关工程人员。压缩包共4个文件包含3个M脚本和1个TXT数据文件覆盖主程序、适应度评估、距离计算等核心模块并配有城市坐标数据可直接运行验证3KB体量轻量易读。资源已有2167人学习实用价值得到验证。通过源码注释可清晰理解PSO初始化、适应度计算、全局最优更新、速度位置更新以及引入混合策略改善早熟收敛的完整流程调整惯性权重、学习因子与混合算子参数可进一步观察算法在不同规模TSP问题上的收敛性与稳定性便于迁移至路径规划、调度优化等实际场景。1. 混合粒子群算法解 TSP为什么标准 PSO 会卡在“半路”旅行商问题TSP在组合优化里是最直白也最折磨人的一类城市一多路径数量爆炸式增长精确算法在 50 个城市以上就很难硬算。工程上更常见的做法是用启发式算法逼近最优解而粒子群优化PSO是其中收敛最快的一类。但如果直接把标准 PSO 的速度-位置更新公式套到 TSP 上会发现粒子在连续空间里飞得再快也映射不成一条合法的城市访问序列——这就是标准 PSO 和离散组合优化之间那道著名的鸿沟。混合粒子群算法HPSO的思路是在 PSO 框架里引入遗传操作或局部搜索把离散路径的“交叉”和“变异”变成粒子更新的手段。这套代码以 TSPLIB 的 eil51 实例为测试数据包含 main.m、fitness.m、dist.m 三个核心文件带注释、可直接运行适合两类人一是刚接触启发式算法、想看懂 PSO 在离散问题上怎么落地的新手二是已经在用遗传算法或模拟退火解 TSP、想对比混合策略收益的熟手。接下来按“建模 → 主循环 → 混合策略 → 调参验证”的顺序拆开讲每个环节都给出可复现的命令和参数。2. TSP 的粒子和适应度建模从连续速度到交换序2.1 城市编码与距离矩阵的两种构建方式TSP 的解是一条闭合环游路径所以每个粒子必须表示成一个离散的城市排列例如[3, 7, 1, 4, 2, 6, 5]表示从城市 3 出发、依次经过 7/1/4/2/6/5、最后回到 3。这个排列就是粒子在搜索空间中的“位置”。速度向量在连续 PSO 里是位移增量但在 TSP 里速度必须设计成“交换序”——一组有限次的交换操作。交换序可以理解为要对当前路径做哪些位置交换才能变成一个新路径。代码包里的 eil51.txt 是 TSPLIB 的 51 城市欧氏距离实例每行是城市编号和二维坐标城市编号不要求连续读取时要按行处理成矩阵。先用load或importdata读入后调用 dist.m 计算距离矩阵。dist.m 的核心逻辑如下function D dist(coords) n size(coords, 1); D zeros(n, n); for i 1:n for j 1:n D(i, j) sqrt(sum((coords(i,:) - coords(j,:)).^2)); end end end这段代码逐对计算城市间的欧氏距离生成一个 n×n 的对称矩阵对角线全为 0。距离矩阵 D 是后续所有适应度计算的基础它只需要在初始化阶段算一次之后每次评估路径长度直接查表不需要重复计算坐标。对于 51 个城市这种方法完全够快如果城市规模到了 500 以上可以考虑用pdist2或分块向量化来加速但可读性会下降。注意 TSPLIB 的 eil51 最优解距离目前公认是 426部分资料写 425.2取决于四舍五入规则后面验证算法时要以这个数为标杆。如果读入的坐标顺序是正确的dist.m 算出的距离矩阵应该让完整环游的最短距离落在 426 附近这能提前验证数据是否读对。2.2 适应度函数与非法路径的惩罚处理适应度函数在 TSP 里就是路径总长度公式为从第一个城市出发依次累加相邻城市距离最后还要加上最后一个城市回到第一个城市的距离。若粒子代表的排列是route长度为 L则路径总距离为D(route(1), route(2)) D(route(2), route(3)) ... D(route(L), route(1))。function fit fitness(route, D) L length(route); dist_sum 0; for i 1:L-1 dist_sum dist_sum D(route(i), route(i1)); end dist_sum dist_sum D(route(L), route(1)); % 回到起点 fit dist_sum; end这个函数把所有相邻距离累计最后闭合回路。注意一个常见误用有人只累加到第 L-1 个城市忘记加回起点导致适应度偏小、粒子“看起来”更优实际上路径不闭合。另一个问题是路径编码中出现重复城市有些变体允许这种现象然后加惩罚项但更干净的做法是在速度更新和交叉操作里强制保持排列的合法性——所有混合策略的操作都必须输出一个“恰好包含每座城市一次”的排列这样适应度函数就无需处理非法路径。这里的适应度是标量且越小越好和标准 PSO 的“最小化问题”直接对应不需要额外转换。3. HPSO 核心循环main.m 的位置更新与记忆机制3.1 粒子群参数与速度定义打开 main.m第一段通常是参数初始化这是整个算法最容易改也最容易改错的地方。代码中常见的参数设置如下nPop 40; % 粒子数量 MaxIter 500; % 最大迭代次数 w 0.8; % 惯性权重 c1 1.5; % 个体学习因子 c2 1.5; % 全局学习因子粒子数量的选择与城市规模相关对于 eil51 这种中小规模实例3050 个粒子足够如果城市规模到 100 以上建议提高到 80100否则搜索覆盖面不够。MaxIter 直接影响运行时间和最终质量500 次迭代配合 40 个粒子在 51 城市上通常能在 10 秒内完成收敛但更稳妥的做法是先跑 300 次迭代看收敛曲线是否已经平台化再决定是否加长。惯性权重 w 控制粒子继承上一时刻“速度”的程度0.8 属于中等偏大的设置有较强的全局搜索能力如果发现收敛后晃动大可以线性递减到 0.4。c1 和 c2 分别控制向个体历史最优和全局历史最优学习的强度经典取法是 c1 c2 1.5如果算法陷入局部最优可以适当调大 c2 但不要超过 2.0。每个粒子除了位置路径排列还需要两个记忆个体最优路径pbest和个体最优长度以及全局最优路径gbest和全局最优长度。初始化时每个粒子随机生成一个全排列并计算适应度将 pbest 初始化为自身gbest 初始化为所有粒子中最优的那个。3.2 交换序与交换基本子的速度更新公式标准 PSO 的速度更新公式是v w*v c1*r1*(pbest - x) c2*r2*(gbest - x)但这个公式在 TSP 里没法直接运算因为粒子位置不是连续向量。混合粒子群算法的做法是把速度定义为交换序的集合即一组“交换哪个位置”的操作列表。两个粒子之间的“差”定义为从粒子 A 变为粒子 B 所需要的最少交换序列。这样速度更新就变成了对当前位置施加一系列交换操作。% 计算 pbest 与当前粒子的差异交换序 diff_pbest path_diff(route, pbest); % 计算 gbest 与当前粒子的差异 diff_gbest path_diff(route, gbest); % 按概率取交换操作 v_new []; for k 1:length(diff_pbest) if rand c1 v_new [v_new, diff_pbest(k)]; end end for k 1:length(diff_gbest) if rand c2 v_new [v_new, diff_gbest(k)]; end end % 应用到当前位置 route apply_swap_seq(route, v_new);这里的path_diff逐个比对两个排列中每个位置的元素如果不一致就找到目标排列中该元素的位置进行一次交换然后记录这个“交换基本子”。apply_swap_seq则依次执行这些交换操作得到新路径。参考粒子速度的惯性部分可以用一个小概率随机扰动替代即每次有一定概率随机交换当前路径中的两个位置相当于保留了探索能力。整个过程并不难但要注意交换序的叠加顺序会影响最终结果所以每次执行交换必须严格按顺序应用不能打乱。参数 c1 和 c2 在这里不再作为物理加速度系数而是作为“接受这次交换操作的概率阈值”。如果某个粒子与 pbest 差异很大diff_pbest长度就长按概率选取后施加的交换就多反之则少。这种机制保持了 PSO 的核心思想向个人最优和全局最优靠拢只是驱动力变成了离散的交换操作。3.3 早熟收敛的判断与重启策略TSP 的适应度景观非常崎岖HPSO 运行到中后期容易遇到所有粒子都聚集到同一个局部最优附近此时 pbest 和 gbest 不再变化交换序的差异趋近于零算法进入死循环。处理方式是在主循环内部加一个停滞计数器如果连续 N 代通常 2050 代gbest 没有改善就触发重启机制。if stall_count 30 % 重新初始化除 gbest 外的所有粒子 for i 1:nPop if i ~ gbest_index route{i} randperm(nCity); end end stall_count 0; end这种重启策略的效果在中等规模 TSP 上非常明显它能破坏粒子群的同质性为搜索引入新的多样性。需要注意的是重启后 pbest 也要重新评估但 gbest 要保留——那可能是目前找到的最好解不应该被丢弃。这是混合粒子群和标准 PSO 在工程实现上一个比较关键的处理差异。4. 混合策略与调参遗传交叉、局部搜索和参数组合4.1 选择、交叉、变异如何嵌入 PSO 主循环混合粒子群算法之所以叫“混合”核心在于把遗传算法GA的操作嵌入 PSO 的循环中而不是把两者前后串联跑两遍。常见做法是每次迭代在速度更新之后按一定概率对粒子的路径执行顺序交叉Order Crossover, OX和逆转变异Inversion Mutation。% 顺序交叉示例取两个位置之间的片段保持相对顺序填充 function child order_crossover(p1, p2) n length(p1); child zeros(1, n); idx sort(randperm(n, 2)); % 随机选两个交叉点 child(idx(1):idx(2)) p1(idx(1):idx(2)); rest p2(~ismember(p2, child(idx(1):idx(2)))); child([1:idx(1)-1, idx(2)1:end]) rest; end这段代码实现的是经典 OX先从父本 p1 中拷贝一段连续片段到子代再从父本 p2 中按顺序填入其余未出现的城市。OX 的好处有两个一是保证子代是一个合法排列不会出现重复或缺失城市二是尽可能保留父本中的相对顺序信息这对路径类问题很重要因为路径的价值不在于某两个城市是不是相邻位置而在于它们的访问先后关系。选择父本时常见做法是以 pbest 和当前粒子作为父本或以 gbest 和当前粒子作为父本这样交叉结果自然融合了全局最优信息与 PSO 的更新方向一致。变异操作则简单得多随机选择路径中的两个位置将中间的段反转。TSP 中经常出现交叉边和自交叉逆转变异能快速消除这种局部混乱且不破坏路径的合法性。嵌入位置放在速度更新之后、适应度评估之前这样每一代的搜索顺序是速度更新 → 交叉 → 变异 → 评估适应度 → 更新 pbest/gbest。这个顺序背后的逻辑是速度更新负责宏观靠拢向已知最优地区移动交叉负责组合已有优良片段变异负责局部扰动跳出局部极值。4.2 参数怎么设惯性权重、学习因子、交叉概率的配合参数之间是相互牵制的。只是简单把 c1 和 c2 调大并不能提升性能关键是让“速度更新产生的交换量”和“遗传操作的数量级”相匹配。下面是一组在 51 城市实例上表现稳定的参数组合参数推荐值作用说明惯性权重 w0.5 ~ 0.9 线性递减前期保持探索后期收敛个体学习因子 c11.2 ~ 1.8控制向 pbest 靠拢的交换量全局学习因子 c21.2 ~ 1.8控制向 gbest 靠拢的交换量交叉概率0.7 ~ 0.9路径片段重组的频率变异概率0.1 ~ 0.3扰动强度过大易破坏好解停滞重启阈值20 ~ 50 代过早重启会打断收敛过晚浪费计算我一般先固定交叉概率 0.8、变异概率 0.2只调 w、c1、c2 三个参数在这组基础上如果仍陷入局部最优再逐步提高变异概率到 0.3。需要特别强调的是变异概率不是越大越好。在 51 城市上变异概率超过 0.4 时粒子几乎变成了随机搜索收敛精度反而下降。对比标准 PSO 和 GA 在 eil51 上的表现可以发现HPSO 的收敛速度通常快于纯 PSO因为交叉引入了组合信息最终解质量优于纯 GA因为 PSO 框架提供了更明确的收敛方向。4.3 常见失败模式收敛太快、不收敛、运算过慢HPSO 跑起来之后常见的失败模式有三种现象和应对措施完全不一样。第一种是前 30 代内 gbest 迅速降到某个值就一动不动这是典型早熟收敛。原因通常是 c2 过大或粒子数太少粒子过早全部涌向 gbest丧失了多样性。应对方式是调大变异概率、提高停滞重启的触发频率或者引入 2-opt 局部搜索。第二种是迭代快结束 gbest 还在明显下降说明收敛不足需要增加迭代次数或改用 w 线性递减让后期更专注于局部精修。第三种是计算速度慢如果 51 城市、40 个粒子跑到 500 代需要几十秒通常是距离矩阵在循环里被重复计算了——每次适应度评估都调用坐标算距离而不是用预先算好的 D 矩阵。效率优化的思路很简单把 dist.m 的结果缓存下来能预计算的绝不进循环。对应地在 main.m 里可以做一个简单的加速优化把适应度评估、交叉、变异全部按粒子循环实现为一个独立函数传入 D 矩阵避免主脚本里重复读取坐标文件。这种改动不改变算法逻辑但比用全局变量或反复 load 数据要快数倍。检查运行时间可以用 Matlab 自带的tic/toc包住主循环输出每次迭代的平均耗时大于 0.1 秒/代就要检查是否出现了低效循环。5. 用 eil51 标准实例验证 HPSO 的收敛质量与稳定性5.1 基准实例与评价指标eil51 是 TSPLIB 中最常用的中等规模基准实例之一最优解为 426四舍五入后。验证 HPSO 是否真正有效不能只看一次运行结果因为启发式算法每次运行的初始粒子都是随机的。正确做法是独立重复运行 10 次统计每次找到的最优路径长度然后计算三个指标最佳值10 次中的最小值、平均值反映算法的期望表现和标准差反映算法的稳定性。最佳值接近 426 说明算法有找到全局最优的潜力平均值与最佳值的差距说明收敛一致性标准差大则说明算法输出不可控需要调参或增加迭代次数。5.2 跑批对比和可视化判断在 Matlab 中可以用一个简单脚本完成 10 次重复实验并绘制收敛曲线的平均值best_history zeros(10, MaxIter); for run 1:10 [gbest_path, gbest_val, history] HPSO_TSP(eil51_coords, nPop, MaxIter); best_history(run, :) history; end mean_curve mean(best_history, 1); plot(mean_curve); xlabel(Iteration); ylabel(Best Distance);history在主循环里记录每一代的 gbest 值最终得到的是 10 条收敛曲线的平均趋势。看这条曲线时重点关注两点曲线下降速度是否在前 100 代内基本进入平台期平台期的均值是否在 430450 区间内。如果在 50 代前就平台化且均值高于 500说明参数组合有严重问题优先检查 c2 是否过大、变异概率是否过低。将 HPSO 与标准 PSO如果实现了和遗传算法跑同一组数据对比时用箱线图展示 10 次实验的分布最直观HPSO 的中位数应明显低于标准 PSO且箱体宽度更窄。这里要提醒一下TSPLIB 中 eil51 的 426 是整数结果但 dist.m 如果不做四舍五入实际最优可能是 425.2 左右所以判断收敛质量时建议以 430 以下作为“合格线”不必强行追求 426。5.3 一个压箱底的小技巧保存历史最优路径用于路径重建最后一公里的检查往往是路径是否合法。很多人在 HPSO 结束后只输出gbest_val和gbest_path但这两个变量只能证明“找到了这个排列”无法直接在图上做几何验证。比较好的做法是在主循环里把每代的 gbest 路径也保存下来运行完后画出来观察是否出现交叉边。TSP 最优解的一个必要条件是路径在几何上不自交在欧氏距离下所以一旦可视化图中出现明显的交叉边说明当前解还没有收敛到足够好——这在几何上很直观可以用来快速判断算法在某次运行中是否失败。figure; plot(coords(gbest_path, 1), coords(gbest_path, 2), o-, LineWidth, 1.5); hold on; plot([coords(gbest_path(1),1), coords(gbest_path(end),1)], ... [coords(gbest_path(1),2), coords(gbest_path(end),2)], r-, LineWidth, 1.5); title(sprintf(Best Route Length: %.2f, gbest_val));可视化验证结束后如果确认解质量稳定可以将gbest_path和gbest_val保存到.mat文件供后续使用save(hpso_eil51_result.mat, gbest_path, gbest_val)。这样跑批时可以把 10 次实验的最优结果统一收集再用sortrows按路径长度排序快速定位哪次运行的gbest_path对应最小值。对于 eil51 这类中小规模实例HPSO 找到 430 以内解是合理的期望如果调参后仍达不到可以进一步叠加 2-opt 局部搜索在每次迭代结束时对 gbest 路径做边交换优化——这也是“混合”的延伸通常能把最后 23 个单位的差距补上。本文还有配套的精品资源点击获取