在运筹学中，运输问题 (Transportation Problem, TP) 的表上作业法通常包含3个步骤：
- 在已知运输表的基础上寻找初始基可行解。通常采用下列3种方法：
1）西北角法（简便但求解效率低）；
2）最小元素法（整体效率适中）；
3）伏格尔（Vogel）法（效率最高，但求解过程繁琐）。 - 通过计算初始基可行解的检验数来判断当前解的最优性。若已为最优解，则直接输出当前解；若非最优解，则转入下一步。判断最优性的方法主要有2种：
1）闭回路法；
2）位势法（对偶变量法）。 - 对当前解进行调整，直到达到最优解。调整时通常采用闭回路法。
本期谈谈第一步：寻找初始基可行解。这里主要讲解伏格尔（Vogel）法的原理和Python实现。
Vogel法简介
伏格尔 (Vogel) 法又称差值法。这种方法考虑到某产地的产品如不能按最小运费就近供应，就考虑次小运费，这里就有一个差额，称为罚数。罚数越大，说明不能按最小运费调运时，运费增加越多。因而对罚数最大处，就应当采用最小运费调运。
Vogel法一般能得到一个比用西北角法和最小元素法所得的初始基本可行解更好的初始基本可行解。使用Vogel法时，先计算出各行各列中的罚数（最小的成本系数cij和次小的成本系数cij之间的差的绝对值），在具有最大罚数的行/列中，选择具有最小的cij的方格来填入基变量值。这样就可以避免将运量分配到该行/列具有次小的成本cij的方格中（如果填入本格会使成本大幅增加），以保证有较小的目标函数值。所以，伏格尔法的基本步骤如下：
- 算出各行各列的罚数，并标出最大罚数（若几个差额同为最大，则可任取其一）；
- 在罚数最大的行或列中的最小元素处填上尽可能大的数（由该行该列对应的资源向量a和需求向量b决定，选择较小的那个值）；
- 划去满足条件（资源已耗尽/需求已满足）的行/列，并修改对应的向量a和b；
- 对未划去的行列重复以上步骤，直到得到一个初始解。
由此可见，Vogel法同最小元素法除在确定供求关系的原则上不同外，其余步骤相同。Vogel法给出的初始基可行解比用最小元素法更接近最优解。
Vogel法应用实例
Example:
假设TP的运输表如下：
使用Vogel法找到的初始基可行解如下（空格位置xij=0）：
在没有退化的情况下，初始基可行解含有 (m+n-1) 个非零变量（m，n分别为成本矩阵的行数和列数，本例中m=3，n=4，基可行解含有6个非零变量）；而在有些运输问题中，初始基可行解的非零变量个数可能少于 (m+n-1)，这是因为可行解中有部分基变量等于0。这种情形称为解的退化。总体而言，运输问题基可行解中非零变量个数不应大于 (m+n-1)。
Vogel法的Python语句
将上例的运输表保存为EXCEL文件，文件名为TP_PPT_Sample1.xlsx：
然后在Python中执行以下代码：
#运输问题求解：使用Vogel逼近法寻找初始基本可行解
import numpy as np
import copy
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)
x=np.zeros(c.shape)
M=pow(10,9)
for factor in c.reshape(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.remove(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.remove(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)
pun.extend(c_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):
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.array(np.nonzero(x)).shape[1]
if varnum!=m+n-1:
print('【注意：问题含有退化解】')
return (cost,x)
path=r'C:\Users\spurs\Desktop\MCM_ICM\Data files\TP_PPT_Sample1.xlsx'
mat=pd.read_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])
运行最终结果：
c:
[[ 3. 11. 3. 10.]
[ 1. 9. 2. 8.]
[ 7. 4. 10. 5.]]
行罚数： [0.0, 1.0, 1.0]
列罚数： [2.0, 5.0, 1.0, 3.0]
罚数向量： [0.0, 1.0, 1.0, 2.0, 5.0, 1.0, 3.0]
最大罚数： 5.0 元素序号： 5
对第 2 列进行操作：
[11. 9. 4.]
最小成本所在行索引： 2
本次迭代后的x矩阵：
[[0. 0. 0. 0.]
[0. 0. 0. 0.]
[0. 6. 0. 0.]]
a: [7. 4. 3.]
b: [3. 0. 5. 6.]
c:
[[3.e+00 1.e+09 3.e+00 1.e+01]
[1.e+00 1.e+09 2.e+00 8.e+00]
[7.e+00 1.e+09 1.e+01 5.e+00]]
【迭代未完成】
------------------------------------------------------------
c:
[[3.e+00 1.e+09 3.e+00 1.e+01]
[1.e+00 1.e+09 2.e+00 8.e+00]
[7.e+00 1.e+09 1.e+01 5.e+00]]
行罚数： [0.0, 1.0, 2.0]
列罚数： [2.0, 0.0, 1.0, 3.0]
罚数向量： [0.0, 1.0, 2.0, 2.0, 0.0, 1.0, 3.0]
最大罚数： 3.0 元素序号： 7
对第 4 列进行操作：
[10. 8. 5.]
最小成本所在行索引： 2
本次迭代后的x矩阵：
[[0. 0. 0. 0.]
[0. 0. 0. 0.]
[0. 6. 0. 3.]]
a: [7. 4. 0.]
b: [3. 0. 5. 3.]
c:
[[3.e+00 1.e+09 3.e+00 1.e+01]
[1.e+00 1.e+09 2.e+00 8.e+00]
[1.e+09 1.e+09 1.e+09 1.e+09]]
【迭代未完成】
------------------------------------------------------------
c:
[[3.e+00 1.e+09 3.e+00 1.e+01]
[1.e+00 1.e+09 2.e+00 8.e+00]
[1.e+09 1.e+09 1.e+09 1.e+09]]
行罚数： [0.0, 1.0, 0.0]
列罚数： [2.0, 0.0, 1.0, 2.0]
罚数向量： [0.0, 1.0, 0.0, 2.0, 0.0, 1.0, 2.0]
最大罚数： 2.0 元素序号： 4
对第 1 列进行操作：
[3.e+00 1.e+00 1.e+09]
最小成本所在行索引： 1
本次迭代后的x矩阵：
[[0. 0. 0. 0.]
[3. 0. 0. 0.]
[0. 6. 0. 3.]]
a: [7. 1. 0.]
b: [0. 0. 5. 3.]
c:
[[1.e+09 1.e+09 3.e+00 1.e+01]
[1.e+09 1.e+09 2.e+00 8.e+00]
[1.e+09 1.e+09 1.e+09 1.e+09]]
【迭代未完成】
------------------------------------------------------------
c:
[[1.e+09 1.e+09 3.e+00 1.e+01]
[1.e+09 1.e+09 2.e+00 8.e+00]
[1.e+09 1.e+09 1.e+09 1.e+09]]
行罚数： [7.0, 6.0, 0.0]
列罚数： [0.0, 0.0, 1.0, 2.0]
罚数向量： [7.0, 6.0, 0.0, 0.0, 0.0, 1.0, 2.0]
最大罚数： 7.0 元素序号： 1
对第 1 行进行操作：
[1.e+09 1.e+09 3.e+00 1.e+01]
最小成本所在列索引： 2
本次迭代后的x矩阵：
[[0. 0. 5. 0.]
[3. 0. 0. 0.]
[0. 6. 0. 3.]]
a: [2. 1. 0.]
b: [0. 0. 0. 3.]
c:
[[1.e+09 1.e+09 1.e+09 1.e+01]
[1.e+09 1.e+09 1.e+09 8.e+00]
[1.e+09 1.e+09 1.e+09 1.e+09]]
【迭代未完成】
------------------------------------------------------------
c:
[[1.e+09 1.e+09 1.e+09 1.e+01]
[1.e+09 1.e+09 1.e+09 8.e+00]
[1.e+09 1.e+09 1.e+09 1.e+09]]
行罚数： [999999990.0, 999999992.0, 0.0]
列罚数： [0.0, 0.0, 0.0, 2.0]
罚数向量： [999999990.0, 999999992.0, 0.0, 0.0, 0.0, 0.0, 2.0]
最大罚数： 999999992.0 元素序号： 2
对第 2 行进行操作：
[1.e+09 1.e+09 1.e+09 8.e+00]
最小成本所在列索引： 3
本次迭代后的x矩阵：
[[0. 0. 5. 0.]
[3. 0. 0. 1.]
[0. 6. 0. 3.]]
a: [2. 0. 0.]
b: [0. 0. 0. 2.]
c:
[[1.e+09 1.e+09 1.e+09 1.e+01]
[1.e+09 1.e+09 1.e+09 1.e+09]
[1.e+09 1.e+09 1.e+09 1.e+09]]
【迭代未完成】
------------------------------------------------------------
c:
[[1.e+09 1.e+09 1.e+09 1.e+01]
[1.e+09 1.e+09 1.e+09 1.e+09]
[1.e+09 1.e+09 1.e+09 1.e+09]]
行罚数： [999999990.0, 0.0, 0.0]
列罚数： [0.0, 0.0, 0.0, 999999990.0]
罚数向量： [999999990.0, 0.0, 0.0, 0.0, 0.0, 0.0, 999999990.0]
最大罚数： 999999990.0 元素序号： 1
对第 1 行进行操作：
[1.e+09 1.e+09 1.e+09 1.e+01]
最小成本所在列索引： 3
本次迭代后的x矩阵：
[[0. 0. 5. 2.]
[3. 0. 0. 1.]
[0. 6. 0. 3.]]
a: [0. 0. 0.]
b: [0. 0. 0. 0.]
c:
[[1.e+09 1.e+09 1.e+09 1.e+09]
[1.e+09 1.e+09 1.e+09 1.e+09]
[1.e+09 1.e+09 1.e+09 1.e+09]]
【迭代完成】
------------------------------------------------------------
>>>供求平衡<<<
>>>初始基本可行解x*：
[[0. 0. 5. 2.]
[3. 0. 0. 1.]
[0. 6. 0. 3.]]
>>>当前总成本： 85.0
下面举一个退化解的例子：
Example:
一运输问题的运输表如下：
※本例其实是指派问题，但是也可以用TP的方法解。关于指派问题的专门解法，笔者（可能）会在将来进行介绍。
将以上表格保存为TP_Text_5.9.xlsx，在Python中求解。程序运行的最终结果为：
>>>供求平衡<<<
>>>初始基本可行解x*：
[[1. 0. 0. 0. 0.]
[0. 0. 1. 0. 0.]
[0. 0. 0. 1. 0.]
[0. 0. 0. 0. 1.]
[0. 1. 0. 0. 0.]]
>>>当前总成本： 0.8999999999999999
【注意：问题含有退化解】
最终得到的初始基可行解只含有5个非零变量，而 (m+n-1)=9，可见发生了退化。
函数使用方法：
TP_vogel (var): 变量var可以是以下3种类型：
- numpy.ndarray：将TP运输表保存为numpy.ndarray类型的矩阵；
- tuple：将 (成本矩阵c, 供给向量a, 需求向量b) 以 (c,a,b) 保存为元组类型；
- list：将 (成本矩阵c, 供给向量a, 需求向量b) 以 [c,a,b] 保存为列表类型。
注：本函数返回的是以成本矩阵和初始基可行解矩阵构成的元组，即 (c,x)。
下一期将介绍表上作业法的第二步：判断可行解最优性的位势法及其Python实现：
感谢阅读，欢迎指正。
扫一扫
8066
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
