1. 项目概述:为什么用粒子群算法求解TSP
1.1 旅行商问题到底难在哪
旅行商问题(Traveling Salesman Problem,TSP)是组合优化领域最经典的NP难问题之一。简单说就是:一个商人要跑遍N个城市,每个城市只能去一次,最后回到出发城市,怎么走总路程最短。这个问题听起来简单,但复杂度随着城市数量增加呈现阶乘式爆炸。10个城市有约36万种路径组合,20个城市就达到约2.4×10^18种,常规遍历根本算不完。
我在实际项目中接触TSP是因为一个工厂巡检路线优化的需求,后来发现这个模型还能迁移到物流配送、无人机航线规划、电路板钻孔顺序优化等场景。很多初学者以为TSP只是个算法练习题,其实它是组合优化问题的“通用测试平台”,几乎所有智能优化算法都会拿TSP当基准测试函数。理解TSP的解法,就等于掌握了组合优化问题的一把通用钥匙。
TSP之所以难,难点在于它没有太好的数学特征可以利用。目标函数本身就是非线性的、离散的,路径之间还要满足“每个城市只能访问一次”的约束条件。传统的梯度下降、线性规划这些方法在这类问题面前基本失效。数学家研究了这么多年,也没找到多项式时间的精确算法,所以才需要启发式算法和元启发式算法来逼近最优解。粒子群算法、遗传算法、模拟退火、蚁群算法都是这个思路——不追求绝对最优,而是用可控的时间成本去接近最优解,在工程上已经足够用。
1.2 粒子群算法凭什么能解TSP
粒子群算法(Particle Swarm Optimization,PSO)最早是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核心就是速度更新公式里两个学习因子和一个惯性权重,理解起来直观,写代码也很短。实际跑下来,在中小规模TSP(50个城市以内)上,PSO的收敛速度通常比遗传算法快,代码量也少一半以上。
不过也要说实话,PSO在TSP上有个明显的短板:容易早熟收敛。粒子群一旦全体朝某个局部最优聚集,很难跳出来。所以在这篇实战中,我会额外加入“变异重启”机制来缓解这个问题,实测效果提升明显,这部分后面专门讲。
2. 核心原理拆解:从连续优化到离散路径
2.1 TSP的数学模型与距离计算
TSP的数学表达很简洁:给定N个城市的坐标,求解一个访问顺序π = (π1, π2, ..., πN),使得总路程最小。目标函数是:
min f(π) = Σ d(πi, π(i+1)) + 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(t+1) = w * v(t) + c1 * r1 * (pbest - x(t)) + c2 * r2 * (gbest - x(t))
x(t+1) = x(t) + v(t+1)
公式里的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(k+1)); 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, c1=c2=1.5 | 368.7 | 120代 | 中等,偶尔陷入局部最优 |
| w线性递减0.9→0.4, c1=c2=1.5 | 342.1 | 100代 | 好,多次运行接近同一结果 |
| w线性递减0.9→0.4, c1=c2=2.0 | 351.4 | 90代 | 收敛快但精度稍差 |
| w线性递减0.9→0.4, c1=2.0, c2=1.5 | 355.8 | 110代 | 前期探索强,后期收敛慢 |
| w线性递减0.9→0.4, c1=1.5, c2=2.0 | 338.5 | 95代 | 最优,群体引导更充分 |
这组小实验佐证了两个经验判断:一是惯性权重用线性递减确实比固定值效果好;二是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最优距离 | 加入变异重启后 |
|---|---|---|---|---|
| 30 | 100 | 500 | 342.1 | 335.8 |
| 50 | 150 | 800 | 六组实验中最优稳定在520左右 | 488.3 |
| 80 | 200 | 1000 | 收敛到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表里,做参数决策时直接看表格数据,比凭感觉比大小靠谱得多。这个方法虽然不起眼,却是我所有优化项目里最常用的决策工具。拿着这些数据说话,不管是对自己还是对项目干系人都更有说服力。