单纯形法的分析和实现

单纯形法是求解一类被称为“线性规划(LP)”问题的通用方法。高二学习的一类含有两个变量的简单线性规划可以通过在平面上作图得出，同理三个变量可以投到空间里计算，但超过四个变量呢？单纯形法即是对这种方法的推广，其过程类似于在n维几何体的顶点上不断尝试，但由于高维几何体过于玄学，大多数资料都使用了更直观的代数角度阐释。

一个线性规划的标准型被定义为：

M a x : \sum i = 1 m c i x i S . T . : a 11 x 1 + a 12 x 2 + \dots + a 1 m x m \leq b 1 a 21 x 1 + a 22 x 2 + \dots + a 2 m x m \leq b 2 \dots a n 1 x 1 + a n 2 x 2 + \dots + a n m x m \leq b n x i \geq 0

由于不等式 x≤y 可以通过引入一个松弛变量 x′≥0，写做 x=y−x′，因而不等式可以被等式替代称为松弛型

用矩阵表示为：

M a x : c X S . T . : A X = B x i \geq 0

而单纯形法的思路是：每次将一个“可以增加值的”变量的值增加到临界点，直到所有变量都不能增加值。例如线性规划：

z = 5 x 1 - 3 x 2 x 1 - x 2 + x 3 = 1 2 x 1 + x 2 + x 4 = 2

(

x3,x4是松弛变量)

则当 当前的非松弛变量x1,x2 取0时z的值为一个可行值。

观察到z=5x1−3x2，因而x1增大可以使得z增大。而(2)(3)式中x1增大的 贡献能力 即等式右边数字与x1的系数之比一定，随便挑一个即可。例如挑(2)：

5 x 1 - 3 x 2 = z (1) x 1 - x 2 + x 3 = 1 (2) 2 x 1 + x 2 + x 4 = 2 (3)

选定x1为 替入变量 ，选定(2)因而x3为替出变量。则用初等行变换消去(1)(3)中的x1：

- 2 5 x 2 + x 3 = 1 - 1 5 z (1) x 1 - x 2 + x 3 = 1 (2) 3 x 2 + x 4 - 2 x 3 = 0 (3)

整理得到：

z = 5 + 2 x 2 - 5 x 3 (1) x 1 = 1 + x 2 - x 3 (2) x 4 = 2 x 3 - 3 x 2 (3)

这时当 当前的非松弛变量x2,x3 都取0时，z = 5。注意到(3)指出 Δx3:Δx2=3:2，代入(1)后显然不能使z更大。因而不会继续转动。

是不是一次就能取得最大值呢？不是。因此我们需要多次转动直到找不到更大的值。事实上，转动次数的上界是指数级的。

代码实现

#include <bits/stdc++.h>
using namespace std;

const double eps = 1e-9;
const double inf = 1e100;

inline int sgn(double d)
{ return abs(d) <= eps ? 0 : (d > 0 ? 1 : -1); }

const int maxn = 1005, maxm = 1005;
int n, m; // 等式数，变量数(不含松弛变量)
int A[maxn][maxm];

void pivot(int l, int e) // 一次转动, l为替出变量，e为替入变量
{
    cout << l << " " << e << endl;
    for (int i = 1; i <= n; i++) if (i != l && sgn(A[i][e])) {
        for (int j = 1; j <= m; j++)
            if (j != e) A[i][j] -= A[l][j] * A[i][e]/A[l][e];
        A[i][e] /= -A[l][e];
    }
    for (int i = 1; i <= m; i++) if (i != e) A[l][i] /= A[l][e]; A[l][e] = 1/A[l][e];
}

double simplex()
{
    for (int e, l; ; ) {
        double t = inf;
        for (e = 1; e < m; e++) if (sgn(A[n][e]) > 0) break;
        if (e == m) return -A[n][m];
        for (int i = 1; i < n; i++)
            if (sgn(A[i][e]) > 0 && t > A[i][m]/A[i][e]) t = A[l=i][m]/A[i][e];
        if (t == inf) return inf;
        pivot(l, e);
    }
}

int main()
{
    cin >> n >> m;
    for (int i = 1; i <= n; i++)
        for (int j = 1; j <= m; j++)
            cin >> A[i][j];
    printf("%g", simplex());
    return 0;
}

把问题表述成线性规划

一个经典的例子是最短路径，详细介绍在《算法导论》上有，在此不再赘述。给出一个利用单纯形法求解最短路径的代码：


#include <bits/stdc++.h>
using namespace std;

const double eps = 1e-9;
const double inf = 1e100;

inline int sgn(double d)
{ return abs(d) <= eps ? 0 : (d > 0 ? 1 : -1); }

const int maxn = 10005, maxm = 1005;
int n, m; // 等式数，变量数(不含松弛变量)
int A[maxn][maxm], B[maxn][maxm];

void pivot(int l, int e) // 一次转动, l为替出变量，e为替入变量
{
    //cout << l << " " << e << endl;
    for (int i = 1; i <= n; i++) if (i != l && sgn(A[i][e])) {
        for (int j = 1; j <= m; j++)
            if (j != e) A[i][j] -= A[l][j] * A[i][e]/A[l][e];
        A[i][e] /= -A[l][e];
    }
    for (int i = 1; i <= m; i++) if (i != e) A[l][i] /= A[l][e]; A[l][e] = 1/A[l][e];
}

double simplex()
{
    for (int e, l; ; ) {
        double t = inf;
        for (e = 1; e < m; e++) if (sgn(A[n][e]) > 0) break;
        if (e == m) return -A[n][m];
        for (int i = 1; i < n; i++)
            if (sgn(A[i][e]) > 0 && t > A[i][m]/A[i][e]) t = A[l=i][m]/A[i][e];
        if (t == inf) return inf;
        pivot(l, e);
        // cout << "Get : " << -A[n][m] << endl;
    }
}

int main()
{
    int v, e, s, t;
    cin >> v >> e >> s;
    memset(A, 0, sizeof A);
    n = 0, m = v+1;
    for (int i = 1; i <= e; i++) {
        int a, b; double c;
        scanf("%d%d%lf", &a, &b, &c);
        A[++n][b] = 1, A[n][a] = -1, A[n][m] = c;
    }
    A[++n][s] = 1; ++n;
    for (int i = 1; i <= n; i++) for (int j = 1; j <= m; j++) B[i][j] = A[i][j];
    for (int i = 1; i <= v; i++) {
        for (int j = 1; j <= n; j++) for (int k = 1; k <= m; k++) A[j][k] = B[j][k];
        A[n][i-1] = 0, A[n][i] = 1;
        double ans = simplex();
        if (ans == inf) printf("2147483647 ");
        else printf("%g ", ans);
    }
    return 0;
}