- Email: [email protected]
An optimization algorithm for solving the rich vehicle routing problem based on Variable Neighborhood Search and Tabu Search metaheuristics
Journal of Computational and Applied Mathematics (
)
–
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
An optimization algorithm for solving the rich vehicle routing problem based on Variable Neighborhood Search and Tabu Search metaheuristics Juan Antonio Sicilia a,∗ , Carlos Quemada b , Beatriz Royo c , David Escuín d a
Escuela de Ingeniería, Universidad Internacional de La Rioja, 26002 Logroño, Spain
b
Virginia Polytechnic Institute and State University in Blacksurg, Virginia, USA
c
Department of Mechanical Engineering, University of Zaragoza, 50018 Zaragoza, Spain
d
Instituto Tecnológico de Aragón, María Luna 1, 50018 Zaragoza, Spain
article
info
Article history: Received 13 October 2014 Received in revised form 9 March 2015 Keywords: Rich vehicle routing problem Optimization algorithm Metaheuristic methods Load balance
abstract This paper presents a novel optimization algorithm that consists of metaheuristic processes to solve the problem of the capillary distribution of goods in major urban areas taking into consideration the features encountered in real life: time windows, capacity constraints, compatibility between orders and vehicles, maximum number of orders per vehicle, orders that depend on the pickup and delivery and not returning to the depot. With the intention of reducing the wide variety of constraints and complexities, known as the Rich Vehicle Routing Problem, this algorithm proposes feasible alternatives in order to achieve the main objective of this research work: the reduction of costs by minimizing distances and reducing the number of vehicles used as long as the service quality to customers is optimum and a load balance among vehicles is maintained. © 2015 Elsevier B.V. All rights reserved.
1. Introduction A vast number of problems from Applied Science including engineering can be brought by means of solving a nonlinear equation using mathematical modeling [1–3]. One of that problems, in concrete the capillary transport of goods problem, is studied in this paper. The logistics of the capillary transport of goods has an essential urban dimension. On account of this, distribution requires efficient systems and the process between warehouses and customers must be efficient and clean. Therefore, it is necessary to optimize urban logistics efficiently and improve the connections between urban and interurban freight transport in order to ensure efficient distribution [4]. In the European Union, over 60% of the population lives in urban areas, generating just under 85% of the EU’s Gross Domestic Product (GDP) [5]. In addition to this, pickup and delivery actions in urban areas represent about 40% of the total cost of transport activities carried out at home. These expenses are further increased by the reduction of stocks and size of goods and by the increase in the number of requests.
∗
Corresponding author. E-mail addresses: [email protected] (J.A. Sicilia), [email protected] (C. Quemada), [email protected] (B. Royo), [email protected] (D. Escuín). http://dx.doi.org/10.1016/j.cam.2015.03.050 0377-0427/© 2015 Elsevier B.V. All rights reserved.
2
J.A. Sicilia et al. / Journal of Computational and Applied Mathematics (
)
–
Due to this, there is great potential for freight transport in these areas. Therefore, it is necessary to find logistics solutions for urban freight distribution that are consistent with the restrictions imposed by the law in order to protect the interests of the public. This paper presents an optimization algorithm that solves a Rich Vehicle Routing Problem (RVRP) and arises from a research project carried out for an important Spanish distribution company. The main goal of this project is to manage its resources in urban areas by reducing costs caused by inefficiency and ineffectiveness as much as possible. Therefore, this research work consists in reducing the constraints and complexities currently encountered in the urban distribution of goods by using the Variable Neighborhood Descent (VND), the General Variable Neighborhood Search (GVNS) and Tabu Search (TS) metaheuristic methods. Most problems in industry, particularly in city logistics, are multi-objective in nature. This makes it difficult to apply classical routing models to real-life problems because these models are developed with the single objective of minimizing the cost of the solution [6]. However, in real life there may be different requirements related to a single tour such as balancing workloads, time, distance, cost, etc. That need to be taken into consideration simply by adding new constraints. This paper is organized as follows. In Section 2, the literature review is presented. Subsequently, Section 3 presents a formal mathematical definition of the RVRP. Section 4 proposes an optimization algorithm to solve this problem using metaheuristic methods. Section 5 shows the computational calculation results of this algorithm and these results are discussed in Section 6. Finally, conclusions are set out in Section 7. 2. Literature review Over the last years, several variants of multiconstrained Vehicle Routing Problems have been studied, forming a class of problems known as Rich Vehicle Routing Problems. A comprehensive and relevant taxonomy for the RVRP literature and an attempt to provide a definition of RVRPs can be found in [7]. Several applications of such problems are close to real world as the work shown by [8] to deal with animal welfare and global constraints to handle production and inventory. Other authors have proposed to refer RVRPs as Multi-Attribute Vehicle Routing Problems (MAVRP). The work presented by [9] considers a multi-attribute vehicle routing problem derived from a real-life milk collection system. This problem is characterized by the presence of a heterogeneous fleet of vehicles, multiple depots, and several resource constraints. [10] claims that the attributes of vehicle routing problems are additional characteristics or constraints that aim to better take into account the specifications of real applications. It then takes a closer look at the concepts of 64 remarkable meta-heuristics, selected objectively for their outstanding performance on 15 classic MAVRP with different attributes. This cross-analysis leads to the identification of winning strategies in designing effective heuristics for MAVRP. This is an important step in the development of general and efficient solution methods for dealing with the large range of vehicle routing variants. The problem RVRP considered in this paper is based on the optimization of the freight distribution in large urban areas taking into consideration the real characteristics presented in companies and trades. Given a fleet of vehicles and a set of customers dispersed over a geographic area, the system must be able to attain a feasible solution that leads to a decrease of the total cost of the problem. The purpose is to reach a solution with the minimum number of vehicles and minimizing the total distance traveled. The model of the problem depends on the following requirements:
• CVRP (Capacitated Vehicle Routing Problem): Vehicles have a maximum limit of capacity and are homogeneous [11]. • VRPTW (Vehicle Routing Problem with Time Windows): Each customer must be serviced in a set time interval [12]. • SDVRP (Site-Dependent Vehicle Routing Problem): Certain orders can only be transported by specific vehicles. Therefore, the fleet must be heterogeneous [13].
• VRPPD (Vehicle Routing Problem with Pickups and Deliveries): Customers may receive and send goods [14]. • VRPB (Vehicle Routing Problem with Backhauls): The goods that leave the depot must be delivered to customers and the goods that are picked up from customers must be transported to the depot [15].
• OVRP (Open Vehicle Routing Problem): Vehicles are not required to return to the depot after visiting the last customer and the route may finish near to the drivers home [16].
• Load balance among routes: It is important to balance the workload among vehicles since the competitiveness of a transport company depends not only on minimizing costs and distances to obtain higher profit margins but also on its ability to treat employees fairly [17]. 3. Mathematical formulation Prior to describing the mathematical formulation, it is important to define what an order is in the context of this paper. An order is a request made by a customer for the transport of goods. There are three different kinds of orders: one-pickup orders, one-delivery orders or pickup–delivery orders at dependent nodes. These nodes have to be passed through in a specific and proper sequence by the same vehicle with the time window imposed by the customer. Each vehicle has a maximum of a day to complete the service. Since a customer might request more than one order, each customer can be visited more than once by different vehicles. In addition to this, there are dependent orders where vehicles have to pick up goods at a specific node and deliver them to another node without the need to return to the warehouse.
J.A. Sicilia et al. / Journal of Computational and Applied Mathematics (
)
–
3
Furthermore, vehicle capacity constraints and compatibilities between orders and vehicles might lead certain vehicles to be excluded from allocating orders or certain orders to be transported only by specific vehicles. This situation requires a heterogeneous fleet in capacity and technical properties. The vehicles may have different capacities and types of loads (normal, refrigerated, isotherm, etc.). Two different location points are assigned to each vehicle. The end node may be the warehouse or the driver’s home and the initial node corresponds to the warehouse regardless of where the vehicle have finished the previous day. Every order has a service time (load/unload time of goods) and a time window in which the operation must be carried out. It is important to point out that this time window does not prevent any vehicle from arriving before the lower time limit and this might cause losses due to the inactivity of a vehicle during a period of time. The intention is for the arrival of each vehicle to be within that time interval. Another important factor that should be taken into account is the load balance among drivers. A maximum number of orders per vehicle is established so a route is carried out independently of its costs. In addition to this, another balance between the distance and the driver actions must be attained in all routes. Since the number of orders might vary each day, the number of customers (geographic nodes) is not fixed and the distances have to be exact in urban areas. A Geographic Information System (GIS) is used to calculate the precise distance from the GPS position of each node to be visited. This distance is only calculated among orders that may go on the same route. Previously, different sets of orders that may be compatible with each other are analyzed. The aim of this work is to create optimized routes that minimize the total cost of the transport operations by minimizing the total distance, reducing the number of vehicles and eliminating inactivity times. Furthermore, a load balance among vehicles must be established, the available resources must be maximized and there must be a commitment to the service quality. One route is defined as an ordered sequence of orders that the same vehicle has to dispatch. Regardless of cost, time is the most important factor to take into account. Therefore, it is necessary to fulfill the pickup and delivery schedule. The RVRP can be described as follows. Let G = (N , A) be a complete graph with N = {1}∪ O ∪{o + 1} the set of nodes and O = {2, . . . , o} the set of orders to be served. The nodes {1} and {o + 1} represent the depot and are added to the network to allow empty tours. The aim is to design a set of routes for a heterogeneous fleet of V vehicles (each has a capacity, a load type and a limited number of orders to be served), servicing a set of O orders. Each arc (i, j) ∈ A has a cost cij . The interval [ai , bi ] denotes the time windows for order i ∈ O and ti represents the time service. A vehicle cannot start servicing node i before ai and after bi . A vehicle can arrive before ai and wait for service. Drivers homes are represented by some nodes and vehicles have to serve them just before returning to the depot. As explained above, it is possible to distinguish three kinds of nodes:
• Deliveries: vehicles depart from depot with goods. • Pickups: vehicles pick up goods at these nodes and transport them to depot. • Pickup-and-delivery dependents: vehicles have to pick up goods at pickup node and deliver them at delivery node. On each route the total quantity of goods picked up and delivered cannot exceed vehicle capacity. Each node i represents both a delivery node, di , and a pickup node, pi . The pickup node of a delivery node i is represented by ri . In addition to this, a node may have a demand to be picked up and delivered. Different nodes i ∈ N may correspond to the same geographical location. The set of available vehicles is represented by V , being Qk the capacity of vehicle k ∈ V , Tk the maximum number of nodes that vehicle k can serve. hk represents the last node to be served (driver’s home) before returning to depot, hk ∈ O. Furthermore, mik denotes the compatibility between node i and vehicle k, ∀i ∈ N, ∀k ∈ V . Node i may be served by vehicle k if and only if mik = 1, mik = 0 otherwise. The RVRP is formulated as follows: Decision variables: xkij =
if arc (i, j) ∈ A is traveled by vehicle k ∈ V otherwise.
1 0
(1)
Objective function: min
cij xkij −
k∈V i∈N j∈N
xk1hk M .
(2)
k∈V
The first term represents the accumulated costs for orders served by vehicles; the second term represents the number of utilized vehicles. Subject to:
xkij = 1 ∀i ∈ O
(3)
j∈N k∈V
i∈N
xkih −
j∈N
xkhj = 0 ∀h ∈ O, ∀k ∈ V
(4)
4
J.A. Sicilia et al. / Journal of Computational and Applied Mathematics (
)
–
xkii = 0 ∀i ∈ V , ∀h ∈ N
(5)
xko+1i = 0 ∀k ∈ V , ∀i ∈ N
(6)
xkhk o+1 = 1 ∀k ∈ V
(7)
(8)
xk1j = 1 ∀k ∈ V
j∈N
ski + ti − skj + cij ≤ (1 − xkij )M
∀k ∈ V , ∀i ∈ N , ∀j ∈ N
sk1 = 0 ∀k ∈ V ai ≤
ski
(9) (10)
= 0 ∀k ∈ V , ∀i ∈ N
(11)
∀k ∈ V , ∀i ∈ N
(12)
bi ≥ ski = 0
mik + mjk ≥ 2xkij
∀k ∈ V , ∀i ∈ N
lki ≤ Qk
∀k ∈ V , ∀i ∈ O, ∀j ∈ O
di xkij = lk1
(13) (14)
∀k ∈ V
(15)
i∈N j∈N
pi xkij = lko+1
∀k ∈ V
(16)
i∈N j∈N
lki − dj + pj − lkj ≤ (1 − xkij )M Tk + 1 ≥
xkij
∀k ∈ V , ∀i ∈ N , ∀j ∈ N
∀k ∈ V
(17) (18)
i∈N j∈N
1−
xkrj h ≤ (1 − xkij )M
∀k ∈ V , ∀i ∈ N , ∀j ∈ O
(19)
h∈N
skri ≤ ski
∀k ∈ V , ∀i ∈ O.
(20)
The domain of decision variables is set by (1). The objective function (2) states that the distance traveled and number of vehicles should be minimized, M being a sufficiently large constant. (3) ensures that one and only one node can be visited after an order. Flow conservation is ensured by constraints (4) and (5) avoids cycles. Constraint (8) states that all vehicles leave the depot while constraint (7) denotes all vehicles returning to depot after servicing the home drivers. Constraint (6) ensures no vehicles leave the virtual depot and vehicles end their routes at node o + 1. Demand constraint is ensured by (14)–(17) with lki to denote the load of the vehicle k at node i after serving node i. Compatibility between vehicles and orders is represented by (13). Time window constraints are ensured by constraints (9)–(12) with ski to denote when the vehicle k starts the service at node i. Constraints (11) and (12) ensure that node i cannot be served before ai and after bi . Constraint (18) ensures that the number of customers served by vehicle k is not greater than the maximum permitted. Constraints (19) and (20) represent precedence constraints for pickup-and-delivery dependents. 4. Optimization algorithm The main difficulty of the RVRP problem is that if the number of elements to be combined is relatively large, the combinatorial result grows exponentially turning it into an NP-hard problem. On account of this, it is necessary to apply optimization techniques in order to obtain a high-quality solution in reasonable computational times. In this section an efficient algorithm to solve the problem is presented. The RVRP is solved by means of a powerful route construction method that starts with the determination of a feasible initial solution and the metaheuristics VND, GVNS and TS based on various elementary neighborhood operators to obtain an improved solution. 4.1. Route construction method At the beginning of this process there are several distribution warehouses, a heterogeneous fleet of vehicles and a set of orders to be distributed. The algorithm will generate as many routes as vehicles are available, taking into consideration the constraints and compatibilities.
J.A. Sicilia et al. / Journal of Computational and Applied Mathematics (
)
–
5
The route construction method consists of assigning orders to vehicles depending on certain parameters and taking into consideration the limits of load balance between routes. This procedure has been implemented randomly to avoid obtaining of the same solutions and to diversify the total space of the search. In the beginning there is no route and the system proceeds to create routes by choosing the first order according to six available criteria. This does not mean that this order is the first one of the final route. It has been selected only to start the insertion process. In this insertion process, orders are inserted until the route is not feasible or until it has reached the maximum number of orders allowed per route. The process finishes when there are no orders dispatched in the routes. This is the initial solution that goes to the local improvements phase. To select the first order i of an empty route, the following insertion criteria might be considered:
• The node i with the biggest value of c1i + ti + cio+1 . That is, the node with the latest arrival time to the depot without • • • •
considering the home driver yet. The node i with the smallest value of bi : earliest close of the time window of the node. The node i with the shortest value of bi − ai : the tightest time window. The node i with the smallest value of cio+1 : the effective distance from the customer to depot. Best average of the above criteria.
The last one is the most common criterion because it includes the rest of criteria. In this case, all criteria are calculated in order to establish a list with orders that obtains the best value in terms of the rest of criteria. The insertion process of the rest of orders is based on the following criteria:
• The largest time margin: it is the extra time available to distribute other orders taking into consideration the time taken to distribute orders that have already been inserted.
• The shortest distance of the route: choice of a nearby customer. The order will be inserted in the position of the route with minimum distance. 4.2. Local improvements and metaheuristics Local improvement operators are methods applied to a set of one or more routes in order to improve the evaluation function. The basic process is made up of movements and exchanges of order segments. It is necessary to apply the operators to improve specific aspects of the problem in the majority of the calculations. These operators are the well-known intraroute improvements applied to a single route and that are computationally very fast, OR-OPT [18] and IOPT [19], and the inter-route improvements that exchange segments of orders between two routes, ICROSS [20] and 2-OPT* [21]. The OROPT operator moves one segment of the route into other position. IOPT means that the sequence is reversed. The CROSS exchanges two (one being potentially empty) customer sequences s1 and s2 from two routes. The customer sequences s1 and s2 can eventually be reversed in this process to produce a slightly larger neighborhood, called I-CROSS. 2-OPT refers to the fact that there is an exchange or swap of a pair of edges. Once the initial solution, s, has been obtained, this is improved through the application of a set of metaheuristics. A sequence of operators with their respective configurations must be specified in order to define a neighborhood search. The metaheuristic process consists of using a neighborhood that exchanges segments. The final state of each operator of a neighborhood determines the subsequent execution of the following neighborhood. Therefore, a good configuration of operators enables the obtaining of one solution or another [22]. 4.2.1. Variable Neighborhood Descent VND is a variant of Variable Neighborhood Search that consists of replacing the current solution iteratively with the result from the local solution if there is an improvement [23]. There is a deterministic change in the structure of neighborhoods defined by different operators every time a local minimum is reached. Initially, the algorithm starts by setting the number of neighborhoods, k (a value of k = 3 is commonly well accepted). For each neighborhood, it is necessary to choice the set of structures defined by the local operators. A neighborhood is a combination of some of the operators described above. Multi-operators are used in order to explore the solution space more extensively. Afterward, the VND algorithm repeats the following process until an improvement is attained. First, the loop begins with k = 1 and the next neighborhood Nk (1) is selected. Nk (k) represents the set of solutions in the neighborhood of s linked by neighborhood relationships (also is usually called search space). If s′ , the best solution in Nk (k), is better than s, then s ←− s′ , otherwise, k = k + 1 and start the loop again until k > kmax = 3. The VND algorithm is shown in Fig. 1. 4.2.2. General Variable Neighborhood Search with Tabu Search GVNS is an extended variant which uses more than one neighborhood in a local search [24] and TS is a metaheuristic that uses the memory of the search process to improve the solutions obtained. TS is a local search strategy where the best solution in the neighborhood of the current solution is selected as the new current solution, even if it leads to an increase in solution cost [25].
6
J.A. Sicilia et al. / Journal of Computational and Applied Mathematics (
)
–
Fig. 1. Variable Neighborhood Descent.
Fig. 2. General Variable Neighborhood Search with Tabu Search.
Due to the fact that the technique of changing the neighborhood structure systematically is simple and very powerful, another way to extend the VNS is incorporating it into other metaheuristics. This may give rise to new methods. The metaheuristic selected is the TS that uses a neighborhood structure to execute up and down movements working with different memory types. The new algorithm (GVNS method with TS) is described in Fig. 2. The initialization of this algorithm is the selection of the set of neighborhood structures defined by the intra-route and inter-route operators both for the TS phase and for the local search. The route construction method enables the algorithm to determine an initial solution s. Then the algorithm repeats the following process until a stopping condition is fulfilled. First of all, the next neighborhood Nk (s) is selected. Afterward, a new solution s′ is generated by using TS and replacing the shaking phase of the basic GVNS algorithm. This has been carried out to eliminate randomness and to diversify the exploration neighborhood, avoiding the decrease in local minimums and allowing the search for better solutions. Subsequently, VND is applied to obtain the new solution s′′ . This new solution s′′ is accepted as the next one if the objective is improved. If no stopping condition is met, the algorithm continues with the next iteration. 5. Computational experiments The computational experiments presented in this section have been extracted from real cases experienced by a large transport company in Spain. This company is in charge of the distribution of goods in urban areas where the distance among different customers is less than 100 km. All the necessary characteristics of the elements that constitute the problem have been extracted from the historical data. The first elements were the orders including the geographical position of customers for pickup, delivery or pickup–delivery dependent operations, the quantity and type of the load, the time required to complete the operation and the time windows in which each customer had to be served. The second element was the heterogeneous fleet with its capacity limit and the type of transportable load. The fulfillment of the time constraints imposed by each customer and the compatibilities between goods and vehicles maintaining a load balance between vehicles (maximum orders per vehicle) are the general constraints imposed for all test problems. In addition to this, it is mandatory to serve all orders in order to maintain quality of service. All vehicles must begin their daily route at the depot but they have no obligation to return to the depot. The drivers home may be the last node of the route being defined as one more order on the route with null quantity. The test problem generation is defined by the number of orders (O), the number of available vehicles (V) and the maximum number of orders per vehicle (L) (see Table 1). Both the service time and the time window of each order are imposed by the customer so they are independent and cannot be changed. The results of the experiments are listed in Tables 2 and 3, obtained by multiple runs of the algorithm. The computing time was less than four minutes in any case. The first four columns define the setting parameters of the experiment. The rest of the columns denotes the efficiency of all the solutions by means of two performance metrics that make reference to the objective function. The first one is the total traveled distance (∆(d)) and the second one is the number of vehicles needed to transport the goods (∆(v)). One of the main conclusions that arise from the results is to examine as the construction method obtains the best solutions for a small volume of orders. This is mainly due to the high number of constraints imposed on the problem initially. This is also confirmed by the previous test carried out with the GLPK tool (a package intended for solving large-scale linear
J.A. Sicilia et al. / Journal of Computational and Applied Mathematics (
)
–
7
Table 1 Parameters of the computational experiments. Problems P
Orders O
Vehicles V
Maximum orders per vehicle L
P1–P16 P17–P46 P47–P70 P71–P94
50 100 200 500
7–8–9–10 14–15–16–17–18 26–27–28–30 64–65–66–68
8–10–12–14 8–10–12–14–16–20 8–10–12–14–16–20 8–10–12–14–16–20
Table 2 Algorithm results of the test problems (P1–P46). Distance is measured in kilometers. Prs
O
V
L
Constructor
∆ (d ) P1 P2 P3 P4 P4 P6 P7 P8 P9 P10 P11 P12 P13 P14 P15 P16
50
P17 P18 P19 P20 P21 P22 P23 P24 P25 P26 P27 P28 P29 P30 P31 P32 P33 P34 P35 P36 P37 P38 P39 P40 P41 P42 P43 P44 P45 P46
100
7
8
9
10
14
15
16
18
20
VND
∆(v)
∆(d)
GVNS with TS
∆(v)
∆ (d )
∆(v)
8 10 12 14 8 10 12 14 8 10 12 14 8 10 12 14
691.84 700.66 670.72 670.72 776.34 785.15 755.22 755.22 841.05 871.32 835.80 835.80 892.43 918.00 918.00 918.00
7 7 7 7 8 8 8 8 9 9 9 9 10 10 10 10
690.04 698.85 668.92 668.92 776.26 785.08 755.14 755.14 840.97 871.25 835.73 835.73 892.36 917.92 917.92 917.92
7 7 7 7 8 8 8 8 9 9 9 9 10 10 10 10
691.84 700.66 670.72 670.72 776.26 785.08 794.24 756.98 840.97 876.47 847.15 835.73 892.36 944.61 945.56 918.00
7 7 7 7 8 8 6 6 9 8 7 9 10 8 8 10
8 10 12 14 16 20 8 10 12 14 16 20 8 10 12 14 16 20 8 10 12 14 16 20 8 10 12 14 16 20
1552.14 1497.62 1509.37 1533.78 1533.78 1533.78 1594.31 1598.67 1613.30 1613.30 1613.30 1613.30 1758.53 1698.94 1672.37 1672.37 1672.37 1672.37 1884.82 1871.18 1845.52 1845.52 1845.52 1845.52 2061.55 2063.54 2037.88 2037.88 2037.88 2037.88
14 14 14 14 14 14 15 15 15 15 15 15 16 16 16 16 16 16 16 16 16 16 16 16 18 18 18 18 18 18
1466.64 1414.18 1430.12 1455.45 1455.45 1455.45 1543.27 1522.30 1538.17 1538.17 1538.17 1538.17 1714.83 1616.12 1626.70 1626.70 1626.70 1626.70 1855.57 1802.65 1801.05 1801.05 1801.05 1801.05 2029.57 1994.90 1993.30 1993.30 1993.30 1993.30
14 14 14 14 14 14 15 15 15 15 15 15 16 16 16 16 16 16 16 16 16 16 16 16 18 18 18 18 18 18
1466.64 1395.64 1452.43 1470.26 1496.56 1460.98 1546.27 1493.21 1552.14 1532.45 1579.28 1548.48 1714.83 1649.32 1635.46 1645.23 1616.73 1633.65 1867.58 1814.47 1794.56 1800.05 1790.08 1807.86 2041.58 1993.72 1986.72 1992.81 2000.68 2024.37
14 12 12 13 13 13 15 12 11 14 14 14 16 13 15 15 15 15 16 14 15 15 14 14 18 15 17 17 17 16
programming (LP) and mixed integer programming (MIP)) by the simplex method that proved to be time-consuming and never-ending. In the experiments with a large number of orders, no solutions were found even after a time of 60 min. The results, shown in Tables 2 and 3, refer to the three methods that compose the algorithm. The first one is the construction method that obtains an initial feasible solution due to the quantity of constraints imposed, the second one is the VND metaheuristic and the last one is the GVNS method with TS.
8
J.A. Sicilia et al. / Journal of Computational and Applied Mathematics (
)
–
Table 3 Algorithm results of the test problems (P47–P94). Distance is measured in kilometers. Prs
O
P47 P48 P49 P50 P51 P52 P53 P54 P55 P56 P57 P58 P59 P60 P61 P62 P63 P64 P65 P66 P67 P68 P69 P70
200
P71 P72 P73 P74 P75 P76 P77 P78 P79 P80 P81 P82 P83 P84 P85 P86 P87 P88 P89 P90 P91 P92 P93 P94
500
V
26
27
28
30
64
65
66
68
L
Constructor
VND
GVNS with TS
∆ (d )
∆(v)
∆ (d )
∆(v)
∆(d)
∆(v)
8 10 12 14 16 20 8 10 12 14 16 20 8 10 12 14 16 20 8 10 12 14 16 20
3269.24 3156.34 2995.92 1990.20 2939.90 2937.52 3259.45 3152.76 3167.86 3163.83 3164.66 3164.66 3406.95 3256.71 3293.72 3324.40 3312.08 3312.08 3458.49 3354.54 3329.70 3424.26 3361.37 3347.08
26 26 26 26 26 26 27 27 27 27 27 27 28 28 28 28 28 28 30 30 30 30 30 30
2981.14 2868.32 1821.34 2818.59 2813.78 2812.26 3005.19 2900.68 2097.10 2910.76 2911.59 2911.59 3137.69 3006.36 3025.08 3024.63 3046.55 3046.55 3283.34 3163.60 3153.73 3169.83 3169.63 3153.96
26 26 26 26 26 26 27 27 27 27 27 27 28 28 28 28 28 28 30 30 30 30 30 30
2982.26 2865.17 2880.28 2833.39 2834.27 2852.44 3005.19 2937.82 2944.53 2975.75 2942.26 2935.87 3137.69 3043.65 3073.58 3047.86 3055.36 3107.04 3283.10 3138.64 3157.61 3168.94 3169.55 3147.87
26 24 24 25 25 23 27 26 26 22 20 21 28 27 27 25 27 20 30 26 29 30 30 25
8 10 12 14 16 20 8 10 12 14 16 20 8 10 12 14 16 20 8 10 12 14 16 20
4784.44 4618.51 4735.81 4786.16 4780.48 4790.88 4756.12 4726.18 4715.73 4830.47 4726.07 4706.99 4837.54 4683.28 4821.26 4986.74 4889.89 4866.25 4937.82 4726.20 4891.40 4839.39 4705.25 4746.83
64 61 61 61 61 61 65 62 62 62 62 62 66 63 63 63 63 63 68 64 65 65 65 65
3886.63 3915.05 3932.03 3875.15 3881.34 3891.85 3953.87 3972.34 3871.10 3918.63 3923.62 3924.17 3991.28 3964.06 3993.59 4005.54 3974.84 3977.15 4113.27 4046.15 4068.07 3981.25 3991.96 4016.21
64 61 61 61 61 61 65 62 62 62 62 62 66 63 63 63 63 63 68 64 65 65 65 65
3866.63 4618.51 3990.26 4786.16 3953.05 3908.25 4756.12 4726.18 3949.81 4830.47 4056.05 4706.99 3990.49 4683.28 4100.81 4986.74 4101.60 4001.04 4937.82 4042.66 4113.26 3993.13 4705.25 4746.83
64 61 57 61 60 60 65 62 59 62 59 62 66 63 59 63 59 61 68 64 63 64 65 65
6. Discussion of results As the number of available vehicles increases, the distance increases since the algorithm balances the freight to the extent possible. For instance, note that in the instance P13 the distance increases by 22% compared to the instance P14 for the three methods. On the other hand, the distance does not rise with respect to the increase in the number of available vehicles and even decreases for most instances with sizes of orders greater than 200 orders especially for metaheuristics methods. The GVNS method with TS in the instance P57 reduces the distance by 2% compared to the instance P47 and in the instance P85 increases the distance by 5% compared to the instance P71. In addition to this, the metaheuristic methods encounter better solutions for large volumes especially in terms of a reduction of used vehicles. This can be seen in some instances. For example, 11 of the 15 available vehicles are used in the instance P25, 21 of 27 are used in P57 and 57 of 64 are used in P73. Concerning the result comparison of the three methods, it is important to mention that the VND metaheuristic is performed to minimize the distance of the obtained solution by the construction method while the GVND method with
J.A. Sicilia et al. / Journal of Computational and Applied Mathematics (
)
–
9
TS prioritizes the reduction in the number of used vehicles although the distance is longer under the condition that the maximum number of orders per vehicle is not exceeded. Furthermore, VND method increases the reduction percentage of the total distance in respect to the construction method for a larger volume of orders. For 50 orders, the results of the construction method and the VND method are similar. However, for 100 orders the VND method reduces the distance by 2%–5%, for 200 reduces the distance by 6%–8% and for 500 orders reduces the distance by 15%–17%. The GVNS method with TS also improves the results of the construction method with regard to the distance as the volume of orders increases. For 50 orders results are similar, but for 100 orders the distance is reduced by 3%–6%, for 200 orders is reduced by 4%–7% and for 500 orders is reduced by 13%–16%. It can also be appreciated that initially the algorithm uses the greatest number of available vehicles, mainly for few orders, since the construction method only reduces the available vehicles for 500 orders or more. 7. Conclusions An optimization algorithm that solves the RVRP problem is presented in this paper. This problem takes into consideration constraints and complexities actually encountered in urban freight distribution which are not considered in the classical problems. Among them, time windows, capacity constraints, compatibility between orders and vehicles, maximum number of orders per vehicle, pickup–delivery dependent orders and not returning to the depot have been treated carefully throughout this paper with the aim of finding a general solution. Although the attraction to this problem has increasingly been growing in recent years, the relatively modest number of publications would seem to indicate that the domain of vehicle routing problems with real-life constraints and complexities is still young. Under this assumption, this work arose from the necessity to deal with it in a research project undertaken for a large Spanish distribution company operating in urban areas. The algorithm is currently integrated into a commercial software tool developed by ALERCE, a Spanish ICT enterprise, under the name of SENDA [26] (Advanced Dynamic Deliveries System) which is used in a daily manner and whose real data have helped to test our algorithm. The main contribution of the paper is twofold: the first one lies on the fact that the algorithm can further powered by adding more constraints other than those presented in the paper since the RVRP centralizes most of the variants of the classical problem; secondly, the algorithm has been properly wrapped in a computer tool that is commercialized by a ICT enterprise with presence in various countries around the world. A large part of the potential of the proposed methods are due to the definition of neighborhood structures which are based on reversing segments of routes and exchanging segments between routes. Furthermore, this algorithm is not only able to solve real-world instances with hundreds of orders efficiently but also computes the solutions within a reasonable time. The methods have been verified by means of several instances extracted from real data provided by a large transport company in Spain. Their performance has been optimized by reducing the total cost of the transport operations. This has been attained by minimizing the total traveled distance, reducing the number of used vehicles and managing the different resources optimally. Future research could aim to use this proposal to solve more general forms of vehicle routing problems which could contain more real-world objectives and constraints. The powerful capacity of this algorithm to find efficient solutions to difficult combinatorial optimization problems should make it useful for solving many other problems in transportation and logistics. Acknowledgments This activity has been partially supported by the Universidad Internacional de La Rioja (UNIR, http://www.unir.net), under the Plan Propio de Investigación, Desarrollo e Innovación [2013–2015]. Research group: Matemática aplicada al mundo real (MAMUR), in part by UNIR Research Support Strategy 2013–2015, under the CYBERSECURITICS.es Research Group [http://research.unir.net] and by the grant MTM2014-52016-C2-1-P. References [1] Á.A. Magreñán, A. Cordero, J.M. Gutiérrez, J.R. Torregrosa, Real qualitative behavior of a fourth-order family of iterative methods by using the convergence plane, Math. Comput. Simul. 105 (2014) 49–61. [2] A. Cordero, L. Feng, Á.A. Magreñán, J.R. Torregrosa, A new fourth-order family for solving nonlinear problems and its dynamics, J. Math. Chem. 53 (3) (2015) 893–910. [3] I.K. Argyros, A. Cordero, Á.A. Magreñán, J.R. Torregrosa, On the convergence of a Damped Secant method with modified right-hand side vector, Appl. Math. Comput. 252 (2015) 315–323. [4] Ministerio de Fomento, Observatorio de transporte de mercancías por carretera, Gobierno de España, Madrid, España, 2011. [5] European Commission, Green Paper: Towards A New Culture for Urban Mobility, European Union, Brussels, Belgium, 2007. [6] P. Ji, K. Chen, The vehicle routing problem: the case of the Hong Kong Postal, Transp. Plan. Technol. 30 (2) (2007) 167–182. [7] R. Lahyani, M. Khemakhem, F. Semet, Rich vehicle routing problems: From a taxonomy to a definition, European J. Oper. Res. 241 (1) (2015) 1–14. [8] J. Oppen, A. Lokketangen, J. Desrosiers, Solving a rich vehicle routing and inventory problem using column generation, Comput. Oper. Res. 37 (7) (2010) 1308–1317. [9] T.G. Crainic, M. Gendreau, W. Rei, A column generation approach for a multi-attribute vehicle routing problem, European J. Oper. Res. 241 (3) (2015) 886–906.
10
J.A. Sicilia et al. / Journal of Computational and Applied Mathematics (
)
–
[10] T. Vidal, T.G. Crainic, M. Gendreau, C. Prins, Heuristics for multi-attribute vehicle routing problems: A survey and synthesis, European J. Oper. Res. 231 (1) (2013) 1–21. [11] M. Gendreau, G. Laporte, J.Y. Potvin, Metaheuristics for the capacitated VRP, in: P. Toth, D. Vigo (Eds.), The Vehicle Routing Problem, SIAM Publishing, Philadelphia, USA, 2002, pp. 129–154. [12] N.A. El-sherbeny, Vehicle routing with time windows: An overview of exact, heuristic and metaheuristic methods, J. King Saud Univ. - Sci. 22 (3) (2010) 123–131. [13] I.M. Chao, T.S. Liou, A new Tabu Search heuristic for the site-dependent vehicle routing problem, in: B. Golden, S. Raghavan, E. Wasil (Eds.), The Next Wave in Computing, Optimization and Decision Technologies, Springer, USA, 2005, pp. 107–119. [14] F.P. Goksal, I. Karaoglan, F. Altiparmak, A hybrid discrete particle swarm optimization for vehicle routing problem with simultaneous pickup and delivery, Comput. Ind. Eng. 65 (2013) 39–53. [15] E.E. Zachariadis, C.T. Kiranoudis, An effective local search approach for the vehicle routing problem with backhauls, Expert Syst. Appl. 39 (3) (2012) 3174–3184. [16] P.P. Repoussis, C.D. Tarantilis, O. Braysy, G. Ioannou, A hybrid evolution strategy for the open vehicle routing problem, Comput. Oper. Res. 37 (3) (2010) 443–455. [17] N. Jozefowiez, F. Semet, E.G. Talbi, An evolutionary algorithm for the vehicle routing problem with route balancing, European J. Oper. Res. 195 (3) (2009) 761–769. [18] I. Or, Travelling salesman-type combinatorial problems and their relation to the logistics of blood banking (Ph.D. Thesis), Department of Industrial Engineering and Management Science, Northwestern University, Evanston, USA, 1976. [19] 0 Bräysy, G. Hasle, W. Dullaert, A multi-start local search algorithm for the vehicle routing problem with time windows, European J. Oper. Res. 159 (2004) 586–605. [20] 0 Bräysy, W. Dullaert, A fast evolutionary metaheuristic for the vehicle routing problem with time windows, Int. J. Artif. Intell. Tools 12 (2) (2003) 153–172. [21] J.Y. Potvin, J.M. Rousseau, An exchange heuristic for routing problems with time windows, J. Oper. Res. Soc. 46 (12) (1995) 1433–1446. [22] P. Hansen, N. Mladenovic, Variable neighborhood search: principles and applications, European J. Oper. Res. 130 (3) (2001) 449–467. [23] M. Tagmouti, M. Gendreau, J.Y. Potvin, A variable neighborhood descent heuristic for arc routing problems with time-dependent service costs, Comput. Ind. Eng. 59 (4) (2010) 954–963. [24] N. Mladenovic, M. Drazic, V. Kovacevic-Vujcic, M. Cangalovic, General variable neighborhood search for the continuous optimization, European J. Oper. Res. 191 (3) (2008) 753–770. [25] E.E. Zachariadis, C.D. Tarantilis, C.T. Kiranoudis, A guided tabu search for the vehicle routing problem with two-dimensional loading constraints, European J. Oper. Res. 195 (3) (2009) 729–743. [26] Visit http://www.alerce.es/es/senda.html for further information about the commercial tool (in Spanish).
)
–
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
An optimization algorithm for solving the rich vehicle routing problem based on Variable Neighborhood Search and Tabu Search metaheuristics Juan Antonio Sicilia a,∗ , Carlos Quemada b , Beatriz Royo c , David Escuín d a
Escuela de Ingeniería, Universidad Internacional de La Rioja, 26002 Logroño, Spain
b
Virginia Polytechnic Institute and State University in Blacksurg, Virginia, USA
c
Department of Mechanical Engineering, University of Zaragoza, 50018 Zaragoza, Spain
d
Instituto Tecnológico de Aragón, María Luna 1, 50018 Zaragoza, Spain
article
info
Article history: Received 13 October 2014 Received in revised form 9 March 2015 Keywords: Rich vehicle routing problem Optimization algorithm Metaheuristic methods Load balance
abstract This paper presents a novel optimization algorithm that consists of metaheuristic processes to solve the problem of the capillary distribution of goods in major urban areas taking into consideration the features encountered in real life: time windows, capacity constraints, compatibility between orders and vehicles, maximum number of orders per vehicle, orders that depend on the pickup and delivery and not returning to the depot. With the intention of reducing the wide variety of constraints and complexities, known as the Rich Vehicle Routing Problem, this algorithm proposes feasible alternatives in order to achieve the main objective of this research work: the reduction of costs by minimizing distances and reducing the number of vehicles used as long as the service quality to customers is optimum and a load balance among vehicles is maintained. © 2015 Elsevier B.V. All rights reserved.
1. Introduction A vast number of problems from Applied Science including engineering can be brought by means of solving a nonlinear equation using mathematical modeling [1–3]. One of that problems, in concrete the capillary transport of goods problem, is studied in this paper. The logistics of the capillary transport of goods has an essential urban dimension. On account of this, distribution requires efficient systems and the process between warehouses and customers must be efficient and clean. Therefore, it is necessary to optimize urban logistics efficiently and improve the connections between urban and interurban freight transport in order to ensure efficient distribution [4]. In the European Union, over 60% of the population lives in urban areas, generating just under 85% of the EU’s Gross Domestic Product (GDP) [5]. In addition to this, pickup and delivery actions in urban areas represent about 40% of the total cost of transport activities carried out at home. These expenses are further increased by the reduction of stocks and size of goods and by the increase in the number of requests.
∗
Corresponding author. E-mail addresses: [email protected] (J.A. Sicilia), [email protected] (C. Quemada), [email protected] (B. Royo), [email protected] (D. Escuín). http://dx.doi.org/10.1016/j.cam.2015.03.050 0377-0427/© 2015 Elsevier B.V. All rights reserved.
2
J.A. Sicilia et al. / Journal of Computational and Applied Mathematics (
)
–
Due to this, there is great potential for freight transport in these areas. Therefore, it is necessary to find logistics solutions for urban freight distribution that are consistent with the restrictions imposed by the law in order to protect the interests of the public. This paper presents an optimization algorithm that solves a Rich Vehicle Routing Problem (RVRP) and arises from a research project carried out for an important Spanish distribution company. The main goal of this project is to manage its resources in urban areas by reducing costs caused by inefficiency and ineffectiveness as much as possible. Therefore, this research work consists in reducing the constraints and complexities currently encountered in the urban distribution of goods by using the Variable Neighborhood Descent (VND), the General Variable Neighborhood Search (GVNS) and Tabu Search (TS) metaheuristic methods. Most problems in industry, particularly in city logistics, are multi-objective in nature. This makes it difficult to apply classical routing models to real-life problems because these models are developed with the single objective of minimizing the cost of the solution [6]. However, in real life there may be different requirements related to a single tour such as balancing workloads, time, distance, cost, etc. That need to be taken into consideration simply by adding new constraints. This paper is organized as follows. In Section 2, the literature review is presented. Subsequently, Section 3 presents a formal mathematical definition of the RVRP. Section 4 proposes an optimization algorithm to solve this problem using metaheuristic methods. Section 5 shows the computational calculation results of this algorithm and these results are discussed in Section 6. Finally, conclusions are set out in Section 7. 2. Literature review Over the last years, several variants of multiconstrained Vehicle Routing Problems have been studied, forming a class of problems known as Rich Vehicle Routing Problems. A comprehensive and relevant taxonomy for the RVRP literature and an attempt to provide a definition of RVRPs can be found in [7]. Several applications of such problems are close to real world as the work shown by [8] to deal with animal welfare and global constraints to handle production and inventory. Other authors have proposed to refer RVRPs as Multi-Attribute Vehicle Routing Problems (MAVRP). The work presented by [9] considers a multi-attribute vehicle routing problem derived from a real-life milk collection system. This problem is characterized by the presence of a heterogeneous fleet of vehicles, multiple depots, and several resource constraints. [10] claims that the attributes of vehicle routing problems are additional characteristics or constraints that aim to better take into account the specifications of real applications. It then takes a closer look at the concepts of 64 remarkable meta-heuristics, selected objectively for their outstanding performance on 15 classic MAVRP with different attributes. This cross-analysis leads to the identification of winning strategies in designing effective heuristics for MAVRP. This is an important step in the development of general and efficient solution methods for dealing with the large range of vehicle routing variants. The problem RVRP considered in this paper is based on the optimization of the freight distribution in large urban areas taking into consideration the real characteristics presented in companies and trades. Given a fleet of vehicles and a set of customers dispersed over a geographic area, the system must be able to attain a feasible solution that leads to a decrease of the total cost of the problem. The purpose is to reach a solution with the minimum number of vehicles and minimizing the total distance traveled. The model of the problem depends on the following requirements:
• CVRP (Capacitated Vehicle Routing Problem): Vehicles have a maximum limit of capacity and are homogeneous [11]. • VRPTW (Vehicle Routing Problem with Time Windows): Each customer must be serviced in a set time interval [12]. • SDVRP (Site-Dependent Vehicle Routing Problem): Certain orders can only be transported by specific vehicles. Therefore, the fleet must be heterogeneous [13].
• VRPPD (Vehicle Routing Problem with Pickups and Deliveries): Customers may receive and send goods [14]. • VRPB (Vehicle Routing Problem with Backhauls): The goods that leave the depot must be delivered to customers and the goods that are picked up from customers must be transported to the depot [15].
• OVRP (Open Vehicle Routing Problem): Vehicles are not required to return to the depot after visiting the last customer and the route may finish near to the drivers home [16].
• Load balance among routes: It is important to balance the workload among vehicles since the competitiveness of a transport company depends not only on minimizing costs and distances to obtain higher profit margins but also on its ability to treat employees fairly [17]. 3. Mathematical formulation Prior to describing the mathematical formulation, it is important to define what an order is in the context of this paper. An order is a request made by a customer for the transport of goods. There are three different kinds of orders: one-pickup orders, one-delivery orders or pickup–delivery orders at dependent nodes. These nodes have to be passed through in a specific and proper sequence by the same vehicle with the time window imposed by the customer. Each vehicle has a maximum of a day to complete the service. Since a customer might request more than one order, each customer can be visited more than once by different vehicles. In addition to this, there are dependent orders where vehicles have to pick up goods at a specific node and deliver them to another node without the need to return to the warehouse.
J.A. Sicilia et al. / Journal of Computational and Applied Mathematics (
)
–
3
Furthermore, vehicle capacity constraints and compatibilities between orders and vehicles might lead certain vehicles to be excluded from allocating orders or certain orders to be transported only by specific vehicles. This situation requires a heterogeneous fleet in capacity and technical properties. The vehicles may have different capacities and types of loads (normal, refrigerated, isotherm, etc.). Two different location points are assigned to each vehicle. The end node may be the warehouse or the driver’s home and the initial node corresponds to the warehouse regardless of where the vehicle have finished the previous day. Every order has a service time (load/unload time of goods) and a time window in which the operation must be carried out. It is important to point out that this time window does not prevent any vehicle from arriving before the lower time limit and this might cause losses due to the inactivity of a vehicle during a period of time. The intention is for the arrival of each vehicle to be within that time interval. Another important factor that should be taken into account is the load balance among drivers. A maximum number of orders per vehicle is established so a route is carried out independently of its costs. In addition to this, another balance between the distance and the driver actions must be attained in all routes. Since the number of orders might vary each day, the number of customers (geographic nodes) is not fixed and the distances have to be exact in urban areas. A Geographic Information System (GIS) is used to calculate the precise distance from the GPS position of each node to be visited. This distance is only calculated among orders that may go on the same route. Previously, different sets of orders that may be compatible with each other are analyzed. The aim of this work is to create optimized routes that minimize the total cost of the transport operations by minimizing the total distance, reducing the number of vehicles and eliminating inactivity times. Furthermore, a load balance among vehicles must be established, the available resources must be maximized and there must be a commitment to the service quality. One route is defined as an ordered sequence of orders that the same vehicle has to dispatch. Regardless of cost, time is the most important factor to take into account. Therefore, it is necessary to fulfill the pickup and delivery schedule. The RVRP can be described as follows. Let G = (N , A) be a complete graph with N = {1}∪ O ∪{o + 1} the set of nodes and O = {2, . . . , o} the set of orders to be served. The nodes {1} and {o + 1} represent the depot and are added to the network to allow empty tours. The aim is to design a set of routes for a heterogeneous fleet of V vehicles (each has a capacity, a load type and a limited number of orders to be served), servicing a set of O orders. Each arc (i, j) ∈ A has a cost cij . The interval [ai , bi ] denotes the time windows for order i ∈ O and ti represents the time service. A vehicle cannot start servicing node i before ai and after bi . A vehicle can arrive before ai and wait for service. Drivers homes are represented by some nodes and vehicles have to serve them just before returning to the depot. As explained above, it is possible to distinguish three kinds of nodes:
• Deliveries: vehicles depart from depot with goods. • Pickups: vehicles pick up goods at these nodes and transport them to depot. • Pickup-and-delivery dependents: vehicles have to pick up goods at pickup node and deliver them at delivery node. On each route the total quantity of goods picked up and delivered cannot exceed vehicle capacity. Each node i represents both a delivery node, di , and a pickup node, pi . The pickup node of a delivery node i is represented by ri . In addition to this, a node may have a demand to be picked up and delivered. Different nodes i ∈ N may correspond to the same geographical location. The set of available vehicles is represented by V , being Qk the capacity of vehicle k ∈ V , Tk the maximum number of nodes that vehicle k can serve. hk represents the last node to be served (driver’s home) before returning to depot, hk ∈ O. Furthermore, mik denotes the compatibility between node i and vehicle k, ∀i ∈ N, ∀k ∈ V . Node i may be served by vehicle k if and only if mik = 1, mik = 0 otherwise. The RVRP is formulated as follows: Decision variables: xkij =
if arc (i, j) ∈ A is traveled by vehicle k ∈ V otherwise.
1 0
(1)
Objective function: min
cij xkij −
k∈V i∈N j∈N
xk1hk M .
(2)
k∈V
The first term represents the accumulated costs for orders served by vehicles; the second term represents the number of utilized vehicles. Subject to:
xkij = 1 ∀i ∈ O
(3)
j∈N k∈V
i∈N
xkih −
j∈N
xkhj = 0 ∀h ∈ O, ∀k ∈ V
(4)
4
J.A. Sicilia et al. / Journal of Computational and Applied Mathematics (
)
–
xkii = 0 ∀i ∈ V , ∀h ∈ N
(5)
xko+1i = 0 ∀k ∈ V , ∀i ∈ N
(6)
xkhk o+1 = 1 ∀k ∈ V
(7)
(8)
xk1j = 1 ∀k ∈ V
j∈N
ski + ti − skj + cij ≤ (1 − xkij )M
∀k ∈ V , ∀i ∈ N , ∀j ∈ N
sk1 = 0 ∀k ∈ V ai ≤
ski
(9) (10)
= 0 ∀k ∈ V , ∀i ∈ N
(11)
∀k ∈ V , ∀i ∈ N
(12)
bi ≥ ski = 0
mik + mjk ≥ 2xkij
∀k ∈ V , ∀i ∈ N
lki ≤ Qk
∀k ∈ V , ∀i ∈ O, ∀j ∈ O
di xkij = lk1
(13) (14)
∀k ∈ V
(15)
i∈N j∈N
pi xkij = lko+1
∀k ∈ V
(16)
i∈N j∈N
lki − dj + pj − lkj ≤ (1 − xkij )M Tk + 1 ≥
xkij
∀k ∈ V , ∀i ∈ N , ∀j ∈ N
∀k ∈ V
(17) (18)
i∈N j∈N
1−
xkrj h ≤ (1 − xkij )M
∀k ∈ V , ∀i ∈ N , ∀j ∈ O
(19)
h∈N
skri ≤ ski
∀k ∈ V , ∀i ∈ O.
(20)
The domain of decision variables is set by (1). The objective function (2) states that the distance traveled and number of vehicles should be minimized, M being a sufficiently large constant. (3) ensures that one and only one node can be visited after an order. Flow conservation is ensured by constraints (4) and (5) avoids cycles. Constraint (8) states that all vehicles leave the depot while constraint (7) denotes all vehicles returning to depot after servicing the home drivers. Constraint (6) ensures no vehicles leave the virtual depot and vehicles end their routes at node o + 1. Demand constraint is ensured by (14)–(17) with lki to denote the load of the vehicle k at node i after serving node i. Compatibility between vehicles and orders is represented by (13). Time window constraints are ensured by constraints (9)–(12) with ski to denote when the vehicle k starts the service at node i. Constraints (11) and (12) ensure that node i cannot be served before ai and after bi . Constraint (18) ensures that the number of customers served by vehicle k is not greater than the maximum permitted. Constraints (19) and (20) represent precedence constraints for pickup-and-delivery dependents. 4. Optimization algorithm The main difficulty of the RVRP problem is that if the number of elements to be combined is relatively large, the combinatorial result grows exponentially turning it into an NP-hard problem. On account of this, it is necessary to apply optimization techniques in order to obtain a high-quality solution in reasonable computational times. In this section an efficient algorithm to solve the problem is presented. The RVRP is solved by means of a powerful route construction method that starts with the determination of a feasible initial solution and the metaheuristics VND, GVNS and TS based on various elementary neighborhood operators to obtain an improved solution. 4.1. Route construction method At the beginning of this process there are several distribution warehouses, a heterogeneous fleet of vehicles and a set of orders to be distributed. The algorithm will generate as many routes as vehicles are available, taking into consideration the constraints and compatibilities.
J.A. Sicilia et al. / Journal of Computational and Applied Mathematics (
)
–
5
The route construction method consists of assigning orders to vehicles depending on certain parameters and taking into consideration the limits of load balance between routes. This procedure has been implemented randomly to avoid obtaining of the same solutions and to diversify the total space of the search. In the beginning there is no route and the system proceeds to create routes by choosing the first order according to six available criteria. This does not mean that this order is the first one of the final route. It has been selected only to start the insertion process. In this insertion process, orders are inserted until the route is not feasible or until it has reached the maximum number of orders allowed per route. The process finishes when there are no orders dispatched in the routes. This is the initial solution that goes to the local improvements phase. To select the first order i of an empty route, the following insertion criteria might be considered:
• The node i with the biggest value of c1i + ti + cio+1 . That is, the node with the latest arrival time to the depot without • • • •
considering the home driver yet. The node i with the smallest value of bi : earliest close of the time window of the node. The node i with the shortest value of bi − ai : the tightest time window. The node i with the smallest value of cio+1 : the effective distance from the customer to depot. Best average of the above criteria.
The last one is the most common criterion because it includes the rest of criteria. In this case, all criteria are calculated in order to establish a list with orders that obtains the best value in terms of the rest of criteria. The insertion process of the rest of orders is based on the following criteria:
• The largest time margin: it is the extra time available to distribute other orders taking into consideration the time taken to distribute orders that have already been inserted.
• The shortest distance of the route: choice of a nearby customer. The order will be inserted in the position of the route with minimum distance. 4.2. Local improvements and metaheuristics Local improvement operators are methods applied to a set of one or more routes in order to improve the evaluation function. The basic process is made up of movements and exchanges of order segments. It is necessary to apply the operators to improve specific aspects of the problem in the majority of the calculations. These operators are the well-known intraroute improvements applied to a single route and that are computationally very fast, OR-OPT [18] and IOPT [19], and the inter-route improvements that exchange segments of orders between two routes, ICROSS [20] and 2-OPT* [21]. The OROPT operator moves one segment of the route into other position. IOPT means that the sequence is reversed. The CROSS exchanges two (one being potentially empty) customer sequences s1 and s2 from two routes. The customer sequences s1 and s2 can eventually be reversed in this process to produce a slightly larger neighborhood, called I-CROSS. 2-OPT refers to the fact that there is an exchange or swap of a pair of edges. Once the initial solution, s, has been obtained, this is improved through the application of a set of metaheuristics. A sequence of operators with their respective configurations must be specified in order to define a neighborhood search. The metaheuristic process consists of using a neighborhood that exchanges segments. The final state of each operator of a neighborhood determines the subsequent execution of the following neighborhood. Therefore, a good configuration of operators enables the obtaining of one solution or another [22]. 4.2.1. Variable Neighborhood Descent VND is a variant of Variable Neighborhood Search that consists of replacing the current solution iteratively with the result from the local solution if there is an improvement [23]. There is a deterministic change in the structure of neighborhoods defined by different operators every time a local minimum is reached. Initially, the algorithm starts by setting the number of neighborhoods, k (a value of k = 3 is commonly well accepted). For each neighborhood, it is necessary to choice the set of structures defined by the local operators. A neighborhood is a combination of some of the operators described above. Multi-operators are used in order to explore the solution space more extensively. Afterward, the VND algorithm repeats the following process until an improvement is attained. First, the loop begins with k = 1 and the next neighborhood Nk (1) is selected. Nk (k) represents the set of solutions in the neighborhood of s linked by neighborhood relationships (also is usually called search space). If s′ , the best solution in Nk (k), is better than s, then s ←− s′ , otherwise, k = k + 1 and start the loop again until k > kmax = 3. The VND algorithm is shown in Fig. 1. 4.2.2. General Variable Neighborhood Search with Tabu Search GVNS is an extended variant which uses more than one neighborhood in a local search [24] and TS is a metaheuristic that uses the memory of the search process to improve the solutions obtained. TS is a local search strategy where the best solution in the neighborhood of the current solution is selected as the new current solution, even if it leads to an increase in solution cost [25].
6
J.A. Sicilia et al. / Journal of Computational and Applied Mathematics (
)
–
Fig. 1. Variable Neighborhood Descent.
Fig. 2. General Variable Neighborhood Search with Tabu Search.
Due to the fact that the technique of changing the neighborhood structure systematically is simple and very powerful, another way to extend the VNS is incorporating it into other metaheuristics. This may give rise to new methods. The metaheuristic selected is the TS that uses a neighborhood structure to execute up and down movements working with different memory types. The new algorithm (GVNS method with TS) is described in Fig. 2. The initialization of this algorithm is the selection of the set of neighborhood structures defined by the intra-route and inter-route operators both for the TS phase and for the local search. The route construction method enables the algorithm to determine an initial solution s. Then the algorithm repeats the following process until a stopping condition is fulfilled. First of all, the next neighborhood Nk (s) is selected. Afterward, a new solution s′ is generated by using TS and replacing the shaking phase of the basic GVNS algorithm. This has been carried out to eliminate randomness and to diversify the exploration neighborhood, avoiding the decrease in local minimums and allowing the search for better solutions. Subsequently, VND is applied to obtain the new solution s′′ . This new solution s′′ is accepted as the next one if the objective is improved. If no stopping condition is met, the algorithm continues with the next iteration. 5. Computational experiments The computational experiments presented in this section have been extracted from real cases experienced by a large transport company in Spain. This company is in charge of the distribution of goods in urban areas where the distance among different customers is less than 100 km. All the necessary characteristics of the elements that constitute the problem have been extracted from the historical data. The first elements were the orders including the geographical position of customers for pickup, delivery or pickup–delivery dependent operations, the quantity and type of the load, the time required to complete the operation and the time windows in which each customer had to be served. The second element was the heterogeneous fleet with its capacity limit and the type of transportable load. The fulfillment of the time constraints imposed by each customer and the compatibilities between goods and vehicles maintaining a load balance between vehicles (maximum orders per vehicle) are the general constraints imposed for all test problems. In addition to this, it is mandatory to serve all orders in order to maintain quality of service. All vehicles must begin their daily route at the depot but they have no obligation to return to the depot. The drivers home may be the last node of the route being defined as one more order on the route with null quantity. The test problem generation is defined by the number of orders (O), the number of available vehicles (V) and the maximum number of orders per vehicle (L) (see Table 1). Both the service time and the time window of each order are imposed by the customer so they are independent and cannot be changed. The results of the experiments are listed in Tables 2 and 3, obtained by multiple runs of the algorithm. The computing time was less than four minutes in any case. The first four columns define the setting parameters of the experiment. The rest of the columns denotes the efficiency of all the solutions by means of two performance metrics that make reference to the objective function. The first one is the total traveled distance (∆(d)) and the second one is the number of vehicles needed to transport the goods (∆(v)). One of the main conclusions that arise from the results is to examine as the construction method obtains the best solutions for a small volume of orders. This is mainly due to the high number of constraints imposed on the problem initially. This is also confirmed by the previous test carried out with the GLPK tool (a package intended for solving large-scale linear
J.A. Sicilia et al. / Journal of Computational and Applied Mathematics (
)
–
7
Table 1 Parameters of the computational experiments. Problems P
Orders O
Vehicles V
Maximum orders per vehicle L
P1–P16 P17–P46 P47–P70 P71–P94
50 100 200 500
7–8–9–10 14–15–16–17–18 26–27–28–30 64–65–66–68
8–10–12–14 8–10–12–14–16–20 8–10–12–14–16–20 8–10–12–14–16–20
Table 2 Algorithm results of the test problems (P1–P46). Distance is measured in kilometers. Prs
O
V
L
Constructor
∆ (d ) P1 P2 P3 P4 P4 P6 P7 P8 P9 P10 P11 P12 P13 P14 P15 P16
50
P17 P18 P19 P20 P21 P22 P23 P24 P25 P26 P27 P28 P29 P30 P31 P32 P33 P34 P35 P36 P37 P38 P39 P40 P41 P42 P43 P44 P45 P46
100
7
8
9
10
14
15
16
18
20
VND
∆(v)
∆(d)
GVNS with TS
∆(v)
∆ (d )
∆(v)
8 10 12 14 8 10 12 14 8 10 12 14 8 10 12 14
691.84 700.66 670.72 670.72 776.34 785.15 755.22 755.22 841.05 871.32 835.80 835.80 892.43 918.00 918.00 918.00
7 7 7 7 8 8 8 8 9 9 9 9 10 10 10 10
690.04 698.85 668.92 668.92 776.26 785.08 755.14 755.14 840.97 871.25 835.73 835.73 892.36 917.92 917.92 917.92
7 7 7 7 8 8 8 8 9 9 9 9 10 10 10 10
691.84 700.66 670.72 670.72 776.26 785.08 794.24 756.98 840.97 876.47 847.15 835.73 892.36 944.61 945.56 918.00
7 7 7 7 8 8 6 6 9 8 7 9 10 8 8 10
8 10 12 14 16 20 8 10 12 14 16 20 8 10 12 14 16 20 8 10 12 14 16 20 8 10 12 14 16 20
1552.14 1497.62 1509.37 1533.78 1533.78 1533.78 1594.31 1598.67 1613.30 1613.30 1613.30 1613.30 1758.53 1698.94 1672.37 1672.37 1672.37 1672.37 1884.82 1871.18 1845.52 1845.52 1845.52 1845.52 2061.55 2063.54 2037.88 2037.88 2037.88 2037.88
14 14 14 14 14 14 15 15 15 15 15 15 16 16 16 16 16 16 16 16 16 16 16 16 18 18 18 18 18 18
1466.64 1414.18 1430.12 1455.45 1455.45 1455.45 1543.27 1522.30 1538.17 1538.17 1538.17 1538.17 1714.83 1616.12 1626.70 1626.70 1626.70 1626.70 1855.57 1802.65 1801.05 1801.05 1801.05 1801.05 2029.57 1994.90 1993.30 1993.30 1993.30 1993.30
14 14 14 14 14 14 15 15 15 15 15 15 16 16 16 16 16 16 16 16 16 16 16 16 18 18 18 18 18 18
1466.64 1395.64 1452.43 1470.26 1496.56 1460.98 1546.27 1493.21 1552.14 1532.45 1579.28 1548.48 1714.83 1649.32 1635.46 1645.23 1616.73 1633.65 1867.58 1814.47 1794.56 1800.05 1790.08 1807.86 2041.58 1993.72 1986.72 1992.81 2000.68 2024.37
14 12 12 13 13 13 15 12 11 14 14 14 16 13 15 15 15 15 16 14 15 15 14 14 18 15 17 17 17 16
programming (LP) and mixed integer programming (MIP)) by the simplex method that proved to be time-consuming and never-ending. In the experiments with a large number of orders, no solutions were found even after a time of 60 min. The results, shown in Tables 2 and 3, refer to the three methods that compose the algorithm. The first one is the construction method that obtains an initial feasible solution due to the quantity of constraints imposed, the second one is the VND metaheuristic and the last one is the GVNS method with TS.
8
J.A. Sicilia et al. / Journal of Computational and Applied Mathematics (
)
–
Table 3 Algorithm results of the test problems (P47–P94). Distance is measured in kilometers. Prs
O
P47 P48 P49 P50 P51 P52 P53 P54 P55 P56 P57 P58 P59 P60 P61 P62 P63 P64 P65 P66 P67 P68 P69 P70
200
P71 P72 P73 P74 P75 P76 P77 P78 P79 P80 P81 P82 P83 P84 P85 P86 P87 P88 P89 P90 P91 P92 P93 P94
500
V
26
27
28
30
64
65
66
68
L
Constructor
VND
GVNS with TS
∆ (d )
∆(v)
∆ (d )
∆(v)
∆(d)
∆(v)
8 10 12 14 16 20 8 10 12 14 16 20 8 10 12 14 16 20 8 10 12 14 16 20
3269.24 3156.34 2995.92 1990.20 2939.90 2937.52 3259.45 3152.76 3167.86 3163.83 3164.66 3164.66 3406.95 3256.71 3293.72 3324.40 3312.08 3312.08 3458.49 3354.54 3329.70 3424.26 3361.37 3347.08
26 26 26 26 26 26 27 27 27 27 27 27 28 28 28 28 28 28 30 30 30 30 30 30
2981.14 2868.32 1821.34 2818.59 2813.78 2812.26 3005.19 2900.68 2097.10 2910.76 2911.59 2911.59 3137.69 3006.36 3025.08 3024.63 3046.55 3046.55 3283.34 3163.60 3153.73 3169.83 3169.63 3153.96
26 26 26 26 26 26 27 27 27 27 27 27 28 28 28 28 28 28 30 30 30 30 30 30
2982.26 2865.17 2880.28 2833.39 2834.27 2852.44 3005.19 2937.82 2944.53 2975.75 2942.26 2935.87 3137.69 3043.65 3073.58 3047.86 3055.36 3107.04 3283.10 3138.64 3157.61 3168.94 3169.55 3147.87
26 24 24 25 25 23 27 26 26 22 20 21 28 27 27 25 27 20 30 26 29 30 30 25
8 10 12 14 16 20 8 10 12 14 16 20 8 10 12 14 16 20 8 10 12 14 16 20
4784.44 4618.51 4735.81 4786.16 4780.48 4790.88 4756.12 4726.18 4715.73 4830.47 4726.07 4706.99 4837.54 4683.28 4821.26 4986.74 4889.89 4866.25 4937.82 4726.20 4891.40 4839.39 4705.25 4746.83
64 61 61 61 61 61 65 62 62 62 62 62 66 63 63 63 63 63 68 64 65 65 65 65
3886.63 3915.05 3932.03 3875.15 3881.34 3891.85 3953.87 3972.34 3871.10 3918.63 3923.62 3924.17 3991.28 3964.06 3993.59 4005.54 3974.84 3977.15 4113.27 4046.15 4068.07 3981.25 3991.96 4016.21
64 61 61 61 61 61 65 62 62 62 62 62 66 63 63 63 63 63 68 64 65 65 65 65
3866.63 4618.51 3990.26 4786.16 3953.05 3908.25 4756.12 4726.18 3949.81 4830.47 4056.05 4706.99 3990.49 4683.28 4100.81 4986.74 4101.60 4001.04 4937.82 4042.66 4113.26 3993.13 4705.25 4746.83
64 61 57 61 60 60 65 62 59 62 59 62 66 63 59 63 59 61 68 64 63 64 65 65
6. Discussion of results As the number of available vehicles increases, the distance increases since the algorithm balances the freight to the extent possible. For instance, note that in the instance P13 the distance increases by 22% compared to the instance P14 for the three methods. On the other hand, the distance does not rise with respect to the increase in the number of available vehicles and even decreases for most instances with sizes of orders greater than 200 orders especially for metaheuristics methods. The GVNS method with TS in the instance P57 reduces the distance by 2% compared to the instance P47 and in the instance P85 increases the distance by 5% compared to the instance P71. In addition to this, the metaheuristic methods encounter better solutions for large volumes especially in terms of a reduction of used vehicles. This can be seen in some instances. For example, 11 of the 15 available vehicles are used in the instance P25, 21 of 27 are used in P57 and 57 of 64 are used in P73. Concerning the result comparison of the three methods, it is important to mention that the VND metaheuristic is performed to minimize the distance of the obtained solution by the construction method while the GVND method with
J.A. Sicilia et al. / Journal of Computational and Applied Mathematics (
)
–
9
TS prioritizes the reduction in the number of used vehicles although the distance is longer under the condition that the maximum number of orders per vehicle is not exceeded. Furthermore, VND method increases the reduction percentage of the total distance in respect to the construction method for a larger volume of orders. For 50 orders, the results of the construction method and the VND method are similar. However, for 100 orders the VND method reduces the distance by 2%–5%, for 200 reduces the distance by 6%–8% and for 500 orders reduces the distance by 15%–17%. The GVNS method with TS also improves the results of the construction method with regard to the distance as the volume of orders increases. For 50 orders results are similar, but for 100 orders the distance is reduced by 3%–6%, for 200 orders is reduced by 4%–7% and for 500 orders is reduced by 13%–16%. It can also be appreciated that initially the algorithm uses the greatest number of available vehicles, mainly for few orders, since the construction method only reduces the available vehicles for 500 orders or more. 7. Conclusions An optimization algorithm that solves the RVRP problem is presented in this paper. This problem takes into consideration constraints and complexities actually encountered in urban freight distribution which are not considered in the classical problems. Among them, time windows, capacity constraints, compatibility between orders and vehicles, maximum number of orders per vehicle, pickup–delivery dependent orders and not returning to the depot have been treated carefully throughout this paper with the aim of finding a general solution. Although the attraction to this problem has increasingly been growing in recent years, the relatively modest number of publications would seem to indicate that the domain of vehicle routing problems with real-life constraints and complexities is still young. Under this assumption, this work arose from the necessity to deal with it in a research project undertaken for a large Spanish distribution company operating in urban areas. The algorithm is currently integrated into a commercial software tool developed by ALERCE, a Spanish ICT enterprise, under the name of SENDA [26] (Advanced Dynamic Deliveries System) which is used in a daily manner and whose real data have helped to test our algorithm. The main contribution of the paper is twofold: the first one lies on the fact that the algorithm can further powered by adding more constraints other than those presented in the paper since the RVRP centralizes most of the variants of the classical problem; secondly, the algorithm has been properly wrapped in a computer tool that is commercialized by a ICT enterprise with presence in various countries around the world. A large part of the potential of the proposed methods are due to the definition of neighborhood structures which are based on reversing segments of routes and exchanging segments between routes. Furthermore, this algorithm is not only able to solve real-world instances with hundreds of orders efficiently but also computes the solutions within a reasonable time. The methods have been verified by means of several instances extracted from real data provided by a large transport company in Spain. Their performance has been optimized by reducing the total cost of the transport operations. This has been attained by minimizing the total traveled distance, reducing the number of used vehicles and managing the different resources optimally. Future research could aim to use this proposal to solve more general forms of vehicle routing problems which could contain more real-world objectives and constraints. The powerful capacity of this algorithm to find efficient solutions to difficult combinatorial optimization problems should make it useful for solving many other problems in transportation and logistics. Acknowledgments This activity has been partially supported by the Universidad Internacional de La Rioja (UNIR, http://www.unir.net), under the Plan Propio de Investigación, Desarrollo e Innovación [2013–2015]. Research group: Matemática aplicada al mundo real (MAMUR), in part by UNIR Research Support Strategy 2013–2015, under the CYBERSECURITICS.es Research Group [http://research.unir.net] and by the grant MTM2014-52016-C2-1-P. References [1] Á.A. Magreñán, A. Cordero, J.M. Gutiérrez, J.R. Torregrosa, Real qualitative behavior of a fourth-order family of iterative methods by using the convergence plane, Math. Comput. Simul. 105 (2014) 49–61. [2] A. Cordero, L. Feng, Á.A. Magreñán, J.R. Torregrosa, A new fourth-order family for solving nonlinear problems and its dynamics, J. Math. Chem. 53 (3) (2015) 893–910. [3] I.K. Argyros, A. Cordero, Á.A. Magreñán, J.R. Torregrosa, On the convergence of a Damped Secant method with modified right-hand side vector, Appl. Math. Comput. 252 (2015) 315–323. [4] Ministerio de Fomento, Observatorio de transporte de mercancías por carretera, Gobierno de España, Madrid, España, 2011. [5] European Commission, Green Paper: Towards A New Culture for Urban Mobility, European Union, Brussels, Belgium, 2007. [6] P. Ji, K. Chen, The vehicle routing problem: the case of the Hong Kong Postal, Transp. Plan. Technol. 30 (2) (2007) 167–182. [7] R. Lahyani, M. Khemakhem, F. Semet, Rich vehicle routing problems: From a taxonomy to a definition, European J. Oper. Res. 241 (1) (2015) 1–14. [8] J. Oppen, A. Lokketangen, J. Desrosiers, Solving a rich vehicle routing and inventory problem using column generation, Comput. Oper. Res. 37 (7) (2010) 1308–1317. [9] T.G. Crainic, M. Gendreau, W. Rei, A column generation approach for a multi-attribute vehicle routing problem, European J. Oper. Res. 241 (3) (2015) 886–906.
10
J.A. Sicilia et al. / Journal of Computational and Applied Mathematics (
)
–
[10] T. Vidal, T.G. Crainic, M. Gendreau, C. Prins, Heuristics for multi-attribute vehicle routing problems: A survey and synthesis, European J. Oper. Res. 231 (1) (2013) 1–21. [11] M. Gendreau, G. Laporte, J.Y. Potvin, Metaheuristics for the capacitated VRP, in: P. Toth, D. Vigo (Eds.), The Vehicle Routing Problem, SIAM Publishing, Philadelphia, USA, 2002, pp. 129–154. [12] N.A. El-sherbeny, Vehicle routing with time windows: An overview of exact, heuristic and metaheuristic methods, J. King Saud Univ. - Sci. 22 (3) (2010) 123–131. [13] I.M. Chao, T.S. Liou, A new Tabu Search heuristic for the site-dependent vehicle routing problem, in: B. Golden, S. Raghavan, E. Wasil (Eds.), The Next Wave in Computing, Optimization and Decision Technologies, Springer, USA, 2005, pp. 107–119. [14] F.P. Goksal, I. Karaoglan, F. Altiparmak, A hybrid discrete particle swarm optimization for vehicle routing problem with simultaneous pickup and delivery, Comput. Ind. Eng. 65 (2013) 39–53. [15] E.E. Zachariadis, C.T. Kiranoudis, An effective local search approach for the vehicle routing problem with backhauls, Expert Syst. Appl. 39 (3) (2012) 3174–3184. [16] P.P. Repoussis, C.D. Tarantilis, O. Braysy, G. Ioannou, A hybrid evolution strategy for the open vehicle routing problem, Comput. Oper. Res. 37 (3) (2010) 443–455. [17] N. Jozefowiez, F. Semet, E.G. Talbi, An evolutionary algorithm for the vehicle routing problem with route balancing, European J. Oper. Res. 195 (3) (2009) 761–769. [18] I. Or, Travelling salesman-type combinatorial problems and their relation to the logistics of blood banking (Ph.D. Thesis), Department of Industrial Engineering and Management Science, Northwestern University, Evanston, USA, 1976. [19] 0 Bräysy, G. Hasle, W. Dullaert, A multi-start local search algorithm for the vehicle routing problem with time windows, European J. Oper. Res. 159 (2004) 586–605. [20] 0 Bräysy, W. Dullaert, A fast evolutionary metaheuristic for the vehicle routing problem with time windows, Int. J. Artif. Intell. Tools 12 (2) (2003) 153–172. [21] J.Y. Potvin, J.M. Rousseau, An exchange heuristic for routing problems with time windows, J. Oper. Res. Soc. 46 (12) (1995) 1433–1446. [22] P. Hansen, N. Mladenovic, Variable neighborhood search: principles and applications, European J. Oper. Res. 130 (3) (2001) 449–467. [23] M. Tagmouti, M. Gendreau, J.Y. Potvin, A variable neighborhood descent heuristic for arc routing problems with time-dependent service costs, Comput. Ind. Eng. 59 (4) (2010) 954–963. [24] N. Mladenovic, M. Drazic, V. Kovacevic-Vujcic, M. Cangalovic, General variable neighborhood search for the continuous optimization, European J. Oper. Res. 191 (3) (2008) 753–770. [25] E.E. Zachariadis, C.D. Tarantilis, C.T. Kiranoudis, A guided tabu search for the vehicle routing problem with two-dimensional loading constraints, European J. Oper. Res. 195 (3) (2009) 729–743. [26] Visit http://www.alerce.es/es/senda.html for further information about the commercial tool (in Spanish).