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 2Bs.t. 2A B ≤ 8A 2B ≤ 7A, 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 的主要参数如下参数含义c目标函数系数一维数组注意是最小化方向A_ub“≤”约束的系数矩阵二维数组每行一个约束b_ub“≤”约束的右侧常数一维数组A_eq“”约束的系数矩阵b_eq“”约束的右侧常数bounds每个变量的上下界默认 (0, None) 就是 [0, ∞)method求解算法建议用 highs 或 highs-ds最容易踩的坑有三个。第一个是忘记 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_ubA_ub, b_ubb_ub, boundsbounds, methodhighs-ds) print(res)结果里 fun 是最小化目标值 -13也就是原问题的最大利润 13 万元x 是 [3, 2]说明 A 产品生产 3 件、B 产品生产 2 件最优。这里用 methodhighs-ds 是想要明确走双单纯形法路线如果版本较老或者不识别换成 methodhighs 即可它内部会自动选择合适算法。你可能觉得这个过程太轻松了确实生产环境里大部分线性规划问题都应该先用这类库函数解一遍。因为它用了大规模稀疏矩阵、预求解、对偶单纯形法甚至内点法手写代码很难在性能上超越它。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 xs.t. A x bx ≥ 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_geNone, b_geNone, A_eqNone, b_eqNone, M1e5, tol1e-8, max_iter200): 大M法单纯形求解器最小化方向。 A_le, b_le: 约束 A_ge, b_ge: 约束 A_eq, b_eq: 约束 返回 (最优x, 最小目标值, 状态) c np.asarray(c, dtypefloat) n len(c) if A_ge is None: A_ge np.zeros((0, n)) else: A_ge np.asarray(A_ge, dtypefloat).reshape(-1, n) if A_eq is None: A_eq np.zeros((0, n)) else: A_eq np.asarray(A_eq, dtypefloat).reshape(-1, n) A_le np.asarray(A_le, dtypefloat).reshape(-1, n) b_le np.asarray(b_le, dtypefloat).reshape(-1) b_ge np.asarray(b_ge, dtypefloat).reshape(-1) b_eq np.asarray(b_eq, dtypefloat).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约束还是 2AB≤8A2B≤7A,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 关键维度对比表从工程视角出发我把两种方式放在一张表里维度库函数 linprog手写大M法求解速度极快大规模稀疏优化较慢稠密矩阵求逆是瓶颈数值稳定性高成熟求解器有预求解和后处理中等M值选取不当会抖动可维护性代码量少升级方便逻辑复杂需要维护数学细节可解释性黑盒只看得到结果每一步都在掌控中适用规模上千变量没问题几十到上百变量还能接受教学价值低高能彻底理解单纯形法这张表不是绝对的。如果你的不等式数量很大但稀疏度很高手写一个基于稀疏矩阵的单纯形法也能做到不错但开发成本高得多。库函数的优势在于它把所有边界情况都处理掉了你不需要关心预求解器、反退化策略、迭代精度这些底层问题。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 - 2x2s.t. 2x1 x2 ≤ 20x1 3x2 ≤ 30-x1 - x2 ≤ -6x1 - x2 0x1, 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_ubA_ub, b_ubb_ub, A_eqA_eq, bounds[(0, None), (0, None)], methodhighs-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.6667x2 ≈ 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。我在代码里用了 tol1e-8 来判断检验数是否非负同时把基变量的检验数手动归零这些都是为了压制浮点误差。如果问题规模更大、数值范围差异大建议改用更高精度的库或做约束矩阵归一化。6.4 速查表报错信息对应原因返回状态最可能原因排查方向singular basis初始基矩阵奇异约束存在线性相关检查是否有重复约束或人工变量添加不完整infeasible约束矛盾或 M 不够大用小模型逐条检查约束区间unbounded缺约束或目标函数方向有误重看变量是否漏了非负边界max_iter退化或循环或问题规模太大增加迭代上限或启用 Bland 规则optimal正常结束与原约束方程做残差验证最后分享一个个人习惯手写算法时每完成一个版本就立刻拿之前用库函数解过的历史案例做回归。哪怕只是一个小数点偏差也说明某个实现细节有问题。这种“双路对照”的做法比单纯对着公式背算法有效得多。后面积累的案例多了手写代码就能反过来成为检验自己是否真正理解运筹学算法的最好镜子。