1. 项目概述为什么用粒子群算法求解TSP1.1 旅行商问题到底难在哪旅行商问题Traveling Salesman ProblemTSP是组合优化领域最经典的NP难问题之一。简单说就是一个商人要跑遍N个城市每个城市只能去一次最后回到出发城市怎么走总路程最短。这个问题听起来简单但复杂度随着城市数量增加呈现阶乘式爆炸。10个城市有约36万种路径组合20个城市就达到约2.4×10^18种常规遍历根本算不完。我在实际项目中接触TSP是因为一个工厂巡检路线优化的需求后来发现这个模型还能迁移到物流配送、无人机航线规划、电路板钻孔顺序优化等场景。很多初学者以为TSP只是个算法练习题其实它是组合优化问题的“通用测试平台”几乎所有智能优化算法都会拿TSP当基准测试函数。理解TSP的解法就等于掌握了组合优化问题的一把通用钥匙。TSP之所以难难点在于它没有太好的数学特征可以利用。目标函数本身就是非线性的、离散的路径之间还要满足“每个城市只能访问一次”的约束条件。传统的梯度下降、线性规划这些方法在这类问题面前基本失效。数学家研究了这么多年也没找到多项式时间的精确算法所以才需要启发式算法和元启发式算法来逼近最优解。粒子群算法、遗传算法、模拟退火、蚁群算法都是这个思路——不追求绝对最优而是用可控的时间成本去接近最优解在工程上已经足够用。1.2 粒子群算法凭什么能解TSP粒子群算法Particle Swarm OptimizationPSO最早是Kennedy和Eberhart在1995年受鸟群觅食行为启发提出的。基本原理不复杂一群粒子在搜索空间里飞行每个粒子都有位置和速度通过追踪“自己找到过的最好位置”和“整个群体找到过的最好位置”来调整飞行方向逐步收敛到全局最优区域。这里有个核心矛盾需要先解决标准粒子群算法是为连续优化问题设计的粒子的位置和速度都是连续实数而TSP的解是离散的城市排列顺序。怎么把两者桥接起来是整个项目最关键的思考点。我采用的方案是“随机键表示法”Random Key Representation。每个粒子的位置是一个长度为N的连续实数向量把这个向量从小到大排序排序后得到的城市索引顺序就是一条完整路径。举个例子有5个城市某个粒子的位置向量是[0.8, 0.3, 0.6, 0.1, 0.9]从小到大对应的城市索引是4、2、3、1、5那么路径就是4→2→3→1→5→4最后回到起点。这个映射方式非常巧妙它把连续空间的任意位置向量都能合法地转换成一个不重复城市的合法路径天然满足TSP“每个城市只走一次”的约束。选择PSO而不是遗传算法是因为PSO结构简单、参数少、收敛速度快。遗传算法需要处理选择、交叉、变异三个算子交叉和变异概率都要靠经验调配调参周期长。而PSO核心就是速度更新公式里两个学习因子和一个惯性权重理解起来直观写代码也很短。实际跑下来在中小规模TSP50个城市以内上PSO的收敛速度通常比遗传算法快代码量也少一半以上。不过也要说实话PSO在TSP上有个明显的短板容易早熟收敛。粒子群一旦全体朝某个局部最优聚集很难跳出来。所以在这篇实战中我会额外加入“变异重启”机制来缓解这个问题实测效果提升明显这部分后面专门讲。2. 核心原理拆解从连续优化到离散路径2.1 TSP的数学模型与距离计算TSP的数学表达很简洁给定N个城市的坐标求解一个访问顺序π (π1, π2, ..., πN)使得总路程最小。目标函数是min f(π) Σ d(πi, π(i1)) d(πN, π1)其中d(i,j)表示城市i到城市j的距离。在“静态欧式TSP”场景下距离直接用欧氏距离公式计算d(i,j) sqrt((xi - xj)^2 (yi - yj)^2)写代码时提前算好距离矩阵存成N×N的矩阵后续每次计算适应度就不用重复开方运算了能省不少时间。这是所有TSP代码实现里最基础的优化点。城市坐标我一般用两种方式生成一是随机生成方便测试不同规模二是用网上公开的TSP标准测试集比如berlin52柏林52个城市、eil51、att48这些。标准测试集的好处是有已知最优解可以验证算法实现的正确性。我下面代码里用随机生成的30个城市方便大家直接跑通想换标准数据集也很容易把坐标替换成数据集里的坐标就行。2.2 随机键表示法让粒子连续运动起来随机键表示法是我在这个项目里最想强调的一个设计它是整个算法能否跑通的关键。大家都知道标准PSO的速度和位置更新公式是v(t1) w * v(t) c1 * r1 * (pbest - x(t)) c2 * r2 * (gbest - x(t))x(t1) x(t) v(t1)公式里的x、v都是连续实数而TSP的解却是一个排列。如果直接把城市编号作为位置两个向量做减法得到的速度根本没有物理意义位置更新后还会出现“重复城市”的非法路径。随机键表示法就是为了解决这个不匹配问题。具体做法是把粒子的位置定义为一个N维连续向量每个维度对应一个城市但这个向量的具体数值大小并不直接代表城市编号而是代表一种“优先级”或者“排序权重”。当我们把向量按数值从小到大排序数值最小的那个维度对应的城市排在路径最前面数值第二小的维度对应的城市排在第二位以此类推。这样任何连续实数向量都能唯一对应一条合法路径。同样的道理粒子之间的“速度”就可以定义为这个连续向量的变化量。粒子位置更新时直接对随机键向量做加减法运算完全不用考虑会不会产生非法解。路径的合法性在排序这一步就已经保证了。这个思路的巧妙之处在于它把离散问题转化成了连续问题可以无缝套用PSO的所有公式不用对算法做额外改造。我第一次看论文时觉得这就是个小技巧实际做下来才意识到这其实是整个方案的地基。没有地基后面的所有迭代逻辑都站不住脚。2.3 适应度函数与速度更新公式的改造点适应度函数就是用来评价一个解好坏的函数。TSP的适应度就是路径总距离但要注意PSO算法习惯上适应度越大越好而TSP追求距离越小越好所以一种做法是直接取路程的倒数作为适应度另一种做法是直接比较路程大小来更新pbest和gbest。我个人偏向后者——不搞花样直接以距离值做比较代码更直观也减少一层计算开销。速度更新公式在这个项目里有几个改造点值得说。第一惯性权重w我用了线性递减策略从0.9递减到0.4。前期权重较大粒子大步探索避免过早聚拢后期权重变小粒子精细搜索提高收敛精度。这个策略几乎不增加计算量却能明显改善结果是性价比极高的一步优化。第二学习因子c1和c2我设置为1.5和1.5。c1代表粒子向自身历史最优学习的程度c2代表向群体最优学习的程度。r1、r2是0到1之间的均匀随机数用于增加搜索的随机性。这两个参数比较常见变化小一般不需要特别调。第三速度上限和位置范围要设置边界约束。速度太大容易震荡发散太小又会让搜索停滞我一般把速度限制在[-1, 1]范围内。由于随机键向量本身不需要有严格边界但为了保持排序的区分度我会把位置范围限制在[-10, 10]或者[-100, 100]之间边界问题在后续的常见问题里会详细展开。3. 完整代码实现与逐段解析3.1 主函数框架参数初始化与城市坐标生成我直接把完整的Matlab代码贴出来然后逐段解析。这个代码结构清晰注释完整在Matlab R2016b及以上版本都能直接运行。%% 粒子群算法解决TSP问题随机键表示法 clc; clear; close all; %% 1. 参数设置 numCity 30; % 城市数量 popSize 100; % 种群规模 maxIter 500; % 最大迭代次数 wStart 0.9; % 惯性权重初始值 wEnd 0.4; % 惯性权重最终值 c1 1.5; % 个体学习因子 c2 1.5; % 群体学习因子 posMin -10; % 粒子位置下界 posMax 10; % 粒子位置上界 velMin -1; % 速度下界 velMax 1; % 速度上界 %% 2. 生成城市坐标并计算距离矩阵 rng(42); % 固定随机种子保证实验结果可复现 cityPos rand(numCity, 2) * 100; % 生成[0,100]区间内的随机坐标 % 计算距离矩阵 distMatrix zeros(numCity, numCity); for i 1:numCity for j 1:numCity if i ~ j distMatrix(i,j) sqrt((cityPos(i,1)-cityPos(j,1))^2 ... (cityPos(i,2)-cityPos(j,2))^2); else distMatrix(i,j) inf; % 对角线设为inf避免自己到自己 end end end这里有几个细节要注意。rng(42)固定随机种子非常关键否则每次运行城市坐标不同、随机初始值不同结果无法横向比较。在城市坐标数据量较小的情况下用双重循环算距离矩阵没问题城市数量上千时就要改为向量化计算不然会卡到怀疑人生。有一回我测试500个城市的算例双重循环直接跑了几分钟才算出距离矩阵后来改成向量化瞬间完成这个优化收益非常大。3.2 核心循环适应度计算、个体最优与全局最优更新%% 3. 初始化粒子群 particlePos rand(popSize, numCity) * (posMax - posMin) posMin; particleVel zeros(popSize, numCity); pbest particlePos; % 个体历史最优位置 pbestDist inf(1, popSize); % 个体历史最优距离 gbest particlePos(1, :); % 全局最优位置 gbestDist inf; % 全局最优距离 %% 4. 主迭代循环 bestHistory zeros(maxIter, 1); % 记录每轮全局最优距离 avgHistory zeros(maxIter, 1); % 记录每轮平均距离 for iter 1:maxIter % 惯性权重线性递减 w wStart - (wStart - wEnd) * (iter / maxIter); % 计算每个粒子的适应度路径总距离 for i 1:popSize [route, dist] decodeRoute(particlePos(i,:), cityPos, distMatrix); % 更新个体历史最优 if dist pbestDist(i) pbestDist(i) dist; pbest(i, :) particlePos(i, :); end % 更新全局最优 if dist gbestDist gbestDist dist; gbest particlePos(i, :); end end % 速度与位置更新 for i 1:popSize r1 rand(1, numCity); r2 rand(1, numCity); particleVel(i, :) w * particleVel(i, :) ... c1 * r1 .* (pbest(i,:) - particlePos(i,:)) ... c2 * r2 .* (gbest - particlePos(i,:)); % 速度边界约束 particleVel(i, :) max(particleVel(i, :), velMin); particleVel(i, :) min(particleVel(i, :), velMax); % 位置更新 particlePos(i, :) particlePos(i, :) particleVel(i, :); % 位置边界约束 particlePos(i, :) max(particlePos(i, :), posMin); particlePos(i, :) min(particlePos(i, :), posMax); end bestHistory(iter) gbestDist; avgHistory(iter) mean(pbestDist); % 每50代打印一次进度 if mod(iter, 50) 0 fprintf(迭代次数: %d, 当前最优距离: %.2f\n, iter, gbestDist); end end解码函数decodeRoute是随机键表示法的核心需要单独定义代码在后面给出。每次循环里所有粒子都要解码一次计算路径距离这是整个迭代过程中计算量最大的部分。如果想优化性能可以考虑向量化距离计算但对中小规模问题影响不大先把逻辑跑通更重要。速度边界那里我用了max和min两次截断相当于把速度限制在[-1,1]区间内。这个细节很重要不加的话粒子速度可能越来越大位置朝边界猛冲后期所有粒子的随机键都跑到上下界附近排序区分度急剧下降算法就失效了。3.3 解码函数与结果可视化%% 5. 解码函数随机键向量 - 城市序列 function [route, dist] decodeRoute(pos, cityPos, distMatrix) [~, idx] sort(pos); % 排序得到城市索引 route idx; % 城市访问顺序 N length(route); dist 0; for k 1:N-1 dist dist distMatrix(route(k), route(k1)); end dist dist distMatrix(route(N), route(1)); % 回到起点 end %% 6. 输出最终结果并绘图 [bestRoute, bestRouteDist] decodeRoute(gbest, cityPos, distMatrix); fprintf(最优路径距离: %.2f\n, bestRouteDist); disp(最优访问顺序:); disp(bestRoute); figure; subplot(1, 2, 1); plot(cityPos(bestRoute, 1), cityPos(bestRoute, 2), k-o, LineWidth, 1.5, MarkerSize, 4); hold on; plot(cityPos(bestRoute(1), 1), cityPos(bestRoute(1), 2), rs, MarkerSize, 8); hold on; plot(cityPos(bestRoute(end), 1), cityPos(bestRoute(end), 2), g^, MarkerSize, 8); grid on; title(最优路径规划); xlabel(X坐标); ylabel(Y坐标); subplot(1, 2, 2); plot(1:maxIter, bestHistory, b-, LineWidth, 1.5); hold on; plot(1:maxIter, avgHistory, r--, LineWidth, 1); legend(全局最优, 种群平均, Location, northeast); xlabel(迭代次数); ylabel(路径距离); title(收敛曲线); grid on;这段代码里sort函数返回两个值第二个输出idx就是排序后的索引也就是城市路径顺序。这里很多人容易忽略一点Matlab的sort默认是升序排列所以得到的路径顺序是“数值小的随机键对应城市排在前面”这并不影响最终结果因为路径的距离只取决于城市之间的相对顺序而不取决于排序方向。绘图部分我习惯分两个子图左边画最优路径用方形标记起点、三角形标记终点方便直观看到路线的起始位置右边画收敛曲线同时显示全局最优和种群平均判断算法是否收敛。4. 参数调优与实验对比4.1 惯性权重与学习因子的影响参数调优是整个PSO实战里最需要耐心的环节。我做了几组对比实验用30个城市固定坐标集每组参数跑10次取平均值结果如下表所示参数组合最优距离平均值收敛代数稳定性评估w常数0.8, c1c21.5368.7120代中等偶尔陷入局部最优w线性递减0.9→0.4, c1c21.5342.1100代好多次运行接近同一结果w线性递减0.9→0.4, c1c22.0351.490代收敛快但精度稍差w线性递减0.9→0.4, c12.0, c21.5355.8110代前期探索强后期收敛慢w线性递减0.9→0.4, c11.5, c22.0338.595代最优群体引导更充分这组小实验佐证了两个经验判断一是惯性权重用线性递减确实比固定值效果好二是c2略大于c1时算法更倾向于群体协作在TSP这种解空间复杂的场景下效果更好。当然这只是30个城市的实验结果城市规模变化后参数最优区间可能会偏移需要按实际情况调整。4.2 种群规模与迭代次数的搭配种群规模和迭代次数的组合本质上是在“每代的搜索广度”和“搜索深度”之间做取舍。种群太小比如20个粒子搜索覆盖面不足很容易陷入局部最优跑500代也跳不出来。种群太大比如500个粒子虽然单次搜索能力强但每代的计算量翻倍迭代次数就得相应减少总时间成本未必划算。我的经验区间是城市数量不超过50时种群规模100、迭代次数500的组合性价比最高。城市数量100以上时建议种群规模200、迭代次数1000起步。有一个判断算法是否“吃饱”的技巧看收敛曲线尾部是否已经平坦。如果连续100代全局最优距离都没有变化说明算法已经收敛再多跑也没有意义。如果时间受限还有一个加速小技巧前300代用较大的惯性权重快速逼近后面200代用小权重精调。这就是模拟退火的思路借用到PSO里很多论文里的“改进PSO”本质就是在做这类操作工程上直接做线性递减就够用了。4.3 不同城市规模下的收敛表现我测试了三种城市规模30个、50个和80个。30个城市时PSO表现得非常轻松前50代快速收敛最后能得到接近最优的路径。50个城市时算法性能开始分化必须配合变异重启机制才能稳定逼近优解。80个城市时单纯用基础PSO已经很难看500代跑完得到的路径明显有交叉说明算法早熟了。城市数量种群规模迭代次数基础PSO最优距离加入变异重启后30100500342.1335.850150800六组实验中最优稳定在520左右488.3802001000收敛到780附近后停滞702.6这个对比说明一个扎心的事实基础PSO在TSP上的表现随着城市数量增加会快速恶化。所以如果打算用PSO处理大规模TSP一定要加机制增强跳出局部最优的能力。变异重启就是其中一种后面常见问题章节我会详细介绍具体实现。4.4 变异重启机制一个低成本高收益的改进针对PSO容易早熟的问题我加了一个非常简单的变异操作每次迭代结束时随机选择一部分粒子把它们的随机键向量打乱重来。实现方式是在主循环末尾加几行判断if rand 0.05 index randi(popSize); particlePos(index, :) rand(1, numCity) * (posMax - posMin) posMin; particleVel(index, :) zeros(1, numCity); pbestDist(index) inf; end这段代码的作用是每代有5%的概率随机挑一个粒子重置。被重置的粒子的历史最优信息也清空相当于让这个粒子“失忆”后重新出发。这个机制看似简单但能有效防止全部粒子都朝同一个局部最优坍缩。实测效果是50个城市场景下最优距离从520左右改善到488提升幅度近6%。代价仅仅是每次迭代多一次判断性能损耗几乎可以忽略。在使用变异重启时建议把重置概率控制在2%到10%之间。太低起不到作用太高则会导致算法一直在“随机搜索”而损失收敛精度。另外重置的对象也不宜太多一次只重置一两个粒子最合适否则会打乱粒子群已经形成的良好聚拢结构。5. 常见问题与踩坑实录5.1 粒子位置更新后城市重复了怎么办这是很多初学者在把PSO套到TSP上时容易卡住的第一道坎。直接拿城市编号当粒子位置来更新必然会出现一个城市被访问两次、另一个城市从未访问的问题。我给的解决方案前面已经说明了就是随机键表示法。这里想再强调一遍把“解”和“编码”分离这是处理离散优化问题的核心思路。粒子位置向量的每个分量不直接对应城市编号而是对应一个排序权重。只要最终解码时用sort排序路径必然合法重复问题从根源上消失。如果你用的是交换序类的方法粒子位置就是城市序列速度定义为交换操作那处理重复路径的方式就完全不同每次更新后要逐一检查非法城市编号并做修复。那个方案在实现上要复杂得多测试下来效果也不比随机键稳定。从工程效率角度我推荐随机键方案。5.2 算法早熟收敛怎么破早熟收敛的表现是迭代没几步收敛曲线就完全走平画出来的路径图有明显交叉。这本质上是粒子群多样性丧失所有粒子都跑到同一个局部最优附近了。我的处理优先级如下第一检查惯性权重是否线性递减固定大权重会导致后期无法精细搜索固定小权重则前期就聚拢。第二加入前面提到的变异重启机制。第三如果条件允许试试把粒子群分成几个子群各自独立搜索每隔若干代交换信息——这种“小生境”方案能进一步提升多样性不过代码复杂度会高一些。有时我也遇到“假收敛”的情况看起来收敛曲线平了但实际是收敛到了比较差的解然后算法卡住不动。这时候我会检查一下种群平均距离和全局最优距离的关系。如果两者非常接近说明群体多样性已经非常差需要强制重置部分粒子如果平均距离还比最优距离大很多说明粒子还在探索中多跑几代说不定还有改善。5.3 结果不稳定、每次运行都不一样PSO是随机算法结果有波动非常正常。但如果你发现每次运行的结果差异特别大可能有两种原因一是没有固定随机种子城市坐标和初始粒子群每次都不一样这必然导致结果漂移。排查方法是在代码开头加rng(固定数字)保证每次运行随机序列一致。二是算法本身稳定性差收敛不到一个可重复的较优解区间这时候需要按5.2节的方法增强全局搜索能力。如果只是想验证算法在某个算例上的表现建议跑多次取平均值和标准差。很少有一次运行就能得出可靠结论的跑10次取均值是比较常规的评估方式。我在项目报告里一般会给出三组数字最优值、平均值、标准差三者共同描述算法的性能表现。5.4 距离计算与边界约束的几个隐蔽坑距离计算有个坑如果忘记将对角线距离设为inf某个城市的随机键恰好排在相邻位置时会出现路径中连续两个相同城市总距离中会包含一个零距离项表面上路径看起来合法但路由本质上是无效的。虽然后面排序解码时连续相同索引出现的概率很低但计算总距离时最好还是过滤一下。边界约束那个问题我已经在前面提过。如果位置上限设得太小比如只有1那么所有随机键都挤在[0,1]区间里排序时数值区分度可能不足最终路径只是不同随机初始值的排序结果算法的搜索能力就废了。把位置范围放宽到[-10,10]或者更大可以给粒子足够的变化空间。Matlab绘图还有一个常见烦恼城市坐标点多时线一多就糊成一团。可以通过设置LineWidth和MarkerSize让路径更醒目也可以把路径图上的城市编号标出来用text函数逐一添加方便检查路径是否存在交叉。5.5 Matlab代码执行效率优化很多人跑PSO时发现代码很慢第一个想到的是把循环改小其实真正的瓶颈往往在解码环节。仔细分析会发现每次迭代要做popSize次解码每次都从sort开始sort的复杂度是O(N log N)总复杂度就是O(maxIter × popSize × N log N)。当N上千时这个开销就非常可观了。优化手段有两个方向一是把解码函数里的距离累加循环改为向量化计算例如先取出路径对应的连续城市索引再使用sum和diag配合提取距离矩阵中的对应元素二是减少不必要的解码次数比如只有当距离可能改善时才解码。不过代码可读性也很重要中小规模直接跑就好不要为了优化牺牲理解成本。另外一个实际经验是Matlab里rand和rng的调用开销不小如果在循环体内反复生成新的RandStream对象会非常拖慢速度。建议每一代只调用一次rng或者在初始化时固定后续循环直接用rand自然生成随机数即可。6. 从代码到项目的经验总结写完这套代码之后我自己最大的一个感受是粒子群算法的“难”不在于算法本身而在于如何把一个实际问题转化成算法能处理的形式。随机键表示法就是一个很好的例子它只有三行代码却解决了连续优化和离散问题之间的鸿沟。很多技术难点想通了原理就简单想不通就在那里绕圈。另一个体会是盲目追求“改进算法”不如先把基础版本调好跑透。网上各种期刊论文里的PSO变体让人眼花缭乱但很多时候基础版PSO加一个简单的变异重启就已经能满足工程需求。初学阶段务必先把标准算法吃透再有针对性地做改进。我见过很多初学者一上来就加十几个改进点最后代码复杂到根本定位不了问题反而事倍功半。如果你打算把这套代码应用到自己的项目我建议的扩展路径是先解决VRP车辆路径问题在TSP的路径编码基础上增加车辆容量约束和时间窗约束再考虑多目标优化比如同时优化总路程和行驶时间最后还可以尝试把PSO和局部搜索算法结合用PSO做全局探索用2-opt做局部精调这种混合策略在100个城市以上的场景下效果会显著提升。最后分享一个小技巧调参时不要每次都从头跑完整代码可以把收敛历史保存下来画对比图。我通常会把每一次跑完的最优距离、平均距离和运行时间记录在一个Excel表里做参数决策时直接看表格数据比凭感觉比大小靠谱得多。这个方法虽然不起眼却是我所有优化项目里最常用的决策工具。拿着这些数据说话不管是对自己还是对项目干系人都更有说服力。 粒子群算法求解TSP旅行商问题：随机键表示法与Matlab实战
1. 项目概述为什么用粒子群算法求解TSP1.1 旅行商问题到底难在哪旅行商问题Traveling Salesman ProblemTSP是组合优化领域最经典的NP难问题之一。简单说就是一个商人要跑遍N个城市每个城市只能去一次最后回到出发城市怎么走总路程最短。这个问题听起来简单但复杂度随着城市数量增加呈现阶乘式爆炸。10个城市有约36万种路径组合20个城市就达到约2.4×10^18种常规遍历根本算不完。我在实际项目中接触TSP是因为一个工厂巡检路线优化的需求后来发现这个模型还能迁移到物流配送、无人机航线规划、电路板钻孔顺序优化等场景。很多初学者以为TSP只是个算法练习题其实它是组合优化问题的“通用测试平台”几乎所有智能优化算法都会拿TSP当基准测试函数。理解TSP的解法就等于掌握了组合优化问题的一把通用钥匙。TSP之所以难难点在于它没有太好的数学特征可以利用。目标函数本身就是非线性的、离散的路径之间还要满足“每个城市只能访问一次”的约束条件。传统的梯度下降、线性规划这些方法在这类问题面前基本失效。数学家研究了这么多年也没找到多项式时间的精确算法所以才需要启发式算法和元启发式算法来逼近最优解。粒子群算法、遗传算法、模拟退火、蚁群算法都是这个思路——不追求绝对最优而是用可控的时间成本去接近最优解在工程上已经足够用。1.2 粒子群算法凭什么能解TSP粒子群算法Particle Swarm OptimizationPSO最早是Kennedy和Eberhart在1995年受鸟群觅食行为启发提出的。基本原理不复杂一群粒子在搜索空间里飞行每个粒子都有位置和速度通过追踪“自己找到过的最好位置”和“整个群体找到过的最好位置”来调整飞行方向逐步收敛到全局最优区域。这里有个核心矛盾需要先解决标准粒子群算法是为连续优化问题设计的粒子的位置和速度都是连续实数而TSP的解是离散的城市排列顺序。怎么把两者桥接起来是整个项目最关键的思考点。我采用的方案是“随机键表示法”Random Key Representation。每个粒子的位置是一个长度为N的连续实数向量把这个向量从小到大排序排序后得到的城市索引顺序就是一条完整路径。举个例子有5个城市某个粒子的位置向量是[0.8, 0.3, 0.6, 0.1, 0.9]从小到大对应的城市索引是4、2、3、1、5那么路径就是4→2→3→1→5→4最后回到起点。这个映射方式非常巧妙它把连续空间的任意位置向量都能合法地转换成一个不重复城市的合法路径天然满足TSP“每个城市只走一次”的约束。选择PSO而不是遗传算法是因为PSO结构简单、参数少、收敛速度快。遗传算法需要处理选择、交叉、变异三个算子交叉和变异概率都要靠经验调配调参周期长。而PSO核心就是速度更新公式里两个学习因子和一个惯性权重理解起来直观写代码也很短。实际跑下来在中小规模TSP50个城市以内上PSO的收敛速度通常比遗传算法快代码量也少一半以上。不过也要说实话PSO在TSP上有个明显的短板容易早熟收敛。粒子群一旦全体朝某个局部最优聚集很难跳出来。所以在这篇实战中我会额外加入“变异重启”机制来缓解这个问题实测效果提升明显这部分后面专门讲。2. 核心原理拆解从连续优化到离散路径2.1 TSP的数学模型与距离计算TSP的数学表达很简洁给定N个城市的坐标求解一个访问顺序π (π1, π2, ..., πN)使得总路程最小。目标函数是min f(π) Σ d(πi, π(i1)) d(πN, π1)其中d(i,j)表示城市i到城市j的距离。在“静态欧式TSP”场景下距离直接用欧氏距离公式计算d(i,j) sqrt((xi - xj)^2 (yi - yj)^2)写代码时提前算好距离矩阵存成N×N的矩阵后续每次计算适应度就不用重复开方运算了能省不少时间。这是所有TSP代码实现里最基础的优化点。城市坐标我一般用两种方式生成一是随机生成方便测试不同规模二是用网上公开的TSP标准测试集比如berlin52柏林52个城市、eil51、att48这些。标准测试集的好处是有已知最优解可以验证算法实现的正确性。我下面代码里用随机生成的30个城市方便大家直接跑通想换标准数据集也很容易把坐标替换成数据集里的坐标就行。2.2 随机键表示法让粒子连续运动起来随机键表示法是我在这个项目里最想强调的一个设计它是整个算法能否跑通的关键。大家都知道标准PSO的速度和位置更新公式是v(t1) w * v(t) c1 * r1 * (pbest - x(t)) c2 * r2 * (gbest - x(t))x(t1) x(t) v(t1)公式里的x、v都是连续实数而TSP的解却是一个排列。如果直接把城市编号作为位置两个向量做减法得到的速度根本没有物理意义位置更新后还会出现“重复城市”的非法路径。随机键表示法就是为了解决这个不匹配问题。具体做法是把粒子的位置定义为一个N维连续向量每个维度对应一个城市但这个向量的具体数值大小并不直接代表城市编号而是代表一种“优先级”或者“排序权重”。当我们把向量按数值从小到大排序数值最小的那个维度对应的城市排在路径最前面数值第二小的维度对应的城市排在第二位以此类推。这样任何连续实数向量都能唯一对应一条合法路径。同样的道理粒子之间的“速度”就可以定义为这个连续向量的变化量。粒子位置更新时直接对随机键向量做加减法运算完全不用考虑会不会产生非法解。路径的合法性在排序这一步就已经保证了。这个思路的巧妙之处在于它把离散问题转化成了连续问题可以无缝套用PSO的所有公式不用对算法做额外改造。我第一次看论文时觉得这就是个小技巧实际做下来才意识到这其实是整个方案的地基。没有地基后面的所有迭代逻辑都站不住脚。2.3 适应度函数与速度更新公式的改造点适应度函数就是用来评价一个解好坏的函数。TSP的适应度就是路径总距离但要注意PSO算法习惯上适应度越大越好而TSP追求距离越小越好所以一种做法是直接取路程的倒数作为适应度另一种做法是直接比较路程大小来更新pbest和gbest。我个人偏向后者——不搞花样直接以距离值做比较代码更直观也减少一层计算开销。速度更新公式在这个项目里有几个改造点值得说。第一惯性权重w我用了线性递减策略从0.9递减到0.4。前期权重较大粒子大步探索避免过早聚拢后期权重变小粒子精细搜索提高收敛精度。这个策略几乎不增加计算量却能明显改善结果是性价比极高的一步优化。第二学习因子c1和c2我设置为1.5和1.5。c1代表粒子向自身历史最优学习的程度c2代表向群体最优学习的程度。r1、r2是0到1之间的均匀随机数用于增加搜索的随机性。这两个参数比较常见变化小一般不需要特别调。第三速度上限和位置范围要设置边界约束。速度太大容易震荡发散太小又会让搜索停滞我一般把速度限制在[-1, 1]范围内。由于随机键向量本身不需要有严格边界但为了保持排序的区分度我会把位置范围限制在[-10, 10]或者[-100, 100]之间边界问题在后续的常见问题里会详细展开。3. 完整代码实现与逐段解析3.1 主函数框架参数初始化与城市坐标生成我直接把完整的Matlab代码贴出来然后逐段解析。这个代码结构清晰注释完整在Matlab R2016b及以上版本都能直接运行。%% 粒子群算法解决TSP问题随机键表示法 clc; clear; close all; %% 1. 参数设置 numCity 30; % 城市数量 popSize 100; % 种群规模 maxIter 500; % 最大迭代次数 wStart 0.9; % 惯性权重初始值 wEnd 0.4; % 惯性权重最终值 c1 1.5; % 个体学习因子 c2 1.5; % 群体学习因子 posMin -10; % 粒子位置下界 posMax 10; % 粒子位置上界 velMin -1; % 速度下界 velMax 1; % 速度上界 %% 2. 生成城市坐标并计算距离矩阵 rng(42); % 固定随机种子保证实验结果可复现 cityPos rand(numCity, 2) * 100; % 生成[0,100]区间内的随机坐标 % 计算距离矩阵 distMatrix zeros(numCity, numCity); for i 1:numCity for j 1:numCity if i ~ j distMatrix(i,j) sqrt((cityPos(i,1)-cityPos(j,1))^2 ... (cityPos(i,2)-cityPos(j,2))^2); else distMatrix(i,j) inf; % 对角线设为inf避免自己到自己 end end end这里有几个细节要注意。rng(42)固定随机种子非常关键否则每次运行城市坐标不同、随机初始值不同结果无法横向比较。在城市坐标数据量较小的情况下用双重循环算距离矩阵没问题城市数量上千时就要改为向量化计算不然会卡到怀疑人生。有一回我测试500个城市的算例双重循环直接跑了几分钟才算出距离矩阵后来改成向量化瞬间完成这个优化收益非常大。3.2 核心循环适应度计算、个体最优与全局最优更新%% 3. 初始化粒子群 particlePos rand(popSize, numCity) * (posMax - posMin) posMin; particleVel zeros(popSize, numCity); pbest particlePos; % 个体历史最优位置 pbestDist inf(1, popSize); % 个体历史最优距离 gbest particlePos(1, :); % 全局最优位置 gbestDist inf; % 全局最优距离 %% 4. 主迭代循环 bestHistory zeros(maxIter, 1); % 记录每轮全局最优距离 avgHistory zeros(maxIter, 1); % 记录每轮平均距离 for iter 1:maxIter % 惯性权重线性递减 w wStart - (wStart - wEnd) * (iter / maxIter); % 计算每个粒子的适应度路径总距离 for i 1:popSize [route, dist] decodeRoute(particlePos(i,:), cityPos, distMatrix); % 更新个体历史最优 if dist pbestDist(i) pbestDist(i) dist; pbest(i, :) particlePos(i, :); end % 更新全局最优 if dist gbestDist gbestDist dist; gbest particlePos(i, :); end end % 速度与位置更新 for i 1:popSize r1 rand(1, numCity); r2 rand(1, numCity); particleVel(i, :) w * particleVel(i, :) ... c1 * r1 .* (pbest(i,:) - particlePos(i,:)) ... c2 * r2 .* (gbest - particlePos(i,:)); % 速度边界约束 particleVel(i, :) max(particleVel(i, :), velMin); particleVel(i, :) min(particleVel(i, :), velMax); % 位置更新 particlePos(i, :) particlePos(i, :) particleVel(i, :); % 位置边界约束 particlePos(i, :) max(particlePos(i, :), posMin); particlePos(i, :) min(particlePos(i, :), posMax); end bestHistory(iter) gbestDist; avgHistory(iter) mean(pbestDist); % 每50代打印一次进度 if mod(iter, 50) 0 fprintf(迭代次数: %d, 当前最优距离: %.2f\n, iter, gbestDist); end end解码函数decodeRoute是随机键表示法的核心需要单独定义代码在后面给出。每次循环里所有粒子都要解码一次计算路径距离这是整个迭代过程中计算量最大的部分。如果想优化性能可以考虑向量化距离计算但对中小规模问题影响不大先把逻辑跑通更重要。速度边界那里我用了max和min两次截断相当于把速度限制在[-1,1]区间内。这个细节很重要不加的话粒子速度可能越来越大位置朝边界猛冲后期所有粒子的随机键都跑到上下界附近排序区分度急剧下降算法就失效了。3.3 解码函数与结果可视化%% 5. 解码函数随机键向量 - 城市序列 function [route, dist] decodeRoute(pos, cityPos, distMatrix) [~, idx] sort(pos); % 排序得到城市索引 route idx; % 城市访问顺序 N length(route); dist 0; for k 1:N-1 dist dist distMatrix(route(k), route(k1)); end dist dist distMatrix(route(N), route(1)); % 回到起点 end %% 6. 输出最终结果并绘图 [bestRoute, bestRouteDist] decodeRoute(gbest, cityPos, distMatrix); fprintf(最优路径距离: %.2f\n, bestRouteDist); disp(最优访问顺序:); disp(bestRoute); figure; subplot(1, 2, 1); plot(cityPos(bestRoute, 1), cityPos(bestRoute, 2), k-o, LineWidth, 1.5, MarkerSize, 4); hold on; plot(cityPos(bestRoute(1), 1), cityPos(bestRoute(1), 2), rs, MarkerSize, 8); hold on; plot(cityPos(bestRoute(end), 1), cityPos(bestRoute(end), 2), g^, MarkerSize, 8); grid on; title(最优路径规划); xlabel(X坐标); ylabel(Y坐标); subplot(1, 2, 2); plot(1:maxIter, bestHistory, b-, LineWidth, 1.5); hold on; plot(1:maxIter, avgHistory, r--, LineWidth, 1); legend(全局最优, 种群平均, Location, northeast); xlabel(迭代次数); ylabel(路径距离); title(收敛曲线); grid on;这段代码里sort函数返回两个值第二个输出idx就是排序后的索引也就是城市路径顺序。这里很多人容易忽略一点Matlab的sort默认是升序排列所以得到的路径顺序是“数值小的随机键对应城市排在前面”这并不影响最终结果因为路径的距离只取决于城市之间的相对顺序而不取决于排序方向。绘图部分我习惯分两个子图左边画最优路径用方形标记起点、三角形标记终点方便直观看到路线的起始位置右边画收敛曲线同时显示全局最优和种群平均判断算法是否收敛。4. 参数调优与实验对比4.1 惯性权重与学习因子的影响参数调优是整个PSO实战里最需要耐心的环节。我做了几组对比实验用30个城市固定坐标集每组参数跑10次取平均值结果如下表所示参数组合最优距离平均值收敛代数稳定性评估w常数0.8, c1c21.5368.7120代中等偶尔陷入局部最优w线性递减0.9→0.4, c1c21.5342.1100代好多次运行接近同一结果w线性递减0.9→0.4, c1c22.0351.490代收敛快但精度稍差w线性递减0.9→0.4, c12.0, c21.5355.8110代前期探索强后期收敛慢w线性递减0.9→0.4, c11.5, c22.0338.595代最优群体引导更充分这组小实验佐证了两个经验判断一是惯性权重用线性递减确实比固定值效果好二是c2略大于c1时算法更倾向于群体协作在TSP这种解空间复杂的场景下效果更好。当然这只是30个城市的实验结果城市规模变化后参数最优区间可能会偏移需要按实际情况调整。4.2 种群规模与迭代次数的搭配种群规模和迭代次数的组合本质上是在“每代的搜索广度”和“搜索深度”之间做取舍。种群太小比如20个粒子搜索覆盖面不足很容易陷入局部最优跑500代也跳不出来。种群太大比如500个粒子虽然单次搜索能力强但每代的计算量翻倍迭代次数就得相应减少总时间成本未必划算。我的经验区间是城市数量不超过50时种群规模100、迭代次数500的组合性价比最高。城市数量100以上时建议种群规模200、迭代次数1000起步。有一个判断算法是否“吃饱”的技巧看收敛曲线尾部是否已经平坦。如果连续100代全局最优距离都没有变化说明算法已经收敛再多跑也没有意义。如果时间受限还有一个加速小技巧前300代用较大的惯性权重快速逼近后面200代用小权重精调。这就是模拟退火的思路借用到PSO里很多论文里的“改进PSO”本质就是在做这类操作工程上直接做线性递减就够用了。4.3 不同城市规模下的收敛表现我测试了三种城市规模30个、50个和80个。30个城市时PSO表现得非常轻松前50代快速收敛最后能得到接近最优的路径。50个城市时算法性能开始分化必须配合变异重启机制才能稳定逼近优解。80个城市时单纯用基础PSO已经很难看500代跑完得到的路径明显有交叉说明算法早熟了。城市数量种群规模迭代次数基础PSO最优距离加入变异重启后30100500342.1335.850150800六组实验中最优稳定在520左右488.3802001000收敛到780附近后停滞702.6这个对比说明一个扎心的事实基础PSO在TSP上的表现随着城市数量增加会快速恶化。所以如果打算用PSO处理大规模TSP一定要加机制增强跳出局部最优的能力。变异重启就是其中一种后面常见问题章节我会详细介绍具体实现。4.4 变异重启机制一个低成本高收益的改进针对PSO容易早熟的问题我加了一个非常简单的变异操作每次迭代结束时随机选择一部分粒子把它们的随机键向量打乱重来。实现方式是在主循环末尾加几行判断if rand 0.05 index randi(popSize); particlePos(index, :) rand(1, numCity) * (posMax - posMin) posMin; particleVel(index, :) zeros(1, numCity); pbestDist(index) inf; end这段代码的作用是每代有5%的概率随机挑一个粒子重置。被重置的粒子的历史最优信息也清空相当于让这个粒子“失忆”后重新出发。这个机制看似简单但能有效防止全部粒子都朝同一个局部最优坍缩。实测效果是50个城市场景下最优距离从520左右改善到488提升幅度近6%。代价仅仅是每次迭代多一次判断性能损耗几乎可以忽略。在使用变异重启时建议把重置概率控制在2%到10%之间。太低起不到作用太高则会导致算法一直在“随机搜索”而损失收敛精度。另外重置的对象也不宜太多一次只重置一两个粒子最合适否则会打乱粒子群已经形成的良好聚拢结构。5. 常见问题与踩坑实录5.1 粒子位置更新后城市重复了怎么办这是很多初学者在把PSO套到TSP上时容易卡住的第一道坎。直接拿城市编号当粒子位置来更新必然会出现一个城市被访问两次、另一个城市从未访问的问题。我给的解决方案前面已经说明了就是随机键表示法。这里想再强调一遍把“解”和“编码”分离这是处理离散优化问题的核心思路。粒子位置向量的每个分量不直接对应城市编号而是对应一个排序权重。只要最终解码时用sort排序路径必然合法重复问题从根源上消失。如果你用的是交换序类的方法粒子位置就是城市序列速度定义为交换操作那处理重复路径的方式就完全不同每次更新后要逐一检查非法城市编号并做修复。那个方案在实现上要复杂得多测试下来效果也不比随机键稳定。从工程效率角度我推荐随机键方案。5.2 算法早熟收敛怎么破早熟收敛的表现是迭代没几步收敛曲线就完全走平画出来的路径图有明显交叉。这本质上是粒子群多样性丧失所有粒子都跑到同一个局部最优附近了。我的处理优先级如下第一检查惯性权重是否线性递减固定大权重会导致后期无法精细搜索固定小权重则前期就聚拢。第二加入前面提到的变异重启机制。第三如果条件允许试试把粒子群分成几个子群各自独立搜索每隔若干代交换信息——这种“小生境”方案能进一步提升多样性不过代码复杂度会高一些。有时我也遇到“假收敛”的情况看起来收敛曲线平了但实际是收敛到了比较差的解然后算法卡住不动。这时候我会检查一下种群平均距离和全局最优距离的关系。如果两者非常接近说明群体多样性已经非常差需要强制重置部分粒子如果平均距离还比最优距离大很多说明粒子还在探索中多跑几代说不定还有改善。5.3 结果不稳定、每次运行都不一样PSO是随机算法结果有波动非常正常。但如果你发现每次运行的结果差异特别大可能有两种原因一是没有固定随机种子城市坐标和初始粒子群每次都不一样这必然导致结果漂移。排查方法是在代码开头加rng(固定数字)保证每次运行随机序列一致。二是算法本身稳定性差收敛不到一个可重复的较优解区间这时候需要按5.2节的方法增强全局搜索能力。如果只是想验证算法在某个算例上的表现建议跑多次取平均值和标准差。很少有一次运行就能得出可靠结论的跑10次取均值是比较常规的评估方式。我在项目报告里一般会给出三组数字最优值、平均值、标准差三者共同描述算法的性能表现。5.4 距离计算与边界约束的几个隐蔽坑距离计算有个坑如果忘记将对角线距离设为inf某个城市的随机键恰好排在相邻位置时会出现路径中连续两个相同城市总距离中会包含一个零距离项表面上路径看起来合法但路由本质上是无效的。虽然后面排序解码时连续相同索引出现的概率很低但计算总距离时最好还是过滤一下。边界约束那个问题我已经在前面提过。如果位置上限设得太小比如只有1那么所有随机键都挤在[0,1]区间里排序时数值区分度可能不足最终路径只是不同随机初始值的排序结果算法的搜索能力就废了。把位置范围放宽到[-10,10]或者更大可以给粒子足够的变化空间。Matlab绘图还有一个常见烦恼城市坐标点多时线一多就糊成一团。可以通过设置LineWidth和MarkerSize让路径更醒目也可以把路径图上的城市编号标出来用text函数逐一添加方便检查路径是否存在交叉。5.5 Matlab代码执行效率优化很多人跑PSO时发现代码很慢第一个想到的是把循环改小其实真正的瓶颈往往在解码环节。仔细分析会发现每次迭代要做popSize次解码每次都从sort开始sort的复杂度是O(N log N)总复杂度就是O(maxIter × popSize × N log N)。当N上千时这个开销就非常可观了。优化手段有两个方向一是把解码函数里的距离累加循环改为向量化计算例如先取出路径对应的连续城市索引再使用sum和diag配合提取距离矩阵中的对应元素二是减少不必要的解码次数比如只有当距离可能改善时才解码。不过代码可读性也很重要中小规模直接跑就好不要为了优化牺牲理解成本。另外一个实际经验是Matlab里rand和rng的调用开销不小如果在循环体内反复生成新的RandStream对象会非常拖慢速度。建议每一代只调用一次rng或者在初始化时固定后续循环直接用rand自然生成随机数即可。6. 从代码到项目的经验总结写完这套代码之后我自己最大的一个感受是粒子群算法的“难”不在于算法本身而在于如何把一个实际问题转化成算法能处理的形式。随机键表示法就是一个很好的例子它只有三行代码却解决了连续优化和离散问题之间的鸿沟。很多技术难点想通了原理就简单想不通就在那里绕圈。另一个体会是盲目追求“改进算法”不如先把基础版本调好跑透。网上各种期刊论文里的PSO变体让人眼花缭乱但很多时候基础版PSO加一个简单的变异重启就已经能满足工程需求。初学阶段务必先把标准算法吃透再有针对性地做改进。我见过很多初学者一上来就加十几个改进点最后代码复杂到根本定位不了问题反而事倍功半。如果你打算把这套代码应用到自己的项目我建议的扩展路径是先解决VRP车辆路径问题在TSP的路径编码基础上增加车辆容量约束和时间窗约束再考虑多目标优化比如同时优化总路程和行驶时间最后还可以尝试把PSO和局部搜索算法结合用PSO做全局探索用2-opt做局部精调这种混合策略在100个城市以上的场景下效果会显著提升。最后分享一个小技巧调参时不要每次都从头跑完整代码可以把收敛历史保存下来画对比图。我通常会把每一次跑完的最优距离、平均距离和运行时间记录在一个Excel表里做参数决策时直接看表格数据比凭感觉比大小靠谱得多。这个方法虽然不起眼却是我所有优化项目里最常用的决策工具。拿着这些数据说话不管是对自己还是对项目干系人都更有说服力。