混合粒子群算法求解TSP的Matlab实现与参数调优

简介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排列空间含义位置向量x1×N城市排列速度v由若干交换对组成的cell数组xv依次执行v中的交换对pbest-x把x改成pbest所需的一组交换对w·v按概率保留速度中的部分交换对c·r·(pbest-x)按概率抽取pbest-x中的交换对加入速度从表里可以看到理论上的纯交换序PSO并不需要引入遗传算子。但实际实现中会有两个麻烦从x到pbest的交换集合并不唯一组合数很大多个交换对叠加后可能产生无效操作。因此行业里的常见做法是用遗传算法的交叉算子来替代“相减”和“速度合成”保留PSO的全局最优引导骨架这就是“混合”二字的来源。2.3 混合策略交叉、变异与局部搜索的组合方式混合粒子群算法不是一个固定算法而是一个算法族。我会根据问题规模在以下三种组合里选择混合类型具体做法适用场景PSO遗传交叉粒子位置与pbest、gbest做顺序交叉再以变异概率做扰动城市数量50以下代码简单PSO2-opt每固定代数对gbest执行2-opt局部搜索收敛快中等规模TSPPSO模拟退火接受差解并用退火温度控制扰动幅度城市多、易早熟的场景这三种可以叠加使用。后续代码里选择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~60maxIter最大迭代次数100~500wStart惯性权重初值0.9wEnd惯性权重终值0.4c1个体学习概率0.7c2全局学习概率0.8pC交叉算子内部概率0.9pM变异概率0.1use2opt是否启用2-opt局部搜索1或03.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超过1rand始终小于它等于每次都强制与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 i1:n if j n continue; end a path(i); b path(i1); c path(j); d path(j1); delta - distMat(a,b) - distMat(c,d) ... distMat(a,c) distMat(b,d); if delta -1e-6 path(i1:j) path(j:-1:i1); % 反转i1到j段 bestLen bestLen delta; improved true; end end end end enddelta是交换前后的路径长度变化量只有负改善才接受。jn时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 参数速查表与默认建议参数推荐范围对结果的影响popSize20~60种群太小容易早熟太大收敛慢maxIter100~500城市数量增加时按倍数增长wStart0.8~1.0初始惯性保留概率大则多样性好wEnd0.3~0.5后期收敛速度过小会快速收敛到局部最优c10.6~0.9向个体历史最优交叉的概率c20.6~0.9向全局最优交叉的概率pC0.7~0.95顺序交叉执行概率pM0.05~0.2变异概率越大越容易跳出局部最优use2opt间隔5~10代间隔太短会压制全局搜索默认值针对30~50城市。如果城市数量增加到100popSize建议调整到60~100maxIter增加到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是交叉触发概率。c10.7、c20.8意味着每次更新有70%的概率与个体历史最优交叉80%的概率与全局最优交叉两次触发相互独立。c2略大于c1体现全局最优更强的引导作用。如果c1和c2都超过0.95粒子会在每一次迭代中都同时与pbest和gbest交叉自身路径被反复拆解收敛曲线容易出现剧烈震荡。pC和c1/c2是嵌套关系。pC控制交叉算子内部是否真正执行c1/c2控制是否发起与某个父本的交叉。例如c10.7、pC0.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开始这是二次探索成功的常见前提。本文还有配套的精品资源点击获取