Using CPLEX and python for finding an exact solution for the CVRP
After seeing the formulation of our main issue which is CVRP ,
Now let’s move to the serious work and implement a good program that allow us to generate an exact solution for the CVRP
First of all CPLEX IBM API is an optimizer solves integer programming problems and very large linear programming problems . It has a modeling layer called Concert that provides interfaces to the C++, C, and Java languages. There is a Python language interface based on the C interface. Additionally, connectors to Microsoft Excel and MATLAB are provided. Finally, a stand-alone Interactive Optimizer executable is provided for debugging and other purposes.
The rest of tools that we will be used is described in the table below
After describing our tools that we need to implement an exact solution for the CVRP we will focus on the development by explaining each iteration below:
- Making random numbers with the package Numpy and Creating an object called “rnd” for generating random numbers.
- Initializing our parameters witch are:
n: Number of customers.
Q: Maximum Capacity of the Vehicle.
N: The set of nodes.
V: The set of nodes + the node 0 which is the depot.
d: A collection that contains the demand of each node. - Initializing a random location for each node on a plot with two dimensions (x=200, y=100).
- Plotting the initial view of our random nodes.
- Initializing our set of arcs A and their cost c (c = The set of distances between each two nodes) . The cost in our case is the distance between each two nodes by calculating the differences between their locations in the plot.
- Calling the CPLEX API, initializing our decisions variables, constraints and printing the solution.
- Lastly, plotting our solution.
Coding the exact solution
#Importing the package numpy as np.import numpy as np# "rnd" is an object that generate random numbers.rnd = np.random#"seed(0)" is a method that reset (every time), the same random set of numbers.rnd.seed(0)# Number of customers.n = 10# Maximum capacity of the vehicle.Q = 15#The set of nodes without the depot.N = [i for i in range(1,n+1)]# The set of nodes + the depot.V = [0]+ N# A collection that contains the demand of each node.d ={i:rnd.randint(1,10)for i in N}# Generating random numbers between (0 and 15) * 200.loc_x = rnd.rand(len(V))*200# Generating random numbers between (0 and 15) * 100.loc_y = rnd.rand(len(V))*100#Importing the package matplotlib.pyplot as plt.import matplotlib.pyplot as plt#Plotting the n nodes without the node 0 (depot) and chose the color blue for each node.plt.scatter(loc_x[1:],loc_y[1:],c='b')# Associating and plotting each demand in the right of each blue node (customer).for i in N:plt.annotate('$d_%d=%d$'%(i,d[i]),(loc_x[i]+2,loc_y[i]))#Ploting the node 0, chosing the red like its color and the square form like a marker.plt.plot(loc_x[0],loc_y[0],c='r' ,marker='s')#Showing the Initial plot.plt.show()#Intializing the set of arcs A.A = [(i,j) for i in V for j in V if i!=j]#Calculating the distance between each node.c= {(i,j):np.hypot(loc_x[i]-loc_x[j],loc_y[i]-loc_y[j]) for i,j in A}#Importing the docplex.mp.model from the CPLEX as Modelfrom docplex.mp.model import Modelmdl = Model('CVRP')#Initializing our binary variable x_i,jx=mdl.binary_var_dict (A,name='x')#Initializing our cumulative demand uu=mdl.continuous_var_dict (N,ub=Q ,name = 'u')#Initializing the objectif functionmdl.minimize(mdl.sum(c[i,j]*x[i,j]for i,j in A))#Initialzing the first constraintmdl.add_constraints(mdl.sum(x[i,j]for j in V if j!=i)==1 for i in N)#Initialzing the second constraintmdl.add_constraints(mdl.sum(x[i,j]for i in V if i!=j)==1 for j in N)#Initialzing the third constraintmdl.add_indicator_constraints_(mdl.indicator_constraint(x[i,j],u[i]+d[j]==u[j])for i,j in A if i!=0 and j!=0)#Initialzing the fourth constraintmdl.add_constraints(u[i]>=d[i] for i in N)#Getting the solutionsolution =mdl.solve(log_output=True)#Printing the solutionprint(solution)#Identifing the active arcs.active_arcs =[ a for a in A if x[a].solution_value> 0.9]plt.scatter(loc_x[1:],loc_y[1:],c='b')for i in N:plt.annotate('$d_%d=%d$'%(i,d[i]),(loc_x[i]+2,loc_y[i]))for i,j in active_arcs :#Coloring the active arcsplt.plot([loc_x[i],loc_x[j]],[loc_y[i],loc_y[j]],c='g',alpha=0.3)plt.plot(loc_x[0],loc_y[0],c='r' ,marker='s')#Plotting the solutionplt.show()
Interpreting and verifying the reliability of the results.
In this part we will verify our constraints if they are satisfied or not form the figure (A plot for the exact solution of the CVRP (10 costumers))that shows that there is no sub-tours so we can conclude that the two first constraints are satisfied with success.
For the rest of the constraints we will use the output below of the program that shows that always the cumulative demand ui + dj = uj and do not exceeds the total capacity Q
The program was tested on 10 costumers,as we have already seen the time elapsed to generate the exact solution equal to 1.00 sec. However when we increased the number of costumers to 11 the program took 114.41 sec (I let you try it by yourselves). We notice that we observe a difference that equal to 113.41 sec so we can conclude and affirm that the complexity of CVRP is exponential.
In a large scale problem the implementation of an exact method could not be a good real solution given the complexity of the CVRP,As we seen in our review of literature the researchers used to work with heuristics and meta-heuristics methods to reduce the complexity in the other hand they find approximative results.
In our future work we will implement an hybrid method CMSA algorthim using the java language on the CVRP by combining this exact method with an heuristic called local search for the purpose to find a best approximative results.
