1. 从“调包”到“手撕”：为什么QP值得你花时间
很多做算法落地的朋友，第一次接触二次规划（Quadratic Programming，简称QP）都是在MATLAB的quadprog或者Python的cvxopt.solvers.qp里填几个矩阵就完事了。模型跑通那一刻觉得“不过如此”，直到某天遇到求解器返回infeasible或者结果飘得离谱，才意识到自己连问题到底是不是凸的、约束有没有写对、求解器内部在干什么都一无所知。我见过太多项目卡在QP这一环：模型预测控制（MPC）里每毫秒要解一次QP，投资组合优化里协方差矩阵稍微病态一点就崩，SVM对偶问题里核矩阵规模一上来内存直接爆掉。这些场景逼着你不能只当调包侠，必须理解QP的结构、求解器的脾气，以及怎么把原始问题改写成数值上更友好的形式。
这篇文章面向的是已经会用求解器、但想搞清楚“为什么这么写”“为什么解不出来”“怎么让它解得更快更稳”的从业者。我会从QP的标准形式讲起，拆解凸性判断、KKT条件、有效集法、内点法的核心逻辑，然后落到实操：怎么选求解器、怎么诊断不可行、怎么处理病态问题、怎么利用稀疏性加速。全文没有花哨的数学炫技，只有我踩过的坑和验证过的方案。如果你正在做MPC、组合优化、SVM、图像去噪或者任何需要最小化二次目标且带约束的任务，这篇内容应该能帮你省下不少调试时间。
2. QP的标准形式与凸性判断：别在非凸问题上浪费时间
2.1 标准形式到底长什么样
二次规划的标准形式通常写成：
minimize (1/2) * x^T * P * x + q^T * x subject to A * x <= b A_eq * x = b_eq lb <= x <= ub其中P是Hessian矩阵，q是线性项系数，A和A_eq分别是不等式和等式约束矩阵。目标函数是二次的，约束是线性的，这就是QP的边界。如果约束里出现x^2或者x*y这种项，那就不是QP了，得归到非线性规划（NLP）去。
我刚开始做MPC的时候，总喜欢把目标函数写成(x - x_ref)^T * Q * (x - x_ref) + u^T * R * u这种形式，然后直接丢给求解器。后来发现，展开成标准形式后P矩阵的结构直接决定了求解器能不能用Cholesky分解，进而决定求解速度。展开过程其实很简单：
(x - x_ref)^T * Q * (x - x_ref) = x^T * Q * x - 2 * x_ref^T * Q * x + const常数项对优化变量没有影响，可以直接扔掉。所以P = 2Q（如果标准形式里带1/2的话），q = -2 * Q * x_ref。这个展开看起来是初中代数，但很多人在这里犯错：忘记乘2，或者把Q和2Q搞混，导致求解结果和预期差一个系数。我的建议是，不管多简单的展开，都写个小脚本用有限差分验证一下梯度，确认P和q没写错。
2.2 凸性判断：P矩阵的半正定性是生死线
QP分凸和非凸，凸QP有全局最优解，非凸QP可能陷入局部最优甚至无界。判断标准就一条：Hessian矩阵P必须是半正定的（所有特征值非负）。如果P是正定的（所有特征值严格大于0），那目标函数是严格凸的，最优解唯一。
实际项目中，P矩阵往往来自物理模型的能量函数或者统计模型的协方差矩阵，理论上应该半正定，但数值上经常出问题。我遇到过好几次：协方差矩阵因为数据点太少或者特征高度相关，出现很小的负特征值（比如 -1e-12），求解器直接报“矩阵不是半正定”。这时候别急着换求解器，先检查P的构造过程。常见的修复手段有：
- 特征值截断：对
P做特征分解，把所有小于tol的特征值置为tol（比如1e-8），然后重组矩阵。这相当于给问题加了一点正则化，牺牲一点精度换数值稳定性。 - 对角加载：直接
P = P + epsilon * I，epsilon取1e-6到1e-8。这个方法简单粗暴，但会改变问题的最优解，需要评估对结果的影响。 - 对称化：有时候
P因为浮点误差不对称，求解器会认为它不正定。强制P = (P + P^T) / 2能解决一部分问题。
注意：对角加载虽然好用，但在MPC里会改变闭环系统的稳定性。如果
P来自LQR的代价矩阵，加载后相当于改变了权重，可能导致控制性能下降。我的经验是，优先用特征值截断，实在不行再加载，并且加载后一定要重新仿真验证。
2.3 等式约束和不等式约束的处理差异
等式约束A_eq * x = b_eq在QP里是“硬约束”，必须严格满足。不等式约束A * x <= b是“软约束”，允许在边界上。这个区别直接影响求解器的设计：等式约束可以用消元法直接降维，不等式约束需要迭代确定哪些是活跃的。
我见过一个典型的坑：在组合优化里，权重和等于1是等式约束，单个权重非负是不等式约束。有人把权重和写成sum(w) <= 1.0001和sum(w) >= 0.9999两个不等式，结果求解器在边界附近来回震荡，收敛极慢。正确做法是保留一个等式约束sum(w) = 1，非负用下界lb = 0处理。等式约束越少，问题维度越低，求解越快。
3. KKT条件与对偶问题：理解求解器在算什么
3.1 KKT条件：最优解的“身份证”
对于凸QP，KKT条件是最优解的充要条件。标准形式的KKT条件包括：
- 平稳性：
P * x + q + A^T * lambda + A_eq^T * nu = 0 - 原始可行性：
A * x <= b，A_eq * x = b_eq - 对偶可行性：
lambda >= 0 - 互补松弛：
lambda_i * (b_i - A_i * x) = 0对所有不等式约束成立
互补松弛是理解有效集法的关键：如果某个不等式约束没有取等号（即A_i * x < b_i），那么对应的拉格朗日乘子lambda_i必须为0；反过来，如果lambda_i > 0，那个约束必须取等号。这意味着最优解处，只有一部分不等式约束是“活跃”的，其余都是“非活跃”的。
我在调试MPC时经常用这个性质来验证结果：求解器返回的lambda向量里，非零元素对应的约束应该正好是那些取等号的。如果发现某个约束明显没取等号但lambda很大，那要么是求解器没收敛，要么是问题本身有数值问题。
3.2 对偶问题：换个角度看QP
QP的对偶问题也是QP，而且有时候对偶形式比原始形式更好解。原始问题：
min (1/2) x^T P x + q^T x s.t. A x <= b拉格朗日函数是L(x, lambda) = (1/2) x^T P x + q^T x + lambda^T (A x - b)。对x求导令其为零，得到x = -P^{-1} (q + A^T lambda)。代回拉格朗日函数，得到对偶问题：
max - (1/2) (q + A^T lambda)^T P^{-1} (q + A^T lambda) - b^T lambda s.t. lambda >= 0对偶问题的变量维度等于不等式约束的个数，而不是原始变量的维度。在SVM里，原始变量是权重向量（维度可能上万），对偶变量是样本对应的乘子（维度等于样本数）。当样本数远小于特征维度时，解对偶问题更划算。这就是SVM用对偶形式的根本原因。
但要注意：对偶问题需要计算P^{-1}，如果P病态或者奇异，对偶形式反而更不稳定。我在做投资组合优化时，协方差矩阵经常接近奇异，这时候直接解原始问题（用内点法）比解对偶问题更稳。选择原始还是对偶，取决于P的条件数和约束数量，没有绝对答案。
3.3 强对偶性与Slater条件
凸QP满足Slater条件时，强对偶性成立，原始问题和对偶问题的最优值相等。Slater条件要求存在一个严格可行点，即满足A x < b的点。如果所有可行点都落在边界上（比如等式约束和不等式约束把可行域压成一条线），Slater条件可能不满足，强对偶性就不保证。
实际中，如果发现原始问题和对偶问题的最优值对不上，先检查Slater条件。我遇到过一次：可行域是一个点（所有约束都是等式），原始问题有唯一解，但对偶问题无解。这种情况下，求解器可能返回奇怪的结果，需要手动处理退化情况。
4. 有效集法 vs 内点法：求解器的两条技术路线
4.1 有效集法：像剥洋葱一样逐层确定活跃约束
有效集法的思路很直观：先猜一组活跃约束（通常从等式约束开始），把不等式约束当成等式来解，然后检查解是否满足所有不等式。如果不满足，就把违反最严重的约束加入活跃集；如果某个活跃约束的乘子为负，就把它移出活跃集。重复这个过程，直到活跃集不再变化。
这个方法的优点是：
- 对于小规模问题（变量少于100个），收敛很快，通常几次迭代就搞定。
- 热启动非常方便：如果连续解一系列相似的QP（比如MPC的每个时间步），可以把上一步的活跃集作为初始猜测，迭代次数能减少一个数量级。
- 返回的活跃集直接告诉你哪些约束在起作用，对分析问题很有帮助。
缺点是：
- 最坏情况下迭代次数可能是指数级的（虽然实际中很少遇到）。
- 对于大规模稀疏问题，活跃集的更新可能破坏稀疏性，导致计算量爆炸。
- 处理退化问题（多个约束同时活跃）时可能循环。
我在MPC里用的qpoases就是有效集法的代表。它的热启动功能让MPC的求解时间从毫秒级降到微秒级。但有一次遇到约束退化，求解器在两组活跃集之间来回跳，最后超时。解决办法是给约束加一点松弛（A x <= b + 1e-8），打破退化。
4.2 内点法：沿着中心路径逼近最优解
内点法（也叫障碍法）的思路完全不同：它把不等式约束通过障碍函数（比如对数障碍）加到目标函数里，把约束问题转化成无约束问题，然后用牛顿法求解。随着障碍参数趋近于0，无约束问题的解趋近于原始问题的最优解。
内点法的优点是：
- 多项式时间复杂度，最坏情况下迭代次数是
O(sqrt(n) * log(1/epsilon))，其中n是变量数，epsilon是精度。 - 对大规模稀疏问题非常友好，因为牛顿法每步只需要解一个稀疏线性系统。
- 不需要猜活跃集，对退化问题不敏感。
缺点是：
- 热启动不如有效集法方便。虽然有一些热启动技术，但效果通常不如有效集法。
- 每次迭代的计算量比有效集法大，对于小规模问题反而更慢。
- 障碍参数趋近于0时，线性系统可能变得病态，需要仔细处理数值精度。
cvxopt的默认QP求解器就是内点法。我在做图像去噪时，问题规模是1000x1000的像素，变量数上百万，但约束是稀疏的。内点法每步解一个稀疏系统，几分钟就能收敛。换成有效集法，活跃集更新会破坏稀疏性，内存直接爆掉。
4.3 怎么选：看问题规模和实时性要求
| 特性 | 有效集法 | 内点法 |
|---|---|---|
| 小规模问题（n<100） | 快，迭代少 | 相对慢 |
| 大规模稀疏问题 | 稀疏性易破坏 | 稀疏性保持好 |
| 热启动 | 非常方便 | 不方便 |
| 退化问题 | 可能循环 | 稳定 |
| 返回活跃集 | 直接返回 | 需要额外计算 |
| 实时控制（MPC） | 首选 | 次选 |
| 组合优化（金融） | 中等规模可用 | 大规模首选 |
我的经验是：MPC、机器人控制这类需要毫秒级响应的场景，优先用有效集法加暖启动；图像处理、统计学习这类离线或准在线场景，用内点法更省心。如果问题规模在100到1000之间，两种方法都可以试，看哪个实际跑得快。
5. 实操中的数值陷阱：从不可行到病态
5.1 不可行诊断：到底是模型错了还是求解器没找到
求解器返回infeasible时，新手第一反应是“模型写错了”，但有时候模型没问题，只是数值上不可行。比如约束x >= 1和x <= 1 - 1e-16，理论上可行域是一个点，但浮点表示下可能判定为空。
诊断不可行的标准方法是解一个松弛问题：给每个约束加一个非负的松弛变量，最小化松弛变量的和。如果最小和大于0，说明原始问题确实不可行；如果最小和等于0但求解器仍报不可行，那就是数值问题。
# 用cvxopt诊断不可行 from cvxopt import matrix, solvers # 原始约束 A*x <= b # 松弛问题: min sum(s) s.t. A*x <= b + s, s >= 0 # 如果最优sum(s) > tol，则不可行我遇到过一次：MPC的约束是x_min <= x <= x_max，但x_min和x_max因为传感器噪声出现了x_min > x_max的情况。求解器报不可行，我查了半天模型，最后发现是传感器数据的问题。所以，先检查约束边界本身是否合理，再怀疑求解器。
5.2 病态问题：条件数太大怎么办
P矩阵的条件数（最大特征值除以最小特征值）如果超过1e10，求解器基本就废了。条件数大的原因是变量之间的尺度差异太大。比如一个变量是电压（量级1e3），另一个是电流（量级1e-3），Hessian矩阵的对角元素相差1e6，条件数直接爆炸。
解决办法是变量缩放：把每个变量除以其典型量级，让所有变量的量级接近1。具体操作是构造一个对角矩阵D，令x = D * y，然后对y求解。缩放后的Hessian是D^T * P * D，条件数通常能降几个数量级。
我在做电路优化时，电压和电流的缩放比例是1000:1，缩放后求解时间从30秒降到0.5秒。缩放矩阵的构造没有固定公式，我的经验是取每个变量在可行域内的典型值或者上下界的中值。
提示：缩放后记得把解变换回原始变量
x = D * y，并且检查约束是否也需要相应缩放。等式约束A_eq * x = b_eq变成A_eq * D * y = b_eq，不等式同理。
5.3 稀疏性利用：别让求解器做无用功
大规模QP的P和A通常是稀疏的，但如果你用密集矩阵存储，内存和计算量都会爆炸。cvxopt支持稀疏矩阵，但需要手动构造。scipy.sparse构造的矩阵可以直接传给cvxopt，但要注意索引格式（CSC vs CSR）和数据类型。
我做过一个实验：10000个变量的QP，P是稀疏的（每行只有10个非零元素）。用密集矩阵存储需要800MB内存，求解时间超过10分钟；用稀疏矩阵存储只需要8MB，求解时间20秒。差距就是这么夸张。
构造稀疏矩阵时，尽量用coo_matrix或lil_matrix逐行构造，然后转成csc_matrix传给求解器。csc_matrix对列操作友好，内点法的Cholesky分解通常按列进行，所以CSC格式最快。
6. 手写一个QP求解器：从理论到代码的最后一公里
6.1 为什么建议你手写一遍
调包能解决80%的问题，但剩下20%的疑难杂症——比如求解器不收敛、结果不符合预期、需要自定义约束——逼着你理解内部机制。手写一个简单的有效集法求解器，不需要多复杂，100行代码就能搞定等式约束QP和简单的不等式约束QP。写完之后，你对KKT条件、活跃集、数值精度的理解会上一个台阶。
6.2 等式约束QP的解析解
等式约束QP有闭式解。问题：
min (1/2) x^T P x + q^T x s.t. A_eq x = b_eqKKT系统是：
[ P A_eq^T ] [ x ] = [ -q ] [ A_eq 0 ] [ nu] [ b_eq ]这是一个线性系统，直接解就行。但要注意：如果A_eq行满秩且P正定，这个系统有唯一解。如果A_eq行不满秩（约束冗余），系统奇异，需要处理。我的做法是用最小二乘或者加一点正则化。
import numpy as np def solve_eq_qp(P, q, A_eq, b_eq): n = P.shape[0] m = A_eq.shape[0] KKT = np.block([[P, A_eq.T], [A_eq, np.zeros((m, m))]]) rhs = np.concatenate([-q, b_eq]) sol = np.linalg.solve(KKT, rhs) return sol[:n], sol[n:]这个函数能解任何等式约束凸QP。我用手写版本验证过cvxopt的结果，在P条件数不大的情况下，两者差异在1e-10以内。
6.3 有效集法的完整实现框架
有效集法的核心循环：
- 初始化工作集
W（通常包含所有等式约束）。 - 解等式约束QP（只考虑
W中的约束）。 - 检查解是否满足所有不等式约束。如果满足，检查
W中不等式约束的乘子是否非负。如果都非负，停止；否则移除乘子为负的约束。 - 如果有不等式约束被违反，把违反最严重的加入
W。 - 重复直到收敛。
这个框架看起来简单，但实现时有几个坑：
- 步长控制：从当前点移动到等式约束QP的解时，不能越过其他约束的边界。需要计算最大步长。
- 退化处理：如果步长为0（当前点已经在某个约束边界上），需要特殊处理，否则死循环。
- 数值容差：判断约束是否活跃、乘子是否为负，都需要容差。容差太小会循环，太大会漏掉约束。
我写的第一版有效集法在退化问题上卡了整整一天，后来加了步长下限1e-12和乘子容差1e-8才稳定。这些经验在教科书里不会写，只有自己踩过才知道。
7. 求解器选型实战：不同场景下的最优选择
7.1 开源求解器对比
| 求解器 | 方法 | 语言 | 热启动 | 稀疏支持 | 适用场景 |
|---|---|---|---|---|---|
| qpoases | 有效集 | C++ | 优秀 | 中等 | MPC、嵌入式 |
| cvxopt | 内点 | Python | 差 | 优秀 | 统计学习、金融 |
| OSQP | 算子分裂 | C | 优秀 | 优秀 | 大规模稀疏、MPC |
| quadprog | 有效集 | MATLAB | 中等 | 中等 | 教学、小规模 |
| Gurobi | 内点+有效集 | C | 优秀 | 优秀 | 商业级、大规模 |
OSQP是近几年比较火的求解器，基于ADMM（交替方向乘子法），对大规模稀疏问题非常友好，而且支持热启动。我在一个MPC项目里用OSQP替换了qpoases，求解时间从5ms降到1ms，代码还更简洁。但OSQP的精度不如内点法，对于需要高精度解的场景（比如金融组合优化），还是得用内点法。
7.2 商业求解器的优势与代价
Gurobi和Mosek是商业求解器的代表，它们的特点是：
- 自动选择算法（内点法、有效集法、混合）。
- 数值稳定性极好，几乎不会因为病态问题崩溃。
- 支持各种约束类型（整数、锥、半定）。
- 提供详细的诊断信息。
代价是贵。一个单机license一年几万到几十万不等。如果项目预算有限，开源求解器完全够用，只是需要更多调试。我的建议是：先用开源求解器跑通，遇到无法解决的数值问题再考虑商业求解器。很多时候，问题出在模型本身而不是求解器。
7.3 嵌入式场景的特殊考量
嵌入式平台（ARM、DSP、FPGA）上跑QP，内存和算力都受限。这时候：
- 优先选有效集法，因为迭代次数少，内存占用小。
- 避免动态内存分配，所有数组预分配。
- 用定点数代替浮点数（如果精度允许）。
- 利用问题的稀疏结构，手动展开矩阵运算。
我在一个ARM Cortex-M7上跑MPC，变量数20，约束数40，用qpoases的嵌入式版本，单次求解0.3ms，完全满足1kHz的控制频率。关键是把所有矩阵都放在静态内存里，避免malloc。
8. 从QP到MPC：一个完整的落地案例
8.1 MPC的QP形式化
模型预测控制（MPC）在每个时间步解一个QP：
min sum_{k=0}^{N-1} (x_k^T Q x_k + u_k^T R u_k) + x_N^T Qf x_N s.t. x_{k+1} = A x_k + B u_k x_min <= x_k <= x_max u_min <= u_k <= u_max把状态和输入堆叠成一个大向量z = [x_0, u_0, x_1, u_1, ..., x_N]，动力学方程变成等式约束，代价函数变成二次型，边界变成不等式约束。这就是一个标准QP。
8.2 热启动与暖启动的区别
热启动（warm start）是用上一次的解作为初始猜测，暖启动（hot start）是用上一次的活跃集作为初始工作集。有效集法对暖启动非常敏感，用对了能减少90%的迭代。OSQP也支持热启动，但效果不如有效集法明显。
我在MPC里用qpoases的暖启动，第一次求解用了15次迭代，之后每次只要2到3次。如果没有暖启动，每次都要15次，控制频率直接减半。
8.3 实时性保障的工程技巧
- 固定迭代次数：不要等求解器收敛，设定最大迭代次数（比如5次），返回当前最优解。MPC对次优解有容忍度，只要闭环稳定就行。
- 提前终止：如果当前解已经满足所有约束且乘子非负，立即返回。
- 降级策略：如果求解失败，用上一次的控制量或者零控制量，保证系统安全。
- 代码生成：用
CVXGEN或者FORCES生成定制化的求解器代码，比通用求解器快一个数量级。
这些技巧在教科书里不会讲，但实际项目中缺一不可。我见过一个无人机项目，因为MPC求解超时导致炸机，后来加了固定迭代次数和降级策略才稳定。
9. 我踩过的三个坑和对应的解决方案
9.1 坑一：P矩阵不对称导致求解器报错
有一次用cvxopt解QP，报错“matrix is not positive definite”。我检查了特征值，全是正的，纳闷了半天。后来发现P矩阵因为浮点误差不对称，P[0,1]和P[1,0]差了1e-16。cvxopt内部用Cholesky分解，要求矩阵严格对称。解决办法很简单：P = (P + P.T) / 2。这个坑我踩了两次，后来养成习惯，构造完P第一件事就是对称化。
9.2 坑二：约束缩放不一致导致不可行
在一个组合优化问题里，约束是sum(w) = 1和w >= 0。我把w缩放了100倍，但忘记缩放等式约束，结果sum(w) = 1变成sum(w_scaled) = 100，求解器报不可行。缩放必须一致：变量缩放后，所有涉及该变量的约束都要相应缩放。等式约束右边乘以缩放因子，不等式约束两边都乘。
9.3 坑三：内点法精度设置过高导致不收敛
cvxopt的默认精度是1e-7，我为了“更精确”改成了1e-12，结果求解器跑了10分钟没收敛。内点法的精度和迭代次数是对数关系，精度提高5个数量级，迭代次数增加约5倍，而且数值误差可能让求解器无法达到那么高的精度。我的经验是：精度设1e-6到1e-8就够了，再高没有意义，反而增加不收敛的风险。
10. 进阶方向：从QP到更一般的凸优化
QP是凸优化的一个特例。如果你已经掌握了QP，下一步可以往这几个方向扩展：
- 二阶锥规划（SOCP）：约束从线性变成二阶锥，能处理鲁棒优化、范数约束等问题。求解器和QP类似，只是每步多了一个锥投影。
- 半定规划（SDP）：变量从向量变成矩阵，约束是矩阵半正定。能处理控制中的Lyapunov稳定性、组合优化中的松弛问题。
- 整数二次规划（MIQP）：变量部分或全部是整数，需要分支定界。求解时间随整数变量数指数增长，适合小规模问题。
这些扩展的求解器（如Mosek、Gurobi）都支持，但底层还是QP的求解逻辑。把QP吃透，这些扩展只是多了一层外壳。
我在实际项目中的体会是：QP的难点不在数学，而在数值和工程。数学推导一下午就能看懂，但数值稳定性、求解器选型、热启动、代码生成这些工程细节，需要反复试错才能掌握。希望这篇内容能帮你少走一些弯路。如果你正在做MPC或者组合优化，遇到求解器不收敛的问题，先检查P矩阵的对称性和条件数，再检查约束缩放，最后看求解器精度设置。这三个地方能解决80%的QP求解问题。
