引言
线性规划（Linear Programming, LP）是运筹学中最基础也是应用最广泛的优化工具，从生产计划到投资组合优化，从物流配送到机器学习中的支持向量机，都能看到它的身影。单纯形法（Simplex Method）由George Dantzig于1947年提出，是求解线性规划问题的经典算法。尽管存在多项式时间的内点法，单纯形法在实际工程中仍然因其高效性和直观性而被广泛使用。本文将用Java完整实现单纯形表算法，并讲解对偶理论的工程意义。
问题建模
标准型线性规划问题定义为：
最大化目标函数：c^T x
约束条件：Ax <= b, x >= 0
其中A是m x n的约束矩阵，b是m维资源向量，c是n维价值向量。
松弛变量引入
为了使用单纯形法，我们需要将不等式约束转化为等式约束。引入松弛变量s >= 0：
Ax + s = b
这样，初始基变量就可以选择所有松弛变量，构成单位矩阵作为初始可行基。
算法原理
单纯形表结构
单纯形表是一个(m+1) x (n+m+1)的矩阵，包含：
- 第0行：目标函数系数（检验数）
- 第1~m行：约束等式系数
- 最后一列：右端项b的值
核心步骤
- 初始基可行解：选择松弛变量作为初始基变量
- 最优性检验：若所有检验数（reduced cost）非正，则当前解最优
- 入基变量选择：选择检验数最大的非基变量入基（保证目标函数严格增加）
- 出基变量选择：使用最小比值法则（Minimum Ratio Test）确定出基变量
- 转轴运算：高斯消元更新单纯形表
- 重复步骤2-5直至最优
对偶理论
每个线性规划问题（原问题）都有一个对应的对偶问题。对偶定理告诉我们：原问题的最优值等于对偶问题的最优值。这一性质在经济学中被称为”影子价格”，在机器学习中对应着SVM的KKT条件。
Java完整实现
import java.util.*;
/**
* 单纯形法求解标准型线性规划问题
* 标准型: max c^T * x, 约束 Ax <= b, x >= 0
*/
public class SimplexMethod {
private static final double EPS = 1e-9;
private static final double INF = 1e18;
private final int m; // 约束数量（行数）
private final int n; // 原始变量数量
private final double[][] a; // 单纯形表: (m+1) x (n+m+1)
private final int[] basis; // 基变量索引，basis[i]表示第i行对应的基变量
private final int[] nonBasis; // 非基变量索引
/**
* 构造单纯形法求解器
* @param A 约束矩阵 m x n
* @param b 右端项向量 m
* @param c 目标函数系数 n
*/
public SimplexMethod(double[][] A, double[] b, double[] c) {
this.m = A.length;
this.n = A[0].length;
// 单纯形表: m个约束 + 1个目标函数行
// 列: n个原始变量 + m个松弛变量 + 1个右端项
this.a = new double[m + 1][n + m + 1];
this.basis = new int[m];
this.nonBasis = new int[n];
// 初始化非基变量（原始变量）
for (int j = 0; j < n; j++) {
nonBasis[j] = j;
}
// 初始化基变量（松弛变量）
for (int i = 0; i < m; i++) {
basis[i] = n + i;
}
// 填充约束行
for (int i = 0; i < m; i++) {
for (int j = 0; j < n; j++) {
a[i][j] = A[i][j];
}
// 松弛变量构成单位矩阵
a[i][n + i] = 1.0;
a[i][n + m] = b[i];
}
// 填充目标函数行（注意：单纯形表中存的是负的检验数）
for (int j = 0; j < n; j++) {
a[m][j] = -c[j];
}
a[m][n + m] = 0.0; // 初始目标函数值
}
/**
* 执行单纯形法迭代求解
* @return 最优目标函数值，若问题无界返回INF
*/
public double solve() {
while (true) {
// 步骤1: 寻找入基变量（检验数最大的非基变量）
int enter = -1;
double maxReducedCost = -EPS;
for (int j = 0; j < n + m; j++) {
// 只考虑非基变量
if (!isBasis(j) && a[m][j] > maxReducedCost) {
maxReducedCost = a[m][j];
enter = j;
}
}
// 所有检验数 <= 0，当前解最优
if (enter == -1) {
break;
}
// 步骤2: 寻找出基变量（最小比值法则）
int leave = -1;
double minRatio = INF;
for (int i = 0; i < m; i++) {
if (a[i][enter] > EPS) {
double ratio = a[i][n + m] / a[i][enter];
if (ratio < minRatio - EPS ||
(Math.abs(ratio - minRatio) < EPS && basis[i] < basis[leave])) {
minRatio = ratio;
leave = i;
}
}
}
// 若没有任何约束限制入基变量，则问题无界
if (leave == -1) {
return INF;
}
// 步骤3: 转轴运算（高斯消元）
pivot(leave, enter);
}
return a[m][n + m];
}
/**
* 转轴运算：以a[leave][enter]为主元进行高斯消元
*/
private void pivot(int leave, int enter) {
double pivotVal = a[leave][enter];
// 主元行归一化
for (int j = 0; j <= n + m; j++) {
a[leave][j] /= pivotVal;
}
// 消去其他行的enter列
for (int i = 0; i <= m; i++) {
if (i != leave && Math.abs(a[i][enter]) > EPS) {
double factor = a[i][enter];
for (int j = 0; j <= n + m; j++) {
a[i][j] -= factor * a[leave][j];
}
}
}
// 更新基变量和非基变量集合
int oldBasis = basis[leave];
basis[leave] = enter;
// 找到enter原来在非基中的位置，替换为oldBasis
for (int j = 0; j < n; j++) {
if (nonBasis[j] == enter) {
nonBasis[j] = oldBasis;
break;
}
}
}
/**
* 判断变量是否在基中
*/
private boolean isBasis(int var) {
for (int b : basis) {
if (b == var) return true;
}
return false;
}
/**
* 获取最优解中各变量的取值
* @return 所有变量（包括松弛变量）的取值数组
*/
public double[] getSolution() {
double[] x = new double[n + m];
Arrays.fill(x, 0.0);
for (int i = 0; i < m; i++) {
x[basis[i]] = a[i][n + m];
}
return x;
}
/**
* 获取原始变量的最优解
* @return 长度为n的解向量
*/
public double[] getOriginalSolution() {
double[] sol = getSolution();
double[] result = new double[n];
System.arraycopy(sol, 0, result, 0, n);
return result;
}
/**
* 打印单纯形表（调试用）
*/
public void printTableau() {
System.out.println("单纯形表：");
for (int i = 0; i <= m; i++) {
System.out.printf("Row %d: ", i);
for (int j = 0; j <= n + m; j++) {
System.out.printf("%8.3f ", a[i][j]);
}
System.out.println();
}
System.out.println("基变量: " + Arrays.toString(basis));
}
public static void main(String[] args) {
// 测试案例: max z = 3x1 + x2 + 2x3
// 约束:
// x1 + x2 + 3x3 <= 30
// 2x1 + 2x2 + 5x3 <= 24
// 4x1 + x2 + 2x3 <= 36
// x1, x2, x3 >= 0
double[][] A = {
{1, 1, 3},
{2, 2, 5},
{4, 1, 2}
};
double[] b = {30, 24, 36};
double[] c = {3, 1, 2};
SimplexMethod simplex = new SimplexMethod(A, b, c);
double optimal = simplex.solve();
System.out.println("最优目标函数值: " + optimal);
double[] solution = simplex.getOriginalSolution();
System.out.println("最优解: x1=" + solution[0] +
", x2=" + solution[1] +
", x3=" + solution[2]);
}
}
算法扩展：两阶段法
上述实现假设原点（所有决策变量为0）是可行解。对于包含”>=”或”=”约束的问题，需要使用两阶段单纯形法：
- 第一阶段：引入人工变量，构造辅助问题最小化人工变量之和，寻找初始可行基
- 第二阶段：以第一阶段得到的可行基为起点，求解原问题
对偶理论的工程意义
| 原问题 | 对偶问题 | 经济学解释 |
|---|---|---|
| max c^T x | min b^T y | 利润最大化 ↔ 成本最小化 |
| 约束Ax <= b | 变量y >= 0 | 资源限制 ↔ 影子价格 |
| 变量x >= 0 | 约束A^T y >= c | 产品产量 ↔ 机会成本 |
在资源分配问题中，对偶变量y_i表示第i种资源的影子价格——即该资源增加一个单位时目标函数的改进量。
复杂度与局限性
| 特性 | 说明 |
|---|---|
| 最坏时间复杂度 | 指数级（Klee-Minty构造型例） |
| 平均时间复杂度 | 多项式级，实际表现极优 |
| 空间复杂度 | O(mn) |
| 数值稳定性 | 需处理浮点精度问题（本实现使用EPS） |
总结
单纯形法不仅是一个算法，更是一种思维方式——通过在可行域的顶点之间移动来寻找最优解。理解单纯形法后，读者可以进一步学习：
- 对偶单纯形法：从对偶可行但原问题不可行的解出发
- 内点法：Karmarkar算法，具有多项式时间复杂度
- 整数线性规划：分支定界法与割平面法
- 网络单纯形法：针对网络流问题的特殊优化
单纯形法的核心思想——基变换与转轴运算——也是许多现代优化算法的基础组件。
