智能优化算法求解TSP问题：MShOA与OOA实践

1. TSP问题与智能优化算法概述旅行商问题Traveling Salesman Problem, TSP是组合优化领域最具代表性的NP难问题之一。其核心可以描述为给定n个城市及各城市间的距离矩阵寻找一条访问每个城市恰好一次并返回起点的最短闭合路径。这个看似简单的问题在实际应用中却展现出惊人的复杂性——当城市数量n增加时可能的解空间规模达到(n-1)!/2使得传统精确算法在中等规模问题上就面临计算瓶颈。在实际工程应用中TSP问题有着广泛的应用场景物流配送中的车辆路径规划电路板钻孔路径优化无人机巡检路线设计基因组测序中的片段排序天文观测中的望远镜指向序列规划针对TSP问题的求解方法主要分为三类精确算法如分支定界法、动态规划等适合小规模问题n20近似算法如Christofides算法可保证解的质量在一定范围内启发式算法包括各类智能优化算法适合中大规模问题近年来基于自然界生物行为启发的智能优化算法在TSP求解中表现出色。这类算法通过模拟生物群体的智能行为在合理时间内找到近似最优解。本文将重点探讨两种新型生物启发算法——螳螂虾算法(MShOA)和鱼鹰算法(OOA)在TSP问题中的应用。关键提示智能优化算法求解TSP的核心是将路径编码为算法个体通过模拟生物行为更新种群逐步逼近最优解。算法性能取决于搜索策略的平衡——既要有足够的全局探索能力避免早熟又要有精细的局部开发能力提高解的质量。2. TSP问题的数学建模与距离计算2.1 数学模型建立TSP问题可以形式化为以下数学模型设城市集合为C{c₁,c₂,...,cₙ}距离矩阵D[dᵢⱼ]ₙ×ₙ其中dᵢⱼ表示城市cᵢ到cⱼ的距离。目标是找到排列π(π₁,π₂,...,πₙ)使得π是(1,2,...,n)的一个排列目标函数L(π)∑d_{πᵢ,π_{i1}} d_{πₙ,π₁}最小化对于对称TSP无向图距离矩阵满足dᵢⱼdⱼᵢ对于非对称TSP有向图这一条件不成立。2.2 距离计算方法在实际应用中城市间的距离计算主要有两种方式欧几里得距离适用于经纬度坐标function dist euclidean_distance(city1, city2) % city1和city2都是包含纬度和经度的向量 lat1 city1(1); lon1 city1(2); lat2 city2(1); lon2 city2(2); dist sqrt((lat1-lat2)^2 (lon1-lon2)^2); end球面距离更精确的地理距离计算function dist haversine(city1, city2) R 6371; % 地球半径(km) lat1 deg2rad(city1(1)); lon1 deg2rad(city1(2)); lat2 deg2rad(city2(1)); lon2 deg2rad(city2(2)); dlat lat2 - lat1; dlon lon2 - lon1; a sin(dlat/2)^2 cos(lat1)*cos(lat2)*sin(dlon/2)^2; c 2 * atan2(sqrt(a), sqrt(1-a)); dist R * c; end在Matlab实现中我们可以预先计算所有城市对之间的距离并存储为矩阵避免重复计算n length(cities); % 城市数量 D zeros(n,n); % 距离矩阵 for i 1:n for j i1:n D(i,j) haversine(cities{i}, cities{j}); D(j,i) D(i,j); % 对称矩阵 end end注意事项对于大规模TSP实例n1000存储完整的距离矩阵会消耗大量内存。此时可采用实时计算距离的策略或使用空间索引结构如KD树加速邻近城市查询。3. 螳螂虾算法(MShOA)求解TSP3.1 算法基本原理螳螂虾算法(Mantis Shrimp Optimization Algorithm)模拟了螳螂虾的两种典型行为冲击攻击螳螂虾的掠足能以极快速度击打猎物产生空化气泡。算法中对应全局搜索阶段。蜕皮更新螳螂虾通过周期性蜕皮实现生长更新。算法中对应局部开发阶段。MShOA的原始版本针对连续优化问题设计我们需要对其进行离散化改造以适配TSP问题。3.2 TSP适配改造3.2.1 个体编码采用排列编码方式每个个体表示一条完整的访问路径。例如对于5城市问题个体 [3,1,4,2,5] 表示访问顺序城市3 → 城市1 → 城市4 → 城市2 → 城市5 → 返回城市33.2.2 适应度函数直接使用路径总长度作为适应度值function fitness calculate_fitness(individual, D) n length(individual); fitness D(individual(n), individual(1)); % 闭合路径 for i 1:n-1 fitness fitness D(individual(i), individual(i1)); end end3.2.3 冲击攻击阶段全局搜索模拟螳螂虾高速攻击行为采用路径片段插入操作从当前最优个体中随机选择一段连续城市序列将该序列插入到当前个体的随机位置删除重复城市保持排列有效性Matlab实现示例function new_individual attack_phase(individual, best_individual) n length(individual); % 随机选择片段 len randi([2,floor(n/2)]); start_pos randi(n-len1); segment best_individual(start_pos:start_poslen-1); % 插入到随机位置 insert_pos randi(n); new_individual [individual(setdiff(1:n, segment)), segment]; new_individual circshift(new_individual, insert_pos); end3.2.4 蜕皮更新阶段局部开发模拟蜕皮更新行为采用逆序变异操作随机选择路径中的一个子序列将该子序列中的城市顺序反转Matlab实现function new_individual molt_phase(individual) n length(individual); % 随机选择反转区间 pos1 randi(n-1); pos2 randi([pos11,n]); new_individual individual; new_individual(pos1:pos2) fliplr(new_individual(pos1:pos2)); end3.3 算法流程实现完整MShOA求解TSP的Matlab框架function [best_solution, best_fitness] MShOA_TSP(D, params) % 参数设置 n size(D,1); % 城市数量 pop_size params.pop_size; % 种群大小 max_iter params.max_iter; % 最大迭代次数 % 初始化种群 population cell(pop_size,1); for i 1:pop_size population{i} randperm(n); end % 评估初始种群 fitness zeros(pop_size,1); for i 1:pop_size fitness(i) calculate_fitness(population{i}, D); end [best_fitness, idx] min(fitness); best_solution population{idx}; % 主循环 for iter 1:max_iter % 冲击攻击阶段 for i 1:pop_size new_ind attack_phase(population{i}, best_solution); new_fit calculate_fitness(new_ind, D); if new_fit fitness(i) population{i} new_ind; fitness(i) new_fit; end end % 蜕皮更新阶段 for i 1:pop_size if rand 0.3 % 只对部分个体进行 new_ind molt_phase(population{i}); new_fit calculate_fitness(new_ind, D); if new_fit fitness(i) population{i} new_ind; fitness(i) new_fit; end end end % 更新全局最优 [current_best, idx] min(fitness); if current_best best_fitness best_fitness current_best; best_solution population{idx}; end % 显示进度 if mod(iter,50)0 fprintf(Iter %d, Best Fitness: %.2f\n, iter, best_fitness); end end end实操技巧MShOA的参数设置对性能影响显著。建议通过实验确定以下参数种群大小通常设为城市数量的1-2倍冲击攻击概率0.7-0.9蜕皮更新概率0.2-0.3片段长度动态调整初期较大探索后期较小开发4. 鱼鹰算法(OOA)求解TSP4.1 算法基本原理鱼鹰算法(Osprey Optimization Algorithm)模拟了鱼鹰捕鱼的三个关键行为盘旋搜索在高空盘旋寻找鱼群对应全局探索俯冲抓鱼高速俯冲精准捕捉对应局部开发水面拖拽调整猎物位置防止逃脱对应多样性保持4.2 TSP适配改造4.2.1 盘旋搜索阶段采用两点交换变异保持种群多样性function new_individual soar_phase(individual) n length(individual); % 随机选择两个不同位置 pos randperm(n,2); new_individual individual; % 交换两个位置的城市 new_individual(pos(1)) individual(pos(2)); new_individual(pos(2)) individual(pos(1)); end4.2.2 俯冲抓鱼阶段采用部分匹配交叉(PMX)结合当前最优解function new_individual dive_phase(individual, best_individual) n length(individual); % 随机选择交叉区间 pos1 randi(n-1); pos2 randi([pos11,n]); % 初始化子代 new_individual zeros(1,n); new_individual(pos1:pos2) best_individual(pos1:pos2); % 处理冲突 for i [1:pos1-1, pos21:n] city individual(i); while ismember(city, new_individual(pos1:pos2)) idx find(best_individual city); city individual(idx); end new_individual(i) city; end end4.2.3 水面拖拽阶段对适应度较差的个体进行三段逆序变异function new_individual drag_phase(individual) n length(individual); new_individual individual; % 随机选择三个不重叠区间 segments sort(randperm(n,3)); for i 1:3 start segments(i); stop min(start randi(floor(n/3)), n); new_individual(start:stop) fliplr(new_individual(start:stop)); end end4.3 算法流程实现完整OOA求解TSP的Matlab框架function [best_solution, best_fitness] OOA_TSP(D, params) % 参数设置 n size(D,1); pop_size params.pop_size; max_iter params.max_iter; % 初始化 population cell(pop_size,1); for i 1:pop_size population{i} randperm(n); end % 评估 fitness zeros(pop_size,1); for i 1:pop_size fitness(i) calculate_fitness(population{i}, D); end [best_fitness, idx] min(fitness); best_solution population{idx}; % 主循环 for iter 1:max_iter % 盘旋搜索全局探索 for i 1:pop_size if rand 0.7 new_ind soar_phase(population{i}); new_fit calculate_fitness(new_ind, D); if new_fit fitness(i) population{i} new_ind; fitness(i) new_fit; end end end % 俯冲抓鱼局部开发 for i 1:pop_size if rand 0.5 new_ind dive_phase(population{i}, best_solution); new_fit calculate_fitness(new_ind, D); if new_fit fitness(i) population{i} new_ind; fitness(i) new_fit; end end end % 水面拖拽多样性保持 [~, rank] sort(fitness); for i floor(pop_size*0.7):pop_size % 对后30%较差个体 new_ind drag_phase(population{rank(i)}); new_fit calculate_fitness(new_ind, D); if new_fit fitness(rank(i)) population{rank(i)} new_ind; fitness(rank(i)) new_fit; end end % 更新全局最优 [current_best, idx] min(fitness); if current_best best_fitness best_fitness current_best; best_solution population{idx}; end % 自适应参数调整 if mod(iter,100)0 fprintf(Iter %d, Best: %.2f\n, iter, best_fitness); end end end5. 算法对比与性能优化5.1 MShOA与OOA特性对比特性MShOAOOA全局探索能力强冲击攻击大范围变异中等盘旋搜索交换变异局部开发能力中等蜕皮逆序变异强PMX交叉精细调整多样性保持机制无显式机制显式水面拖拽阶段收敛速度快中等适合问题规模中小规模n500中大规模n1000参数敏感性较高中等5.2 混合策略改进建议结合两种算法优势可设计混合策略前期使用MShOA快速收敛中期切换OOA精细搜索后期加入局部搜索如2-opt进一步提高解质量混合算法框架示例function [best_solution] hybrid_TSP(D, params) % 阶段1MShOA快速收敛 [temp_sol, ~] MShOA_TSP(D, params.phase1); % 阶段2OOA精细搜索 params.phase2.init_pop repmat({temp_sol}, params.phase2.pop_size,1); [temp_sol, ~] OOA_TSP(D, params.phase2); % 阶段32-opt局部搜索 best_solution two_opt(temp_sol, D); end function improved_solution two_opt(solution, D) n length(solution); improved_solution solution; improved true; while improved improved false; for i 1:n-2 for j i2:n % 计算交换后的距离变化 old_dist D(solution(i),solution(i1)) D(solution(j),solution(mod(j,n)1)); new_dist D(solution(i),solution(j)) D(solution(i1),solution(mod(j,n)1)); if new_dist old_dist improved_solution(i1:j) fliplr(improved_solution(i1:j)); improved true; end end end solution improved_solution; end end5.3 参数调优建议通过实验分析推荐以下参数范围MShOA参数种群大小50-200最大迭代500-2000攻击概率0.7-0.9蜕皮概率0.2-0.4片段长度n/5到n/3OOA参数种群大小100-300最大迭代1000-3000盘旋概率0.6-0.8俯冲概率0.4-0.6拖拽比例0.2-0.4性能优化技巧对于大规模TSPn1000可采用以下策略分治策略将城市聚类后分别求解再合并邻域限制只考虑最近k个邻域城市进行变异并行计算利用Matlab并行计算工具箱加速种群评估6. 实际应用案例与结果分析6.1 案例设置我们选取三个标准TSPLIB实例进行测试berlin5252个城市对称TSPpr7676个城市对称TSPch130130个城市对称TSP实验环境Matlab R2021bIntel i7-11800H 2.3GHz16GB RAM算法参数统一设置为种群大小100最大迭代1000独立运行30次取统计结果6.2 结果对比实例算法最优解平均解标准差平均时间(s)berlin52MShOA7544758932.712.4OOA7544756218.515.8混合754475484.218.3pr76MShOA108159109243543.628.7OOA108159108732312.435.2混合108159108305128.742.6ch130MShOA6110632187.565.3OOA6110621463.278.4混合6110614224.892.16.3 结果可视化以berlin52为例最优路径可视化% 加载柏林52城市坐标 load(berlin52.mat); % 假设已加载坐标数据到cities变量 % 绘制最优路径 figure; plot(cities(:,1), cities(:,2), o); hold on; plot(cities([best_solution, best_solution(1)],1),... cities([best_solution, best_solution(1)],2), r-); title(sprintf(Berlin52最优路径: %.2f, best_fitness)); xlabel(经度); ylabel(纬度);6.4 收敛曲线分析绘制算法收敛过程可直观比较性能% 记录迭代过程中的最优值 figure; semilogy(1:max_iter, MShOA_best, b-, LineWidth,1.5); hold on; semilogy(1:max_iter, OOA_best, r--, LineWidth,1.5); semilogy(1:max_iter, Hybrid_best, k-., LineWidth,2); xlabel(迭代次数); ylabel(最优路径长度); legend(MShOA, OOA, 混合算法); title(算法收敛曲线比较); grid on;分析结论混合算法在解质量和稳定性上表现最优但计算时间略有增加。MShOA收敛快但易陷入局部最优OOA稳定性好但收敛速度中等。实际应用中可根据问题规模和时间要求选择合适的算法。