最优化方法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{c}$、$\mathbf{x}\in {\text{R}}^n$，$\mathbf{A}\in {\text{R}}^{m\times n}$，$m\le n$且rank$\mathbf{A}=m$，$\mathbf{b}\in {\text{R}}^m$且$\mathbf{b}\ge \mathbf{o}$。设$\mathbf{A}=({\mathbf{\alpha }}_1,{\mathbf{\alpha }}_2\cdots ,{\mathbf{\alpha }}_n)$。$\mathbf{B}$由$\mathbf{A}$中$m$个线性无关列向量构成，满足${\mathbf{B}}^{-1}\mathbf{b}\ge \mathbf{b}$，称为标准型线性规划(1)的一个基矩阵，其列向量下标集记为$B_{ind}$。$\mathbf{A}$中去掉$\mathbf{B}$剩下部分记为$\mathbf{N}$，称为(1)的对应$\mathbf{B}$的非基矩阵，其列向量下标集记为$N_{ind}$。于是线性方程组$\mathbf{Ax}=\mathbf{b}$可写为$(\mathbf{B},\mathbf{N})\left(\begin{array}{c}{\mathbf{x}}_{\mathbf{B}}\\ {\mathbf{x}}_{\mathbf{N}}\end{array}\right)=\mathbf{b}$，即
$$\mathbf{B}{\mathbf{x}}_{\mathbf{B}}+\mathbf{N}{\mathbf{x}}_{\mathbf{N}}=\mathbf{b},$$
亦即
$${\mathbf{x}}_{\mathbf{B}}={\mathbf{B}}^{-1}\mathbf{b}-{\mathbf{B}}^{-1}\mathbf{N}{\mathbf{x}}_{\mathbf{N}}.$$
令自由未知量${\mathbf{x}}_{\mathbf{N}}=\mathbf{o}$得到$\mathbf{Ax}=\mathbf{b}$的解${\mathbf{x}}_0=\left(\begin{array}{c}{\mathbf{x}}_{\mathbf{B}}\\ {\mathbf{x}}_{\mathbf{N}}\end{array}\right)=\left(\begin{array}{c}{\mathbf{B}}^{-1}\mathbf{b}\\ \mathbf{o}\end{array}\right)$称为问题(1)的一个初始基本可行解。对初始基本可行解，对应的目标函数值$f({\mathbf{x}}_0)={\mathbf{c}}^{\top }{\mathbf{x}}_0={\mathbf{c}}_{\mathbf{B}}^{\top }{\mathbf{x}}_{\mathbf{B}}+{\mathbf{c}}_{\mathbf{N}}^{\top }{\mathbf{x}}_{\mathbf{N}}={\mathbf{c}}_{\mathbf{B}}^{\top }{\mathbf{x}}_{\mathbf{B}}={\mathbf{c}}_{\mathbf{B}}^{\top }{\mathbf{B}}^{-1}\mathbf{b}$。对任意的${\mathbf{x}}_{\mathbf{N}}\ge \mathbf{o}$，$\mathbf{x}=\left(\begin{array}{c}{\mathbf{x}}_{\mathbf{B}}\\ {\mathbf{x}}_{\mathbf{N}}\end{array}\right)$，考虑目标函数
$$\begin{array}{rl}f(\mathbf{x})& ={\mathbf{c}}^{\top }\mathbf{x}={\mathbf{c}}_{\mathbf{B}}^{\top }{\mathbf{x}}_{\mathbf{B}}+{\mathbf{c}}_{\mathbf{N}}^{\top }{\mathbf{x}}_{\mathbf{N}}\\ & ={\mathbf{c}}_{\mathbf{B}}^{\top }({\mathbf{B}}^{-1}\mathbf{b}-{\mathbf{B}}^{-1}\mathbf{N}{\mathbf{x}}_{\mathbf{N}})+{\mathbf{c}}_{\mathbf{N}}^{\top }{\mathbf{x}}_{\mathbf{N}}\\ & ={\mathbf{c}}_{\mathbf{B}}^{\top }{\mathbf{B}}^{-1}\mathbf{b}-({\mathbf{c}}_{\mathbf{B}}^{\top }{\mathbf{B}}^{-1}\mathbf{N}{\mathbf{x}}_{\mathbf{N}}-{\mathbf{c}}_{\mathbf{N}}^{\top }{\mathbf{x}}_{\mathbf{N}})\\ & =f({\mathbf{x}}_0)-({\mathbf{c}}_{\mathbf{B}}^{\top }{\mathbf{B}}^{-1}\mathbf{N}-{\mathbf{c}}_{\mathbf{N}}^{\top }){\mathbf{x}}_{\mathbf{N}}.\end{array}$$
令$\mathbf{z}={\mathbf{c}}_{\mathbf{B}}^{\top }{\mathbf{B}}^{-1}\mathbf{N}$，考虑向量$\mathbf{z}-{\mathbf{c}}_{\mathbf{N}}^{\top }$。只要存在$k\in N_{ind}$，使得$z_k-c_k>0$，选取其中下标最小者
e = min ⁡ k ∈ N i n d { k ∣ z k − c k > 0 } e=\min_{k\in N_{ind}}\{k|z_k-c_k>0\} e=k∈Nind​min​{k∣zk​−ck​>0}
解等式约束条件形成的线性方程组，并令自由未知量${\mathbf{x}}_{\mathbf{N}}=\left(\begin{array}{c}0\\ ⋮\\ x_e\\ ⋮\\ 0\end{array}\right)$，即
$$x_k=\{\begin{array}{ll}{x}_e& k\in N_{ind},k=e\\ 0& k\in N_{ind},k\ne e\end{array}$$
其中$x_e\ge 0$。得
$${\mathbf{x}}_{\mathbf{B}}={\mathbf{B}}^{-1}\mathbf{b}-{\mathbf{B}}^{-1}\mathbf{N}{\mathbf{x}}_{\mathbf{N}}={\mathbf{B}}^{-1}\mathbf{b}-{\mathbf{B}}^{-1}{\mathbf{\alpha }}_e{x}_e=\bar{\mathbf{b}}-x_e{\mathbf{y}}_e$$
其中，$\bar{\mathbf{b}}={\mathbf{B}}^{-1}\mathbf{b}=\left(\begin{array}{c}{\bar{b}}_1\\ {\bar{b}}_2\\ ⋮\\ {\bar{b}}_m\end{array}\right)$，${\mathbf{y}}_e={\mathbf{B}}^{-1}{\mathbf{\alpha }}_e=\left(\begin{array}{c}{y}_{e_1}\\ y_{e_2}\\ ⋮\\ y_{e_m}\end{array}\right)$。若${\mathbf{y}}_e\le \mathbf{o}$，则问题（1）为无界问题。设存在$i\in Bind$，使得$y_{e_i}>0$，设
$$l=\min \left\{j|\frac{{\bar{b}}_j}{y_{e_j}}=\underset{i\in B_{ind}}{\min }\left\{\frac{{\bar{b}}_i}{y_{e_i}}|y_{e_i}>0\right\}\right\}$$
即取$l$为 min ⁡ i ∈ B i n d { b ˉ i y e i ∣ y e i > 0 } \min\limits_{i\in B_{ind}}\{\frac{\bar{b}_i}{y_{e_i}}|y_{e_i}>0\} i∈Bind​min​{yei​​bˉi​​∣yei​​>0}中的最小下标（ min ⁡ i ∈ B i n d { b ˉ i y e i ∣ y e i > 0 } \min\limits_{i\in B_{ind}}\{\frac{\bar{b}_i}{y_{e_i}}|y_{e_i}>0\} i∈Bind​min​{yei​​bˉi​​∣yei​​>0}可能有多于1个的等值最小者）。将$\mathbf{B}$中的${\mathbf{\alpha }}_l$与$\mathbf{N}$中的${\mathbf{\alpha }}_e$进行交换，即
{ B = ( B − { α l } ) ∪ { α e } N = ( N − { α e } ) ∪ { α l } \begin{cases} \boldsymbol{B}=(\boldsymbol{B}-\{\boldsymbol{\alpha}_l\})\cup\{\boldsymbol{\alpha}_e\}\\ \boldsymbol{N}=(\boldsymbol{N}-\{\boldsymbol{\alpha}_e\})\cup\{\boldsymbol{\alpha}_l\} \end{cases} {B=(B−{αl​})∪{αe​}N=(N−{αe​})∪{αl​}​
称为问题（1）的一次轴转操作。可以证明，轴转后$\mathbf{B}$仍然是问题（1）的一个基矩阵，且对应的基础可行解${\mathbf{x}}_1$的目标函数值$f({\mathbf{x}}_1)<f({\mathbf{x}}_0)$。轴转操作中入基向量下标$e$和离基限量下标$l$的选取规则
{ e = min ⁡ k ∈ N i n d { k ∣ z k − c k > 0 } l = min ⁡ { j ∣ b ˉ j y e j = min ⁡ i ∈ B i n d { b ˉ i y e i ∣ y e i > 0 } } \begin{cases} e=\min\limits_{k\in N_{ind}}\{k|z_k-c_k>0\}\\ l=\min\left\{j|\frac{\bar{b}_j}{y_{e_j}}=\min\limits_{i\in B_{ind}}\left\{\frac{\bar{b}_i}{y_{e_i}}|y_{e_i}>0\right\}\right\} \end{cases} ⎩ ⎨ ⎧​e=k∈Nind​min​{k∣zk​−ck​>0}l=min{j∣yej​​bˉj​​=i∈Bind​min​{yei​​bˉi​​∣yei​​>0}}​
称为Bland法则。下列代码对问题（1）的初始基矩阵$\mathbf{B}$实现轴转过程。

import numpy as np
def pivot(A, b, c, Bind, Nind, z):
    m, _ = A.shape										#A的行数
    B = A[:, Bind]										#基矩阵
    B1 = np.linalg.inv(B)								#基矩阵的逆
    e = np.min(np.where((z - c[Nind]) > 1e-10)[0])		#计算入基向量下标
    ye = np.matmul(B1, A[:, Nind[e]]).flatten()
    infinite = False									#初始化无界标志
    index = np.where(ye > 1e-10)[0]						#ye中正值元素下标
    if index.size = = 0:								#若ye所有元素不超过0
        infinite = True
        return infinite, None, None, None
    b1 = np.matmul(B1, b.reshape(m, 1)).flatten()
    minl = np.min(b1[index] / ye[index])				#计算离基向量下标
    l = min(index[np.where(b1[index] / ye[index] == minl)[0]])
    Bind[l], Nind[e] = Nind[e], Bind[l]					#轴转
    B = A[:, Bind]										#基矩阵
    B1 = np.linalg.inv(B)								#基矩阵的逆
    w = np.matmul(c[Bind], B1)
    z = np.matmul(w, A[:, Nind]).flatten()				#构造判别式中的z
    return infinite, Bind, Nind, z

函数pivot的参数A，b，c分别表示待解标准型线性规划的等式约束系数矩阵$\mathbf{A}$，常数向量$\mathbf{b}$和目标函数系数向量$\mathbf{c}$。参数Bind和Nind分别表示初始基矩阵$\mathbf{B}$和非基矩阵$\mathbf{N}$的列向量在$\mathbf{A}$中的下标数组。
函数体中第4行，用Bind作为列标，读取A中列向量组成基矩阵$\mathbf{B}$，赋予B。第5行调用Numpy.linalg包的inv函数，计算${\mathbf{B}}^{-1}$赋予B1。第6行按Bland规则计算入基向量下标e。第7行利用e计算向量${\mathbf{y}}_e={\mathbf{B}}^{-1}{\mathbf{\alpha }}_e$赋予ye。第9行用Numpy的where函数计算ye中大于零的元素下标，赋予数组index。第10~12行的if语句，根据ye中元素是否全非正，设置无界标志infinite为True（infinite在第8行初始化为False），并返回infinite及3个空集。第13行计算向量$\bar{\mathbf{b}}={\mathbf{B}}^{-1}\mathbf{b}$赋予b1。第14~15行按Bland规则计算立即向量下标赋予l。第16通过交换Bind和Nind中元素的方式行执行轴转操作
{ B = ( B − { α l } ) ∪ { α e } N = ( N − { α e } ) ∪ { α l } . \begin{cases} \boldsymbol{B}=(\boldsymbol{B}-\{\boldsymbol{\alpha}_l\})\cup\{\boldsymbol{\alpha}_e\}\\ \boldsymbol{N}=(\boldsymbol{N}-\{\boldsymbol{\alpha}_e\})\cup\{\boldsymbol{\alpha}_l\} \end{cases}. {B=(B−{αl​})∪{αe​}N=(N−{αe​})∪{αl​}​.
第17、18行分别计算轴转后的基矩阵及其逆阵。第19~20行计算向量$\mathbf{z}={\mathbf{c}}_{\mathbf{B}}^{\top }{\mathbf{B}}^{-1}\mathbf{N}$，赋予z。第21行返回infinite, Bind，Nind和z。
例1：标准型线性规划
$$\{\begin{array}{l}\text{minimize}-x_1-5x_2\\ \text{s.t.}5x_1+6x_2+x_3=30\\ 3x_1+2x_2+x_4=12\\ x_1,x_2,x_3,x_4\ge 0\end{array}$$
目标函数$f(\mathbf{x})=-x_1-5x_2={\mathbf{c}}^{\top }\mathbf{x}$，其中$\mathbf{x}\in {\text{R}}^4$，${\mathbf{c}}^{\top }=(-1,-5,0,0)$。等式约束系数矩阵和常数向量
$$\mathbf{A}=\left(\begin{array}{cccc}5& 6& 1& 0\\ 3& 2& 0& 1\end{array}\right),\mathbf{b}=\left(\begin{array}{c}30\\ 12\end{array}\right).$$
取基矩阵${\mathbf{B}}_1=({\mathbf{\alpha }}_3,{\mathbf{\alpha }}_4)=\mathbf{I}$，非基矩阵${\mathbf{N}}_1=({\mathbf{\alpha }}_1,{\mathbf{\alpha }}_2)$。对应的基本可行解${\mathbf{x}}_1=\left(\begin{array}{c}0\\ 0\\ 30\\ 12\end{array}\right)$，${\mathbf{x}}_{\mathbf{B}}=\left(\begin{array}{c}30\\ 12\end{array}\right)$，${\mathbf{x}}_{{\mathbf{N}}_1}=\mathbf{o}$。${\mathbf{c}}_{{\mathbf{B}}_1}^{\top }=(0,0)$，${\mathbf{c}}_{{\mathbf{N}}_1}^{\top }=(-1,-5)$，目标函数值$f({\mathbf{x}}_0)={\mathbf{c}}_{{\mathbf{B}}_1}^{\top }{\mathbf{x}}_{{\mathbf{B}}_1}+{\mathbf{c}}_{{\mathbf{N}}_1}^{\top }{\mathbf{x}}_{{\mathbf{N}}_1}=0$。下列代码对基矩阵${\mathbf{B}}_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([[5, 6, 1, 0],							#设置系数矩阵
              [3, 2, 0, 1]])
b = np.array([30, 12])								#设置常数向量
c = np.array([-1, -5, 0, 0])						#设置目标函数系数向量
Bind = np.array([2, 3])								#初始基矩阵列向量下标集
Nind = np.setdiff1d(np.arange(4), Bind)				#初始非基矩阵列向量下标集
B = A[:, Bind]										#初始基矩阵
B1 = np.linalg.inv(B)								#逆矩阵
w = np.matmul(c[Bind], B1)
z = np.matmul(w, A[:, Nind]).flatten()				#计算向量z
_, Bind, Nind, _ = pivot(A, b, c, Bind, Nind, z)	#轴转
print(Bind, Nind)
B = A[:, Bind]										#新的基矩阵
B1 = np.linalg.inv(B)
xB = np.matmul(B1, b.reshape(b.size,1)).flatten()	#基变量解
x = np.zeros(4)										#初始化可行解
x[Bind] = xB										#填充基变量值
f = np.dot(c, x)									#目标函数值
print(x, f)

代码中第2~4行设置输出数据为有理数（分数）格式。第5~8行设置本题标准型线性规划的等式约束系数矩阵，常数向量，目标函数系数向量A、b、c。第9、10行分别设置初始基矩阵和非基矩阵的列向量下标。13~14行计算判别向量中的z。第15行调用pivot函数计算轴转，将返回值中Bind，Nind赋予同名变量。17、18行重新计算基矩阵B及其逆矩阵B1。第19行计算基变量解，20~21行计算新的基本可行解x，22行计算新的基本可行解处函数值。运行程序，输出

[2 0] [3 1]
[4 0 10 0] -4.0

第1行表示新的基矩阵由${\mathbf{\alpha }}_3,{\mathbf{\alpha }}_1$构成（注意Python数组下标从0开始编码），非基矩阵由${\mathbf{\alpha }}_4,{\mathbf{\alpha }}_2$构成。第2行表示对应B的新的基本可行解为${\mathbf{x}}_1=(4,0,10,0{)}^{\top }$，对应的目标函数值为$f({\mathbf{x}}_1)=-4<0=f({\mathbf{x}}_0)$。
写博不易，敬请支持：
如果阅读本文于您有所获，敬请点赞、评论、收藏，谢谢大家的支持！