# Installation solvers
!pip install cplex
!pip install docplex
from docplex.mp.model import Model
import matplotlib.pyplot as plt
import numpy as np

# Define a problem:
Nodes    = ["depot", "node1", "node2", "node3", "node4", "node5"];
Nodes_1  = ["node1", "node2", "node3", "node4", "node5"];
Vehicles = ["v1", "v2"];
depots   = ["depot"]

I   = range(len(Nodes))
J   = range(len(Nodes))
I_1 = range(len(Nodes_1))
J_1 = range(len(Nodes_1))
P   = range(len(Nodes))
K   = range(len(Vehicles))
I_2 = range(len(depots))
J_2 = range(len(depots))


capacity = [5, 5]
demand   = [0, 1, 1, 1, 1, 1]
cost = [
[1000,13,11,12,11,12],
[13,1000,10,13,14,12],
[11,10,1000,12,13,15],
[12,13,12,1000,15,15],
[11,14,13,15,1000,13],
[12,12,15,15,13,1000]
]

# Create a CPLEX model:
mdl = Model('CVRP')

# Define arcs and capacities:
x = {}
for i in I:
  for j in J:
    for k in K:
      x[(i,j,k)] = mdl.binary_var(name='x_{0}_{1}_{2}'.format(i,j,k))

u = {}
for i in I_1:
  u[(i)] = mdl.continuous_var(lb= demand[i], ub=5, name='u_{0}'.format(i))

# Define objective function:
mdl.minimize(mdl.sum(x[(i,j,k)]*cost[i1][j1] for i1, i in enumerate(I) for j1, j in enumerate(J)));

# Add constraints:
for j in J:
  if(j>0):
    mdl.add_constraint(mdl.sum(x[(i,j,k)] for i in I for k in K if i!=j) == 1, "c1")

for i in I:
  if(i>0):
    mdl.add_constraint(mdl.sum(x[(i,j,k)] for j in J for k in K if i!=j) == 1, "c2")

for p in P:
  for k in K:
    mdl.add_constraint(mdl.sum(x[(i,p,k)] for i in I) -  mdl.sum(x[(p,j,k)] for j in J) == 0, "c3")

for k, k1 in enumerate(K):
  mdl.add_constraint(mdl.sum(demand[i1]*x[(i,j,k)] for i,i1 in enumerate(I) for j in J) <= capacity[k1], "c4")

for k in K:
   mdl.add_constraint(mdl.sum(x[(0,j,k)] for j in J) == 1, "c5")

for k in K:
   mdl.add_constraint(mdl.sum(x[(i,0,k)] for i in I) == 1, "c6")

# The first type
for i in I:
  if(i>0):
    for j in J:
      if(j>0):
        for k in K:
          if(i!=j):
            mdl.add_constraint(u[i] - u[j] + (5*x[(i,j,k)]) <= (5-demand[i]), "c6")

# The second type
for i,i1 in enumerate(I_1):
  for j in J_1:
    for k in K:
      if(i!=j):
        mdl.add_constraint(u[i] - u[j] + (5*x[(i,j,k)]) <= (5-demand[i1]), "c6")

# The solve section
mdl.print_information()
print("------------------------------------")
m = mdl.solve(log_output=True)
print("------------------------------------")

x_solution = {}
for i in I:
  for j in J:
    for k in K:
      if(x[(i,j,k)].solution_value > 0.5):
        x_solution[(i,j,k)] = x[(i,j,k)].solution_value
print(x_solution)
print("------------------------------------")

# print(solution)
objective = m.get_objective_value()
print(objective)
print(m.solve_status) # Returns if the solution is Optimal or just Feasible