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