模拟退火算法MATLAB实战:从排产优化到参数标定
简介这份资源面向希望系统掌握模拟退火算法并落地到MATLAB建模的算法学习者与工程人员围绕全局优化这一核心问题展开。内容从算法原理切入讲解初始温度、冷却因子与终止温度三大参数的作用并借助Metropolis准则说明如何以一定概率接受恶化解从而跳出局部最优。压缩包共10个文件以9个m脚本和1个mp4视频为主脚本承载算法实现与案例代码视频用于演示讲解整体约365.64MB。资源涵盖旅行商问题、整数规划等典型建模案例配套迭代解变化与目标函数随温度变化的可视化结果并给出参数调整思路帮助读者理解算法动态行为。已有137人学习适合想提升MATLAB复杂优化问题求解能力、积累实战经验的读者参考。1. 从一次车间排产翻车说起模拟退火算法到底能解决什么去年帮一家做线材裁切的厂子做排产优化需求很朴素给定一批长度不同的原料和一批订单怎么切最省料。我一开始用贪心跑出来废料率 11%车间主任看了一眼说“你这还不如老师傅手排”。换成模拟退火算法重写同一批数据废料率压到 4.3%而且约束一改比如某规格必须优先满足改个目标函数就能重跑。这件事让我彻底改变了对模拟退火算法Simulated AnnealingSA的看法——它不是教科书里那个“爬山会卡局部最优”的玩具而是在 MATLAB 里几十行就能落地、对目标函数几乎无要求的通用优化器。这篇笔记面向三类人正在做数学建模、需要快速出结果的在校同学手上有排产/参数标定/路径规划这类组合优化问题、但不想上商业求解器的工程师以及被粒子群、遗传算法调参折磨过、想找一个“参数少、逻辑透明”的替代方案的人。我会把 SA 的建模思路、MATLAB 实现骨架、参数怎么设、以及我踩过的坑一次讲清你照着改目标函数就能用在自己的问题上。2. 模拟退火算法的建模逻辑为什么它比爬山法更值得投入2.1 从金属退火到 Metropolis 准则模拟退火的核心思想借自冶金高温时原子活跃、能跳出局部低能态缓慢降温后系统趋于全局能量最低。映射到优化问题上目标函数值就是能量解就是原子构型温度 T 是控制探索范围的关键参数。算法每一步做三件事在当前解附近随机扰动生成新解计算目标函数差值 ΔE按 Metropolis 准则决定是否接受。这个准则就是 SA 区别于爬山法的全部秘密ΔE 0新解更优无条件接受。ΔE 0新解更差以概率 exp(-ΔE / T) 接受。温度高时exp(-ΔE/T) 接近 1差解也大概率被接受算法在大范围游走温度低时差解几乎被拒绝算法退化为局部搜索。正是这个“以一定概率接受坏解”的机制让 SA 有机会翻出局部最优的坑。很多人第一次实现 SA 时把接受概率写错比如漏了负号或用了 ΔE/T 直接比较结果要么不收敛要么退化成随机游走这是最常见的翻车点。2.2 为什么在 MATLAB 里选 SA 而不是粒子群数学建模和工程现场我一般这样选连续、低维、目标函数光滑的问题粒子群PSO收敛快但一旦问题变成离散组合优化排产、TSP、装箱、线材裁切PSO 的“速度”概念很难定义编码和修复约束会写到你怀疑人生。SA 对解的形式没有要求只要你能定义“邻域扰动”和“目标函数”两件事它就能跑。维度模拟退火 SA粒子群 PSO解的形式任意离散/连续/混合主要连续需要调的参数初始温度、降温系数、终止温度、链长惯性权重、学习因子、种群规模约束处理罚函数或邻域内修复罚函数为主跳出局部最优概率接受机制天然支持依赖种群多样性MATLAB 实现难度低核心 30 行中需处理速度和边界我一般会先上 SA 把问题跑通、拿到一个可信的基线再考虑要不要换更复杂的算法。先能跑再谈快这是工程上的顺序。2.3 一个最小可跑的 SA 骨架下面这段代码是我每次开新问题时的模板目标函数换成你自己的即可。这里用经典的 Rastrigin 函数做演示它有大量局部极小正好检验 SA 的跳出能力。% sa_min_template.m —— 模拟退火最小骨架 clear; clc; rng(42); % 固定随机种子保证结果可复现 % 1. 目标函数换成你的问题 objfun (x) 20 x(1)^2 x(2)^2 - 10*(cos(2*pi*x(1)) cos(2*pi*x(2))); % 2. 参数设置 T0 100; % 初始温度要足够高让初始接受率在 0.7~0.9 Tend 1e-3; % 终止温度 alpha 0.95; % 降温系数越接近 1 降温越慢、越稳 L 200; % 每个温度下的迭代链长Markov 链长度 lb [-5.12 -5.12]; ub [ 5.12 5.12]; % 3. 初始解 x lb (ub - lb) .* rand(1,2); best_x x; best_f objfun(x); T T0; while T Tend for k 1:L % 邻域扰动高斯扰动步长随温度缩放 x_new x randn(1,2) .* (0.1 * sqrt(T)); x_new min(max(x_new, lb), ub); % 边界裁剪 f_new objfun(x_new); dE f_new - objfun(x); if dE 0 || rand exp(-dE / T) % Metropolis 准则 x x_new; if f_new best_f best_f f_new; best_x x_new; end end end T T * alpha; % 降温 end fprintf(最优解: (%.4f, %.4f), 最优值: %.6f\n, best_x(1), best_x(2), best_f);逻辑说明外层while控制降温内层for是每个温度下的马尔可夫链链长 L 决定该温度下探索是否充分。dE 0 || rand exp(-dE/T)这一行就是 Metropolis 准则注意exp(-dE/T)里 dE 为正、T 为正结果落在 (0,1)和rand比较才成立。参数说明T0太小会导致一开始就困在局部alpha取 0.9~0.99越小降温越快但越容易早熟L一般取问题维度的 10~100 倍。邻域步长0.1*sqrt(T)是我常用的自适应写法温度高时步子大、温度低时步子小比固定步长稳得多。3. 把 SA 用到真实建模案例线材裁切与参数标定3.1 线材裁切问题的 SA 建模回到开头那个线材裁切案例。问题描述原料长度 L6000mm订单需要若干种长度和数量求切割方案使废料最少。这是一个一维下料问题属于典型的整数组合优化。建模的关键是解的编码和邻域扰动。我用的编码是每个订单分配一个“切割模式编号”模式表预先枚举所有可行的切割组合。扰动就是随机改变某个订单的模式编号。% cutting_stock_sa.m —— 线材裁切模拟退火 rng(1); L 6000; % 原料长度 demand [1200 800 1500 600 900]; % 各规格需求长度 qty [10 20 8 15 12]; % 各规格需求数量 % 1. 枚举可行切割模式每种模式是一组各规格的数量 patterns []; for a 0:floor(L/demand(1)) for b 0:floor((L-a*demand(1))/demand(2)) for c 0:floor((L-a*demand(1)-b*demand(2))/demand(3)) for d 0:floor((L-a*demand(1)-b*demand(2)-c*demand(3))/demand(4)) e floor((L-a*demand(1)-b*demand(2)-c*demand(3)-d*demand(4))/demand(5)); if abcde 0 patterns [patterns; a b c d e]; end end end end end np size(patterns,1); % 2. 目标函数给定每种模式的使用次数算总废料 % 这里用“模式使用次数”作为解需满足需求约束罚函数处理 objfun (cnt) sum(cnt) * L - sum(min(sum(cnt .* patterns,1), qty) .* demand) ... 1e4 * sum(max(0, qty - sum(cnt .* patterns,1))); % 3. SA 主循环 T 500; Tend 1e-2; alpha 0.97; Lk 100; cnt randi([0 5], 1, np); % 初始解每种模式随机用几次 best_cnt cnt; best_f objfun(cnt); while T Tend for k 1:Lk cnt_new cnt; idx randi(np); cnt_new(idx) max(0, cnt_new(idx) randi([-2 2])); % 扰动 dE objfun(cnt_new) - objfun(cnt); if dE 0 || rand exp(-dE/T) cnt cnt_new; if objfun(cnt) best_f best_f objfun(cnt); best_cnt cnt; end end end T T * alpha; end fprintf(最少废料: %.1f mm\n, best_f);逻辑说明patterns枚举了所有不超原料长度的切割组合解向量cnt表示每种模式用几次。目标函数里第一项是总用料第二项是实际满足的需求长度第三项是罚函数——需求没满足就加一个大惩罚。罚系数 1e4 要远大于废料量级否则算法会为了省料而牺牲需求满足这是下料问题的经典坑。参数说明T500对应废料量级毫米级初始温度要让接受率合理alpha0.97降温偏慢因为模式空间大扰动randi([-2 2])允许一次改多次加快收敛。3.2 连续参数标定SA 拟合非线性模型SA 不只做离散问题。我做过一个储能衰减建模的活需要拟合一个双指数衰减模型的三个参数用最小二乘做目标函数SA 直接搜参数空间。相比lsqcurvefitSA 的好处是不依赖初值对多峰的目标函数也能找到全局附近。% param_fit_sa.m —— 用 SA 拟合双指数衰减模型 t 0:0.5:50; y_true 80*exp(-0.05*t) 20*exp(-0.5*t) randn(size(t))*0.3; % 含噪观测 objfun (p) sum((p(1)*exp(-p(2)*t) p(3)*exp(-p(4)*t) - y_true).^2); lb [50 0.001 5 0.1]; ub [100 0.2 40 1.0]; T 1e4; Tend 1e-1; alpha 0.96; Lk 300; p lb (ub-lb).*rand(1,4); best_p p; best_f objfun(p); while T Tend for k 1:Lk p_new p randn(1,4).*[2 0.005 2 0.02]*sqrt(T/1e4); p_new min(max(p_new, lb), ub); dE objfun(p_new) - objfun(p); if dE 0 || rand exp(-dE/T) p p_new; if objfun(p) best_f best_f objfun(p); best_p p; end end end T T * alpha; end fprintf(拟合参数: A1%.2f k1%.4f A2%.2f k2%.4f\n, best_p);逻辑说明目标函数是残差平方和四个参数各有物理量纲所以扰动步长按参数尺度分别缩放[2 0.005 2 0.02]不能统一用 randn。这是连续问题 SA 的关键细节——各维尺度差异大时必须归一化或分维设步长否则小量纲的参数几乎不动。参数说明T1e4对应残差平方和的量级Lk300比离散问题大因为连续空间需要更多采样边界lb/ub既是物理约束也是防止参数跑飞的手段。4. 避坑与排查SA 落地时最容易翻车的 5 个地方4.1 现象算法跑完结果还不如贪心原因初始温度太低或降温太快算法一开始就退化成局部搜索根本没探索。解决先做一次“接受率标定”——在初始温度下跑 100 步统计接受比例目标是把初始接受率调到 0.7~0.9。如果低于 0.5把 T0 乘 2 再试。我一般写个小脚本自动扫 T0比拍脑袋靠谱。4.2 现象结果每次跑都不一样没法写进报告原因SA 是随机算法没固定种子。解决rng(42)固定随机流报告里注明种子。但要注意固定种子不等于结果稳定如果参数本身不合理换个种子结果会差很多。判断标准是换 5 个种子跑最优值波动在 1% 以内才算参数调好了。4.3 现象约束总被违反罚函数压不住原因罚系数太小或者邻域扰动直接产生不可行解且没修复。解决罚系数取目标函数典型值的 100~1000 倍更好的做法是在邻域扰动后加一步“修复”比如下料问题里把超需求的模式次数截断。罚函数是最后手段能在邻域内修复就别用罚。4.4 现象跑到后期几乎不动但离最优还差很远原因降温系数 alpha 太小或者链长 L 太短某个温度下还没探索充分就降温了。解决把 alpha 调到 0.98~0.99L 调到维度的 50~100 倍。代价是运行时间变长但 SA 本来就慢工出细活。如果时间紧宁可减少外层降温次数、增大 L也不要反过来。4.5 现象内存爆了或跑了一晚上没结束原因目标函数里反复枚举或重复计算或者终止条件写成了死循环。解决把能预计算的量比如模式表、距离矩阵在循环外算好终止条件除了温度再加一个“最大迭代次数”兜底。我一般会加iter max_iter的硬中断防止参数写错导致无限跑。5. 进阶技巧自适应降温与并行重启5.1 自适应降温让温度自己决定快慢固定 alpha 的问题是简单问题降温太慢浪费时间难问题降温太快早熟。我后来改用基于接受率的自适应降温每个温度链结束后统计接受率如果高于 0.6 就多降一点低于 0.2 就少降一点。% 自适应降温片段 accept_rate accepted / Lk; if accept_rate 0.6 T T * 0.90; % 探索充分加速降温 elseif accept_rate 0.2 T T * 0.99; % 快冻住了放慢 else T T * 0.95; end这个改动让同一套参数能适配不同难度的问题省去了反复手调 alpha 的功夫。代价是要多维护一个accepted计数器逻辑上多几行。5.2 并行重启用多核换稳定性SA 单次运行有随机性工程上我常用并行重启同时跑 8 个不同种子的 SA取最优。MATLAB 里用parfor一行就能并行化。% 并行重启 n_runs 8; results zeros(n_runs, 1); parfor i 1:n_runs rng(i); % 每个 worker 不同种子 results(i) run_sa_once(); % 封装好的单次 SA end [best_val, idx] min(results); fprintf(8 次重启最优: %.4f (第 %d 次)\n, best_val, idx);注意parfor里不能用全局随机流必须每个 worker 单独rng(i)否则会报错或结果不可复现。并行重启的收益是稳定性单次可能翻车8 次里总有一次靠谱。代价是 8 倍算力但对中小规模问题完全可接受。5.3 怎么判断 SA 结果可信我一般做两件事验证一是换种子跑 5 次看最优值波动二是和贪心/枚举小规模实例对比确认 SA 在小规模上能拿到已知最优。如果小规模都对不上说明目标函数或邻域设计有问题别急着上大规模。这套验证习惯帮我省过好几次“以为算法牛、其实是目标函数写错”的尴尬。最后说个我自己的习惯每次用 SA我都会先把目标函数单独拿出来用几个手算的已知解验证一遍确认它算出来的值符合预期再丢进 SA 循环。这个“后悔药”花不了五分钟但能挡掉八成以上的翻车。希望帮到你。本文还有配套的精品资源点击获取