[ simplex algorithm ] 線性規劃 - 單純形法

線性規劃簡介

線性規劃是在線性約束條件下尋找某個目標線性函數的最大、最小值，例如：
$$\begin{array}{lrl}     \min& -3x_1+2x_2 & \\     \mathrm{subject~to}         &4x_1-6x_2&\le 5\\         &5x_1-2x_2&\ge 7\\         &x_1&\ge 0.     \end{array}$$

在最佳化問題中，線性規劃是非常重要的領域，很多的演算法問題都可以轉換成線性規劃問題來解決，例如：最大網路流、最大匹配等。因此研究如何系統化的解決線性規劃問題是非常重要的事情

線性規劃標準形

所有的線性規劃問題都可以轉換成下列標準型：

$$\begin{array}{lrl}     \max &\sum_{j=1}^n c_j \times x_j &\\     \mathrm{subject~to}         &\sum_{j=1}^n A_{i,j} \times x_j &\le b_j ,\; \forall i=1\sim m\\         &x_j& \geq 0, \; \forall  j = 1\sim n     \end{array}$$

一般來說大多數解線性規劃的演算法都會要求將問題轉換成標準形，因此務必要理解其轉換方法。轉換成標準形的步驟如下：

1. 將最小化改成最大化：
   只要改變目標函數的正負號，即可將最小化問題改成標準型設定的最大化問題，也就是說 $$\min \sum_{j=1}^n c_j \times x_j$$ 等價於 $$\max \sum_{j=1}^n -c_j \times x_j$$
2. 去除等式：
   如果約束條件中存在等式，將其轉化為兩個分別為大於等於及小於等於的不等式取代，例如：
   $$x_1+x_2=5$$ 等價於 $$x_1+x_2 \le 5\\x_1+x_2 \ge 5$$
3. 去除大於等於：
   如果約束條件中存在大於等於約束，將約束兩邊取負即可，例如：
   $$x_1+x_2 \ge 5$$ 等價於 $$-x_1-x_2 \le -5$$
4. 去除自由變數：
   如果未知數 $x_i$ 沒有非負約束，我們說 $x_i$ 是自由變數。加入新變量 $x_i'$，並用 $x_i-x_i'$替換原來的變量 $x_i$，並增加$x_i,x_i' \ge 0$的條件即可

標準化範例：

$$\begin{array}{lrl}     \min& -x_1+2x_2 & \\     \mathrm{subject~to}         &x_1+x_2&= 5\\         &x_1-x_2&\le 2\\         &x_2&\ge 0.     \end{array}$$
首先將最小化改成最大化：
$$\begin{array}{lrl}     \max& x_1-2x_2 & \\     \mathrm{subject~to}         &x_1+x_2&= 5\\         &x_1-x_2&\le 2\\         &x_2&\ge 0.     \end{array}$$
接著去除等式：
$$\begin{array}{lrl}     \max& x_1-2x_2 & \\     \mathrm{subject~to}         &x_1+x_2&\le 5\\         &x_1+x_2&\ge 5\\         &x_1-x_2&\le 2\\         &x_2&\ge 0.     \end{array}$$
去除大於等於：
$$\begin{array}{lrl}     \max& x_1-2x_2 & \\     \mathrm{subject~to}         &x_1+x_2&\le 5\\         &-x_1-x_2&\le -5\\         &x_1-x_2&\le 2\\         &x_2&\ge 0.     \end{array}$$
最後去除自由變數，將$x_1$用$x_1-x_3$取代：
$$\begin{array}{lrl}     \max& x_1-2x_2-x_3 & \\     \mathrm{subject~to}         &x_1+x_2-x_3&\le 5\\         &-x_1-x_2+x_3&\le -5\\         &x_1-x_2-x_3&\le 2\\         &x_1,x_2,x_3&\ge 0.     \end{array}$$

單純形法

輸入input

本人實作的單純形法其輸入是一個被稱為單純形表的 $m+1 \times n+1$ 矩陣$A$，標準形中每個項目的值都可以對應到矩陣$A$上：
$$\begin{array}{lrl}     \max &\sum_{j=1}^n A_{0,j} \times x_j &\\     \mathrm{subject~to}         &\sum_{j=1}^n A_{i,j} \times x_j &\le A_{i,0} ,\; \forall i=1\sim m\\         &x_j& \geq 0, \; \forall  j = 1\sim n     \end{array}$$
可以發現：

1. $A_{0,1} \sim A_{0,n}$對應標準形中的$c_1 \sim c_n$
2. $A_{i,0}$對應標準形中的$b_i ,\; \forall i=1\sim m$

為了方便計算請設$A_{0,0}=0$

上面標準化範例寫成單純形表為以下形式：

$$\begin{bmatrix} 0& 1 &-2 &-1 \\ 5 &1& 1 &-1 \\ -5 &-1 &-1 &1 \\ 2 &1 &-1& -1 \end{bmatrix}$$

輸出output

輸出為一個$size=n+1$的vector $ans$，其中：

1. $ans[0]$為線性規劃的解，也就是最大值
2. $ans[1] \sim ans[n]$為變數$x_1 \sim x_n$的其中一組滿足條件的答案
3. 如果無解或是解無限大則$ans$的$size=0$
4. 注意如果答案是 $0$ 由於浮點數誤差的關係可能會產生 $``-0"$ 這種東西，要小心處理

上面標準化範例的解如下：

$$\{0.5,\; 3.5,\; 1.5,\; 0\}$$

意即在 $x_1=3.5,\;x_2=1.5,\;x_3=0$ 的情況下會有最大值 $x_1-2x_2-x_3=0.5$

程式碼

使用C++11

	#include<vector>
	#include<algorithm>
	using namespace std;
	template<class VDB>
	VDB simplex(int m,int n,vector<VDB> a){
	vector<int> left(m+1), up(n+1);
	iota(left.begin(), left.end(), n);
	iota(up.begin(), up.end(), 0);
	auto pivot = [&](int x, int y){
	swap(left[x], up[y]);
	auto k = a[x][y]; a[x][y] = 1;
	vector<int> pos;
	for(int j = 0; j <= n; ++j){
	a[x][j] /= k;
	if(a[x][j] != 0) pos.push_back(j);
	}
	for(int i = 0; i <= m; ++i){
	if(a[i][y]==0 || i == x) continue;
	k = a[i][y], a[i][y] = 0;
	for(int j : pos) a[i][j] -= k*a[x][j];
	}
	};
	for(int x,y;;){
	for(int i=x=1; i <= m; ++i)
	if(a[i][0]<a[x][0]) x = i;
	if(a[x][0]>=0) break;
	for(int j=y=1; j <= n; ++j)
	if(a[x][j]<a[x][y]) y = j;
	if(a[x][y]>=0) return VDB();
	pivot(x, y); // infeasible
	}
	for(int x,y;;){
	for(int j=y=1; j <= n; ++j)
	if(a[0][j] > a[0][y]) y = j;
	if(a[0][y]<=0) break;
	x = -1;
	for(int i=1; i<=m; ++i) if(a[i][y] > 0)
	if(x == -1 || a[i][0]/a[i][y] < a[x][0]/a[x][y])
	x = i;
	if(x == -1) return VDB();
	pivot(x, y); // unbounded
	}
	VDB ans(n + 1);
	for(int i = 1; i <= m; ++i)
	if(left[i] <= n) ans[left[i]] = a[i][0];
	ans[0] = -a[0][0];
	return ans;
	}

view raw simplex.h hosted with ❤ by GitHub

	#include<cstdio>
	#include"simplex.cpp"
	int main(){
	int m = 3, n = 3;
	vector<vector<double>> a=
	{
	{ 0, 1, -2, -1},
	{ 5, 1, 1, -1},
	{-5, -1, -1, 1},
	{ 2, 1, -1, -1}
	};
	auto res=simplex(m,n,a);
	if(res.size()==n+1){
	for(int i=0;i<=n;++i) printf("%f ",res[i]);puts("");
	}
	return 0;
	}

view raw test.cpp hosted with ❤ by GitHub

原理

有空時會補上，先讓我好好休息吧！終於可以正常睡覺了