二次序列规划（SQP）是求解带约束非线性优化问题的高效迭代算法。文章解析了 SQP 的核心思想，包括基于 KKT 条件的搜索方向构建、QP 子问题建模及 Hessian 矩阵近似技术（如 BFGS、L-BFGS）。内容涵盖初始点选择策略、收敛性判据、数值稳定性增强技巧，并结合模型预测控制、结构拓扑优化等场景展示实际应用。通过实例代码演示了如何利用自动微分和可行性恢复方法提升优化效率与精度。
二次序列规划（SQP）是一种用于求解带约束非线性优化问题的高效迭代算法，广泛应用于工程、经济、物理和科学计算等领域。该方法通过在每一步迭代中构建目标函数的二次近似和约束的一阶线性化，将原问题分解为一系列二次规划子问题，逐步逼近最优解。本文深入解析 SQP 的核心思想、关键步骤及 Hessian 矩阵近似技术，并介绍 BFGS、L-BFGS 等常用改进策略，结合实例展示其实际应用效果。
在处理非线性优化问题时，目标函数像一座起伏不定的山脉，约束条件如同密布的围栏。序列二次规划（Sequential Quadratic Programming, SQP）通过每步重新绘制局部地图，用数学方式确定最优搜索方向，避免陷入局部极小值。
我们面对的问题通常形式如下：
$$ \min_x f(x) \quad \text{s.t.} \quad c_i(x) = 0,\ c_j(x) \geq 0 $$
一旦 $f(x)$ 或 $c(x)$ 是非线性的，问题变得棘手。这类问题广泛存在于机器学习超参调优、机器人轨迹规划、金融投资组合构建等领域。它们的共同特点是：可行域可能非凸、存在多个局部极小值、且边界曲率大。传统方法如梯度下降或牛顿法，在这些场景下要么收敛慢，要么根本找不到好解。
SQP 的核心思想是'化整为零'：在脚下画一个小圆圈，假设这片区域是平滑的，判断朝哪个方向下坡最快，然后迈出一步，重复这个过程。即局部近似 + 迭代逼近。
在无约束优化中，我们只需要让梯度为零即可。但在有约束的情况下，一个点要想成为候选最优解，必须满足更严格的 Karush-Kuhn-Tucker（KKT）条件。
| 条件类型 | 数学表达 | 物理意义 |
|---|---|---|
| 平稳性 | $\nabla_x \mathcal{L}(x,\lambda,\mu)=0$ | 梯度力与约束反作用力平衡 |
| 原始可行性 | $c_i(x)=0$, $c_j(x)\geq0$ | 解必须位于可行域内 |
| 对偶可行性 | $\mu_j \geq 0$ | 不等式约束乘子非负 |
| 互补松弛 | $\mu_j c_j(x) = 0$ | 要么约束不起作用，要么紧致 |
其中，拉格朗日函数定义为： $$ \mathcal{L}(x, \lambda, \mu) = f(x) + \sum_{i \in \mathcal{E}} \lambda_i c_i(x) - \sum_{j \in \mathcal{I}} \mu_j c_j(x) $$
这里的 $\lambda_i$ 和 $\mu_j$ 就像是'影子价格'。
import numpy as np
def lagrangian_gradient(x, lambdas, mus, f_grad, c_eq_grads, c_ineq_grads):
"""
计算拉格朗日函数在某点的梯度
参数说明：
x: 当前变量向量 (n,)
lambdas: 等式约束乘子向量 (m_eq,)
mus: 不等式约束乘子向量 (m_ineq,)
f_grad: 目标函数梯度函数 handle -> R^n
c_eq_grads: 等式约束梯度列表 [R^n]
c_ineq_grads: 不等式约束梯度列表 [R^n]
"""
grad_L = f_grad(x)
for i, g in enumerate(c_eq_grads):
grad_L += lambdas[i] * g(x)
for j, h in enumerate(c_ineq_grads):
grad_L -= mus[j] * h(x)
return grad_L
这段代码检查当前点是否接近 KKT 点。如果 lagrangian_gradient 接近零向量，说明系统达到了力学上的平衡状态。
具体做法是在当前点 $x_k$ 处构造一个 QP 子问题：
$$ \begin{aligned} \min_d &\quad \nabla f(x_k)^T d + \frac{1}{2} d^T B_k d \ \text{s.t.}&\quad c_i(x_k) + \nabla c_i(x_k)^T d = 0 \ &\quad c_j(x_k) + \nabla c_j(x_k)^T d \geq 0 \end{aligned} $$
这里面藏着三个关键思想：
整个流程可以用一张图概括：
flowchart TB
A[初始化 x₀, λ₀, μ₀, B₀] --> B{满足终止条件?}
B -- 否 --> C[构建 QP 子问题]
C --> D[求解 QP 得到 dₖ, νₖ]
D --> E[线搜索求 αₖ]
E --> F[更新 xₖ₊₁ = xₖ + αₖdₖ]
F --> G[更新 λₖ₊₁, μₖ₊₁ ← νₖ]
G --> H[更新 Bₖ→Bₖ₊₁ via BFGS]
H --> B
B -- 是 --> I[输出最优解]
这是一个闭环反馈系统，QP 子问题的对偶变量 $\nu_k$ 直接被拿来更新拉格朗日乘子。
考虑这个问题： $$ \min x_1^2 + x_2^2 \quad \text{s.t.}\ x_1^2 + x_2^2 \leq 1 $$ 最优解显然是原点 $(0,0)$。但如果选初始点为边界上的 $(1,0)$，可能会因为忽略圆的曲率导致不可行。
| 初始点位置 | 可行性 | 局部近似误差 | 典型风险 |
|---|---|---|---|
| 内部中心 | 高 | 低 | 收敛快，稳定性好 ✅ |
| 边界平坦区 | 中 | 中 | 步长受限，易振荡 ⚠️ |
| 边界曲率大区 | 低 | 高 | 快速失可行，发散风险 ❌ |
| 外部轻微违规 | 视情况 | 高 | 子问题不可行 🔁 |
如果初始猜测严重违反约束，先做一个'可行性恢复阶段'： $$ \min_{x,s} \sum_i s_i \quad \text{s.t.}\ g_i(x) \leq s_i,\ h_j(x)=0,\ s_i \geq 0 $$
import numpy as np
from scipy.optimize import minimize
def feasibility_restoration(g_funcs, h_funcs, x0):
n = len(x0)
m = len(g_funcs)
p = len(h_funcs)
def objective(z):
x = z[:n]
s = z[n:n+m]
return np.sum(s)
def constraint_eq(z):
x = z[:n]
return [h(x) for h in h_funcs]
def constraint_ineq(z):
x = z[:n]
s = z[n:n+m]
return [s[i] - g_funcs[i](x) for i in range(m)]
z0 = np.concatenate([x0, np.ones(m)])
cons = [
{'type': 'eq', 'fun': lambda z: constraint_eq(z)},
{'type': 'ineq', 'fun': lambda z: constraint_ineq(z)}
]
bounds = [(None, None)] * n + [(0, None)] * m
res = minimize(objective, z0, method='SLSQP', constraints=cons, bounds=bounds)
if res.success:
return res.x[:n]
else:
raise RuntimeError("Failed to find feasible point")
SQP 本质是局部优化器，单一起点很容易被困住。可以使用拉丁超立方采样生成几十个分散的初始点，分别运行 SQP，最后选出全局表现最佳的那个。
很多人觉得'反正只是近似，干嘛不用梯度法？'但忽略二阶信息的代价是惨重的。有了 Hessian 信息，就知道应该'提前转弯'，沿着谷底一路滑下去。
数学上，这就是二阶泰勒展开的魅力： $$ f(x_k + d) \approx f(x_k) + \nabla f(x_k)^T d + \frac{1}{2} d^T \nabla^2 f(x_k) d $$
去掉常数项，QP 子问题的目标函数就成了： $$ \min_d \quad \nabla f(x_k)^T d + \frac{1}{2} d^T B_k d $$
这里 $B_k$ 可以是精确 Hessian，也可以是拟牛顿近似（如 BFGS）。只要它正定，就能保证 QP 有唯一解。
import numpy as np
def quadratic_objective_approximation(xk, dk, grad_f, hess_approx):
linear_term = grad_f.T @ dk
quadratic_term = 0.5 * dk.T @ hess_approx @ dk
return linear_term + quadratic_term
现代做法是用自动微分（AD），比如 JAX、CasADi 或 PyTorch。它们能在一次反向传播中精确计算所有偏导。
import jax
import jax.numpy as jnp
def objective(x):
return jnp.sum(x**2) + jnp.sin(x[0]) * jnp.cos(x[1])
grad_objective = jax.grad(objective)
x = jnp.array([1.0, 2.0])
print("f(x) =", objective(x))
print("∇f(x) =", grad_objective(x))
即使用了 AD，也不能盲目相信线性化模型。我们可以定义一致性指标： $$ \gamma_k = \frac{|c(x_k + p_k)|}{|\tilde{c}(x_k + p_k)|} $$
如果误差持续增长超过 5%，就得警惕了——也许该缩小信任域半径，或者触发重线性化机制。
| 方法 | 收敛特性 | 计算成本 | 适用场景 |
|---|---|---|---|
| 精确 Hessian | 快速二阶收敛 | 高（需解析或 AD） | 小中规模、高精度需求 ✅ |
| BFGS 近似 | 超线性收敛 | 中等（存储完整矩阵） | 一般非线性问题 🔄 |
| L-BFGS 近似 | 接近超线性 | 低（仅存向量历史） | 大规模问题 📈 |
| 单位阵初始化 | 线性收敛 | 最低 | 初期粗略搜索 🔧 |
大多数情况下，推荐用 BFGS 或 L-BFGS。
def bfgs_update(B, s, y):
sy = s @ y
if sy <= 1e-8:
raise ValueError("Curvature condition violated: y^T s <= 0")
Bs = B @ s
outer1 = np.outer(Bs, Bs) / (Bs @ s)
outer2 = np.outer(y, y) / sy
return B - outer1 + outer2
当内存撑不住完整的 $n \times n$ 矩阵时，L-BFGS 登场了。它只保留最近 $m=5\sim20$ 次的 $(s_k, y_k)$ 对。
def lbfgs_direction(grad, s_list, y_list, rho_list, gamma0=1.0):
q = grad.copy()
alpha = np.zeros(len(s_list))
# 向后传递
for i in reversed(range(len(s_list))):
alpha[i] = rho_list[i] * (s_list[i] @ q)
q -= alpha[i] * y_list[i]
# 缩放
z = gamma0 * q
# 向前传递
for i in range(len(s_list)):
beta = rho_list[i] * (y_list[i] @ z)
z += s_list[i] * (alpha[i] - beta)
return -z
在机器人轨迹跟踪中，每毫秒都要解一个带动力学约束的 NLP 问题。SQP 凭借其快速局部收敛能力，能在几毫秒内完成求解。
% MATLAB 伪代码：MPC 中的 SQP 调用
for t = 1:T_horizon
obj = sum((x(t) - x_ref).^2 + u(t)^2);
subject_to:
x(t+1) == A*x(t) + B*u(t); % 动力学约束
u_min <= u(t) <= u_max;
end
sol = sqp_solver(obj, constraints, x0);
apply_control(sol.u(1));
飞机机翼减重设计中，既要满足应力限制，又要防止局部屈曲。这些约束高度非线性，传统方法极易发散。SQP 结合伴随法计算梯度，稳定性和效率双双在线。
在电池等效电路建模中，需联合估计十几个电化学参数。通过 SQP 优化，RMSE 可降低 60% 以上，极大提升 SOC 估算精度。
| 参数 | 初始猜测 | 估计结果 |
|---|---|---|
| R0 | 0.05 Ω | 0.048 Ω |
| R1 | 0.03 Ω | 0.031 Ω |
| C1 | 1000 F | 980 F |
| … | … | … |
别只盯着目标函数下降！要综合判断：
def is_converged(x_prev, x_curr, f_prev, f_curr, kkt_res, max_iter, iter_count):
dx = np.linalg.norm(x_curr - x_prev)
df = abs(f_curr - f_prev)
return (dx < 1e-5) and (df < 1e-7) and (kkt_res < 1e-6) and (iter_count < max_iter)
随着 GPU 加速和 ADMM 框架的发展，分布式 SQP 正在兴起。例如在电网调度中，各节点本地运行子问题，仅交换少量协调变量，既保护隐私又提升效率。
回顾全文，你会发现 SQP 的强大之处不在于某个单一技术，而在于它把多个精巧的思想拧成了一股绳：
所以，无论你是做自动驾驶、智能制造，还是量化交易，只要遇到带约束的连续优化问题，不妨试试 SQP。
微信公众号「极客日志」，在微信中扫描左侧二维码关注。展示文案：极客日志 zeeklog
使用加密算法（如AES、TripleDES、Rabbit或RC4）加密和解密文本明文。 在线工具，加密/解密文本在线工具，online
解析常见 curl 参数并生成 fetch、axios、PHP curl 或 Python requests 示例代码。 在线工具，curl 转代码在线工具，online
将字符串编码和解码为其 Base64 格式表示形式即可。 在线工具，Base64 字符串编码/解码在线工具，online
将字符串、文件或图像转换为其 Base64 表示形式。 在线工具，Base64 文件转换器在线工具，online
将 Markdown（GFM）转为 HTML 片段，浏览器内 marked 解析；与 HTML转Markdown 互为补充。 在线工具，Markdown转HTML在线工具，online
将 HTML 片段转为 GitHub Flavored Markdown，支持标题、列表、链接、代码块与表格等；浏览器内处理，可链接预填。 在线工具，HTML转Markdown在线工具，online