CVRP subtour elimination: Cannot convert MLinExpr to gurobipy.LinExpr

Hello,

I am trying to implement the Capacitated Vehicle Routing Problem (CVRP) with gurobipy, based on this TSP example:

https://colab.research.google.com/github/Gurobi/modeling-examples/blob/master/traveling_salesman/tsp_gcl.ipynb 

I changed the variable type of x from tupledict (m.addVars) to mvar (m,addMVars), and I think this creates a problem when I implement the subtour elimination. 

The code:

model.cbLazy(gp.quicksum(x[i,j] for i, j in combinations(tour_elem, 2)) - (len(tour_elem)-1) <= 0)

raises an error, "TypeError: Cannot convert MLinExpr to gurobipy.LinExpr".

Traceback (most recent call last):
  File "src\gurobipy\callback.pxi", line 209, in gurobipy.CallbackClass.callback
  File "mycode.py", line 92, in subtourelim
    model.cbLazy(gp.quicksum(x[i,j] for i, j in combinations(tour_elem, 2)) - (len(tour_elem)-1) <= 0)
  File "src\gurobipy\model.pxi", line 6594, in gurobipy.Model.cbLazy
  File "src\gurobipy\model.pxi", line 6599, in gurobipy.Model.__cbCutOrLazy
TypeError: Cannot convert MLinExpr to gurobipy.LinExpr
Traceback (most recent call last):
  File "src\gurobipy\callback.pxi", line 209, in gurobipy.CallbackClass.callback

If this error occurred because m.cbLazy function only accepts LinExpr, is there any other way I can use MVars? Or should I use tupledict to use m.cbLazy function?

I attach my code and the CVRP ILP below.

import gurobipy as gp
import numpy as np
from itertools import combinations


def parse():
    """Returns a CVRP problem (a distance matrix, customer demands, and vehicle capacity)
    node 0 indicates the depot.

    :returns cmat : distance matrix (list of list, 11x11)
    :returns demand : customer demands (list of length 11)
    :returns tcap : vehicle capacity

    """
    lines=[[0.07080287595563761,0.8150640110845127,20],
    [0.8332237129635587,0.5759695086915411,5],
    [0.12091926778092488,0.21621155351811572,3],
    [0.7905205792134123,0.5169049580812282,4],
    [0.20318755190207172,0.9270361039593276,5],
    [0.814237813716852,0.5102894891274262,1],
    [0.2485747097938361,0.5440808890128546,7],
    [0.38636268653868444,0.10643108406133994,4],
    [0.07961575423644518,0.3554731208235127,3],
    [0.4373921928144826,0.8962394300910876,5],
    [0.9653080450372932,0.9827630650166035,2]]

    x,y,d=lines[0]
    coordx = [x]
    coordy = [y]
    demand = [0]
    tcap = d

    for l in lines[1:]:
        x,y,d=l
        coordx+=[x]
        coordy+=[y]
        demand+=[d]

    cmat=[]
    for i in range(11):
        xi,yi=coordx[i],coordy[i]
        cline=[]
        for j in range(11):
            xj,yj = coordx[j],coordy[j]
            dx,dy = xi-xj,yi-yj
            ce = (dx*dx+dy*dy) ** 0.5
            cline+=[ce]
        cmat+=[cline]

    return cmat, demand, tcap

"""VRP1: ILP ACVRP
V = {0,1,2,...,n}
minimize sum_{i,j in V} c_ij x_ij
subject to
sum_{i in V} x_ij=1  for j in V-{0}
sum_{j in V} x_ij=1  for i in V-{0}
sum_{i in V} x_i0=k
sum_{j in V} x_0j=k
x_ij binary
(subtour elimination)
sum_{i not in s, j in s} x_ij >= r(S) for all S subset of V-{0}, S!={}
or sum_{i, j in S} x_ij <= |S|-r(S) for all S subset of V-{0}, S!={}
"""

cmat, demand, tcap = parse()

# Create a new model
m = gp.Model()

x = m.addMVar((11,11), vtype='B', name="x")
c = np.array(cmat)  # 2d distance matrix
d = np.array(demand)  # 1d demand vector

# Set objective function
m.setObjective((c*x).sum(), gp.GRB.MINIMIZE)
m.addConstrs(x[i,:].sum() == 1 for i in range(1,11))
m.addConstrs(x[:,i].sum() == 1 for i in range(1,11))
m.addConstr(x[0,:].sum() == x[:,0].sum())

# Callback - use lazy constraints to eliminate sub-tours
def subtourelim(model, where):
    if where == gp.GRB.Callback.MIPSOL:
        # make a list of edges selected in the solution
        vals = model.cbGetSolution(x)
        selected = [[j for j in range(11) if vals[i, j] > 0.5] for i in range(11)]
        # find the shortest cycle in the selected edge list
        tour_elem = subtour(selected)
        if len(tour_elem)>0:
            # add subtour elimination constr. for every pair of cities in subtour
            # doesn't work because of the incompatibility of mlinexpr and linexpr??
            model.cbLazy(gp.quicksum(x[i,j] for i, j in combinations(tour_elem, 2)) - (len(tour_elem)-1) <= 0)

def subtour(edges):
    """Check unvisited nodes from the depot and return them."""
    unvisited = [0]+[1]*10
    q=[0]
    while q:
        current=q.pop()
        for j in edges[current]:
            if unvisited[j]:
                unvisited[j]=0
                q+=[j]
    cycle_elems=[]
    for i,e in enumerate(unvisited):
        if e==1:
            cycle_elems+=[i]
    # todo: optimization (find minimum subtour) and capacity constraints
    return cycle_elems


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

# Solve it!
m.optimize()

print(f"Optimal objective value: {m.objVal}")

Thank you!