SQP序列二次规划法：非线性约束优化问题的原理与Python实践

项目标题用序列二次规划法SQP解决非线性优化问题 相关热词SQP、序列二次规划法、非线性优化问题 摘要描述从数学原理到工程实践拆解序列二次规划法SQP如何求解带约束的非线性优化问题并给出可复现的Python实现与调试经验。真正让我下决心把SQPSequential Quadratic Programming序列二次规划法彻底搞明白是因为几年前做结构参数优化时遇到一个很现实的问题目标函数是非线性的约束条件里有非线性不等式还夹杂着耦合的等式约束。当时图省事直接用罚函数法结果罚因子调了一个星期要么约束严重违约要么迭代收敛慢到怀疑人生换成遗传算法吧收敛速度又太拉胯每次结果还在最优解附近抖来抖去。后来老老实实把SQP的原理和实现啃了一遍才发现这类带约束的非线性优化问题本来就有非常成熟的解法。这篇文章不是我复述教科书而是把从学SQP、写SQP到调SQP整个过程的经验整理出来。内容会覆盖SQP的数学原理、QP子问题的构造逻辑、Python实现以及大量平时踩坑才换来的调试经验。无论你是在做轨迹规划、结构优化、经济调度还是单纯被课程作业卡住了这篇文章应该都能帮到你。我会尽量讲得通俗一点让没有系统学过优化的开发者也能跟上节奏。1. 非线性优化问题为什么难先搞清楚我们在面对什么1.1 一个典型的带约束非线性问题长什么样先看一个大家都会遇到的数学形式min f(x) s.t. g_i(x) 0, i 1..m h_j(x) 0, j 1..p这里的x是决策变量向量f(x)是目标函数g_i(x) 0是不等式约束h_j(x) 0是等式约束。所谓“非线性优化问题”就是f、g、h三个函数里面至少有一个不是线性的。注意这里说的非线性是广义的比如目标函数带平方项、约束里是半径的平方和、或者函数里有指数、对数、三角函数都算。举个例子假设我要在平面上找一个点让它到点(3,4)的距离最短同时这个点必须落在以原点为圆心、半径为3的圆内还要求它的横纵坐标之和不小于2。写成数学表达式就是min (x1-3)^2 (x2-4)^2 s.t. x1^2 x2^2 9 x1 x2 2这就是一个非常典型的带约束非线性优化问题。看起来人畜无害但让程序自动求解事情就没那么简单了。首先目标函数的最优解在(3,4)可这个点根本不在可行域里直接做无约束优化必然会跑到约束外面去。其次约束边界是一个圆弧不是一条直线这导致可行域的边界是弯曲的很多线性规划的方法根本没法用。1.2 无约束方法直接杀过去会发生什么很多人遇到优化问题第一反应就是梯度下降。梯度下降的思路很简单沿着负梯度方向走直到梯度为0或者目标函数不再下降。问题是梯度下降完全没有约束的概念它只关心目标函数怎么走下降最快根本不理会你要求的“必须在圆内”“必须满足某个等式”这些边界条件。如果你硬把约束加到目标函数里比如用罚函数法把违反约束的程度加进去min f(x) ρ * [Σ max(0, g_i(x))^2 Σ h_j(x)^2]这样确实能让解逼近可行域但麻烦接踵而至。罚因子ρ太小最终解可能严重违反约束ρ太大目标函数和罚项的数量级差出好几个零导致Hessian矩阵病态普通优化算法根本走不动。我当年调罚因子调到怀疑人生就是这个原因。乘子法augmented Lagrangian比纯罚函数好一些但依然存在每次迭代都要处理不等式约束的active set识别问题实现复杂度不低。说白了无约束优化方法处理带约束问题的思路都是“绕路”把一个约束问题强行变成无约束问题。而SQP的思路完全不同正面硬刚在每一次迭代中把原问题近似成一个带约束的二次规划问题然后去解这个更简单的子问题。1.3 SQP的核心思路每次迭代解一个局部二次规划SQP的基本想法一句话就能讲清楚在当前迭代点附近把目标函数做二阶泰勒展开把约束函数做一阶泰勒展开即线性化得到一个带线性约束的二次规划Quadratic Programming, QP子问题然后求解这个QP子问题得到搜索方向再通过线搜索确定步长迭代更新直到收敛。为什么非要近似成二次规划因为带线性约束的二次规划已经有了极其成熟的求解方法无论是有效集active set法还是内点法都能快速稳定地求解。而且二次模型比一阶线性模型更接近光滑非线性函数的真实形状配合拟牛顿法对二阶信息的逼近SQP理论上可以达到超线性甚至二次收敛速度。我用一个类比帮你建立直觉你在爬一座山头但山下画了禁行区域。SQP的做法是——每走一步都站在原地感受一下脚下的斜坡曲率搭一个局部“滑梯”二次模型再结合禁行区域的边界方向线性化约束算出下一步往哪滑、滑多远。走几步重新搭一次滑梯最终就能到达禁行区域内的最低点。这个“边走边搭滑梯”的策略就是序列二次规划法名字里“序列”二字的由来。2. SQP原理拆解从KKT条件到QP子问题2.1 Lagrange函数与KKT条件约束优化的路标要理解SQP必须先理解KKTKarush-Kuhn-Tucker条件。KKT条件是非线性约束优化问题的一阶必要条件它告诉我们如果一个点是局部最优解那么在这个点上目标函数的梯度、约束函数的梯度以及拉格朗日乘子之间必须满足一组平衡关系。先写出原问题的Lagrange函数L(x, λ, μ) f(x) Σ λ_i * g_i(x) Σ μ_j * h_j(x)其中λ_i 0是不等式约束对应的拉格朗日乘子μ_j是等式约束对应的乘子可正可负。KKT条件包含四部分梯度平衡∇_x L(x, λ, μ) 0即目标函数梯度与约束梯度的加权和互相抵消原始可行g_i(x) 0h_j(x) 0即解必须在可行域内对偶可行λ_i 0互补松弛λ_i * g_i(x) 0即不等式约束要么不起作用g_i(x) 0此时乘子为0要么正压在边界上g_i(x) 0乘子可以大于0。严格说KKT条件需要在约束规范如LICQ条件成立时才是必要条件但工程里大多数光滑问题都满足。SQP的终止判据就是让迭代点逐步逼近KKT条件。你可以把KKT条件想象成停车场里的划线车位每辆车都必须同时满足“不能越线”“不能压线”“方向和坡度都合适”这三个要求才算停到位。2.2 在迭代点构造二次规划子问题假设现在已经有了一个迭代点x_k我们要在这里构造QP子问题。核心思路是在x_k附近做泰勒展开目标函数保留到二阶项f(x) ≈ f(x_k) ∇f(x_k)^T d (1/2) * d^T * B_k * d约束函数保留到一阶项g_i(x) ≈ g_i(x_k) ∇g_i(x_k)^T dh_j(x) ≈ h_j(x_k) ∇h_j(x_k)^T d这里的d是决策变量在当前迭代点的位移即x x_k d。B_k是Lagrange函数关于x的Hessian矩阵或它的近似。替换到原问题后就得到QP子问题min_d ∇f(x_k)^T d (1/2) * d^T * B_k * d s.t. g_i(x_k) ∇g_i(x_k)^T d 0 h_j(x_k) ∇h_j(x_k)^T d 0注意目标函数里的常数项f(x_k)在求解最小值时不影响d的选择所以通常省略但如果是用来计算函数值别忘了加回来。为什么不把约束也展开到二阶因为如果约束也做二阶近似子问题就变成了带二次约束的二次规划QCQP求解难度立刻上升一个量级。一阶线性化虽然听起来粗糙但SQP依靠迭代更新每轮都会在fresh点重新线性化最终收敛时误差会趋于0这正是SQP能保持超线性收敛效率的关键原因。2.3 用BFGS近似Hessian不要傻傻算二阶导QP子问题里那个B_k到底怎么来理论上B_k应该是Lagrange函数对x的Hessian矩阵也就是要计算所有二阶偏导数。但实际工程中二阶导数的解析式往往非常复杂数值微分又慢又不稳定。所以SQP的工程实现基本都用拟牛顿法来近似Hessian其中BFGSBroyden-Fletcher-Goldfarb-Shanno公式是最常用的。BFGS的更新公式长这样s_k x_{k1} - x_k y_k ∇_x L(x_{k1}, λ_{k1}, μ_{k1}) - ∇_x L(x_k, λ_{k1}, μ_{k1}) B_{k1} B_k - (B_k * s_k * s_k^T * B_k) / (s_k^T * B_k * s_k) (y_k * y_k^T) / (y_k^T * s_k)这个公式看起来很吓人但它的本质其实很简单根据相邻两次迭代中梯度信息的变化去修正当前对二阶曲率的估计。如果实际函数是二次的BFGS经过有限次迭代就能精确还原Hessian对非二次函数它也能给出一个合理的正定近似。这里有几个非常关键的经验点注意一是初始B_0通常设为单位矩阵I随着迭代逐步修正二是必须保证更新后B_{k1}正定。实际代码中通常还要加一个阻尼处理damped BFGS避免y_k^T s_k太小甚至为负导致更新破坏正定性。如果B_k不是正定的QP子问题本身可能变成非凸二次规划求解难度陡增。这也是为什么很多SQP实现会在y^T s过小的时候跳过这次更新或者改用单位矩阵做混合。2.4 线搜索与L1精确罚函数保证每一步都在前进有了搜索方向d接下来还要确定走多远。这里的细节往往决定算法是否稳定如果不做线搜索直接x_{k1} x_k d可能一步就冲出可行域甚至导致目标函数值上升。SQP常用的做法是引入一个“评价函数”merit function将目标函数和约束违反程度统一成一个标量然后在搜索方向上进行一维线搜索。最常见的评价函数是L1精确罚函数φ(x; ρ) f(x) ρ * [Σ max(0, g_i(x)) Σ |h_j(x)|]其中ρ是罚参数。为什么叫“精确”因为只要ρ足够大这个无约束问题的极小点就和原带约束问题的极小点一致不需要让ρ趋向无穷大——这跟罚函数法完全不一样正是它好用之处。然后用Armijo回溯准则来确定步长αφ(x_k α*d_k; ρ) φ(x_k; ρ) - c * α * ||d_k||^2c通常取1e-4这样的小值。如果步长为1时条件不满足就把步长乘0.5继续试。我实测下来这个回溯法虽然有点“笨”但非常稳定很适合工程场景。另一种思路是信赖域trust region方法给子问题加上一个对||d||的半径约束也能保证全局收敛。两种策略各有优劣SciPy的SLSQP其实更接近“序列最小二乘规划”的变体但整体框架依然是SQP那一套。3. 用Python实现SQP核心代码与验证3.1 先讲一个关键点QP子问题谁来解写SQP代码前你要先想明白一件事外层SQP负责组织迭代但每一步里面的QP子问题到底由谁来解如果你自己实现一个完整的QP求解器比如带有效集法至少还要加两百行代码。这显然不是一篇博文能承载的重点。所以我的建议是直接调用成熟的二次规划求解器。在Python生态里scipy.optimize.minimize的SLSQP方法本身就能处理带约束的优化但它更像一个完整的SQP实现如果你想手动搭建外层SQP框架底层QP子问题可以使用qpsolvers库、cvxopt或者顺手用scipy.optimize.minimize去解一个线性约束二次目标的小问题虽然这有点“套娃”但作为教学演示完全够用。先别急着追求工业级写法重点是看清SQP的外层骨架。3.2 快速上手直接用SciPy的SLSQP解非线性约束问题先说最省事的路径。如果你不想手动实现SQP直接用scipy.optimize.minimize(methodSLSQP)就够了。下面这段代码求解的就是前面提到的圆内距离问题import numpy as np from scipy.optimize import minimize def objective(x): return (x[0] - 3.0)**2 (x[1] - 4.0)**2 def obj_grad(x): return np.array([2.0 * (x[0] - 3.0), 2.0 * (x[1] - 4.0)]) def circle_constraint(x): return 9.0 - x[0]**2 - x[1]**2 # SciPy要求 0 def circle_grad(x): return np.array([-2.0 * x[0], -2.0 * x[1]]) def line_constraint(x): return x[0] x[1] - 2.0 # SciPy要求 0 def line_grad(x): return np.array([1.0, 1.0]) cons [ {type: ineq, fun: circle_constraint, jac: circle_grad}, {type: ineq, fun: line_constraint, jac: line_grad}, ] x0 np.array([1.0, 1.0]) # 选择一个可行的初始点 res minimize(objective, x0, methodSLSQP, jacobj_grad, constraintscons, options{ftol: 1e-12, maxiter: 200}) print(res.x) # 预期约 [1.8, 2.4] print(res.fun) # 预期约 4.0 print(res.success, res.message)运行结果应当是x ≈ [1.8, 2.4]目标函数值4.0successTrue。你可以验证一下点(1.8, 2.4)在半径3的圆边界上到(3,4)的距离平方是(1.2)^2 (1.6)^2 4同时x1 x2 4.2 2完全满足约束。这里有一个极其常见的坑我必须重点提醒注意SciPy里ineq约束的fun(x)返回值必须大于等于0即fun(x) 0表示约束成立。如果你写成x[0]**2 x[1]**2 - 9 0这种直接照抄数学式的形式算法会把约束方向理解反结果要么报错要么解出一个完全错误的值。我见过太多人在这里翻车了。3.3 手写一个简化版SQP看清外层骨架为了让你真正理解SQP是怎么迭代的我再写一个教学用的手写版。这个版本只展示外层SQP框架构造QP子问题、线搜索、BFGS更新。核心逻辑是完整的但子问题的求解我直接调用了scipy.optimize.minimize来处理这样代码才能控制在一屏之内。我用一个线性约束的例子来验证这样简化版中的BFGS更新逻辑完全严格约束x1 x2 2无约束最优解(3,4)违反约束真正的最优解应该落在直线x1 x2 2上。import numpy as np from scipy.optimize import minimize def f(x): return (x[0] - 3.0)**2 (x[1] - 4.0)**2 def gf(x): return np.array([2.0 * (x[0] - 3.0), 2.0 * (x[1] - 4.0)]) # 约束表示g(x) 0这里 g x1 x2 - 2 def g(x): return x[0] x[1] - 2.0 def gg(x): return np.array([1.0, 1.0]) def solve_qp(xk, Bk): 解线性约束的二次规划子问题返回搜索方向 d def qp_obj(d): return gf(xk) d 0.5 * d Bk d cons [ {type: ineq, fun: lambda d, xkxk: -(g(xk) gg(xk) d)} ] res minimize(qp_obj, np.zeros(2), methodSLSQP, constraintscons) return res.x def l1_merit(xk, rho): return f(xk) rho * max(0.0, g(xk)) def sqp_basic(x0, max_iter100, tol1e-8): x np.array(x0, dtypefloat) n len(x) B np.eye(n) rho 10.0 for _ in range(max_iter): d solve_qp(x, B) if np.linalg.norm(d) tol: break alpha 1.0 phi0 l1_merit(x, rho) while l1_merit(x alpha * d, rho) phi0 - 1e-4 * alpha * np.linalg.norm(d)**2: alpha * 0.5 s alpha * d x_new x s # 约束梯度为常数直接用 f 的梯度差做BFGS在这种场景下严格成立 y gf(x_new) - gf(x) if y s 1e-12: B B - np.outer(B s, s B) / (s B s) np.outer(y, y) / (y s) x x_new print(fiter{_:3d} x{x} f{f(x):.6f}) return x x_opt sqp_basic(np.array([0.0, 0.0])) print(opt:, x_opt, f:, f(x_opt))跑完以后你应该能看到迭代过程稳定收敛到x ≈ [0.5, 1.5]目标函数值约12.5。这个例子的解析解也可以用拉格朗日乘子法手推出来在直线x1 x2 2上最小化(x1-3)^2 (x2-4)^2代入x2 2 - x1后求导可得x1 0.5x2 1.5。这段代码我特意做了简化有几处要向你说明白第一QP子问题我用scipy的SLSQP去解虽然递归调用了SQP实现但因为子问题的目标是二次、约束是线性SLSQP处理这种子问题毫无压力第二BFGS更新只用了f的梯度差这是因为本例约束梯度是常数完整的SQP应该用Lagrange函数的梯度差第三罚参数rho我是手工固定为10生产环境应该根据约束违背程度动态调整。简单说这个版本是“教学骨架”不是“工业实现”真上生产就用3.2节的代码或者成熟优化库。3.4 生产环境怎么选手工实现还是调用现成库看完上面的代码你应该已经能分清SQP的两件事外层迭代逻辑线性化、二次近似、线搜索和内层QP求解。实际工程开发中绝大多数情况下你不需要手动实现任何一个环节。SciPy的SLSQP、trust-constr以及专业库IPOPT、NLopt都提供了经过大量测试的SQP或内点法实现性能和稳定性远超我们手写的版本。那为什么还要理解SQP原理我的体会是当你面对一个具体优化问题调参调到崩溃、算法给出错误解的时候只有理解了原理你才能判断问题到底出在“约束方向写反”还是“Hessian近似失效”还是“初始点不可行”。优化算法不是黑盒你不能指望不掌握原理就能用好它。4. SQP能解决哪些实际问题工程应用与算法选型4.1 典型应用场景盘点SQP在工程界的应用范围非常广我挑几个典型领域列出来领域典型优化问题为什么适合SQP航空航天飞行器轨迹优化、火箭入轨方案状态变量与控制变量强耦合动力学约束高度非线性机械结构轻量化设计应力、位移约束有限元响应与设计参数之间是非线性隐式关系化工过程反应器操作条件优化、物料平衡化学反应平衡方程复杂常伴随非线性等式约束电力系统经济调度、最优潮流计算电网功率方程是典型非线性潮流约束金融量化带约束的组合优化、风险预算配置目标函数常为二次风险项约束是线性/凸约束在这些场景里SQP之所以受欢迎是因为它能够同时处理等式约束和不等式约束对初值有一定容忍度收敛速度在中小规模问题上非常可观。特别是当你的问题有几十个决策变量、十几个非线性约束时SQP往往能在几十次迭代内给出满足要求的解。4.2 一个带约束的参数拟合实战案例说一个我做过的例子。当时需要标定某个材料本构模型的参数要求拟合曲线不仅总体误差小还要保证某些物理量不越界。目标函数是模型输出与实验数据的残差平方和约束条件是仿真得到的应力峰值和应变范围不能超阈值。这类问题的共同点是目标函数和约束函数都要通过一个数值模拟器计算梯度都得靠前向差分或者伴随方法。用SciPy的SLSQP建模时要注意的关键接口是约束函数的定义约束条件是“应力峰值 - 阈值 0”在代码里要写成threshold - stress_peak(x) 0。我经常看到有人把不等号方向写反导致算法一开始就认为约束已经满足解出来以后邦邦响——全是坑。对这种黑箱仿真模型我强烈建议你尽量提供解析梯度或半解析梯度至少也要用中心差分代替默认的前向差分。因为这个问题的目标函数和约束函数本身就有数值噪声再用大差分离散步长去算梯度SQP相当于闭着眼睛爬山很容易震荡。4.3 SQP、内点法、罚函数法、遗传算法怎么选每次我讲SQP都有人问那内点法呢遗传算法呢我到底该用哪个这里我把主流几类算法放在一张表里对比算法是否需要导数约束处理收敛速度全局搜索能力适用规模SQP需要线性化QP子问题超线性弱落入局部最优中小规模几十到几百变量内点法需要障碍项/路径跟踪超线性~二次弱但大尺度稳定大规模上千变量罚函数法需要罚项并入目标较慢弱约束简单的小问题遗传算法不需要罚/修复/解码慢较强非光滑、离散、黑箱我的个人选择逻辑很简单如果问题光滑、变量规模在几百以内、约束数量不多优先试SQP规模很大、约束很多直接上内点法IPOPT目标和约束有大量离散变量或不可导别再纠结SQP老老实实选择遗传算法这类启发式算法。SQP不是万能的但它是很多工程问题最值得先试的解法。5. 调试SQP的常见坑与排查技巧实录5.1 约束方向写反最普遍的低级错误这个坑真的值得反复强调SciPy里ineq约束的fun(x)必须返回0的值而eq约束的fun(x)必须返回0。也就是说如果你的约束是g(x) 0建模时要写成-g(x) 0如果是h(x) 0直接写成h(x)。我曾经在某个项目里定义了一个不等式约束x[0]**2 x[1]**2 - 9 0程序跑完successTrue结果解出来的点在圆外十几公里。排查了两天才发现是约束方向写反算法早就把可行域理解成“圆外区域”了。所以拿到结果第一步不要信success先手动检查约束是否真的满足这是一个好习惯。5.2 初始点不可行SLSQP直接摆烂怎么办SLSQP的初始点既可以在可行域内也可以在可行域外它一般能自己拉回到可行域。但有时会遇到Inequality constraints incompatible或Positive directional derivative for linesearch这类报错前者通常是初始点离可行域太远或者约束之间本身就互相冲突后者往往是当前点在某个约束边界上搜索方向却又指向了约束外部。我的经验是能给出可行初始点就尽量给可行初始点这个比什么参数都重要。如果求一个可行点本身就很难可以先放弃部分约束先做一次无约束或小罚因子的优化得到一个比较合理的点再作为初始点或者用连续化方法把约束从宽松版本逐步收紧。另一个技巧是把options里的maxiter调大一些比如从默认100改到1000因为SLSQP在前几次迭代可能主要是在往可行域里爬真正优化目标函数反而在后期。5.3 数值梯度把你的SQP带沟里如果你不给jacscipy.optimize.minimize内部会用有限差分自动算梯度。这个方法在实验室小规模问题上还凑合但工程问题里一旦目标函数带噪声、或者变量尺度差异很大前向差分的步长就很容易踩雷步长太小舍入误差爆发步长太大截断误差大。我建议只要函数能写出解析导数就手写梯度传进去。实在写不出来用中心差分而且要自定义一个合理的步长比如eps 1e-6 * (1 abs(x))。另外我习惯在正式跑优化前先随机抽几个点用解析梯度和数值梯度对一下最大误差超过1e-4就说明梯度写错了别急着进下一步。5.4 不收敛、震荡先查变量的尺度和收敛容差SQP对变量尺度很敏感。假设你的变量一个是x1在0.001量级一个是x2在10000量级目标函数几乎只由x2主导QP子问题里Hessian矩阵的条件数会非常大收敛很慢甚至震荡。解决办法很简单对变量做归一化让所有决策变量都在0~1或者-1~1附近优化完成后再映射回原始尺度。这一步对SQP效果提升非常明显我几乎每次都做。收敛容差也要按实际问题来调。ftol1e-12是我常用的值但也别无脑设置太小目标函数本身噪声水平可能在1e-4你却要求迭代到1e-12那算法会卡在噪声里反复挣扎。先跑一遍看目标函数下降曲线再根据自己的精度需求设一个合理的容差。5.5 多初始点试算对付非凸问题的土办法SQP本质上是一个局部优化算法它找到的最优解极大依赖于初始点。对非凸问题一个初始点可能收敛到A局部最优另一个初始点收敛到B局部最优。很多人拿SQP跑一次就下结论“解不对”其实只是初值没选对。我的土办法是用随机采样生成多组初始点比如20组分别跑SQP最后取目标函数值最小的解。如果目标函数和约束都足够光滑这个“多起点SQP”方法在中小规模问题上效果出奇地好经常能逼近全局最优解。更讲究一点的做法是先用遗传算法粗略搜索一遍把得到的较好解作为SQP的初始点做精修两个阶段配合速度和精度都能兼顾。一些个人体会我自己的体会是SQP是一把非常趁手的瑞士军刀但它不是万能的。如果问题规模特别大、约束特别多内点法可能更稳如果目标函数和约束都不光滑就别死磕SQP老老实实先用启发式算法找好初始点或者对函数做光滑近似。但回到最开头那个场景——中小规模的带约束非线性优化问题SQP绝对是我工具箱里最先拿出来用的方法之一。最后再分享一个小建议调试SQP的时候先别急着怀疑算法按顺序检查三件事——约束方向有没有写反梯度算得对不对变量量纲有没有归一化。我踩过的坑里九成以上都能归到这三类。搞定了这些SQP通常都能给你一个满意的答案。