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

(mathcal{Definition})

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

[operatorname{maximize}~~~~z=sum_{i=1}^nc_ix_i\operatorname{s.t.}egin{cases} sum_{j=1}^na_{ij}x_jle b_i&i=1,2,cdots,m\ x_ige0&i=1,2,cdots,nend{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.}egin{cases}sum_{j=1}^na_{ij}x_j=b_i&i=1,2,dots,m\x_ige0&i=1,2,dots,nend{cases} ]

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

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

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

[operatorname{maximize}~~~~z=oldsymbol{c}^Toldsymbol x\operatorname{s.t.}egin{cases} oldsymbol{Ax}=oldsymbol b\ x_ige0end{cases}\ ]

其中，

[oldsymbol c=egin{pmatrix}c_1\c_2\vdots\c_nend{pmatrix}\oldsymbol A=egin{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} ]

并称 (oldsymbol A) 为约束矩阵，(oldsymbol b) 为资源向量，(oldsymbol c) 为价值向量，(oldsymbol 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.}egin{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_6ge0 end{cases} ]

为例，不难看出初始解可为 (oldsymbol x=egin{pmatrix}0,0,0,30,24,36end{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 ot=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_ige 0)。更为精准的算法见学长的博客。
2. 转轴使得 (c_jle0)。取 (j=argmax{c_j|c_j>0})，找出 (i=argmin{frac{b_i}{a_{ij}}|a_{ij}>0})，最后转轴 (x_j-x_{i+n})。如此仍能保持 (b_ige0)。

  无解当且仅当找不到初始解，这显而易见；无界当且仅当第二步找不到 (i)，因为此时任意增大 (c_j>0) 的非基变量 (x_j)，只会使 (z) 和 (oldsymbol 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
", 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 ? ' ' : '
' );
			}
		}
	}
	return 0;
}