Uses Large Neighborhood Search for optimizing Vehicle Routing Problem with Time Windows. This algorithm can easily be improved.

VRPTWLNS.py

import csv
import pulp            #This contains LP Modeling and solver
import time            #This is used to see how much time for solving
import random

#######################################################################
#######################################################################



NumberOfPoints = 10  # number of points/locations in the data

k = 3  # number of vehicles

BigM = 50

StartTime = 12          #12 means 6am. 14 means 8am. half hour is 1 time unit.

EndTime = 44            #44 means 22:00 or 10pm.

ObjectiveWeightForVehicles = [10,15,20] #length of this list should be equal to k (number of vehicles)

IterationLimit = 100

#######################################################################
#######################################################################

assert len(ObjectiveWeightForVehicles) == k  #checking if that condition satisfies.

print()
print("Solving the problem of %d vehicles and %d locations... Please wait" % (k,NumberOfPoints))
print()



#------------------------------------------------------------------------------------------------------
# GETTING & PREPARING DATA
#------------------------------------------------------------------------------------------------------


vehicleList = range(k)
Points = range(1,NumberOfPoints+1)
Nodes = range(0,NumberOfPoints+2) #there are 52 nodes as 0 and numberofpoints+1 are entry and exit depots respectively.



with open('travel_times.csv') as csvfile:
    data = list(csv.reader(csvfile))


travelTimes = {}
for i in Points:
    for j in Points:
        if data[j][i].isdigit() and i!=j:
            travelTimes[(i,j)] = float(int(data[j][i])/10)


#Addition of 0 and NumberOfPoints+1 points as they are depots in this model.
travelTimes[(0,NumberOfPoints+1)] = 0
for i in Points:
    travelTimes[(0,i)] = 0
    travelTimes[(i,(NumberOfPoints+1))] = 0


#Addition of symmetrical travel times
for i in range(2,NumberOfPoints+1):
    for j in range(1,i):
        travelTimes[(i,j)] = travelTimes[(j,i)]

travelFrom,travelTo = zip(*list(travelTimes.keys()))


#-----------------------------------------------------------------------------

with open('location_data.csv') as csvfile1:
    data1 = list(csv.reader(csvfile1))


LocationData = dict([(point,[]) for point in Points])

for j in [1,6,7,9]:
    for i in Points:
        LocationData[i].append(data1[i][j])

CapacityOfVehicle = float(data1[1][4])






#------------------------------------------------------------------------------------------------------
# LP MODELING
#------------------------------------------------------------------------------------------------------


VRP = pulp.LpProblem("Assign_FromPointId_ToPointId_VehicleId, StartTime(S)_PointId_VehicleId", pulp.LpMinimize)

DecisionVar = pulp.LpVariable.dicts("Assign",(Nodes,Nodes,vehicleList),0,1,pulp.LpBinary)  #Binary

TimeDecisionVar = pulp.LpVariable.dicts("Start",(Nodes,vehicleList),0,None, pulp.LpContinuous)
for point in Points:
    for vehicle in vehicleList:
        TimeDecisionVar[point][vehicle].bounds(float(LocationData[point][1]), float(LocationData[point][2]))



VRP += 1 #for getting feasible solution first


for point in Points:
    VRP += pulp.lpSum([DecisionVar[point][node][vehicle] for node in Nodes
                       for vehicle in vehicleList if point!=node]) == 1, "eachCustomerOnce"+str(point)
    #constraint1 (2.2 in mip model.png)



for vehicle in vehicleList:
    VRP += pulp.lpSum([float(LocationData[point][0]) * DecisionVar[point][node][vehicle] for point in Points
                       for node in Nodes if point!=node]) <= float(CapacityOfVehicle), "Capacity"+str(vehicle)
    #constraint2 (2.3 in mip model.png)



for vehicle in vehicleList:
    VRP += pulp.lpSum([DecisionVar[0][node][vehicle] for node in
                       Nodes if node!=0]) == 1, "entryDepotConnection"+str(vehicle)
    #constraint3 (2.4 in mip model.png)
    #there are 21 nodes as 0 and NumberOfPoints+1 are entry and exit depots respectively.



for vehicle in vehicleList:
    VRP += pulp.lpSum([DecisionVar[node][NumberOfPoints+1][vehicle] for node in Nodes
                       if node!=NumberOfPoints+1]) == 1, "exitDepotConnection"+str(vehicle)
    #constraint4 (2.6 in mip model.png)



for point in Points:
    for vehicle in vehicleList:
        VRP += pulp.lpSum([DecisionVar[node][point][vehicle] - DecisionVar[point][node][vehicle] for node in
                           Nodes if node!=point]) == 0, "forTrip"+str(point)+'k'+str(vehicle)
    #constraint5 (2.5 in mip model.png)



for point in Points:
    for vehicle in vehicleList:
        VRP += DecisionVar[point][0][vehicle] == 0, "ToEntryDepot"+str(point)+str(vehicle)
        VRP += DecisionVar[NumberOfPoints+1][point][vehicle] == 0, "FromExitDepot"+str(point)+str(vehicle)
    #constraint6 - For constraint5 to work correctly.



for vehicle in vehicleList:
    for node1 in Nodes:
        if node1 == NumberOfPoints+1:
            continue
        else:
            for node2 in Nodes:
                if node2 == 0:
                    continue
                else:
                    if (node1 != node2):
                        if (node1 == 0):
                            VRP += TimeDecisionVar[node1][vehicle] - TimeDecisionVar[node2][vehicle]\
                   + BigM * DecisionVar[node1][node2][vehicle] <= float(BigM - travelTimes[node1,node2]), \
                               "timewindow"+str(vehicle)+'p'+str(node1)+'p'+str(node2)
                        else:
                            VRP += TimeDecisionVar[node1][vehicle] - TimeDecisionVar[node2][vehicle] \
                                   + BigM * DecisionVar[node1][node2][vehicle] <= float(
                                BigM - travelTimes[node1, node2] - float(LocationData[node1][3])), \
                                   "timewindow" + str(vehicle) + 'p' + str(node1) + 'p' + str(node2)



for vehicle in vehicleList:
    VRP += TimeDecisionVar[NumberOfPoints+1][vehicle] <= EndTime, "EndTime"+str(vehicle)
    VRP += TimeDecisionVar[0][vehicle] >= StartTime, "StartTime"+str(vehicle)


VRP.writeLP("VRPLNSFEASIBLE.lp")


tic = time.time()


VRP.solve() #print("Status:", pulp.LpStatus[VRP.status])


Obj = pulp.LpConstraintVar(name="obj", e=sum([ObjectiveWeightForVehicles[vehicle] * DecisionVar[From][To][vehicle] for From,To in zip(travelFrom,travelTo)
            for vehicle in vehicleList]))

VRP.setObjective(Obj)

VRP.writeLP("VRPLNSOPTIMAL.lp")


bestObjectiveValue = float(pulp.value(VRP.objective))
print("****Feasible solution found, current best objective value: ",float(pulp.value(VRP.objective)))
print()

decVarName = []
bestSolution = []

count=0
for v in VRP.variables():
    if v.name[0] == "A":  #print(v.name,"and",v.varValue)
        decVarName.append(v.name)
        bestSolution.append(v.varValue)
        count +=1 #print(count)


i = 0
while True:
    if i == IterationLimit:
        break
    remainSame = []
    remainSame = random.sample(range(count),int((count*2)/3))   #print(len(remainSame))

    VRP1 = VRP.deepcopy()


    for vehicle in vehicleList:
       for node1 in Nodes:
           for node2 in Nodes:
               for index in remainSame:
                   if DecisionVar[node1][node2][vehicle].name == decVarName[index]:
                       VRP1 += DecisionVar[node1][node2][vehicle] == bestSolution[index], "LNSSetToCurrentBest"+str(index)


    #VRP1.writeLP(str(i)+"VRPONELNS.lp")
    VRP1.solve()

    if (pulp.LpStatus[VRP1.status] == "Optimal") and (float((pulp.value(VRP1.objective))) < bestObjectiveValue):
        print("*** New best solution found")
        bestObjectiveValue = float(pulp.value(VRP1.objective))
        bestSolution = []
        test=[]
        for v in VRP1.variables():
            if v.name[0] == "A":
                bestSolution.append(v.varValue)
        print("Updating best objective value to ", bestObjectiveValue)
        print()

    del VRP1
    print("LNS Iteration: ",i)
    i +=1


print()
print("****LNS is done.")
print()
print("Best solution is: ")
FinalSolution = zip(decVarName,bestSolution)
for x,y in FinalSolution:
    if y != 0:
        print(x,"=",y)




toc = time.time()



#for vehicle in vehicleList:
#    for node1 in Nodes:
#        for node2 in Nodes:
#            if DecisionVar[node1][node2][vehicle].value() == 1.0:
#                print("Point:"+ str(node1)+ " to Point:" + str(node2) + " by vehicle:" + str(vehicle))

#print()
#print("Total cost/Objective value:",pulp.value(VRP.objective))
#print()
#print("**See VRPOutput.csv and VRPMinVehicle.lp for more details.**")
print()
print("Time taken for solving is " + str((toc-tic)) + " sec")