简介：压缩包内为基于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*。
一个关键性质是：对于固定的有限r，x*(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 量级 | 迭代点位置特征 |
|---|---|---|
| 1 | 1e-1 左右 | 明显在外域，目标函数占主导 |
| 1e2 | 1e-2 左右 | 靠近边界，罚项开始起作用 |
| 1e4 | 1e-4 左右 | 基本贴住边界，x 接近最优 |
| 1e6 | 1e-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的倍增节奏、违约度判据和内层病态问题看懂，再往那边走就顺了。
本文还有配套的精品资源，点击获取
