模拟退火算法原理与Matlab实现：以旅行商问题为例

1. 从一个“退火”的物理过程说起如果你在搜索引擎里输入“模拟退火算法”大概率会看到一堆关于“退火”的物理过程描述金属加热到高温然后缓慢冷却原子从高能无序状态逐渐趋于低能稳定状态。这个比喻很经典但第一次接触时我总觉得它离我们写代码解决实际问题有点远。直到有一次我面对一个复杂的排产优化问题试遍了贪心、遗传算法效果都不理想不是陷入局部最优解就是收敛太慢才真正回过头来琢磨这个听起来有点“玄学”的算法。简单来说模拟退火算法Simulated Annealing, SA是一种启发式随机搜索算法它的核心思想就是模仿物理退火过程来寻找一个复杂问题的近似全局最优解。为什么是“近似”因为对于很多NP-hard问题比如我们马上要讲的旅行商问题在有限时间内找到绝对的最优解几乎不可能我们能做的是找到一个“足够好”的解。模拟退火厉害的地方在于它允许在搜索过程中“暂时接受”一个更差的解。这个特性让它有概率跳出局部最优的“陷阱”向着全局最优的山峰攀登。今天我们就以经典的旅行商问题Traveling Salesman Problem, TSP作为战场用Matlab这把熟悉的工具手把手实现一遍模拟退火算法。你会发现它的核心代码可能比你想象的要简洁但其中关于温度调度、邻域搜索、接受准则的每一个参数设置都藏着无数“坑”和经验。这篇文章我会把我从理论到实践再到调参优化过程中踩过的雷、总结的技巧毫无保留地分享给你。2. 旅行商问题一个看似简单却让计算机头疼的经典模型在深入算法之前我们必须先彻底理解我们要打败的“BOSS”——旅行商问题。它的描述简单得令人发指一个旅行商要去N个城市推销商品每个城市只能去一次最后回到起点如何规划路线才能使总路程最短就这么一个问题却成为了组合优化领域的“标杆”。它的解空间有多大呢对于N个城市理论上存在 (N-1)!/2 条不同的哈密顿回路。这个数字增长得有多恐怖5个城市有12条可能路线10个城市就有181440条而30个城市这个数字达到了约 4.4×10^30即使用每秒计算万亿次的超级计算机穷举一遍也需要上百亿年。这就是所谓的“组合爆炸”也是TSP被归为NP-hard问题的原因。在Matlab里我们如何描述一个TSP问题实例最直接的就是一个坐标矩阵。比如我们随机生成10个城市的坐标% 生成10个城市的二维坐标 (范围在[0, 100]内) num_city 10; city_coordinate 100 * rand(num_city, 2);有了坐标任意两个城市i和j之间的距离就可以用欧几里得距离公式计算distance sqrt((x_i - x_j)^2 (y_i - y_j)^2)。我们可以预先计算一个距离矩阵避免在算法中重复计算这是提升效率的关键一步。% 计算距离矩阵 dist_matrix zeros(num_city, num_city); for i 1:num_city for j i1:num_city d norm(city_coordinate(i, :) - city_coordinate(j, :)); dist_matrix(i, j) d; dist_matrix(j, i) d; % 对称矩阵 end end问题的解就是一个城市的访问序列例如[1, 3, 5, 2, 4, 6, 7, 9, 8, 10, 1]最后回到起点。我们的目标函数或称代价函数、能量函数就是这条路线所有相邻城市距离之和。在模拟退火中这个函数值对应着系统的“能量”我们的目标就是让能量最低。注意在实际项目中距离矩阵可能来自真实的道路网络非欧几里得距离或者包含其他代价如时间、费用。算法框架完全通用你只需要替换距离计算部分即可。3. 模拟退火核心三要素温度、扰动与接受理解了问题现在让我们拆解模拟退火这个“黑盒”。它之所以有效依赖于三个核心机制的协同工作我把它比作一个“智能探险家”的登山策略温度 (Temperature)这是算法的“控制阀”。高温时探险家精力充沛系统能量高愿意尝试各种方向甚至走下坡路接受差解以探索更广阔的地形。随着温度降低他越来越“保守”只愿意接受小幅度的上坡或直接下坡最终在低温下稳定在某个低点局部最优或全局最优。温度从初始高温T0按照某个冷却进度表 (Cooling Schedule)逐渐降低至终止温度T_end。状态产生函数 (邻域搜索)这决定了探险家“下一步怎么走”。在TSP中从一个当前路线当前状态如何产生一条新路线新状态常用的邻域操作有交换 (Swap)随机选择两个城市交换它们在序列中的位置。逆转 (Reverse/2-opt)随机选择一段子路径将其顺序完全颠倒。这是TSP中非常高效的一种局部优化操作。插入 (Insert)随机选择一个城市将其插入到序列的另一个随机位置。 这些操作生成了当前解的“邻居”。操作的选择直接影响搜索的效率和效果。Metropolis接受准则这是算法的“灵魂”决定了是否用新解替换旧解。它用一个简单的概率公式表示P exp(-ΔE / T)其中ΔE E_new - E_old是新旧解的目标函数值之差在TSP中就是路径长度的变化。如果ΔE 0新解更优路径更短一定接受。如果ΔE 0新解更差则以概率P接受这个更差的解。 这个准则完美体现了“退火”思想温度T高时即使ΔE很大变得更差很多P也可能较大算法有较大可能“犯错”跳出温度T很低时P趋近于0算法几乎只接受优化解趋于稳定。把这三要素串起来算法的基本骨架就清晰了初始化生成一个初始解如随机路线设定初始温度T0。外循环降温过程当温度T T_end时重复内循环马尔可夫链长度在每个温度下进行L次尝试。对当前解施加一次邻域扰动产生新解。计算目标函数变化ΔE。根据 Metropolis 准则决定是否接受新解。按照冷却进度表降低温度T。输出最终找到的最优解。4. 手把手实现从零编写Matlab代码理论说得再多不如一行代码。我们直接上干货实现一个基础的模拟退火算法求解TSP。我会在关键步骤加上详细注释并分享我调试时的心得。4.1 基础框架搭建首先我们定义问题的规模和参数。%% 1. 问题定义与参数设置 clear; clc; % 城市数量 num_city 20; % 随机生成城市坐标 city_coordinate 100 * rand(num_city, 2); % 模拟退火算法参数 T0 1000; % 初始温度 T_end 1e-8; % 终止温度 alpha 0.99; % 温度衰减系数 (每次乘以alpha) L 100 * num_city; % 马尔可夫链长度每个温度下的迭代次数 % 计算距离矩阵 dist_matrix pdist2(city_coordinate, city_coordinate);这里用了pdist2函数快速计算距离矩阵比自己写循环更简洁高效。参数设置是门艺术T0通常设置为一个能使初始接受概率对差解较高的值。一个经验法则是让初始时exp(-ΔE_avg / T0)接近 1ΔE_avg是随机扰动产生的平均能量差。可以简单设一个较大的数如1000, 10000再调整。alpha通常在0.9到0.999之间。越大降温越慢搜索越细致但耗时越长。L通常与问题规模成正比。100*N是一个常用的经验起点。4.2 核心迭代过程接下来是算法的核心循环。%% 2. 初始化 % 生成初始解随机排列 current_route randperm(num_city); current_distance calculate_distance(current_route, dist_matrix); % 记录最优解 best_route current_route; best_distance current_distance; % 记录迭代过程用于绘图 iter 1; distance_history zeros(1, 10000); % 预分配记录最优距离历史 temperature_history zeros(1, 10000); % 记录温度历史 T T0; while T T_end for i 1:L % 2.1 产生新解邻域搜索 new_route generate_new_route(current_route); new_distance calculate_distance(new_route, dist_matrix); % 2.2 计算能量差 delta_E new_distance - current_distance; % 2.3 Metropolis接受准则 if delta_E 0 % 新解更好直接接受 current_route new_route; current_distance new_distance; % 更新全局最优 if current_distance best_distance best_route current_route; best_distance current_distance; end else % 新解更差以一定概率接受 P exp(-delta_E / T); if rand() P current_route new_route; current_distance new_distance; end % 注意接受差解时不更新全局最优解 end end % 记录数据 distance_history(iter) best_distance; temperature_history(iter) T; iter iter 1; % 2.4 降温 T T * alpha; % 可以添加一个提前终止条件比如连续若干代最优解未改进 end % 截断记录数组 distance_history distance_history(1:iter-1); temperature_history temperature_history(1:iter-1);这里有两个关键的子函数需要实现calculate_distance和generate_new_route。4.3 关键子函数实现计算路径距离函数function total_dist calculate_distance(route, dist_matrix) % 计算一条闭合路径的总距离 n length(route); total_dist 0; for i 1:n-1 total_dist total_dist dist_matrix(route(i), route(i1)); end % 从最后一个城市回到起点 total_dist total_dist dist_matrix(route(n), route(1)); end邻域搜索函数这里采用交换和逆转两种操作function new_route generate_new_route(old_route) % 以一定概率选择不同的邻域操作增加搜索多样性 new_route old_route; n length(old_route); % 操作1交换两个随机城市的位置 (概率0.5) % 操作2逆转一段随机子路径 (概率0.5) if rand() 0.5 % 交换操作 idx randperm(n, 2); % 随机选择两个不同的索引 new_route(idx(1)) old_route(idx(2)); new_route(idx(2)) old_route(idx(1)); else % 逆转操作 (2-opt的一部分) idx sort(randperm(n, 2)); % 随机选择起点和终点并排序 new_route(idx(1):idx(2)) old_route(idx(2):-1:idx(1)); end end实操心得generate_new_route是影响算法性能的关键。纯随机交换虽然简单但效率不高。逆转2-opt操作对于TSP特别有效因为它能直接消除路径中的交叉是局部搜索的利器。在实际代码中我常常会给逆转操作更高的权重比如0.7或者采用更复杂的自适应策略。4.4 可视化与结果分析代码跑完了不看看结果怎么行我们用图形来直观展示算法的收敛过程和最终路线。%% 3. 结果可视化 figure(Position, [100, 100, 1200, 400]); % 子图1优化过程收敛曲线 subplot(1, 3, 1); plot(1:length(distance_history), distance_history, b-, LineWidth, 1.5); xlabel(迭代次数 (外循环)); ylabel(最优路径长度); title(模拟退火优化过程收敛曲线); grid on; % 子图2温度下降曲线 subplot(1, 3, 2); semilogy(1:length(temperature_history), temperature_history, r-, LineWidth, 1.5); % 对数坐标看温度下降 xlabel(迭代次数 (外循环)); ylabel(温度 (对数坐标)); title(温度下降曲线); grid on; % 子图3最优路径可视化 subplot(1, 3, 3); % 将最优路径构造成闭合循环用于绘图 plot_route [best_route, best_route(1)]; plot(city_coordinate(plot_route, 1), city_coordinate(plot_route, 2), ko-, ... LineWidth, 1.5, MarkerSize, 8, MarkerFaceColor, y); hold on; text(city_coordinate(:,1), city_coordinate(:,2), num2str((1:num_city)), ... HorizontalAlignment, center, VerticalAlignment, bottom); xlabel(X坐标); ylabel(Y坐标); title([最优路径图 (总距离: , num2str(best_distance, %.2f), )]); axis equal; grid on; fprintf(初始随机路径长度: %.2f\n, calculate_distance(randperm(num_city), dist_matrix)); fprintf(模拟退火找到的最优路径长度: %.2f\n, best_distance);运行这段代码你会看到三张图一张显示最优路径长度如何随着迭代下降初期可能波动后期趋于平稳一张显示温度如何指数衰减最后一张直观地画出找到的最优旅行路线。5. 参数调优从“能用”到“好用”的关键一跃写完基础版本你可能发现结果有时好有时坏这就是参数调优的战场了。模拟退火被戏称为“炼丹”就是因为参数敏感。下面是我总结的几个核心调优点5.1 初始温度 T0 的设置T0不能拍脑袋定。一个经典的方法是进行一个初始接收率测试随机产生大量如1000次的邻域扰动计算ΔE的绝对值平均值ΔE_avg。设定一个期望的初始接受概率P0例如0.8意味着初始时80%的差解会被接受。根据公式T0 -ΔE_avg / ln(P0)反推初始温度。% 初始温度估计示例 P0 0.8; test_num 1000; delta_E_list zeros(1, test_num); init_route randperm(num_city); init_dist calculate_distance(init_route, dist_matrix); for i 1:test_num test_route generate_new_route(init_route); test_dist calculate_distance(test_route, dist_matrix); delta_E_list(i) abs(test_dist - init_dist); end delta_E_avg mean(delta_E_list); T0_estimated -delta_E_avg / log(P0); fprintf(估计的初始温度 T0: %.2f\n, T0_estimated);5.2 马尔可夫链长度 L 的确定L代表每个温度下搜索的充分程度。太短搜索不充分太长计算耗时。一个自适应策略是在每个温度下直到系统在该温度下达到“准平衡态”再降温。一个可操作的判断是连续接受或拒绝一定次数如10*N的新解后就认为平衡了。我们的简单实现中固定L100*N是一个折中方案对于中小规模问题N100通常够用。5.3 降温策略的选择我们用的是最简单的指数降温T_{k1} alpha * T_k。它简单但降温速度先快后慢。还有其他策略经典模拟退火T_k T0 / (1 k) 其中k是迭代次数。降温较快。快速模拟退火T_k T0 / (1 k)。有理论上的收敛保证。自适应降温根据当前解的接受率动态调整alpha。如果接受率高说明还没充分搜索慢点降接受率低说明已接近稳定快点降。% 自适应降温示例 (集成到主循环中) accept_rate 0; % 需要在内循环中统计接受新解的次数 % ... 内循环结束后 ... current_accept_rate accept_count / L; if current_accept_rate 0.6 alpha 0.99; % 接受率高慢点降温 elseif current_accept_rate 0.4 alpha 0.95; % 接受率低快点降温 else alpha 0.97; % 保持 end5.4 终止条件除了温度低于T_end还可以结合最大迭代次数防止无限循环。最优解连续未改进次数如果最优解连续N代如50代都没有更新可以认为已收敛。温度与能量稳定温度已很低且当前能量在多次迭代中变化极小。6. 进阶技巧与性能优化当城市数量增加到几百上千时基础版本的效率会成为瓶颈。以下是一些提升性能和效果的进阶技巧6.1 增量计算避免重复计算整个路径长度在邻域操作中路径长度变化ΔE通常只涉及被改动的那一小段路径。例如对于交换操作swap(i, j)总距离的变化只与城市i, j及其前后邻居有关无需重新计算整个路径。这能极大加速内循环。function delta_E calc_delta_E_swap(old_route, dist_matrix, i, j) % 计算交换城市i和j位置后路径长度的变化增量计算 n length(old_route); % 处理索引的循环问题前驱和后继 pre_i old_route(mod(i-2, n) 1); % i的前一个城市 suc_i old_route(mod(i, n) 1); % i的后一个城市 pre_j old_route(mod(j-2, n) 1); % j的前一个城市 suc_j old_route(mod(j, n) 1); % j的后一个城市 % 旧边 (pre_i-i), (i-suc_i), (pre_j-j), (j-suc_j) % 新边 (pre_i-j), (j-suc_i), (pre_j-i), (i-suc_j) % 注意当i和j相邻时有些边会重合需要特殊处理。这里假设了i和j不相邻的通用情况。 % 为了简化这里给出一个更通用的思路直接计算受影响的那一小段路径的旧长度和新长度作差。 % 一个更稳健但复杂的方法是先执行交换得到新路径然后计算受影响片段的变化。 % 对于初学者在性能要求不高时完整重算更安全。此处代码从略示意思想。 end6.2 多种邻域操作混合与自适应选择不要只使用一种邻域操作。混合使用交换、逆转、插入等并让算法自适应选择表现好的操作。可以记录每种操作在过去一段时间内成功改进解的次数动态调整其被选中的概率。6.3 并行化与多起点搜索模拟退火的内循环迭代是相互独立的非常适合并行计算。你可以使用Matlab的parfor循环来并行执行内循环的部分迭代。此外从多个不同的初始解开始运行多个独立的模拟退火进程最后取最优结果可以有效避免对单一初始解的依赖。6.4 与局部搜索算法结合混合策略模拟退火擅长全局探索局部搜索如2-opt, 3-opt擅长局部挖掘。一个常见的策略是在模拟退火的每个温度下或者当找到一个新解时立刻对其执行一次快速的局部搜索比如几次2-opt移动将其推到最近的局部最优然后再继续退火过程。这种“模拟退火局部搜索”的混合策略往往能取得更好的效果。7. 避坑指南那些我踩过的雷温度下降过快alpha太小这是新手最常见的问题。温度骤降算法还没来得及充分探索就迅速“冻结”结果和简单的局部搜索差不多很容易陷入局部最优。对策使用较大的alpha如0.995或者采用自适应降温策略并监控接受率。马尔可夫链长度 L 不足在每个温度下只尝试了几次就降温系统远未达到平衡。对策将L设置为与问题规模相关如100*N或者实现“直至平衡”的终止条件。邻域操作设计不当如果邻域操作产生的变化太小搜索效率低下变化太大接受率会很低。对于TSP逆转2-opt操作通常是比单纯交换更有效的邻域。随机数种子模拟退火是随机算法每次运行结果可能不同。为了结果可复现在调试阶段固定随机数种子rng(1)。但在最终评估时应多次运行取统计结果。忽略可视化与中间输出闷头跑程序出了问题不知道在哪。一定要把迭代过程的关键指标当前最优解、当前温度、接受率输出或绘图。收敛曲线能告诉你算法是否在正常工作温度曲线能帮你判断降温速度是否合适。过早放弃对于复杂问题模拟退火可能需要较长时间才能找到优质解。看到初期结果不好就否定算法或参数可能为时过早。耐心观察收敛趋势。模拟退火算法就像一位有经验的探险家它不追求每一步都向上而是懂得在适当的时候“以退为进”从而有更大的机会找到最高的山峰。通过Matlab实现TSP求解我们不仅学会了一个算法更掌握了一种解决复杂优化问题的通用思维框架。记住没有一套参数能通吃所有问题耐心调试、细心观察、结合问题特性设计邻域操作才是用好这把“瑞士军刀”的关键。