C++实现单纯形法求解线性规划：完整指南与工程实践

简介：C++实现的单纯形算法计算程序，是一份面向运筹学课程设计、算法爱好者及编程初学者的可运行源码包。它用C++语言实现了线性规划领域经典的单纯形法，围绕标准型转换、单纯形表构建、基变量与非基变量迭代替换等核心步骤展开，程序结构清晰，能够帮助使用者快速理解线性规划求解流程。

包体紧凑，共3个文件，包含2个C++源文件和1个头文件，分别承担矩阵运算、核心算法与主程序调用等模块，压缩包仅2KB，便于阅读和二次修改。目前已有741人学习/下载，适合正在学习运筹学、数值计算或准备相关课程设计的开发者参考。

通过这份资源，读者可以掌握单纯形法从建模到编码的完整落地思路，学习如何借助STL容器存储矩阵、处理约束条件并实现迭代收敛判断，同时了解基本的异常处理和调试方法，为后续扩展更复杂的优化算法打下基础。 最近接了个小项目，需要在不依赖第三方优化库的前提下实现线性规划求解。犹豫了半天，最后还是决定用C++把单纯形算法完整实现一遍。前前后后花了差不多两天，踩了不少数值坑，也把两阶段法、退化处理这些细节都补全了。这个程序能从标准型线性规划出发，自动求出最优解，也能正确判断无界解和无可行解，实测下来和MATLAB的linprog结果完全一致。

这篇文章就把整个实现过程拆开讲清楚。不管你是正在学运筹学、准备算法面试，还是工作中突然要处理线性规划问题，这篇都能给你一个可以抄作业的参考版本。

1. 项目概述：一个能落地的单纯形法求解器

1.1 这个程序到底解决什么问题

单纯形算法解决的是线性规划问题，就是在一堆线性不等式约束下，求一个线性目标函数的最大值或最小值。典型的场景包括生产计划（有限原料下怎么排产利润最高）、运输调度（多个产地和销地之间怎么调运成本最低）、资源分配（预算有限时怎么投放效果最好）。

这个程序做的事情很直接：你输入约束矩阵、右侧常数向量和目标函数系数，它输出最优解和目标函数值。如果问题无界或者无可行解，它也会明确告诉你，而不是卡死或者给个莫名其妙的结果。

实现过程中最麻烦的其实不是算法主流程，而是各种边界情况：没有初始可行基怎么办、遇到退化循环怎么处理、浮点误差怎么控制。这些我在后面都会逐个展开。

1.2 为什么选C++而不是MATLAB或Python

这个问题我纠结过。Python有scipy.optimize.linprog，MATLAB有linprog，都是现成的，几行代码就出结果。但这次的需求有两个特点：一是要嵌入到已有的C++系统中，不能引入Python运行时；二是数据量不算小，迭代过程中要做大量的矩阵行变换，C++的性能优势很明显。

从学习角度看，自己实现一遍单纯形法，对算法的理解深度完全不是调用API能比的。面试时被问到“单纯形法的入基出基规则”“怎么处理退化”这类问题时，亲手写过的人明显能答到点子上。这就是那些C++八股文里反复强调的底层理解。

当然我也要说实话：如果只是做一次性数据分析，直接用Python完全没问题。这个C++版本适合的场景是学习、教学、嵌入生产系统，以及需要完全掌控中间过程的场景。

2. 算法原理与单纯形表设计

2.1 线性规划标准型与松弛变量

单纯形法要求问题必须先转化成标准形式：

最小化 z = c^T x，满足 Ax = b（b >= 0），且 x >= 0。

也就是说，所有不等式约束都要变成等式约束，所有变量都要是非负的。这里就需要引入松弛变量和剩余变量。

大于等于号约束，左边减掉一个非负剩余变量变成等式；小于等于号约束，左边加上一个非负松弛变量变成等式。比如 2x1 + x2 <= 4，变成 2x1 + x2 + s1 = 4，其中s1就是松弛变量，它在目标函数里的系数是0。

这个转化是后面所有计算的基础，代码里必须处理得干净。我实现时单独写了一个预处理函数，把用户输入的任意不等式约束统一转换成标准型，同时记录哪些列是人工添加的变量列，方便后面两阶段法使用。

2.2 为什么用“单纯形表”而不是矩阵形式

单纯形法有两种经典的实现路径。一种是基于矩阵分解的 revised simplex method，每次迭代只更新基矩阵的逆；另一种就是这里用的 tableau 形式，把所有系数、右端项、检验数放在一张表里做行变换。

tableau 形式的核心就是维护一张增广矩阵：

最后一行是检验数行（reduced cost）。对于最小化问题，只要还有检验数为负，就说明目标函数还能继续下降，需要继续迭代。

我选tableau形式的原因很实际：代码简单直观，调试时可以直接把表格打出来看中间状态，教学时也容易讲清楚。矩阵形式虽然在大规模问题上迭代更快，但实现复杂度高不少，对于这个项目来说属于过度设计。只有当问题规模上千行上万列时，才值得考虑revised simplex和稀疏矩阵存储。

3. 核心实现：C++代码怎么写才不翻车

3.1 类结构设计与存储方案

代码结构上我封装了一个 Simplex 类，内部核心成员就四个：单纯形表、基变量索引、求解状态、一个静态的EPS常量。

class Simplex { public: enum class Status { OPTIMAL, UNBOUNDED, INFEASIBLE, RUNNING }; Simplex(const std::vector<std::vector<double>>& A, const std::vector<double>& b, const std::vector<double>& c, bool maximize); bool solve(); std::vector<double> solution() const; double objectiveValue() const; Status status() const { return status_; } private: std::vector<std::vector<double>> tableau_; std::vector<int> baseIndex_; // 基变量在原始列中的索引 Status status_ = Status::RUNNING; static constexpr double EPS = 1e-9; };

存储直接用 vector<vector > 嵌套二维数组。有人可能会提动态二维数组或者 unique_ptr，但在这个场景下完全没必要。vector 自动管理内存，连续访问快，而且写起来不容易出内存错误。我之前见过有人为了“性能”手写裸指针二维数组，结果析构函数写错了内存泄漏，反而是最不划算的选择。记住：程序正确性永远是第一位的，性能优化在明确瓶颈之后再做。

3.2 入基、出基与pivot操作

单纯形迭代的核心就是三步：找入基变量、找出基变量、做pivot行变换。

入基变量的选择，最小化问题里通常选检验数最小（负得最多）的那一列，这叫Dantzig规则，收敛通常最快。

int Simplex::enteringColumn() const { int enter = -1; double minReducedCost = -EPS; for (int j = 0; j < static_cast<int>(tableau_[0].size()) - 1; ++j) { if (tableau_.back()[j] < minReducedCost) { minReducedCost = tableau_.back()[j]; enter = j; } } return enter; }

需要注意的是比较条件用了< -EPS而不是< 0。这是数值稳定性的大坑之一：浮点运算会让理论上应该为0的检验数变成 -1e-14 这种微小负数，如果不设阈值过滤，程序会在最优解附近反复做无意义的入基出基，甚至死循环。

出基变量的选择遵循最小比值原则。对每一行，如果入基列系数 a_ij > EPS，计算 RHS / a_ij，选比值最小的那一行的基变量出基。

int Simplex::leavingRow(int enter) const { int leave = -1; double minRatio = std::numeric_limits<double>::infinity(); for (int i = 0; i < static_cast<int>(tableau_.size()) - 1; ++i) { double a = tableau_[i][enter]; if (a > EPS) { double ratio = tableau_[i].back() / a; if (ratio < minRatio) { minRatio = ratio; leave = i; } } } return leave; }

如果找不到任何一行满足 a_ij > EPS，说明入基变量可以无限增大，目标函数无下界，直接判定问题无界。

pivot操作就是标准的高斯消元：先把主元行除以主元，让主元变成1，然后消去其他所有行（包括检验数行）中该列的系数。

void Simplex::pivot(int row, int col) { double pivotVal = tableau_[row][col]; for (auto& val : tableau_[row]) { val /= pivotVal; } for (int i = 0; i < static_cast<int>(tableau_.size()); ++i) { if (i == row) continue; double factor = tableau_[i][col]; if (std::abs(factor) < EPS) continue; for (int j = 0; j < static_cast<int>(tableau_[i].size()); ++j) { tableau_[i][j] -= factor * tableau_[row][j]; } } baseIndex_[row] = col; }

这一块坑很多。首先是累加误差，每做一次pivot所有元素都会更新一遍，几十轮迭代后初始输入的信息会混入大量浮点噪声。所以EPS阈值不是摆设，是保证程序能在有限步内终止的关键。其次，如果表格里出现了数值很小的负数（理论应为0但实际是 -1e-12），后续的入基判断就会出错，所以在pivot里对小于EPS的系数直接跳过是必要的。

solve() 主循环很简单：

bool Simplex::solve() { while (true) { int enter = enteringColumn(); if (enter == -1) { status_ = Status::OPTIMAL; return true; } int leave = leavingRow(enter); if (leave == -1) { status_ = Status::UNBOUNDED; return false; } pivot(leave, enter); } }

实际工程里我会加一个最大迭代次数保护（比如约束列数的100倍），防止因数值问题或者退化导致理论上的死循环在实践里真的发生。虽然理论上有Bland规则能保证终止，但程序卡死一次的代价远大于多加一个计数器。

3.3 两阶段法：没有初始可行基怎么办

标准型线性规划要求 Ax = b 且 b >= 0，同时初始基变量存在。但现实输入往往没有天然的单位矩阵作为初始基。比如约束里有大于等于号时，引入剩余变量后系数是 -1，没法直接当基变量。

解决办法是用两阶段法。第一阶段引入人工变量，让每个约束都有一个系数为1的人工变量，然后最小化所有人工变量之和。如果最优值大于EPS，说明原问题无可行解；如果最优值为0，人工变量全部变成0，就可以进入第二阶段，去掉人工变量列，用原目标函数继续迭代。

这里有个关键实现细节：第一阶段目标函数的检验数行不能直接塞进tableau里，因为人工变量的初始基变量对应的检验数必须为0，否则单纯形法初始解就不满足最优性条件。正确做法是先把人工变量的目标系数设为1，然后通过行变换把基变量（人工变量）在目标行里的系数消成0。

我一开始就是在这里翻的车——直接把人工变量的目标系数写了进去，结果第一阶段迭代方向完全错误，解出来的人工变量和不为0，却判断成有可行解。浪费了好几个小时调试。后来打印每一步的单纯形表才发现问题。

两阶段法和给约束加M倍人工变量惩罚的大M法相比，优势是避免了选择M值的困扰。M太小，可能惩罚不够导致人工变量残留；M太大，数值计算中会放大浮点误差。我建议都写两阶段法，更稳健。

4. 数值稳定性与退化循环那些坑

4.1 EPS不能随便设，浮点比较要命

这个项目里我踩的最深的就是浮点比较问题。C++里直接判断两个double是否相等就是灾难，单纯形表里的每一个数都是经过多轮行变换计算出来的，误差在几千万分之一级别非常正常。

我最终的EPS取值是1e-9。这个值不是拍脑袋定的，参考了数值线性代数教材的常见实践，也结合了我自己的测试：在最优解附近，检验数误差通常小于1e-10。EPS太小会在边界处失效，EPS太大会把真正需要入基的变量错误过滤掉，导致次优解。

具体到代码里有三个地方必须用EPS：

1. 入基判断：检验数 < -EPS 才算不合格，< 0 会让你在最优附近空转
2. 出基判断：系数 > EPS 才参与最小比值计算，否则会把 a=1e-12 的数值当成有效主元，产生巨大的比值
3. pivot消元：|factor| < EPS 直接跳过，减少无意义的浮点运算

我用一个生产计划算例验证过：如果不加EPS或者把EPS设成0，程序有时会多迭代两三轮，有时会直接卡死。加上EPS之后一切正常。数值计算里，用阈值隔离浮点噪声，这属于基本功。

4.2 退化循环与Bland规则

退化是指多个基变量同时取0的情况，此时最小比值计算结果为0，出基之后基变量集合可能没变，目标函数值也不变。严重时会出现循环，就是一系列入基出基操作后回到之前的某个基，永远无法到达最优解。

经典教材里的Beale例子就是为了说明这个问题而构造的，会让Dantzig规则直接陷入死循环。虽然实际业务问题里退化循环极少见，但既然做通用求解器，就必须处理。

理论上有Bland规则保证终止：入基时选择检验数为负的最小小标号变量，出基时在比值最小的候选者里选择最小标号。这个规则牺牲了收敛速度，但保证了不会循环。

我的做法是做一个运行开关，默认用Dantzig规则保证收敛速度，同时用一个迭代计数器。当迭代次数超过某个阈值（我设的是约束数量乘100），自动切换成Bland规则。这个策略既保证了大问题上的性能，又在极端情况下保底。如果你追求代码简单，也可以全部用Bland规则，对大部分小规模问题速度差异其实不大。

5. 测试样例与调试实录

5.1 教科书算例：生产计划问题

我用经典的例子做基准测试：

最大化 z = 3x1 + 2x2

约束： 2x1 + x2 <= 4 x1 + x2 <= 3 x1, x2 >= 0

这个例子的手工最优解是 x1=1, x2=2，目标值7。程序输出结果：

迭代1：进入变量 x1，离开变量 s1 迭代2：进入变量 x2，离开变量 s2 最优解：x1 = 1.000000, x2 = 2.000000 目标函数值：7.000000

和手算完全一致。这里有个验证技巧：每次迭代我都把当前单纯形表打印出来，手动检查检验数行和RHS列是否合理。调试数值程序的时候，这种可视化输出比GDB断点好用得多，能直接看到数据变化趋势。

我还用这个例子做了输入鲁棒性测试：把约束顺序打乱、把目标函数系数从int转成double、把右侧常数写成带小数位的浮点数，程序输出结果都一样。这说明初始化和pivot操作对输入顺序不敏感，是我要的效果。

5.2 无界解与无可行解的判定测试

无界解的测试我构造了一个简单例子：最大化 x1 + x2，约束 -x1 + x2 <= 1，x1、x2 >= 0。这个问题的可行域向x1正方向无限延伸，目标函数可以无限增大。程序正确输出了UNBOUNDED，并且能指出是无界解而不是数值故障。

无可行解的测试用这个例子：最小化 x1 + x2，约束 x1 + x2 <= 1 和 x1 + x2 >= 2，x1、x2 >=0。显然两个约束互相矛盾。第一阶段结束时人工变量的和不等于0（大于EPS），程序给出了INFEASIBLE判定。

这两个边界情况是考验两阶段法实现是否正确的关键。很多网上抄来的代码在这两种输入下要么直接崩掉，要么给出错误结果。特别是无可行解的判定，经常有人忽略了第一阶段目标值必须等于0这个条件。

测试时我还对照了Python的scipy.optimize.linprog，随机生成了几百个小规模线性规划问题，两边结果做比较。整个过程帮我抓到了两个隐藏bug：一个是最小化/最大化方向没处理对，另一个是存在退化问题的输入上迭代次数异常增长。这种“自己实现和成熟库对照”的测试思路，我认为很值得推荐。

6. 常见问题排查与性能优化方向

6.1 常见问题速查表

我把实际调试中最容易踩的坑整理成一个速查表，方便大家快速定位问题：

除了这些，我还要提醒一个C++工程上的细节：我最初把EPS作为普通常量放在类内部，后来发现不同算例对精度的要求其实不同，大数值的输入需要更大的EPS，小数值的输入需要更小的EPS。这个方法本身没有做自适应，但至少把EPS设计成了一个参数，方便调用方根据自己的数据规模调整。

6.2 性能优化与扩展思路

单纯形算法理论上是指数级的，但实际工程中因为大规模线性规划几乎都是稀疏的，配合合理的选主元策略，收敛速度通常很好。这个项目用vector<vector >存储，对几百行几百列的问题完全够用。

如果以后要应对更大规模的问题，有几个方向可以做：一是改成稀疏矩阵存储，只保存非零元素，pivot时只更新受影响的行；二是用revised simplex，维护基矩阵的逆而不是整个表格，配合LU分解，数值稳定性更好，迭代速度也更快；三是引入列生成，对列数极多但实际用不到的问题，动态生成列而不是一开始全部加载。

从功能扩展角度，同一个框架可以加不少东西。灵敏度分析就是很自然的下一步：最优解出来了，约束松弛变量的值可以直接解读成影子价格，对应的是资源每增加一个单位目标函数的改善量。整数规划分支定界也常常用单纯形法作为下层求解器，C++版本可以很方便地嵌入进去。

7. 最后分享一点个人体会

写这个程序最大的收获，不是学会了单纯形法，而是深刻理解了数值计算中“理论正确”和“工程正确”的巨大差距。教科书上的算法流程我早就背熟了，但写出一个真正能在各种边界输入下稳定运行的版本，完全是另一回事：EPS设多少、什么时候跳过一次消元、退化要不要处理，这些细节教科书不会告诉你，只能靠测试和踩坑来积累。

如果你也要自己实现这个算法，我的建议是先从2x2的手算例子开始，每一步都打印单纯形表，亲眼看到数据怎么变化。等基本的迭代流程跑通了，再加两阶段法，加退化处理，加无界无解判定。每加一层，就用和成熟库对照的随机测试做回归验证。按照这个顺序来，两天时间完全够用。

本文还有配套的精品资源，点击获取