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 * Iepsilon取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条件最优解的“身份证”对于凸QPKKT条件是最优解的充要条件。标准形式的KKT条件包括平稳性P * x q A^T * lambda A_eq^T * nu 0原始可行性A * x bA_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 对偶问题换个角度看QPQP的对偶问题也是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 怎么选看问题规模和实时性要求特性有效集法内点法小规模问题n100快迭代少相对慢大规模稀疏问题稀疏性易破坏稀疏性保持好热启动非常方便不方便退化问题可能循环稳定返回活跃集直接返回需要额外计算实时控制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-3Hessian矩阵的对角元素相差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个变量的QPP是稀疏的每行只有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优秀优秀大规模稀疏、MPCquadprog有效集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在每个时间步解一个QPmin sum_{k0}^{N-1} (x_k^T Q x_k u_k^T R u_k) x_N^T Qf x_N s.t. x_{k1} 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求解问题。