别再死记硬背了！用Python手搓一个单纯形法求解器，理解每一步迭代

用Python手搓单纯形法求解器：从理论到代码的深度实践

为什么我们需要自己实现单纯形法？

在运筹学和优化领域，单纯形法就像是一把瑞士军刀——它可能不是最快的工具，但却是最可靠、最通用的线性规划求解方法之一。许多商业软件如Excel Solver、MATLAB都内置了这个算法，但作为学习者和实践者，亲自动手实现它能带来三大独特价值：

1. 透视算法本质：通过代码实现，你能真正理解检验数计算、基变换等关键步骤的数学含义
2. 调试思维训练：处理退化、循环等边界情况能大幅提升你的算法调试能力
3. 定制化扩展：自己实现的求解器可以轻松集成到特定业务场景中，不受商业软件限制

下面我们从一个实际的生产优化问题出发，逐步构建完整的单纯形法求解器。

1. 问题建模与标准化

案例：家具厂生产优化

假设一家家具厂生产桌子和椅子，每张桌子利润为3元，每把椅子利润为4元。生产受以下限制：

• 木材供应：2单位/桌 + 1单位/椅 ≤ 40
• 人工工时：1单位/桌 + 3单位/椅 ≤ 30

数学建模：

max Z = 3x₁ + 4x₂ s.t.: 2x₁ + x₂ ≤ 40 x₁ + 3x₂ ≤ 30 x₁, x₂ ≥ 0

标准化转换

单纯形法要求约束条件为等式形式，我们引入松弛变量：

# 原始约束 2x₁ + x₂ + s₁ = 40 x₁ + 3x₂ + s₂ = 30 # 目标函数调整 max Z = 3x₁ + 4x₂ + 0s₁ + 0s₂

Python实现标准化：

def standard_form(c, A, b): """将线性规划转化为标准形 参数： c: 目标函数系数 [n,] A: 约束矩阵 [m,n] b: 约束右侧常数 [m,] 返回： 标准形系数 (c', A', b') """ m, n = A.shape # 添加松弛变量 c_slack = np.zeros(m) A_slack = np.eye(m) return ( np.hstack([c, c_slack]), # 新目标系数 np.hstack([A, A_slack]), # 新约束矩阵 b.copy() # 约束右侧不变 )

2. 单纯形表数据结构设计

高效的单纯形表实现是算法核心。我们使用NumPy数组存储，并设计专门的类来管理：

class SimplexTable: def __init__(self, c, A, b): self.c = c # 目标函数系数 self.A = A # 约束矩阵 self.b = b # 右侧常数 self.basis = [...] # 当前基变量索引 self.z = 0 # 当前目标值 def compute_reduced_cost(self): """计算所有非基变量的检验数""" pass def pivot(self, entering, leaving): """执行旋转操作""" pass def get_solution(self): """提取当前解""" pass

关键数据结构示例：

3. 算法核心步骤实现

3.1 初始可行解确定

def initialize(self): # 检查是否需要两阶段法 if np.all(self.b >= 0): # 简单情况：松弛变量直接形成可行基 self.basis = list(range(len(self.c)-self.A.shape[0], len(self.c))) else: # 需要两阶段法处理 self.phase1()

3.2 检验数计算与入基选择

def select_entering(self): """选择入基变量——最大检验数规则""" reduced_costs = self.compute_reduced_cost() entering = np.argmax(reduced_costs) return entering if reduced_costs[entering] > 0 else None

3.3 θ比值测试与出基选择

def select_leaving(self, entering): """通过最小比值测试选择出基变量""" ratios = [] for i in range(len(self.b)): if self.A[i, entering] > 0: ratios.append(self.b[i] / self.A[i, entering]) else: ratios.append(np.inf) leaving = np.argmin(ratios) return leaving if ratios[leaving] < np.inf else None

3.4 基变换（旋转操作）

def pivot(self, entering, leaving): """执行旋转操作""" pivot_val = self.A[leaving, entering] # 更新主元行 self.A[leaving] /= pivot_val self.b[leaving] /= pivot_val # 更新其他行 for i in range(len(self.b)): if i != leaving and self.A[i, entering] != 0: factor = self.A[i, entering] self.A[i] -= factor * self.A[leaving] self.b[i] -= factor * self.b[leaving] # 更新基变量记录 self.basis[leaving] = entering

4. 完整算法流程

将上述步骤组合成完整算法：

def solve(self): self.initialize() while True: entering = self.select_entering() if entering is None: # 最优解条件 break leaving = self.select_leaving(entering) if leaving is None: # 无界情况 raise ValueError("Problem is unbounded") self.pivot(entering, leaving) return self.get_solution()

5. 可视化迭代过程

理解单纯形法的几何意义至关重要。我们可以用Matplotlib展示解在可行域顶点间的移动：

import matplotlib.pyplot as plt def plot_iteration(path, A, b): """绘制单纯形法的迭代路径""" fig, ax = plt.subplots() # 绘制约束条件 x = np.linspace(0, 25, 100) for i in range(A.shape[0]): y = (b[i] - A[i,0]*x) / A[i,1] ax.plot(x, y, label=f'约束{i+1}') # 绘制可行域顶点 vertices = compute_vertices(A, b) ax.scatter(*zip(*vertices), c='red') # 绘制迭代路径 path = np.array(path) ax.plot(path[:,0], path[:,1], 'bo-') ax.set_xlim(0, 25) ax.set_ylim(0, 15) ax.legend() plt.show()

6. 处理特殊情况的技巧

6.1 退化与循环

当基变量取值为0时可能出现退化，导致算法循环。我们采用Bland规则：

def select_entering_bland(self): """Bland规则：选择下标最小的正检验数变量""" reduced_costs = self.compute_reduced_cost() for j in range(len(reduced_costs)): if reduced_costs[j] > 0: return j return None

6.2 初始可行解获取

当原点不是可行解时，需要使用两阶段法：

def phase1(self): """第一阶段：构造辅助问题寻找初始可行解""" # 添加人工变量构建辅助问题 aux_c = np.zeros(self.A.shape[1]) # 原变量系数为0 aux_c = np.hstack([aux_c, np.ones(self.A.shape[0])]) # 人工变量系数为1 aux_A = np.hstack([self.A, np.eye(self.A.shape[0])]) aux_table = SimplexTable(aux_c, aux_A, self.b.copy()) # 解辅助问题 aux_table.solve() if not np.isclose(aux_table.z, 0): raise ValueError("Problem is infeasible") # 转换到原问题 self.basis = ...

7. 性能优化技巧

工业级实现需要考虑数值稳定性与计算效率：

def revised_simplex(self): """修正单纯形法：只存储基逆矩阵""" Binv = np.eye(len(self.basis)) # 基逆矩阵 while True: # 计算对偶变量 y = self.c[self.basis] @ Binv # 选择入基变量 entering = self.select_entering_revised(y) # 计算入基方向 d = Binv @ self.A[:, entering] # 选择出基变量 leaving = self.select_leaving_revised(d) # 更新基逆矩阵 Binv = self.update_basis_inverse(Binv, entering, leaving)

8. 完整代码架构

最终我们的求解器包含以下模块：

simplex_solver/ │── core/ │ ├── tableau.py # 单纯形表实现 │ ├── phase1.py # 第一阶段处理 │ └── revised.py # 修正单纯形法 │── utils/ │ ├── io.py # 模型读取/输出 │ └── visualization.py# 求解过程可视化 │── examples/ # 示例问题 └── tests/ # 单元测试

使用示例：

from simplex_solver import SimplexSolver # 构建模型 c = np.array([3, 4]) # 目标函数 A = np.array([[2, 1], [1, 3]]) # 约束矩阵 b = np.array([40, 30]) # 右侧常数 # 求解 solver = SimplexSolver(c, A, b) solution = solver.solve() print(f"最优解: x1={solution.x[0]:.2f}, x2={solution.x[1]:.2f}") print(f"最优值: {solution.z:.2f}")

9. 验证与测试

确保算法正确性的关键测试案例：

import unittest class TestSimplex(unittest.TestCase): def test_unbounded(self): c = np.array([1, 1]) A = np.array([[-1, 1], [-1, -1]]) b = np.array([-1, -2]) with self.assertRaises(ValueError): SimplexSolver(c, A, b).solve() def test_optimal(self): c = np.array([3, 2]) A = np.array([[1, 2], [1, 1]]) b = np.array([6, 4]) solver = SimplexSolver(c, A, b) self.assertAlmostEqual(solver.solve().z, 8)

10. 扩展应用

单纯形法不仅适用于标准线性规划，还可扩展应用于：

1. 整数规划松弛：作为分支定界法的底层求解器
2. 灵敏度分析：研究参数变化对解的影响
3. 对偶理论：通过单纯形乘子获取对偶解

class SensitivityAnalyzer: def __init__(self, solver): self.solver = solver def analyze(self, param_range): """分析目标系数变化对解的影响""" results = [] for delta in param_range: modified_c = self.solver.c + delta results.append(SimplexSolver(modified_c, self.solver.A, self.solver.b).solve()) return results

通过这个完整的实现过程，你不仅掌握了单纯形法的数学原理，还获得了将其转化为高效代码的实践经验。这种"通过编码理解算法"的方式，往往比单纯的理论学习更能带来深刻认知。