单纯形法
单纯形法是针对求解线性规划问题的一个算法，这个名称里的'单纯形'是代数拓扑里的一个概念，可以简单将'单纯形'理解为一个凸集，标准的线性规划问题可以表示为:
min（or max） f(x)=cx
s.t. Ax=b
x>=0,b>=0
以上形式称为线性规划标准型，使用单纯型法时，如果约束条件含有不等式时需新增变量（松弛变量、人工变量）转化为标准型，min. f(x)=cx指求函数最小值（也可以是求最大值），x是一个Rn维向量代表有n个变量，线性规划问题主要是面向实际问题，x变量可以代表距离、成本、价格、数量等，线性规划问题中要求x大于等于0，c同样是一个Rn维向量，这样cx实际上就是一个线性函数f(x);s.t.代表subject to代表服从于意思，这里是指变量x需要满足的约束条件，A是一个Rm*n维矩阵，代表有m个等式约束。下面是一个约束是不等式的情形:
min -4x1-x2
s.t. -x1+2x2<=4
2x1+3x2<=12
x1-x2<=3
x1,x2>=0
求解上面这个问题只要初中数学知识即可，具体可以使用代数法或几何的方法轻松得到，考虑到实际问题当中变量x是多维的，约束条件也会比示例多的多，这就需要一个一劳永逸的算法能通过计算机来获得正解，单纯形法就是这样的一个算法。单纯形法最早由 George Dantzig于1947年提出，单纯形法对于求解线性规划问题是具有跨时代意义的，其实不仅仅是针对线性规划，非线性规划问题在求解的过程中也大量依赖单纯形法，George Dantzig本人也被称为'线性规划之父'，他是好莱坞电影《心灵捕手》的原型。
一、凸集及其性质介绍
1.1 凸集概念
单纯形可以理解为一个凸集，介绍单纯形法之前有必要先来了解一下凸集概念及其性质，在很多机器学习算法中都会穿插凸集的概念，凸集上求最优解是凸优化的一个分支，凸集对于简化问题是有重要意义的。凸集可以用图形表示为一个没有洞的联通区域，如下图所示：
凸集可以这样描述:用一条直线连接集合里两个元素，这条连线上的所有元素都在这个集合里，这个集合称为凸集；用代数公式表达是，对于任意x ,y
S,都有αx+(1-α)y
S（α>0）,则S为凸集:
根据定义下面这个图形所定义的区域不是凸集：
回过头来看线性规划的约束条件，将满足约束条件Ax=b的集合称为S,也称S为可行区域，对于其中任意x ,y
S有:
Ax=b,Ay=b
x和y之间任意一个元素可以表示为 v= αx+(1-α)y,计算Av=A[αx+(1-α)y]=αb+b-αb=b,显然有v
S，所以说线性规划中的约束条件描述的是一个凸集。
1.2 凸集的极点（顶点）
从凸集代数定义公式可以发现，凸集中任何一个点都可以用另外两个点线性组合而成，凸集中有一些特殊点用公式表达:x= αx1+(1-α)x2 x1 ,x2S,如果有且只有x=x1=x2,则称x为极点，或者形象的称为凸集的顶点；从图像上来看顶点不在凸集两个不同点连线上，如下图所示：
这里的x1,x2,x3,x4,x5都是顶点,而x不是顶点，特殊的，圆形上的圆周上每点都是顶点，所以说凸集上顶点可以是有限的，也可以是无限多的。在线性规划问题中，我们也把上面凸集称为可行区域，上图x1,x2,x3,x4,x5称为可行区域的顶点，x称为可行点。
对于凸集还有其他一些概念，虽然在单纯形法中没有用到，这里可以顺带提一下，凸集的区域范围可以有界或者是无界的，当凸集有界时，线性规划目标函数是把两个Rn维向量做内积运算得到一个实数R,所以这里本质上是一个有界线性泛函过程，如果没有约束条件的话，显然x向量与c向量呈180度相反方向时可得到最小值，c在这里可以理解为目标函数的梯度。
当凸集无界时时，需要引入方向的概念,设x是Rn维向量，x沿着一个方向延伸得到新的一个点v：v=x+αd x
S;α>0;d>0，如果新的点v
S,则称d为可行方向。从线性代数角度来看，方向作为一个向量，可以由其他的方向线性合成，如果凸集一个方向不能由其他方向合成则称这个方向为极方向，这句话也可以反过来理解，凸集的所有可行方向都可以由极方向线性表示，极方向与极点有着对应关系。凸集的可行方向对于理解'可行区域求最优解'、'KKT条件'等非线性规划有着重要的意义。
二、单纯形法原理
2.1 可行解、最优解、基本可行解
根据上面的描述，凸集中每个元素都满足约束条件，凸集中每个元素称为目标函数的可行解，显然可行解的数量是非常多的，可行解带入目标函数后一般得不到函数的最小值或最大值；如果可行解带入目标函数后可得到最值（最大或最小值），那么这个可行解就成为最优解。接下来介绍基本可行解，基本可行解是实现单纯形法的关键，还是以上面的问题为例，我们新增三变量x3,x4,x5使得约束条件变为等式标准型：
min -4x1-x2+0x3+0x4+0x5
s.t. -x1+2x2+x3+0x4+0x5=4
2x1+3x2 +0x3+x4+0x5=12
x1-x2 +0x3+0x4+x5=3
x1,x2,x3,x4,x5>=0
这里引入的x3,x4,x5称为松弛变量，松弛变量同样是大于等于0的变量，这时原问题常数向量c=[-4,-1,0,0,0],由于引入松弛变量系数都是0，所以新的目标函数和原目标函数是完全等价的，这时约束条件系数是矩阵A=
秩为3，A表示的是一个三维空间，如果将后面三列单位矩阵作为这个三维空间的三个基（也可以理解为三维空间三个坐标轴），A前面2列所代表向量可以通过单位矩阵线性表示出来，当然也可以将矩阵A中任意三列作为三维空间的基，将A中列分成两类，一类列是线性空间的基用B表示，一类列不是基用N表示，这样约束条件可写成矩阵形式：
上面公式中将与基B相乘的变量x3,x4,x5称为基变量，基变量一般用xB表示；与基N相乘的变量x1,x2称为非基变量一般用xN表示，每一个变量都与B或N中一个列对应，上面的公式也可以写为:
如果令xN=0,则有
,本例中B是R3*3满秩矩阵，即B有逆矩阵B-1,可得
。约束条件中要求x>=0,如果此时恰好有xB=B-1b>=0也就满足了所有约束条件，x=（xN,xB）=（0,0，xB）也是目标函数的一个可行解，由于此时非基变量都为0，称此可行解为基本可行解,由此可见，基本可行解中'基'字指此时可行解基变量不全为0，而非基变量皆为0，请注意（0,0，xB）是基本可行解，xB是基本可行解中非0的分量集合，xB本身并不是基本可行解。
从线性空间的角度来理解线性规划问题，目标函数在新增松弛变量后x
R5，即目标函数中的x是一个来自5维空间的向量，Ax=b的含义是通过线性映射将一个5维空间的向量映射为一个固定的三维向量 b,当然从函数角度来看这个映射过程是多对一的关系，即有多个5维的向量线性映射后变为一个3维向量b,这些5维向量组成的集合之前已经介绍过，是一个凸集，目标函数的含义此时可以表述为，在5维空间组成的凸集中找到一个或多个向量，与一个固定向量c的内积最小。
2.2 基本可行解、凸集顶点（极点）
基本可行解与凸集顶点是等价的，而凸集的顶点在线性规划问题中与最优解又有着密切的联系，首先来证明一个重要的定理：凸集中元素x是基本可行解的充分必要条件为x是凸集的顶点。这个定理是实现单纯型法的关键，下面来详细证明一下：
先来证明这个定理必要性:即基本可行解一定是凸集的顶点，设x是基本可行解，根据之前分析可以知道向量x中必有前k个分量大于0，其余分量等于0，这个k个分量对应的是矩阵k个线性无关的列（基），设可行区S内有x1,x2,x1≠x2,根据凸集的定义设基本可行解x=αx1+(1-α)x2,由于基本可行解x前k个分量大于0而其余的分量等于0，又由于α选择的任意性，所以当i>=k+1时，必有x1[i]=x2[i]=0 (x1[i]代表x1第i个分量)。又因为Ax1=Ax2=b，即A(x1-x2)=0,已经知道x1,x2的分量从k+1开始都是0，此时就可以得到
，这里Ai是矩阵中几个线性不相关的列，前面已经知道x1≠x2，则Ai的系数x1[i]-x2[i]必然有不等于0的情形，而如果出现不等于0情况说明这些列Ai是线性相关的，这显然与假设矛盾，所以此时可以得到一个结论，没有两个不同的点x1,x2使得有基本可行解x=αx1+(1-α)x2存在，要使等式成立必须x=x1=x2,满足这个条件时x显然是凸集顶点（极点）。
接下来证明充分性：x是凸集的顶点则x一定是基本可行解。设顶点x有k个分量大于0，其余等于0，由于约束条件要求x分量必须大于等于0，这样分类显然是成立的。 如果这个k个分量对应矩阵的列是线性无关的，那么顶点x就是基本可行解，这里我们可以用反证法，先假设k个分量对应矩阵列是线性相关的，即存在一个分量不全为0的向量
使得有等式：
，对于顶点x而言，前k分量也有相应的等式:
,联合这两个等式同时引入一个实数b，b>0,得到下面两组等式：
这样就从顶点x构造出2个新的向量x1,x2，通过上面公式可以发现，x1、x2也是凸集中元素,上面两个等式可以得到x1=x-ab,x2=x+ab，这里顶点x是向量，由于向量a是不全为0的向量且b是正实数，必然有x1≠x2；通过x1,x2可以反过来得到x的表达式x=x1/2+x2/2,因为x1≠x2这样等式显然与x是顶点的定义矛盾，此时就推翻最初假设得到新的结论：k个分量对应矩阵列是线性不相关的，这同时证明了凸集的顶点则x一定是基本可行解。
2.3 凸集顶点与最优解
线性规划的目标函数是一个线性函数f(x)=cx,这个函数如果没有约束条件的束缚，其值域是正负无穷的，线性函数在定义域上等值线是相互之间平行的向量，等值线不会像非线性规划那样范围逐渐变小，进而变成一个点，这也导致了线性规划问题最值有可能是多个，加入约束条件后，目标函数最值只能在一个凸集中选择，线性规划的等值线运动过程如下图所示：
其中虚线代表目标函数f(x)=cx的等值线，ABCD所包含区域即为满足约束条件的凸集，箭头的方向代表目标函数的梯度，在线性规划问题中这个梯度就是向量c。当虚线沿着箭头方向平移时代表目标函数值不断变大，如果是求最小值则沿着箭头反方向移动，当然这个平移是有限制的，即不能离开凸集所在区域ABCD，最终目标函数一定在顶点A处取得最小值。
我们可以得到这样一个结论，线性规划问题的最值一定在凸集的顶点处。通过之前的分析已经知道，凸集的顶点是满足约束条件的基本可行解，单纯形法的核心思想可以归纳为：找到每一个基本可行解，代入目标函数后计算函数值取其最大或最小值即可。单纯形法上从代数角度是寻找约束条件的每一个基本可行解，从几何意义上来说是遍历凸集的每一个顶点，根据算法的特性有时也称为转轴法。
三、单纯形法的实现
3.1 入基、出基
单纯形法过程就是不断找出基本可行解的过程，基本可行解有一个显著的特征：基变量不为0，非基变量都为0。利用这个特性我们能够很快的找出每一个基本可行解。还是之前的标准型为例：
min -4x1-x2+0x3+0x4+0x5
s.t. -x1+2x2+x3+0x4+0x5=4
2x1+3x2 +0x3+x4+0x5=12
x1-x2 +0x3+0x4+x5=3
x1,x2,x3,x4,x5>=0
约束条件系数矩阵A是一个秩为3的矩阵，代表A是一个有3个基的3维线性空间，约束条件可以描述为A矩阵作用于一个5维向量x1,x2,x3,x4,x5后得到一个三维向量b=(4,12,3),这说明A矩阵是一个线性映射过程，将原来有5个基线性空间线性映射为一个3个基的线性空间，原来5维空间有2个基在线性映射后变成了0，这当然不是偶然的，因为线性映射有一个特性：零空间的维度（原来空间被映射为0的个数）+值域的维度（映射后空间的维度，本例是3）=A的列数（原来空间的维度），这个定理可以通过反证法来证明，这里不详细介绍。
基本可行解的特性决定了在原来5维空间中，每个顶点的形式一定是有3个分量大于0，另外2个等于0，所以顶点形式按排列组合则最多有
种，这与线性映射的机制完全吻合，我们不妨设原来5维空间有5个标准正交基：
假设一个在5维空间凸集中有一个顶点x形式为 (x1,x2,x3,0,0)，写成向量与坐标相乘的形式：
经过线性映射后原线性空间的3个5维基向量被映射为3个3维基向量，另外有2个基被映射为0向量：
线性映射后新向量也可以写成基变量与单位矩阵相乘的形式：
以上线性映射过程中产生了两个数学模型，一个是凸集空间，一个是线性方程模型，5维空间中顶点与线性方程中基本可行解一一对应，这个总结中一一对应非常重要，一一对应从函数角度来说是满射、单射，是可逆的，凸集顶点x:(x1,x2,x3,0,0)形式决定线性映射的方式，从而决定了xB；同样的也可以通过调整xB逆向得到5维空间中不同的顶点，xB是约束条件中基本可行解中大于0的分量，也就是说不断调整基本可行解就可以得到相应的5维空间中不同的顶点，这就像一系列联动的齿轮轴组，调节一端就可以影响另一端，这样一来就能很好理解单纯形法为什么也叫转轴法了。
说明基本可行解时说过，xB是在矩阵A中选择三个列作为单位矩阵得到的，同时提到，也可以选择矩阵其他的列作为3维空间的基，选择不同列作为3维空间基即可产生不同的xB，从而得到5维空间不同的顶点。新选择列称为入基，被替代的列称出基，这种替换的可以有许多种方式，使用暴力方法逐个迭代每一个顶点也可以实现，而当求解最小值问题时，如果借助梯度的概念，在得到一个顶点后，切换下一个顶点方向是函数变小的方向无疑会节省很多时间。
约束条件系数矩阵A中将作为基的列集合记作B，不是基的列集合记作N,则有A=(N，B),同样把x分量中基变量的集合记作xB，非基变量集合记作xN,约束条件可以写成： (N，B)（xN，xB）=b,得到等式 NxN+BxB=b,矩阵B是基矩阵且可逆，得到:
xB=B-1b - B-1NxN ⑴
当xN=0 时显然有xB=B-1b ，如果此时有B-1b>0则称xB是初始基本可行解,目标函数中同样将向量c也分成两部分，基变量对应系数记为cB, 非基变量对应的系数记为cN,这样就得到目标函数另一种形式：
f(x)=cBxB+cNxN ⑵
带入公式（1）后得到目标函数另一种形式:
如果此时有初始基本可行解，可得到f(x)=cBxB=cBB-1b。对上面的函数求导得：
，
大于0时导数小于0，说明在取初始基本解后目标函数还有变小的空间，这时可以选择矩阵非基集合中一列来代替现有基，即选择入基，如何选择入基呢?函数的微分形式有：
越大则导数
越小,目标函数变小的幅度就越大，将此作为选择入基的标准，遍历每个非基列pj,由于每次都只一个列作为入基，xN其他分量依然为0，
简化为cBB-1pj-cj,引入每一个待选列的检验量tj=cBB-1pj-cj, 如果检验量大于0，就选择这其中检验量最大的列作为入基就能保证目标函数变小，如果待选列的检验量都小于等于0则代表此时已经得到最小值。
入基已经准备就绪如何选择出基呢，即如何选择被替代的基列？在选择入基时其实还有一个隐含的问题，入基被选择后，新的基变量值是多少？这些可以从公式（1）中得到答案。
根据以上分析过程完成开篇线性规划问题，制成python代码如下：
import numpy as np
BASEINDEX=2
#列变换实现单纯形法
def simplex_ColTranslate(c ,A,b,flagstable):
B=A[:,BASEINDEX:]
Binv=np.linalg.inv(B)
xb=np.dot(Binv ,b)
cb=c[:,BASEINDEX:]
f=np.dot(cb,xb )
Inbaseindex=-1
OutBaseIndex = -1
tempvalue=0
for i in range(BASEINDEX):
v=np.squeeze(np.dot(np.dot(cb,Binv),A[:,i])-c[:,i],0)
if tempvalue<=v:
#找出入基
tempvalue=v
Inbaseindex=i
if tempvalue==0:
#找到最优解
return xb,f,c,flagstable
else:
#入基出基替换，同时交换系数c
y=np.dot(Binv,A[:,Inbaseindex])
minivalue=None
for i in range(y.shape[0]):
if y[i]>0:
x=xb[i]/y[ i]
#找出最小值，保证xb>=0
if minivalue is None or minivalue>x :
minivalue=x
OutBaseIndex=i+BASEINDEX
for i in range(A.shape[1]):
if i==OutBaseIndex:
for j in range(A.shape[0]):
tmp=A[j,i].copy()
A[j,i]=A[j, Inbaseindex]
A[j, Inbaseindex]=tmp
tempc= np.squeeze( c[:,OutBaseIndex],0).copy()
c[:,OutBaseIndex]=c[:,Inbaseindex]
c[:,Inbaseindex]=tempc
flagstable[OutBaseIndex-BASEINDEX]=Inbaseindex
return simplex_ColTranslate(c, A, b,flagstable)
def useColTranslate():
BASEINDEX = 2
c = np.array([[-4, -1, 0, 0, 0]],dtype=np.float)
A = np.array([[-1, 2, 1, 0, 0], [2, 3, 0, 1, 0], [1, -1, 0, 0, 1]],dtype=np.float)
b = np.array([[4], [12], [3]],dtype=np.float)
# 记录下标变动情况
flagstable = {}
for i in range(c.shape[1] - BASEINDEX):
flagstable[i] = BASEINDEX + i
x, f, c, flags = simplex_ColTranslate(c, A, b, flagstable)
print('最优解:%.2f' % (f[0, 0]))
for j in range(x.shape[0]):
print('x%d=%.1f' % (flags[j] + 1, x[j, 0]))
if __name__ == '__main__':
useColTranslate()运行结果如下：
3.2 单纯形法的表格方法
上面代码中设置了一个常量BASEINDEX用以指示A矩阵中基矩阵的开始索引，此后通过检验量将入基替换出基，将新的列置换出去，同时更新目标函数中的c向量的顺序，将cB与cN中分量同时完成置换，c中分量与x中分量如影随行，上面算法本质上是通过矩阵列变换实现单纯形法，好处是容易理解，与数学推导过程完全一致，弊端是打乱了原来矩阵列的顺序，不易于二次利用其中的参数，如在接下来介绍的两阶段法中需要二次利用矩阵A，列变换的处理方式为后期处理带来很大的难度。
行变换的方法也称单纯形法表格法，行变换不是将入基与出基对调，而是通过消元法将入基变成单位向量，在这个消原过程中，将初始参数b向量放到A矩阵右侧一起参与消元法，加入b后的矩阵称为增广矩阵，表格法选择检验量以及找到最值的过程都是与列变换一致的。初等行变换相当于矩阵A左乘一个矩阵，初等列变换相当于矩阵A右乘一个矩阵，两种变换殊途同归，都是得到一个与xB维度一样的向量，即基本可行解的非零部分，接下来举一个使用表格法的例子：
引入松弛变量x5,x6得到下面的标准型：
接下来建立单纯形表如下图所示：
表的最左侧代表这个阶段的基变量，此时基变量是x4, x5,x6,中间是约束条件矩阵系数A,系数上方的变量名与下面的列关联；右侧是标准型中b向量，最下面一行左边是检验量，右边单位格处是目前的函数值，这个单元格值等于cBB-1b，即基本可行解代入函数后得到的函数值。观察一下检查量，x2这一列检查量大于0且值最大，这时x2需将会升级为基变量，再观察与x2关联的非基向量(1,-1,2),其中只有系数1,2大于0，min{10/1,4/2}=2,这时就找到了x6，与x6关联的基向量将会变成非基向量，我们称2为本次检验的主元，主元是接下来进行行变换乘法因子。
接下来需要通过行变换把x2关联的非基向量(1,-1,2)变成单位向量(0,0,1)，方法是首先主元所在的第三行除以主元，将主元2变为1，然后将调整后第三行加上第二行使得-1变为0，再将第三行乘以-1加上第一行，这样x2关联的列向量(1,-1,2)就变成了(0,0,1)，而此时x6所关联的列不再是单位向量形式，代表已经变成了非基向量。注意这个行变换过程中，b也同时参与运算，这样单纯形表就变为下面的形式：
这时左侧基变量变x6被x2取代，右侧b变成了（8,10,2），由于xB=B-1b,而B始终是单位矩阵，这里x4=8,x5=10,x2=2,
在这一步中重复上面的操作，将x3与x5置换，并通过行变换使向量(0,2,-2)变为（0,1,0），得到下面的表格：
此时最下面一行的检验量都小于等于0代表已经找到最小值，最小值为-19，基变量为x2=12,x3=5,x4=8,这里强调一点，这样的结果代表在原来凸集中顶点是（0,12,5,8,0,0），而对应的基本可行解中非0部分xB的形式是(8,5,12),他们之间非0部分值顺序有可能是不一样的，这是线性映射的结果，根据行变换的算法制成python程序如下：
import numpy as np
#从矩阵得到单位基矩阵
def getbaseIndexsfromMatrix(A):
index=[]
for i in range(A.shape[1]):
if(np.dot(A[:,i],A[:,i].T)==1.0 and np.max(A[:,i])==1.0 ):
index.append(i)
return index
#从基变量表得到单位基矩阵
def getbaseIndexsfromVariables(flags):
index=[]
for key in flags:
index.append(flags[key])
return index
#行变换-单纯形表格法,min=True求最小值
def simplex_RowTranslate(c ,A,b,flagstable,min=True):
#获得单位基矩阵
index=getbaseIndexsfromVariables(flagstable)
cb = c[:, index]
B=A[:,index]
Binv=np.linalg.inv(B)
tempvalue,Inbaseindex,OutBaseIndex = 0,-1,-1
checkvalue=np.zeros(A.shape[1] )
for i in range(A.shape[1]):
if i not in index:
v = np.squeeze(np.dot(np.dot(cb, Binv), A[:, i]) - c[:, i], 0)
checkvalue[i] = v
# 找出入基
if min:
if tempvalue < v:
tempvalue = v
Inbaseindex = i
else:
if tempvalue > v:
tempvalue = v
Inbaseindex = i
if tempvalue<=0 and min :
return b, c, A,flagstable
elif tempvalue>=0 and min==False:
return b, c, A, flagstable
else:
minivalue,privot,proitrow = None ,None,None
for r in flagstable:
if A[r,Inbaseindex] > 0:
x = b[r,0] /A[r,Inbaseindex]
if minivalue is None or minivalue>x :
minivalue=x
OutBaseIndex=flagstable[r]
#所在行，接下来要从该行开始做行变换
proitrow=r
#主元
privot=A[r,Inbaseindex]
A[proitrow,:]=A[proitrow,:]/privot
b[proitrow, 0] = b[proitrow, 0] / privot
for i in range(A.shape[0]):
if i!=proitrow and A[i,Inbaseindex]!=0:
r = -A[i, Inbaseindex]
b[i, 0] = b[i, 0] + r * b[proitrow, 0]
for j in range(A.shape[1]):
A[i, j]=A[i, j]+r*A[proitrow, j]
for key in flagstable:
if flagstable[key]==OutBaseIndex:
flagstable[key] = Inbaseindex
r2 = -checkvalue[Inbaseindex]
tempvalue=0
for i in range(checkvalue.__len__() ):
checkvalue[i] = r2 * A[proitrow, i] + checkvalue[i]
if min:
if tempvalue <= checkvalue[i]:
tempvalue = checkvalue[i]
else:
if tempvalue >= checkvalue[i]:
tempvalue = checkvalue[i]
if tempvalue<=0 and min:
return b,c, A, flagstable
elif tempvalue >= 0 and min == False:
return b,c, A, flagstable
else:
return simplex_RowTranslate(c, A, b, flagstable,min)
def useRowTranslate( ):
c = np.array([[1, -2, 1, 0, 0,0]], dtype=np.float)
A = np.array([[1, 1, -2, 1, 0,0], [2, -1, 4, 0, 1,0], [-1, 2, -4, 0,0, 1]], dtype=np.float)
b = np.array([[10], [8], [4]], dtype=np.float)
# 基变量表
flagstable = {}
index = getbaseIndexsfromMatrix(A)
for i in range(index.__len__()):
flagstable[i] = index[i]
x, c, AA, flags = simplex_RowTranslate(c, A, b, flagstable)
index = getbaseIndexsfromVariables(flagstable)
cb = c[:, index]
B = AA[:, index]
Binv = np.linalg.inv(B)
f = np.dot(np.dot(cb, Binv), b)
print('最优解:%.2f' % (f[0,0] ))
for j in range(x.shape[0]):
print('x%d=%.1f' % (flags[j] + 1, x[j, 0]))
if __name__ == '__main__':
useRowTranslate( )运行结果如下图，可以发现结果与之前运算的结果一致的：
通过比较可以发现行变换的表格法实用性更强，由于每次行变换新的入基总是以单位向量形式出现，此时不再需要指定矩阵A中单位矩阵起始位置，程序代码自动识别基变量以及基列，具体实现表格法的函数是simplex_RowTranslate，其中有一个参数min,当min为真是求解线性规划最小值，min为假时求函数最大值。
四、大M法、两阶段法
4.1 人工变量的引入
上面的例子中在引入松弛变量变成标准型后，系数矩阵A自动有一个单位矩阵，单位矩阵的阶与条件个数一致，A秩为m的Rm*n的矩阵，即A是满秩的。有时加入松弛变量后并不能得到单位矩阵，遇到这种情形时自然不能得到初始基本可行解，比如下面这个例子:
加入松弛变量x3,x4,x5后变为等式的标准型：
上面的标准型矩阵系数A没有单位矩阵，之前算法在这里无效了，如果此时再引入两个变量x6,x7,同样的要求x6,x7>=0，这时就可以得到标准型了：
上面的标准型中x5,x6,x7是基变量,A是R3*7的满秩矩阵，x6,x7被称为人工变量，针对本例是求最大值目标函数可以写为:
观察这个目标函数可以发现松弛变量与人工变量的不同，松弛变量在函数中系数为0，单纯形法处理后松弛变量本身可以不为0，这对求解目标函数没有任何影响，而人工变量前面有系数M，且M>0,这就要求人工变最终必须为0，否则新的目标函数最值与原目标函数是不相等的，M有时也叫做惩罚系数，意味着如果人工变量不为0原函数将永远不能取得最大值。