title: three-index vehicle flow formulation CVRP with docplex and python error
body: <p>I am trying to code the three-index cvrp with CPLEX and python. I have found numerous codes that implement the <a href="https://github.com/industrial-ucn/jupyter-examples/blob/master/optimization/cvrp-cplex.ipynb" rel="nofollow noreferrer">two-index formulation</a> but have not been able to find the extended version. Following the github code, I am using the generalized miller-tucker-zemlin subtour constraint (see picture 1.37-1.38), but I keep getting unfeasible solutions, regardless of the number of vehicles I choose.</p>
  <p><a href="https://i.stack.imgur.com/ewjuD.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/ewjuD.png" alt="enter image description here" /></a></p>
  <p>Below is my code:</p>
  <pre><code>import numpy as np
  import matplotlib.pyplot as plt
  from docplex.mp.model import Model
  
  rnd = np.random
  rnd.seed(0)
  
  n = 10
  Q = 20
  N = [i for i in range(1, n+1)]
  V = [0] + N
  q = {i: rnd.randint(1, 10) for i in N}
  K = range(1,4)
  loc_x = rnd.rand(len(V))*200
  loc_y = rnd.rand(len(V))*100
  
  A = [(i, j, k) for i in V for j in V if i != j for k in K]
  Y = [(i, k) for i in V for k in K]
  U = [(i, k) for i in N for k in K]
  # c = {(i, j): np.hypot(loc_x[i]-loc_x[j], loc_y[i]-loc_y[j]) for i, j, k in A}
  c = {(i, j): np.hypot(loc_x[i]-loc_x[j], loc_y[i]-loc_y[j]) for i, j, k in A}
  
  mdl = Model('CVRP')
  
  x = mdl.binary_var_dict(A, name='x')
  y = mdl.binary_var_dict(Y, name='y')
  u = mdl.continuous_var_dict(U, ub=Q, name='u')
  
  mdl.minimize(mdl.sum(c[i, j]*x[i, j, k] for i, j, k in A))
  
  # constraint 1
  #for i in N:
  mdl.add_constraints(mdl.sum(y[i, k] for k in K) == 1 for i in N)
  
  # constraint 2
  mdl.add_constraint(mdl.sum(y[0, k] for k in K) == len(K))
  
  # constraint 3
  mdl.add_constraints(mdl.sum(x[i, j, k] for j in V if i != j) - mdl.sum(x[j, i, k] for j in V if i != j) - y[i, k] == 0 for i in V for k in K)
  
  # constraints 4 and 5 subtour elimination and capacity
  mdl.add_indicator_constraints(mdl.indicator_constraint(x[i, j, k], u[i, k]+q[j] == u[j, k]) for i, j, k in A if i != 0 and j != 0)
  mdl.add_constraints(u[i, k] &gt;= q[i] for i in N for k in K)
  
  mdl.parameters.timelimit = 30
  solution = mdl.solve(log_output=True)
  </code></pre>
  <p>Are there any other ways to code this formulation in CPLEX? I do need the three index formulation to add cost to vehicles.</p>