Note -「线性规划」学习笔记

\(\mathcal{Definition}\)

  线性规划（Linear Programming, LP）形式上是对如下问题的描述：

\[\operatorname{maximize}~~~~z=\sum_{i=1}^nc_ix_i\\\operatorname{s.t.}\begin{cases} \sum_{j=1}^na_{ij}x_j\le b_i&i=1,2,\cdots,m\\ x_i\ge0&i=1,2,\cdots,n\end{cases} \]

其中，\(\operatorname{maximize} z\) 即最大化 \(z\) 的值，亦可作 \(\operatorname{minimum} z\)（最小化），我们称其为目标函数。\(\operatorname{s.t.}\)（subject to，服从）称为约束条件。像这样的问题形式是线性规划的标准型。对于以下问题：

\[\operatorname{maximize}~~~~z=\sum_{i=1}^nc_ix_i\\\operatorname{s.t.}\begin{cases}\sum_{j=1}^na_{ij}x_j=b_i&i=1,2,\dots,m\\x_i\ge0&i=1,2,\dots,n\end{cases} \]

则是所谓松弛型，因为我们可以通过引入松弛变量将标准型转化为松弛型：

\[\operatorname{maximize}~~~~z=\sum_{i=1}^nc_ix_i\\\operatorname{s.t.}\begin{cases} x_{i+n}=b_i-\sum_{j=1}^na_{ij}x_j&i=1,2,\cdots,m\\ x_i\ge0&i=1,2,\cdots,m+n\end{cases} \]

所以数学考试大方程题的第一问也是线性规划。此时，原有的变量 \(x_1,x_2,\cdots,x_n\) 称为非基变量，松弛用的变量 \(x_{n+1},x_{n+2},\cdots,x_{n+m}\)（\(m\) 为约束个数）称为基变量。更简便地，我们用矩阵来描述标准形式：

\[\operatorname{maximize}~~~~z=\boldsymbol{c}^T\boldsymbol x\\\operatorname{s.t.}\begin{cases} \boldsymbol{Ax}=\boldsymbol b\\ x_i\ge0\end{cases}\\ \]

其中，

\[\boldsymbol c=\begin{pmatrix}c_1\\c_2\\\vdots\\c_n\end{pmatrix}\\\boldsymbol A=\begin{pmatrix} a_{11}&a_{12}&\cdots&a_{1n}\\ a_{21}&a_{22}&\cdots&a_{2n}\\ \vdots&\vdots&\ddots&\vdots\\ a_{m1}&a_{m2}&\cdots&a_{mn}\end{pmatrix} \]

并称 \(\boldsymbol A\) 为约束矩阵，\(\boldsymbol b\) 为资源向量，\(\boldsymbol c\) 为价值向量，\(\boldsymbol x\) 为决策价值变量向量。

\(\mathcal{Algorithm}\)

  说了那么多定义，接下来我们来设计一个算法求解线性规划问题。算法的基本思路为：

  初始找出一个满足所有约束的初始解，通过这个解不断的更新，找到更优的解，直到找不到位置（类似爬山算法）。

  可以发现这样一定能找到最优解，因为如果把求解的 \(X\) 当做 \(n\) 维空间的一个点，约束条件一定在 \(n\) 维空间中框出一个 \(n\) 维凸包，一 个解为凸包的一个顶点。凸面体从哪一个发现看都是单峰的，用爬山的算法是肯定能找到最优解的。

  以

\[\operatorname{maximize}~~~~z=3x_1+x_2+2x_3+0x_4+0x_5+0x_6\\ \operatorname{s.t.}\begin{cases} x_4=30-(x_1+x_2+3x_3)\\ x_5=24-(2x_1+2x_2+5x_3)\\ x_6=36-(4x_1+x_2+2x_3)\\ x_1,x_2,\cdots,x_6\ge0 \end{cases} \]

为例，不难看出初始解可为 \(\boldsymbol x=\begin{pmatrix}0,0,0,30,24,36\end{pmatrix}\)。接着变换第三个式子：

\[x_6=36-(4x_1+x_2+2x_3)\\ \Rightarrow x_1=9-\frac{1}4x_2-\frac{1}2x_3-\frac{1}4x_6 \]

代入 \(z\)，得到：

\[z=27+\frac{3}4x_2+\frac{3}2x_3-\frac{3}4x_6 \]

  如此，交换 \(z\) 表达式中的一个非基变量和一个基变量的操作称为转轴（pivot）。形式地，设交换非基变量 \(x_u\) 与基变量 \(x_{v+n}\)：

\[x_{v+n}=b_v-\sum_{i=1}^na_{vi}x_i\\\Rightarrow~~~~x_u=\frac{b_v-x_{v+n}-\sum_{i=1}^n[i\not=u]a_{vi}x_i}{a_{vu}} \]

上例中，转轴 \(x_1-x_6\)，再将新的等式依次代入其他等式和目标函数，得到了 \(z=27+\frac{3}4x_2+\frac{3}2x_3-\frac{3}4x_6\)。如果我们通过若干次转轴使得 \(z\) 表达式中所有变量系数非正且 \(b_i\) 非负，就能得到此时的常数项为 \(z\) 的最值。总结一下，我们的算法分为两步：

1. 找初始解。一种比较好写的随机算法：随机取 \(b_i<0\)，再随机取 \(a_{ij}<0\)，转轴 \(x_j-x_{i+n}\)。不存在 \(b_i\) 时结束。虽然这样仅保证 \(b_i\ge 0\)。更为精准的算法见学长的博客。
2. 转轴使得 \(c_j\le0\)。取 \(j=\arg\max\{c_j|c_j>0\}\)，找出 \(i=\arg\min\{\frac{b_i}{a_{ij}}|a_{ij}>0\}\)，最后转轴 \(x_j-x_{i+n}\)。如此仍能保持 \(b_i\ge0\)。

  无解当且仅当找不到初始解，这显而易见；无界当且仅当第二步找不到 \(i\)，因为此时任意增大 \(c_j>0\) 的非基变量 \(x_j\)，只会使 \(z\) 和 \(\boldsymbol b\) 中某些元素增大而不可能减少，一定能够保持合法性。

\(\mathcal{Example}\)

  UOJ #179 \(97\) 分代码：

/* Clearink */

#include <ctime>
#include <cstdio>
#include <random>

#define rep( i, l, r ) for ( int i = l, rpbound##i = r; i <= rpbound##i; ++i )
#define per( i, r, l ) for ( int i = r, rpbound##i = l; i >= rpbound##i; --i )

const int MAXN = 50;
const double EPS = 1e-9;
int n, m, type, id[MAXN + 5], ref[MAXN + 5];
double a[MAXN + 5][MAXN + 5];

inline void iswp ( int& a, int& b ) { a ^= b ^= a ^= b; }
inline double dabs ( const double a ) { return a < 0 ? -a : a; }
inline int dcmp ( const double a, const double b ) {
	return dabs ( a - b ) < EPS ? 0 : ( a < b ? -1 : 1 );
}

inline bool halfp () {
	static std::mt19937 rnd ( time ( 0 ) ^ 20120712 );
	return rnd () & 1;
}

inline void pivot ( const int r, const int c ) {
	iswp ( id[r + n], id[c] );
	double tmp = -a[r][c]; a[r][c] = -1;
	rep ( i, 0, n ) a[r][i] /= tmp;
	rep ( i, 0, m ) if ( i != r && dcmp ( a[i][c], 0 ) ) {
		tmp = a[i][c], a[i][c] = 0;
		rep ( j, 0, n ) a[i][j] += tmp * a[r][j];
	}
}

inline int simplex () {
	rep ( i, 1, n ) id[i] = i;
	while ( true ) {
		int i = 0, j = 0;
		rep ( k, 1, m ) {
			if ( !~dcmp ( a[k][0], 0 ) && ( !i || halfp () ) ) {
				i = k;
			}
		}
		if ( !i ) break;
		rep ( k, 1, n ) {
			if ( dcmp ( a[i][k], 0 ) == 1 && ( !j || halfp () ) ) {
				j = k;
			}
		}
		if ( !j ) return 2;
		pivot ( i, j );
	}
	while ( true ) {
		int i = 0, j = 0;
		double mx = 0, mn = 1e9;
		rep ( k, 1, n ) if ( dcmp ( a[0][k], mx ) == 1 ) mx = a[0][j = k];
		if ( !j ) break;
		rep ( k, 1, m ) {
			if ( !~dcmp ( a[k][j], 0 ) && dcmp ( -a[k][0] / a[k][j], mn ) == -1 ) {
				mn = -a[i = k][0] / a[k][j];
			}
		}
		if ( !i ) return 1;
		pivot ( i, j );
	}
	return 0;
}

int main () {
	scanf ( "%d %d %d", &n, &m, &type );
	rep ( i, 1, n ) scanf ( "%lf", &a[0][i] );
	rep ( i, 1, m ) {
		rep ( j, 1, n ) scanf ( "%lf", &a[i][j] ), a[i][j] *= -1;
		scanf ( "%lf", &a[i][0] );
	}
	int res = simplex ();
	if ( res == 2 ) puts ( "Infeasible" );
	else if ( res ) puts ( "Unbounded" );
	else {
		printf ( "%.12f\n", a[0][0] );
		if ( type ) {
			rep ( i, 1, m ) ref[id[i + n]] = i;
			rep ( i, 1, n ) {
				printf ( "%.12f%c", ref[i] ? a[ref[i]][0] : 0, i ^ n ? ' ' : '\n' );
			}
		}
	}
	return 0;
}