Gurobi to solve Vehicle Routing Problem

[Trich]

Im trying to solve a classic Vehicle Routing Problem using Gurobi and Python and I was just wondering if anyone had any created code that I could reference when trying to solve it. I have code to a basic traveling salesman problem (just one vehicle), but i need to change this code a bit to allow for multiple vehicles and a demand at each place. If anyone has any ideas how to change the code listed below it would be greatly appreciated!


import numpy
from numpy import *
from gurobipy import *

# 1) Read data file (distances)

# TSP_data_xxx.csv
# CAUTION:  NEED TO STRIP ALL HEADER INFO FIRST!
c = numpy.genfromtxt('TSP_data_30_python.csv', dtype=None, delimiter=',')
[rows,cols] = c.shape

numCities = rows    # which also equals columns
Cities = range(1,numCities+1)

# print Cities

   
# 2) GUROBI
# Model
m = Model("TSP")


# a) Decision Variable Definitions and Objective Function:
decvarx = []
for r1 in range(0,numCities+1):
    decvarx.append([])
    decvarx[r1] = []
    for r2 in range(0,numCities+1):
        decvarx[r1].append([])
        decvarx[r1][r2] = []
       
print size(decvarx)
           
for k in Cities:
    # print k
    for l in Cities:
        # print l
        decvarx[k][l] = m.addVar(lb=0, ub=1, obj=float(c[k-1][l-1]), vtype=GRB.BINARY, name="x.%d.%d" % (k,l))
       
decvaru = []
decvaru.append([])
for k in Cities:
    # print k
    decvaru.append([])
    decvaru[k] = m.addVar(lb=0, ub=GRB.INFINITY, obj=0, vtype=GRB.CONTINUOUS, name="u.%d" % (k))

# The objective is to minimize the total travel distance.
m.modelSense = GRB.MINIMIZE

# Update model to integrate new variables
m.update()

# Constraint (1)
for j in Cities:
    # print j
    m.addConstr(quicksum(decvarx[i][j] for i in Cities) == 1, "Constr.1.%d" % j)

# Constraint (2)
for i in Cities:
    m.addConstr(quicksum(decvarx[i][j] for j in Cities) == 1, "Constr.2.%d" % i)

# Constraint (3)
for i in range(2,numCities+1):
    for j in range(2,numCities+1):
        if (i != j):
            exprl = LinExpr()
            exprl.clear()
            exprl.add(decvaru[i], 1.0)
            exprl.add(decvaru[j], -1.0)
            exprl.add(decvarx[i][j], numCities)
            exprr = LinExpr()
            exprr.clear()
            exprr.addConstant(numCities - 1)
            m.addConstr(exprl, GRB.LESS_EQUAL, exprr, "Constr.3.%d.%d" % (i, j))
                       
# Solve
m.optimize()

# Print solution
print '\nTOTAL DISTANCE:', m.objVal
print 'SOLUTION:'
for i in Cities:
    for j in Cities:
        if decvarx[i][j].x == 1.0:
            print 'x[',i,'][',j,'] = 1'

exit()