Python实战：用有效集法解决不等式约束二次规划问题（附完整代码）

Python实战用有效集法解决不等式约束二次规划问题附完整代码二次规划问题在金融投资组合优化、工程控制系统设计等领域应用广泛。当问题带有不等式约束时有效集法Active Set Method因其稳定性和高效性成为首选算法之一。本文将手把手带您实现该算法的Python版本从理论到代码落地解决实际工程问题。1. 理解不等式约束二次规划的核心二次规划(QP)的标准形式可表示为\min \frac{1}{2}x^TQx c^Tx \\ s.t. \quad Ax \leq b \\ \quad Ex d其中Q是正定矩阵A和E分别表示不等式和等式约束矩阵。有效集法的精髓在于识别哪些不等式约束在最优解处会激活即变为等式约束。这相当于在迭代过程中动态调整约束条件初始时选择部分约束作为有效集每次迭代求解等式约束QP子问题根据结果调整有效集成员实际项目中常遇到的典型场景包括投资组合权重限制单资产占比不超过15%机器人运动规划的关节角度约束供应链优化的库存容量限制2. 算法实现的关键步骤拆解2.1 初始化可行解寻找初始可行点是算法第一步。简单问题可人工指定复杂情况可采用两阶段法def find_initial_point(A, b, ENone, dNone): # 构造辅助线性规划问题 from scipy.optimize import linprog res linprog(cnp.ones(A.shape[1]), A_ubA, b_ubb, A_eqE, b_eqd) if not res.success: raise ValueError(无可行解!) return res.x2.2 等式约束QP求解器有效集法的核心子问题是求解等式约束QP。这里采用变量消去法实现def solve_equality_qp(Q, c, A, b): 求解 min 0.5xQx cx s.t. Axb 返回最优解和拉格朗日乘子 n Q.shape[0] m A.shape[0] # 构造KKT系统 KKT np.block([ [Q, A.T], [A, np.zeros((m,m))] ]) rhs np.concatenate([-c, b]) solution np.linalg.solve(KKT, rhs) x solution[:n] lam solution[n:] return x, lam2.3 有效集更新逻辑根据迭代结果动态调整有效集是算法核心情况判断条件处理方式最优解在有效集内x_k x* 且 λ ≥ 0找到全局最优需移除约束x_k x* 但存在λ 0移除对应约束需添加约束x* 在可行域内添加新激活约束需限制步长x* 超出可行域计算最大可行步长对应Python实现def update_active_set(x_new, x_curr, lam, active_set, A, b): # 情况1找到最优解 if np.allclose(x_new, x_curr) and np.all(lam -1e-6): return x_new, active_set, True # 情况2需要移除约束 elif np.allclose(x_new, x_curr): min_idx np.argmin(lam) new_active [i for i in active_set if i ! min_idx] return x_new, new_active, False # 情况3需要添加约束 elif is_feasible(x_new, A, b): new_active get_active_constraints(x_new, A, b) return x_new, new_active, False # 情况4需要限制步长 else: alpha compute_max_step(x_curr, x_new, A, b, active_set) x_updated x_curr alpha * (x_new - x_curr) new_active get_active_constraints(x_updated, A, b) return x_updated, new_active, False3. 完整算法实现与性能优化将各模块组合成完整算法def active_set_qp(Q, c, A, b, ENone, dNone, max_iter100): # 初始化 x find_initial_point(A, b, E, d) active_set get_active_constraints(x, A, b) converged False for _ in range(max_iter): # 构造当前有效约束 A_active A[active_set] b_active b[active_set] # 求解等式约束子问题 x_new, lam solve_equality_qp(Q, c, A_active, b_active) # 更新有效集 x, active_set, converged update_active_set( x_new, x, lam, active_set, A, b) if converged: break return x性能优化技巧使用稀疏矩阵存储Q和A当n1000时采用热启动策略warm start加速后续求解对KKT系统使用Cholesky分解而非直接求逆4. 实战案例投资组合优化假设需要优化包含5只股票的投资组合# 收益率协方差矩阵 Q np.array([ [0.2, 0.05, 0.03, 0.01, 0.02], [0.05, 0.3, 0.02, 0.01, 0.03], [0.03, 0.02, 0.15, 0.02, 0.01], [0.01, 0.01, 0.02, 0.1, 0.01], [0.02, 0.03, 0.01, 0.01, 0.25] ]) # 预期收益率求最小方差组合 c np.zeros(5) # 约束条件权重和为1单只股票不超过30% A np.vstack([ np.ones(5), # 权重和1 -np.ones(5), # 权重和1等式约束拆解 np.eye(5), # 每只股票≥0 -np.eye(5) # 每只股票≤0.3 ]) b np.array([1, -1] [0]*5 [0.3]*5) # 求解 weights active_set_qp(Q, c, A, b) print(最优权重分配:, np.round(weights, 3))典型输出结果最优权重分配: [0.3 0.12 0.3 0.18 0.1]调试建议检查初始点可行性监控有效集变化情况验证KKT条件是否满足当算法震荡时适当增加收敛容差5. 与其他QP求解器的对比下表比较了不同QP求解方法的特性方法优点缺点适用场景有效集法迭代次数可预测需要可行初始点中小规模问题内点法理论复杂度低需要精细参数调节大规模问题梯度投影内存需求低收敛速度慢超大规模问题CVXOPT接口简单依赖第三方库快速原型开发实际项目中当约束数量超过1000时建议考虑内点法实现。但对于大多数工程问题有效集法在稳定性和实现复杂度上有着良好平衡。6. 常见问题与解决方案Q1算法陷入无限循环怎么办检查约束线性独立性确保没有冗余约束。可添加迭代次数限制和收敛判断if np.linalg.norm(x_new - x) 1e-6: breakQ2如何处理数值不稳定问题对Q矩阵进行正则化添加小量对角元素使用更稳定的线性求解器如QR分解适当缩放问题参数Q3如何加速算法收敛采用更智能的初始点选择策略在有效集变化时保留部分矩阵分解结果实现warm start机制在量化投资的实际应用中我们发现对协方差矩阵进行适当的收缩估计Ledoit-Wolf方法可以显著提升优化结果的稳定性。另一个实用技巧是对权重施加小幅L2正则化避免极端配置。