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;
}