简介外点法是求解约束非线性规划问题的常用优化算法尤其适合处理可行域复杂、难以直接投影的带约束场景。这份Matlab程序实例面向学习运筹优化、人工智能算法以及需要解决约束建模问题的研究人员与工程师演示了如何构造惩罚函数将约束并入目标通过逐步调整惩罚因子逼近最优解。压缩包共9个文件全部为.m脚本包括目标函数、等式与不等式约束、梯度与雅可比矩阵计算、newton法辅助模块以及用于统一调度的主程序各文件职责划分清楚覆盖从目标函数定义、约束构建、惩罚函数生成到迭代求解与停止判断的关键环节可对照源码逐行推演外点法的完整计算流程。压缩包整体仅3KB内容精简易读几乎没有冗余文件。目前已有674人学习下载。通过学习该实例读者能够掌握惩罚因子的选取与调整策略、约束违反程度的判断方法以及调用fmincon等优化工具箱函数实现求解的具体写法同时还能理解基于梯度的搜索方向更新与收敛停止准则并可将外点法迁移到其他工程优化问题中是一份简洁且具有较强参考价值的算法实现代码。1. 外点法matlab程序实例为什么值得从零手写手写一个外点法外点罚函数法求解器是理解约束优化最直接的路径。与内点法始终把迭代点锁在可行域内不同外点法从不可行点出发靠不断增大的罚因子把解“拉”进可行域。这个“先违反、后收敛”的思路让它可以处理那些无法轻易给出可行初始点的工程问题也很适合与 matlab 优化工具箱配合主循环自己写内层无约束求解交给 fminunc或完全用梯度下降实现。对于正在学优化方法的学生以及需要自定义约束逻辑的工程师来说这套“外层增罚、内层无约束优化”的程序实例都能在半小时内改成可调、可验证的最小实现。什么情况下值得放弃 fmincon 转手写外点法一是目标函数和约束规模小、但结构特殊需要精确控制惩罚方式二是想观察罚因子变化如何影响解的轨迹三是手头没有优化工具箱只能靠基础 MATLAB 函数完成求解。下面从罚函数的数学形式切入给出两版可运行代码最后落到多约束与乘子法改进。2. 外点罚函数法的数学模型与收敛性前提2.1 罚函数怎么写从不等式到 max 算子约束优化问题的一般形式是$$\min_{x\in\mathbb{R}^n} f(x), \quad \text{s.t.}\quad g_i(x)\le 0\ (i1,\dots,m),\quad h_j(x)0\ (j1,\dots,l)$$外点法把约束“吸收”进目标函数构造辅助函数$$P(x,M)f(x)M\sum_{i1}^{m}\left[\max\left(0,\ g_i(x)\right)\right]^{2}M\sum_{j1}^{l}\left[h_j(x)\right]^{2}$$其中 M 为罚因子。对不等式约束max(0, g_i(x)) 保证了只有违反约束时才产生惩罚g_i(x)≤0 时该项恒为 0g_i(x)0 时按违反程度的平方放大。等式约束则始终以平方形式参与惩罚无论从哪一侧偏离 h_j(x)0都会被罚项拉回。选择二次而不是一次惩罚核心原因是可微性。只要 f、g、h 自身足够光滑平方后的罚函数对 x 连续可微梯度与 Hessian 都存在fminunc 或自写梯度法都能直接使用。代价是可行域边界上 max 算子引入一阶不可导点实际计算中拟牛顿法通常能容忍如果追求更强光滑性可以改用三次或更高次惩罚但问题病态程度会同步上升。2.2 罚因子序列为什么必须趋于正无穷固定某个有限 M 时P(x, M) 的最优点 x*(M) 通常并不满足原始约束。目标函数会尽量把解往无约束最优方向推惩罚项则把解往可行域拉双方在约束边界外侧达成平衡。要让 x*(M) 真正回到可行域就必须让惩罚项的相对权重压倒目标函数也就是令$$M_0 M_1 M_2 \dots \to \infty$$收敛性理论上要求三点目标函数与约束函数在最优解邻域内连续罚因子序列严格递增且发散每一轮内层无约束优化都要解得足够精确否则误差会跨轮累积。这些条件满足时x*(M_k) 的任意极限点都是原问题的局部最优解加上凸性假设后可以升级为全局最优解。需要注意“罚因子越大越好”是个常见误解。M 过大会让 P 的 Hessian 变得病态条件数随 M 线性增长内层求解器反而更难收敛。实际工程里通常取 M01~10按 5~10 倍递增而不是一步推到 1e6。2.3 外点与内点的本质区别迭代点的可行性外点法得名于迭代点始终落在可行域之外靠惩罚逐步逼近边界。与之相对的内点法屏障法通过在目标函数里加入屏障项把迭代点限制在可行域内部。两者在使用体验上的差异非常明显对比项外点法内点法/屏障法初始点要求任意点不需要可行必须是严格可行内点迭代点位置从不可行逐渐逼近边界始终在可行域内部约束类型适配等式与不等式都很自然不等式方便等式需额外处理中途“解”的可用性不可行不能直接作为工程解近似可行可提前停机实现成本几十行 MATLAB 代码足够需处理屏障参数与边界有界性因为外点法几乎不挑剔初始点工程上很多场景都在用它。比如模型预测控制里每一步热启动点来自上一时刻的次优解不保证可行外点法可以直接接管。2.4 一个可以解析验证的标准算例选一个能用手算验证的例子后续所有程序都以它为准$$\min\ f(x)x_1^{2}x_2^{2} \quad \text{s.t.}\quad g(x)2-x_1-x_2\le 0$$直观上最优解是直线 x1x22 上离原点最近的点即 x*(1, 1)目标值 f*2。对应的罚函数为$$P(x,M)x_1^{2}x_2^{2}M(2-x_1-x_2)^{2}$$注意第二项只在 x1x22 时非零。对 x1、x2 求偏导并令其为零得到 x1x22M/(12M)。把不同 M 代入可以得到一条完全确定的收敛轨迹Mx1x2f(x)约束 g(x)00.00000.00002.000010.66670.88890.6667100.95241.81410.09521000.99501.98010.00991e40.999951.999900.000101e60.99999951.9999990.000001罚因子每增加两个数量级约束违背大约缩小两个数量级目标函数从下方逼近 2。这条解析轨迹就是后面程序的“标准答案”数值结果落在同一趋势上就说明外循环和内层求解器都没写错。3. 外点法MATLAB程序实例调用fminunc的内外层循环3.1 程序骨架与职责划分一个完整的外点法程序建议按三层职责分开写问题定义层目标函数、不等式约束、等式约束的函数句柄外循环层负责更新罚因子 M、调用内层求解、检查停机条件内层求解层把当前 M 下的 P(x, M) 当作无约束优化问题处理。这种划分的最大好处是换问题时不用动外循环只改第一层换内层求解器时不用碰问题定义。在 MATLAB 里用函数句柄加元胞数组就能把多约束问题描述得足够干净。下面的实例对应 2.4 节标准算例文件名可命名为 outerpoint_demo.m与 outerpointmethod 这个主题词的检索习惯保持一致。3.2 可直接运行的代码% outerpoint_demo.m % 外点法求解: % min f(x) x1^2 x2^2 % s.t. g(x) 2 - x1 - x2 0 clear; clc; % --- 问题定义 --- f_obj (x) x(1)^2 x(2)^2; g_ineq (x) 2 - x(1) - x(2); % g 0 % 罚函数: P f M * max(0, g)^2 penalty (x, M) f_obj(x) M * max(0, g_ineq(x))^2; % --- 外层参数 --- M0 1; M M0; % 罚因子初值 gamma 10; % 罚因子增长倍数 tol_v 1e-6; % 约束违背容差 tol_x 1e-4; % 相邻代理解距离容差 max_outer 50; % 外层最大迭代数 x_prev [0; 0]; % 初始点故意不可行 viol inf; for k 1:max_outer % 内层无约束优化: fminunc 求解 P(x, M) opts optimoptions(fminunc, ... Algorithm, quasi-newton, ... Display, off, ... MaxIterations, 500, ... MaxFunctionEvaluations, 2000, ... StepTolerance, 1e-10, ... OptimalityTolerance, 1e-8); [x_opt, ~, exitflag] fminunc((x) penalty(x, M), x_prev, opts); viol max(0, g_ineq(x_opt)); fprintf(k%2d M%7.1e x(%9.6f, %9.6f) f%.8f viol%.2e\n, ... k, M, x_opt(1), x_opt(2), f_obj(x_opt), viol); % 停机: 约束已满足且前后两次解基本不动 if viol tol_v norm(x_opt - x_prev) tol_x break; end x_prev x_opt; % 热启动 M M * gamma; % 惩罚加重 end fprintf(结果: x*(%.8f, %.8f), f*%.8f, g%.2e\n, ... x_opt(1), x_opt(2), f_obj(x_opt), viol);这段代码在笔者的环境中跑出的轨迹与 2.4 节理论表一致大约 12~18 次外循环后满足停机条件具体次数会随 MATLAB 版本和机器精度略有浮动。exitflag 为 1 表示内层无约束优化正常收敛若为 0说明迭代次数或函数评估次数触顶需要放大对应的 Max 选项而不要怀疑外点法本身的逻辑。3.3 罚因子初值、增长倍数、容差三个参数的设置逻辑外点法参数不多但每个都直接影响迭代轨迹参数建议取值设置逻辑M01 ~ 10太小则前几轮约束几乎不受惩罚白做无约束优化太大则初始问题就病态gamma5 ~ 10增长太慢外循环多太快则相邻两轮问题差异大热启动失效tol_v1e-5 ~ 1e-7决定最终可行性应比约束尺度本身小两个数量级以上tol_x1e-4 左右防止 M 已经很大但解还在缓慢移动时无限循环M 的增长倍数 gamma 值得单独说明。若 gamma 过小比如 1.1外循环需要几十次才能把 M 从 1 推到 100而相邻两轮的最优解变化极微小整体上等于在做大量冗余的无约束优化。若 gamma 过大比如 100M 从 1 跳到 100 后惩罚项突然主导上一轮的最优点作为本轮初值离新问题的最优点太远fminunc 需要额外迭代去追赶外循环次数少了但总成本没降。5~10 是实践中最稳健的区间。3.4 热启动策略为什么每次内层迭代都用上次的最优解外点法每一轮对应不同的 M下一轮的无约束问题与上一轮相比只是惩罚权重变化最优解只移动一小段。如果把初始点固定不变每轮都从零开始优化前几轮还能勉强运行M 变大后每一步都要重新跨越整个解空间计算量成倍上涨。把上一轮的 x_opt 当作本轮初值就是热启动。它对内层算法的迭代次数影响非常大尤其 M 变大后P 的等高面被拉长成狭长形状从远处冷启动很容易在陡峭的惩罚槽里来回震荡。热启动配合拟牛顿法时fminunc 还能复用近似的 Hessian 信息收敛速度会明显优于冷启动。4. 不依赖工具箱的梯度下降实现与收敛性对比4.1 自写最速下降配回溯线搜索的梯度框架如果手头机器没有 matlab 优化工具箱许可证fminunc 不可用。只要罚函数可微自己写一个最速下降法配合回溯线搜索就足够支撑外点法的内层求解。核心函数如下function [x, iter] steepest_descent(fun, grad, x0, max_iter, tol) % fun : 目标函数句柄 % grad : 梯度函数句柄 % x0 : 起始列向量 % max_iter : 最大迭代数 % tol : 梯度无穷范数停机阈值 alpha0 1.0; % 初始步长 rho 0.5; % 步长衰减率 c1 1e-4; % Armijo 充分下降系数 x x0(:); for iter 1:max_iter g grad(x); if norm(g, inf) tol return; end d -g; % 最速下降方向 alpha alpha0; f_now fun(x); % 回溯: 只要不满足 Armijo 条件就缩小步长 while fun(x alpha * d) f_now c1 * alpha * (g * d) alpha rho * alpha; end x x alpha * d; end end这里使用无穷范数停机是因为罚函数在边界附近的梯度分量可能一个极大、另一个极小只用欧氏范数容易误判。Armijo 条件中的 c1 取 1e-4 是数值优化里的经典配置rho 取 0.5 让步长搜索速度与稳定性保持平衡alpha0 取 1.0 则是因为最速下降方向配合充分下降条件时单位步长在罚函数光滑区域通常是可行起点。4.2 解析梯度与数值梯度的验证标准算例下罚函数为 Px1^2x2^2M(2-x1-x2)^2解析梯度为 ∇P(2x1-2M(2-x1-x2), 2x2-2M(2-x1-x2))。用中心差分验证解析式是否正确M 10; x [0.3; 0.8]; g_analytic 2*x - 2*M*(2 - x(1) - x(2))*[1; 1]; % 注意 g0 时才成立 g_numeric zeros(2,1); eps_g 1e-6; for i 1:2 e zeros(2,1); e(i) eps_g; g_numeric(i) (penalty(xe,M) - penalty(x-e,M)) / (2*eps_g); end disp([g_analytic g_numeric]);这一小段验证值得每次改完罚函数都跑一次。解析梯度写错时数值梯度对比能直接指出差在哪个分量。对含 max 的罚函数要特别注意x 恰在约束边界 g(x)0 时max 算子不可导解析梯度与数值梯度都会有偏差且偏差比远离边界处略大属于正常现象。4.3 不同罚因子增长倍数的收敛行为对比把上面的最速下降嵌入外点法在相同容差下改变 gamma可以得到如下趋势数值在同一台机器上观察得到绝对值随版本浮动相对趋势稳定gamma外循环次数总梯度调用约最终违约束量240 以上800 以上1e-6 附近520 左右400 左右1e-61015~18350 左右1e-65010 左右1000 以上1e-6gamma 从 2 提到 10总成本明显下降gamma 提到 50 以后每轮无约束优化难度陡增总梯度调用次数反而上升。这说明参数并非越大越好也印证了前面章节对 gamma 的分析。真实工程中可以先取 gamma10 跑一轮再根据每轮内层迭代是否触顶来微调。4.4 排错内层求解不收敛时先查什么最常遇到的现象有三个fminunc 返回 exitflag0或最优性残差长时间不下降先看 MaxFunctionEvaluations 是否触顶。罚函数含 max 时边界非光滑点会让拟牛顿法反复试探函数评估次数消耗很快。M 增大后回溯线搜索把步长压到极小这是最速下降在病态二次函数上的典型表现。换成 BFGS 或共轭梯度或者把 gamma 调回 5。终解的约束违背比 tol_v 大一个量级多数情况是外层循环提前被 tol_x 停机。应检查 tol_x 是否设得比 tol_v 更严格或直接改成仅靠 viol 与 M 上限作为停机条件。5. 多约束与等式约束的扩展用乘子法弥补外点法的不足5.1 把问题定义改成约束元胞数组实际工程很少只有单个不等式约束。把约束句柄收进元胞数组罚函数模块不需要改结构g_cell { (x) -x(1), ... % x1 0 写成 -x1 0 (x) x(2) - x(1)^2 }; % x2 - x1^2 0 h_cell { (x) x(1) x(2) - 1 }; % 等式约束 function P penalty_all(x, M, f_obj, g_cell, h_cell) P f_obj(x); for i 1:numel(g_cell) P P M * max(0, g_cell{i}(x))^2; % 逐项独立惩罚 end for j 1:numel(h_cell) P P M * h_cell{j}(x)^2; end end逐项独立计算惩罚不要先对所有约束求和再取 max否则某些约束被满足时会被另一些违反的约束“带伤”。这种写法也方便在日志里单独记录每个约束的违背量。5.2 乘子法改进的更新式外点法要把罚因子推到很大才能满足可行性而乘子法在罚项中加入拉格朗日乘子的线性项让有限罚因子也能达到高精度。不等式约束的增广拉格朗日函数常用形式为$$L(x,\lambda,\mu)f(x)\frac{1}{2\mu}\sum_{i}\left[\max\left(0,\ \lambda_i\mu g_i(x)\right)\right]^2-\frac{1}{2\mu}\sum_{i}\lambda_i^2$$乘子更新式为$$\lambda_i \leftarrow \max\left(0,\ \lambda_i\mu g_i(x)\right)$$实现时只需在外循环末尾更新乘子罚因子按 2~5 倍缓慢增长收敛速度通常比纯外点法快一到两个数量级。与第 3 章程序相比改动量只有几行非常划算。5.3 用参考解验证扩展程序把 5.1 的约束与目标 f(x1-2)^2(x2-1)^2 组合可行域由 x1x21 与 x2≤x1^2 相交而成该问题没有直观的闭式解。验证方法是先用 fmincon 的 sqp 算法跑出参考解再让外点法程序从同一初始点出发两步都收敛后比较目标值与约束违背量差异在 1e-5 以内即视为实现正确。还可以把每轮 x_opt 和 viol 存入数组绘制 viol 随 M 变化的 semilogy 曲线观察曲线是否呈线性下降趋势。此时再回头对照 2.4 节的解析轨迹外点法、增广拉格朗日与 fmincon 三条路径的目标值应一致到 1e-5 以内构成完整的交叉验证闭环。本文还有配套的精品资源点击获取