Comparative Study of Algorithms Metaheuristics Based Applied to the Solution of the Capacitated Vehicle Routing Problem

3.1 Savings based heuristics

There are several types of heuristic methods to solve the VRP, which are addressed extensively in Braekers et al. [3], trying to generate broad, nonexclusive categories, where the best-known heuristics are located to solve this problem; among them one of the most used and popular algorithms is the one of Clarke & Wright, and that has had contributions from different authors [7].

This algorithm is based on successively combining subtours until a Hamiltonian cycle is obtained, of which the subtours have a common node or vertex called base or initial.

The method can be described as follows:

• Having a solution of two different routes 0…i0and 0j…0can be joined by creating a new route 0…ij…0.

• The distance savings obtained by the union is

In Eq. (1) sijis the savings on the total distance traveled if the two routes 0…i0and 0j…0are joined.

• An initial solution is started in this algorithm, and the unions that give greater savings are made as long as they do not violate the restrictions of the problem.

• When the maximum saving is negative, the combinations of the routes will increase the distance traveled, but the amount of routes in the solution will decrease; depending on the characteristics of each problem, this can generate circular or radial routes that can be avoided by placing a reference value λ, which penalizes the union routes with distant customers. Saving is proposed as

In Eq. (2) sijrepresents penalized savings with the weight λin the total distance traveled if the two routes 0…i0and 0j…0are joined, which prevents when possible, merging routes with nodes far apart.

3.1.2 Example of application of the savings algorithm

A company wants to solve the problem of routing and design of the fleet of its product xto its 10 customers in the city and has a homogeneous fleet of trucks with capacity for 100 units of x product, with locations and demand shown in Table 1.

Table 1.

Cardinal coordinates of storage and customers with their demand.

There are details in the Cartesian coordinates of the warehouse and each customer with the demand, while in Figure 1, the position of each customer and the warehouse is shown.

1. Step 1: The matrix of Euclidean distances between each pair of nodes is calculated: Table 2 shows the distance matrix between all customers along with the warehouse. This matrix is symmetrical, that is, it has the same distance to go from client ito client jand vice versa. Point 0 has been considered for the warehouse (whs).

Cij	whs	c1	c2	c3	c4	c5	c6	c7	c8	c9	c10
whs	—	25.46	19.80	19.65	7.07	25.00	15.81	23.60	29.15	24.19	11.05
c1	25.46	—	5.66	23.02	26.42	33.06	40.22	37.22	54.20	6.71	33.62
c2	19.80	5.66	—	19.65	21.02	29.61	34.67	33.00	48.60	7.28	28.18
c3	19.65	23.02	19.65	—	26.00	42.20	34.00	14.87	39.45	17.12	18.97
c4	7.07	26.42	21.02	26.00	—	18.03	14.14	30.61	31.62	26.93	17.09
c5	25.00	33.06	29.61	42.20	18.03	—	25.00	48.60	45.00	36.88	34.89
c6	15.81	40.22	34.67	34.00	14.14	25.00	—	32.20	20.00	39.81	17.09
c7	23.60	37.22	33.00	14.87	30.61	48.60	32.20	—	29.61	31.78	15.26
c8	29.15	54.20	48.60	39.45	31.62	45.00	20.00	29.61	—	51.62	21.26
c9	24.19	6.71	7.28	17.12	26.93	36.88	39.81	31.78	51.62	—	30.48
c10	11.05	33.62	28.18	18.97	17.09	34.89	17.09	15.26	21.26	30.48	—

Table 2.

Distance matrix.

1. Step 2: Once the distance matrix is obtained, the savings are calculated. For the savings matrix, no row or column is placed for the warehouse.

   For example, the savings between customer 1 and customer 2 is

   s12=c1bdg+cbdg2−c12E3
   s12=25.46+19.80−5.66=39E4

   In Table 3 all the savings are shown; in the same way the matrix is ​​symmetric.

Sij	c1	c2	c3	c4	c5	c6	c7	c8	c9	c10
c1	—	39.60	22.08	6.11	17.40	1.04	11.84	0.41	42.93	2.89
c2		—	19.80	5.85	15.18	0.94	10.40	0.35	36.71	2.67
c3			—	0.72	2.44	1.46	28.38	9.36	26.72	11.72
c4				—	14.04	8.74	0.06	4.60	4.33	1.03
c5					—	15.81	0.00	9.15	12.31	1.16
c6						—	7.21	24.97	0.19	9.77
c7							—	23.14	16.01	19.38
c8								—	1.72	18.94
c9									—	4.75
c10										—

Table 3.

Savings matrix.

1. Step 3: The route 0i0is built for each client. In Figure 2 each route is shown from the warehouse to each customer and back to the warehouse.

1. Step 4: Savings are organized from the highest to lowest.

   In the list of savings to choose, only the savings that can be chosen are considered.

   When the list is prepared with all the savings, those savings that one or both clients have already considered in a previous route are discarded.

2. Step 5: To assemble the routes, the restrictions are considered; for this example the only restriction is the capacity of the truck that does not exceed 100 units of the product x.

   For the first route, the highest savings are chosen and placed in the form 0i0; in this case 0is the warehouse (whs), as the savings are chosen to add them to the existing route or create a new one; the demands of each client are added, and the route is closed when the sum of the demands is equal to or less than the capacity of the truck 100 units, but when adding one more client, the demand exceeds the capacity, and you can no longer choose that customer.

   The composition of the routes is displayed step by step in Figure 3, and the complete route diagram is shown in Figure 4. Below is the composition of the routes with the demands.

Route 1: whs, c1, c9, whs	Route demand 1: 17 + 14 = 31
Route 1: whs, c2, c1, c9, whs	Route demand 1: 31 + 25 = 56
Route 2: whs, c3, c7, whs	Route demand 2: 10 + 28 = 38
Route 1: whs, c2, c1, c9, c3, c7, whs	Route demand 1: 56 + 38 = 94
Route 2: whs, c6, c8, whs	Route demand 2: 15 + 36 = 51
Route 2: whs, c6, c8, c10, whs	Route demand 2: 51 + 24 = 75
Route 2: whs, c5, c6, c8, c10, whs	Route demand 2: 75 + 20 = 95
Route 3: whs, c4, whs	Route demand 3: 29

Consequently, Clarke & Wright algorithm determines a solution for the routing problem in which the distance traveled is 204.20 units in length.

3.2 Heuristic method of assigning first, routing after

Sweep heuristics are the best-known method of assigning first, routing later.

This method is solved in two phases. First, groups of customers called clusters are created considering the capacity constraints of the vehicles, and second for each cluster, a route is generated that visits all customers.

In sweep heuristic, clusters are created by turning a half-straight in the central tank from the horizontal counterclockwise; after that the customers are incorporated into the mentioned group until the maximum capacity restriction of the vehicles is met.

This heuristic is used to find solutions to geographical problems, that is to say in which the nodes or vertices correspond to a point in the plane. It is assumed that the location of each client ican be represented through its polar coordinates riθihaving a single central deposit. It defines

3.2.2 Example of sweep heuristics

We will take the example of the savings algorithm.

1. Step 1: Formula (5) is used to obtain the polar coordinate table, where θiis expressed in radians and riis the directed distance. Table 4 shows the polar coordinates for each client i. The change of polar coordinates is displayed (Figure 5).

	ri	θi
c1	25.46	−2.36
c2	19.80	−2.36
c3	19.65	2.88
c4	7.07	−0.79
c5	25.00	−0.93
c6	15.81	0.32
c7	23.60	2.21
c8	29.15	1.03
c9	24.19	−2.62
c10	11.05	1.66

Table 4.

Polar coordinates for each i.

1. Step 2: It is sorted by θifrom least to greatest, as seen in Table 5.

	ri	θi
c9	24.19	−2.62
c2	19.80	−2.36
c1	25.46	−2.36
c5	25.00	−0.93
c4	7.07	−0.79
c6	15.81	0.32
c8	29.15	1.03
c10	11.05	1.66
c7	23.60	2.21
c3	19.65	2.88

Table 5.

Ascending ordering of customers.

1. Step 3: To elaborate the routes, it is done in two phases, the first one where the clients are grouped by the sweep method and the second one where a TSP is resolved (step 4).

   For the sweep method, the angles from the smallest to the largest are chosen, and it moves counterclockwise.

   As can be seen in Table 5, customers are already sorted in ascending order by their angular polar coordinate, and customers are chosen until they fail to comply with the capacity restriction of the truck that is 100 units of product x. Considering this the routes are as follows:

Route 1: whs, c9, c2, c1, c5, whs	Route demand 1: 14+25+17+20 = 76
Route 2: whs, c4, c6, c8, whs	Route demand 2: 29+15+36 = 80
Route 3: whs, c10, c7, c3, whs	Route demand 3: 24+28+10 = 62

In Figure 6, the sweeps are visualized starting with client 9 that has the greatest angle, thus grouping them in zones in this case by colors and within each one for their best distance. The sweep groups customers do not violate the restriction of the truck.

In Figure 7, the solution is shown with three routes before the TSP is applied.

1. Step 4: In the second phase to each route already generated in the first, it is resolved by TSP, for this case with the nearest node.

   The routes are as follows:

Route 1: whs, c9, c1, c2, c5, whs	Route demand 1: 14+17+25+20 = 76
Route 2: whs, c4, c6, c8, whs	Route demand 2: 29+15+36 = 80
Route 3: whs, c10, c7, c3, whs	Route demand 3: 24+28+10 = 62

Thus, the sweep algorithm has a local solution for the routing problem in which the distance traveled is 222.36 units in length.

In Figure 8, the route diagram is displayed.

1. Step 5: Repeat step 2 where the customers already ordered from Table 5, continue to rotate the position until the first returns to be first, and for each rotation, steps 3 and 4 are performed, and after all iterations, the best is selected.

   Below is the iteration that had the best result. Table 6 shows the fifth iteration of nine where you start with client six (Figure 9).

	ri	θi
c6	15.81	0.32
c8	29.15	1.03
c10	11.05	1.66
c7	23.60	2.21
c3	19.65	2.88
c9	24.19	−2.62
c1	25.46	−2.36
c2	19.80	−2.36
c5	25.00	−0.93
c4	7.07	−0.79

Table 6.

Customer ordering—fifth iteration.

After performing step 3 in Figure 10, the grouping of customers is appreciated to not violate the restriction of the truck’s capacity.

In step 4, each grouping is resolved with a TSP, and the following routes are obtained:

Route 1: whs, c6, c8, c10, whs	Route demand 1: 15+36+24 = 75
Route 2: whs, c7, c3, c9, c1, c2, whs	Route demand 2: 28+10+14+25+17 = 94
Route 3: whs, c5, c4, whs	Route demand 3: 20+29 = 49

In Figure 10, the diagram of the routes of the fifth iteration is displayed, which obtained the best response.

The sweep algorithm determines a solution for the routing problem in which the distance traveled is 205.96 units in length, that is, a solution of lower quality than that obtained by the Clarke & Wright algorithm with 204.20 units in length.

4.1 GRASP

GRASP methods had their origins at the end of the 1980s in order to find a solution to problems of joint coverings, and in 1995 by Feo and Resende, this metaheuristic is of general purpose [8].

The word GRASP comes from the acronym of greedy randomized adaptive search procedures that would be something like search procedures based on voracious adaptive random functions.

GRASP has a multistart process in which each step has a construction and an improvement phase. In the construction phase, the constructive heuristic process obtains a good initial solution, which is improved in the second phase by a local search algorithm. The best of all solutions examined is saved as the final result.

There are many implementations of GRASP metaheuristics, including variants and hybridizations with other procedures such as variable neighborhood search or path relinking, with which this metaheuristic has proven to work very well in practice as demonstrated in Marti and Sandoya [9]. A simple scheme to represent the operation of this algorithm is as follows:

While (stop condition)

Construction phase:

• Choose a list of candidate elements.

• Have a restricted list with the best candidates.

• Select an item randomly from the restricted list.

Improvement phase:

• Perform a local search process based on the solution built until it can no longer be improved.

Update:

• If the solution obtained improves to the best stored, update it.

In the construction phase, a possible solution is built iteratively, considering an element in each step. In each iteration the choice of the next element to be added to the partial solution is determined by a greedy function, which examines the benefit of adding each of the elements according to the objective function and choosing the best one.

This metaheuristic works with a restricted list of the best candidates, which makes the best candidate randomly selected for each iteration of the construction phase.

In the improvement phase, the results that are obtained from the construction phase are not usually local optimal; therefore, a local search procedure is applied as post-processing to perfect the solution obtained.

Performing several iterations is a way of sampling the solution space.

4.2 Simulated annealing

The simulated annealing metaheuristics was introduced in the 1950s by Metropolis Hastings to be used in the field of statistical thermodynamics simulating cooling processes of a material.

In 1983 the method was refocused to solve combinatorial optimization problems of great complexity by Scott Kirkpatrick, C. Daniel Gelatt and Mario P. Vecchi, and independently in 1985 by Vlado Cerny. For its implementation ease, this metaheuristic had a great boom in the 1980s.

Simulated annealing has its procedure based on local search by environments that is characterized by an acceptance criterion of neighboring solutions that are adapted throughout its execution.

A temperature variable is used, T, that determines the extent to which neighboring solutions, worse than the currentn,can be accepted. This variable is about starting it with a high value, which is called the initial temperature, T0, which generates a high probability of accepting a nonimprovement movement. In each iteration the temperature decreases through a temperature cooling mechanism, α, having a smaller probability until approaching the optimal solution and reaching a final temperature, Tf. Costs also decrease as the temperature decreases, making it increasingly difficult to accept bad movements in search of the solution.

In each iteration a specific number of neighbors is generated, which can be fixed for the entire execution or depend on each iteration.

Each time a neighbor is generated, the acceptance criterion is applied to see if it replaces the current solution:

• If the neighboring solution is better than the current one, it is automatically accepted, as it would be done in the classic local search.

• If the neighboring solution is worse than the current one, there is still a chance that the neighbor will replace the current solution. This allows the algorithm to exit from local optimum, in which the classic local search would be trapped.

This model is given by the following structure:

Take an initial solution x

Take an initial temperature T

While (not frozen)

Perform Ltimes

1. Take x’from Nx

2. d = f(x′) – f(x)

3. If d<0dox=x′

4. Ifd>0dox=x′with p=e−d/T

5. Take action of the cooling mechanism T=rT

The following parameters are determined:

1. Initial temperature: it is established by doing a series of tests to reach a certain fraction of accepted movements.

2. Cooling speed r.

3. LengthLthat must be proportional to the expected size of Nx.

4. When the cooling sequence ends, it is frozen cont=cont+1when a temperature is completed and the percentage of movements accepted is less than MinPercent.cont=0when the best stored solution is improved.

5.1 Test instances

The cases to be evaluated are divided into three groups classified by the type of client with 10 examples each. Next, some tables will be presented, which have the name of instance, the truck’s capacity in columncap, the number of customers in column n, the number of vehicles to be used in columnk, and the optimal solution in column opt.

1. Clustered clients, as shown in Table 7, belong to the Augerat B set in 1995 [10] and specify that the coordinates are points between [0,100] × [0,100] in the grid that are chosen to form neighborhood groups (NC) closest, where k≤NC−1. The demands have a uniform distribution U130; however, n/10was multiplied by3.

2. Random clients, as shown in Table 8, belong to the Augerat set A in 1995 and specify that the coordinates are points between [0,100] × [0,100] placed randomly. The demands have a uniform distribution U130; however, n/10was multiplied by 3.

3. Clustered and random clients, as shown in Table 9, belong to the Augerat set X in 1995 and specify that the coordinates are points between [0,1000] x [0,1000] that are grouped and placed randomly, where kis the minimum feasible number of vehicles.

6.1 Test cases resolved by heuristics

Clarke & Wright heuristics have better quality solutions, solving problems where customers with a small nare grouped.

Table 10 shows that for the group of clients with gathered positions, the gap is 3.63%; for the positions of random clients, the gap is 5.18%; and in less effective way for customers with grouped and random positions, it has a gap of 6.55%.

It also compares the number of vehicles kthat were obtained when solving each case against the optimal solution, and it was obtained that 12 cases had a vehicle more than the optimum B-n41-k6, B-n45-k6, A-n33-k6, A-n34-k5, A-n38-k5, A-n44-k6, X-n106-k14, X-n115-k10, X-n134-k13, and X-n139-k10, and two cases had three more than the optimal vehicles X-n101-k25 and X-n125-k30.

Sweep heuristics are more effective in solving problems where customers with a small nare grouped together. However, the difference in the average gap between random and grouped customers is short.

Table 11 shows that for the group of customers with grouped positions, the gap is 8.68%; for random customer positions; the gap is 8.85%; and in a less effective way for customers with grouped and random positions, it has a gap of 17.00%.

It also compares the number of vehicle numbers kthat were obtained when solving each case against the optimal solution, and it was obtained that seven cases had a vehicle more than the optimal B-n45-k6, A-n38-k5, A-n39-k5, X-n115-k10, X-n129-k18, X-n134-k13, and X-n139-k10; one case had five more vehicles than the optimum X-n101-k25, and one case had six more vehicles than the optimal X-n125-k30.

In each group of clients, the sweep heuristic obtained better answers than the Clarke & Wright heuristics with 30% in the group of clients with a grouped position, 30% in the group of clients with a random position, and 20% in the group of clients with grouped and random position. In other words, Clarke & Wright heuristics are superior with 70% in the first two groups of clients and with 80% in the last group of clients.

A comparison among the values of the Distance traveled in the solution found by the heuristic, the optimal solution and the GAP for each one of the considered test instances is shown in Table 12.

On the other hand, Table 13 shows a summary of the minimum, maximum and average gap for each of the three classes of problems considered: Clustered, Random and Clustered, and Random.

6.2 Test cases resolved through metaheuristics

The GRASP metaheuristics are based on a previous solution for which Clarke & Wright heuristic responses were selected since their responses are of better quality than the sweep heuristics.

The following parameters were considered for its implementation:

• α = 0.5

• Number of iterations = 10,000

GRASP’s metaheuristics are more effective in solving problems where customers with a small nare grouped together.

Table 14 shows that for the clients with grouped positions; the gap is 3.09%; for the positions of random clients, the gap is 4.38%; and less effectively for customers with grouped and random positions, it has a gap of 5.97%.

On average metaheuristic GRASP based is better than Clarke & Wright heuristics by 0.53%, 0.77%, and 0.38% in the solutions of the positions of the grouped, random, and grouped and random clients, the group of clients with random positions being those that obtained a greater improvement in the quality of the solutions.

The simulated annealing metaheuristics start from a previous solution for which Clarke & Wright heuristic responses were selected since their responses are of better quality than the sweep heuristics.

The following parameters were considered for its implementation:

• Current temperature = 250

• Final temperature = 10

• Cooling coefficient = 0.8

• Number of iterations = 10,000

The simulated annealing metaheuristic is more effective in solving problems where customers with a small nare grouped together.

Table 15 and 17 shows that for clients with grouped positions, the gap is 3.13%; for the positions of random clients, the gap is 4.41%; and less effectively for clients with grouped and random positions, it has a gap of 6.26%.

On average simulated annealing heuristics based is better than Clarke & Wright heuristics by 0.52%, 0.71%, and 0.21% in the solutions of grouped, random, and grouped and random clients’ positions, the group of clients with random positions being those that obtained a greater improvement in the solutions quality.

Within each group of clients, the simulated annealing metaheuristics obtained better answers than the GRASP metaheuristics with 30% in the group of clients with a grouped position, 50% in the group of clients with a random position, and 40% in the group of clients with grouped and random position. That is, the GRASP metaheuristic is superior with 70% in the first group and with 60% in the third group of clients and is equal with 50% in the second group of clients.

A comparison between the minimum and maximum gap within each group is established in Table 16 test results with clustered clients have better results. Therefore, obtaining a minimum gap of 0%, that is, in the case of B-n34-k5, the GRASP metaheuristic obtained the optimal solution (Table 17).