很多人一听到非线性优化、约束优化，第一反应就是打开MATLAB直接调fmincon。这个习惯没错，fmincon内部确实有一套非常成熟的算法框架，对大多数常规问题都能给出满意的解。但如果你做的是算法研究、嵌入式落地，或者需要把优化过程嵌到自研的仿真平台上，你迟早会碰到一个绕不开的问题：底层算法到底是怎么运转的？它的收敛性判断、子问题求解、步长修正每一步在做什么？这时候，手头有一套自己写的、能跑能看的SQP源码，价值就完全不一样了。
我这次整理了一套自编的序列二次规划（SQP）法求解非线性优化问题的MATLAB源码，并且配套了35个难度递进的测试算例，从简单的无约束凸函数到带非线性等式和不等式约束的工程优化题，每个算例都跑通了，附上了收敛过程曲线和最终的数值结果。这篇文章就围绕这套项目展开，把SQP的算法拆解、MATLAB实现上的关键细节、示例设计思路以及我在调试这35个例子过程中踩过的坑，一次性说清楚。
1. 从问题到方案：为什么偏偏选了SQP
1.1 非线性优化问题，难点到底在哪里
非线性优化问题从数学形式上就比线性规划麻烦得多，它的目标函数和约束条件至少有一个是非线性的。这意味着可行域可能是一个弯曲的、不规则的几何体，最优解可能在边界上，也可能在某个凸凹交错的谷底。经典的单纯形法、梯度法在面对非线性约束时常常施展不开，因为可行方向的计算本身就变得极其复杂。
举个例子，如果你要优化一个机械臂的关节角度，让末端执行器的定位误差最小，同时保证每个关节的力矩不超过电机极限，这就是一个典型的非线性约束优化问题。目标函数是末端误差的范数，约束条件里的关节力矩又和角度、角速度之间存在着强非线性关系。直接枚举角度组合不现实，一阶梯度法又很难在约束边界上精确地停住，所以必须依赖更高级的迭代算法。
1.2 SQP的核心思想：把难啃的骨头切成小块
SQP全称Sequential Quadratic Programming，核心策略非常朴素却极其有效：把一个复杂的非线性优化问题，在每一个迭代点处用一个二次规划子问题来局部逼近，然后反复求解这些相对简单的子问题，让迭代点逐步逼近原问题的KKT点。
为什么是二次规划？因为二次规划（QP）的理论和算法已经非常成熟，等式约束、不等式约束都能处理，求解效率高、数值稳定性好。SQP把一个一般的非线性规划问题在迭代点处进行Taylor展开，目标函数取到二阶项，约束函数取到一阶项，就得到了一个带线性约束的二次规划模型。求解这个QP子问题得到搜索方向，再通过适当的一维搜索确定步长，就完成了一次迭代。这个过程重复下去，本质上是在“用一系列好解的问题去逼近一个难解的问题”，而不是试图一步登天直接算出全局最优点。
1.3 和MATLAB自带优化工具箱相比，自编的价值在哪里
MATLAB的fmincon当然很强，它是商业级的实现，经过了无数工程案例的验证。但正因为它是封闭的，你只能看到输入输出，中间发生了什么对你是一个黑盒。我自编SQP的动机有几个：
- 算法教学和原理验证的需要，比如对比不同Hessian矩阵更新策略对收敛速度的影响，
fmincon做不到这么细粒度。 - 定制化的需要，比如在某些工业场景里，需要对变量加特殊的逻辑映射，或者需要每迭代一步就输出一批中间计算数据，
fmincon的输出控制没那么方便。 - 防止对工具箱的依赖，换一套没有优化工具箱的MATLAB环境，自编代码照样能跑。
当然，自编也有自编的代价，需要自己处理数值微分、矩阵正定性修正、线搜索策略等一大堆细节。但正是这些细节，构成了对一个优化算法最深层的理解。
2. 算法落地：SQP实现的几个关键环节
2.1 整体框架：迭代循环里到底发生了什么
标准的SQP迭代框架可以概括为以下几步。首先，给定初始点 (x_0) 和初始Hessian矩阵估计 (H_0)，通常用单位阵或有限差分近似。然后在每个迭代点 (x_k) 处，把原问题转化为QP子问题：
[ \min_{d} \quad \frac{1}{2} d^T H_k d + \nabla f(x_k)^T d ]
[ \text{s.t.} \quad \nabla g_i(x_k)^T d + g_i(x_k) \le 0, \quad i \in \mathcal{I} ]
[ \quad \nabla h_j(x_k)^T d + h_j(x_k) = 0, \quad j \in \mathcal{E} ]
其中 (g_i) 是不等式约束，(h_j) 是等式约束。这个子问题的解 (d_k) 就是搜索方向。接着，把搜索方向代入原目标函数和约束中，做一维搜索求步长 (\alpha_k)，得到 (x_{k+1} = x_k + \alpha_k d_k)。最后，用BFGS公式或者SRI公式对Hessian矩阵进行修正，更新 (H_{k+1})，进入下一轮迭代。当梯度范数、约束违反量和步长都小于阈值时，迭代终止。
2.2 QP子问题求解的实现取舍
QP子问题的求解是SQP的核心之一。我在这套代码里采用的策略是：先用MATLAB内置的quadprog求解QP子问题，然后对quadprog返回的无解或数值异常情况做后备处理。这样做的理由很实际——自编一套完整有效的active-set QP求解器工作量巨大，而且数值稳定性很难超过MATLAB官方的实现。但整个SQP的主逻辑、Hessian更新、线搜索、收敛判断都是自己写的，所以核心算法依然是SQP，而没有把整个问题丢给fmincon。
这里有个非常容易踩的坑：quadprog默认求解的是严格凸二次规划，如果Hessian矩阵不正定，它会直接报错。所以必须对Hessian矩阵做修正。我采用的修正方案是：检查每个特征值，小于阈值就把特征值强制拉到一个较小的正数上，比如(10^{-8})，然后重组矩阵。这个修正也叫“Hessian正则化”，看似简单，实际对稳定性提升很大。
2.3 一维搜索、罚函数与步长修正
获得了QP子问题的解 (d_k) 之后，直接令 (\alpha_k = 1) 未必稳妥，因为子问题只是局部近似，大步长可能导致目标函数上升或者约束严重违反。这时需要一维搜索。标准做法是构造一个“效益函数”或“罚函数”，在方向 (d_k) 上对步长进行一维极小化。
我采用的是经典的 (L_1) 罚函数形式：
[ \Phi(x; \mu) = f(x) + \mu \sum \max(0, g_i(x)) + \mu \sum |h_j(x)| ]
其中罚因子 (\mu) 的选取有一点学问：如果太小，约束违反得不到充分惩罚；如果太大，会过于限制步长，导致收敛变慢。在我的代码里，罚因子的初始值设为试探值，每次迭代时若约束违反量没有下降，就按比例放大 (\mu)。这一步是实际调试中最容易出问题的地方之一，后面会在常见问题里详细展开。
2.4 梯度信息与数值微分的选择
SQP需要目标函数和约束函数的梯度信息。对于理论测试算例，当然可以用解析梯度，这样收敛快、精度高。但为了适用性更广，我的代码默认采用中心差分格式来自动计算梯度，即每个变量加一个小扰动 (\epsilon) 和减去一个小扰动，计算函数值之差除以 (2\epsilon)。这个做法的优点是无需为每个新问题手动推导梯度表达式，缺点是计算量变大，而且对 (\epsilon) 的取值很敏感。(\epsilon) 太小会因浮点舍入误差而失真，太大会因为泰勒展开截断误差而失真。我在这35个示例里普遍使用 (\epsilon = 10^{-6})，对多数连续可微测试函数效果都不错；如果函数梯度变化非常剧烈，需要适当调大到 (10^{-5})。
2.5 收敛判据的三重门
收敛判据决定了算法什么时候“觉得够了”。我设置了三个条件，必须同时满足才判定收敛：
- 梯度投影范数：( abla f(x_k)) 在可行方向上的投影小于阈值 (10^{-6})。
- 约束违反量：所有等式约束的绝对值和不等式约束的正部最大值小于 (10^{-6})。
- 迭代步长：(|d_k|) 小于 (10^{-8})。
三个条件同时满足，说明当前点既满足“一阶必要性”，又“站在可行域内”，还“动不了了”。这三个条件缺一不可，单看梯度不看约束，很容易把点停在不可行区域。
3. 35个示例的设计套路：从简单到硬核的递进
3.1 示例分类与难度梯度安排
35个示例不是随便凑数的，而是按照“变量数量、约束类型、非线性强度、特殊构造”四个维度做了系统排布。大致可以分成以下几类：
| 类型 | 编号范围 | 特点 | 典型算例 |
|---|---|---|---|
| 无约束凸优化 | 例1-5 | 目标函数为简单的二次或四次函数 | Rosenbrock函数 |
| 带边界约束 | 例6-10 | 变量有上下界限制 | 边界盒约束下的多项式优化 |
| 线性约束非线性目标 | 例11-15 | 约束为线性等式/不等式 | 投资组合优化 |
| 非线性等式约束 | 例16-20 | 约束包含非线性方程 | 多体机构装配误差最小化 |
| 混合非线性约束 | 例21-30 | 同时包含等式和不等式约束 | 压力容器设计、齿轮传动比优化 |
| 工程基准题 | 例31-35 | 经典工程设计测试题 | Himmelblau问题、焊接梁设计等 |
这个梯度安排的最大价值是：如果你逐例跑下来，你能直观地感受到问题复杂度上升时，SQP算法在迭代次数、Hessian修正频率、步长选择策略等方面的行为变化。这是任何教科书都给不了的直觉经验。
3.2 经典题型的构造逻辑：目标函数、约束条件怎么搭
以经典的Rosenbrock函数为例，目标函数为 (f(x) = 100(x_2 - x_1^2)^2 + (1 - x_1)^2)。这个函数的等高线沿着一个窄而弯曲的“香蕉形”山谷延伸，最优点是 ((1,1))。由于山谷曲率变化很大，一阶梯度法很容易在山谷里震荡，但SQP配BFGS修正后，能够很快识别出弯曲方向，十几步内就收敛。
而像压力容器设计这样的工程题，目标函数是材料成本，决策变量是容器半径、壁厚等，约束条件包括应力极限、几何尺寸下限等。这类问题高度非线性，而且不等式约束很多，初始点往往不可行，这正好检验SQP处理不可行初始点的能力。
3.3 工程风的算例：为什么拿“经典基准题”当试金石
在工程设计领域，有一批流传了几十年的经典测试问题，比如Welded Beam Design（焊接梁设计）、Pressure Vessel Design（压力容器设计）、Three-Bar Truss（三杆桁架设计）。这些问题的数学模型并不复杂，但约束条件中充斥着大量非线性不等式，变量之间的耦合关系非常强，而且最优点经常位于多个约束边界的交汇处。
我用这套自编SQP跑这些经典题时，重点观察的是：算法能否从多个不同的初始点出发最终收敛到同一最优解附近。比如压力容器设计，经验最优解位于整数变量和连续变量混合的边界上，SQP如果不处理整数变量，会在连续空间里给出一个参考最优值，配合舍入策略可以作为工程初步设计值。这类结果对于做结构优化和机械设计的人是很实用的参考。
3.4 不同初始点、不同步长策略的结果对比
同一道题，初始点取不同位置，SQP的表现差异可能非常大。我在示例里对部分题目专门做了初始点敏感性分析，比如把初始点分别取在可行域内部、边界外很远、以及高度非凸区域，然后对比迭代路径。结果显示，SQP对初始点位置的敏感性明显低于普通梯度法，但在高度非凸问题中，不同初始点仍可能收敛到不同局部最优点。这提醒我们，SQP本质上是局部优化算法，全局搜索需要搭配多起点策略或者全局优化框架。
3.5 检查单调性、极值、可行性等细节
每跑完一个例子，我不仅看最终优化结果，还会检查三个衍生指标：目标函数随迭代次数的下降曲线是否单调、KKT残差是否随着迭代不断缩小、每个约束在最优点的裕量是多少。这些检查看起来零碎，实际非常有用。比如目标函数单调性是判断线搜索是否有效的重要依据，如果目标函数曲线出现明显上跳且后线没有恢复，基本可以断定罚因子或线搜索策略有bug。
4. 代码结构、输入输出设计与结果可视化
4.1 主函数、子函数、示例脚本的模块划分
整套代码模块划分成三层：
- 核心层：主程序
sqp_solver.m，负责迭代流程控制，内部调用来解决QP子问题。 - 功能层：
compute_objective.m、compute_constraints.m、compute_gradient.m、bfgs_update.m、line_search.m、hessian_modify.m等子函数，分别封装求解流程中的独立环节。 - 示例层：35个独立的示例脚本，分别命名为
example_01_rosenbrock.m、example_02_quadratic_simple.m等，每个脚本定义自己独特的目标函数和约束函数，调用核心层完成求解并绘制结果图。
这种分层的最大好处是可替换性极强。想测试一个新的优化问题，只需要新写一个示例脚本，定义好目标函数和约束函数，不需要改动核心层任何代码。想更换Hessian更新策略，比如从BFGS换成DFP，只需要修改bfgs_update.m一个文件。
4.2 输入参数约定与收敛判据
核心求解器的函数签名如下：
function [x_opt, f_opt, info] = sqp_solver(fun, con, x0, options)其中fun是目标函数指针，con是约束函数指针，约定返回两个结构体字段：等式约束值向量ceq和不等式约束值向量cineq。约束的编写约定为：不等式约束以 (g(x) \le 0) 的形式给出，等式约束为 (h(x) = 0)。这个约定如果不统一，算法计算结果一定会出错。
options结构体里主要字段包括：最大迭代次数、梯度计算步长、收敛阈值、罚因子初始值、是否显示迭代过程。全部有默认值，一般用户不需要逐一配置。
options = struct(... 'max_iter', 500, ... 'tol_f', 1e-8, ... 'tol_con', 1e-6, ... 'eps_diff', 1e-6, ... 'mu_init', 10, ... 'verbose', true);4.3 结果展示：表格与收敛曲线
每个示例运行结束后会输出三样东西：
- 终端表格：包含迭代次数、目标函数值、最大约束违反、步长、梯度范数等信息。信息量很大，方便做逐轮分析。
- 收敛曲线图：以迭代次数为横轴，目标函数值和约束违反量分别画两条对数坐标曲线，观察收敛趋势。
- 最优点截图：将最终结果与理论最优解（有解析解时）做对照，输出绝对误差和相对误差。
用对数坐标查看约束违反量下降趋势是很好的习惯。正常情况下，约束违反量应当呈锯齿状但总体下降，最终跌破阈值。如果约束违反量在某几次迭代中暴涨，说明罚因子更新或者线搜索存在隐患，需要回头排查。
5. 常见问题与排查技巧实录
5.1 初始点在不可行域，迭代直接发散怎么办
这是跑非线性约束优化时最常见的情景。初始点不在可行域内部，而是落在某个约束的违反地带，SQP直接把子问题搭在一个不可行点上，QP子问题本身就可能无解。
我的处理办法是分两步走：先做可行性修复，暂时把目标函数优先级降低，先用最小二乘方式求解约束违反量最小化问题，把迭代点拉回可行域附近，再切换到正常SQP流程。这个技术也叫“弹性约束法”，虽然牺牲了少量计算时间，但显著提升了算法的鲁棒性，尤其适合工程场景，因为工程初始参数往往是拍脑袋设的，不会恰好满足全部约束。
5.2 数值微分误差导致收敛慢或卡在非最优点
中心差分虽然实现简单，但在目标函数是分段连续或者包含绝对值项时，会出现差分失真。比如某个函数在某点附近有剧烈尖峰，差分步长 (\epsilon) 跨过了尖峰，算出的梯度方向完全反了。
排查这类问题的方法是：把渐变和突变分开测试，对每个变量单独绘制目标函数剖面图，观察函数在局部是否平滑、差分步长是否跨过间断点。如果发现这种情况，解决办法是改用解析梯度，或者在函数里特殊处理不可导点，比如用平滑近似替代绝对值函数。
5.3 不等式约束数量偏多时，QP子问题求解失败
不等式约束多了之后，QP子问题的active set会变得很复杂，quadprog有时会返回“遇到数值问题”的警告。我遇到这类情况，首先检查约束函数的量纲是否统一。很多工程问题里，约束值差了好几个数量级，比如一个约束数值在(10^3)量级，另一个在(10^{-3})量级，合到一起后大规模数值差异会导致矩阵条件数变差。
统一量纲的方式很简单，在约束函数计算返回值时做缩放处理，让所有约束数量级落在0.1到10之间。这个预处理往往能立竿见影地解决QP子问题求解失败的问题。
5.4 Hessian修正频率过高，算法退化成梯度法
Hessian矩阵修正本身是好事，但如果每一次迭代的特征值修正都出现在对角线上，说明BFGS更新质量极差，Hessian矩阵离正定越来越远，算法实际退化成带步长的梯度法。
出现这种现象，通常是因为线搜索做得不充分，导致BFGS更新公式里的差值向量 (y) 和 (s) 不满足曲率条件。解决办法是加强线搜索，让步长满足Armijo条件和Wolfe条件，确保每次迭代更新都是“有益”的更新。
6. 操作实录：把一个具体算例完整跑通的流程
这一节拿35个示例里的第21个示例来说，它是一个带两个非线性等式约束、五个不等式约束的测试题，目标函数是一个四变量多项式。完整流程如下：
第一步，在示例脚本里定义目标函数：
function f = obj_fun(x) f = x(1)^4 + 3*x(2)^2 + x(3)^3 + (x(4) - 1)^2; end第二步，定义约束函数：
function [ceq, cineq] = con_fun(x) ceq = [x(1)^2 + x(2)^2 + x(3)^2 + x(4)^2 - 10; x(2)*x(3) - x(1)*x(4) - 2]; cineq = [x(1)^2 + x(4)^2 - 3; 4*x(2) - 2*x(3) - 1; x(3)^2 + x(4) - 5; 0.5 - x(1); x(2) - 2]; end第三步，在主脚本中调用求解器：
x0 = [2; 2; 2; 2]; options.verbose = true; [x_opt, f_opt, info] = sqp_solver(@obj_fun, @con_fun, x0, options);跑完之后观察终端输出，重点看每一步的约束违反量的变化。我实测下来，这个例子约28步收敛，目标函数从初始点的约58下降到最优点附近的约9.3，最大约束违反量从初始的10.6降到低于(10^{-7})，迭代曲线整体没有出现明显的发散抖动。
这个操作流程也展示了最核心的使用模式：任何新问题，只需替换目标函数和约束函数，即可套用这套流程完成优化求解。
7. 经验复盘：这套35示例项目带给我最深的三点体会
第一点是中间过程的可视化比最终结果更重要。如果一个优化算法只告诉你答案，那你永远无法判断它有多可靠。我在35个示例中强制输出每一轮迭代的梯度、约束违反、步长和Hessian修正次数，才有了足够多的数据去对比不同策略之间的优劣。
第二点是罚因子更新策略是SQP的隐形瓶颈。很多人关注Hessian更新和QP子问题求解，但忽略了罚因子这个看似外围的参数。实际调试发现，罚因子选得不好，要么迭代后期在约束边界附近反复跳，要么前期搜索被约束严重压制走不动。这套代码里我采用的是自适应罚因子，虽然代码只有十几行，却是让35个示例全部跑通的重要一环。
第三点是对局部最优解要有清醒认识。SQP很强大，但它不保证全局最优，35个示例里有好几个多解问题。跑这类问题时，我从多个初始点出发，对比收敛到不同局部最优解的目标函数值，最后人为筛选出最小者。要做工程决策时，这一步骤绝对不能省。
如果你也想把优化算法吃透，我建议你亲手把SQP框架从零搭建一遍，不要满足于调工具箱接口。能把35个不同难度的问题真正跑通，确保每一个的收敛曲线、约束违反量下降趋势都符合预期，你对于非线性优化算法的理解会比只看十遍教材都更扎实。遇到跑不通的案例，不要急着怀疑算法，先检查梯度的数值计算、约束的量纲、Hessian的正定性，这三点修好了，大部分问题都能解决。
