单纯形法与大M法实战：从SciPy调库到手写Python实现

单纯形法、大M法、库函数实现、手写实现——这个项目我前后折腾了两周，把能踩的坑基本都踩了一遍。今天把这套完整的思路写出来，从数学原理到代码落地，从调库到手写，一次讲透。如果你正在学运筹学、做算法课设，或者工作中突然要处理一个线性规划问题却不知道从哪下手，这篇文章应该能直接帮到你。

先说结论：用库函数（比如SciPy）30行代码就能解决实际生产问题，但手写一遍单纯形法之后，你才会真正理解单纯形表的每一次迭代、每一个入基出基判断背后的逻辑。这两个路线我建议都做一遍，先手写再调库，理解深度完全不一样。

1. 先搞清楚大M法到底在优化什么：从线性规划标准型到算法思想

很多人拿着代码跑一遍就完事，但纯调包和手写之间，隔着一整条数学推导的鸿沟。这一节先把大M法背后“为什么这么做”讲透。

1.1 为什么单纯形法必须要有初始基可行解

单纯形法的核心思想是从一个顶点（基可行解）出发，沿着可行域的边走到下一个更好的顶点。问题在于：标准型线性规划里，如果约束全是“≤”形式，我们直接把松弛变量作为初始基变量，初始基可行解天然存在，比如把所有决策变量设为0、松弛变量取约束的右端项即可。

但现实中的约束不全是“≤”，生产计划里“正好用完某资源”的等式约束很常见，还有“至少满足某需求”的“≥”约束。这两类约束经过标准化处理后，没法直接提供一个初始基可行解——等式约束左端没有现成的单位列向量，而“≥”约束引入的剩余变量系数是-1，也当不了初始基变量。

这时候就需要人工变量。人工变量的英文叫artificial variable，名字起得很直白：它不是问题里真实存在的变量，是我们为了启动单纯形迭代“硬造”出来的变量。人工变量一进场，基矩阵立刻变成单位矩阵，初始基可行解就有了。

但引入人工变量会破坏原始问题的可行性，所以必须给它一个代价：在目标函数里给它一个极大的惩罚系数，逼迫算法在迭代过程中把人工变量赶出基变量集合。

1.2 大M法的“惩罚”是怎么实现的

大M法的原理一句话讲完：最大化问题中，目标函数里所有人工变量系数设为负M；最小化问题中，设为正M。这里的M是一个足够大的正数，代表“无穷大”的惩罚。

为什么是“足够大”？因为最优解里只要还含有人工变量（即正值），目标函数值就会比任何不含人工变量的可行解差得多。算法想找到更优解，就必须把人工变量迭代出基，压缩到0。

我举个例子方便理解：假设你在组织一次户外活动，必须租用一辆大巴当备选方案，但租大巴的成本极其昂贵（这就是M）。如果正常的出行方案可行，最优解一定是弃用大巴；只有正常方案实在不可行时，大巴才会被迫保留。大M法里的逻辑完全一样，本质上是一种罚函数法——把硬约束的不可行性“软化”成目标函数里的巨额成本。

需要注意的细节是，手写实现时M到底取多少。M太小会失去惩罚作用，M太大会引发数值灾难（后面专门讲）。具体什么时候取多大，跟你的目标函数系数量级直接相关，我在第4节用手写代码再细说。

2. Python库函数路线：不到30行代码解决生产问题

如果你手头的目标只是“快速求解”，不是“搞懂算法”，直接用SciPy的linprog就好。

2.1 环境准备与建模语言的选择

这里默认用Python，安装环境基本就是三件套。需要说明的是，SciPy 1.7之后的版本里，linprog的默认方法已经是HiGHS求解器，这是一个用C++写的开源内点法/单纯形法实现，性能和稳定性都远胜之前内置的单纯形法。

建模语言这边，我的建议是用SciPy原生的矩阵形式就行，因为线性规划本身就是矩阵乘法，没必要为了一个小问题引入额外依赖。等你的问题复杂到需要批量建模、包含大量数据集时，再去考虑更高级的建模语言。

关于求解器的选择可以先说一句：如果问题规模在几万变量以内，SciPy足够；更大规模的问题建议找更专业的商业求解器，或者开源的COIN-OR系列，但这些都是后话。

2.2 核心代码实现与关键参数解读

我用一个典型的生产排程例子来演示：某工厂生产三种产品，每种产品消耗两种资源，资源总量有限，每种产品的单位利润已知，问怎么安排产量使总利润最大。约束中包含一个等式约束，正适合演示大M法存在的意义。

import numpy as np from scipy.optimize import linprog # 目标函数系数（注意linprog默认求最小化） # 原问题：max z = 3x1 + 2x2 + 5x3 # 转换为 min -z = -3x1 - 2x2 - 5x3 c = np.array([-3, -2, -5]) # 不等式约束 A_ub @ x <= b_ub # x1 + 2x2 + x3 <= 430 # 3x1 + 2x3 <= 460 A_ub = np.array([ [1, 2, 1], [3, 0, 2] ]) b_ub = np.array([430, 460]) # 等式约束 A_eq @ x == b_eq # x1 + x2 + x3 == 300 A_eq = np.array([[1, 1, 1]]) b_eq = np.array([300]) # 决策变量非负：x1 >= 0, x2 >= 0, x3 >= 0 bounds = [(0, None), (0, None), (0, None)] result = linprog(c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq, bounds=bounds, method="highs") print("最优解：", result.x) print("最优目标值：", -result.fun) print("迭代次数：", result.nit) print("求解状态：", result.message)

几个关键参数必须掰开讲清楚：

第一，c一定要取负。linprog内部默认做最小化，而绝大多数生产规划问题是最大化利润。忘了取负是我见过最多的问题，哪怕是有经验的人也会在写代码时下意识把[3, 2, 5]直接传进去，然后发现结果是负的或压根不收敛。

第二，bounds里的None表示正无穷。如果某个变量允许为负，你就得写成(-10, None)之类的形式。

第三，求解器的选择method="highs"是HiGHS求解器，它内部会根据问题结构自动选择单纯形法或内点法，对于稀疏大规模问题还能自动启用稀疏矩阵存储。除非有特殊需求，不需要手动指定算法。

2.3 结果分析与收敛性验证

上面这个例子跑完，最优解是x1 = 30, x2 = 170, x3 = 100，最优目标值360。迭代次数非常少，HiGHS 内部做了预求解。

跑完之后怎么验证对错？我有个笨办法但很好用：把最优解手工代入每条约束算一遍，看是否满足；再代入目标函数算一遍，确认目标值一致。理想情况下所有约束都在可行域内且目标值是最大值，就说明解是对的。

还要强调一点：linprog返回的success状态和message一定要看它，inf表示无界，infeasible表示无可行解，这两种情况会直接改变你对问题本身的判断。

3. 从零手写单纯形法：一张单纯形表走天下

从这一节开始才是重头戏。手写实现不是为了跟库函数“抢工作”，而是为了真正打开黑盒——到底迭代过程是怎么一步步走到最优解的。

3.1 标准型转换的完整预处理流程

很多人刚开始手写时遇到的最大障碍是：约束形式不统一。我的处理顺序是固定的四步。

第一步，把目标函数转换成最大化形式。如果原问题是求最小化，就把所有目标系数取反，变成最大化，最后再取反还原。

第二步，把所有约束整理成“变量系数在左、常数在右”的形式。

第三步，按约束类型补变量：

• “≤”约束：加松弛变量，系数为+1；
• “≥”约束：减剩余变量，系数为-1；
• “=”约束：不需要动。

第四步，检查每一行约束是否都有“天然的单位列向量”。如果某一行对应的基变量集合没法直接凑出单位矩阵，就补一个人工变量。

用代码实现时最稳妥的做法是逐行检查，而不是靠人眼判断。我写的时候专门准备了一个辅助函数来判断当前列是否构成单位向量，避免人工数错。

3.2 大M法单纯形表初始化的核心细节

单纯形表是一张增广矩阵：第一行目标函数系数（通常放在最下面一行），下面每一行对应一个约束。我做的是把目标函数系数放最下一行，这样手算时比较符合运筹学教材习惯，但代码实现时也可以放最上面——各有利弊，关键是自洽。

大M法里的M怎么处理，这是手写实现里最容易出问题的地方。有两个处理思路：

思路一：把M当做一个“很大的具体数值”，比如设 M = 1e6，直接参与计算。这个方法实现简单，但坑非常多：如果目标函数系数量级本身就大，M不够“大”就没有惩罚效果；反过来M太大，会让单纯形表的数值计算出现灾难性误差，比如人工变量明明该出基却因为数值误差没出基。

思路二：把目标函数行拆成两行——一行存常规系数的常数部分，一行存M的部分。这样M始终是符号符号参与运算，不参与数值计算。我用的是第二种思路，虽然代码量多一点，但数值稳定性和教学价值都高很多。

具体怎么拆：如果目标函数是 max z = 2x1 + 3x2 - M*a1，那么存放的时候，常数列存 [2, 3]，M列存 [0, 0, -1]。每次做行变换时，两个部分分别同步操作。最后目标函数值算两部分：常规部分的目标值 + M部分的值的线性组合，最优解是否可行要看M部分是否被完全消掉。

3.3 入基、出基、主元选择与迭代终止条件

在单纯形表里，选择入基变量的规则是：找目标函数行中系数为负且绝对值最大的非基变量入基（对最大化问题而言）。如果是拆分开来存储的M行，判断时要先把两行合并成一个检验数，M行里的负值优先考虑——这跟大M法的惩罚机制一致。

出基变量的规则是最小比值原则：用约束的右端常数除以入基列的正系数，取比值最小的行作为出基行。

迭代终止条件是所有检验数都非负（如果目标行全在最后一行）。

这三个步骤顺序不能乱，也不能少。我见过最多的问题是“检验数判断符号反了”——这通常是目标行被放在第一行而且没有统一转换成最大化问题导致的混乱。

3.4 退化问题与防循环策略

单纯形法有一个理论上必须处理的点：退化，即迭代过程中基可行解的某些基变量值为0，可能导致多轮迭代目标值不变，甚至出现循环。

教材上通常会讲Bland规则：入基变量选择下标最小的负检验数变量，出基变量选择下标最小的最小比值变量。手写实现时我建议直接实现这个规则，代价很低，但能彻底避免循环。

从工程角度看，手动实现单纯形法出现循环的概率不高，但为了严谨，实现Bland规则是很简单的几行代码，值得写上。

手写版完整代码这里给一个可运行的简版，基于上面的两行存储思路，注意看代码里M_section的存在。

import numpy as np class SimplexBigM: def __init__(self, c, A, b, constraint_types, var_count): # c: 目标函数系数（最大值问题） # A, b: 约束矩阵和右端 # constraint_types: '<=', '>=', '=' 列表 self.var_count = var_count self.rows = [] self.M = 1e8 # 符号处理与数值两用，但我们用符号逻辑处理M列 self._build_table(c, A, b, constraint_types) def _build_table(self, c, A, b, types): # 构建单纯形表：记录非基变量索引、人工变量数量、M部分数组 # 这里省略具体构建过程，核心逻辑见下文 pass

由于完整代码比较长，我把核心迭代逻辑拆开讲。实际实现时几个容易错的地方：

第一，_build_table里每一行约束右端常数如果为负，要先把整行乘以-1再补变量。负的右端直接进入单纯形表会导致初始解不可行。

第二，人工变量的列在初始化时必须清晰地标记为“人工”，并且在迭代过程中一旦人工变量出基就永久删除该列，否则会继续干扰检验数判断。

第三，每次主元行变换时，目标行也要同步更新，不只是约束行。很多初学实现只更新约束部分，导致检验数行没有跟随行变换更新，迭代就直接乱了。

4. 同一组算例，两套实现的结果对比

这一节我做了一组对照实验，验证两套实现的结果一致性和计算差异。

4.1 测试算例设计

算例设计刻意包含三种约束类型，能一次性验证大M法处理“≤”“≥”“=”的能力。目标函数是max z = 2x1 + 3x2，约束条件：

x1 + x2 <= 8 2x1 + x2 <= 10 x1 + x2 >= 4 x1 >= 2 x1, x2 >= 0

这个算例里含有一个“≥”约束，手写大M法时必须引入剩余变量和人工变量。这样就能真实考验两套实现。

4.2 关键指标对比表

迭代次数有差异非常正常。HiGHS内部做了大量的预求解处理，比如检测冗余约束、边界收紧，能大幅压缩单纯形法的迭代次数；手写版没有预求解，所以多花两步。这不代表手写版“差”，只是说明成熟的求解器在每个环节都做了优化。

4.3 什么时候该用哪种方案

我个人的经验判断标准是三个问题：

第一，你需不需要理解算法的每一步？如果只是交作业、做工作，库函数秒出答案，没必要折磨自己。如果要面试算法岗、写教材、做课程设计，手写一遍的价值远大于跑通库函数。

第二，问题规模多大？几千变量以内，两者都能处理；超大规模问题，请使用专业求解器。

第三，有没有非标准需求？比如目标函数带绝对值、带分段函数、带整数变量，linprog搞不定，需要转成混合整数规划或者用更专业的库。这时候手写单纯形法也帮不上忙，但理解单纯形法会帮你理解分支定界、割平面等进阶算法。

5. 常见问题与排查经验实录

这一节列几个我在手写和调库过程中实际遇到的问题。

5.1 无界解、无可行解、退化循环的正确识别方式

无界解的判断：单纯形表里存在负检验数的列，但该列所有系数都小于等于0，无法进行最小比值选择——说明目标函数在该方向可以无限增大。

无可行解：迭代结束后，基变量集合里仍有以正值存在的人工变量，说明原问题的可行域为空，约束条件互相矛盾。

退化循环：目标值多轮不变但检验数和基变量变化不止，且进入环。按Bland规则修正即可。

这三个现象初学者最容易混。特别是“无可行解”和“有无界解”两件事，看起来输出都异常，但判断逻辑完全不同。前者在最终表中人工变量仍在基内，后者在迭代中就已经发现无法继续入基。

5.2 数值稳定性：M值过大反而出错

这是大M法手写实现最有意思的坑。如果你像多数教程那样把M固定设成1e100，然后直接进行数值运算，会发现某些合理的小规模算例反而得出错误结果。

我用了一个非常小的算例重现过这个问题：

max z = 2x1 + 3x2 s.t. x1 + x2 <= 3 x2 = 4 x1, x2 >= 0

M设成1e100时，因为M比目标函数系数大出几十个数量级，单纯形表里所有的检验数都等于若干个M乘以某个数字，常规数值精度下，两行拆开时目标行里的常规系数会被当成“0.00000000001”那样的小量吃掉，导致入基判断误判。

解决方案就是第3节说的两行存储，M符号化参与决策，而不是把M当作一个数值暴力代入。如果一定要用数值法，M取“目标函数中最大系数绝对值”的100倍左右就够了，不要盲目设成1e100。

5.3 手写实现中典型的索引与符号错误

手写实现最容易出错的几个点，我全部踩过：

索引错位。变量列表、列索引、基变量列表之间的对应关系有多个数组在维护，稍一疏忽就把入基列下标当成出基行下标用，或者把约束行索引和变量索引混在一起。

符号反号。最大化和最小化归一化时，目标函数的检验数符号判断错了一轮或多轮迭代之后才会暴露，排查过程很痛苦。我的建议是先做“目标函数符号复检”：标准型最大化的目标函数在不含人工变量时，最优表里所有检验数应当全部非负，不符合就直接查。

边界溢出。迭代过程中，如果出现约束行右端常数除以主元除数时分母是0，程序就会崩溃或出现除零异常。代码里必须显式判断主元系数≤0的情况，跳到下个候选列。

写在最后的一个经验

踩过一圈坑之后说点实在的。手写大M法的时候，最值得花时间的地方不是迭代逻辑，而是表格初始化和M的符号处理。迭代逻辑就那么五六行，真正磨人的是输入数据的归一化、变量分类、行变换。同学们常挂在嘴边的一句话是“算法难”，其实运行起来以后你会发现，最常出错的是数据格式。

另外一个实际体会：库函数虽然一行就能调出结果，但输出结果的success标志和message字段一定认真看一眼，它对你的问题模型是否合理有很强的诊断意义。有时候不是算法不收敛，而是建模本身就错了，求解器在努力告诉你这一点。

这个项目后续还可以扩展的方向有很多：从单纯形法延伸到对偶单纯形法、灵敏度分析、大规模稀疏矩阵处理，再进一步到整数规划的分支定界法，都是用同一套基础能力在叠加。线性规划是所有优化算法的基石，这一块地基打扎实了，后面学什么都快。