本文介绍了如何使用Python实现运输问题(TP)算法，包括North-West-Corner和Minimum-Cell-Cost的初始化方法，以及SteppingStone和Modified Distribution的迭代算法。通过实例展示了如何处理不平衡问题并提供完整性和退化处理。最后在main.py中演示了完整的算法应用和成本计算。
基于前面几篇文章，我用python把TP问题算法实现一下
首先定义TPAlgorithm类，它的构造函数中输入supply、demand和cost参数，可以选择用North-West-Corner或者Minimum-Cell-Cost的初始化算法，用SteppingStone或者MODI的迭代算法；handleUnbalance方法用于在实际计算前进行平衡化处理；
import numpy as np
from NorthWestCorner import NorthWestCorner
from MinCellCost import MinCellCost
from SteppingStone import SteppingStone
from ModifiedDistribution import ModifiedDistribution
class TPAlgorithm:
def __init__(self, supply, demand, cost, initMethod='NorthWestCorner', iterMethod='SteppingStone'):
'''
supply: [a_1,a_2,...,a_m]
demand: [b_1,b_2,...,b_n]
cost: [
[c_11, c_12, ... , c_1n],
[c_21, c_22, ... , c_2n],
...
[c_m1, c_m2, ... , c_mn]
]
initMethod: 'NorthWestCorner' or 'MinCellCost'
iterMethod: 'SteppingStone' or 'ModifiedDistribution'
'''
self.supply=supply
self.demand=demand
self.cost=cost
self.initMethod=initMethod
self.iterMethod=iterMethod
self.handleUnbalance()
def handleUnbalance(self):
'''
Check whether the problem is unbalanced, if so, extend the tableau
'''
sum_supply=np.sum(self.supply)
sum_demand=np.sum(self.demand)
m=self.supply.shape[0]
n=self.demand.shape[0]
if(sum_supply<sum_demand):#If total demand is larger than total supply, add a dummy row
self.supply=np.append(self.supply,[sum_demand-sum_supply])
dummy_row=np.zeros((1,n))
self.cost=np.vstack((self.cost,dummy_row))
elif(sum_supply>sum_demand):#If total demand is less than total supply, add a dummy column
self.demand=np.append(self.demand,[sum_supply-sum_demand])
dummy_column=np.zeros((m,1))
self.cost=np.hstack((self.cost,dummy_column))
def solve(self):
'''
return: [
[x_11, x_12, ... , x_1n],
[x_21, x_22, ... , x_2n],
...
[x_m1, x_m2, ... , x_mn]
]
'''
initMethodAlgorithm=None
if self.initMethod=='NorthWestCorner':
initMethodAlgorithm=NorthWestCorner(self.supply, self.demand, self.cost)
elif self.initMethod=='MinCellCost':
initMethodAlgorithm=MinCellCost(self.supply, self.demand, self.cost)
initSolution=initMethodAlgorithm.calculate()
iterAlgorithm=None
if self.iterMethod=='SteppingStone':
iterAlgorithm=SteppingStone(self.supply,self.demand,self.cost,initSolution)
if self.iterMethod=='ModifiedDistribution':
iterAlgorithm=ModifiedDistribution(self.supply,self.demand,self.cost,initSolution)
return iterAlgorithm.calculate()
def solutionCost(self,solution):
m=self.supply.shape[0]
n=self.demand.shape[0]
totalCost=0
for i in range(0,m):
for j in range(0,n):
cell_value=0 if np.isnan(solution[i,j]) else solution[i,j]
totalCost=totalCost+(self.cost[i,j])*cell_value
return totalCost
接着实现两种初始化算法，因为两种算法的接口类似，所以定义了一个统一的InitMethodInterface接口类：
import abc
import numpy as np
class InitMethodInterface(metaclass=abc.ABCMeta):
'''
The interface for northwest-corner, minimum-cell-cost, and Vogel's approximation methods
'''
def __init__(self, supply, demand, cost):
'''
supply: [a_1,a_2,...,a_m]
demand: [b_1,b_2,...,b_n]
cost: [
[c_11, c_12, ... , c_1n],
[c_21, c_22, ... , c_2n],
...
[c_m1, c_m2, ... , c_mn]
]
'''
self.supply=supply
self.demand=demand
self.cost=cost
@abc.abstractmethod
def calculate(self):
pass
def isComplete(self,solution):
demand_sums=np.nansum(solution,axis=0)
for inx,sum_Value in enumerate(demand_sums):
if sum_Value!=self.demand[inx]:
return False
supply_sums=np.nansum(solution,axis=1)
for inx,sum_Value in enumerate(supply_sums):
if sum_Value!=self.supply[inx]:
return False
return True
import numpy as np
from InitMethodInterface import InitMethodInterface
'''
Use North-West-Corner method to calculate Basic Feasible Solution
'''
class NorthWestCorner(InitMethodInterface):
def calculate(self):
'''
step 1 : Allocate as much as possible to the cell in the upper left-hand corner, subject to the
supply and demand constraints
step 2 : Allocate as much as possible to the next adjacent feasible cell
step 3 : Repeat step 2 until all rim requirements have been met.
'''
m=self.supply.shape[0]
n=self.demand.shape[0]
solution=np.full((m,n),np.nan)
solution[0,0]=min(self.supply[0],self.demand[0])
index_i=0
index_j=0
while self.isComplete(solution) is not True:
index_j=index_j+1
if np.nansum(solution,1)[index_i]==self.supply[index_i]:
index_j=index_j-1
index_i=index_i+1
temp=self.demand[index_j]
for i in range(0,index_i):
cell_value=0 if np.isnan(solution[i,index_j]) else solution[i,index_j]
temp=temp-cell_value
solution[index_i][index_j]=min(temp,self.supply[index_j])
else:
temp=self.supply[index_i]
for j in range(0,index_j):
cell_value=0 if np.isnan(solution[index_i,j]) else solution[index_i,j]
temp=temp-cell_value
solution[index_i][index_j]=min(temp,self.demand[index_i])
return solution
import numpy as np
from InitMethodInterface import InitMethodInterface
'''
Use Minimum-Cell-Cost method to calculate Basic Feasible Solution
'''
class MinCellCost(InitMethodInterface):
def minCostCell(self,solution):
'''
find the cell with minimized cost, available and feasible
'''
m=self.supply.shape[0]
n=self.demand.shape[0]
min_i=0
min_j=0
min_cost=float('inf')
for i in range(0,m):
for j in range(0,n):
if np.isnan(solution[i,j]) is False:
continue
demand_sum=np.nansum(solution,axis=0)[j]
supply_sum=np.nansum(solution,axis=1)[i]
if demand_sum==self.demand[j] or supply_sum==self.supply[i]:
continue
if min_cost==float('inf'):
min_i=i
min_j=j
min_cost=self.cost[i,j]
else:
if min_cost>self.cost[i,j]:
min_i=i
min_j=j
min_cost=self.cost[i,j]
return (min_i,min_j)
def calculate(self):
'''
step 1 : Allocate as much as possible to the feasible cell with the minimum transportation
cost, and adjust the rim requirements
step 2 : Repeat step 1 until all rim requirements have been met.
'''
m=self.supply.shape[0]
n=self.demand.shape[0]
solution=np.full((m,n),np.nan)
while self.isComplete(solution) is not True:
(index_i,index_j)=self.minCostCell(solution)
available_demand=self.demand[index_j]
for i in range(0,m):
if i==index_i:
continue
cell_value=0 if np.isnan(solution[i,index_j]) else solution[i,index_j]
available_demand=available_demand-cell_value
available_supply=self.supply[index_i]
for j in range(0,n):
if j==index_j:
continue
cell_value=0 if np.isnan(solution[index_i,j]) else solution[index_i,j]
available_supply=available_supply-cell_value
solution[index_i,index_j]=min(available_demand,available_supply)
return solution
SteppingStone类实现了Stepping-Stone算法；它有几个方法比较复杂，findLoop方法用于寻找Loop，具体实现上采用递归的方式；handleDegeneracy方法会检查是否存在退化，如果存在则持续处理直到跳出退化；calculate里实现了算法的主流程，每次循环前首先调用handleDegeneracy方法处理退化
import numpy as np
class SteppingStone:
def __init__(self, supply, demand, cost, initSolution):
self.supply=supply
self.demand=demand
self.cost=cost
self.initSolution=initSolution
def pathIsClose(self, path):
'''
Check whether the path is closed
'''
if len(path)>3 and ( path[0][0]==path[-1][0] or path[0][1]==path[-1][1]):
return True
return False
def getAvailableRowCells(self, solution, path):
lastCellInPath=path[-1]
m=solution.shape[0]
n=solution.shape[1]
numPathThisRow=0
result=[]
rowInxs=[k for k in range(lastCellInPath[1]+1,n)]
rowInxs.extend([k for k in range(lastCellInPath[1]-1,-1,-1)])
for j in rowInxs:
if (lastCellInPath[0], j) in path:
return []
if np.isnan(solution[lastCellInPath[0],j]):
continue
result.append((lastCellInPath[0], j))
return result
def getAvailableColumnCells(self, solution, path):
lastCellInPath=path[-1]
m=solution.shape[0]
n=solution.shape[1]
numPathThisColumn=0
result=[]
colInxs=[k for k in range(lastCellInPath[0]-1,-1,-1)]
colInxs.extend([k for k in range(lastCellInPath[0]+1,m)])
for i in colInxs:
if (i,lastCellInPath[1]) in path:
return []
if np.isnan(solution[i,lastCellInPath[1]]):
continue
result.append((i,lastCellInPath[1]))
return result
def findLoop(self, solution, path, nextCell):
'''
Find a loop by recursion
'''
currentPath=path
if len(path)==1 and nextCell in path:
pass
else:
currentPath.append(nextCell)
if self.pathIsClose(currentPath) is True:
return currentPath
availableRowCells=self.getAvailableRowCells(solution, currentPath)
for cell in availableRowCells:
loop=self.findLoop(solution, currentPath, cell)
if loop is not None:
return loop
availableColumnCells=self.getAvailableColumnCells(solution, currentPath)
for cell in availableColumnCells:
loop=self.findLoop(solution, currentPath, cell)
if loop is not None:
return loop
currentPath.remove(currentPath[-1])
return None
def costChangeThisLoop(self,loop):
sum=0
for inx,cell in enumerate(loop):
sum+=(1 if inx%2==0 else -1)*self.cost[cell[0],cell[1]]
return sum
def existLoop(self, solution):
m=solution.shape[0]
n=solution.shape[1]
for i in range(0,m):
for j in range(0,n):
if np.isnan(solution[i,j]):
loop=self.findLoop(solution, [(i,j)], (i,j))
if loop is not None:
return True
return False
def handleDegeneracy(self, solution):
m=solution.shape[0]
n=solution.shape[1]
while self.isDegeneracy(solution) is True:
done=False
for i in range(0,m):
if done:
break
for j in range(0,n):
if np.isnan(solution[i,j]) == False:
continue
solution[i,j]=0 # Allocate zero to this cell, check whether loop exists
if self.existLoop(solution):
done=True
break
else:
solution[i,j]=np.nan
return solution
def calculate(self):
'''
1. Determine the stepping-stone paths and cost changes for each empty cell in the tableau.
2. Allocate as much as possible to the empty cell with the greatest net decrease in cost.
3. Repeat steps 1 and 2 until all empty cells have positive cost changes that indicate an
optimal solution.
'''
solution=self.initSolution
m=solution.shape[0]
n=solution.shape[1]
index=1
while True:
self.handleDegeneracy(solution)
minChange=float('inf')
minChangeLoop=None
for i in range(0,m):
for j in range(0,n):
if np.isnan(solution[i,j]):
loop=self.findLoop(solution, [(i,j)], (i,j))
if loop is None:
continue
change=self.costChangeThisLoop(loop)
if change<minChange:
minChange=change
minChangeLoop=loop
print('######################### iteration '+str(index)+" #############################")
print('Current solution : ')
print(solution)
index+=1
if minChange>=0:
return solution
solution=self.reAllocate(solution,minChangeLoop)
def reAllocate(self, solution, loop):
maxPlusValue=float('inf')
for inx,cell in enumerate(loop):
if inx%2==0:
continue
if maxPlusValue>solution[cell[0],cell[1]]:
maxPlusValue=solution[cell[0],cell[1]]
for inx,cell in enumerate(loop):
if inx%2==0:
cell_value=0 if np.isnan(solution[cell[0],cell[1]]) else solution[cell[0],cell[1]]
solution[cell[0],cell[1]]=cell_value+maxPlusValue
else:
solution[cell[0],cell[1]]-=maxPlusValue
solution[cell[0],cell[1]]=np.nan if solution[cell[0],cell[1]]==0 else solution[cell[0],cell[1]]
return solution
def isDegeneracy(self, solution):
m=solution.shape[0]
n=solution.shape[1]
num_nonnan=m*n-np.isnan(solution).sum()
if num_nonnan==m+n-1:
return False
else:
return True
ModifiedDistribution类实现了MODI算法，因为该算法是Stepping-Stone算法的变体，会复用到很多SteppingStone类的方法，所以直接从SteppingStone类继承；相比于SteppingStone类，主要添加了calculateUV方法，用于uuu和vvv的计算
import numpy as np
from SteppingStone import SteppingStone
class ModifiedDistribution(SteppingStone):
'''
Modified Distribution Method can be cosidered extending from Stepping Stone Method
'''
def calculateUV(self, solution):
'''
Calculate the values of extra u and v
'''
u=np.full((self.supply.shape[0],1),np.nan)# Extra Column to store u values
v=np.full((self.demand.shape[0],1),np.nan)# Extra Row to store v values
isnan_u=np.isnan(u)
isnan_v=np.isnan(v)
u[0]=0
m=solution.shape[0]
n=solution.shape[1]
while True in isnan_u or True in isnan_v:
for i in range(0,m):
for j in range(0,n):
if np.isnan(solution[i,j]):
continue
if np.isnan(u[i]) and np.isnan(v[j]):
continue
if np.isnan(u[i]) and np.isnan(v[j]) == False:
u[i]=self.cost[i,j]-v[j]
if np.isnan(u[i]) == False and np.isnan(v[j]):
v[j]=self.cost[i,j]-u[i]
isnan_u=np.isnan(u)
isnan_v=np.isnan(v)
return (u,v)
def calculateCostChange(self, u, v, solution):
'''
Calculate the cost change of each empty cell
'''
minChange=float('inf')
min_i=0
mim_j=0
m=solution.shape[0]
n=solution.shape[1]
for i in range(0,m):
for j in range(0,n):
if np.isnan(solution[i,j]) is False:
continue
k=self.cost[i,j]-u[i]-v[j]
if k<minChange:
minChange=k
min_i=i
min_j=j
return (minChange,min_i,min_j)
def calculate(self):
'''
1. Develop an initial solution using one of the three methods available.
2. Compute u_i and v_j values for each row and column by applying the formula u_i+v_j=c_ij
to each cell that has an allocation.
3. Compute the cost change, k_ij , for each empty cell using c_ij-u_i-v_j=k_ij
4. Allocate as much as possible to the empty cell that will result in the greatest net
decrease in cost (most negative k_ij). Allocate according to the stepping-stone path for
the selected cell.
5. Repeat steps 2 through 4 until all k_ij values are positive or zero.
'''
solution=self.initSolution
m=solution.shape[0]
n=solution.shape[1]
index=0
while True:
self.handleDegeneracy(solution)
(u,v)=self.calculateUV(solution)
(minChange,i,j)=self.calculateCostChange(u,v,solution)
if minChange>=0:
return solution
loop=self.findLoop(solution, [(i,j)], (i,j))
solution=self.reAllocate(solution,loop)
print('######################### iteration '+str(index)+" #############################")
print('Current solution : ')
print(solution)
index+=1
在main.py文件中，进行算法测试：
import numpy as np
from NorthWestCorner import NorthWestCorner
from MinCellCost import MinCellCost
from TPAlogirthm import TPAlgorithm
from SteppingStone import SteppingStone
from ModifiedDistribution import ModifiedDistribution
supply=np.array([150,175,275])
demand=np.array([200,100,300])
cost=np.array([
[6,8,10],
[7,11,11],
[4,5,12]
])
tpAlgorithm=TPAlgorithm(supply,demand,cost,'MinCellCost','ModifiedDistribution')
solution=tpAlgorithm.solve()
print('Final result : ')
print(solution)
print('Optimal cost : ')
print(tpAlgorithm.solutionCost(solution))
扫一扫
329
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
