I am trying to code the three-index cvrp with CPLEX and python. I have found numerous codes that implement the two-index formulation 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.
Below is my 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] >= q[i] for i in N for k in K)
mdl.parameters.timelimit = 30
solution = mdl.solve(log_output=True)
Are there any other ways to code this formulation in CPLEX? I do need the three index formulation to add cost to vehicles.