第 6 章：单纯形法原理 —— 代数上的"爬山"算法¶
场景： 图解法告诉我们最优解在可行域的顶点，但如果有 100 个变量，顶点有上百万个，如何系统地找到正确的那个？ 单纯形法 （Simplex Method）由 George Dantzig 于 1947 年发明，它像智能导航一样，从一个顶点走到相邻的更好顶点，逐步逼近最优解。
6.1 单纯形法的核心思想¶
核心比喻：爬山找峰顶
想象你在浓雾中爬山（可行域内），你能做的是：
- 站在当前点（当前顶点）
- 环顾四周，找出相邻的更高的点（相邻顶点）
- 走过去（基变换）
- 重复直到四周没有更高的点
这就是单纯形法—— 每一步都严格改善目标值 ，最终停在最高点。
6.2 单纯形法的关键概念¶
基、基变量与非基变量¶
在标准形式中，从 \(n\) 个变量和 \(m\) 个等式约束中选取 \(m\) 个作为 基变量 （basic variables），其余 \(n-m\) 个为 非基变量 （non-basic variables）。
"""
奶茶店问题的标准形式：
min z' = -5x1 - 3x2 + 0*s1 + 0*s2
s.t. 2x1 + x2 + s1 = 10 (s1 是松弛变量)
x1 + 2x2 - s2 = 8 (s2 是剩余变量)
x1, x2, s1, s2 >= 0
初始可行基：B = [s1, s2]（单位矩阵）
初始基变量：s1 = 10, s2 = 8
初始非基变量：x1 = 0, x2 = 0
初始目标值：z' = 0, 即 z = 0
"""
print("单纯形法演示：奶茶店问题")
print("=" * 60)
print("\n初始基: B = [s1, s2]（单位矩阵，容易求逆）")
print("初始解: x1=0, x2=0, s1=10, s2=8")
print("初始目标值: z = 5×0 + 3×0 = 0")
print("\n这个解是什么意思？")
print(" - 不生产任何奶茶，销量为 0")
print(" - 所有原料都没有使用")
print(" - 利润为 0——显然不是最优，还有提升空间")
渲染效果：
单纯形法演示：奶茶店问题
============================================================
初始基: B = [s1, s2]（单位矩阵，容易求逆）
初始解: x1=0, x2=0, s1=10, s2=8
初始目标值: z = 5×0 + 3×0 = 0
6.3 最优性检验：检验数¶
非基变量的 检验数 （reduced cost）\(\sigma_j\) 告诉我们：将该变量从 0 增加 1 单位时，目标函数值会如何变化。
\(\sigma_j = c_j - \mathbf{y}^T \mathbf{a}_j\)
其中 \(\mathbf{y}^T = \mathbf{c}_B^T B^{-1}\) 是 影子价格向量 。
def compute_reduced_cost():
"""计算各变量的检验数（reduced cost）"""
# 当前基: B = [s1, s2] = 单位矩阵
B_inv = np.eye(2)
c_B = np.array([0, 0]) # 松弛/剩余变量的目标系数都是 0
y = np.dot(c_B, B_inv) # y = [0, 0]
print("当前基变量成本 c_B = [0, 0]")
print(f"影子价格 y = c_B * B^(-1) = {y}")
# 各变量的系数列
a_x1 = np.array([2, 1])
a_x2 = np.array([1, 2])
a_s1 = np.array([1, 0])
a_s2 = np.array([0, -1])
# 检验数 = c_j - y^T * a_j
sigma_x1 = 5 - np.dot(y, a_x1) # c1 = 5
sigma_x2 = 3 - np.dot(y, a_x2) # c2 = 3
print(f"\n检验数计算:")
print(f" σ_x1 = c1 - y^T*a1 = 5 - {y}·{a_x1} = {sigma_x1}")
print(f" σ_x2 = c2 - y^T*a2 = 3 - {y}·{a_x2} = {sigma_x2}")
print(f"\n解释: x1 从 0→1，z 增加 {sigma_x1} 元！")
print(f" x2 从 0→1，z 增加 {sigma_x2} 元！")
print(f"\n→ 两个检验数都是负的（对于 min 问题），说明当前解不是最优")
print(f"→ 应该让 x1 和 x2 进入基变量（生产奶茶）")
compute_reduced_cost()
渲染效果：
当前基变量成本 c_B = [0, 0]
影子价格 y = c_B * B^(-1) = [0. 0.]
检验数计算:
σ_x1 = c1 - y^T*a1 = 5 - [0. 0.]·[2 1] = 5.0
σ_x2 = c2 - y^T*a2 = 3 - [0. 0.]·[1 2] = 3.0
解释: x1 从 0→1，z 增加 5 元！
x2 从 0→1，z 增加 3 元！
→ 两个检验数都是负的（对于 min 问题），说明当前解不是最优
→ 应该让 x1 和 x2 进入基变量（生产奶茶）
检验数的直观理解
- 检验数 > 0（对 max 问题） ：让这个变量增加能提高利润，值得"投资"
- 检验数 < 0（对 max 问题） ：让这个变量增加会降低利润，应该"止损"
6.4 基变换：换入变量与换出变量¶
换入变量的选择¶
选取 负检验数最大 （对 max 问题）的非基变量作为换入变量——这个变量"最有价值"。
对于奶茶店问题（max 形式）： - \(\sigma_{x1} = 5\)（最大） - \(\sigma_{x2} = 3\)
→ 选 \(x_1\) 作为 换入变量 （每增加 1 单位利润增加最多）
换出变量的选择：最小比值原则¶
换入变量增加时，某些基变量会减少（因为总资源有限）。为了保持可行性，我们选择 减少到零最早 的那个基变量作为换出变量。
def determine_pivot():
"""
确定换入换出变量
换入: x1 (σ_x1 = 5 最大)
当前基: s1 = 10 - 2x1, s2 = 8 - x1
x1 增加时，s1 和 s2 都会减少
"""
print("换入变量: x1（检验数 5 最大）")
print("\n换出变量的最小比值分析:")
# 基变量表达式（从约束方程中解出）
print(" s1 = 10 - 2*x1 (当前值 10)")
print(" s2 = 8 - x1 (当前值 8)")
print("\n让 x1 增加，基变量如何变化:")
print(" 当 x1 = 5 时: s1 = 10 - 2*5 = 0 (s1 第一个到达 0)")
print(" 当 x1 = 8 时: s2 = 8 - 8 = 0 (s2 到达 0)")
print("\n → x1 最多增加到 5（受 s1 限制）")
print(" → s1 被换出，s2 继续作为基变量")
print("\n旋转元 (Pivot Element): a_11 = 2（s1 行的 x1 系数）")
determine_pivot()
渲染效果：
换入变量: x1（检验数 5 最大）
换出变量的最小比值分析:
s1 = 10 - 2*x1 (当前值 10)
s2 = 8 - x1 (当前值 8)
让 x1 增加，基变量如何变化:
当 x1 = 5 时: s1 = 10 - 2*5 = 0 (s1 第一个到达 0)
当 x1 = 8 时: s2 = 8 - 8 = 0 (s2 到达 0)
→ x1 最多增加到 5（受 s1 限制）
→ s1 被换出，s2 继续作为基变量
旋转元 (Pivot Element): a_11 = 2（s1 行的 x1 系数）
6.5 旋转运算（Pivot Operation）¶
用初等行变换把旋转元（pivot element）变成 1，同列其他元素变成 0：
def pivot_operation():
"""
旋转运算：把 x1 换入，s1 换出
旋转元: a_11 = 2
旋转步骤：
1. 旋转行除以 pivot → pivot 变成 1
2. 其他行减去 pivot 列的倍数 → 其他行的 pivot 列变成 0
"""
print("旋转运算（从初始表到第一次迭代）")
print("=" * 55)
# 初始单纯形表（增广矩阵）
print("\n初始表（旋转前）:")
print(" x1 x2 s1 s2 | b")
print("s1: 2 1 1 0 | 10")
print("s2: 1 2 0 1 | 8")
print("z: -5 -3 0 0 | 0")
print("\n旋转运算（以 a_11=2 为旋转元）:")
print("\n步骤1: s1 行 ÷ 2")
print(" x1 x2 s1 s2 | b")
print("s1: 1 0.5 0.5 0 | 5 ← 新 s1（实际是 x1）")
print("s2: 1 2 0 1 | 8")
print("z: -5 -3 0 0 | 0")
print("\n步骤2: s2 行 - s1 行×1")
print(" x1 x2 s1 s2 | b")
print("s1: 1 0.5 0.5 0 | 5")
print("s2: 0 1.5 -0.5 1 | 3 ← 新 s2")
print("z: -5 -3 0 0 | 0")
print("\n步骤3: z 行 + s1 行×5")
print(" x1 x2 s1 s2 | b")
print("s1: 1 0.5 0.5 0 | 5 ← x1 = 5")
print("s2: 0 1.5 -0.5 1 | 3 ← s2 = 3")
print("z: 0 -0.5 2.5 0 | 25 ← z = 25")
print("\n新基: B = [x1, s2]")
print("新解: x1=5, x2=0, s1=0, s2=3")
print("目标值: z = 25")
pivot_operation()
渲染效果：
旋转运算（从初始表到第一次迭代）
=======================================================
步骤1: s1 行 ÷ 2
...
新基: B = [x1, s2]
新解: x1=5, x2=0, s1=0, s2=3
目标值: z = 25
6.6 第二次迭代：继续改进¶
从新基 \(B = [x_1, s_2]\) 出发，继续检验最优性：
def second_iteration():
"""
第二次迭代
当前基: x1, s2
当前解: x1=5, x2=0, s1=0, s2=3
目标值: z=25
"""
print("第二次迭代：检验最优性")
print("=" * 50)
# B = [x1, s2] 的逆
# 从上一步结果反推
# x1 + 0.5*s1 = 5 - 0.5*s2?
# 实际上 B = [[1, 0], [0, 1.5]]?
print("当前基: B = [x1, s2]")
print("基变量: x1=5, s2=3")
print("非基变量: x2=0, s1=0")
# 检验数
print("\n计算检验数（对 min 问题）:")
# s1 的检验数 > 0？→ 继续迭代
# x2 的检验数: -0.5 < 0 → 不需要进基
print(" σ_s1 = 2.5 > 0 → s1 可以进入基（减少约束浪费）")
print(" σ_x2 = -0.5 < 0 → x2 不值得进入")
print("\n→ 非基变量 s1 的检验数为正（对 min 问题）")
print("→ s1 进基，增加 s1 会让 z' 变大（z 变小）")
print("→ 需要继续迭代！")
print("\n最小比值分析:")
print(" s1 行: 无约束（s1 前系数 0.5）")
print(" s2 行: (3 - 0) / 0.5 = 6")
print(" → s2 出基，s1 进基")
print("\n第二次旋转后:")
print(" x1 x2 s1 s2 | b")
print("x1: 1 2 0 -1 | 4 ← 新 x1")
print("s1: 0 3 -1 2 | 6 ← 新 s1")
print("z: 0 4 0 5 | 26 ← z = 26")
print("\n新解: x1=4, x2=2, s1=0, s2=0")
print("目标值: z = 26")
second_iteration()
6.7 最优性检验与收敛¶
def optimality_check():
"""最优性检验"""
print("最优性检验")
print("=" * 50)
print("\n当前解: x1=4, x2=2, z=26")
print("\n所有非基变量的检验数（min 问题要求 ≤ 0）:")
print(" σ_x2 = -0.5 ≤ 0 ✓")
print(" σ_s2 = 0.5 ≤ 0 ✓")
print("\n→ 所有检验数 ≤ 0，当前解是最优解！")
print("\n最终结果:")
print(" 每天生产 4 杯 A 奶茶，2 杯 B 奶茶")
print(" 牛奶消耗: 2×4 + 1×2 = 10 单位（刚好用完）")
print(" 茶叶消耗: 1×4 + 2×2 = 8 单位（刚好用完）")
print(" 最大利润: 5×4 + 3×2 = 26 元")
optimality_check()
渲染效果：
最优性检验
==================================================
当前解: x1=4, x2=2, z=26
所有非基变量的检验数（min 问题要求 ≤ 0）:
σ_x2 = -0.5 ≤ 0 ✓
σ_s2 = 0.5 ≤ 0 ✓
→ 所有检验数 ≤ 0，当前解是最优解！
6.8 单纯形法的收敛性¶
为什么单纯形法一定能在有限步内收敛？
- 可行域顶点个数有限（最多 \(C_n^m\) 个）
- 每次迭代目标函数值严格改善
- 不可能循环（除非退化情况）
- 所以一定能在有限步内找到最优解（或判断无界）
要点总结¶
- 单纯形法 = 从一个可行顶点沿边走到更好的相邻顶点
- 初始可行基通常选单位矩阵（松弛/剩余变量）
- 检验数 \(\sigma_j = c_j - y^T a_j\) 决定是否继续迭代
- 换入变量：负检验数最大（max 问题）/ 正检验数最小（min 问题）
- 换出变量：最小比值原则 \(\theta = \min(b_i / a_{ik})\)
- 旋转运算：用初等行变换更新基
- 所有检验数 \(\le 0\)（min）或 \(\ge 0\)（max）→ 最优解
课后练习¶
- 手算练习 ：用单纯形法求解： \(\max \quad z = 2x_1 + x_2\)
约束：\(x_1 + 2x_2 \le 6\), \(x_1 + x_2 \le 4\), \(x_1, x_2 \ge 0\)
-
Python 验证 ：用 scipy.optimize.linprog 验证上面的结果。
-
思考题 ：如果目标函数是 \(z = x_1 + 1000x_2\)，单纯形法会先让哪个变量进基？这反映了算法的什么特点？
下一章预告： 第 6 章建立了单纯形法的直觉。第 7 章将介绍 表格单纯形法 （更系统的手工计算工具）、 人工变量法 （处理非单位矩阵初始基）和 两阶段法 。