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

用Python从零构建单纯形法求解器原理剖析与工业级实现运筹优化领域中线性规划问题如同夜空中的北极星为复杂决策问题指引方向。而单纯形法则是寻找这颗北极星最经典的罗盘——它不仅是教科书中的算法典范更是工业级优化系统的核心引擎。本文将带您深入算法内核用Python实现一个支持完整迭代可视化与边界条件处理的单纯形法求解器。1. 单纯形法的核心思想与数学表达单纯形法的本质是在多维空间中沿着可行域的边进行智能爬山。想象一个多面体每个顶点代表一个基本可行解算法的工作就是从一个顶点出发沿着使目标函数增长最快的边移动到相邻顶点。算法数学表述 给定标准型线性规划问题maximize cᵀx subject to Ax ≤ b, x ≥ 0其中x ∈ ℝⁿ是决策变量向量A ∈ ℝ^(m×n)是约束矩阵b ∈ ℝ^m是资源向量c ∈ ℝⁿ是目标函数系数通过引入松弛变量s将不等式转化为等式Ax Is b x, s ≥ 0此时基变量由初始的松弛变量构成对应初始可行解x0。2. Python实现框架设计我们采用面向对象方式构建求解器核心类结构如下class SimplexSolver: def __init__(self, c, A, b): self.c np.array(c) # 目标函数系数 self.A np.array(A) # 约束矩阵 self.b np.array(b) # 右端项 self.m, self.n A.shape # 约束数、决策变量数 self.tableau None # 单纯形表 self.basis None # 当前基变量索引 self.history [] # 迭代历史记录 def _initialize(self): 构建初始单纯形表 self.tableau np.zeros((self.m1, self.nself.m1)) self.tableau[:-1, :self.n] self.A self.tableau[:-1, self.n:self.nself.m] np.eye(self.m) self.tableau[:-1, -1] self.b self.tableau[-1, :self.n] -self.c self.basis list(range(self.n, self.nself.m)) def solve(self, max_iter100): self._initialize() for _ in range(max_iter): if self._check_optimality(): break entering self._select_entering() leaving self._select_leaving(entering) self._pivot(entering, leaving) self.history.append(self._current_solution()) return self._get_solution()关键数据结构说明属性类型描述tableaunumpy.ndarray(m1)×(nm1)的单纯形表矩阵basisList[int]当前基变量在表中的列索引集合historyList[dict]记录每次迭代后的基变量和目标值3. 迭代过程的核心操作3.1 最优性检验与入基变量选择最优性检验通过检查目标行系数实现def _check_optimality(self): 检查当前解是否最优 return all(x 0 for x in self.tableau[-1, :-1]) def _select_entering(self): 选择入基变量Bland规则避免循环 for j in range(self.n self.m): if self.tableau[-1, j] 0 and j not in self.basis: return j return None提示实际工业实现中常采用Dantzig规则选择负系数最大的变量入基但Bland规则能保证算法终止3.2 出基变量选择与主元运算出基选择通过最小比值测试实现def _select_leaving(self, entering): 执行最小比值测试选择出基变量 ratios [] for i in range(self.m): if self.tableau[i, entering] 0: ratio self.tableau[i, -1] / self.tableau[i, entering] ratios.append((ratio, i)) if not ratios: raise ValueError(Problem is unbounded) return min(ratios)[1] def _pivot(self, entering, leaving): 执行主元运算 pivot_val self.tableau[leaving, entering] self.tableau[leaving, :] / pivot_val for i in range(self.m 1): if i ! leaving: self.tableau[i, :] - self.tableau[i, entering] * self.tableau[leaving, :] self.basis[leaving] entering4. 边界条件处理与算法强化工业级实现需要处理各种边界情况4.1 退化情形检测当最小比值测试出现多个相同值时算法进入退化状态def _check_degeneracy(self): 检测当前解是否退化 basic_vars self.tableau[:-1, -1] return any(abs(x) 1e-10 for x in basic_vars)4.2 两阶段法处理人工变量对于非标准形式问题需要引入两阶段法class TwoPhaseSimplex(SimplexSolver): def __init__(self, c, A, b): super().__init__(c, A, b) self.artificial [] # 记录人工变量索引 def _add_artificial_vars(self): 添加人工变量构建辅助问题 num_artificial sum(1 for bi in self.b if bi 0) self.artificial list(range(self.nself.m, self.nself.mnum_artificial)) # 扩展单纯形表... def phase1(self): 第一阶段最小化人工变量和 self._add_artificial_vars() # 修改目标函数为最小化人工变量 # 执行标准单纯形法... if abs(self.tableau[-1, -1]) 1e-6: raise ValueError(Infeasible problem) def phase2(self): 第二阶段原始问题求解 # 移除人工变量相关列 # 恢复原始目标函数 super().solve()5. 可视化迭代过程使用Matplotlib实现迭代过程动画def plot_iteration(self, iteration): 绘制指定迭代时的单纯形表状态 fig, ax plt.subplots(figsize(12, 6)) data self.history[iteration][tableau] columns [fx{i} for i in range(self.n)] [fs{i} for i in range(self.m)] [RHS] rows [f约束{i} for i in range(self.m)] [目标] ax.axis(off) table ax.table(cellTextnp.round(data, 2), colLabelscolumns, rowLabelsrows, loccenter, cellLoccenter) # 高亮显示基变量 for j in self.history[iteration][basis]: table[(0, j)].set_facecolor(#FFDDC1) plt.title(f迭代 {iteration1}: 目标值 {self.history[iteration][objective]:.2f}) plt.show()典型迭代过程输出示例迭代基变量目标值备注0s1, s20初始解1x1, s212第一次旋转2x1, x236达到最优解6. 性能优化与工业实践6.1 稀疏矩阵处理大规模问题中约束矩阵通常非常稀疏from scipy.sparse import csr_matrix class SparseSimplexSolver(SimplexSolver): def _initialize(self): 使用稀疏矩阵初始化 self.A_sparse csr_matrix(self.A) # 构建稀疏形式的单纯形表... def _pivot(self, entering, leaving): 稀疏矩阵的旋转操作 # 专用稀疏矩阵运算...6.2 数值稳定性增强添加扰动处理防止数值误差累积def _pivot(self, entering, leaving): 带数值稳定的主元运算 pivot_val self.tableau[leaving, entering] if abs(pivot_val) 1e-10: # 添加随机扰动打破退化 self.tableau[leaving, entering] np.random.uniform(1e-8, 1e-7) pivot_val self.tableau[leaving, entering] # 继续正常旋转...7. 完整代码实现与测试案例完整求解器代码包含以下核心文件结构simplex_solver/ ├── core.py # 基础单纯形法实现 ├── two_phase.py # 两阶段法扩展 ├── sparse.py # 稀疏矩阵版本 ├── viz.py # 可视化工具 └── tests/ # 测试案例 ├── test_wyndor.py # 经典Wyndor案例 └── test_degen.py # 退化情形测试示例测试案例Wyndor Glass问题def test_wyndor(): Wyndor Glass公司案例测试 c [3, 5] # 目标函数系数 A [[1, 0], # 约束矩阵 [0, 2], [3, 2]] b [4, 12, 18] # 右端项 solver SimplexSolver(c, A, b) solution solver.solve() assert np.allclose(solution[x], [2, 6]) assert abs(solution[objective] - 36) 1e-6 print(测试通过最优解, solution)在实现过程中有几个关键点需要特别注意主元选择时的数值精度处理循环检测与预防机制内存效率与大规模问题支持与其他优化库如SciPy的结果对比验证