线性规划单纯形法、大M法与两阶段法：完整求解程序实现与避坑指南-北京尧图网络

简介这份资源面向学习最优化方法、运筹学或线性规划课程的高校学生与工程技术人员聚焦单纯形法求解中的大M法与两阶段法提供可直接运行的MATLAB程序实现。压缩包共3个文件均为.m脚本整体约4KB分别承担主程序入口、两阶段单纯形法求解以及大M法与单纯形法核心算法代码注释详尽便于对照教材理解算法流程。目前已有2547人学习下载说明其在教学与自学场景中具有一定参考价值。读者可借助主函数输入约束方程与目标函数自动调用两阶段法或大M法完成规划求解从而直观比较两种方法在人工变量处理与迭代过程上的差异同时可结合注释梳理单纯形表更新、检验数判断与最优解输出等关键环节为课程作业、算法复现与后续扩展提供可复用的脚本基础。1. 从一份“最优化方法”作业说起线性规划单纯形法到底怎么跑起来很多同学第一次接触最优化方法都是在“线性规划”这一章被单纯形法卡住的。课本上把单纯形表讲得头头是道可一到编程实现面对大 M 法和两阶段法就不知道从哪下手。这份资源就是一份完整的线性规划求解程序实现把单纯形法、大 M 法和两阶段法三种经典算法都写成了可运行的代码。它适合正在做运筹学课程设计的学生也适合需要把线性规划求解器嵌到自己项目里的开发者。拿到手之后你不需要从零推导公式直接看代码就能明白每一步迭代在干什么改几个参数就能解自己的问题。2. 项目结构与算法选型三种方法怎么分工2.1 目录结构一览这份资源解压后是一个典型的课程设计工程目录结构不复杂但该有的模块都有。常见做法是按算法拆文件再配一个统一的入口和测试用例。LinearProgramming/ ├── simplex.py # 标准单纯形法核心实现 ├── big_m.py # 大M法求解器 ├── two_phase.py # 两阶段法求解器 ├── utils.py # 矩阵运算、表格打印等工具函数 ├── main.py # 命令行入口读取问题并调度算法 ├── test_cases/ # 测试用例每个用例一个文本文件 │ ├── case1.txt │ ├── case2.txt │ └── case3.txt └── README.md # 使用说明这个划分方式的好处是职责清晰simplex.py只负责处理标准型约束全是等式、右端项非负、有初始可行基的迭代big_m.py和two_phase.py负责把一般形式转化成标准型再调用核心迭代。你如果要加对偶单纯形法只需要新增一个文件不用动已有的逻辑。2.2 三种算法的适用场景标准单纯形法要求约束已经是等式且右端项非负现实中很少直接满足。大 M 法和两阶段法就是用来处理“初始可行基不好找”的情况。方法核心思路适用场景注意点标准单纯形法从已知可行基出发迭代约束为等式且右端项非负需要手动给出初始基大 M 法引入人工变量目标函数加惩罚项 M约束含 ≥ 或 M 取值要足够大否则可能误判两阶段法第一阶段最小化人工变量之和第二阶段解原问题同上数值更稳定多一次完整迭代代码稍长我一般会优先用两阶段法因为大 M 法里那个 M 到底取多大是个玄学问题取小了可能把人工变量留在基里取大了又容易引起数值误差。两阶段法虽然多写几十行但结果更可靠。2.3 数据输入格式测试用例文件采用简单的文本格式第一行是目标函数系数第二行是约束数量后面每行是一条约束。常见做法是用空格分隔数字用、、表示不等号。# case1.txt 2 3 0 0 0 3 1 1 4 1 3 6 1 0 2第一行2 3 0 0 0表示目标函数max z 2x1 3x2后面补零是因为要转成标准型后变量数会变。第二行3表示有三条约束。这种格式的好处是直观改起来方便你甚至可以直接在文本编辑器里手写一个新用例。3. 核心模块拆解单纯形表迭代的代码实现3.1 标准单纯形法的迭代主循环simplex.py是整个项目的心脏它接收一个标准型问题返回最优解和目标函数值。核心逻辑就是不断找入基变量、出基变量直到没有正检验数为止。import numpy as np def simplex(c, A, b, basis): 标准单纯形法 c: 目标函数系数已转为求 max 的标准型 A: 约束系数矩阵形状 m x n b: 右端项向量长度 m basis: 初始基变量的下标列表长度 m 返回: (最优值, 最优解向量, 迭代次数) m, n A.shape # 构造单纯形表最后一行放检验数 tableau np.zeros((m 1, n 1)) tableau[:m, :n] A tableau[:m, n] b tableau[m, :n] -np.array(c) # 检验数行求 max 时取负 iteration 0 while True: iteration 1 # 1. 找入基变量检验数行中最小负值对应的列 reduced tableau[m, :n] if np.all(reduced -1e-9): break # 没有负检验数已达最优 enter int(np.argmin(reduced)) # 2. 找出基变量比值检验 ratios [] for i in range(m): if tableau[i, enter] 1e-9: ratios.append(tableau[i, n] / tableau[i, enter]) else: ratios.append(np.inf) if all(r np.inf for r in ratios): raise ValueError(问题无界) leave int(np.argmin(ratios)) # 3. 高斯消元把入基列变成单位向量 pivot tableau[leave, enter] tableau[leave, :] / pivot for i in range(m 1): if i ! leave: tableau[i, :] - tableau[i, enter] * tableau[leave, :] basis[leave] enter # 提取最优解 x np.zeros(n) for i in range(m): x[basis[i]] tableau[i, n] return tableau[m, n], x, iteration这段代码有几个关键点。第一检验数行取负是因为我们统一按求最大值处理如果原问题是求最小值在调用前先取反即可。第二比值检验里用1e-9而不是零是为了避免浮点误差导致除零。第三入基变量选的是最小负检验数对应的列这是 Dantzig 规则虽然收敛慢但实现简单。如果你要改成 Bland 规则避免循环把argmin换成第一个负检验数的下标就行。参数basis是初始基变量的下标列表比如[3, 4, 5]表示松弛变量x4, x5, x6构成初始基。这个列表在迭代中会被原地修改所以调用方要注意是否需要保留原始副本。3.2 大 M 法的惩罚项处理大 M 法的核心是在目标函数里给人工变量加上-M的系数求 max 时迫使它们在最优解中为零。big_m.py负责构造这个扩展问题。def build_big_m_problem(c, A, b, constraints): 将一般线性规划转为大M法标准型 c: 原始目标函数系数 A: 原始约束系数矩阵 b: 原始右端项 constraints: 每条约束的类型列表如 [, , ] 返回: (新目标系数, 新系数矩阵, 新右端项, 初始基下标) m, n A.shape new_c list(c) new_A [list(row) for row in A] new_b list(b) basis [] artificial_indices [] for i, typ in enumerate(constraints): if typ : # 加松弛变量系数为 1 for j in range(m): new_A[j].append(1 if j i else 0) new_c.append(0) basis.append(len(new_c) - 1) elif typ : # 减剩余变量系数为 -1 for j in range(m): new_A[j].append(-1 if j i else 0) new_c.append(0) # 再加人工变量系数为 1 for j in range(m): new_A[j].append(1 if j i else 0) new_c.append(-1e6) # 大M惩罚 basis.append(len(new_c) - 1) artificial_indices.append(len(new_c) - 1) else: # # 直接加人工变量 for j in range(m): new_A[j].append(1 if j i else 0) new_c.append(-1e6) basis.append(len(new_c) - 1) artificial_indices.append(len(new_c) - 1) return new_c, np.array(new_A), np.array(new_b), basis, artificial_indices这里M取的是1e6是个经验值。如果问题系数本身很大这个 M 可能不够大导致人工变量没被赶出去。我一般会在求解后检查一下人工变量是否为零如果不为零就提示用户调大 M。另外注意型约束先减剩余变量再加人工变量顺序不能反否则初始基的单位向量构造不出来。3.3 两阶段法的两轮求解两阶段法把问题拆成两个标准型第一阶段最小化人工变量之和第二阶段用第一阶段得到的基去解原问题。two_phase.py的代码结构比大 M 法清晰因为不需要纠结 M 的大小。def two_phase_solve(c, A, b, constraints): 两阶段法求解一般线性规划返回: (最优值, 最优解, 总迭代次数) # 第一阶段构造辅助问题目标是最小化人工变量之和 m, n A.shape c1 [0] * n A1 [list(row) for row in A] basis [] artificial [] for i, typ in enumerate(constraints): if typ : for j in range(m): A1[j].append(1 if j i else 0) c1.append(0) basis.append(len(c1) - 1) elif typ : for j in range(m): A1[j].append(-1 if j i else 0) c1.append(0) for j in range(m): A1[j].append(1 if j i else 0) c1.append(1) # 第一阶段目标最小化人工变量 basis.append(len(c1) - 1) artificial.append(len(c1) - 1) else: for j in range(m): A1[j].append(1 if j i else 0) c1.append(1) basis.append(len(c1) - 1) artificial.append(len(c1) - 1) # 第一阶段求 min转成 max 取负 val1, x1, iter1 simplex([-v for v in c1], np.array(A1), np.array(b), basis[:]) if abs(val1) 1e-6: raise ValueError(问题无可行解) # 第二阶段去掉人工变量列用第一阶段的基继续迭代 keep [j for j in range(len(c1)) if j not in artificial] A2 np.array(A1)[:, keep] c2 list(c) [0] * (len(keep) - n) # 重新映射基下标 new_basis [] for b_idx in basis: if b_idx in keep: new_basis.append(keep.index(b_idx)) else: # 人工变量在基里但值为零需要换基这里简化处理 new_basis.append(0) val2, x2, iter2 simplex(c2, A2, np.array(b), new_basis) return val2, x2[:n], iter1 iter2第一阶段的目标函数系数里人工变量为 1其余为 0然后取负转成求 max。如果第一阶段最优值不为零说明原问题无可行解。第二阶段要去掉人工变量列同时把基下标重新映射到新矩阵的列上。这里有个坑如果人工变量在第一阶段结束时仍在基里但值为零直接删列会导致基不可用常见做法是找一个非人工变量替换它或者保留该列但强制其值为零。代码里做了简化处理实际使用时建议加一个换基步骤。4. 避坑与排查单纯形法实现中的常见翻车点4.1 迭代死循环现象程序跑了几千次迭代还没停检验数一直在负值和零之间跳。原因单纯形法在退化情况下可能循环Dantzig 规则不保证有限终止。解决改用 Bland 规则即入基选第一个负检验数对应的变量出基选比值相同情况下下标最小的。虽然收敛慢但能保证不循环。4.2 大 M 法结果里人工变量没消掉现象最优解里人工变量不为零目标函数值明显偏小。原因M 取值不够大或者问题本身无可行解。解决先检查第一阶段是否可行如果可行再调大 M比如从1e6调到1e8。但要注意 M 太大会导致数值不稳定两阶段法是更稳妥的选择。4.3 浮点误差导致检验数判断错误现象明明应该继续迭代程序却认为已经最优或者相反。原因浮点运算累积误差检验数在零附近震荡。解决所有和零比较的地方都加容差比如 -1e-9而不是 0。比值检验里分母也要判断是否大于1e-9。4.4 输入格式解析失败现象读取测试用例时报ValueError或数组形状不匹配。原因文本文件里有多余空行、制表符或者不等号写成了而不是。解决在main.py里加一层清洗用line.strip().split()统一处理空白不等号只取第一个字符。我一般会先打印解析后的矩阵确认一遍再求解。4.5 无界问题误判为有解现象程序返回了一个很大的目标函数值但实际问题是发散的。原因入基列所有系数都小于等于零时比值检验全为无穷应该报无界。解决在比值检验后加判断如果所有比值都是inf直接抛出“问题无界”异常不要继续迭代。5. 二次开发与扩展从能跑到好用5.1 加一个灵敏度分析模块单纯形法解出最优解只是第一步实际项目里经常要看目标系数或右端项变化时最优解怎么变。你可以在simplex.py返回时把最终的单纯形表也带出来然后写一个sensitivity.py计算影子价格和允许变化范围。def sensitivity_analysis(tableau, basis, m, n): 从最终单纯形表提取灵敏度信息 tableau: 最终单纯形表 basis: 最终基变量下标返回: 影子价格列表、目标系数允许范围 shadow_prices tableau[m, n:] # 松弛变量对应的检验数 ranges [] for i in range(m): # 简化处理只输出当前基变量的系数 ranges.append({ var: basis[i], value: tableau[i, n], coefficient: tableau[m, basis[i]] }) return shadow_prices, ranges影子价格就是松弛变量在最终表中的检验数经济含义是约束右端项每增加一个单位目标函数的增量。这个功能对做资源分配类问题的同学特别有用。5.2 换成 Bland 规则避免循环把simplex函数里的入基选择改一行就行# 原代码enter int(np.argmin(reduced)) # 改为 Bland 规则 enter None for j in range(n): if reduced[j] -1e-9: enter j break出基选择也要相应调整在比值相同的情况下选下标最小的基变量。这个改动会让迭代次数变多但再也不会出现死循环。我一般会在求解前加一个开关默认用 Dantzig如果迭代超过 1000 次自动切到 Bland。5.3 用对偶单纯形法处理新增约束实际项目里经常遇到“已经解好了突然加一条约束”的情况。从头再解一遍太浪费对偶单纯形法可以从原最优基出发快速调整到新最优。实现思路是把新约束加进最终表如果右端项为负就用对偶单纯形法迭代。这部分代码可以单独放一个dual_simplex.py复用已有的高斯消元逻辑。5.4 验证方法用标准测试集对拍改完代码怎么确认没改坏我习惯准备一组已知答案的小规模问题每次改动后跑一遍对拍。比如经典的 Beale 问题、Klee-Minty 立方体问题这些都有公开的最优解。另外可以用scipy.optimize.linprog做交叉验证虽然它用的是内点法但最优值应该一致。from scipy.optimize import linprog def verify_with_scipy(c, A, b): 用 scipy 对拍验证 res linprog(c[-v for v in c], A_ubA, b_ubb, methodhighs) return -res.fun if res.success else None注意linprog默认求最小值所以目标系数要取负。如果两个结果差在1e-6以内基本可以认为实现正确。5.5 一个我踩过的坑有一次我把两阶段法第二阶段的目标系数直接用了原始c忘了补松弛变量对应的零结果程序跑出来最优解完全不对。排查了半天才发现是数组长度不匹配simplex函数里c和A的列数对不上但 numpy 广播没报错只是静默算错了。从那以后我每次调用simplex之前都强制走一遍assert len(c) A.shape[1]这个习惯帮我省了很多后悔药。希望你在改代码的时候也加上类似的检查别让静默错误浪费一晚上。本文还有配套的精品资源点击获取

线性规划单纯形法、大M法与两阶段法：完整求解程序实现与避坑指南

相关推荐

JavaScript核心教程：从变量、DOM到异步编程的实战指南

Azure图标集实战：SVG命名、四种接入、主题与避坑体检

Python requests响应处理全指南：正文、状态码、Cookie、耗时、响应头

更多资讯

PHP+XML生成家谱树：从模型到渲染的完整实践

阿里Coding Plan实测：OpenClaw接入TaoToken的API Key配置指南

《孤胆枪手2026》安装与性能优化指南

Java SE教学资源管理系统：JDBC+Swing毕设开发全解析

Java Swing人事管理系统源码：Swing+JDBC实现桌面端CRUD全链路

E430刷BIOS 2.5：版本确认、U盘离线刷写与避坑

今日更新

Verilog实现FPGA密码锁：矩阵键盘与状态机设计实战

数据服务监控运维实战：从指标体系到告警治理与故障排查

FreeRTOS任务通知实战详解：内核原理、API与避坑指南

本周精选

金融AI落地指南：从调研报告到技术选型与工程实践

大模型上下文成本与注意力优化技术深度解析——为何O(n²)平方复杂度无法被常规优化消除｜TaoToken 统一 Key 通道实测

基于GMM的语音识别系统原理与MATLAB实现详解

本月热门

我发现了一个新思路：用 Remotion + Claude Code 像写代码一样自动化生成短视频

Windows下 Codex 中 Chrome 和 Computer Use 插件不可用问题排查及解决参考方式：TaoToken 统一 Key 配置与验证

2026 大模型集体涨价：用 Python 做企业 Token 成本测算与选型避坑（附配置）

读完文章，想聊聊您的网站？