别再死记硬背单纯形法了！用C++手写一个两阶段求解器，从原理到代码一次搞定
从零实现两阶段单纯形法：C++代码与数学原理深度解析
在优化问题的求解中，线性规划占据着核心地位。而单纯形法作为解决线性规划问题的经典算法，其重要性不言而喻。但很多学习者在掌握理论后，往往面临一个困境：如何将这些数学步骤转化为实际可运行的代码？本文将彻底解决这个问题，带你从数学原理出发，用C++完整实现一个两阶段单纯形法求解器。
1. 单纯形法的数学基础与实现框架
单纯形法的核心思想是在可行域的顶点间移动，逐步逼近最优解。要实现这一算法，我们需要先理解几个关键概念：
- 基本变量(Basic variables)：当前解中非零的变量
- 非基本变量(Non-basic variables)：当前解中为零的变量
- 单纯形表(Simplex tableau)：包含所有约束条件和目标函数的矩阵表示
两阶段法的特殊之处在于：
- 第一阶段通过引入人工变量寻找初始可行解
- 第二阶段在找到的可行解基础上优化原始目标函数
下面是我们C++实现的核心类结构：
class LinearProgram { public: LinearProgram(); ~LinearProgram(); void solve(); private: // 核心算法方法 int enter(int objrow); // 选择入基变量 int leave(int col); // 选择离基变量 void pivot(int row, int col); // 转轴变换 int phase1(); // 第一阶段求解 int phase2(); // 第二阶段求解 // 辅助方法 void swapbasic(int row, int col); void Make2DArray(double** &a, int m, int n); // 数据成员 int m, n; // 约束数和变量数 double** a; // 单纯形表 int* basic, *nonbasic; // 基本和非基本变量索引 };2. 数据结构设计与初始化
实现单纯形法的第一步是设计合适的数据结构来存储单纯形表和相关变量。我们使用二维数组a来表示整个单纯形表，其中：
- 第0行存储目标函数系数
- 第1到m行存储约束条件
- 第m+1行用于第一阶段的人工目标函数
初始化过程需要处理三种约束类型：
- ≤约束：添加松弛变量
- =约束：直接保留
- ≥约束：添加剩余变量和人工变量
LinearProgram::LinearProgram() { // 读取输入参数 cin >> minmax; // 1为max，-1为min cin >> m >> n; // 约束数和变量数 cin >> m1 >> m2 >> m3; // ≤, =, ≥约束的数量 // 分配内存 Make2DArray(a, m+2, n1+1); basic = new int[m+2]; nonbasic = new int[n1+1]; // 初始化单纯形表 for(int i=0; i<=m+1; i++) for(int j=0; j<=n1; j++) a[i][j] = 0.0; // 处理≥约束，添加人工变量 for(int i=m-m3+1, j=n+1; i<=m; i++, j++) { a[i][j] = -1.0; a[m+1][j] = -1.0; // 人工变量在辅助目标函数中的系数 } }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; }3.2 离基变量选择
离基变量的选择需要保证解始终可行，我们使用最小比值测试：
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.3 转轴变换
转轴变换是单纯形法的核心操作，它实现了基变量的替换：
void LinearProgram::pivot(int row, int col) { // 标准化主元行 for(int j=0; j<=n1; j++) { if(j != col) a[row][j] /= a[row][col]; } a[row][col] = 1.0 / a[row][col]; // 更新其他行 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); // 更新基变量索引 }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. 数值稳定性与工程实践
在实际实现中，我们需要特别注意数值稳定性问题：
- 浮点数比较：使用
DBL_EPSILON作为零的阈值 - Bland法则：防止算法陷入循环
- 内存管理：正确释放动态分配的数组
LinearProgram::~LinearProgram() { delete[] basic; delete[] nonbasic; Delet2DArray(a, m+2, n1+1); } void LinearProgram::Delet2DArray(double** &a, int m, int n) { for(int i=0; i<=m; i++) delete[] a[i]; delete[] a; }6. 使用示例与结果输出
完整的求解流程封装在solve()方法中：
void LinearProgram::solve() { error = compute(); switch(error) { case 0: output(); break; case 1: cout << "输入数据错误" << endl; break; case 2: cout << "无界解" << endl; break; case 3: cout << "无可行解" << endl; } } void LinearProgram::output() { cout << "最优值：" << -minmax * a[0][0] << endl; cout << "最优解：" << endl; for(int j=1; j<=n; j++) { cout << "x" << j << " = "; bool found = false; for(int i=1; i<=m; i++) { if(basic[i] == j) { cout << a[i][0]; found = true; break; } } if(!found) cout << "0.0"; cout << endl; } }这个实现完整展示了两阶段单纯形法的所有关键环节，从数学原理到实际代码的转化过程。通过这个案例，我们不仅理解了算法的运作机制，还掌握了如何将数学算法转化为健壮的软件实现。