考虑线性规划的标准形式：
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+Im?y=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=[zB?zN??].
考虑到 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?=[xB?yB??],zN?=[xN?yN??],
其中 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=[zB?zN??]=?????xB?yB?xN?yN???????.
接下来考虑两种情形：
情形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=???xB?xN?yN?????
此时，原问题的基变量为 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=?????xB?yB?xN?yN???????
这种情况下，不能直接把 ( 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?[Ik?0?R1?0?].
最后 m ? k m-k m?k??? 行是冗余的。删除它们，然后把 x B x_B xB? 作为原问题的基变量。
算法实现
import numpy as npfrom simplex_ad import SimplexAD # 之前实现的单纯形算法class Simplex2P(object):def __init__(self, c, A, b):self._c = cself._A = Aself._b = bself._m = len(self._A)self._n = len(self._c)self._status = None# The following variables are for phase 1self._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 = Nonedef _construct_phase1_instance(self):c1 = [0] * self._n + [1] * self._mI = np.eye(self._m)A = self._AA1 = np.array(np.bmat("A I"))b1 = self._bv0 = [i + self._n for i in range(self._m)]return c1, A1, b1, v0def _initialize(self):# solve phase1 problemsim = SimplexAD(*self._ins1).solve(print_info=False)# check feasibilityif 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 basicbasic_indices2 = basic_non_art_indices + basic_art_indicesself._basic_vars = [x[1] for x in basic_indices2]# row swappingB_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_N1def _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) @ NB_inv_N1 = self._swapping_columns_and_rows(B_inv_N)# The number of non basic, non artificial variablesk = 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 existrr = 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 instancereserved_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 ins1c1, A1, b1, v0 = self._ins1A1 = 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_numreturn [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 variablesin_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]# replacementfor j in range(self._m):if self._basic_vars[j] == out_var:self._basic_vars[j] = in_varbreakfor 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 + jself._R12[:, in_var] = np.array([1 if i == row_ind else 0 for i in range(self._m)])breakdef _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_indk = len(self._R12[row_id])for j in range(k):if abs(self._R12[row_id][j]) > 1e-6:return jreturn Nonedef 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.（点此下载）