Using Column Generation for Optimizing Vehicle Routing Problem with Time Windows. The Objective here is only to minimize the number of vehicles used. Uses PuLP library in python.

0VRPTWColGen.py

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


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



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

k = 1  # 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 = [1] #length of this list should be equal to k (number of vehicles)



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

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


#------------------------------------------------------------------------------------------------------
# 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]:
    for i in Points:
        LocationData[i].append(data1[i][j])

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







#------------------------------------------------------------------------------------------------------
# MASTER PROBLEM MIP MODELING
#------------------------------------------------------------------------------------------------------


Paths = list(range(NumberOfPoints+1))

MASVRP = pulp.LpProblem("Master Problem VRP", pulp.LpMinimize)

MDecisionVar = []

for path in Paths:
    MDecisionVar.append(pulp.LpVariable("Path"+str(path),0,1,pulp.LpContinuous))


Obj = pulp.LpConstraintVar(name="OBJ",e=sum([MDecisionVar[path] for path in Paths]))
MASVRP.setObjective(Obj)


Customer = []
for path in Paths:
    Customer.append(pulp.LpConstraintVar(name="Customer"+str(path),sense=pulp.LpConstraintEQ,rhs=1,e=MDecisionVar[path]))
    MASVRP += Customer[path]





#------------------------------------------------------------------------------------------------------
# SUBPROBLEM MIP MODELING
#------------------------------------------------------------------------------------------------------


SUBVRP = pulp.LpProblem("Sub Problem VRP", pulp.LpMinimize)


ConstantCost = pulp.LpVariable("ConstantCost",1,1)
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]))



for vehicle in vehicleList:
    SUBVRP += 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:
    SUBVRP += 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:
    SUBVRP += 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:
        SUBVRP += 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:
        SUBVRP += DecisionVar[point][0][vehicle] == 0, "ToEntryDepot"+str(point)+str(vehicle)
        SUBVRP += 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):
                        SUBVRP += 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)



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





tic = time.time()




#------------------------------------------------------------------------------------------------------
# COLUMN GENERATION ITERATIONS
#------------------------------------------------------------------------------------------------------



i = 0
while True:
    MASVRP.writeLP(str(i)+"Master.lp")
    MASVRP.solve()
    price = {}
    price.clear()

    for path in range(0,NumberOfPoints+1):
        price[path] = float(MASVRP.constraints["Customer"+str(path)].pi)

    print("Dual values: ",price)
    print()

    SUBVRP += sum([ConstantCost-[price[From] * DecisionVar[From][To][vehicle] for From in Nodes for To in Nodes
         for vehicle in vehicleList if not(From == To or From == NumberOfPoints+1)]])
    # Objective

    SUBVRP.solve()

    SUBVRP.writeLP(str(i)+"SUBVRP.lp")

    if (pulp.value(SUBVRP.objective) > -1) or i==50:
        break

    expression = Obj
    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))
                    if not (node1 == NumberOfPoints+1):
                        expression += Customer[node1]
    print()
    print("Master LP problem objective value: ",pulp.value(MASVRP.objective))
    print("Subproblem Objective value: ",pulp.value(SUBVRP.objective))
    print()

   # print(expression)
   # print()

    MDecisionVar.append(pulp.LpVariable("Path"+str(len(Paths)),0,1,pulp.LpContinuous,e=expression))

    Paths.append(len(Paths))

    print("*************************************************************")
    print("Column Generation Iteration: ",i)
    print()
    i+=1




print("##################*********************######################")
print("****** Column Generation Done. Setting all variables to integers in the Master Problem and solving it... ******")
print()
for variables in MASVRP.variables():
    variables.cat = pulp.LpInteger


MASVRP.writeLP("MasterInteger.lp")

MASVRP.solve()

toc = time.time()


for path in Paths:
    if MDecisionVar[path].value() == 1.0:
        print("Path selected is: ",str(path))

print()
print("Objective value:",pulp.value(MASVRP.objective))


print("Time taken for solving is " + str((toc-tic)) + " sec")
print()