Matlab模拟退火与fmincon混合优化:全局搜索+局部精修实战 写这篇东西的起因是我前阵子在一个参数标定项目里被Matlab的fmincon折磨得不行目标函数是仿真系统一跑就是几分钟的黑盒输出带回声、有噪声还带着一堆局部极小值。初始点稍微偏一点结果就稳稳陷在某个盆地不出来。后来把模拟退火算法Simulated Annealing和优化工具箱组合起来用先用全局搜索打底再用fmincon精修效果完全不一样。这篇就是我对这套方案的完整复盘包含工具箱选型、simulannealbnd核心参数、三个能直接跑的案例以及踩坑后的调参思路。1. 为什么需要模拟退火这类全局优化算法1.1 基于梯度的优化器并不万能先聊点实际的。大部分人在Matlab里接触优化第一个想到的就是fmincon。它的原理是基于梯度信息沿着最速下降方向或者拟牛顿方向一步步爬到邻近的极小点。对凸函数、光滑函数、初始点比较靠谱的问题它又快又准确实是工程利器。但工程里的优化问题往往没这么讲武德。目标函数可能是仿真结果直接算出来的连表达式都拿不到可能存在噪声数值微分出来的梯度全是乱跳更常见的是多峰函数到处都是局部极小值。这时候fmincon本质上是“站在山腰往山下走”它只能看到脚下的坡看不到远处的谷。你换个初值它就可能爬到另一个完全不同的局部坑里。模拟退火解决的就是这个痛点。它不依赖梯度不要求函数光滑也不会看见局部极小就一头扎进去。核心思路特别朴素先到处乱逛逛得差不多了再慢慢收网。实际用下来它就是给fmincon这类局部优化器打前站的最佳搭档。1.2 模拟退火的核心机制允许“坏解”出现模拟退火这个名字来源于冶金学里的退火工艺金属加热到高温后慢慢冷却原子在高温下剧烈运动冷却过程中逐渐排列成低能量的规则晶格。算法就是把目标函数当成能量函数把一组参数当成系统状态通过模拟“加热—冷却”的过程来找能量最低的状态。关键机制叫Metropolis准则。每次迭代随机扰动当前解得到一个新解。如果新解比当前解好直接接受如果新解更差也不是直接扔掉而是以一定概率接受它P exp(-(ΔE) / T)其中ΔE是新解和旧解的能量差T是当前温度。温度越高接受差解的概率越大温度越低这个概率越小。也就是说算法在前期“体温”高的时候允许翻山越岭往坏地方走借此逃出局部陷阱后期“体温”降下来就老老实实做局部精修。很多初学者不理解为什么要接受差解。你想想登山的场景你想从A山走到更高的B山中间隔着一条山谷如果只允许往上走你永远到不了B。接受差解就是允许你先下到山谷再爬到B。这就是全局搜索和局部搜索的本质区别。1.3 全局优化工具箱的定位Matlab里有一个常见的误会把“优化工具箱”当成所有优化函数的统称。其实官方划分很清晰Optimization Toolbox主要是基于梯度的局部优化和线性、二次规划fmincon、fminunc、linprog、quadprog都在这。Global Optimization Toolbox专门做全局优化模拟退火对应的simulannealbnd、遗传算法ga、粒子群particleswarm都在这里。这两个工具箱一般是配合使用的。模拟退火负责大范围找好区域fmincon负责在好区域里精确收敛。别一门心思只盯着一个函数最实用的方案往往是“全局粗搜 局部精修”的混合策略。2. Matlab优化工具箱全景与模拟退火的选型2.1 常用求解器横向对比先给一张表把几个主流求解器的特性摆出来方便你选型。求解器搜索机制是否全局约束支持典型适用场景fmincon梯度下降/拟牛顿否局部边界、线性、非线性约束凸问题、初值靠近最优解、有精确梯度fminunc梯度下降/拟牛顿否局部无约束光滑无约束问题simulannealbnd模拟退火是只原生支持边界约束黑盒、多峰、不可导、中等维度ga遗传算法是支持约束但处理非线性约束较麻烦离散变量、混合整数、复杂约束particleswarm粒子群是主要支持边界约束连续变量、计算量可控、多峰搜索这张表最关键的信息是fmincon求的是“局部最优”simulannealbnd求的是“尽量靠近全局最优”。如果问题本身是凸的别杀鸡用牛刀直接fmincon就行。如果问题肉眼可见地多峰、黑盒、不可导就别在fmincon的初值上反复赌运气了直接上模拟退火。2.2 什么场景优先选模拟退火我自己的判断标准很直接满足下面两三条就优先考虑SA目标函数没有解析表达式只能通过仿真、实验或者外部程序返回值。目标函数存在多个极值点而且你不知道哪个是全局。函数不可导、不连续或者对参数变化很敏感。变量维度不算太高一般几十维以内太高的话SA收敛会非常慢。你愿意用几分钟甚至几小时的全局搜索时间换一个稳定的最终方案。遇到过很多项目目标函数跑一次要好几分钟这种情况下你不可能让SA像无头苍蝇一样乱撞几万次。我的做法是先用低精度的代理模型或者降采样版本跑SA确定最优参数大概在哪个区域再用高精度仿真在那个小范围里做局部优化。这个思路比单纯堆迭代次数管用得多。2.3 一个必须澄清的误区工具箱归属去网上搜“Matlab模拟退火”很多人会直接输入simulannealbnd结果报错“Undefined function”。原因不是函数名记错了而是没有安装Global Optimization Toolbox。检查工具箱是否存在可以用license(test, GADS_Toolbox)如果返回1说明有全局优化工具箱如果返回0说明只有基础的Optimization Toolbox。另外还可以在Matlab里运行ver看输出列表里有没有Global Optimization Toolbox。这个细节特别容易坑到新手。我见过有人折腾了半天代码结果发现是工具箱都没装。如果你所在的公司或者学校只有Optimization Toolbox也别灰心模拟退火算法本身不复杂自己写一个核心循环也就几十行后面第4节我会给你一个能跑的手写版本。3. simulannealbnd核心参数与实操细节3.1 目标函数怎么写给SA才不踩坑simulannealbnd对目标函数的要求比fmincon宽松很多但有几个基础规则必须遵守。第一输出必须是标量第二尽量别返回NaN或Inf第三输入x的维度要和初始点x0保持一致。一个典型的匿名函数写法fun (x) (x(1)-2)^2 (x(2)1)^2 5*sin(3*x(1)); x0 [0, 0]; lb [-10, -10]; ub [10, 10]; [x, fval] simulannealbnd(fun, x0, lb, ub);这里的x是行向量还是列向量取决于你的x0。如果你传递的是行向量SA内部会把x当行向量处理目标函数里要注意x(1)、x(2)的顺序别写反。另外一个经验是如果目标函数内部有报错风险比如根号下负数、对数里出现0最好在函数开头加保护逻辑返回一个大数值而不是让程序崩溃。SA会在搜索过程中频繁试探边界和极端位置目标函数必须足够健壮。3.2 边界约束与非线性约束处理simulannealbnd原生支持的是边界约束lb和ub这个很简单。没有下界或上界时传[]就行[x, fval] simulannealbnd(fun, x0, [], ub);但实际工程问题经常有非线性约束比如某个指标不能超过阈值、某种资源总和不能超限。SA本身不直接支持非线性约束我通常用外点惩罚函数来处理。基本思路是把约束违反量放大后加到目标函数上让算法在搜索时自动避开不可行区域。惩罚函数的标准形式是cost_pen cost_original lambda * max(0, g(x))^2其中g(x) 0是不等式约束lambda是惩罚系数。lambda太小约束会被无视lambda太大会扭曲目标函数形状让SA的收敛变得很困难。我一般从1e3开始试逐步增大到1e5量级直到最优解的约束违反量小到可以忽略。3.3 温度、冷却率、迭代次数怎么调模拟退火的效果很大程度取决于温度调度。温度高前期探索范围大温度低后期容易陷入局部。降温太快跟没退火一样降温太慢计算时间会爆炸。在simulannealbnd里这些因素分散在几个参数中参数作用经验值InitialTemperature初始温度决定前期接受差解的概率根据目标函数值量级调整常用10~1000TemperatureFcn温度下降函数temperatureexp、temperaturefast默认常用指数降温ReannealInterval多少代后重新升温再降温防止过早收敛100~1000MaxIter最大迭代次数500~10000MaxFunEvals最大目标函数调用次数通常比MaxIter大几倍TolFun目标函数变化小于该值则停止1e-6~1e-8如果你完全不知道初始温度怎么设有一个相对靠谱的校准方法先固定一个初始解做若干次随机扰动统计目标函数增量的平均值Δ_mean然后设定初始接受率P0比如0.8用公式反推T0 ≈ -Δ_mean / ln(P0)意思是刚开始有80%的概率接受那些变差的解。这样算法前期才有足够的探索能力。很多默认参数跑不出好结果问题往往就出在初始温度和目标函数值量级不匹配。3.4 saoptimset关键选项速查老版本用saoptimset设置参数新版本也可以继续用。我给一个常用的配置模板options saoptimset( ... Display, iter, ... MaxIter, 3000, ... MaxFunEvals, 30000, ... InitialTemperature, 100, ... TemperatureFcn, temperatureexp, ... ReannealInterval, 200, ... TolFun, 1e-8); [x, fval] simulannealbnd(fun, x0, lb, ub, options);如果用的是较新的Matlab版本也可以改成optimoptions方式参数名基本一致options optimoptions(simulannealbnd, ... Display, iter, ... MaxIterations, 3000, ... InitialTemperature, 100);有一个细节MaxIter和MaxFunEvals是两套完全不同的终止条件。MaxIter是迭代轮数每一轮可能调用多次目标函数MaxFunEvals是目标函数调用次数上限。如果你的目标函数跑一次特别贵更应该限制MaxFunEvals而不是MaxIter。4. 三个可复现案例从单目标到组合优化4.1 案例A多峰函数寻优和fmincon对比第一个案例用经典的Rastrigin函数做测试。这个函数的特点是全局最小值在原点但周围布满了密密麻麻的局部极小点专门用来考验全局优化算法。clear; clc; rng(1); fun (x) sum(x.^2 - 10*cos(2*pi*x) 10); lb [-5, -5]; ub [5, 5]; x0 [3, 4]; % 模拟退火 options saoptimset(Display, off, MaxIter, 5000, ... InitialTemperature, 50, TemperatureFcn, temperatureexp); [x_sa, f_sa] simulannealbnd(fun, x0, lb, ub, options); % fmincon做对比 opts_fmin optimoptions(fmincon, Display, off, Algorithm, sqp); [x_f, f_f] fmincon(fun, x0, [], [], [], [], lb, ub, [], opts_fmin); fprintf(SA: x[%.4f %.4f], fval%.6f\n, x_sa(1), x_sa(2), f_sa); fprintf(fmincon: x[%.4f %.4f], fval%.6f\n, x_f(1), x_f(2), f_f);在我本地的测试里SA基本能稳定找到接近(0,0)的位置fmincon从(3,4)出发则容易停在一个局部极小点函数值明显更大。这个案例想说明的不是fmincon不行而是它确实缺乏“跨过山谷”的能力。如果换成凸函数fmincon的表现反而会更好所以选求解器一定要看问题形状。4.2 案例B带非线性约束的工程参数优化第二个案例更贴近工程目标函数很简单但有一条非线性约束。问题定义如下min (x1-3)^2 (x2-4)^2约束条件x1^2 x2^2 16这个问题的解析最优解大概在(3, sqrt(7))附近也就是约束边界上。先用惩罚函数把约束带进目标函数跑SA再用fmincon在SA的结果上精修。clear; clc; rng(2); obj (x) (x(1)-3)^2 (x(2)-4)^2; con (x) x(1)^2 x(2)^2 - 16; lambda 1e4; penObj (x) obj(x) lambda * max(0, con(x))^2; lb [-10, -10]; ub [10, 10]; x0 [0, 0]; options saoptimset(Display, off, MaxIter, 4000, ... InitialTemperature, 100); [x_sa, f_sa] simulannealbnd(penObj, x0, lb, ub, options); % 用fmincon做约束精修 nonlcon (x) deal(con(x), []); [x_refine, f_refine] fmincon(obj, x_sa, [], [], [], [], lb, ub, nonlcon); fprintf(SA粗搜索: x[%.4f %.4f], 约束值%.4f\n, x_sa(1), x_sa(2), con(x_sa)); fprintf(fmincon精修: x[%.4f %.4f], 约束值%.4f\n, x_refine(1), x_refine(2), con(x_refine));这里要注意fmincon精修时使用的是原始目标函数不再带惩罚项而是把约束正式传给nonlcon。这样SA负责找到边界附近的可行区域fmincon负责把解精确落到边界上两者分工明确。我带过一个项目最初只靠惩罚函数跑SA结果总是差那么一点后来改成这种“SA找区域fmincon找精确解”的两段式精度瞬间就上去了。4.3 案例C手写SA求解TSP自定义邻域是关键simulannealbnd并不是为离散组合优化设计的比如旅行商问题TSP它的变量是城市访问顺序是置换编码不是连续向量。强行套simulannealbnd会很别扭。这种情况下我更推荐直接用Matlab手写一个SA核心循环代码量也不大。TSP的SA思路用随机排列表示一条路径。邻域操作采用“交换两个城市的位置”或者“反转一段路径”。新路径长度更短就接受更长则以exp(-Δ/T)概率接受。温度按指数冷却T T * alpha。核心代码框架nCity 30; rng(3); cityPos rand(nCity, 2) * 100; distMat pdist2(cityPos, cityPos); path randperm(nCity); path [path, path(1)]; % 首尾相连 cost pathCost(path, distMat); T 1000; alpha 0.995; while T 1e-3 newPath swapTwo(path); newCost pathCost(newPath, distMat); delta newCost - cost; if delta 0 || exp(-delta / T) rand path newPath; cost newCost; end T T * alpha; end fprintf(最终路径长度: %.2f\n, cost); function len pathCost(p, D) len 0; for i 1:length(p)-1 len len D(p(i), p(i1)); end end function p2 swapTwo(p) p2 p; idx randperm(length(p)-2, 2) 1; % 不交换首尾 p2(idx) p2(fliplr(idx)); end这个案例的价值在于当你面对的不是标准连续优化问题而是特殊编码的组合优化问题时理解算法原理远比赛会背工具箱接口更重要。SA的核心就三件事解的表达、邻域扰动、温度调度。5. 常见问题与调优技巧实录5.1 结果重复性差怎么处理很多用户第一次跑SA会惊讶地发现每次结果都不一样。这是正常现象因为SA本质是随机算法。但如果每次结果天差地别说明算法没有充分收敛。我处理重复性问题有一个标准流程先固定随机种子rng(1)保证每次实验可复现然后增加MaxIter和MaxFunEvals最后如果计算资源允许多次运行SA取目标函数值最小的解作为最终结果。这个“多跑几次取最优”的策略虽然简单但在实际工程里非常实用。5.2 SA收敛太慢除了等待还能做什么如果跑了几千次迭代还在原地转大概率是温度策略不适合当前问题。常用的加速手段降低InitialTemperature减少前期无效的乱逛。把TemperatureFcn换成temperaturefast降温更激进。增大ReannealInterval减少重新升温的次数让算法更专心收敛。检查目标函数是否有大量平坦区域如果有SA在平坦区里会瞎转这种问题更适合加一点局部梯度信息。我踩过的一个坑是为了追求全局最优把初始温度设得特别高结果前一半计算时间全在随机乱跳真正的局部精修反而不够。5.3 惩罚函数把结果带偏了惩罚函数几乎是处理非线性约束的标配但它有个副作用惩罚系数太大目标函数表面会变得特别陡峭SA在边界附近很难稳定惩罚系数太小约束会被无视。我的建议是不要一次性定死惩罚系数而是跑一轮SA检查最终解的约束违反量。如果违反量大就把lambda乘以10再跑如果结果已经在可行域内但精度不够再用fmincon精修。实际项目中我经常用SA惩罚函数找到边界区域然后靠fmincon去满足严格的约束。5.4 SAfmincon混合优化模板综合前面的经验我目前最常用的优化流程是这样的先用模拟退火做全局搜索得到一个候补解。检查候补解的目标函数值和约束违反量。以候补解为初值调用fmincon做局部精修。如果fmincon结果不符合预期回到第1步换一个随机种子再跑一次。这段逻辑可以封装成一个函数核心代码如下function [x_best, f_best] hybrid_optimize(objFun, x0, lb, ub, nonlcon) saOptions saoptimset(Display, off, MaxIter, 3000); [x_sa, ~] simulannealbnd(objFun, x0, lb, ub, saOptions); fminOptions optimoptions(fmincon, Display, off); if nargin 5 nonlcon []; end [x_best, f_best] fmincon(objFun, x_sa, [], [], [], [], lb, ub, nonlcon, fminOptions); end这种混合策略几乎成了我做中型优化项目的默认起手式。它既保留了SA的全局搜索能力又利用了fmincon在局部的高精度收敛两个工具箱才算真正打出了配合。5.5 问题速查表现象常见原因解决办法每次运行结果都不同随机算法未充分收敛固定rng、增加迭代次数、多次运行取最优收敛到明显很差的解初始温度过低或降温过快提高InitialTemperature更换temperatureexp搜索过程极其缓慢初始温度过高、退火区间过长降低初始温度、增大ReannealInterval约束始终不满足惩罚系数太小增大lambda或改用fmincon处理约束目标函数返回NaN搜索到无效区域在目标函数内加保护逻辑返回大数结果停留在初始点附近邻域扰动幅度太小检查x0向量范围调整初始温度6. 我的实操心得与后续扩展6.1 踩坑后总结的习惯我现在做任何一个优化项目都会先花时间把目标函数的行为摸一遍。能画图就画图不能画图就随机采样若干组参数看一眼目标函数的量级和粗糙程度。这个习惯帮我省了无数调参时间。第二个习惯是永远不要把SA的结果当最终答案。SA擅长的是“找到好区域”而不是“给出极高精度的解”。除非你的问题对精度要求不高否则建议在SA结果后面再接一个局部优化器。哪怕是fminunc这种无约束的局部求解器也够用。第三个习惯是记录每一次SA运行的轨迹。Matlab里可以通过自定义输出函数或绘图函数把每次迭代的当前最优值记录下来后期分析收敛过程非常有用。但自定义绘图函数在新版本里有时兼容性一般我直接保存变量画图自己来。6.2 还可以怎么扩展模拟退火的灵活性在于温度调度和邻域操作都可以自定义。我在一个项目里给温度函数加入“阶段性升温”策略每跑固定轮数后如果检测到连续多轮最优值没有变化就自动升温重新搜索。效果比固定降温曲线好不少。如果你对Matlab的全局优化工具箱已经很熟还可以进一步把SA与其他方法组合先用SA做变量筛选再用线性回归或响应面模型精修或者在SA的目标函数中嵌入仿真模型实现真正的黑盒优化。只要理解了“全局粗搜局部精修”这个核心思想工具箱就只是工具真正解决问题的是你的调度策略。模拟退火不是万能的但在我处理过的多峰、黑盒、不可导问题里它几乎是最稳的兜底方案。配合Matlab优化工具箱里的局部求解器这一套组合拳值得成为你优化武器库里的常备选项。