二次序列规划（SQP）算法详解与实战

这段代码检查当前点是否接近 KKT 点。如果 lagrangian_gradient 接近零向量，说明系统达到了力学上的平衡状态。

SQP 如何一步步逼近 KKT 点？

具体做法是在当前点 $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} $$

这里面藏着三个关键思想：

1. 二阶泰勒展开：保留目标函数的曲率信息，避免纯梯度法的锯齿震荡；
2. 一阶线性化：把非线性约束变成线段，便于 QP 求解；
3. 方向更新机制：解出搜索方向 $d_k$ 后，通过步长控制推进新点 $x_{k+1} = x_k + \alpha_k d_k$。

整个流程可以用一张图概括：

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，最后选出全局表现最佳的那个。

QP 子问题建模

目标函数的二次逼近

很多人觉得'反正只是近似，干嘛不用梯度法？'但忽略二阶信息的代价是惨重的。有了 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

Jacobian 计算：自动微分才是王道

现代做法是用自动微分（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 近似策略

精确 Hessian vs 拟牛顿法

大多数情况下，推荐用 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

L-BFGS：大规模问题的秘密武器

当内存撑不住完整的 $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

工程实战：SQP 在真实世界的应用场景

场景一：模型预测控制（MPC）

在机器人轨迹跟踪中，每毫秒都要解一个带动力学约束的 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 估算精度。

收敛性保障与性能调优

多维度收敛判据

别只盯着目标函数下降！要综合判断：

• 梯度范数 $ |\nabla \mathcal{L}| < 10^{-6} $
• 变量增量 $ |Δx| < 10^{-5} $
• 目标变化 $ |Δf| < 10^{-7} $
• 约束违反度 $ |c(x)| < 10^{-6} $

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)

数值稳定性增强技巧

• 正则化 Hessian：加个小扰动 $\mu I$，防止奇异；
• 信任域限制：不让步子迈太大，避免模型失效；
• 混合求解器切换：QP 失败时自动降级到活动集法或投影梯度法。

未来展望：并行化与分布式 SQP

随着 GPU 加速和 ADMM 框架的发展，分布式 SQP 正在兴起。例如在电网调度中，各节点本地运行子问题，仅交换少量协调变量，既保护隐私又提升效率。

结语

回顾全文，你会发现 SQP 的强大之处不在于某个单一技术，而在于它把多个精巧的思想拧成了一股绳：

• 用 KKT 条件锚定理论基础，
• 用 QP 子问题实现局部精确建模，
• 用拟牛顿法平衡效率与精度，
• 用模块化设计支撑工程落地。

所以，无论你是做自动驾驶、智能制造，还是量化交易，只要遇到带约束的连续优化问题，不妨试试 SQP。