简介这是一份讲解拉格朗日乘子法及其MATLAB实现的入门资源适合理工科学生、科研人员以及需要求解约束优化问题的工程师参考。资料从拉格朗日函数与KKT条件的基本原理讲起结合优化理论说明其在SVM、资源分配等场景中的应用并重点展示如何使用fmincon函数完成数值求解。资源共3个文件包含2个可直接运行的MATLAB脚本m文件和1个说明文档docx压缩包仅11KB内容精炼、便于快速上手。目前已有1204人学习下载。通过学习这份资料读者可以理解乘子法的数学推导过程掌握fmincon的调用方法与参数设置并借助示例代码迁移到自己的优化问题中是兼顾原理与实战的实用工具包。1. 拉格朗日乘子法遇上fmincon约束优化不再靠调罚函数“拉格朗日乘子法”这个名字做优化的人几乎天天见但真要在一个带约束的非线性问题上手写代码大多数人第一反应还是罚函数法把约束变成惩罚项然后调那个权重参数。这个参数就是个玄学调小了约束不满足调大了目标函数变得病态收敛慢得让人想砸键盘。fmincon 是 MATLAB 自带的约束优化求解器它内部用的核心思想正是拉格朗日乘子法——通过构造拉格朗日函数、盯住 KKT 条件把“有约束的最优化”转成一系列“无约束的迭代”最终同时给出最优解和乘子。这篇笔记面向两类人一类是正在做结构优化、参数辨识、经济调度被非线性约束折磨的工程师另一类是刚开始用 fmincon、想知道参数怎么设、结果怎么验证的科研人员。读完你能得到一套可以直接抄走的求解流程以及调试时最该看的几个坑。2. 从罚函数到拉格朗日乘子为什么KKT条件是fmincon的“裁判”2.1 罚函数法的问题惩罚因子那根拐杖很多工程师在遇到带约束的优化问题时第一反应是把约束罚进目标函数。比如一个典型问题 min f(x) s.t. c(x) ≤ 0用罚函数法会写成 min f(x) ρ·max(c(x), 0)²。这个形式直白代码也好写但有一个绕不开的毛病ρ 的取值没有通用规则。ρ 取小了最终解离可行域很远约束违例在工程上根本不能接受ρ 取大了目标函数在约束边界附近会变得极陡Hessian 矩阵的条件数随 ρ 增大而快速恶化梯度法在边界附近基本原地踏步。这里有个数量级上的直觉当目标函数值在 10³ 量级、约束违例在 10⁻² 量级时ρ 到底取 10⁴ 还是 10⁶结果可能差很远而且换一个初始点原先好用的 ρ 又失效了。更麻烦的是罚函数法给出的解总是不严格满足约束的它只是在惩罚项和目标函数之间找一个折中你永远说不清这个解是“真最优”还是“罚出来的”。这就是大家常说的“罚函数调参翻车”不是数学不对而是工程上没法稳定使用。拉格朗日乘子法的思路完全不同。它不把惩罚因子推向无穷而是把约束“原样”放进拉格朗日函数里通过迭代更新乘子 λ让解在满足约束的同时逼近最优。用直白的话说罚函数法是靠一根越拧越紧的拐杖把解逼到可行域乘子法则是让约束和乘子互相纠偏。fmincon 的 interior-point 算法和 SQP 算法骨子里都是这一套先写出拉格朗日函数再对 KKT 条件做迭代求解而不是去猜某个惩罚权重。2.2 拉格朗日函数与KKT条件手算一遍就懂给定一般形式的非线性规划问题min f(x)s.t. c_i(x) ≤ 0i 1…mh_j(x) 0j 1…p拉格朗日函数定义为L(x, λ, μ) f(x) Σ λ_i·c_i(x) Σ μ_j·h_j(x)如果 x* 是局部最优解并且约束满足正则性条件通常指约束梯度线性无关即 LICQ那么一定存在乘子 λ* ≥ 0 和 μ*使得∇f(x*) Σ λ_i·∇c_i(x*) Σ μ_j·∇h_j(x*) 0λ_i·c_i(x*) 0对所有 i 成立第二组式子叫互补松弛条件。它有一个非常直观的工程解释一个不等式约束要么它不活跃此时 λ_i 0约束形同虚设要么它被“压死”在边界上此时 c_i(x*) 0乘子才能非零。fmincon 的迭代本质就是在找一个满足上述条件的点每次迭代都会输出一个一阶最优性度量optimality measure这个度量就是上面第一个式子的残差范数。举一个能用手算验证的例子后面第 3 章的代码也会复用min f(x) (x₁−2)² (x₂−1)²s.t. x₁² x₂² ≤ 1目标梯度是 [2x₁−4, 2x₂−2]约束梯度是 [2x₁, 2x₂]。由于最优解一定落在单位圆边界上互补松弛条件退化为约束取等号。联立梯度方程和 x₁²x₂²1可以解出 x* (2/√5, 1/√5) ≈ (0.8944, 0.4472)乘子 λ* √5−1 ≈ 1.236目标函数值 f* 6−2√5 ≈ 1.5279。这个手算结果就是后面验证 fmincon 输出的标尺。2.3 乘子值本身是有用的“副产品”fmincon 在返回最优解 x 的同时会返回一个 lambda 结构体这就是拉格朗日乘子的数值估计。对不等式约束lambda.ineqnonlin 是非负的对等式约束lambda.eqnonlin 可正可负。很多使用者把这个结构体忽略掉只拿 x 和 fval这是很亏的乘子的大小直接告诉你哪个约束在“卡脖子”。乘子越大说明这个约束稍微放松一点点目标函数就能改善很多这就是所谓的“影子价格”。需要提醒的是乘子的符号和尺度。fmincon 内部把不等式统一写成 c(x) ≤ 0 的标准形式所以返回的乘子非负如果你手算时习惯把约束写成 g(x) ≥ 0那乘子会差一个符号。另外目标函数整体放大 10 倍乘子也会放大 10 倍所以拿两个不同量纲的问题去比较乘子绝对值没有意义必须先归一化。3. 用fmincon跑通拉格朗日乘子法最小可复现的约束优化代码3.1 目标函数与约束的MATLAB描述先跑出第一个结果这一节直接进入能复现的代码。用第 2 章手算过的那个问题目标函数和约束各写成一个函数文件。这是 fmincon 的标准接口所有非线性问题都按这个套路来。% objfun.m function f objfun(x) % 目标函数min (x1-2)^2 (x2-1)^2 f (x(1)-2)^2 (x(2)-1)^2; end% confun.m function [c, ceq] confun(x) % 不等式约束 x1^2x2^2 1写成 c(x) 0 的标准形式 c x(1)^2 x(2)^2 - 1; ceq []; % 本例没有等式约束 end目标函数返回一个标量约束函数返回两个输出c 是不等式约束向量每个元素必须 ≤ 0ceq 是等式约束向量每个元素必须等于 0。如果没有等式约束ceq 返回空数组。注意约束函数的顺序不能换fmincon 内部依赖这个固定接口。接下来是主脚本调用 fmincon% solve_lagrange.m x0 [0.5; 0.5]; % 初始点尽量靠近可行域 options optimoptions(fmincon, ... Algorithm, interior-point, ... Display, iter, ... SpecifyObjectiveGradient, false, ... SpecifyConstraintGradient, false); [x, fval, exitflag, output, lambda] fmincon(objfun, x0, ... [], [], [], [], [], [], confun, options); fprintf(x [%.6f, %.6f]\n, x(1), x(2)); fprintf(fval %.6f\n, fval); fprintf(exitflag %d\n, exitflag); fprintf(lambda.ineqnonlin %.6f\n, lambda.ineqnonlin);fmincon 的函数签名依次是目标函数句柄、初始点 x0、线性不等式矩阵 A 和向量 b、线性等式矩阵 Aeq 和向量 beq、下界 lb、上界 ub、非线性约束函数。本例全部线性约束和边界都没有所以用空数组占位。非线性约束必须放在最后一位写成函数句柄。运行后x 应该约等于 (0.8944, 0.4472)fval 约等于 1.5279lambda.ineqnonlin 约等于 1.2361exitflag 为 1表示一阶最优性条件满足。这个结果和第 2 节的手算完全对得上说明求解器没有跑偏。3.2 给fmincon喂解析梯度收敛速度和精度都能上一个台阶上面的运行使用的是有限差分梯度fmincon 内部对每个变量做扰动来估计导数。问题是有限差分的精度受限于步长选择步长太大会引入截断误差太小会引入舍入误差一般只能把“最优性”做到 1e-4 到 1e-5 量级。如果你的问题对精度要求高或者变量维数大应该把解析梯度写进去。下面给出带梯度的目标函数和约束函数版本% objfun_grad.m function [f, g] objfun_grad(x) f (x(1)-2)^2 (x(2)-1)^2; if nargout 1 % 目标梯度是列向量维度与 x 相同 g [2*x(1)-4; 2*x(2)-2]; end end% confun_grad.m function [c, ceq, gradc, gradceq] confun_grad(x) c x(1)^2 x(2)^2 - 1; ceq []; if nargout 2 % 约束雅可比矩阵每列对应一个不等式约束的偏导数 gradc [2*x(1); 2*x(2)]; gradceq []; end end这里有一个非常关键的约定目标函数的梯度是列向量维度是 n×1非线性不等式约束的梯度是一个 n×m 矩阵n 是变量个数m 是不等式约束个数每一列对应一个约束函数对所有变量的偏导。如果只有一个不等式约束gradc 就是一个 n×1 列向量也就是上面代码的写法。如果写了多个约束比如 c [c1; c2]那么 gradc 要写成 [∇c1, ∇c2]即第一列是 c1 的偏导第二列是 c2 的偏导这个顺序不能错。写完之后在 options 里打开两个开关options optimoptions(fmincon, ... Algorithm, interior-point, ... SpecifyObjectiveGradient, true, ... SpecifyConstraintGradient, true);打开开关后fmincon 就不再调用有限差分直接用你给的梯度去构造 KKT 条件。这时 OptimalityTolerance 的设置才有意义——否则数值差分误差本身就有 1e-5 量级你却把容差设成 1e-10纯属自己骗自己不仅迭代不收敛还会在边界附近来回震荡。3.3 梯度写没写错用中心差分做一次体检解析梯度写错是 fmincon 调参中最常见的事故而且症状常常很隐蔽不是立刻报错而是收敛变慢、结果偏移、甚至直接发散。我一般会先跑一段中心差分对比代码成本极低能省下后面几个小时的排查时间。% check_grad.m x [0.4; 0.8]; % 任取一个不在最优点的检验点 h 1e-6; % 中心差分步长 % 解析梯度 [f_parsed, g_parsed] objfun_grad(x); % 数值梯度 g_num zeros(2, 1); for i 1:2 xp x; xp(i) xp(i) h; xm x; xm(i) xm(i) - h; g_num(i) (objfun(xp) - objfun(xm)) / (2 * h); end disp(解析梯度 数值梯度); disp([g_parsed, g_num]);如果两列数值差异超过 1e-4说明解析梯度大概率写错了。常见的错误有两类一是把目标梯度写成了行向量二是约束雅可比的行列顺序搞反了。梯度校验通过后再把解析梯度用于正式求解这时候 fmincon 的迭代步数通常会明显减少收敛精度能压到 1e-8 以下。3.4 验证KKT条件光看exitflag不够exitflag1 只代表求解器内部判定“一阶最优性条件满足”但这个结论建立在它自己的梯度估计上。如果用了解析梯度可以做一次独立的 KKT 残差检查把求解结果拿到解析公式里算一遍% verify_kkt.m [x, fval, ~, ~, lambda] fmincon(objfun_grad, x0, ... [], [], [], [], [], [], confun_grad, options); grad_f [2*x(1)-4; 2*x(2)-2]; grad_c [2*x(1); 2*x(2)]; % KKT 梯度平衡条件 residual grad_f lambda.ineqnonlin * grad_c; fprintf(KKT residual norm %e\n, norm(residual)); % 互补松弛条件检查 fprintf(c(x) %.6e, lambda %.6f\n, x(1)^2x(2)^2-1, lambda.ineqnonlin);如果 residual 的范数在 1e-6 量级而 c(x) 在 -1e-8 到 0 之间说明解确实站在约束边界上并且乘子符号和大小都合理。这个验证过程不仅是在检验 fmincon更是在检验你对问题的建模是否正确。4. fmincon算法与参数调优内点法、SQP和五个必调选项4.1 Algorithm怎么选interior-point、sqp还是active-setfmincon 的 Algorithm 选项经常让人困惑因为名字里没有“拉格朗日乘子法”。实际上内点法和 SQP 都属于乘子法框架区别在于各自如何构造迭代子问题。算法适用场景特点典型坑interior-point大规模问题、不等式约束多用障碍函数把不等式约束“软化”迭代稳定可处理稀疏 Hessian初始点在可行域外时可能多花一个恢复阶段sqp中规模问题、约束梯度便宜每步解一个二次规划子问题收敛快边界跟踪直接每步代价高几十个变量时明显变慢active-set小规模问题、想观察乘子变化过程迭代路径可解释接近经典乘子法的操作方式对初始点敏感约束多时容易抖动trust-region-reflective只有边界约束或线性等式约束需要提供梯度不能直接处理一般非线性不等式约束很多人误用它去解不等式约束问题直接报错选型口诀可以这么说问题变量少、约束个数在个位、希望结果稳定可复现选 sqp问题规模大、不等式约束几十上百条、只关心最终结果选 interior-point做算法对比研究、想看清楚每一步哪个约束被激活选 active-set。这个选择可以先用你的真实数据分别跑一遍对比迭代步数和总耗时再定不丢人。4.2 五个必调参数它们各自管什么fmincon 的参数面板里有十几个选项但实际工程中真正需要动的就五个。其余保持默认往往比乱调更健康。参数默认值作用建议TolFun1e-6目标函数两次迭代的变化容差常规问题 1e-7需要高精度时 1e-9但必须配合解析梯度TolCon1e-6约束违反容差衡量“不可行度”结构安全类问题设 1e-8约束必须严格满足OptimalityTolerance1e-6KKT 一阶最优性残差容差不要小于数值差分能支持的量级有限差分时建议 1e-5 到 1e-6StepTolerance1e-10迭代步长容差判断“走不动了”一般不用动如果 exitflag2 且约束未满足说明卡在不可行点MaxIterations / MaxFunctionEvaluations与问题默认相关迭代上限和函数计算上限复杂问题常需要放大到 1000~3000但不建议盲目放大TolCon 和 OptimalityTolerance 是最容易被混淆的一对。TolCon 管的是“解离可行域多远”OptimalityTolerance 管的是“解离真最优多远”。举一个实际场景你的解约束违例只有 1e-7满足 TolCon1e-6但目标函数还有 0.5 的改善空间原因是 OptimalityTolerance 设得太大求解器提前收工了。反过来目标函数收敛得很好但约束违例 1e-5在安全评估里照样不合格。所以每次求解完这两个量都要打印出来看。4.3 读懂Displayiter的输出NLP constraint violation是硬指标设置 Displayiter 后控制台会打印一张迭代表。每一行依次是迭代步数、函数计算次数、目标函数值、NLP 约束违例量、一阶最优性度量、步长。这张表里最需要盯的是第二列往后的趋势目标函数是否在连续下降约束违例是否在朝 0 收敛。如果看到约束违例一直在 1e-2 附近抖动说明求解器正在可行域边界上“找缝”如果约束违例稳步下降但目标函数始终不动大概率是约束条件跟目标函数方向矛盾需要检查建模而不是调参数。exitflag 的速查也有用1 表示一阶最优性满足0 表示迭代步数或计算次数超限-2 表示问题不可行2 表示步长小于 StepTolerance 但约束违例不满足。看到 exitflag0 时先把 MaxIterations 提到 2000 再试一次如果还是 0问题几乎可以确定出在约束本身而不是参数上。5. 拉格朗日乘子法fmincon避坑清单五个常见翻车现场与排查5.1 现象exitflag0控制台显示迭代一步都没动我遇到过好几次这种情况目标函数写得没问题约束也没问题但 fmincon 从第一步开始就不收敛迭代次数一直跳到上限。原因多半是初始点离可行域太远。interior-point 算法有一个“恢复阶段”专门处理不可行初始点但如果约束是非线性且初始点严重违例这个阶段可能反复失败。解决思路是先做一个“可行性恢复”单独求解一个只求可行、不管目标的最优化问题以它的解作为新的初始点。% feasible_start.m function [c, ceq] feasible_confun(x) % 原约束不变 c x(1)^2 x(2)^2 - 1; ceq []; end % 目标函数临时改成“惩罚所有约束违例” obj_feas (x) sum(max(confun_constraint(x), 0).^2);更简单的方法是直接换初始点或者先用 sqp 算法跑一次它对不可行初始点的容忍度往往更高再用返回的 x 作为 interior-point 的初始点。这个技巧在工程问题里非常实用。5.2 现象开了解析梯度后结果反而变差甚至发散这是个特别容易让人怀疑人生的场景不给梯度时结果还行给了梯度反而崩了于是有人得出“解析梯度没用”的结论。真相几乎总是梯度写错。原因就出在符号约定和矩阵形状上。目标梯度写成行向量是重灾区约束雅可比的行列顺序搞反也很常见。尤其是多个不等式约束时gradc 的每一列对应一个约束而不是每一行。解决方法是跑一次 3.3 节的中心差分对比把解析梯度和数值梯度打印出来逐项看。血泪经验梯度错了fmincon 的 KKT 条件就建立在错误信息上它不但不会更快反而会朝着错误方向迭代。这不是求解器的锅是人的锅。5.3 现象手算的拉格朗日乘子与lambda.ineqnonlin对不上做完手算发现差个负号或者数值差了好几倍这是非常常见的核对事故。原因之一是符号约定。fmincon 把不等式约束写成 c(x) ≤ 0乘子取非负而你手算时可能习惯写成 x₁²x₂² ≥ 1 这种形式乘子自然反号。原因之二是尺度问题目标函数放大 10 倍乘子跟着放大 10 倍这不是 bug是数学性质。解决方法是先把所有约束统一改写成标准形式 c(x) ≤ 0再手算。数值对不上时优先检查是不是单位换算比如约束用的是兆帕还是帕其次再检查梯度方向。记住乘子不是固定的物理常数它依赖你建模时的尺度核对的唯一依据是 KKT 方程本身而不是直觉。5.4 现象约束条件在数量级上差距太大怎么调都别扭比如一个约束的数值是 1e6另一个是 1e-3目标函数又是 1e2。这种情况下 fmincon 的 KKT 系统很容易病态迭代在几个约束之间反复横跳。原因在数学上很明确KKT 系统里各项的量级差太多数值求解的舍入误差被放大。这不是 fmincon 的问题而是模型的问题。解决方法是做变量和约束的归一化。工程上常见的做法是把所有变量缩放到 0.1~10 区间约束也做无量纲化处理。比如应力约束写成“当前应力/许用应力 ≤ 1”而不是“当前应力 ≤ 许用应力”。归一化之后乘子的大小才有了跨问题比较的意义求解器的稳定性也会明显提升。5.5 现象某个约束乘子为零但你觉得它明明在起作用这可能是拉格朗日乘子法应用中最容易误判的一个坑。你看到 x* 在某条约束边界上但 lambda.ineqnonlin 对应分量是 0于是怀疑求解器出了问题。原因通常是约束冗余即当前点有多条约束同时活跃但它们对应的梯度线性相关导致 LICQ 条件不满足KKT 乘子存在但不唯一。求解器会给其中一条约束分配乘子另一条给 0。也有可能是约束在该点不可微比如绝对值约束或 max 型约束。解决方法是检查约束雅可比矩阵在解处的秩。如果秩小于活跃约束个数说明有冗余约束删掉一条再求解。对不可微约束需要做平滑近似比如用 sqrt(c²ε) 代替 |c|ε 取 1e-4 量级。乘子为零不一定说明约束没用更准确的说法是“这条约束在当前点提供的信息已经被其他约束覆盖了”。6. 用拉格朗日乘子做敏感性分析验证最优解的一把尺子6.1 用乘子预测“放松约束带来的收益”fmincon 返回的 lambda 结构体不是摆设它可以直接用来做敏感性分析。对第 3 章的示例问题不等式约束是 x₁²x₂² ≤ 1。如果我把约束右端项从 1 改成 1.001可行域稍微放大一点最优目标函数值会下降多少根据 KKT 理论目标函数对约束右端项的偏导数约等于乘子的负值。本例 λ ≈ 1.236所以可行域放大 0.001 后目标函数值大约下降 1.236×0.001 ≈ 0.00124。可以写一个脚本验证这个预测。% sensitivity.m % 把约束半径平方作为参数传入 function [c, ceq] confun_radius(x, r_sq) c x(1)^2 x(2)^2 - r_sq; ceq []; end % 用原半径求一次 [x1, fval1] run_with_radius(1.0); lambda1 get_lambda(1.0); % 用放大后的半径再求一次 delta 0.001; [x2, fval2] run_with_radius(1.0 delta); % 比较预测与实际 predicted fval1 - lambda1 * delta; fprintf(预测下降: %.6f\n, fval1 - predicted); fprintf(实际下降: %.6f\n, fval1 - fval2);如果预测值和实际值在同一个量级说明乘子不仅准而且整个 KKT 系统是自洽的。这个验证方法可以在每次求解完之后当作例行检查既能发现梯度写错、约束建模错误也能帮你理解当前设计到底被哪条约束卡住。6.2 多约束问题里如何快速判断“瓶颈约束”当问题有几十条约束时逐个去看约束违例非常低效。我的习惯是直接打印所有活跃约束的乘子按乘子绝对值排序排在前面的一定是当前优化中最“值钱”的约束。比如结构减重问题里如果应力约束的乘子远大于位移约束那说明继续减重的瓶颈在应力改材料比改截面更划算。用乘子做敏感性分析还有一个好处它把“定性判断”变成“定量判断”。不需要反复重跑优化来试只要看乘子就知道哪个约束值得去改善。这个方法我在多个项目里用过确实能省不少计算时间。早期我做优化时只盯着 x 和 fval从不看 lambda结果换一组初值就翻车后来养成了习惯每次求解完必打印 lambda再做一次约束扰动验证核完才敢拿结果去出图。希望这个习惯对你也有帮助。本文还有配套的精品资源点击获取 
长河君 · 2026/10/26 7:45:26 · 阅读 1.2k 拉格朗日乘子法与fmincon：约束优化从罚函数到KKT条件实战
简介这是一份讲解拉格朗日乘子法及其MATLAB实现的入门资源适合理工科学生、科研人员以及需要求解约束优化问题的工程师参考。资料从拉格朗日函数与KKT条件的基本原理讲起结合优化理论说明其在SVM、资源分配等场景中的应用并重点展示如何使用fmincon函数完成数值求解。资源共3个文件包含2个可直接运行的MATLAB脚本m文件和1个说明文档docx压缩包仅11KB内容精炼、便于快速上手。目前已有1204人学习下载。通过学习这份资料读者可以理解乘子法的数学推导过程掌握fmincon的调用方法与参数设置并借助示例代码迁移到自己的优化问题中是兼顾原理与实战的实用工具包。1. 拉格朗日乘子法遇上fmincon约束优化不再靠调罚函数“拉格朗日乘子法”这个名字做优化的人几乎天天见但真要在一个带约束的非线性问题上手写代码大多数人第一反应还是罚函数法把约束变成惩罚项然后调那个权重参数。这个参数就是个玄学调小了约束不满足调大了目标函数变得病态收敛慢得让人想砸键盘。fmincon 是 MATLAB 自带的约束优化求解器它内部用的核心思想正是拉格朗日乘子法——通过构造拉格朗日函数、盯住 KKT 条件把“有约束的最优化”转成一系列“无约束的迭代”最终同时给出最优解和乘子。这篇笔记面向两类人一类是正在做结构优化、参数辨识、经济调度被非线性约束折磨的工程师另一类是刚开始用 fmincon、想知道参数怎么设、结果怎么验证的科研人员。读完你能得到一套可以直接抄走的求解流程以及调试时最该看的几个坑。2. 从罚函数到拉格朗日乘子为什么KKT条件是fmincon的“裁判”2.1 罚函数法的问题惩罚因子那根拐杖很多工程师在遇到带约束的优化问题时第一反应是把约束罚进目标函数。比如一个典型问题 min f(x) s.t. c(x) ≤ 0用罚函数法会写成 min f(x) ρ·max(c(x), 0)²。这个形式直白代码也好写但有一个绕不开的毛病ρ 的取值没有通用规则。ρ 取小了最终解离可行域很远约束违例在工程上根本不能接受ρ 取大了目标函数在约束边界附近会变得极陡Hessian 矩阵的条件数随 ρ 增大而快速恶化梯度法在边界附近基本原地踏步。这里有个数量级上的直觉当目标函数值在 10³ 量级、约束违例在 10⁻² 量级时ρ 到底取 10⁴ 还是 10⁶结果可能差很远而且换一个初始点原先好用的 ρ 又失效了。更麻烦的是罚函数法给出的解总是不严格满足约束的它只是在惩罚项和目标函数之间找一个折中你永远说不清这个解是“真最优”还是“罚出来的”。这就是大家常说的“罚函数调参翻车”不是数学不对而是工程上没法稳定使用。拉格朗日乘子法的思路完全不同。它不把惩罚因子推向无穷而是把约束“原样”放进拉格朗日函数里通过迭代更新乘子 λ让解在满足约束的同时逼近最优。用直白的话说罚函数法是靠一根越拧越紧的拐杖把解逼到可行域乘子法则是让约束和乘子互相纠偏。fmincon 的 interior-point 算法和 SQP 算法骨子里都是这一套先写出拉格朗日函数再对 KKT 条件做迭代求解而不是去猜某个惩罚权重。2.2 拉格朗日函数与KKT条件手算一遍就懂给定一般形式的非线性规划问题min f(x)s.t. c_i(x) ≤ 0i 1…mh_j(x) 0j 1…p拉格朗日函数定义为L(x, λ, μ) f(x) Σ λ_i·c_i(x) Σ μ_j·h_j(x)如果 x* 是局部最优解并且约束满足正则性条件通常指约束梯度线性无关即 LICQ那么一定存在乘子 λ* ≥ 0 和 μ*使得∇f(x*) Σ λ_i·∇c_i(x*) Σ μ_j·∇h_j(x*) 0λ_i·c_i(x*) 0对所有 i 成立第二组式子叫互补松弛条件。它有一个非常直观的工程解释一个不等式约束要么它不活跃此时 λ_i 0约束形同虚设要么它被“压死”在边界上此时 c_i(x*) 0乘子才能非零。fmincon 的迭代本质就是在找一个满足上述条件的点每次迭代都会输出一个一阶最优性度量optimality measure这个度量就是上面第一个式子的残差范数。举一个能用手算验证的例子后面第 3 章的代码也会复用min f(x) (x₁−2)² (x₂−1)²s.t. x₁² x₂² ≤ 1目标梯度是 [2x₁−4, 2x₂−2]约束梯度是 [2x₁, 2x₂]。由于最优解一定落在单位圆边界上互补松弛条件退化为约束取等号。联立梯度方程和 x₁²x₂²1可以解出 x* (2/√5, 1/√5) ≈ (0.8944, 0.4472)乘子 λ* √5−1 ≈ 1.236目标函数值 f* 6−2√5 ≈ 1.5279。这个手算结果就是后面验证 fmincon 输出的标尺。2.3 乘子值本身是有用的“副产品”fmincon 在返回最优解 x 的同时会返回一个 lambda 结构体这就是拉格朗日乘子的数值估计。对不等式约束lambda.ineqnonlin 是非负的对等式约束lambda.eqnonlin 可正可负。很多使用者把这个结构体忽略掉只拿 x 和 fval这是很亏的乘子的大小直接告诉你哪个约束在“卡脖子”。乘子越大说明这个约束稍微放松一点点目标函数就能改善很多这就是所谓的“影子价格”。需要提醒的是乘子的符号和尺度。fmincon 内部把不等式统一写成 c(x) ≤ 0 的标准形式所以返回的乘子非负如果你手算时习惯把约束写成 g(x) ≥ 0那乘子会差一个符号。另外目标函数整体放大 10 倍乘子也会放大 10 倍所以拿两个不同量纲的问题去比较乘子绝对值没有意义必须先归一化。3. 用fmincon跑通拉格朗日乘子法最小可复现的约束优化代码3.1 目标函数与约束的MATLAB描述先跑出第一个结果这一节直接进入能复现的代码。用第 2 章手算过的那个问题目标函数和约束各写成一个函数文件。这是 fmincon 的标准接口所有非线性问题都按这个套路来。% objfun.m function f objfun(x) % 目标函数min (x1-2)^2 (x2-1)^2 f (x(1)-2)^2 (x(2)-1)^2; end% confun.m function [c, ceq] confun(x) % 不等式约束 x1^2x2^2 1写成 c(x) 0 的标准形式 c x(1)^2 x(2)^2 - 1; ceq []; % 本例没有等式约束 end目标函数返回一个标量约束函数返回两个输出c 是不等式约束向量每个元素必须 ≤ 0ceq 是等式约束向量每个元素必须等于 0。如果没有等式约束ceq 返回空数组。注意约束函数的顺序不能换fmincon 内部依赖这个固定接口。接下来是主脚本调用 fmincon% solve_lagrange.m x0 [0.5; 0.5]; % 初始点尽量靠近可行域 options optimoptions(fmincon, ... Algorithm, interior-point, ... Display, iter, ... SpecifyObjectiveGradient, false, ... SpecifyConstraintGradient, false); [x, fval, exitflag, output, lambda] fmincon(objfun, x0, ... [], [], [], [], [], [], confun, options); fprintf(x [%.6f, %.6f]\n, x(1), x(2)); fprintf(fval %.6f\n, fval); fprintf(exitflag %d\n, exitflag); fprintf(lambda.ineqnonlin %.6f\n, lambda.ineqnonlin);fmincon 的函数签名依次是目标函数句柄、初始点 x0、线性不等式矩阵 A 和向量 b、线性等式矩阵 Aeq 和向量 beq、下界 lb、上界 ub、非线性约束函数。本例全部线性约束和边界都没有所以用空数组占位。非线性约束必须放在最后一位写成函数句柄。运行后x 应该约等于 (0.8944, 0.4472)fval 约等于 1.5279lambda.ineqnonlin 约等于 1.2361exitflag 为 1表示一阶最优性条件满足。这个结果和第 2 节的手算完全对得上说明求解器没有跑偏。3.2 给fmincon喂解析梯度收敛速度和精度都能上一个台阶上面的运行使用的是有限差分梯度fmincon 内部对每个变量做扰动来估计导数。问题是有限差分的精度受限于步长选择步长太大会引入截断误差太小会引入舍入误差一般只能把“最优性”做到 1e-4 到 1e-5 量级。如果你的问题对精度要求高或者变量维数大应该把解析梯度写进去。下面给出带梯度的目标函数和约束函数版本% objfun_grad.m function [f, g] objfun_grad(x) f (x(1)-2)^2 (x(2)-1)^2; if nargout 1 % 目标梯度是列向量维度与 x 相同 g [2*x(1)-4; 2*x(2)-2]; end end% confun_grad.m function [c, ceq, gradc, gradceq] confun_grad(x) c x(1)^2 x(2)^2 - 1; ceq []; if nargout 2 % 约束雅可比矩阵每列对应一个不等式约束的偏导数 gradc [2*x(1); 2*x(2)]; gradceq []; end end这里有一个非常关键的约定目标函数的梯度是列向量维度是 n×1非线性不等式约束的梯度是一个 n×m 矩阵n 是变量个数m 是不等式约束个数每一列对应一个约束函数对所有变量的偏导。如果只有一个不等式约束gradc 就是一个 n×1 列向量也就是上面代码的写法。如果写了多个约束比如 c [c1; c2]那么 gradc 要写成 [∇c1, ∇c2]即第一列是 c1 的偏导第二列是 c2 的偏导这个顺序不能错。写完之后在 options 里打开两个开关options optimoptions(fmincon, ... Algorithm, interior-point, ... SpecifyObjectiveGradient, true, ... SpecifyConstraintGradient, true);打开开关后fmincon 就不再调用有限差分直接用你给的梯度去构造 KKT 条件。这时 OptimalityTolerance 的设置才有意义——否则数值差分误差本身就有 1e-5 量级你却把容差设成 1e-10纯属自己骗自己不仅迭代不收敛还会在边界附近来回震荡。3.3 梯度写没写错用中心差分做一次体检解析梯度写错是 fmincon 调参中最常见的事故而且症状常常很隐蔽不是立刻报错而是收敛变慢、结果偏移、甚至直接发散。我一般会先跑一段中心差分对比代码成本极低能省下后面几个小时的排查时间。% check_grad.m x [0.4; 0.8]; % 任取一个不在最优点的检验点 h 1e-6; % 中心差分步长 % 解析梯度 [f_parsed, g_parsed] objfun_grad(x); % 数值梯度 g_num zeros(2, 1); for i 1:2 xp x; xp(i) xp(i) h; xm x; xm(i) xm(i) - h; g_num(i) (objfun(xp) - objfun(xm)) / (2 * h); end disp(解析梯度 数值梯度); disp([g_parsed, g_num]);如果两列数值差异超过 1e-4说明解析梯度大概率写错了。常见的错误有两类一是把目标梯度写成了行向量二是约束雅可比的行列顺序搞反了。梯度校验通过后再把解析梯度用于正式求解这时候 fmincon 的迭代步数通常会明显减少收敛精度能压到 1e-8 以下。3.4 验证KKT条件光看exitflag不够exitflag1 只代表求解器内部判定“一阶最优性条件满足”但这个结论建立在它自己的梯度估计上。如果用了解析梯度可以做一次独立的 KKT 残差检查把求解结果拿到解析公式里算一遍% verify_kkt.m [x, fval, ~, ~, lambda] fmincon(objfun_grad, x0, ... [], [], [], [], [], [], confun_grad, options); grad_f [2*x(1)-4; 2*x(2)-2]; grad_c [2*x(1); 2*x(2)]; % KKT 梯度平衡条件 residual grad_f lambda.ineqnonlin * grad_c; fprintf(KKT residual norm %e\n, norm(residual)); % 互补松弛条件检查 fprintf(c(x) %.6e, lambda %.6f\n, x(1)^2x(2)^2-1, lambda.ineqnonlin);如果 residual 的范数在 1e-6 量级而 c(x) 在 -1e-8 到 0 之间说明解确实站在约束边界上并且乘子符号和大小都合理。这个验证过程不仅是在检验 fmincon更是在检验你对问题的建模是否正确。4. fmincon算法与参数调优内点法、SQP和五个必调选项4.1 Algorithm怎么选interior-point、sqp还是active-setfmincon 的 Algorithm 选项经常让人困惑因为名字里没有“拉格朗日乘子法”。实际上内点法和 SQP 都属于乘子法框架区别在于各自如何构造迭代子问题。算法适用场景特点典型坑interior-point大规模问题、不等式约束多用障碍函数把不等式约束“软化”迭代稳定可处理稀疏 Hessian初始点在可行域外时可能多花一个恢复阶段sqp中规模问题、约束梯度便宜每步解一个二次规划子问题收敛快边界跟踪直接每步代价高几十个变量时明显变慢active-set小规模问题、想观察乘子变化过程迭代路径可解释接近经典乘子法的操作方式对初始点敏感约束多时容易抖动trust-region-reflective只有边界约束或线性等式约束需要提供梯度不能直接处理一般非线性不等式约束很多人误用它去解不等式约束问题直接报错选型口诀可以这么说问题变量少、约束个数在个位、希望结果稳定可复现选 sqp问题规模大、不等式约束几十上百条、只关心最终结果选 interior-point做算法对比研究、想看清楚每一步哪个约束被激活选 active-set。这个选择可以先用你的真实数据分别跑一遍对比迭代步数和总耗时再定不丢人。4.2 五个必调参数它们各自管什么fmincon 的参数面板里有十几个选项但实际工程中真正需要动的就五个。其余保持默认往往比乱调更健康。参数默认值作用建议TolFun1e-6目标函数两次迭代的变化容差常规问题 1e-7需要高精度时 1e-9但必须配合解析梯度TolCon1e-6约束违反容差衡量“不可行度”结构安全类问题设 1e-8约束必须严格满足OptimalityTolerance1e-6KKT 一阶最优性残差容差不要小于数值差分能支持的量级有限差分时建议 1e-5 到 1e-6StepTolerance1e-10迭代步长容差判断“走不动了”一般不用动如果 exitflag2 且约束未满足说明卡在不可行点MaxIterations / MaxFunctionEvaluations与问题默认相关迭代上限和函数计算上限复杂问题常需要放大到 1000~3000但不建议盲目放大TolCon 和 OptimalityTolerance 是最容易被混淆的一对。TolCon 管的是“解离可行域多远”OptimalityTolerance 管的是“解离真最优多远”。举一个实际场景你的解约束违例只有 1e-7满足 TolCon1e-6但目标函数还有 0.5 的改善空间原因是 OptimalityTolerance 设得太大求解器提前收工了。反过来目标函数收敛得很好但约束违例 1e-5在安全评估里照样不合格。所以每次求解完这两个量都要打印出来看。4.3 读懂Displayiter的输出NLP constraint violation是硬指标设置 Displayiter 后控制台会打印一张迭代表。每一行依次是迭代步数、函数计算次数、目标函数值、NLP 约束违例量、一阶最优性度量、步长。这张表里最需要盯的是第二列往后的趋势目标函数是否在连续下降约束违例是否在朝 0 收敛。如果看到约束违例一直在 1e-2 附近抖动说明求解器正在可行域边界上“找缝”如果约束违例稳步下降但目标函数始终不动大概率是约束条件跟目标函数方向矛盾需要检查建模而不是调参数。exitflag 的速查也有用1 表示一阶最优性满足0 表示迭代步数或计算次数超限-2 表示问题不可行2 表示步长小于 StepTolerance 但约束违例不满足。看到 exitflag0 时先把 MaxIterations 提到 2000 再试一次如果还是 0问题几乎可以确定出在约束本身而不是参数上。5. 拉格朗日乘子法fmincon避坑清单五个常见翻车现场与排查5.1 现象exitflag0控制台显示迭代一步都没动我遇到过好几次这种情况目标函数写得没问题约束也没问题但 fmincon 从第一步开始就不收敛迭代次数一直跳到上限。原因多半是初始点离可行域太远。interior-point 算法有一个“恢复阶段”专门处理不可行初始点但如果约束是非线性且初始点严重违例这个阶段可能反复失败。解决思路是先做一个“可行性恢复”单独求解一个只求可行、不管目标的最优化问题以它的解作为新的初始点。% feasible_start.m function [c, ceq] feasible_confun(x) % 原约束不变 c x(1)^2 x(2)^2 - 1; ceq []; end % 目标函数临时改成“惩罚所有约束违例” obj_feas (x) sum(max(confun_constraint(x), 0).^2);更简单的方法是直接换初始点或者先用 sqp 算法跑一次它对不可行初始点的容忍度往往更高再用返回的 x 作为 interior-point 的初始点。这个技巧在工程问题里非常实用。5.2 现象开了解析梯度后结果反而变差甚至发散这是个特别容易让人怀疑人生的场景不给梯度时结果还行给了梯度反而崩了于是有人得出“解析梯度没用”的结论。真相几乎总是梯度写错。原因就出在符号约定和矩阵形状上。目标梯度写成行向量是重灾区约束雅可比的行列顺序搞反也很常见。尤其是多个不等式约束时gradc 的每一列对应一个约束而不是每一行。解决方法是跑一次 3.3 节的中心差分对比把解析梯度和数值梯度打印出来逐项看。血泪经验梯度错了fmincon 的 KKT 条件就建立在错误信息上它不但不会更快反而会朝着错误方向迭代。这不是求解器的锅是人的锅。5.3 现象手算的拉格朗日乘子与lambda.ineqnonlin对不上做完手算发现差个负号或者数值差了好几倍这是非常常见的核对事故。原因之一是符号约定。fmincon 把不等式约束写成 c(x) ≤ 0乘子取非负而你手算时可能习惯写成 x₁²x₂² ≥ 1 这种形式乘子自然反号。原因之二是尺度问题目标函数放大 10 倍乘子跟着放大 10 倍这不是 bug是数学性质。解决方法是先把所有约束统一改写成标准形式 c(x) ≤ 0再手算。数值对不上时优先检查是不是单位换算比如约束用的是兆帕还是帕其次再检查梯度方向。记住乘子不是固定的物理常数它依赖你建模时的尺度核对的唯一依据是 KKT 方程本身而不是直觉。5.4 现象约束条件在数量级上差距太大怎么调都别扭比如一个约束的数值是 1e6另一个是 1e-3目标函数又是 1e2。这种情况下 fmincon 的 KKT 系统很容易病态迭代在几个约束之间反复横跳。原因在数学上很明确KKT 系统里各项的量级差太多数值求解的舍入误差被放大。这不是 fmincon 的问题而是模型的问题。解决方法是做变量和约束的归一化。工程上常见的做法是把所有变量缩放到 0.1~10 区间约束也做无量纲化处理。比如应力约束写成“当前应力/许用应力 ≤ 1”而不是“当前应力 ≤ 许用应力”。归一化之后乘子的大小才有了跨问题比较的意义求解器的稳定性也会明显提升。5.5 现象某个约束乘子为零但你觉得它明明在起作用这可能是拉格朗日乘子法应用中最容易误判的一个坑。你看到 x* 在某条约束边界上但 lambda.ineqnonlin 对应分量是 0于是怀疑求解器出了问题。原因通常是约束冗余即当前点有多条约束同时活跃但它们对应的梯度线性相关导致 LICQ 条件不满足KKT 乘子存在但不唯一。求解器会给其中一条约束分配乘子另一条给 0。也有可能是约束在该点不可微比如绝对值约束或 max 型约束。解决方法是检查约束雅可比矩阵在解处的秩。如果秩小于活跃约束个数说明有冗余约束删掉一条再求解。对不可微约束需要做平滑近似比如用 sqrt(c²ε) 代替 |c|ε 取 1e-4 量级。乘子为零不一定说明约束没用更准确的说法是“这条约束在当前点提供的信息已经被其他约束覆盖了”。6. 用拉格朗日乘子做敏感性分析验证最优解的一把尺子6.1 用乘子预测“放松约束带来的收益”fmincon 返回的 lambda 结构体不是摆设它可以直接用来做敏感性分析。对第 3 章的示例问题不等式约束是 x₁²x₂² ≤ 1。如果我把约束右端项从 1 改成 1.001可行域稍微放大一点最优目标函数值会下降多少根据 KKT 理论目标函数对约束右端项的偏导数约等于乘子的负值。本例 λ ≈ 1.236所以可行域放大 0.001 后目标函数值大约下降 1.236×0.001 ≈ 0.00124。可以写一个脚本验证这个预测。% sensitivity.m % 把约束半径平方作为参数传入 function [c, ceq] confun_radius(x, r_sq) c x(1)^2 x(2)^2 - r_sq; ceq []; end % 用原半径求一次 [x1, fval1] run_with_radius(1.0); lambda1 get_lambda(1.0); % 用放大后的半径再求一次 delta 0.001; [x2, fval2] run_with_radius(1.0 delta); % 比较预测与实际 predicted fval1 - lambda1 * delta; fprintf(预测下降: %.6f\n, fval1 - predicted); fprintf(实际下降: %.6f\n, fval1 - fval2);如果预测值和实际值在同一个量级说明乘子不仅准而且整个 KKT 系统是自洽的。这个验证方法可以在每次求解完之后当作例行检查既能发现梯度写错、约束建模错误也能帮你理解当前设计到底被哪条约束卡住。6.2 多约束问题里如何快速判断“瓶颈约束”当问题有几十条约束时逐个去看约束违例非常低效。我的习惯是直接打印所有活跃约束的乘子按乘子绝对值排序排在前面的一定是当前优化中最“值钱”的约束。比如结构减重问题里如果应力约束的乘子远大于位移约束那说明继续减重的瓶颈在应力改材料比改截面更划算。用乘子做敏感性分析还有一个好处它把“定性判断”变成“定量判断”。不需要反复重跑优化来试只要看乘子就知道哪个约束值得去改善。这个方法我在多个项目里用过确实能省不少计算时间。早期我做优化时只盯着 x 和 fval从不看 lambda结果换一组初值就翻车后来养成了习惯每次求解完必打印 lambda再做一次约束扰动验证核完才敢拿结果去出图。希望这个习惯对你也有帮助。本文还有配套的精品资源点击获取