Correctly formulate a constraint(Distance constraint on TSP)

Hello,

I am using a python IDE with gurobi to model and solve a problem for an assignment. The problem, although not stated seems to be an OVRP(with some extra constraints). The code i have right now is:

from itertools import combinations
import gurobipy as gp
from gurobipy import GRB

def find_length(list):
    length = 0
    for i in range(len(list)):
        length += distances[list[i][0]][list[i][1]]
    return length


def find_tours(list):
    tours = []
    for start in list:
        this_tour = []
        begin = start[0]
        end = start[1]
        if begin == 0:

            this_tour.append(start)
            current = start
            while next_route_exists(current, list):
                current = find_next_route(current, list)
                this_tour.append(current)

            tours.append(this_tour)
    return tours


def find_next_route(current_route, list):
    end = current_route[1]
    for tuple in list:
        if tuple[0] == end:
            list.remove(tuple)
            return tuple


def next_route_exists(current_route, list):
    end = current_route[1]
    for tuple in list:
        if tuple[0] == end:
            return True
    return False


def subtourelim(model, where):
    if where == GRB.Callback.MIPSOL:
        # make a list of edges selected in the solution
        vals = model.cbGetSolution(model._vars)
        selected = gp.tuplelist((i, j) for i, j in model._vars.keys() if vals[i, j] > 0.5)

        tours = find_tours(selected)
        for tour in tours:
            if find_length(tour) > 400:
                # add subtour elimination constr. for every pair of cities in tour
                model.cbLazy(gp.quicksum(distances[i][j] for i, j in tour) <= 400)


# Distances between cities

distances = [[0, 180, 240, 85, 285, 205, 235, 255, 155, 120, 230, 340, 220, 160, 240],  # Α
             [0, 0, 255, 150, 125, 100, 175, 235, 25, 65, 95, 210, 55, 135, 215],  # Β1
             [0, 0, 0, 160, 235, 160, 105, 45, 245, 255, 195, 225, 310, 115, 50],  # Β2
             [0, 0, 0, 0, 215, 125, 150, 170, 130, 110, 160, 260, 200, 75, 150],  # Β3
             [0, 0, 0, 0, 0, 90, 130, 195, 150, 190, 50, 85, 155, 155, 185],  # Β4
             [0, 0, 0, 0, 0, 0, 70, 135, 110, 135, 40, 135, 155, 65, 120],  # Β5
             [0, 0, 0, 0, 0, 0, 0, 70, 175, 200, 95, 125, 225, 70, 55],  # Β6
             [0, 0, 0, 0, 0, 0, 0, 0, 235, 250, 165, 180, 290, 110, 20],  # Β7
             [0, 0, 0, 0, 0, 0, 0, 0, 0, 40, 110, 225, 70, 130, 215],  # Β8
             [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 145, 265, 100, 135, 230],  # Β9
             [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 115, 135, 105, 150],  # Β10
             [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 240, 185, 175],  # Β11
             [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 190, 270],  # Β12
             [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 90],  # Β13
             [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]  # B14

n = 15
for i in range(n):
    for j in range(n):
        if j > i:
            distances[j][i] = distances[i][j]

dist = {(i, j): distances[i][j] for i in range(n) for j in range(n)}
m = gp.Model()

vars = m.addVars(dist.keys(), obj=dist, vtype=GRB.BINARY, name='route')

m.addConstrs(vars.sum(i, '*') <= 1 for i in range(1, n))  # a city might be the departure
m.addConstrs(vars.sum('*', i) == 1 for i in range(1, n))  # every city must be the destination

m.addConstr(vars.sum(0, '*') >= 1, "start")  # could have more than 1 vehicles
m.addConstr(vars.sum('*', 0) == 0)  # can't go back to start
m.addConstr(vars.sum(0, '*') <= 8, "start")  # can't have more than 8 vehicle

m.addConstrs(vars[i, i] == 0 for i in range(1, n))  # can't go to  current city from current city
m.addConstrs(vars[i, j] + vars[j, i] <= 1 for i in range(1, n) for j in range(1, n) if  j != i)  # eliminate double edges between cities

route_km = vars.prod(dist)
m.addConstr(route_km <= 400 * vars.sum(0, '*'))  # total km driven




# total vehicles used

total_vehicles = vars.sum(0, '*')
m.setObjective(total_vehicles, GRB.MINIMIZE)

# optimize model

m._vars = vars
m.Params.lazyConstraints = 1
m.optimize(subtourelim)
#m.optimize()

# get edges

vals = m.getAttr('x', vars)
selected = gp.tuplelist((i, j) for i, j in vals.keys() if vals[i, j] > 0.5)
print('')
tours = find_tours(selected)

i = 0
total=0
for tour in tours:
    i += 1
    print('Tour No ' + str(i) + ':')
    print(tour)
    print('is ' + str(find_length(tour)) + ' km long')
    print('')
    total+=find_length(tour)
print('Total : ' + str(total) + ' km')

The problem i am having right now is that, since i am new to gurobi, i can't express a constraint i want to add properly and that leads to a wrong solution coming up. The problem is with this constraint:

route_km = vars.prod(dist)
m.addConstr(route_km <= 400 * vars.sum(0, '*'))  # total km driven

which is wrong. I wanted to add a constraint that each truck can drive up to 400 km. Once that distance is driven or can't be driven if another "city" is chosen, a new truck has to be deployed(total number of trucks is 8, as stated in the constraints). If this constraint is formulated i believe i can eliminate this lazy constraint as well:

tours = find_tours(selected)
for tour in tours:
    if find_length(tour) > 400:
        # add subtour elimination constr. for every pair of cities in tour
        model.cbLazy(gp.quicksum(distances[i][j] for i, j in tour) <= 400)

Thank you in advance for any help,

George