【Python 实现】运输问题的表上作业法（⼆）：利⽤位势法判断当前解的最优性
上期谈到⽤Python实现求解运输问题 (Transportation Problem, TP) 的表上作业法的第⼀步——利⽤Vogel法寻初始基可⾏解：这期来讲讲到初始基可⾏解之后怎样判断当前解是否是最优解。如果当前解已达到最优，那么⽆需再进⾏操作；如果当前解⾮最优，那么还要对当前解进⾏调整以达到最优。调整解的操作放到下⼀期再讲，本期先谈谈怎样判断最优性。
⽂章⽬录
通常判断TP问题解的最优性有两种⽅法：闭回路法和位势法。考虑到两种⽅法代码的复杂程度，我们重点讲讲更容易实现的位势法。位势法（Potentials）⼜称对偶变量法，是利⽤线性规划的对偶原理计算解的检验数，从⽽通过检验数判断最优性的⽅法。关于对偶原理，请参考运筹学相关书籍，这⾥不进⾏阐述，⽽只对位势法原理进⾏简单介绍。
位势法（对偶变量法）
位势法应⽤实例
下⾯通过⼀个运输问题的例题进⼀步说明位势法原理。
Example：
已知该TP问题的运输（成本）表如下（将中间3x4的成本表记为表1）：
1.
通过最⼩元素法到⼀组初始基可⾏解：
2. 把初始基可⾏解中⾮零元素对应位置的成本系数剥离出来，并添加上位势（图中u -u 、v -v
），记为表2：
3. 构建位势⽅程组（需要注意的是，TP问题的可⾏解中⾮零变量个数为m+n-1，其中m、n分别是成本矩阵的⾏数和列数（本例中m=3,n=4），所以位势⽅程组应该含有m+n-1个⽅程；但是位势的个数为m+n，即⽅程组有m+n个未知量，多于⽅程的个数，可知
⽅程组是解不出来的，因此我们随便令某个位势=0，即可求解这个⽅程组，并且此操作对最后检验数⽆影响）：
4. 求出位势后，将u +v
填⼊成本表中对应解变量为0的位置，记为表3：
5. 表1-表3（σ=c -(u +v
)，对应位置元素相减），得到检验数表σ：
6. 判断最优性：由于本例是极⼩化问题，如果检验数表中存在σ<0，说明当前解未达到最优，需要对解进⾏进⼀步的调整；否则，当前解达到最优（本例中当前初始基可⾏解并⾮最优解，因为σ=-1<0）。如果采⽤的检验数σ=(u +v )-c ，则全部反向判断即可。位势法的Python 语句
同样使⽤上⾯的例题，但在寻初始基可⾏解时运⽤Vogel法，编制Python程序进⾏求解。
⾸先将运输表写⼊EXCEL表格并保存，⽂件名为TP_PPT_Sample1.xlsx：
然后执⾏以下代码：
#Vogel 法寻初始基可⾏解+位势法判断解的最优性
import  numpy as  np
import  copy
1314i j ij ij i j ij 24ij i j ij
import pandas as pd
def TP_split_matrix(mat):
c=mat[:-1,:-1]
a=mat[:-1,-1]
b=mat[-1,:-1]
return(c,a,b)
def TP_vogel(var):#Vogel法代码，变量var可以是以numpy.ndarray保存的运输表，或以tuple或list保存的(成本矩阵，供给向量，需求向量) import numpy
typevar1=type(var)==numpy.ndarray
typevar2=type(var)==tuple
typevar3=type(var)==list
if typevar1==False and typevar2==False and typevar3==False:
print('>>>⾮法变量<<<')
(cost,x)=(None,None)
else:
if typevar1==True:
[c,a,b]=TP_split_matrix(var)
elif typevar2==True or typevar3==True:
[c,a,b]=var
cost=copy.deepcopy(c)
s(c.shape)
M=pow(10,9)
for factor shape(1,-1)[0]:
while int(factor)!=M:
if np.all(c==M):
break
else:
print('c:\n',c)
#获取⾏/列最⼩值数组
row_mini1=[]
row_mini2=[]
for row in range(c.shape[0]):
Row=list(c[row,:])
row_min=min(Row)
row_mini1.append(row_min)
row_2nd_min=min(Row)
row_mini2.append(row_2nd_min)
#print(row_mini1,'\n',row_mini2)
r_pun=[row_mini2[i]-row_mini1[i]for i in range(len(row_mini1))]
print('⾏罚数：',r_pun)
#计算列罚数
col_mini1=[]
col_mini2=[]
for col in range(c.shape[1]):
Col=list(c[:,col])
col_min=min(Col)
col_mini1.append(col_min)
col_2nd_min=min(Col)
col_mini2.append(col_2nd_min)
c_pun=[col_mini2[i]-col_mini1[i]for i in range(len(col_mini1))]
print('列罚数：',c_pun)
pun=copy.deepcopy(r_pun)
print('罚数向量：',pun)
max_pun=max(pun)
max_pun_index=pun.index(max(pun))
max_pun_num=max_pun_index+1
print('最⼤罚数：',max_pun,'元素序号：',max_pun_num)
if max_pun_num<=len(r_pun):
row_num=max_pun_num
print('对第',row_num,'⾏进⾏操作：')
row_index=row_num-1
catch_row=c[row_index,:]
print(catch_row)
min_cost_colindex=int(np.argwhere(catch_row==min(catch_row)))
print('最⼩成本所在列索引：',min_cost_colindex)
if a[row_index]<=b[min_cost_colindex]:
x[row_index,min_cost_colindex]=a[row_index]
c1=copy.deepcopy(c)
c1[row_index,:]=[M]*c1.shape[1]
b[min_cost_colindex]-=a[row_index]
a[row_index]-=a[row_index]
else:
x[row_index,min_cost_colindex]=b[min_cost_colindex]
c1=copy.deepcopy(c)
c1[:,min_cost_colindex]=[M]*c1.shape[0]
a[row_index]-=b[min_cost_colindex]
b[min_cost_colindex]-=b[min_cost_colindex]
else:
col_num=max_pun_num-len(r_pun)
col_index=col_num-1
print('对第',col_num,'列进⾏操作：')
catch_col=c[:,col_index]
print(catch_col)
#寻最⼤罚数所在⾏/列的最⼩成本系数
min_cost_rowindex=int(np.argwhere(catch_col==min(catch_col)))
print('最⼩成本所在⾏索引：',min_cost_rowindex)
#计算将该位置应填⼊x矩阵的数值（a,b中较⼩值）
if a[min_cost_rowindex]<=b[col_index]:
x[min_cost_rowindex,col_index]=a[min_cost_rowindex]
c1=copy.deepcopy(c)
c1[min_cost_rowindex,:]=[M]*c1.shape[1]
b[col_index]-=a[min_cost_rowindex]
a[min_cost_rowindex]-=a[min_cost_rowindex]
else:
x[min_cost_rowindex,col_index]=b[col_index]
#填⼊后删除已满⾜/耗尽资源系数的⾏/列，得到剩余的成本矩阵，并改写资源系数                            c1=copy.deepcopy(c)
c1[:,col_index]=[M]*c1.shape[0]
a[min_cost_rowindex]-=b[col_index]
b[col_index]-=b[col_index]
c=c1
print('本次迭代后的x矩阵：\n',x)
print('a:',a)
print('b:',b)
print('c:\n',c)
if np.all(c==M):
print('【迭代完成】')
print('-'*60)
else:
print('【迭代未完成】')
print('-'*60)
total_cost=np.sum(np.multiply(x,cost))
if np.all(a==0):能运行python的软件
if np.all(b==0):
print('>>>供求平衡<<<')
else:
print('>>>供不应求，需求⽅有余量<<<')
elif np.all(b==0):
print('>>>供⼤于求，供给⽅有余量<<<')
else:
print('>>>⽆法到初始基可⾏解<<<')
print('>>>初始基本可⾏解x*：\n',x)
print('>>>当前总成本：',total_cost)
[m,n]=x.shape
varnum=np.o(x)).shape[1]
if varnum!=m+n-1:
print('【注意：问题含有退化解】')
return(cost,x)
def create_c_nonzeros(c,x):
import numpy as np
import copy
nonzeros=copy.deepcopy(x)
for i in range(x.shape[0]):
for j in range(x.shape[1]):
if x[i,j]!=0:
nonzeros[i,j]=1
#print(nonzeros)
c_nonzeros=np.multiply(c,nonzeros)
return c_nonzeros
def if_dsquare(a,b):
print('a:',a.shape,'\n','b:',b.shape)
correct='>>>位势⽅程组可解<<<'
error='>>>位势⽅程组不可解，请检查基变量个数是否等于(m+n-1)及位势未知量个数是否等于(m+n)<<<' if len(a.shape)==2:
if len(b.shape)==2:
if a.shape[0]==a.shape[1]and a.shape==b.shape:
print(correct)
if_dsquare=True
else:
print(error)
if_dsquare=False
elif len(b.shape)==1and b.shape[0]!=0:
if a.shape[0]==a.shape[1]and a.shape[0]==b.shape[0]:
print(correct)
if_dsquare=True
else:
print(error)
if_dsquare=False
else:
print(error)
if_dsquare=False
elif len(a.shape)==1:
if len(b.shape)==2:
if b.shape[0]==b.shape[1]and a.shape[0]==b.shape[0]:
print('>>>位势⽅程组系数矩阵与⽅程组值向量位置错误<<<')
if_dsquare='True but not solvable'
else:
print(error)
if_dsquare=False
elif len(b.shape)==1:
print(error)
if_dsquare=False
else:
print(error)
if_dsquare=False
else:
print(error)
if_dsquare=False
return if_dsquare
def TP_potential(cost,x):
[m,n]=x.shape
varnum=np.o(x)).shape[1]
if varnum!=m+n-1:
sigma=None
print('【问题含有退化解，暂时⽆法判断最优性】')
else:
#print(c_nonzeros.shape)
c_nonzeros=create_c_nonzeros(c,x)
cc_nonzeros=np.o(c_nonzeros))
A=[]
B=[]
length=c_nonzeros.shape[0]+c_nonzeros.shape[1]
s((1,length))[0]
for i in range(cc_nonzeros.shape[1]):
s((1,length))[0]
zeros[cc_nonzeros[0,i]]=1
zeros[cc_nonzeros[1,i]+c_nonzeros.shape[0]]=1
A.append(zeros)
B.append(c_nonzeros[cc_nonzeros[0,i],cc_nonzeros[1,i]])
l=[1]
for j in range(length-1):
l.append(0)#补充⼀个x1=0的⽅程以满⾜求解条件
A.append(l)
B.append(0)
#print(A)
#print(B)
A=np.array(A)
B=np.array(B)
judge=if_dsquare(A,B)
if judge==True:
sol=np.linalg.solve(A,B)#求解条件：A的⾏数（⽅程个数)=A的列数(变量个数）=B的个数（⽅程结果个数）才能解#print(sol)
var=[]#创建位势名称数组
for i in range(c_nonzeros.shape[0]):
var.append('u'+str(i+1))
for j in range(c_nonzeros.shape[1]):
var.append('v'+str(j+1))
#print(var)
solset=dict(zip(var,sol))
print('>>>当前位势:\n',solset)
u=[]
v=[]
[m,n]=c_nonzeros.shape
zero_anspose(np.argwhere(c_nonzeros==0))
for i in range(m):
u.append(sol[i])
for j in range(n):
v.append(sol[j+m])
for k in range(zero_places.shape[1]):
c_nonzeros[zero_places[0,k],zero_places[1,k]]=u[zero_places[0,k]]+v[zero_places[1,k]] #print(c_nonzeros)
sigma=cost-c_nonzeros
print('>>>检验表σ：\n',sigma)
if np.all(sigma>=0):
print('>>>TP已达到最优解<<<')
else:
print('>>>TP未达到最优解<<<')
else:
sigma=None
print('>>>位势⽅程组不可解<<<')
return sigma
path=r'C:\Users\spurs\Desktop\MCM_ICM\Data files\TP_PPT_Sample1.xlsx'
ad_excel(path,header=None).values
#c=np.array([[3,11,3,10],[1,9,2,8],[7,4,10,5]])
#a=np.array([7,4,9])
#b=np.array([3,6,5,6])
[c,x]=TP_vogel(mat)
#[c,x]=TP_vogel([c,a,b])
sigma=TP_potential(c,x)

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。