简介这份资源是一份基于Matlab的模拟退火算法入门源码包面向正在学习启发式优化算法、需要快速上手实现模拟退火求解组合优化问题的科研人员与工程师。模拟退火算法源于固体退火过程适合处理旅行商、装载、图着色等具有多个局部最优解的NP难问题包内通过Matlab脚本展示了状态生成、能量计算、Boltzmann接受准则与温度调度等核心环节可帮助读者理解算法流程并在此基础上改造自己的目标函数。资源共2个文件均为m脚本包含一个可运行的主程序和一个核心算法函数压缩包整体仅817B轻量精简便于对照学习与二次开发。目前已有555人学习下载适合具备一定Matlab基础、希望从代码层面快速掌握模拟退火实现思路的初学者。通过阅读这两个脚本读者可以直观看到工具箱调用方式、参数初始化逻辑以及结果输出的组织方法为后续将模拟退火应用到实际优化问题中提供可直接参考的模板。 遇到一个没有解析梯度的非线性优化问题大多数人第一反应是跑fmincon跑几次发现结果每次不同才想起模拟退火算法。Matlab 里的模拟退火工具通常有两种形态全局优化工具箱里封装好的simulannealbnd以及网上流传、面向教学和二次开发的模拟退火算法工具箱源码。后者解压后往往是一组.m文件目标函数、邻域扰动、温度调度都被拆开看起来能直接套用但真正的门槛在初温、调度方式和邻域扰动三者的配合。下面按“原理—最小实现—源码拆解—用接受率调参”的顺序讲适合要快速把手头问题接进模拟退火的人也适合想读懂这类源码再改的 Matlab 用户。1. 模拟退火算法Matlab工具箱不是有simulannealbnd就够了2. 模拟退火的收敛机制与Matlab工具箱选型先看Metropolis再看参数2.1 Metropolis准则退火算法能逃出局部最优的唯一开关模拟退火算法本身不复杂核心就一条以概率接受比当前解更差的新解。设当前目标值为f(x)新解目标值为f(x_new)能量差ΔE f(x_new) - f(x)。当ΔE 0时新解必被接受当ΔE 0时接受概率为exp(-ΔE/T)。这里T是当前温度也是整个工具箱里最关键的参数。温度高时exp(-ΔE/T)接近 1算法基本在随机游走温度低时这个概率趋近于 0算法退化成爬山法。模拟退火能跳出局部最优就是靠温度还没降下来时接受劣解。如果初始温度给得太低算法从头到尾都在“爬山”那和fmincon的差别只剩下收敛慢如果降温太快则等价于快速淬火随机性还没发挥作用就已经冻住。理解这一条后面看源码才知道哪些参数能碰、哪些不能动。2.2 Matlab官方优化工具箱与源码自建怎么选Matlab 里跑模拟退火算法最快的方式是全局优化工具箱的simulannealbndoptions saoptimset(InitialTemperature, 100, ... TemperatureFcn, temperatureexp, ... MaxIterations, 3000); [x, fval] simulannealbnd(rastrigin, x0, lb, ub, options);这条命令能跑通但很多时候只是“能跑”。saoptimset暴露出来的调度函数有限邻域扰动也是内置的对具体问题做定制时反而受限制。压缩包里常见的教学型源码走的是另一条路主循环自己写接受判断自己写温度调度自己写可定制程度高很多。对比项官方 simulannealbnd自建/教学源码上手成本一条命令默认参数可跑25 个文件需要自己组织邻域扰动内置随机方向扰动可针对问题写 2-opt、反转、变异温度调度几种温度函数可选完全可控可做回火约束表达边界约束友好罚函数或自定义接受机制排错成本封装深报错绕每步可打印容易定位我一般拿到这类源码包第一件事不是直接调参数而是把它拆成“目标函数、邻域扰动、接受判断、温度调度”四块再逐个替换成自己的实现。官方工具箱适合验证一个想法源码才适合改成一个真正合身的最优化工具。2.3 初温、末温、降温系数对收敛行为的影响三个参数决定 90% 的收敛行为。初温T0决定算法最初能接受多差的解末温T_end决定算法什么时候停止降温系数alpha决定温度下降曲线。最常见的调度是指数降温也就是每个外循环把温度乘上一个固定比例T_new max(T * alpha, T_end);alpha取 0.9 时温度从 100 降到 1e-8 大约需要 210 轮取 0.99 时需要 2200 轮左右。这不是“多算一会儿”的问题而是直接影响邻域搜索步长。若源码里把扰动幅度写成T * randn(size(x))那么温度降得越快搜索半径缩得越急很多区域还没来得及探索就已经被冻结了。初温的设定没有解析公式工程上常用做法是先让解做几十次随机游走统计目标值变化的平均值取它的 520 倍作为T0。我在第 5 章会给具体代码。末温相对简单只要保证温度低到exp(-ΔE/T)对最大可接受的劣解都不敏感即可通常取 1e-6 到 1e-8 已经足够。3. 用Matlab从零实现最小模拟退火算法求解器可运行的源码3.1 主函数框架与边界处理先给一个可直接保存成.m文件的最小实现。这个版本适合连续优化问题变量用列向量表示边界采用反弹策略。完整源码如下function [bestx, bestf, trace] sa_minimize(fun, x0, lb, ub, opts) % 模拟退火最小化求解器 % 输入: % fun : 目标函数句柄输入列向量 x返回标量 f % x0 : 初始解列向量 % lb,ub: 变量下界与上界 % opts : 结构体 % opts.T0 初始温度 % opts.T_end 终止温度 % opts.alpha 降温系数 % opts.max_inner 每个温度下抽样次数 % opts.max_outer 最大降温轮数 x x0(:); fx fun(x); bestx x; bestf fx; T opts.T0; trace zeros(opts.max_outer, 3); for k 1:opts.max_outer n_accept 0; for j 1:opts.max_inner % 邻域扰动步长与当前温度挂钩 x_new x T * randn(size(x)); % 越界反弹 low x_new lb; high x_new ub; x_new(low) 2 * lb(low) - x_new(low); x_new(high) 2 * ub(high) - x_new(high); f_new fun(x_new); delta f_new - fx; % Metropolis 接受判断 if delta 0 || rand() exp(-delta / T) x x_new; fx f_new; n_accept n_accept 1; if fx bestf bestf fx; bestx x; end end end % 记录温度、历史最优值、接受率 trace(k, :) [T, bestf, n_accept / opts.max_inner]; % 指数降温 T max(T * opts.alpha, opts.T_end); if T opts.T_end break; end end % 截断未用到的 trace 行 trace trace(1:k, :); end这个实现刻意保持“单文件”方便读懂流程。x_new x T * randn(size(x))是连续问题最常见的扰动方式标准差直接取当前温度高温时大步探索低温时小步精调。这种耦合写法简单但第 4 章会说明它的副作用。越界反弹用2 * bound - value而不是裁切是为了尽量保留扰动脉冲的方向信息避免大量点堆积在边界上。3.2 接受判断与降温循环接受判断是这段源码里最值得细看的地方。标准写法rand() exp(-delta / T)在delta很大或T很小时exp参数接近-inf计算结果大概率下溢为 0此时接受判断没问题但当delta/T在 30 到 700 之间时浮点计算会损失精度。更稳定的写法是把不等式两边取对数变成delta / T -log(rand())这个留到 4.4 展开。外层for循环控制温度内层for控制每个温度下的马尔可夫链长度。max_inner取 50200 是常见区间。取值太小时每个温度下只试了几个点温度就降了算法偏向随机爬山取值太大收敛轮数增加但对问题维度较高时更稳。工程上我会先用 100 起步观察接受率再调整。3.3 跑一版Rastrigin测试并输出收敛轨迹用 Rastrigin 函数做冒烟测试这个函数在多维空间里布满了局部极小点是检验模拟退火算法的标准题目。调用代码如下rastrigin (x) sum(x(:).^2 - 10*cos(2*pi*x(:)) 10); opts struct(T0, 100, T_end, 1e-8, alpha, 0.9, ... max_inner, 80, max_outer, 300); [bx, bf, tr] sa_minimize(rastrigin, [2.5; -3.1], ... [-5.12; -5.12], [5.12; 5.12], opts); fprintf(最优解 x[%.4f, %.4f] 目标值 %.6f\n, bx(1), bx(2), bf);alpha 0.9时温度从 100 降到 1e-8 约需 210 轮max_outer 300留有足够余量。Rastrigin 二维问题的全局最优点是[0, 0]目标值 0。如果运行结果反复停在其他位置优先怀疑初温过低或max_inner太小而不是代码有 bug。4. 模拟退火工具箱源码结构拆解目标接口、邻域扰动与调度策略4.1 工具箱应该拆成哪些文件一个能用于多类问题的模拟退火算法工具箱很少把所有逻辑塞进一个文件。拿到网上流传的“模拟退火算法工具箱源码”压缩包时我一般会先把它重构成下面这种结构再逐文件对照学习sa_toolbox/ ├── sa_minimize.m % 主入口控制外层降温与内层抽样 ├── sa_init.m % 解析 opts填充默认参数 ├── sa_neighbor.m % 邻域扰动函数 ├── sa_accept.m % Metropolis 接受判断 ├── sa_cooling.m % 温度调度函数 ├── sa_restart.m % 回火重启策略 └── demo_rastrigin.m % 示例脚本sa_minimize.m调用其他模块而不是自己包办所有逻辑目的是让使用者能只替换sa_neighbor.m就把连续优化改成组合优化。源码可读性很大程度上取决于目标函数接口是否统一。常见约定是目标函数句柄接收一个列向量返回一个标量。工具箱内部不管具体问题只负责按当前解产生新解、计算目标值、按接受率更新温度。4.2 邻域扰动连续问题与组合问题的两种实现连续问题的扰动在第 3 章已经出现。组合问题则完全不同以旅行商问题为例最常用的 2-opt 反转的伪代码风格实现如下function new_s sa_neighbor_tsp(s) n numel(s); i randi([1, n - 1]); len randi([1, n - i]); j i len; new_s s; new_s(i:j) new_s(j:-1:i); % 反转区间保持路径闭环 end这段代码里的i和j是随机选取的两个断点new_s(i:j) new_s(j:-1:i)把区间内访问顺序反转。这种邻域扰动对 TSP 天然有效因为它改变了路径交叉而不破坏路径的闭环结构。对比第 3 章连续问题的x T * randn(size(x))你会发现源码质量的差别主要在这里扰动方式是否匹配问题结构直接决定收敛速度和最终质量。4.3 温度调度的三种升级写法指数降温是默认方案但不是唯一方案。教学源码包通常只演示指数调度实际做二次开发时我会在sa_cooling.m里加入线性降温和周期回火两种模式function T sa_cooling(k, opts) % 模式1: 指数降温 T opts.T0 * opts.alpha^k; % 模式2: 周期回火每隔 reanneal_interval 轮升温一次 if opts.reanneal_interval 0 mod(k, opts.reanneal_interval) 0 T T opts.T_rise; end end回火的操作是每reanneal_interval轮把温度加上T_rise让算法重新获得跳出局部最优的能力。它和“重启”不一样回火保留了当前最优解的信息只是短暂允许更大范围的扰动。T_rise一般取当前温度的 0.10.5 倍取太大等于重新初始化取太小没有效果。这个技巧在官方simulannealbnd里对应内置的ReannealInterval自建源码时更好控制。4.4 源码里最容易改错的三个细节第一温度与步长强耦合。x_new x T * randn(size(x))虽然直观但低温期步长变得极小算法在最后阶段几乎只在当前点附近打转。改进做法是把扰动步长sigma独立出来或让sigma按接受率自动调整接受率太高就增大步长太低就缩小。第二exp(-delta/T)的数值稳定性。上一章提到的对数重写版本如下function ok sa_accept(delta, T) if delta 0 ok true; else % 对 rand exp(-delta/T) 两边取对数避免 exp 下溢 ok (delta / T) -log(rand()); end end逻辑上它与rand() exp(-delta/T)完全等价但不会因为exp下溢导致接受概率被人为截断成 0。第三历史最优解与当前解混用。模拟退火允许接受劣解所以算法结束时的x大概率不是历史最优。工具箱源码里必须在接受新解的同时记录单独的bestx和bestf最后返回的是bestx而不是当前解。这个错误在初版实现里非常常见诊断方法是观察trace的第二列是否单调不增如果出现回退就是记录逻辑混了。5. 参数是否有问题看接受率比看收敛曲线更快5.1 用单次运行诊断初始温度sa_minimize的导出trace第三列是每轮接受率。接受率是判断初温是否合理的第一个信号第一轮接受率如果低于 0.5说明T0太小算法初始阶段就已经开始“爬山”如果高于 0.99说明温度过高前 1/3 的轮次基本浪费在纯随机游走上。正常区间是第一轮 0.70.95末轮低于 0.01。不想盲调T0的话先用随机游走估算目标值变化量级比反复跑完整退火快得多dE zeros(200, 1); for i 1:200 x_try x0 randn(size(x0)); dE(i) abs(fun(x_try) - fun(x0)); end T0 10 * mean(dE);这段代码不跑退火只统计邻域扰动的目标值平均变化再把T0设为该均值的 10 倍。对新问题尤其有效可以一次性把初温设定到一个可用量级后续再按接受率修正。5.2 回火重试与收敛曲线的配合连续多轮接受率为 0 且历史最优值不再变化通常是陷在一个较深的局部极小。常见的处理是给sa_restart加回火逻辑检测到连续 30 轮没有任何接受时把当前温度恢复到T0 * 0.3重新加热升温而不是整体重启。和整体重启相比它保住了已找到的bestx同时给当前解一个重新探索的机会。调试时同时看温度和最优值的双纵轴图一段代码就能判断退火是否进行得合理figure; yyaxis left; semilogy(tr(:,1)); ylabel(Temperature); yyaxis right; plot(tr(:,2)); ylabel(Best f); xlabel(outer iteration);左轴温度指数下降右轴最优值应呈阶梯状下降而不是直线。如果bestf曲线在最后 20% 轮次还在大幅变动说明max_outer不够或降温系数偏大如果早就平了但接受率还在 0.3 以上说明温度没降到位。把接受率输出放在每轮结束连续三轮低于 0.05 时把温度恢复到当前值的四倍再跑二十轮比重新全局初始化省时间。本文还有配套的精品资源点击获取