最优化方法Python计算：求解一般形式线性规划问题的两阶段方法

利用博文《最优化方法Python计算：线性规划的标准化》定义的standardizing函数、《最优化方法Python计算：标准型线性规划的单纯形算法》定义的simplex函数、《最优化方法Python计算：标准型线性规划的辅助问题》定义的prepro函数和《最优化方法Python计算：辅助问题最优解的预后处理》定义的prognostic函数装配成如下的求解一般形式线性规划
$$\{\begin{array}{l}\text{minimize}{\mathbf{c}}^{\top }\mathbf{x}\\ \text{s.t. }{\mathbf{A}}_{iq}\mathbf{x}\le {\mathbf{b}}_{iq}\\ {\mathbf{A}}_{eq}\mathbf{x}={\mathbf{b}}_{eq}\\ \mathbf{x}\ge \mathbf{o}\end{array}(1)$$
的linprog函数

import numpy as np
def linprog(c, Aiq=None, biq=None, Aeq=None, beq=None):
    if isinstance(Aiq, np.ndarray):
        nx = Aiq.shape[1]
    else:
        nx = Aeq.shape[1]
    A, b = standardizing(Aiq, biq, Aeq, beq)
    E, Bind1 = prepro(A)
    m, n = A.shape
    c = np.concatenate((c, np.zeros(n - c.size)), axis = 0)
    n1 = E.shape[1]
    A1 = np.hstack((A, E))
    Nind1 = np.setdiff1d(np.arange(n + n1), Bind1)
    if n1 > 0:
        c1 = np.concatenate((np.zeros(n), np.ones(n1)), axis = 0)
        Bindst = simplex(A1, b, c1, Bind1, Nind1)
        Bst = A1[:, Bindst]
        Bst1 = np.linalg.inv(Bst)
        xBst = np.matmul(Bst1, b.reshape(m, 1)).flatten()
        art = np.where(Bindst >= n)[0]
        if art.size > 0 and not(np.abs(xBst[art]) < 1e-10).all():
            print('不可行问题...', end=',')
            Bind = None
        else:
            A, Bindst, b=prognostic(A, E, Bindst, b)
            Nindst = np.setdiff1d(np.arange(n), Bindst)
            Bind = simplex(A, b, c, Bindst, Nindst)
    else:
        Bind = simplex(A, b, c, Bind1, Nind1)
    if isinstance(Bind, np.ndarray):
        x = np.zeros(n)
        B = A[:, Bind]
        B1 = np.linalg.inv(B)
        xB = np.matmul(B1, b.reshape(b.size, 1)).flatten()
        x = np.zeros(n)
        x[Bind] = xB
        if n > nx:
            x=(x[:nx],x[nx:])
    else:
        x = None
    return x

程序的第2~41行定义的linprog函数接受的参数参数c表示目标函数的系数向量，Aiq，biq表示不等式约束的系数矩阵和常数向量。Aeq，beq表示等式约束的系数矩阵和常数向量。
函数体中，第3~6行的if-else语句记录原问题决策变量数于nx。第7行调用standardizing函数将输入的线性规划问题标准化，返回标准型等式约束的系数矩阵A和常数向量b。第8行调用prepro函数构造原问题标准型的辅助问题，返回值赋予E和Bind1表示添加的人工变量的系数矩阵和辅助问题的初始基矩阵。第14~29行的if-else语句根据是标准型问题否添加了人工变量，决定一次直接解决原问题（第29行的操作）还是通过两个阶段完成计算（第15~27行的操作）。对有人工变量的情形，第16行调用simplex函数，完成第一阶段求解辅助问题的最优解。返回值赋予表示辅助问题最优解对应的基矩阵${\mathbf{B}}^{\ast }$向量下标集的Bindst。第17~20行根据Bindst算得问题辅助问题的最优解中对应人工变量非零值的片段art。第21~27行的内嵌if-else语句根据art是否为空判定标准型原问题是否可行。第22~23行处理不可行情形，第25~27行处理可行情形。其中，第25行调用prognostic函数对第一阶段算得的${\mathbf{B}}^{\ast }$以及标准型的$\mathbf{A}$与$\mathbf{b}$作预后处理。第27行用重构的标准型问题的系数矩阵A，常数向量b即其初始基矩阵列向量下标Bindst，再次调用simplex，完成第二阶段计算。第30~40行的if-else语句根据当前表示标准型问题最终基矩阵是否为空集，计算最终解x。若基矩阵为空集，意味着原问题不可行或无界，第40行将x置为空。否则第31~38行利用Bind构造标准型问题的最优解x。其中，第37~38行的if语句对具有松弛变量的情形，将最优解分解成原问题变量部分（x[:nx]）和松弛变量部分（x[nx:]）。利用linprog函数可求解各种形式的线性规划问题。
综合案例： 一个农民承包了6块耕地，共300亩。准备播种小麦、玉米、蔬菜、瓜果这4种农产品。每个地块由于土质地形不一样，单位面积不同作物的收益是不一样的，如下表所示（本案例数据来自https://opt.aliyun.com/platform/case}）。

	地块1	地块2	地块3	地块4	地块5	地块6
小麦	500	550	630	1000	800	700
玉米	800	700	600	950	90	930
蔬菜	1200	1040	980	860	880	780
瓜果	1000	960	840	650	600	700

6块地的面积分别为42，56，44，39，60和59亩。四种作物分别最多种植76，88，40和96亩。问题是如何安排种植计划，即在每块地中如何安排4种植物种植面积，可以得到最大的总收益。
设第$i$中植物在第$j$块地的种植面积为$s_{ij}$，$i=1,2,3,4$，$j=1,2,3,4,5,6$。根据上表，记
$$\mathbf{P}=(p_{ij}{)}_{4\times 6}=\left(\begin{array}{cccccc}500& 550& 630& 1000& 800& 700\\ 800& 700& 600& 950& 90& 930\\ 1200& 1040& 980& 860& 880& 780\\ 1000& 960& 840& 650& 600& 700\end{array}\right),$$
则总收益为
$$f=\sum _{i=1}^4\sum _{j=1}^6{p}_{ij}{s}_{ij}.$$
受地块大小的限制，有
∑ i = 1 4 s i 1 ≤ 42 , ∑ i = 1 4 s i 2 ≤ 46 , ∑