单纯形法（simplex algorithm）适用于求解如下的含 n 个变量和 m 个约束条件的线性规划标准型问题：
\begin{aligned}
\max&\quad Z = \bm c^{\text{T}} \bm x \\
\text{s.t.}&\quad \bm A \bm x \le \bm b \\
&\quad \bm x \ge \bm 0.
\end{aligned} \tag{1}
称满足式 (1) 中的两个约束条件的 \bm x 为式 (1) 的一个 可行解（feasible solution）．
从几何直观上讲，\bm A \bm x \le \bm b 和 \bm x \ge \bm 0 围成了一个 \mathbb R^n 中的凸集，最优解（optimal solution）（若存在）一定可以在凸集的顶点处取得．单纯形法以凸集的一个顶点作为初始解，每次迭代从当前解所在顶点的相邻顶点中选择一个更优的作为新的解，直到所有相邻顶点都不更优为止．单纯形法的大致流程如下：
对于任意一个凸集顶点处对应的解（称为 顶点解（corner-point solution）），其应能使得 \bm A \bm x \le \bm b 和 \bm x \ge \bm 0 这 n + m 个不等式中的 n 个取到等号．若 \bm A \bm x \le \bm b 中的某个不等式取到等号，则对应的松弛变量为 0；若 \bm x \ge \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) 改写如下：
证明. 先证明每次迭代后 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 \ge \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 \ge \bm 0 不成立时，我们无法找到平凡的初始可行解，需引入两阶段法（two-phase method）以求解．它首先作一个辅助线性规划问题以求得一初始可行解，然后在此基础上求解原问题．
/// 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
}
::::
作为本文的结尾，我们给出线性规划基本定理（fundamental theorem of linear programming）的证明．
参考文献：
清华大学工业工程系《运筹学（确定性方法）》（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