资讯详情

C++手搓单纯形法:从建模到求解的完整计算程序

📅 2026/10/23 13:05:26 | 华诺云谱 👁 阅读
C++手搓单纯形法:从建模到求解的完整计算程序
简介这份资源是面向算法学习者与C开发者的单纯形法完整实现程序用于求解线性规划问题的最优解。单纯形法由丹齐格于1947年提出通过迭代在基本可行解之间转换逐步逼近目标函数最优值是运筹学与数值计算中的经典算法。资源包共3个文件包含2个cpp源文件与1个h头文件压缩包约2KB代码围绕单纯形表、基变量与非基变量替换、迭代终止条件等核心环节展开结构紧凑便于直接阅读与二次开发。目前已有742人学习下载适合作为课程实验、算法练习或工程入门的参考实现。读者可从中理解线性规划标准形式的转换思路、单纯形表的构造方式以及C中矩阵与向量的组织技巧并借鉴异常处理与退出条件检查等工程细节为后续性能优化和测试调试打下基础。1. 用 C 手搓单纯形法从建模到求解的完整计算程序很多做运筹学作业或线性规划求解的开发者第一次接触单纯形法都是在纸上推演表格一旦变量数超过五个手工迭代基本就崩了。这份 C 实现的单纯形算法计算程序解决的正是这个痛点把标准型线性规划问题丢进去程序自动完成基变换、入基出基判断、迭代收敛最后输出最优解和目标函数值。它适合两类人一类是正在学运筹学、需要验证手工计算结果的学生另一类是想把线性规划求解嵌进自己工程里的 C 开发者拿过去改改输入输出就能用。整个项目不依赖第三方数学库纯标准 C 实现编译即跑这也是它最实在的地方。2. 项目结构与技术选型为什么用纯 C 而不是调库2.1 目录组织与文件职责拿到这份源码第一件事是看目录。这类单纯形法程序的结构通常很紧凑常见做法是拆成三个部分主入口、算法核心、以及输入输出辅助。我一般会按下面的方式组织读者对照自己手上的源码可以快速定位simplex/ ├── main.cpp // 程序入口负责读取问题、调用求解、打印结果 ├── simplex.h // 单纯形算法类声明定义数据结构与接口 ├── simplex.cpp // 算法核心实现迭代逻辑都在这里 └── README.md // 编译说明与输入格式约定main.cpp只做三件事构造问题、调用SimplexSolver::solve()、输出结果。真正的计算全部压在simplex.cpp里。这种拆分的好处是你想把求解器嵌到别的项目只需要拿走simplex.h和simplex.cppmain.cpp直接丢掉换自己的调用代码。2.2 数据结构设计表格法还是修正单纯形单纯形法有两种主流实现路线表格法Tableau和修正单纯形法Revised Simplex。这份程序用的是表格法原因是它直观、易调试适合教学和中小规模问题。核心数据结构是一张二维表常见做法是用std::vectorstd::vectordouble存整张单纯形表行代表约束列代表变量加松弛变量最后一列存右端项 b。// simplex.h class SimplexSolver { public: // 构造时传入目标函数系数、约束系数矩阵、右端项 SimplexSolver(const std::vectordouble c, const std::vectorstd::vectordouble A, const std::vectordouble b); bool solve(); // 执行迭代返回是否找到最优解 double getOptimalValue() const; // 取最优目标函数值 std::vectordouble getSolution() const; // 取最优解向量 private: int m, n; // 约束数、原始变量数 std::vectorstd::vectordouble tableau; // 单纯形表 std::vectorint basis; // 基变量索引 void buildTableau(); // 初始化表格加入松弛变量 int findEnteringVariable(); // 入基变量检验数为正且最大 int findLeavingVariable(int enterCol); // 出基变量最小比值检验 void pivot(int row, int col); // 高斯消元做基变换 };这里几个参数要说明白m是约束个数n是决策变量个数。tableau的列数是n m 1多出来的m列是松弛变量最后 1 列是右端项。basis记录每一行当前哪个变量在基里初始时全是松弛变量。findEnteringVariable返回检验数最大的列最大化问题findLeavingVariable用最小比值法则避免可行解跑出可行域。2.3 为什么不用第三方库有开发者会问为什么不直接上某优化库原因很实际这份程序的目标是教学和轻量嵌入。第三方库动辄几十兆编译链复杂而单纯形法本身逻辑清晰纯 C 实现不到三百行。对于变量数在几十以内的线性规划问题表格法的性能完全够用。当然如果你的问题规模上千那确实该换修正单纯形法或专业求解器这是后话。3. 核心算法拆解一次迭代到底发生了什么3.1 初始化把标准型变成单纯形表单纯形法的第一步是把问题转成标准型。假设原始问题是最大化c^T x约束是Ax ≤ bx ≥ 0。程序通过加入松弛变量把不等式变成等式构造初始基可行解。这一步在buildTableau()里完成。void SimplexSolver::buildTableau() { // 表格列数 决策变量 n 松弛变量 m 右端项 1 tableau.assign(m, std::vectordouble(n m 1, 0.0)); for (int i 0; i m; i) { for (int j 0; j n; j) { tableau[i][j] A[i][j]; // 填入约束系数 } tableau[i][n i] 1.0; // 第 i 行松弛变量系数为 1 tableau[i][n m] b[i]; // 右端项 } // 目标函数行检验数行放在最后一行初始为 -c std::vectordouble objRow(n m 1, 0.0); for (int j 0; j n; j) { objRow[j] -c[j]; // 最大化问题取负便于统一判断 } tableau.push_back(objRow); // 初始基全是松弛变量 basis.resize(m); for (int i 0; i m; i) basis[i] n i; }逻辑说明tableau比约束多一行最后一行是检验数行。把-c放进去是因为最大化问题里检验数为负表示还有改进空间迭代到全非负就最优。basis初始化为松弛变量索引对应单位矩阵保证初始解可行。参数上A[i][j]是第 i 个约束里第 j 个变量的系数b[i]必须非负否则初始基不可行需要两阶段法这份程序默认输入满足b ≥ 0。3.2 入基与出基最小比值法则的代码落地迭代的核心是选谁进基、谁出基。入基选检验数最负的那一列对应最大化改进最快出基用最小比值法则保证右端项在消元后仍然非负。int SimplexSolver::findEnteringVariable() { int enterCol -1; double minVal -1e-9; // 容差避免浮点误差误判 int lastRow m; // 检验数行索引 for (int j 0; j n m; j) { if (tableau[lastRow][j] minVal) { minVal tableau[lastRow][j]; enterCol j; } } return enterCol; // 返回 -1 表示已最优 } int SimplexSolver::findLeavingVariable(int enterCol) { int leaveRow -1; double minRatio 1e18; for (int i 0; i m; i) { double coef tableau[i][enterCol]; if (coef 1e-9) { // 只考虑正系数否则比值无意义 double ratio tableau[i][n m] / coef; if (ratio minRatio) { minRatio ratio; leaveRow i; } } } return leaveRow; // 返回 -1 表示问题无界 }参数说明1e-9是浮点比较容差单纯形法迭代多了容易积累误差没有容差会出现死循环或误判。findEnteringVariable返回 -1 时说明所有检验数非负当前解就是最优解。findLeavingVariable返回 -1 时说明入基列所有系数非正问题无界目标函数可以无限增大。这两个返回值在solve()里分别对应「最优」和「无界」两种终止状态。3.3 基变换pivot 消元的实现细节选定入基列和出基行后做高斯消元把入基列变成单位向量。这一步是数值计算最容易翻车的地方主元太小会导致误差爆炸。void SimplexSolver::pivot(int row, int col) { double pivotVal tableau[row][col]; // 归一化主元行 for (int j 0; j n m; j) { tableau[row][j] / pivotVal; } // 其他行消元包括检验数行 for (int i 0; i m; i) { if (i row) continue; double factor tableau[i][col]; if (std::fabs(factor) 1e-12) continue; for (int j 0; j n m; j) { tableau[i][j] - factor * tableau[row][j]; } } basis[row] col; // 更新基变量索引 }逻辑说明先把主元行整体除以主元让主元变成 1然后对其他每一行减去factor倍的主元行把该列其他元素清零。注意循环上界是i m因为检验数行也要消元。basis[row] col更新基索引保证后续能正确还原解向量。主元判断用1e-12是防止除以极小数实际工程里如果主元绝对值小于阈值通常要触发重新选主元或报数值不稳定。3.4 主循环与终止条件把上面几个函数串起来solve()就是一个 while 循环找入基、找出基、pivot直到最优或无界。bool SimplexSolver::solve() { buildTableau(); const int MAX_ITER 10000; // 防止死循环的后悔药 int iter 0; while (iter MAX_ITER) { int enterCol findEnteringVariable(); if (enterCol -1) return true; // 最优 int leaveRow findLeavingVariable(enterCol); if (leaveRow -1) return false; // 无界 pivot(leaveRow, enterCol); } return false; // 迭代超限视为失败 }MAX_ITER是血泪经验浮点误差可能让检验数在 0 附近反复横跳没有迭代上限程序可能永远不退出。设成 10000 对中小规模问题足够真超了说明输入有问题或需要换数值更稳的算法。4. 编译运行与输入格式从源码到可执行文件4.1 编译命令与依赖这份程序没有外部依赖只要有 C11 以上的编译器就能编。Linux 或 macOS 下g -stdc11 -O2 -o simplex main.cpp simplex.cpp ./simplexWindows 下用 MinGW 或 Visual Studio 命令行同理把源文件一起编译即可。-O2开优化对单纯形法这种循环密集的代码有实际提升。如果编译报fabs未定义检查有没有#include cmath这是新手最常踩的坑。4.2 输入格式约定程序读取问题的方式通常有两种硬编码在main.cpp里或者从标准输入读。常见做法是前者方便调试。一个典型输入长这样// main.cpp 中的问题定义 std::vectordouble c {3, 5}; // 目标函数系数 std::vectorstd::vectordouble A { {1, 0}, // 约束1: x1 4 {0, 2}, // 约束2: 2*x2 12 {3, 2} // 约束3: 3*x1 2*x2 18 }; std::vectordouble b {4, 12, 18}; // 右端项 SimplexSolver solver(c, A, b); if (solver.solve()) { std::cout 最优值: solver.getOptimalValue() std::endl; auto sol solver.getSolution(); for (size_t i 0; i sol.size(); i) std::cout x i 1 sol[i] std::endl; }参数说明c是目标函数系数长度等于决策变量数A是约束系数矩阵行数等于约束数列数等于变量数b是右端项必须非负。这个例子是经典的线性规划入门题最优解应该是x1 2, x2 6目标值 36。跑出来对不上先检查b有没有负数再检查约束方向是不是都统一成了≤。4.3 结果还原从表格里取出解向量迭代结束后解藏在tableau的右端项列里但只有基变量有值非基变量全是 0。getSolution()要做的是按basis索引把值填回原变量位置。std::vectordouble SimplexSolver::getSolution() const { std::vectordouble x(n, 0.0); // 决策变量初始为 0 for (int i 0; i m; i) { int var basis[i]; if (var n) { // 只关心原始决策变量 x[var] tableau[i][n m]; } } return x; }逻辑说明遍历每一行看基变量索引var如果它落在原始变量范围[0, n)内就把该行右端项赋给对应变量。松弛变量的值不需要输出。getOptimalValue()更简单直接把检验数行最后一列取出来注意符号——因为初始化时目标函数取了负这里要取反。5. 避坑与常见问题迭代不收敛、结果对不上怎么办5.1 现象程序跑出无界但问题明明有最优解原因通常是约束方向没统一。单纯形法标准型要求所有约束是≤且b ≥ 0。如果原始问题里有≥或等式约束直接塞进去会导致初始基不可行程序误判无界。解决办法是先把问题转成标准型≥约束两边乘 -1 变≤但这样b可能变负需要引入人工变量走两阶段法。这份程序默认不做两阶段所以输入前必须自己保证b非负。5.2 现象迭代次数很多但检验数一直在 0 附近震荡这是浮点误差导致的退化循环。原因在于主元太小消元后误差被放大。解决方式有两个一是加容差把1e-9调大一点比如1e-7二是加 Bland 规则入基时选下标最小的负检验数列而不是选最负的能理论上避免循环。我一般会在findEnteringVariable里加个开关退化严重时切 Bland 规则。5.3 现象最优值和手工算的对不上差一点点先别怀疑算法九成是浮点精度问题。检查输出时有没有做四舍五入0.9999999和1.0在程序里是两个数。常见做法是在输出前统一处理小于1e-6的当成 0接近整数的保留整数。另外检查目标函数值有没有取反初始化时-c进去取结果时忘了取反符号就反了。5.4 现象变量数一多程序变慢甚至卡死表格法的复杂度是 O(m*n) 每轮迭代迭代轮数又和变量数相关规模上去确实会慢。如果变量数超过几百建议换修正单纯形法只维护基矩阵的逆内存和计算量都小一个量级。这份程序定位是教学和中小规模硬扛大规模问题不是它的设计目标。5.5 现象编译通过但运行崩溃提示 vector 越界多半是A的行列数和c、b对不上。A的行数必须等于b的长度列数必须等于c的长度。构造SimplexSolver时最好加个断言检查assert(A.size() b.size())和assert(A[0].size() c.size())把错误挡在构造阶段比运行到一半崩溃好排查得多。6. 二次开发改造从教学程序到工程可用的求解器这份程序拿来直接跑没问题但要嵌进真实工程还得做几处改造。第一处是输入输出解耦。现在问题硬编码在main.cpp改成从文件或标准输入读格式用简单的文本协议比如第一行放m n接着c再接着A和b。这样你的上层业务代码可以用任意语言生成问题文件C 程序只负责求解。第二处是错误处理。现在的solve()只返回 bool信息量太少。我一般会改成返回枚举OPTIMAL、UNBOUNDED、INFEASIBLE、MAX_ITER_REACHED再配一个getErrorMessage()。这样调用方能区分「问题本身无解」和「程序没算出来」排查方向完全不同。第三处是数值稳定性。把主元阈值、检验数容差做成可配置参数退化时自动切 Bland 规则。下面是一个改造后的接口示意enum class SolveStatus { OPTIMAL, UNBOUNDED, INFEASIBLE, MAX_ITER }; struct SolverConfig { double pivotTol 1e-9; // 主元最小阈值 double objTol 1e-7; // 检验数容差 int maxIter 10000; // 迭代上限 bool useBland false; // 是否强制 Bland 规则 }; class SimplexSolver { public: SimplexSolver(const std::vectordouble c, const std::vectorstd::vectordouble A, const std::vectordouble b, const SolverConfig cfg {}); SolveStatus solve(); // ... 其余接口不变 };改造完怎么验证我的习惯是准备一组已知答案的测试用例覆盖最优、无界、退化三种情况每次改完跑一遍。特别是退化用例构造一个检验数反复为 0 的问题看程序能不能在有限步内退出。这一步别省单纯形法的数值坑基本都藏在退化场景里。从那以后我每次改求解器核心都强制走一遍「已知答案用例 退化用例 大规模随机用例」三件套少一件都不敢提交。希望这份拆解帮到你拿到源码先跑通那个x12, x26的例子再动手改输入输出路子就顺了。本文还有配套的精品资源点击获取
📝

华诺云谱内容团队

资深建站顾问 · 行业研究员

10年+企业数字化服务经验,专注智能建站、SEO优化与品牌营销,持续输出建站技巧、行业洞察与营销干货,已帮助5000+企业实现数字化增长。

你可能需要的服务

订阅华诺云谱资讯周报

每周一封,精选建站技巧、SEO与营销干货,直达邮箱。已有 8,000+ 企业主订阅,助你少走弯路。

↑