單純形算法從一個基本可行解出發，朝着目标函數值下降的方向疊代，直到最優。從對偶的角度來看，原問題目标函數下降的方向，就是對偶問題的對偶解可行的方向，當對偶解可行時，目标函數達到最優。
本文介紹對偶單純形法。它的思路是從對偶可行解出發，朝着原問題可行的方向疊代，直到原問題可行，于是得到最優解。
對偶可行解
考慮線性規劃問題：
min c T x s.t. A x = b x ≥ 0 \begin{aligned} \min~ & c^T x \\ \text{s.t.}~ & Ax=b\\ & x\geq 0 \end{aligned} min s.t. cTxAx=bx≥0
和它的對偶問題：
max b T y s.t. A T y ≤ c \begin{aligned} \max~ & b^T y \\ \text{s.t.}~ & A^T y \leq c \end{aligned} max s.t. bTyATy≤c
其中 c , x ∈ R n c, x \in \mathbb{R}^n c,x∈Rn， y ∈ R m y\in\mathbb{R}^m y∈Rm， A ∈ R m × n A\in\mathbb{R}^{m\times n} A∈Rm×n， b ∈ R m ≥ 0 b\in \mathbb{R}^m \geq \mathbf{0} b∈Rm≥0。
令 B , N B, N B,N 分别代表基矩陣和非基矩陣。對偶問題可以寫成如下形式：
max b T y s.t. B T y ≤ c B N T y ≤ c N \begin{aligned} \max~ & b^T y \\ \text{s.t.}~ & B^T y \leq c_B\\ & N^T y \leq c_N \end{aligned} max s.t. bTyBTy≤cBNTy≤cN
回顧原問題的基本解 x x x，稱為 原始解，是這樣的形式：
x = [ x B x N ] = [ B − 1 b 0 ] . x = \begin{bmatrix} x_B\\ x_N \end{bmatrix} = \begin{bmatrix} B^{-1}b\\ \mathbf{0} \end{bmatrix}. x=[xBxN]=[B−1b0].
當 B − 1 b ≥ 0 B^{-1}b \geq 0 B−1b≥0 時， x ≥ 0 x \geq 0 x≥0，它就是基本可行解。
再看對偶問題，它的基本解 y y y，稱為 對偶解，是這樣的形式：
y = ( B T ) − 1 c B . y = (B^T)^{-1}c_B. y=(BT)−1cB.
檢查對偶問題的兩個限制：第一個 B T y 0 ≤ c B B^Ty_0\leq c_B BTy0≤cB 自然成立；第二個限制 c N − N T y ≥ 0 c_N - N^T y \geq 0 cN−NTy≥0，即
c N − N T y = c N − N T ( B T ) − 1 c B = c N − N T ( B − 1 ) T c B = c N − ( B − 1 N ) T c B ≥ 0. \begin{aligned} c_N - N^T y &= c_N - N^T (B^T)^{-1} c_B\\ & = c_N - N^T (B^{-1})^T c_B\\ &= c_N - (B^{-1}N)^T c_B \\ & \geq 0. \end{aligned} cN−NTy=cN−NT(BT)−1cB=cN−NT(B−1)TcB=cN−(B−1N)TcB≥0.
令 μ : = c N − ( B − 1 N ) T c B \mu := c_N - (B^{-1}N)^T c_B μ:=cN−(B−1N)TcB，如果 μ ≥ 0 \mu \geq 0 μ≥0，那麼 y y y 就是對偶問題的基本可行解。
回顧單純形算法， μ \mu μ 實際上是原問題的 Reduced Cost，即原問題目标函數關于 x N x_N xN 的導數。是以，對偶可行意味着原始“最優”，即原問題的目标函數值關于 x N x_N xN 無法降低。如果原始解 x x x 也是可行的，即 B − 1 b ≥ 0 B^{-1}b \geq 0 B−1b≥0，那麼 x x x 和 y y y 分别為原問題和對偶問題的最優解。
如何找到初始的對偶可行解？本文不做介紹，詳情見文末的參考文獻[1]。
如何疊代
對偶可行保證原始解的最優性，但可能損失可行性。疊代的思路就讓原始解滿足最優性的同時，也變得可行。
已知基矩陣 B B B，我們介紹出入基的規則。
已知對偶可行解 y = ( B T ) − 1 c B y = (B^T)^{-1}c_B y=(BT)−1cB。令 b ~ : = B − 1 b \tilde{b}:= B^{-1}b b~:=B−1b，如果 b ~ ≥ 0 \tilde{b}\geq 0 b~≥0，那麼原始解可行， x x x 和 y y y 即為原問題和對偶問題的最優解。
否則存在分量 b ~ i < 0 \tilde{b}_i < 0 b~i<0 。為了讓 x x x 變得可行，自然的想法是讓 B B B 的第 i i i 行對應的變量 x B i x_{B_i} xBi 出基，是以 x B i x_{B_i} xBi 是出基變量。
那麼入基變量如何計算？
先把對偶問題換一種寫法。令 y ~ = B T y \tilde{y} = B^T y y~=BTy，于是 y = ( B T ) − 1 y ~ y=(B^T)^{-1}\tilde{y} y=(BT)−1y~。代入對偶問題，得到下面的等價形式：
max b ~ T y ~ s.t. y ~ ≤ c B ( B − 1 N ) T y ~ ≤ c N \begin{aligned} \max~ & \tilde{b}^T \tilde{y} \\ \text{s.t.}~ & \tilde{y} \leq c_B\\ & (B^{-1}N)^T \tilde{y} \leq c_N \end{aligned} max s.t. b~Ty~y~≤cB(B−1N)Ty~≤cN
前面假設 b ~ i < 0 \tilde{b}_i < 0 b~i<0，那麼減少 y ~ i \tilde{y}_i y~i 可以增加目标函數值，同時需要保證
( B − 1 N ) T y ~ ≤ c N . (B^{-1}N)^T \tilde{y} \leq c_N. (B−1N)Ty~≤cN.
令 J J J 代表非基變量的下标。 ∀ j ∈ J \forall j\in J ∀j∈J，令 a ~ j = B − 1 a j \tilde{a}_j = B^{-1}a_j a~j=B−1aj，其中 a j a_j aj 代表 A A A 的第 j j j 列。上面的限制條件可以寫成：
a ~ j T y ~ ≤ c j , ∀ j ∈ J . \tilde{a}_j^T \tilde{y} \leq c_j,\quad \forall j\in J. a~jTy~≤cj,∀j∈J.
令 y ~ i \tilde{y}_i y~i 減小 δ \delta δ，其中 δ > 0 \delta > 0 δ>0，我們有
a ~ j T y ~ − a ~ i j δ ≤ c j , ∀ j ∈ J . \tilde{a}_j^T \tilde{y} - \tilde{a}_{ij}\delta \leq c_j,\quad \forall j\in J. a~jTy~−a~ijδ≤cj,∀j∈J.
如果 a ~ i j > 0 \tilde{a}_{ij} > 0 a~ij>0， ∀ j ∈ J \forall j\in J ∀j∈J，那麼 δ \delta δ 可以為無窮大，對偶問題的最優目标函數值為 + ∞ +\infty +∞，那麼原問題無可行解。
假設存在 j ∈ J j\in J j∈J，使得 a ~ i j < 0 \tilde{a}_{ij} < 0 a~ij<0。為了滿足限制條件，我們有
δ : = min { c j − a ~ j T y ~ − a ~ i j and a ~ i j < 0 , ∀ j ∈ J } . \delta := \min \left\{\frac{c_j - \tilde{a}^T_j \tilde{y}}{-\tilde{a}_{ij}} \text{ and } \tilde{a}_{ij} < 0, ~\forall j\in J\right\}. δ:=min{−a~ijcj−a~jTy~ and a~ij<0, ∀j∈J}.
注意到
a ~ j T y ~ = ( B − 1 a j ) T B T y = a j T ( B − 1 ) T B T y = a j T y , \begin{aligned} \tilde{a}^T_j \tilde{y} & = (B^{-1}a_j)^T B^Ty \\ & = a_j^T (B^{-1})^TB^T y\\ & = a^T_j y, \end{aligned} a~jTy~=(B−1aj)TBTy=ajT(B−1)TBTy=ajTy,
我們有
δ = min { c j − a j T y − a ~ i j and a ~ i j < 0 , ∀ j ∈ J } . \delta = \min \left\{\frac{c_j - a^T_j y}{-\tilde{a}_{ij}} \text{ and } \tilde{a}_{ij} < 0, ~\forall j\in J\right\}. δ=min{−a~ijcj−ajTy and a~ij<0, ∀j∈J}.
上式稱為 Minimum Ratio Test，取得最小值對應的下标 j j j 即為入基變量 x j x_j xj 的下标。
算法描述
第0步：輸入對偶可行的基矩陣 B B B。
第1步：判斷原始解是否可行。如果 b ~ ≥ 0 \tilde{b} \geq 0 b~≥0，目前是最優解，算法停止。
第2步：計算入基變量和出基變量。如果存在 b ~ i < 0 \tilde{b}_i < 0 b~i<0，那麼 x i x_i xi 是出基變量。如果存在 j ∈ J j \in J j∈J 使得 a ~ i j < 0 \tilde{a}_{ij} < 0 a~ij<0，根據 Minimum Ratio Test， 找到入基變量 x j x_j xj。
第3步：判斷問題是否無界。如果 a ~ i j ≥ 0 \tilde{a}_{ij} \geq 0 a~ij≥0， ∀ i \forall i ∀i，則對偶問題無界，原問題無可行解，算法停止。
第4步：執行出入基操作，更新基矩陣 B B B，然後執行第1步。
算法實作
下面我們用Python來實作對偶單純形算法。
先定義算法的輸入和輸出。
class DualSimplex(object):
"""
對偶單純形算法（基本版）。
Note:
1、系數矩陣滿秩。
2、未處理退化情形。
3、輸入對偶可行解（對應的列）。
"""
def __init__(self, c, A, b, v1):
"""
:param c: n * 1 vector
:param A: m * n matrix
:param b: m * 1 vector
:param v1: dual feasible solution, list of column indices
"""
# 輸入
self._c = np.array(c)
self._A = np.array(A)
self._b = np.array(b)
self._basic_vars = v0
self._m = len(A)
self._n = len(c)
self._non_basic_vars = self._init_non_basic_vars()
# 輔助變量
self._iter_num = 0
self._B_inv = None
self._N_bar = None # N_bar = B^{-1}N
self._y = None # dual solution (i.e., shadow price)
# 輸出
self._sol = None # primal solution
self._obj_dual = None # dual objective
self._status = None
接下來要實作對偶單純形算法
DualSimplex.solve()
，思路如下。
class DualSimplex(object):
# ...
# 其它函數省略……
def solve(self):
self._iter_num = 0 # 記錄疊代次數
self._check_init_solution() # 檢查初始的對偶解是否可行
self._update_solutions()
self._update_obj()
self._print_info()
while not self._is_optimal(): # 判斷是否最優或者不可行
if self._status == "INFEASIBLE":
break
self._pivot() # 疊代（選主元入基，執行Minimum Ratio Test，然後出基）
self._update_solutions()
self._update_obj()
self._iter_num += 1
if self._status != 'INFEASIBLE':
self._status = 'OPTIMAL'
上面的關鍵步驟是實作每次疊代的出入基操作，即
SimplexA._pivot()
。
class DualSimplex(object):
# ...
# 其它函數省略……
def _pivot(self):
# 出基變量 x_i
i = int(np.argmin(self._sol))
# 出基變量 x_i 在 self._basic_vars 中的index
i_ind = None
for k in range(self._m):
if self._basic_vars[k] == i:
i_ind = k
break
# 判斷原問題是否可行
# 對偶問題無界，則原問題不可行
if not self._check_feasibility(i_ind):
self._status = 'INFEASIBLE'
return
# 入基變量 x_j 在 self._basic_vars 中的index
j_ind = self._minimum_ratio_test(i_ind)
# 入基變量 x_j
j = self._non_basic_vars[j_ind]
# update basic vars
for k in range(self._m):
if self._basic_vars[k] == i:
self._basic_vars[k] = j
break
# update non basic vars
self._non_basic_vars[j_ind] = i
def _check_feasibility(self, row_ind):
N = np.array([self._A[:, j] for j in self._non_basic_vars]).transpose()
self._N_bar = self._B_inv @ N
for x in self._N_bar[row_ind]:
if x < 0:
return True
return False
def _minimum_ratio_test(self, ind_out):
N = np.array([self._A[:, j] for j in self._non_basic_vars]).transpose()
c_N = np.array([self._c[j] for j in self._non_basic_vars])
c_bar = c_N - self._y @ N
self._N_bar = self._B_inv @ N
a_bar = self._N_bar[ind_out] * -1
ratios = list(map(lambda c, a: c/a if a > 0 else np.infty, c_bar, a_bar))
return int(np.argmin(ratios))
完整代碼
結語
本文介紹了對偶單純形法，但是還有兩個問題沒有解決：一個是對偶可行解的初始化；另一個處理退化情形。解決思路與單純形法的處理思路類似，本文不做介紹，詳情可以看下面的參考文獻[1]。
我們已經知道了單純形算法，為什麼還需要對偶單純形算法？它有什麼樣的應用場景？
具體來說，有如下兩點：
1、對偶單純形法在實際中表現不錯。不僅如此，原問題如果退化，其對偶問題可能是正常的。
2、對偶單純形法可應用于整數規劃的求解。分支定界是求解整數規劃問題的經典方法，每個分支需要求解一個線性規劃子問題。子問題相比原來問題，增加了新的限制，一般導緻原始解不可行，卻不會違背對偶可行。于是可以用之前的對偶解作為子問題的初始可行解。
參考文獻
[1] Mihai Banciu. Dual Simplex (lecture notes). Bucknell University, 2010.
[2] David P. Williamson. ORIE 6300 Mathematical Programming I (lecture notes). 2014.
![平衡優化算法（EO）結合粘液黴菌算法（SMA）（Matlab代碼實作）[圖]](data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACwAAAAAAQABAAACAkQBADs=)