单纯形法是一种求解线性规划问题的经典算法，其核心思想是通过在可行域顶点间的迭代移动，逐步逼近并最终找到使目标函数最优的顶点 [ref_1]。其算法实现通常包含几个关键步骤：将问题转换为标准形式、构建单纯形表、迭代寻找最优解（包括确定入基和出基变量、进行枢轴运算等）[ref_3][ref_4]。 以下是一个遵循上述原理，能够处理标准最大化问题（目标函数为 `max z = c^T x`，约束为 `Ax <= b, x >= 0`，并假设 `b >= 0`）的Python实现。该实现包含了构建初始单纯形表、迭代优化以及结果提取的全过程 [ref_2][ref_5]。 ```python import numpy as np class SimplexSolver: """ 单纯形法求解器，用于求解标准形式的线性规划问题： 最大化: c^T * x 约束条件: A * x <= b, x >= 0 假设 b >= 0，以便构造初始可行基。 """ def __init__(self, c, A, b): """ 初始化单纯形求解器。 参数: c (list/np.array): 目标函数系数向量 (长度为 n)。 A (list/np.array): 约束系数矩阵 (m x n)。 b (list/np.array): 约束右侧常数向量 (长度为 m)，要求非负。 """ self.c = np.array(c, dtype=float).flatten() # 目标函数系数 self.A = np.array(A, dtype=float) # 约束矩阵 self.b = np.array(b, dtype=float).flatten() # 右侧常数 self.m, self.n = self.A.shape # m个约束，n个原始变量 # 构建初始单纯形表 (m+1 行， n+m+1 列) # 最后m列为松弛变量，最后一列为右侧常数(RHS) self.tableau = np.zeros((self.m + 1, self.n + self.m + 1)) # 填充约束部分: [A | I | b] self.tableau[:self.m, :self.n] = self.A self.tableau[:self.m, self.n:self.n+self.m] = np.eye(self.m) self.tableau[:self.m, -1] = self.b # 填充目标函数行: [-c | 0 | 0] self.tableau[-1, :self.n] = -self.c # 基变量索引 (初始为松弛变量的索引) self.basis = list(range(self.n, self.n + self.m)) def solve(self, max_iterations=1000): """ 执行单纯形算法主循环。 参数: max_iterations (int): 最大迭代次数，防止无限循环。 返回: dict: 包含最优值、最优解、状态等信息。 """ iteration = 0 while iteration < max_iterations: iteration += 1 # 1. 检查最优性：目标函数行（最后一行）的前n+m列（检验数）是否都非负 # 如果所有检验数都 >= 0，则当前解为最优解 [ref_1] if np.all(self.tableau[-1, :-1] >= -1e-10): # 考虑数值误差 break # 2. 选择入基变量：选择检验数最小的列（最负的，Dantzig规则）[ref_3] # 入基变量将进入基，替换一个现有的基变量 entering = np.argmin(self.tableau[-1, :-1]) # 3. 检查无界性：如果入基列的所有系数都 <= 0，则问题无界 [ref_4] if np.all(self.tableau[:-1, entering] <= 1e-10): return { 'status': 'unbounded', 'message': '问题无界，目标函数值可无限增大。' } # 4. 选择出基变量：使用最小比值测试 (Minimum Ratio Test) [ref_1][ref_5] # 计算比值：RHS / 入基列的正系数，选择比值最小的行 ratios = [] for i in range(self.m): if self.tableau[i, entering] > 1e-10: # 只考虑正系数 ratio = self.tableau[i, -1] / self.tableau[i, entering] ratios.append((ratio, i)) else: ratios.append((np.inf, i)) # 非正系数，比值为无穷大 # 找到比值最小的行，该行对应的基变量将出基 min_ratio, leaving_row = min(ratios, key=lambda x: x[0]) # 5. 更新基变量索引：用入基变量替换出基变量 leaving_var = self.basis[leaving_row] self.basis[leaving_row] = entering # 6. 执行枢轴运算 (Pivot Operation)，使入基列变为单位向量 [ref_2] pivot = self.tableau[leaving_row, entering] # 枢轴行除以枢轴元素 self.tableau[leaving_row, :] /= pivot # 其他行减去（该行入基列系数 * 新的枢轴行） for i in range(self.m + 1): if i != leaving_row: factor = self.tableau[i, entering] self.tableau[i, :] -= factor * self.tableau[leaving_row, :] if iteration >= max_iterations: return { 'status': 'max_iterations_exceeded', 'message': f'达到最大迭代次数 {max_iterations}，未找到最优解。' } # 提取结果 return self._get_solution() def _get_solution(self): """ 从最终的单纯形表中提取最优解和最优值。 """ solution = np.zeros(self.n + self.m) # 包含原始变量和松弛变量 # 基变量的值等于对应行的RHS值 for i, var_index in enumerate(self.basis): if var_index < len(solution): solution[var_index] = self.tableau[i, -1] # 原始变量的最优解 x_optimal = solution[:self.n] # 最优目标函数值 (单纯形表最后一行最后一列，符号取反) z_optimal = self.tableau[-1, -1] return { 'status': 'optimal', 'optimal_value': z_optimal, 'solution': x_optimal, 'basis': self.basis, 'tableau': self.tableau } def print_tableau(self, title="当前单纯形表"): """ 打印当前单纯形表，用于调试。 """ print(f"\n{title}:") print(self.tableau) print(f"基变量索引: {self.basis}") # ==================== 使用示例 ==================== if __name__ == "__main__": # 示例问题：最大化 z = 3x1 + 2x2 # 约束: # 2x1 + x2 <= 100 # x1 + x2 <= 80 # x1 <= 40 # x1, x2 >= 0 # 期望最优解: x1=20, x2=60, z=3*20+2*60=180 [ref_3] # 目标函数系数 (max 3x1 + 2x2 -> 输入系数 [3, 2]) c = [3, 2] # 约束矩阵 A (3个约束) A = [ [2, 1], [1, 1], [1, 0] ] # 约束右侧常数 b b = [100, 80, 40] # 创建求解器实例并求解 solver = SimplexSolver(c, A, b) print("初始单纯形表:") solver.print_tableau("初始表") result = solver.solve() print("\n" + "="*50) print("求解结果:") print(f"状态: {result['status']}") if result['status'] == 'optimal': print(f"最优目标函数值: {result['optimal_value']:.2f}") for i, val in enumerate(result['solution']): print(f" x{i+1} = {val:.2f}") # 验证结果 z_calc = np.dot(c, result['solution']) print(f"验证计算的目标值: {z_calc:.2f}") print(f"基变量索引 (对应单纯形表中的列): {result['basis']}") ``` 上述代码的核心逻辑与步骤如下表所示： | 步骤 | 关键操作 | 目的与说明 | | :--- | :--- | :--- | | **1. 初始化** | 构建 `(m+1) x (n+m+1)` 的单纯形表，将目标函数行设为 `-c`，约束部分加入单位矩阵的松弛变量 [ref_2]。 | 将不等式约束 `Ax <= b` 转化为等式 `Ax + s = b`，并构造一个初始基可行解（松弛变量为基变量）。 | | **2. 最优性检验** | 检查目标函数行（最后一行）前 `n+m` 个检验数是否全部 `>= 0` [ref_1]。 | 若全部非负，则当前基可行解即为最优解，算法终止。 | | **3. 选择入基变量** | 选取检验数中**最小（最负）** 的列索引 `entering`（Dantzig规则）[ref_3]。 | 该变量进入基将最大程度地提升目标函数值。 | | **4. 无界性检验** | 检查入基列 `tableau[:-1, entering]` 的所有系数是否 `<= 0`。 | 若全是非正数，则目标函数可无限增大，问题无界。 | | **5. 选择出基变量** | 对每一行 `i`，计算比值 `rhs_i / a_i_entering`（仅当 `a_i_entering > 0`），选取**最小非负比值**对应的行 `leaving_row` [ref_5]。 | 保证迭代后所有变量仍非负，并确定哪个基变量应离开基。 | | **6. 枢轴运算** | 1. 将 `leaving_row` 行除以枢轴元素 `pivot`。<br>2. 对其他所有行，减去该行入基列系数与新的枢轴行的乘积。 | 通过行变换，使入基列变为单位向量（该列仅在 `leaving_row` 行为1，其余为0），完成基的替换 [ref_2]。 | | **7. 迭代与结果提取** | 重复步骤2-6直至满足最优条件。最后，基变量的值等于其对应行的RHS值，非基变量值为0；最优值位于表格右下角 [ref_4]。 | 获得最优解 `x_optimal` 和最优目标值 `z_optimal`。 | **重要注意事项与扩展**： 1. **标准形式**：上述代码假设问题是**最大化**且约束为“小于等于”。对于最小化问题或包含“大于等于”约束，需引入人工变量或使用两阶段法、大M法进行预处理 [ref_3]。 2. **初始可行基**：代码依赖 `b >= 0` 的假设，通过松弛变量直接获得初始可行基。若 `b` 包含负数，则需要使用两阶段单纯形法来寻找初始基可行解 [ref_6]。 3. **数值稳定性**：代码中使用了 `1e-10` 作为零的容差，以处理浮点数计算误差。 4. **对偶单纯形法**：当初始解不可行但对偶解可行时（常见于添加新约束后），可以使用对偶单纯形法，它从对偶可行的解开始，逐步恢复原始可行性直至达到最优 [ref_6]。其Python实现逻辑与本例类似，但选择出基和入基变量的规则恰好相反。 5. **实际应用**：对于大规模问题，推荐使用更专业的优化库，如 `scipy.optimize.linprog`（内部可能使用修订单纯形法）或 `PuLP`、`CVXOPT` 等。这里的实现主要用于教育目的，以清晰展示单纯形法的核心迭代机制 [ref_1][ref_5]。