粒子群算法求解TSP组合优化问题：Matlab实现与2-opt局部搜索优化

粒子群算法到底是只能解连续优化，还是也能啃组合优化问题？这个问题困扰了我挺长时间。直到我拿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，这是个阶乘级的天文数字。n=20的时候，理论上要检查的路线数量已经是4万亿（实际上是约1.2万亿），哪怕计算机一秒钟检查100万条路线，也要跑好几年。n=50或者n=100的时候，这个数字直接把穷举判了死刑。

精确算法当然存在，比如动态规划可以做到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]范围内的城市坐标点。接下来要算距离矩阵distMat，distMat(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

双重循环在n=20时无所谓，但如果城市数上到几百，建议用向量化写法一行搞定：

d2 = (cityXY(:,1) - cityXY(:,1)').^2 + (cityXY(:,2) - cityXY(:,2)').^2; distMat = sqrt(d2);

两种写法效果完全一样，选哪种取决于你是否在意性能。注意distMat必须是方阵，对角线元素是0，这关系到后面计算回路总长度时最后一步“回起点”的距离。

2.2 参数怎么定：粒子数、迭代次数、w与c1/c2的初值

PSO有四组核心参数直接决定搜索效果，我会在表格里列出来并说明调整逻辑。

惯性权重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(i+1)); 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 = i+1: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个城市连接第j+1个城市；如果把i到j这一段反转，那么变成第i-1个城市连接第j个城市，第i个城市连接第j+1个城市。如果后者总长更短，反转就值得做。循环条件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和PSO+2-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再观察算法表现。方向对了，后面的路会顺很多。