简介增广拉格朗日算法是求解带约束非线性优化问题的经典迭代方法它结合拉格朗日乘子法与惩罚函数法的优点通过引入惩罚项和增广乘子逐步逼近最优解被广泛应用于图像处理、机器学习、信号处理等领域。这份MATLAB源码针对这类需求编写特别适合需要处理非凸、非光滑或大规模优化问题的研究者与工程师可作为直接运行的参考实现。源码基于近端交替线性化最小化PALM框架在每一轮迭代中先求解线性化局部最小问题以更新决策变量再根据约束偏差更新拉格朗日乘子并动态增大惩罚系数以强化约束满足同时包含初始化、收敛判断等完整逻辑结构清晰、便于二次开发与教学演示。资源包非常精简仅包含1个m脚本文件总大小只有2KB可以无缝集成到已有项目中。目前已有3092人学习下载足以说明其实用性是学习增广拉格朗日算法和PALM方法的优质样例。1. 增广拉格朗日算法是约束优化的稳定主力MATLAB 里比罚函数更值得先掌握做约束优化时很多人第一反应是罚函数把等式约束乘一个大系数塞进目标函数权重从 10 调到 1e6解倒是接近可行域了海森矩阵却越来越病态迭代速度肉眼可见地变慢。增广拉格朗日算法Augmented Lagrangian Method也叫增广拉格朗日乘子法正好补上这一环它保留罚项的稳定结构但通过乘子迭代让罚系数保持有限就能精确收敛到可行点。对中等规模非线性受约束问题它实现简单、对初值不敏感思想上也与 MATLAB 优化工具箱内置算法同源。适合想自己控制迭代过程、做算法对比或写科研与课程代码的工程师和学生。2. 乘子法原理先说透拉格朗日函数、罚函数病态与 L_ρ 的构造增广拉格朗日不是天外来客它是在拉格朗日乘子法和罚函数法之间架起的一座桥。先把两侧概念对齐再理解中间多了什么后面的 MATLAB 代码才有依据。2.1 从拉格朗日乘子法与 KKT 条件说起等式约束的精确解在哪考虑光滑问题 min f(x)约束 c(x)0。定义拉格朗日函数 L(x,λ)f(x)λᵀc(x)λ 是维数等于约束个数的乘子向量。在最优点处KKT 条件给出两个方程∇f(x)∇c(x)ᵀλ0 与 c(x)0。前者是驻点条件后者是约束可行性联立解出的是局部最优的候选点还需要再配合二阶充分条件才能确认。直接用牛顿法解这个方程组小规模可行工程里却常受挫于两个现实一是 f、c 的解析梯度和海森往往拿不到仿真模型通常只能给出函数值二是牛顿迭代收敛半径窄初始点不佳就会发散。所以实用策略是把问题拆成两层内层做无约束极小化外层更新乘子。这正是增广拉格朗日算法的主框架也是后面 MATLAB 代码的骨架。2.2 纯罚函数的病态为什么 ρ 不能无限大罚函数法把问题改写成 min f(x)(ρ/2)||c(x)||²。理论保证 ρ→∞ 时罚函数的最优点趋于原问题最优点但工程上ρ→∞是个危险信号。ρ 变大意味着目标函数在约束面法向上越来越陡Hessian 条件数近似正比于 ρ梯度法步长被压小拟牛顿法也需要更多迭代才能在内层站稳。以 min x² s.t. x1 为例罚函数子问题解析解是 x*ρ/(ρ2)约束残差为 2/(ρ2)。要把残差压到 1e-4ρ 至少要取到上万此时子问题已经明显病态。增广拉格朗日固定 ρ2靠乘子外循环约 25 轮就能把同一个残差压到 2e-8ρ 自始至终不用变大。下表是同一问题下两类方法的数值对照第 4 章的 MATLAB 实验会复现这个趋势。求解方式参数设置最优解 x*约束残差备注罚函数ρ10.3336.7e-1残差过大罚函数ρ1000.9802.0e-2残差勉强可用罚函数ρ100000.99982.0e-4条件数开始恶化增广拉格朗日ρ2外环 25 轮1.0000~2e-8ρ 保持有限这张表里增广拉格朗日那一行约束满足不是靠加大 ρ 逼出来的而是靠外层乘子更新把残余逐轮顶回零。ρ 在这里只负责给内层子问题提供合适曲率这正是它和纯罚函数的分水岭。2.3 从 L 到 L_ρ增广项为什么能稳定迭代增广拉格朗日函数写为 L_ρ(x,λ)f(x)λᵀc(x)(ρ/2)||c(x)||²。对比纯拉格朗日它多了最后的二次惩罚项对比纯罚函数它多了乘子线性项。两个性质让它特别适合迭代。其一对任意固定的 ρ0L_ρ 的驻点就是原问题的 KKT 候选点不会因 ρ 有限而系统性偏离可行方向。其二由内层驻点条件 ∇f(x_k)∇c(x_k)ᵀ(λ_kρc(x_k))0可以直接读出乘子估计 λ_kρc(x_k)于是乘子更新公式 λ_{k1}λ_kρc(x_k) 几乎是自然写出来的。把这段逻辑译成白话内层解完一个子问题后约束残差 c(x_k) 会残留在拉格朗日梯度里乘子更新就是把这部分残余收进λ下一轮子问题就在修正方向上继续。fmincon 的内点法和 SQP 内部也有类似修正步但增广拉格朗日把它显式放在外层每轮做了什么都能打印检查这是它长期活跃在科研代码和课程作业里的原因。3. 用 MATLAB 手写增广拉格朗日乘子法等式与不等式约束的最小可跑代码先说工具箱边界下面代码用到的 fminunc 属于 MATLAB 优化工具箱如果只有基础 MATLAB把内层求解换成 fminsearch 同样能跑通只是收敛慢一两个数量级。这个差别本身就是自写增广拉格朗日的理由之一——外层框架可以自由挂接不同的内层求解器。3.1 单等式约束问题的完整 MATLAB 实现代码处理 min fx1²x2²约束 x1x2−20解析解是 x(1,1)。故意从不可行点出发打印残差和乘子序列观察收敛过程。% augmented_lagrangian_eq.m % min f(x) x1^2 x2^2 s.t. c(x) x1 x2 - 2 0 % fminunc 需要 MATLAB 优化工具箱基础 MATLAB 可换成 fminsearch f (x) x(1)^2 x(2)^2; ceq (x) x(1) x(2) - 2; x0 [-1; 2]; % 初始点故意不取可行点 lambda 0; % 乘子初始值单等式约束所以是标量 rho 2; % 罚系数保持有限即可 tol_c 1e-8; % 约束残差停机阈值 history []; % 记录残差序列 for k 1:30 % 内层固定 lambda 与 rho做一次无约束极小化 Lag (z) f(z) lambda * ceq(z) 0.5 * rho * ceq(z)^2; opts optimoptions(fminunc, Display, off, ... Algorithm, quasi-newton); x fminunc(Lag, x0, opts); x0 x; % 热启动下一轮从本轮解出发 cval ceq(x); history(k, :) [k, abs(cval)]; fprintf(iter%2d x(%8.5f,%8.5f) c(x)%10.3e lambda%8.4f\n, ... k, x(1), x(2), cval, lambda); if abs(cval) tol_c break; end % 乘子更新区别于罚函数法的核心一行 lambda lambda rho * cval; end figure; semilogy(history(:,1), history(:,2), o-); xlabel(外层迭代 k); ylabel(|c(x)|); grid on;代码分三层理解。外层 for 循环每轮做三件事把当前 λ、ρ 封进匿名函数 Lag用 fminunc 求一次无约束极小化检查残差并更新乘子。内层用 quasi-newton 算法对初始点不敏感但每次调用要做有限差分梯度所以让 x 在上轮结果上继续迭代即热启动能显著减少内层工作量。history 记录每轮残差结尾的 semilogy 画出来是近似直线下降对应线性收敛速率。参数含义要分清rho 既是罚系数也决定乘子更新步长初始值取 1~10 较稳妥参考 c(x) 的典型量级lambda 取 0 即可知道对偶解近似值时可以取它热启动tol_c 是外层停机阈值常用 1e-6 到 1e-8。想压到 1e-10 就必须同时收紧内层 fminunc 的容差否则会被内层精度顶住。3.2 多等式约束与向量化写法约束变多时乘子 λ 和残差 c 都变成向量更新公式仍是 λλρc只是变成逐维推进。最容易犯的错是把增广项写成矩阵而忘了求和。% 两个等式约束c1 x1 x2 - 2; c2 x1 - x2 ceq_vec (x) [x(1) x(2) - 2; x(1) - x(2)]; lambda zeros(2, 1); % 乘子向量化 rho 2; % 增广项写作 0.5*rho*sum(ceq_vec(z).^2) Lag (z) f(z) lambda. * ceq_vec(z) 0.5 * rho * sum(ceq_vec(z).^2); % 乘子更新向量形式逐维向前推进 x fminunc(Lag, x0, opts); lambda lambda rho * ceq_vec(x);提示匿名函数里 ceq_vec(z) 被调用了两次约束函数计算昂贵时应先存成局部变量。MATLAB 匿名函数没有惰性求值这是约束次数上百时最常见的隐藏开销。3.3 不等式约束用 max 形式更新与截断一步完成对不等式 g(x)≤0乘子必须非负不能直接套等式公式。标准写法是定义增广拉格朗日函数 Φ(x,λ,ρ)f(x)(1/(2ρ))Σᵢ[max(0, λᵢρgᵢ(x))²−λᵢ²]乘子更新同步改为 λᵢ←max(0, λᵢρgᵢ(x))。% 不等式约束 g(x) x(1) - 0.8 0 的实现片段 g (x) x(1) - 0.8; rho 2; lambda 0; % 内层子问题max 形式自动判断约束是否活跃 LagI (z) f(z) (1/(2*rho)) * ... (max(0, lambda rho * g(z))^2 - lambda^2); x fminunc(LagI, x0, opts); % 乘子更新max 截断保证非负 lambda max(0, lambda rho * g(x));与等式约束版的差别只在 max 截断处。调试时可以借用这个性质不活跃约束的乘子会停在 0正好对应 KKT 互补松弛条件。混合等式与不等式约束时分别维护两个乘子向量内层目标把两类增广项加总即可。4. MATLAB 优化工具箱中的拉格朗日相关算法与参数怎么选从 R2016a 到最近版本fminunc、fmincon 的调用接口保持稳定下面内容可以直接对照你自己的 MATLAB 优化工具箱环境。4.1 fmincon 三种算法与增广拉格朗日的关系装了优化工具箱后最常用的入口是 fmincon。optimoptions 里有 sqp、interior-point、active-set 三种主要算法。它们都不是字面上的增广拉格朗日实现但内部都依赖子问题迭代 乘子类变量修正这一框架。算法适用规模每步做什么与增广拉格朗日的关系interior-point大中规模解带障碍函数的 KKT 系统障碍参数下降与乘子更新交替进行sqp中规模每步解一个 QP 子问题用拉格朗日函数海森近似含乘子修正active-set中小规模维护活跃约束集合迭代每步解等式约束子问题类似外层乘子思想选择建议目标和约束都光滑、变量几百个以内时sqp 通常最省事约束多或达到百万变量时interior-point 更稳已知活跃集或要反复热启动时active-set 合适。大规模内存受限时可以给 interior-point 设HessianApproximationlbfgs避免显式组装海森矩阵。既然有 fmincon为什么还要自写增广拉格朗日常见场景有三个目标函数来自仿真模型、梯度噪声大外层乘子慢更新能平滑噪声需要把子问题定向到特定内层算法或硬件论文或课程实验要输出乘子序列和残差曲线。这些场景下工具箱的封装反而碍事。4.2 关键参数表ρ、容差与停机条件的推荐值自写增广拉格朗日的参数只有五六个但每个都直接影响收敛。下面是我们常用的初始值参考。参数作用推荐初始值需要调整的信号rho子问题曲率与乘子步长1~10多轮外循环后残差不下降lambda乘子初值0 或对偶估计已知对偶解时取估计值加速tol_c约束残差停机阈值1e-6~1e-8残差达不到说明内层没解稳内层容差fminunc 停机精度1e-6~1e-8残差卡在固定量级rho 增长倍率外环放大系数2~5只在前一指标触发时使用注意ρ 只在约束残差下降不及预期时才放大不是每轮都放大。逐轮放大 ρ 等于退化成纯罚函数法。判断指标是相邻两轮残差比值若残差没有缩小到上一轮的 0.3~0.5 倍再把 ρ 乘 2~5。4.3 一个对照实验罚函数 vs 增广拉格朗日用同一个问题 min x² s.t. x1 做实验两者的差别一目了然代码可以直接复制到 MATLAB 里跑。% compare_penalty_vs_auglag.m f (x) x.^2; c (x) x - 1; fprintf(--- 纯罚函数提高 rho ---\n); for rho [1 10 100 1e4] xp fminbnd((x) f(x) 0.5*rho*c(x)^2, -1, 3); fprintf(rho%7g x%.6f 残差%.2e\n, rho, xp, abs(c(xp))); end fprintf(--- 增广拉格朗日固定 rho2更新乘子 ---\n); rho 2; lambda 0; x 0; for k 1:20 x fminbnd((x) f(x) lambda*c(x) 0.5*rho*c(x)^2, -1, 3); lambda lambda rho * (x - 1); fprintf(k%2d x%.8f 残差%.2e lambda%8.4f\n, ... k, x, abs(x - 1), lambda); end预期输出里罚函数要把 ρ 从 1 拉到 1e4残差才降到 2e-4增广拉格朗日固定 ρ2约 20 轮外迭代残差就低于 1e-6乘子 λ 同时收敛到 −2。一维问题用 fminbnd 是因为黄金分割法比有限差分更稳定这里的重点不是内层选型而是观察外层乘子更新如何在不放大 ρ 的情况下把残差压下去。搭出这个最小对照后再往多约束、非凸问题上迁移排查逻辑错误会容易得多。5. 增广拉格朗日收敛性验证、热启动与三个具体坑5.1 用 KKT 残差验证真收敛只看约束残差会误判。解满足 c(x)0 但不一定是极小点尤其目标函数非凸或初始点较远时。正确做法是同时检查驻点残差与可行残差KKT 残差 max(||∇f(x)∇c(x)ᵀλ||∞, ||c(x)||∞)。% KKT 残差检查接续第 3.1 节的例子 grad_f [2*x(1); 2*x(2)]; grad_c [1; 1]; stationarity norm(grad_f lambda * grad_c, Inf); feasibility abs(ceq(x)); kkt_res max(stationarity, feasibility); fprintf(KKT 残差%.3e (stationarity%.3e, feasibility%.3e)\n, ... kkt_res, stationarity, feasibility);判断标准两个分量都低于外层容差时才可以认为收敛。stationarity 降不下去而 feasibility 正常通常是内层子问题没解准反过来则是 ρ 或外环迭代不足。5.2 热启动把上一轮的 λ 与 x 直接传给下一轮热启动是增广拉格朗日对比罚函数的天然优势。进一步可以加乘子预热头两轮用较小的 ρ 跑让 λ 先进入合理量级再调到目标 ρ。非线性约束下 λ 初值不好会导致前几轮残差抖动这个预热能明显平滑过程。5.3 ρ 初值、内层精度与约束量纲三个容易踩的坑第一个坑是 ρ 初始值过大。ρ 很大会让内层子问题立即贴住可行域但乘子更新量 ρc(x) 也大外环容易在解附近来回摆动表现为残差不单调下降、λ 符号反复翻转。对策是先设小 ρ 跑两轮再逐步提上去。第二个坑是内层精度与外层容差不匹配。fminunc 默认停机容差在 1e-6 量级外层 tol_c 设 1e-10 时最后一轮会卡在内层极限上。要压到 1e-8 以下需要同步调StepTolerance和OptimalityTolerance或提供解析梯度改用 trust-region。第三个坑是多约束共用一个 ρ。位移与应力这类量纲不同的约束差几个数量级时量纲大的会在子问题里压倒一切。处理办法是给每个约束维护独立 ρᵢ只在残差下降慢的那一维上放大 ρᵢ相当于给乘子迭代加了自适应预条件。这个细节在参数辨识、形貌优化等混合量纲问题里出现频率很高建议从第一版代码就把 ρ 写成向量而不是标量。把这里的乘子框架再进一步对 x 做分裂就走向 ADMM但在引入 ADMM 之前先用上面的最小代码把 λ 序列和残差曲线画出来你会比直接调用工具箱更清楚每一步在修正什么。本文还有配套的精品资源点击获取 增广拉格朗日算法原理与MATLAB实现：从罚函数到乘子迭代
简介增广拉格朗日算法是求解带约束非线性优化问题的经典迭代方法它结合拉格朗日乘子法与惩罚函数法的优点通过引入惩罚项和增广乘子逐步逼近最优解被广泛应用于图像处理、机器学习、信号处理等领域。这份MATLAB源码针对这类需求编写特别适合需要处理非凸、非光滑或大规模优化问题的研究者与工程师可作为直接运行的参考实现。源码基于近端交替线性化最小化PALM框架在每一轮迭代中先求解线性化局部最小问题以更新决策变量再根据约束偏差更新拉格朗日乘子并动态增大惩罚系数以强化约束满足同时包含初始化、收敛判断等完整逻辑结构清晰、便于二次开发与教学演示。资源包非常精简仅包含1个m脚本文件总大小只有2KB可以无缝集成到已有项目中。目前已有3092人学习下载足以说明其实用性是学习增广拉格朗日算法和PALM方法的优质样例。1. 增广拉格朗日算法是约束优化的稳定主力MATLAB 里比罚函数更值得先掌握做约束优化时很多人第一反应是罚函数把等式约束乘一个大系数塞进目标函数权重从 10 调到 1e6解倒是接近可行域了海森矩阵却越来越病态迭代速度肉眼可见地变慢。增广拉格朗日算法Augmented Lagrangian Method也叫增广拉格朗日乘子法正好补上这一环它保留罚项的稳定结构但通过乘子迭代让罚系数保持有限就能精确收敛到可行点。对中等规模非线性受约束问题它实现简单、对初值不敏感思想上也与 MATLAB 优化工具箱内置算法同源。适合想自己控制迭代过程、做算法对比或写科研与课程代码的工程师和学生。2. 乘子法原理先说透拉格朗日函数、罚函数病态与 L_ρ 的构造增广拉格朗日不是天外来客它是在拉格朗日乘子法和罚函数法之间架起的一座桥。先把两侧概念对齐再理解中间多了什么后面的 MATLAB 代码才有依据。2.1 从拉格朗日乘子法与 KKT 条件说起等式约束的精确解在哪考虑光滑问题 min f(x)约束 c(x)0。定义拉格朗日函数 L(x,λ)f(x)λᵀc(x)λ 是维数等于约束个数的乘子向量。在最优点处KKT 条件给出两个方程∇f(x)∇c(x)ᵀλ0 与 c(x)0。前者是驻点条件后者是约束可行性联立解出的是局部最优的候选点还需要再配合二阶充分条件才能确认。直接用牛顿法解这个方程组小规模可行工程里却常受挫于两个现实一是 f、c 的解析梯度和海森往往拿不到仿真模型通常只能给出函数值二是牛顿迭代收敛半径窄初始点不佳就会发散。所以实用策略是把问题拆成两层内层做无约束极小化外层更新乘子。这正是增广拉格朗日算法的主框架也是后面 MATLAB 代码的骨架。2.2 纯罚函数的病态为什么 ρ 不能无限大罚函数法把问题改写成 min f(x)(ρ/2)||c(x)||²。理论保证 ρ→∞ 时罚函数的最优点趋于原问题最优点但工程上ρ→∞是个危险信号。ρ 变大意味着目标函数在约束面法向上越来越陡Hessian 条件数近似正比于 ρ梯度法步长被压小拟牛顿法也需要更多迭代才能在内层站稳。以 min x² s.t. x1 为例罚函数子问题解析解是 x*ρ/(ρ2)约束残差为 2/(ρ2)。要把残差压到 1e-4ρ 至少要取到上万此时子问题已经明显病态。增广拉格朗日固定 ρ2靠乘子外循环约 25 轮就能把同一个残差压到 2e-8ρ 自始至终不用变大。下表是同一问题下两类方法的数值对照第 4 章的 MATLAB 实验会复现这个趋势。求解方式参数设置最优解 x*约束残差备注罚函数ρ10.3336.7e-1残差过大罚函数ρ1000.9802.0e-2残差勉强可用罚函数ρ100000.99982.0e-4条件数开始恶化增广拉格朗日ρ2外环 25 轮1.0000~2e-8ρ 保持有限这张表里增广拉格朗日那一行约束满足不是靠加大 ρ 逼出来的而是靠外层乘子更新把残余逐轮顶回零。ρ 在这里只负责给内层子问题提供合适曲率这正是它和纯罚函数的分水岭。2.3 从 L 到 L_ρ增广项为什么能稳定迭代增广拉格朗日函数写为 L_ρ(x,λ)f(x)λᵀc(x)(ρ/2)||c(x)||²。对比纯拉格朗日它多了最后的二次惩罚项对比纯罚函数它多了乘子线性项。两个性质让它特别适合迭代。其一对任意固定的 ρ0L_ρ 的驻点就是原问题的 KKT 候选点不会因 ρ 有限而系统性偏离可行方向。其二由内层驻点条件 ∇f(x_k)∇c(x_k)ᵀ(λ_kρc(x_k))0可以直接读出乘子估计 λ_kρc(x_k)于是乘子更新公式 λ_{k1}λ_kρc(x_k) 几乎是自然写出来的。把这段逻辑译成白话内层解完一个子问题后约束残差 c(x_k) 会残留在拉格朗日梯度里乘子更新就是把这部分残余收进λ下一轮子问题就在修正方向上继续。fmincon 的内点法和 SQP 内部也有类似修正步但增广拉格朗日把它显式放在外层每轮做了什么都能打印检查这是它长期活跃在科研代码和课程作业里的原因。3. 用 MATLAB 手写增广拉格朗日乘子法等式与不等式约束的最小可跑代码先说工具箱边界下面代码用到的 fminunc 属于 MATLAB 优化工具箱如果只有基础 MATLAB把内层求解换成 fminsearch 同样能跑通只是收敛慢一两个数量级。这个差别本身就是自写增广拉格朗日的理由之一——外层框架可以自由挂接不同的内层求解器。3.1 单等式约束问题的完整 MATLAB 实现代码处理 min fx1²x2²约束 x1x2−20解析解是 x(1,1)。故意从不可行点出发打印残差和乘子序列观察收敛过程。% augmented_lagrangian_eq.m % min f(x) x1^2 x2^2 s.t. c(x) x1 x2 - 2 0 % fminunc 需要 MATLAB 优化工具箱基础 MATLAB 可换成 fminsearch f (x) x(1)^2 x(2)^2; ceq (x) x(1) x(2) - 2; x0 [-1; 2]; % 初始点故意不取可行点 lambda 0; % 乘子初始值单等式约束所以是标量 rho 2; % 罚系数保持有限即可 tol_c 1e-8; % 约束残差停机阈值 history []; % 记录残差序列 for k 1:30 % 内层固定 lambda 与 rho做一次无约束极小化 Lag (z) f(z) lambda * ceq(z) 0.5 * rho * ceq(z)^2; opts optimoptions(fminunc, Display, off, ... Algorithm, quasi-newton); x fminunc(Lag, x0, opts); x0 x; % 热启动下一轮从本轮解出发 cval ceq(x); history(k, :) [k, abs(cval)]; fprintf(iter%2d x(%8.5f,%8.5f) c(x)%10.3e lambda%8.4f\n, ... k, x(1), x(2), cval, lambda); if abs(cval) tol_c break; end % 乘子更新区别于罚函数法的核心一行 lambda lambda rho * cval; end figure; semilogy(history(:,1), history(:,2), o-); xlabel(外层迭代 k); ylabel(|c(x)|); grid on;代码分三层理解。外层 for 循环每轮做三件事把当前 λ、ρ 封进匿名函数 Lag用 fminunc 求一次无约束极小化检查残差并更新乘子。内层用 quasi-newton 算法对初始点不敏感但每次调用要做有限差分梯度所以让 x 在上轮结果上继续迭代即热启动能显著减少内层工作量。history 记录每轮残差结尾的 semilogy 画出来是近似直线下降对应线性收敛速率。参数含义要分清rho 既是罚系数也决定乘子更新步长初始值取 1~10 较稳妥参考 c(x) 的典型量级lambda 取 0 即可知道对偶解近似值时可以取它热启动tol_c 是外层停机阈值常用 1e-6 到 1e-8。想压到 1e-10 就必须同时收紧内层 fminunc 的容差否则会被内层精度顶住。3.2 多等式约束与向量化写法约束变多时乘子 λ 和残差 c 都变成向量更新公式仍是 λλρc只是变成逐维推进。最容易犯的错是把增广项写成矩阵而忘了求和。% 两个等式约束c1 x1 x2 - 2; c2 x1 - x2 ceq_vec (x) [x(1) x(2) - 2; x(1) - x(2)]; lambda zeros(2, 1); % 乘子向量化 rho 2; % 增广项写作 0.5*rho*sum(ceq_vec(z).^2) Lag (z) f(z) lambda. * ceq_vec(z) 0.5 * rho * sum(ceq_vec(z).^2); % 乘子更新向量形式逐维向前推进 x fminunc(Lag, x0, opts); lambda lambda rho * ceq_vec(x);提示匿名函数里 ceq_vec(z) 被调用了两次约束函数计算昂贵时应先存成局部变量。MATLAB 匿名函数没有惰性求值这是约束次数上百时最常见的隐藏开销。3.3 不等式约束用 max 形式更新与截断一步完成对不等式 g(x)≤0乘子必须非负不能直接套等式公式。标准写法是定义增广拉格朗日函数 Φ(x,λ,ρ)f(x)(1/(2ρ))Σᵢ[max(0, λᵢρgᵢ(x))²−λᵢ²]乘子更新同步改为 λᵢ←max(0, λᵢρgᵢ(x))。% 不等式约束 g(x) x(1) - 0.8 0 的实现片段 g (x) x(1) - 0.8; rho 2; lambda 0; % 内层子问题max 形式自动判断约束是否活跃 LagI (z) f(z) (1/(2*rho)) * ... (max(0, lambda rho * g(z))^2 - lambda^2); x fminunc(LagI, x0, opts); % 乘子更新max 截断保证非负 lambda max(0, lambda rho * g(x));与等式约束版的差别只在 max 截断处。调试时可以借用这个性质不活跃约束的乘子会停在 0正好对应 KKT 互补松弛条件。混合等式与不等式约束时分别维护两个乘子向量内层目标把两类增广项加总即可。4. MATLAB 优化工具箱中的拉格朗日相关算法与参数怎么选从 R2016a 到最近版本fminunc、fmincon 的调用接口保持稳定下面内容可以直接对照你自己的 MATLAB 优化工具箱环境。4.1 fmincon 三种算法与增广拉格朗日的关系装了优化工具箱后最常用的入口是 fmincon。optimoptions 里有 sqp、interior-point、active-set 三种主要算法。它们都不是字面上的增广拉格朗日实现但内部都依赖子问题迭代 乘子类变量修正这一框架。算法适用规模每步做什么与增广拉格朗日的关系interior-point大中规模解带障碍函数的 KKT 系统障碍参数下降与乘子更新交替进行sqp中规模每步解一个 QP 子问题用拉格朗日函数海森近似含乘子修正active-set中小规模维护活跃约束集合迭代每步解等式约束子问题类似外层乘子思想选择建议目标和约束都光滑、变量几百个以内时sqp 通常最省事约束多或达到百万变量时interior-point 更稳已知活跃集或要反复热启动时active-set 合适。大规模内存受限时可以给 interior-point 设HessianApproximationlbfgs避免显式组装海森矩阵。既然有 fmincon为什么还要自写增广拉格朗日常见场景有三个目标函数来自仿真模型、梯度噪声大外层乘子慢更新能平滑噪声需要把子问题定向到特定内层算法或硬件论文或课程实验要输出乘子序列和残差曲线。这些场景下工具箱的封装反而碍事。4.2 关键参数表ρ、容差与停机条件的推荐值自写增广拉格朗日的参数只有五六个但每个都直接影响收敛。下面是我们常用的初始值参考。参数作用推荐初始值需要调整的信号rho子问题曲率与乘子步长1~10多轮外循环后残差不下降lambda乘子初值0 或对偶估计已知对偶解时取估计值加速tol_c约束残差停机阈值1e-6~1e-8残差达不到说明内层没解稳内层容差fminunc 停机精度1e-6~1e-8残差卡在固定量级rho 增长倍率外环放大系数2~5只在前一指标触发时使用注意ρ 只在约束残差下降不及预期时才放大不是每轮都放大。逐轮放大 ρ 等于退化成纯罚函数法。判断指标是相邻两轮残差比值若残差没有缩小到上一轮的 0.3~0.5 倍再把 ρ 乘 2~5。4.3 一个对照实验罚函数 vs 增广拉格朗日用同一个问题 min x² s.t. x1 做实验两者的差别一目了然代码可以直接复制到 MATLAB 里跑。% compare_penalty_vs_auglag.m f (x) x.^2; c (x) x - 1; fprintf(--- 纯罚函数提高 rho ---\n); for rho [1 10 100 1e4] xp fminbnd((x) f(x) 0.5*rho*c(x)^2, -1, 3); fprintf(rho%7g x%.6f 残差%.2e\n, rho, xp, abs(c(xp))); end fprintf(--- 增广拉格朗日固定 rho2更新乘子 ---\n); rho 2; lambda 0; x 0; for k 1:20 x fminbnd((x) f(x) lambda*c(x) 0.5*rho*c(x)^2, -1, 3); lambda lambda rho * (x - 1); fprintf(k%2d x%.8f 残差%.2e lambda%8.4f\n, ... k, x, abs(x - 1), lambda); end预期输出里罚函数要把 ρ 从 1 拉到 1e4残差才降到 2e-4增广拉格朗日固定 ρ2约 20 轮外迭代残差就低于 1e-6乘子 λ 同时收敛到 −2。一维问题用 fminbnd 是因为黄金分割法比有限差分更稳定这里的重点不是内层选型而是观察外层乘子更新如何在不放大 ρ 的情况下把残差压下去。搭出这个最小对照后再往多约束、非凸问题上迁移排查逻辑错误会容易得多。5. 增广拉格朗日收敛性验证、热启动与三个具体坑5.1 用 KKT 残差验证真收敛只看约束残差会误判。解满足 c(x)0 但不一定是极小点尤其目标函数非凸或初始点较远时。正确做法是同时检查驻点残差与可行残差KKT 残差 max(||∇f(x)∇c(x)ᵀλ||∞, ||c(x)||∞)。% KKT 残差检查接续第 3.1 节的例子 grad_f [2*x(1); 2*x(2)]; grad_c [1; 1]; stationarity norm(grad_f lambda * grad_c, Inf); feasibility abs(ceq(x)); kkt_res max(stationarity, feasibility); fprintf(KKT 残差%.3e (stationarity%.3e, feasibility%.3e)\n, ... kkt_res, stationarity, feasibility);判断标准两个分量都低于外层容差时才可以认为收敛。stationarity 降不下去而 feasibility 正常通常是内层子问题没解准反过来则是 ρ 或外环迭代不足。5.2 热启动把上一轮的 λ 与 x 直接传给下一轮热启动是增广拉格朗日对比罚函数的天然优势。进一步可以加乘子预热头两轮用较小的 ρ 跑让 λ 先进入合理量级再调到目标 ρ。非线性约束下 λ 初值不好会导致前几轮残差抖动这个预热能明显平滑过程。5.3 ρ 初值、内层精度与约束量纲三个容易踩的坑第一个坑是 ρ 初始值过大。ρ 很大会让内层子问题立即贴住可行域但乘子更新量 ρc(x) 也大外环容易在解附近来回摆动表现为残差不单调下降、λ 符号反复翻转。对策是先设小 ρ 跑两轮再逐步提上去。第二个坑是内层精度与外层容差不匹配。fminunc 默认停机容差在 1e-6 量级外层 tol_c 设 1e-10 时最后一轮会卡在内层极限上。要压到 1e-8 以下需要同步调StepTolerance和OptimalityTolerance或提供解析梯度改用 trust-region。第三个坑是多约束共用一个 ρ。位移与应力这类量纲不同的约束差几个数量级时量纲大的会在子问题里压倒一切。处理办法是给每个约束维护独立 ρᵢ只在残差下降慢的那一维上放大 ρᵢ相当于给乘子迭代加了自适应预条件。这个细节在参数辨识、形貌优化等混合量纲问题里出现频率很高建议从第一版代码就把 ρ 写成向量而不是标量。把这里的乘子框架再进一步对 x 做分裂就走向 ADMM但在引入 ADMM 之前先用上面的最小代码把 λ 序列和残差曲线画出来你会比直接调用工具箱更清楚每一步在修正什么。本文还有配套的精品资源点击获取