單純形法（Simplex Method）是解決線性規劃問題的一種高效且廣泛使用的演算法。由喬治·丹齊克（George Dantzig）在20世紀40年代提出，這一方法透過系統地檢查可行解空間的極點，從而找到最優解。由於其計算效率高，單純形法迅速成為線性規劃問題中最重要和最常用的演算法之一。它的應用範圍廣泛，能夠有效解決實際中的大規模最佳化問題，因此在現代工業和經濟管理中扮演著關鍵角色。
| 初始單純形表 | 迭代後單純形表 |
|---|---|
 |
 |
一、單純形表的消元法建構
線性規劃的標準型為
其中： \(\boldsymbol{A}\) 是 \(m \times n\) 矩陣， \(m \leqslant n\) 且秩 \(r(\boldsymbol{A})=m\) ，即 \(\boldsymbol{A}\) 中至少有一個 \(m \times m\) 滿秩子矩陣。
1.1 基本可行解的代數消元(最優解判斷條件)
不妨設 \(\boldsymbol{B}=\left(\begin{array}{llll}\boldsymbol{P}_1 & \boldsymbol{P}_2 & \cdots & \boldsymbol{P}_m\end{array}\right)\) 是線性規劃的一個基, 將有關矩陣和向量分塊，記 \(\boldsymbol{A}=(\boldsymbol{B}, N), C=\left(\boldsymbol{C}_B, \boldsymbol{C}_N\right), X=\left(\boldsymbol{X}_B, \boldsymbol{X}_N\right)^{\mathrm{T}}\) ，為求線性規劃的基本最優解，先要求線性規劃的基本可行解，這樣就需將約束中的基變數用非基變數表示出來。
用 \(\boldsymbol{B}^{-1}\) 左乘約束方程 \(\boldsymbol{A} \boldsymbol{X}=\boldsymbol{b}\) 的兩端，得
即
其中 \(\boldsymbol{E}\) 是單位矩陣。整理, 得
將其代入目標函式中，得
即
非基變數 \(x_j\) 前面的係數 \(c_j-\boldsymbol{C}_B \boldsymbol{B}^{-1} \boldsymbol{P}_j\) 稱為變數 \(x_j\) 的檢驗數，表示該變數增加或減少一個單位所引起目標函式值的變化，它可以用來判斷將 \(x_j\) 變成基變數後能否改進目標函式值，以後記為 \(\sigma_j=c_j-\boldsymbol{C}_B \boldsymbol{B}^{-1} \boldsymbol{P}_j\) 。
最優解判定定理：對某基本可行解 \(\boldsymbol{X}_B=\boldsymbol{B}^{-1} \boldsymbol{b}\) 其他 \(\boldsymbol{x}_N=0\), 若所有的 \(\sigma_j=\boldsymbol{c}_{\boldsymbol{j}}\) \(-\boldsymbol{C}_B \boldsymbol{B}^{-1} \boldsymbol{P}_j \leqslant 0\) ，則該解為最優解。
1.2 基本可行解的矩陣消元
考慮線性方程組 \(\left\{\begin{array}{l}\boldsymbol{A} \boldsymbol{X}=\boldsymbol{b} \\ \boldsymbol{z}=\boldsymbol{C} \boldsymbol{X}\end{array}\right.\), 其變數為 \(\left[\begin{array}{l}\boldsymbol{X} \\ \boldsymbol{z}\end{array}\right]\), 為便於求解, 整理得方程組 \(\left\{\begin{array}{c}0 \cdot z+\boldsymbol{A} \boldsymbol{X}=\boldsymbol{b} \\ -z+\boldsymbol{C X}=\boldsymbol{0}\end{array}\right.\), 其增廣矩陣見下表, 應用高斯消元法, 求解線性方程組的解。
| 常數項 | z | XB | XN | |
|---|---|---|---|---|
| b | 0 | B | N | （1-1） |
| 0 | -1 | \(C_B\) | \(C_N\) | （1-2） |
用 \(\boldsymbol{B}^{-1}\) 左乘方程組(1-1)的兩端, 將 \(\boldsymbol{X}_B\) 的係數化為單位矩陣, 得下表。
| 常數項 | z | $$X_B$$ | $$X_N$$ | |
|---|---|---|---|---|
| \(B^{-1} b\) | 0 | \(B^{-1}B=E\) | \(B^{-1}N\) | （1-3） |
| 0 | -1 | \(C_B\) | \(C_N\) |
將方程組（1-3）左乘 \(-\boldsymbol{C}_B\) 加到方程（1-2）兩邊，消元化簡得下表。
| 常數項 | z | \(X_B\) | \(X_N\) | |
|---|---|---|---|---|
| \(B^{-1} b\) | 0 | $$B^{-1}B=E$$ | $$B^{-1}N$$ | |
| $$-C_BB^{-1}b$$ | -1 | $$C_B - C_BB^{-1}B = 0$$ | $$C_N - C_BB^{-1}N $$ | （1-4） |
注意式(1-4)中基向量 \(\boldsymbol{X}_B\) 、非基向量 \(\boldsymbol{X}_N\) 的係數，它們形式相似，當前者 \(\boldsymbol{C}_B-\boldsymbol{C}_B \boldsymbol{B}^{-1} \boldsymbol{B}=0\)時，後者即為檢驗數 \(\boldsymbol{C}_{\boldsymbol{N}}-\boldsymbol{C}_B \boldsymbol{B}^{-1} \boldsymbol{N}\) 的矩陣形式。也就是說，用消元法將目標函式中基變數的係數化為零的同時，就會得到非基變數的檢驗數。事實上， \(\boldsymbol{C}_B-\boldsymbol{C}_B \boldsymbol{B}^{-1} \boldsymbol{B}\) 也可看成基向量的檢驗數, 應用消元法後, 當基向量的檢驗數化為零時, 非基變數的係數 \(\boldsymbol{C}_{\boldsymbol{N}}-\boldsymbol{C}_B \boldsymbol{B}^{-1} \boldsymbol{N}\) 就是其檢驗數。
從上述論述可知, 計算檢驗數的方法有兩種：一是透過公式 \(\sigma_j=c_j-\boldsymbol{C}_B \boldsymbol{B}^{-1} \boldsymbol{P}_j\) 計算得出；二是應用消元法，將目標函式中基向量的係數都化為零時，非基向量的係數就是檢驗數。這也是單純表結構的出處。
1.3 單純形法=逆矩陣法
單純形法的核心思想是透過變換基本變數和非基本變數，找到一個新的基本可行解，並透過比較檢驗數判斷該解是否為最優解。在單純形法中，每一個基本可行解都對應著一個基矩陣的逆矩陣。基矩陣是從約束條件係數矩陣中選取的列所構成的矩陣，它的逆矩陣在每次迭代中都要更新。
這個逆矩陣的求解是單純形法消元過程的核心，因為它直接影響到新解的生成和目標函式的最佳化。具體來說，當我們選擇一個進入基的變數並排除一個離開基的變數時，這相當於對基矩陣進行了一次列替換操作。為了保持計算的有效性，我們需要快速更新逆矩陣。這通常透過一個修正的高斯消元法來完成，使得新基的逆矩陣可以透過舊基的逆矩陣迅速計算出來。
單純形的每次迭代，都可以使用逆矩陣來計算檢驗數，從而判斷是否可以進一步最佳化目標函式。檢驗數的計算過程依賴於逆矩陣，因為檢驗數用於評估當前解的可行性及其對目標函式的影響。如果所有檢驗數都滿足最優性條件，那麼當前解就是最優解；否則，我們需要透過逆矩陣的進一步操作繼續迭代。
例1：已知初始單純形表和最終單純形表, 試求解以下問題。
（1）在初始單純形表中找出最優基 \(\boldsymbol{B}\) ，在最終單純形表裡找出 \(\boldsymbol{B}^{-1}\) 。
（2）完成最終單純形表。(3) 給出最優解與最優值。
| $$C_B$$ | $$X_B$$ | $$B^{-1}b$$ | $$x_1$$ | $$x_2$$ | $$x_3$$ | $$x_4$$ | $$x_5$$ | $$x_6$$ | |
|---|---|---|---|---|---|---|---|---|---|
| 初始表 | 0 | \(x_4\) | 60 | 3 | 1 | 1 | 1 | 0 | 0 |
| 0 | \(x_5\) | 10 | 1 | -1 | 2 | 0 | 1 | 0 | |
| 0 | \(x_6\) | 20 | 1 | 1 | -1 | 0 | 0 | 1 | |
| \(σ_j\) | 2 | -1 | 1 | 0 | 0 | 0 | |||
| 最終表 | \(x_4\) | -1 | -2 | ||||||
| \(x_1\) | 1/2 | 1/2 | |||||||
| \(x_2\) | -1/2 | 1/2 | |||||||
| \(σ_j\) |
確定最優基和\(B^{-1}\)
- 最優基 \(\boldsymbol{B}\)：從初始基變數\(x_4\)、\(x_5\)、\(x_6\) 透過單純形法操作，最終基變數變為\(x_4\)、\(x_1\)、\(x_2\)
- \(B^{-1}\)：透過最終單純形表的\(x_4\)、\(x_5\)、\(x_6\)列得到 \(B^{-1}\)：
計算\(B^{-1}\)
- 右端項\(b = \begin{bmatrix} 60 \quad 10 \quad 20 \end{bmatrix}\)
計算最終單純形表中的列向量
-\(x_3\)列向量計算：
- 初始表中\(x_3\)列為\([1, 2, -1]^T\)
因此\(x_3\)列為\(\begin{bmatrix} 1 \quad 0.5 \quad -1.5 \end{bmatrix}^T\)。
計算檢驗數\(σ_j\)
目標函式為\(z = 2x_1 -x_2+x_3\)
- 基變數的成本系數向量\(C_B = [0, 2, -1]\)
| 變數 | $$P_j$$ | $$σ_j = c_j-C_B B^{-1}P_j - $$ | 計算過程 |
|---|---|---|---|
| $$x_1 $$ | $$ \begin{bmatrix} 0 \quad 1 \quad 0 \end{bmatrix} $$ | 0 | $$ 2-[0, 2, -1] \cdot \begin{bmatrix} 0 \ 1 \ 0 \end{bmatrix} = 0 $$ |
| $$x_2$$ | $$ \begin{bmatrix} 0 \quad 0 \quad 1 \end{bmatrix} $$ | 0 | $$-1-[0, 2, -1] \cdot \begin{bmatrix} 0 \ 0 \ 1 \end{bmatrix} = 0 $$ |
| $$ x_3 $$ | $$ \begin{bmatrix} 1 \quad 0.5 \quad -1.5 \end{bmatrix} $$ | -1.5 | $$1-[0, 2, -1] \cdot \begin{bmatrix} 1 \ 0.5 \ -1.5 \end{bmatrix} = -1.5 $$ |
| $$x_4 $$ | $$ \begin{bmatrix} 1 \quad 0 \quad 0 \end{bmatrix} $$ | 0 | $$ 0-[0, 2, -1] \cdot \begin{bmatrix} 1 \ 0 \ 0 \end{bmatrix} = 0 $$ |
| $$ x_5$$ | $$ \begin{bmatrix} -1 \quad 0.5 \quad -0.5 \end{bmatrix} $$ | -1.5 | $$0-[0, 2, -1] \cdot \begin{bmatrix} -1 \ 0.5 \ -0.5 \end{bmatrix} = -1.5 $$ |
| $$ x_6 $$ | $$ \begin{bmatrix} -2 \quad 0.5 \quad 0.5 \end{bmatrix} $$ | -0.5 | $$0- [0, 2, -1] \cdot \begin{bmatrix} -2 \ 0.5 \ 0.5 \end{bmatrix} = -0.5 $$ |
最終單純形表
| $$C_B$$ | $$X_B$$ | $$B^{-1}b$$ | $$x_1$$ | $$x_2$$ | $$x_3$$ | $$x_4$$ | $$x_5$$ | $$x_6$$ | |
|---|---|---|---|---|---|---|---|---|---|
| 最終表 | 0 | $$x_4$$ | 10 | 0 | 0 | 1 | 1 | -1 | -2 |
| 3 | $$x_1$$ | 15 | 1 | 0 | 0.5 | 0 | 1/2 | 1/2 | |
| 2 | $$x_2$$ | 5 | 0 | 1 | -1.5 | 0 | -1/2 | 1/2 | |
| $$σ_j$$ | 0 | 0 | -1.5 | 0 | -1.5 | -0.5 |
最優解與最優值
- 最優解：\(X^*=\left[\begin{array}{llllll}15 & 5 & 0 & 10 & 0 & 0\end{array}\right]^T\)。
- 最優值：\(z^*=25\)
二、單純形的迭代步驟
例1：對於線性規劃問題
其中， \(x_3\) 和 \(x_4\) 是松馳變數。
決策變數：\(\mathbf{x}=\left[\begin{array}{l}x_1& x_2 & x_3 & x_4\end{array}\right]^T\) 是決策變數向量。
目標函式: Max \(\mathbf{c}^{\top} \mathbf{x}\) 其中， \(\mathbf{c}=\left[\begin{array}{c}2\quad 3\quad 0 \quad0\end{array}\right]^T\) 為目標函式的係數向量。
約束條件: \(\mathbf{A} \mathbf{x}=\mathbf{b}\) ，其中 \(\mathbf{A}=\left[\begin{array}{cccc}4 & 1 & 1 & 0 \\ -1 & 1 & 0 & 1\end{array}\right]\) 是約束係數矩陣。
\(\mathbf{b}=\left[\begin{array}{c}16 & 6\end{array}\right]^T\) 是約束右端項向量。因此，標準型的線性規劃問題可以表示為：
初始時令 \(x_N=[0,0]\)
**第一輪迭代: **此時 \(c_N-c_B B^{-1} N=[2,3]\) 因此選擇 \(x_2\) 作為入基變數更為高效，且 \(b-N x_N=\left[\begin{array}{c}16-4 x_1-x_2 \\ 6+x_1-x_2\end{array}\right]>=0\) 故 \(x_2=6\) 並令 \(x_4=0, x_2\) 入基, \(x_4\) 出基。經過此輪迭代後，各變數如下
**第二輪迭代: **此時 \(c_N-c_B B^{-1} N=[5,-3]\) ，選擇 \(x_1\) 作為入基變數更為高效，且
\(b-N x_N=\left[\begin{array}{c}16-4 x_1 \\ 6+x_1\end{array}\right]>=0\) 故 \(x_1=4\) 並令 \(x_3=0, x_1\) 入基， \(x_3\) 出基。經過此輪迭代後，各變數如下
**第三輪迭代: **此時 \(c_N-c_B B^{-1} N=[-1,-2]\) 都小於 0 ，達到收斂條件，此時令 \(x_N=[0,0]\) 解得 \(x_B=[2,8]\) 最小值為28。
**例2：*求解下面線性模型
初始單純形表
| \(x_1\) | \(x_2\) | \(x_3\) | \(x_4\) | \(x_5\) | \(b\) | \(\theta\) | |
|---|---|---|---|---|---|---|---|
| 目標函式 | 70 | 30 | 0 | 0 | 0 | ||
| 約束 1 | 3 | 9 | 1 | 0 | 0 | 540 | |
| 約束 2 | 5 | 5 | 0 | 1 | 0 | 450 | |
| 約束 3 | 9 | 3 | 0 | 0 | 1 | 720 |
判斷當前頂點是否是最優解
對於最大化問題,若當前目標函式中的非基變數的係數小於等於0時，則所得的解為最優解。而在當前的例子中，非基變數的係數分別為70和30，意味著在可行域內隨著非基變數\(x_1\)和\(x_2\)的增大，目標函式就會繼續增大，因此當前的解不是最優解。故判斷當前所得的解是否為最優解時，只需判斷目標函式中非基變數的係數是否小於等於0。
進基和出基變數
變數的出基與入基，在幾何影像上表現為頂點的變化。入基的規則為選擇使目標函式z變化最快的非基變數入基，即選擇係數最大且為正數的非基變數入基，故在本例中選擇\(x_1\)入基。出基的規則則需要引入一個新的量\(\theta\)，\(\theta=b/a_i\)（\(a_i\)為非基變數係數，\(a_i\)的選擇的是剛剛入基的非基變數的係數），選擇最小的\(\theta\)出基。
| \(x_1\) | \(x_2\) | \(x_3\) | \(x_4\) | \(x_5\) | \(b\) | \(\theta\) | |
|---|---|---|---|---|---|---|---|
| 目標函式 | 70 | 30 | 0 | 0 | 0 | ||
| 約束 1 | 3 | 9 | 1 | 0 | 0 | 540 | 180 |
| 約束 2 | 5 | 5 | 0 | 1 | 0 | 450 | 90 |
| 約束 3 | 9 | 3 | 0 | 0 | 1 | 720 | 80 |
第一次迭代，經過\(x_1\)入基與\(x_5\)出基的運算後，得到的結果如下所示。
| \(x_1\) | \(x_2\) | \(x_3\) | \(x_4\) | \(x_5\) | \(b\) | \(\theta\) | |
|---|---|---|---|---|---|---|---|
| \(Z\) | 0 | 20/3 | 0 | 0 | -70/9 | 5600 | |
| \(x_3\) | 0 | 8 | 1 | 0 | -1/3 | 300 | |
| \(x_4\) | 0 | 10/3 | 0 | 1 | -5/9 | 50 | |
| \(x_1\) | 1 | 1/3 | 0 | 0 | 1/9 | 80 |
令非基變數的值為0，則得到第一次迭代的解，如下所示
判斷該解是否為最優解(重複第二個步驟，直到是最優解為止)，顯然由於非基變數\(x_2\)的係數大於0，故當前位置還不是最優解，再次重複上面步驟。
三、練習
| Food | Energy（能量） | Protein（蛋白質） | Calcium（鈣） | Price |
|---|---|---|---|---|
| Oatmeal（燕麥） | 110 | 4 | 2 | 3 |
| Whole milk（全奶） | 160 | 8 | 285 | 9 |
| Cherry pie（草莓派） | 420 | 4 | 22 | 20 |
| Pork with beans（豬肉） | 260 | 14 | 80 | 19 |
決策變數：
- \(x_1\)：購買燕麥（Oatmeal）的數量
- \(x_2\)：購買全奶（Whole milk）的數量
- \(x_3\)：購買草莓派（Cherry pie）的數量
- \(x_4\)：購買豬肉（Pork with beans）的數量
目標函式：
最小化總價格：
約束條件：
- 能量需求：至少2000單位能量
- 蛋白質需求：至少55單位蛋白質
- 鈣需求：至少800單位鈣
- 非負性約束
from pulp import LpProblem, LpMinimize, LpVariable, lpSum, value, LpStatus
# 定義問題
problem = LpProblem("Minimize Cost", LpMinimize)
# 定義變數
x1 = LpVariable('x1', lowBound=0, cat='Continuous')
x2 = LpVariable('x2', lowBound=0, cat='Continuous')
x3 = LpVariable('x3', lowBound=0, cat='Continuous')
x4 = LpVariable('x4', lowBound=0, cat='Continuous')
# 目標函式
problem += 3*x1 + 9*x2 + 20*x3 + 19*x4, "Total Cost"
# 約束條件
problem += 110*x1 + 160*x2 + 420*x3 + 260*x4 >= 2000, "Energy"
problem += 4*x1 + 8*x2 + 4*x3 + 14*x4 >= 55, "Protein"
problem += 2*x1 + 285*x2 + 22*x3 + 80*x4 >= 800, "Calcium"
# 求解
problem.solve()
# 輸出結果
print(f"Status: {LpStatus[problem.status]}")
print(f"Oatmeal (燕麥) units: {value(x1)}")
print(f"Whole milk (全奶) units: {value(x2)}")
print(f"Cherry pie (草莓派) units: {value(x3)}")
print(f"Pork with beans (豬肉) units: {value(x4)}")
print(f"Total minimum cost: {value(problem.objective)}")
Status: Optimal
Oatmeal (燕麥) units: 14.24428
Whole milk (全奶) units: 2.7070577
Cherry pie (草莓派) units: 0.0
Pork with beans (豬肉) units: 0.0
Total minimum cost: 67.09635929999999
總結
單純形法是線性規劃問題中的一種經典演算法，它透過逐步最佳化，找到能夠最大化或最小化目標函式的可行解。儘管單純形法在最壞情況下可能需要指數級時間，但在實際應用中，單純形法通常能高效地找到最優解，因此它被廣泛應用於各種線性規劃問題中，例如生產計劃、資源分配、運輸最佳化等。
單純形法的直觀性在於從一個頂點開始沿邊界移動到另一個頂點，不斷提高目標函式值，直到達到最優解。這種方法簡單且易於理解，對大多數線性規劃問題都能有效求解。然而，單純形法也存在一些侷限性，例如在處理退化問題時可能會出現迴圈，導致演算法陷入無限迴圈的困境。為了克服這些侷限性，研究者們提出了多種改進方案，如反週期規則來防止迴圈、對偶單純形法來處理不可行的初始解等。此外，內點法作為一種替代演算法，透過從可行域的內部逼近最優解，提供了不同於單純形法的求解思路，並且在某些情況下表現出更好的最壞情況效能。透過這些改進和新方法的引入，線性規劃問題的求解在理論和實踐上都得到了極大的豐富和發展，使得我們能夠解決更復雜、更大規模的問題。這些進步不僅提升了演算法的效率，也擴充了線性規劃在各個領域的應用範圍。單純形法及其改進方法的持續發展，確保了線性規劃在最佳化和決策問題中的核心地位。
參考文獻
- https://zhuanlan.zhihu.com/p/672071565
- python最最佳化演算法實戰---線性規劃之單純形法
在給定的營業時間內（14小時），啤酒廠需要制定一個生產計劃，以最大化其收益。生產生啤需要1小時，生產黑啤需要2小時。生啤的售價為20美元/箱，黑啤的售價為30美元/箱。每日銷售上限為100箱。