宣告
本文為本人原創，轉載請註明出處。本文僅發表在部落格園，作者LightningStar。
問題描述
所有的線性規劃問題都可以歸約到標準型的問題，規約過程比較簡單且已經超出本文範圍，不再描述，可以參考擴充閱讀部分。下面直接給出線性規劃標準型描述。
標準型描述
線性規劃問題標準型的矩陣描述[1]：
目標：
約束：
我們的最終目標在約束條件下就是找到一個解\(\mathbf{x}\)，使得\(z\)最大。注意我們的描述，我們的解是找到一組值\(\mathbf{x}\)使得\(z\)最大，這組值才是問題的一個解，\(z\)得最大值究竟是多少並不是問題的解。
在後文中，粗體小寫字母一般表示向量，粗體大寫字母一般表示矩陣，小寫字母表示標量。
鬆弛型
在用單純型法求解線性規劃問題之前，必須先把線性規劃問題轉換成增廣矩陣形式。增廣矩陣形式引入非負鬆弛變數將不等式約束變成等式約束。
引入鬆弛變數\(\mathbf{\hat{x}} = \mathbf{b} - \mathbf{Ax}\)，則有：
將上述公式表示為矩陣形式則為：
我們引入的鬆弛變數\(\mathbf{\hat{x}}\)又稱基本變數，\(\mathbf{x}\)又稱非基本變數。
單純型演算法
一個例子
為便於理解和描述，我們通過一個例子來講解迭代過程：
劃為矩陣表示為：
令$\mathbf{y}=[x_1, x_2, x_3, \hat{x_1}, \hat{x_2}, \hat{x_3}]^T, \mathbf{d}=[\mathbf{c}, 0]^T, z = \mathbf{d^Ty}, \mathbf{\hat{A}=[A, I]} $，則有:
為方便求解\(z\)的最大值，我們可以設計如下的增廣矩陣[2]，通過對增廣矩陣的迭代計算可以得到\(z\)的最大值：迭代結束時增廣矩陣右上角的值的相反數。
下面開始對增廣矩陣進行迭代：
-
原線性規劃問題的一個初始解是\(\mathbf{x=0}, \mathbf{\hat{x}=b}\)，即\(\mathbf{y_0} = [0, 0, 0, 30, 24, 36]^T\)，初始\(\mathbf{d_0}=[3, 1, 2, 0, 0, 0]^T\)，\(z=\mathbf{d_0^T}\mathbf{y_0}=0\)
-
由\(\mathbf{d}\)可知，\(y_0\)的收益最大，因此選擇增大\(y_0\)以獲取更大收益。判斷依據是\(max(\mathbf{d}) == d_0\)，並且\(d_0 = 3 \gt 0\)
-
下面判斷\(y_0\)最大可以是多少。取\(\mathbf{b} ./ \hat{A}[:,0] = [30, 12, 9]\)中的最小正整數，即\(y_0 = 9\)
-
依據高斯消元法[3]，將增廣矩陣第4行作為基礎向量，將第4行作為基礎向量的依據是\(\mathbf{b} ./ \hat{A}[:,0]\)的最小值就在增廣矩陣的第4行。將增廣矩陣中其他行的\(y_0\)的係數化為0，結果為
-
下面開始新一輪的迭代過程，\(max(\mathbf{d}) == d_2\)，並且\(d_2 = 0.5 \gt 0\)，因此選擇增大\(y_2\)
-
取\(\mathbf{b} ./ \hat{A}[:,2] = [8.4, 1.5, 18]\)中的最小正整數，即\(y_2 = 1.5\)
-
取增廣矩陣的第3行作為基本向量，對增廣矩陣運用高斯消元法將\(y_2\)的其他行的係數劃為0得
-
下面開始新一輪的迭代過程，\(max(\mathbf{d}) == d_1\)，並且\(d_1 = 0.0625 \gt 0\)，因此選擇增大\(y_1\)
-
取\(\mathbf{b} ./ \hat{A}[:,1] = [-92, 4, 132]\)中的最小正整數，即\(y_1 = 4\)
-
取增廣矩陣的第3行作為基本向量，對增廣矩陣運用高斯消元法得
- \(max(\mathbf{d}) == d_1\)，並且\(d_0 = 0\)，因此迭代結束。最大值\(z=28\)(增廣矩陣右上角的值的相反數)
到目前為止，我們已經求得了標準型問題中\(z\)的最大值，但是還沒有給出一個解。我們僅僅知道如何求出\(z\)的最大值，但是什麼樣的\(\mathbf{x}\)會使得\(z\)取得最大值呢？這比知道\(z\)的最大值更重要。
現在觀察一下我們已知的一些資訊，已知\(z=28\)，已知\(y_1 = 4\)。在迭代過程中我們似乎也求得了\(y_0 = 9\)和\(y_2=1.5\)，但是實際上這是不對的。因為只有最後一次迭代的結果是準確的，而在迭代過程中得到的只是中間結果，因此我們只知道\(z=28, y_1 = 4\)。另外還有增廣矩陣。在本文開頭我們有公式：
現在我們將已知的值帶入上述公式，得到:
通過解方程（本文不涉及如何解方程）可以得到一個可行解為：
又已知原始\(\mathbf{c} = [3, 1, 2]^T\)，得\(z = \mathbf{c^T} \mathbf{x} = 28\)。
演算法過程
問題描述
為了防止讀者忘記我們要解決的問題，這裡再囉嗦一下，我們要解決的是線性規劃問題，並將所有的線性規劃問題都歸約到標準型上。因此最終問題變成了對標準型的求解。在上文中我們已經通過了一個例子來介紹如何單純形演算法的演算過程，並如何通過迭代的結果求得一個解。下面我們來將這個過程用演算法的形式表示出來，但是這個演算法僅包含迭代過程，至於如何通過迭代出來的結果求得解，則不是本文關心的內容。
演算法的輸入與輸出
這裡我們來搞清楚演算法的輸入和輸出。我們在問題中已知的是\(\mathbf{c^T}\)和矩陣\(\mathbf{A}\)，以及\(\mathbf{b}\)。因此這些已知值就是演算法的輸入。而演算法的輸出則是迭代的最後結果\(z\)的值和\((i, x_i)\)。\((i, x_i)\)是一個元組，其中\(i\)是\(\mathbf{x}\)中的第\(i\)個元素，而\(x_i\)是\(\mathbf{x}\)中的第\(i\)個元素的值（下標從0開始索引）。
演算法python實現
python的程式碼可以當作偽碼去閱讀，這裡直接給出python的實現過程。
def solve(c, A, b):
NUM_NON_BASIC_VARIABLES = c.shape[0]
NUM_BASIC_VARIABLES = b.shape[0]
# z = -d[-1]
d = np.hstack((c, np.zeros(NUM_BASIC_VARIABLES + 1)))
# 初始化增廣矩陣
_A = np.hstack((A, np.identity(NUM_BASIC_VARIABLES)))
A_hat = np.c_[_A, b]
_step = 0
last_update_x_inx = -1
last_update_x = 0
while True:
i = np.nanargmax(d[:-1])
if d[i] <= 0:
break
# 利用高斯消元法求解
_res = A_hat[:, -1] / A_hat[:, i]
# 將忽略 divided by zero 的結果，係數小於等於0的也不能考慮在內
j = np.where(_res > 0, _res, np.inf).argmin()
if _res[j] <= 0: # 係數小於等於0的會違反了 >= 0 的基本約束條件
break
last_update_x_inx = i
last_update_x = _res[j]
# 下面計算y中除了y[i]之外的值
# 1.運用高斯消元法
A_hat[j, :] = A_hat[j, :] / A_hat[j, i] # A_hat[j,i] = 1
# for _row in range(A_hat.shape[0]):
# if _row != j:
# A_hat[_row,:] = A_hat[_row,:] - A_hat[_row,i] * A_hat[j,:]
# 下面四行等價於上述的for迴圈
_tmp = np.copy(A_hat[j, :])
_A = np.outer(A_hat[:, i], _tmp) # 列向量乘以行向量
A_hat -= _A
A_hat[j, :] = _tmp
d = d - d[i] * A_hat[j, :]
# 列印中間過程
_step += 1
# print('step:', _step)
# print('d = ', d)
# print('A_hat = ', A_hat)
# print('z = ', -d[-1])
z = -d[-1]
if last_update_x_inx == -1:
return None
return (z, (last_update_x_inx, last_update_x)) # return z