简介面向最优化方法与运筹学课程的 MATLAB 程序包聚焦线性规划单纯形法、大M法与两阶段法的程序实现适用于正在学习最优化算法、准备运筹学考试或需要做课程设计的本科生与研究生。整份资料共 3 个 m 文件覆盖 main.m、twophase.m、MySimplex_method.m 三个模块压缩包仅 4KB轻量便于直接打开阅读。其中 MySimplex_method.m 同时融合大M法和单纯形法求解twophase.m 采用两阶段策略并通过调用 MySimplex_method.m 完成单纯形迭代main 函数则提供约束方程和目标函数的输入入口自动驱动整体求解。程序注释详细、逻辑清晰不仅展示了两阶段法与大M法之间的内在联系还帮助读者理解人工变量、基变量与最优性判断等关键步骤适合对照代码边读边调试也可作为算法对比实验的起点。目前已有 2541 人学习下载对希望快速上手 MATLAB 线性规划编程、想在短时间看清算法实现细节的学习者具有不错的参考价值。 那个压缩包传到课程平台之前我盯着程序输出的z 2.000000看了很久才把文件名从final改成final_v2。这门课叫最优化理论与方法研一第一学期作业要求手写实现线性规划里的三个经典算法单纯形法、大M法和两阶段法。如果你也正在为这份作业挠头或者刚接触最优化想知道单纯形法怎么从公式层面落地成能跑的代码这篇文章应该能帮你省下不少调试时间。我会从原理、代码结构、算法差异到调试经验一起讲清楚尽量把那些教材不写、老师不讲的细节补出来。1. 为什么大M法和两阶段法总是成对出现1.1 单纯形法能跑起来的隐含前提单纯形法解决的是线性规划标准形的问题min c^T x s.t. Ax b x ≥ 0这里最关键的一个隐含前提是模型里必须有一个现成的初始基可行解。什么情况下有现成的就是约束条件全是“≤”且右端项 b 非负的时候。比如约束是2x1 3x2 ≤ 6你加一个松弛变量 x3变成2x1 3x2 x3 6那么 x3 6 就是一个现成的基变量所有的 b 都大于等于零初始单纯形表直接就能填。但实际情况没这么理想。题目里一旦出现“≥”约束或者“”约束松弛变量那一套就失效了。“≥”约束需要减去一个剩余变量但这会让系数变成 -1没法构成单位矩阵等号约束更直接根本没有可加的松弛变量位置。这时候你就找不到初始基可行解单纯形法连第一脚都迈不出去。这就是我第一版程序写出来却在几个算例上直接报错的原因。我那时候没想明白以为单纯形法对所有线性规划是通用的结果在约束条件稍复杂的情况下初始单纯形表里根本没有单位矩阵算法还没开始迭代就已经处于非法状态。1.2 “凑”出初始基的两条技术路线为了给算法一个合法的起点教材里引入了人工变量的概念。人工变量不是模型原本的变量纯粹是为了凑出单位矩阵加的“假变量”。问题来了你加了假变量进去如果不管它最终解很可能停留在人工变量不等于零的位置那个解对原问题毫无意义。于是必须想办法让人工变量在迭代结束后变成零。大M法和两阶段法就是两种“逼人工变量离基”的手段一个靠惩罚一个靠分阶段。大M法的思路很暴力在目标函数里给人工变量赋一个极大的惩罚系数 M对最小化问题就是 M对最大化问题就是 -M。因为 M 是一个非常大的正数只要人工变量还在基里、还大于零目标函数值就会非常差单纯形法为了保证目标函数最优会本能地先把人工变量踢出基变量队伍。理论上迭代到最优时人工变量应该等于零这样得到的解就是原问题的可行解。两阶段法的思路则是分步处理。第一阶段先不看原来的目标函数单独设定一个新的目标函数“最小化所有人 工变量之和”在这个目标下强行搜索可行解。如果最后这个和是零说明原问题确实存在可行解且人工变量已经被赶出基了如果和大于零说明原问题根本无可行解直接判死刑。第二阶段再把原目标函数装回去在第一阶段结束的单纯形表基础上继续优化。课程把这两个方法放在一起教正是因为它们在解决同一个问题没有初始基时怎么办。但对写代码的人来说这两个方法实现难度并不相同。我第一次写的时候以为大M法更简单实际跑起来才发现两阶段法在某些场景下反而更稳下面我会详细讲为什么。2. 程序骨架单纯形表该怎么在代码里落地2.1 数据结构设计是第一个分水岭我见过很多同学的实现是用三个二维数组分别存约束矩阵 A、右端项 b、目标函数系数 c再单独维护基变量列表。这样做不是不行但逻辑绕来绕去容易在换基的时候把自己绕晕。我更推荐用一张完整的单纯形表把所有信息都放在同一个二维数组里行和列的含义固定下来# tableau 布局 # 行 0..m-1: 约束行col 0..n-1 是变量系数col n 是右端项 b # 行 m: 目标函数行的检验数加入 M 或阶段切换时会变用 Python 列表表达就是tableau [[0.0] * (nvars 1) for _ in range(nconstraints 1)]。基变量单独用一个整数列表 base_vars 记录比如base_vars[i]表示第 i 行约束对应的基变量是哪个原始变量索引。这样程序的结构非常清晰每次换基只需要改 base_vars 里的一个元素然后对 tableau 做一次高斯消去。核心的 pivot 操作也就是换基运算本质就是线性代数里的初等行变换。伪代码可以写成这样def pivot(tableau, row, col): # 先把主元所在行归一化 pivot_val tableau[row][col] tableau[row][:] [v / pivot_val for v in tableau[row]] # 再对其他所有行做消去把 col 列变成单位向量 for r in range(len(tableau)): if r row: continue factor tableau[r][col] if factor ! 0.0: tableau[r][:] [tv - factor * rv for rv, tv in zip(tableau[row], tableau[r])]这段代码看起来简单但我第一次写的时候漏了一个细节目标函数行也要参与消去。很多初学版本的实现只对约束行做消去而目标函数行是用“检验数”公式单独计算的这样两者容易脱节。统一用全表行变换处理虽然数据上多算了几次但逻辑一致性高很多后期调试不容易出错。2.2 入基出基规则与退化震荡的防卫单纯形法的迭代分三步找入基变量、找出基变量、做 pivot。入基变量的选择规则很简单对于最大化问题选检验数最大的正数那一列对于最小化问题选检验数最小的负数那一列。这里要注意很多教材的表述不一致有的默认最大化有的默认最小化程序里最好统一先转换成最小化标准形或者在每个函数注释里写清楚当前表是哪种形式不然几轮迭代下来符号就乱了。出基变量用的是最小比值原则。对每一行如果入基列系数 a_ic 0计算 b_i / a_ic取其中最小的那一行作为出基行。这样做是为了保证换基之后右端项仍然非负从而维持可行性。我调试时遇到的第一个诡异现象是程序在某些算例上会无限循环变量不断入基出基目标函数值却怎么都不变。这就是退化震荡问题。退化是指某个基变量取值为零导致多个行的最小比值相等此时出基行选择不唯一。如果每次选的都是同一个方向可能绕一圈又回到原点永远迭代不完。解决退化震荡的经典方法是 Bland 规则入基时选检验数满足条件且编号最小的列出基时选满足最小比值且编号最小的行。虽然 Bland 规则在实际算例中可能会让迭代步数变多但它能严格保证算法终止。我后来直接在代码里默认启用 Bland 规则的一个简化版本宁可多迭代几轮也不想程序当场死循环。3. 大M法与两阶段法在代码里的核心差异3.1 M值怎么选以及大M法的真正痛点大M法在代码实现上确实省事在目标函数行中人工变量对应的位置直接填上惩罚系数即可。如果是最小化问题且需要加人工变量把那个位置设成 M如果是最大化问题就设成 -M。然后正常走单纯形迭代就行不需要额外判断阶段切换。但问题出在 M 的取值上。我第一次图省事取了 M 100000结果在某个包含较多约束的算例上目标函数值在迭代前期疯狂跳动最终精度还差得离谱。原因很简单M 太大和正常的系数放在同一个浮点数体系里计算会出现严重的数值抵消。比如 M 与人工变量系数的乘积把原问题目标函数里那些较小的系数完全淹没了单纯形表的判断几乎只依赖 M 的符号数值误差一累积连检验数的正负都可能判断错。我的经验是M 不要拍脑袋取极大值而是根据模型里c和A的数据规模来定。一般取M max(abs(c)) * 10 max(abs(A)) * 10或者干脆M 1000左右先用小规模算例验证再根据误差放大。而且 M 并非越大越好这一点是教材里不会专门强调的。另一个大M法的隐藏问题是如果问题本身是无界的或者无可行解大M法的单纯形表会出现一些极端状态。比如无界解时入基列所有 a_ic ≤ 0找不到出基变量此时程序如果没有显式判断就会数组越界或者死循环我在验证用例时专门给这个分支加了提前退出。3.2 两阶段法第一阶段结束后的“换挡”细节两阶段法的代码比大M法绕但数值稳定性明显更好。第一阶段的目标函数是所有人工变量之和这在代码里很好实现额外加一个 m 个单位矩阵的列然后把目标函数行写成人 工变量的系数都是 1其他变量系数都是 0。迭代结束后需要检查目标函数值是否等于 0。这里有一个很多人第一次写都会掉的坑第一阶段结束后单纯形表里人工变量列虽然已经变成单位矩阵的一部分但基变量里可能还残留着某些人工变量。如果残留在基里且最终目标值为零说明有冗余约束这时候人工变量列理论上可以删掉但在代码里操作数组不是“删列”这么简单因为你还要同时维护 base_vars 的索引。阶段切换的正确姿势是这样的先判断第一阶段目标函数值的绝对值是否小于某个容差比如 1e-9如果可行把当前表中所有人工变量列标记为“忽略”或者物理删除这些列并同步压缩所有变量索引把目标函数行恢复成原问题的 c 向量同时重新计算当前基对应的检验数行然后再进入单纯形法的主循环。第 3 步是必须做的因为第一阶段结束时目标函数行里存储的还是人工变量之和的系数直接接着迭代就等于还在优化错误的目标。我一开始就是少了这步第二阶段跑了四五轮之后才发现目标值完全不对检查了很久才意识到目标行没有重 建。对比下来如果让我给学生作业推荐我会说大M法适合写起来快、算例规模小的场景两阶段法适合作为你真正想长期使用的求解器骨架。因为它不依赖那个很玄学的 M 值数值行为更加可控这也是很多商业求解器在预处理阶段会采用类似策略的原因。4. 调试实录最容易翻车的几个位置4.1 浮点数比较与 -0.0000 问题单纯形法几乎所有判断都是浮点数判断检验数是否为负、右端项是否为正、目标值是否为零。直接用 0是新手最容易犯的错因为单纯形表在迭代过程中会产生大量形如1e-16的浮点噪声。原本该是零的地方经过几次行变换后变成-3.8e-15这在数学上就是零但程序不知道。我后来的处理方式是设置一个全局容差EPS 1e-9所有判断都通过abs(val) EPS来判定。输出的时候也要防一手直接用 Python 的print(f{val:.6f})会把-0.000000这种值原样打出来看着像 bug其实是浮点数符号位的锅。在打印前统一加一个判断如果绝对值小于 EPS 就直接置为 0.0观感好很多排查问题时也少了很多干扰。4.2 人工变量没删干净导致的异常结果两阶段法第一阶段结束后如果只是把人工变量列的系数设成零而不删列第二阶段迭代时程序还会去算那些列的检验数一旦入基规则扫描到了人工变量列可能会把人 工变量重新拉回基里得出的结果自然就错了。这个问题在代码里很隐蔽因为单纯形表不会报错目标函数值看起来也正常但最优解和手写答案对不上。我当时花了差不多一个晚上才发现是这个问题。建议在设计数据结构的时候就给每列加一个active标记第一阶段结束后把所有人工变量列设成activeFalse入基变量扫描时跳过它们。这个标记比真正删列更安全因为你不需要频繁改数组长度索引也不容易错。4.3 退化案例被同一个循环折磨到怀疑人生前面说过退化震荡的问题。我在一个经典的退化算例上程序连续迭代了六十多轮还没停。后来加上了 Bland 规则才稳定收敛。一个小技巧调试阶段可以在每次迭代后打印当前基变量列表和目标函数值。当基变量列表重复出现时说明已经循环了。我在代码里加了一个“历史状态集合”每次新迭代前检查当前基变量元组是否已经出现过如果出现过就直接中断并提示退化。这虽然只是调试用不是严谨的数学处理但能很快帮你定位是不是进入了循环不至于傻等几十分钟。4.4 输入数据标准化先完成再谈优化这个问题很隐蔽但很致命如果你输入的模型不是标准形比如某个变量没有写非负限制、某个约束是x1 x2 10但你忘了加剩余变量程序会静默地给出一个你完全无法解释的结果。我的建议是主程序里写一个standardize函数先完成以下步骤再进入求解器把 min/max 统一如果是 max 问题把 c 取负转成 min自由变量拆成两个非负变量的差或者做变量替换对每个约束补松弛变量、剩余变量、人工变量并记录哪些是人工变量检查 b 向量如果有负数行把整行乘以 -1。这个标准化函数写完可以一劳永逸后面跑任何算例都安全了。千万不要把这三板斧写死在具体算例里否则换个题目又要改代码。5. 验证算例怎么确定程序没写错5.1 一个混合约束的测试用例我用的验证算例是教材上的经典案例因为手算答案可查能直接对拍max z 3x1 - x2 - x3 s.t. x1 - 2x2 x3 ≤ 11 -4x1 x2 2x3 ≥ 3 -2x1 x3 1 x1, x2, x3 ≥ 0注意这个算例里同时出现了“≤”、“≥”、“”三种约束完美覆盖了需要加人工变量的场景。第一个约束加松弛变量 x4第二个约束减剩余变量 x5 再加人工变量 x6第三个约束直接加人工变量 x7。手算的最终结果是x14, x21, x39目标值z2。你的程序只要跑出这个结果基本可以判定核心逻辑没问题。5.2 两阶段法跑通后的迭代表现我程序里加了verbose日志两阶段法跑这个算例时第一阶段的目标值从初始的2.0左右一路降到一个接近于0的数此时人工变量 x6、x7 都已经离基。第二阶段换回原目标后大约再迭代 4 轮左右检验数行全部 ≤ 0程序输出optimal solution found。对照最终单纯形表的输出基变量是x14.000, x21.000, x39.000z 值是2.000000。大M法在同一个算例上也能收敛到同样结果不过中间过程的迭代表格数值波动明显更大而且迭代次数多了两轮。如果你自己的程序这两种方法都能给出这个答案并且输出的小数和手算一致恭喜你这个程序可以拿去交作业了。再分享一个交叉验证的小技巧同一个算例先用小M值跑一遍再用大M值跑一遍两次结果是同一个才说明没被数值误差带偏。如果两次结果出现毫厘级偏差优先去检查浮点容差的设置而不是怀疑算法本身。5.3 边界条件的测试比常规算例更重要确保程序能正确解出最优解还不够你还得测两个反直觉的边界条件无界解和无可行解。无界解比如max z x1约束只要x1 ≤ 100但 x2 自由且目标里没它……需要构造一个入基列没有正值的场景。我在代码里专门加了判断如果入基列的所有约束系数都 ≤ 0就输出unbounded并终止。无可行解只需要用两阶段法人为构造一个第一阶段目标值不为零的算例比如两个互相矛盾的约束。我的程序会把这种情况判定为infeasible并退出而不是继续用错误的基变量迭代下去。这两个分支如果没写程序遇到这些情况会直接抛异常或给出荒谬结果让整个作业的完成度大打折扣。6. 代码扩展从交作业到真正可用的小求解器拿到这份代码之后我其实不建议你真的拿它去解大规模线性规划因为单纯形法的朴素实现效率有限。但你可以在这个基础上做不少扩展让你的作业代码看起来不只是“实现了”而是“有工程意识”。我比较推荐的扩展方向有三个加入对偶单纯形法对偶单纯形法在处理“初始不可行但检验数全满足条件”的问题时非常高效比如某些整数规划的松弛问题。实现思路和单纯形法几乎是镜像的替换入基出基逻辑即可。加入预处理自动删除全零行、合并冗余约束、检测自由变量。让用户输入原始模型就能跑而不是必须先手动转成标准形这对实际使用体验的提升非常明显。输出更丰富的统计信息比如迭代次数、退化检测次数、计算耗时。这些信息在分析复杂算例时很有用也是课程汇报里一个不错的加分亮点。我后来还给这个程序写了一个简单的文本解析器能直接读取形如自由格式的 LP 文件自动完成标准化、加入工变量等全部流程。从那之后这个程序就不只是课程作业了我遇到小型线性规划问题都会直接扔给它跑方便程度超过预期。最后再分享一个我踩过的坑如果是在 Windows 环境下用 C 语言写这套算法注意浮点数在 Release 模式下的优化选项可能会改变某些运算顺序导致结果和 Debug 模式不一致碰到这种玄学情况先把浮点模型设为严格模式再检查有没有未初始化的数组。用 Python 的话倒是没这个烦恼就是大算例速度会慢一些不过交作业完全够用了。本文还有配套的精品资源点击获取 线性规划单纯形法、大M法与两阶段法代码实现及调试全解析
简介面向最优化方法与运筹学课程的 MATLAB 程序包聚焦线性规划单纯形法、大M法与两阶段法的程序实现适用于正在学习最优化算法、准备运筹学考试或需要做课程设计的本科生与研究生。整份资料共 3 个 m 文件覆盖 main.m、twophase.m、MySimplex_method.m 三个模块压缩包仅 4KB轻量便于直接打开阅读。其中 MySimplex_method.m 同时融合大M法和单纯形法求解twophase.m 采用两阶段策略并通过调用 MySimplex_method.m 完成单纯形迭代main 函数则提供约束方程和目标函数的输入入口自动驱动整体求解。程序注释详细、逻辑清晰不仅展示了两阶段法与大M法之间的内在联系还帮助读者理解人工变量、基变量与最优性判断等关键步骤适合对照代码边读边调试也可作为算法对比实验的起点。目前已有 2541 人学习下载对希望快速上手 MATLAB 线性规划编程、想在短时间看清算法实现细节的学习者具有不错的参考价值。 那个压缩包传到课程平台之前我盯着程序输出的z 2.000000看了很久才把文件名从final改成final_v2。这门课叫最优化理论与方法研一第一学期作业要求手写实现线性规划里的三个经典算法单纯形法、大M法和两阶段法。如果你也正在为这份作业挠头或者刚接触最优化想知道单纯形法怎么从公式层面落地成能跑的代码这篇文章应该能帮你省下不少调试时间。我会从原理、代码结构、算法差异到调试经验一起讲清楚尽量把那些教材不写、老师不讲的细节补出来。1. 为什么大M法和两阶段法总是成对出现1.1 单纯形法能跑起来的隐含前提单纯形法解决的是线性规划标准形的问题min c^T x s.t. Ax b x ≥ 0这里最关键的一个隐含前提是模型里必须有一个现成的初始基可行解。什么情况下有现成的就是约束条件全是“≤”且右端项 b 非负的时候。比如约束是2x1 3x2 ≤ 6你加一个松弛变量 x3变成2x1 3x2 x3 6那么 x3 6 就是一个现成的基变量所有的 b 都大于等于零初始单纯形表直接就能填。但实际情况没这么理想。题目里一旦出现“≥”约束或者“”约束松弛变量那一套就失效了。“≥”约束需要减去一个剩余变量但这会让系数变成 -1没法构成单位矩阵等号约束更直接根本没有可加的松弛变量位置。这时候你就找不到初始基可行解单纯形法连第一脚都迈不出去。这就是我第一版程序写出来却在几个算例上直接报错的原因。我那时候没想明白以为单纯形法对所有线性规划是通用的结果在约束条件稍复杂的情况下初始单纯形表里根本没有单位矩阵算法还没开始迭代就已经处于非法状态。1.2 “凑”出初始基的两条技术路线为了给算法一个合法的起点教材里引入了人工变量的概念。人工变量不是模型原本的变量纯粹是为了凑出单位矩阵加的“假变量”。问题来了你加了假变量进去如果不管它最终解很可能停留在人工变量不等于零的位置那个解对原问题毫无意义。于是必须想办法让人工变量在迭代结束后变成零。大M法和两阶段法就是两种“逼人工变量离基”的手段一个靠惩罚一个靠分阶段。大M法的思路很暴力在目标函数里给人工变量赋一个极大的惩罚系数 M对最小化问题就是 M对最大化问题就是 -M。因为 M 是一个非常大的正数只要人工变量还在基里、还大于零目标函数值就会非常差单纯形法为了保证目标函数最优会本能地先把人工变量踢出基变量队伍。理论上迭代到最优时人工变量应该等于零这样得到的解就是原问题的可行解。两阶段法的思路则是分步处理。第一阶段先不看原来的目标函数单独设定一个新的目标函数“最小化所有人 工变量之和”在这个目标下强行搜索可行解。如果最后这个和是零说明原问题确实存在可行解且人工变量已经被赶出基了如果和大于零说明原问题根本无可行解直接判死刑。第二阶段再把原目标函数装回去在第一阶段结束的单纯形表基础上继续优化。课程把这两个方法放在一起教正是因为它们在解决同一个问题没有初始基时怎么办。但对写代码的人来说这两个方法实现难度并不相同。我第一次写的时候以为大M法更简单实际跑起来才发现两阶段法在某些场景下反而更稳下面我会详细讲为什么。2. 程序骨架单纯形表该怎么在代码里落地2.1 数据结构设计是第一个分水岭我见过很多同学的实现是用三个二维数组分别存约束矩阵 A、右端项 b、目标函数系数 c再单独维护基变量列表。这样做不是不行但逻辑绕来绕去容易在换基的时候把自己绕晕。我更推荐用一张完整的单纯形表把所有信息都放在同一个二维数组里行和列的含义固定下来# tableau 布局 # 行 0..m-1: 约束行col 0..n-1 是变量系数col n 是右端项 b # 行 m: 目标函数行的检验数加入 M 或阶段切换时会变用 Python 列表表达就是tableau [[0.0] * (nvars 1) for _ in range(nconstraints 1)]。基变量单独用一个整数列表 base_vars 记录比如base_vars[i]表示第 i 行约束对应的基变量是哪个原始变量索引。这样程序的结构非常清晰每次换基只需要改 base_vars 里的一个元素然后对 tableau 做一次高斯消去。核心的 pivot 操作也就是换基运算本质就是线性代数里的初等行变换。伪代码可以写成这样def pivot(tableau, row, col): # 先把主元所在行归一化 pivot_val tableau[row][col] tableau[row][:] [v / pivot_val for v in tableau[row]] # 再对其他所有行做消去把 col 列变成单位向量 for r in range(len(tableau)): if r row: continue factor tableau[r][col] if factor ! 0.0: tableau[r][:] [tv - factor * rv for rv, tv in zip(tableau[row], tableau[r])]这段代码看起来简单但我第一次写的时候漏了一个细节目标函数行也要参与消去。很多初学版本的实现只对约束行做消去而目标函数行是用“检验数”公式单独计算的这样两者容易脱节。统一用全表行变换处理虽然数据上多算了几次但逻辑一致性高很多后期调试不容易出错。2.2 入基出基规则与退化震荡的防卫单纯形法的迭代分三步找入基变量、找出基变量、做 pivot。入基变量的选择规则很简单对于最大化问题选检验数最大的正数那一列对于最小化问题选检验数最小的负数那一列。这里要注意很多教材的表述不一致有的默认最大化有的默认最小化程序里最好统一先转换成最小化标准形或者在每个函数注释里写清楚当前表是哪种形式不然几轮迭代下来符号就乱了。出基变量用的是最小比值原则。对每一行如果入基列系数 a_ic 0计算 b_i / a_ic取其中最小的那一行作为出基行。这样做是为了保证换基之后右端项仍然非负从而维持可行性。我调试时遇到的第一个诡异现象是程序在某些算例上会无限循环变量不断入基出基目标函数值却怎么都不变。这就是退化震荡问题。退化是指某个基变量取值为零导致多个行的最小比值相等此时出基行选择不唯一。如果每次选的都是同一个方向可能绕一圈又回到原点永远迭代不完。解决退化震荡的经典方法是 Bland 规则入基时选检验数满足条件且编号最小的列出基时选满足最小比值且编号最小的行。虽然 Bland 规则在实际算例中可能会让迭代步数变多但它能严格保证算法终止。我后来直接在代码里默认启用 Bland 规则的一个简化版本宁可多迭代几轮也不想程序当场死循环。3. 大M法与两阶段法在代码里的核心差异3.1 M值怎么选以及大M法的真正痛点大M法在代码实现上确实省事在目标函数行中人工变量对应的位置直接填上惩罚系数即可。如果是最小化问题且需要加人工变量把那个位置设成 M如果是最大化问题就设成 -M。然后正常走单纯形迭代就行不需要额外判断阶段切换。但问题出在 M 的取值上。我第一次图省事取了 M 100000结果在某个包含较多约束的算例上目标函数值在迭代前期疯狂跳动最终精度还差得离谱。原因很简单M 太大和正常的系数放在同一个浮点数体系里计算会出现严重的数值抵消。比如 M 与人工变量系数的乘积把原问题目标函数里那些较小的系数完全淹没了单纯形表的判断几乎只依赖 M 的符号数值误差一累积连检验数的正负都可能判断错。我的经验是M 不要拍脑袋取极大值而是根据模型里c和A的数据规模来定。一般取M max(abs(c)) * 10 max(abs(A)) * 10或者干脆M 1000左右先用小规模算例验证再根据误差放大。而且 M 并非越大越好这一点是教材里不会专门强调的。另一个大M法的隐藏问题是如果问题本身是无界的或者无可行解大M法的单纯形表会出现一些极端状态。比如无界解时入基列所有 a_ic ≤ 0找不到出基变量此时程序如果没有显式判断就会数组越界或者死循环我在验证用例时专门给这个分支加了提前退出。3.2 两阶段法第一阶段结束后的“换挡”细节两阶段法的代码比大M法绕但数值稳定性明显更好。第一阶段的目标函数是所有人工变量之和这在代码里很好实现额外加一个 m 个单位矩阵的列然后把目标函数行写成人 工变量的系数都是 1其他变量系数都是 0。迭代结束后需要检查目标函数值是否等于 0。这里有一个很多人第一次写都会掉的坑第一阶段结束后单纯形表里人工变量列虽然已经变成单位矩阵的一部分但基变量里可能还残留着某些人工变量。如果残留在基里且最终目标值为零说明有冗余约束这时候人工变量列理论上可以删掉但在代码里操作数组不是“删列”这么简单因为你还要同时维护 base_vars 的索引。阶段切换的正确姿势是这样的先判断第一阶段目标函数值的绝对值是否小于某个容差比如 1e-9如果可行把当前表中所有人工变量列标记为“忽略”或者物理删除这些列并同步压缩所有变量索引把目标函数行恢复成原问题的 c 向量同时重新计算当前基对应的检验数行然后再进入单纯形法的主循环。第 3 步是必须做的因为第一阶段结束时目标函数行里存储的还是人工变量之和的系数直接接着迭代就等于还在优化错误的目标。我一开始就是少了这步第二阶段跑了四五轮之后才发现目标值完全不对检查了很久才意识到目标行没有重 建。对比下来如果让我给学生作业推荐我会说大M法适合写起来快、算例规模小的场景两阶段法适合作为你真正想长期使用的求解器骨架。因为它不依赖那个很玄学的 M 值数值行为更加可控这也是很多商业求解器在预处理阶段会采用类似策略的原因。4. 调试实录最容易翻车的几个位置4.1 浮点数比较与 -0.0000 问题单纯形法几乎所有判断都是浮点数判断检验数是否为负、右端项是否为正、目标值是否为零。直接用 0是新手最容易犯的错因为单纯形表在迭代过程中会产生大量形如1e-16的浮点噪声。原本该是零的地方经过几次行变换后变成-3.8e-15这在数学上就是零但程序不知道。我后来的处理方式是设置一个全局容差EPS 1e-9所有判断都通过abs(val) EPS来判定。输出的时候也要防一手直接用 Python 的print(f{val:.6f})会把-0.000000这种值原样打出来看着像 bug其实是浮点数符号位的锅。在打印前统一加一个判断如果绝对值小于 EPS 就直接置为 0.0观感好很多排查问题时也少了很多干扰。4.2 人工变量没删干净导致的异常结果两阶段法第一阶段结束后如果只是把人工变量列的系数设成零而不删列第二阶段迭代时程序还会去算那些列的检验数一旦入基规则扫描到了人工变量列可能会把人 工变量重新拉回基里得出的结果自然就错了。这个问题在代码里很隐蔽因为单纯形表不会报错目标函数值看起来也正常但最优解和手写答案对不上。我当时花了差不多一个晚上才发现是这个问题。建议在设计数据结构的时候就给每列加一个active标记第一阶段结束后把所有人工变量列设成activeFalse入基变量扫描时跳过它们。这个标记比真正删列更安全因为你不需要频繁改数组长度索引也不容易错。4.3 退化案例被同一个循环折磨到怀疑人生前面说过退化震荡的问题。我在一个经典的退化算例上程序连续迭代了六十多轮还没停。后来加上了 Bland 规则才稳定收敛。一个小技巧调试阶段可以在每次迭代后打印当前基变量列表和目标函数值。当基变量列表重复出现时说明已经循环了。我在代码里加了一个“历史状态集合”每次新迭代前检查当前基变量元组是否已经出现过如果出现过就直接中断并提示退化。这虽然只是调试用不是严谨的数学处理但能很快帮你定位是不是进入了循环不至于傻等几十分钟。4.4 输入数据标准化先完成再谈优化这个问题很隐蔽但很致命如果你输入的模型不是标准形比如某个变量没有写非负限制、某个约束是x1 x2 10但你忘了加剩余变量程序会静默地给出一个你完全无法解释的结果。我的建议是主程序里写一个standardize函数先完成以下步骤再进入求解器把 min/max 统一如果是 max 问题把 c 取负转成 min自由变量拆成两个非负变量的差或者做变量替换对每个约束补松弛变量、剩余变量、人工变量并记录哪些是人工变量检查 b 向量如果有负数行把整行乘以 -1。这个标准化函数写完可以一劳永逸后面跑任何算例都安全了。千万不要把这三板斧写死在具体算例里否则换个题目又要改代码。5. 验证算例怎么确定程序没写错5.1 一个混合约束的测试用例我用的验证算例是教材上的经典案例因为手算答案可查能直接对拍max z 3x1 - x2 - x3 s.t. x1 - 2x2 x3 ≤ 11 -4x1 x2 2x3 ≥ 3 -2x1 x3 1 x1, x2, x3 ≥ 0注意这个算例里同时出现了“≤”、“≥”、“”三种约束完美覆盖了需要加人工变量的场景。第一个约束加松弛变量 x4第二个约束减剩余变量 x5 再加人工变量 x6第三个约束直接加人工变量 x7。手算的最终结果是x14, x21, x39目标值z2。你的程序只要跑出这个结果基本可以判定核心逻辑没问题。5.2 两阶段法跑通后的迭代表现我程序里加了verbose日志两阶段法跑这个算例时第一阶段的目标值从初始的2.0左右一路降到一个接近于0的数此时人工变量 x6、x7 都已经离基。第二阶段换回原目标后大约再迭代 4 轮左右检验数行全部 ≤ 0程序输出optimal solution found。对照最终单纯形表的输出基变量是x14.000, x21.000, x39.000z 值是2.000000。大M法在同一个算例上也能收敛到同样结果不过中间过程的迭代表格数值波动明显更大而且迭代次数多了两轮。如果你自己的程序这两种方法都能给出这个答案并且输出的小数和手算一致恭喜你这个程序可以拿去交作业了。再分享一个交叉验证的小技巧同一个算例先用小M值跑一遍再用大M值跑一遍两次结果是同一个才说明没被数值误差带偏。如果两次结果出现毫厘级偏差优先去检查浮点容差的设置而不是怀疑算法本身。5.3 边界条件的测试比常规算例更重要确保程序能正确解出最优解还不够你还得测两个反直觉的边界条件无界解和无可行解。无界解比如max z x1约束只要x1 ≤ 100但 x2 自由且目标里没它……需要构造一个入基列没有正值的场景。我在代码里专门加了判断如果入基列的所有约束系数都 ≤ 0就输出unbounded并终止。无可行解只需要用两阶段法人为构造一个第一阶段目标值不为零的算例比如两个互相矛盾的约束。我的程序会把这种情况判定为infeasible并退出而不是继续用错误的基变量迭代下去。这两个分支如果没写程序遇到这些情况会直接抛异常或给出荒谬结果让整个作业的完成度大打折扣。6. 代码扩展从交作业到真正可用的小求解器拿到这份代码之后我其实不建议你真的拿它去解大规模线性规划因为单纯形法的朴素实现效率有限。但你可以在这个基础上做不少扩展让你的作业代码看起来不只是“实现了”而是“有工程意识”。我比较推荐的扩展方向有三个加入对偶单纯形法对偶单纯形法在处理“初始不可行但检验数全满足条件”的问题时非常高效比如某些整数规划的松弛问题。实现思路和单纯形法几乎是镜像的替换入基出基逻辑即可。加入预处理自动删除全零行、合并冗余约束、检测自由变量。让用户输入原始模型就能跑而不是必须先手动转成标准形这对实际使用体验的提升非常明显。输出更丰富的统计信息比如迭代次数、退化检测次数、计算耗时。这些信息在分析复杂算例时很有用也是课程汇报里一个不错的加分亮点。我后来还给这个程序写了一个简单的文本解析器能直接读取形如自由格式的 LP 文件自动完成标准化、加入工变量等全部流程。从那之后这个程序就不只是课程作业了我遇到小型线性规划问题都会直接扔给它跑方便程度超过预期。最后再分享一个我踩过的坑如果是在 Windows 环境下用 C 语言写这套算法注意浮点数在 Release 模式下的优化选项可能会改变某些运算顺序导致结果和 Debug 模式不一致碰到这种玄学情况先把浮点模型设为严格模式再检查有没有未初始化的数组。用 Python 的话倒是没这个烦恼就是大算例速度会慢一些不过交作业完全够用了。本文还有配套的精品资源点击获取