企业官网建设流程全解析

单纯形法是线性规划最经典的求解方法，而大M法是单纯形法里处理“初始可行基不存在”这一难题的最直接补丁。我这次用Python把大M法从零手写了一遍，又把同一个算例丢给scipy.optimize.linprog去解，两边结果对照着看，把算法原理、代码细节、工程选型一次说透。这篇文章适合正在学运筹学、需要手撕算法的学生，也适合那些在项目里只想快速求解线性规划、但又被各种约束形式搞得头疼的工程师。看懂之后，你既能在作业里讲清楚大M法的每一步，也能在真实项目里准确调用现成工具。

1. 大M法解决的是什么问题

1.1 标准型与“随手能用的初始可行基”

线性规划的标准型长这样：最小化目标函数，所有约束写成等式，所有变量非负。数学上就是 min c^T x，满足 Ax = b，且 x >= 0。单纯形法的思路是沿着可行域的顶点不断跳跃，直到找到最优顶点。整个过程有一个前提：你得先有一个初始顶点，也就是一组初始基可行解。

如果约束全是“小于等于”形式，而且右端项 b 非负，事情就很简单。每个不等式加一个松弛变量，这些松弛变量天然构成一个单位矩阵，直接就是初始基。可现实中的约束哪有这么听话？只要出现“大于等于”或者“等于”的约束，你加完松弛变量后会得到系数为 -1 的列，它不能作为正的基变量，那就没有现成的单位矩阵可用了。

这个时候标准做法有两类：一类是两阶段法，先最小化人工变量之和；另一类就是大M法，在目标函数里直接给人工变量挂上巨额惩罚系数。本质是一样的，大M法用起来更直观，代码实现也更紧凑。

1.2 人工变量与M惩罚：把“不可行”拉回“可行”

人工变量这名字很形象。它不是原问题里真实存在的变量，而是我们强行加进去“凑数”的。它的作用就是填满单位矩阵的缺口，让单纯形法有一个可以起步的初始基。

问题在于，人工变量一旦加上，约束条件就被“放松”了，求解出来的点可能根本不在原问题的可行域里。为了逼着算法把这些人工变量赶出基，我们在目标函数里给它们加一个巨大的正系数 M。因为是最小化问题，算法每让一个人工变量留在基里，就要付出 M 的代价，只要 M 足够大，最优解里人工变量就会被压到零，从而回到原问题的可行域。

打个比方：人工变量像是临时工，让项目能马上开工，但工资特别高，项目一结束你恨不得立刻把人请走。M 就是这个高工资。

这里有个关键：M 必须“足够大”，但到底多大算足够？这在实际数值计算里是个很微妙的问题，后面我单独讲。

1.3 一个能当测试基准的算例

为了把两种方案都跑通，我选了一个非常典型的算例：

min z = 2x1 + 3x2

约束条件：

• x1 + x2 >= 2
• x1 + 2x2 >= 3
• x1, x2 >= 0

为什么说它典型？因为两条约束全是“大于等于”，直接加松弛变量的话，松弛变量前面的系数是 -1，没法当正基变量用，所以必然要引入人工变量。这个算例的最优解可以手算验证：x1 = 1，x2 = 1，最优目标值 z = 5。拿它当测试基准，无论手写代码还是 scipy 求解，结果对不对一眼就能看出来。

先把不等式改成等式。第一条约束加剩余变量 s1，第二条加剩余变量 s2：

• x1 + x2 - s1 = 2
• x1 + 2x2 - s2 = 3

于是标准型里的系数矩阵 A 就是 [[1, 1, -1, 0], [1, 2, 0, -1]]，目标系数 c 是 [2, 3, 0, 0]，右端项 b 是 [2, 3]。注意，这里的 s1、s2 是非负变量，但它们的列是 -1，所以初始基还是得靠人工变量。

2. 手写大M法：从单纯形表到Python代码

2.1 数据结构设计：别用二维表，用基变量集合

教材里的单纯形法通常画一张单纯形表，一行一行更新。Python实现时，我建议不要维护整张表，而是维护基变量索引集合和非基变量索引集合。原因是每次迭代只需要访问基矩阵、目标系数和右端项，不需要把整张表搬来搬去，逻辑更清晰，也更容易排查问题。

我的核心数据结构是：A_full 是“原系数矩阵 + 单位矩阵”的增广矩阵，c_full 是“原目标系数 + M”的增广系数。basis 存基变量在增广矩阵里的列索引，nonbasis 存非基变量索引。初始状态，每个约束行补一个人工变量，所以 basis 初始化为最后一组单位列，nonbasis 则是其余变量。

这样设计的好处是，代码里的每一次换基，本质上就是 basis 里换一个元素，同时把被换出的变量丢进 nonbasis。整个求解过程就是在这两个集合之间做“交换”，非常直观。

2.2 三个关键步骤：检验数、最小比值、换基

单纯形法的每次迭代可以压缩成三个步骤。

第一步，计算非基变量的检验数。检验数表示“非基变量每增加一个单位，目标函数会变化多少”。在最小化问题里，如果所有非基变量的检验数都大于等于零，说明再增加任何一个非基变量都无法让目标更小，当前解就是最优解。如果有负的检验数，就选最负的那个作为进基变量。

第二步，确定离基变量。进基变量确定后，它要从当前值 0 开始增大，但其他约束不能突破，所以要看这个新方向上的每个正分量的比值，取最小的那个作为限制。这个规则叫最小比值原则，数学上就是让基变量中第一个被“压到零”的变量出基。

第三步，换基。把进基变量加入 basis，把离基变量从 basis 里移除，完成一轮旋转。理论上新基对应的目标值会比上一轮更小，只要没有退化和循环，算法就会收敛。

三个步骤里，最容易出错的是检验数计算。很多人会记错符号，尤其是用最小化形式写的时候。这里统一约定：检验数 r_j = c_j - c_B^T B^{-1} A_j，最小化问题选 r_j 最小的负值进基，所有 r_j >= 0 时停止。这个约定贯穿我后面的所有代码。

2.3 完整代码：一个可以直接跑的大M法类

下面是我写的一个精简但完整的实现：

import numpy as np class BigMSimplex: def __init__(self, c, A, b, M=1e6): self.c = np.array(c, dtype=float) self.A = np.array(A, dtype=float) self.b = np.array(b, dtype=float) self.M = float(M) self.m, self.n = self.A.shape # 每一行补一个人工变量，形成 [A | I]，确保初始基是单位矩阵 self.A_full = np.hstack([self.A, np.eye(self.m)]) self.c_full = np.concatenate([self.c, np.full(self.m, self.M)]) # 初始基变量：所有人工变量 self.basis = list(range(self.n, self.n + self.m)) self.nonbasis = list(range(self.n)) self.eps = 1e-9 def solve(self, max_iter=1000): for _ in range(max_iter): B = self.A_full[:, self.basis] cB = self.c_full[self.basis] try: B_inv = np.linalg.inv(B) except np.linalg.LinAlgError: return None, "基矩阵奇异，请检查约束是否线性相关" xB = B_inv @ self.b if np.any(xB < -self.eps): return None, "基变量出现负值，数值异常" # 检验数 = c_j - c_B^T B^{-1} A_j y = cB @ B_inv entering = None min_rc = -self.eps for j in self.nonbasis: rc = self.c_full[j] - y @ self.A_full[:, j] if rc < min_rc: min_rc = rc entering = j if entering is None: return self._build_result(xB) # 最小比值原则确定离基变量 d = B_inv @ self.A_full[:, entering] leaving = None theta = np.inf for i, di in enumerate(d): if di > self.eps: ratio = xB[i] / di if ratio < theta: theta = ratio leaving = self.basis[i] if leaving is None: return None, "无界解：进基方向不存在正分量" # 换基：离基变量出，进基变量入 pos = self.basis.index(leaving) self.basis[pos] = entering self.nonbasis.remove(entering) self.nonbasis.append(leaving) return None, "达到最大迭代次数，疑似退化循环" def _build_result(self, xB): x_full = np.zeros(self.n + self.m) for idx, val in zip(self.basis, xB): x_full[idx] = val # 人工变量仍大于阈值，说明原问题无可行解 art_vars = range(self.n, self.n + self.m) for idx in art_vars: if x_full[idx] > 1e-6: return None, "无可行解：人工变量未出基" return { "x": x_full[:self.n], "z": self.c_full @ x_full, "basis": list(self.basis), }

使用这个类有一个前置动作：你必须先把所有约束手动转换成等式标准型，并且保证 b 全为正。比如上面的算例，我已经把“大于等于”约束减掉了剩余变量，所以传入的 A 是 [[1, 1, -1, 0], [1, 2, 0, -1]]。如果 b 里有负数，约束两边先乘 -1 再处理，否则初始基变量会出现负值，整个算法就废了。

2.4 跑一下算例：手算的每一步都有代码对应

把算例喂进去：

model = BigMSimplex( c=[2, 3, 0, 0], A=[[1, 1, -1, 0], [1, 2, 0, -1]], b=[2, 3], M=1e6 ) res, msg = model.solve() print(res) print(msg)

输出结果是：

{'x': array([1., 1., 0., 0.]), 'z': 5.0, 'basis': [0, 1]}

x 的前两个值就是 x1 = 1、x2 = 1，目标值 z = 5，跟手算完全一致。basis 是 [0, 1]，说明最优基里两个原变量都在，人工变量已经全部被赶出去了，这意味着解落在原问题的可行域内。

建议你自己手算一遍这个算例，你会发现迭代过程中目标值的变化轨迹：初始目标值高得离谱，因为人工变量在基里背着巨额 M；然后每一步人工变量逐渐被挤出基，目标值快速下降；最后人工变量全部出基，得到干净的可行解。这个过程只要跑通一遍，大M法就算真正入门了。

3. scipy.optimize.linprog：三行代码搞定，但要小心符号

3.1 linprog的输入约定与约束转换

手写代码适合理解原理，真到了工程里，我只推荐直接用现成库。scipy.optimize.linprog 是 SciPy 里求解线性规划的标准接口，用法极简，但符号约定必须先搞清楚。

linprog 默认求解的是最小化问题，约束分两部分：不等式约束统一写成 A_ub @ x <= b_ub，等式约束用 A_eq @ x = b_eq。变量的下界上界由 bounds 参数控制，默认是变量非负，但我建议永远显式写清楚。

“大于等于”约束怎么处理？两边乘以 -1，变成“小于等于”。比如 x1 + x2 >= 2，就写成 -x1 - x2 <= -2。很多人第一次用都会在这里翻车，忘了转符号，然后对着一个迷惑的解百思不得其解。

3.2 用同一算例验证

上面的算例用 linprog 求解，代码长这样：

from scipy.optimize import linprog c = [2, 3] # x1 + x2 >= 2 -> -x1 - x2 <= -2 # x1 + 2x2 >= 3 -> -x1 - 2x2 <= -3 A_ub = [[-1, -1], [-1, -2]] b_ub = [-2, -3] bounds = [(0, None), (0, None)] res = linprog(c, A_ub=A_ub, b_ub=b_ub, bounds=bounds, method='highs') print(res.x) print(res.fun) print(res.success) print(res.message)

输出：

[1. 1.] 5.0 True Optimization terminated successfully.

结果一致：x1 = 1，x2 = 1，最优值 5.0。这说明手写代码和 scipy 在高层面是等价的——都求出了同一个最优解。

这里有个体验上的细节：手写大M法的类需要你先手工转换等式标准型，再传增广矩阵；而 linprog 直接吃原始不等式约束，省掉了不少体力活。scipy 内部会自动处理初始基的问题，它本质上用的不是大M法的大数惩罚策略，而是更稳健的两阶段法，所以不会遇到大M法那种 M 值调参的麻烦。

3.3 method参数怎么选

linprog 的 method 参数在新版 SciPy 里有几个选项：highs、highs-ds、highs-ipm、interior-point、simplex。如果不上点心，很容易被这一堆名字绕晕。

我的建议是：直接用默认的 method='highs'。这个参数对应 HiGHS 求解器，是当前版本默认选项，综合表现最稳。它在内部会根据问题规模自动选择单纯形法、对偶单纯形法或者内点法，而且支持稀疏矩阵，处理几千个变量的问题毫无压力。

旧版的 method='simplex' 是两阶段单纯形法的实现，适合小规模教学验证，但在新版文档里已经有废弃倾向，不建议在新的代码里依赖它。如果你遇到大规模稀疏问题，可以试试 highs-ipm，那是内点法路线，某些场景下比单纯形法快。

还有一个调试小技巧：linprog 支持 options 参数，比如 options={'disp': True} 可以输出迭代日志。配合 callback 函数还能拿到每一步的迭代点，这在排查“为什么结果不对”的时候非常有用。

4. 两种方案深度对比：教学、调试、上生产该怎么选

4.1 手写大M法 vs scipy：一张表看明白

我整理了一张对比表，尽量把两种方案的差异说清楚：

这个表格基本反映我的真实结论：手写大M法和 scipy 不是同一个定位的东西。前者是教学工具、学习工具，后者是工业工具。你要是拿手写大M法去求解一个 5000 变量的生产问题，大概率跑到天荒地老，而且数值上还不见得收敛。

4.2 我实际使用的选型路线

我在实际项目里的经验是分场景走：

第一类场景，教学和作业。目标是理解单纯形法和大M法的每一步，那必须自己手写，而且要刻意不用 scipy，把旋转过程打印出来，一行一行看懂。

第二类场景，原型验证。需求是快速知道一个线性规划问题的最优解长什么样，比如做一个资源分配模型的demo，直接用 linprog，五分钟出结果，重点放在业务建模而不是算法细节。

第三类场景，生产环境。如果问题规模很大或者约束性质特殊，我会优先考虑 linprog 的 highs，它已经能覆盖绝大多数需求；如果问题里还有整数变量，那就是混合整数规划，linprog 就不够用了，得换专门的求解器。

一个容易被忽略的点是：手写大M法的代码虽然简单，但它隐含了一个假设——约束矩阵的行必须线性无关。如果约束互相矛盾或者冗余，基矩阵求逆会直接失败，而这类问题在工程数据里非常常见。linprog 内部的预处理会把冗余约束、矛盾约束先筛一遍，这层保护手写代码是没有的。

5. 常见问题与排查技巧实录

5.1 结果为nan或浮点异常

手写大M法最常见的问题就是结果出现 nan。原因通常出在 M 太大上。当 M 取 1e10 甚至更大时，目标函数系数里出现 1e10 和 3 这样的量级差异，单纯形表里的检验数会变成 1e10 减去某个数，浮点精度根本接不住这种动态范围。

我调试时遇到过检验数算出来是 -9999999.999 和 0.000001 这种数字，看着好像没问题，但一次旋转之后，基变量值全乱了。解决思路是把 M 降到 1e4 到 1e6 这个区间，够用就行，不要一味求大。

如果你发现结果是 nan，优先检查三件事：M 是否过大、b 是否有负数、约束矩阵是否线性相关。这三点占了九成以上的问题。

5.2 M值过大或过小的连锁反应

M 值太小会怎样？假设某个问题里人工变量带来的惩罚不够大，算法可能宁可让一个人工变量留在基里，也不愿意付出增加某个变量的代价，最后给出的“最优解”里人工变量是个正数，原问题根本没有可行解。

M 值太大又会怎样？浮点数在极端量级下会丢失精度，大M法的检验数本质上是一个“大数减小数”，一旦两个数都大得离谱，差值就不准了，迭代方向可能选错，甚至直接不收敛。

我在实践中找 M 的经验是：先看目标函数系数和约束系数的量级，取一个比它们大三到四个数量级的数，但又不要超过 1e7。比如系数都在 10 以内，M 取 1e5 到 1e6 就很合适。做完之后一定要检查人工变量是否真的是零，别只依赖代码里那个 1e-6 的阈值，最稳的办法是代入原约束手工验证。

5.3 最优解不唯一怎么判断

单纯形法告诉我“所有非基变量检验数都大于等于零”就停下来了，但如果你算出来的检验数里有零，那就说明可能存在多个最优解。比如目标函数刚好平行于某条约束边界，整条边界上的点都是最优。

这个现象在手写代码里很容易被忽略，因为代码只要找到一个最优解就返回了。实际业务里，最优解不唯一有时候是好事，表示方案有弹性；有时候是坏事，因为你不确定拿到的解是不是最符合业务偏好的那个。想看所有最优解，可以给一个零检验数的非基变量尝试进基，比如加一个很小的扰动，看目标值是否不变。linprog 本身不会告诉你解的多样性，你只能手动分析。

5.4 退化与循环：不是大M的专利

退化是指换基之后目标函数值没有下降，这在实践中比想象中常见。如果连续多轮目标值不变，算法可能在几个基之间反复横跳，这就是循环。经典教材里会提到 Bland 规则来防止循环，简单理解就是“进基和离基都选最小下标”，代价是收敛可能变慢。

我的手写代码里没有实现 Bland 规则，而是加了最大迭代次数上限，就是为了防止陷入死循环。出现这种情况时，我的处理方式有两个：要么改用 Bland 规则重跑，要么加一个极小的随机扰动把退化打破。生产环境里更彻底的办法是直接用 linprog，它内部有完善的防循环机制。

5.5 无界解与无可行解的识别

无界解的特征比较明显：选中的进基变量方向上的系数全是非正的，也就是说无论这个变量增到多大，都没有任何约束拦住它，目标函数可以一路降到负无穷。我的代码里对应的是“leaving is None”这个分支，这时候就要检查模型是不是少了约束。

无可行解的情况更隐蔽。大M法里，如果迭代结束但人工变量还在基里且值大于零，说明原问题根本找不到满足全部约束的解。我们代码里 _build_result 最后一段就是干这个的。你可以把算例改成一个矛盾的约束，比如同时要求 x >= 5 和 x <= 3，试试代码会不会正确报出“无可行解”，这个实验对理解人工变量的意义很有帮助。

个人在实际操作中的体会是：手写一遍大M法，最大的收获不是能默写代码，而是理解了为什么工程工具里几乎没有用纯大数惩罚的。数值计算最怕量级失衡，两阶段法之所以成为主流，就是因为它绕开了大M数，把问题拆成了两个干净的小问题。如果你是初学，我非常建议把大M法和两阶段法对照着写一遍；如果你只是要用线性规划，那 linprog 的 highs 就是当前最省心的选择，把符号和约束转换搞对，剩下的交给库就好。