企业官网建设流程全解析

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

线性规划作为运筹学中的经典工具，在资源分配、生产计划、投资组合等场景中有着广泛应用。虽然市面上有成熟的求解器如GLPK、CPLEX，但理解算法底层实现对于开发者而言，不仅能提升问题建模能力，还能在特殊需求时进行定制优化。本文将带你用C++从零实现两阶段单纯形算法，通过面向对象设计封装完整求解过程。

1. 线性规划与单纯形算法基础

线性规划问题的标准形式可表示为：

最大化 cᵀx 约束条件 Ax ≤ b x ≥ 0

其中x为决策变量向量，c为目标函数系数，A为约束矩阵，b为资源限制向量。

单纯形算法的核心思想是通过在可行域的顶点间移动，逐步逼近最优解。其数学本质是：

• 将不等式约束通过松弛变量转化为等式
• 构建单纯形表表示当前基本可行解
• 通过转轴运算在基变量与非基变量间交换

两阶段法的特殊之处在于：

1. 第一阶段：引入人工变量构造辅助问题，找到初始可行解
2. 第二阶段：用找到的可行解求解原问题

// 典型问题输入格式示例 1 // 1:max, -1:min 4 4 // 约束数m=4, 变量数n=4 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 // x1+x2+x3+x4 = 9 0 1 -1 2 1 // x2-x3+2x4 ≥1 1 1 3 -1 // 目标函数 z = x1+x2+3x3-x4

2. 类架构设计与核心成员

我们设计的LinearProgram类需要处理三大核心任务：

• 问题输入与初始化
• 两阶段求解过程
• 结果输出与验证

2.1 数据成员规划

class LinearProgram { private: int m, n; // 约束总数和变量数 int m1, m2, m3; // ≤/=≥约束计数 int n1, n2; // 扩展变量后的维度 double **a; // 单纯形表二维数组 int *basic; // 基变量下标 int *nonbasic; // 非基变量下标 double minmax; // 1=max, -1=min int error; // 错误代码 };

关键设计考量：

• 动态二维数组：使用Make2DArray/Delet2DArray管理内存
• 灵活的变量标记：basic/nonbasic数组跟踪变量状态
• 错误处理机制：通过error代码标识无解、无界等情况

2.2 核心方法分解

3. 算法实现关键细节

3.1 入基与离基选择

入基变量选择遵循Bland法则避免循环：

int LinearProgram::enter(int objrow) { double temp = DBL_EPSILON; int col = 0; for(int j=1; j<=n1; 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 i=1; i<=m; 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 转轴运算实现

转轴运算包含三个关键步骤：

1. 归一化主元行
2. 消去其他行的主元列
3. 更新基变量记录

void LinearProgram::pivot(int row, int col) { // 1. 归一化主元行 for(int j=0; j<=n1; j++) { if(j != col) a[row][j] /= a[row][col]; } a[row][col] = 1.0/a[row][col]; // 2. 消去其他行 for(int i=0; i<=m+1; i++) { if(i != row) { for(int j=0; j<=n1; 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(m+1); // 使用辅助目标函数 if(error > 0) return error; // 检查人工变量是否全部消去 for(int i=1; i<=m; i++) { if(basic[i] > n2) { if(a[i][0] > DBL_EPSILON) return 3; // 无可行解 for(int j=1; j<=n1; 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 i=row+1; i<=m; 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 i=0; i<m; i++) { a[i] = new double[n](); // 值初始化 } } void LinearProgram::Delet2DArray(double** &a, int m, int n) { for(int i=0; i<m; 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 =3x1+5x2

输出结果示例：

最优值：34.0000 最优解： x1 = 3.0000 x2 = 5.0000

7. 进阶扩展方向

对于希望进一步优化的开发者，可以考虑：

1. 稀疏矩阵存储：对大规模问题使用CSR/CSC格式
2. 并行化处理：将矩阵运算转为BLAS实现
3. 预处理技术：检测冗余约束、固定变量等
4. 对偶单纯形法：处理约束条件变化的情况

// 稀疏矩阵存储示例 struct SparseMatrix { vector<double> values; vector<int> col_indices; vector<int> row_ptr; };

实现过程中最易出错的是转轴运算的索引处理，建议使用断言检查数组边界：

void LinearProgram::pivot(int row, int col) { assert(row >=1 && row <=m); assert(col >=1 && col <=n1); // ...其余代码 }