简介这份资源面向学习最优化方法、线性规划课程的高校学生与工程技术人员聚焦单纯形法求解中的大M法与两阶段法编程实现帮助读者理解算法原理并完成课程作业或实际建模求解。压缩包共3个文件均为MATLAB的m脚本整体约4KB体量轻便包含主程序入口、两阶段法函数与大M法单纯形函数代码注释详尽、逻辑通俗易懂。主程序负责输入约束方程与目标函数自动调用两阶段法或大M法完成规划求解其中两阶段法的单纯形迭代通过调用大M法函数实现便于对照两种方法的异同。目前已有2547人学习下载读者可借此掌握单纯形表构建、人工变量处理、迭代终止判断等关键环节并在此基础上修改参数求解自己的线性规划问题是入门与复习单纯形法的实用参考。1. 单纯形法为什么需要人造变量从标准型到初始基可行解的鸿沟线性规划的标准型看起来人畜无害目标函数求极大约束全是等式变量非负。教科书上写得很轻松——找到初始基可行解然后迭代换基直到检验数全部非正。但真正动手写代码的人会撞上一堵墙约束里那些大于等于和等式根本凑不出一组现成的单位基向量。没有单位基单纯形表连第一行都填不出来。这就是大M法和两阶段法存在的理由。它们不是单纯形法的高级变体而是专门用来跨越初始基可行解不存在这道鸿沟的工程手段。很多人学最优化方法时把这两个方法当成两套独立算法背下来其实它们解决的是同一个问题只是代价函数的设计哲学不同。搞懂这一点才算真正理解了单纯形法的完整链路。这份资源把大M法和两阶段法的程序实现放在一起恰好覆盖了从数学模型到可运行代码的关键环节。下面从为什么需要人造变量开始一层层拆到代码里的具体处理。2. 人造变量的本质为什么不能直接对不等式做消元2.1 松弛变量与剩余变量的分工先把约束类型理清楚。标准型要求所有约束是等式所以不等式要转化小于等于约束左边加一个松弛变量系数 1直接充当基变量大于等于约束左边减一个剩余变量系数 -1但 -1 不能当基变量基变量要求对应列是单位向量等式约束没有任何多余变量可以借用问题就出在后两种。剩余变量系数是 -1等式约束连个多余的变量都没有。单纯形法的第一步是构造初始单纯形表要求基矩阵是单位阵。如果约束里凑不出单位阵迭代就无法启动。常见的实现思路是对每个缺基变量的约束人为添加一个变量让它在该约束中系数为 1在其他约束中系数为 0。这样所有添加的变量就构成了一组单位基。这些变量就是人造变量。但人造变量是凭空造出来的它们并不属于原问题。如果最终解里人造变量不为零说明原问题无可行解。所以人造变量必须被赶出基或者被迫取零值。大M法和两阶段法的区别就在于用什么手段逼人造变量归零。2.2 大M法的惩罚逻辑大M法的思路很直接在目标函数里给人造变量加一个巨大的惩罚系数 M。如果是求极大值人造变量的系数就是 -M求极小值则相反。M 是一个理论上无穷大的正数。这样做的效果是只要人造变量还留在基里且取正值目标函数就会被扣掉一个巨大的量单纯形法在迭代中会优先把它换出去。代码实现时M 不能真的取无穷大通常取一个远大于所有原始系数的数比如 10^6 或 10^8。# 大M法构造初始单纯形表的核心逻辑 def build_big_m_table(c, A, b, constraint_types, M1e6): c: 原始目标函数系数求极大 A: 约束系数矩阵 b: 右端项 constraint_types: 每个约束的类型 、、 M: 惩罚系数 m, n A.shape total_vars n artificial_indices [] for i, ctype in enumerate(constraint_types): if ctype : # 加松弛变量系数 1可作基变量 slack np.zeros((m, 1)) slack[i] 1 A np.hstack([A, slack]) c np.append(c, 0) total_vars 1 elif ctype : # 加剩余变量系数 -1不能作基变量 surplus np.zeros((m, 1)) surplus[i] -1 A np.hstack([A, surplus]) c np.append(c, 0) total_vars 1 # 再加人造变量系数 1 art np.zeros((m, 1)) art[i] 1 A np.hstack([A, art]) c np.append(c, -M) # 求极大惩罚为 -M artificial_indices.append(total_vars) total_vars 1 else: # art np.zeros((m, 1)) art[i] 1 A np.hstack([A, art]) c np.append(c, -M) artificial_indices.append(total_vars) total_vars 1 return A, b, c, artificial_indices这段代码的关键在于松弛变量系数 1 可以直接当基变量剩余变量系数 -1 不行所以大于等于约束需要同时加剩余变量和人造变量。等式约束只需要加人造变量。人造变量在目标函数里的系数是 -M这个 -M 会随着迭代进入检验数计算驱动人造变量出基。M 的取值是个工程问题。取太小惩罚力度不够人造变量可能赖在基里不走取太大浮点运算时容易和其他系数产生数量级差异导致数值不稳定。常见的做法是取原始系数最大绝对值的 1000 到 10000 倍。2.3 两阶段法的分步策略两阶段法不依赖 M 这个魔法数字而是把问题拆成两个阶段第一阶段构造一个辅助问题目标函数是极小化所有人造变量之和。如果辅助问题的最优值为零说明人造变量可以全部取零原问题有可行解如果最优值大于零原问题无可行解。第二阶段在第一阶段得到的基可行解基础上去掉人造变量列恢复原始目标函数继续单纯形迭代。# 两阶段法第一阶段极小化人造变量之和 def phase_one(A, b, artificial_indices): 构造辅助目标函数min sum(artificial_vars) 等价于 max -sum(artificial_vars) n A.shape[1] c_phase1 np.zeros(n) for idx in artificial_indices: c_phase1[idx] -1 # 求极大化 -sum等价于极小化 sum # 用单纯形法求解辅助问题 result simplex(A, b, c_phase1, basisartificial_indices) if abs(result[objective]) 1e-8: return None # 原问题无可行解 # 检查人造变量是否还在基里若在且值为零尝试换出 return result两阶段法的优势在于不需要选 M。M 选得不好会导致数值问题而两阶段法把可行性判断和最优性优化彻底分开逻辑更清晰。代价是要跑两遍单纯形计算量略大。从设计取舍看大M法代码更紧凑适合手写实现和教学演示两阶段法数值更稳健适合放进通用求解器。这份资源把两者都实现了对比着看能理解为什么工业级求解器更倾向两阶段法。3. 单纯形表的数据结构基变量、检验数与非基变量怎么组织3.1 单纯形表的矩阵表示单纯形法的每次迭代本质上是一次矩阵变换。把约束写成 $Ax b$把 $A$ 按基变量和非基变量分块$A [B | N]$其中 $B$ 是基矩阵$N$ 是非基矩阵。对应的变量也分成 $x_B$ 和 $x_N$。约束变为 $B x_B N x_N b$所以 $x_B B^{-1}b - B^{-1}N x_N$。代入目标函数 $z c_B^T x_B c_N^T x_N$得到$$z c_B^T B^{-1}b (c_N^T - c_B^T B^{-1}N) x_N$$令非基变量 $x_N 0$得到基可行解 $x_B B^{-1}b$目标值 $z c_B^T B^{-1}b$。检验数就是 $\sigma_N c_N^T - c_B^T B^{-1}N$。代码里不需要每次都求 $B^{-1}$而是用高斯-若尔当消元维护整张单纯形表。表的结构是基变量$x_1$$x_2$...$x_n$右端项$x_{B1}$$a_{11}$$a_{12}$...$a_{1n}$$b_1$$x_{B2}$$a_{21}$$a_{22}$...$a_{2n}$$b_2$..................检验数$\sigma_1$$\sigma_2$...$\sigma_n$$z$每次迭代选一个检验数为正的非基变量入基求极大时再按最小比值法则选一个基变量出基然后做一次初等行变换把入基变量对应的列变成单位向量。3.2 入基与出基的判定细节入基变量的选择直接影响收敛速度。最简单的规则是选检验数最大的那个叫 Dantzig 规则。但 Dantzig 规则容易在退化问题上循环。常见的改进是 Bland 规则选下标最小的正检验数变量入基同时出基也选下标最小的保证不循环但收敛慢。出基判定用最小比值法则对每个系数 $a_{ik} 0$ 的行计算 $b_i / a_{ik}$选最小的那个对应的基变量出基。如果所有 $a_{ik} \leq 0$说明问题无界。def simplex_iteration(table, basis, c, maximizeTrue): table: 单纯形表最后一行为检验数 basis: 基变量下标列表 c: 目标函数系数 m, n table.shape n - 1 # 最后一列是右端项 # 1. 选入基变量找最大正检验数求极大 reduced_costs table[-1, :n] if maximize: entering np.argmax(reduced_costs) if reduced_costs[entering] 1e-8: return None # 已最优 else: entering np.argmin(reduced_costs) if reduced_costs[entering] -1e-8: return None # 2. 最小比值法则选出基变量 ratios [] for i in range(m - 1): if table[i, entering] 1e-8: ratios.append(table[i, -1] / table[i, entering]) else: ratios.append(np.inf) if min(ratios) np.inf: return unbounded # 无界解 leaving np.argmin(ratios) # 3. 初等行变换把 entering 列变成单位向量 pivot table[leaving, entering] table[leaving, :] / pivot for i in range(m): if i ! leaving: table[i, :] - table[i, entering] * table[leaving, :] basis[leaving] entering return table, basis这段代码里有个容易翻车的点浮点数比较不能用 0必须用容差。检验数判断最优时用 1e-8比值判断时用 1e-8。容差取太小会误判取太大会丢精度。常见的做法是取 1e-9 到 1e-6 之间具体看问题规模。3.3 退化与循环的处理退化是指基可行解里出现零分量导致最小比值出现多个相同的最小值。退化本身不致命但可能引发循环——基变量组合重复出现迭代永远不终止。处理循环的常见手段有三种Bland 规则、字典序法、扰动法。Bland 规则实现最简单但收敛慢字典序法在出基判定时引入额外的比较维度实现复杂但效率好扰动法给右端项加微小扰动破坏退化条件实现也不复杂。代码里如果只做教学演示用 Bland 规则就够了。如果要做通用求解器建议用字典序法。这份资源里的实现如果没处理循环遇到退化问题可能会卡住这是使用时需要留意的。4. 从数学公式到可运行代码大M法与两阶段法的实现差异4.1 大M法的完整迭代流程大M法的代码结构可以拆成四步构造初始表、迭代换基、判断终止、提取解。构造初始表时需要根据约束类型决定加松弛变量、剩余变量还是人造变量。这一步的代码在 2.2 节已经展示过。迭代换基就是反复调用 3.2 节的simplex_iteration。终止条件有两个检验数全部非正最优或者入基列所有系数非正无界。提取解的时候要注意如果最优解里人造变量不为零说明原问题无可行解。这个判断不能省否则会把一个无可行解的问题当成有解输出。def solve_big_m(c, A, b, constraint_types, M1e6): A, b, c_ext, art_indices build_big_m_table(c, A, b, constraint_types, M) m, n A.shape # 初始基变量松弛变量和人造变量 basis [] for i, ctype in enumerate(constraint_types): if ctype : basis.append(original_n i) # 松弛变量下标 else: basis.append(art_indices.pop(0)) # 构造单纯形表 table np.hstack([A, b.reshape(-1, 1)]) # 计算初始检验数sigma c - c_B * B^{-1} * A c_B c_ext[basis] table[-1, :] c_ext - c_B table[:-1, :] # 迭代 while True: result simplex_iteration(table, basis, c_ext) if result is None: break if result unbounded: return {status: unbounded} # 检查人造变量是否为零 solution np.zeros(n) for i, b_idx in enumerate(basis): solution[b_idx] table[i, -1] for idx in art_indices: if abs(solution[idx]) 1e-8: return {status: infeasible} return {status: optimal, solution: solution[:original_n], objective: table[-1, -1]}这段代码里有个细节初始检验数的计算不能直接抄目标函数系数因为基变量的检验数必须是零。需要用 $c_B$ 做一次修正。很多人手写单纯形法时在这里翻车迭代半天发现检验数不对其实是初始表没构造对。4.2 两阶段法的阶段衔接两阶段法的代码比大M法多一个阶段切换的逻辑。第一阶段结束后如果辅助问题最优值为零需要做两件事去掉人造变量列恢复原始目标函数重新计算检验数。去掉人造变量列时要注意如果某个人造变量还在基里但值为零不能直接删列否则基矩阵会缺一个单位向量。常见的处理是找一个非基的非人造变量通过行变换把它换进基再删掉人造变量列。如果找不到这样的变量说明该约束是冗余的可以直接删掉整行。def solve_two_phase(c, A, b, constraint_types): # 第一阶段 A1, b1, c1, art_indices build_phase_one_table(A, b, constraint_types) result1 simplex_solve(A1, b1, c1, basisart_indices) if abs(result1[objective]) 1e-8: return {status: infeasible} # 阶段切换去掉人造变量列恢复原始目标 keep_cols [i for i in range(A1.shape[1]) if i not in art_indices] A2 A1[:, keep_cols] c2 np.zeros(A2.shape[1]) c2[:len(c)] c # 原始变量系数 # 用第一阶段的基作为第二阶段的初始基 basis2 [keep_cols.index(b) for b in result1[basis] if b in keep_cols] # 第二阶段迭代 result2 simplex_solve(A2, b1, c2, basisbasis2) return result2阶段切换时最容易出问题的地方是基变量的映射。第一阶段的基变量下标是原始变量、松弛变量、剩余变量、人造变量的混合下标去掉人造变量列后下标需要重新映射。如果映射错了第二阶段的初始基就不是可行基迭代会从错误的地方开始。4.3 两种方法的数值表现对比对比维度大M法两阶段法代码复杂度低一次构造高需要阶段切换数值稳定性依赖 M 的取值较稳定迭代次数通常较少通常较多无可行解判断检查人造变量是否为零检查第一阶段最优值适用场景教学演示、小规模问题通用求解器、大规模问题从工程角度看大M法的 M 是个魔法数字取大了数值不稳定取小了惩罚不够。两阶段法没有这个问题但代码量更大。如果只是手写一个求解器验证算法大M法更快如果要处理实际数据两阶段法更可靠。5. 避坑与排查单纯形法实现中最容易翻车的五个地方5.1 检验数符号搞反导致迭代方向错误现象程序迭代几次后目标函数值越来越小或者直接报无界但问题明明有有限最优解。原因求极大和求极小的检验数判定方向相反。求极大时选正检验数入基求极小时选负检验数入基。如果目标函数是求极小但代码里用了求极大的判定逻辑迭代方向就反了。解决在代码里显式区分 maximize 和 minimize或者在构造目标函数时统一转成求极大。转法很简单min z 等价于 max -z把目标函数系数全部取反即可。5.2 最小比值法则漏判无界现象程序在某次迭代中找不到出基变量但没有报无界而是继续迭代或者崩溃。原因最小比值法则要求只对入基列系数为正的行计算比值。如果所有系数都非正说明入基变量可以无限增大而不违反约束问题无界。代码里如果没判断这个条件argmin会返回一个无效下标。解决在计算比值前先检查入基列是否有正系数。如果没有直接返回无界。代码里用if min(ratios) np.inf或者if all(a 1e-8 for a in column)来判断。5.3 浮点误差累积导致检验数误判现象程序在接近最优时反复迭代检验数在零附近震荡永远不满足终止条件。原因浮点运算有误差理论上应该为零的检验数可能变成 1e-15 或者 -1e-15。如果终止条件用sigma 0严格判断就会因为一个极小的正数继续迭代。解决所有浮点比较都加容差。检验数判断用sigma 1e-8比值判断用a 1e-8。容差不能太大否则会丢失精度也不能太小否则起不到作用。1e-8 是个常用的经验值。5.4 人造变量未完全出基导致误判无可行解现象大M法求解一个明明有可行解的问题程序却报无可行解。原因人造变量在最优解里可能还留在基里但值为零。如果代码只检查人造变量是否在基里而不检查其值是否为零就会误判。或者 M 取得太小人造变量没有被完全赶出基值不为零。解决判断无可行解的条件是人造变量的值大于容差而不是人造变量是否在基里。同时检查 M 的取值是否足够大通常取原始系数最大绝对值的 1000 倍以上。5.5 阶段切换时基变量映射错误现象两阶段法第一阶段正常结束第二阶段一开始就报无界或者迭代到错误的最优解。原因第一阶段结束后基变量列表里可能包含人造变量。去掉人造变量列时如果直接删列而不调整基变量下标基矩阵就不再是单位阵第二阶段的初始基不可行。解决阶段切换时先检查基里是否有人造变量。如果有且值为零找一个非基的非人造变量换进去如果找不到说明该约束冗余删掉整行。然后再删人造变量列重新映射基变量下标。6. 设计取舍为什么工业求解器很少用大M法大M法和两阶段法在教科书里是并列的但翻开源码级的求解器实现会发现两阶段法或者它的变体才是主流。这不是偶然而是数值计算领域几十年踩坑的结果。M 的选取本质上是一个没有正确答案的问题。取 10^6遇到系数本身就是 10^5 量级的问题惩罚力度不够取 10^12浮点数的有效位数只有 15 到 16 位M 和其他系数的加减运算会直接丢精度。更麻烦的是M 的影响会渗透到整个单纯形表的每一次行变换中误差不断累积最终可能导致基变量判定错误。我一般会这么理解大M法是用一个人造的大数去模拟无穷大但计算机里的浮点数没有无穷大只有有限精度这个模拟从一开始就是有损的。两阶段法把可行性判断和最优性优化拆开第一阶段只关心人造变量能不能归零不涉及原始目标函数。这样做的代价是多跑一遍单纯形但换来的是数值稳定性。对于大规模问题多跑一遍的代价远小于数值不稳定导致的错误结果。工业求解器还会在第一阶段用专门的可行性松弛技术进一步降低数值风险。从代码实现的角度看两阶段法的阶段切换逻辑确实更复杂需要处理基变量映射、冗余约束删除、人造变量换出等边界情况。但这些复杂度是一次性的——写对了就对了。大M法的复杂度是参数性的——M 取多少永远需要根据问题调整换个数据集可能就要重新调。对于需要反复求解不同问题的场景两阶段法的总成本更低。这份资源把两种方法都实现了价值在于能对比着看。先跑大M法理解惩罚机制再跑两阶段法理解可行性判断的独立性最后对比两者在同一组测试数据上的迭代次数和解的质量。如果发现大M法在某些数据上结果异常大概率是 M 的取值问题换成两阶段法就能验证。从那以后我每次实现单纯形法都会先用两阶段法跑一遍确认可行性和最优解再用大M法做交叉验证。如果两者结果不一致优先信两阶段法然后回头查大M法的 M 是不是取小了。这个习惯帮我省了很多调试时间。希望帮到你。本文还有配套的精品资源点击获取