最优化方法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}.$$
受地块大小的限制，有
$$\begin{array}{l}\sum _{i=1}^4{s}_{i1}\le 42,\\ \sum _{i=1}^4{s}_{i2}\le 46,\\ \sum _{i=1}^4{s}_{i3}\le 44,\\ \sum _{i=1}^4{s}_{i4}\le 39,\\ \sum _{i=1}^4{s}_{i5}\le 60,\\ \sum _{i=1}^4{s}_{i6}\le 59.\end{array}$$
又由于每种作物的最大种植面积的限制，有
$$\begin{array}{l}\sum _{j=1}^6{s}_{1j}\le 76,\\ \sum _{j=1}^6{s}_{2j}\le 88,\\ \sum _{j=1}^6{s}_{3j}\le 40,\\ \sum _{j=1}^6{s}_{4j}\le 96.\end{array}$$
为便于将问题表示成线性规划模型，令$k=6(i-1)+j$（$i=1,2,3,4$，$j=1,2,3,4,5,6$），则$s_{ij}=x_k$，$p_{ij}=c_k$，$k=1,2,\cdots ,24$。于是，本问题模型化为线性规划：
$$\{\begin{array}{l}\text{maximize}\sum _{k=1}^{24}{c}_k{x}_k\\ \text{s.t. }\sum _{i=1}^4{x}_{6(i-1)+1}\le 42,\\ \sum _{i=1}^4{x}_{6(i-1)+2}\le 46,\\ \sum _{i=1}^4{x}_{6(i-1)+3}\le 44,\\ \sum _{i=1}^4{x}_{6(i-1)+4}\le 39,\\ \sum _{i=1}^4{x}_{6(i-1)+5}\le 60,\\ \sum _{i=1}^4{x}_{6(i-1)+6}\le 59,\\ \sum _{j=1}^6{x}_j\le 76,\\ \sum _{j=1}^6{x}_{6+j}\le 88,\\ \sum _{j=1}^6{x}_{12+j}\le 40,\\ \sum _{j=1}^6{x}_{18+j}\le 96.\end{array}$$
记$\mathbf{c}=\left(\begin{array}{c}{c}_1\\ ⋮\\ c_{24}\end{array}\right)$，$\mathbf{x}=\left(\begin{array}{c}{x}_1\\ ⋮\\ x_{24}\end{array}\right)$，
$${\mathbf{A}}_{iq}=\left(\begin{array}{cccccccccccccccccccccccc}1& 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 1\\ 1& 1& 1& 1& 1& 1& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 1& 1& 1& 1& 1& 1& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 1& 1& 1& 1& 1& 1& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 1& 1& 1& 1& 1& 1\end{array}\right),$$
$${\mathbf{b}}_{iq}^{\top }=(42,46,44,39,60,59,76,88,40,96).$$
则线性规划的矩阵形式为
$$\{\begin{array}{l}\text{minmize}-{\mathbf{c}}^{\top }\mathbf{x}\\ \text{s.t.}{\mathbf{A}}_{iq}\mathbf{x}\le {\mathbf{b}}_{iq}\\ \mathbf{x}\ge \mathbf{o}\end{array}.$$
下列代码完成本案例的求解计算。

import numpy as np										#导入numpy
from fractions import Fraction as F						#设置输出格式
np.set_printoptions(formatter = {'all':lambda x:
                               str(F(x).limit_denominator())})
P = np.array([[500, 550, 630, 1000, 800, 700],			#单位面积收益矩阵
              [800, 700, 600, 950, 90, 930],
              [1200, 1040, 980, 860, 880, 280],
              [1000, 960, 840, 650, 600, 700]])
c = -P.flatten()										#目标函数系数向量
Up = np.eye(6)
Lo = np.zeros((4, 6))
Lo[0] = np.ones(6)
for i in range(3):										#构造系数矩阵组成块
    Up = np.hstack((Up, np.eye(6)))
    T = np.zeros((4, 6))
    T[i + 1] = np.ones(6)
    Lo = np.hstack((Lo,T))
Ai = np.vstack((Up, Lo))								#不等式约束系数矩阵
bi = np.array([42, 46, 44, 39, 60, 59, 76, 88, 40, 96])	#常数向量
linprog(c, Ai, bi)

程序的第5~8行将表格所列数据构造矩阵P。第9行按行优先顺序将P扁平化为并取负值赋予表示目标函数系数向量的c。注意到问题的不等式约束矩阵${\mathbf{A}}_{iq}$的结构：前6行由4个6阶单位阵并列而成。后4行由4个$4\times 6$的矩阵并列而成，每一个均有一行全为1，其余元素为0。元素为1的行位于第$i$行，$i=1,2,3,4$。按此观察，程序的第10行将Up初始化为6阶单位阵，第11~12行将Lo初始化为第1行元素为1，其余元素全为0的$4\times 6$矩阵。第13~17行的for循环，将Up和Lo分别组成系数矩阵${\mathbf{A}}_{iq}$的前6行和后4行。第18行将Up和Lo纵向叠加为系数矩阵，赋予Ai。第19行设置不等式约束的常数向量${\mathbf{b}}_{iq}$为bi。运行程序，输出

(array([0, 0, 0, 16, 60, 0, 0, 0, 0, 23, 0, 59,
        40, 0, 0, 0, 0, 0, 2, 46, 44, 0, 0, 0]),
 array([0, 0, 0, 0, 0, 0, 0, 6, 0, 4]))

意为问题的最优解为${\mathbf{x}}_0^{\top }=(0,0,0,16,60,0,0,0,0,23,0,59,40,0,0,0,0,0,2,46,44,0,0,0)$，松弛变量对应值为$(x_{25},x_{26},x_{27},x_{28},x_{29},x_{30},x_{31},x_{32},x_{33},x_{34})=(0,0,0,0,0,0,0,6,0,4)$。即小麦在4号地块种植16亩，第5号地块种植60亩；玉米在4号地块种23亩，6号地块种59亩；蔬菜在1号地块种40亩；瓜果在1号地块种2亩，2号地块种46亩，3号地块种44亩可使收益最大。松弛变量$x_{32}=6$意味着玉米实际种植面积距最大计划种植面积少6亩。相仿地，瓜果实际种植面积比最大计划种植面积少$x_{34}=4$亩。
写博不易，敬请支持：
如果阅读本文于您有所获，敬请点赞、评论、收藏，谢谢大家的支持！