企业官网建设流程全解析

前阵子一个排产优化的需求找过来：变量不到十个，约束也不复杂，但既有不等式又有等式，整个模型的形态很“不规则”。朋友问我说，到底该直接用现成的库函数一行求解，还是要自己写单纯形法？我的回答是：两条路都得会。这也是今天这篇东西要系统讲清楚的题——线性规划的 Python 实现，重点放在单纯形法里的大M法，分别说透不依赖任何库函数的手写版本，和调用标准库函数的版本，各自的原理、代码、优缺点和适用边界。适合正在学运筹学、准备面试手撕算法、或者第一次在工程里用线性规划建模的读者。

1. 两条实现路线背后真正要解决的问题

1.1 库函数解决的是“算得快、算得稳”

工程里面对一个线性规划模型，我的第一反应永远是先看有没有成熟求解器。Python 生态里最常用的就是 SciPy 里的linprog，底层有 HiGHS 这样的高性能求解器支撑，纯 C++ 实现，数值稳定性远好于自己造的轮子。几十个变量、几百个约束这种规模，库函数基本是毫秒级出结果，而且对退化、数值病态、约束冗余都有统一处理。

但库函数有个前提：你必须先把模型整理成它认得的格式。约定是极小化目标，约束拆成小于等于、等式、变量边界三块。很多新手拿一个极大化、带大于等于约束的问题直接往linprog里塞，返回一个错误状态，就以为库不靠谱，其实是建模入口错了。所以库函数不是帮你不懂标准形，而是帮你省略从标准形到单纯形表的过程。

1.2 手写实现解决的是“看得懂、控得住”

手写单纯形法的价值不在生产环境，在于它逼你把整个算法内部走一遍。单纯形法的核心动作就三个：确定入基变量、确定出基变量、做一次高斯主元消去。但真正到你写代码的时候，你会不断卡在细节上：初始基可行解去哪找？等式约束怎么处理？人工变量的 M 应该多大？检验数符号方向为什么总是搞反？

这些问题调库时根本不会暴露。一旦你完整手写过大M法，再用linprog就会自然理解它内部在做什么。比如method='highs'默认方法为什么输出status=0还是status=2，为什么有时候会有“退化”提示，这些都能从算法内部逻辑找到依据。

1.3 我平时选择路线的三个判断条件

第一，问题规模。几十个变量以下、模型还要反复修改调试，我用linprog先跑通；但如果要复现课本算法、做教学演示，或者要截取迭代中间过程的单纯形表，我只能自己写。

第二，结果可信度。刚搭的模型有低级建模错误时，手写版反而更像“照妖镜”。它会给出奇怪的中间值，逼你去检查约束方向；库函数救不了这种错误。

第三，环境限制。有些部署环境装不了 SciPy，或者项目规范禁止引入不必要的第三方依赖，这时候手写一个纯 Python 版本就是兜底方案。当然它只适合小规模，大规模还是要靠成熟求解器。

2. 大M法为什么能处理“没有初始可行解”的问题

2.1 标准形是所有单纯形法的起点

几乎所有单纯形法教材的第一步，都是先把一般线性规划化为标准形：

$$ \min \ c^T x $$

$$ s.t. \ A x = b, \ x \ge 0 $$

只有等式约束，只有非负变量，目标统一极小化。这一步做不干净，后面全乱。遇到小于等于约束，就加一个非负松弛变量；遇到大于等于约束，就减一个非负剩余变量；遇到等式约束，就什么也不加，保持不动。

问题是，加了剩余变量之后，约束矩阵里未必能立刻找到一个现成的单位矩阵作为初始基。比如一个等式约束或者一个大于等于约束，没有松弛变量可以直接放进基里。单纯形法需要一个初始基可行解，没有它，迭代根本起不了步。这就是大M法登场的原因。

2.2 人工变量与大M惩罚：先加再赶

思路挺直观：给那些“缺基”的约束行，硬塞一个只属于该行的人工变量，让它先充当初始基变量，让算法能跑起来。但人工变量是编出来的，不是原模型的一部分，所以必须在目标函数里给每个工变量配上极大的惩罚系数 M。

如果是极小化问题，那就是目标函数加上 $M \sum x_a$。这个 M 足够大，大到哪怕原目标因此损失一丁点，也会逼迫单纯形法在迭代过程中把人工变量赶出基。迭代结束之后，人工变量如果不是基变量，说明它已经被降到零，等于没存在过；如果人工变量还在基里而且是正数，说明原问题根本没有可行解。

这里最容易被忽略的一点是：M 并不是数学上的“无穷大”，它只是一个数值上足够大的惩罚。M 取多少合适，本身就是实战里的一个坑。我后面专门讲。

2.3 一个能讲清楚全流程的测试模型

为了后面两套实现方便对拍，我固定用下面这个测试模型：

原问题为极小化：

$$ \min \ z = -3x_1 + x_2 + x_3 $$

约束：

$$ x_1 - 2x_2 + x_3 \le 11 $$

$$ -4x_1 + x_2 + 2x_3 \ge 3 $$

$$ -2x_1 + x_3 = 1 $$

$$ x_1, x_2, x_3 \ge 0 $$

这个问题既有小于等于、又有大于等于、还有等式约束，一个模型就覆盖了大M法需要的所有机制。实际约束不等式不太好猜最优解，因而很适合用来验证算法正确性。我知道的最优结果是 $(x_1, x_2, x_3) = (4, 1, 9)$，目标值 $z = -2$，后面两套实现都必须收敛到这个值。

3. 不借助任何库函数，纯Python手写大M法

3.1 单纯形表的数据结构约定

手写版我选择最基础的数据结构：用二维列表tableau存单纯形表，约束行在上面，目标函数行放最下面。每个约束行先放各变量的系数，最后一列放右端常数，目标行的最后一列放当前目标值。

所有代码只用 Python 标准库，不用 NumPy，不调用任何优化库。这样做的目的很纯粹：让每一个数学动作都显式写在循环里，而不是藏在矩阵运算里。你在代码里看到的tableau[i][j]，对应书上的 $a_{ij}$；看到的last_row[j]，对应的是检验数。

3.2 自动标准化：松弛变量、剩余变量和人工变量怎么进表

我写的函数会接收原始模型，自动完成标准化：

• 约束类型是<=，添加一个松弛变量，系数为+1；
• 约束类型是>=，先添加一个剩余变量，系数为-1，再添加一个工变量，系数为+1；
• 约束类型是=，只添加工变量，系数为+1。

这些人工变量会被记录在一个集合里，目标函数扩展到对应长度后，把这些列的目标系数全部改成 M。初始基也同步建立：松弛变量和人工变量所在的列就是初始基。

这里有一个很容易错的地方：人工变量进目标函数时，所带的 M 惩罚系数需要和它在目标行里的系数保持一致。但单纯形表要求基变量在目标行里的系数必须为 0，所以目标行要先对所有基变量列做一次消元，把基变量在目标行中的值清零。不消的话，检验数根本是错的，第一步入基选列就乱掉。

3.3 迭代主循环：检验数、入基、出基、主元消去

单纯形法的迭代部分可以拆成四步：

1. 计算检验数。目标行经过基变量列清零后，非基变量列的数值就是对应检验数。
2. 选入基变量。极小化问题里，选择检验数负得最厉害的非基变量列。
3. 选出基变量。对每一行，用右端常数除以入基列系数，只保留系数为正的行，取最小比值对应的行为出基行。
4. 主元消去。先归一化出基行，再把它下面所有行（包括目标行）在入基列的系数全部消成 0。

如果所有检验数都不小于 0，迭代结束，当前基就是最优基。如果在选入基变量时没有可选列，说明无界；如果在选入基变量后没有任何行满足正系数比值条件，同样是无界。最后要检查人工变量是否还在基里，如果在且取值为正，则原问题无可行解。

3.4 完整代码和运行结果

下面是我实际能跑的完整实现：

def big_m_simplex(A, b, c, cons_types, M=1e6, eps=1e-9): m = len(A) n = len(c) ext_A = [row[:] for row in A] basis = [] artificial_cols = set() for i, typ in enumerate(cons_types): if typ == '<=': col = len(ext_A[0]) for r in range(m): ext_A[r].append(1 if r == i else 0) basis.append(col) elif typ == '>=': col_surplus = len(ext_A[0]) for r in range(m): ext_A[r].append(-1 if r == i else 0) col_art = len(ext_A[0]) for r in range(m): ext_A[r].append(1 if r == i else 0) basis.append(col_art) artificial_cols.add(col_art) else: col_art = len(ext_A[0]) for r in range(m): ext_A[r].append(1 if r == i else 0) basis.append(col_art) artificial_cols.add(col_art) ext_cols = len(ext_A[0]) c_ext = list(c) + [0] * (ext_cols - n) for j in artificial_cols: c_ext[j] = M tableau = [row[:] + [b[i]] for i, row in enumerate(ext_A)] tableau.append(c_ext[:] + [0]) for row_idx, var_idx in enumerate(basis): coeff = tableau[-1][var_idx] if abs(coeff) > eps: for j in range(ext_cols + 1): tableau[-1][j] -= coeff * tableau[row_idx][j] while True: last_row = tableau[-1] enter = None min_cost = -eps for j in range(ext_cols): if j in basis: continue if last_row[j] < min_cost: min_cost = last_row[j] enter = j if enter is None: break leave = None min_ratio = float('inf') for i in range(m): a = tableau[i][enter] if a > eps: ratio = tableau[i][-1] / a if ratio < min_ratio - eps: min_ratio = ratio leave = i if leave is None: return {'status': 'unbounded'} pivot = tableau[leave][enter] for j in range(ext_cols + 1): tableau[leave][j] /= pivot for i in range(m + 1): if i == leave: continue factor = tableau[i][enter] if abs(factor) > eps: for j in range(ext_cols + 1): tableau[i][j] -= factor * tableau[leave][j] basis[leave] = enter for row_idx, var_idx in enumerate(basis): if var_idx in artificial_cols and abs(tableau[row_idx][-1]) > eps: return {'status': 'infeasible'} x = [0.0] * ext_cols for row_idx, var_idx in enumerate(basis): x[var_idx] = tableau[row_idx][-1] z = -tableau[-1][-1] return {'status': 'optimal', 'x': x[:n], 'z': z}

调用时这样写：

A = [ [1, -2, 1], [-4, 1, 2], [-2, 0, 1], ] b = [11, 3, 1] c = [-3, 1, 1] cons_types = ['<=', '>=', '='] result = big_m_simplex(A, b, c, cons_types) print(result)

运行结果是：

{'status': 'optimal', 'x': [4.0, 1.0, 9.0], 'z': -2.0}

和处理论最优解完全一致。整个过程只用到了列表切片、循环和基本四则运算，没有任何第三方库。我给这段代码留了一个eps=1e-9的容差，实际数值计算里浮点误差不可能完全为零，判断时不要用严格相等。

4. 用库函数实现：十行代码跑通同一问题

4.1 linprog 接口基本约定

SciPy 的linprog是长期维护的线性规划接口，极小化形式约定如下：

$$ \min \ c^T x $$

约束分别是 $A_{ub} x \le b_{ub}$、$A_{eq} x = b_{eq}$、$x$ 的边界 bounds。极大化目标要自己转成极小化，把 c 整体取负即可；大于等于约束也建议统一转成小于等于，否则容易在边界条件上绕晕。

转我们那个模型时，要注意第二条大于等于约束：

$$ -4x_1 + x_2 + 2x_3 \ge 3 $$

等价于：

$$ 4x_1 - x_2 - 2x_3 \le -3 $$

所以这一行放进A_ub时系数全部取相反数，右端也取相反数。等式约束直接放第二组参数。变量边界统一是(0, None)，也就是非负且没有上界。

4.2 method 参数怎么选

现在的linprog默认method='highs'，这是 HiGHS 求解器的入口，内部既包含单纯形法也包含内点法，性能好、接口稳定。如果为了和本文主题呼应，想看单纯形法效果，可以在旧版本里用method='simplex'；但在新版 SciPy 里我更推荐直接highs，不需要也不建议手动切到老式方法。

一个常见误区是：highs默认解决的是极小化标准问题，如果直接拿极大化模型套进去，出来的目标值方向会反。你只需要记得，所有原问题的目标系数都在 c 里给定，linprog返回的fun永远是这套极小化约定下的目标值。

4.3 示例代码和结果

from scipy.optimize import linprog c = [-3, 1, 1] A_ub = [ [1, -2, 1], [4, -1, -2], ] b_ub = [11, -3] A_eq = [ [-2, 0, 1], ] b_eq = [1] bounds = [(0, None)] * 3 res = linprog(c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq, bounds=bounds, method='highs') print(res)

输出关键字段：

res.x = [4. 1. 9.] res.fun = -2.0 res.status = 0

status=0表示成功找到最优解。res.x和手写版结果一致，目标值都是-2。到这里，两套实现已经对上了。

5. 两套实现对拍：结果、差异和各自的适用边界

5.1 结果对拍：目标值一致还不够

目标值和最优解一致，只是第一层验证。真正有价值的对拍还要看：无可行解时都返回什么状态，无界时返回什么状态，遇到退化时谁先卡住。手写版我故意保留了返回状态字段，方便你在不同模型上做交叉验证。

我实际跑过多组小模型，手写版和linprog的最优目标值几乎都能对上，除非 M 设置过大导致浮点误差积累，或者约束矩阵本身是病态的。碰到对不上的情况，先检查你自己模型的方向有没有写反，不要一上来就怀疑算法。

5.2 六个维度的对比

这张表不是想说哪个更强，而是两条路线解决的是不同层面的问题。手写版把算法过程完全摊开，库函数把算法效率推到极致。

5.3 为什么说两者是互补而不是替代

有一次我调试一个带十几个等式约束的排班模型，linprog总是报不可行。我把同样的数据喂给手写大M法，发现它返回的也是不可行，而且能通过人工变量所在的基位置，反推是哪一行约束把可行域直接卡没了。库函数只能给你一个状态码，手写版却能告诉你问题出在哪个约束附近。

反过来，模型一旦到了几百个变量，手写版跑完一次优化可能比库函数慢几十倍，而且 M 取值稍有问题就会出现诡异的迭代震荡。这种规模下我几乎不会用自研轮子，直接用highs求结果。互补关系就体现在这里：小模型用手写版理解问题，大模型用库函数拿结果。

6. 实战中会踩的坑和我的缓解习惯

6.1 M 值不是越大越好

刚接触大M法时，我一度觉得 M 直接取1e8好了，反正是个惩罚。实际跑下来发现，M 太大会让目标函数里的惩罚数值过大，单纯形表在做主元消去时产生严重的数值误差，甚至让本应离开基的人工变量因为浮点判断失误而留在基里。M 太小又起不到惩罚效果，人工变量可能留在最优解里。人经验是：先看原目标系数量级，M 取它最大值的一百到一千倍，测试问题时我用1e6够了；如果结果出现奇怪的小数残差，再把 M 降到1e4到1e5试一次。工程上如果 M 的取值敏感，说明这个问题本身更适合用两阶段单纯形法，大M法只适合快速试算。

6.2 检验数与出基规则的符号约定

这是手写代码最容易翻车的地方。同样是检验数，有的教材用 $c_j - z_j$，有的用 $z_j - c_j$，极大化和极小化的入基方向又不一样，一个符号没对齐，迭代就永远不收敛。我的解决办法是在代码注释里固定一套约定：极小化问题，目标行先消去基变量系数，然后选择目标行中负得最多的非基变量列入基；所有检验数非负即停止。你自己写的时候也必须在注释里写清这一句，不然三天后再看代码保证忘。

6.3 退化循环与无界无解判定

退化现象就是多个出基行最小比值相同，导致某次迭代换了基变量但目标值没下降，算法可能多绕几轮甚至循环。正规做法是加 Bland 规则，也就是在入基、出基选择出现平局时，一律选最小列号、最小行号，能有效避免循环。我手写代码里用“最小比值时优先保留非人工变量行”来缓解人工变量引起的退化问题，步骤虽小，但效果明显。无界判定则很简单：入基变量选出来之后，如果所有约束行在入基列上的系数都小于等于零，意味着目标可以无限下降，直接返回无界。

6.4 linprog 版本和参数带来的坑

SciPy 版本更新很快，老教程里的method='interior-point'、method='simplex'在新版本里可能报警告甚至不可用，所以不要照抄旧代码。最好统一method='highs'。另外bounds参数默认是每个变量非负，但如果不显式写，某些版本的行为会让人误判负值变量；我建议无论每个变量是什么约束，都一个个写清楚。还有个高频错误是A_ub行数变量从第二个约束开始复制粘贴时错位，导致维度对不上，linprog会直接报错，这种报错反而是好事，能逼你在建模阶段发现问题。

最后再分享一个我自己的习惯：求解任何新的模型，我都会同时跑一遍手写版和库函数版，哪怕手写版只用来画单纯形表前几步。这一步很像“双重记账”，能自动暴露建模时的方向错误。等模型稳定了，再去追求性能，把手写版扔一边，专心依赖highs。这套组合打法，帮我少走了很多弯路。