For understanding basics
of column generation please refer to my previous post http://practicalappliedanalytics.blogspot.in/2017/09/column-generation-column.html
The informal formulation
for VRPTW (Desaulniers, G.,
Desrosiers, J. and Solomon, M.M., Column generation, Chapter 3,Springer, 2005,
DOI: 10.1007/b135457) is as follows:
The restricted master problem for path generation is:
The sub-problem with
reduced cost factor in the objective is:
Objective of sub-problem
is
Assumptions:
1.
Homogeneous fleet of
vehicles
2.
Starting and ending
depot are same
3.
Cost is assumed to be
proportional to distance
4.
Travel time is
calculated with assuming fixed vehicle speed in all lanes
5.
The solution of the
model might require further branching on vehicles or flow variables to get
integer solution. Here column generated model is build
6. Upload 'VRPTWData.csv'
file for extracting input parameters. This data is based on Solomon instances.
(Link:https://gist.github.com/adarsh70/f4ecb197ce0bc2252bd8c01f6d5978f1)
#Vehicle Routing Problem with time windows (Adapted from Solomon)
from pulp import*
import pandas as pd
#Data for the model both in list and dataframe
aggrdat = pd.DataFrame.from_csv('VRPTWData.csv',index_col = ['CustNo'])
aggrdatrec = pd.DataFrame.to_records(aggrdat)
aggrdatser = pd.Series(aggrdatrec)
aggrdatlis = aggrdatser.tolist()
nodesdf = aggrdat.iloc[:,[]]
custdf = aggrdat.iloc[1:,[]]
demanddf = aggrdat.iloc[1:,[2]]
timewindowdf = aggrdat.iloc[0:,[3,4,5]]
nodes = []
for i,x in enumerate(aggrdatlis):
for j,y in enumerate(x):
if j == 0:
nodes.append(y)
demand = []
for i,x in enumerate(aggrdatlis):
if i != 0:
demand.extend([([x[0],x[3]])])
timewindow = []
for i,x in enumerate(aggrdatlis):
timewindow.extend([[x[0],x[4],x[5],x[6]]])
cust = []
for i,x in enumerate(nodes):
if i > 0:
cust.append(x)
depot = []
for i,x in enumerate(nodes):
if i == 0:
depot.append(i)
loccoord = []
for i,x in enumerate(aggrdatlis):
loccoord.extend([[x[0],x[1],x[2]]])
#Calculate distance and time based on assumed vehicle speed of 5 units per time period
def dist(a,b):
dist = []
dm = (((a[1] - b[1])**2) + ((a[2] - b[2])**2))**(1/2)
dist.extend([a[0],b[0],dm])
return dist
def time(a,b):
time = []
dm = ((((a[1] - b[1])**2) + ((a[2] - b[2])**2))**(1/2))/5
time.extend([a[0],b[0],dm])
return time
arcs = []
for i,x in enumerate(loccoord):
for j,y in enumerate(loccoord):
arcs.extend([[x[0],y[0]]])
arccost = []
for i,x in enumerate(loccoord):
for j,y in enumerate(loccoord):
arccost.append(dist(x,y))
arctime = []
for i,x in enumerate(loccoord):
for j,y in enumerate(loccoord):
arctime.append(time(x,y))
M = []
for i,x in enumerate(arctime):
s = timewindow[x[0]-1][2] - timewindow[x[1]-1][1] + x[2]
M.extend([[x[0],x[1],s]])
arccostdf = pd.DataFrame(arccost,columns=['fromnode','tonode','cost'])
arccostdf = arccostdf.set_index(['fromnode','tonode'])
arcsdf = pd.DataFrame(arcs,columns=['fromnode','tonode'])
arctimedf = pd.DataFrame(arctime,columns=['fromnode','tonode','time'])
arctimedf = arctimedf.set_index(['fromnode','tonode'])
Mdf = pd.DataFrame(M,columns=['fromnode','tonode','Mvalue'])
Mdf = Mdf.set_index(['fromnode','tonode'])
#Master Problem relaxed for calculating duals
class MasterProblem:
def __init__(self,nodes,arcs,cust,depot,demand,timewindow,arccost,arctime,vehcap,genpaths,arccostdf,demanddf,arcsdf,nodesdf,custdf,arctimedf,timewindowdf,Mdf,problemname):
self.nodes = nodes
self.arcs = arcs
self.cust = cust
self.depot = depot
self.demand = demand
self.timewindow = timewindow
self.arccost = arccost
self.arctime = arctime
self.vehcap = vehcap
self.genpaths = genpaths
self.arccostdf = arccostdf
self.demanddf = demanddf
self.arcsdf = arcsdf
self.nodesdf = nodesdf
self.custdf = custdf
self.arctimedf = arctimedf
self.timewindowdf = timewindowdf
self.Mdf = Mdf
self.problemname = problemname
self.prob = LpProblem(problemname,LpMinimize)
self.obj = LpConstraintVar("obj")
self.prob.setObjective(self.obj)
self.PathVars = []
self.constraintList = []
for i,x in enumerate(cust):
var = LpConstraintVar("Cust"+str(i),LpConstraintEQ,1)
self.constraintList.append(var)
self.prob += var
ct = []
for i,x in enumerate(self.genpaths):
temp = []
for j,y in enumerate(x):
if y>0:
for a,b in enumerate(self.arcs):
if a == j:
temp.append(b[0])
temp.append(b[1])
ct.extend([temp])
custcount = []
temp = []
for i,x in enumerate(ct):
for a,b in enumerate(self.cust):
temp.append(x.count(b))
custcount.extend([temp])
temp = []
for i,x in enumerate(self.genpaths):
pathcost = 0
for j,y in enumerate(x):
if y>0:
pathcost += (self.arccost[j][2])
var = LpVariable("Path"+str(i),0,None,LpContinuous,lpSum(self.obj * pathcost + [self.constraintList[a] * b for a,b in enumerate(custcount[i])]))
self.PathVars.append(var)
def solve(self):
self.prob.writeLP("Masterprob.lp")
self.prob.solve()
return [self.prob.constraints[i].pi for i in self.prob.constraints]
def addPath(self,path):
self.genpaths.append(path)
temp = []
for j,y in enumerate(path):
if y>0:
for a,b in enumerate(self.arcs):
if a == j:
temp.append(b[0])
temp.append(b[1])
custcount = []
for i,x in enumerate(self.cust):
custcount.append(temp.count(x))
pathcost = 0
for j,y in enumerate(path):
if y>0:
pathcost += (self.arccost[j][2])
var = LpVariable("Path"+str(len(self.genpaths)),0,None,LpContinuous,lpSum(self.obj * pathcost + [self.constraintList[a] * b for a,b in enumerate(custcount)]))
self.PathVars.append(var)
def startSubProblem(self,duals):
newSubProb = SubProblem(duals,self.arcs,self.nodes,self.vehcap,self.nodesdf,self.custdf,self.arccostdf,self.demanddf,self.arcsdf,self.timewindowdf,self.Mdf,self.arctimedf)
path=newSubProb.returnPath()
return path
def setRelaxed(self,relaxed):
if relaxed==False:
for var in self.prob.variables():
var.cat = LpBinary
def getObjective(self):
return value(self.prob.objective)
def getUsedPaths(self):
usedPathList=[]
for i,x in enumerate(self.PathVars):
if value(x)>0:
usedPathList.append((value(x),self.genpaths[i]))
return usedPathList
#Sub Problem for finding new path variables based on reduced cost
class SubProblem:
def __init__(self,duals,arcs,nodes,vehcap,nodesdf,custdf,arccostdf,demanddf,arcsdf,timewindowdf,Mdf,arctimedf):
self.subprob = LpProblem("Sub Problem solver",LpMinimize)
self.varList = LpVariable.dicts("Arc_variable", ((i,j) for i,j in arcs),cat='Binary')
self.subprob += lpSum([self.varList[i,j] * (arccostdf.loc[(i,j),'cost'] - duals[i - 2]) for i,j in arccostdf.index if i > 1]
+ [self.varList[i,j] * arccostdf.loc[(i,j),'cost'] for i,j in arccostdf.index if i == 1])
self.subprob += lpSum([self.varList[i,j] * demanddf.loc[(i),'Demand'] for i,j in arcs if i != 1]) <= vehcap
self.subprob += lpSum([self.varList[i,j] for i,j in arcs if i == 1]) == 1
self.subprob += lpSum([self.varList[i,j] for i,j in arcs if j == 1]) == 1
for i in nodes:
for j in nodes:
if i == j:
self.subprob += self.varList[i,j] == 0
for c in cust:
self.subprob += (lpSum([self.varList[i,j] for i,j in arcs if j == c]) == lpSum([self.varList[i,j] for i,j in arcs if i == c]))
self.twvarList = LpVariable.dicts("Time_variable", ((i) for i in nodes),cat='Continuous')
for i in nodes:
self.subprob += self.twvarList[i] >= timewindowdf.loc[(i),'ReadyTime']
self.subprob += self.twvarList[i] <= timewindowdf.loc[(i),'DueDate']
for i in nodes:
for j in nodes:
self.subprob += (lpSum([self.twvarList[i] - self.twvarList[j] + Mdf.loc[(i,j),'Mvalue'] * self.varList[i,j]]) <= lpSum([Mdf.loc[(i,j),'Mvalue'] - arctimedf.loc[(i,j),'time']]))
self.subprob.writeLP('subprob.lp')
self.subprob.solve()
def returnPath(self):
path=False
if value(self.subprob.objective) < 0:
path=[]
for i,v in enumerate(self.subprob.variables()):
if i < len(arcs):
path.append(v.varValue)
return path
#Generate initial paths for each customer individually (i.e. Depot-Cust-Depot) and run iteration on Master and Sub-Problem to get the solution
def conjugatepos(arcpos,arcs):
for i,x in enumerate(arcs):
if i == arcpos:
first = x[1]
second = x[0]
for i,x in enumerate(arcs):
if x[0] == first and x[1] == second:
pos = i
return pos
genpaths=[]
print ("Generating start paths:")
## generate simple start paths
for x in range(len(cust)):
temp=[0.0 for y in range(len(arcs))]
genpaths.append(temp)
for i,x in enumerate(genpaths):
genpaths[i][i+1] = 1.0
conjpos = conjugatepos(i+1,arcs)
genpaths[i][conjpos] = 1.0
print (genpaths[i])
print ("\n\nTrying to solve problem")
CGprob=MasterProblem(nodes,arcs,cust,depot,demand,timewindow,arccost,arctime,20,genpaths,arccostdf,demanddf,arcsdf,nodesdf,custdf,arctimedf,timewindowdf,Mdf,'Vehicle Routing Problem with time windows')
relaxed=True
while relaxed==True:
duals=CGprob.solve()
newPath=CGprob.startSubProblem(duals)
print ('New path: %s' % newPath)
if newPath:
CGprob.addPath(newPath)
else:
CGprob.setRelaxed(False)
CGprob.solve()
relaxed=False
print ("\n\nSolution: %s Optimal path cost is" % CGprob.getObjective())
t=CGprob.getUsedPaths()
for i,x in enumerate(t):
print ("Path %s: selected %s times %s" % (i,x[0],x[1]))
No comments:
Post a Comment