对标准型线性规划
{
minimize
c
⊤
x
s.t.
A
x
=
b
x
≥
o
(
1
)
\begin{cases} \text{minimize}\quad\quad\boldsymbol{c}^\top\boldsymbol{x}\\ \text{s.t.\ \ \ \ }\quad\quad\quad\boldsymbol{Ax}=\boldsymbol{b}\\ \quad\quad\quad\quad\quad\quad\boldsymbol{x}\geq\boldsymbol{o} \end{cases}\quad\quad(1)
⎩
⎨
⎧minimizec⊤xs.t. Ax=bx≥o(1)
其中
A
∈
R
m
×
n
\boldsymbol{A}\in\text{R}^{m\times n}
A∈Rm×n,且rank
A
=
m
≤
n
\boldsymbol{A}=m\leq n
A=m≤n,
b
≥
o
\boldsymbol{b}\geq\boldsymbol{o}
b≥o。假定
A
\boldsymbol{A}
A中存在部分列向量
α
I
1
,
⋯
,
α
I
m
1
\boldsymbol{\alpha}_{I_1},\cdots,\boldsymbol{\alpha}_{I_{m_1}}
αI1,⋯,αIm1构成单位阵
I
m
\boldsymbol{I}_m
Im的一部分。其中,
m
i
≤
m
m_i\leq m
mi≤m。记
n
1
=
m
−
m
1
n_1=m-m_1
n1=m−m1,
E
\boldsymbol{E}
E为
I
m
\boldsymbol{I}_m
Im中去掉
α
I
1
,
⋯
,
α
I
m
1
\boldsymbol{\alpha}_{I_1},\cdots,\boldsymbol{\alpha}_{I_{m_1}}
αI1,⋯,αIm1剩下的部分,则
E
∈
R
m
×
n
1
\boldsymbol{E}\in\text{R}^{m\times n_1}
E∈Rm×n1。引入人工变量
x
a
∈
R
n
1
=
(
x
a
1
⋮
x
a
n
1
)
\boldsymbol{x}_a\in\text{R}^{n_1}=\begin{pmatrix}x_{a_1}\\\vdots\\x_{a_{n_1}}\end{pmatrix}
xa∈Rn1=
xa1⋮xan1
,令
A
′
=
(
A
,
E
)
\boldsymbol{A}'=(\boldsymbol{A},\boldsymbol{E})
A′=(A,E),构造线性规划
{
minimize
e
⊤
x
a
s.t.
A
′
(
x
x
a
)
=
A
x
+
E
x
a
=
b
x
,
x
a
≥
o
(
2
)
\begin{cases} \text{minimize}\quad\quad\boldsymbol{e}^\top\boldsymbol{x}_a\\ \text{s.t.\ \ \ \ }\quad\quad\quad\boldsymbol{A}'\begin{pmatrix}\boldsymbol{x}\\\boldsymbol{x}_a\end{pmatrix} =\boldsymbol{Ax}+\boldsymbol{Ex}_a =\boldsymbol{b}\\ \quad\quad\quad\quad\quad\quad\boldsymbol{x},\boldsymbol{x}_a\geq\boldsymbol{o} \end{cases}\quad\quad(2)
⎩
⎨
⎧minimizee⊤xas.t. A′(xxa)=Ax+Exa=bx,xa≥o(2)
其中,
e
∈
R
n
\boldsymbol{e}\in\text{R}^n
e∈Rn,所有元素均为1。(2)称为线性规划(1)的辅助线性规划。辅助问题总有一个明显的初始基矩阵:
A
′
=
(
A
,
E
)
\boldsymbol{A}'=(\boldsymbol{A},\boldsymbol{E})
A′=(A,E)的最后
m
m
m列构成的单位阵
I
m
\boldsymbol{I}_m
Im。
例1 标准型线性规划
{
minimize
−
2
x
1
+
x
2
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
≥
0
.
\begin{cases} \text{minimize}\quad\quad -2x_1+x_2\\ \text{s.t.\ \ }\quad\quad\quad\quad x_1+x_2-x_3=2\\ \quad\quad\quad\quad\quad\quad x_1-x_2-x_4=1\\ \quad\quad\quad\quad\quad\quad x_1+x_5=3\\ \quad\quad\quad\quad\quad\quad x_1,x_2,x_3,x_4,x_5\geq0 \end{cases}.
⎩
⎨
⎧minimize−2x1+x2s.t. x1+x2−x3=2x1−x2−x4=1x1+x5=3x1,x2,x3,x4,x5≥0.
其等式约束系数矩阵和常数向量
A
=
(
1
1
−
1
0
0
1
−
1
0
−
1
0
1
0
0
0
1
)
,
b
=
(
2
1
3
)
.
\boldsymbol{A}=\begin{pmatrix}1&1&-1&0&0\\1&-1&0&-1&0\\1&0&0&0&1\end{pmatrix},\boldsymbol{b}=\begin{pmatrix}2\\1\\3\end{pmatrix} .
A=
1111−10−1000−10001
,b=
213
.
其中,
α
5
\boldsymbol{\alpha}_5
α5构成单位阵中一部分。添加人工变量
x
6
,
x
7
x_6,x_7
x6,x7,构造其辅助线性规划
{
minimize
x
6
+
x
7
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
≥
0
.
\begin{cases} \text{minimize}\quad\quad x_6+x_7\\ \text{s.t.\ \ }\quad\quad\quad\quad x_1+x_5=3\\ \quad\quad\quad\quad\quad\quad x_1+x_2-x_3+x_6=2\\ \quad\quad\quad\quad\quad\quad x_1-x_2-x_4+x_7=1\\ \quad\quad\quad\quad\quad\quad x_1,x_2,x_3,x_4,x_5,x_6,x_7\geq0 \end{cases}.
⎩
⎨
⎧minimizex6+x7s.t. x1+x5=3x1+x2−x3+x6=2x1−x2−x4+x7=1x1,x2,x3,x4,x5,x6,x7≥0.
其等式约束系数矩阵为
A
′
=
(
A
,
E
)
=
(
1
1
−
1
0
0
1
0
1
−
1
0
−
1
0
0
1
1
0
0
0
1
0
0
)
.
\boldsymbol{A}'=(\boldsymbol{A},\boldsymbol{E})=\begin{pmatrix}1&1&-1&0&0&1&0\\1&-1&0&-1&0&0&1\\1&0&0&0&1&0&0\end{pmatrix}.
A′=(A,E)=
1111−10−1000−10001100010
.
其中,
E
=
(
α
6
,
α
7
)
=
(
1
0
0
1
0
0
)
\boldsymbol{E}=(\boldsymbol{\alpha}_6,\boldsymbol{\alpha}_7)=\begin{pmatrix}1&0\\0&1\\0&0\end{pmatrix}
E=(α6,α7)=
100010
。显然,
(
α
6
,
α
7
,
α
5
)
=
I
3
(\boldsymbol{\alpha}_6,\boldsymbol{\alpha}_7,\boldsymbol{\alpha}_5)=\boldsymbol{I}_3
(α6,α7,α5)=I3成为辅助线性规划的一个基矩阵。
下列代码定义构造标准型线性规划(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表示标准型问题的系数矩阵
A
\boldsymbol{A}
A。函数体内第12、13行分别将表示人工变量系数矩阵的E和单位向量在
(
A
,
E
)
(\boldsymbol{A},\boldsymbol{E})
(A,E)中位置下标的序列的pos初始化为空集和
m
m
m个
−
1
-1
−1的数组。第14行定义单位阵e,其中的每一个行向量e[i]构成一个单位向量。第17~25行的while循环确定每个单位向量e[i]在
(
A
,
E
)
(\boldsymbol{A},\boldsymbol{E})
(A,E)中位置(或在A中,或追加于E),记录于pos。其中,第18行调用的find_e函数在矩阵A中查找与单位向量e[i]。该函数是定义于第3~10行内部函数。若e[i]为A中的第j列,返回j,否则返回
−
1
-1
−1。本函数最终返回E和pos。
例2 用prepro构造
{
minimize
−
2
x
1
+
x
2
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
≥
0
.
\begin{cases} \text{minimize}\quad\quad -2x_1+x_2\\ \text{s.t.\ \ }\quad\quad\quad\quad x_1+x_2-x_3=2\\ \quad\quad\quad\quad\quad\quad x_1-x_2-x_4=1\\ \quad\quad\quad\quad\quad\quad x_1+x_5=3\\ \quad\quad\quad\quad\quad\quad x_1,x_2,x_3,x_4,x_5\geq0 \end{cases}.
⎩
⎨
⎧minimize−2x1+x2s.t. x1+x2−x3=2x1−x2−x4=1x1+x5=3x1,x2,x3,x4,x5≥0.
的辅助问题的等式约束系数矩阵。
解: 根据例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]]
意为辅助问题中人工变量的系数矩阵
E
=
(
1
0
0
1
0
0
)
\boldsymbol{E}=\begin{pmatrix}1&0\\0&1\\0&0\end{pmatrix}
E=
100010
,初始基矩阵列向量为
(
α
6
,
α
7
,
α
5
)
(\boldsymbol{\alpha}_6,\boldsymbol{\alpha}_7,\boldsymbol{\alpha}_5)
(α6,α7,α5)。辅助问题等式约束矩阵
A
′
=
(
A
,
E
)
=
(
1
1
−
1
0
0
1
0
1
−
1
0
−
1
0
0
1
1
0
0
0
1
0
0
)
\boldsymbol{A}'=(\boldsymbol{A},\boldsymbol{E})=\begin{pmatrix}1&1&-1&0&0&1&0\\1&-1&0&-1&0&0&1\\1&0&0&0&1&0&0\end{pmatrix}
A′=(A,E)=
1111−10−1000−10001100010
。与例1的计算结果一致。
辅助问题(2)由于其目标函数
e
⊤
x
a
≥
0
\boldsymbol{e}^\top\boldsymbol{x}_a\geq0
e⊤xa≥0有下界,故必有最优解。辅助问题(2)与原问题(1)有如下的关系
定理1 标准型线性规划(1)有可行解,当且仅当其辅助线性规划(2)的最优解为
x
a
∗
=
o
\boldsymbol{x}_a^{*}=\boldsymbol{o}
xa∗=o。
对问题(1)运用单纯形算法可算得其辅助问题(2)最优解
x
∗
=
(
x
x
a
)
\boldsymbol{x}^{*}=\begin{pmatrix}\boldsymbol{x}\\\boldsymbol{x}_a\end{pmatrix}
x∗=(xxa),求解过程称为第一阶段。按定理1,若
x
a
=
o
\boldsymbol{x}_a=\boldsymbol{o}
xa=o则原问题(1)可行。并且,若
x
∗
\boldsymbol{x}^{*}
x∗中的原问题决策变量部分
x
\boldsymbol{x}
x恰巧是是(1)的基本可行解,即第一阶段计算所得的辅助问题的基矩阵
B
∗
\boldsymbol{B}^{*}
B∗恰包含在原问题(1)的系数矩阵
A
\boldsymbol{A}
A中,则以此为起点,可运用单纯形算法求解问题(1),这一过程称为第二阶段。
例3 标准型线性规划
{
minimize
−
2
x
1
+
x
2
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
≥
0
.
\begin{cases} \text{minimize}\quad\quad -2x_1+x_2\\ \text{s.t.\ \ }\quad\quad\quad\quad x_1+x_2-x_3=2\\ \quad\quad\quad\quad\quad\quad x_1-x_2-x_4=1\\ \quad\quad\quad\quad\quad\quad x_1+x_5=3\\ \quad\quad\quad\quad\quad\quad x_1,x_2,x_3,x_4,x_5\geq0 \end{cases}.
⎩
⎨
⎧minimize−2x1+x2s.t. x1+x2−x3=2x1−x2−x4=1x1+x5=3x1,x2,x3,x4,x5≥0.
其目标函数系数向量
c
=
(
−
2
1
0
0
0
)
\boldsymbol{c}=\begin{pmatrix}-2\\1\\0\\0\\0\end{pmatrix}
c=
−21000
,等式约束系数矩阵和常数向量
A
=
(
1
1
−
1
0
0
1
−
1
0
−
1
0
1
0
0
0
1
)
,
b
=
(
2
1
3
)
.
\boldsymbol{A}=\begin{pmatrix}1&1&-1&0&0\\1&-1&0&-1&0\\1&0&0&0&1\end{pmatrix},\boldsymbol{b}=\begin{pmatrix}2\\1\\3\end{pmatrix} .
A=
1111−10−1000−10001
,b=
213
.
下列代码求解该问题。
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
x_6=x_7=0
x6=x7=0,故原问题可行。辅助问题最优解对应基矩阵列向量下标为2,1,5均含于原问题系数矩阵中,故可利用它直接求解原问题。
写博不易,敬请支持:
如果阅读本文于您有所获,敬请点赞、评论、收藏,谢谢大家的支持!