简介：Matlab混合粒子群算法（HPSO）求解TSP的完整代码实例，面向智能优化算法初学者、Matlab开发者、运筹优化课程设计等场景。算法在标准粒子群基础上引入遗传操作或局部搜索，以更有效地逼近旅行商问题的最短路径，也为与其他启发式方法对比和调参练习提供了直接可用的基础工程。压缩包共4个文件，包含3个m文件和1个txt文件：m文件分别负责主程序、路径长度适应度计算与距离矩阵生成，txt文件提供城市坐标测试数据，整体压缩包仅3KB，解压后即可直接运行。已有2167人浏览/学习。代码对关键步骤均有详细注释，覆盖参数初始化、适应度计算、个体与全局最优更新、混合策略以及终止判断等环节，能帮助读者逐段对照算法原理理解实现细节，也方便在此基础上调整惯性权重、学习因子或替换局部搜索方式，灵活扩展为其他组合优化问题的求解模板。
1. 混合粒子群算法求解TSP：从连续搜索到排列空间
当你在Matlab里用randperm生成路径，再用for循环计算总距离时，很快会撞上一堵墙：标准粒子群算法的速度和位置公式全部建立在实数向量空间上，而TSP问题要求输出的是合法的城市排列。排列不是向量，对城市编号做加减乘除得到的大概率不是可行路径。混合粒子群算法就是为这类离散排列问题准备的：它在粒子群框架里引入交叉、变异和2-opt局部搜索，用遗传算子模拟速度合成，再保留PSO的全局协作引导能力。这里给出一套带中文注释的Matlab代码实例，讲清楚每个参数怎么设、每段代码在解决什么，并讨论收敛验证和加速技巧。适合不依赖matlab优化工具箱、想手写优化算法并理解内部机制的工程师。
2. 混合粒子群算法的核心机制：编码、适应度与交换序运算
动手写代码之前，必须先把“位置”和“速度”在排列空间里的意义定义清楚。如果这一步不落地，后续的混合策略只会变成随机搜索。这一章先从TSP场景解释三个基础概念：排列编码、路径适应度、交换序，再给出三种常见混合策略，为Matlab实现铺路。
2.1 排列编码与适应度函数：从城市坐标到路径长度
TSP问题的标准输入是一个城市坐标矩阵，输出是一条周游所有城市并回到起点的闭环路径。在Matlab中，我把路径编码为1×N的整数排列，排列的每一位对应一个城市编号。假设路径表示从城市出发，按顺序经过后续城市，最后回到。这种编码简洁，适应度函数就是路径上的欧氏距离总和。
计算路径长度最直接的写法是循环累加，但Matlab更适合用向量化方式一次性求和。下面的函数同时处理了“回到起点”这一步：
function total = calcPathLength(path, distMatrix) % 计算闭合环路长度，输入distMatrix为n x n距离矩阵 idx = [path, path(1)]; % 末尾补起点，形成闭环 row = idx(1:end-1); col = idx(2:end); linIdx = sub2ind(size(distMatrix), row, col); total = sum(distMatrix(linIdx)); end代码逻辑很直接：把路径变成相邻城市对，用sub2ind把行列下标映射到距离矩阵的线性索引，最后求和。这里把distMatrix单独抽出来作为参数，意味着你可以传入欧氏距离、曼哈顿距离甚至任意业务自定义权重矩阵。需要注意的是，如果城市坐标的量纲差异很大，不要直接计算欧氏距离，先对坐标做归一化，否则距离矩阵会被某一个坐标维度主导，这个细节会在最后一章展开。
2.2 交换序：在排列空间上定义“速度”
标准PSO速度更新公式里有三个关键操作：保留惯性、向个体历史最优靠拢、向全局历史最优靠拢。在连续空间里，这些操作都是基于向量加减和实数系数相乘。排列空间没有线性结构，所以理论上有一种做法叫交换序。
举例说明。假设当前排列，而某个粒子的历史最优排列。这两个排列的差异可以用一个交换对来表示：交换的第1位和第2位，就得到。这个交换对就是“从当前排列到目标排列的速度”。如果差异不止一处，就用一组交换对表示。位置更新就是依次执行这一组交换对。这样一来，PSO的速度更新可以映射成下面这张表：
| 标准PSO概念 | TSP排列空间含义 |
|---|---|
| 位置向量x | 1×N城市排列 |
| 速度v | 由若干交换对组成的cell数组 |
| x+v | 依次执行v中的交换对 |
| pbest-x | 把x改成pbest所需的一组交换对 |
| w·v | 按概率保留速度中的部分交换对 |
| c·r·(pbest-x) | 按概率抽取pbest-x中的交换对加入速度 |
从表里可以看到，理论上的纯交换序PSO并不需要引入遗传算子。但实际实现中会有两个麻烦：从x到pbest的交换集合并不唯一，组合数很大；多个交换对叠加后可能产生无效操作。因此行业里的常见做法是用遗传算法的交叉算子来替代“相减”和“速度合成”，保留PSO的全局最优引导骨架，这就是“混合”二字的来源。
2.3 混合策略：交叉、变异与局部搜索的组合方式
混合粒子群算法不是一个固定算法，而是一个算法族。我会根据问题规模在以下三种组合里选择：
| 混合类型 | 具体做法 | 适用场景 |
|---|---|---|
| PSO+遗传交叉 | 粒子位置与pbest、gbest做顺序交叉，再以变异概率做扰动 | 城市数量50以下，代码简单 |
| PSO+2-opt | 每固定代数对gbest执行2-opt局部搜索 | 收敛快，中等规模TSP |
| PSO+模拟退火 | 接受差解并用退火温度控制扰动幅度 | 城市多、易早熟的场景 |
这三种可以叠加使用。后续代码里选择PSO+顺序交叉+2-opt的组合：交叉算子负责全局搜索方向的传递，2-opt负责把最终插队结果打磨成更短路径。这个组合不追求理论纯度，但在Matlab里实现清晰，调参也容易预测。选择交叉算子时，顺序交叉能保留父本中一段区间的相对顺序，比部分映射交叉更容易保留环路结构，对TSP更友好。
3. 带注释的Matlab代码实例：主循环与混合策略实现
这一章直接给可运行的代码。所有函数都放进同一个文件hybrid_PSO_TSP.m里，Matlab R2016b之后支持在函数文件末尾追加局部函数，读者只需按顺序复制保存，不需要额外添加路径。
3.1 主程序文件：参数接口与距离矩阵计算
主函数参数设计成结构体，便于后续批量实验。代码第一段负责读参数、计算距离矩阵：
function [bestPath, bestLen] = hybrid_PSO_TSP(cityXY, params) % 混合粒子群算法求解TSP，带中文注释 % 输入： % cityXY - n x 2 矩阵，第i行存放第i个城市坐标 % params - 结构体，字段如下： % popSize : 种群数量 % maxIter : 最大迭代代数 % wStart : 惯性权重起始值 % wEnd : 惯性权重结束值 % c1 : 个体学习概率 % c2 : 全局学习概率 % pC : 交叉概率 % pM : 变异概率 % use2opt : 是否启用2-opt局部搜索 % 输出： % bestPath - 最短路径城市排列 % bestLen - 最短路径长度 n = size(cityXY, 1); dx = cityXY(:,1) - cityXY(:,1)'; % 两两城市x坐标差矩阵 dy = cityXY(:,2) - cityXY(:,2)'; % 两两城市y坐标差矩阵 distMat = sqrt(dx.^2 + dy.^2); % 欧氏距离矩阵这里没有用matlab优化工具箱的pdist2，因为pdist2属于统计和机器学习工具箱，手写坐标差矩阵可以让基础版Matlab直接运行。dx是n×n矩阵，第i行第j列表示城市i与城市j的x坐标差，配合dy的平方和开根号得到完整距离矩阵。
params字段与默认建议值可以对照下面这张表：
| params字段 | 语义 | 建议取值 |
|---|---|---|
| popSize | 种群大小 | 20~60 |
| maxIter | 最大迭代次数 | 100~500 |
| wStart | 惯性权重初值 | 0.9 |
| wEnd | 惯性权重终值 | 0.4 |
| c1 | 个体学习概率 | 0.7 |
| c2 | 全局学习概率 | 0.8 |
| pC | 交叉算子内部概率 | 0.9 |
| pM | 变异概率 | 0.1 |
| use2opt | 是否启用2-opt局部搜索 | 1或0 |
3.2 种群初始化与个体历史最优
接着初始化粒子种群。每个粒子用randperm生成一个随机排列。速度不单独保存，因为混合策略里用交叉算子代替了显式速度：
pop = cell(params.popSize, 1); for i = 1:params.popSize pop{i} = randperm(n); end pbestPop = pop; % 个体历史最优位置 pbestLen = inf(params.popSize, 1); % 个体历史最优距离 for i = 1:params.popSize pbestLen(i) = calcPathLength(pop{i}, distMat); end [gbestLen, gidx] = min(pbestLen); % 全局最优距离 gbestPop = pbestPop{gidx}; % 全局最优路径 bestHistory = zeros(params.maxIter, 1);把个体最优直接复制一份作为pbestPop，内存开销很小，因为每组排列只是1×n整数向量。calcPathLength函数在第2章已经定义，完整文件末尾需要再放一份。
3.3 主循环：惯性权重递减与粒子更新
主循环里每一代做四件事：更新惯性权重、更新每个粒子、更新个体最优与全局最优、周期执行2-opt。代码如下：
for iter = 1:params.maxIter w = params.wStart - (params.wStart - params.wEnd) * iter / params.maxIter; for i = 1:params.popSize new = psoUpdate(pop{i}, pbestPop{i}, gbestPop, ... w, params.c1, params.c2, params.pC, params.pM); pop{i} = new; len = calcPathLength(new, distMat); if len < pbestLen(i) pbestLen(i) = len; pbestPop{i} = new; end if len < gbestLen gbestLen = len; gbestPop = new; end end if params.use2opt && mod(iter, 5) == 0 [gbestPop, gbestLen] = twoOptMove(gbestPop, distMat, gbestLen); end bestHistory(iter) = gbestLen; end bestPath = gbestPop; bestLen = gbestLen; endw按线性策略从wStart衰减到wEnd。前期w接近0.9，粒子保留较多自身路径样本，种群维持多样性；后期w接近0.4，粒子更容易接受pbest与gbest引导的交叉，收敛速度加快。mod(iter,5)==0让2-opt每五代执行一次，而不是每代执行，降低计算量并防止过早陷入局部最优。bestHistory保存每一代全局最优长度，用于绘制收敛曲线或判断早熟。
3.4 粒子更新函数：交叉、变异与惯性保留
psoUpdate是混合粒子群的核心更新函数。它用随机数把更新过程分成三部分：惯性保留、向个体最优点交叉、向全局最优点交叉。这里是完整的核心代码：
function new = psoUpdate(x, pbest, gbest, w, c1, c2, pC, pM) % 混合更新：用顺序交叉和交换变异实现PSO的速度合成 r1 = rand; if r1 < w new = x; % 保留惯性路径 else new = swapMutate(x); % 惯性不足时先做一次交换变异 end if rand < c1 new = orderCrossover(new, pbest, pC); % 向个体历史最优学习 end if rand < c2 new = orderCrossover(new, gbest, pC); % 向全局最优学习 end if rand < pM new = swapMutate(new); % 最终变异扰动 end end这里把c1和c2当作交叉触发概率，而不是连续PSO中的系数，所以取值范围建议控制在0.6~0.9。如果c1超过1，rand始终小于它，等于每次都强制与pbest交叉，会破坏粒子的自身结构。w作为惯性保留概率，当惯性分支被触发但后续又发生交叉时，原路径也会被修改，这正是速度合成与遗传操作结合的效果。
3.5 顺序交叉、交换变异与2-opt局部搜索
顺序交叉是TSP中很常用的交叉算子，它从第二个父本选取一个子区间，按顺序填补到第一个父本中，保留两个父本的顺序特性：
function child = orderCrossover(parent, donor, pC) % parent与donor都是1 x n城市排列，pC为交叉概率 n = length(parent); if rand > pC child = parent; return; end point1 = randi(n-2); % 随机交叉起点 point2 = point1 + randi(n-point1); % 随机交叉终点 child = -ones(1,n); % 先用-1占位 child(point1:point2) = parent(point1:point2); % 用donor中未出现的城市按顺序填充其余空位 j = point2 + 1; for i = 1:n k = mod(point2 + i - 1, n) + 1; % 从donor交叉点后开始循环 if ~ismember(donor(k), child) if j > n j = 1; end child(j) = donor(k); j = j + 1; end end endchild保留parent在[point1, point2]区间的片段，剩余位置用donor中未出现的城市按顺序填充，保证子代是合法排列。ismember在n较小时性能可接受，城市超过200个时建议改成逻辑数组记录已使用城市。
交换变异很简单：
function x = swapMutate(x) % 交换变异：随机交换两个位置的城市编号 idx = randperm(length(x), 2); x([idx(1) idx(2)]) = x([idx(2) idx(1)]); end2-opt局部搜索用来做路径精化，对路径中两条不相邻的边进行断开重连：
function [path, bestLen] = twoOptMove(path, distMat, bestLen) % 2-opt：尝试反转一段子路径，若路径变短则接受 n = length(path); improved = true; while improved improved = false; for i = 1:n-2 for j = i+1:n if j == n continue; end a = path(i); b = path(i+1); c = path(j); d = path(j+1); delta = - distMat(a,b) - distMat(c,d) ... + distMat(a,c) + distMat(b,d); if delta < -1e-6 path(i+1:j) = path(j:-1:i+1); % 反转i+1到j段 bestLen = bestLen + delta; improved = true; end end end end enddelta是交换前后的路径长度变化量，只有负改善才接受。j==n时continue，代表不处理包含回到起点那条边的反转，这是简化实现，对30~100城市的问题足够稳定。需要把calcPathLength、psoUpdate、orderCrossover、swapMutate、twoOptMove这些局部函数按顺序放在hybrid_PSO_TSP.m文件末尾，主函数结束时需要end，每个局部函数也要有对应的end。
使用示例：
rng(1); cityXY = 100 * rand(30, 2); % 30个城市 params = struct('popSize',40,'maxIter',200,'wStart',0.9,'wEnd',0.4,... 'c1',0.7,'c2',0.8,'pC',0.9,'pM',0.1,'use2opt',1); [bestPath, bestLen] = hybrid_PSO_TSP(cityXY, params); disp(bestLen);运行后可以看到bestLen随迭代下降。如果始终不下降，优先检查距离矩阵是否对称，再检查orderCrossover生成的子代是否合法。
4. 参数设定与收敛性分析：惯性权重、学习因子与局部搜索强度
很多读者把代码跑通后的第一件事就是改参数。混合粒子群算法可调参数多，随意组合很容易出现“跑很久不如随机搜索”的现象。这一章从工程角度给出调参方向，并解释每个参数为什么会影响收敛。
4.1 参数速查表与默认建议
| 参数 | 推荐范围 | 对结果的影响 |
|---|---|---|
| popSize | 20~60 | 种群太小容易早熟，太大收敛慢 |
| maxIter | 100~500 | 城市数量增加时按倍数增长 |
| wStart | 0.8~1.0 | 初始惯性保留概率，大则多样性好 |
| wEnd | 0.3~0.5 | 后期收敛速度，过小会快速收敛到局部最优 |
| c1 | 0.6~0.9 | 向个体历史最优交叉的概率 |
| c2 | 0.6~0.9 | 向全局最优交叉的概率 |
| pC | 0.7~0.95 | 顺序交叉执行概率 |
| pM | 0.05~0.2 | 变异概率，越大越容易跳出局部最优 |
| use2opt间隔 | 5~10代 | 间隔太短会压制全局搜索 |
默认值针对30~50城市。如果城市数量增加到100，popSize建议调整到60~100，maxIter增加到300~1000。
4.2 惯性权重w的线性递减策略对收敛的影响
代码里使用wStart到wEnd的线性衰减：
w = params.wStart - (params.wStart - params.wEnd) * iter / params.maxIter;这个策略在连续PSO中被广泛验证，在混合PSO里同样有效。迭代初期，w接近0.9，粒子保留较多自身路径，种群拥有足够探索空间；迭代后期，w接近0.4，粒子更愿意接受pbest与gbest引导的交叉，快速向当前最佳区域收敛。
这里有一个容易踩的坑：wEnd不要设为0。如果后期w变成0，惯性分支完全消失，粒子一旦被gbest同质化，就再也没有路径保留机制，种群会迅速收敛到局部最优。即使pM不为0，靠少量变异很难找回多样性。wEnd建议不低于0.3。
4.3 学习因子c1/c2与交叉概率pC的配合原则
在本文的代码中，c1和c2是交叉触发概率。c1=0.7、c2=0.8意味着每次更新有70%的概率与个体历史最优交叉，80%的概率与全局最优交叉，两次触发相互独立。c2略大于c1，体现全局最优更强的引导作用。如果c1和c2都超过0.95，粒子会在每一次迭代中都同时与pbest和gbest交叉，自身路径被反复拆解，收敛曲线容易出现剧烈震荡。
pC和c1/c2是嵌套关系。pC控制交叉算子内部是否真正执行，c1/c2控制是否发起与某个父本的交叉。例如c1=0.7、pC=0.9，实际发生有效个体交叉的概率约0.63。如果交叉太频繁，可以优先降低c2，让粒子更多依靠自身历史信息修正路径，而不是被全局最优过早拉拢。
4.4 局部搜索强度与早熟判定的工程技巧
2-opt局部搜索的计算量是O(n²)，每5代执行一次。在30城市问题上，加入2-opt后总耗时增加不到2倍，但最终结果通常能缩短5%~15%。每代都做2-opt并不好，因为局部搜索过强会让种群迅速集中到gbest附近，其他粒子提供的多样性被浪费。间隔5~10代是比较稳妥的折衷。
早熟判定可以这样实现：记录bestHistory连续多少代没有变化，超过阈值就执行粒子重启：
COUNTER_LIMIT = 30; if iter > COUNTER_LIMIT && bestHistory(iter) == bestHistory(iter-COUNTER_LIMIT) [~, gidx] = min(pbestLen); % 当前gbest所在粒子 for i = 1:params.popSize if i ~= gidx pop{i} = swapMutate(pop{i}); end end end这段代码放在主循环内部，通过min(pbestLen)动态找到当前全局最优所在粒子，并在重启时跳过它。交换变异能打散非最优粒子的路径结构，保留全局最优位置不被破坏。需要注意，重启后w仍按原线性衰减策略继续下降，如果重启发生在迭代后期，w已经很低，重新探索效果有限。更复杂做法是重启时把iter拉回较小的值，或者临时重置w，这属于迭代调度范畴。
5. 进阶：归一化距离矩阵、多次验证与退火式重启
5.1 距离矩阵归一化：避免量纲偏置
如果城市坐标来自经纬度、像素坐标或业务数据，x和y的量纲可能相差很大，直接算欧氏距离会让距离矩阵被较大的坐标维度主导。常见做法是在计算距离矩阵之前，对cityXY做Min-Max归一化到[0,1]区间。注意不要对距离矩阵本身归一化，而是要先把坐标归一化，再算距离，否则城市间的相对位置关系会被压缩变形。
cityXY = (cityXY - min(cityXY)) ./ (max(cityXY) - min(cityXY));这里的max和min作用于整个矩阵，把坐标压缩到[0,1]。如果城市分布有长尾异常值，可以先裁剪最多5%的极值再做归一化，避免某一个极端坐标把其他城市之间的距离全部压到几乎为0。
5.2 用多组随机初始化验证算法稳定性
混合粒子群算法是随机算法，单次运行结果可能是运气成分。至少做20次独立实验，统计最优长度的均值和标准差，同时记录bestHistory的收敛曲线。如果均值距离已知最优解较远，不要直接加倍popSize，先观察收敛曲线是否早期进入平台。平台期出现得太早，说明w下降过快或c2过高，优先把c2调低0.05重新实验。
5.3 模拟退火式重启的第二种实现
第4章的早熟重启只做交换变异，这里给出更接近模拟退火的做法：当bestHistory连续40代无改善时，用温度T控制粒子打乱程度。温度随迭代次数降低，后期打乱概率也降低：
if iter > 40 && bestHistory(iter) == bestHistory(iter-40) T = 0.6 * (1 - iter / params.maxIter) + 0.1; for i = 1:params.popSize if rand < T for k = 1:ceil(n/10) pop{i} = swapMutate(pop{i}); end end end endT从0.7左右线性衰减到0.1，表示后期即使判定早熟，也保留10%的扰动概率。这种重启比直接重置w更温和，不会把已经获得的好路径片段全部打散。重启后的w建议设在0.9而不是从0.4开始，这是二次探索成功的常见前提。
本文还有配套的精品资源，点击获取
