运筹系列2：线性规划两阶段法python代码

提示：本文参考了scipy的linprog源码，对源码感兴趣的小伙伴可以直接去读源码，注释真的是非常详尽了，比代码都长。

1. 补充问题

上一节中的代码在运行时还有很多细节没有处理，这里补充两个比较重要的情况：

1. 存在等式约束
   如果有等式约束，那么就没法通过添加松弛变量直接给出初始可行解，需要用大M法或者两阶段法求解。计算机求解一般使用二阶段法，即首先给等式约束条件添加人工变量，使得问题有一个初始可行解。二阶段法是第一阶段最小化人工变量的和使其等于0，而大M法则是给人工变量添加一个非常大的系数M使其等于0。
   注意人工变量和松弛变量/剩余变量不一样。松弛变量/剩余变量是用于将不等式变为等式（标准形），而人工变量则是在等式约束中额外添加一个变量，这个变量最后一定要等于0才行，否则就是没有初始可行解，不能进入第二阶段。
   在第一阶段，每个不等式约束对应一个松弛变量，每个等式约束对应一个人工变量。另外，b<0的不等式约束也需要添加人工变量，因为对应的松弛变量等于b不能作为基变量。第一阶段需要添加一个额外的目标函数：$\min \Sigma x_{av}$，$x_{av}$指的是所有的人工变量。
2. RHS <0
   如果right hand side（也就是b）小于0，则这一行也无法给出初始可行解。此时的一般做法是左右都乘以-1，添加松弛变量变为等式约束，然后添加人工变量得到初始可行解。
3. 没有有限最优解（目标函数无界）
   对应的情况是：目标函数中有某个系数大于0，对应变量在约束条件中的系数全都≤0，问题结束，没有最优解。

首先假设问题除了不等式约束，还有等式约束：
求解问题为：
min $cx$
s.t.
$A_{ub}x\le b_{ub}$
$A_{eq}x=b_{eq}$

2. 基于scipy的python算法实现

下面代码的核心步骤是构造两阶段的标准单纯形矩阵，具体求解函数可参考运筹系列1，这里我们直接调用scipy的_solve_simplex函数。

首先定义一些变量：
T：用于求解的矩阵。第二阶段比第一阶段要多一行。
basis：基向量的下标集合，通过solve_simplex函数改变
maxiter：最大迭代次数
tol：求解精度
OptimizeResult：格式输出函数
solve_simplex：单纯形法求解函数。

完整的求解步骤如下：

from scipy.optimize._linprog_simplex import _solve_simplex # 求解标准单纯形矩阵的函数
from scipy.optimize import OptimizeResult # 格式化输出结果的函数
def linprog(c, A_ub=None, b_ub=None, A_eq=None, b_eq=None, maxiter=1000, tol=1.0E-12):
    cc = np.asarray(c)
    f0 = 0
    n = len(c)
    mub = len(b_ub)
    meq = len(b_eq)
    m = mub + meq  # 变量总数
    n_slack = mub  # 不等号添加松弛变量
    n_artificial = meq + np.count_nonzero(b_ub < 0)  # 等号/没有初始可行解的不等号添加人工变量
    Aub_rows, Aub_cols = A_ub.shape
    Aeq_rows, Aeq_cols = A_eq.shape
    slcount = 0
    avcount = 0
    status = 0
    messages = {0: "Optimization terminated successfully.",
                1: "Iteration limit reached.",
                2: "Optimization failed. Unable to find a feasible"
                   " starting point.",
                3: "Optimization failed. The problem appears to be unbounded.",
                4: "Optimization failed. Singular matrix encountered."}
                
    # 建立第一阶段单纯型表T。
    T = np.zeros([m + 2, n + n_slack + n_artificial + 1])  # 添加两行和一列
    T[-2, :n] = cc  # T倒数第二行用来记系数
    T[-2, -1] = f0
    b = T[:-2, -1]  # T倒数第一列用来记b
    if meq > 0:  # 给T赋值
        T[:meq, :n] = A_eq
        b[:meq] = b_eq
    if mub > 0:
        T[meq:meq + mub, :n] = A_ub
        b[meq:meq + mub] = b_ub
        np.fill_diagonal(T[meq:m, n:n + n_slack], 1) # 剩余变量对应系数1

    basis = np.zeros(m, dtype=int)
    r_artificial = np.zeros(n_artificial, dtype=int)
    for i in range(m):
        if i < meq or b[i] < 0:
            basis[i] = n + n_slack + avcount
            r_artificial[avcount] = i
            avcount += 1
            if b[i] < 0:
                b[i] *= -1
                T[i, :-1] *= -1
            T[i, basis[i]] = 1
            T[-1, basis[i]] = 1 # T的最后一行是新增加的约束条件
        else:
            basis[i] = n + slcount
            slcount += 1
    for r in r_artificial:
        T[-1, :] = T[-1, :] - T[r, :] # 将人工变量的1全部变为0，T转为标准形。至此T构造完毕。

    # 第一阶段求解
    nit1, status = _solve_simplex(T, n, basis)

    # 如果第一阶段的目标函数为0，则删除人工变量，进入第二阶段
    if abs(T[-1, -1]) <= tol:
        T = T[:-1, :]
        T = np.delete(T, np.s_[n + n_slack:n + n_slack + n_artificial], 1)
    else:
        status = 2

    if status != 0:
        message = messages[status]
        return OptimizeResult(x=np.nan, fun=-T[-1, -1], nit=nit1, status=status, message=message, success=False)
    nit2, status = _solve_simplex(T, n, basis, maxiter=maxiter - nit1, phase=2, tol=tol, nit0=nit1)

    if status != 0:
        message = messages[status]
        return OptimizeResult(x=np.nan, fun=-T[-1, -1], nit=nit2, status=status, message=message, success=False)
    solution = np.zeros(n + n_slack + n_artificial)
    solution[basis[:m]] = T[:m, -1]
    x = solution[:n]
    slack = solution[n:n + n_slack]
    obj = -T[-1, -1]
    return OptimizeResult(x=x, fun=obj, nit=int(nit2), status=status, slack=slack, message=messages[status],
                          success=(status == 0))

3. 例子说明

下面是例子：

import numpy as np
c=np.array([2,3,-5])
A_ub=np.array([[-2,5,-1]])
B_ub=np.array([-10])
A_eq=np.array([[1,1,1]])
B_eq=np.array([7])
res=linprog(-c,A_ub,B_ub,A_eq,B_eq)
print(res)

对应 $\max 2x_1+3x_2-5x_3$
s.t.
$-2x_1+5x_2-x_3\le -10$
$x_1+x_2+x_3=7$

输出为

     fun: -14.571428571428571
 message: 'Optimization terminated successfully.'
     nit: 2
   slack: array([0.])
  status: 0
 success: True
       x: array([6.42857143, 0.57142857, 0.        ])

下面具体说明一下中间过程：
首先，将max转为min问题，需要将c取负数，最后的结果相应的也取负数即可。
其次，转变为标准形，添加了2个slack variable，添加了1个artificial variable（-10小于0，因此需要一个人工变量）。
第一阶段的问题需要添加一个额外的约束，如下图：

此时矩阵T为：

调用标准单纯形法进行求解，最后的T如下：

其中T[-1,-1]代表第一阶段的最优值为0，可以继续进行第二阶段的求解。删除最后一行，以及标号为4、5的两列，求解后如下表：

因为只有2个约束条件，因此基变量也只有2个，剩下的两个变量=0。

4. 代码跟踪

我们如果用运筹系列1的代码进行跟踪，第一阶段：
初始化：

做一个简单的转换，将最后添加的两个人工变量变成基变量：

第一次pivoting：

第二次pivoting：

接下来进入第二阶段，首先剔除人工变量所在的列，然后将最后一行替换为原来的c:

进行消元，生成初始可行解：

然后将d，s送入求解器，直接就得到最优解：

参考文档

[1] Dantzig, George B., Linear programming and extensions. Rand Corporation Research Study Princeton Univ. Press, Princeton, NJ, 1963
[2] Hillier, S.H. and Lieberman, G.J. (1995), “Introduction to Mathematical Programming”, McGraw-Hill, Chapter 4.
[3] Bland, Robert G. New finite pivoting rules for the simplex method. Mathematics of Operations Research (2), 1977: pp. 103-107.
[4] Andersen, Erling D., and Knud D. Andersen. “The MOSEK interior point optimizer for linear programming: an implementation of the homogeneous algorithm.” High performance optimization. Springer US, 2000. 197-232.
[5] Andersen, Erling D. “Finding all linearly dependent rows in large-scale linear programming.” Optimization Methods and Software 6.3 (1995): 219-227.
[6] Freund, Robert M. “Primal-Dual Interior-Point Methods for Linear Programming based on Newton’s Method.” Unpublished Course Notes, March 2004. Available 2/25/2017 at https://ocw.mit.edu/courses/sloan-school-of-management/15-084j-nonlinear-programming-spring-2004/lecture-notes/lec14_int_pt_mthd.pdf
[7] Fourer, Robert. “Solving Linear Programs by Interior-Point Methods.” Unpublished Course Notes, August 26, 2005. Available 2/25/2017 at http://www.4er.org/CourseNotes/Book%20B/B-III.pdf
[8] Andersen, Erling D., and Knud D. Andersen. “Presolving in linear programming.” Mathematical Programming 71.2 (1995): 221-245.
[9] Bertsimas, Dimitris, and J. Tsitsiklis. “Introduction to linear programming.” Athena Scientific 1 (1997): 997.
[10] Andersen, Erling D., et al. Implementation of interior point methods for large scale linear programming. HEC/Universite de Geneve, 1996.