“容易理解的”单纯形法
单纯形法
松弛型线性规划
线性规划标准形式：
\[\begin{aligned}
\text{maximize }& \sum_{i=1}^n c_i x_i\\
\text{s. t. }& \sum_{i=1}^n a_{j, i} x_i \le b_j, \forall j = 1, \dots, m\\
\text{and }& x_i \ge 0, \forall i = 1, \dots, n
\end{aligned}
\]
我们可以将上述问题改写，使得只有等号：
\[\begin{aligned}
\text{maximize }& \sum_{i=1}^n c_i x_i\\
\text{s. t. }& x_{j + n} = b_j - \sum_{i=1}^n a_{i, j} x_i, \forall j = 1, \dots, m\\
\text{and }& x_i \ge 0, \forall i = 1, \dots, n + m\\
\end{aligned}
\]
称 \(x_i, i \in 1, \dots, n\) 为基变量，其余的 \(m\) 个为非基变量。
转轴
定义一次转轴 \(\text{pivot}(j, i)\) 表示将基变量 \(x_i\) 用非基变量 \(x_{j + n}\) 改写，之后 \(x_i\) 将变为非基变量，\(x_{j + n}\) 变为基变量。更具体地，\(\text{pivot}(j, i)\) 分以下两个步骤：
-
改写 \(x_{j + n}\) 的等式：
\[\begin{aligned} x_{j + n} &= b_j - \sum_{i=1}^n a_{i, j} x_i\\ \implies x_i &= b_j - \sum_{\substack{k=1\\k \ne i}}^n\frac{a_{k, j}}{a_{i, j}} x_k - \frac{1}{a_{i, j}} x_{j + n} \end{aligned} \] -
对于其他的等式及 \(\sum c_ix_i\)，将其中的 \(x_i\) 用上一步的等式替换。
进行 \(\text{pivot}\) 操作之后，基变量的数量不会变化。注意，基变量满足 \(x_i = 0\)，但满足 \(x_i = 0\) 不一定是基变量。
爬山
首先，解集是凸的。并且，容易发现存在最优解满足其中 \(n\) 个 \(x_i = 0\)。因此，我们可以不断进行 \(\text{pivot}\) 操作，直到这 \(n\) 个基变量能组成一组最优解。
假设当前有一组可行解 \(\lang x_1, \dots, x_n\rang = \lang 0, \dots, 0 \rang\)，并且所有 \(b_i\) 非负。我们要尝试只对其中某个 \(x_i\) 进行增大。步骤如下：
- 找到一个 \(c_i \gt 0\)。如果不存在，说明增大任意 \(x_i\) 都不会使总和变大，即此时已取到最优解。
- 找到一个 \(a_{i, j} \gt 0\)。如果不存在，说明 \(x_i\) 可以取任意大值，即答案为 \(+\infty\)。
- 找到 \(a_{i, j} \gt 0\) 中 \(\frac{b_i}{a_{i, j}}\) 最小的那一个 \(j\)，然后执行 \(\text{pivot}(j, i)\)。在这之后，其他的 \(b'_k \gets b_k - a_{k, j} \times\frac{b_i}{a_{i, j}} \ge 0\)。
最后一步实际上让总和增加了 \(\frac{c_jb_i}{a_{i, j}}\)，并且得到了新的一组可行解 \(\lang x_1, \dots, x_{i - 1}, x_{j + n}, x_{i + 1}, \dots, x_n\rang = \lang 0, \dots, 0 \rang\)。
只要不断进行爬山操作，就一定能找到最优解。最坏操作次数为 \(O(\binom{n + m}{n})\)，但在随机数据下通常为 \(O(m)\)。
处理负权
普遍情况下 \(b_i\) 并非非负。这个时候也不能使用解 \(\lang x_1, \dots, x_n\rang = \lang 0, \dots, 0 \rang\)。我们可以引入一个变量 \(x_0\)，并且求解下列线性规划：
\[\begin{aligned}
\text{maximize }& -x_0\\
\text{s. t. }& x_{j + n} = b_j - \sum_{i=1}^n a_{i, j} x_i + x_0, \forall j = 1, \dots, m\\
\text{and }& x_i \ge 0, \forall i = 0, \dots, n + m
\end{aligned}
\]
若该线性规划的最优解为 \(0\)（即 \(x_0 = 0\)），那么我们就求得了一组可行解。否则，原线性规划无解。
我们找到最小的 \(b_j\)，执行 \(\text{pivot}(j, 0)\)，这样所有的 \(b'_k = b_k - b_i \ge 0\)。之后就可以用爬山算法求解了。
一次 \(\text{pivot}\) 的复杂度为 \(O(nm)\)，随机数据下期望复杂度为 \(O(nm^2)\)。
#include <algorithm>
#include <iostream>
#include <iomanip>
#include <numeric>
#include <limits>
using namespace std;
#define MAXN 205
#define eps 1e-8
using ld = double;
template <typename T>
constexpr T INF = numeric_limits<T>::max() / 2;
int n, m;
ld a[MAXN][MAXN], c[MAXN];
int id[MAXN];
int out[MAXN];
void pivot(int r, int b, int op) // x_b will be written as a function of x_{r + n ( + op)} (rewrite the r-th equation)
{
swap(id[r + n + op], id[b]);
ld coef = -a[r][b]; // swapping the coefficients
a[r][b] = -1;
for (int i = 0; i <= n + op; i++)
{
a[r][i] /= coef;
}
for (int i = 0; i <= m; i++)
{
if (i != r && abs(a[i][b]) > eps) // replacing x_b
{
coef = a[i][b];
a[i][b] = 0;
for (int j = 0; j <= n + op; j++)
{
a[i][j] += coef * a[r][j];
}
}
}
}
int main()
{
cin >> n >> m;
for (int i = 1; i <= n; i++) // reading c[i]
{
cin >> a[0][i];
}
for (int i = 1; i <= m; i++) // reading b[i], a[i][j]
{
for (int j = 1; j <= n; j++)
{
cin >> a[i][j];
a[i][j] = -a[i][j];
}
cin >> a[i][0]; // [0] is the constant term
}
iota(id + 1, id + n + m + 2, 1); // x_{id[i]} is on the left side of the i-th equation
do { // Step 1. find the primitive valid root
ld w = -eps;
int i = 0;
for (int k = 1; k <= m; k++)
{
if (a[k][0] < w)
{
w = a[i = k][0];
}
}
if (!i) // `forall x = 0` is a valid root
{
break;
}
copy(a[0] + 1, a[0] + n + 1, c + 1); // backup c
fill(a[0] + 1, a[0] + n + 1, 0);
for (int j = 1; j <= m; j++) // introduce the auxiliary variable x_{n + 1}
{
a[j][n + 1] = 1;
}
a[0][n + 1] = -1;
pivot(i, n + 1, 1); // after this stage, forall b[i] >= 0
for (; ; ) // (Step 2. finding the maximal c * x)
{
int j = 0;
w = eps;
for (int k = 1; k <= n + 1; k++) // finding the optimal increasing option, a.k.a., the max c[j]
{
if (a[0][k] > w)
{
w = a[0][j = k];
}
}
if (!j) // the current value is optimal
{
break;
}
w = INF<ld>;
i = 0;
for (int k = 1; k <= m; k++) // finding the min b[i] / a[i][j], when a[i][j] > 0 (a[i][j] is stored as -a[i][j])
{
if (-a[k][j] > eps && (!i || a[k][0] / -a[k][j] < w))
{
w = a[i = k][0] / -a[k][j];
}
}
if (!i) { } // -x_{n + 1} can be infinity, which is a contradiction (won't happen)
pivot(i, j, 1);
}
if (abs(a[0][0]) > eps) // valid root doesn't exist
{
cout << "Infeasible\n";
return 0;
}
for (i = 1; i <= m; i++) // pivot x_{n + 1} back to a base variable
{
if (id[i + n + 1] == n + 1)
{
for (int j = 1; j <= n + 1; j++)
{
if (abs(a[i][j]) > eps)
{
pivot(i, j, 1);
break;
}
}
break;
}
}
for (i = 1; i <= n; i++) // move x_{n + 1} back to its original position
{
if (id[i] == n + 1) // swap column i, n + 1
{
swap(id[i], id[n + 1]);
for (int j = 0; j <= m; j++)
{
swap(a[j][i], a[j][n + 1]);
}
}
}
for (i = n + 1; i <= n + m; i++) // remove the occurence of x_{n + 1}
{
id[i] = id[i + 1];
}
for (i = 1; i <= n + m; i++)
{
id[i] -= (id[i] > n);
}
for (i = 0; i <= n; i++) // restore c
{
a[0][i] = 0;
}
for (i = 0; i <= m; i++)
{
a[i][n + 1] = 0;
}
for (i = 1; i <= n; i++)
{
if (id[i] <= n) // direct contribution
{
a[0][i] += c[id[i]];
}
}
for (i = 1; i <= m; i++)
{
if (id[i + n] <= n) // indirect contribution, which can be rewritten into parts from the base variables
{
for (int j = 0; j <= n; j++)
{
a[0][j] += c[id[i + n]] * a[i][j];
}
}
}
} while (false);
for (; ; ) // Step 2. finding the maximal c * x
{
int i = 0, j = 0;
ld w = eps;
for (int k = 1; k <= n; k++)
{
if (a[0][k] > w)
{
w = a[0][j = k];
}
}
if (!j)
{
break;
}
w = INF<ld>;
for (int k = 1; k <= m; k++)
{
if (-a[k][j] > eps && (!i || a[k][0] / -a[k][j] < w))
{
w = a[i = k][0] / -a[k][j];
}
}
if (!i) // the sought answer can be infinity
{
cout << "Unbounded\n";
return 0;
}
pivot(i, j, 0);
}
// Step 3. print the answer
cout << fixed << setprecision(10) << a[0][0] << '\n';
for (int i = n + 1; i <= n + m; i++) // only non-base variables are non-zero
{
if (id[i] <= n)
{
out[id[i]] = i - n;
}
}
for (int i = 1; i <= n; i++)
{
if (!out[i]) // base variable
{
cout << "0 ";
}
else
{
cout << a[out[i]][0] << ' ';
}
}
cout << '\n';
return 0;
}
浙公网安备 33010602011771号