1. 从资源分配这个老问题说起线性规划这个词听起来像是课本里的东西但如果你做过任何跟资源分配沾边的事情你其实已经在跟它打交道了。工厂排产、物流调度、投资组合配比、甚至食堂配菜控制成本本质上都是同一类问题在一堆约束条件里找到一个让目标最优的解。我最早接触线性规划是在做一个排班系统的时候。当时的需求很朴素——给定一批员工、一批班次、每个人可工作的时间段和成本怎么排班能让总人力成本最低同时每个班次都有人覆盖。我一开始想用贪心结果发现贪心在某些边界情况下会翻车后来想用搜索状态空间又大到离谱。直到有人提醒我这不就是个线性规划问题吗我才意识到自己绕了一大圈其实站在一个非常成熟的数学工具门口。而线性规划里最经典、最核心的求解算法就是单纯形法Simplex。它诞生于上世纪四十年代由 George Dantzig 提出至今仍然是很多求解器的底层引擎之一。你可能会觉得奇怪都这么多年了难道没有更快的算法吗有比如内点法在某些场景下确实更快但单纯形法依然是实践中最常用的方法之一原因很简单——它在绝大多数实际问题上的表现非常稳而且它能给出一个顶点解这对很多业务场景来说比内部解更有意义。这篇文章我想做的事情很明确把单纯形法从课本里的公式拉回到能上手用的工具。我会讲清楚它到底在干什么、为什么这么干、怎么一步步算出来以及在实际写代码时有哪些坑。不管你是刚学运筹学的学生还是需要在项目里落地优化算法的工程师都能从里面拿到能直接用的东西。2. 线性规划的标准形态与单纯形法的切入点2.1 为什么所有线性规划都能写成同一个样子线性规划问题有各种长相有的求最大值有的求最小值有的约束是小于等于有的是大于等于还有的是等式变量有的要求非负有的没有限制。如果每个问题都单独处理算法会变得极其复杂。所以运筹学里有一个约定俗成的做法把所有线性规划问题统一转换成标准型。标准型通常长这样目标函数最大化 ( z c_1x_1 c_2x_2 \dots c_nx_n )约束条件所有约束都是等式 ( a_{i1}x_1 \dots a_{in}x_n b_i )且 ( b_i \geq 0 )变量约束所有 ( x_j \geq 0 )这个转换过程本身不复杂但它是理解单纯形法的前提。我见过不少人直接跳到单纯形表的操作步骤结果因为没搞懂标准型遇到大于等于约束时完全不知道松弛变量和剩余变量的区别算着算着就乱了。转换规则其实就三条目标函数统一为最大化如果原问题是求最小化直接取负号变成最大化最后结果再取负回来。不等式变等式小于等于约束加一个松弛变量大于等于约束减一个剩余变量等式约束不动。右边常数保证非负如果某个约束右边是负数两边同时乘 -1注意不等号方向要翻转。举个例子。假设有个问题最小化 z 3x1 2x2 约束 x1 x2 4 2x1 x2 6 x1, x2 0转换成标准型最大化 z -3x1 - 2x2 约束 x1 x2 - s1 4 2x1 x2 s2 6 x1, x2, s1, s2 0这里的 ( s_1 ) 是剩余变量( s_2 ) 是松弛变量。它们不是凑数的而是有实际含义的——松弛变量代表资源没用完的部分剩余变量代表超出最低要求的量。理解这一点后面看单纯形表里的数字就不会觉得是凭空冒出来的。2.2 可行域是个多面体最优解一定在顶点上标准型搞定之后下一步是理解线性规划的几何意义。所有约束条件在 n 维空间里划出一个区域叫做可行域。因为约束都是线性的所以可行域是一个凸多面体。这里有一个非常关键的定理如果线性规划有最优解那么至少存在一个最优解落在可行域的顶点上。这个定理是单纯形法的立足点。为什么因为可行域里的点有无穷多个你不可能一个个试。但顶点是有限的——对于一个有 m 个约束、n 个变量的标准型问题顶点的数量最多是 ( \binom{n}{m} ) 个。虽然这个数也可能很大但至少它是有限的而且可以通过系统化的方式去搜索。单纯形法的核心思想就是从一个顶点出发沿着可行域的边走到相邻的顶点每次移动都让目标函数值变得更好直到没有更好的邻居为止。这时候你就到了最优顶点。我第一次真正理解这个思路的时候觉得它特别像爬山——你在一个多面体的表面上走每一步都往更高的方向走走到一个局部最高的点。但因为可行域是凸的局部最高就是全局最高。这个性质是线性规划之所以好解的根本原因。2.3 基变量、非基变量与基本可行解要把几何直觉翻译成代数操作需要引入几个概念。在标准型里假设有 m 个约束、n 个变量n m。我们可以把变量分成两组基变量选 m 个变量它们的系数矩阵构成一个可逆的方阵。非基变量剩下的 n - m 个变量。令所有非基变量等于 0然后解出基变量的值得到的解叫做基本解。如果这个基本解里所有变量都非负那它就是一个基本可行解对应可行域的一个顶点。这个对应关系是单纯形法的灵魂。你在单纯形表里做的所有行变换本质上就是在不同的基之间切换从一个基本可行解跳到另一个基本可行解。我刚开始学的时候最大的困惑是为什么非基变量要设成 0后来想明白了非基变量设成 0意味着你站在某个顶点上那些没有进入基的变量暂时不参与你只用基变量来满足约束。当你把某个非基变量换入基时就相当于沿着一条边走到了相邻顶点。3. 单纯形法的完整计算流程拆解3.1 初始基本可行解的构造大M法与两阶段法单纯形法需要一个起点也就是一个初始基本可行解。对于所有约束都是小于等于且右边非负的问题松弛变量天然构成一个单位矩阵直接就能当初始基。但一旦出现大于等于或等于约束就没有这么幸运了。这时候有两种常用处理方式大M法给每个需要人工变量的约束引入一个人工变量并在目标函数里给它一个极大的惩罚系数 M最大化问题里是 -M。这样算法会优先把人工变量赶出基因为留着它们会让目标函数变得极差。大M法的优点是写起来简单缺点是 M 的取值很敏感——太小了惩罚不够太大了会导致数值不稳定。两阶段法第一阶段先最小化人工变量之和如果能把它们都降到 0说明找到了一个可行基第二阶段再在这个基上优化原目标函数。两阶段法在数值上更稳健实际求解器里用得更多。我在自己写求解器的时候一开始用的是大M法结果在一个系数跨度很大的问题上出现了精度问题——M 取 1e9 的时候浮点运算的误差把有效数字吃掉了。后来换成两阶段法问题就消失了。所以如果你要自己实现我建议直接上两阶段法别图省事。3.2 单纯形表的迭代入基、出基与最小比值检验有了初始基本可行解之后就进入迭代循环。每一轮迭代做三件事选入基变量看目标函数行也叫 z 行里非基变量的系数。如果是最大化问题选系数为正且最大的那个或者按 Bland 规则选最小的下标避免循环。这个变量进入基意味着我们沿着让目标函数增长最快的方向走。选出基变量对每一行用右边的常数除以该行入基变量的系数得到比值。在系数为正的行里选比值最小的那一行对应的基变量出基。这一步叫最小比值检验目的是保证新解仍然可行——如果选错了基变量会变成负数就跑到可行域外面去了。做行变换用高斯消元把入基变量所在的列变成单位向量同时更新 z 行。这一步是纯代数操作但手工算的时候最容易出错。我用一个具体例子走一遍。假设标准型是最大化 z 3x1 2x2 约束 x1 x2 s1 4 x1 3x2 s2 6 x1, x2, s1, s2 0初始基是 ( s_1, s_2 )初始基本可行解是 ( (0, 0, 4, 6) )z 0。初始单纯形表基变量x1x2s1s2右端项s111104s213016z-3-2000注意 z 行写的是目标函数系数的相反数这是单纯形表的惯例方便做行变换。第一轮x1 的系数 -3 最小对应原目标系数最大选 x1 入基。比值检验4/1 46/1 6最小是 4所以 s1 出基。行变换后基变量x1x2s1s2右端项x111104s202-112z013012第二轮x2 的系数是 1正数选 x2 入基。比值检验4/1 42/2 1最小是 1所以 s2 出基。行变换后基变量x1x2s1s2右端项x1101.5-0.53x201-0.50.51z003.50.513z 行所有系数都非负了说明没有改进空间迭代结束。最优解是 ( x_1 3, x_2 1 )最大值 z 13。这个例子很简单但它完整展示了单纯形法的每一步。手工算的时候我建议每一步都检查两件事右端项是否全部非负可行性z 行是否还有负系数最优性。这两个检查能帮你抓住大部分计算错误。3.3 退化、循环与Bland规则理论上单纯形法可能遇到一个问题退化。所谓退化就是某个基本可行解里有基变量等于 0。这时候最小比值检验可能出现多个相同的最小比值选不同的出基变量会导致不同的迭代路径极端情况下可能陷入循环——一直在几个基之间打转永远到不了最优解。循环在实际问题中很少见但不是不可能。学术界构造过专门触发循环的例子比如 Beale 例子。为了避免循环可以用Bland 规则入基时选下标最小的正系数变量出基时在比值相同的行里选下标最小的基变量。Bland 规则能保证算法终止但代价是迭代次数可能变多。我在实际项目里的做法是默认用 Dantzig 规则选最大正系数同时设置一个迭代上限比如 10000 次。如果超了就切换到 Bland 规则重新跑。这样在绝大多数情况下享受 Dantzig 的快速收敛极端情况下也有兜底。4. 从手算到代码实现单纯形法的关键决策4.1 数据结构选择表格法还是修订单纯形法如果你只是教学或者处理小规模问题直接用表格法就够了——用一个二维数组存单纯形表每次迭代做行变换。代码直观调试方便。但如果你要处理几百上千个变量的问题表格法就力不从心了。原因有两个一是每次迭代都更新整个表计算量大二是浮点误差会累积导致数值不稳定。这时候需要修订单纯形法。它不存整个表只存基矩阵的逆 ( B^{-1} )每次迭代通过更新 ( B^{-1} ) 来计算需要的列。这样内存占用小而且可以用 LU 分解等数值技巧提高稳定性。我自己的经验是变量数在 100 以内表格法完全够用超过 500就该考虑修订单纯形法或者直接调用成熟求解器了。自己从头写一个工业级求解器投入产出比并不高。4.2 数值稳定性浮点误差是怎么毁掉结果的单纯形法对浮点误差很敏感。一个典型场景是某个变量理论上应该等于 0但因为误差变成了 1e-12。如果这个变量在后续迭代里被选入基误差会被放大最终导致结果完全错误。常见的应对手段有相对误差阈值判断一个数是否为零时不用绝对等于 0而是看它是否小于某个相对阈值比如 1e-9 乘以该行最大元素的绝对值。重新计算基矩阵的逆不要一直用更新公式每隔若干次迭代重新从原始矩阵算一次 ( B^{-1} )把累积误差清掉。使用有理数运算对于小规模问题可以用分数精确计算完全避免浮点误差。Python 的 fractions 模块就能做这件事但速度会慢很多。我在一个投资组合优化的小工具里用过有理数版本变量数只有 20 多个跑起来完全没问题而且结果干净漂亮没有任何应该是 0 但显示 1e-16的尴尬。4.3 用Python实现一个最小可用的单纯形法下面是一个基于表格法的简化实现支持标准型的最大化问题。代码不长但覆盖了核心逻辑。import numpy as np def simplex(c, A, b, max_iter1000): c: 目标函数系数 (n,) A: 约束系数矩阵 (m, n) b: 右端项 (m,) 返回: (最优值, 最优解, 是否成功) m, n A.shape # 构造初始单纯形表 tableau np.zeros((m 1, n m 1)) tableau[:m, :n] A tableau[:m, n:nm] np.eye(m) tableau[:m, -1] b tableau[m, :n] -c # z 行 basis list(range(n, n m)) # 初始基变量下标 for _ in range(max_iter): # 检查最优性 z_row tableau[m, :-1] if np.all(z_row -1e-9): break # 选入基变量 enter np.argmin(z_row) # 最小比值检验 ratios [] for i in range(m): if tableau[i, enter] 1e-9: ratios.append(tableau[i, -1] / tableau[i, enter]) else: ratios.append(np.inf) if all(r np.inf for r in ratios): return None, None, False # 无界 leave int(np.argmin(ratios)) # 行变换 pivot tableau[leave, enter] tableau[leave, :] / pivot for i in range(m 1): if i ! leave: factor tableau[i, enter] tableau[i, :] - factor * tableau[leave, :] basis[leave] enter # 提取解 x np.zeros(n) for i, var in enumerate(basis): if var n: x[var] tableau[i, -1] return tableau[m, -1], x, True这段代码有几个地方值得说明。第一判断最优性时用了-1e-9而不是0这是为了容忍浮点误差。第二最小比值检验里只考虑系数大于1e-9的行避免除以零或负数。第三如果所有比值都是无穷大说明问题无界直接返回失败。你可以用前面的例子测试一下c np.array([3, 2]) A np.array([[1, 1], [1, 3]]) b np.array([4, 6]) val, x, ok simplex(c, A, b) print(val, x, ok) # 应该输出 13.0, [3. 1.], True这个实现只支持小于等于约束实际使用中你需要先做标准型转换或者扩展代码支持人工变量。但作为理解算法逻辑的起点它已经足够了。5. 单纯形法在实际场景中的表现与边界5.1 它擅长什么稀疏约束与顶点解单纯形法在稀疏约束矩阵的问题上表现特别好。所谓稀疏就是大部分系数是 0。现实中的很多问题都是稀疏的——比如网络流问题每个节点只和少数几个节点相连再比如排班问题每个员工只涉及少数几个班次。稀疏性带来的好处是每次迭代实际需要计算的列很少修订单纯形法可以利用这一点大幅加速。这也是为什么很多商业求解器在处理大规模稀疏问题时默认还是用单纯形法而不是内点法。另一个优势是顶点解。内点法给出的是可行域内部的一个点虽然目标值最优但可能所有变量都是非零的小数。而单纯形法给出的顶点解很多变量天然就是 0。在排班、选址、选品这类场景里0 或非 0的决策比每个都分一点更符合业务直觉。5.2 它不擅长什么大规模稠密问题与数值病态问题单纯形法最怕的是稠密矩阵。如果约束矩阵里几乎没有 0每次迭代的计算量会急剧上升。变量数上千、约束数上千的稠密问题单纯形法可能会跑很久。另一个软肋是数值病态。如果矩阵的条件数很大浮点误差会非常严重甚至导致算法给出错误的最优解。这时候内点法通常更稳健因为它不依赖基矩阵的逆。我在一个图像处理相关的优化问题里遇到过这种情况约束矩阵的条件数接近 1e12单纯形法跑出来的结果和预期差了十万八千里。换成内点法之后结果就正常了。所以选算法不能只看哪个更有名得看问题本身的特点。5.3 什么时候该用现成求解器什么时候自己写我的建议很直接教学、原型验证、变量数小于 50自己写能加深理解调试也方便。生产环境、变量数超过 100用成熟求解器比如开源的 HiGHS、GLPK或者商业的 Gurobi、CPLEX。它们在数值稳定性、预处理、并行化上的积累不是个人短时间能追上的。特殊结构问题如果问题有特殊结构比如网络流、指派问题可以用专门的算法往往比通用单纯形法快几个数量级。自己写求解器最大的价值是知道里面在发生什么。当你用 Gurobi 跑出一个奇怪结果时如果你理解单纯形法的原理就能判断是模型写错了、数值出问题了还是问题本身无界。这种判断力是调包调不出来的。6. 几个我踩过的坑和对应的处理方式6.1 忘记处理无界和无可行解线性规划不一定有最优解。可能无界目标函数可以无限增大也可能无可行解约束互相矛盾。我早期写的代码只处理了有最优解的情况遇到无界问题就死循环遇到无可行解就返回一堆垃圾数字。正确的做法是在算法里显式检测这两种情况无界入基变量选定后如果所有行的系数都小于等于 0说明这个方向可以无限延伸问题无界。无可行解两阶段法的第一阶段结束时如果人工变量之和大于 0说明原问题无可行解。这两个检测一定要加否则你的求解器在真实数据上会表现得非常不可靠。6.2 目标函数方向的混淆最大化还是最小化这个看起来是小事但我在项目里见过不止一次因为方向搞反而导致结果完全相反的案例。尤其是当代码里同时处理两种问题时很容易在某个分支忘了取负号。我的做法是在入口处统一转换成最大化所有内部逻辑只处理最大化。出口处如果需要再转换回原始方向。这样只需要在一个地方处理方向问题减少出错概率。6.3 约束右端项为负时的处理标准型要求右端项非负。如果原始问题里某个约束的右端项是负数需要两边乘 -1。但乘 -1 之后不等号方向会翻转这时候松弛变量和剩余变量的选择也要跟着变。我踩过的坑是只记得乘 -1忘了翻转不等号结果把大于等于当成了小于等于加错了变量类型。后来我养成了一个习惯转换标准型时每一步都在纸上写清楚转换完再对照检查一遍。这个习惯帮我省了很多调试时间。6.4 迭代次数异常时的排查思路如果你的单纯形法实现跑了很多次迭代还没收敛可以按这个顺序排查检查最小比值检验是否正确处理了系数为 0 或负数的情况。检查是否有退化导致的循环尝试切换到 Bland 规则。检查浮点误差是否导致某个本该为 0 的变量被反复选入基。检查问题本身是否无界但无界检测没触发。我遇到过一次迭代不收敛的情况最后发现是最小比值检验里用了 0而不是 1e-9导致系数为 0 的行也被纳入比值计算除出了无穷大选错了出基变量。改成严格大于一个小阈值之后问题就解决了。7. 单纯形法给我的几点启发单纯形法最让我佩服的地方不是它有多快而是它把一个看似无从下手的连续优化问题转化成了在有限个顶点上的系统性搜索。这个转化本身就是一种非常漂亮的思维方式当你面对无限多的可能性时先找到结构把无限变成有限。另一个启发是关于够用就好。内点法在理论上有多项式时间复杂度单纯形法在最坏情况下是指数级的。但在实践中单纯形法往往更快。这提醒我理论上的最优和实际中的最优是两回事做工程决策时不能只看理论指标。最后一点是关于理解深度。我用过很多求解器但真正让我对线性规划有信心的是自己从头实现过一遍单纯形法。知道每一步在干什么知道哪里可能出错知道怎么排查这种底气是调包给不了的。如果你正在学运筹学或者优化算法我强烈建议你至少手写一次单纯形法哪怕只处理两三个变量的问题。写完之后你对线性规划的理解会上一个台阶。
