Simplex Algorithm Done Right

registerGen · 2026-03-28 22:43:19 · 算法·理论

本文同步发表于博客园．本文假设读者具有基本的线性代数知识．

摘要 本文尝试用线性代数的语言叙述单纯形法，证明其正确性并导出线性规划基本定理．本文同时给出了单纯形法的 Rust 语言实现．

单纯形法 (simplex algorithm) 适用于求解如下的含 n 个变量和 m 个约束条件的线性规划标准型问题：

称满足式 (1) 中两个约束条件的 \bm x 为式 (1) 的一个 可行解 (feasible solution)．

从几何直观上讲，\bm A \bm x \leqslant \bm b 和 \bm x \geqslant \bm 0 围成了一个 \mathbb R^n 中的凸集，最优解 (optimal solution)（若存在）一定可以在凸集的顶点处取得（定理 10）．单纯形法以凸集的一个顶点作为初始解，每次迭代从当前解所在顶点的相邻顶点中选择一个更优的作为新的解，直到所有相邻顶点都不更优为止．单纯形法的大致流程如下：

1. 选定一组初始可行解；
2. 不断迭代，直到找到最优解或检测到解无界为止：
   2a. 判断当前解是否为最优解，如是，则退出迭代；
   2b. 选定解移动的方向；
   2c. 选定解移动的距离，若为无穷大，则解无界，退出迭代；
   2d. 计算移动后的解，并将其作为新的当前解，回到步骤 2a．

下面我们从代数的角度分析这一算法．首先，我们添加 m 个 松弛变量 (slack variable) \bm s \geqslant \bm 0，将 \bm A \bm x \leqslant \bm b 中的不等号变为等号：

并记 \bar{\bm A} = \begin{pmatrix} \bm A & \bm I_m \end{pmatrix}，\bar{\bm x} = \begin{pmatrix}\bm x^{\text{T}} & \bm s^{\text{T}}\end{pmatrix}^{\text{T}}，\bar{\bm c} = \begin{pmatrix} \bm c^{\text{T}} & \bm 0^{\text{T}}\end{pmatrix}^{\text{T}}，式 (1) 便化为

对于任意一个凸集顶点处对应的解（称为 顶点解 (corner-point solution)），其应能使得 \bm A \bm x \leqslant \bm b 和 \bm x \geqslant \bm 0 这 n + m 个不等式中的 n 个取到等号．若 \bm A \bm x \leqslant \bm b 中的某个不等式取到等号，则对应的松弛变量为 0；若 \bm x \geqslant \bm 0 中的某个不等式取到等号，则对应的变量为 0．总之，对任意顶点解 \bar{\bm x}，\bar{\bm x} 有 n 个分量为 0．我们将这 n 个分量称作 非基变量 (non-basic variable)，其余 m 个称作 基变量 (basic variable)．用 \mathcal B 和 \mathcal N 表示基变量和非基变量的下标集合，\bm x_{\mathcal B} 和 \bm x_{\mathcal N} 表示 \bar{\bm x} 中基变量和非基变量组成的列向量，\bm c_{\mathcal B} 和 \bm c_{\mathcal N} 表示 \bar{\bm c} 中基变量和非基变量对应的系数组成的列向量，\bm B 和 \bm N 表示 \bar{\bm A} 中基变量和非基变量所在列组成的矩阵．使用这些记号，我们将式 (3) 改写如下：

我们假设 \bm B 总是可逆的（从后面的讨论中可知在整个过程中这一假设始终成立），这样 \bm B 便是 \mathbb R^m 中的一组基．

对于每次迭代，我们首先选定一组基 \bm B，在式 (4) 中用约束条件从目标函数中消去 \bm x_{\mathcal B} 得

在式 (5) 中令 \bm x_{\mathcal N} = \bm 0 即得基为 \bm B 时的解

现在我们尝试寻找一个更优的顶点解．对于一个方向 \bm d = \begin{pmatrix} \bm d_{\mathcal B}^{\text{T}} & \bm d_{\mathcal N}^{\text{T}} \end{pmatrix}^{\text{T}} \in \mathbb R^{n+m}，令当前解 \bar{\bm x} 向 \bm d 这一方向移动一个充分小的标量 \tau \gt 0 得到新解 \bar{\bm x}' = \bar{\bm x} + \tau \bm d．首先 \bar{\bm x} 需满足式 (3) 中的约束条件 \bar{\bm A} \bar{\bm x}' = \bm b 和 \bar{\bm x}' \geqslant \bm 0．对前者，由 \bar{\bm A} \bar{\bm x} = \bm b 知

这是一个充要条件．而对后者，由于 \tau 是小量，\bar{\bm x}' \geqslant \bm 0 的充要条件是

其次 \bar{\bm x}' 要使得 Z' = \bar{\bm c}^{\text{T}} \bar{\bm x}' \gt Z，而 Z = \bar{\bm c}^{\text{T}} \bar{\bm x}，故充要条件为

当我们从一个顶点移动到另一个顶点时，\bm A \bm x \leqslant \bm b 和 \bm x \geqslant \bm 0 这 n + m 个不等式中，原来的 n 个取到等号的不等式中仍有 n - 1 个取到等号，其对应的非基变量的值仍为 0；而剩下的那个取不到等号，其对应的非基变量的值就从 0 变成了一个正数．由此，我们便得到了 n 个可能的方向——对于每个 k \in \mathcal N，当 \bar x_k 变成正数时，其对应的方向 \bm d^{(k)} 就有 d^{(k)}_k = 1 和 d^{(k)}_j = 0（对所有 j \in \mathcal N \setminus \{k\}）．将其代入式 (7)，同时将 \bar{\bm A} 改写为 \begin{pmatrix} \bm B & \bm N \end{pmatrix} 并记 \bm a_k 为 \bm N 中 \bar x_k 所在的列，我们有 \bar{\bm A} \bm d = \bm B \bm d^{(k)}_{\mathcal B} + \bm N \bm d^{(k)}_{\mathcal N} = \bm B \bm d^{(k)}_{\mathcal B} + \bm a_k = \bm 0，解得

若所有基变量取值均不为 0（称为 非退化 (non-degenerate) 情形），即

则 \bm d^{(k)} 一定满足式 (8)．式 (14) 会自动处理退化情形，详见后文．又记 c_k 为 \bar{\bm c} 的第 k 个分量，如有

则 \bm d^{(k)} 满足式 (9)，也即 \bm d^{(k)} 是一个可使目标函数 Z 增大的方向．

我们有如下的定理 1．注意式 (13) 的左边正是式 (5) 中 \bm x_{\mathcal N} 的系数，于是我们可在解得式 (5) 的同时判断当前解是否为最优解．

定理 1 (最优性判别) 对于从式 (3) 的 \bar{\bm A} 中选定的一组基 \bm B，若

-\bm c_{\mathcal B}^{\text{T}} \bm B^{-1} \bm N + \bm c_{\mathcal N}^{\text{T}} \leqslant \bm 0^{\text{T}}, \tag{13}

则式 (6) 给出了式 (3) 的一组最优解，因此我们无法找到使解更优的方向．

证明 设 \hat{\bar{\bm x}} 为式 (3) 的任意一组可行解，\hat{\bar{\bm x}} = \begin{pmatrix} \hat{\bm x}_{\mathcal B}^{\text{T}} & \hat{\bm x}_{\mathcal N}^{\text{T}} \end{pmatrix}^{\text{T}}，注意到 \bm B, \bm N, \bm c_{\mathcal B}^{\text{T}}, \bm c_{\mathcal N}^{\text{T}}, \bm b 是固定的，于是类似式 (5) 有 \hat Z = \bm c_{\mathcal B}^{\text{T}} \bm B^{-1} \bm b + (\bm c_{\mathcal N}^{\text{T}} - \bm c_{\mathcal B}^{\text{T}} \bm B^{-1} \bm N) \hat{\bm x}_{\mathcal N}．与式 (6) 比较得 \hat Z - Z = (\bm c_{\mathcal N}^{\text{T}} - \bm c_{\mathcal B}^{\text{T}} \bm B^{-1} \bm N) \hat{\bm x}_{\mathcal N}，由式 (13) 和 \hat{\bm x}_{\mathcal N} \geqslant \bm 0 知 \hat Z - Z \leqslant 0，即 Z 是一组最优解．\blacksquare

现在我们讨论当方向 \bm d^{(k)} 确定后，最多能沿它移动多远的距离．设移动后的解为 \bar{\bm x}^{(k)} = \bar{\bm x} + t \bm d^{(k)}，其中 t \geqslant 0 为移动的距离．容易证明 \bar{\bm x}^{(k)} 满足 \bar{\bm A} \bar{\bm x}^{(k)} = \bm b 且 \bar{\bm x}^{(k)} 是优于 \bar{\bm x} 的解（式 (7) 和式 (9) 不依赖于 \tau 是小量）．由 \bm d^{(k)}_{\mathcal N} 的构造过程可知 \bm d^{(k)}_{\mathcal N} \geqslant \bm 0，故 \bm x_{\mathcal N} + t \bm d_{\mathcal N}^{(k)} \geqslant \bm 0 始终成立，因此欲使 \bar{\bm x}^{(k)} \geqslant \bm 0，需成立 \bm x_{\mathcal B} + t \bm d_{\mathcal B}^{(k)} = \bm x_{\mathcal B} - t \bm B^{-1} \bm a_k \geqslant \bm 0，即

其中 \bm e_j 为第 j 个分量为 1，其余分量为 0 的列向量．特别地，若所有 \bm e_j^{\text{T}} \bm B^{-1} \bm a_k 均非正，则 t 可取任意正数，解无界．

定理 2 (无界性判别) 对于从式 (3) 的 \bar{\bm A} 中选定的一组基 \bm B 和一个满足式 (12) 的 k \in \mathcal N，若 \bm B^{-1} \bm a_k \leqslant \bm 0，则解无界．

证明 作射线 \bar{\bm x}(t) = \bar{\bm x} + t \bm d^{(k)}，其中 t \geqslant 0．由上文讨论知 \bar{\bm A} \bar{\bm x}(t) = \bm b 和 \bm x_{\mathcal N}(t) = \bm x_{\mathcal N} + t \bm d_{\mathcal N}^{(k)} \geqslant \bm 0，而由 \bm B^{-1} \bm a_k \leqslant \bm 0 知 \bm x_{\mathcal B}(t) = \bm x_{\mathcal B} - t \bm B^{-1} \bm a_k \geqslant \bm 0 对任意 t \geqslant 0 均成立，故 \bar{\bm x}(t)（t \geqslant 0）都是式 (3) 的可行解．设 \bar{\bm x}(t) 对应的目标函数值为 Z(t)，则有 Z(t) = \bar{\bm c}^{\text{T}} \bar{\bm x}(t) = \bar{\bm c}^{\text{T}} \bar{\bm x} + t \bar{\bm c}^{\text{T}} \bm d^{(k)} = Z + t \bar{\bm c}^{\text{T}} \bm d^{(k)}．由式 (12)，\bar{\bm c}^{\text{T}} \bm d^{(k)} \gt 0，故当 t \to +\infty 时 Z(t) \to +\infty，即解无界．\blacksquare

当非退化（式 (11)）时，式 (14) 右边为正数，t 有正解．令 t 等于式 (14) 右边的值，j 为使 t 取到该值时对应的下标，则 \bar x^{(k)}_j = 0，\bar x^{(k)}_k \gt 0（因为 d^{(k)}_k = 1 \gt 0），对任意 l \in \mathcal N 且 l \ne k 有 \bar x^{(k)}_l = 0（因为 \bar x_l = d^{(k)}_l = 0），故 \bar{\bm x}^{(k)} 有 n 个分量为 0，是一个顶点解．我们可以取 \mathcal B^{(k)} = (\mathcal B \setminus \{j\}) \cup \{k\} 作为解 \bar{\bm x}^{(k)} 对应的基变量下标集合．

若退化，则式 (14) 右边可能为正数或 0．如果是正数则可按照非退化情况处理，否则只能有 t = 0，即 \bar{\bm x}^{(k)} = \bar{\bm x}，而 j 仍需从 \mathcal B 中满足 \bar x_j = 0 且 \bm e_j^{\text{T}} \bm B^{-1} \bm a_k \gt 0 的元素中选择（显然这样的 j 存在）．同样地，\bar{\bm x}^{(k)} 是顶点解，取 \mathcal B^{(k)} = (\mathcal B \setminus \{j\}) \cup \{k\}．注意到，t = 0 时 \bar x_j = 0，而由式 (10) 有 d^{(k)}_j = -\bm e_j^{\text{T}} \bm B^{-1} \bm a_k \lt 0，故 \bm d^{(k)} 不满足式 (8)，易知其逆命题亦成立，这也就是为什么式 (14) 会自动处理退化情形——此时解不移动而基改变．

不论如何，在上述过程中我们都将一个变量从基中剔除并将另一变量加入基中．称 \bar x_j 为这轮迭代的 离基变量 (leaving basic variable)，\bar x_k 为这轮迭代的 进基变量 (entering basic variable)．在一轮迭代中，我们更新了当前解 \bar{\bm x} 和对应的基 \bm B，这样便可进入下一轮迭代．

现在我们来证明若 \bm B 初始时可逆，则它始终可逆．注意定理 3 并不依赖于式 (11)，因而也适用于退化情形．

定理 3 对于一次迭代，设 \mathcal B 为迭代开始时的基变量下标集合，\mathcal B' = (\mathcal B \setminus \{j\}) \cup \{k\}，\bm B 和 \bm B' 分别为对应的基．若 \bm B 可逆，则 \bm B' 可逆．

证明 由 \bm B' 的构造过程得

回忆线性代数中的结论 \det(\bm A + \bm u \bm v^{\text{T}}) = (1 + \bm v^{\text{T}} \bm A^{-1} \bm u) \det \bm A（当 \bm A 可逆时），将其应用到式 (15) 得

由 j 的构造过程以及式 (14) 知 \bm e_j^{\text{T}} \bm B^{-1} \bm a_k \ne 0，又 \bm B 可逆，故式 (16) 给出 \det \bm B' \ne 0，即 \bm B' 可逆．\blacksquare

下面我们介绍这一算法的具体实现．求逆矩阵是计算量很大的工作，而事实上我们也不会真正求出 \bm B^{-1}，而是通过初等行变换进行消元以得到同样的效果．我们将用矩阵的语言具体描述这一算法．我们先讨论 \bm b \geqslant \bm 0 的情况．

首先在式 (1) 中将目标函数改成 Z 并添加一个约束条件 Z - \bm c^{\text{T}} \bm x = 0，于是目标函数的信息便“嵌入”了矩阵之中．结合式 (2)，我们将式 (3) 改写如下：

由于 \bm b \geqslant \bm 0，我们可将所有松弛变量作为初始的 \mathcal B，于是初始时 \bm B = \bm I，\bm x_{\mathcal B} = \bm s = \bm b \geqslant \bm 0．

每次迭代开始时，我们保证矩阵经过列重排（等效地，代码中记录了每一行中基变量所在的列，这样便无需进行列重排）后具有式 (21) 的形式（特别地，初始时式 (17) 亦符合此形式）．令 \bm x_{\mathcal N} = \bm 0，得到基为 \bm B 时的解

若 \bm c_{\mathcal B}^{\text{T}} \bm B^{-1} \bm N - \bm c_{\mathcal N}^{\text{T}} \geqslant \bm 0，由定理 1，式 (18) 为最优解，退出迭代；否则继续寻找更优解．

此时我们需选择进基变量和离基变量，常见的选择方法为：

• 选择 \bm c_{\mathcal B}^{\text{T}} \bm B^{-1} \bm N - \bm c_{\mathcal N}^{\text{T}} 中绝对值最大的负值分量对应的变量作为进基变量，若有多个则任取一个；
• 选择使得式 (14) 取得等号时对应的变量作为离基变量，若不存在则说明解无界，退出迭代，若有多个则任取一个．注意由式 (21) 可得式 (14) 中 \bm x_{\mathcal B} = \bm B^{-1} \bm b，\bm B^{-1} \bm a_k 即为 \bm B^{-1} \bm N 中 \bar x_k 所在的列．

这种选择方法在非退化情况下表现良好，但在退化情况下算法可能不会终止．Bland 规则可解决这一问题（见定理 5），但效率低下，因此实际应用中一般不会使用 Bland 规则．

作如下的变量替换：用新的基变量和新的非基变量所在列组成的矩阵替换掉原来的 \bm B 和 \bm N，用 \begin{pmatrix} z & \bm b^{\text{T}} \end{pmatrix}^{\text{T}} 替换掉原来的式 (21) 右边．经过列重排后，矩阵如下：

由定理 3 知 \bm B 可逆．作初等行变换，即在式 (19) 两边左乘恰当的矩阵将 -\bm c_{\mathcal B}^{\text{T}} 元消去（对应求解式 (5) 过程中消去 \bm x_{\mathcal B} 这一操作）：

化简，得

于是我们便可进入下一轮迭代．

现在我们来证明单纯形法的正确性．

定理 4 (单纯形法正确性) 若每次迭代均成立式 (11)（即不退化），则单纯形法在有限步内终止，且终止时或报告原问题解无界，或给出原问题的一组最优解．

证明 先证明每次迭代后 Z 严格增大．若迭代尚未终止，即式 (13) 不满足，则存在 k \in \mathcal N 使得 k 满足式 (12)．我们选择这样的一个 \bar x_k 作为进基变量，并由式 (14) 计算移动距离 t 和新解 \bar{\bm x}^{(k)} = \bar{\bm x} + t \bm d^{(k)} 或判断出解无界（定理 2）．由 \bm x_{\mathcal B} \gt \bm 0 知 t \gt 0，并由式 (9) 知 \bar{\bm c}^{\text{T}} \bm d^{(k)} \gt 0．于是设 \bar{\bm x}^{(k)} 对应的目标函数值为 Z^{(k)}，则有 Z^{(k)} = \bar{\bm c}^{\text{T}} (\bar{\bm x} + t \bm d^{(k)}) = Z + t \bar{\bm c}^{\text{T}} \bm d^{(k)} \gt Z，即 Z 严格增大．

再证明不存在被重复访问的基．每组基由其下标集合 \mathcal B \in \{1, 2, \ldots, n + m\} 唯一确定，而 |\mathcal B| = m，故不同的基的总数至多为 n + m \choose m，是一个有限数．Z 在迭代过程中严格增大，而每组基 B 由式 (6) 唯一确定一个 Z 值，故同一组基不可能在迭代序列中出现两次．

因此，算法至多访问 n + m \choose m 组不同的基，即算法一定在至多 n + m \choose m 次迭代后终止．终止时或判断出解无界，或满足式 (13) ——即找到最优解（定理 1）．\blacksquare

定理 4 的证明过程提示我们，单纯形法的时间复杂度在最坏情况是指数级的．

定理 5 Bland 规则可保证在退化时算法也能有限步内终止．

证明 见 Bland (1977)．\blacksquare

::::info[\bm b \geqslant \bm 0 时单纯形法的 Rust 代码]

#[derive(Debug, Clone)]
struct LPData {
    /// Number of variables, including slacks.
    n: usize,
    /// Number of constraints.
    m: usize,
    /// Constraints coefficients and right-hand sides.
    cons: AugMat,
    /// Negated objective coefficients.
    obj: Vec<f64>,
    /// Right-hand side of the objective function.
    rhs: f64,
    /// bv[i] is the index of the basic variable for the i-th constraint.
    bv: Vec<usize>,
}

impl LPData {
    fn new(n: usize, m: usize, cons: AugMat, obj: Vec<f64>, rhs: f64, bv: Vec<usize>) -> Self {
        assert!(n >= 1 && m >= 1);
        assert_eq!(cons.m, m);
        assert_eq!(cons.n, n);
        assert_eq!(obj.len(), n);
        assert_eq!(bv.len(), m);
        LPData {
            n,
            m,
            cons,
            obj,
            rhs,
            bv,
        }
    }

    fn pivot_col(&self) -> Option<usize> {
        Some(
            self.obj
                .iter()
                .enumerate()
                .min_by(|x, y| x.1.partial_cmp(y.1).unwrap())
                .filter(|x| *x.1 < -EPS)?
                .0,
        )
    }

    fn pivot_row(&self, col: usize) -> Option<usize> {
        Some(
            (0..self.m)
                .filter_map(|i| {
                    if self.cons.data[i][col] > EPS {
                        Some((i, self.cons.aug[i] / self.cons.data[i][col]))
                    } else {
                        None
                    }
                })
                .min_by(|x, y| x.1.partial_cmp(&y.1).unwrap())?
                .0,
        )
    }

    fn change_basis(&mut self, row: usize, col: usize) {
        let pivot = self.cons.data[row][col];
        let obj_col = self.obj[col];

        self.bv[row] = col;

        self.cons.scale(row, 1.0 / pivot);

        let scalars = (0..self.m)
            .map(|i| -self.cons.data[i][col])
            .collect::<Vec<_>>();
        scalars.into_iter().enumerate().for_each(|(i, x)| {
            if i != row {
                self.cons.add_to(i, row, x)
            }
        });

        self.cons
            .add_to_vec(&mut self.obj, &mut self.rhs, row, -obj_col);
    }

    fn run(&mut self) -> Option<()> {
        while let Some(col) = self.pivot_col() {
            self.change_basis(self.pivot_row(col)?, col);
        }
        Some(())
    }

    fn extract(&self) -> LPResult {
        let mut sln = vec![0.0; self.n];
        self.bv.iter().enumerate().for_each(|(i, &j)| {
            sln[j] = self.cons.aug[i];
        });
        LPResult::Optimal { ans: self.rhs, sln }
    }

    fn optimize(mut self) -> LPResult {
        match self.run() {
            Some(()) => self.extract(),
            None => LPResult::Unbounded,
        }
    }
}

::::

当 \bm b \geqslant \bm 0 不成立时，我们无法找到平凡的初始可行解，需引入两阶段法 (two-phase method) 以求解．它首先作一个辅助线性规划问题以求得一初始可行解，然后在此基础上求解原问题．

首先进行预处理：将式 (3) 中所有满足 b_i \lt 0 的约束条件两边同乘 -1．预处理后有 \bm b \geqslant \bm 0．

第一阶段：引入 m 个 人工变量 (artificial variable) \bm r \geqslant \bm 0（等效地，代码中只对 b_i \lt 0 的约束条件引入人工变量），求解如下线性规划问题：

这里 \bm b \geqslant \bm 0，因此可以取所有人工变量作为初始的 \mathcal B，于是初始时 \bm B = \bm I，\bm x_{\mathcal B} = \bm r = \bm b \geqslant \bm 0．在此基础上，我们便可对式 (22) 施 \bm b \geqslant \bm 0 时的单纯形法．

定理 6 式 (3) 有可行解当且仅当式 (22) 的最优解 W^* = 0．

证明 若式 (3) 有可行解 \bar{\bm x}^*，则 \begin{pmatrix} (\bar{\bm x}^*)^{\text{T}} & \bm 0^{\text{T}}\end{pmatrix}^{\text{T}} 是式 (22) 的可行解，对应的目标函数值为 0，另一方面 \bm r \geqslant \bm 0，故 W = -\bm 1^{\text{T}} \bm r \leqslant 0，于是 W^* = 0．反之，若 W^* = 0，则最优解 \begin{pmatrix} (\bar{\bm x}^*)^{\text{T}} & (\bm r^*)^{\text{T}}\end{pmatrix}^{\text{T}} 满足 -\bm 1^{\text{T}} \bm r^* = 0，由 \bm r^* \geqslant \bm 0 知 \bm r^* = \bm 0，从而 \bar{\bm A} \bar{\bm x}^* = \bm b．又因为 \bar{\bm x}^* \geqslant \bm 0，故 \bar{\bm x}^* 是式 (3) 的一个可行解．\blacksquare

若 W^* \lt 0，由定理 6 知原问题无可行解，算法终止；否则原问题有可行解，算法继续．但这时的基 \mathcal B 中可能仍含有人工变量（注意此时所有人工变量取值为 0，故此为退化情形）．对于基中每个人工变量 r_i，若存在 k \in \mathcal N 使得 \bm e_i^{\text{T}} \bm B^{-1} \bm a_k \ne 0，则以 \bar x_k 为进基变量、\bar r_i 为离基变量执行一次基变换，将人工变量剔除出基．否则对所有 k \in \mathcal N 都有 \bm e_i^{\text{T}} \bm B^{-1} \bm a_k = 0，即式 (21) 中的 \bm B^{-1} \bm N 的第 i 行全为 0，而经过初等行变换后这一行只含变量 r_i，系数为 1，故这一约束条件化为 r_i = 0，是一个冗余约束，可将其直接删除，并相应地修改 m 的值．经过上述处理后，我们得到一组完全由原变量构成的基，便可前往第二阶段．

第二阶段：以处理后的已不含人工变量的基作为初始基，恢复原目标函数 Z = \bar{\bm c}^{\text{T}} \bar{\bm x}，删去人工变量，对式 (3) 施 \bm b \geqslant \bm 0 时的单纯形法．

下面给出两阶段法正确性的完整陈述．

定理 7 (两阶段法正确性) 若每次迭代均成立式 (11)（即不退化），则两阶段法在有限步内终止，且终止时或报告原问题无可行解，或报告原问题解无界，或给出原问题的一组最优解．若引入 Bland 规则，则前文对退化情况也适用．

证明 第一阶段对式 (22) 运行单纯形法．式 (22) 的目标函数有上界 0，解不可能无界．由定理 4（或定理 5）知第一阶段在有限步内终止．若 W^* \lt 0，由定理 6 知原问题无可行解．若 W^* = 0，经剔除人工变量的处理后，我们得到原问题的一组解和对应的基．第二阶段对式 (3) 运行单纯形法，由定理 4（或定理 5）知亦在有限步内终止，终止时或由定理 2 判定解无界，或由定理 1 给出最优解．\blacksquare

::::info[两阶段法的 Rust 代码]

/// Solves the linear program:
///   maximize obj^T x
///   subject to cons.data[i]^T x <= cons.aug[i] for all i.
/// where n is the number of variables, m is the number of constraints.
pub fn solve_lp(n: usize, m: usize, cons: AugMat, obj: Vec<f64>) -> LPResult {
    let arts = cons.aug.iter().filter(|&&x| x < -EPS).count();
    let mut cons1 = AugMat::new(m, n + m + arts);
    let mut bv1 = vec![0usize; m];
    let mut art_idx = n + m;

    for (i, bv) in bv1.iter_mut().enumerate() {
        if cons.aug[i] < -EPS {
            (0..n).for_each(|j| cons1.data[i][j] = -cons.data[i][j]);
            cons1.aug[i] = -cons.aug[i];
            cons1.data[i][n + i] = -1.0;
            cons1.data[i][art_idx] = 1.0;
            *bv = art_idx;
            art_idx += 1;
        } else {
            (0..n).for_each(|j| cons1.data[i][j] = cons.data[i][j]);
            cons1.aug[i] = cons.aug[i];
            cons1.data[i][n + i] = 1.0;
            *bv = n + i;
        }
    }

    let mut obj1 = std::iter::repeat_n(0.0, n + m)
        .chain(std::iter::repeat_n(1.0, arts))
        .collect::<Vec<_>>();
    let mut rhs1 = 0.0;

    for i in 0..m {
        if cons.aug[i] < -EPS {
            cons1.add_to_vec(&mut obj1, &mut rhs1, i, -1.0);
        }
    }

    let mut lp1 = LPData::new(n + m + arts, m, cons1, obj1, rhs1, bv1);
    lp1.run().expect("Phase 1 should not be unbounded");
    if lp1.rhs < -EPS {
        return LPResult::Infeasible;
    }

    for i in 0..m {
        if lp1.bv[i] >= n + m {
            if let Some(col) = (0..n + m).find(|&j| lp1.cons.data[i][j].abs() > EPS) {
                lp1.change_basis(i, col);
            } else {
                // This constraint is redundant, we can ignore it.
            }
        }
    }

    let mut cons2 = AugMat::new(m, n + m);
    for i in 0..m {
        cons2.data[i].copy_from_slice(&lp1.cons.data[i][..n + m]);
        cons2.aug[i] = lp1.cons.aug[i];
    }

    let mut obj2 = obj
        .iter()
        .map(|x| -*x)
        .chain(std::iter::repeat_n(0.0, m))
        .collect::<Vec<_>>();
    let mut rhs2 = 0.0;

    for i in 0..m {
        let col = lp1.bv[i];
        if col < n + m {
            let x = -obj2[col];
            lp1.cons.add_to_vec(&mut obj2, &mut rhs2, i, x);
        }
    }

    let lp2 = LPData::new(n + m, m, cons2, obj2, rhs2, lp1.bv);

    let mut result = lp2.optimize();
    if let LPResult::Optimal { sln, .. } = &mut result {
        sln.truncate(n);
    }
    result
}

::::

于是我们便可给出大 M 法的正确性陈述．

定理 9 (大 M 法正确性) 若每次迭代均成立式 (11)（即不退化），则大 M 法在有限步内终止，且终止时或报告原问题无可行解，或报告原问题解无界，或给出原问题的一组最优解．若引入 Bland 规则，则前文对退化情况也适用．

证明 对式 (23) 施 \bm b \geqslant \bm 0 时的单纯形法，由定理 4（或定理 5）知算法在有限步内终止，终止时算法或判定原问题无可行解，或判定原问题解无界，或给出原问题的一组最优解（定理 8）．\blacksquare

不难发现，两阶段法和大 M 法本质上是同一算法的两种不同实现方式．两者的逻辑完全相同——先最小化人工变量的和以求得一个可行解，再在此基础上最大化原目标函数以求得最优解（或判定解无界）．大 M 法实际上是将两阶段法的两个阶段合为一个问题进行求解．

作为本文的结尾，我们给出线性规划基本定理 (fundamental theorem of linear programming) 的证明．

定理 10 (线性规划基本定理) 对于式 (1)，若其有可行解，则其必有可行顶点解；若其有最优解，则其必有顶点解使得它最优．

证明 由前文知式 (1) 的可行解集与式 (3) 的可行解集存在一一对应关系，且对应的目标函数值相同．以下就式 (3) 进行论述．

若式 (3) 有可行解，由定理 6 知式 (22) 的最优解 W^* = 0，由定理 7，引入 Bland 规则的两阶段法的第一阶段必可在有限步内终止．由 W^* = 0，第一阶段不会报告原问题无可行解．此后经剔除人工变量的处理后，我们得到原问题的一组顶点解和对应的基．

若式 (3) 有最优解，由定理 7 知引入 Bland 规则的两阶段法一定会在有限步内终止，并得到三种结果中的某一种，下面我们排除算法报告原问题无可行解或解无界的两种情况．假设算法报告无可行解，由于式 (3) 有最优解蕴含式 (3) 有可行解，故可由上文可得出矛盾．假设算法报告解无界，则由定理 2 可知存在 \bar{\bm x}(t) = \bar{\bm x} + t \bm d^{(k)} 使得对应的 Z(t) \to +\infty，与式 (3) 有（有限）最优解矛盾．综上，算法会给出一满足式 (13) 的顶点解，由定理 1 知该顶点解最优．\blacksquare

利用大 M 法（定理 9）亦可证明定理 10，留给读者．

值得注意的是，本文对单纯形法的论述中，我们并未使用任何关于线性规划问题的几何性质的论断，而是完全基于算法的代数性质进行分析的．因此，定理 10 的证明也不依赖于任何几何性质（例如可行域的凸性），这和常规的证明方法不同．

参考文献：
清华大学工业工程系《运筹学（确定性方法）》（2025-2026 春季学期）课程讲义
Introduction to Operations Research (11th Edition) by F. S. Hillier and G. J. Lieberman, 2021
New Finite Pivoting Rules for the Simplex Method by R. Bland, 1977