1. 模型

常見的線性規劃模型如下：
max z=cx$z=cx$
s.t. Ax=b$Ax=b$

2. 求解步驟

假設B是基變數集合，通過矩陣的線性變換，基變數可由非基變量表示：
x′i=ci+Σj∉Bmi,jx′j,(i∈B)$x_i^{\prime }=c_i+{\mathrm{\Sigma }}_{j\notin B}{m}_{i,j}{x}_j^{\prime },(i\in B)$
目標函式z也可以完全由非基變量表示：
z=z0+Σj∉Bcjx′j$z=z_0+{\mathrm{\Sigma }}_{j\notin B}{c}_j{x}_j^{\prime }$
當達到最優解時，目標函式中所有的係數c≤0，這樣非基變數等於0時，目標函式可以取到最大值。以此為目標，每次將最大的正係數max{cj$c_j$}對應的非基變數替換為基變數，同時將min{bj/ai,j$b_j/a_{i,j}$}對應的基變數替換為非基變數

3. python演算法實現

這裡假設原問題都是小於等於約束，這樣新增鬆弛變數之後，問題一定有初始可行解；同時假設問題存在有限最優解。特殊情況將在下一節進行處理。程式碼為：

import numpy as np

def pivot():
    l = list(d[0][:-2])
    jnum = l.index(max(l)) #轉入編號
    m = []
    for i in range(bn):
        if d[i][jnum] == 0:
            m.append(0.)
        else
 
:
            m.append(d[i][-1]/d[i][jnum])
    inum = m.index(min([x for x in m[1:] if x!=0]))  #轉出下標
    s[inum-1] = jnum
    r = d[inum][jnum]
    d[inum] /= r
    for i in [x for x in range(bn) if x !=inum]:
        r = d[i][jnum]
        d[i] -= r * d[inum]

def solve():
    flag = True
    while flag:
        if
 
 max(list(d[0][:-1])) <= 0: #直至所有係數小於等於0
            flag = False
        else:
            pivot()

def printSol():
    for i in range(cn - 1):
        if i in s:
            print("x"+str(i)+"=%.2f" % d[s.index(i)+1][-1])
        else:
            print("x"+str(i)+"=0.00")
    print("objective is %.2f"%(-d[0][-1]))

呼叫的例子：

d = np.loadtxt("data.txt", dtype=np.float)
(bn,cn) = d.shape
s = list(range(cn-bn,cn-1)) #基變數列表
solve()
printSol()

data.txt檔案中的內容為：

1 14 6 0 0 0 0 0

1 1 1 1 0 0 0 4

1 0 0 0 1 0 0 2

0 0 1 0 0 1 0 3

0 3 1 0 0 0 1 6

代表的求解模型是：

max z=x0+14∗x1+6∗x2$z=x_0+14\ast x_1+6\ast x_2$

s.t.

x0+x1+x2≤4$x_0+x_1+x_2\le 4$

x0≤2$x_0\le 2$

x2≤3$x_2\le 3$

3∗x1+x2≤6$3\ast x_1+x_2\le 6$

執行後輸出結果為：

x0=0.00

x1=1.00

x2=3.00

x3=0.00

x4=2.00

x5=0.00

x6=0.00

objective is 32.00

4. 寫後感

將simplex用程式碼寫出來，才覺得以前糾結那麼久的問題原來那麼簡單。兩三行程式碼能說清楚的事，何必寫一堆看得人眼花繚亂的數學公式呢。
另外，線性規劃還有一些很基礎的理論要掌握好：
1. 極點和極方向的理論，這個是單純型法的理論基礎。可以參考這裡
2. 對偶理論，這個在以後經常會用到。