# 从零到一：用Python构建你的第一个单纯形法求解器 如果你曾经尝试过手动求解一个超过三个变量的线性规划问题，大概率会陷入表格和数字的海洋，感觉既繁琐又容易出错。这正是运筹学中**单纯形法**的魅力所在——它提供了一套系统化的迭代规则，能将复杂的优化问题转化为一系列可计算的步骤。而今天，我们不再满足于纸笔演算，而是要亲手用Python将这套精妙的数学逻辑“翻译”成代码，打造一个属于你自己的、能自动寻找最优解的求解引擎。 这篇文章面向的是那些已经熟悉Python基础语法，但对运筹学或线性规划算法实现感到好奇的开发者。我们不会止步于调用现成的`scipy.optimize.linprog`，而是要深入算法的“黑箱”，从零开始，一步步构建核心的迭代过程。你会看到数学中的“入基”、“出基”、“检验数”如何变成程序中的变量和条件判断，理解算法每一步的决策逻辑，并最终获得一个可以处理标准形式线性规划问题的、功能完整的Python类。这不仅是学习一个算法，更是一次绝佳的“将理论转化为实践”的思维训练。 ## 1. 理解战场：线性规划与单纯形法核心概念 在开始编码之前，我们必须清晰地定义战场。一个线性规划问题，通常是在一组线性等式或不等式的约束下，寻找一个线性目标函数的最大值或最小值。为了应用单纯形法，我们首先需要将问题转化为**标准形式**。 **标准形式**有三个关键特征： 1. 目标函数为**最大化**（Maximize）。 2. 所有约束条件均为**等式**。 3. 所有决策变量均为**非负**。 例如，一个生产计划问题：“生产产品A每件利润3元，耗时2小时；产品B每件利润4元，耗时1小时。总工时40小时，且产品B的生产受限于某种原料，最多关联生产30件（简化模型）。”其数学模型和标准形式转化如下： * **原始模型**： * 最大化利润：`Z = 3*x1 + 4*x2` * 约束： * `2*x1 + x2 <= 40` （工时约束） * `x1 + 3*x2 <= 30` （原料关联约束） * `x1, x2 >= 0` * **转化为标准形式**： 我们需要引入**松弛变量** `x3` 和 `x4`，将不等式变为等式。 * 最大化：`Z = 3*x1 + 4*x2 + 0*x3 + 0*x4` * 约束： * `2*x1 + x2 + x3 = 40` * `x1 + 3*x2 + x4 = 30` * `x1, x2, x3, x4 >= 0` 这里，`x3`和`x4`就是松弛变量，它们代表了未被利用的工时和原料“余量”，在目标函数中系数为0，表示它们不直接产生利润。 > 提示：将最小化问题转为最大化，只需将目标函数系数乘以-1。对于“>=”约束，则需要引入“剩余变量”和“人工变量”，这属于两阶段法或大M法的范畴，我们将在后续进阶部分简要讨论。 单纯形法的核心思想是在**基可行解**的顶点之间“跳跃”，每次跳跃都让目标函数值变得更好，直到找到最优解。这涉及到几个核心操作，我们将用一张表来厘清它们在数学和程序中的对应关系： | 数学概念 | 程序中的对应物 | 核心作用 | | :--- | :--- | :--- | | **基变量 (Basic Variables)** | 一组索引，指向当前解中不为零的变量。 | 构成当前可行解的“骨架”，其数量等于约束等式的个数。 | | **非基变量 (Non-basic Variables)** | 除基变量索引外的其他变量索引，其值固定为0。 | 等待进入基变量集合、可能改善目标函数的候选者。 | | **单纯形表 (Simplex Tableau)** | 一个增广矩阵，通常用二维数组（如NumPy矩阵）表示。 | 存储所有系数（A）、约束右端项（b）、目标函数系数（c）以及当前解的信息，是迭代运算的载体。 | | **检验数 (Reduced Cost, σ)** | 对于最大化问题，计算 `c_j - (当前基变量系数向量与对应A矩阵列的內积)`。 | **入基变量选择的依据**。对于最大化问题，任何正检验数的非基变量入基，都有可能提升目标函数值。 | | **θ比率 (Theta Ratio)** | 对每个约束行，用当前解的右端项(`b_i`)除以入基变量在该行的正系数(`A[i, enter_idx]`)。 | **出基变量选择的依据**。选择最小非负θ值对应的基变量出基，能保证迭代后所有变量仍满足非负约束。 | | **枢轴元 (Pivot Element)** | 入基变量列与出基变量行交叉点的系数。 | 执行**行变换**的基准点，通过它将该列化为单位向量，从而完成基的替换。 | 理解了这张“地图”，我们就可以开始搭建程序的“脚手架”了。 ## 2. 搭建框架：设计SimplexSolver类 好的程序结构始于清晰的设计。我们将创建一个 `SimplexSolver` 类，它封装单纯形法求解的全部状态和操作。类的初始化需要接收标准形式的输入。 ```python import numpy as np class SimplexSolver: """ 单纯形法求解器（仅处理标准最大化形式，约束为等式，变量非负）。 假设初始可行基由松弛变量构成。 """ def __init__(self, c, A, b): """ 初始化求解器。 参数: c (list/np.array): 目标函数系数向量，长度 n (决策变量+松弛变量)。 A (list/np.array): 约束系数矩阵，形状 (m, n)。 b (list/np.array): 约束右端项向量，长度 m，应全部非负。 """ # 转换为NumPy数组以便计算 self.c = np.array(c, dtype=float).flatten() # 目标函数系数 (1 x n) self.A = np.array(A, dtype=float) # 约束矩阵 (m x n) self.b = np.array(b, dtype=float).flatten() # 右端项 (m,) self.m, self.n = self.A.shape # m个约束，n个变量 # 验证输入维度 assert len(self.c) == self.n, "c的长度必须等于A的列数(n)" assert len(self.b) == self.m, "b的长度必须等于A的行数(m)" assert (self.b >= 0).all(), "初始右端项b应全部非负（或通过预处理保证）" # 构建初始单纯形表（增广矩阵） # 表格结构: [ A | b ] # [ c | 0 ] (实际存储时，c行是检验数计算的基础) self.tableau = np.zeros((self.m + 1, self.n + 1)) self.tableau[:-1, :-1] = self.A self.tableau[:-1, -1] = self.b # 最后一行前n列存储的是原始c的负值，便于检验数计算，具体见`_compute_tableau`方法 # 初始基变量假设：后m个变量是松弛变量，构成初始基。 # 这是基于问题已通过添加松弛变量转化为标准形式的前提。 self.basic_vars = list(range(self.n - self.m, self.n)) # 基变量索引列表 self.non_basic_vars = [i for i in range(self.n) if i not in self.basic_vars] # 当前目标函数值 self.obj_value = 0.0 # 记录迭代次数 self.iterations = 0 # 状态: 'optimal', 'unbounded', 'infeasible', 'processing' self.status = 'processing' # 初始化表格（计算初始检验数行） self._compute_tableau() ``` 这个 `__init__` 方法完成了数据的接收、验证和基础结构的搭建。其中，`_compute_tableau` 是一个关键的内部方法，它负责根据当前的基变量集合，计算单纯形表的最后一行（检验数行）和当前的目标函数值。让我们来实现它： ```python def _compute_tableau(self): """计算当前基下的单纯形表（主要是检验数行和目标值）。""" # 提取基变量对应的系数矩阵部分 B 和其在目标函数中的系数 c_B B = self.A[:, self.basic_vars] c_B = self.c[self.basic_vars] # 计算基的逆（实际计算中，我们通过行变换更新整个表，而非显式求逆） # 这里更准确的做法是：当前表格的基变量列已经通过之前的行变换变成了单位矩阵。 # 我们直接利用更新后的表格计算检验数。 # 计算检验数 sigma_j = c_j - c_B * B^{-1} * A_j # 对于当前表格，基变量列是单位阵，c_B * B^{-1} 就是 c_B 本身（因为B^{-1}作用后是单位向量）。 # 更通用的计算方式：y = c_B (在当前表结构下，由于基列是单位阵，y就是c_B) y = c_B # 计算所有变量的检验数 for j in range(self.n): if j in self.basic_vars: # 基变量的检验数应为0（理论上） self.tableau[-1, j] = 0 else: # 非基变量检验数: c_j - y * A[:, j] self.tableau[-1, j] = self.c[j] - np.dot(y, self.A[:, j]) # 计算当前目标函数值: c_B * b self.obj_value = np.dot(c_B, self.b) self.tableau[-1, -1] = -self.obj_value # 通常表中存储的是 -Z，这里我们存储 -obj_value ``` > 注意：在教科书的单纯形表表示中，最后一行通常存储 `c_j - z_j`（即检验数）和当前的 `-Z` 值。我们的实现遵循了这一惯例，将目标函数值的相反数放在表格右下角。 ## 3. 核心迭代：入基、出基与枢轴运算 有了初始化的表格，单纯形法就进入了一个循环：检查是否最优，如果不是，则选择入基和出基变量，执行枢轴运算（行变换）更新表格。这是算法跳动的心脏。 ### 3.1 最优性判定与入基变量选择 对于最大化问题，最优解的一个判定条件是：**所有非基变量的检验数都小于或等于0**。如果有正检验数，说明让对应的非基变量“入场”（从0增加）能提升目标函数。 ```python def _get_entering_variable(self): """根据最大检验数规则选择入基变量。返回索引，若无正检验数则返回-1。""" # 最后一行前n列是检验数 sigma_row = self.tableau[-1, :-1] # 只考虑非基变量（尽管基变量检验数理论为0，但计算误差可能导致非零） # 简化处理：直接找最大正检验数 max_sigma = 0 enter_idx = -1 for j in self.non_basic_vars: if self.tableau[-1, j] > max_sigma: max_sigma = self.tableau[-1, j] enter_idx = j return enter_idx ``` ### 3.2 无界性判断与出基变量选择 选择了入基变量 `enter_idx` 后，我们需要决定哪个现有的基变量要“离场”。这通过计算θ比率来完成。如果入基变量在所有约束中的系数都非正，则问题可能是无界的（目标函数可以无限增大）。 ```python def _get_leaving_variable(self, enter_idx): """ 根据最小比值规则选择出基变量。 参数: enter_idx: 入基变量的索引。 返回: 出基变量在`basic_vars`列表中的位置索引，若无可行比值则返回-1（无界）。 """ # 获取入基变量在约束中的系数列（前m行） pivot_column = self.tableau[:-1, enter_idx] # 获取当前右端项（前m行，最后一列） rhs = self.tableau[:-1, -1] min_ratio = float('inf') leave_pos = -1 # 记录在basic_vars列表中的位置 for i in range(self.m): if pivot_column[i] > 1e-10: # 避免除零和数值误差，只考虑正系数 ratio = rhs[i] / pivot_column[i] if ratio >= 0 and ratio < min_ratio: # 比值应非负 min_ratio = ratio leave_pos = i # 第i个约束行对应的基变量 return leave_pos # 如果为-1，表示无界 ``` ### 3.3 枢轴运算：高斯-约当消元 确定了入基变量列 (`enter_idx`) 和出基变量行 (`leave_pos`) 后，交叉点的元素就是**枢轴元**。我们需要通过行变换，将枢轴元变为1，并将其所在列的其他元素（包括检验数行）变为0。这本质上是执行一次高斯-约当消元，用入基变量替换出基变量，形成新的基。 ```python def _pivot(self, enter_idx, leave_pos): """执行枢轴运算，更新单纯形表。""" # leave_pos 是出基变量在 basic_vars 列表中的索引 leave_idx = self.basic_vars[leave_pos] # 出基变量的全局索引 # 1. 更新基变量和非基变量列表 self.basic_vars[leave_pos] = enter_idx self.non_basic_vars.remove(enter_idx) self.non_basic_vars.append(leave_idx) # 2. 对枢轴行进行归一化：使枢轴元变为1 pivot_element = self.tableau[leave_pos, enter_idx] self.tableau[leave_pos, :] /= pivot_element # 3. 对其他行（包括检验数行）进行消元：使入基变量列变为单位向量 for i in range(self.m + 1): if i != leave_pos: # 跳过枢轴行本身 factor = self.tableau[i, enter_idx] if abs(factor) > 1e-10: # 如果该行在此列系数不为零 self.tableau[i, :] -= factor * self.tableau[leave_pos, :] # 4. 由于行变换，右端项b和目标值已自动更新，但需要同步更新内部状态 self.b = self.tableau[:-1, -1].copy() self.obj_value = -self.tableau[-1, -1] # 表中存储的是 -Z self.iterations += 1 ``` ## 4. 组装与测试：让求解器运转起来 现在，我们将所有部件组装到主循环中，并添加一个友好的求解接口。 ```python def solve(self, max_iterations=1000): """执行单纯形法主循环。""" print(f"开始单纯形法求解，初始基变量: {self.basic_vars}") print(f"初始目标函数值: {self.obj_value:.2f}") while self.iterations < max_iterations: # 步骤1: 检查最优性，选择入基变量 enter_idx = self._get_entering_variable() if enter_idx == -1: self.status = 'optimal' print(f"在第 {self.iterations} 次迭代后达到最优解。") break # 步骤2: 选择出基变量 leave_pos = self._get_leaving_variable(enter_idx) if leave_pos == -1: self.status = 'unbounded' print("问题无界（目标函数可无限增大）。") break # 步骤3: 执行枢轴运算 print(f"迭代 {self.iterations+1}: 入基 x{enter_idx+1}, 出基 x{self.basic_vars[leave_pos]+1}") self._pivot(enter_idx, leave_pos) print(f" 当前目标值: {self.obj_value:.2f}, 基变量: {[idx+1 for idx in self.basic_vars]}") else: print(f"达到最大迭代次数 {max_iterations}，可能未收敛。") self.status = 'max_iterations_exceeded' return self.status, self.obj_value, self.get_solution() def get_solution(self): """获取当前解向量。""" solution = np.zeros(self.n) for i, var_idx in enumerate(self.basic_vars): solution[var_idx] = self.tableau[i, -1] # 基变量的值等于对应行的右端项 # 非基变量值默认为0 return solution ``` 是时候用我们开头的例子来测试一下了！ ```python # 定义问题：最大化 Z = 3*x1 + 4*x2 # 约束: # 2*x1 + x2 <= 40 -> 引入松弛变量 x3 # x1 + 3*x2 <= 30 -> 引入松弛变量 x4 # 标准形式系数: c = [3, 4, 0, 0] # 目标函数系数 (x1, x2, x3, x4) A = [ [2, 1, 1, 0], # 第一个约束 [1, 3, 0, 1] # 第二个约束 ] b = [40, 30] # 右端项 # 创建求解器并求解 solver = SimplexSolver(c, A, b) status, opt_value, solution = solver.solve() print("\n" + "="*50) print("求解结果:") print(f"状态: {status}") print(f"最优目标函数值: {opt_value}") print(f"最优解: x1 = {solution[0]:.2f}, x2 = {solution[1]:.2f}") print(f"松弛变量: x3 = {solution[2]:.2f}, x4 = {solution[3]:.2f}") print("="*50) ``` 运行这段代码，你会看到类似以下的输出，它清晰地展示了算法的迭代过程： ``` 开始单纯形法求解，初始基变量: [2, 3] 初始目标函数值: 0.00 迭代 1: 入基 x2, 出基 x4 当前目标值: 40.00, 基变量: [3, 2] 迭代 2: 入基 x1, 出基 x3 当前目标值: 70.00, 基变量: [1, 2] 在第 2 次迭代后达到最优解。 ================================================== 求解结果: 状态: optimal 最优目标函数值: 70.0 最优解: x1 = 18.00, x2 = 4.00 松弛变量: x3 = 0.00, x4 = 0.00 ================================================== ``` 看，我们的求解器正确地找到了最优解：生产18件产品A和4件产品B，可以获得最大利润70元，并且工时和原料约束都被完全利用（松弛变量为0）。这个迭代过程——`x2`入基替换`x4`，然后`x1`入基替换`x3`——与手动计算的单纯形表迭代完全一致，但现在是由你的代码自动完成的！ ## 5. 超越基础：处理更复杂的情况与优化 我们构建的求解器虽然能处理标准形式的松弛变量初始基，但现实问题往往更复杂。一个健壮的求解器还需要考虑以下几点： **1. 初始可行基的寻找（两阶段法）：** 很多问题在添加松弛变量后，松弛变量的系数可能不是单位阵（例如处理“>=”约束时引入了人工变量），或者右端项`b`有负数。这时初始解可能不可行。**两阶段法**是标准解决方案： * **第一阶段**：引入人工变量，构造一个辅助问题，其目标是最小化人工变量之和。用单纯形法求解这个辅助问题。如果最优值为0，则找到了原问题的一个可行基；否则原问题无可行解。 * **第二阶段**：用第一阶段找到的可行基作为起点，移除人工变量，恢复原目标函数，继续用单纯形法求解。 **2. 避免循环与退化：** 理论上，单纯形法可能在某些退化情况下（基变量值为0）陷入循环，无限迭代。**布兰德规则**是一种避免循环的确定性规则：在入基和出基变量选择出现多个候选时，总是选择下标最小的变量。 **3. 数值稳定性优化：** * **枢轴元选择**：在入基变量列有多个正系数时，选择θ比率最小的行出基，但有时选择绝对值较大的枢轴元能提高数值稳定性。 * **矩阵更新**：对于大规模问题，显式维护和更新基矩阵`B`的逆矩阵（乘积形式或LU分解形式），比操作整个单纯形表更高效、更稳定。这就是**修正单纯形法**的核心思想。 **4. 添加对偶单纯形法：** 当初始解是“对偶可行”（检验数行满足最优性条件）但原始不可行（右端项有负数）时，对偶单纯形法非常有效。它通过在行之间迭代，逐步消除不可行性，同时保持对偶可行性（检验数非正）。这在处理增加新约束或灵敏度分析时特别有用。 实现这些进阶功能会大大增加代码的复杂度，但也正是这些挑战让运筹优化软件（如CPLEX, Gurobi）变得强大。作为练习，你可以尝试为你的求解器添加两阶段法的第一阶段。这里是一个简化的思路框架： ```python def phase_one(self, A, b): """两阶段法的第一阶段，寻找初始可行基。""" m, n_orig = A.shape # 为每个约束添加人工变量 # 构造辅助问题：minimize sum of artificial vars # 等价于 maximize -sum of artificial vars # 具体实现需要扩展A矩阵，构造新的c向量... # 调用simplex求解... # 判断最优值是否为0，并提取可行基... pass ``` 亲手实现一个单纯形法求解器，最大的收获不是代码本身，而是对线性优化核心逻辑的深刻理解。你不再将`linprog()`视为魔法，而是明白了每一次“求解”背后进行的矩阵变换和逻辑决策。下次当你面对资源分配、排产计划或投资组合优化问题时，你不仅知道如何建模，更清楚求解器会如何一步步为你找到最优方案。这种从数学到代码的贯通感，是单纯调用库函数无法替代的。