OSQP: Operator Splitting Quadratic Program

📖 1. 概述

OSQP（Operator Splitting Quadratic Program）是一种一阶算子分裂求解器，专门用于求解凸二次规划问题。它基于 ADMM（Alternating Direction Method of Multipliers，交替方向乘子法）框架，将复杂的 QP 问题分解为一系列简单的子问题迭代求解。

📏 2. 二次规划问题

2.1 标准形式

OSQP 求解的 QP 问题标准形式为：

其中各变量含义如下：

2.2 直观示例

考虑一个简单的二维 QP 问题：

对应 OSQP 标准形式：

🔄 3. ADMM 算法

3.1 基本原理

ADMM（Alternating Direction Method of Multipliers）是解决可分离凸优化问题的一阶方法，结合了对偶分解和增广拉格朗日方法的优点。

考虑如下可分离问题：

ADMM 的增广拉格朗日函数为：

其中 \(\rho > 0\) 为惩罚参数，\(y\) 为对偶变量（拉格朗日乘子）。

3.2 分裂策略

ADMM 将原问题分裂为三个交替更新步骤：

3.3 收敛性分析

🔧 4. QP 到 ADMM 重构

OSQP 的核心创新在于将标准 QP 问题重构为 ADMM 可高效求解的形式。

4.1 第一步：变量分裂

引入辅助变量 \(z\)，将原问题等价转化为：

其中 \(\delta_{[l,u]}(z)\) 是区间 \([l,u]\) 的指示函数：

4.2 第二步：增广拉格朗日

对约束 \(Ax = z\) 引入对偶变量 \(y\)，构造增广拉格朗日函数：

增广项 \(\frac{\rho}{2}\|Ax - z\|_2^2\) 的作用：

• 保证子问题严格凸，提高数值稳定性
• 使收敛不依赖于原始函数的凸性
• \(\rho\) 越大，对约束违反的惩罚越重

4.3 第三步：算子分裂

将增广拉格朗日函数按变量 \(x\) 和 \(z\) 分裂：

🧩 5. 求解流程

基于上述重构，OSQP 的 ADMM 迭代步骤如下：

5.1 x 更新（线性方程组求解）

这是 x 更新的核心公式。令 \(\mathcal{H} = P + \rho A^T A\)（常称为组合 Hessian），则：

5.2 z 更新（区间投影）

投影操作为逐元素裁剪：

5.3 对偶变量更新

这是简单的梯度上升步，对偶变量 \(\rho^{-1}y\) 收敛到约束 \(Ax = z\) 的拉格朗日乘子。

5.4 终止条件

OSQP 同时监控原残差和对偶残差：

5.5 完整算法流程

💻 6. 代码实现

6.1 Python 接口示例

import osqp
import numpy as np
from scipy import sparse

# 定义 QP 问题数据
P = sparse.csc_matrix(np.array([[2.0, 0.0], [0.0, 1.0]]))
q = np.array([1.0, 1.0])

A = sparse.csc_matrix(np.array([[1.0, 0.0], [0.0, 1.0]]))
l = np.array([0.0, 0.0])
u = np.array([2.0, 3.0])

# 创建求解器并求解
solver = osqp.OSQP()
solver.setup(P, q, A, l, u, alpha=1.6, eps_abs=1e-6, eps_rel=1e-6)
result = solver.solve()

print(f"最优解: x = {result.x}")
print(f"最优值: p = {result.info.obj_val}")
print(f"迭代次数: {result.info.iterations}")

6.2 手动实现 OSQP 核心

import numpy as np
from scipy import sparse
from scipy.sparse.linalg import spsolve

def osqp_solve(P, q, A, l, u, rho=1.0, max_iter=4000,
               eps_abs=1e-4, eps_rel=1e-4, alpha=1.6):
    n = P.shape[0]
    m = A.shape[0]

    # 预处理: 计算 H = P + rho * A^T A
    H = P + rho * A.T @ A
    H = H.tocsc()

    # 初始化
    x = np.zeros(n)
    z = np.zeros(m)
    y = np.zeros(m)

    for k in range(max_iter):
        # x 更新: H x = -q + A^T(rho*z - y)
        b = -q + A.T @ (rho * z - y)
        x = spsolve(H, b)

        # z 更新: 区间投影
        Ax = A @ x
        z_old = z.copy()
        z = np.clip(Ax + y / rho, l, u)

        # 对偶变量更新
        y = y + rho * (Ax - z)

        # 计算残差
        primal_res = np.linalg.norm(Ax - z)
        dual_res = np.linalg.norm(rho * A.T @ (z - z_old))

        # 检查终止条件
        eps_pri = eps_abs + eps_rel * max(np.linalg.norm(Ax), np.linalg.norm(z))
        eps_dual = eps_abs + eps_rel * np.linalg.norm(rho * A.T @ z)
        if primal_res < eps_pri and dual_res < eps_dual:
            print(f"收敛于 {k+1} 次迭代")
            break

    return x, y

6.3 C 语言核心循环

void osqp_solve(Workspace *work) {
    int n = work->n;
    int m = work->m;

    for (int iter = 0; iter < work->settings->max_iter; iter++) {
        // x 更新
        vec_add_scaled(work->x, work->z, work->y, 1.0, -1.0/rho, m);
        mat_vec_mul(work->A, work->x, work->tmp, m, n);
        vec_add_scaled(work->rhs, work->q, work->y, -1.0, 1.0, n);
        // 求解 H * x_new = rhs (使用预计算的 Cholesky 分解)
        ldl_solve(work->chol, work->x, work->rhs);

        // z 更新 (区间投影)
        mat_vec_mul(work->A, work->x, work->Ax, m, n);
        for (int i = 0; i < m; i++) {
            work->z[i] = proj(work->Ax[i] + work->y[i]/rho,
                              work->l[i], work->u[i]);
        }

        // 对偶变量更新
        vec_add_scaled(work->y, work->Ax, work->z, 1.0, -1.0, m);

        // 残差检查
        if (check_convergence(work)) break;
    }
}

float proj(float val, float lo, float hi) {
    return fmax(lo, fmin(hi, val));
}

🚀 7. 应用场景

# MPC 问题: 每步求解 QP
# min sum_{k=0}^{N-1} (x_k' Q x_k + u_k' R u_k)
# s.t. x_{k+1} = A x_k + B u_k
#      x_min <= x_k <= x_max
#      u_min <= u_k <= u_max

# 将 MPC 问题转化为标准 QP 形式
# P: 分块对角矩阵 (Q, R)
# A: 状态转移约束矩阵
# l, u: 变量上下界
solver = osqp.OSQP()
solver.setup(P, q, A, l, u, warm_start=True)
# 每个控制周期只需很少的迭代次数（热启动）
result = solver.solve()

📊 8. 对比

❓ 9. 常见问题

9.1 为什么 OSQP 不需要矩阵分解？

9.2 \(\rho\) 和 \(\alpha\) 如何选择？

9.3 如何处理非凸 QP？

9.4 热启动的原理？

热启动（Warm Start）利用上一次 QP 的解作为当前 QP 的初始点：

在 MPC 等连续时间步问题中，相邻 QP 的解通常接近，热启动可将迭代次数从数百次减少到 10-50 次。

9.5 OSQP 与 cvxpy 的关系？

CVXPY 是一个建模语言，自动将问题转化为标准形式并选择后端求解器：

📚 10. 参考文献

1. Stellato, B., et al. (2020). OSQP: An Operator Splitting Solver for Quadratic Programs. Mathematical Programming Computation, 12, 637-672.
2. Boyd, S., et al. (2011). Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers. Foundations and Trends in Machine Learning, 3(1), 1-122.
3. Goldfarb, D., et al. (2017). Operator Splitting Methods for Large Scale Convex Optimization. IMA Journal of Numerical Analysis.
4. Parikh, N., & Boyd, S. (2014). Proximal Algorithms. Foundations and Trends in Optimization, 1(3), 127-239.
5. osqp.org: OSQP 官方文档与教程.