考虑线性规划的标准形式(更多介绍参考《线性规划的标准形》):
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
其中
c
,
x
∈
R
n
c, x \in \mathbb{R}^n
c,x∈Rn,
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。
单纯形算法的思路是从多面体的一个顶点出发,然后沿着降低目标的方向,迭代到另一个顶点,直到目标值无法降低,于是得到最优解。
基本可行解
把
A
A
A 拆成两个部分:
A
=
[
B
N
]
A = [B \quad N]
A=[BN]
其中
B
∈
R
m
×
m
B\in \mathbb{R}^{m\times m}
B∈Rm×m,
N
∈
R
m
×
(
n
−
m
)
N \in \mathbb{R}^{m \times (n-m)}
N∈Rm×(n−m)。
相应地,把
c
c
c 和
x
x
x 分别拆成两个部分,即
c
=
[
c
B
c
N
]
x
=
[
x
B
x
N
]
.
c = \begin{bmatrix} c_B \\ c_N\end{bmatrix} \quad x = \begin{bmatrix} x_B \\ x_N \end{bmatrix}.
c=[cBcN]x=[xBxN].
注意:
c
B
,
x
B
c_B, x_B
cB,xB 的下标 与
B
B
B 的列分别对应。
假设 A A A 满秩,我们可以交换列的位置保证 B B B 可逆,所以 B B B 的列向量是一组 基。 B B B 称为 基矩阵, N N N 称为 非基矩阵。 把 x B x_B xB 对应的变量称为 基变量, x N x_N xN 对应的变量称为 非基变量。
这样一来,约束条件可以写成
[
B
N
]
[
x
B
x
N
]
=
b
⇔
B
x
B
+
N
x
N
=
b
.
[B \quad N] \begin{bmatrix} x_B\\ x_N \end{bmatrix} = b \Leftrightarrow Bx_B + N x_N = b.
[BN][xBxN]=b⇔BxB+NxN=b.
我们有
x
B
=
B
−
1
b
−
B
−
1
N
x
N
.
x_B = B^{-1}b - B^{-1}Nx_N.
xB=B−1b−B−1NxN.
令
x
N
=
0
x_N = 0
xN=0,得到 基本解
x
0
=
[
B
−
1
b
0
]
.
x_0 = \begin{bmatrix} B^{-1}b\\ 0 \end{bmatrix}.
x0=[B−1b0].
换句话说,基本解
x
0
x_0
x0 满足约束条件
A
x
0
=
b
Ax_0 = b
Ax0=b,但有可能不可行。如果
x
0
≥
0
x_0 \geq \mathbf{0}
x0≥0,那么称它为 基本可行解。换句话说,基本可行解是多面体的一个顶点(证明略)。
假设已知一个初始的基本可行解,接下来我们要想办法迭代到另一个基本可行解,并且使得目标函数值下降。
迭代条件
要想目标函数降低,我们自然想到梯度下降法。令 f ( x ) f(x) f(x) 代表目标函数,梯度 ∇ f ( x ) \nabla f(x) ∇f(x) 就是函数值增加的方向。所以应该沿着负梯度的方向移动 x x x。
把
x
B
x_B
xB 代入目标函数
f
(
x
)
:
=
c
T
x
f(x):=c^Tx
f(x):=cTx 可得
f
(
x
)
=
c
T
x
=
c
B
T
x
B
+
c
N
T
x
N
=
c
B
T
B
−
1
b
+
(
c
N
T
−
c
B
T
B
−
1
N
)
x
N
.
\begin{aligned} f(x) & = c^Tx \\ & = c_B^Tx_B + c_N^T x_N \\ & = c_B^T B^{-1}b + ( c_N^T - c_B^T B^{-1}N) x_N. \end{aligned}
f(x)=cTx=cBTxB+cNTxN=cBTB−1b+(cNT−cBTB−1N)xN.
对
f
(
x
)
f(x)
f(x) 关于
b
b
b 求导, 得到 Shadow Price (也称为 对偶变量):
λ
T
:
=
c
B
T
B
−
1
.
\lambda^T := c_B^T B^{-1}.
λT:=cBTB−1.
对
f
(
x
)
f(x)
f(x) 关于
x
N
x_N
xN 求导(求梯度),得到 Reduced Cost :
μ
T
:
=
c
N
T
−
c
B
T
B
−
1
N
=
c
N
T
−
λ
T
N
.
\mu^T := c_N^T - c_B^TB^{-1}N = c_N^T - \lambda^T N.
μT:=cNT−cBTB−1N=cNT−λTN.
令
J
J
J 代表非基变量的下标。
∀
j
∈
J
\forall j\in J
∀j∈J,
x
j
x_j
xj 每增加一个单位,目标函数就增加
μ
j
\mu_j
μj。
如果存在 μ j < 0 \mu_j < 0 μj<0,我们只要把 x j x_j xj 从 0 增加 δ \delta δ,那么目标函数就可以降低 − μ j δ -\mu_j \delta −μjδ。
当 μ j ≥ 0 \mu_j \geq 0 μj≥0, ∀ j ∈ J \forall j\in J ∀j∈J,这意味着目标函数值无法降低,此时达到最优解(证明略)。
如何迭代
假设存在 j ∈ J j\in J j∈J 使得 μ j < 0 \mu_j < 0 μj<0,于是我们想增加 x j x_j xj 使得目标函数降低。但是问题来了,增加多少可以保证 x x x 可行?
回顾基本可行解的定义
x
=
[
x
B
x
N
]
=
[
B
−
1
b
−
B
−
1
N
x
N
x
N
]
.
x = \begin{bmatrix} x_B\\ x_N \end{bmatrix} = \begin{bmatrix} B^{-1}b - B^{-1}Nx_N \\ x_N \end{bmatrix}.
x=[xBxN]=[B−1b−B−1NxNxN].
注意到
x
N
=
0
x_N = 0
xN=0,我们只要保证
x
B
≥
0
x_B \geq 0
xB≥0。接下来把
x
B
x_B
xB 换一种写法,令
b
~
=
B
−
1
b
\tilde{b} = B^{-1}b
b~=B−1b,
a
~
j
=
B
−
1
a
j
\tilde{a}_j = B^{-1}a_j
a~j=B−1aj,我们有
x
B
=
[
x
B
1
x
B
2
⋮
x
B
m
]
=
[
b
~
1
b
~
1
⋮
b
~
m
]
−
[
a
~
j
1
a
~
j
2
⋮
a
~
j
m
]
x
j
=
[
b
~
1
−
a
~
j
1
x
j
b
~
2
−
a
~
j
2
x
j
⋮
b
~
m
−
a
~
j
m
x
j
]
.
x_B = \begin{bmatrix} x_{B_1}\\ x_{B_2}\\ \vdots\\ x_{B_m}\end{bmatrix} = \begin{bmatrix} \tilde{b}_1\\ \tilde{b}_1\\ \vdots\\ \tilde{b}_m \end{bmatrix} - \begin{bmatrix} \tilde{a}_{j_1}\\ \tilde{a}_{j_2}\\ \vdots\\ \tilde{a}_{j_m} \end{bmatrix}x_j = \begin{bmatrix} \tilde{b}_1 - \tilde{a}_{j_1}x_j\\ \tilde{b}_2 - \tilde{a}_{j_2}x_j\\ \vdots\\ \tilde{b}_m - \tilde{a}_{j_m}x_j \end{bmatrix}.
xB=⎣⎢⎢⎢⎡xB1xB2⋮xBm⎦⎥⎥⎥⎤=⎣⎢⎢⎢⎡b~1b~1⋮b~m⎦⎥⎥⎥⎤−⎣⎢⎢⎢⎡a~j1a~j2⋮a~jm⎦⎥⎥⎥⎤xj=⎣⎢⎢⎢⎡b~1−a~j1xjb~2−a~j2xj⋮b~m−a~jmxj⎦⎥⎥⎥⎤.
现在要增加
x
j
x_j
xj, 那么只要保证
x
B
x_B
xB 的每个分量非负即可。
注意到
x
0
x_0
x0 是基本可行解,所以
b
~
=
B
−
1
b
≥
0
\tilde{b} = B^{-1}b \geq 0
b~=B−1b≥0。容易验证,下面的取值(称之为 Minimum Ratio Test)可以保证
x
B
≥
0
x_B \geq 0
xB≥0。
x
j
:
=
min
{
b
~
i
a
~
j
i
and
a
~
j
i
>
0
,
i
=
1
,
2
,
…
,
m
}
.
x_j := \min \left\{ \frac{\tilde{b}_i}{\tilde{a}_{j_i}} \text{ and } \tilde{a}_{j_i}>0, \quad i=1,2,\ldots, m\right\}.
xj:=min{a~jib~i and a~ji>0,i=1,2,…,m}.
这样一来,
x
j
≥
0
x_j \geq 0
xj≥0,相应地,另一个变量
x
B
i
=
b
~
i
−
a
~
j
i
=
0
x_{B_i} = \tilde{b}_i - \tilde{a}_{j_i} = 0
xBi=b~i−a~ji=0, 其中
i
i
i 是上式达到最小值的下标。换句话说,此时
x
j
x_j
xj 成为基变量(入基),而
x
B
i
x_{B_i}
xBi 成为非基变量(出基)。可以证明,新得到的解仍然是一个基本可行解。
注意:如果 a ~ j i ≤ 0 \tilde{a}_{j_i} \leq 0 a~ji≤0, ∀ i \forall i ∀i ,这意味着 x j x_j xj 可以无限大,最优目标函数值为 − ∞ -\infty −∞,于是最优解不存在。
算法描述
第0步:输入基本可行解对应的基矩阵 B B B。
第1步:判断当前的解是否最优。如果 μ ≥ 0 \mu \geq 0 μ≥0,当前是最优解,算法停止。
第2步:计算入基变量和出基变量。如果存在 j ∈ J j\in J j∈J 使得 μ j < 0 \mu_j < 0 μj<0,那么 x j x_j xj 是入基变量。执行 Minimum Ratio Test, 找到出基变量 x i x_i xi。
第3步:判断问题是否无界。如果 a ~ j i ≤ 0 \tilde{a}_{j_i} \leq 0 a~ji≤0, ∀ i \forall i ∀i,则代表无界,算法停止。
第4步:执行出入基操作,更新基矩阵 B B B,然后执行第1步。
算法实现
下面我们用Python来实现单纯形算法。
先定义算法的输入和输出。
class SimplexA(object):
"""
单纯形算法(基本版)。
Note:
1、系数矩阵满秩。
2、未处理退化情形。
3、输入基本可行解(对应的列)。
"""
def __init__(self, c, A, b, v0):
"""
:param c: n * 1 vector
:param A: m * n matrix
:param b: m * 1 vector
:param v0: basic variables, list of variable indices
注意:v0是 B 的列下标。x0 = B^{-1}b 即为基本可行解(需要保证x0非负)。
"""
# 输入
self._c = np.array(c)
self._A = np.array(A)
self._b = np.array(b)
self._basic_vars = v0 # basic variables
self._m = len(A)
self._n = len(c)
self._non_basic_vars = self._init_non_basic_vars() # non basic variables
# 辅助变量
self._iter_num = 0
self._B_inv = None # inverse of B
self._lambda = None # shadow price
self._mu = None # reduced cost
# 输出
self._obj = None # objective function value
self._sol = None # solution
self._status = None
接下来要实现单纯形算法 SimplexA.solve()
,思路如下。
class SimplexA(object):
# ...
# 其它函数省略……
def solve(self):
self._iter_num = 0 # 记录迭代次数
self._check_init_solution() # 检查初始基本解是否可行
self._update_reduced_cost()
self._update_obj()
self._update_solution()
while not self._is_optimal(): # 判断是否最优或者无界
if self._status == "UNBOUNDED":
break
self._pivot() # 迭代(选主元入基,执行Minimum Ratio Test,然后出基)
self._update_reduced_cost() # 更新Reduced Cost: mu = c_N - lambda * N
self._update_obj()
self._update_solution()
self._iter_num += 1
if self._status != "UNBOUNDED":
self._status = 'OPTIMAL'
print("Done >> status: {}".format(self._status))
return self
目标函数值和Reduced Cost的计算可以直接套公式。关键是实现每次迭代的出入基操作,即 SimplexA._pivot()
。
class SimplexA(object):
# ...
# 其它函数省略……
def _pivot(self):
""" 选主元,入基和出基。
"""
j_ind = np.argmin(self._mu)
# 入基变量 x_j
j = self._non_basic_vars[j_ind]
# 出基变量 x_i
i = self._minimum_ratio_test(j)
if i is None:
self._status = 'UNBOUNDED'
return
# 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 _minimum_ratio_test(self, j):
""" Minimum Ratio Test.
给定入基的非基变量,返回出基的基变量。
:param j: 入基变量 x_j 的下标 j
:return: 出基变量 x_i 的下标 i
"""
b_bar = np.dot(self._B_inv, self._b)
a_in = np.dot(self._B_inv, self._A[:, j])
ratios = list(map(lambda b, a: b/a if a > 1e-6 else np.infty, b_bar, a_in))
i_ind = np.argmin(ratios)
if ratios[i_ind] != np.infty:
return self._basic_vars[i_ind]
else:
return None
遗留问题
需要注意的是,前面描述的算法并没有完全解决问题。比如:初始的基本可行解如何找到?矩阵 A A A 不是满秩如何处理?
此外,迭代过程有没有可能出现死循环?比如初始解为 x 0 x_0 x0, 经过一系列迭代又回到 x 0 x_0 x0。遗憾的是,这种情况有可能发生,称为 退化情形。
矩阵 A A A 不满秩时意味着有冗余的约束。如果已知初始的基本可行解,于是已知 A A A 的基,那么去掉冗余的约束即可。综上所述,还有两个问题有待解决:
- 计算一个初始基本可行解;
- 处理退化情形。
我们将在后续的教程中介绍处理的方法。
参考文献
[1] Christopher Griffin. Linear Programming: Penn State Math 484 Lecture Notes. Version 1.8.3. Chapter 5 (点此下载)