别再死记硬背单纯形表了！用Python手搓一个单纯形法求解器（附完整代码）

用Python实现单纯形法求解器告别死记硬背的表格计算在运筹学和优化算法学习中单纯形法是一个让人又爱又恨的存在——它强大到能解决各种线性规划问题但繁琐的表格计算和迭代步骤常常让学习者陷入机械记忆的困境。今天我们将用Python从头构建一个单纯形法求解器通过代码理解算法本质让抽象的数学原理变成可运行的工程实践。1. 线性规划与单纯形法基础线性规划Linear Programming是数学优化的重要分支用于在线性约束条件下寻找目标函数的极值。其标准形式可表示为maximize: c^T * x subject to: A * x b x 0单纯形法的核心思想是通过顶点迭代在可行域的顶点间移动逐步逼近最优解。与图解法相比它的优势在于系统性适用于任意维度的线性规划问题高效性平均情况下能在多项式时间内找到解灵活性能处理等式约束、不等式约束等多种形式传统手工计算需要反复绘制和更新单纯形表而我们将用Python实现这一过程的自动化。2. 算法实现的关键步骤2.1 标准化处理任何线性规划问题都需要先转化为标准形式def standard_form(c, A, b): 将线性规划转化为标准形式 # 添加松弛变量 m, n A.shape slack np.eye(m) new_A np.hstack([A, slack]) new_c np.hstack([c, np.zeros(m)]) return new_c, new_A, b转换规则目标函数统一为最大化约束条件全部转化为等式所有变量非负2.2 初始基可行解构建我们采用松弛变量法构建初始基def initialize_simplex(c, A, b): # 构建初始表格 m, n A.shape tableau np.zeros((m1, nm1)) tableau[:-1, :n] A tableau[:-1, n:nm] np.eye(m) # 松弛变量 tableau[:-1, -1] b tableau[-1, :n] -c # 目标函数行 # 基变量索引 basis list(range(n, nm)) return tableau, basis2.3 迭代优化过程核心迭代逻辑包括三个关键操作def simplex_iterate(tableau, basis): # 选择入基变量检验数最大的非基变量 entering np.argmax(tableau[-1, :-1]) # 计算比率并选择出基变量 ratios tableau[:-1, -1] / tableau[:-1, entering] leaving np.argmin(np.where(ratios 0, ratios, np.inf)) # 高斯消元枢轴运算 pivot tableau[leaving, entering] tableau[leaving, :] / pivot for i in range(len(tableau)): if i ! leaving: tableau[i, :] - tableau[i, entering] * tableau[leaving, :] # 更新基变量集合 basis[leaving] entering return tableau, basis3. 完整求解器实现整合各模块后的完整求解器import numpy as np class SimplexSolver: def __init__(self, c, A, b): self.c c self.A A self.b b self.tableau None self.basis None self.history [] # 记录迭代过程 def solve(self, max_iter100): # 标准化问题 c_std, A_std, b_std self.standard_form() # 初始化单纯形表 self.initialize_tableau(c_std, A_std, b_std) # 迭代优化 for _ in range(max_iter): if np.all(self.tableau[-1, :-1] 0): # 最优解条件 break self.tableau, self.basis self.iterate() self.history.append(self.tableau.copy()) return self.extract_solution() def standard_form(self): # 实现标准化转换 pass def initialize_tableau(self): # 构建初始表格 pass def iterate(self): # 执行单次迭代 pass def extract_solution(self): # 从最终表格提取解 x np.zeros(self.tableau.shape[1]-1) for i, var in enumerate(self.basis): if var len(x): x[var] self.tableau[i, -1] return x4. 边界情况处理完善的求解器需要处理各种特殊情况4.1 无界解检测def check_unbounded(tableau, entering): ratios tableau[:-1, -1] / tableau[:-1, entering] return np.all(ratios 0)4.2 退化处理当基变量取值为0时可能出现循环我们采用Bland规则def select_entering(tableau): # Bland规则选择下标最小的正检验数变量 for j in range(len(tableau[-1])-1): if tableau[-1, j] 0: return j return None4.3 两阶段法处理人工变量对于无明显可行解的问题def two_phase_solve(c, A, b): # 第一阶段最小化人工变量 phase1_c np.zeros(A.shape[1]) phase1_c np.hstack([phase1_c, np.ones(A.shape[0])]) # 人工变量 # 第二阶段用第一阶段结果作为初始解 pass5. 可视化迭代过程为增强理解我们可以可视化单纯形表的迭代def plot_iteration(history, iteration): plt.figure(figsize(12, 6)) sns.heatmap(history[iteration], annotTrue, fmt.2f) plt.title(fSimplex Tableau - Iteration {iteration1}) plt.show()典型输出示例Iteration 1: [[ 2. 1. 1. 0. 40.] [ 1. 3. 0. 1. 30.] [-3. -4. 0. 0. 0.]] Iteration 2: [[ 1.67 0. 1. -0.33 30. ] [ 0.33 1. 0. 0.33 10. ] [-1.67 0. 0. 1.33 40. ]] Optimal Solution: x1 18.0, x2 4.0, Value 30.06. 工程实践建议在实际应用中还需考虑数值稳定性使用QR分解等数值稳定方法处理大型稀疏矩阵预处理通过缩放和行变换改善条件数并行计算对大规模问题采用分布式计算与商业求解器集成如CPLEX、Gurobi的Python接口# 示例使用Gurobi验证结果 import gurobipy as gp model gp.Model() x model.addVars(2, namex) model.setObjective(3*x[0] 4*x[1], sensegp.GRB.MAXIMIZE) model.addConstr(2*x[0] x[1] 40) model.addConstr(x[0] 3*x[1] 30) model.optimize()通过这个项目我们不仅掌握了单纯形法的实现更重要的是理解了其背后的数学原理。代码已开源在GitHub包含完整的单元测试和示例问题读者可以自行扩展对偶单纯形法、内点法等高级功能。