Multi Depot Vehicle Routing Problem
with Inventory Constraints
A.J. uit het Broek
s1992503
Abstract
We introduce a new problem called the Multi Depot Vehicle Routing Problem with Inventory con-straints (MDVRP-I). In this problem, vehicles are stationed at several depots, a deterministic set of customers has to be satis ed and each depot has a known inventory level, i.e., the maximum demand served from the depot. In addition, we study the extension in which customers are allowed to be served by more than one vehicle, known as split deliveries. Furthermore, we are the rst to include serving times into the multi depot vehicle routing problem with split deliveries and to elaborate on the e ciency of the corresponding integer linear program. The Multi Depot Vehicle Routing Problem (MDVRP) and the Split Delivery Vehicle Routing Problem (SDVRP) are two special cases of the new problem. A hybrid genetic algorithm is developed that can solve the MDVRP-I and the two special cases. The performance of the heuristic on the two special cases is tested on benchmark instances from the literature, and the results are highly competitive to the current state-of-the-art heuristics. Many best-known solutions were equaled, and for the SDVRP benchmark, new improved solutions were found for 8 out of 48 instances. A new benchmark is created to test the performance of the heuristic on the MDVRP-I as introduced in this paper. The heuristic nds optimal solutions for all instances up to 25 customers that are solved by the integer linear program. Finally, we conducted experiments on the allocation of inventory among depots, and showed that it can have a signi cant in uence on the routing.
1 Introduction
It is increasingly common for customers to order products online. Many retailers try to bene t from this development by allowing online ordering, thereby augmenting their current sales in physical stores. In the highly competitive market, delivery times are of increasing importance since customers are becoming more demanding about service levels. Beside price, they base their decisions on dierent distribution performance criteria like delivery speed, delivery time windows and insights into the status of their purchase's delivery. Customers tend to order their products from rms that have the best delivery services, and excellent e-ful llment seems essential for retailers to bene t from the growing online market. To increase delivery speed and reduce delivery costs, rms can deliver products from a retail store closer to the customer instead of delivering it from a central depot. A secondary bene t is that lost sales might be avoided by maintaining local inventory at the retail store when there is no inventory left at the central depot. Obviously, to be able to deliver a product from a certain depot, or store, this depot must have su cient stock on hand. Especially for retail stores, this may not be the case, and the inventory level will play an important role in the problem.
The well-known Vehicle Routing Problem (VRP), introduced by Dantzig and Ramser (1959), has been studied for decades and is widely applicable in many elds, e.g., transportations, production planning, and telecommunications. In the VRP, all vehicles are dispatched from a single depot, while in the Multi-Depot VRP (MDVRP) vehicles are stationed at several depots. Due to its wide applicability, many articles, literature reviews and books have been written about the VRP and its extensions. Reviewing the current literature reveals that some problem classes have received less attention than others. Examples of less studied extensions are the MDVRP and the Multi-Depot Split Delivery VRP (MDSDVRP), which extends the MDVRP by the possibility that customers will be served by more than one vehicle. Another observation is that most papers consider a speci c problem and are not able to handle a larger class of problems. Furthermore, there is no literature about deterministic depot inventories for the MDVRP scenario and its extensions.
In this paper, a new problem called the Multi Depot Vehicle Routing Problem with Inventory constraints (MDVRP-I) is introduced. In this problem, vehicles are stationed at several depots, a deterministic set of customers has to be satis ed and each depot has a known inventory level. The aim is to minimise the distance travelled by the vehicles while each depot obeys its inventory level, i.e., the maximum demand served from the depot. When customers can be served by more than one vehicle, the problem is called the Multi Depot Split Delivery Vehicle Routing Problem with Inventory constraints (MDSDVRP-I). The MDVRP and the SDVRp are two special cases of the new problem which could be obtained by using an in nite inventory level. Moreover, we identify characteristic properties of the problem with respect to the depot inventory.
A Hybrid Genetic Algorithm (HGA) is proposed that can solve a broad class of problems, including the VRP, MDVRP, and the MDSDVRP-I, with dierent settings like maximum driving times and serving times. Furthermore, the Integer Linear Program (ILP) is studied to obtain exact solutions, and its performance is improved by cuts.
In Section 2 the relevant literature is discussed, and Section 3 introduces the notation and the formal problem de nition as well as the ILP with some possible cuts. A detailed explanation of the HGA is given in Section 4, while the performance of the HGA is given in Section 5. Section 6 concludes the paper, and Section 7 suggest ideas for further research.
2 Literature Review
The literature for the MDVRP starts in 1970s, the rst literature about the Split Delivery VRP (SDVRP) dates from the end of the 1980s, and the MDSDVRP that combines both problems is rst discussed in 2011. This section will review these problems and will end with the genetic algorithm.
of limited size are solved by exact algorithms, and many heuristics have been proposed. The rst eective heuristic was the tabu search given by Cordeau et al. (1997) and later the simulated annealing method described by Lim and Zhu (2006). In the following years, the adaptive large neighborhood search (ALNS) of Ropke and Pisinger (2006) was the state-of-the-art heuristic for solving the MDVRP.
In the past decade, these heuristics have been extended by evolutionary approaches like genetic algorithms. The rst approaches were not competitive to the existing literature, see Thangiah and Salhi (2001) and Ombuki-Berman and Hanshar (2009). More recent papers such as Ho et al. (2008) gave competitive results, and the method discussed by Vidal et al. (2012) is superior to the previous state-of-the-art methods, which indicates that an evolutionary approach is a promising method to tackle the MDVRP. Shortly after this progression, the evolutionary approaches were already being used for extensions of the MDVRP as in Roodbergen and Kolman (2016). A recent literature review for the MDVRP is given by Montoya-Torres et al. (2015).
The idea of obtaining some savings by split deliveries was introduced by Dror and Trudeau (1989) and further developed by Dror and Trudeau (1990), who also described some fundamental properties of optimal solutions of the SDVRP. One of these properties states that there always exists an optimal solution in which no two routes have more than one split customer in common. Although the SDVRP is a relaxation of the VRP, Dror et al. (1994) show that obtaining the optimal solution is harder. Furthermore, the eect of split delivery is most present for problems with high demand customers, i.e., customers with a demand higher than 10% of the vehicle capacity. Later on, the SDVRP is extended with serving times. Mostly, it is assumed that the serving time of customers is independent of the quantities that are delivered by a vehicle; see Gendreau et al. (2006), and Desaulniers (2010). More recently, Vacca and Salani (2011) proposed making the serving time and the delivered quantity dependent.
Several exact methods and heuristics have been proposed to solve the SDVRP. Archetti and Speranza (2006), and Aleman and Hill (2010) solved the SDVRP with a tabu search algorithm and Jin et al. (2007) used a two-stage algorithm. A constructive heuristic together with an adaptive memory algorithm are given by Aleman et al. (2010), and the new constructive heuristic of IV and Cavalier (2012) can quickly create initial feasible solutions for the SDVRP. Most recently, a new state-of-the-art heuristic has been described by Silva et al. (2015), who introduced an iterative local search. Several papers improve the LP relaxation by cuts; see Archetti et al. (2011), Vacca and Salani (2011) and Archetti et al. (2014). All these approaches are aimed speci cally at the SDVRP cannot handle other extensions of the VRP, whereas our proposed HGA can solve other extensions of the VRP too. A more extended review is given by Archetti and Speranza (2008) and Gulczynski et al. (2008).
The MDSDVRP was rst proposed by Gulczynski et al. (2011), who used a two-stage heuristic based on integer programming that rst assigns customers to depots and then solves the SDVRP per depot. Ray et al. (2014) designed a new heuristic that takes into account the geographical placement of the depots. To the best of our knowledge, these are the only two pieces of literature discussing the MDSDVRP, neither of which discusses the e ciency of the ILP or takes into account serving times.
The idea of a Genetic Algorithm was rst described by Holland (1975) and belongs to the class of adaptive stochastic optimisation algorithms. This class is a subset of the evolutionary algorithms, which mimic the natural process from the biology eld. An extensive literature review focused on GAs for the MDVRP is given by Karakati£ and Podgorelec (2015), who conclude that GAs are well applicable for the MDVRP. To the best of our knowledge, there exists no literature on a VRP variant with inventory restrictions at the depots. One eld that combines distribution with inventory management over time is called the Inventory Routing Problem (IRP). The IRP considers multiple periods in which each customer has a deterministic demand for each period. The main characteristic of the IRP is that each customer has a maximum local storage capacity and the depots have an in nite amount of inventory. Hence, the IRP literature assumes in nite inventory levels at the depots as well.
classes. Moreover, the literature about the MDSDVRP is scarce, and there is no literature that studies the e ciency of the ILP formulation for the problem or takes into account serving times. To ll those gaps, we introduce inventory constraints and propose a HGA that is able to solve all these problems, including serving times.
3 Problem Formulation
This section formally states the MDSDVRP-I and introduces the notation. The second half introduces the Integer Linear Program (ILP) and provides some possible cuts.
The VRP consists of a set of customers Vc and a set Vd with a single depot. Let G = (V, A) be a
complete directed graph where V = Vd∪ Vcand the sets Vdand Vcare mutually exclusive and collectively
exhaustive. The arc set A ⊆ V × V is assumed to be complete and of Euclidean distance denoted by cij: N×N → R. A homogeneous eet of m ∈ N vehicles with capacity Q ∈ N and maximum working time
T ∈ R+is stationed at the depot. Let K be the set of m vehicles for a depot. Furthermore, each customer
has a nonnegative demand qi∈ N and a nonnegative service duration τi∈ R+. The duration of a single
route is the sum of the total driving time and the total service time to serve all customers included in the route. The objective is to nd a set of vehicle routes visiting all customers which minimizes the total travel time and obeys the capacity and route duration constraint.
The MDVRP extends the VRP to a setting with ηd ∈ N identical depots where each depot has a
homogeneous eet K of m ∈ N vehicles with capacity Q ∈ N and maximum working time T ∈ R+. For
ease of notation, de ne ˜K ≡ (K × Vd)as the set of all vehicles in the instance. The MDSDVRP adds the
possibility of serving a customer by several vehicles, which approach is called split deliveries. This allows for customers with a demand higher than the vehicle capacity and may result in lower objective values. The MDVRP-I shares all assumptions and settings of the MDVRP. In addition, there is a maximum available inventory at each depot, denoted by Ij ∈ N for j ∈ Vd. When split deliveries are allowed the
problem is called MDSDVRP-I. The opportunity for split deliveries allows for customers with a demand higher than the highest depot inventory.
3.1 Integer Linear Program
Let xijk`be a binary variable that indicates whether arc (i, j) ∈ A is used by vehicle k ∈ K which departs
from depot ` ∈ Vd. Let cij be the costs of traversing arc (i, j). Next, let yik` be a non-negative integer
variable that indicates the amount vehicle k, which departs from depot `, serves to customer i. The ILP model for the MDSDVRP-I is formulated as,
X j∈V xijk` = 0 ∀` ∈ Vd, i ∈ Vd\{`}, k ∈ K (8) X j∈V xjik`− X j∈V xijk` = 0 ∀i ∈ V, k ∈ K, ` ∈ Vd (9) 2 ≤ uik`≤ |Vc| + 1 ∀i ∈ Vc, k ∈ K, ` ∈ Vd (10)
uik`− ujk`− 1 ≥ |Vc| · (xijk`− 1) ∀i ∈ Vc, j ∈ Vc\{i}, k ∈ K, ` ∈ Vd (11)
yik`∈ N ∀i ∈ Vc, k ∈ K, ` ∈ Vd (12)
uik`∈ N ∀i ∈ Vc, k ∈ K, ` ∈ Vd (13)
xijk` ∈ {0, 1} ∀(i, j) ∈ V, k ∈ K, ` ∈ Vd (14)
Objective (1) minimizes the total transport costs. Constraint (2) implies that each customer is served completely, and constraints (3), (4), and (5) enforce limits on the inventory usage, the vehicle capacity and the route duration. Constraint (6) ensures that a customer only receives products from a vehicle that visits the customer. Constraints (7) and (8) make sure that each vehicle only leaves a depot once and each vehicle can only leave its corresponding depot. The ow conservation is modelled by constraint (9), and constraints (10) and (11) eliminate all subtours.
The previous ILP notation is based on two other ILP models. The rst one is a formulation for the MDVRP which takes into account serving times and vehicle capacities, see Vidal et al. (2012). The second one includes split deliveries into the MDVRP but has no serving times, see Ray et al. (2014). We have combined these and added inventory constraints.
3.2 Cuts for the ILP Formulation
This section provides three possible cuts to improve the lower bound given by de LP-relaxation, which can improve the performance of exact methods to solve the ILP. First, the subtour elimination constraints are considered. Second, a theorem that holds for the SDVRP is shown to hold for the MDSDVRP-I. Third, a cut speci c for the SDVRP is rewritten to a cut for the MDSDVRP-I. Fourth, a new cut is introduced based on the given theorem.
The given ILP formulation uses the basic variant of the MillerTuckerZemlin (MTZ) based constraints to eliminate subtours, as introduced by Miller et al. (1960). However, for the asymmetric-TSP this notation leads to a relative deviation between the LP relaxation bounds and the corresponding optimal values of 7.39%, as shown in a comparison by Oncan et al. (2009). The lifted variant described by Sherali and Driscoll (2002) only leads to a deviation of 4.39%. A formulation of the lifted MTZ constraints for the multiple travelling salesman problem, which is a VRP without capacity constraint, is given by Kara and Bektas (2006). These constraints are rewritten to our formulation as follows,
uik`− ujk`+ |Vc| · xijk`+ (|Vc| − 2) · xjik`≤ |Vc| − 1 ∀i ∈ Vc, j ∈ Vc\{i}, k ∈ K, ` ∈ Vd
The sub-tour elimination constraints might be further improved by the multi-commodity ow notation introduced by Sherali et al. (2006) which only has a deviation of 0.43%, according to Oncan et al. (2009). For the single depot SDVRP there are several de nitions and theorems that can be used to improve the ILP formulation. First the de nition of a k-split cycle is given, followed by the theorems for the SDVRP and the MDSDVRP-I.
De nition Consider a set C = {ii, i2, . . . , ik}of customers and suppose there exist k routes r1, r2, . . . , rk
with k ≥ 2, such that rwcontains customers iwand iw+1, where w = 1, . . . , k−1 and rkcontains customers
i1 and ik. Such a scenario is called a k-split cycle.
Theorem 1. If the distance matrix satis es the triangle inequality, then there exists an optimal solution to the SDVRP in which no two vehicles have more than one split customer in common.
Lemma 1. If the costs satisfy the triangle inequality, then there exists an optimal solution to the SDVRP which does not contain k-split cycles, where k ≥ 2.
Lemma 2. If the costs satisfy the triangle inequality, then there exists an optimal solution to the SDVRP such that each edge between two customers is traversed at most once.
Lemma 3. If the costs satisfy the triangle inequality, then there exists an optimal solution to the SDVRP such that the total number of splits is lower than the number of routes.
Theorem 1 is proven by Dror and Trudeau (1989), Lemma 1 is shown by Dror and Trudeau (1989), Lemma 2 is shown by Desaulniers (2010), and Lemma 3 is shown by Archetti et al. (2008). Note that the lemmas follow directly from Theorem 1. Furthermore, Theorem 1 is easily extended to the multi depot setting since the original proof does not use properties of the depot. It follows that the lemmas hold for the MDSDVRP as well. We can extend this to the case with inventory constraints as shown by the following theorem.
Theorem 2. If the distance matrix satis es the triangle inequality, then there exists an optimal solution to the MDSDVRP-I in which no two vehicles have more than one split customer in common.
Proof. Suppose we have a feasible solution with two customers a and b who receive all their demand qa
and qbin a split by the same two vehicles, called vehicle 1 and 2. Let q1a and q2abe de ned as the amount
delivered to customer a by vehicle 1 and 2, respectively. q1b and q2b are de ned similarly for customer b.
From the triangle inequality, it immediately follows that the driving costs incurred by a vehicle decrease if the number of customer visits of that vehicle are decreased by one. One has to prove that one of the four deliveries, q1a, q1b, q2a, or q2b, can be removed without violating any constraint.
The vehicle capacity constraint is already considered by Dror and Trudeau (1989), and the time needed by a vehicle only decreases if customers are removed from a tour; hence, the maximum driving time constraint is not violated. It is left to be proven that the inventory constraints of the depots are not violated.
Let qmin= min{q1a, q1b, q2a, q2b}, and without loss of generality, assume qmin= q1a. Vehicles 1 and 2 can
exchange some loading such that vehicle 2 takes over the loading of customer a, which equals an amount of q1a, and vehicle 1 takes over the same amount from customer b. The new delivery quantities for the
vehicles are, qnew1a = q1a− q1a = 0 qnew1b = q1b+ q1a qnew2a = q2a+ q1a qnew2b = q2b− q1a which implies, qnew1a + qnew1b = q1a+ q1b qnew2a + qnew2b = q2a+ q2b
It follows that vehicle 1 only serves customer b and thus one split delivery is eliminated, while the total loading of both vehicles stays unchanged. Hence, the inventory constraint of both depots is satis ed. Lemma 4. If the costs satisfy the triangle inequality, then there exists an optimal solution to the MDSDVRP-I which does not contain k-split cycles, where k ≥ 2.
Lemma 5. If the costs satisfy the triangle inequality, then there exists an optimal solution to the MDSDVRP-I such that each edge between two customers is traversed at most once.
Proof. Follows directly from Theorem 2
Lemma 6. If the costs satisfy the triangle inequality, then there exists an optimal solution to the MDSDVRP-I such that the total number of splits is lower than the number of routes.
Proof. Follows directly from Theorem 2
Desaulniers (2010) describes a valid inequality for the SDVRP to increase the quality of the LP relaxation. The inequality is derived from Lemma 2 and is proven to be eective in the single commodity ow formulation. The inequality is rewritten to the MDSDVRP-I case with aid of Lemma 5 and equals,
X k∈K X `∈Vd xijk`+ X k∈K X `∈Vd
xjik` ≤ 1 ∀i ∈ Vc, j ∈ Vc\{i}
Theorem 2 can be used to de ne a new cut that avoids routes with more than one split delivery in common. First one has to identify all combinations of customers visited by vehicle k1`1 which is accomplished by
P
i∈Vcxij1k1`1 +
P
i∈Vcxij2k1`1 ≥ 2 where (j1, j2) ∈ Vc× Vc and j1 6= j2. All other vehicles k2`2 might
only visit one of the customers (j1, j2) ∈ Vc× Vcwhich is modeled as Pi∈Vcxij1k2`2+
P
i∈Vcxij2k2`2 ≤ 1.
The two constraints are combined as, X i∈Vc xij1k1`1+ X i∈Vc xij2k1`1 ! + X i∈Vc xij1k2`2+ X i∈Vc xij2k2`2 ! ≤ 3 ∀j1∈ Vc, j2∈ Vc\{j1} ∀(k1, `1) ∈ ˜K (15) ∀(k2, `2) ∈ ˜K\{(k1, `1)}
Recall the de nition ˜K ≡ (K × Vd). According to Lemma 4 the previous cut could be generalized to avoid
k-split cycles for k ≥ 2.
4 The Genetic Algorithm
Parent selection via tournament Crossover Inherit mix Inherit parent 2 Inherit parent 1 Make complete Education Individual Evaluation Create 4µ individuals to fill population Maximum number of childs without improvement Max number Iterations or time limit Stop Yes Kill all but 1/3 µ of both subpopulations
Update penalty weights evaluation Feasible
Yes probability PrepRepair with
Multiple of 100 iterations Yes No No No Start Survivor Selection Yes No
Figure 1: Flow chart of the overall design of the proposed HGA.
a xed number of ospring solutions has been created. The general overview of the HGA is graphically shown in Figure 1.
The search space and data structure are explained in Section 4.1, and the evaluation of individuals is described in Section 4.2. Section 4.3 gives the crossover procedures followed by the educational phase in Section 4.4. The population initialisation, the selection procedure, the repair phase and the diversify procedure are considered in detail in Section 4.5.
4.1 Search Space and Solution Representation
Exploration of both feasible and infeasible solutions may lead to better quality solutions, as shown by Glover and Hao (2011). Therefore, the search space S is de ned as the set of feasible and infeasible solutions, where infeasible solutions are obtained by relaxing the vehicle capacity, the maximum driving time and the inventory restriction. The number of vehicles per depot is not relaxed since preliminary analysis shows it is challenging to decrease the number of vehicles once they are used.
Let R(σ) be the set of routes in the solution σ ∈ S, let Rj(σ)denote the subset of routes departing from
depot j ∈ Vd and let Vcr be de ned as the subset of customers visited by route r ∈ R(σ). Each route r
departs from a depot vr
0∈ Vd, visits a sequence of nrcustomer locations vr1, v2r, . . . , vrnr ∈ Vc and returns
to the originating depot vr
nr+1 = v
r
0. A route r is described by its driving time cr(σ) = P nr
i=0cvr ivi+1r ,
its total duration τr(σ) = cr(σ) +Pi∈Vr
cτi, and its load qr(σ) =
P
i∈Vr cq
r
i(σ)where qri(σ)denotes the
amount vehicle r delivers to customer i.
Since infeasible solutions are allowed, there is a penalty for exceeding the depot inventory, the vehicle capacity, and the maximum driving time which equals ωI, ωQ, and ωD, respectively. Moreover, there is a
penalty scaling parameter ω used in the repair phase of the HGA, see Section 4.5. The total costs φr(σ)
The costs φ(σ) of an individual σ equals the routing costs plus the penalty for depot inventory excess, φ(σ) = X r∈R(σ) φr(σ) + ω · ωI· X j∈Vd X r∈Rj(σ) qr(σ) − Ij + .
A solution σ ∈ S is represented by four chromosomes. The rst chromosome is a vector of routes R(σ) where each route r ∈ R(σ) is represented by a sub-chromosome, see Figure 2c. Each route stores additional information like the load of the vehicle qr(σ), the duration τr(σ), the total length cr(σ), the
total penalized costs φr(σ), the initial inventory of the depot Ir(σ), and whether the route itself is feasible.
The route pattern contains a city index i and the delivery amount qr
i(σ)of each visit. The chromosome
has a route for each allowed vehicle even if the vehicle is unused.
The second chromosome is the depot assignment chromosome, denoted as χi(σ), and stores the depot
assignment for each customer i ∈ Vc, see Figure 2b. Third, there is a chromosome to store the remaining
inventory Ij(σ)of depot j ∈ Vd. The last chromosome stores for each split delivery k a set with the route
and visit indices, denoted as ϕk(σ). Moreover, an individual σ stores additional information like the total
length, the total penalty including the inventory penalty, and whether it is feasible with respect to the inventory limitation, the capacity constraints and the maximum driving time constraints.
The data structure allows the crossover procedure to easily copy routes to the child and enables each operator to ll unused vehicles. During education, the inventory chromosome enables a fast feasibility check whereas the split chromosome is used to nd all split deliveries in O(1). The depot assignment chromosome is used to nd the dierence between two solutions in O(ηc), as explained in Section 4.2.
4.2 Individual Evaluation
The individual evaluation determines the quality of a solution with respect to the other solutions in the population. The evaluation leads to a so-called biased tness of the individual that is used to select parents for the crossover phase (Section 4.3) and to determine which individuals will live on in the next generation (Section 4.5), these are called survivors. Diversity of the population is essential to allow for a broad search perspective. Therefore, the biased tness is based on both the total costs and the diversity contribution of a solution. Taking into account the diversity contribution of an individual may avoid that the population converge to a local optimum in an early stage.
The diversity contribution ∆(σ) of an individual σ is de ned as the average distance to its nclose ∈ N
closest neighbors. The distances between two individuals σ1and σ2is expressed as the dierence in depot
assignment, see equation (16). The closest neighbor for σ1 is the individual with the smallest distance
δ(σ1, σ2)and the set of the nclose closest neighbors is called Nclose.
δ(σ1, σ2) = 1 ηc X i∈Vc Iχi(σ1)6=χi(σ2), (16)
where ηc is the number of customers and I is an indicator function indicating whether the depot
assign-ment χi(σ)of customer i is dierent in σ1 and σ2. In other words, the distance between two solutions is
expressed as the number of customers that have a dierent depot assignment between the solutions. The diversity contribution ∆(σ1)is now stated as,
∆(σ1) = 1 nclose X σ2∈Nclose δ(σ1, σ2). (17)
The nal evaluation value, or biased tness, is based on Roodbergen and Kolman (2016) and Vidal et al. (2012). It is expressed as ranks of the total costs and the diversity contribution of an individual,
a b c d e f g h i 3 6 6 3 1 3 7 6 2 3 2
(a) Possible solution for MDSDVRP-I exam-ple with a split delivery for customer d
Depot i Ii(σ) a 5 b 1 Split k ϕk(σ) 1 {(1, 3), (3, 2)} Customer i χi(σ) c a d (a, b) e a f b g b h b i b
(b) Inventory chromosome (left top), split chromosome (right top), and depot assignment chromosome (bottom)
R(σ) r = 1 r = 2 r = 3 r = 4 c1(σ) = 18 q1(σ) = 5 a c, 2 d, 1 e, 2 a c2(σ) = 0 q2(σ) = 0 a a c3(σ) = 17 q3(σ) = 5 b d, 1 f, 2 g, 2 b c4(σ) = 7 q4(σ) = 4 b h, 2 i, 2 b
(c) Pattern chromosome. For clarity only ci(σ)and qi(σ)are shown.
Figure 2: The chromosome representation for an arbitrarily feasible MDSDVRP-I solution. The example has two depots with both two vehicles with capacity ten. Next, each customer demands two items and the inventory of each depots is ten items.
where ξtc ∈ R and ξdc ∈ R are dynamic weight parameters, see Roodbergen and Kolman (2016), and
nelite ∈ N is the number of elite individuals. Elite individuals are individuals with the best objective
value and will survive the selection procedure for sure. Notice that the best individual has the lowest evaluation value. Both the ranks and the biased tness are continuously updated after each new created ospring. The dynamic weight parameters are de ned such that the weight of the cost increases and the weight of the diversity contribution decreases during the execution of the algorithm. These parameters are de ned similar as in Roodbergen and Kolman (2016),
ξtc= 0.5 + α · ITactual ITmax and, ξdc= 0.5 − α · ITactual ITmax
where ITactual is the current iteration, ITmaxis the maximum number of iterations the HGA considers
4.3 Reproduction Procedures
A tournament selection as given in Miller and Goldberg (1996) is used to select parents for recombination. For a tournament we randomly select individuals from the population. The individual with the best biased tness wins the tournament and will be a parent for the next child. The size of the tournament is a way to control the selection pressure and is introduced as a parameter, called nparents. Large tournaments
favor tter individuals, which will increase the convergence rate of the HGA. Note that both feasible and infeasible individuals may win the tournament. In this way, we will especially search at the borders of feasibility since the high quality solutions mostly are found there.
We propose a new crossover based on the periodic crossover with insertions (PIX) described by Vidal et al. (2012). The crossover may transmit good groups of customers' assignments for a depot while enabling depot-customer recombinations. The aim is to have a versatile crossover that allows for a wide exploration of the search space and small re nements of high quality solutions. The rst is ensured by the possibility that a child inherits material from both parents in an approximately equal proportion, and the second is ensured by the possibility that a child inherits mostly from one parent.
Special attention goes to the inheritance of split deliveries as the crossover must be able to inherit split deliveries, including between dierent depots, while it avoids an excessive increase of split deliveries. The crossover uses an auxiliary data-structure that indicates how much unserved demand each customer in the child solution has, called ˆqi. At the start of the crossover phase, we have ˆqi= qi, ∀i ∈ Vc.
First, the depots are randomly divided into three mutually exclusive sets, called Λ1, Λmix and Λ2. The
depots in Λ1 are inherited from parent 1, the depots in Λmixare partially inherited from the rst parent
and the depots in Λ2are inherited from the second parent. For each depot in Λ1we choose to inherit all
corresponding routes, and for the depots in Λmixa random number of corresponding routes is inherited.
Split deliveries are inherited as well and are not restored to complete visits. Note if a part of a split customer i is served by a depot in Λ1∪ Λmixand the remaining demand is served by a depot in Λ2, we
have 0 < ˆqi ≤ qi. In this scenario, the split customer i is not fully inherited from the rst parent. If all
parts of split delivery i were served by depots in Λ1∪ Λmix, it is most likely that the split delivery is
completely inherited, which implies ˆqi= 0, however this is not necessary.
The next step is to inherit the routes corresponding to the depots in Λ2, which will be inherited from
the second parent. First, the customers who are already, partially or completely, inherited from the rst parent are removed from the routes. Then all routes are inherited by the child.
At the end of the crossover, it might occur that there are some customers with unsatis ed demand. The unsatis ed demand is inserted with the cheapest insertion heuristic. The heuristic takes into account the penalty costs and is allowed to insert the demand as a split over two routes. This can eliminate split deliveries that were present in one of the parents. When split deliveries are not allowed, the unsatis ed demand is inserted as a complete visit, and penalty costs are not taken into account.
Preliminary analysis shows that the proposed crossover is eective for a multi depot setting but lacks performance for a single depot instance. Therefore, the crossover works on a vehicle level instead of a depot level when the instance has only one depot. All other steps stay unchanged.
4.4 Education
The crossover phase is followed by an educational phase that improves the quality of the ospring solution. The operators can handle moves between dierent depots, which contrasts with Vidal et al. (2012). The rst operators are the well-known VRP operators, whereas the other operators are dedicated to optimising split deliveries.
Let u, v ∈ Vc be customers and let r(u) denote the route containing customer u. The partial route
between the customers u1, u2∈ {u1, u2 ∈ Vc× Vc | r(u1) = r(u2)} is denoted as (u1, u2). Next, let x be
the successor of u in r(u) and y be the successor of v in r(v). Note that the successors can be customers as well as depots. The operators require u 6= y and v 6= x, and when a successor location, i.e. x or y, is used by an operator this has to be a customer location.
The rst operator is relocate. It removes customer u and inserts it after customer v. The dierence in capacity and driving time can both be calculated in constant time. The dierence in inventory is only taken into account if vr(u)
0 6= v
r(v)
0 and due to the inventory chromosome the inventory penalty is
calculated in constant time too. There are two additional variants of the relocate operator. The rst one removes (u, x) and inserts the visits after v, the second one removes (u, x) and inserts the visits in reversed order after v. The possible moves with the relocate operators are shown in Figure 3.
The fourth operator, called swap, interchanges two customer visits u and v. There are three additional swap operators that swap (u, x) with v, (u, x) with v while reversing (u, x), and (u, x) with (v, y), respectively. A single move of these operators runs in O(1), thus to nd a rst improvement takes O(η2
c). The possible moves with the swap operators are shown in Figure 4.
Before relocate r(u) r(v) u v After relocate u v
(a) Single relocate
Before relocate r(u) r(v) u x v After relocate u x v (b) Multi relocate Before relocate r(u) r(v) u x v After relocate x u v
(c) Multi relocate with reverse Figure 3: All possible relocate moves. These moves do not require r(u) 6= r(v).
Before swap r(u) r(v) u v After swap u v (a) Single swap
Before swap r(u) r(v) u x v After swap u x v (b) Multi swap I Before swap r(u) r(v) u x v After swap x u v
(c) Multi swap I with reverse
Before swap r(u) r(v) u x v y After swap u x v y (d) Multi swap II Figure 4: All possible swap moves. These moves do not require r(u) 6= r(v).
The eighth operator is 2-opt which has three dierent implementations. The rst implementation is an intra-route 2-opt move that replaces (u, x) and (v, y) by (u, v) and (x, y). The second one is an inter-route move and replaces (u, x) and (v, y) by (u, v) and (x, y), the third one is an inter-inter-route move as well and replaces (u, x) and (v, y) by (u, y) and (x, v). The intra-route 2-opt requires r(u) = r(v) while the inter-route moves require r(u) 6= r(v) and vr(u)
0 = v
r(v)
0 . The eect of the intra-route move is found
in O(1) while the inter-route moves need O(ηc)time, the total running time of the moves is O(ηc2) and
O(η3
c), respectively. The possible 2-opt moves are shown in Figure 5.
Before 2-opt r(u) r(v) u x v y After 2-opt u x v y
(a) Intra-route 2-opt
Before 2-opt r(u) r(v) u x v y After 2-opt u x v y (b) Inter-route 2-opt I Before 2-opt r(u) r(v) u x v y After 2-opt u x v y (c) Inter-route 2-opt II
Figure 5: All possible 2-OPT moves. The intra-route move requires r(u) = r(v) whereas the inter-route moves require r(u) 6= r(v) where both tours start at the same depot.
The last operator is the splitReinsertion which is an adaptation of the operator described by Silva et al. (2015), see Algorithm 1. The parameters for the operator are a customer u and a given forbidden route ¯
r, by default ¯r = r(u). The forbidden route is used to avoid a intra-route relocate move.
First, the dierence in length and penalty are calculated for removing customer u from the current route. Second, de ne L = ∅ as a set of routes and each route r ∈ R(σ)\{¯r} for which an insertion is feasible, i.e., q(r) < Q and Ivr
0 > 0, is added to the set L. The operator is terminated if the residual capacity
of the vehicles in the set L is insu cient to have a feasible split insertion of customer u. Next, for each route r ∈ L the cheapest insert location for customer u is found, the cost of this insertion is called br.
A small modi cation of the greedy knapsack heuristic is used to determine the routes in which customer u is reinserted. Let ar = min{Q − q(r), Ivr
0}, then the routes in L are sorted in increasing order for
br/ar. Lastly, customer u is inserted in the rst route in L, the next route in L is considered as long as
the demand of customer u is not satis ed. Moreover, for each route it is checked whether the condition Ivr
0 > 0still holds. Notice that this operator may result in a inter-route relocate.
Instead of the greedy knapsack heuristic one can use an e cient dynamic programming. However, as described by Boudia et al. (2007) the greedy knapsack heuristic is approximately 30% faster and gives an optimal solution in 66% of the cases. Furthermore, for their instances there is no signi cant dierence in the nal solution quality between the two methods.
Algorithm 1 splitReinsert
Input: A complete graph G = (V, A), solution σ ∈ S, customer u ∈ Vc and a forbidden route ¯r ∈ R(σ)
Output: Same or improved solution
1: L ← A ← U ← ∅ 2: sumResidual ← 0 3: for r ∈ R(σ)\{¯r} do 4: ar← min{Q − q(r), Ivr 0} 5: if ar> 0then 6: L ← L ∪ {r} 7: A ← A ∪ ar 8: sumResidual ← sumResidual + ar 9: end if 10: end for 11: 12: if sumResidual < qi then 13: return σ 14: end if 15: 16: for r ∈ L do
17: Let ur be the cheapest insert costs for customer i in route r
18: U ← U ∪ ur
19: end for
4.5 Population Management
The HGA uses two separated subpopulations for feasible and infeasible solutions that are independently managed. Both have a minimum size µ and a maximum size µ + λ. After the education of a new created ospring, it is immediately inserted into the corresponding subpopulation. The remainder of this section explains the initialisation of the population, the diversi cation procedure, the repair phase, the updating procedure for the penalty weight parameters and the survivor selection.
As initialisation 4µ individuals are created by randomly assigning customers to depots and then creating random routes for the given depot assignment. Each initial individual is educated, but the operators are not allowed to make moves between depots. In this education each operator is at most considered ninitEduc times. Afterwards, the individual is inserted into the corresponding subpopulation according to
its feasibility. When a subpopulation reaches its maximum size, the survivor selection procedure reduces the size to µ individuals, as described in the subsequent sections. Note, at the end of the initialisation it may occur that one subpopulation has less than µ individuals.
The diversi cation procedure is activated when the best individual is not improved in the last itdiv
iterations. The procedure kills all but the ttest µ/3 individuals of both sub-populations. Next, 4µ new individuals are created as in the initialisation phase.
When an ospring solution is feasible after it has education it is called naturally feasible, otherwise it is called naturally infeasible. According to its natural feasibility the ospring solution is added to its corresponding subpopulation. Each naturally infeasible solution is repaired with probability Prep and if
the repair phase succeeds the solution is added to the feasible subpopulation while the original infeasible solution is not removed from the infeasible subpopulation. The repair phase re-educates the individual with ω = 10, when this results in another infeasible solution the individual is re-educated with ω = 100. Due to this temporarily compelling increase of penalties for infeasibility the education phase is biased towards feasible solutions.
The penalty parameters are dynamically updated during runtime to guide the HGA search along the boundary of feasibility. The initial values are ω = ωI = ωD = 1.0 and ωQ = ¯c/¯q, where ¯c denotes the
average distance between customer locations and ¯q equals the average demand. The desired proportion of naturally feasible individuals is denoted as ξref while ξI, ξD, and ξQ denote the realized proportion of
naturally feasible individuals in the last hundred ospring individuals with respect to the depot inventory, the maximum driving time and the vehicle capacity, respectively. The parameters are updated every 100 iterations as follows, where β ∈ {I, Q, D}:
If ξβ≤ ξref− 0.05, then ωβ = ωβ· 1.20
If ξβ ≥ ξref+ 0.05, then ωβ = ωβ· 0.85
The nal procedure is the survivor selection that determines the µ individuals that will survive to the next generation and eliminates the others. To avoid premature convergence, it is essential to ensure di-versity of both subpopulations, however, to evolve to good quality solutions, the most promising solutions should survive. The survivor selection favors diversity by using the biased tness function while the best individuals in terms of costs, called elite individuals, will survive for sure. The number of elite individuals is introduced as a parameter called nelite.
5 Numerical Analysis
This section describes the parameter calibration, and next, the performance of the HGA is analysed. A random instance set, see Section 5.1, is used for the parameter calibration whereas general benchmark sets from the literature are used to test the performance of the algorithm under dierent settings. Furthermore, the HGA is benchmarked versus the ILP model given in Section 3.1.
The ILP model is coded in CPLEX 12.6 whereas the HGA is coded in C++14 which is compiled with "gcc -march=core-avx2 std=c++14 -g - to -O2". All tests are performed on a 64bit Intel i5-4670K 3.40GHz CPU with 8GB RAM. Moreover, CPLEX 12.6 is allowed to run in parallel (up to four threads) and use MIP-pooling techniques, the HGA is only allowed to run sequential on one core.
5.1 Instances
The instances have ηccustomers with ηc∈ {25, 50, 100, 150, 200, 250}, and ηddepots with ηd ∈ {1, 2, 4, 6}.
Inspired by Cordeau et al. (1997) we create the test instance set as follows. First, randomly locate ηd
depots in the [−50, 50]2 square according to a discrete uniform distribution. Then randomly locate a
customer i in the [−100, 100]2 square according to a discrete uniform distribution. Find the distance to
the closest depot, called ˜`i, and discard the customer with probability 1 − e−0.05˜`i. This step is repeated
until there are ηccustomers. The serving time τiand the demand qiare randomly chosen with a discrete
uniform distribution in the interval [1, 25]. The exponential probability to discard customers leads to instances with customers clustered around depots. Let Ψ = Pi∈Vcqi denote the total demand in an
instance. The total inventory is equally allocated among the depots and the inventory per depot equals dψ · Ψ
ηde, where ψ indicates the inventory tightness.
The next settings are the capacity and the maximum driving time of the vehicles. We aim for solutions with ten to twenty customers per vehicle where some instances are restricted by the capacity and others by the time constraint. Three dierent settings for the vehicle capacity and the maximum driving time are used, see Table 1. The number of vehicles per depot is determined as dηc
t·Ψe, where t is the expected
average amount of customers per vehicle.
1 2 3
t 10 15 20 T 325 450 485 Q 150 200 250
Table 1: Settings for the random instances. t is the expected average number of customers in a vehicle, T is the maximum driving time of a vehicle and Q equals the capacity of a vehicle.
5.2 Parameter Calibration
The performance of the HGA depends on a correlated set of parameters. A set of parameter values is obtained by training the HGA on a restricted set of 144 training instances, see Section 5.1. Like Vidal et al. (2012), the calibration is repeated for each problem class. The parameters to calibrate and the considered ranges, which are based on values from existing literature, are shown in Table 2.
Preliminary analysis shows that the parameters ncloseand nelite should depend on the population size µ
while nitDivshould depend on the maximum number of iterations ITmax. Therefore, three new parameters
are introduced. Let nc be de ned such that nclose= nc · µ, let el be de ned such that nelite= el · µ and
let nd be de ned such that nitDiv = nd · ITmax.
The training starts with the nal parameter values of Vidal et al. (2012), ninitEduc = 50, nparents = 2
and ITmax = 50, 000. Then each individual parameter is increased and decreased, and the parameter
Parameter Range SDVRP MDVRP MDSDVRP-I µ [5, 200] 25 25 25 λ [1, 200] 85 85 85 nparents [2, 10] 2 2 2 nc [0, 0.25] 0.20 0.20 0.20 el [0, 1] 0.45 0.45 0.45 nd [0, 0.5] 0.30 0.30 0.30 Prep [0, 1] 0.60 0.60 0.60 ξref [0, 1] 0.30 0.30 0.30 ξtc [0, 1] 0.50 0.50 0.50 ξdc [0, 1] 0.50 0.50 0.50 ninitEduc [0, 1000] 50 100 75 α [0, 0.5] 0.10 0.10 0.10
Table 2: Results of the parameter calibration for the dierent problem classes. The range indicates the values that are considered in the calibration whereas the nal values are shown in the subsequent columns.
5.3 Benchmark
Two main aspects of the MDSDVRP-I are the multi depot setting and the possibility to have split deliveries. Benchmark instances from the literature are used to determine how the HGA performs on these two aspects. For each benchmark, the HGA is compared to the current state-of-the-art heuristic. All instances are solved ten times, the minimum and average objectives are reported and the time denotes the running time up to the last improvement. Furthermore, the average deviation from the Best Known Solution (BKS) is reported for all benchmarks. Objective values equal to the BKS are indicated in boldface, while objective values lower than the BKS are boldfaced and underlined.
The rst benchmark consist of instances for the SDVRP designed by Belenguer et al. (2000). The set consists of 25 problems where each vehicle has a limited capacity and no maximum driving time. The rst 11 instances are modi cations from the standard TSPLIB Reinelt (1991) whereas the others are generated as follows. In each instance the vehicles have a capacity of Q = 160 and the coordinates in S51D*, S76D*, and S101D* are equal to the coordinates in eil51, eil76, and eil101, respectively. The demand of the customers is generated as in Dror and Trudeau (1989). The instances are used for two dierent scenarios, namely one with a limited eet size and one with an unlimited eet size.
The results for the limited and unlimited eet size are compared with the results of Aleman and Hill (2010), Archetti et al. (2011), and Silva et al. (2015), see Table 3 and Table 4. For the unlimited eet the average gap between the BKS and the HGA is 0.13% while for the state-of-the-art method it is 0.15%, for the minimum of ten runs these values equals -0.01% and +0.01%, respectively. In addition, the HGA found new best known solutions for four instances and for several instances objectives are found that are lower than the current state-of-the-art heuristic. Moreover, the solutions' standard deviation per instance is on average 0.07% which indicates that the results of the HGA are stable. It is concluded that the HGA gives strongly competitive results to the specialized iterated local search described by Silva et al. (2015). The average running time of the heuristic is short, on average 2.42 minutes are needed to nd the nal solution, and applicable for many real-life decisions. The increase in running time relative to Silva et al. (2015) is explained by the redundant steps made to handle inventory and a multi depot setting.
Aleman Archetti Silva HGA
Problem BKS min. min. min. avg. min. avg. time (min)
eil22 375.28 375.28 - 375.28 375.28 375.28 375.28 0.0 eil23 568.56 569.75 - 568.56 568.56 568.56 568.56 0.0 eil30 505.01 505.01 - 505.01 505.01 505.01 505.01 0.0 eil33 837.06 843.64 - 837.06 837.06 837.06 837.06 0.0 eil51 524.61 527.67 527.98 524.61 524.61 524.61 524.61 0.0 eilA76 823.89 853.20 832.71 823.89 824.92 823.89 823.89 0.1 eilB76 1009.04 1034.21 1047.17 1009.04 1012.07 1009.04 1010.18 3.1 eilC76 738.67 761.55 765.30 738.67 739.89 738.67 738.67 0.3 eilD76 687.60 695.96 705.26 687.60 689.36 686.70 687.06 1.4 eilA101 826.14 844.21 854.03 826.14 826.58 826.14 827.44 3.1 eilB101 1076.26 1112.15 1119.17 1076.26 1079.15 1076.07 1078.70 7.1 S51D1 459.50 468.79 459.50 459.50 459.50 459.50 459.50 0.0 S51D2 708.42 718.69 717.18 709.29 709.49 708.42 708.63 1.0 S51D3 947.97 969.75 960.38 948.06 950.12 947.97 947.97 0.3 S51D4 1561.34 1628.20 1569.92 1562.01 1563.29 1563.18 1565.36 1.1 S51D5 1333.67 1362.19 1339.38 1333.67 1333.67 1333.67 1334.09 1.6 S51D6 2169.10 2236.16 2182.13 2169.10 2177.78 2169.10 2177.78 2.3 S76D1 598.94 613.70 633.83 598.94 598.94 598.94 598.94 0.0 S76D2 1087.40 1128.15 1104.56 1087.40 1089.45 1087.40 1087.97 4.2 S76D3 1427.81 1472.92 1435.10 1427.86 1429.26 1425.73 1428.48 4.4 S76D4 2079.76 2180.13 2106.64 2079.76 2081.16 2082.78 2088.10 5.3 S101D1 726.59 749.19 791.11 726.59 728.45 726.59 726.59 1.0 S101D2 1378.43 1409.03 1426.15 1378.43 1386.03 1367.63 1382.38 12.2 S101D3 1874.81 1947.62 1911.40 1874.81 1880.62 1878.56 1883.70 10.9 S101D5 2790.70 2910.71 2824.16 2791.22 2795.36 2801.89 2815.15 11.4 Avg Gap to BKS +2.42 % +2.28 % +0.01 % +0.15 % -0.01% +0.13%
Avg time 1.95 min - 0.02 min 2.83 min
Table 3: Results on Belenguer et al. (2000) instances of the SDVRP with an unlimited eet size. `Aleman' refers to Aleman and Hill (2010), `Archetti' refers to Archetti et al. (2011), and `Silva' refers to Silva et al. (2015).
Aleman Archetti Silva HGA
Problem BKS min. min. min. avg. min. avg. time (min)
eil22 375.28 375.28 - 375.28 375.28 375.28 375.28 0.0 eil23 568.56 569.75 - 568.56 568.56 568.56 568.56 0.0 eil30 512.72 512.72 - 512.72 512.72 512.72 512.72 0.0 eil33 837.06 853.10 - 837.06 837.06 837.06 837.06 0.0 eil51 524.61 524.61 527.98 524.61 524.61 524.61 524.61 0.0 eilA76 823.89 851.24 829.58 823.89 825.22 823.89 823.89 0.1 eilB76 1009.04 1059.57 1023.23 1009.04 1011.19 1009.04 1010.61 3.4 eilC76 738.67 753.29 761.98 738.67 739.83 738.67 738.67 0.4 eilD76 687.60 699.35 709.51 687.60 688.37 686.70 686.89 1.7 eilA101 826.14 852.74 851.39 826.14 826.26 826.14 826.89 1.3 eilB101 1076.26 1139.27 1118.93 1076.26 1078.58 1076.26 1078.90 10.9 S51D1 459.50 471.92 464.78 459.50 459.50 459.50 459.50 0.0 S51D2 708.42 731.01 719.82 708.42 709.54 708.42 708.54 1.3 S51D3 947.97 1001.22 955.56 948.01 949.96 947.97 947.97 0.9 S51D4 1561.01 1680.66 1614.36 1561.01 1563.25 1576,18 1582.69 3.7 S51D5 1333.67 1389.40 1343.07 1333.67 1333.85 1334,92 1339.69 1.9 S51D6 2174.71 2218.23 2179.65 2169.10 2174.71 2185,14 2190.41 4.2 S76D1 598.94 606.47 634.93 598.94 598.98 598.94 598.94 0.1 S76D2 1087.99 1143.36 1129.06 1087.99 1089.69 1087.40 1088,57 5.8 S76D3 1427.81 1490.08 1442.08 1427.81 1429.01 1425.73 1429.43 6.2 S76D4 2079.76 2173.61 2100.15 2079.76 2080.76 2083.32 2083.90 6.4 S101D1 726.59 749.19 788.33 726.59 728.44 726.59 726.59 0.9 S101D2 1378.19 1443.44 1417.87 1383.35 1386.45 1367.63 1382.88 11.9 S101D3 1876.97 1988.78 1906.28 1876.97 1881.26 1885.76 1889.48 13.0 S101D5 2792.01 2984.48 2896.19 2792.01 2795.73 2883.79 2910.39 15.9 Avg Gap to BKS +3.40% +2.52 % +0.02% +0.13% +0.19% +0.36%
Avg time 2.93 min - 0.03 min 3.52 min
The second benchmark uses 33 instances for the MDVRP designed by Cordeau et al. (1997) where the number of customers vary from 50 till 417 and there are at most 6 depots. The number of vehicles per depot varies between 2 and 14 and the vehicles have a maximum capacity and driving time. For some instances an exact solution is found by Baldacci and Mingozzi (2009). These instances are marked by an asterisk, see Table 5.
The results are compared with the results obtained by Vidal et al. (2012). Their heuristic is the current-state-of-the-art heuristic and nds the BKS on all instances. The average gap between the BKS and the proposed HGA equals 0.11%. The solutions' standard deviation per instances is on average 0.05% which indicates that the HGA gives stable results for the MDVRP. The average running time of the heuristic is 27.66 minutes which is reasonable for the MDVRP. Furthermore, for the instances where the minimum is unequal to the BKS it is seen that the HGA is not converged. If the HGA runs as an overnight process the BKS is found on all instances. It is concluded that the HGA gives competitive solutions to the current state-of-the-art heuristic in reasonable time.
Vidal et al. (2012) HGA
BKS min. avg. min. avg. time (min)
*P01 576.87 576.87 576.87 576.87 576.87 0.0 *P02 473.53 473.53 473.53 473.53 473.53 0.0 P03 641.19 641.19 641.19 641.19 641.19 0.0 P04 1001.04 1001.04 1001.23 1001.04 1001.31 5.9 P05 750.03 750.03 750.03 750.03 750.03 1.8 *P06 876.50 876.50 876.50 876.50 877.28 2.2 *P07 881.97 881.97 884.43 881.97 882.56 7.6 P08 4372.78 4372.78 4397.42 4391.58 4409.12 139.4 P09 3858.66 3858.66 3868.59 3863.25 3873.39 96.8 P10 3631.11 3631.11 3636.08 3631.11 3640.63 83.1 P11 3546.06 3546.06 3548.25 3546.06 3547.60 68.9 *P12 1318.95 1318.95 1318.95 1318.95 1318.95 0.0 P13 1318.95 1318.95 1318.95 1318.95 1318.95 0.0 P14 1360.12 1360.12 1360.12 1360.12 1360.12 0.0 P15 2505.42 2505.42 2505.42 2505.42 2505.42 1.3 P16 2572.23 2572.23 2572.23 2572.23 2572.23 0.0 P17 2709.09 2709.09 2709.09 2709.09 2709.09 0.0 P18 3702.85 3702.85 3702.85 3702.85 3702.85 14.5 P19 3827.06 3827.06 3827.06 3827.06 3827.06 1.6 P20 4058.07 4058.07 4058.07 4058.07 4058.07 0.9 P21 5474.84 5474.84 5476.41 5474.84 5476.01 94.6 P22 5702.16 5702.16 5702.16 5702.16 5702.16 19.8 P23 6078.75 6078.75 6078.75 6078.75 6078.75 13.9 pr01 861.32 861.32 861.32 861.32 861.32 0.0 pr02 1307.34 1307.34 1307.34 1307.34 1307.34 1.2 pr03 1803.80 1803.80 1803.80 1803.80 1803.80 8.9 pr04 2058.31 2058.31 2059.36 2058.31 2058.52 7.4 pr05 2331.20 2331.20 2340.29 2337.84 2343.78 78.8 pr06 2676.30 2676.30 2681.93 2677.64 2684.68 107.2 pr07 1089.56 1089.56 1089.56 1089.56 1089.56 0.0 pr08 1664.85 1664.85 1665.05 1664.85 1664.85 1.1 pr09 2133.20 2133.20 2134.17 2133.20 2134.29 36.7 pr10 2868.26 2868.26 2886.59 2878.63 2899.01 118.9 Avg Gap to BKS +0.00 % +0.08 % +0.04 % +0.11 %
Avg time 4.42 min 27.66 min
5.4 Results MDSDVRP-I Instances
This section shows the performance of the HGA for the MDSDVRP-I instances. We will examine the impact of inventory levels at the depots, the in uence of the allocation skewness of the inventory and the relevance of the possibility for split deliveries. For this analysis, we created a new set of instances as explained in Section 5.1. However, to get useful results from the ILP model, the instances are smaller, ηc∈ {5, 10, . . . , 40}. For each setting given in Table 1 there are three random instances; thus, there are
nine instances for each combination of ηc and ηd. The inventory in an instance might be reallocated
without changing the coordinates of the customers and depots. We expect that a tighter inventory will consistently lead to higher objective values whereas the skewness may result in higher or lower objective values, depending on the allocation of customers.
The rst analysis considers the in uence of inventory levels. The instances are solved for two dierent inventory tightness values, namely ψ = 1.25 and ψ = 1.05. For ψ = 1.05, there is just enough inventory to have feasible solutions. The results of the ILP model are shown in Table 6. It is seen that the ILP e ciently solves the small instances, whereas medium instances, i.e., 25 customers and more, are not solved within ve days. Instances in which the optimal solution contains split deliveries were signi cantly more time-consuming than instances without a split delivery (p-value: 0.000).
For the HGA the same parameter settings as found in the calibration section are used except for the number of iterations, which is set to 10.000 since the instances are smaller than those used for the calibration. The results are given as an average result of ten runs, see Table 6. Already for the smallest instances the HGA is faster than the ILP, furthermore the runtime of the HGA stays consistently low for most instance sets. For instances with 20 customers or more, the HGA outperforms the ILP, both in nding feasible solutions and solution quality. Moreover, the HGA nds the optimal solution for all instances solved by the ILP and for no instance the ILP gives a better solution than the HGA. It is concluded that the HGA gives high quality solutions for the MDSDVRP-I.
Table 6 shows that tight inventory levels result in higher objective values and longer runtimes for the HGA. Furthermore, it is seen that large instances are less aected by tight inventory levels than small instances, which is explained by the number of split deliveries. For ψ = 1.25, there are on average 0.37 split deliveries for the rst 100 instances, while for ψ = 1.05, there are on average 0.90 split deliveries. For large instances, these numbers equal 0.07 and 0.16, respectively. It follows that large instances have more possibilities for avoiding split deliveries and thus are less disadvantaged by tight inventory levels. Figure 6a shows the average objective increase compared to no inventory constraint for the instances in sets 4, 12 and 20 for dierent values for the inventory tightness. These instance sets have 5, 15 and 25 customers, respectively. First, notice that tight inventory levels mainly aects small instances, in which routing costs can increase by almost 150%. Second, the graph shows some stair-wise steps that are caused by eliminating split deliveries. Third, small instances need higher inventory levels before the consequences of the inventory constraints fade out. It is shown that the in uence of inventory to routing should not be underrated since extra inventory can result in signi cant savings on the routing part. Furthermore, an underestimate of demand is especially costly for rms with only a few delivery locations.
The second analysis considers the eect of split deliveries. For small instances, i.e., less than 25 customers, split deliveries were necessary to obtain feasible solutions. Instances with 25 customers or more are feasible without the possibility for split deliveries, although the routing cost can be decreased up to 10% by introducing split deliveries. Furthermore, savings due to split deliveries become larger if the inventory levels of the depots decrease. For some instances, the savings could be enormous: Suppose there are three depots and three customers. The rst two depots, called A and B, are located close to the customers while the third one is located far away. Depot A and B together have enough inventory to satisfy all customers, but no depot has enough inventory to satisfy two of the customers. When split deliveries are allowed, the demand of all customers can be satis ed by depot A and B, whereas we need depot C, which is located far away, if split deliveries are not allowed.
1.0 1.5 2.0 2.5 3.0 0 50 100 150 Inventory tightness (ψ) Cost increase (%) 5 customers 15 customers 25 customers
(a) Relative increase in routing costs compared to a scenario with unlimited inventory.
0.5 1.0 1.5 2.0 0 20 40 60 80 100 Inventory skewness (θ) Cost increase (%) 5 customers 15 customers 25 customers
(b) Relative increase in routing costs compared to the inventory allocation with the lowest routing costs Figure 6: In uence of inventory tightness ψ, and inventory skewness θ.
there is one central depot, say a distribution center (DC), while the other depots are smaller depots, say stores. The same instances are used, and the depot closest to the (0,0) coordinate becomes the DC and has an inventory level of (ψ − θ) · Ψ, with 0 ≤ θ ≤ ψ, whereas the inventory of the stores equals θ · Ψ
ηd−1.
Thus, ψ indicates the inventory tightness, and θ gives the skewness of the inventory. As can be seen, the total inventory among the depots is independent of θ.
1 5 1 262.09 262.09 0.00 262.09 0.00 9/9 0.00 0.00 262.09 262.09 0.00 262.09 0.00 9/9 0.00 0.00 0.00 2 5 2 223.42 223.42 0.00 223.42 0.00 9/9 0.00 0.00 241.37 241.37 0.00 241.37 0.00 9/9 0.00 0.00 8.03 3 5 4 266.25 266.25 0.00 266.25 0.00 9/9 0.00 0.00 336.76 336.76 0.00 336.76 0.00 9/9 0.00 0.00 26.48 4 5 6 525.87 525.87 0.00 525.87 0.00 9/9 0.00 0.00 667.73 667.73 0.01 667.73 0.00 9/9 0.00 0.00 26.98 5 10 1 331.43 331.43 0.01 331.43 0.00 9/9 0.00 0.00 331.43 331.43 0.00 331.43 0.00 9/9 0.00 0.00 0.00 6 10 2 401.29 401.29 0.08 401.29 0.00 9/9 0.00 0.00 426.90 426.90 0.04 426.90 0.00 9/9 0.00 0.00 6.38 7 10 4 557.35 590.00 24.85 590.00 0.00 7/7 5.86 0.00 596.44 653.47 514.28 653.43 0.00 6/6 9.56 -0.00 10.75 8 10 6 485.50 491.84 486.89 491.84 0.00 8/8 1.31 0.00 537.95 570.10 659.98 570.10 0.00 6/6 5.98 0.00 15.91 9 15 1 426.36 426.36 0.13 426.36 0.00 9/9 0.00 0.00 426.36 426.36 0.09 426.36 0.00 9/9 0.00 0.00 0.00 10 15 2 382.20 382.20 0.09 382.20 0.00 9/9 0.00 0.00 384.12 384.12 0.12 384.12 0.00 9/9 0.00 0.00 0.50 11 15 4 528.40 545.93 3659.40 545.42 0.00 7/7 3.32 0.00 542.01 574.26 1150.50 572.13 0.00 6/6 5.56 -0.37 4.90 12 15 6 637.18 672.74 2829.13 672.74 0.00 4/4 5.58 0.00 666.94 765.35 1653.58 765.13 0.00 3/3 14.72 -0.03 13.73 13 20 1 623.02 623.02 122.43 623.02 0.00 9/9 0.00 0.00 623.02 623.02 122.80 623.02 0.00 9/9 0.00 0.00 0.00 14 20 2 548.88 ∞ 2521.22 567.95 0.00 7/7 3.47 -∞ 559.00 574.98 723.17 574.98 0.00 7/7 2.86 0.00 1.24 15 20 4 538.76 546.80 715.73 546.80 0.02 7/7 1.49 0.00 553.22 ∞ 2450.75 596.32 0.00 3/3 7.79 -∞ 9.06 16 20 6 542.38 714.42 1545.23 701.57 0.14 0/0 29.35 -1.80 556.36 ∞ 3923.84 773.85 0.18 0/0 39.09 -∞ 10.30 17 25 1 572.77 ∞ 348.49 581.81 0.00 6/6 1.58 -∞ 547.38 604.13 147.59 581.81 0.00 6/6 6.29 -3.69 0.00 18 25 2 539.84 ∞ 111.98 575.80 0.00 6/6 6.66 -∞ - - - 588.50 0.00 - - - 2.21 19 25 4 595.52 ∞ 3496.83 643.18 0.00 0/0 8.00 -∞ - - - 666.35 0.00 - - - 3.60 20 25 6 - - - 762.11 1.31 - - - 836.62 2.69 - - - 9.78 21 30 1 - - - 660.86 0.01 - - - 660.86 0.39 - - - 0.00 22 30 2 - - - 659.48 0.01 - - - 667.44 0.00 - - - 1.21 23 30 4 - - - 749.82 0.00 - - - 784.06 0.23 - - - 4.57 24 30 6 - - - 805.97 4.75 - - - 877.22 4.20 - - - 8.84 25 35 1 - - - 738.31 0.00 - - - 738.31 0.00 - - - 0.00 26 35 2 - - - 777.43 0.00 - - - 797.09 0.02 - - - 2.26 27 35 4 - - - 812.33 0.00 - - - 836.40 0.03 - - - 2.96 28 35 6 - - - 952.00 2.21 - - - 1015.96 2.56 - - - 6.72 29 40 1 651.00 ∞ 6691.17 762.56 0.00 0/0 17.14 -∞ - - - 762.56 0.00 - - - 0.00 30 40 2 - - - 894.69 0.04 - - - 899.67 0.04 - - - 0.56 31 40 4 - - - 804.30 0.00 - - - 816.60 0.03 - - - 1.53 32 40 6 - - - 917.92 0.07 - - - 975.65 5.85 - - - 6.29
Table 6: Performance on the MDSDVRP-I Instances. Each instance set contains nine random instances and each number is the average of ten runs, except for the results of the ILP. Column ηc indicates the number of customers and column ηd gives the number of depots. LB and UB show the average lower bound and upper bound found by the
6 Conclusion
We introduced inventory constraints for the MDVRP that can have a signi cant impact to the optimal solution. In this problem, vehicles are stationed at several depots, a deterministic set of customers has to be satis ed and each depot has a known inventory level. The aim is to minimise the distance travelled by the vehicles while each depot obeys its inventory level, i.e., the maximum demand served from the depot. Two special cases of the introduced problem are the MDVRP and the SDVRP.
Suppose there is a distribution center and several stores with a local inventory. A rm can decrease its routing costs by using the local inventory at the stores, however, it is not obvious that these stores will always have enough inventory to satisfy all demand. In a setting like this, it is not realistic to assume that each depot can satisfy all demand. To solve the problem one should take into account the limited inventory of the stores, which is achieved by adding inventory constraints.
The problem is formulated as an integer linear program and some cuts are provided to improve the performance. Since the problem is NP-hard, a hybrid genetic algorithm is proposed that is capable of solving a large class of problems including the two special cases. The heuristic is able to handle dierent constraints like a maximum driving time, a vehicle capacity and a serving time per customer.
Using benchmark instances from the literature shows that the heuristic gives highly competitive solutions for the SDVRP in a short time interval. For 8 out of 48 instances, a new best known solution is found and many other best known solutions were equaled. For the MDVRP the heuristic nds competitive solutions in reasonable time, for most instances the best known solution is found. New test instances are created for the MDSDVRP-I, and instances up to 25 customers are solved by the ILP. For all instances the HGA nds the same or better solutions than the ILP in a short running time.
Some characteristics of the problem we introduced are described. First, inventory tightness among the depots can have a signi cant in uence on the routing costs. For small instances, i.e., 15 customers, routing costs can increase by more than 50% compared to no inventory constraints. For larger instances, this eect becomes smaller but is still present. Second, it is seen that the geographically most central depot should have slightly more inventory than the other depots. The third characteristic states that the need for split deliveries increases due to inventory constraints. For small instances, split deliveries were necessary to have feasible solutions, while for larger instances, savings up to 10% can be obtained.
7 Further Research
There is ample of scope for further research. A straight forward extension is to add time windows or a setting with multiple days. A conceptual newer extension is to add the possibility to reallocate inventory among depots for lower transportation costs, the lower costs could be explained by larger vehicles. The properties that indicate whether a split might be fruitful described by Archetti et al. (2008) might be used to improve the local search operators for split deliveries. Furthermore, the benchmark of Belenguer et al. (2000) shows that the HGA performs better for the unlimited eet variant than for the limited eet variant while for many instances the unlimited eet solution is also feasible for the limited eet variant. It might help to allow infeasible solutions with respect to the number of vehicles.
The current distance measure is designed for a multi depot setting, as a result the dierence between any solutions of a single depot instance is always zero. The solution quality might improve when a dierent distance measure is used for the single depot setting or the runtime will decrease by approximately 25% by removing the distance measure for this case. Furthermore, for instances that allow split deliveries almost 15% of the time is spend to the splitReinsertion operator. One may test whether the solution quality is aected when this operator is applied less often to gain a performance bene t.
Aleman, R.E. and R Hill (2010). A tabu search with vocabulary building approach for the vehicle routing problem with split demands. International Journal of Metaheuristics 1 (1), 5580.
Aleman, R.E., X Zhang, and R Hill (2010). An adaptive memory algorithm for the split delivery vehicle routing problem. Journal of Heuristics 16 (3), 441473.
Archetti, C., N Bianchessie, and M Speranza (2014). Branch-and-cut algorithms for the split delivery vehicle routing problem. European Journal of Operational Research 238 (3), 685698.
Archetti, C., M Savelsbergh, and M Speranza (2008). To split or not to split: That is the question. Transportation Research Part E: Logistics and Transportation Review 44 (1), 114123.
Archetti, C., G Speranza, and N Bianchessie (2011). A column generation approach for the split delivery vehicle routing problem. Networks 58 (4), 241254.
Archetti, C. and M Speranza (2006). A tabu search algorithm for the split delivery vehicle routing problem. Transportation Science 40 (1), 6473.
Archetti, C. and M Speranza (2008). The split delivery vehicle routing problem: A survey. In The vehicle routing problem: Latest advances and new challenges, pp. 103122. Springer.
Archetti, C., M Speranza, and M Savelsbergh (2008). An optimization-based heuristic for the split delivery vehicle routing problem. Transportation Science 42 (1), 2231.
Baldacci, R. and A Mingozzi (2009). A uni ed exact method for solving dierent classes of vehicle routing problems. Mathematical Programming 120 (2), 347380.
Belenguer, J. M., M. C Martinez, and E Mota (2000). A lower bound for the split delivery vehicle routing problem. Operational Research 48 (5), 801810.
Boudia, M., C Prins, and M Reghioui (2007). An eective memetic algorithm with population man-agement for the split delivery vehicle routing problem. In Thomas Bartz-Beielstein, María José Blesa Aguilera, Christian Blum, Boris Naujoks, Andrea Roli, Günter Rudolph, and Michael Sampels (Eds.), Hybrid Metaheuristics, pp. 1630. Springer.
Cordeau, J.F., M Gendreau, and G Laporte (1997). A tabu search heuristic for periodic and multi-depot vehicle routing problems. Networks 30 (2), 105119.
Dantzig, G.B. and J Ramser (1959). The truck dispatching problem. Management science 6 (1), 8091. Desaulniers, G. (2010). Branch-and-price-and-cut for the split delivery vehicle routing with time windows.
Operations Research 58 (1), 179192.
Dror, M., G Laporte, and P Trudeau (1994). Vehicle routing with split deliveries. Discrete Applied Mathematics 50 (3), 239254.
Dror, M. and P Trudeau (1989). Savings by split delivery routing. Transportation Science 23 (2), 141145. Dror, M. and P Trudeau (1990). Split delivery routing. Naval Research Logistics 37 (3), 383402. Gendreau, M., P Dejax, D Feillet, and C Gueguen (2006). Vehicle routing with time windows and split
deliveries. Technical Paper 851.
Glover, F. and J.-K Hao (2011). The case for strategic oscillation. Annals of Operations Research 183 (1), 163173.
Gulczynski, D.J., B Golden, and E Wasil (2008). Recent developments in modeling and solving the split delivery vehicle routing problem. Tutorials in Operations Research, 170180.
