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

这段代码虽然短，但它干的是'验明正身'的活儿：检查当前点是否接近 KKT 点。如果 lagrangian_gradient 接近零向量，说明系统达到了力学上的平衡状态，离最优解不远了！

不过要注意，KKT 只是必要条件，不是充分条件。除非问题是凸的，否则你找到的可能是鞍点或者局部极小。这就像是找到了一个山谷底部，但不能确定是不是全球最低点。

SQP 如何一步步逼近 KKT 点？

既然我们没法一次性搞定全局最优，那就换个思路：每一步都尝试让当前估计点更靠近 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)$，会发生什么？

• 在该点处，约束梯度是 $(2,0)$，线性化后变成 $2(x\_1 - 1) \leq 0 \Rightarrow x\_1 \leq 1$
• 这个近似忽略了圆的曲率，允许沿 $x\_2$ 方向自由移动
• 但由于真实约束是圆形，稍大一点的步长就会导致不可行

结果就是：算法在边界附近反复试探、收缩步长、甚至发散。明明离最优解很近，却因为'地图画歪了'而寸步难行。

所以啊，初始点不仅是起点，更是决定算法命运的第一推手。

如何获得高质量的初始解？

方法一：从松弛问题出发

如果你手头的初始猜测严重违反约束，别硬上！先做一个'可行性恢复阶段'：

$$ \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 $$

这个辅助问题的目标是尽可能减少总的约束违反量。一旦找到一个勉强可行的点，再把它丢给主 SQP 流程。

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)]

    # 构造变量：[x; s]
    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] # 返回修复后的 x
    else:
        raise RuntimeError("Failed to find feasible point")

这段代码简直就是'急救包'——哪怕你的初始点完全不可行，也能把它拉回正轨。

方法二：多起点策略，打破局部最优魔咒

SQP 本质是局部优化器，单一起点很容易被困住。怎么办？简单粗暴但有效：多跑几次，挑最好的那次结果。

你可以用拉丁超立方采样生成几十个分散的初始点，分别运行 SQP，最后选出全局表现最佳的那个。

from multiprocessing import Pool
import numpy as np

def multi_start_sqp(objective, grad_obj, constraints, bounds, n_starts=20):
    results = []
    dim = len(bounds)

    def worker(seed):
        np.random.seed(seed)
        x0 = np.array([np.random.uniform(low, high) for low, high in bounds])
        try:
            # 先做可行性修复
            x0_feas = feasibility_restoration(
                [lambda x: con['fun'](x) for con in constraints if con['type']=='ineq'],
                [lambda x: con['fun'](x) for con in constraints if con['type']=='eq'],
                x0
            )
            # 运行主 SQP
            res = minimize(objective, x0_feas, jac=grad_obj, method='SLSQP', bounds=bounds, constraints=constraints)
            return res if res.success else None
        except:
            return None

    with Pool() as pool:
        results = pool.map(worker, range(n_starts))
    valid_results = [r for r in results if r is not None]
    if valid_results:
        best = min(valid_results, key=lambda r: r.fun)
        return best
    else:
        raise RuntimeError("No successful run among all starts.")

实验表明，即使在 1030 维空间，使用 2050 次重启就能大幅提升找到全局优解的概率。这招在实际工程中特别管用，尤其是在缺乏先验知识的时候。

QP 子问题建模：如何打造可靠的'局部地图'

目标函数的二次逼近：为什么要保留曲率？

很多人觉得'反正只是近似，干嘛不用梯度法？'但忽略二阶信息的代价是惨重的。

想象你要穿过一条弯曲的山谷。如果只看梯度方向，你会不断左右摇摆（zig-zag），效率极低。而有了 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

别小看这一行矩阵乘法，它背后可是决定了你能否顺利穿越