在 Codeforces 有人问了我两个优化问题，一个是非线性规划（具体说是凸优化问题）。一般来说非线性规划没有什么具体的算法，但是凸优化，可以转化成凸包，然后转换成线段（如果是二维的）上的最值问题就搞定了。另一个是线性规划（所有线性规划其实也是特殊的凸优化），线性规划有著名的单纯形算法，7 年前就学过，但一直没写过代码。趁这次机会重新学习一下单纯形算法，并给出代码。参考教材：运筹学 第三版 清华大学出版社，另外 OI-wiki 上讲的简洁清晰但是不够全面。

Gauss 消元法

之前一直不写这个模板的原因：可能无解，可能唯一解，可能无穷多个解，double 有判断，很烦。

求解 $Ax = b$，如果无解就输出空向量，否则输出（某一个）答案向量，无穷解的话随便输出一个。

浮点数版

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
72
73
74
75
76
77
78
79
80
#include <bits/stdc++.h>
#define watch(x) std::cout << (#x) << " is " << (x) << std::endl
using LL = long long;
using pii = std::pair<int, int>;
using pll = std::pair<LL, LL>;

std::vector<double> Gauss(std::vector<std::vector<double>> A, std::vector<double> b) {
	int n = A.size(), m = A[0].size();
	std::vector<double> x(m);
	std::vector<int> p(m);
	std::iota(p.begin(), p.end(), 0);
	const double eps = 1e-12;
	auto findNonZero = [&](int i) { // 实际上找最大的比较好
		for (int row = i; row < n; ++row) if (fabs(A[row][i]) > eps) return row;
		return n;
	};
	auto triangleGauss = [&](int sz) { // A[i][i] = 1
		std::vector<double> x(sz);
		for (int i = sz - 1; i >=0; --i) {
			x[i] = b[i];
			for (int row = 0; row < i; ++row) b[row] -= A[row][i] * x[i];
		}
		x.resize(A[0].size());
		return x;
	};
	int sz = n;
	for (int i = 0, row; i < n; ++i) {
		while (i < m) {
			row = findNonZero(i);
			if (row != n) break;
			for (int j = 0; j < n; ++j) A[j][i] = A[j][m - 1];
			std::swap(p[i], p[--m]);
		}
		if (i == m) {
			for (int row = m; row < n; ++row) if (fabs(b[row])) {
				// std::cout << "\nNo answer\n";
				return std::vector<double>();
			}
			sz = i;
			break;
		}
		if (row != i) {
			std::swap(A[row], A[i]);
			std::swap(b[row], b[i]);
		}
		b[i] /= A[i][i];
		for (int j = m - 1; j >= i; --j) A[i][j] /= A[i][i];
		for (int row = i + 1; row < n; ++row) {
			b[row] -= A[row][i] * b[i];
			for (int j = m - 1; j >= i; --j) {
				A[row][j] -= A[row][i] * A[i][j];
			}
		}
	}
	// if (sz != A[0].size()) std::cout << "\nInfinite answer\n";
	auto xt = triangleGauss(sz);
	for (int t = 0; t < A[0].size(); ++t) x[p[t]] = xt[t];
	return x;
}

int main() {
	// freopen("in","r",stdin);
	std::ios::sync_with_stdio(false);
	std::cin.tie(nullptr);
	int cas = 1;
	std::cin >> cas;
	while (cas--) {
		int n, m;
		std::cin >> n >> m;
		std::vector<std::vector<double>> a(n, std::vector<double>(m, 0));
		for (auto &x : a) {
			for (auto &i : x) std::cin >> i;
		}
		std::vector<double> b(n);
		for (auto &x : b) std::cin >> x;
		auto x = Gauss(a, b);
		for (auto t : x) std::cout << t << "\n";
	}
	return 0;
}

有限域版

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
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
#include <bits/stdc++.h>
#define watch(x) std::cout << (#x) << " is " << (x) << std::endl
using LL = long long;
using pii = std::pair<int, int>;
using pll = std::pair<LL, LL>;

std::vector<LL> Gauss(std::vector<std::vector<LL>> A, std::vector<LL> b) {
	int n = A.size(), m = A[0].size();
	std::vector<LL> x(m);
	std::vector<int> p(m);
	std::iota(p.begin(), p.end(), 0);
	const LL M = 998244353;
	std::function<LL(LL)> inv = [&](LL a) -> LL {
		return a == 1 ? 1 : (M - M / a) * inv(M % a) % M; 
	};
	auto sub = [](LL &x, LL y) {
		(x -= y) < 0 && (x += M);
	};
	auto triangleGauss = [&](int sz) { // A[i][i] = 1
		std::vector<LL> x(sz);
		for (int i = sz - 1; i >=0; --i) {
			x[i] = (b[i] + M) % M;
			for (int row = 0; row < i; ++row) sub(b[row], A[row][i] * x[i] % M);
		}
		x.resize(A[0].size());
		return x;
	};
	auto findNonZero = [&](int i) {
		for (int row = i; row < n; ++row) if (A[row][i]) return row;
		return n;
	};
	int sz = n;
	for (int i = 0, row; i < n; ++i) {
		while (i < m) {
			row = findNonZero(i);
			if (row != n) break;
			for (int j = 0; j < n; ++j) A[j][i] = A[j][m - 1];
			std::swap(p[i], p[--m]);
		}
		if (i == m) {
			for (int row = m; row < n; ++row) if (b[row]) {
				// std::cout << "\nNo answer\n";
				return std::vector<LL>();
			}
			sz = i;
			break;
		}
		if (row != i) {
			std::swap(A[i], A[row]);
			std::swap(b[i], b[row]);
		}
		LL inva = inv(A[i][i]);
		(b[i] *= inva) %= M;
		for (int j = m - 1; j >= i; --j) (A[i][j] *= inva) %= M;
		for (int row = i + 1; row < n; ++row) {
			sub(b[row], A[row][i] * b[i] % M);
			for (int j = m - 1; j >= i; --j) {
				sub(A[row][j], A[row][i] * A[i][j] % M);
			}
		}
	}
	// if (sz != A[0].size()) std::cout << "\nInfinite answer\n";
	auto xt = triangleGauss(sz);
	for (int t = 0; t < A[0].size(); ++t) x[p[t]] = xt[t];
	return x;
}

int main() {
	// freopen("in","r",stdin);
	std::ios::sync_with_stdio(false);
	std::cin.tie(nullptr);
	int cas = 1;
	std::cin >> cas;
	while (cas--) {
		int n, m;
		std::cin >> n >> m;
		std::vector<std::vector<LL>> a(n, std::vector<LL>(m, 0));
		for (auto &x : a) {
			for (auto &i : x) std::cin >> i;
		}
		std::vector<LL> b(n);
		for (auto &x : b) std::cin >> x;
		auto x = Gauss(a, b);
		for (auto t : x) std::cout << t << "\n";
	}
	return 0;
}

上面做法是先化成上三角再求，其实也可以先直接化成对角的，两种都挺好，就不过不改了。

单纯形法

首先无论是什么样的线性规划问题，都先化成标准形式：$\max z = \sum c_i x_i$，其中 $Ax = b, x_i \geq 0$（这里输入的 $b$ 必然都是非负，否则显然无可行解），并且保证 $A$ 的左边是一个单位阵。简单的说通过 “大 M 法” 使得 ($b_1, \cdots, b_n, 0, \cdots, 0$) 是可以可行解。

  • 求极小值通过 $c$ 取负号解决
  • $\geq b$ 通过加一个变量变成等号
  • $\leq b$ 通过减一个变量变成等号
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
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
#include <bits/stdc++.h>
#define watch(x) std::cout << (#x) << " is " << (x) << std::endl
using LL = long long;
using pii = std::pair<int, int>;
using pll = std::pair<LL, LL>;

using VD = std::vector<double>;
const double eps = 1e-10;
const double inf = 1e10;
// make sure that A = (I, A') and b >= 0, compute max cx
VD simplexCore(VD c, std::vector<VD> A, VD b) {
	int n = A.size(), m = c.size();
	std::vector<int> p(m);
	std::iota(p.begin(), p.end(), 0);
	for (int i = 0; i < n; ++i) A[i].emplace_back(b[i]);
	c.emplace_back(0);
	A.emplace_back(c);
	for (int j = n; j <= m; ++j) {
		for (int i = 0; i < n; ++i) {
			A[n][j] -= A[n][i] * A[i][j];
		}
	}
	auto check = [&]() -> bool {
		for (int j = n; j < m; ++j) if (A[n][j] > eps) {
			bool flag = false;
			for (int i = 0; i < n; ++i) if (A[i][j] > eps) {
				flag = true;
				break;
			}
			if (!flag) return false;
		}
		return true;
	};
	while (1) {
		int ch = std::max_element(A[n].begin() + n, A[n].begin() + m) - A[n].begin(), hc;
		if (A[n][ch] < eps) break;
		assert(check()); // otherwise unbounded, no max solution 
		double theta = DBL_MAX;
		for (int i = 0; i < n; ++i) if (A[i][ch] > eps && A[i].back() / A[i][ch] < theta) {
			theta = A[i].back() / A[i][ch];
			hc = i;
		}
		std::swap(p[ch], p[hc]);
		double tmp = 1 / A[hc][ch];
		for (int j = n; j <= m; ++j) A[hc][j] *= tmp;
		for (int i = 0; i <= n; ++i) if (i != hc) {
			for (int j = n; j <= m; ++j) if (j != ch) {
				A[i][j] -= A[i][ch] * A[hc][j];
			}
		}
		for (int i = 0; i <= n; ++i) A[i][ch] *= -tmp;
		A[hc][ch] = tmp;
	}
	VD x(m);
	for (int i = 0; i < n; ++i) x[p[i]] = A[i].back();
	watch(-A.back().back()); // max_val
	return x; // point Corresponds to max_val
}
// compute max cx, with Aqx = bq and Alq x <= blq, end of 0 can be ommit in A and Aq
VD simplex(VD c, std::vector<VD> Aq, VD bq, std::vector<VD> Alq, VD blq) {
	assert(Aq.size() == bq.size());
	assert(Alq.size() == blq.size());
	int n = Aq.size() + Alq.size();
	int m = c.size();
	for (int i = 0; i < bq.size(); ++i) if (bq[i] < -eps) {
		for (auto &x : Aq[i]) x = -x;
		bq[i] = -bq[i];
	}
	for (int i = 0; i < blq.size(); ++i) if (blq[i] < -eps) {
		for (auto &x : Alq[i]) x = -x;
		++m;
	}
	std::vector<VD> A(n, VD(n + m));
	VD f(n + m), b(n);
	int now = n + c.size();
	for (int i = 0; i < n; ++i) A[i][i] = 1;
	for (int i = 0; i < Aq.size(); ++i) {
		for (int j = 0; j < Aq[i].size(); ++j) A[i][n + j] = Aq[i][j];
		b[i] = bq[i];
		f[i] = -inf;
	}
	for (int i = 0; i < Alq.size(); ++i) {
		for (int j = 0; j < Alq[i].size(); ++j) A[i + Aq.size()][n + j] = Alq[i][j];
		if (blq[i] < -eps) {
			A[i + Aq.size()][now++] = -1;
			f[i + Aq.size()] = -inf;
		}
		b[i + Aq.size()] = fabs(blq[i]);
	}
	for (int i = 0; i < c.size(); ++i) f[n + i] = c[i];
	auto x = simplexCore(f, A, b);
	return VD(x.begin() + n, x.begin() + n + c.size());
}
int main() {
	// freopen("in","r",stdin);
	std::ios::sync_with_stdio(false);
	std::cin.tie(nullptr);
	int cas = 1;
	std::cin >> cas;
	while (cas--) {
		int nlq, nq, m;
		std::cin >> nlq >> nq >> m;
		std::vector<std::vector<double>> Alq(nlq, std::vector<double>(m, 0));
		std::vector<std::vector<double>> Aq(nq, std::vector<double>(m, 0));
		std::vector<double> blq(nlq), bq(nq), c(m);
		for (auto &x : c) std::cin >> x;
		for (auto &x : Alq) for (auto &i : x) std::cin >> i;
		for (auto &x : blq) std::cin >> x;
		for (auto &x : Aq) for (auto &i : x) std::cin >> i;
		for (auto &x : bq) std::cin >> x;
		auto x = simplex(c, Aq, bq, Alq, blq);
		for (auto t : x) std::cout << t << "\n";
	}
	return 0;
}
/* 输入样例：（和 matlab linprog 样例一致）
1
6 1 2
1 0.33333333

1 1
1 0.25
1 -1
-0.25 -1
-1 -1
-1 1 
2 1 2 1 -1 2

1 0.25
0.5
*/

SagMath 线性规划文档 以及 Matlab linprog 文档

对偶理论

对偶理论

例题：605C，利用对偶理论之后，此题两种做法，一种是利用三分法搞定，当然了要特别注意精度问题，另一种利用半平面的交（可以查看 HDU 模板）。

以下内容以后有空再补吧

整数规划

混合型规划