用Python实现单纯形法：从数学原理到代码实战
第一次接触单纯形法时，我被那些表格和迭代规则弄得晕头转向。直到有一天，我决定用Python从头实现这个算法，才发现原来每一步迭代背后都有清晰的几何意义和代数逻辑。本文将带你用NumPy一步步构建单纯形法求解器，通过可视化观察算法如何在可行域的顶点间"跳跃"，最终找到最优解。
1. 单纯形法核心思想解析
单纯形法的美妙之处在于它将几何直觉与代数计算完美结合。想象一个多维空间中的多面体（可行域），每个顶点代表一个基本可行解。算法从一个顶点出发，沿着边移动到相邻的顶点，每次移动都使目标函数值改善，直到找到最优解。
关键几何性质：
- 最优解必定出现在顶点处（凸多面体的极值点）
- 两个顶点相邻当且仅当它们共享n-1个活跃约束
- 沿着边移动相当于交换一个基变量和非基变量
用Python表示这些概念非常直观。我们先定义线性规划的标准形式：
import numpy as np class LinearProgram: def __init__(self, c, A, b): """ 标准形式线性规划: max c^T x s.t. Ax <= b x >= 0 """ self.c = np.array(c) # 目标函数系数 self.A = np.array(A) # 约束矩阵 self.b = np.array(b) # 约束右端项 self.m, self.n = A.shape # 约束数, 变量数2. 单纯形表与初始基构建
将不等式转化为等式需要引入松弛变量。对于m个约束，我们添加m个松弛变量，形成扩展问题：
def initialize_simplex_tableau(self): # 添加松弛变量 slack_vars = np.eye(self.m) self.A_ext = np.hstack([self.A, slack_vars]) self.c_ext = np.hstack([self.c, np.zeros(self.m)]) # 构建初始单纯形表 tableau = np.zeros((self.m+1, self.n+self.m+1)) tableau[:-1, :self.n+self.m] = self.A_ext tableau[:-1, -1] = self.b tableau[-1, :self.n] = -self.c # 目标函数行 return tableau初始基通常选择松弛变量，这对应于原点(0,...,0)。表格中基变量的系数矩阵为单位阵，这是单纯形法启动的关键。
单纯形表结构示例：
| 基变量 | x1 | x2 | s1 | s2 | s3 | 解 |
|---|---|---|---|---|---|---|
| s1 | 1 | 0 | 1 | 0 | 0 | 4 |
| s2 | 0 | 2 | 0 | 1 | 0 | 12 |
| s3 | 3 | 2 | 0 | 0 | 1 | 18 |
| Z | -3 | -5 | 0 | 0 | 0 | 0 |
3. 迭代过程实现
单纯形法的每次迭代包含三个关键步骤：选择入基变量、确定出基变量、高斯消元更新表格。
def simplex_iteration(self, tableau, basis): # 最优性检验：选择负的最小的检验数 reduced_costs = tableau[-1, :-1] entering = np.argmin(reduced_costs) if reduced_costs[entering] >= 0: return True # 达到最优 # 比率测试确定出基变量 pivot_col = tableau[:-1, entering] ratios = tableau[:-1, -1] / pivot_col ratios[pivot_col <= 0] = np.inf # 忽略非正分母 leaving = np.argmin(ratios) if ratios[leaving] == np.inf: raise ValueError("问题无界") # 高斯消元 pivot = tableau[leaving, entering] tableau[leaving, :] /= pivot for i in range(self.m+1): if i != leaving: tableau[i, :] -= tableau[i, entering] * tableau[leaving, :] # 更新基 basis[leaving] = entering return False关键实现细节：
- 使用Bland规则避免循环（当有多个候选时选择下标最小的变量）
- 处理退化情况（比率为0）
- 无界解检测（所有约束系数非正）
4. 可视化迭代路径
为了直观理解算法的运行过程，我们可以将二维问题的迭代路径可视化：
import matplotlib.pyplot as plt def plot_feasible_region(A, b, vertices, path=None): # 绘制可行域和迭代路径 fig, ax = plt.subplots() x = np.linspace(0, 10, 400) 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 = np.array(vertices) ax.plot(vertices[:,0], vertices[:,1], 'ro') if path is not None: path = np.array(path) ax.plot(path[:,0], path[:,1], 'bo-') ax.set_xlim(0, max(vertices[:,0])+1) ax.set_ylim(0, max(vertices[:,1])+1) plt.legend() plt.show()对于Wyndor Glass公司的例子，你会看到算法如何从(0,0)出发，经过(0,6)最终到达最优解(2,6)。
5. 处理特殊情况的技巧
实际实现中需要考虑各种边界情况：
退化处理：
# 在比率测试中添加扰动避免循环 ratios = tableau[:-1, -1] / (pivot_col + 1e-10)两阶段法处理人工变量：
def phase_one(self): # 添加人工变量构建辅助问题 art_vars = np.eye(self.m) self.A_art = np.hstack([self.A_ext, art_vars]) c_art = np.hstack([np.zeros(self.n+self.m), np.ones(self.m)]) # 第一阶段求解 tableau = self.initialize_phase_one_tableau() basis = list(range(self.n+self.m, self.n+self.m+self.m)) while not self.simplex_iteration(tableau, basis): pass if not np.isclose(tableau[-1,-1], 0): raise ValueError("问题无可行解") return tableau, basis灵敏度分析实现：
def sensitivity_analysis(self, tableau, basis): # 计算影子价格 shadow_prices = tableau[-1, self.n:self.n+self.m] # 计算允许的变化范围 basis_inv = tableau[:-1, self.n:self.n+self.m] delta_b = np.linalg.solve(basis_inv, np.eye(self.m)) return { 'shadow_prices': shadow_prices, 'allowable_changes': delta_b }6. 完整算法实现与测试
将上述组件组合成完整的单纯形法求解器：
class SimplexSolver: def __init__(self, c, A, b): self.lp = LinearProgram(c, A, b) def solve(self): # 两阶段法处理 tableau, basis = self.lp.phase_one() # 第二阶段 tableau = self.lp.initialize_phase_two(tableau, basis) while not self.lp.simplex_iteration(tableau, basis): pass # 提取解 solution = np.zeros(self.lp.n + self.lp.m) for i, var in enumerate(basis): solution[var] = tableau[i, -1] return { 'optimal_value': tableau[-1, -1], 'solution': solution[:self.lp.n], 'basis': basis, 'tableau': tableau }测试案例：Wyndor Glass公司问题
c = [3, 5] A = [[1, 0], [0, 2], [3, 2]] b = [4, 12, 18] solver = SimplexSolver(c, A, b) result = solver.solve() print(f"最优值: {result['optimal_value']}") print(f"最优解: x1={result['solution'][0]}, x2={result['solution'][1]}")7. 性能优化与工程实践
在实际应用中，单纯形法的实现需要考虑数值稳定性和计算效率：
稀疏矩阵优化：
from scipy.sparse import csr_matrix class SparseSimplexTableau: def __init__(self, data, indices, indptr, shape): self.data = np.array(data) self.indices = np.array(indices) self.indptr = np.array(indptr) self.shape = shape def pivot(self, row, col): # 实现稀疏矩阵的高斯消元 pass修正单纯形法：
def revised_simplex(self): # 只存储基逆矩阵而非完整表格 basis_inv = np.eye(self.m) while True: # 计算简约成本 c_B = self.c_ext[basis] y = c_B @ basis_inv reduced_costs = self.c_ext - y @ self.A_ext # 选择入基变量 entering = np.argmin(reduced_costs) if reduced_costs[entering] >= -1e-8: break # 计算入基方向 d = basis_inv @ self.A_ext[:, entering] # 比率测试 x_B = basis_inv @ self.b ratios = x_B / d ratios[d <= 0] = np.inf leaving = np.argmin(ratios) # 更新基逆 E = np.eye(self.m) E[:, leaving] = -d / d[leaving] E[leaving, leaving] = 1 / d[leaving] basis_inv = E @ basis_inv # 更新基 basis[leaving] = entering实现单纯形法最令人兴奋的时刻是看到它成功找到最优解的那一刻。记得第一次我的代码正确求解出Wyndor问题时，我立刻画出了迭代路径——看到算法如何沿着可行域的边缘一步步爬升到最优顶点，那种直观的理解是单纯看数学推导无法获得的。
