简介面向需要处理带约束最优化问题的数学工程与机器学习学习者这份MATLAB实践资料围绕拉格朗日乘子法原理与实现展开既讲清拉格朗日函数、拉格朗日方程以及KKT条件等核心概念也演示使用fmincon函数求解带约束非线性优化问题的完整流程。压缩包共3个文件包含2个m源码脚本和1个docx说明文档整体仅11KB轻量易用适合快速下载与本地查阅。其中m脚本为可运行的示例程序docx则重点分析不同初始点对优化结果的影响帮助读者理解乘子法求解的关键细节。已有1201人学习下载适合正在学习优化理论或需要在MATLAB中快速搭建求解代码的读者。结合原理讲解与上机实践可掌握拉格朗日乘子法从公式推导到代码落地的思路并能迁移至支持向量机、资源分配等典型场景。1. 拉格朗日乘子法不是理论玩具约束优化里它是 fmincon 的底层逻辑做结构优化的人常碰这个问题想让重量最小但应力不能超限做经济调度的人也一样想让成本最低但发电功率必须平衡。这类带约束的最小化问题背后都站着同一个数学工具——拉格朗日乘子法。而 MATLAB 里被用得最多的约束优化求解器 fmincon本质上就是在替你迭代求解拉格朗日乘子满足的 KKT 条件。很多人把 fmincon 当黑匣子用了一年却不知道它内部在算什么出了问题也不知道从哪里排查。这篇笔记把拉格朗日乘子法从手写推导讲到 fmincon 落地再把乘子迭代、参数设置和常见翻车场景串起来帮你在工程里真正敢用、会用这套东西。2. 手写拉格朗日乘子法与 KKT 条件先搞懂 fmincon 在替你算什么2.1 从等式约束出发拉格朗日函数的极值条件拉格朗日乘子法解决的是带等式约束的优化问题。问题写成min f(x)满足 h(x) 0。核心思想是把约束“吸收”进目标函数构造一个新的无约束函数L(x, λ) f(x) λᵀ h(x)这里的 λ 就是拉格朗日乘子。把约束以加权形式加到目标函数里之后原问题的最优解必然满足新函数对所有变量和乘子的偏导数为零。这个条件把约束优化问题转化成了方程组求解问题直观上也说得通在最优点处目标函数的梯度必须与约束曲面的法向平行否则沿约束曲面移动还能继续降低目标函数值。λ 就是那个度量梯度之间比例关系的系数。举个能用手算的例子。min f x₁² x₂²约束 h x₁ x₂ - 1 0。构造拉格朗日函数 L x₁² x₂² λ(x₁ x₂ - 1)对 x₁、x₂、λ 分别求偏导并令其为零∂L/∂x₁ 2x₁ λ 0∂L/∂x₂ 2x₂ λ 0∂L/∂λ x₁ x₂ - 1 0由前两个式子得到 x₁ x₂ -λ/2代入约束得 -λ 1即 λ -1因此 x₁ x₂ 0.5最优目标值 f 0.5。这个结果从几何上看也很清楚在直线 x₁ x₂ 1 上离原点最近的点就是 (0.5, 0.5)。拉格朗日乘子在这里为 -1表示约束每放松一个单位最优目标值大致下降 1 个单位这就是乘子的“对偶价格”含义后面验证 fmincon 结果时还会用到。2.2 不等式约束与 KKT 条件互补松弛的直观理解实际工程问题里约束大多是不等式比如应力不能超过许用值、流量不能低于下限。问题变成 min f(x)满足 g(x) ≤ 0。这时单靠拉格朗日函数求偏导不够需要补充一套条件就是 KKT 条件。KKT 条件包含四条原问题可行约束必须满足、梯度平衡拉格朗日函数对 x 的偏导为零、对偶可行不等式对应的乘子 λ ≥ 0、互补松弛λᵢgᵢ(x) 0。互补松弛是最容易忽略的一条。它说的是乘子与对应的约束函数至少有一个为零约束不起作用时乘子为零乘子非零时约束必然处于边界。这个性质让拉格朗日乘子法能自动识别哪些约束“紧”、哪些约束“松”。比如一个储层压力约束如果离边界很远它的乘子逼近 0求解器就知道这个约束暂时不用管一旦迭代接近边界乘子开始变大把解“推”回可行域内。fmincon 的 interior-point 算法内部就是在解 KKT 系统的变形。它把不等式约束通过障碍函数对数屏障转成一系列等式约束的子问题再对每个子问题求解拉格朗日函数的驻点。所以你在 fmincon 里设置约束时它其实是在求解一个包含乘子的方程组而不是简单地“把不满足约束的解剔除”。理解这层含义后面看 lambda 输出、调容差参数时就不会一头雾水。2.3 用符号计算验证手写 KKT一段能跑的 MATLAB 代码理论推完直接用 MATLAB 符号工具箱验证 2.1 的例子。这段代码把拉格朗日函数写出来求梯度并解方程组和手算结果对照。syms x1 x2 lambda f x1^2 x2^2; h x1 x2 - 1; L f lambda * h; gradL gradient(L, [x1, x2, lambda]); sol solve(gradL(1) 0, gradL(2) 0, gradL(3) 0, ... [x1, x2, lambda], Real, true); disp(double([sol.x1, sol.x2, sol.lambda]));逻辑说明gradient 函数对符号表达式求梯度得到三个方程分别对应 ∂L/∂x₁、∂L/∂x₂、∂L/∂λsolve 用符号方式解方程组Real, true 只保留实数解避免复数干扰。输出结果应为 0.5、0.5、-1与手算一致。参数说明syms 声明符号变量时lambda 在 MATLAB 里不是保留字可以放心用gradient 的第二个参数传变量向量返回的 gradL 是列向量顺序与变量顺序一致。如果约束变成非线性比如 h x₁² x₂² - 1solve 可能解不出显式根这时换成 vpasolve 做数值求解更稳妥。syms x1 x2 lambda f (x1 - 1)^2 (x2 - 2)^2; h x1^2 x2^2 - 1; L f lambda * h; gradL gradient(L, [x1, x2, lambda]); sol vpasolve(gradL(1) 0, gradL(2) 0, gradL(3) 0, ... [x1, x2, lambda]); disp(double([sol.x1, sol.x2, sol.lambda]));vpasolve 用数值方法在默认搜索范围内找一组解。注意非线性方程的驻点可能不止一个vpasolve 只返回一个解要穷举全部驻点需要带初值反复调用或者画等值线看。这一步也是后面用 fmincon 时的参照系fmincon 同样只能保证找到局部最优解多个驻点时结果依赖初始点。3. 用 fmincon 跑通带约束最小化从目标函数到求解器的完整落地3.1 最小可运行案例匿名函数还是函数文件fmincon 的标准调用形式是 [x, fval, exitflag, output, lambda] fmincon(fun, x0, A, b, Aeq, beq, lb, ub, nonlcon, options)。fun 是目标函数可以用匿名函数写在脚本里也可以单独写成 function 文件。二者选型原则很简单目标函数只有一两行表达式时用匿名函数省去建文件的开销目标函数长、需要调用外部数据或者要提供梯度时用函数文件更清晰。拿 2.1 的例子做个最小实现。约束只有等式 h x₁ x₂ - 1写成线性形式 Aeq·x beq其中 Aeq [1, 1]beq 1。目标函数用匿名函数完整代码如下。fun (x) x(1)^2 x(2)^2; x0 [0, 0]; Aeq [1, 1]; beq 1; [x, fval, exitflag, output, lambda] fmincon(fun, x0, [], [], Aeq, beq); disp(x); disp(fval); disp(lambda.eqlin);运行结果x [0.5000, 0.5000]fval 0.5000lambda.eqlin -1.0000。这个输出和手写拉格朗日乘子法得到的乘子完全一致。注意 A、b、lb、ub 这些不需要的输入用 [] 占位这是 fmincon 调用里最容易写错的地方——占位符少一个后面所有参数整体错位报错还往往指向奇怪的维度问题。3.2 fmincon 六个输入分组的含义线性约束、边界与非线性约束fmincon 把约束分成四组每组在内部走不同的处理路径。A、b 描述线性不等式 A·x ≤ bAeq、beq 描述线性等式 Aeq·x beqlb、ub 是变量硬边界nonlcon 描述非线性不等式 c(x) ≤ 0 和非线性等式 ceq(x) 0。分类的意义在于线性约束可以提前约简减少求解器变量维度边界约束在算法里处理方式与一般约束不同interior-point 对 lb、ub 的障碍函数处理效率更高非线性约束则每轮迭代都要重新计算函数值和梯度最贵。实际工程里常见的错误是把线性约束写进 nonlcon。比如变量 x₁ x₂ ≤ 5 明明可以写成 A [1, 1], b 5有人却写成非线性约束函数。后果是每次迭代都多付一次函数求值的开销而且数值梯度计算误差会让收敛变慢。判断标准很简单约束里有没有变量之间的非线性运算乘积、指数、三角函数没有就是线性约束必须放进 A、b 或 Aeq、beq。非线性约束的例子给 2.1 的例子里加一个不等式 x₁·x₂ ≥ 2这是一个非线性不等式写成标准形式 c(x) ≤ 0 就是 2 - x₁·x₂ ≤ 0。function [c, ceq] mycon(x) c 2 - x(1) * x(2); ceq []; end主程序里通过 nonlcon 参数传入 mycon。注意函数返回两个输出c 是行向量形式的非线性不等式组ceq 是等式组没有等式约束时 ceq 必须返回空数组不能省略第二个输出。省略了 MATLAB 会报错提示输出参数个数不足。3.3 必调的四个参数Algorithm、Tolerance、MaxIterations、Displayoptions 用 optimoptions 创建。这里给出最常调整的四个参数及典型取值。参数典型取值作用Algorithminterior-point / sqp / active-set决定 KKT 系统的求解策略ConstraintTolerance1e-6默认约束满足程度影响可行性OptimalityTolerance1e-6默认KKT 条件的梯度残差容忍度MaxIterations400默认防死循环非凸问题加大Displayiter / final控制命令行输出粒度Algorithm 的选择直接影响拉格朗日乘子法的落地方式。interior-point 适合大规模问题、约束较多的情况默认推荐能处理稀疏大规模矩阵sqp 适合中小规模、约束函数计算昂贵的情况每步收敛扎实但内存消耗比 interior-point 大active-set 是早期算法对初值敏感新代码不建议用。碰到非光滑目标函数时可以用 sqp 替代 interior-point 试一下有时会显著改善。Tolerance 参数是避坑关键。ConstraintTolerance 设太严比如 1e-10求解器可能反复迭代都无法满足约束最后报迭代终止但约束不满足OptimalityTolerance 设太松解停在远离真正极值点的地方。工程上一般先保持默认 1e-6计算完成后看 exitflag 和 output.constrviolation 判断是否需要调整。options optimoptions(fmincon, ... Algorithm, sqp, ... ConstraintTolerance, 1e-6, ... OptimalityTolerance, 1e-6, ... MaxIterations, 800, ... Display, iter); [x, fval, exitflag, output, lambda] fmincon(fun, x0, [], [], Aeq, beq, ... [], [], [], options);参数说明optimoptions 的第一个参数必须写 fmincon 指定求解器类型因为不同求解器支持的选项不一样写错会直接报错。Display 设为 iter 可以看到每轮迭代的目标函数下降量和约束违反度排查收敛异常时非常有用。MaxIterations 是非凸工程问题里最常调大的参数默认 400 对简单问题足够但对强非线性约束的问题经常不够用。4. 乘子法家族与罚函数什么时候该绕过 fmincon 自己写迭代4.1 罚函数为什么会被淘汰病态问题的来源在增广拉格朗日法普及之前外点罚函数法是处理约束的常见手段。做法是把问题改写成 min f(x) ρ·h(x)²ρ 是一个很大的正数逼迫解靠近可行域。看起来简单但 ρ 一大Hessian 矩阵条件数随 ρ 增长梯度下降和拟牛顿法的收敛速度急剧恶化。也就是说为了让约束严格满足而把 ρ 调大反而让无约束优化的子问题变得病态迭代步数暴涨甚至发散。用一个数值例子说明。min f x₁² x₂²约束 h x₁ x₂ - 1 0罚函数写成 Φ x₁² x₂² ρ(x₁ x₂ - 1)²。解析求 ∇Φ 0得到 x₁ x₂ ρ/(1 2ρ)ρ → ∞ 时逼近 0.5但 ρ 1000 时计算出的数值解和真实解仍有一段距离。而且沿最速下降方向Hessian 的特征值分别是 2 和 2 2ρ条件数约 1 ρρ 1e6 时条件数 1e6标准 BFGS 直接失去精度。这就是罚函数的本质缺陷精确性靠大 ρ 换数值稳定性随 ρ 恶化。4.2 增广拉格朗日乘子法罚项与乘子迭代的配合增广拉格朗日乘子法解决了罚函数的病态问题。它在拉格朗日函数后面加一个二次罚项L_A(x, λ, ρ) f(x) λᵀh(x) (ρ/2)‖h(x)‖²关键变化是求解 L_A 的驻点后用公式 λ ← λ ρ·h(x) 更新乘子再进行下一轮。直观理解二次罚项先帮解靠近可行域乘子项再逐步修正让最终解精确落在约束曲面上。ρ 不需要无限增大保持适中的数值就能达到精度Hessian 条件数被控制住子问题求解稳定得多。乘子更新公式的推导很直接。设第 k 轮解为 x_kKKT 条件要求 ∇f(x_k) λ∇h(x_k) 0。对 L_A 求梯度令其为零∇f(x_k) λ∇h(x_k) ρh(x_k)∇h(x_k) 0。对比两式λ 的修正量就是 ρ·h(x_k)所以 λ_{k1} λ_k ρ·h(x_k)。这是一个一阶迭代收敛速率取决于 ρ 的选择和问题的非线性程度。4.3 手写乘子法循环收敛判断与参数选择用增广拉格朗日乘子法解 2.1 的例子代码展示完整循环。内层用 fminunc 求解无约束子问题外层更新乘子直到约束违反度小于容差。% 增广拉格朗日乘子法求解 min x1^2x2^2 s.t. x1x2-10 % 内层用 fminunc外层更新乘子 f (x) x(1)^2 x(2)^2; h (x) x(1) x(2) - 1; rho 1; % 罚参数初始值 lambda 0; % 拉格朗日乘子初始值 x [0, 0]; % 变量初始值 tol 1e-8; % 约束违反度容差 for k 1:50 % 构造增广拉格朗日函数作为内层目标 LA (xx) f(xx) lambda * h(xx) (rho/2) * h(xx)^2; optionsIn optimoptions(fminunc, Display, off, Algorithm, quasi-newton); x fminunc(LA, x, optionsIn); % 检查约束违反度决定退出还是继续 hval h(x); if abs(hval) tol fprintf(k%d, x(%.8f, %.8f), lambda%.8f, h%.2e\n, k, x(1), x(2), lambda, hval); break; end % 乘子更新 lambda lambda rho * hval; fprintf(k%d, x(%.8f, %.8f), lambda%.8f, h%.2e\n, k, x(1), x(2), lambda, hval); end逻辑说明内层 fminunc 求的是当前 λ 和 ρ 下的增广拉格朗日函数驻点外层用 hval 判断约束是否满足不满足就按 λ ← λ ρ·h 更新乘子。由于问题简单通常在 20 轮以内收敛x 趋近 (0.5, 0.5)λ 趋近 -1与 fmincon 的 lambda.eqlin 输出一致。参数说明ρ 初始值选 1更新过程中不改变它实际工程里 ρ 可以在每轮乘子更新不明显时按 ρ ← min(2ρ, ρ_max) 增长ρ_max 一般取 1e6。λ 初始值取 0 即可取非零初值可以加速但需要问题背景支撑。退出条件用约束违反度这是增广拉格朗日法的天然停机准则比看相邻两轮 x 的差值更可靠。什么时候值得绕过 fmincon 自己写乘子法迭代两种情况。第一种你的问题里乘子有明确的物理意义比如经济调度中的电价影子价格、结构优化中的应力约束乘子代表灵敏度你希望精确控制乘子的收敛路径。第二种fmincon 对你的问题反复失败而你能提供一个好的 λ 初值和 ρ 更新策略。除此之外fmincon 内部封装的增广拉格朗日变体通常比自己写的实现更鲁棒不要为了“用乘子法”而用乘子法。5. 拉格朗日乘子法避坑指南五个让求解器翻车的常见问题5.1 现象fmincon 返回 complex value目标函数算不下去现象是 exitflag 为负fval 返回复数或者报错“Objective function returned a complex value”。原因通常是目标函数或约束里带 sqrt、log、x^0.5 这类运算而迭代过程中变量被试探到负域。解决分两层能用边界约束挡住的在 lb、ub 里限制变量为正挡不住的在目标函数里加保护逻辑例如用 x.^2 替代 sqrt 的平方形式。这类问题在参数估计和力学计算里很常见尤其分母上带变量的表达式极易在乘子迭代的试探步里踩到奇异点。5.2 现象不同初始点得到不同结果乘子值也跟着变现象是同一个约束优化问题换 x0 后最优解变了lambda 输出也变。原因是非凸问题的 KKT 条件存在多个局部解interior-point 和 sqp 都是局部优化算法只能保证收敛到某个驻点。解决手段有两种用 MultiStart 配合 fmincon 做多初值扫描或者用 GlobalSearch对乘子比较敏感的问题可以从乘子的物理含义出发选初值——比如知道乘子大概量级把 x0 取在让它接近约束边界的位置上。多起点扫描的代价是计算量翻倍但换来对“局部最优还是全局最优”的判断力工程上值得。5.3 现象等式约束始终满足不了报“Equation solved but... ”现象是 exitflag 为 1 但 output.constrviolation 明显大于 0或者提示等式约束不满足。原因多为 ConstraintTolerance 设得比 1e-8 还严而等式约束本身是强非线性、初始点离约束曲面太远迭代步长不足以精确着陆。解决路径先把 ConstraintTolerance 放宽到 1e-6 或 1e-5 看可行性是否恢复如果还不行给 nonlcon 的等式约束做一个比例缩放把 ceq 除以一个特征尺度让约束函数值在 1 的量级避免数值消去。另一个隐蔽坑等式约束最好写成两个不等式形式 c ≤ 0 和 -c ≤ 0 吗不要这样做这会破坏约束雅可比结构直接用 ceq 返回最稳。5.4 现象手写增广拉格朗日循环发散或震荡现象是自己写的乘子迭代里h 值先小后大或者 λ 在某个值附近来回跳。原因多是 ρ 更新太快乘子修正步过长导致内层 fminunc 的子问题每次都从很差的初值启动另一个可能是内层容差太松子问题没收敛就更新 λ误差一路累积。解决ρ 增长率控制在每轮 1.5 到 2 倍之间不要超过 10内层 fminunc 的 OptimalityTolerance 设为 1e-10 或更严保证子问题的解足够准每轮输出 h 和 λ观察震荡模式再决定调整方向。5.5 现象lambda 输出符号与手写 KKT 对不上现象是手算拉格朗日乘子得到 λ -1fmincon 的 lambda.eqlin 恰好也是 -1但换成不等式约束时符号总和自己推的差一个负号。原因在于 fmincon 的乘子符号约定和常见教科书写法不同。fmincon 文档里对不等式约束 g(x) ≤ 0返回的 lambda.ineqnonlin ≥ 0而对等式约束lambda.eqlin 可正可负。自己写 KKT 条件时如果用 h(x) 0 作为等式乘子符号取决于梯度平衡方程的写法。解决统一以 fmincon 的输出为基准手写验证时把符号约定代入别花时间质疑求解器。具体验证方法是把 KKT 残差打印出来见下一章。6. 用拉格朗日乘子法验证 fmincon 结果梯度检验与对偶间隙的实操技巧fmincon 返回的 lambda 结构体不是摆设它是验证解是否真正满足 KKT 条件的钥匙。拿到解 x 和乘子后手工计算目标函数梯度和约束梯度组装 KKT 残差∇f(x) Σλᵢ∇gᵢ(x) Σμⱼ∇hⱼ(x) 的范数应当接近求解器的 OptimalityTolerance。工程上我习惯用有限差分梯度对比 fmincon 输出的 grad 字段能快速发现目标函数写错、变量顺序不一致这类低级错误。一个实用技巧开启 SpecifyObjectiveGradient 和 SpecifyConstraintGradient 选项把自己的解析梯度传给 fmincon同时把 CheckGradients 设为 true让求解器在第一步帮你校验梯度。这个选项平时关闭但每当更换目标函数或重构约束时跑一次能省掉大量排查时间。更隐蔽的做法是观察对偶间隙如果问题是凸的fmincon 求解结束时 fval 与原问题对偶问题的目标值应当几乎相等。对偶间隙大就意味着原始问题非凸或者数值精度不够此时乘子值只能当参考不能直接物理解读。最后说个血泪经验不要一上来就信 lambda 数值。任何自动微分和有限差分梯度都有误差乘子的灵敏度分析结论要配合约束边界的小扰动验证。我一般会在最优解旁边做一次 pull 测试——把约束边界人为偏移 1%重算最优值用差分近似乘子的物理解释与 lambda 对比。这套习惯坚持下来拉格朗日乘子法就从一个理论名词变成了手上最趁手的工程工具希望帮到你。本文还有配套的精品资源点击获取
