单纯形法

提出单纯形的思路

　　我们知道，线性规划(LP)问题如果有最优解，必可在某个极点（基本可行解）上达到。一个直观的想法是：对于LP问题，找出所有的基本可行解，然后逐个比较，即枚举法。但是事实上，时间开销会非常大，假设原问题中有n个变量，m个约束条件，则时间开销为$C^{m}_{n}$，而$C^{m}_{n}$会随着m，n的增大而迅速地增大。显然，这不可行。

　　换种思路：若从某一基本可行解出发，每次总是寻求比上一个更“好”的基本可行解，直至找到最优解。这样就可以大大减少计算量，其实这样的思想应用在非常多的地方，比如梯度下降等等。

　　如果我们要使用这种迭代的求解方法，我们需要解决下面的三个问题:

• (1) 如何判断当前的基本可行解是最优解（迭代终点）；
• (2) 如何寻找下一个改善的基本可行解（迭代关系式）；
• (3) 如何得到一个初始的基本可行解（迭代起点）；

　　美国数学家丹奇格（G.B.Dantzig）解决了上面的三个问题，也就得到了一种快速求解LP问题的方法，我们称为单纯形法。

通过一个例子进行分析

　　举个例子。

问题描述

　　设有某LP问题如下

　　$max z=cx\ s.t.
left{ egin{gathered}
Ax=b hfill \
xgeq 0 hfill \
end{gathered} ight.$

　　其中，$c=egin{bmatrix}5&2&0&0&0end{bmatrix}$，$x=egin{bmatrix}x_1\x_2\x_3\x_4\x_5end{bmatrix}$，$A=egin{bmatrix}30&20&1&&\5&1&&1&\1&&&&1\ end{bmatrix}$，$B=egin{bmatrix}160\15\4end{bmatrix}$。

解：　　

　　令$B^{(0)}=(p_3,p_4,p_5)=egin{bmatrix}1&&\&1&\&&1end{bmatrix}$为基矩阵，即$x_3,x_4,x_5$为基变量，$x_1,x_2$为非基变量。

　　使用非基变量表示基变量，有$left{ egin{gathered} 
x_3=160-30x_1-20x_2 hfill \ 
x_4=15-5x_1-x_2 hfill \ 
x_5=4-x_1 hfill \ 
end{gathered} ight.$。因此，目标函数的非基表示为$z=0+5x_1+2x_2$。

　　令$x_1=0,x_2=0$，则得到一个基本可行解为$x^{(0)}=(0,0,160,15,4)^T$，目标函数$z^{(0)}=0$，显然这不是最优解。我们需要换一组基变量（即从一个顶点移动到另一个顶点）。于是，我们选取$(x_1,x_2)$中的一个变量（换入基变量），$(x_3,x_4,x_5)$中的一个变量（换出基变量），进行换入换出。现在，假设换入基变量为$x_1$，下面确定换出基变量。(在这里我们称进行基变换的两个变量为换入基变量和换出基变量，也可以称为进基变量和离基变量)

　　由于$x_2$仍为非基变量，故$x_2$仍取零值。因此，有$left{ egin{gathered} 
x_3=160-30x_1 hfill \ 
x_4=15-5x_1 hfill \ 
x_5=4-x_1 hfill \ 
end{gathered} ight.$由于有$xgeq 0$，因此$x_1=min{160/30,15/5,4/1}=3$，故$x_4$为换出基变量，这样得到一个新的基本可行解$x^{(1)}=(3,0,70,0,1)$，相应的基矩阵$B^{(1)}=(p_3,p_1,p_5)$，基变量为$x_3,x_1,x_5$，非基变量为$x_2,x_4$，对应的目标函数值$z^{(1)}=15>z^{(0)}=0$。因此，$x^{(1)}$是比$x^{(0)}$改善的基本可行解。

　　下面分析$x^{(1)}$是否是最优解。在上面我们得到目标函数值$z=0+5x_1+2x_2$，根据约束方程组，我们可以得到$x_1=3-frac{1}{5}x_2-frac{1}{5}x_4$，代入目标函数式，我们得到$z=15+x_2-x_4$。观察左式，我们发现非基变量$x_2$前的系数为正，增加$x_2$，目标函数值仍可以增大，或者说，让$x_2$做换入基变量，进行换入，我们有可能使目标函数值增大。因此，$x^{(1)}$不是最优解。

　　使用相同的方法，以$x_2$为换入基变量，我们得到$x_3$为换出基变量，得到新的基本可行解$x^{(2)}=(2,5,0,0,2)^T$，目标函数$z=20-frac{1}{14}x_3-frac{4}{7}x_4$，在左式中，我们可以看出非基变量$x_3$或$x_4$由零值增加，只会使得目标函数值$z$减少，故$x^{(3)}$是最优解。

　　下面归纳一下单纯形法的基本解题步骤。

单纯形法的基本过程

　　第一步：构造一个初始的基本可行解。

　　第二步：判断当前基本可行解是否为最优解。

　　第三步：若当前解不是最优解，则要进行基变换迭代到下一个基本可行解。

最优性检验

　　这一部分，我们来判断当前解是否是最优解。我们将目标函数化为如左的形式$z=z_0+sum_{j=m+1}^{n}{sigma_jcdot x_j}$。我们称$sigma_j$为非基变量$x_j$的检验数。有以下定理

　　若关于非基变量的所有检验数$sigma_jleq 0$成立，则当前基本可行解$x^*$就是最优解。

基变换

　　我们通过基变换迭代到下一个改善了的基本可行解。问题在于，我们如何选择换入基变量和换出基变量。

• (1) 换入基变量的选择   对于$z=z_0+sum_{j=m+1}^{n}{sigma_jcdot x_j}$，若有两个以上的$sigma>0$，则选最大$sigma$。
• (2) 换出基变量的选择   假设换入基变量为$x_k$。在上一轮迭代中的多个基变量中，增长$x_k$最先达到0的基变量为换出基变量。

无穷多个最优解及无界解的判定

　　实际上，有些LP问题无解或者会有无穷多个解。这两类情况有些特殊，需要在迭代中进行额外的判断，防止程序陷入死循环或者出错。

• (1) 无穷多个最优解的判定   如果当前基本可行解$x^{(k)}$的所有非基变量的检验数$sigma_i$满足$sigma_ileq 0$，且其中一个$sigma_j=0$，则该LP问题有无穷多个最优解。
• (2) 无界解的判定   如果当前基本可行解$x^{(k)}$，其中一个非基变量$x_i$的检验数$sigma_i>0$，且$x_i$对应的系数列向量$p_i=(a_1,a_2,...,a_m)$中，所有分量$a_jleq 0$，则该线性规划问题具有无界解（或者称无最优解）。这个可能比较费解，换种说法，可能会容易理解一点：某次迭代中，对于换入基变量$x_i$，无法找到一个换出基变量，则此问题无界。

单纯形法的python实现

　　在编程的过程中，我们会思考一些具体的细节，也许会有新发现哦。

　　

#python 3
import numpy as np
from itertools import combinations

# 实体类　Solution
# 控制类　Simplex

class Solution:
	def __init__(self):
		pass

	def set_para(self,A,b,z):
		# A m*n
		# b m*1
		# z 1*(m+1)
		self.A=A
		self.b=b
		self.z=z
		
		self.m,self.n=A.shape
		self.x_index=[i for i in range(self.n)]

	def get_init_solution(self):
		for JB in combinations(range(self.n),self.m):
			if self._is_solution(JB):
				JB,JN=self._rearrange(JB)
				self._set_JB_JN(JB,JN)
				return True
		return False

	def _is_solution(self,JB):
		B=np.hstack([self.A[:,i] for i in JB])
		# B的行列式的值不为0，则B是一个可逆矩阵
		if np.linalg.det(B):
			return True
		return False

	def _rearrange(self,JB):
		JN=[i for i in range(self.n) if i not in JB]
		B=np.hstack([self.A[:,i] for i in JB])
		N=np.hstack([self.A[:,i] for i in JN])
		self.z=[self.z[0,i] for i in JB]+[self.z[0,i] for i in JN]+[self.z[0,-1]]
		self.z=np.matrix([self.z],dtype=float)
		# 将B转化为单位矩阵
		# 即(B|N)x=b -> (I|B'N)x=B'b
		self.A=np.hstack((np.eye(self.m),B.I*N))
		self.b=np.dot(B.I,self.b)
		# 相应的重新定义JN，JB
		self.x_index=list(JB)+list(JN)
		JB=[i for i in range(self.m)]
		JN=[i for i in range(self.m,self.n)]
		# 在目标函数中，使用非基变量替代基变量
		for i in range(self.m):
			_change=np.zeros((1,self.n+1))
			_change[0,:self.n]=self.A[i,:]
			_change[0,-1]=-self.b[i,0]
			self.z-=_change*self.z[0,i]
		return JB,JN

	def _set_JB_JN(self,JB,JN):
		self.JB=JB
		self.JN=JN

	def is_best(self):
		best,inf_solution=True,False
		for i in self.JN:
			sigma=self.z[0,i]
			if sigma>0:
				best=False
			elif sigma==0:
				inf_solution=True
		return best,inf_solution

	def get_inVar(self):
		greatest_sigma=0
		for i in self.JN:
			sigma=self.z[0,i]
			if greatest_sigma<sigma:
				greatest_sigma=sigma
				inVar=i
		return inVar

	def get_outVar(self,inVar):
		min_ratio=self.b[self.JB[0],0]/self.A[self.JB[0],inVar]
		outVar=self.JB[0]
		flag=False
		for i in self.JB[1:]:
			k=self.A[i,inVar]
			if k>0:
				flag=True
				_tmp=self.b[i,0]/k
				if _tmp<min_ratio:
					min_ratio=_tmp
					outVar=i
		if flag==False:
			return None
		return outVar

	def in_and_out(self,inVar,outVar):
		self.A[:,[inVar,outVar]]=self.A[:,[outVar,inVar]]
		self.x_index[outVar],self.x_index[inVar]=self.x_index[inVar],self.x_index[outVar]
		self.z[0,inVar],self.z[0,outVar]=self.z[0,outVar],self.z[0,inVar]
		B=np.hstack([self.A[:,i] for i in self.JB])
		N=np.hstack([self.A[:,i] for i in self.JN])
		# 将B转化为单位矩阵
		# 即(B|N)x=b -> (I|B'N)x=B'b
		self.A=np.hstack((np.eye(self.m),B.I*N))
		self.b=np.dot(B.I,self.b)
		# 在目标函数中，使用非基变量替代基变量
		
		for i in range(self.m):
			_change=np.zeros((1,self.n+1))
			_change[0,:self.n]=self.A[i,:]
			_change[0,-1]=-self.b[i,0]
			self.z-=_change*self.z[0,i]
	def getX(self):
		x=[0]*self.n
		for i in self.JB:
			x[self.x_index[i]]=self.b[i,0]
		return x
		
class Simplex:
	def __init__(self):
		self.solution=Solution()

		# 0 正常，尚未得到最优解，继续迭代
		# 1 无解，无界解
		# 2 达到最优解
		# 3 问题有无数个最优解
		self.status=0

	def set_para(self,A,b,z):
		# A,b,z 需以矩阵的形式输入
		self.solution.set_para(A,b,z)
	
	def output_result(self):
		self._main()
		if self.status==1:
			print("此问题无界")
		elif self.status==2:
			print("此问题有一个最优解")
		elif self.status==3:
			print("此问题有无穷多个最优解")
	
	def _main(self):
		# 获得初始可行解
		self._get_init_solution()
		if self.status==1:
			return 

		while True:
			print ("--------------------")
			print ("z:",self.solution.z[0,-1])
			print ("x:",self.solution.getX())
			# 最优性检验
			self._is_best()
			if self.status in (2,3):
				return 

			# 换入换出
			self._mainloop()
			if self.status in (1,2):
				return

	def _get_init_solution(self):
		if self.solution.get_init_solution():
			self.status=0
		else:
			self.status=1

	def _is_best(self):
		best,inf_solution=self.solution.is_best()
		if best==True and inf_solution==False:
			self.status=2
		elif best==True and inf_solution==True:
			self.status=3
		else:
			self.status=0

	def _mainloop(self):
		inVar =self.solution.get_inVar()
		outVar=self.solution.get_outVar(inVar)
		# 未找到换出基变量，此问题有无界解
		if outVar==None:
			self.status=1
			return
		self.solution.in_and_out(inVar,outVar)

if __name__=="__main__":
	s=Simplex()
	A=np.matrix([[30,20,1,0,0],
				 [ 5, 1,0,1,0],
				 [ 1, 0,0,0,1]])

	b=np.matrix([[160,15,4]]).T
	#  			  sigma,...,z0
	z=np.matrix([[5,2,0,0,0, 0]])
	s.set_para(A,b,z)
	s.output_result()

单纯形表法及其python实现

　　上面了说了一下单纯形的原理和程序实现，下面我们看一下单纯形表。

表格格式

　　设线性规划问题

　　　　$max z=cx\ s.t.
left{ egin{gathered} 
Ax=b hfill \ 
xgeq 0 hfill \ 
end{gathered} ight.$

　　其中，$A=egin{bmatrix}1&&&&a_{1,m+1}&cdots&a_{1,n}\&1&&&a_{2,m+1}&cdots&a_{2,n}\&&ddots&&&&\&&&1&a_{m,m+1}&cdots&a_{m,n} end{bmatrix}$,$c=egin{bmatrix}c_1&c_2&cdots &c_nend{bmatrix}$，

　　　　$x=egin{bmatrix}x_1&x_2&cdots&x_nend{bmatrix}^T$，$B=egin{bmatrix}b_1&b_2&cdots&b_mend{bmatrix}^T$。

　　我们建立如下表格

　　

　　其中，

• (1) 第1行是价值系数，标出了决策变量$x_j$的价值系数$c_j(j=1,2,cdots,n)$。
• (2) 第2行是标识行，标出表中主体各列的含义。
• (3) 最后1行是检验数行，除$(-z_0)$是表示当前解的目标函数$z$的负值外，其余各元素均为对应决策变量$x_j$的检验数$sigma_j(j=0,1,cdots,n)$（基变量的检验数为0）。
• (4) 第1列是$c_B$列，标出基变量的价值系数。
• (5) 第2列是$x_B$列，标出当前基变量的名称。
• (6) 第3列是右端项，前$m$个元素是当前基本可行解的基变量的取值。第$(m+1)$个元素是$(-z)$的取值。
• (7) 其余个列标出了约束方程组中决策变量$x_j$的系数列向量$p_j(j=0,1,cdots,n)$（因为$z$永远不会被换出，因此表中省去了$z$的系数列）。
• (8) 最后一列留做最小比值准则求各个比值时填数据用。(最小比值准则，即选择换出基变量的准则)。

具体步骤

• (1) 检验当前基本可行解是否为最优解？

　　观察单纯形表的检验数行，若所有的$sigma_jleq 0$，则停止计算，已得到最优解，否则进行下一步。

• (2) 检验是否为无界解？

　　在$sigma_j>0(j∈J_N)$，若有一个$sigma_{m+t}>0$，而在单纯形表中$sigma_{m+t}$所在列的其他元素，即$p_{m+t}$列的所有分量$a_{i,m+t}leq 0(i=1,2,cdots,m)$，则该问题无最优解，停止计算，否则进入下一步。

• (3) 选择换入基变量。

　　由换入基变量选择准则：$max{sigma_j}(sigma_j>0,j∈J_N) = sigma_{m+t}$，选择$x_{m+t}$为换入基变量，相应$p_{m+t}$为换入基向量。称$p_{m+t}$所在列为主列。

• (4) 选择换出基变量。

　　由换出基变量的最小比值准则：

　　　　$ heta=min_i{frac{b_i}{a_{i,m+t}}|a_{i,m+t}>0}=frac{b_l}{a_{l,m+t}}$

　　则称第$l$行为主行，与主行所对应的基变量$x_l$为换出基变量。

　　在求最小比值时，可将每一个比值$frac{b_i}{a_{i,m+t}}(a_{i,m+t}>0)$纪在单纯形表的最后一列($ heta$列)的对应位置上，然后从中选出最小值。

• (5) 基变换。

　　将可行基由$(p_1,cdots,p_l,cdots,p_m)$变换为$(p_1,cdots,p_{l-1},p_{m+t},p_{l+1},cdots,p_m)$，且将主列$p_{m+t}$化为单位列向量$e_l$即

　　　　$p_{m+t}=egin{bmatrix}a_{1,m+t}\a_{2,m+t}\ vdots \a_{m,m+t} end{bmatrix} underset{Longrightarrow}{ 化为 } p_l=egin{bmatrix}0\ vdots\1\ vdots \0end{bmatrix}leftarrow 第l个分量$

• (6) 回到(1)，对新解做最优性检验。

　　看起来，是不是特别复杂，有没有感觉单纯表很麻烦。其实，单纯形表明确了同时大大简化了，单纯形的整个计算过程。下面以一个例子进行分析。

例子

　　仍然使用最开始的例子，即某问题的数学模型标准形式如下：

$max{z}=5x_1+2x_2+0x_3+0x_4+0x_5$;

$s.t.left { egin{gathered} 
30x_1+20x_2+x3=160 hfill \ 
5x_1+x_2+x_4=15 hfill \ 
x_1+x_5=4 hfill \ 
x_jgeq0,j=1,2,3,4,5 hfill \ 
end{gathered} ight.$

解　　首先，我们建立初始的单纯形表，如下

• 建表　　初始可行基$B^{(0)}=(p_3,p_4,p_5)$，基变量为$x_3,x_4,x_5$，非基变量为$x_1,x_2$。
• 判优　　由于$sigma_1=5>0,sigma_2=2>0$，故当前解不是最优解；同时判断得本题也不是无最优解问题。
• 选择换入换出基变量　　由于$sigma_1>sigma_2$，所以选择$x_1$为换入基变量。根据最小比值准则：$frac{160}{30},frac{15}{5},frac{4}{1}$，可见最小比值为3，选择$x_4$为换出基变量，$a_{2,1}$为主元素。
• 换入换出　　使用初等行变换，将主元素$[5]$化为1，将主列$p_1$的其他各元素化为0，得到下表。（具体来说，首先将主元素那一行(主行)的所有元素除以5，然后主行上一行减去主行*30，主行下一行减去主行*1，此外检验数那一行减去主行*5）

 　

• 判优　　由于$sigma_2=1>0$，故当前解不是最优解；同时判断得本题也不是无最优解问题。
• 选择换入换出基变量　　由于$sigma_2=1>0$，所以选择$x_2$为换入基变量。根据最小比值准则：$frac{70}{14},frac{3}{1/5}$，可见最小比值为3，选择$x_3$为换出基变量，$a_{1,2}$为主元素。
• 换入换出　　使用初等行变换，将主元素$[14]$化为1，将主列$p_2$的其他各元素化为0，得到下表。

　　

• 判优　　由于所有的检验数都小于0，故当前解是最优解。因此，$z$的最大值为20，此时$x_1=2,x_2=5,x_5=2$。

python实现

import numpy as np
# 实体类　Table
# 控制类　Simplex

class Table:
	def __init__(self):
		pass

	def set_para(self,A,b,c,base,z0):
		"""
		输入LP必须已经化为标准形式
		"""
		self.A=A
		self.b=b
		self.c=c 
		self.z0=z0
		self.base=base
		self.m,self.n=self.A.shape


	def build(self):
		self.table=np.zeros((self.m+1,self.n+1))
		self.table[:-1,:1]=self.b.T
		self.table[-1 ,0]=self.z0
		self.table[:-1,1:]=self.A
		self.table[-1, 1:]=c
		self.baseVar=base
		

	def is_best(self):
		for sigma_index in range(self.n):
			if sigma_index not in self.baseVar:
				sigma=self.table[-1,1+sigma_index]
				if sigma>0:
					return False
		return True

	def is_no_solution(self):
		for sigma_index in range(self.n):
			if sigma_index not in self.baseVar:
				sigma=self.table[-1,1+sigma_index]
				if sigma>0:
					no_solution_flag=True
					for a in self.table[:-1,1+sigma_index]:
						if a>0:
							no_solution_flag=False
					if no_solution_flag==True:
						return True
		return False

	def get_inVar(self):
		max_sigma=0
		inVar=None
		for sigma_index in range(self.n):
			if sigma_index not in self.baseVar:
				sigma=self.table[-1,1+sigma_index]
				if sigma>max_sigma:
					max_sigma=sigma
					inVar=sigma_index
		return inVar

	def get_outVar(self,inVar):
		rates=[]
		for nobaseVar in range(self.m):
			a=self.table[nobaseVar,inVar+1]
			b=self.table[nobaseVar,     0 ]
			if a>0:
				rate=b/a
				rates.append((rate,nobaseVar))
		return min(rates)[1]
		
	def in_out(self,inVar,outVar):
		a=self.table[outVar,inVar+1]
		self.table[outVar,:]/=a 
		for i in range(self.m+1):
			if i != outVar:
				self.table[i,:]-=self.table[outVar,:]*self.table[i,inVar+1]
		self.baseVar[outVar]=inVar

	def show(self):
		print ('基变量/取值：',self.baseVar,end='/')
		print (self.table[:-1,0])
		print ("单纯形表")
		for i in range(self.m+1):
			for j in range(self.n+1):
				print ('%6.2f'%self.table[i,j],end=' ')
			print ()
		print ()

class Simplex:
	def __init__(self):
		self.table=Table()

		# 0 正常，尚未得到最优解，继续迭代
		# 1 无解，无界解
		# 2 达到最优解
		self.status=0
		self.inVar=None
		self.outVar=None

	def set_para(self,A,b,c,base,z0=0):
		self.table.set_para(A,b,c,base,z0)

	def output_result(self):
		self._main()
		if self.status==1:
			print("此问题无界")
		elif self.status==2:
			print("此问题有一个最优解")
		elif self.status==3:
			print("此问题有无穷多个最优解")

	def _main(self):
		self._build_table()
		while 1:
			self.table.show()
			if self._is_best() or self._is_no_solution():
				return

			self._get_inVar()
			self._get_outVar()
			self._in_out()

	def _build_table(self):
		self.table.build()

	def _is_best(self):
		if self.table.is_best():
			self.status=2
			return True
		return False

	def _is_no_solution(self):
		if self.table.is_no_solution():
			self.status=1
			return True
		return False

	def _get_inVar(self):
		self.inVar=self.table.get_inVar()

	def _get_outVar(self):
		self.outVar=self.table.get_outVar(self.inVar)

	def _in_out(self):
		self.table.in_out(self.inVar,self.outVar)

if __name__=="__main__":
	s=Simplex()
	A=np.matrix([[30,20,1,0,0],
				 [ 5, 1,0,1,0],
				 [ 1, 0,0,0,1]])

	b=np.matrix([[160,15,4]])
	c=np.matrix([[5,2,0,0,0]])
	base=[2,3,4]
	s.set_para(A,b,c,base)
	s.output_result()

 人工变量及其处理办法

　　需要注意的是，上面的单纯形表格法的基本假设是，输入的A中包含一个单位矩阵(也就是说，初始可行基为单位矩阵)。往往，我们需要添加人工变量来达到这个条件。例如

　　$A=egin{bmatrix}30&20&1&1&2\5&1&2&1&3\1&2&3&4&1\ end{bmatrix}$，对于A这个矩阵，如果我们通过初等行变换将其转化为一个包含单位矩阵的矩阵，这是比较麻烦的；为了解决这个问题，我们通过添加人工变量$x_6,x_7,x_8$的方法，将A变成$egin{bmatrix}30&20&1&1&2&1&&\5&1&2&1&3&&1&\1&2&3&4&1&&1\ end{bmatrix}$，这样A中就含有一个单位阵。我们可以以$x_6,x_7,x_8$为基变量，开始迭代。但是添加了人工变量后的问题与原问题已经不等价了。只有当最优解中，人工变量都取0值时，才可以认为两个问题的最优解是相当的。

　　为了解决这个问题，两种方法被引入。

　1.大M法

　　假设原目标函数为$max{z}=sum{c_jx_j}$，添加人工变量之后，我们将目标函数修改为$max{z}=sum{c_jx_j}-Mx_{n+1}-Mx_{n+2}-cdots-Mx_{n+m}$。其中，$M$是个很大的正数。为了对目标函数实现最大化，人工变量就会被迅速换出去。

　2.两阶段法

 　　当线性规划问题添加了人工变量之后，我们可以将问题拆成两个问题。

　　第一阶段求解第一个线性规划问题，目标函数为$min{w}=sum{x_{n+i}}$，即目标是对所有人工变量之和求最小。需要说明的是，当求得的最优解中，至少有一个人工变量不为0值，说明原问题无可行解，不需要继续进行第二阶段的计算。

　　第二阶段，目标函数为添加人工变量之前的目标函数。

1