最优化方法Python计算：标准型线性规划的辅助问题

对标准型线性规划
$$\{\begin{array}{l}\text{minimize}{\mathbf{c}}^{\top }\mathbf{x}\\ \text{s.t. }\mathbf{Ax}=\mathbf{b}\\ \mathbf{x}\ge \mathbf{o}\end{array}(1)$$
其中$\mathbf{A}\in {\text{R}}^{m\times n}$，且rank$\mathbf{A}=m\le n$，$\mathbf{b}\ge \mathbf{o}$。假定$\mathbf{A}$中存在部分列向量${\mathbf{\alpha }}_{I_1},\cdots ,{\mathbf{\alpha }}_{I_{m_1}}$构成单位阵${\mathbf{I}}_m$的一部分。其中，$m_i\le m$。记$n_1=m-m_1$，$\mathbf{E}$为${\mathbf{I}}_m$中去掉${\mathbf{\alpha }}_{I_1},\cdots ,{\mathbf{\alpha }}_{I_{m_1}}$剩下的部分，则$\mathbf{E}\in {\text{R}}^{m\times n_1}$。引入人工变量${\mathbf{x}}_a\in {\text{R}}^{n_1}=\left(\begin{array}{c}{x}_{a_1}\\ ⋮\\ x_{a_{n_1}}\end{array}\right)$，令${\mathbf{A}}^{\prime }=(\mathbf{A},\mathbf{E})$，构造线性规划
$$\{\begin{array}{l}\text{minimize}{\mathbf{e}}^{\top }{\mathbf{x}}_a\\ \text{s.t. }{\mathbf{A}}^{\prime }\left(\begin{array}{c}\mathbf{x}\\ {\mathbf{x}}_a\end{array}\right)=\mathbf{Ax}+{\mathbf{Ex}}_a=\mathbf{b}\\ \mathbf{x},{\mathbf{x}}_a\ge \mathbf{o}\end{array}(2)$$
其中，$\mathbf{e}\in {\text{R}}^n$，所有元素均为1。(2)称为线性规划(1)的辅助线性规划。辅助问题总有一个明显的初始基矩阵：${\mathbf{A}}^{\prime }=(\mathbf{A},\mathbf{E})$的最后$m$列构成的单位阵${\mathbf{I}}_m$。
例1 标准型线性规划
$$\{\begin{array}{l}\text{minimize}-2x_1+x_2\\ \text{s.t. }{x}_1+x_2-x_3=2\\ x_1-x_2-x_4=1\\ x_1+x_5=3\\ x_1,x_2,x_3,x_4,x_5\ge 0\end{array}.$$
其等式约束系数矩阵和常数向量
$$\mathbf{A}=\left(\begin{array}{ccccc}1& 1& -1& 0& 0\\ 1& -1& 0& -1& 0\\ 1& 0& 0& 0& 1\end{array}\right),\mathbf{b}=\left(\begin{array}{c}2\\ 1\\ 3\end{array}\right).$$
其中，${\mathbf{\alpha }}_5$构成单位阵中一部分。添加人工变量$x_6,x_7$，构造其辅助线性规划
$$\{\begin{array}{l}\text{minimize}{x}_6+x_7\\ \text{s.t. }{x}_1+x_5=3\\ x_1+x_2-x_3+x_6=2\\ x_1-x_2-x_4+x_7=1\\ x_1,x_2,x_3,x_4,x_5,x_6,x_7\ge 0\end{array}.$$
其等式约束系数矩阵为
$${\mathbf{A}}^{\prime }=(\mathbf{A},\mathbf{E})=\left(\begin{array}{ccccccc}1& 1& -1& 0& 0& 1& 0\\ 1& -1& 0& -1& 0& 0& 1\\ 1& 0& 0& 0& 1& 0& 0\end{array}\right).$$
其中，$\mathbf{E}=({\mathbf{\alpha }}_6,{\mathbf{\alpha }}_7)=\left(\begin{array}{cc}1& 0\\ 0& 1\\ 0& 0\end{array}\right)$。显然，$({\mathbf{\alpha }}_6,{\mathbf{\alpha }}_7,{\mathbf{\alpha }}_5)={\mathbf{I}}_3$成为辅助线性规划的一个基矩阵。
下列代码定义构造标准型线性规划（1）的辅助问题（2）的Python函数。

import numpy as np									#导入numpy
def prepro(A):
    def find_e(A, ei):								#查找A中单位向量位置
        n = A.shape[1]
        j = 0
        while j < n:
            if(abs(A.T[j] - ei) < 1e-10).all():
                return j
            j += 1
        return -1
    m, n = A.shape
    pos = np.array([-1] * m)						#单位向量位置下标序列
    E = np.array([[]for i in range(m)])				#人工变量系数矩阵
    e = np.eye(m)									#单位阵
    i = 0
    k = 0
    while i < m:									#对每一个单位向量
        j = find_e(A, e[i])							#查找在A中位置
        if j >= 0:									#若在A中出现
            pos[i] = j								#记录位置下标
        else:										#未在A中
            E = np.hstack((E, e[i].reshape(m, 1)))	#追加到E
            pos[i] = n + k							#记录位置
            k += 1
        i += 1
    return E, pos

程序的第2~26行定义的prepro函数实现构造标准型线性规划辅助问题的过程。其唯一的参数A表示标准型问题的系数矩阵$\mathbf{A}$。函数体内第12、13行分别将表示人工变量系数矩阵的E和单位向量在$(\mathbf{A},\mathbf{E})$中位置下标的序列的pos初始化为空集和$m$个$-1$的数组。第14行定义单位阵e，其中的每一个行向量e[i]构成一个单位向量。第17~25行的while循环确定每个单位向量e[i]在$(\mathbf{A},\mathbf{E})$中位置（或在A中，或追加于E），记录于pos。其中，第18行调用的find_e函数在矩阵A中查找与单位向量e[i]。该函数是定义于第3~10行内部函数。若e[i]为A中的第j列，返回j，否则返回$-1$。本函数最终返回E和pos。
例2 用prepro构造
$$\{\begin{array}{l}\text{minimize}-2x_1+x_2\\ \text{s.t. }{x}_1+x_2-x_3=2\\ x_1-x_2-x_4=1\\ x_1+x_5=3\\ x_1,x_2,x_3,x_4,x_5\ge 0\end{array}.$$
的辅助问题的等式约束系数矩阵。
解： 根据例1的数据，下列程序完成计算

import numpy as np#导入numpy
from fractions import Fraction as F
np.set_printoptions(formatter={'all':lambda x:#设置输出格式
                               str(F(x).limit_denominator())})
A = np.array([[1, 1, -1, 0, 0],	#原问题系数矩阵
              [1, -1, 0, -1, 0],
              [1, 0, 0, 0, 1]])
E, pos = prepro(A)				#引入人工变量
print(E, pos)
A1 = np.hstack((A,E))
print(A1)

程序的第2~4行设置数据输出格式。第5~7行设置原问题等式约束系数矩阵A。第8行调用prepro函数，传递A计算辅助系统中人工变量的系数矩阵E和辅助问题的初始基矩阵列向量下标集pos。第10行构造辅助问题等式约束系数矩阵。运行程序，输出

[[1 0]
 [0 1]
 [0 0]] [5 6 4]
[[1 1 -1 0 0 1 0]
 [1 -1 0 -1 0 0 1]
 [1 0 0 0 1 0 0]]

意为辅助问题中人工变量的系数矩阵$\mathbf{E}=\left(\begin{array}{cc}1& 0\\ 0& 1\\ 0& 0\end{array}\right)$，初始基矩阵列向量为$({\mathbf{\alpha }}_6,{\mathbf{\alpha }}_7,{\mathbf{\alpha }}_5)$。辅助问题等式约束矩阵${\mathbf{A}}^{\prime }=(\mathbf{A},\mathbf{E})=\left(\begin{array}{ccccccc}1& 1& -1& 0& 0& 1& 0\\ 1& -1& 0& -1& 0& 0& 1\\ 1& 0& 0& 0& 1& 0& 0\end{array}\right)$。与例1的计算结果一致。
辅助问题（2）由于其目标函数${\mathbf{e}}^{\top }{\mathbf{x}}_a\ge 0$有下界，故必有最优解。辅助问题（2）与原问题（1）有如下的关系
定理1 标准型线性规划(1)有可行解，当且仅当其辅助线性规划(2)的最优解为${\mathbf{x}}_a^{\ast }=\mathbf{o}$。
对问题（1）运用单纯形算法可算得其辅助问题(2)最优解${\mathbf{x}}^{\ast }=\left(\begin{array}{c}\mathbf{x}\\ {\mathbf{x}}_a\end{array}\right)$，求解过程称为第一阶段。按定理1，若${\mathbf{x}}_a=\mathbf{o}$则原问题(1)可行。并且，若${\mathbf{x}}^{\ast }$中的原问题决策变量部分$\mathbf{x}$恰巧是是(1)的基本可行解，即第一阶段计算所得的辅助问题的基矩阵${\mathbf{B}}^{\ast }$恰包含在原问题（1）的系数矩阵$\mathbf{A}$中，则以此为起点，可运用单纯形算法求解问题(1)，这一过程称为第二阶段。
例3 标准型线性规划
$$\{\begin{array}{l}\text{minimize}-2x_1+x_2\\ \text{s.t. }{x}_1+x_2-x_3=2\\ x_1-x_2-x_4=1\\ x_1+x_5=3\\ x_1,x_2,x_3,x_4,x_5\ge 0\end{array}.$$
其目标函数系数向量$\mathbf{c}=\left(\begin{array}{c}-2\\ 1\\ 0\\ 0\\ 0\end{array}\right)$，等式约束系数矩阵和常数向量
$$\mathbf{A}=\left(\begin{array}{ccccc}1& 1& -1& 0& 0\\ 1& -1& 0& -1& 0\\ 1& 0& 0& 0& 1\end{array}\right),\mathbf{b}=\left(\begin{array}{c}2\\ 1\\ 3\end{array}\right).$$
下列代码求解该问题。

import numpy as np													#导入numpy
from fractions import Fraction as F									#设置输出格式
np.set_printoptions(formatter={'all':lambda x:
                               str(F(x).limit_denominator())})
A = np.array([[1, 1, -1, 0, 0],										#原问题数据设置
              [1, -1, 0, -1, 0],
              [1, 0, 0, 0, 1]])
b = np.array([2, 1, 3])
c = np.array([-1, 2, 0, 0, 0])
E, pos = prepro(A)													#构造辅助问题系数矩阵
A1 = np.hstack((A, E))
Bind = pos															#构造辅助问题初始基矩阵
Nind = np.setdiff1d(np.arange(5+E.shape[1]),Bind)
c1 = np.concatenate((np.zeros(5), np.ones(E.shape[1])), axis = 0)	#辅助问题目标函数系数
Bindst = simplex(A1, b, c1, Bind, Nind)								#求解辅助问题最优解
Nindst = np.setdiff1d(np.arange(5), Bindst)
Bst = A1[:, Bindst]													#辅助问题最优解基矩阵
Bst1 = np.linalg.inv(B)												#逆阵
xBst = np.matmul(Bst1,b.reshape(3, 1)).flatten()					#对应基矩阵最优解部分
x = np.zeros(5 + E.shape[1])
x[Bindst] = xBst													#辅助问题最优解
print('辅助问题最优解x=%s'%x)
print('辅助问题最优解对应的基矩阵向量下标：%s'%Bindst)					#辅助问题最优解基矩阵包含于原问题系数矩阵
Bind = simplex(A, b, c, Bindst, Nindst)								#求解原问题最优解
B = A[:, Bind]														#当前基矩阵
B1 = np.linalg.inv(B)												#逆阵
xB = np.matmul(B1,b.reshape(3, 1)).flatten()						#对应基矩阵最优解部分
x = np.zeros(5)
x[Bind] = xB														#原问题最优解
print('原问题最优解x=%s'%x)

程序的前14行与例2中的代码一样调用prepro构造原问题辅助问题。第15行调用博文《最优化方法Python计算：标准型线性规划的单纯形算法》定义的simplex函数求解辅助问题，完成第一阶段计算。第17~22行计算并输出辅助问题最优解。根据第15所得辅助问题最优解对应的基矩阵（第23行输出）包含于原问题的系数矩阵中（见下列输出），第18行调用simplex求解原问题，完成第二阶段计算。运行程序，输出

辅助问题最优解x=[3/2 1/2 0 0 3/2 0 0]
辅助问题最优解对应的基矩阵向量下标：[1 0 4]
原问题最优解x=[3 0 1 2 0]

注意辅助问题最优解中，人工变量$x_6=x_7=0$，故原问题可行。辅助问题最优解对应基矩阵列向量下标为2，1，5均含于原问题系数矩阵中，故可利用它直接求解原问题。
写博不易，敬请支持：
如果阅读本文于您有所获，敬请点赞、评论、收藏，谢谢大家的支持！