考虑线性规划的标准形式:
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。
回顾单纯形算法:从一个基本可行解开始,沿着目标函数下降的方向,迭代到新的基本可行解,直到目标函数达到最优。
那么如何找到一个初始的基本可行解?本文介绍一个方法:两阶段法。
两阶段法
两阶段法的思路如下:
1、构造辅助变量,把原问题转化成“新问题”,使得新问题的基本可行解是显而易见的。
2、用单纯形算法求解新问题,从而得到原问题的一个基本可行解。
3、再次使用单纯形算法,求解原问题。
下面介绍新问题的定义。定义辅助变量 y ∈ R m y \in \mathbb{R}^m y∈Rm, e = ( 1 , 1 , … , 1 ) T ∈ R n e = (1, 1, \ldots, 1)^T \in \mathbb{R}^n e=(1,1,…,1)T∈Rn, I m I_m Im 是 m m m 维单位矩阵。
考虑如下问题(下文称它为新问题):
min
e
T
y
s.t.
A
x
+
I
m
y
=
b
x
,
y
≥
0
\begin{aligned} \min~ & e^T y\\ \text{s.t.}~ & Ax + I_m y =b\\ & x, y\geq 0 \end{aligned}
min s.t. eTyAx+Imy=bx,y≥0
原问题有可行解等价于新问题的最优目标函数值为零(证明略)。
注意到 I m I_m Im 是新问题的基,用单纯形算法计算新问题,得到最优解 x ∗ , y ∗ x^*, y^* x∗,y∗。如果原问题可行,那么 y ∗ = 0 y^*=\mathbf{0} y∗=0,意味着 A x ∗ = b A x^* = b Ax∗=b。换句话说, x ∗ x^* x∗ 是原问题的可行解。
不过问题来了, x ∗ x^* x∗ 虽然是原问题的可行解,但是不能直接当作单纯形算法的输入,因为新问题和原问题的基变量不一定相同。接下来我们讨论如何计算原问题的基变量。
基变量
为了简化记号,把新问题的最优解记作
z
z
z,即
z
=
(
x
1
,
…
,
x
n
,
y
1
,
…
,
y
m
)
T
.
z = (x_1,\ldots,x_n,y_1,\ldots,y_m)^T.
z=(x1,…,xn,y1,…,ym)T.
令
z
B
,
z
N
z_B, z_N
zB,zN 分别代表基变量和非基变量。于是
z
z
z 也可以写成
z
=
[
z
B
z
N
]
.
z = \begin{bmatrix} z_B\\ z_N \end{bmatrix}.
z=[zBzN].
考虑到
z
B
,
z
N
z_B, z_N
zB,zN 有可能包含
x
x
x 和
y
y
y 中的分量,令
z
B
=
[
x
B
y
B
]
,
z
N
=
[
x
N
y
N
]
,
z_B = \begin{bmatrix} x_B\\ y_B \end{bmatrix}, \quad z_N = \begin{bmatrix} x_N\\ y_N \end{bmatrix},
zB=[xByB],zN=[xNyN],
其中
x
B
,
y
B
x_B, y_B
xB,yB 是
x
,
y
x, y
x,y 对应的基变量,
x
N
,
y
N
x_N, y_N
xN,yN 是
x
,
y
x, y
x,y 对应的非基变量。
于是,
z
z
z 可以进一步写成
z
=
[
z
B
z
N
]
=
[
x
B
y
B
x
N
y
N
]
.
z = \begin{bmatrix} z_B \\ z_N \end{bmatrix} = \begin{bmatrix} x_B\\ y_B\\ x_N\\ y_N \end{bmatrix}.
z=[zBzN]=⎣⎢⎢⎡xByBxNyN⎦⎥⎥⎤.
接下来考虑两种情形:
情形1:基变量
z
B
z_B
zB 不包含
y
y
y 。
z
=
[
x
B
x
N
y
N
]
z = \begin{bmatrix} x_B\\ x_N\\ y_N \end{bmatrix}
z=⎣⎡xBxNyN⎦⎤
此时,原问题的基变量为
x
B
x_B
xB。
情形2:基变量
z
B
z_B
zB 包含
y
y
y 的分量。
z
=
[
x
B
y
B
x
N
y
N
]
z = \begin{bmatrix} x_B\\ y_B\\ x_N\\ y_N \end{bmatrix}
z=⎣⎢⎢⎡xByBxNyN⎦⎥⎥⎤
这种情况下,不能直接把
(
x
B
,
x
N
)
(x_B, x_N)
(xB,xN) 作为原问题的初始基本可行解,因为
x
B
x_B
xB 的维度小于
m
m
m 。
能不能把 y B y_B yB 与 x N x_N xN 对应的列交换?
答案是可以的。因为
x
N
=
0
x_N = \mathbf{0}
xN=0,选择某一列入基不会降低目标函数值。
具体来说,交换行和列的顺序,先把最优解和对应的系数矩阵写成如下形式。
交换的思路是把
I
m
−
k
I_{m-k}
Im−k 与
R
2
R_2
R2 中非零元素所在的行列,按列进行交换。
考虑矩阵
R
2
R_2
R2 的第
i
i
i 行,找到一个非零元所在的列
j
j
j, 然后交换
y
B
i
y_{B_i}
yBi 和
x
N
j
x_{N_j}
xNj。如果
j
j
j 不存在,那么
R
2
R_2
R2 的第
i
i
i 行是
0
\mathbf{0}
0,这意味着
y
B
i
y_{B_i}
yBi 对应的约束是冗余的。
从 i = 0 i = 0 i=0 开始,执行上述步骤,最终得到的结果是:
(1) R 2 R_2 R2 被转换成 I m − k I_{m-k} Im−k 或者
(2)执行到某一步时 R 2 R_2 R2 中的元素全为0导致无法继续交换。
第一种情况意味着 y B y_B yB 全部被替换,此时 z B = x B z_B = x_B zB=xB,即原问题的基变量。
第二种情况意味着原问题的系数矩阵不满秩:
A
∼
[
I
k
R
1
0
0
]
.
A \sim \begin{bmatrix} I_k & R_1 \\ \mathbf{0} & \mathbf{0} \end{bmatrix}.
A∼[Ik0R10].
最后
m
−
k
m-k
m−k 行是冗余的。删除它们,然后把
x
B
x_B
xB 作为原问题的基变量。
算法实现
import numpy as np
from simplex_ad import SimplexAD # 之前实现的单纯形算法
class Simplex2P(object):
def __init__(self, c, A, b):
self._c = c
self._A = A
self._b = b
self._m = len(self._A)
self._n = len(self._c)
self._status = None
# The following variables are for phase 1
self._ins1 = self._construct_phase1_instance()
self._basic_vars = []
self._basic_art_vars = []
self._basic_art_vars_num = 0
# 处理退化情形的中间变量
# R = B^{-1}*N.
# R12 对应 R 的前 k 列,其中 k 代表 basic non artificial variable 的个数
self._R12 = None
def _construct_phase1_instance(self):
c1 = [0] * self._n + [1] * self._m
I = np.eye(self._m)
A = self._A
A1 = np.array(np.bmat("A I"))
b1 = self._b
v0 = [i + self._n for i in range(self._m)]
return c1, A1, b1, v0
def _initialize(self):
# solve phase1 problem
sim = SimplexAD(*self._ins1).solve(print_info=False)
# check feasibility
if abs(sim.get_objective()) > 1e-6 or sim.get_status() != 'OPTIMAL':
self._status = 'INFEASIBLE'
return
# If feasible, consider the following two cases.
# Case1(normal): artificial vars not in basic vars.
self._basic_vars = sim.get_basic_vars()
# Case2(degenerate): artificial vars in basic vars.
# Then exchange artificial columns with non-basic columns.
self._basic_art_vars = [i for i in self._basic_vars if i >= self._n]
self._basic_art_vars_num = len(self._basic_art_vars)
if self._basic_art_vars_num > 0:
self._resolve_degeneracy()
def _swapping_columns_and_rows(self, B_inv_N):
basic_indices = list(zip(range(self._m), self._basic_vars))
basic_non_art_indices = [item for item in basic_indices if item[1] < self._n]
basic_art_indices = [item for item in basic_indices if item[1] >= self._n]
# column swapping
# reorder basic indices w.r.t. non basic
basic_indices2 = basic_non_art_indices + basic_art_indices
self._basic_vars = [x[1] for x in basic_indices2]
# row swapping
B_inv_N1 = np.zeros(B_inv_N.shape)
ind2 = [x[0] for x in basic_indices2]
for i in range(self._m):
B_inv_N1[i, :] = B_inv_N[ind2[i], :]
return B_inv_N1
def _format_R12(self):
B = self._get_submatrix_of_A1(self._basic_vars)
non_basic_vars = list(set(range(self._n + self._m)) - set(self._basic_vars))
N = self._get_submatrix_of_A1(non_basic_vars)
B_inv_N = np.linalg.inv(B) @ N
B_inv_N1 = self._swapping_columns_and_rows(B_inv_N)
# The number of non basic, non artificial variables
k = self._n - (self._m - self._basic_art_vars_num)
self._R12 = B_inv_N1[:, :k]
def _resolve_degeneracy(self):
self._format_R12()
# replace artificial basic columns with non-basic columns.
self._replace_artificial_basis()
# remove redundant rows if there exist
rr = self._get_redundant_rows()
if len(rr) > 0:
self._remove_redundant_rows(rr)
# remove redundant basis w.r.t. redundant rows.
reserved_indices = set(range(len(self._basic_vars))) - set(rr)
self._basic_vars = [self._basic_vars[i] for i in reserved_indices]
def _remove_redundant_rows(self, redundant_rows):
# remove from original instance
reserved_rows = set(range(self._m)) - set(redundant_rows)
self._A = np.array([self._A[i] for i in reserved_rows])
self._b = np.array([self._b[i] for i in reserved_rows])
self._m = len(self._A)
# remove from ins1
c1, A1, b1, v0 = self._ins1
A1 = np.array([A1[i] for i in reserved_rows])
b1 = np.array([b1[i] for i in reserved_rows])
self._ins1 = (c1, A1, b1, v0)
def _get_redundant_rows(self):
# return zero-vector indices w.r.t. R2.
k = len(self._basic_vars) - self._basic_art_vars_num
return [i for i in range(k, self._m) if sum(np.absolute(self._R12[i])) < 1e-6]
def _get_submatrix_of_A1(self, cols):
""" Return the sub matrix of A1 w.r.t. column indices.
"""
mat = []
A = self._ins1[1]
for j in cols:
mat.append(A[:, j])
return np.array(mat).transpose()
def _replace_artificial_basis(self):
for i in range(self._basic_art_vars_num):
# "in" variable from non basic variables
in_var = self._get_in_var(self._basic_art_vars_num, i)
if in_var is None:
continue
# "out" variable from basic variables (which is artificial)
out_var = self._basic_art_vars[i]
# replacement
for j in range(self._m):
if self._basic_vars[j] == out_var:
self._basic_vars[j] = in_var
break
for j in range(self._basic_art_vars_num):
if self._basic_art_vars[j] == out_var:
row_ind = self._m - self._basic_art_vars_num + j
self._R12[:, in_var] = np.array([1 if i == row_ind else 0 for i in range(self._m)])
break
def _get_in_var(self, basic_art_vars_num, art_ind):
"""
:param basic_art_vars_num: the total number of "basic" artificial variables
:param art_ind: the index of the artificial variable in "basic_art_vars"
:return: "in" variable from non basic variables
"""
row_id = self._m - basic_art_vars_num + art_ind
k = len(self._R12[row_id])
for j in range(k):
if abs(self._R12[row_id][j]) > 1e-6:
return j
return None
def solve(self):
self._initialize()
if self._status == 'INFEASIBLE':
print('>> INFEASIBLE.')
else:
print("Init: feasible.")
sim = SimplexAD(self._c, self._A, self._b, self._basic_vars).solve()
self._status = sim.get_status()
return self
参考文献
[1] Christopher Griffin. Linear Programming: Penn State Math 484 Lecture Notes. Version 1.8.3. Chapter 6.(点此下载)