粒子群算法到底是只能解连续优化还是也能啃组合优化问题这个问题困扰了我挺长时间。直到我拿Matlab把粒子群跑在旅行商问题TSP上才发现思路一旦打开代码量甚至比遗传算法还少效果也相当能打。这篇就把整套实现过程、完整代码和我在调参过程中踩过的坑一次讲清楚适合正在学智能算法、准备课程大作业或者工作中遇到路径排程问题、想用Matlab快速验证算法的朋友。我当时面对TSP的第一个想法是粒子群更新公式里全是连续实数的加减乘而TSP要输出一个不重复的城市访问顺序这两者怎么对上后来我用的是“随机键编码”也就是让粒子继续在连续空间飞解码时靠排序得到路径。整套流程跑通之后你会发现所谓离散优化很多时候只是换了个编码方式粒子群的核心迭代机制一点不用改。1. 粒子群算法与TSP的适配逻辑打破“连续优化只解连续题”的刻板印象1.1 粒子群算法的核心机制三股力如何推动搜索粒子群算法的思想其实非常朴素你想象一群鸟在一片区域里找食物每只鸟都不知道食物在哪但它们能记住自己经过的最好位置也能看到群体里其他鸟发出信号于是每只鸟下一次飞行方向由三股力共同决定。第一股力是“惯性”也就是它当前飞行的速度代表了对上一步状态的保持第二股力是“个体认知”它会被自己历史最优位置吸引代表个体经验第三股力是“社会认知”它会被整个群体当前最优位置吸引代表群体分享信息。三股力加权合成后粒子就完成了从当前位置到下一位置的移动。用公式写就是v(i) w * v(i) c1 * r1 * (pbest(i) - x(i)) c2 * r2 * (gbest - x(i)) x(i) x(i) v(i)其中w是惯性权重c1和c2是学习因子r1和r2是[0,1]之间的随机数。这个公式是粒子群算法的灵魂也是后面所有代码的主干。理解了这个机制再看TSP问题你要思考的就只有一个点粒子位置x到底表示什么。1.2 TSP为什么难组合爆炸的问题描述TSP问题的描述看起来极短给定n个城市和两两之间的距离找一条经过每个城市恰好一次、最后回到起点的最短回路。听起来简单但它是典型的NP-hard问题。为什么难因为城市访问顺序的候选数量是(n-1)!/2这是个阶乘级的天文数字。n20的时候理论上要检查的路线数量已经是4万亿实际上是约1.2万亿哪怕计算机一秒钟检查100万条路线也要跑好几年。n50或者n100的时候这个数字直接把穷举判了死刑。精确算法当然存在比如动态规划可以做到O(n²·2ⁿ)但n超过30之后内存和时间都会爆炸。因此实际工程和学术研究中大家更常用的是启发式算法或元启发式算法粒子群、遗传、模拟退火、蚁群都属于后者。粒子群的特点是群体并行搜索、参数少、实现简单再加上粒子记住个体最优和全局最优的特性天然适合去逼近这类大规模组合优化问题的近似最优解。1.3 “连续粒子”如何变成“一条路径”随机键编码思想这里就是全文最关键的地方。TSP的解是“城市排列”本质上是离散的、不重复的序列粒子群默认的位置x是连续实数向量直接套用会产生非法解。解决办法就是随机键编码。随机键编码的做法非常聪明粒子位置x保持一个n维连续实数向量解码时对这个向量做一次排序排序后的索引就作为城市的访问顺序。举个例子。假设有5个城市一个粒子的位置是x [0.28, 0.91, 0.05, 0.62, 0.33]按值从小到大排序对应索引是[3, 1, 5, 4, 2]那么这条路径就是先访问城市3再访问城市1再访问城市5再访问城市4最后访问城市2然后回到城市3。这个索引序列永远不会重复因为它就是对1到n的排列天然合法。随机键编码的好处在于PSO的连续空间更新公式完全不用改粒子的位置还是老老实实的实数飞行也在连续空间进行唯一多出来的步骤就是每次评估之前先sort一下解码成路径。这就把离散优化问题转换成了连续空间上的搜索问题逻辑上完全自洽。有人可能会问那为什么不直接用交换序列、把速度定义成“一系列交换操作”来更新路径这个方案也存在但实现起来要处理交换序列的拼接、截断、去重复杂度高得多而且容易写出隐藏bug。随机键编码是“换一种空间”交换序列是“在离散空间硬建模”。新手我强烈推荐随机键编码先跑通再谈进阶。2. 开始编码前先把三张表搭好城市、距离矩阵与PSO参数2.1 城市坐标与距离矩阵不要用双重循环写的太笨我用20个城市的随机坐标来演示。城市数太少体现不出算法效果太多又会让演示迭代时间变长20是个刚好能看懂的规模。numCities 20; cityXY 100 * rand(numCities, 2);这样就生成了20个[0,100]×[0,100]范围内的城市坐标点。接下来要算距离矩阵distMatdistMat(i,j)表示城市i到城市j的欧氏距离。这里给一个中规中矩的写法distMat zeros(numCities, numCities); for i 1:numCities for j 1:numCities distMat(i, j) sqrt((cityXY(i,1) - cityXY(j,1))^2 ... (cityXY(i,2) - cityXY(j,2))^2); end end双重循环在n20时无所谓但如果城市数上到几百建议用向量化写法一行搞定d2 (cityXY(:,1) - cityXY(:,1)).^2 (cityXY(:,2) - cityXY(:,2)).^2; distMat sqrt(d2);两种写法效果完全一样选哪种取决于你是否在意性能。注意distMat必须是方阵对角线元素是0这关系到后面计算回路总长度时最后一步“回起点”的距离。2.2 参数怎么定粒子数、迭代次数、w与c1/c2的初值PSO有四组核心参数直接决定搜索效果我会在表格里列出来并说明调整逻辑。参数典型范围我的建议说明粒子数numParticles30~10050越多探索能力越强但每代计算量线性增加最大迭代次数maxIter100~500200小规模问题200代基本够大规模优先加迭代惯性权重w0.4~0.90.9大权重利于全局探索可线性递减实现先探索后精修学习因子c1、c21.5~2.52.0c1管个体经验c2管群体经验两者相等是保守做法惯性权重w是第一个值得细说的参数。w越大粒子越倾向维持当前速度不容易被个体和群体带偏全局探索更强w越小粒子越容易被最优位置吸引局部开发更强。常用的进阶做法是让w在迭代过程中从0.9线性衰减到0.4前期大范围找后期精细挖。基础版本里我固定用0.9代码更简单效果也过得去。c1和c2通常取2.0左右让个体认知和社会认知大致平衡。如果你发现算法很容易早熟可以把c2调小一点让粒子多摸索一会儿如果你发现收敛太慢可以适当加大c2让群体最优的吸引力更强。但调参有个基本原则一次只动一个参数记录结果别一上来就乱调。还有一个小细节容易被忽略Matlab的rand每次运行结果都不同这会让实验结果无法复现。教学演示时我建议先用rng(42)固定随机种子方便对照结果真正求解问题时再放开并在多次运行中取最优值或平均值。rng(42); % 固定随机种子让结果可复现 numParticles 50; maxIter 200; w 0.9; c1 2.0; c2 2.0;2.3 评价函数路径长度怎么算目标函数是路径总长度也就是按照路径顺序把相邻城市距离累加起来最后加上“从最后一个城市回到起点”的闭路距离。function totalDist evalRoute(route, distMat) n numel(route); totalDist 0; for i 1:n-1 totalDist totalDist distMat(route(i), route(i1)); end totalDist totalDist distMat(route(n), route(1)); end这段代码非常简单但恰恰是最容易写错的地方。实际调试中我看到很多人漏掉最后一句totalDist totalDist distMat(route(n), route(1))导致评价函数少算了一段回路距离算法却在朝着错误目标优化结果画出来的路径首尾不闭合。遇到结果不对劲时先用小规模数据手工验算评价函数比如4个城市去算一条已知路线的人工距离和函数输出对一下。3. 核心代码完整拆解跑通PSO-TSP的每一行都说明白3.1 主脚本骨架初始化、迭代、绘图我把完整主脚本放在这里然后逐段拆解因为很多细节不逐行说明你复制代码跑通是一回事真正理解又是另一回事。% pso_tsp_demo.m clear; clc; close all; % 1. 生成20个城市的随机坐标 numCities 20; cityXY 100 * rand(numCities, 2); % 2. 距离矩阵 distMat zeros(numCities, numCities); for i 1:numCities for j 1:numCities distMat(i, j) sqrt((cityXY(i,1) - cityXY(j,1))^2 ... (cityXY(i,2) - cityXY(j,2))^2); end end % 3. PSO参数 numParticles 50; maxIter 200; w 0.9; c1 2.0; c2 2.0; % 4. 初始化 particle rand(numParticles, numCities); velocity zeros(numParticles, numCities); pbest particle; pbestVal inf(numParticles, 1); gbestVal inf; for i 1:numParticles route decodeRoute(particle(i, :)); pbestVal(i) evalRoute(route, distMat); if pbestVal(i) gbestVal gbestVal pbestVal(i); gbest particle(i, :); gbestRoute route; end end bestHist zeros(maxIter, 1); % 5. 迭代 for iter 1:maxIter for i 1:numParticles r1 rand(1, numCities); r2 rand(1, numCities); velocity(i, :) w * velocity(i, :) ... c1 * r1 .* (pbest(i, :) - particle(i, :)) ... c2 * r2 .* (gbest - particle(i, :)); particle(i, :) particle(i, :) velocity(i, :); route decodeRoute(particle(i, :)); val evalRoute(route, distMat); if val pbestVal(i) pbestVal(i) val; pbest(i, :) particle(i, :); end if val gbestVal gbestVal val; gbest particle(i, :); gbestRoute route; end end bestHist(iter) gbestVal; end % 6. 绘图 figure(Position, [100 100 1000 420]); subplot(1, 2, 1); plot(cityXY(:, 1), cityXY(:, 2), ko, MarkerFaceColor, k); hold on; plot(cityXY(gbestRoute, 1), cityXY(gbestRoute, 2), r-, LineWidth, 1.6); title(最优路径); xlabel(x); ylabel(y); grid on; subplot(1, 2, 2); plot(bestHist, b-, LineWidth, 1.5); title(收敛曲线); xlabel(迭代次数); ylabel(最优路径长度); grid on;初始化部分最关键的一点particle rand(numParticles, numCities)生成的是连续实数矩阵也就是50个粒子每个粒子是20维实数向量。这些实数向量本身不是路径只是路径的“编码”。pbest和pbestVal分别保存每个粒子个体历史最优的位置和对应的路径长度gbest和gbestVal保存全局最优的位置和长度。迭代循环是核心。内层for i对每个粒子做三件事先按标准公式更新速度其中r1和r2是1×20的随机向量能让每个维度产生不同的随机扰动再按位置公式更新粒子位置然后解码、评估检查是否需要更新pbest和gbest。这里要注意particle(i,:)更新后有可能超出初始范围导致数值非常大或非常小但这在随机键编码下完全没问题因为sort只关心元素的相对大小不关心绝对数值。这点我第一次就想歪了总觉得粒子飞出边界是不是出bug了其实排序解码让这个问题自然消失了。3.2 解码与评估函数为什么简单到只有几行解码函数用了Matlab的sort返回值技巧。sort(x)返回两个值第一个是排序后的数组第二个是排序时各元素在原数组中的下标。我们要的就是这个下标function route decodeRoute(x) [~, route] sort(x); end这行代码是整个随机键方案的精髓也是我写完整块代码之后最想画红线的地方。sort之后得到的是1到n的一个排列它天然满足TSP“每个城市只能访问一次”的约束。不需要做去重不需要修复非法解这就是编码设计带来的红利。评估函数上一节已经给过。这两段加在一起不到20行但把“连续粒子”和“离散路径”两个世界完整对接上了。3.3 初始化陷阱为什么一开始不要用randperm有个非常自然的念头既然最终路径是城市排列那我初始化的时候直接用randperm生成一组随机排列省得再sort一次。这种想法我理解但一定不要这么干。原因在于PSO的速度-位置更新公式是在连续实数空间定义的。如果用randperm初始化粒子位置这个位置本身就是整数排列下一轮更新时velocity加上去之后粒子位置会变成带小数的向量比如[1.2, 3.7, 2.1, 4.9]你无法直接把它解释成一个合法路径。强行取整会出现重复城市造成非法解算法直接失灵。随机键编码之所以成立是因为粒子的连续位置完全不需要解释为路径路径只是排序的副产物。编码空间和解空间分离算法才能顺畅运行。所以初始化必须用rand不能用randperm。这个逻辑想通了整篇文章的核心也就抓住了。4. 跑通之后的结果判读与三个真实踩坑记录4.1 第一次运行收敛曲线和路线图怎么看代码跑通之后你会看到左右两张图。右边是收敛曲线横轴迭代次数纵轴最优路径长度。典型的收敛曲线是前30至50代快速下降后面逐渐平缓。如果曲线在中后期基本不再变化说明群体已经收敛如果直到最后一代还在明显下降说明maxIter给少了要继续跑。左边是路径图。最优路线会是一个闭合回路把20个城市都用线段连起来好的路线没有交叉形状比较“舒展”。如果出现很多交叉线段大概率就是还没收敛或者粒子数太少导致搜索能力不足这时候可以加迭代次数、加大粒子数或者引入局部搜索这个后面细说。还要验证一下结果合理性。你可以随机生成100条路径算平均长度和算法找到的长度做对比。粒子群找到的结果一般明显好于随机路径平均值如果在你的城市规模下结果和随机路径差不多那大概率参数有问题或者是代码某处写错了。4.2 坑一把粒子位置直接当路径速度更新把解变得“不合法”这个坑我踩得最惨发生在我第一次尝试“简化”的时候。我用randperm生成初始位置然后直接对位置向量做速度更新结果粒子位置变成了带小数的随机数解码时一堆重复城市得到的所谓路径根本不能算路径。后来我才意识到我混淆了“粒子空间”和“解空间”。粒子群更新永远发生在编码空间也就是连续实数空间TSP路径只是解码空间里的一个投影。想让算法工作你必须保证编码空间数学结构完整至于解码后的排列是升序还是降序、代表起点还是终点都只是约定。牢记这一点就能避开很多莫名其妙的bug。4.3 坑二距离矩阵写成非对称或漏掉回程距离另一个隐蔽的坑在距离矩阵。TSP的欧氏距离是对称的distMat(i,j)distMat(j,i)对角线为0。但有些人在手动构造数据时把distMat写成了只填对角线一侧的下三角矩阵或者坐标差算错符号导致距离矩阵不对称。PSO不会立刻报错它只会沿着错误距离优化最后给出的路线看上去“神秘地短”。排查方法很简单算一下distMat是否对称用isequal(distMat, distMat)检查再检查对角线是否全为0。另外evalRoute里漏掉回起点那段距离是最常见的路径长度会少算一段导致算法认为某个解很好实际画出来路径是断开的。每次修改评价函数后用小型人工数据验证一下整个过程不过十几秒但能省掉后面好几小时的困惑。4.4 坑三不看随机特性就断言结果越好粒子群是随机搜索算法同样的代码、不同的随机种子跑出来的结果是有波动的。我见过有人在对比实验里只跑一次就宣称A算法比B算法好这在方法学上是站不住脚的。一个稳妥的做法是独立运行20次记录每次的最优路径长度算平均值、标准差、最优值和最差值。平均值反映算法的一般水平标准差反映稳定性。如果你要发小论文或者做项目汇报可以顺便画个箱线图比单次结果有说服力得多。Matlab写个外层for循环跑20次就行把gbestVal记成数组再分析。numRuns 20; results zeros(numRuns, 1); for run 1:numRuns % 这里放入整个PSO流程结束后把gbestVal存入results(run) end disp([平均: , num2str(mean(results))]); disp([标准差: , num2str(std(results))]);这也是为什么教学演示时我推荐先rng(42)固定种子至少读者复制代码跑出来的数字是一样的能验证自己是否成功复现。5. 从能跑到跑得好给基础PSO加一个2-opt局部搜索5.1 基础PSO为什么在TSP上容易“差点意思”随机键编码的PSO有全局搜索能力这是它的优点但它也有一个明显短板局部精修能力弱。粒子在连续空间里“飞”排序解码后的路径在离散空间里的微小变化可能对应连续空间里一次不小的位置变动这种非线性会让算法很难对路径做精细调整。实战中常见的问题是最优路径基本成型但总有一两条边是交叉或次优的看起来就差那么一点。单纯加大迭代次数往往改善有限因为粒子群在收敛后期群体多样性严重下降大家挤在同一个局部最优附近谁也跳不出来。这时候最适合往算法里加一个局部搜索算子。行业通用做法是2-opt专门用来消除路径中的交叉边和次优段。2-opt的思想非常直观取路径中的一段把它反转如果反转后的总路径更短就接受这个新路径反复执行直到没有改进为止。这就像你出门一看地图发现两段路交叉了你直接把其中一段掉个头路径立刻顺很多。5.2 2-opt的Matlab实现与接入方式2-opt函数实现如下function route twoOpt(route, distMat) n numel(route); improved true; while improved improved false; for i 2:n-1 for j i1:n current distMat(route(i-1), route(i)) ... distMat(route(j), route(mod(j, n)1)); changed distMat(route(i-1), route(j)) ... distMat(route(i), route(mod(j, n)1)); if changed current - 1e-10 route(i:j) route(j:-1:i); improved true; end end end end end理解这个函数抓住一个关键点它只比较两段边的总长度。原始路径上第i-1个城市连接第i个城市第j个城市连接第j1个城市如果把i到j这一段反转那么变成第i-1个城市连接第j个城市第i个城市连接第j1个城市。如果后者总长更短反转就值得做。循环条件improved让整个过程持续到找不到任何改进操作为止。在主循环末尾接入2-opt记得更新全局最优bestHist(iter) gbestVal; newRoute twoOpt(gbestRoute, distMat); newVal evalRoute(newRoute, distMat); if newVal gbestVal gbestRoute newRoute; gbestVal newVal; end bestHist(iter) gbestVal;这里有个细节要说明一下2-opt之后gbestRoute和gbestVal变了但gbest这个连续位置向量没有同步改变。因为2-opt是在解空间做的局部精修我们并不需要把这个精修后的路径再映射回连续空间下次粒子更新的社会认知依然指向编码空间里的gbest位置。这样设计的好处是保持随机键空间的数学一致性缺点是gbestVal和gbest编码不再完美对应。实际不影响算法运行因为gbestVal只是用来输出最优值社会认知看的是gbest编码位置而不是gbestVal。我在代码注释里专门标了这一句防止后来看代码的人困惑。5.3 加了2-opt之后实验对比和使用建议在相同城市数据和相同参数下基础PSO和PSO2-opt的结果一般会有可感知的差距。以20个城市为例基础PSO跑200代最优路径可能在320到350之间加2-opt之后能压到300到320。城市数越多2-opt带来的改善越明显因为大规模TSP里交叉边和次优段出现的概率更高。不过也要提醒一句2-opt不是万能的它本质上是局部搜索只能把已有路径“捋直”不能像PSO那样在全局范围重新搜索。所以合理的搭配是PSO负责全局搜索、把解的“小区”找对2-opt负责局部精装修把小区里的“户型”调整到最好。这也是很多学术论文里“混合算法”的思路来源。如果你的数据规模上百个城市还可以考虑把2-opt作用到每个粒子的pbest上而不仅仅是作用于gbest。代价是每代计算量增加不少但搜索质量也会进一步提升。先用全局最优做2-opt如果不够再往每个粒子铺开这是我自己的调优顺序。个人体会是粒子群解TSP更像“搭桥”编码设计是桥墩参数调优是桥面2-opt这样的局部搜索就是护栏。桥墩打歪了后面再怎么修饰都白搭但桥面铺好了加上护栏才真正能通车。这篇代码跑通之后你可以尝试把惯性权重改成线性递减、把满足条件的2-opt接入粒子个体、或者把城市数加到50再观察算法表现。方向对了后面的路会顺很多。