简介压缩包内为基于Matlab实现的外点法程序实例通过源代码直观演示外点法求解含约束非线性规划问题的完整流程面向相关课程学生、科研人员与优化算法初学者也适合中高级研究者快速复现算法。资源共9个文件全部为M脚本体积仅3KB脚本按模块拆分分别承担目标函数定义、等式与不等式约束、梯度与雅可比矩阵计算、牛顿迭代求解等任务结构清晰便于逐模块对照运行和修改。目前已有674人学习/下载。外点法的核心思想是通过惩罚因子将约束纳入目标函数从不可行域逐步逼近最优解这套程序正好演示了惩罚函数如何随迭代调整、收敛阈值如何判定以及fmincon等优化工具箱的调用流程。在理解代码后替换目标函数与约束条件即可将外点法迁移到工程优化、路径规划或人工智能模型调参等实际场景中。1. 外点法把约束揉进目标函数的 MATLAB 入门路线约束优化里有个反直觉的结论最优解往往不在可行域内部舒服地待着而是被某个不等式或者等式约束压到边界上。处理这类问题最直白的手段不是给 MATLAB 装一堆求解器而是把约束“罚”进目标函数让迭代点从禁止区外面一点点靠近边界——这就是外点法也叫外部罚函数法。它思路短、代码省适合快速验证约束优化的想法也是理解增广拉格朗日乘子法的必经一站。这篇博文就从罚函数构造讲起给出一套不依赖优化工具箱也能跑通的外点法 MATLAB 程序实例把惩罚因子、约束违约度、内层无约束优化这几个关键旋钮逐一拧开。2. 外点法数学模型与惩罚函数的构造2.1 罚函数从哪来一次罚、二次罚与 max(0, ·)²考虑标准约束优化问题min f(x) s.t. g_i(x) ≤ 0, i 1,...,m h_j(x) 0, j 1,...,n外点法的核心是把上述问题改写成无约束优化构造增广目标函数P(x, r) f(x) r/2 * ( Σ max(0, g_i(x))² Σ h_j(x)² )其中r 0是惩罚因子。注意这里对不等式约束只惩罚“违反”的部分当g_i(x) ≤ 0时罚项为 0一旦越界max(0, g_i(x))²开始贡献正的惩罚。等式约束则不管正负h_j(x)²一律计入。为什么用二次罚而不是一次罚一个直接原因max(0, g)²在g 0处导数存在且连续左导数和右导数都为 0梯度表达式干净而一次罚max(0, g)在交界处不可微内层无约束优化会在边界附近出现锯齿状的震荡。系数写成r/2也纯粹是约定这样求导后系数正好抵消梯度表达式更清爽。% 罚函数 P 和它的梯度 g (x) 3 - x(1) - x(2); grad_g (x) [-1; -1]; h (x) x(1) - x(2); grad_h (x) [1; -1]; P (x, r) f(x) r/2 * ( max(0, g(x)).^2 h(x).^2 ); gradP (x, r) grad_f(x) r * ( max(0, g(x)) * grad_g(x) h(x) * grad_h(x) );这段代码先定义了约束函数和梯度再拼出罚函数。max(0, g(x))返回的是违反量本身如果g(x)已经是负数违反量为 0罚项不生效。这种逐项叠加的方式非常方便扩展到几十个约束只需要在求和循环里累加即可。2.2 从不可行域逼近的收敛逻辑为什么 r 必须趋于无穷罚函数的直观行为是r越小目标函数f(x)的话语权越大迭代点可能停留在不可行域深处r越大罚项权重越高无约束极小点会被拖向可行域边界。理论上当r → ∞时罚函数的无约束最优解x*(r)收敛到原问题的最优解x*。一个关键性质是对于固定的有限rx*(r)通常落在外侧——它为了压低目标函数宁可付出一点约束违约的代价也不完全回到可行域。这正是“外点”名字的来源。比如下面这个例子min f(x) (x1-1)² (x2-2)² s.t. g(x) 3 - x1 - x2 ≤ 0 h(x) x1 - x2 0解析解是x* (1.5, 1.5)f* 0.5约束g取等号最优解在边界上。若取初始点x0 (0, 0)它既不满足g也不满足h完全在可行域外面迭代路径会沿着“外域”逐步贴向边界点。收敛逻辑可以用一组不等式直观说明。罚函数的无约束最优解满足∇f r( ... ) 0观察等式约束部分r * h(x)在收敛时趋于有限值这个值其实就是拉格朗日乘子λ的近似。它告诉我们两件事一是h(x)的违约量大致按1/r级别衰减二是迭代到后期乘子近似稳定罚函数问题与原问题在最优解处逐渐重合。所以程序里惩罚因子不能固定在某个值必须每轮乘一个mu 1递增。2.3 收敛判据约束违约度与相邻迭代差有了罚函数还必须定义“什么时候算收敛”。单纯看相邻两次罚函数值的变化不够可靠因为r很大时罚函数值被罚项主导数值上可能变化极小但约束还没对齐。常见的判据有两个约束违约度viol sqrt( Σ max(0,g_i(x))² Σ h_j(x)² )它直接度量当前点离可行域有多远相邻迭代点位移norm(x_k - x_{k-1})反映无约束优化是否已经稳定在某个不动点上。viol sqrt( max(0, g(x)).^2 h(x).^2 ); if viol 1e-7 norm(x - x_old) tol break; end这两个条件需要同时满足。只卡viol可能出现一种情况r很大但内层优化没收敛点还没走到罚函数极小点就误判成功只卡位移又可能在可行域深处提前退出因为目标函数那边梯度太平坦点几乎不动但约束还差得远。实践中我把违约度的阈值设得比位移阈值更严格通常取viol 1e-7、tol 1e-8。3. 外点法 MATLAB 程序实例完整脚本与逐段说明3.1 完整可运行脚本 outerpoint_demo.m下面这个脚本可以直接保存运行不依赖 Optimization Toolbox内层用最速下降加 Armijo 回溯线搜索。解的问题就是上一节那个带一个不等式和一个等式约束的小例子读者跑完可以把目标、约束替换成自己的问题。%% outerpoint_demo.m —— 外点法 MATLAB 程序实例 % min f(x) (x1-1)^2 (x2-2)^2 % s.t. g(x) 3-x1-x2 0, h(x) x1-x2 0 % 解析解: x* (1.5, 1.5), f* 0.5 % 目标函数及其梯度 f (x) (x(1)-1)^2 (x(2)-2)^2; grad_f (x) [2*(x(1)-1); 2*(x(2)-2)]; % 约束函数及其梯度 g (x) 3 - x(1) - x(2); grad_g (x) [-1; -1]; h (x) x(1) - x(2); grad_h (x) [1; -1]; % 罚函数r/2 * ( max(0,g)^2 h^2 ) P (x, r) f(x) r/2 * ( max(0, g(x)).^2 h(x).^2 ); gradP (x, r) grad_f(x) r * ( max(0, g(x)) * grad_g(x) h(x) * grad_h(x) ); % 参数设置 x0 [0; 0]; % 初始点位于不可行域外点 r 1; % 初始惩罚因子 mu 10; % 惩罚因子放大倍数 tol 1e-8; % 迭代位移容差 maxOut 50; % 外层最大循环次数 x_old x0 1; % 保证第一轮可以进入循环 history zeros(0, 4); % [r, 到真解距离, 违约度, fval] for k 1:maxOut % ---------- 内层无约束优化最速下降 Armijo ---------- for it 1:2000 d -gradP(x0, r); % 负梯度方向 if norm(d) 1e-12 break; % 梯度已为零内层收敛 end alpha 1; % 初始步长 c 0.4; % Armijo 条件常数 while P(x0 alpha*d, r) P(x0, r) c * alpha * (gradP(x0, r) * d) alpha alpha * 0.5; end x0 x0 alpha * d; end % ---------- 收敛判据与记录 ---------- viol sqrt( max(0, g(x0))^2 h(x0)^2 ); dist norm(x0 - [1.5; 1.5]); % 与解析解的距离仅本例可用 history [history; r, dist, viol, f(x0)]; if viol 1e-7 norm(x0 - x_old) tol fprintf(第%d轮收敛: r%.2e, x(%.8f, %.8f), f%.8e, viol%.2e\n, ... k, r, x0(1), x0(2), f(x0), viol); break; end x_old x0; r mu * r; % 惩罚因子倍增 end % 打印每轮记录 disp( r ||x-x*|| viol fval); disp(history);外层跑 50 轮内层跑 2000 步对二维问题只算罚函数和梯度毫秒级完成。运算量不是瓶颈重点是把外点法的循环骨架跑通内层在固定r下求无约束极小点外层判断约束违约和位移不达标就放大r继续。3.2 代码逐段拆解梯度计算与 Armijo 线搜索脚本里最容易被替换的是内层优化器。这里刻意没用fminunc因为很多机器上的 MATLAB 不带 Optimization Toolbox而最速下降二十行就能写完逻辑透明。负梯度方向d是下降方向但要配合“走多远”的策略否则梯度下降要么振荡要么慢如蜗牛。Armijo 回溯线搜索的思想很朴素先试alpha 1如果步长太大导致函数值下降不足反复减半直到满足P(x alpha*d) ≤ P(x) c * alpha * gradP * d其中c ∈ (0, 0.5)控制接受的下降量我取 0.4偏宽松避免一轮里回溯太多次数。这个条件保证每一步目标值都有足够下降程序里嵌在while循环里直到不等式成立才更新x0。gradP 表达式为什么没有二阶导因为二次罚项max(0,g)^2对x求导后是2 * max(0,g) * grad_g再乘上前面的r/2正好消掉系数 2变成r * max(0,g) * grad_g。对等式约束同理。如果读者把罚项系数写成了r而不是r/2梯度也要相应乘以 2这里最容易出错。3.3 运行结果与典型参数表运行脚本后每轮r对应的违约度大致按1/r的节奏下降。我没有列出精确输出因为不同 MATLAB 版本和浮点环境会有细微差异但量级关系稳定可参照下表核对程序行为惩罚因子 r约束违约度 viol 量级迭代点位置特征11e-1 左右明显在外域目标函数占主导1e21e-2 左右靠近边界罚项开始起作用1e41e-4 左右基本贴住边界x 接近最优1e61e-6 左右违约度低于大多数默认容差违约度与r近似反比是二次罚函数的典型性质。如果发现违约度下降速度远慢于这个规律先怀疑内层优化没收敛再怀疑约束函数数值量纲差异过大。3.4 扩展多约束、非光滑约束与 fminunc 替身真实问题往往不止两个约束。扩展方式是让g和h变成函数句柄数组罚项用循环累加g_list {...}; % 多个不等式约束 P_pen 0; for i 1:length(g_list) P_pen P_pen r/2 * max(0, g_list{i}(x))^2; endgradP也要同步叠加max(0, g_i) * grad_g_i。量纲差异大时给每个约束配独立的权重r_i比如位移约束和角度约束数值差几个数量级统一用一个r会导致小量纲约束几乎不被惩罚迭代点长期违反它。如果机器装有 Optimization Toolbox内层可以直接换用拟牛顿法替换掉最速下降循环options optimoptions(fminunc, Display, off, Algorithm, quasi-newton); x0 fminunc((t) P(t, r), x0, options);fminunc稳定、收敛快但对无约束问题底层实现也是迭代法和外点法循环嵌套没有问题。4. 外点法调参技巧与最优性验证4.1 r0 与 mu 的影响看约束违约度曲线初始惩罚因子r0和放大倍数mu是外点法最核心的参数直接决定收敛速度和数值稳定性。r0太小前几轮基本在优化无约束目标白白浪费迭代r0太大罚函数从一开始就病态梯度方向振荡。常见策略是r0 1起步观察第一轮违约度如果已经小于1e-3说明起点太靠近可行域可以适当调小r0看趋势。mu取 5 到 20 之间比较实用。mu 10是默认选择每轮违约度下降约一个数量级便于观察收敛曲线。mu偏大如 100能少跑几轮外层但惩罚因子从 1 跳到 100 时内层无约束优化的初始点突然变得“非常不可行”需要额外迭代才能跟上变化mu偏小如 2每轮进展太慢外层轮数翻倍。建议在调试阶段把每轮的viol打出来画出对数坐标下违约度随轮数的曲线斜率大致稳定的话参数的节奏就在合理区间。4.2 病态 Hessian 与内层优化器选择二次罚函数有个绕不开的代价r增大时罚函数在边界法方向的 Hessian 特征值以r量级增长条件数趋近1 r。最速下降法在线性条件下的收敛速度与条件数直接挂钩条件数增大后梯度方向几乎正交于指向最优解的方向迭代路径会走之字形。这时把内层换成拟牛顿法或者共轭梯度法效果立竿见影。% 用 fminunc 替换内层后外层判据不变 options optimoptions(fminunc, Display, off, ... Algorithm, quasi-newton, SpecifyObjectiveGradient, false);也可以继续用最速下降但放宽内层收敛标准——即内层不需要完全收敛跑固定步数就退出外层判断因为下一轮r变大后当前点本来就会被继续修正。这种“不完全内层优化”配合mu增大有时反而比每轮都死磕到高精度更省总耗时。4.3 用 KKT 残差和 fmincon 复核最优解程序说收敛还不够得用独立手段验证。一阶最优性条件要求存在拉格朗日乘子使得∇f(x*) Σ λ_i * ∇g_i(x*) Σ μ_j * ∇h_j(x*) 0 λ_i ≥ 0, λ_i * g_i(x*) 0外点法的罚项在收敛时给出乘子近似不等式乘子取λ_i ≈ r * max(0, g_i(x))等式乘子取μ_j ≈ r * h_j(x)。验证代码片段lambda r * max(0, g(x0)); % 不等式乘子近似 mu_eq r * h(x0); % 等式乘子近似 KKT_grad grad_f(x0) lambda * grad_g(x0) mu_eq * grad_h(x0); fprintf(KKT 梯度残差: %.3e\n, norm(KKT_grad));如果残差大于1e-5说明当前点和真正的最优之间还有距离。交叉验证最直接的方法是调用fmincon解原问题对比外点法的结果。两个解在1e-4量级内一致才能放心把罚函数换进自己的代码里。外点法虽然古老却是理解增广拉格朗日乘子法的基石。乘子法的改进思路正是在罚项后面再加上一项λ * h(x)的线性补偿让有限惩罚因子也能得到精确约束满足。把这里r的倍增节奏、违约度判据和内层病态问题看懂再往那边走就顺了。本文还有配套的精品资源点击获取