用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.m1, self.nself.m1)) tableau[:-1, :self.nself.m] self.A_ext tableau[:-1, -1] self.b tableau[-1, :self.n] -self.c # 目标函数行 return tableau初始基通常选择松弛变量这对应于原点(0,...,0)。表格中基变量的系数矩阵为单位阵这是单纯形法启动的关键。单纯形表结构示例基变量x1x2s1s2s3解s1101004s20201012s33200118Z-3-500003. 迭代过程实现单纯形法的每次迭代包含三个关键步骤选择入基变量、确定出基变量、高斯消元更新表格。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.m1): 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, pathNone): # 绘制可行域和迭代路径 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, labelf约束{i1}) # 标记顶点和路径 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.nself.m), np.ones(self.m)]) # 第一阶段求解 tableau self.initialize_phase_one_tableau() basis list(range(self.nself.m, self.nself.mself.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.nself.m] # 计算允许的变化范围 basis_inv tableau[:-1, self.n:self.nself.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问题时我立刻画出了迭代路径——看到算法如何沿着可行域的边缘一步步爬升到最优顶点那种直观的理解是单纯看数学推导无法获得的。
别再死记硬背表格了！用Python手把手实现单纯形法，理解每一步迭代的底层逻辑
用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.m1, self.nself.m1)) tableau[:-1, :self.nself.m] self.A_ext tableau[:-1, -1] self.b tableau[-1, :self.n] -self.c # 目标函数行 return tableau初始基通常选择松弛变量这对应于原点(0,...,0)。表格中基变量的系数矩阵为单位阵这是单纯形法启动的关键。单纯形表结构示例基变量x1x2s1s2s3解s1101004s20201012s33200118Z-3-500003. 迭代过程实现单纯形法的每次迭代包含三个关键步骤选择入基变量、确定出基变量、高斯消元更新表格。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.m1): 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, pathNone): # 绘制可行域和迭代路径 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, labelf约束{i1}) # 标记顶点和路径 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.nself.m), np.ones(self.m)]) # 第一阶段求解 tableau self.initialize_phase_one_tableau() basis list(range(self.nself.m, self.nself.mself.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.nself.m] # 计算允许的变化范围 basis_inv tableau[:-1, self.n:self.nself.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问题时我立刻画出了迭代路径——看到算法如何沿着可行域的边缘一步步爬升到最优顶点那种直观的理解是单纯看数学推导无法获得的。
