从零实现两阶段单纯形算法：一个C++类搞定线性规划求解（附完整代码）

从零实现两阶段单纯形算法一个C类搞定线性规划求解附完整代码线性规划作为运筹学中的经典工具在资源分配、生产计划、投资组合等场景中有着广泛应用。虽然市面上有成熟的求解器如GLPK、CPLEX但理解算法底层实现对于开发者而言不仅能提升问题建模能力还能在特殊需求时进行定制优化。本文将带你用C从零实现两阶段单纯形算法通过面向对象设计封装完整求解过程。1. 线性规划与单纯形算法基础线性规划问题的标准形式可表示为最大化 cᵀx 约束条件 Ax ≤ b x ≥ 0其中x为决策变量向量c为目标函数系数A为约束矩阵b为资源限制向量。单纯形算法的核心思想是通过在可行域的顶点间移动逐步逼近最优解。其数学本质是将不等式约束通过松弛变量转化为等式构建单纯形表表示当前基本可行解通过转轴运算在基变量与非基变量间交换两阶段法的特殊之处在于第一阶段引入人工变量构造辅助问题找到初始可行解第二阶段用找到的可行解求解原问题// 典型问题输入格式示例 1 // 1:max, -1:min 4 4 // 约束数m4, 变量数n4 2 1 1 // ≤约束2个, 约束1个, ≥约束1个 1 0 2 0 18 // x1 2x3 ≤ 18 0 2 0 -7 0 // 2x2 -7x4 ≤ 0 1 1 1 1 9 // x1x2x3x4 9 0 1 -1 2 1 // x2-x32x4 ≥1 1 1 3 -1 // 目标函数 z x1x23x3-x42. 类架构设计与核心成员我们设计的LinearProgram类需要处理三大核心任务问题输入与初始化两阶段求解过程结果输出与验证2.1 数据成员规划class LinearProgram { private: int m, n; // 约束总数和变量数 int m1, m2, m3; // ≤/≥约束计数 int n1, n2; // 扩展变量后的维度 double **a; // 单纯形表二维数组 int *basic; // 基变量下标 int *nonbasic; // 非基变量下标 double minmax; // 1max, -1min int error; // 错误代码 };关键设计考量动态二维数组使用Make2DArray/Delet2DArray管理内存灵活的变量标记basic/nonbasic数组跟踪变量状态错误处理机制通过error代码标识无解、无界等情况2.2 核心方法分解方法名功能描述数学对应enter()选择入基变量最大正检验数规则leave()选择离基变量最小比值测试pivot()执行转轴运算高斯-约当消元phase1()第一阶段求解人工变量消去phase2()第二阶段求解原始目标函数优化3. 算法实现关键细节3.1 入基与离基选择入基变量选择遵循Bland法则避免循环int LinearProgram::enter(int objrow) { double temp DBL_EPSILON; int col 0; for(int j1; jn1; j) { if(nonbasic[j]n2 a[objrow][j]temp) { col j; temp a[objrow][j]; break; // Bland规则 } } return col; }离基变量通过最小比值测试确定int LinearProgram::leave(int col) { double temp DBL_MAX; int row 0; for(int i1; im; i) { double val a[i][col]; if(val DBL_EPSILON) { val a[i][0]/val; if(val temp) { row i; temp val; } } } return row; }3.2 转轴运算实现转轴运算包含三个关键步骤归一化主元行消去其他行的主元列更新基变量记录void LinearProgram::pivot(int row, int col) { // 1. 归一化主元行 for(int j0; jn1; j) { if(j ! col) a[row][j] / a[row][col]; } a[row][col] 1.0/a[row][col]; // 2. 消去其他行 for(int i0; im1; i) { if(i ! row) { for(int j0; jn1; j) { if(j ! col) { a[i][j] - a[i][col]*a[row][j]; if(fabs(a[i][j]) DBL_EPSILON) a[i][j] 0.0; } } a[i][col] -a[i][col]*a[row][col]; } } swapbasic(row, col); }注意浮点数比较使用DBL_EPSILON阈值避免精度误差导致错误判断4. 两阶段法完整流程4.1 第一阶段构造辅助问题int LinearProgram::phase1() { error simplex(m1); // 使用辅助目标函数 if(error 0) return error; // 检查人工变量是否全部消去 for(int i1; im; i) { if(basic[i] n2) { if(a[i][0] DBL_EPSILON) return 3; // 无可行解 for(int j1; jn1; j) { if(fabs(a[i][j]) DBL_EPSILON) { pivot(i,j); break; } } } } return 0; }4.2 第二阶段求解原问题int LinearProgram::phase2() { return simplex(0); // 使用原始目标函数 } int LinearProgram::compute() { if(error 0) return error; if(m ! m1) { // 存在非≤约束时需要第一阶段 error phase1(); if(error 0) return error; } return phase2(); }5. 工程实践中的优化技巧在实际编码中我们还需要处理以下关键问题5.1 数值稳定性处理主元缩放避免极小主元导致数值不稳定定期矩阵清理将接近零的元素显式设为0迭代次数限制防止病态问题无限循环// 改进后的pivot操作片段 if(fabs(a[row][col]) 1e-10) { // 寻找同列其他非零主元 for(int irow1; im; i) { if(fabs(a[i][col]) 1e-8) { row_swap(row, i); break; } } }5.2 内存管理策略使用RAII模式管理动态数组void LinearProgram::Make2DArray(double** a, int m, int n) { a new double*[m]; for(int i0; im; i) { a[i] new double[n](); // 值初始化 } } void LinearProgram::Delet2DArray(double** a, int m, int n) { for(int i0; im; i) { delete[] a[i]; } delete[] a; a nullptr; }6. 完整代码实现与测试最终的类实现包含以下关键文件结构LinearProgram/ ├── include/ │ └── LinearProgram.h // 类声明 ├── src/ │ └── LinearProgram.cpp // 类实现 └── test/ └── main.cpp // 测试用例典型测试案例// 生产计划问题 1 // 最大化 3 2 // 3约束2变量 2 0 1 // 2个≤, 1个≥ 2 1 12 // 2x1 x2 ≤12 1 1 8 // x1 x2 ≤8 1 0 3 // x1 ≥3 3 5 // z 3x15x2输出结果示例最优值34.0000 最优解 x1 3.0000 x2 5.00007. 进阶扩展方向对于希望进一步优化的开发者可以考虑稀疏矩阵存储对大规模问题使用CSR/CSC格式并行化处理将矩阵运算转为BLAS实现预处理技术检测冗余约束、固定变量等对偶单纯形法处理约束条件变化的情况// 稀疏矩阵存储示例 struct SparseMatrix { vectordouble values; vectorint col_indices; vectorint row_ptr; };实现过程中最易出错的是转轴运算的索引处理建议使用断言检查数组边界void LinearProgram::pivot(int row, int col) { assert(row 1 row m); assert(col 1 col n1); // ...其余代码 }