# 线性规划实战：两阶段单纯形法从入门到精通（附Python代码示例） 线性规划是优化领域的一块基石，它试图在满足一系列线性等式或不等式约束的条件下，最大化或最小化一个线性目标函数。单纯形法作为求解线性规划问题的经典算法，其核心思想是在可行域的顶点（基本可行解）之间“行走”，直到找到最优解。然而，对于许多实际问题，我们很难直接找到一个初始的顶点作为起点。这就好比在一片未知的荒野中寻找最高点，你不仅需要一张地图（目标函数），还需要一个明确的出发营地（初始基本可行解）。两阶段单纯形法正是为了解决这个“营地搭建”问题而诞生的精巧策略。 本文将从一个实践者的视角，深入剖析两阶段单纯形法的原理、步骤与实现细节。我们不仅会探讨其背后的数学逻辑，更会通过完整的Python代码实现，将理论转化为可运行、可调试的工具。无论你是希望巩固算法理解的学生，还是需要在项目中应用线性规划的技术人员，这篇文章都将带你从“知道”走向“精通”，并为你提供一套可直接复用的代码框架。 ## 1. 理解两阶段法的必要性：从“无家可归”到“安营扎寨” 在标准的单纯形法中，我们从一个已知的基本可行解开始迭代。但现实中的线性规划模型，其约束条件往往不直接提供一个现成的、所有变量非负的“单位矩阵”作为初始基。例如，考虑以下问题： **最大化：** Z = 3x₁ + 5x₂ **约束条件：** 1. x₁ + x₂ ≥ 4 2. 2x₁ + x₂ ≤ 10 3. x₁, x₂ ≥ 0 为了应用单纯形法，我们首先需要将其转化为标准形式（等式约束，右端项非负）。引入松弛变量s₁和剩余变量s₂后，约束变为： 1. x₁ + x₂ - s₁ = 4 (s₁ ≥ 0) 2. 2x₁ + x₂ + s₂ = 10 (s₂ ≥ 0) 现在，我们观察系数矩阵。第二个约束中，s₂的系数是1，可以作为一个初始基变量。但第一个约束中，没有现成的、系数为1且只出现在该约束中的变量（x₁, x₂, s₁的系数都不是[1,0]或[0,1]这样的单位向量）。我们无法直接找到一个由松弛变量构成的初始基。这就是所谓的“**初始基本可行解不易获得**”的典型情况。 > 注意：当约束是“≤”且右端项非负时，引入的松弛变量天然构成初始基。但“≥”或“=”约束破坏了这一便利。 为了解决这个问题，两阶段法引入了一个巧妙的“脚手架”——**人工变量**。其核心思想分为两个阶段： * **第一阶段（搭建脚手架）**：构造一个辅助的线性规划问题，其目标是**最小化所有人工变量的和**。通过求解这个辅助问题，我们判断原问题是否可行（目标和是否为0），并同时为原问题“雕刻”出一个初始的基本可行解。 * **第二阶段（拆除脚手架，建造主体）**：利用第一阶段找到的可行解作为起点，移除非必需的人工变量，转而求解原始的目标函数，直至找到最优解。 这种方法将寻找起点的困难，转化为求解另一个更简单的线性规划问题，在理论上保证了算法的完备性。 ## 2. 第一阶段：构造与求解辅助问题 第一阶段的目标是判断可行性并找到一个初始顶点。让我们用数学语言和代码来具体化这个过程。 ### 2.1 引入人工变量与构造辅助问题 对于每个没有明显初始基变量的约束，我们引入一个非负的**人工变量**。以上述问题为例，第一个约束需要引入人工变量a₁。 原约束系统变为： 1. x₁ + x₂ - s₁ + a₁ = 4 2. 2x₁ + x₂ + s₂ = 10 3. 所有变量 ≥ 0 现在，我们可以轻松地构造一个初始基：令基变量为 `[a₁, s₂]`，非基变量 `[x₁, x₂, s₁]` 设为0。这给出了一个初始基本解：`a₁=4, s₂=10, 其他=0`。这个解满足所有约束，但它不是原问题的可行解，因为包含了本不该存在的人工变量a₁。 为了“驱逐”这些人工变量，我们构造一个**辅助问题（Phase I Problem）**： * **目标函数**：最小化 W = a₁ + ... + a_m (所有人工变量之和)。我们希望将它们降到0。 * **约束条件**：包含人工变量后的新等式约束。 * **变量**：包含所有原变量、松弛/剩余变量以及人工变量。 如果原问题是可行的，那么必然存在一个解使得所有人工变量都为0，从而辅助问题的最优目标值 `W* = 0`。如果 `W* > 0`，则说明我们无法在不违反原约束的情况下将所有人工变量归零，即原问题**无可行解**。 ### 2.2 Python实现：第一阶段单纯形表初始化 我们将使用经典的单纯形表形式来实现。首先，定义问题的数据结构。 ```python import numpy as np class TwoPhaseSimplex: def __init__(self, c, A, b, constraint_types): """ 初始化线性规划问题。 参数: c: 目标函数系数向量 (n,)， 最大化目标。 A: 约束系数矩阵 (m, n) b: 右端项向量 (m,)，应已转换为非负。 constraint_types: 约束类型列表，元素为 '<='， '>='， 或 '=' """ self.c = np.array(c, dtype=float) self.A_orig = np.array(A, dtype=float) self.b = np.array(b, dtype=float) self.constraint_types = constraint_types self.m, self.n = A.shape self.num_artificial = 0 self.artificial_indices = [] ``` 接下来，是标准化的关键函数，它负责添加松弛变量、剩余变量和人工变量，并构建第一阶段的单纯形表。 ```python def _standardize(self): """将问题转化为标准形式，并构建第一阶段单纯形表.""" # 首先，处理 '<=' 和 '>=' 约束，添加松弛/剩余变量 slack_surplus_vars = [] current_var_count = self.n # 创建扩展的系数矩阵和成本向量 A_extended = self.A_orig.copy() c_extended = list(self.c) + [0] * self.m # 为松弛/剩余变量预留空间，初始成本为0 phase1_c = [] # 第一阶段目标函数系数 # 第一阶段：原变量和松弛/剩余变量的成本为0 phase1_c = [0] * (self.n + self.m) for i, constr_type in enumerate(self.constraint_types): if constr_type == '<=': # 添加松弛变量，系数为1 new_col = np.zeros((self.m, 1)) new_col[i, 0] = 1 A_extended = np.hstack([A_extended, new_col]) slack_surplus_vars.append(current_var_count) current_var_count += 1 # 第一阶段成本为0 phase1_c.append(0) elif constr_type == '>=': # 添加剩余变量，系数为-1 new_col = np.zeros((self.m, 1)) new_col[i, 0] = -1 A_extended = np.hstack([A_extended, new_col]) slack_surplus_vars.append(current_var_count) current_var_count += 1 # 第一阶段成本为0 phase1_c.append(0) # 对于'>='约束，还需要添加人工变量 art_col = np.zeros((self.m, 1)) art_col[i, 0] = 1 A_extended = np.hstack([A_extended, art_col]) self.artificial_indices.append(current_var_count) phase1_c.append(1) # 人工变量在第一阶段成本为1 self.num_artificial += 1 current_var_count += 1 elif constr_type == '=': # 等号约束直接添加人工变量 art_col = np.zeros((self.m, 1)) art_col[i, 0] = 1 A_extended = np.hstack([A_extended, art_col]) self.artificial_indices.append(current_var_count) phase1_c.append(1) # 人工变量在第一阶段成本为1 self.num_artificial += 1 current_var_count += 1 # 构建第一阶段单纯形表 # 表结构: [A_extended | b] # [phase1_c | 0] (成本行) self.tableau = np.zeros((self.m + 1, current_var_count + 1)) self.tableau[:-1, :-1] = A_extended self.tableau[:-1, -1] = self.b self.tableau[-1, :-1] = phase1_c # 初始基变量：人工变量和松弛变量（来自'<='约束） self.basic_vars = [] for i, constr_type in enumerate(self.constraint_types): if constr_type in ['>=', '=']: # 基变量是人工变量 self.basic_vars.append(self.artificial_indices.pop(0)) else: # '<=' 约束 # 基变量是松弛变量，找到它 self.basic_vars.append(slack_surplus_vars.pop(0)) # 由于成本行对应的是人工变量（成本1），而基变量是人工变量，需要将成本行中基变量对应的列消元为0 self._update_cost_row_for_basis(phase1_c) ``` `_update_cost_row_for_basis` 函数用于确保单纯形表成本行中，所有基变量对应的列系数为0，这是单纯形表的标准形式。 ```python def _update_cost_row_for_basis(self, original_cost): """更新成本行，使得基变量对应的列为0。""" for basic_var in self.basic_vars: if original_cost[basic_var] != 0: # 找到基变量所在的行 row_idx = self.basic_vars.index(basic_var) pivot_row = self.tableau[row_idx, :] # 成本行减去 pivot_row * cost_coeff cost_coeff = original_cost[basic_var] self.tableau[-1, :] -= cost_coeff * pivot_row ``` ### 2.3 迭代求解与可行性判断 构建好第一阶段单纯形表后，我们使用标准的单纯形法进行迭代，但目标是最小化人工变量之和。迭代步骤包括选择进基变量（检验数为正，因为我们是最小化问题，通常我们转化为最大化 -W 来处理，但这里我们直接处理最小化）、选择出基变量（最小比值检验）和主元旋转。 ```python def _simplex_iteration(self, phase): """执行一次单纯形迭代。phase=1 或 2。""" # 对于第一阶段（最小化），我们看成本行（最后一行）是否有正数（因为标准单纯形表针对最大化）。 # 通常处理方式是：最小化问题转化为最大化其负值。在我们的表结构中，最后一行存储的是 -z + ... = 0 的系数。 # 因此，对于最小化，如果最后一行（成本行）在非基变量列有正数，则目标函数值还能下降。 cost_row = self.tableau[-1, :-1] if phase == 1: # 第一阶段：最小化人工变量和。检查是否有正检验数（忽略人工变量列？） # 更稳健的做法：只考虑非基变量中的非人工变量 non_basic = [j for j in range(self.tableau.shape[1] - 1) if j not in self.basic_vars] # 在非基变量中，找出检验数（成本行系数）最大的正数（因为我们是最大化 -W，等价于最小化 W） entering_candidates = [(cost_row[j], j) for j in non_basic if cost_row[j] > 1e-10] else: # 第二阶段：最大化原目标。检查是否有正检验数（标准最大化问题） entering_candidates = [(cost_row[j], j) for j in range(self.tableau.shape[1] - 1) if cost_row[j] > 1e-10 and j not in self.basic_vars] if not entering_candidates: return True # 达到最优 # 选择进基变量：最大检验数规则（Bland规则可防止循环，此处简化） entering_var = max(entering_candidates, key=lambda x: x[0])[1] # 最小比值检验确定出基变量 ratios = [] for i in range(self.m): if self.tableau[i, entering_var] > 1e-10: # 避免除零或负数 ratio = self.tableau[i, -1] / self.tableau[i, entering_var] ratios.append((ratio, i)) if not ratios: print("问题无界！") return None # 无界 leaving_row = min(ratios, key=lambda x: x[0])[1] leaving_var = self.basic_vars[leaving_row] # 主元旋转 pivot = self.tableau[leaving_row, entering_var] self.tableau[leaving_row, :] /= pivot for i in range(self.m + 1): if i != leaving_row: factor = self.tableau[i, entering_var] self.tableau[i, :] -= factor * self.tableau[leaving_row, :] # 更新基变量记录 self.basic_vars[leaving_row] = entering_var return False # 未达到最优，继续迭代 ``` 第一阶段的主循环会反复调用 `_simplex_iteration(phase=1)`，直到达到最优。然后，我们检查最优值（单纯形表右下角的值，即 -W*）。如果 `W*` 非常接近0（考虑数值误差），则说明原问题可行，可以进入第二阶段；否则，报告无可行解。 ## 3. 第二阶段：求解原始目标函数 第一阶段成功后，我们得到了一个不含人工变量（或人工变量值为0）的基本可行解。现在，我们需要“切换赛道”，开始优化原始的目标函数。 ### 3.1 过渡与表格重构 首先，我们需要从单纯形表中移除人工变量（如果它们还在的话），并将最后一行（成本行）替换为原始目标函数的系数。但更优雅的做法是，直接在当前表的基础上，将成本行更新为原目标函数的系数，并再次利用 `_update_cost_row_for_basis` 函数，确保基变量对应的成本列系数为0。 ```python def _phase2_setup(self): """准备第二阶段：将成本行替换为原目标函数系数，并重新计算.""" # 创建扩展的原目标函数成本向量，长度与当前表的变量数一致 extended_c = np.zeros(self.tableau.shape[1] - 1) # 前n个是原变量 extended_c[:self.n] = self.c # 松弛/剩余变量的成本为0，人工变量的成本也为0（但它们应该被踢出基） # 注意：人工变量可能还在表中，但我们不关心它们的成本。 # 将表的最后一行替换为扩展的成本向量（但符号取反，因为我们的表是 -z + c^T x = 0 的形式） self.tableau[-1, :-1] = -extended_c # 单纯形表标准形式是 -z + c x = 0 self.tableau[-1, -1] = 0 # 现在，需要将成本行中基变量对应的列消元为0 # 我们需要知道当前基变量对应的原成本系数 basis_cost_coeffs = [] for basic_var in self.basic_vars: if basic_var < len(extended_c): basis_cost_coeffs.append(extended_c[basic_var]) else: # 如果是人工变量（理论上第一阶段后值应为0，且可能已是非基变量），成本为0 basis_cost_coeffs.append(0.0) # 手动进行行变换，将基变量列的成本消为0 for i, basic_var in enumerate(self.basic_vars): cost_coeff = basis_cost_coeffs[i] if abs(cost_coeff) > 1e-10: self.tableau[-1, :] -= cost_coeff * self.tableau[i, :] ``` ### 3.2 第二阶段迭代与结果提取 设置好第二阶段表格后，我们再次进入单纯形迭代循环，但这次调用 `_simplex_iteration(phase=2)`。迭代结束后，单纯形表就包含了原问题的最优解信息。 我们需要从最终的单纯形表中提取解和最优值。基变量的值等于右端项（b列）对应行的值。非基变量的值为0。最优目标值位于表格的右下角（对于我们的表结构，它是 `-z` 的值，所以需要取反）。 ```python def solve(self): """执行两阶段单纯形法求解。""" print("=== 第一阶段开始 ===") self._standardize() self._print_tableau("初始第一阶段表") iteration = 0 while iteration < 1000: # 防止无限循环 optimal = self._simplex_iteration(phase=1) if optimal is None: print("第一阶段：问题无界（异常情况）。") return None if optimal: break iteration += 1 # self._print_tableau(f"第一阶段迭代 {iteration} 后") # 调试用 W_star = -self.tableau[-1, -1] # 第一阶段目标值 W* print(f"第一阶段结束。W* = {W_star:.6f}") if abs(W_star) > 1e-6: print(f"原问题无可行解 (W* = {W_star} > 0).") return None print("\n=== 第二阶段开始 ===") self._phase2_setup() self._print_tableau("初始第二阶段表") iteration = 0 while iteration < 1000: optimal = self._simplex_iteration(phase=2) if optimal is None: print("第二阶段：原问题无界。") return None if optimal: break iteration += 1 # 提取结果 solution = np.zeros(self.tableau.shape[1] - 1) for i, var_idx in enumerate(self.basic_vars): if var_idx < len(solution): solution[var_idx] = self.tableau[i, -1] opt_value = -self.tableau[-1, -1] # 原问题最优值 print(f"\n求解完成。最优值: {opt_value:.4f}") print("最优解 (原变量):", solution[:self.n]) return opt_value, solution[:self.n] ``` ## 4. 实战演练、调试技巧与性能考量 让我们用一个完整的例子来测试我们的实现，并讨论一些关键的实践要点。 ### 4.1 完整示例与代码测试 考虑本文开头提出的问题： 最大化 Z = 3x₁ + 5x₂ 约束： x₁ + x₂ ≥ 4 2x₁ + x₂ ≤ 10 x₁, x₂ ≥ 0 ```python if __name__ == "__main__": # 定义问题：最大化 3x1 + 5x2 c = [3, 5] A = [ [1, 1], # x1 + x2 >= 4 [2, 1] # 2x1 + x2 <= 10 ] b = [4, 10] constraint_types = ['>=', '<='] solver = TwoPhaseSimplex(c, A, b, constraint_types) result = solver.solve() if result: opt_val, opt_sol = result print(f"\n验证：") x1, x2 = opt_sol print(f" 目标函数值: 3*{x1:.2f} + 5*{x2:.2f} = {3*x1+5*x2:.2f}") print(f" 约束1 (>=4): {x1 + x2:.2f}") print(f" 约束2 (<=10): {2*x1 + x2:.2f}") ``` 运行这段代码，你应该能看到类似以下的输出，表明算法成功找到了最优解 `x₁=0, x₂=10`，最大值为50。 ``` === 第一阶段开始 === 第一阶段结束。W* = 0.000000 === 第二阶段开始 === 求解完成。最优值: 50.0000 最优解 (原变量): [ 0. 10.] 验证： 目标函数值: 3*0.00 + 5*10.00 = 50.00 约束1 (>=4): 10.00 约束2 (<=10): 10.00 ``` ### 4.2 常见陷阱与调试策略 在实现两阶段法时，以下几个坑点需要特别注意： 1. **数值稳定性**：浮点数计算会引入舍入误差。在判断检验数是否为0、比值是否为负时，应使用一个很小的容差（如 `1e-10`），而不是直接与0比较。 2. **退化与循环**：当最小比值检验出现平局时，可能导致迭代在几个基之间循环，永远无法达到最优。虽然理论上存在避免循环的规则（如Bland规则），但在实际中，由于数值误差，完全退化导致的无限循环较少见，但选择出基变量时添加一个确定性规则（如选择下标最小的变量）是良好的实践。 3. **人工变量仍在基中**：第一阶段结束时，可能出现人工变量仍在基中但其值为0的情况（退化）。在进入第二阶段前，必须通过额外的行变换将其替换出基，否则第二阶段单纯形表可能不可行。我们的代码通过 `_phase2_setup` 中的成本行变换来处理，但更严谨的做法是检查基变量列表，如果含有人工变量，尝试用非人工的非基变量进行主元旋转将其替换。 4. **无界解判断**：在选择进基变量后，如果对应列的所有系数都小于或等于0，则问题无界。代码中 `if not ratios:` 部分对此进行了处理。 为了便于调试，实现一个打印单纯形表的辅助方法非常有用。 ```python def _print_tableau(self, msg=""): """打印当前单纯形表（调试用）。""" print(f"\n--- {msg} ---") m, n = self.tableau.shape for i in range(m-1): row_str = " | ".join([f"{self.tableau[i, j]:7.2f}" for j in range(n-1)]) print(f"B[{self.basic_vars[i]:2d}]: [{row_str} | {self.tableau[i, -1]:7.2f}]") print("-" * (8*(n-1) + 10)) cost_row_str = " | ".join([f"{self.tableau[-1, j]:7.2f}" for j in range(n-1)]) print(f"Cost : [{cost_row_str} | {self.tableau[-1, -1]:7.2f}]") ``` ### 4.3 性能优化与扩展思考 我们实现的是一种教科书式的单纯形法。对于大规模问题，其性能可能成为瓶颈。以下是一些优化和扩展方向： * **稀疏矩阵存储**：实际问题的约束矩阵通常非常稀疏。使用 `scipy.sparse` 等库的稀疏矩阵数据结构可以极大减少内存占用和计算量。 * **修正单纯形法**：我们实现的是**表格单纯形法**，它每一步都操作整个矩阵，计算量大。**修正单纯形法**只存储和操作基矩阵的逆，在变量远多于约束时效率更高。其核心步骤是： 1. 计算基矩阵的逆 B⁻¹。 2. 计算单纯形乘子 π = c_Bᵀ B⁻¹。 3. 计算非基变量的检验数 σ_N = c_Nᵀ - π N。 4. 选择进基变量。 5. 计算进基列在原始空间的表示 a = B⁻¹ A_j。 6. 最小比值检验确定出基变量。 7. 更新基矩阵的逆（使用秩1更新，如乘积形式逆或LU更新）。 * **与大模型求解器的对比**：对于生产环境，更推荐使用成熟的优化库，如 `PuLP`、`CVXOPT` 或商业求解器 `Gurobi`、`CPLEX`。它们集成了更高效的算法（如内点法）、预处理技术、并行计算和鲁棒的数值处理。自己实现单纯形法的价值在于教学、定制化需求或理解算法本质。 * **处理各种边界情况**：我们的示例代码是一个简化版本。一个健壮的实现还需要处理： * 变量无界（free variables）。 * 所有约束右端项非负的预处理。 * 更复杂的退化处理策略。 两阶段单纯形法将寻找可行起点的艺术与单纯形迭代的科学相结合，是线性规划求解器中不可或缺的组成部分。通过亲手实现它，你不仅能深刻理解单纯形法如何从“无”到“有”地开始工作，还能洞察到数值计算中那些微妙而重要的细节。当你下次调用 `scipy.optimize.linprog` 并瞬间得到答案时，或许会想起这个在顶点间耐心行走、先搭脚手架再建高楼的经典算法。