Seventh Mist

On Simplex Method

线性规划是一类非常常见的问题, 许多最优化问题都可以规约到线性规划.

~其实, 单纯形是骗分利器……网络流不会建模咋办? 单纯形强上呀.~

标准形式

下面设未知数个数为nn, 约束条件个数(AA的行数)为mm.

Maximize

cTx\displaystyle \mathbf{c}^T \mathbf{x}

s.t.

Axb\displaystyle A \mathbf{x} \le \mathbf{b} x0\displaystyle \mathbf{x} \ge \mathbf{0}

松弛形式

Maximize

cTx\displaystyle \mathbf{c}^T \mathbf{x}

s.t.

y=bAx0\displaystyle \mathbf{y} = \mathbf{b} - A \mathbf{x} \ge 0

在实践中我们一般不区分x\mathbf{x}y\mathbf{y}, 而把x\mathbf{x}写作x1n\mathbf{x}_{1 \ldots n}, 把y\mathbf{y}写作xn+1n+m\mathbf{x}_{n + 1\ldots n + m}.

单纯形法

一些定义

称松弛型中, 写在等式左边的为基(本)变量(如上面的y\mathbf{y}), 其余的为非基(本)变量(如上面的x\mathbf{x}).

算法流程

先假定我们找到了一组初始解x\mathbf{x}, 然后尝试通过迭代得到最优解.

我们进行如下的过程:

  1. 首先找到任意一项满足ci>0\mathbf{c}_i \> 0ii;
  2. 然后我们枚举所有的约束, 然后找到在满足非负约束的条件下xi\mathbf{x}_i至多能增加多少, 设约束矩阵的第jj行给出的约束导致xi\mathbf{x}_i至多能增加dd;
  3. xi\mathbf{x}_i为主元反写第jj行, 并代入到目标函数和约束矩阵的其他行中, 从而使得xi\mathbf{x}_i成为了基变量, 而原来第jj行的基变量则变成了非基变量;
  4. 回到1, 直到找不到满足条件的ii.

算法的正确性以及复杂度证明可以参考CLRS.

虽然单纯形的最坏复杂度是指数级的, 但在OI中基本上只要存得下, 就能跑得出.

找一组可行的初始解

如何找到一组可行的初始解呢? 我们可以尝试增加一个变量yy, 然后做一个新的线性规划:

Maximize

z\displaystyle -z

s.t.

z0\displaystyle z \ge 0 x0\displaystyle \mathbf{x} \ge \mathbf{0} y=bAx+z0\displaystyle \mathbf{y} = \mathbf{b} - A \mathbf{x} + z \ge \mathbf{0}

如果z=0-z = 0, 那么原问题有解, 否则无解.

而解这个线性规划还是比较简单的.

判断是否无界

如果在算法的步骤2中无法找到xi\mathbf{x}_i增长的一个上界, 那么这个线性规划问题就是无界的.

代码实现

[lang: cpp]
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
#include <cmath>
#include <cstdio>
#include <algorithm>
const int N = 50 + 10;
const double eps = 1e-9;
int n, m;
double a[N][N], b[N], c[N];
int row[N], col[N];
void pivot(int x, int y, int ext = -1) {
std::swap(row[x], col[y]);
double temp = -1. / a[x][y];
a[x][y] = -1.;
for (int i = 0; i <= n; ++i) a[x][i] *= temp;
for (int i = 1; i <= m + 1 || i <= ext; ++i) {
if (i == x || fabs(a[i][y]) < eps || (i > m + 1 && i != ext)) continue;
temp = a[i][y];
a[i][y] = 0.;
for (int j = 0; j <= n; ++j) a[i][j] += temp * a[x][j];
}
}
bool simplex(int ext = -1) {
while (1) {
int p = 0;
for (int i = 1; i <= n; ++i) if (a[m + 1][i] > eps && (!p || col[i] < col[p])) p = i;
if (!p) break;
int k = 0;
double cur;
for (int i = 1; i <= m; ++i) {
if (a[i][p] > -eps) continue;
double temp = -a[i][0] / a[i][p];
if (!k || temp + eps < cur || (fabs(temp - cur) < eps && row[i] < row[k])) cur = temp, k = i;
}
if (k) pivot(k, p, ext); else return false;
}
return true;
}
int main() {
int type;
scanf("%d%d%d", &n, &m, &type);
++n;
for (int i = 1; i < n; ++i) scanf("%lf", &a[m + 3][i]);
for (int i = 1; i <= m; ++i) {
for (int j = 1; j < n; ++j) scanf("%lf", &a[i][j]), a[i][j] = -a[i][j];
scanf("%lf", &a[i][0]);
}
for (int i = 1; i <= m; ++i) a[i][n] = 1.;
for (int i = 1; i <= n; ++i) col[i] = i;
for (int i = 1; i <= m + 2; ++i) row[i] = n + i;
a[m + 1][n] = -1.;
int k = 1;
for (int i = 2; i <= m; ++i) if (a[i][0] < a[k][0]) k = i;
pivot(k, n, m + 3);
simplex(m + 3);
if (a[m + 1][0] < -eps) return puts("Infeasible"), 0;
for (int i = 1; i <= n; ++i) a[m + 2][i] = -a[m + 1][i];
m += 2;
if (!simplex()) return puts("Unbounded"), 0;
printf("%.12f\n", a[m + 1][0]);
if (type) {
static int pos[N];
for (int i = 1; i <= m; ++i) if (row[i] < n) pos[row[i]] = i;
for (int i = 1; i < n; ++i) printf("%.12f ", pos[i] ? a[pos[i]][0] : 0.);
}
return 0;
}