尧图网站建设 
写代码的人应该都有过这种经历目标函数明明是个漂亮的凸二次函数求个导、令它等于零几秒钟就能写完解析解的代码结果一旦加上约束解出来的点直接跑到不可行域外面去了。我在做小车MPC轨迹跟踪的时候就被这个问题卡过大半天——二次规划、约束条件、状态限幅都在算出来的控制量却让执行器饱和整条轨迹全是毛刺。后来把思路放到积极集法active set method上一步步看它怎么拆解约束、怎么切换工作集才算真正把二次规划的问题看透。这篇文章把二次规划、积极集法的核心原理、手算推演、可运行代码以及工程落地中常见的问题完整讲一遍适合刚碰带约束二次规划的学生也适合已经会调现成求解器、但想搞清楚黑盒内部到底在迭代什么的工程师。1. 先弄明白二次规划到底难在哪1.1 标准形式二次目标加线性约束二次规划的标准形式可以写成这样[ \min_x ; \frac{1}{2} x^T H x g^T x ] [ s.t. ; A x \le b ]如果问题里还有等式约束就再加一项 (A_{eq} x b_{eq})。这里 (H) 是Hessian矩阵凸二次规划要求 (H) 至少半正定也就是说目标函数是个开口向上的碗碗底是唯一的严格正定时。(g) 是线性项系数约束全部是线性的。这种形式在工程里出现频率高得惊人。最小二乘回归加约束就是一个典型例子比如 (y X\beta) 的回归问题加上 (\beta \ge 0) 或者 (\sum \beta_i 1) 之后目标函数展开就是 (0.5 \beta^T (X^T X) \beta - (X^T y)^T \beta const)其中 (H X^T X) 天然对称半正定(g -X^T y)约束是线性不等式。做投资组合优化时目标是让组合方差最小(H) 是协方差矩阵约束是权重和为1、权重非负这也是一个典型的凸二次规划。另一个我接触最多的是模型预测控制MPC每个控制周期都要解一个二次规划目标函数里含状态偏差和控制增量的惩罚约束是执行机构饱和幅度、状态量上下限本质上就是在每个采样时刻求一组最优控制增量。非线性规划领域里的SQP方法也是把非线性问题在每次迭代点处展开解一个二次规划子问题来得到搜索方向。可以说只要优化问题里出现二次目标 线性约束二次规划就出现了。1.2 求导失灵不等式约束带来的组合难题无约束二次规划的解法太简单了目标函数对 (x) 求导并令其等于零得到 (Hx g 0)一步得到解析解 (x^* -H^{-1}g)。如果这个解恰好落在可行域内那问题确实结束了。但绝大多数时候它不会那么听话。我早年犯过的错误就是拿到带约束的问题后先算无约束解再手动把越界的变量掰回边界上。简单问题这么糊弄还能勉强对上结果稍微复杂一点就崩。原因在于不等式约束的最优解并不满足 (\nabla f 0)它满足的是KKT条件[ Hx g A^T \lambda 0 ] [ \lambda_i \ge 0, \quad \lambda_i (b_i - a_i^T x) 0 ]其中 (\lambda) 是拉格朗日乘子(a_i^T) 是第 (i) 条约束的系数向量。互补条件 (\lambda_i (b_i - a_i^T x) 0) 的含义是如果第 (i) 条约束在最优解处没有触到边界即 (b_i - a_i^T x 0)那么它对应的乘子必须为零反过来如果约束生效了目标函数极值点上就一定会有该约束对应的推力这个推力就是乘子。难就难在哪些约束会生效这件事事先不知道。(m) 条约束里任意组合都可能成为有效约束这是一个 (2^m) 级的组合选择问题。无约束求导解决不了它因为它本质上不是连续优化问题里沿梯度走到极值的问题而是要在离散的约束激活模式里找出正确的一套。用个生活类比出门旅行收拾行李你不可能把所有季节的衣服都背上只能先按最可能的天气打包路上天气变了再往包里加衣服或减衣服。积极集法干的就是这个事。1.3 积极集法的破局思路主动猜出有效约束积极集法的核心思路很朴素既然最优解处只有一部分约束生效那我就维护一个我认为生效的约束集合叫工作集working set。在这个工作集下把不等式约束当成等式约束来处理解一个等式约束的二次规划子问题。解完之后做两个检查检查新点是否满足所有约束检查工作集约束对应的乘子是否非负。如果出了问题就调整工作集——要么加一条刚触界的约束要么踢掉一条乘子为负的约束——然后重新解子问题。每次迭代只调整工作集里的一个或少数几个约束逐步逼近真正的最优激活模式。这避免了直接面对 (2^m) 的组合爆炸问题理论上只要不退化有限步就能收敛。而且它有一个非常好的工程性质如果前后两个QP问题的约束激活模式差得不多上一次解完的工作集可以直接拿来当这次的初始工作集省掉大量重复迭代。这个性质在MPC里几乎是决定性的优点后面我会专门展开。2. 积极集法的两个核心判据触界与乘子2.1 工作集(active set)到底是什么在任意一个可行点 (x) 上把所有满足 (a_i^T x b_i) 的约束收集起来叫这个点的有效约束集。有效约束就是在当前点触到边界的约束。积极集法在迭代过程中维护的工作集本质上就是当前迭代点处认为应该生效的约束集合。为什么这个集合这么重要因为在最优点处真正起作用的只有这些约束。假设你已经知道了最终的最优有效约束集合 (W^*)那么原问题就退化成一个等式约束问题[ \min_x ; \frac{1}{2} x^T H x g^T x, \quad s.t. ; a_i^T x b_i, ; i \in W^* ]等式约束问题就好解多了直接拉格朗日乘子法转线性方程组。所以积极集法整个过程可以概括成不断猜测 (W)求解对应的等式约束子问题再用KKT条件验证猜测是否正确不正确就修正。2.2 等式约束子问题是怎么被解出来的假设当前迭代点是 (x_k)工作集是 (W)。我不直接解以 (x) 为变量的等式约束问题而是换成求解一个搜索方向 (p)[ \min_p ; \frac{1}{2} p^T H p (Hx_k g)^T p ] [ s.t. ; a_i^T p 0, ; i \in W ]这个转换的意义在于工作集内的约束在当前点已经被钉住了所以沿搜索方向走时这些约束必须保持等式成立即 (a_i^T (x_k \alpha p) b_i)等价于 (a_i^T p 0)。目标函数里 (p) 的线性项系数是当前点处的梯度 (Hx_k g)。对这个问题写拉格朗日函数对 (p) 和乘子 (\lambda) 求导置零会得到一个线性方程组[ \begin{bmatrix} H A_W^T \ A_W 0 \end{bmatrix} \begin{bmatrix} p \ \lambda \end{bmatrix}\begin{bmatrix} -(Hx_k g) \ 0 \end{bmatrix} ]其中 (A_W) 是工作集约束系数拼成的矩阵每一行是一条约束。只要 (H) 正定、工作集约束线性无关这个KKT矩阵就是非奇异的直接用高斯消去法或Cholesky分解就能解。这也是为什么约束线性无关在积极集法里是个前提条件实际工程中遇到约束线性相关时必须做处理否则矩阵奇异一步就报错。每次迭代的核心计算量就集中在这个线性方程组的求解上。2.3 什么时候加约束什么时候踢约束解出 (p) 之后分两种情况。第一种情况是 (p) 足够接近零向量。这说明在当前工作集下目标函数已经无法再改进了当前点是一个局部最优候选点。此时需要检查工作集内约束的拉格朗日乘子 (\lambda)。如果所有乘子都大于等于零数值上允许一定容差KKT条件全部满足当前点就是全局最优解。但如果出现了某个 (\lambda_j 0)说明第 (j) 条约束在当前点拉住了目标函数就像一只拖后腿的手把目标困在了一个并不是真正最小的位置。释放这条约束目标函数还能继续下降。所以要把乘子最负的那条约束从工作集里踢掉然后重新求解子问题。第二种情况是 (p) 不为零。(p) 是当前工作集下让目标下降最快的可行方向沿着它走一步目标一定下降。能走多远上限是1因为子问题沿这个方向的精确最优步长就是1再远目标就开始反弹。但在到达步长1之前很可能先碰到某条不属于工作集的约束的边界不能再往前走了。对所有不在工作集的约束算一下最大可行步长[ \alpha_{\max} \min_{i \notin W, \ a_i^T p 0} \frac{b_i - a_i^T x_k}{a_i^T p} ](a_i^T p 0) 表示这条约束的边界正在被逼近分母为正如果 (a_i^T p \le 0)表示方向在远离这条边界不需要限制。最终步长取 (\alpha \min(1, \alpha_{\max}))。如果 (\alpha) 被某条约束卡在小于1的位置这条约束就是新触界的拦截者要走过去就得把它加进工作集如果 (\alpha 1)那就直接走到子问题最优点工作集暂时不变继续下一轮迭代。3. 一个四步走通的手算例子3.1 问题、初始点与工作集光看算法流程还是抽象我拿一个二维问题完整推演一遍全部手算都能验算[ \min_x ; f(x) \frac{1}{2}(x_1^2 x_2^2) - x_1 - x_2 ]约束 [ x_1 x_2 \le 1, \quad x_1 \ge 0, \quad x_2 \ge 0 ]矩阵形式写出来就是 (H I)单位阵(g (-1, -1))。无约束极小点在 ((1,1))显然不满足 (x_1 x_2 \le 1)所以最优解一定落在约束边界上。几何上看可行域是三角形区域目标等值线是以 ((1,1)) 为圆心的同心圆可行域内离 ((1,1)) 最近的点就是 ((0.5, 0.5))这个点就是我们最终期望算到的结果。取初始点 (x_0 (0, 0.5))。这个点处 (x_1 0)所以约束 (x_1 \ge 0) 是生效的初始工作集设为 (W {x_1 \ge 0})。3.2 完整的迭代过程记录我按每一轮迭代把关键量列成表方便对照。轮次当前点 (x)工作集 (W)梯度 (Hxg)方向 (p)步长 (\alpha)动作0(0, 0.5){(x_1 \ge 0)}(-1, -0.5)(0, 0.5)1.0触碰到 (x_1x_2 \le 1)加入新约束1(0, 1){(x_1 \ge 0), (x_1x_2 \le 1)}(-1, 0)(0, 0)—乘子出现负值 (-1)踢掉 (x_1 \ge 0)2(0, 1){(x_1x_2 \le 1)}(-1, 0)(0.5, -0.5)1.0无新约束直接走到子问题最优点3(0.5, 0.5){(x_1x_2 \le 1)}(-0.5, -0.5)(0, 0)—乘子 (\lambda0.5 0)最优第0轮算出来 (p (0, 0.5))意思是保持 (x_10) 不动只把 (x_2) 往上推。它先碰到的约束是 (x_1x_2 \le 1)因为从 (x_20.5) 出发沿该方向走到边界要走的步长正好是1所以实际走了单位步长到达 ((0,1))并把这条约束加进工作集。这一步直观上也很好理解既然无约束最优点在 ((1,1))在当前点最想做的就是增大两个变量被边界拦住是很自然的事。第1轮是最关键的一轮。工作集里现在有 (x_10) 和 (x_1x_21) 两条约束当前点 ((0,1)) 是可行域的角点。解子问题得到 (p (0,0))说明在被这两条等式约束钉死的状态下任何方向都无法改善目标。但检查乘子时发现约束 (x_1 \ge 0) 对应的乘子 (\lambda -1)小于零。这个负乘子的几何含义等下一节细说操作上是把 (x_1 \ge 0) 这条约束从工作集里移除。第2轮有意思的地方在于移除约束后工作集里只剩 (x_1x_2 \le 1)从 ((0,1)) 出发的方向变成 ((0.5, -0.5))。这正好是沿着斜边边界 (x_1x_21) 往下滑的方向。步长计算时(x_2 \ge 0) 这条约束允许走2步但单位步长1更小所以 (\alpha 1)新点落在 ((0.5, 0.5))没有触碰新边界工作集保持不变。第3轮验证收尾在 ((0.5, 0.5)) 处再解一次子问题得到 (p(0,0))乘子 (\lambda 0.5 0)KKT条件全部满足计算结束。3.3 每一步背后在干什么第1轮的负乘子是最容易困惑的地方。直观解释是这样的当前点 ((0,1)) 同时被两条约束钉住其中 (x_10) 这条约束其实在帮倒忙。如果把这条约束松绑允许 (x_1) 稍微增大一点目标函数沿哪个方向能下降沿着斜边 (x_1x_21) 朝 ((0.5,0.5)) 移动目标从 (f(0,1) -0.5) 降到 (f(0.5,0.5) -0.75)确实还有改进空间。乘子为负就是算法用数学语言告诉你这条约束在当前点不是支撑最优解而是限制你得过头了。第2轮的方向 ((0.5,-0.5)) 也有清晰几何意义。((0.5,-0.5)) 正是无约束最优点 ((1,1)) 到斜边垂线的投影方向所以沿这个方向走单位步长恰好落到最优投影点。而且这轮里虽然 (x_2 \ge 0) 约束还允许继续走但单位步长已经到极限因为子问题里沿 (p) 方向的最优步长就是1再往下走目标函数会开始上升所以停在这里是合理的。还有一点值得体会第1轮和第3轮都得到 (p0)但一轮判定继续迭代、一轮判定收敛差异完全来自乘子符号。这说明积极集法里求解子问题只是半程检查乘子才是判断是否到终点的关键。初学者最容易漏掉乘子检查结果在 (p0) 时直接认为已经最优然后拿着错误答案去调试怎么都想不通。4. 手写一个能跑的Python实现4.1 核心子程序解KKT系统前面的推导已经说明每次迭代的核心是解一个KKT线性方程组。我用NumPy写一个最小可用的实现不求性能最优但逻辑清晰、可以直接改来用。import numpy as np def solve_kkt(H, grad, A_w): 解等式约束子问题: min 0.5 * p H p grad p s.t. A_w p 0 返回 (p, lambda) n H.shape[0] m A_w.shape[0] KKT np.block([ [H, A_w.T], [A_w, np.zeros((m, m))] ]) rhs np.concatenate([-grad, np.zeros(m)]) sol np.linalg.solve(KKT, rhs) return sol[:n], sol[n:]这个函数的输入是当前迭代点处的梯度 (grad Hx g) 和工作集约束矩阵 (A_w)。返回值中 (p) 是搜索方向(\lambda) 是对应工作集约束的拉格朗日乘子。注意到我直接用np.linalg.solve要求KKT矩阵非奇异这就是前面说的约束线性无关条件。严格来说在约束工作集里可能混入冗余约束导致奇异完善的实现要先做秩检查这里先不展开。4.2 主循环步长、触界、乘子检查主循环按之前描述的算法流程组织解子问题、判断 (p) 是否为零、计算步长、更新工作集。def active_set_qp(H, g, A, b, x0, W0None, tol1e-9, max_iter200): x x0.copy() m A.shape[0] # 初始工作集默认取所有在当前点触界的约束 if W0 is None: W0 [i for i in range(m) if abs(A[i] x - b[i]) tol] W list(W0) for _ in range(max_iter): grad H x g A_w A[W, :] p, lam solve_kkt(H, grad, A_w) # 情况1: 方向接近零检查乘子 if np.linalg.norm(p) tol: if len(W) 0 and lam.min() -tol: j_local int(np.argmin(lam)) j W.pop(j_local) continue else: return x, W # 情况2: 方向非零计算最大可行步长 alpha 1.0 for i in range(m): if i in W: continue den A[i] p if den tol: cand (b[i] - A[i] x) / den alpha min(alpha, cand) x_new x alpha * p # 新点处触界的约束加入工作集 touched [ i for i in range(m) if i not in W and abs(A[i] x_new - b[i]) tol ] for i in touched: if i not in W: W.append(i) x x_new raise RuntimeError(达到最大迭代次数可能出现了退化或数值问题)这里有一个实现细节值得说明我在计算步长时对所有不在工作集的约束都检查一遍分母 (a_i^T p)只有为正时才可能构成限制。同时在更新点之后统一检查哪些约束在新点处触界一次性加入工作集这样即使步长为1且有约束在终点恰好被碰到也不会漏掉。4.3 用刚才的例子验证代码用第3节手算的例子做测试H np.array([[1.0, 0.0], [0.0, 1.0]]) g np.array([-1.0, -1.0]) # 约束顺序x1x21, x10, x20 A np.array([[1.0, 1.0], [-1.0, 0.0], [0.0, -1.0]]) b np.array([1.0, 0.0, 0.0]) x0 np.array([0.0, 0.5]) x_opt, W_opt active_set_qp(H, g, A, b, x0, W0[1]) print(最优解:, x_opt) print(最终工作集:, W_opt)在我的环境里运行输出是最优解: [0.5 0.5] 最终工作集: [0]恰好对应手算推演的结果(x^* (0.5, 0.5))最终生效的约束只有 (x_1x_2 \le 1)。代码和手算互相验证说明整个实现逻辑是自洽的。5. 工程落地的实战经验5.1 热启动是积极集法的杀手锏真正在工程里用积极集法最大的优势不是它原理简单而是热启动warm start能力。在做MPC时相邻两个控制周期求解的QP问题在结构上几乎一样只是目标函数里的参考轨迹和测量值变了而实际的约束激活模式往往变化很小。上一周期算完得到的最优解和最优点处的有效约束集合直接作为下一周期的初始点 (x_0) 和初始工作集 (W_0)通常只需要几次迭代就能收敛。这个性质对嵌入式系统非常友好每次迭代的耗时波动小最坏情况可控不像某些算法每轮迭代内部还有不定次数的回溯搜索。我在一个算力受限的控制器上把QP求解时间从几毫秒压到几百微秒靠的就是热启动加提前把KKT矩阵的分解结果缓存起来。相比之下内点法虽然迭代次数与规模关系不大但每次迭代都要解一个更大的线性系统而且暖启动的优势没有积极集法这么直接。5.2 退化循环与数值容差积极集法最经典的理论问题叫退化degeneracy表现出来就是算法在一个点附近反复加约束、踢约束陷入死循环。最常见的原因是多个约束在同一个点同时触界或者约束之间存在近似线性相关导致乘子判定时出现平局工作集在几条约束之间反复横跳。我在写数值实验时第一次遇到这个问题卡在循环里出不来打印工作集轨迹才发现某条约束被反复加进来又被踢出去。工程上的应对手段通常有几个层次。第一约束预处理删掉重复或近似重复的约束把变量和约束都做缩放让数值在同一个量级。第二判断规则上做文章比如乘子为负时总选最负的那个或者加一个很小的随机扰动打破平局。第三设置最大迭代次数的兜底一旦超限就触发一个降级策略。另外(H) 只是半正定而非严格正定时KKT矩阵可能奇异常见的处理是给Hessian加一个很小的正则化项 (H \epsilon I)让问题变成严格凸代价是解的精度会受影响但稳定很多。5.3 什么时候该换内点法或ADMM积极集法不是万能的。它的迭代次数和约束切换次数强相关如果问题规模大、约束模式变动剧烈或者矩阵结构本身适合做稀疏分解内点法和ADMM往往更有优势。我整理了一个简单的选型对比直接给出结论性的经验维度积极集法内点法ADMM类如OSQP典型规模中小规模稠密问题大规模稀疏问题大规模稀疏问题热启动能力极强工作集可直接复用一般复用受限中等可复用残差信息迭代次数与约束切换次数相关基本固定几十轮几十到几百轮取决于精度精度高直接解线性系统高精度需要更多迭代中精度受容差限制使用场景MPC、SQP内层、小规模组合优化通用大规模凸优化大规模MPC、稀疏QP比如在金融组合优化里资产数量几千、协方差矩阵稀疏度一般用内点法或OSQP更常见。而传统MPC领域因为每个周期的问题小且热启动收益巨大qpOASES这类积极集法求解器反而更流行。OSQP虽然基于ADMM但在大规模稀疏场景下极快也是现代MPC的热门选择。5.4 对接现成库与两个调试技巧如果不想自己实现常见的现成库按接口易用度排序quadprog是Goldfarb-Idnani型积极集法适合中小规模稠密问题接口极简qpOASES专门为MPC的在线优化设计支持热启动和动态问题OSQP基于ADMM适合大规模稀疏问题scipy.optimize.minimize里的SLSQP和trust-constr也能解但性能一般适合原型验证。前提是确认库的底层算法和你对问题的预期匹配。调试积极集法有两个技巧我屡试不爽。第一个是打印每一轮的工作集、(p)、(\lambda)、(\alpha)特别是观察约束索引的序列是否在反复出现同一组索引——如果在基本就是退化循环。第二个是调试前先把约束做归一化处理把每条约束的 (a_i) 除以它的范数这样步长和乘子的数值尺度一致容差判断才不容易出偏差。我见过很多看起来约束没问题但求解器老报不可行的案例最后查出来是变量单位不一致毫米和米混用数值差了一千倍导致容差判断全部失效。我一直觉得积极集法虽然是个老方法但它把约束激活模式的离散决策问题讲得最透彻。你一旦手写过一次理解了KKT条件里互补性条件到底在说什么后面再用内点法、ADMM或者其他任何先进求解器遇到原理性问题时心里都有底。放在现在的嵌入式MPC项目里我依然首选带热启动的积极集法原因只有一个它在实际问题上就是稳定、快速、可控。 二次规划与积极集法：原理、实现与工程实战
写代码的人应该都有过这种经历目标函数明明是个漂亮的凸二次函数求个导、令它等于零几秒钟就能写完解析解的代码结果一旦加上约束解出来的点直接跑到不可行域外面去了。我在做小车MPC轨迹跟踪的时候就被这个问题卡过大半天——二次规划、约束条件、状态限幅都在算出来的控制量却让执行器饱和整条轨迹全是毛刺。后来把思路放到积极集法active set method上一步步看它怎么拆解约束、怎么切换工作集才算真正把二次规划的问题看透。这篇文章把二次规划、积极集法的核心原理、手算推演、可运行代码以及工程落地中常见的问题完整讲一遍适合刚碰带约束二次规划的学生也适合已经会调现成求解器、但想搞清楚黑盒内部到底在迭代什么的工程师。1. 先弄明白二次规划到底难在哪1.1 标准形式二次目标加线性约束二次规划的标准形式可以写成这样[ \min_x ; \frac{1}{2} x^T H x g^T x ] [ s.t. ; A x \le b ]如果问题里还有等式约束就再加一项 (A_{eq} x b_{eq})。这里 (H) 是Hessian矩阵凸二次规划要求 (H) 至少半正定也就是说目标函数是个开口向上的碗碗底是唯一的严格正定时。(g) 是线性项系数约束全部是线性的。这种形式在工程里出现频率高得惊人。最小二乘回归加约束就是一个典型例子比如 (y X\beta) 的回归问题加上 (\beta \ge 0) 或者 (\sum \beta_i 1) 之后目标函数展开就是 (0.5 \beta^T (X^T X) \beta - (X^T y)^T \beta const)其中 (H X^T X) 天然对称半正定(g -X^T y)约束是线性不等式。做投资组合优化时目标是让组合方差最小(H) 是协方差矩阵约束是权重和为1、权重非负这也是一个典型的凸二次规划。另一个我接触最多的是模型预测控制MPC每个控制周期都要解一个二次规划目标函数里含状态偏差和控制增量的惩罚约束是执行机构饱和幅度、状态量上下限本质上就是在每个采样时刻求一组最优控制增量。非线性规划领域里的SQP方法也是把非线性问题在每次迭代点处展开解一个二次规划子问题来得到搜索方向。可以说只要优化问题里出现二次目标 线性约束二次规划就出现了。1.2 求导失灵不等式约束带来的组合难题无约束二次规划的解法太简单了目标函数对 (x) 求导并令其等于零得到 (Hx g 0)一步得到解析解 (x^* -H^{-1}g)。如果这个解恰好落在可行域内那问题确实结束了。但绝大多数时候它不会那么听话。我早年犯过的错误就是拿到带约束的问题后先算无约束解再手动把越界的变量掰回边界上。简单问题这么糊弄还能勉强对上结果稍微复杂一点就崩。原因在于不等式约束的最优解并不满足 (\nabla f 0)它满足的是KKT条件[ Hx g A^T \lambda 0 ] [ \lambda_i \ge 0, \quad \lambda_i (b_i - a_i^T x) 0 ]其中 (\lambda) 是拉格朗日乘子(a_i^T) 是第 (i) 条约束的系数向量。互补条件 (\lambda_i (b_i - a_i^T x) 0) 的含义是如果第 (i) 条约束在最优解处没有触到边界即 (b_i - a_i^T x 0)那么它对应的乘子必须为零反过来如果约束生效了目标函数极值点上就一定会有该约束对应的推力这个推力就是乘子。难就难在哪些约束会生效这件事事先不知道。(m) 条约束里任意组合都可能成为有效约束这是一个 (2^m) 级的组合选择问题。无约束求导解决不了它因为它本质上不是连续优化问题里沿梯度走到极值的问题而是要在离散的约束激活模式里找出正确的一套。用个生活类比出门旅行收拾行李你不可能把所有季节的衣服都背上只能先按最可能的天气打包路上天气变了再往包里加衣服或减衣服。积极集法干的就是这个事。1.3 积极集法的破局思路主动猜出有效约束积极集法的核心思路很朴素既然最优解处只有一部分约束生效那我就维护一个我认为生效的约束集合叫工作集working set。在这个工作集下把不等式约束当成等式约束来处理解一个等式约束的二次规划子问题。解完之后做两个检查检查新点是否满足所有约束检查工作集约束对应的乘子是否非负。如果出了问题就调整工作集——要么加一条刚触界的约束要么踢掉一条乘子为负的约束——然后重新解子问题。每次迭代只调整工作集里的一个或少数几个约束逐步逼近真正的最优激活模式。这避免了直接面对 (2^m) 的组合爆炸问题理论上只要不退化有限步就能收敛。而且它有一个非常好的工程性质如果前后两个QP问题的约束激活模式差得不多上一次解完的工作集可以直接拿来当这次的初始工作集省掉大量重复迭代。这个性质在MPC里几乎是决定性的优点后面我会专门展开。2. 积极集法的两个核心判据触界与乘子2.1 工作集(active set)到底是什么在任意一个可行点 (x) 上把所有满足 (a_i^T x b_i) 的约束收集起来叫这个点的有效约束集。有效约束就是在当前点触到边界的约束。积极集法在迭代过程中维护的工作集本质上就是当前迭代点处认为应该生效的约束集合。为什么这个集合这么重要因为在最优点处真正起作用的只有这些约束。假设你已经知道了最终的最优有效约束集合 (W^*)那么原问题就退化成一个等式约束问题[ \min_x ; \frac{1}{2} x^T H x g^T x, \quad s.t. ; a_i^T x b_i, ; i \in W^* ]等式约束问题就好解多了直接拉格朗日乘子法转线性方程组。所以积极集法整个过程可以概括成不断猜测 (W)求解对应的等式约束子问题再用KKT条件验证猜测是否正确不正确就修正。2.2 等式约束子问题是怎么被解出来的假设当前迭代点是 (x_k)工作集是 (W)。我不直接解以 (x) 为变量的等式约束问题而是换成求解一个搜索方向 (p)[ \min_p ; \frac{1}{2} p^T H p (Hx_k g)^T p ] [ s.t. ; a_i^T p 0, ; i \in W ]这个转换的意义在于工作集内的约束在当前点已经被钉住了所以沿搜索方向走时这些约束必须保持等式成立即 (a_i^T (x_k \alpha p) b_i)等价于 (a_i^T p 0)。目标函数里 (p) 的线性项系数是当前点处的梯度 (Hx_k g)。对这个问题写拉格朗日函数对 (p) 和乘子 (\lambda) 求导置零会得到一个线性方程组[ \begin{bmatrix} H A_W^T \ A_W 0 \end{bmatrix} \begin{bmatrix} p \ \lambda \end{bmatrix}\begin{bmatrix} -(Hx_k g) \ 0 \end{bmatrix} ]其中 (A_W) 是工作集约束系数拼成的矩阵每一行是一条约束。只要 (H) 正定、工作集约束线性无关这个KKT矩阵就是非奇异的直接用高斯消去法或Cholesky分解就能解。这也是为什么约束线性无关在积极集法里是个前提条件实际工程中遇到约束线性相关时必须做处理否则矩阵奇异一步就报错。每次迭代的核心计算量就集中在这个线性方程组的求解上。2.3 什么时候加约束什么时候踢约束解出 (p) 之后分两种情况。第一种情况是 (p) 足够接近零向量。这说明在当前工作集下目标函数已经无法再改进了当前点是一个局部最优候选点。此时需要检查工作集内约束的拉格朗日乘子 (\lambda)。如果所有乘子都大于等于零数值上允许一定容差KKT条件全部满足当前点就是全局最优解。但如果出现了某个 (\lambda_j 0)说明第 (j) 条约束在当前点拉住了目标函数就像一只拖后腿的手把目标困在了一个并不是真正最小的位置。释放这条约束目标函数还能继续下降。所以要把乘子最负的那条约束从工作集里踢掉然后重新求解子问题。第二种情况是 (p) 不为零。(p) 是当前工作集下让目标下降最快的可行方向沿着它走一步目标一定下降。能走多远上限是1因为子问题沿这个方向的精确最优步长就是1再远目标就开始反弹。但在到达步长1之前很可能先碰到某条不属于工作集的约束的边界不能再往前走了。对所有不在工作集的约束算一下最大可行步长[ \alpha_{\max} \min_{i \notin W, \ a_i^T p 0} \frac{b_i - a_i^T x_k}{a_i^T p} ](a_i^T p 0) 表示这条约束的边界正在被逼近分母为正如果 (a_i^T p \le 0)表示方向在远离这条边界不需要限制。最终步长取 (\alpha \min(1, \alpha_{\max}))。如果 (\alpha) 被某条约束卡在小于1的位置这条约束就是新触界的拦截者要走过去就得把它加进工作集如果 (\alpha 1)那就直接走到子问题最优点工作集暂时不变继续下一轮迭代。3. 一个四步走通的手算例子3.1 问题、初始点与工作集光看算法流程还是抽象我拿一个二维问题完整推演一遍全部手算都能验算[ \min_x ; f(x) \frac{1}{2}(x_1^2 x_2^2) - x_1 - x_2 ]约束 [ x_1 x_2 \le 1, \quad x_1 \ge 0, \quad x_2 \ge 0 ]矩阵形式写出来就是 (H I)单位阵(g (-1, -1))。无约束极小点在 ((1,1))显然不满足 (x_1 x_2 \le 1)所以最优解一定落在约束边界上。几何上看可行域是三角形区域目标等值线是以 ((1,1)) 为圆心的同心圆可行域内离 ((1,1)) 最近的点就是 ((0.5, 0.5))这个点就是我们最终期望算到的结果。取初始点 (x_0 (0, 0.5))。这个点处 (x_1 0)所以约束 (x_1 \ge 0) 是生效的初始工作集设为 (W {x_1 \ge 0})。3.2 完整的迭代过程记录我按每一轮迭代把关键量列成表方便对照。轮次当前点 (x)工作集 (W)梯度 (Hxg)方向 (p)步长 (\alpha)动作0(0, 0.5){(x_1 \ge 0)}(-1, -0.5)(0, 0.5)1.0触碰到 (x_1x_2 \le 1)加入新约束1(0, 1){(x_1 \ge 0), (x_1x_2 \le 1)}(-1, 0)(0, 0)—乘子出现负值 (-1)踢掉 (x_1 \ge 0)2(0, 1){(x_1x_2 \le 1)}(-1, 0)(0.5, -0.5)1.0无新约束直接走到子问题最优点3(0.5, 0.5){(x_1x_2 \le 1)}(-0.5, -0.5)(0, 0)—乘子 (\lambda0.5 0)最优第0轮算出来 (p (0, 0.5))意思是保持 (x_10) 不动只把 (x_2) 往上推。它先碰到的约束是 (x_1x_2 \le 1)因为从 (x_20.5) 出发沿该方向走到边界要走的步长正好是1所以实际走了单位步长到达 ((0,1))并把这条约束加进工作集。这一步直观上也很好理解既然无约束最优点在 ((1,1))在当前点最想做的就是增大两个变量被边界拦住是很自然的事。第1轮是最关键的一轮。工作集里现在有 (x_10) 和 (x_1x_21) 两条约束当前点 ((0,1)) 是可行域的角点。解子问题得到 (p (0,0))说明在被这两条等式约束钉死的状态下任何方向都无法改善目标。但检查乘子时发现约束 (x_1 \ge 0) 对应的乘子 (\lambda -1)小于零。这个负乘子的几何含义等下一节细说操作上是把 (x_1 \ge 0) 这条约束从工作集里移除。第2轮有意思的地方在于移除约束后工作集里只剩 (x_1x_2 \le 1)从 ((0,1)) 出发的方向变成 ((0.5, -0.5))。这正好是沿着斜边边界 (x_1x_21) 往下滑的方向。步长计算时(x_2 \ge 0) 这条约束允许走2步但单位步长1更小所以 (\alpha 1)新点落在 ((0.5, 0.5))没有触碰新边界工作集保持不变。第3轮验证收尾在 ((0.5, 0.5)) 处再解一次子问题得到 (p(0,0))乘子 (\lambda 0.5 0)KKT条件全部满足计算结束。3.3 每一步背后在干什么第1轮的负乘子是最容易困惑的地方。直观解释是这样的当前点 ((0,1)) 同时被两条约束钉住其中 (x_10) 这条约束其实在帮倒忙。如果把这条约束松绑允许 (x_1) 稍微增大一点目标函数沿哪个方向能下降沿着斜边 (x_1x_21) 朝 ((0.5,0.5)) 移动目标从 (f(0,1) -0.5) 降到 (f(0.5,0.5) -0.75)确实还有改进空间。乘子为负就是算法用数学语言告诉你这条约束在当前点不是支撑最优解而是限制你得过头了。第2轮的方向 ((0.5,-0.5)) 也有清晰几何意义。((0.5,-0.5)) 正是无约束最优点 ((1,1)) 到斜边垂线的投影方向所以沿这个方向走单位步长恰好落到最优投影点。而且这轮里虽然 (x_2 \ge 0) 约束还允许继续走但单位步长已经到极限因为子问题里沿 (p) 方向的最优步长就是1再往下走目标函数会开始上升所以停在这里是合理的。还有一点值得体会第1轮和第3轮都得到 (p0)但一轮判定继续迭代、一轮判定收敛差异完全来自乘子符号。这说明积极集法里求解子问题只是半程检查乘子才是判断是否到终点的关键。初学者最容易漏掉乘子检查结果在 (p0) 时直接认为已经最优然后拿着错误答案去调试怎么都想不通。4. 手写一个能跑的Python实现4.1 核心子程序解KKT系统前面的推导已经说明每次迭代的核心是解一个KKT线性方程组。我用NumPy写一个最小可用的实现不求性能最优但逻辑清晰、可以直接改来用。import numpy as np def solve_kkt(H, grad, A_w): 解等式约束子问题: min 0.5 * p H p grad p s.t. A_w p 0 返回 (p, lambda) n H.shape[0] m A_w.shape[0] KKT np.block([ [H, A_w.T], [A_w, np.zeros((m, m))] ]) rhs np.concatenate([-grad, np.zeros(m)]) sol np.linalg.solve(KKT, rhs) return sol[:n], sol[n:]这个函数的输入是当前迭代点处的梯度 (grad Hx g) 和工作集约束矩阵 (A_w)。返回值中 (p) 是搜索方向(\lambda) 是对应工作集约束的拉格朗日乘子。注意到我直接用np.linalg.solve要求KKT矩阵非奇异这就是前面说的约束线性无关条件。严格来说在约束工作集里可能混入冗余约束导致奇异完善的实现要先做秩检查这里先不展开。4.2 主循环步长、触界、乘子检查主循环按之前描述的算法流程组织解子问题、判断 (p) 是否为零、计算步长、更新工作集。def active_set_qp(H, g, A, b, x0, W0None, tol1e-9, max_iter200): x x0.copy() m A.shape[0] # 初始工作集默认取所有在当前点触界的约束 if W0 is None: W0 [i for i in range(m) if abs(A[i] x - b[i]) tol] W list(W0) for _ in range(max_iter): grad H x g A_w A[W, :] p, lam solve_kkt(H, grad, A_w) # 情况1: 方向接近零检查乘子 if np.linalg.norm(p) tol: if len(W) 0 and lam.min() -tol: j_local int(np.argmin(lam)) j W.pop(j_local) continue else: return x, W # 情况2: 方向非零计算最大可行步长 alpha 1.0 for i in range(m): if i in W: continue den A[i] p if den tol: cand (b[i] - A[i] x) / den alpha min(alpha, cand) x_new x alpha * p # 新点处触界的约束加入工作集 touched [ i for i in range(m) if i not in W and abs(A[i] x_new - b[i]) tol ] for i in touched: if i not in W: W.append(i) x x_new raise RuntimeError(达到最大迭代次数可能出现了退化或数值问题)这里有一个实现细节值得说明我在计算步长时对所有不在工作集的约束都检查一遍分母 (a_i^T p)只有为正时才可能构成限制。同时在更新点之后统一检查哪些约束在新点处触界一次性加入工作集这样即使步长为1且有约束在终点恰好被碰到也不会漏掉。4.3 用刚才的例子验证代码用第3节手算的例子做测试H np.array([[1.0, 0.0], [0.0, 1.0]]) g np.array([-1.0, -1.0]) # 约束顺序x1x21, x10, x20 A np.array([[1.0, 1.0], [-1.0, 0.0], [0.0, -1.0]]) b np.array([1.0, 0.0, 0.0]) x0 np.array([0.0, 0.5]) x_opt, W_opt active_set_qp(H, g, A, b, x0, W0[1]) print(最优解:, x_opt) print(最终工作集:, W_opt)在我的环境里运行输出是最优解: [0.5 0.5] 最终工作集: [0]恰好对应手算推演的结果(x^* (0.5, 0.5))最终生效的约束只有 (x_1x_2 \le 1)。代码和手算互相验证说明整个实现逻辑是自洽的。5. 工程落地的实战经验5.1 热启动是积极集法的杀手锏真正在工程里用积极集法最大的优势不是它原理简单而是热启动warm start能力。在做MPC时相邻两个控制周期求解的QP问题在结构上几乎一样只是目标函数里的参考轨迹和测量值变了而实际的约束激活模式往往变化很小。上一周期算完得到的最优解和最优点处的有效约束集合直接作为下一周期的初始点 (x_0) 和初始工作集 (W_0)通常只需要几次迭代就能收敛。这个性质对嵌入式系统非常友好每次迭代的耗时波动小最坏情况可控不像某些算法每轮迭代内部还有不定次数的回溯搜索。我在一个算力受限的控制器上把QP求解时间从几毫秒压到几百微秒靠的就是热启动加提前把KKT矩阵的分解结果缓存起来。相比之下内点法虽然迭代次数与规模关系不大但每次迭代都要解一个更大的线性系统而且暖启动的优势没有积极集法这么直接。5.2 退化循环与数值容差积极集法最经典的理论问题叫退化degeneracy表现出来就是算法在一个点附近反复加约束、踢约束陷入死循环。最常见的原因是多个约束在同一个点同时触界或者约束之间存在近似线性相关导致乘子判定时出现平局工作集在几条约束之间反复横跳。我在写数值实验时第一次遇到这个问题卡在循环里出不来打印工作集轨迹才发现某条约束被反复加进来又被踢出去。工程上的应对手段通常有几个层次。第一约束预处理删掉重复或近似重复的约束把变量和约束都做缩放让数值在同一个量级。第二判断规则上做文章比如乘子为负时总选最负的那个或者加一个很小的随机扰动打破平局。第三设置最大迭代次数的兜底一旦超限就触发一个降级策略。另外(H) 只是半正定而非严格正定时KKT矩阵可能奇异常见的处理是给Hessian加一个很小的正则化项 (H \epsilon I)让问题变成严格凸代价是解的精度会受影响但稳定很多。5.3 什么时候该换内点法或ADMM积极集法不是万能的。它的迭代次数和约束切换次数强相关如果问题规模大、约束模式变动剧烈或者矩阵结构本身适合做稀疏分解内点法和ADMM往往更有优势。我整理了一个简单的选型对比直接给出结论性的经验维度积极集法内点法ADMM类如OSQP典型规模中小规模稠密问题大规模稀疏问题大规模稀疏问题热启动能力极强工作集可直接复用一般复用受限中等可复用残差信息迭代次数与约束切换次数相关基本固定几十轮几十到几百轮取决于精度精度高直接解线性系统高精度需要更多迭代中精度受容差限制使用场景MPC、SQP内层、小规模组合优化通用大规模凸优化大规模MPC、稀疏QP比如在金融组合优化里资产数量几千、协方差矩阵稀疏度一般用内点法或OSQP更常见。而传统MPC领域因为每个周期的问题小且热启动收益巨大qpOASES这类积极集法求解器反而更流行。OSQP虽然基于ADMM但在大规模稀疏场景下极快也是现代MPC的热门选择。5.4 对接现成库与两个调试技巧如果不想自己实现常见的现成库按接口易用度排序quadprog是Goldfarb-Idnani型积极集法适合中小规模稠密问题接口极简qpOASES专门为MPC的在线优化设计支持热启动和动态问题OSQP基于ADMM适合大规模稀疏问题scipy.optimize.minimize里的SLSQP和trust-constr也能解但性能一般适合原型验证。前提是确认库的底层算法和你对问题的预期匹配。调试积极集法有两个技巧我屡试不爽。第一个是打印每一轮的工作集、(p)、(\lambda)、(\alpha)特别是观察约束索引的序列是否在反复出现同一组索引——如果在基本就是退化循环。第二个是调试前先把约束做归一化处理把每条约束的 (a_i) 除以它的范数这样步长和乘子的数值尺度一致容差判断才不容易出偏差。我见过很多看起来约束没问题但求解器老报不可行的案例最后查出来是变量单位不一致毫米和米混用数值差了一千倍导致容差判断全部失效。我一直觉得积极集法虽然是个老方法但它把约束激活模式的离散决策问题讲得最透彻。你一旦手写过一次理解了KKT条件里互补性条件到底在说什么后面再用内点法、ADMM或者其他任何先进求解器遇到原理性问题时心里都有底。放在现在的嵌入式MPC项目里我依然首选带热启动的积极集法原因只有一个它在实际问题上就是稳定、快速、可控。