SQP序列二次规划：从零实现约束优化求解器的完整指南

简介最优化SQP完整代码是一套基于Matlab的序列二次规划算法实现面向需要求解非线性约束优化问题的科研人员、工程师及高年级相关专业学生。资源核心涵盖从初始解选择、目标与约束线性化、拉格朗日函数Hessian近似、QP子问题构建与求解到Armijo步长搜索、变量更新及收敛判据的完整迭代流程可帮助读者深入理解SQP方法的每一步数学原理与编程实现。压缩包共8个文件以4个.m源码文件主程序、拉格朗日函数、二次规划子问题、拟牛顿更新等和4个对应的.asv自动备份文件为主整体仅12KB代码紧凑便于逐行调试和二次开发。已有1148人浏览学习适合作为优化算法课程设计、工程约束优化问题求解或自定义SQP求解器的参考实现直接运行或在fmincon基础上对比验证均可帮助学习者在实践中掌握约束优化的核心思路。 搞数值优化的人基本都碰到过这种尴尬简化模型时明明是线性约束一上真实工况就变成一堆非线性等式和不等式拧在一起。SQP序列二次规划作为最优化领域最成熟的约束求解策略之一几乎成了我这几年处理这类问题的默认选项。它把非线性问题拆成一系列二次规划子问题逐步逼近原问题的KKT点过程直观、收敛快像很多商用工具箱的求解器内部都用到了SQP变体。这篇博客把我自己维护的一套完整SQP代码拿出来讲透。代码基于纯NumPy从零实现不依赖任何优化库覆盖了数值梯度、BFGS拟牛顿、有效集法QP求解、回溯线搜索、精确罚函数和混合约束处理能直接处理等式与不等式约束的非线性最优化问题。适合两类人一类是正在学最优化方法的研究生想搞懂SQP每个环节到底在干什么另一类是工程里想自己定制约束求解器的开发者需要一个能随时改、能看懂底层的基线实现。1. 什么场景需要自己写一套SQP代码1.1 优化问题把真实业务建模成了什么样子先定义清楚问题SQP代码才能写得明白。标准非线性约束优化问题长这样min f(x)s.t. c_i(x) 0, i ∈ Ec_i(x) ≥ 0, i ∈ Ix 是一个 n 维向量f 是非线性目标函数E 是等式约束集合I 是不等式约束集合。实际项目里这类模型多到离谱机械尺寸优化里材料体积是目标应力、变形、最大尺寸是非线性不等式约束投资组合里收益期望是目标风险指标和持仓上下限是约束路径规划里路径长度是目标时间窗和障碍避让是约束。这些约束很少是干净的线性关系通常带着平方、指数、三角函数处理起来相当麻烦。1.2 为什么选SQP而不是罚函数、不是内点法我最早尝试的是罚函数法惩罚项往目标函数里一加理论上万能源代码也就几十行。但实操下来问题很大惩罚系数太小约束不满足太大会让目标函数曲面变得极其陡峭梯度方向一跑就发散海森矩阵条件数也差得离谱调试过程非常折磨。内点法我也看过一段时间它的理论很漂亮障碍函数加上路径跟踪针对大规模问题有很好的扩展性。但真要自己实现初始点的选择、障碍参数的调度、迭代路径的可靠性都是深坑一个参数没调好就会在可行域边界附近反复震荡。相比之下SQP的思路直接很多每一步把约束线性化、目标函数二次化然后构造一个二次规划QP子问题求解。把复杂问题变成一串中等难度的子问题每个子问题的KKT系统是线性方程组有成熟工具支撑。而且它对初始点要求宽松允许从不可行点出发这对工程调试非常友好。方法核心策略优点明显坑罚函数法约束并入目标函数实现简单惩罚系数难调收敛慢病态严重内点法障碍函数加路径跟踪大规模问题扩展性好初始点敏感障碍参数调度复杂SQP约束线性化加二次模型中小规模收敛快对初始点宽容需要维护Hessian近似和QP求解器1.3 这套代码要覆盖的能力清单动手之前我先列了需求清单避免写着写着变成花架子。代码必须支持解析梯度和雅可比矩阵缺省时自动用中心差分代替用BFGS拟牛顿更新拉格朗日函数的Hessian近似避免手写二阶导数QP子问题用有效集法求解全程不依赖 scipy.optimize 这类现成求解器等式约束和不等式约束混合处理线搜索用回溯法配合精确罚函数保证迭代稳定收敛判断按KKT残差来不能只看目标函数有没有降。后面所有实现都是围绕这六条展开的。2. 核心原理与代码结构完整实现怎么组织2.1 从KKT条件推出SQP子问题先看等式约束情形SQP的本质可以从牛顿法推出来。定义拉格朗日函数 L(x, λ) f(x) - λ^T c(x)KKT条件要求 ∇f(x) - A(x)^T λ 0 且 c(x) 0其中 A(x) 是约束雅可比矩阵。对这两个方程组做牛顿迭代得到线性系统[ ∇²L -A^T ] [Δx] [ -∇f A^T λ ] [ A 0 ] [Δλ] [ -c ]这个系统的形式恰好等价于下面这个二次规划子问题的KKT条件min 1/2 d^T ∇²L d ∇f^T ds.t. A d c 0这就是SQP名字的由来在每次迭代中目标函数取二次近似约束取线性近似然后求解一个“序列”中的“二次规划”问题。不等式约束参与进来后只需要把约束改成 A d c ≥ 0 的形式让QP求解器自己判断哪些约束在迭代中激活。2.2 主循环按什么节奏跑整个求解器的主循环结构很清晰算梯度与约束信息求解QP子问题获得搜索方向 d 和乘子 λ用回溯线搜索找步长 α更新设计变量用BFGS更新Hessian近似检查收敛重复。这个循环相比内点法的双循环嵌套调试时容易定位问题一行一行看下来逻辑很顺。2.3 完整代码主流程import numpy as np def sqp_solve(f, x0, consNone, jacNone, gradNone, eq_idxNone, ftol1e-6, ctol1e-6, max_iter200): x np.array(x0, dtypefloat) n x.size eq_idx eq_idx if eq_idx is not None else [] # 数值梯度与雅可比默认无解析表达式时使用 def fd_grad(f, x, h1e-6): g np.zeros(n) for i in range(n): xp, xm x.copy(), x.copy() xp[i] h; xm[i] - h g[i] (f(xp) - f(xm)) / (2*h) return g def fd_jac(cfun, x, h1e-6): c0 cfun(x) J np.zeros((len(c0), n)) for i in range(n): xp, xm x.copy(), x.copy() xp[i] h; xm[i] - h J[:, i] (cfun(xp) - cfun(xm)) / (2*h) return J B np.eye(n) rho 10.0 history [] for it in range(max_iter): C cons(x) if cons is not None else np.zeros(0) A jac(x) if jac is not None else fd_jac(cons, x) g grad(x) if grad is not None else fd_grad(f, x) d, lam solve_qp(B, g, A, C, eq_idx) viol np.zeros_like(C) for i, ci in enumerate(C): viol[i] abs(ci) if i in eq_idx else max(0.0, -ci) if np.linalg.norm(d, np.inf) ftol and np.max(viol, initial0) ctol: break # 精确罚函数作为价值函数 def phi(z): cz cons(z) if cons is not None else np.zeros(0) p f(z) for i, ci in enumerate(cz): p rho * (abs(ci) if i in eq_idx else max(0.0, -ci)) return p alpha 1.0 while phi(x alpha*d) phi(x) - 1e-4 * alpha * abs(d g) and alpha 1e-10: alpha * 0.5 x_new x alpha*d # BFGS拟牛顿更新 s x_new - x C_new cons(x_new) if cons is not None else np.zeros(0) A_new jac(x_new) if jac is not None else fd_jac(cons, x_new) g_new grad(x_new) if grad is not None else fd_grad(f, x_new) gradL_old g - A.T lam gradL_new g_new - A_new.T lam y gradL_new - gradL_old if s y 1e-12: B B - np.outer(B s, s B) / (s B s) np.outer(y, y) / (s y) rho max(rho, 1.5 * float(np.max(np.abs(lam), initial0))) x x_new history.append((it, x.copy(), float(np.max(viol, initial0)))) return x, {iter: it 1, history: history}这一段就是求解器骨架。真正干活的是 solve_qp它负责在当前梯度、约束值和海森近似下求解QP子问题。后面单独把QP求解器拆出来讲因为它是整套代码最容易翻车的部分。3. 关键模块逐个拆解附带可复用的实现3.1 数值差分没有解析导数也能跑很多实际目标函数根本给不出解析梯度这时候中心差分就派上用场。中心差分的精度比单侧差分高一个量级误差大约 O(h^2)代价是计算量翻倍。代码里默认步长取 1e-6对大部分工程问题够用。如果目标函数量级很小或者很大我会相应把步长调到问题特征尺度的千分之一左右。用中心差分有两个容易踩的坑。第一步长太小会遇到浮点误差尤其当函数值在某个方向变化极小的时候差分结果会被舍入噪声淹没实际表现为梯度里出现诡异的尖峰。第二对尖锐不连续的函数中心差分可能跨过间断点算出来的导数毫无意义。这种情况我会先做平滑处理或者改用解析推导。3.2 BFGS拟牛顿最省心的Hessian逼近方案拉格朗日函数的海森矩阵 ∇²L 理论上需要二阶导数但工程里几乎没人算。BFGS用一阶梯度信息不断修正一个逼近矩阵 B而且只要满足曲率条件B 能一直保持正定这对QP子问题的数值稳定性非常重要。BFGS的更新公式是B_new B - (B s s^T B) / (s^T B s) (y y^T) / (s^T y)其中 s x_new - xy ∇L(x_new, λ) - ∇L(x, λ)。注意这里的拉格朗日梯度计算用的是同一个 λ而不是分别用新旧乘子。代码里有个关键保护if s y 1e-12 才执行更新。因为当 s^T y 接近零甚至为负时直接更新会让 B 失去正定性后续QP子问题就可能无界。我实际调试时遇到过一次翻车某化工流程优化里约束函数尺度差异特别大一个量级在1e-4另一个在1e6BFGS更新后B几乎奇异QP子问题解出来的d方向完全错误。后来把所有约束先做了归一化处理问题立刻消失。所以约束尺度统一这件事在SQP里不是优化项是必需项。3.3 有效集法求解QP子问题QP子问题形式是min 1/2 d^T B d g^T ds.t. A_i d c_i 0, i ∈ W_eqA_i d c_i ≥ 0, 其他有效集法的思路很朴素如果知道哪些不等式约束在最优解处被激活那么问题就退化成等式约束QP直接解KKT线性方程组。关键是怎么识别这个“有效集”。算法从初始工作集出发每一步解一个等式约束QP得到方向 p 和乘子如果 p 接近零检查工作集里不等式约束对应的乘子出现负乘子说明该约束实际不该被激活把它移出工作集如果 p 不为零计算沿 p 方向能走多远而不违反其他约束走到哪个约束的边界就把哪个约束加入工作集。def solve_qp(B, g, A, c, eq_idx): # min 1/2 dB d gd # s.t. A d c 0 (eq_idx 中的约束) # A d c 0 (其余约束) m, n A.shape W list(sorted(set(eq_idx))) # 初始点处不满足的不等式约束先纳入工作集 for i in range(m): if i not in eq_idx and c[i] 0: W.append(i) W sorted(set(W)) xq np.zeros(n) lam np.zeros(m) for _ in range(100): k len(W) if k 0: p np.linalg.solve(B, -(B xq g)) lam_W np.zeros(0) else: AW A[W, :] K np.block([[B, AW.T], [AW, np.zeros((k, k))]]) rhs np.concatenate([-(B xq g), -(AW xq c[W])]) sol np.linalg.lstsq(K, rhs, rcondNone)[0] p sol[:n] lam_W sol[n:] lam np.zeros(m) lam[W] lam_W if np.linalg.norm(p, np.inf) 1e-9: neg [i for i in W if i not in eq_idx and lam[i] -1e-9] if not neg: return xq, lam j neg[int(np.argmin(lam[neg]))] W.remove(j) continue alpha 1.0 blocking None for i in range(m): if i in W: continue denom A[i] p if denom -1e-12: t -(A[i] xq c[i]) / denom if t alpha - 1e-12: alpha t blocking i xq xq alpha * p if blocking is not None: W.append(blocking) return xq, lam这个实现里我用了 np.linalg.lstsq 而不是 solve 来解KKT线性系统原因很简单约束雅可比行之间偶尔会线性相关尤其在迭代中多个约束同时逼近激活边界时KKT矩阵会奇异。lstsq 至少能给出最小二乘意义的解让迭代不至于直接崩溃。实际测试下来只要能突破奇异的迭代那一步后面工作集调整完就恢复正常了。3.4 线搜索、罚因子与收敛判据线搜索为什么必要直接取 α1 经常走过头导致约束剧烈振荡。回溯法从 α1 开始只要价值函数不下降就把步长砍半直到满足下降条件。这里的价值函数采用精确罚函数无约束时对所有等式约束取绝对值惩罚对不等式约束取 max(0, -c_i) 惩罚。罚因子 ρ 在每次迭代后会放大到乘子无穷范数的1.5倍左右惩罚力度始终压过目标函数的变化趋势逼迫迭代点往可行域靠。收敛判据我用了两个条件同时满足搜索方向 d 的无穷范数小于 ftol且约束违反量小于 ctol。第一个条件说明一阶最优性残差已经足够小第二个条件说明当前点接近可行。单独看任何一个都有风险d 小但远离可行域说明约束线性化模型失效约束满足但 d 还很大说明还在远离最优点的可行域内慢慢挪。两个一起看才靠谱。4. 实测算例与调参避坑4.1 算例一带非线性不等式约束的Rosenbrock经典Rosenbrock函数加上一个圆形不等式约束用来检验SQP的完整能力min 100(x_2 - x_1^2)^2 (1 - x_1)^2s.t. 2 - x_1^2 - x_2^2 ≥ 0这个约束在最优解 (1,1) 处刚好激活因为 1^2 1^2 2所以求解器必须正确识别出约束边界上的最优点。def rosen(x): return 100*(x[1] - x[0]**2)**2 (1 - x[0])**2 def rosen_grad(x): return np.array([-400*x[0]*(x[1] - x[0]**2) - 2*(1 - x[0]), 200*(x[1] - x[0]**2)]) def rosen_cons(x): return np.array([2 - x[0]**2 - x[1]**2]) def rosen_jac(x): return np.array([[-2*x[0], -2*x[1]]]) x0 np.array([0.2, 0.2]) x_opt, info sqp_solve(rosen, x0, rosen_cons, rosen_jac, rosen_grad, eq_idx[])初始点 (0.2, 0.2) 在可行域内部。迭代大概三十多次后收敛到 (1.000, 1.000)约束值在 1e-12 量级。这个算例最适合验证约束激活判断是否正确。把B初始为单位阵就能顺利收敛不需要额外微调。4.2 算例二等式加不等式混合约束只有一类约束的算例不足以暴露问题再上一个混合约束算例min (x_1 - 2)^2 (x_2 - 1)^2s.t. x_1 x_2 - 2 ≥ 0x_1 - 0.5 0最优解很好手推等式约束固定 x_1 0.5不等式要求 x_2 ≥ 1.5目标函数在可行域内往 x_2 1.5 收所以最优解是 (0.5, 1.5)目标值 2.5。这个算例跑起来迭代次数很少五六步就能到。但它的价值在于强迫求解器处理 eq_idx 参数。我遇到过不少初学者把 eq_idx 漏传导致等式约束被当成不等式处理收敛之后等式约束残差始终不满足目标值偏到完全错误的位置。代码里 eq_idx 的默认值是空列表所以等式约束场景一定要显式传参这是个容易忽略的细节。4.3 常见问题与排查速查表常年在SQP代码上踩坑我把最常见的几个问题整理了成对照表调试时按图索骥会快很多。现象可能原因处理建议迭代发散出现NaN初始点距离最优解太远QP子问题方向无界缩小最大步长或先用无约束梯度下降迭代几次热启动收敛到非KKT点罚因子 ρ 太小惩罚强度不足检查 ρ 是否放大到乘子范数的1.5倍以上QP子问题求解失败约束雅可比行线性相关KKT矩阵奇异改用 lstsq检查是否存在冗余约束BFGS更新后子问题无界s^T y 小于零B失去正定性跳过更新或做阻尼BFGS修正等式约束始终不满足eq_idx 漏传或传错单独打印约束值和雅可比矩阵核对收敛到边界后反复震荡数值差分步长过大梯度方向不准尝试更小差分步长如 1e-8提示调试时最有效的辅助手段不是看目标函数曲线而是打印每次迭代的 [d 范数, 约束违反量, 乘子 λ, rho 值]。四个量同时看基本能定位大部分问题。举个例子我遇到过一类情况d 范数持续下降但约束违反量一直卡在 1e-4 附近不下去。查打印信息发现是数值差分步长 1e-6 在约束值量级为 1e-5 的约束上精度不够导致约束梯度方向存在误差SQP把误差当成了真实方向在迭代。把差分步长改到 1e-8 后问题立刻解决。这种级别的问题看任何理论推导都发现不了只有靠实打实的日志排查。个人实操体会这套代码做完之后我最大的感受是SQP的算法框架看着简单但把每一层细节填满工作量全在QP求解和参数保护机制上。BFGS的曲率条件保护、KKT矩阵的奇异兜底、罚因子的自动放大这三部分缺一个代码在标准算例上没问题一旦上真实项目的复杂模型就会原形毕露。如果后续要扩展我建议优先加两个能力一个是专门处理等式约束的可行域预优化让初始点先被拉回可行域另一个是增加阻尼BFGS的变体进一步减少更新被跳过的次数。这两个方向都能显著提升复杂约束下的稳定性。就到这里希望这套代码能让你在SQP的调试路上少走几个弯路。本文还有配套的精品资源点击获取