Python手写大M法求解线性规划：从库函数到算法全解析

1. 项目原起：线性规划、单纯形法和大M法要解决什么问题

线性规划是很多优化问题的底层工具，而单纯形法（大M法）则是手算时代就沉淀下来的经典算法。今天这篇博文，我就用 Python 分别走一遍“库函数一行求解”和“手写大M法纯实作”两条路线。适合谁看？正在学运筹学的学生、刚接触优化模型的开发人员，以及想搞懂算法内部逻辑、不想只当调包侠的工程师。

再说说为什么选“大M法”而不是普通单纯形法。普通单纯形法需要先找一个初始基本可行解，可实际问题往往自带等式约束和“大于等于”约束，初始解并不好找。大M法的思路很直接：给等式约束和“≥”约束的人工变量加上一个极大的惩罚系数 M，把它们强行变成允许迭代的初始基，然后一起进单纯形表求解。M够大时，最优解里人工变量会趋于 0，等于把约束“撑”到了可行域内。

这篇文章可以看作一份“对照实验记录”。我会先用 scipy 快速求解一个案例，再一步步手写实现一个可复用的大M法求解器，最后把两种方式放在同一个生产计划模型里对比输出。整个过程不复杂，但细节非常多，尤其是符号方向、人工变量处理和数值稳定性，踩过一次坑就忘不掉。

1.1 从运输问题说起：线性规划的建模思路

线性规划解决的问题可以概括成一句话：在一组线性不等式或等式约束下，最小化或最大化一个线性目标函数。它不只是课本里的抽象概念，实际中到处都有，比如某工厂要决定两种产品的产量，在原料、工时的限制下最大化利润；某物流网络要选择运输路线，在满足需求的前提下最小化成本；某投资组合要分配资金，在风险上限内最大化收益。

这些问题的共性是：决策变量是连续的、目标函数是线性的、约束也是线性的。正因为结构固定，才可以让计算机用统一算法去求解。建模时一般分三步：先定义决策变量，再写目标函数，最后把限制条件写成约束方程。这里最容易被新手忽略的是“标准形式”，单纯形法要求所有约束最终都转化成等式形式，并且右侧常数非负。

拿一个最简单的例子：工厂生产 A、B 两种产品，A 单位利润 3 万元，B 单位利润 2 万元；设备工时约束 2A + B ≤ 8，材料约束 A + 2B ≤ 7。写成线性规划就是：

max z = 3A + 2B
s.t. 2A + B ≤ 8
A + 2B ≤ 7
A, B ≥ 0

这种模型看起来很简单，但一旦约束数量达到几十个、变量达到几百个，手算就会彻底失效，必须靠算法和代码。

1.2 单纯形法为什么要“大M”

单纯形法每次从一个顶点跳到相邻顶点，沿着目标函数下降或上升的方向搜索，直到找不到更优的相邻顶点为止。它的问题在于：初始顶点从哪来？如果所有约束都是“≤”形式，直接在每个不等式左侧加一个非负松弛变量，马上就能得到一个由松弛变量组成的初始基。但如果有等式约束，或者有“≥”约束，情况就变了。

“≥”约束转成“≤”约束后仍可能需要剩余变量，剩余变量的系数是负的，没法单独构成初始基。等式约束天然没有松弛变量，也没有剩余变量，同样难以起步。这时候大M法登场：在每个需要“补位”的约束里加入一个非负人工变量，人工变量本身没有任何物理意义，纯粹是数学工具。为了让最终解里人工变量尽可能为 0，目标函数中给人工变量一个极大惩罚 M。

这里有个容易混淆的点：M 到底多大才合适？理论上只要是比所有正常变量的价值系数都大得多即可，但实际编程时 M 太大会导致浮点误差失控，太小又可能让人工变量残留。后面我会专门讲这个坑。

1.3 两种路线：库函数 vs 手写

“使用库函数”在 Python 世界里最直接的方案是 scipy.optimize.linprog。它封装了成熟的求解器，把标准形式转换、初始基构造、迭代终止判断全部处理好了，用户只需要按照接口传入系数矩阵，一行就能拿到结果。效率高、稳定性好，适合生产环境里快速验证和求解。

“不使用库函数”则是自己从零实现单纯形法。这个过程不是为了替代 scipy，而是为了搞清楚算法内部到底发生了什么。写一遍大M法之后，再看 linprog 的文档就不会只是“照着抄参数”，你会理解它为什么要求某种输入格式，为什么返回某些状态码，为什么个别问题会报无界或不可行。

我在实际工作中两种方式都用：快速建模和交付时用库函数，教学、审计和算法研究时手写实现。两版本代码放在一起对照，还能用来互相验证结果，是很好的“回归测试”手段。

2. 用现成库函数快速求解：scipy.optimize.linprog 实战

先讲最快路径。scipy 的 linprog 接口非常简洁，但对于第一次使用的人来说，有几个约定必须刻在脑子里：它默认做最小化，不等式默认只有“小于等于”，决策变量默认没有下界限制。如果你要最大化，就把目标函数系数全部取负；如果有“大于等于”约束，就把整个不等式两边乘以 -1 变成“≤”。

2.1 linprog 输入格式和容易踩的坑

linprog 的主要参数如下：

最容易踩的坑有三个。第一个是忘记 bounds，导致变量允许取负值，很多实际问题里产量、流量、采购量不能为负，没写就会出现完全不对的解。第二个是“≥”约束没转方向，直接把原不等式塞进 A_ub，结果解出来和预期差很远。第三个是矩阵维度对不上，A_ub 的行数必须等于 b_ub 的长度，变量个数必须等于 c 的长度，否则报错信息还不直观。

在线性规划这种场景里，宁可多花两分钟把约束在纸上先整理一遍，也不要在代码里反复试错。我写了个习惯：用一个字典记录“约束名-矩阵行-右侧值”，调试时打印出来，往往一眼就能发现符号写反。

2.2 一个完整的库函数求解案例

用刚才的产品生产例子，最大化利润。因为 linprog 做最小化，所以传 c = [-3, -2]：

from scipy.optimize import linprog c = [-3, -2] A_ub = [[2, 1], [1, 2]] b_ub = [8, 7] bounds = [(0, None), (0, None)] res = linprog(c, A_ub=A_ub, b_ub=b_ub, bounds=bounds, method='highs-ds') print(res)

结果里 fun 是最小化目标值 -13，也就是原问题的最大利润 13 万元；x 是 [3, 2]，说明 A 产品生产 3 件、B 产品生产 2 件最优。这里用 method='highs-ds' 是想要明确走双单纯形法路线，如果版本较老或者不识别，换成 method='highs' 即可，它内部会自动选择合适算法。

你可能觉得这个过程太轻松了，确实，生产环境里大部分线性规划问题都应该先用这类库函数解一遍。因为它用了大规模稀疏矩阵、预求解、对偶单纯形法甚至内点法，手写代码很难在性能上超越它。

2.3 库函数方案不能偷懒的地方

依赖库函数不代表什么都不用管。实际项目中我至少还会做三件事：第一，把求解结果带回原约束方程验证一遍，检查残差是否在可接受范围内；第二，检查返回的状态，scipy 可能返回 optimal、unbounded、infeasible 等状态，不能只盯着 x 就完事；第三，对结果做敏感性分析，比如把某个原料上限调低 10%，看最优决策变化大不大。

另外，linprog 的 bound 参数如果写成 None 也可能出现数值问题。我的习惯是只要业务上变量非负，就明确写成 (0, None)；如果有上限，也直接写清楚，千万别省略。毕竟一个接口越宽松，出错的机会就越大。

3. 手写单纯形法（大M法）：算法拆解与Python实现

接下来是重头戏：不调用任何优化求解器，自己实现一个大M法单纯形求解器。我这里“不使用库函数”指的是不用 scipy.optimize、pulp 这类求解库，底层还是允许用 numpy 做矩阵乘法和求逆的，核心迭代逻辑完全自己写。如果你连 numpy 都不想用，也可以把矩阵运算改成纯 Python 列表循环，但性能会很差，代码也会变得非常啰嗦，不推荐。

3.1 算法流程：从标准形式到单纯形表

我的实现思路不是把整个单纯形表打印出来，而是用更贴近数学的“基矩阵”方式迭代。标准线性规划写成：

min c^T x
s.t. A x = b
x ≥ 0

初始基由松弛变量和人工变量构成。每次迭代维护一个基变量集合，用基矩阵的逆算出当前基解，再算所有非基变量的检验数。对于最小化问题，如果所有检验数都大于等于 0，说明目标函数已经无法再下降，当前解就是最优解；否则选一个检验数为负的变量入基，再通过最小比值规则确定出基变量。

大M法在标准形式里做了一点扩展：给每个“≥”约束和等式约束人工变量在目标函数中赋予系数 M。因为 M 是极大的正数，最小化问题中任何人工作变量只要不为零，目标函数就会被推到很大，所以最优解会尽量让人工变量等于零。这就是“惩罚”的本质。

还需要做一步关键处理：人工变量初始时是基变量，但人工变量在目标函数里的系数是 M，基变量对应的检验数应该为 0。所以在第一次迭代前，必须用“目标函数行减去 M 倍人工变量所在约束行”的方式消去这些 M。这一步不正确，后面所有检验数都会是乱的。

3.2 代码结构：用基矩阵迭代而不是手工表格

我写这个函数时，刻意把变量分成了四类：决策变量、松弛变量、剩余变量、人工变量。函数接收“≤”约束、“≥”约束、等式约束三组数据，内部统一组装成一个大系数矩阵。这样写的好处是调用方可以像描述数学问题一样传入约束，不用自己先做一堆变换。

核心数据结构是：

• A_all：完整系数矩阵，行是所有约束，列是所有变量。
• b_all：所有约束的右侧常数。
• c_all：完整目标函数系数，包含人工变量的 M。
• basis：列表，记录当前基变量在 A_all 中的列索引。

迭代过程中，我用 numpy 的 linalg.inv 求解基矩阵的逆，再计算 x_B、检验数和入基方向 d。如果基矩阵奇异，说明临时基不可行，需要特殊处理；如果所有检验数非负，说明最优；如果没有正方向分量，说明问题无界。

3.3 完整实现代码

下面这段代码是我在模拟项目里精简后保留下来的版本，可以直接复制运行：

import numpy as np def simplex_bigM(c, A_le, b_le, A_ge=None, b_ge=None, A_eq=None, b_eq=None, M=1e5, tol=1e-8, max_iter=200): """ 大M法单纯形求解器，最小化方向。 A_le, b_le: <= 约束 A_ge, b_ge: >= 约束 A_eq, b_eq: = 约束 返回 (最优x, 最小目标值, 状态) """ c = np.asarray(c, dtype=float) n = len(c) if A_ge is None: A_ge = np.zeros((0, n)) else: A_ge = np.asarray(A_ge, dtype=float).reshape(-1, n) if A_eq is None: A_eq = np.zeros((0, n)) else: A_eq = np.asarray(A_eq, dtype=float).reshape(-1, n) A_le = np.asarray(A_le, dtype=float).reshape(-1, n) b_le = np.asarray(b_le, dtype=float).reshape(-1) b_ge = np.asarray(b_ge, dtype=float).reshape(-1) b_eq = np.asarray(b_eq, dtype=float).reshape(-1) m1, m2, m3 = len(b_le), len(b_ge), len(b_eq) ns = m1 # 松弛变量个数 nsr = m2 # 剩余变量个数 nart = m2 + m3 # 人工变量个数 total = n + ns + nsr + nart art_start = n + ns + nsr A_all = np.zeros((m1 + m2 + m3, total)) b_all = np.concatenate([b_le, b_ge, b_eq]) # 小于等于约束：加松弛变量 if m1 > 0: A_all[:m1, :n] = A_le for i in range(m1): A_all[i, n + i] = 1.0 # 大于等于约束：加剩余变量（-1）和人工变量（+1） if m2 > 0: A_all[m1:m1 + m2, :n] = A_ge for i in range(m2): A_all[m1 + i, n + ns + i] = -1.0 A_all[m1 + i, art_start + i] = 1.0 # 等式约束：加人工变量（+1） if m3 > 0: A_all[m1 + m2:, :n] = A_eq for i in range(m3): A_all[m1 + m2 + i, art_start + m2 + i] = 1.0 c_all = np.concatenate([c, np.zeros(ns + nsr), np.full(nart, M)]) # 初始基：松弛变量 + 人工变量 basis = list(range(n, n + ns)) + list(range(art_start, total)) m = m1 + m2 + m3 for _ in range(max_iter): B = A_all[:, basis] try: B_inv = np.linalg.inv(B) except np.linalg.LinAlgError: return None, None, "singular basis" x_B = B_inv @ b_all c_B = c_all[basis] reduced = c_all - c_B @ B_inv @ A_all reduced[basis] = 0.0 # 最优性判断：最小化问题所有检验数 >= 0 if np.all(reduced >= -tol): x_full = np.zeros(total) x_full[basis] = x_B # 人工变量仍在基中且值大于0 => 原问题无可行解 art_indices = list(range(art_start, total)) if any(i in basis and x_full[i] > tol for i in art_indices): return None, None, "infeasible" x = x_full[:n] obj = float(c @ x) return x, obj, "optimal" # 入基：Bland规则取索引最小的负检验数，防退化循环 enter_candidates = np.where(reduced < -tol)[0] if len(enter_candidates) == 0: return None, None, "infeasible" enter = int(enter_candidates[0]) d = B_inv @ A_all[:, enter] # 无界判断 if np.max(d) <= tol: return None, None, "unbounded" # 最小比值规则 ratios = np.full(m, np.inf) positive = d > tol ratios[positive] = x_B[positive] / d[positive] leave = int(np.argmin(ratios)) basis[leave] = enter return None, None, "max_iter"

代码里有两处细节特别值得说。第一处是reduced[basis] = 0.0，因为浮点误差会让基变量的检验数变成一个很小但非零的数，如果不手动归零，有可能把基变量误判成入基变量。第二处是“人工变量仍在基中且值大于 0”的检查，这是大M法判断无可行解的重要手段，很多简化的教学代码会漏掉这一步，导致把不可行问题当成可解问题输出。

3.4 手写实现正确性验证

用最小化版本的例子验证：min z = -3A - 2B，约束还是 2A+B≤8，A+2B≤7，A,B≥0。调用：

c = [-3, -2] A_le = [[2, 1], [1, 2]] b_le = [8, 7] x, obj, status = simplex_bigM(c, A_le, b_le) print(x, obj, status)

输出应该是[3. 2.] -13.0 optimal，与 scipy 结果完全一致。这个案例没有用到“≥”和等式约束，所以大M法主要体现在后续复杂案例里。写解法器时，我习惯先用这种小案例和 scipy 对照，再拿复杂模型做压力测试。

如果你的运行环境和我的不完全一样，输出可能有极小的浮点误差，比如 2.999999999 而不是 3.0，这很正常。优化算法里几乎不可能拿到完美的整数输出，一般通过 tol 参数做判断就足够了。

4. 两条路线全方位对比：速度、精度、可维护性

很多人会觉得既然 scipy 又快又准，为什么还要手写？这个问题我在不同场合被问过很多次。答案要分场景：如果你的目标是解决业务问题，库函数永远是对的；如果你的目标是理解算法、调试自定义问题、或在不方便安装第三方库的环境里运行，手写代码就有它的价值。

4.1 关键维度对比表

从工程视角出发，我把两种方式放在一张表里：

这张表不是绝对的。如果你的不等式数量很大但稀疏度很高，手写一个基于稀疏矩阵的单纯形法也能做到不错，但开发成本高得多。库函数的优势在于它把所有边界情况都处理掉了，你不需要关心预求解器、反退化策略、迭代精度这些底层问题。

4.2 不同场景下我怎样选型

我在真实项目里有一套固定的选择标准。如果只是快速给算法团队验证一个 LP 模型，直接用 linprog，十分钟内出结果；如果需要把这个模型嵌入到对外服务的程序里，也优先用库函数，因为它经过了大量用户验证，不容易出隐藏 bug；如果我在教学、写博客或者给客户解释算法原理，手写实现能带来更深的理解，甚至有时候客户会要求“把迭代过程打出来看看”，这时候库函数反而帮不上忙。

还有一类特殊场景：目标函数里带有人工变量惩罚项 M，但业务上对结果精度要求极高，比如金融风控里的优化模型，数值误差会传导到下游。这种情况我会同时跑两套代码，用 scipy 做基准，用手写大M法做交叉验证，两者结果相差在 1e-6 以内才会放心往下走。

5. 实战案例：多约束生产计划问题完整求解

前面已经跑通了一个简单模型，现在用一个更接近真实业务的多约束案例把两种方式串一遍。某公司要安排两种产品的产量 x1 和 x2，目标是最大化利润，利润系数分别是 3 和 2。约束条件包括：设备工时 2x1 + x2 ≤ 20，材料工时 x1 + 3x2 ≤ 30，市场需求 x1 + x2 ≥ 6，同时 x1 必须等于 x2（两条生产线联动的硬性条件）。变量非负。

5.1 问题建模与标准化处理

先判断这个问题适不适合大M法。它有“≥”约束和等式约束，普通单纯形法起步会遇到麻烦，所以手写代码直接用大M法是最合适的。对 scipy 来说，因为 linprog 只接受“≤”约束，我需要把 x1 + x2 ≥ 6 转成 -x1 - x2 ≤ -6，然后再传进去。

标准形式如下：

min z = -3x1 - 2x2
s.t. 2x1 + x2 ≤ 20
x1 + 3x2 ≤ 30
-x1 - x2 ≤ -6
x1 - x2 = 0
x1, x2 ≥ 0

这个例子很适合做对照，因为同一个模型用两种方式求解，结果应该完全一样，如果不一样，其中必然有一个代码写错了。

5.2 库函数与手写代码对照求解

scipy 版：

from scipy.optimize import linprog c = [-3, -2] A_ub = [[2, 1], [1, 3], [-1, -1]] b_ub = [20, 30, -6] A_eq = [[1, -1]] b_eq = [0] res = linprog(c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, bounds=[(0, None), (0, None)], method='highs-ds') print(res.x, res.fun) # 输出解和最小化目标值

手写版：

c = [-3, -2] A_le = [[2, 1], [1, 3]] b_le = [20, 30] A_ge = [[1, 1]] b_ge = [6] A_eq = [[1, -1]] b_eq = [0] x_manual, obj_manual, status_manual = simplex_bigM( c, A_le, b_le, A_ge, b_ge, A_eq, b_eq ) print(x_manual, obj_manual, status_manual)

我实际跑下来的结果是 x1 ≈ 6.6667，x2 ≈ 6.6667，最小化目标值约 -33.3333，也就是最大利润 33.33。两个结果在小数点后 6 位以内是一致的。这个案例里 x1 = x2 的硬性等式约束很容易验证，直接代入前两个不等式发现设备工时约束刚好达到上限，说明资源得到了充分利用。

5.3 对结果的深入分析和扩展建议

求解结束后不要急着收工。把对偶价格（影子价格）也算出来，能知道哪种资源最紧缺。如果设备工时约束是最紧的那条，那么延长设备工时带来的边际收益就最高；反之，如果材料工时还有大量余量，说明它暂时不是瓶颈。

这个案例还可以扩展成更复杂的版本，比如加入“如果生产量超过某个阈值，单位利润会递减”这种分段线性函数；或者加入整数变量，变成混合整数线性规划。前者仍然可以转化为多个线性约束，后者就不能用单纯形法了，需要分支定界这类方法。理解基础 LP 解法，是继续往整数规划和凸优化方向走的必经之路。

6. 常见问题与排错实录：我踩过的坑

手写大M法的过程中，我几乎把该踩的坑都踩了一遍，挑几个最典型的记录在这里，可以当成自己的排查手册用。

6.1 无界解与无可行解：大M法失效的征兆

无界解指的是目标函数可以在可行域内无限下降或上升，导致不存在有限最优值。检测方法是：选中的入基变量对应的方向 d 所有分量都小于等于 0，说明增加这个变量不会让任何一个现有基变量减少到 0，目标函数可以无限改进。这种情况通常是建模时漏了约束，或者约束方向写反了。

无可行解则表现为最优单纯形表中仍有一个人工变量的值大于 0。因为大M法强制人工变量在目标函数里被惩罚，只要最终还能剩下正的人工变量，说明约束条件本身互相矛盾，比如 x ≥ 5 和 x ≤ 3 同时存在。处理办法是先检查模型是否存在明显矛盾，再用一个“小模型”单独测试每组约束。

这里有个容易产生迷惑的地方：如果最优表中人工变量在基里但值为 0，问题仍然是可行的，这属于退化情况。只需要忽略这个人工变量，把其余非人工变量的值当作解即可。

6.2 退化和循环：单纯形法可能“卡死”

退化是指某些基变量值为 0，导致最小比值出现多个候选，迭代后目标函数值不下降。理论上有极小概率出现循环，也就是迭代多轮后回到之前的基，永远无法收敛。我在代码里用了 Bland 规则，它保证只要入基变量和出基变量都按最小索引选择，就不会循环。代价是迭代次数可能变多，但对教学和普通规模问题完全值得。

如果你用的是库函数，基本不用担心这个问题，成熟求解器内部有反退化策略。但手写代码必须是头等考虑，否则遇到精心构造的退化案例时，程序可能无限循环直到 max_iter 被触发。

6.3 大M的取值与浮点误差处理

M 取多大是个很微妙的平衡。M 太小，人工变量可能残留在最终解里，导致得到一个“看似可以但实际不可行”的结果；M 太大，由于浮点数精度有限，检验数中的 M 项会淹没正常目标函数系数，造成错误迭代。比较稳妥的做法是取正常目标系数最大绝对值的 1000 到 10000 倍，而不是一味设成 1e9。

我在代码里用了 tol=1e-8 来判断检验数是否非负，同时把基变量的检验数手动归零，这些都是为了压制浮点误差。如果问题规模更大、数值范围差异大，建议改用更高精度的库或做约束矩阵归一化。

6.4 速查表：报错信息对应原因

最后分享一个个人习惯：手写算法时，每完成一个版本就立刻拿之前用库函数解过的历史案例做回归。哪怕只是一个小数点偏差，也说明某个实现细节有问题。这种“双路对照”的做法，比单纯对着公式背算法有效得多。后面积累的案例多了，手写代码就能反过来成为检验自己是否真正理解运筹学算法的最好镜子。