粒子群算法PSO跑TSP这种组合优化问题在Matlab里实现起来其实比很多人想象中简单但难点不在“写代码”而在“怎么把连续优化的那一套逻辑翻译成排列问题能用的版本”。很多教程一上来就甩一大段代码代码能跑可你想改一个参数都不知道从哪里下手更别提自己复现和改进。这篇文章就拿旅行商问题TSP当靶子给出一套完整的Matlab粒子群算法实现把每个模块的来龙去脉都拆开讲清楚包括位置怎么编码、速度怎么定义、交换序列怎么生成、主循环怎么更新。代码全部用Matlab基础函数写不需要优化工具箱适合课程设计、建模比赛、毕业设计里做智能算法对比实验的读者直接抄作业也适合想真正搞懂PSO离散化原理的人逐行研究。先说一个很多人问的问题标准粒子群算法的位置和速度都是连续实数更新公式就是“位置加速度”可TSP的解是一个城市的排列比如1→5→3→2→4→1你总不能把5.7个城市算进去吧。所以要做的事情其实是重新定义三样东西什么叫两个解之间的“差”什么叫“速度”什么叫“位置加上速度”。我下面用的方案是交换序列法这也是目前解决排列类问题最常用的离散PSO做法之一。理解了这一层后面看代码就是顺水推舟的事。1. 为什么TSP适合用粒子群来实验1.1 TSP的难度到底在哪里TSP本身定义非常简单给定N个城市坐标找一条从某个城市出发、经过所有城市恰好一次、最后回到出发城市的最短闭合回路。难点在于城市数量一上来暴力枚举根本扛不住。10个城市的排列是10!也就是3628800种听起来还凑合到了20个城市就是2.43e18种搜索空间超级计算机也顶不住。而这个问题的复杂度本质上是组合爆炸属于NP-hard级别的优化问题。所以在实际工程和研究里正常的思路不是找绝对最优解而是在合理时间内找到一个“足够好”的解。这就给各种启发式算法提供了用武之地粒子群算法就是其中之一。它不像穷举那样把整个搜索空间翻一遍而是通过一群粒子彼此协作、互相借鉴经验在解空间里不断朝着更优的区域移动。虽然不保证全局最优但通常能在很短时间内收敛到接近最优的路径这对工程上完全够用。1.2 连续版PSO为什么不能直接搬过来标准PSO里每个粒子有两个核心属性位置X和速度V。更新公式是V_new w * V c1 * r1 * (pbest - X) c2 * r2 * (gbest - X)X_new X V_new关键在于pbest - X和gbest - X这两个表达式。在连续空间里两个向量相减得到的是另一个向量表示“我离个体最优还差多少、离全局最优还差多少”然后速度带着粒子往这个方向飞。但TSP的解是一个排列两个排列做减法没有数学意义1→2→3和3→2→1相减等于多少没人知道。所以必须换一种语言来刻画“差异”。这就是交换序列出场的时机。两个排列之间的差完全可以定义为一组“交换对”的集合我交换哪两个位置上的城市就能从当前排列一步一步变成目标排列。比如当前路径是[1 2 3 4 5]目标路径是[3 2 1 4 5]那一个可行的交换序列就是“交换位置1和位置3上的城市”这一对操作就能完成转换。理解了这个离散PSO的核心思路就清晰了把pbest - X看成“把当前排列变成个体历史最优所需的交换对集合”把gbest - X看成“把当前排列变成全局最优所需的交换对集合”。速度就是一组待执行的交换对位置更新就是把这些交换对依次施加到当前排列上。1.3 一个形象的类比快递员的路线修正为了让你更有体感可以把粒子想象成一个跑固定区域的快递员。他每天送完货会记下一条自己走过的最短路线这就是pbest。整个配送站有一个全公司最好的路线时不时张贴在墙上这就是gbest。第二天出发前他心里琢磨三件事昨天那条路线我已经跑熟了保留一部分直觉惯性权重我自己的历史最优路线里有一段很顺值得改过来个体学习公司墙上的最优路线里有一段明显更好也值得借鉴全局学习。他不会把整条路线推翻重来而是在当前路线基础上把几个路段的顺序交换调整一下。这个“调整路段”的动作在代码里就是交换对一系列交换对拼起来就是粒子的速度。所以离散PSO的本质并没有偏离标准PSO只是把“矢量加减”换成了“排列之间的交换距离”算法骨架完全一样。2. Matlab完整代码分模块解析下面这段代码可以直接复制到Matlab里运行建议保存成tsp_pso.m脚本。我先把完整代码整体放出来然后一个模块一个模块拆开讲它干了什么以及每个关键变量到底在表达什么。clear; clc; close all; rng(1); % 1. 数据准备随机生成31个城市坐标 numCity 31; city 100 * rand(numCity, 2); dist zeros(numCity, numCity); for i 1:numCity for j 1:numCity dist(i, j) sqrt((city(i,1)-city(j,1))^2 (city(i,2)-city(j,2))^2); end end % 2. 粒子群参数设置 nPop 60; % 种群大小 maxIter 300; % 迭代次数 w 0.9; % 初始惯性权重 wEnd 0.4; % 结束惯性权重 c1 1.2; % 个体学习概率 c2 1.6; % 全局学习概率 % 3. 定义粒子结构体并初始化 empty_particle.Position []; empty_particle.Velocity []; empty_particle.Cost []; empty_particle.Pbest []; empty_particle.PbestCost []; particle repmat(empty_particle, nPop, 1); gbest.Position []; gbest.Cost inf; for i 1:nPop particle(i).Position randperm(numCity); particle(i).Velocity []; particle(i).Cost tourLength(particle(i).Position, dist); particle(i).Pbest particle(i).Position; particle(i).PbestCost particle(i).Cost; if particle(i).Cost gbest.Cost gbest.Position particle(i).Position; gbest.Cost particle(i).Cost; end end % 4. 主循环 bestHistory zeros(maxIter, 1); for it 1:maxIter w w - (it / maxIter) * (w - wEnd); for i 1:nPop seqPB calcSwapSeq(particle(i).Pbest, particle(i).Position); seqGB calcSwapSeq(gbest.Position, particle(i).Position); newVel []; if ~isempty(particle(i).Velocity) keepIdx rand(size(particle(i).Velocity, 1), 1) w; newVel particle(i).Velocity(keepIdx, :); end if ~isempty(seqPB) keepIdx rand(size(seqPB, 1), 1) c1; newVel [newVel; seqPB(keepIdx, :)]; end if ~isempty(seqGB) keepIdx rand(size(seqGB, 1), 1) c2; newVel [newVel; seqGB(keepIdx, :)]; end % 限制速度长度避免路径被过度打乱 if size(newVel, 1) floor(numCity / 2) newVel newVel(1:floor(numCity / 2), :); end pos particle(i).Position; for s 1:size(newVel, 1) pos swapPositions(pos, newVel(s, 1), newVel(s, 2)); end particle(i).Position pos; particle(i).Velocity newVel; particle(i).Cost tourLength(particle(i).Position, dist); if particle(i).Cost particle(i).PbestCost particle(i).Pbest particle(i).Position; particle(i).PbestCost particle(i).Cost; end if particle(i).Cost gbest.Cost gbest.Position particle(i).Position; gbest.Cost particle(i).Cost; end end bestHistory(it) gbest.Cost; fprintf(Iter %d, best cost %.2f\n, it, gbest.Cost); end % 5. 结果可视化 subplot(1, 2, 1); plot(city(gbest.Position, 1), city(gbest.Position, 2), o-, LineWidth, 1.5); title(Best TSP Tour); grid on; subplot(1, 2, 2); semilogy(bestHistory, LineWidth, 1.5); xlabel(Iteration); ylabel(Best Cost); title(Convergence Curve); grid on;还需要两个辅助函数保存到同一个目录下即可。function L tourLength(route, dist) L 0; n length(route); for i 1:n-1 L L dist(route(i), route(i1)); end L L dist(route(n), route(1)); endfunction seq calcSwapSeq(target, current) seq []; cur current; n length(cur); for j 1:n if cur(j) ~ target(j) idx find(cur target(j), 1); seq [seq; j, idx]; cur([j, idx]) cur([idx, j]); end end endfunction newPos swapPositions(pos, i, j) newPos pos; newPos([i, j]) pos([j, i]); end2.1 数据准备距离矩阵为什么要提前算好代码里先随机生成了31个城市的坐标这一步直接决定后面实验的可复现性。我用rng(1)固定了随机种子这是很多人写Matlab算法容易忽略的一个小细节。如果不固定你每次运行都会得到完全不同的城市分布那对比实验结果就没什么意义了因为你不知道这个结果到底是因为算法好还是因为你这次随机到一组特别简单的点。固定种子之后大家可以一起复现同一个算例调试和讨论都方便很多。距离矩阵dist是整个算法里使用频率最高的数据每个粒子每轮迭代都要计算路径长度而路径长度本质上就是一堆距离之和。如果每算一次都临时去开平方、求距离几百次迭代下来会浪费大量时间。所以正确做法是一开始就把所有城市两两之间的距离算好存成一个N×N的矩阵后面查表就行。dist(i,j)表示城市i到城市j的欧氏距离主对角线上是0。2.2 粒子初始化和路径长度计算每个粒子的Position就是一条路线用randperm(numCity)生成一个1到31的随机排列。注意这个排列是有方向性的[1 2 3]和[2 1 3]在路径长度上不一定相等所以位置实际上还隐含了起点。不过TSP的闭合特性意味着任何一个城市都可以作为起点所以dist(route(n), route(1))把最后一个城市和第一个城市连接起来保证路径是闭合的。tourLength函数就是干这个的把路径上相邻城市的距离累加起来再加上最后一个城市回到第一个城市的距离。这里有个容易出错的地方很多人第一次写会漏掉最后一段回程导致适应度值偏小。我在调代码的时候踩过这个坑路径看起来没问题收敛曲线也正常但算出来的“最短路径”其实就是一条不闭合的折线最后补上回程那一段才对了。初始化的另一层意义是为全局最优gbest提供一个起点。我先把gbest.Cost设成inf然后遍历所有粒子只要发现某个粒子的Cost小于当前gbest.Cost就更新gbest。这样一轮下来gbest保证是初始种群里的最优路线。2.3 交换序列函数两个排列的“差”怎么求calcSwapSeq这个函数是整个离散PSO的灵魂。传入两个参数target和current返回一个两列的矩阵seq每一行表示一次交换操作第一列和第二列是要交换的两个位置。核心逻辑并不复杂从左到右扫一遍current如果第j个位置上的城市和target不一样就找到target(j)这个城市在current的哪个位置把这两个位置上的城市交换然后记录下这次交换。举个例子target是[3 2 1 4 5]current是[1 2 3 4 5]。第一轮j1cur(1)1不等于target(1)3在cur里找到3的位置是3记录交换对(1,3)交换以后cur变成[3 2 1 4 5]继续往后扫2、1、4、5全部对上了结束。所以返回的交换序列只有一行。把这个交换对作用到最初的current上就能得到target。这就是“两个排列之间的差”的数学表达。这个函数的时间复杂度在最坏情况下是O(n^2)因为find每一次都是线性搜索。31、50个城市完全没感觉但如果做到几百个城市这个函数会变成性能瓶颈可以考虑用映射表记录每个城市当前所在的位置把复杂度降到O(n)。文章后面还会再提这个问题。2.4 速度更新惯性、个体学习、全局学习如何体现主循环里最重要的部分就是速度的合成。我把速度定义成一个交换对列表所以三部分贡献都体现在“往这个列表里追加交换对”上。保留旧速度的操作是取出原有Velocity按照概率w随机保留部分交换对pbest方向的贡献是取seqPB中按概率c1保留的行gbest方向的贡献是取seqGB中按概率c2保留的行。三个来源合并后就得到当前粒子的新速度。看到这里可能有读者想问c1和c2在标准PSO里是加速系数为什么到了这里变成了概率这确实是离散PSO实现里最常见的变体。严格来说标准的交换序列式PSO还可以给每个交换对分配一个“速度收益搜索算子的保留概率”但为了代码简洁和便于理解我更倾向于直接用概率采样来决定“这一批交换对是否执行”。这样做的好处是逻辑直观缺点是对c1、c2数值的敏感度比连续PSO高一点。后面参数调优部分还会细说一般c1取0.8到1.5之间c2取1.2到2.0之间运行效果都比较稳。速度生成了以后还要执行“位置加速度”这一步。我把当前位置pos复制一份然后依次把newVel里的每一行交换对作用到路上。这里要特别强调一下交换对里的两个数字是位置下标不是城市编号。比如交换对是(1,5)意思是把当前路径的第1个位置和第5个位置上的城市互换。在实际调试中我见过不少新手把城市编号当位置下标结果路径里出现重复城市或者整个路径顺序越变越乱最后收敛曲线直接崩掉。为了防止粒子速度里的交换对过多把路径拆得七零八落我在更新后加了一个限制如果交换对数量超过城市数量的一半就截断到floor(numCity/2)。这个操作是我做实验的时候发现很有必要的因为三部分速度合并之后交换对数量有时候会膨胀得很厉害有些粒子一轮就要执行十几次交换等于把整条路线完全洗牌反而找不到更优解。加入这个上限后粒子更新的幅度更有节制收敛也更稳。2.5 收敛曲线的可视化与结果判断绘图部分我用了两个子图。左边把gbest.Position按顺序连接起来画出最终最优路径。这里有个绘图函数的小技巧plot(city(gbest.Position, 1), city(gbest.Position, 2), o-)里用gbest.Position作为下标来索引城市坐标这样Matlab会自动按照路径顺序连接各个点不需要手动去重排坐标矩阵。路径是不是合理、有没有交叉边一眼就能看出来。右边画收敛曲线我用的是semilogy纵轴取对数。因为路径长度从初始值几千一路降到几百如果用线性坐标前面几轮的大数值会把后面的细节压得看不清用对数坐标能更清楚看到后期收敛过程。如果你发现收敛曲线是一条斜线持续往下走说明算法还在稳定改进可以放心加大迭代次数。如果曲线很快变成一条平线大概率是粒子群已经聚集到了某个局部最优。bestHistory这个数组是调试算法最重要的工具我在主循环里每轮记录一次gbest.Cost。很多同学写算法只看最终结果不看中间过程。但实际调参的时候曲线形状能告诉你很多信息到底是收敛太慢、还是收敛太快、还是根本不收敛。后面第4节问题排查全都是围绕这条曲线展开的。3. 参数怎么调从原理到经验值3.1 参数速查表粒子群算法在TSP上的表现参数影响非常直接。我整理了一张表对应代码里的变量名给出推荐范围和我的实测感受。参数常用范围作用注意事项nPop40-100种群规模决定并行搜索的广度城市数越多种群越大但太大计算量猛增maxIter200-1000迭代轮数先用小迭代数跑通再放大量级w0.9→0.4惯性权重保留旧速度交换对的比例建议线性递减不要固定在一个值c10.8-1.5向个体最优学习的概率太大会让每个粒子沉迷自身路径c21.2-2.0向全局最优学习的概率太大容易早熟收敛到局部最优3.2 为什么惯性权重要线性递减代码里我把w从0.9线性递减到0.4这是很多文献里都验证过的经验做法。早期粒子需要保持较高的运动活性尽可能覆盖更广的搜索区域所以保留更多旧速度让路线调整的幅度大一些。到了后期粒子群已经摸到比较有希望的区域了这时候如果还在大幅度乱跳会破坏好不容易积累下来的结构所以要把交换对保留比例降下来让粒子把更多精力放在精细打磨局部路径上。对比实验也能明显看到差距如果w固定为0.9整个算法后期收敛非常慢最优解的路径图上偶尔还会残留交叉边如果w固定为0.4前期搜索能力太弱粒子群很快就抱成一团最后解的质量很差。线性递减实际上是“先广后精”的策略在TSP这种搜索空间巨大的问题里尤其适用。3.3 加速系数对搜索行为的影响c1和c2在我的实现里表现为“保留交换对的概率”所以它们对搜索行为的影响非常直观。c1越大粒子越倾向于把自己的路线往pbest方向改c2越大粒子越倾向于往gbest方向改。理想情况下两者应该平衡让粒子既借鉴全局最优又保留自己的个性。我个人的调参习惯是c2比c1略大比如c11.2、c21.6因为TSP问题中存在大量局部最优全局最优信息相对更有引导价值。但也不能差太多。如果你的收敛曲线在第50轮之前就趴平了而且最终路径图上交叉边很多基本可以断定是c2设置过大粒子群被某个局部最优强吸引住了。这时候把c2降到1.2左右把c1提到1.5左右让粒子多探索一些不同路径往往能救回来。调参的时候一定要记住一个原则一次只改一个参数。很多同学喜欢同时动w、c1、c2和种群大小结果效果变差了根本不知道是谁的锅。每次只调一个变量记录下收敛曲线和解的质量再决定下一步怎么走这才是科学的实验方法。4. 常见问题与排查实录4.1 收敛曲线不动了怎么办这是最常遇到的问题也是新手最容易误判的情况。我做实验的时候一开始就遇到过前50轮曲线快速下降后面200轮纹丝不动仿佛算法已经收敛了。但仔细看路径图还有一些明显的交叉边这说明算法陷入了局部最优。如果你遇到这种情况首先检查收敛曲线的纵轴是不是到了对数坐标下非常小的量级。如果路径长度还在几百的量级说明改进空间还很大需要调整参数。有效的手段有三个一是把惯性权重下限wEnd调高一点比如从0.4调到0.5让后期粒子还保留一定活跃度二是把c1调大一些让粒子不要完全被全局最优裹挟三是在更新速度时加入“变异机制”也就是在newVel里随机插入若干组随机生成的交换对模仿遗传算法的变异操作增加跳出局部最优的概率。这部分我写得很直接因为我自己就是靠第三个方法在31城算例上把结果又推进了一截。4.2 路径上出现重复城市怎么定位如果你在调试代码的时候发现某条路径里有重复城市同时少了某个城市那问题几乎可以锁定在交换操作的实现上。对应到代码里要么是swapPositions这个函数的参数传错了要么是calcSwapSeq里find(cur target(j), 1)没找到目标城市。后者在正常逻辑下不会发生因为target和cur是由同一组城市排列构成的每个城市必然存在。前者就很容易出错具体表现是你在计算newVel时把城市编号当成了位置下标结果交换的是“城市编号为1和城市编号为5”的路径点而不是“第1个位置和第5个位置”。定位这种问题有个特别有效的小实验把种群大小设为1迭代次数设为1把粒子初始化成一条已知路径比如1到N的顺序排列然后在主循环里打印更新前后Position。如果更新后路径漏了城市或者多了城市直接用脚本看是哪一行代码造成了破坏。这种单步调试比盯着代码看半天快得多。我写第一版交换序列PSO的时候就是这么把问题抓出来的。4.3 TSP、VRP和静态欧式TSP别搞混了搜索热词里同时出现了TSP和VRP说明这两个概念确实容易被放一起讨论顺手澄清一下。TSP是最经典的单一旅行商回路问题一辆车出发、经过所有客户点、回到起点没有容量限制不需要考虑多车协同。VRP是全称“车辆路径问题”本质是TSP的推广多辆车共享一组客户点每辆车都有容量上限甚至时间窗约束目标是让总行驶距离最短或者使用的车辆数最少。VRP和TSP的算法思路有很大交集但每次粒子更新一条路线已经不够了必须同时维护多条路线并且处理客户点分配和路径顺序两重优化复杂度上了一个台阶。至于静态欧式TSP指的是所有城市坐标固定不变距离用欧氏距离计算。对应的是动态TSP城市坐标会随时间变化或者路网距离不是直线距离。写论文或者做报告的时候这些前提条件一定要写清楚否则别人复现你的实验根本对不上结果。本文代码默认就是静态欧式TSP这也是最经典的入门假设。4.4 规模扩大之后要注意什么31个城市跑300次迭代在Matlab里基本上几秒到十几秒能跑完完全没有压力。但是如果你把城市数量提高到100、200情况就不一样了。这个时候有两个地方会成为瓶颈一是每次calcSwapSeq都要O(n^2)时间扫描二是每轮迭代每个粒子都要把速度里的交换对全部执行一遍速度累计起来很可观。我的建议是把种群大小适当提高到80到120迭代次数可以保持300到500因为大规模TSP靠一味增加迭代次数收益会越来越小。同时可以考虑在calcSwapSeq里加入一个“城市位置映射表”比如用一个长度为n的数组posOfCity保存每个城市当前所在位置交换城市时同步更新映射表查找目标城市位置的时间就从O(n)降到O(1)。这个优化写起来也不复杂但对大规模算例帮助很明显。再进一步可以在每轮迭代后对gbest执行一次2-opt局部搜索把交叉的边解掉再更新历史最优这就是混合算法了效果比单纯调参数提升更快。5. 扩展想法和一点私人经验代码本身到这里已经能完整跑通了。如果你只想交作业到这里就够了。但如果你是想把这个实验写进论文或者作为毕设的一部分我建议再往前走两步。一是大胆引入2-opt局部搜索。TSP问题天然适合2-opt把路径里的两条边断开然后反向重连如果总长度变短就接受。这个操作实现起来也就二三十行Matlab代码但效果立竿见影。把局部搜索加在gbest更新之后每一轮对全局最优做几次2-opt扰动可以明显弥补PSO局部搜索能力弱的短板。很多论文里说的“混合粒子群算法”其实就是在PSO框架里嵌入这类局部搜索算子。二是跑多个随机种子做统计分析。很多课设报告只用一组随机种子跑一遍然后拍脑袋说“算法收敛到了XX”这种结论其实站不住脚。正确做法是选5个或10个不同的随机种子每个都固定下来分别跑完整算法记录每次的最优路径长度、收敛代数、运行时间最后算平均值和标准差。这样写进报告里才是有说服力的实验对比。说回我个人实操的体会。粒子群算法在TSP上的表现很依赖参数和随机种子不同算例下同一组参数可能一个效果好一个效果差。所以我写这类算法实验的时候一定会把随机种子、参数配置、城市规模这些信息全部记录在注释里保证任何一次实验都能复现。这个习惯后来帮我排除过很多“这次怎么跑出来的结果跟上次不一样”的困惑。调参的时候我习惯先把w的初值和终值定下来再调c1和c2的比例关系最后再动种群大小和迭代次数这样每次只动一个变量问题定位快很多。最后再分享一个小技巧调试阶段把城市数量改成8个左右迭代次数改成100这样整个程序一秒内就能跑完你就能在电脑前反复折腾交换序列、速度限制这些细节等逻辑完全跑通之后再改回大算例正式实验效率高得多心也不会累。