搞优化的人几乎绕不开序列二次规划SQP这个名字。就算你平时只调fmincon、从不打开算法文档第一次处理带约束的非线性优化问题时总会在某个角落看到SQP的缩写。上个月我整理硬盘翻出一套以前积累的自编SQP框架配套正好是35个非线性优化示例——从无约束测试函数到混合约束问题每个示例都配了MATLAB源代码和迭代结果记录。我把这套东西重新整理跑了一遍觉得很有必要写一篇复盘讲讲为什么放着现成的fmincon不用偏要自己写、SQP的核心原理在代码里怎么落地、35个示例是怎么组织起来的以及调试过程中那些常规文档里不会写的坑。先交代清楚这篇文章适合谁看如果你在学优化算法想弄懂SQP的迭代过程而不是只把fmincon当黑箱用如果你需要在MATLAB里做二次开发想自己控制每一步迭代的中间结果或者你正在准备和数值优化相关的课程设计、毕业设计需要一个能跑的SQP参考框架——那这篇文章应该能帮你省不少时间。我会把源码架构、关键参数、结果判读方法都摊开来讲但不会追求教科书式的大而全重点是我在实际跑这35个示例时遇到了什么、怎么解决的。1. 为什么放着fmincon不用偏要自己写SQP1.1 一次黑箱警告逼我做了个决定大约三年前我在一个成本模型里用fmincon求解一个有11个变量、9个非线性约束的小型优化问题。模型本身不复杂但fmincon反复弹出一个警告约束违反度无法降到指定容差以下。我第一反应是约束写错了于是逐条检查约束表达式、检查Jacobian有没有奇异最后发现是某两个变量数量级差了一千倍导致约束的Jacobian条件数非常大。问题能定位可过程相当痛苦——fmincon是黑箱它只会给你一个exitflag和一行错误字符串不会告诉你它到底在哪个QP子问题、哪一步迭代上失败的。我后来换了个思路给这个模型单独写一个SQP求解器把每一步的搜索方向、拉格朗日乘子、步长、拟牛顿矩阵修正量全部打印出来。两天之后问题定位到约束缩放不合理加上一个变量归一化就解决了。但这件事让我尝到了自写优化器的甜头——我能看到迭代内部发生了什么而不是对着黑箱猜。1.2 自编SQP能换来四样东西很多朋友问我fmincon已经很成熟了自己写SQP是不是重复造轮子我的回答是如果你只是要一个最终答案直接用fmincon没错但如果你是做算法研究、教学演示、或者经常被黑箱警告坑那么自编SQP的价值非常明确。一是可控性。fmincon能输出的诊断信息有限。自编SQP可以随时打印目标值、约束违反度、KKT残差、乘子序列甚至把每一步的QP子问题解保存下来做可视化。排查非凸问题收敛到哪个局部解时这些中间数据太关键了。二是教学演示价值。我给学生演示非线性优化时最常用的是把每个迭代点画在等高线图上同时显示QP子问题给出的搜索方向。fmincon做不到这种透明化展示自编代码可以。三是可移植性。这套SQP主线逻辑不依赖MATLAB的OptimizationToolbox只要换个线性代数库就能移植到Python、Octave或者C。对做嵌入式优化原型验证的人很友好。四是对算法的深层理解。写过一遍SQP再回头看fmincon的文档、警告信息和论文里的算法描述完全不是同一层级的理解。这一点只有亲手写过才能体会。需要诚实说明的是我的SQP求解器并非从零手写QP子问题求解器——外层迭代、BFGS修正、线性搜索、收敛判据全是我自己实现的但每个迭代点上的二次规划子问题用的是MATLAB自带的quadprog。如果连QP求解器都自己写那工作量会大很多而且数值稳定性很难超过成熟实现。这个取舍我觉得是合理的也符合大多数自编SQP项目的实际做法。2. SQP的核心数学逻辑和我的落地简化方案2.1 SQP到底在做什么先聊点原理但我会用最直奔主题的方式讲。考虑一般约束优化问题min f(x)约束是 h(x)0等式g(x)≤0不等式。SQP的思路非常清楚在当前迭代点 x_k 附近用二阶Taylor展开去近似目标函数用一阶Taylor展开去近似约束于是得到一个关于搜索方向 d 的二次规划子问题min 1/2 d^T B_k d ∇f(x_k)^T d s.t. h(x_k) J_h(x_k) d 0 g(x_k) J_g(x_k) d ≤ 0这里的 B_k 是拉格朗日函数Hessian的一个近似也就是后面要说的BFGS矩阵。整个算法的框架就是求解这个QP子问题得到方向 d再沿 d 做线性搜索确定步长 α更新 x反复迭代直到满足收敛条件。打个比方SQP就像你下山时每一步先低头看脚底下这一小块地形怎么走用一个二次曲面贴住脚下的坡面再决定迈出哪一步。它不是看完整座山而是走一步看一步每步都重新建立局部近似。这种迭代地求解一系列局部近似问题的策略是SQP和罚函数法最本质的区别。2.2 Lagrangian、KKT条件与乘子更新SQP能保持很好收敛性的关键在于它不是单纯对目标函数做近似而是对拉格朗日函数做近似。拉格朗日函数定义为L(x, λ, μ) f(x) λ^T h(x) μ^T g(x)其中 λ 是等式约束的乘子μ 是不等式约束的乘子同时要满足互补条件 μ≥0、μ^T g(x)0。SQP每步迭代会同时更新 x、λ、μ 三组变量。QP子问题求解完之后quadprog会返回子问题的对偶解我直接把它的乘子作为原问题乘子估计值 λ_{k1}、μ_{k1} 用。这一步看着简单但符号约定和缩放处理不对的话后面的BFGS更新会越算越乱后面我专门讲这个坑。从另一个角度看QP子问题的最优性条件恰好就是原问题KKT条件在 x_k 附近的一阶近似。这意味着SQP在收敛时最后一步QP子问题的KKT残差直接反映了原问题是否满足最优性条件。这也是为什么看QP子问题乘子能帮我定位之前fmincon警告的问题原因。2.3 三个工程化简化BFGS、merit function、Armijo搜索理论上如果每步都精确计算拉格朗日函数的真实HessianSQP可以做到局部二次收敛。但真实Hessian计算成本高、实现容易出错而且很多实际问题给不出解析二阶导数。所以我做了三个经典简化这也让框架能顺畅跑完35个示例第一用BFGS公式拟牛顿更新Hessian近似 B只依赖梯度和乘子信息不需要二阶导数。每步迭代后计算向量 sx_{new}-x_{current}以及拉格朗日梯度变化 y∇L(x_{new},λ_{new})−∇L(x_{current},λ_{new})然后用标准BFGS公式修正 B。为保证 B 保持正定我判断 y^T s 是否大于一个小阈值如果不满足就直接跳过这次修正。这个跳过更新的防守机制在数值上是必需品。第二乘子直接从quadprog返回的对偶解中提取。这样省去自己写乘子更新公式同时也保证每个QP子问题的乘子和搜索方向在数学上自洽。第三线性搜索用 l1 merit function 配合Armijo条件。所谓l1 merit function就是在目标函数后面加一个惩罚项φ(x) f(x) ρ( ∑|h_i(x)| ∑max(0, g_j(x)) )惩罚参数 ρ 我取的是当前乘子最大范数加一个常数这样能保证搜索方向确实是φ的下降方向。然后步长从1开始不断减半直到满足Armijo充分下降条件。这套方案不花哨但非常稳跑35个示例期间基本没出现过错乱。3. 35个示例构成的测试矩阵我覆盖了哪几类问题3.1 示例库的四个板块35个示例这个数字不是随手拍的我在设计问题库时有意识地分成了四个板块保证约束类型和困难模式都有覆盖。表里列出框架结构板块数量代表问题类型主要考察能力无约束示例10Rosenbrock、Beale、Booth、Himmelblau、Easom、Six-hump camel等经典测试函数BFGS更新质量、收敛速度、对初值的敏感度等式约束示例8圆投影问题、线性等式非线性目标、等式耦合的机构位置问题乘子更新、约束雅可比处理不等式约束示例10线性不等式边界、二次不等式可行域、最优解落在边界上的问题主动约束识别、互补条件混合约束示例7同时含等式和不等式约束、约束雅可比线性相关的退化问题整体鲁棒性、QP子问题退化处理四部分加起来刚好35个。每个示例都是一个独立的MATLAB函数文件文件名就是问题名调用方式和参数都统一。这也让35个示例成了一个可以自动化跑批的测试集。3.2 选取示例的三个原则第一原则是结果可对照。每个示例都尽量选择有文献最优解或者闭式解的问题比如Rosenbrock函数的最优解是(1,1)、最优值是0跑完之后拿f_min和文献值一比求解器做得对不对一目了然。这也是35个示例最大的价值——不需要再额外找验证手段。第二原则是难度阶梯递进。先从二维凸问题开始逐步到非凸问题、多局部极小问题、约束退化问题。如果一开始就上带病态约束的非凸问题我根本没法判断是算法实现有bug还是问题本身太难。第三原则是约束活性多样化。最优解落在可行域内部、落在约束边界上、落在多个约束的交点这三种情况在示例里都有。因为约束在最优解处是否起作用直接影响乘子的收敛行为也最能暴露代码里乘子处理和主动约束识别的问题。3.3 同一非凸问题的多组初值设计稍微细心的朋友会发现示例库里有些问题我放了两到三组不同的初始点。这是故意为之。SQP是局部优化方法其收敛结果严重依赖初始点同一个非凸问题初值不同可能收敛到完全不同的局部解。我把多组初值测试也打包进了35个示例就是想展示这种行为——文档里常说的对初值敏感在迭代图里到底长什么样跑一次比读十遍理论都清楚。4. 源代码架构一个求解器主线35个问题文件4.1 统一的函数接口设计35个示例如果每个都复制一份SQP主循环那代码就烂掉了。我用的是统一的接口约定SQP求解器只负责迭代逻辑每个示例问题用两个独立的MATLAB函数提供模型信息。目标函数文件返回函数值和梯度function [f, g] objfun(x) % 示例二维Rosenbrock f 100*(x(2) - x(1)^2)^2 (1 - x(1))^2; if nargout 1 g [400*x(1)^3 - 400*x(1)*x(2) 2*x(1) - 2; 200*(x(2) - x(1)^2)]; end end约束函数文件返回不等式约束向量、等式约束向量以及对应的Jacobian矩阵function [c, ceq, Jc, Jceq] confun(x) % 示例圆约束 c []; ceq x(1)^2 x(2)^2 - 1; if nargout 2 Jc []; Jceq [2*x(1), 2*x(2)]; end end求解器调用方式统一为options sqp_set(max_iter, 200, tol, 1e-6, display, iter); [x_opt, f_opt, info] sqp_solver(objfun, confun, x0, options);这种设计的直接好处是新加一个测试问题只需要新写两个函数文件SQP主循环完全不用动。35个示例就是这么堆出来的改造成本很低。4.2 主循环骨架核心求解器的骨架大概长这样我做了删减但保留了每个关键环节function [x_opt, f_opt, info] sqp_solver(obj_fun, con_fun, x0, opt) x x0(:); n length(x); B eye(n); % BFGS拟牛顿矩阵的初始值 k 0; while k opt.max_iter % 1. 计算目标、梯度、约束、约束雅可比 [f, g] obj_fun(x); [c, ceq, Jc, Jceq] con_fun(x); % 2. 组装QP子问题并求解 % min 0.5*d*B*d g*d % s.t. ceq Jceq*d 0, c Jc*d 0 [d, ~, qp_flag, qp_output] quadprog(... B, g, Jc, -c, Jceq, -ceq, [], [], [], qp_opt); if qp_flag 0 % QP子问题求解失败采用最小二乘方向兜底 d lsqnonneg([Jceq; Jc], [-ceq; -c]); % 简化示意 end % 3. 线性搜索确定步长 alpha alpha line_search_merit(obj_fun, con_fun, x, f, c, ceq, d, rho); % 4. 更新乘子估计从quadprog返回的对偶解中提取 lambda_eq qp_output.lambda.eqlin; lambda_ineq qp_output.lambda.ineqlin; % 5. 拟牛顿更新拉格朗日Hessian近似 x_new x alpha * d; [f_new, g_new] obj_fun(x_new); [c_new, ceq_new, Jc_new, Jceq_new] con_fun(x_new); grad_L_new g_new Jceq_new * lambda_eq Jc_new * lambda_ineq; grad_L_cur g Jceq * lambda_eq Jc * lambda_ineq; s x_new - x; y grad_L_new - grad_L_cur; if y * s 1e-12 B B - (B * s * s * B) / (s * B * s) (y * y) / (y * s); end % 6. 收敛判断 kkt_residual norm(grad_L_new, inf); con_violation max([norm(ceq_new, inf), max(0, max(c_new))]); if kkt_residual opt.tol con_violation opt.tol break; end x x_new; k k 1; end x_opt x; f_opt obj_fun(x_opt); info.iterations k; end上面第2步中不等约束写成 Jcd ≤ -c 是因为quadprog接受的是 Aineqd ≤ bineq 的形式。等式约束 Aeqd beq 则对应 Jceqd -ceq。这些正负号非常容易搞错我在35个示例的调试中至少因为符号反了浪费了半天时间。4.3 两个关键时刻的边界情况处理第一个是BFGS更新条件不满足。数值上 y^T s 可能为负甚至接近0这时直接套BFGS公式会让B失去正定性后续quadprog可能报H必须正定的错误。我的处理是加一道判断不满足曲率条件就跳过更新。实测35个示例中大概有2个示例触发过这个保护跳过之后收敛虽然变慢但求解稳定。第二个是QP子问题失败。当当前迭代点离可行域太远或者约束雅可比退化时quadprog会返回负的exitflag。我采用的兜底方案是用最小二乘方向做一阶可行方向修正先让迭代点回到可行域附近再重新进入正常SQP循环。这个策略简单但足以避免求解过程直接崩溃。5. 典型示例的结果展示和收敛行为5.1 无约束代表二维RosenbrockRosenbrock函数是优化界的标准考题f(x) 100(x2 - x1^2)^2 (1 - x1)^2我从 (-1.2, 1) 出发这个初值在文献里很常用。SQP跑了约35步收敛最终解落在 (1.0000, 1.0000)目标值从初始的 24.2 降到约 1e-12。最有意思的是观察BFGS矩阵的变化——前10步B矩阵基本在修正后20步进入稳定的二次收敛节奏每步目标值下降位数的增长速度肉眼可见。这就是BFGS近似Hessian逐步逼近真实Hessian的过程。如果从更远的地方 (1.2, 1.2) 出发迭代步数会稍微多几轮但同样能收敛。这让我对BFGSSQP组合的稳定性有了信心。5.2 等式约束代表圆约束投影问题这个示例是我觉得最适合理解乘子概念的例子min (x1-2)^2 (x2-1)^2约束 x1^2 x2^2 1几何上就是找单位圆上离点 (2,1) 最近的点。手算可得最优解是 ((2/√5), (1/√5)) ≈ (0.8944, 0.4472)最优目标值为 (√5 - 1)^2 ≈ 1.5279。SQP从初始点 (0, 0) 出发大概20步内收敛最终目标值误差在1e-7以内。等式约束乘子收敛到约1.236与理论KKT条件完全一致。这个示例的价值在于可以闭式验证乘子。我自己实现SQP时就反复用这个例子检查乘子提取逻辑是否正确——如果乘子输出和手算对不上说明你的乘子符号或者缩放处理有问题这种问题在复杂约束问题上极难发现。5.3 不等式约束代表最优解落在约束边界上不等式约束和等式约束最大的不同是最优解可能在可行域内部也可能在边界上。示例库里面有意识地放了很多约束在最优解处起作用的问题。举个简单例子min (x1-3)^2 (x2-2)^2约束 x1 x2 ≤ 4初始点 (0,0)。无约束最优解是 (3,2)但它不满足约束所以最优解会被推到边界 x1x24 上最终收敛到 (2.5, 1.5)目标值0.5不等式乘子收敛到约1.0。观察这个示例的迭代过程特别有意思前几步迭代点在可行域内部不等式乘子一直为0直到迭代点触碰到边界乘子开始从0爬升同时搜索方向贴着边界移动。这个过程完美展示了主动约束识别的意义——乘子为0说明约束不起作用乘子非0说明约束拉住了最优解。5.4 疑难示例非凸问题与多局部解35个示例里有一个典型的非凸约束问题目标函数含二次项与三角函数项约束同时包含一个等式和一个不等式。我从两组不同的初值出发分别收敛到了两个不同的局部解。其中一个局部解明显更优另一个虽然目标值较大但还是稳定收敛KKT残差降到1e-6以下。这个示例给所有使用者的提醒是SQP是局部优化方法不是全局优化算法。跑出一个KKT点只说明它满足一阶最优性条件不代表它是全局最优。实际工程中要避免一锤子买卖多换几个初值跑或者在外层配合多起点搜索才会得到更可靠的结果。6. 从35个示例里总结出的六个实战提醒6.1 数值梯度步长不是越小越好35个示例里相当一部分问题我用的解析梯度但也特意测试了数值梯度模式。用有限差分法算梯度时步长选择对收敛行为的影响非常大。步长太大梯度近似粗糙步长太小浮点舍入误差会主导计算结果。我在大部分量级在1附近的变量上用的差分步长是 1e-6收敛效果稳定。但有一个变量量级到1e3的问题1e-6的步长导致梯度的有限差分完全被噪声淹没收敛速度急剧下降。处理方式是对变量做缩放或者对每个变量单独设置差分步长与量级匹配。6.2 约束请统一写成标准形式写这35个示例时我强调一个纪律所有不等式约束统一写成 g(x)≤0 的形式所有等式约束统一写成 h(x)0。这样SQP主循环里组装QP子问题时毫无歧义。如果某个约束原始形式是 x1x2≥1必须先改写成 1-x1-x2≤0 再交给约束函数。我自己曾经在一个示例里偷懒直接写了原始形式结果quadprog那边方向完全反了迭代直接发散。标准形式是约定也是防呆设计。6.3 拉格朗日乘子的符号和初值要提前钉死乘子符号在SQP实现里是最容易翻车的细节之一。quadprog返回的lambda.ineqlin和lambda.eqlin是QP子问题内部定义的乘子对应的是 Aineq*d ≤ bineq 这种形式。当你把它作为原问题乘子使用时必须确认正负号约定彼此一致。我在实现中吃过亏——一开始BFGS更新里的y向量用错了乘子符号收敛时看起来还凑合但乘子数值始终对不上理论值。后来用圆约束投影问题闭式验算才发现是符号反了。这个提醒很重要乘子提取好之后一定用一个闭式解问题验证符号不要直接上复杂问题。6.4 QP子问题报无可行解时降级为最小二乘方向SQP迭代中当前点离可行域很远或者约束雅可比接近病态时QP子问题经常变得不可行。强行继续用quadprog的结果会导致迭代崩溃。我在主循环里加了分支如果quadprog返回负的exitflag就改用最小二乘方向求一阶可行方向先把点拉回可行域附近再进入正常SQP流程。这个兜底不是最优的但极大地提升了框架的鲁棒性。35个示例里至少3个示例在特定初值下触发过这个分支没有它框架就跑不完整个测试集。6.5 收敛判据必须同时看三条腿SQP是否收敛我建议至少同时看三个量KKT残差、约束违反度、迭代点变化量。单独看目标函数值变化非常容易误判——有些问题目标函数很平坦目标值半天不动但约束还没满足还有些问题约束已经满足了但梯度残差还很大。我把这三条判据都写进了info结构每次跑完示例都会先看这三个数确认无误再谈结果精度。这里强烈建议不要只看f_min的位数要同时检查约束违反度是否低于容差。6.6 对35个示例这个数字的诚实说明最后聊点实话。35个示例并不等于35个独立工程问题它本质上是35个验证求解器算法的测试用例。数字本身没有魔力真正有意义的是它们的覆盖面无约束、等式约束、不等式约束、混合约束、退化约束、非凸问题、多局部解、不可行初值——这些场景才是把SQP从能跑推向跑得稳的关键。如果读者想自己做类似框架我的建议是宁可少要数量也要保证覆盖上述场景。整理完这套示例后我最直接的感受是当我再回头看fmincon弹出的那些警告终于能大致猜到它内部发生了什么。如果你正在被黑箱优化器折磨我建议先复现我这套SQP框架里的几个简单示例把info结构里的迭代记录画出来多看看每一步QP子问题给出的搜索方向长什么样。看到第10个示例的时候对SQP应该会有一种原来如此的顿悟感。这个方向上的深水区还有很多——比如稀疏QP、filter方法、非光滑约束处理都可以在这套骨架上继续扩展但那是下一步的话题了。