单纯形法入门笔记（目前无法通过 UOJ 模板题）
太难了啊。
考虑线性规划的标准型
注意到在 \(n\) 维空间内，每条限制相当于一个 \(n - 1\) 维的平面去切割这个空间，这是凸的，凸的空间和凸的空间的交仍然是凸的。
这是一个很强的性质，也就是说，此时直接进行爬山算法就是正确的（容易想到如果爬不出最大值，解空间就不是凸的）。
prework
具体来说，先转化为松弛型，就是添加 \(m\) 个元 \(x_{n + 1}, x_{n + 2}, \dots, x_{n + m}\)：
记 \(z = \sum\limits_{i = 1} ^ n c_i x_i\)，也就是目标最大值。记初始时 \(S = \{1, 2, \dots, n\}, \ N = \{n + 1, n + 2, \dots, n + m\}\)，以及上面的所有系数为 \(f\)，这样我们就有了 \(m\) 条等式和目标最大值的式子：
下面，我们将 \(S\) 中的元称之为基本元，\(N\) 中的元称之为非基本元。
爬山算法需要一个初始点，我们先默认所有 \(b_i\ge 0\)，这样就有了一组初始解 \(x_1 = x_2 = \dots = x_n = 0\)。在算法过程中，我们始终默认 \(\forall i\in N, \ x_i = 0\)，并时刻保证 \(\forall i\in S, \ f_{i, 0}\ge 0\)。
转动（pivot）
单纯形法的具体爬山方法是进行转动（pivot）。一次转动中：
-
选择一个非基本元 \(x_i\) 和基本元 \(x_j\)，并且满足 \(f_{i, j} \not = 0\)。
-
将式子 \(x_i = f_{i, 0} + \sum\limits_{k \in N} f_{i, k} x_k\) 改写为 \(-f_{i, j} x_j = f_{i, 0} - x_i + \sum\limits_{k \in N \land k \not = j} f_{i, k} x_k\)。
-
将式子除以 \(-f_{i, j}\)，将 \(i\) 从 \(N\) 中移到 \(S\)，将 \(j\) 从 \(S\) 中移到 \(N\)，然后整理得 \(x_j = f_{j, 0} + \sum\limits_{k \in N} f_{j, k} x_k\)。
-
对于 \(p \in S \land p \not = j\) 以及 \(p = 0\)，将 \(f_{p, i}\) 展开。
这样就实现了基本元和非基本元之间的交换。
爬山
当我们发现存在 \(i\in N\) 满足 \(f_{0, i} > 0\) 时，我们希望能尽可能增大 \(x_i\) 的值，于是考虑将 \(x_i\) 通过 pivot 操作换到非基本元。
找到这样一个 \(i\) 后寻找一个能让 \(x_i\) 增大的限制最紧的 \(j(j \in S)\)。具体来讲，对于式子 \(x_j = f_{j, 0} + \sum\limits_{i \in S} f_{j, i} x_i\)，如果 \(f_{j, i}\ge 0\)，那么 \(x_i\) 的增大是没有限制的；如果 \(f_{j, i} < 0\)，那么 \(x_i\) 最大取到 \(-\dfrac {f_{j, 0}} {f_{j, i}}\)。我们找到这个分式最小的 \(j\)，如果所有 \(j\) 都满足 \(f_{j, i}\ge 0\) 说明答案是无界的（可以取到 \(+\infty\)）。
找到这样的 \(j\) 后，将两者 pivot 一下。由于 \(f_{j, i} < 0\)，所以 pivot 后 \(f_{i, 0}\) 一定非负，且 \(f_{0, 0}\) 即答案会增大。
不断重复此过程，即可得到最优解，答案为 \(f_{0, 0}\)。
初始可行解
如果初始存在整数 \(i\in [1, m]\) 满足 \(b_i < 0\)，我们还需要先找一组可行解才能开始爬山。
考虑加入一个人工变量 \(x_{n + m + 1}\)，将 \(n + m + 1\) 加入 \(N\)，最大化 \(z' = -x_{n + m + 1}\)，每条条式子改为 \(x_{i + n} = f_{i + n, 0} + \sum\limits_{j = 1} ^ n f_{i, j} x_j + x_{n + m + 1}\)。
然后找到满足 \(f_{j, 0}\) 最小的 \(j(j\in [n + 1, n + m])\)，此时 \(b_j < 0\)，然后 pivot 一下 \(j\) 和 \(n + m + 1\)。
pivot 之后 \(f_{n + m + 1, 0}\) 就变成了 \(-f_{j, 0}\)，这是一个非负数。然后对于其他 \(k\in [n + 1, n + m] \land k\not = j\)，将 \(f_{k, n + m + 1} = 1\) 展开后 \(f_{k, 0}\) 会加上 \(f_{j, 0}\)，也会变成一个非负数。
先求出这个线性规划的最优解，如果 \(z' = -x_{n + m + 1}\) 取不到 \(0\)，说明无解。
否则，如果 \(n + m + 1 \in N\)，在对应式子中随便找一个 \(i\in N \land f_{n + m + 1, i} \not = 0\)，然后 pivot 一下 \(n + m + 1\) 和 \(i\)。注意这里的正确性，因为此时一定有 \(f_{n + m + 1} = 0\)。
然后，在所有式子中删去 \(x_{n + m + 1}\)，最大化 \(z = \sum\limits_{i = 1} ^ n c_ix_i\) 再求一次线性规划最优解。
Code（无法通过 UOJ 模板题）
点击查看代码
#include <bits/stdc++.h>
#define ll long long
#define LL long long
#define uLL unsigned LL
#define pb push_back
#define pir pair <ll, ll>
#define mkp make_pair
#define fi first
#define se second
#define i128 __int128
using namespace std;
template <class T>
void rd(T &x) {
char ch; bool f = 0;
while(!isdigit(ch = getchar()))
if(ch == '-') f = 1;
x = ch - '0';
while(isdigit(ch = getchar()))
x = (x << 1) + (x << 3) + ch - '0';
if(f) x = -x;
}
const ll maxn = 1e4 + 10, inf = 2e9 + 5, mod = 998244353, iv = mod - mod / 2;
const LL INF = 1e18 + 5;
ll power(ll a, ll b = mod - 2, ll p = mod) {
ll s = 1;
while(b) {
if(b & 1) s = 1ll * s * a % p;
a = 1ll * a * a % p, b >>= 1;
} return s;
}
template <class T1, class T2>
void add(T1 &x, const T2 y) { x = x + y >= mod? x + y - mod : x + y; }
template <class T1, class T2>
void sub(T1 &x, const T2 y) { x = x < y? x + mod - y : x - y; }
template <class T1, class T2>
ll pls(const T1 x, const T2 y) { return x + y >= mod? x + y - mod : x + y; }
template <class T1, class T2>
ll mus(const T1 x, const T2 y) { return x < y? x + mod - y : x - y; }
template <class T1, class T2>
void chkmin(T1 &a, const T2 b) { a = a < b? a : b; }
template <class T1, class T2>
void chkmax(T1 &a, const T2 b) { a = a < b? b : a; }
const double eps = 1e-9;
ll n, m, c[maxn], a[maxn][maxn], b[maxn], id1[maxn], id2[maxn];
double f[maxn][maxn], res[maxn];
ll N, M, op;
void pivot(ll r, ll c) {
double k = -f[r][c];
f[r][c] = -1, swap(id1[r], id2[c]);
for(ll i = 0; i <= N; i++) f[r][i] /= k;
for(ll i = 0; i <= M; i++) {
if(i == r) continue;
double w = f[i][c]; f[i][c] = 0;
for(ll j = 0; j <= N; j++) f[i][j] += w * f[r][j];
}
}
void simplex() {
ll c = 0, r = 0;
for(ll i = 1; i <= N; i++)
if(f[0][i] > eps) { c = i; break; }
if(!c) return;
for(ll i = 1; i <= M; i++)
if(f[i][c] < -eps && (!r || f[r][0] / (-f[r][c]) > f[i][0] / (-f[i][c])))
r = i;
if(!r) puts("Unbounded"), exit(0);
pivot(r, c);
return simplex();
}
void print() {
printf("z = %lf", f[0][0]);
for(ll i = 1; i <= N; i++) printf(" + %lfx%d", f[0][i], id2[i]);
puts("");
for(ll i = 1; i <= M; i++) {
printf("x%d = %lf", id1[i], f[i][0]);
for(ll j = 1; j <= N; j++)
printf(" + %lfx%d", f[i][j], id2[j]);
puts("");
}
}
int main() {
rd(n), rd(m), op = 1;
for(ll i = 1; i <= n; i++) rd(c[i]);
for(ll i = 1; i <= m; i++) {
for(ll j = 1; j <= n; j++) rd(a[i][j]);
rd(b[i]);
}
for(ll i = 1; i <= n; i++) id2[i] = i;
id2[n + 1] = n + m + 1, f[0][n + 1] = -1;
N = n + 1, M = m;
ll pos = 0;
for(ll i = 1; i <= m; i++) {
for(ll j = 1; j <= n; j++) f[i][j] = -a[i][j];
f[i][n + 1] = 1;
f[i][0] = b[i], id1[i] = i + n;
if(b[i] < 0 && (!pos || b[pos] > b[i])) pos = i;
}
if(pos) {
pivot(pos, n + 1);
simplex();
if(f[0][0] < -eps) return puts("Infeasible"), 0;
for(ll i = 1; i <= M; i++)
if(id1[i] == n + m + 1) {
assert(abs(f[i][0]) < eps);
for(ll j = 1; j <= N; j++)
if(abs(f[i][j]) > eps) {
pivot(i, j);
break;
}
break;
}
pos = 1;
while(id2[pos] <= n + m) ++pos;
for(ll i = 0; i <= M; i++) swap(f[i][pos], f[i][n + 1]);
swap(id2[pos], id2[n + 1]);
}
for(ll i = 0; i <= n; i++) f[0][i] = 0;
for(ll i = 1; i <= n; i++)
if(id2[i] <= n) f[0][i] = c[id2[i]];
for(ll i = 1; i <= m; i++)
if(id1[i] <= n) {
for(ll j = 0; j <= n; j++)
f[0][j] += f[i][j] * c[id1[i]];
}
N = n;
simplex();
printf("%.9lf\n", f[0][0]);
if(!op) return 0;
for(ll i = 1; i <= m; i++)
if(id1[i] <= n) res[id1[i]] = f[i][0];
for(ll i = 1; i <= n; i++) printf("%.9lf ", res[i]); puts("");
return 0;
}
浙公网安备 33010602011771号