17.3 二次规划求解器 (QP Solver)
相关章节
17.3.1 理论基础
二次规划问题定义
二次规划(Quadratic Programming, QP)是一类特殊的数学优化问题，其目标函数是二次的，约束条件是线性的：
\min_{x} \quad & q(x) = x^T D x + c^T x + c_0 \\ \text{s.t.} \quad & Ax \leq b \\ & l \leq x \leq u \end{aligned}$$ 其中： - $x \in \mathbb{R}^n$ 是决策变量 - $D \in \mathbb{R}^{n \times n}$ 是Hessian矩阵（对称，通常要求半正定） - $c \in \mathbb{R}^n$ 是线性项系数 - $c_0 \in \mathbb{R}$ 是常数项 - $A \in \mathbb{R}^{m \times n}$ 是约束矩阵 - $b \in \mathbb{R}^m$ 是约束右端项 - $l, u \in \mathbb{R}^n$ 是变量上下界 ### 线性规划作为特例 当 $D = 0$ 时，QP退化为线性规划(LP)： $$\min_{x} \{ c^T x + c_0 \mid Ax \leq b, l \leq x \leq u \}$$ ### 凸性与最优性 **凸QP**：当 $D$ 是半正定矩阵时，QP是凸优化问题，具有以下性质： - 任何局部最优解都是全局最优解 - KKT条件是充分必要条件 - 存在多项式时间算法 **非凸QP**：当 $D$ 不定时，问题可能是NP难的。 ### 单纯形法原理 CGAL的QP求解器基于**改进的单纯形法**（Revised Simplex Method），核心思想： 1. **基变量与非基变量**：将变量分为基变量 $x_B$（在边界内）和非基变量 $x_N$（在边界上） 2. **基本可行解**：由 $m$ 个线性无关的约束确定一个顶点解 3. **转轴操作**：沿可行方向移动，从一个顶点转移到相邻的更好顶点 4. **最优性条件**：当不存在改进方向时达到最优 ### KKT条件 对于凸QP，KKT条件给出最优解的充要条件： $$\begin{aligned} 2Dx + c + A^T \lambda + \mu_l - \mu_u &= 0 \\ Ax - s &= b \\ x - l &\geq 0, \quad u - x \geq 0 \\ s &\geq 0 \\ \lambda^T s &= 0, \quad \mu_l^T(x-l) = 0, \quad \mu_u^T(u-x) = 0 \end{aligned}$$ 其中 $\lambda$ 是不等式约束的乘子，$\mu_l, \mu_u$ 是边界约束的乘子，$s$ 是松弛变量。 ## 17.3.2 架构分析 ### 模块结构 ``` QP_solver/ ├── include/CGAL/ │ ├── QP_models.h # 问题建模 │ ├── QP_functions.h # 求解函数 │ ├── QP_solution.h # 解的表示 │ ├── QP_options.h # 求解选项 │ └── QP_solver/ │ ├── QP_solver.h # 主求解器类 │ ├── QP_solver_impl.h # 实现细节 │ ├── QP_basis_inverse.h # 基逆矩阵 │ ├── QP_pricing_strategy.h # 定价策略基类 │ ├── QP_full_exact_pricing.h # 全精确定价 │ ├── QP_partial_exact_pricing.h # 部分精确定价 │ ├── QP_full_filtered_pricing.h # 全过滤定价 │ ├── QP_partial_filtered_pricing.h │ ├── QP_exact_bland_pricing.h # Bland规则 │ ├── Unbounded_direction.h # 无界方向 │ └── basic.h, functors.h # 基础设施 └── examples/ ├── first_qp.cpp # 基础示例 ├── first_lp.cpp # 线性规划 ├── first_nonnegative_qp.cpp # 非负QP └── investment.cpp # 投资组合优化 ``` ### 核心类层次 ``` Quadratic_program<NT> # 稀疏QP模型 ├── Quadratic_program_from_mps<NT> # MPS格式导入 ├── Quadratic_program_from_iterators # 迭代器构造 ├── Linear_program_from_iterators # LP特化 ├── Nonnegative_quadratic_program_from_iterators └── Nonnegative_linear_program_from_iterators QP_solver<Q, ET, Tags> # 主求解器 ├── QP_basis_inverse<ET, Is_linear> # 基逆管理 └── QP_pricing_strategy<Q, ET, Tags> # 定价策略 Quadratic_program_solution<ET> # 解对象 ``` ### 模板参数设计 ```cpp // Q: 问题类型（模型） // ET: 精确数类型（如Gmpz, MP_Float） // Tags: 编译时标签 // - Is_linear: Tag_true(LP) / Tag_false(QP) // - Is_nonnegative: Tag_true(标准形) / Tag_false(一般形) template < typename Q, typename ET, typename Tags > class QP_solver; ``` ### 问题建模接口 **稀疏模型**（推荐用于大多数应用）： ```cpp template <typename NT_> class Quadratic_program { public: typedef NT_ NT; // 设置约束矩阵 void set_a(int variable, int constraint, const NT& val); void set_b(int constraint, const NT& val); void set_r(int constraint, CGAL::Comparison_result val); // <=, =, >= // 设置边界 void set_l(int variable, bool is_finite, const NT& val = NT(0)); void set_u(int variable, bool is_finite, const NT& val = NT(0)); // 设置目标函数 void set_c(int variable, const NT& val); // 线性项 void set_d(int i, int j, const NT& val); // 二次项（i >= j） void set_c0(const NT& val); // 常数项 }; ``` **迭代器模型**（用于已有数据的高效包装）： ```cpp template < typename A_it, // 约束矩阵迭代器（列优先） typename B_it, // 右端项迭代器 typename R_it, // 关系迭代器 typename FL_it, // 下界有限性迭代器 typename L_it, // 下界迭代器 typename FU_it, // 上界有限性迭代器 typename U_it, // 上界迭代器 typename D_it, // 二次矩阵迭代器（行优先） typename C_it // 目标函数迭代器 > class Quadratic_program_from_iterators; ``` ## 17.3.3 实现细节 ### 求解器状态机 ```cpp // 来自 QP_solver.h enum Quadratic_program_status { QP_OPTIMAL, // 找到最优解 QP_INFEASIBLE, // 问题不可行 QP_UNBOUNDED, // 问题无界 QP_ERROR // 数值错误 }; class QP_solver { int m_phase; // 1: Phase I, 2: Phase II, 3: 完成 Quadratic_program_status m_status; int m_pivots; // 转轴次数 // 阶段I：寻找初始可行解 // 阶段II：优化目标函数 }; ``` ### 基逆矩阵管理 ```cpp // 来自 QP_basis_inverse.h template <class ET, class Is_linear> class QP_basis_inverse { // 存储基逆矩阵 M_B^{-1} // 对于LP：仅存储约束部分 // 对于QP：存储扩展矩阵 [A_C, 2D_{B_O,B_O}] void update(int entering, int leaving, ...); // 秩一更新 void transition(...); // 从Phase I到Phase II的转换 }; ``` ### 定价策略 定价策略决定选择哪个非基变量进入基： ```cpp // 来自 QP_pricing_strategy.h template <typename Q, typename ET, typename Tags> class QP_pricing_strategy { public: virtual int pricing(int& direction, ET& mu, ET& nu) = 0; virtual void leaving_basis(int i) = 0; }; // 实现类： // - QP_full_exact_pricing: 检查所有非基变量（精确） // - QP_partial_exact_pricing: 检查部分变量（更快） // - QP_full_filtered_pricing: 使用过滤的近似检查 // - QP_exact_bland_pricing: Bland规则（防循环） ``` ### 核心求解循环 ```cpp // 来自 QP_solver_impl.h Quadratic_program_status QP_solver::solve() { // Phase I: 寻找初始可行解 while (phase() == 1) { pivot_step(); } if (status() == QP_INFEASIBLE) { return QP_INFEASIBLE; } // 转换到Phase II transition(); // Phase II: 优化 while (phase() == 2) { pivot_step(); } return status(); } void QP_solver::pivot_step() { // 1. 定价：选择进入变量 pricing(); // 2. 比率测试：选择离开变量 ratio_test_init(); ratio_test_1(); if (needs_ratio_test_2()) { ratio_test_2(); } // 3. 更新：执行转轴 update_1(); update_2(); } ``` ### 比率测试 比率测试确定沿进入变量方向的步长： ```cpp // 来自 QP_solver_impl.h void QP_solver::ratio_test_1() { // 计算 q_x = M_B^{-1} * A_{C,j} // 确定最大可行步长 // 对于每个基变量，检查边界约束 for (each basic variable i) { if (q_i > 0) { // 检查下界 t_max = min(t_max, (x_i - l_i) / q_i); } else if (q_i < 0) { // 检查上界 t_max = min(t_max, (u_i - x_i) / (-q_i)); } } } ``` ### 解的提取 ```cpp // 来自 QP_solution.h template <class ET> class Quadratic_program_solution { public: // 解状态 bool is_optimal() const; bool is_infeasible() const; bool is_unbounded() const; // 目标函数值 ET objective_value() const; // 变量值（迭代器访问） Variable_value_iterator variable_values_begin() const; Variable_value_iterator variable_values_end() const; // 解的验证 template <typename QP> bool solves_quadratic_program(const QP& qp) const; }; ``` ## 17.3.4 使用示例 ### 基础QP求解 ```cpp #include <CGAL/QP_models.h> #include <CGAL/QP_functions.h> #include <CGAL/Gmpz.h> // 精确整数类型 #include <iostream> #include <assert> typedef CGAL::Gmpz ET; // 精确数类型 typedef CGAL::Quadratic_program<int> Program; typedef CGAL::Quadratic_program_solution<ET> Solution; int main() { // 创建非负QP模型，默认约束为 Ax <= b Program qp(CGAL::SMALLER, // 关系：<= true, // 下界有限 0, // 下界值 false, // 上界有限 0); // 上界值（未使用） // 变量索引 const int X = 0; const int Y = 1; // 设置约束 // 约束0: x + y <= 7 qp.set_a(X, 0, 1); qp.set_a(Y, 0, 1); qp.set_b(0, 7); // 约束1: -x + 2y <= 4 qp.set_a(X, 1, -1); qp.set_a(Y, 1, 2); qp.set_b(1, 4); // 设置上界 y <= 4 qp.set_u(Y, true, 4); // 设置二次目标函数: 2x^2 + 8y^2 - 32y + 64 // 注意：set_d设置的是2D矩阵的元素 qp.set_d(X, X, 2); // 2 * D_xx = 2 => D_xx = 1 qp.set_d(Y, Y, 8); // 2 * D_yy = 8 => D_yy = 4 qp.set_c(Y, -32); // 线性项 qp.set_c0(64); // 常数项 // 求解 Solution s = CGAL::solve_quadratic_program(qp, ET()); // 验证解 assert(s.solves_quadratic_program(qp)); // 输出结果 if (s.is_optimal()) { std::cout << "Optimal solution found\!" << std::endl; std::cout << "Objective value: " << s.objective_value() << std::endl; // 获取变量值 Solution::Variable_value_iterator it = s.variable_values_begin(); std::cout << "x = " << *it++ << std::endl; std::cout << "y = " << *it << std::endl; } return 0; } ``` ### 线性规划示例 ```cpp #include <CGAL/QP_models.h> #include <CGAL/QP_functions.h> #include <CGAL/Gmpz.h> typedef CGAL::Gmpz ET; int main() { // 创建LP（无二次项） CGAL::Quadratic_program<int> lp(CGAL::SMALLER, true, 0, false, 0); // 最大化 3x + 4y，转换为最小化 -3x - 4y lp.set_c(0, -3); // x的系数 lp.set_c(1, -4); // y的系数 // 约束 lp.set_a(0, 0, 1); // x lp.set_a(1, 0, 1); // y lp.set_b(0, 10); // x + y <= 10 lp.set_a(0, 1, 2); lp.set_a(1, 1, -1); lp.set_b(1, 4); // 2x - y <= 4 // 求解 CGAL::Quadratic_program_solution<ET> s = CGAL::solve_linear_program(lp, ET()); if (s.is_optimal()) { std::cout << "Optimal value: " << s.objective_value() << std::endl; } return 0; } ``` ### 从MPS文件导入 ```cpp #include <CGAL/QP_models.h> #include <CGAL/QP_functions.h> #include <fstream> typedef CGAL::Gmpz ET; int main() { // 从MPS文件读取问题 std::ifstream in("problem.mps"); CGAL::Quadratic_program_from_mps<int> qp(in); if (\!qp.is_valid()) { std::cerr << "Error reading MPS file: " << qp.get_error() << std::endl; return 1; } std::cout << "Problem name: " << qp.get_problem_name() << std::endl; std::cout << "Variables: " << qp.get_n() << std::endl; std::cout << "Constraints: " << qp.get_m() << std::endl; // 求解 CGAL::Quadratic_program_solution<ET> s = CGAL::solve_quadratic_program(qp, ET()); return 0; } ``` ### 迭代器接口（零拷贝） ```cpp #include <CGAL/QP_models.h> #include <CGAL/QP_functions.h> #include <vector> typedef CGAL::Gmpz ET; int main() { // 已有数据存储 std::vector<std::vector<int>> A = {{1, 1}, {-1, 2}}; std::vector<int> b = {7, 4}; std::vector<CGAL::Comparison_result> r = {CGAL::SMALLER, CGAL::SMALLER}; std::vector<bool> fl = {true, true}; std::vector<int> l = {0, 0}; std::vector<bool> fu = {false, true}; std::vector<int> u = {0, 4}; std::vector<std::vector<int>> D = {{2, 0}, {0, 8}}; std::vector<int> c = {0, -32}; // 使用迭代器构造（无数据拷贝） auto qp = CGAL::make_quadratic_program_from_iterators( 2, // 变量数 2, // 约束数 A.begin(), b.begin(), r.begin(), fl.begin(), l.begin(), fu.begin(), u.begin(), D.begin(), c.begin(), 64 // c0 ); CGAL::Quadratic_program_solution<ET> s = CGAL::solve_quadratic_program(qp, ET()); return 0; } ``` ### 投资组合优化 ```cpp #include <CGAL/QP_models.h> #include <CGAL/QP_functions.h> #include <vector> #include <math> typedef CGAL::Gmpz ET; // 马科维茨投资组合优化 // 最小化风险（方差）同时达到目标收益 int main() { const int n = 3; // 资产数量 // 预期收益率 std::vector<double> expected_return = {0.1, 0.15, 0.2}; // 协方差矩阵（对称） std::vector<std::vector<double>> cov = { {0.04, 0.02, 0.01}, {0.02, 0.09, 0.03}, {0.01, 0.03, 0.16} }; double target_return = 0.12; // 创建QP：变量为投资比例 CGAL::Quadratic_program<double> qp(CGAL::EQUAL, true, 0, false, 0); // 设置二次项（风险）：2 * D = 协方差矩阵 for (int i = 0; i < n; ++i) { for (int j = 0; j <= i; ++j) { qp.set_d(i, j, 2 * cov[i][j]); // 注意：set_d设置2D的元素 } } // 约束1：投资比例之和为1 for (int i = 0; i < n; ++i) { qp.set_a(i, 0, 1); } qp.set_b(0, 1); qp.set_r(0, CGAL::EQUAL); // 约束2：预期收益等于目标 for (int i = 0; i < n; ++i) { qp.set_a(i, 1, expected_return[i]); } qp.set_b(1, target_return); qp.set_r(1, CGAL::EQUAL); // 求解 CGAL::Quadratic_program_solution<ET> s = CGAL::solve_quadratic_program(qp, ET()); if (s.is_optimal()) { std::cout << "Minimum risk portfolio:" << std::endl; auto it = s.variable_values_begin(); for (int i = 0; i < n; ++i) { std::cout << "Asset " << i << ": " << CGAL::to_double(*it++) << std::endl; } std::cout << "Risk (variance): " << CGAL::to_double(s.objective_value()) << std::endl; } return 0; } ``` ## 17.3.5 复杂度分析 ### 时间复杂度 | 操作 | 最坏情况 | 平均情况 | 说明 | |------|----------|----------|------| | Phase I | $O(m^3 \cdot n)$ | $O(m^2 \cdot n)$ | 寻找可行解 | | Phase II | $O(2^n \cdot m^3)$ | $O(n \cdot m^3)$ | 单纯形迭代 | | 每次转轴 | $O(m^2)$ | $O(m)$ | 基逆更新 | | 定价 | $O(n \cdot m)$ | $O(n)$ | 取决于策略 | 其中： - $n$：变量数 - $m$：约束数 ### 空间复杂度 - 问题数据：$O(n^2 + nm)$（二次矩阵 + 约束矩阵） - 基逆矩阵：$O(m^2)$ - 工作向量：$O(n + m)$ - **总体**：$O(n^2 + nm)$ ### 数值稳定性 1. **精确算术**：使用ET（如Gmpz）进行关键计算，避免舍入误差 2. **过滤定价**：先用浮点快速筛选，再用精确计算验证 3. **Bland规则**：防止循环，保证终止性 4. **缩放**：内部自动处理数值缩放 ### 性能优化建议 1. **选择合适的定价策略**： - 小规模问题：全精确定价 - 大规模问题：部分过滤定价 2. **利用问题结构**： - 非负变量使用`Nonnegative_quadratic_program` - 线性问题使用`solve_linear_program` 3. **稀疏性**：使用`Quadratic_program`稀疏接口 ## 17.3.6 与其他优化库的比较 | 特性 | CGAL QP | GLPK | CPLEX | Gurobi | |------|---------|------|-------|--------| | 精确算术 | 是 | 否 | 否 | 否 | | 开源 | 是 | 是 | 否 | 否 | | 浮点过滤 | 是 | 否 | 是 | 是 | | 问题规模 | 中小 | 大 | 很大 | 很大 | | 几何集成 | 优秀 | 无 | 无 | 无 | CGAL QP求解器的优势在于： - **精确性**：保证结果的正确性 - **几何集成**：与CGAL其他模块无缝协作 - **模板设计**：灵活支持各种数类型 ## 17.3.7 几何应用 ### 最近点问题 ```cpp // 在凸多面体上寻找距离给定点最近的点 // 可转化为QP问题 ``` ### 最小包围球 ```cpp // 通过QP对偶问题求解最小包围球 // 与CGAL::Min_sphere_d结果一致 ``` ### 碰撞检测 ```cpp // 分离轴定理可转化为QP可行性问题 ``` ### 网格变形能量最小化 ```cpp // 基于物理的变形可建模为QP // 能量 = 二次形变能 + 线性约束 ```