Combined location and routing problems for drug distribution

Discrete Applied Mathematics 165 (2014) 130–145

Contents lists available at ScienceDirect

Discrete Applied Mathematics journal homepage: www.elsevier.com/locate/dam

Combined location and routing problems for drug distribution Alberto Ceselli, Giovanni Righini, Emanuele Tresoldi ∗ Dipartimento di Tecnologie dell’Informazione, Università degli Studi di Milano, Italy

article

info

Article history: Received 16 December 2011 Received in revised form 5 June 2013 Accepted 22 July 2013 Available online 12 August 2013 Keywords: Location-routing Column generation Branch-and-price

abstract We present a model for the optimization of logistics operations in emergency health care systems; in particular, we study the problem of efficient distribution of vaccines or drugs through the simultaneous and coordinated use of distribution centers and vehicles. We devise an exact algorithm based on column generation with three different types of columns and branch-and-bound. The pricing subproblems are solved through advanced dynamic programming techniques. In order to strengthen the dual bounds, we adapt two families of cuts from the literature and we introduce a new one. Our framework also includes primal heuristics and ad-hoc branching rules. An experimental campaign on realistic data proves our method to be effective and flexible. © 2013 Elsevier B.V. All rights reserved.

1. Introduction Mathematical programming models and algorithms have been successfully used for decades to optimize operations in distribution logistics: typical examples concern freight carriers, mail services and on-demand pick-up and delivery services. A more recent field of investigation concerns the application of operations research techniques to the optimization of logistics operations in health care systems and emergency management. These sectors are characterized by a larger dependency on human behavior (which is often unpredictable), the need for fairness in service provision (which is not an issue in industrial logistics) and the lack of reliable historical data (because of the uniqueness of the events considered, especially in the case of emergencies). In particular in this paper we discuss the application of mathematical programming techniques to the distribution of drugs or vaccines in the case of an emergency. Such situations are hardly foreseeable and almost impossible to avoid; hence it is crucial to be able to set up an efficient response system to provide the required services in an optimized way when needed. Optimization in this case is such an essential component of a good response system that the U.S. Centers for Diseases Control and Prevention in 2007 evaluated ‘‘the timeliness of distributing antibiotics to the general public as an effective measure against a release of anthrax’’ [6]. Operations research techniques have already been applied to emergency drug distribution: in [20] the authors present a stochastic model for a Vehicle Routing Problem (VRP) arising in large-scale bio-terrorism emergency, in which the maximum number of citizens must be reached within a specified time limit using a heterogeneous fleet of vehicles. Such a model is then reformulated as a deterministic VRP and solved with a tabu search algorithm. In this paper we expand the model introduced in [20], also including an additional distribution strategy based on distribution centers. This allows us to consider more realistic situations in which multiple distribution channels can be exploited at the same time. In particular, we explore the option of reaching citizens in two ways: either by delivering drugs at their homes with a heterogeneous fleet of vehicles, or by establishing distribution centers where the citizens go by their own means to receive treatments or drugs. Hence a combined Location and Routing Problem (LRP) arises. Location-routing

∗

Corresponding author. Tel.: +39 036360562; fax: +39 363360803. E-mail addresses: [email protected] (A. Ceselli), [email protected] (G. Righini), [email protected] (E. Tresoldi).

0166-218X/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.dam.2013.07.016

A. Ceselli et al. / Discrete Applied Mathematics 165 (2014) 130–145

131

is a lively research area [18]; in particular, recent exact algorithms like [3,4,8] try to cope with refined handling of limited resources. Problems with up to 200 customers and 14 potential depots for the so-called Capacitated LRP can be solved to proven optimality, even if several hours of computation might be needed. Our problem shares some features also with the Multi Depot Heterogeneous VRP [5], in which the starting point of each vehicle is part of the decision process, and for which state-of-the-art exact algorithms can solve problems with hundred customers in a few hours. Finally, for what concerns the objective of the optimization, our problem is similar to the Team Orienteering Problem (TOP) [22], being concerned in reaching the largest amount of demand with limited resources: TOP instances involving up to one hundred customers can be tackled by state-of-the-art algorithms. To the best of our knowledge our double channel distribution problem has never been addressed before. We present an exact algorithm, which is based on dynamic column and cut generation and branch-and-bound, where the pricing subproblems are solved through advanced dynamic programming techniques. The paper is organized as follows: in Section 2 we describe the problem and we point out its main features. A mathematical model is given in Section 3 and our solution approach is detailed in Section 4. In Sections 5 and 6 the results of our experimental campaign and some conclusions are reported. 2. Problem description The problem we address fits into the following scenario: the population living in a given region has been partitioned into small groups and a corresponding suitable set of delivery sites, one for each group, has been identified. These sites represent the demand points, i.e. the destinations where the drugs have to be brought. We assume that once a delivery site is visited, all people assigned to it will get the drugs within a very short time. A limited number of the delivery sites can be selected as facilities where the drugs can be stored; these selected facilities are the origins of the distribution. When a facility is selected it can employ two different strategies to deliver the drugs to the sites assigned to it. Distribution center strategy. Facilities using this strategy are called Distribution Centers (DCs). They are characterized by a capacity and a distribution range (Fig. 1(a)). The capacity is an estimate of the overall demand that the facility can satisfy. The range defines the maximum distance within which a DC can provide its services to the delivery sites assigned to it. Several distribution methods can lie behind this strategy, as the use of many small vehicles, one for every site, self-service points, or the use of the existing postal network for last-mile services. We remark that when the capacity of a DC is less than the overall demand falling within its range, we assume that it is possible to decide which delivery sites are assigned to the DC and which of them are not. Routing strategy. This strategy, illustrated in Fig. 1(b), is the classical VRP-like distribution method. When a facility uses this strategy it is called depot. With every depot is associated a fleet of heterogeneous vehicles with different capacities. Every vehicle can perform a single route starting from the depot, visiting some delivery sites and coming back to the same depot. Decision variables. Fig. 1(c) shows the big picture of the logistics system with mixed distribution strategy, where the problem is to choose which facilities are to be selected and what is the most appropriate distribution strategy for each of them. Optimal routes must also be computed for the vehicles. Constraints. The following constraints must be respected. I No delivery site can be served multiple times by a single DC or visited more than once by a single route. Multiple assignments to different distribution centers and multiple visits by different routes are allowed but they do not increase the number of citizens served. II No split delivery: when a delivery site is assigned to a DC or visited by a vehicle, its demand is entirely served. III If a facility is used it can be either a depot or a DC, but not both at the same time. IV The total number of DCs that can be opened is limited by a budget constraint. V The overall demand of the delivery sites assigned to a DC must not exceed the capacity of the DC. VI All delivery sites assigned to a DC must be within its range. VII The number of available vehicles of each type is limited. VIII Each route starts from a depot. IX The overall demand of the sites visited by a vehicle must not exceed the capacity of the vehicle. X Every delivery site in a route must be visited before the deadline. The deadline represents a time limit within which the pharmacological treatments are most effective. It is worth noting that is not necessary for a vehicle to be back at the depot within the time limit. Objective function. We want to maximize the demand served within the deadline. We are aware that many other objective functions can and should be considered in the design of optimized emergency systems, e.g. equity in service provision. However we decided to concentrate our analysis on this objective for two main reasons: first of all, in the case of a pandemic outbreak it is of paramount importance to reach the largest part of the population as quickly as possible; this carries an advantage also to those who do not receive any treatment because it limits the spread of the illness. The second reason is that this is the objective function considered by Shen et al. [20] and therefore it allows for a comparison of the two algorithmic approaches.

132

A. Ceselli et al. / Discrete Applied Mathematics 165 (2014) 130–145

(a) Distribution Center Strategy: the dotted circle represents the range of the DC; P is the capacity of the DC; d is the demand of the sites. Black sites are those assigned to the DC.

(b) Routing Strategy: a depot and three vehicles: f is the capacity of the vehicle; d is the demand of the sites. Black sites are those visited by the vehicles.

(c) Mixed Strategy: combination of DC-based and vehicle-based strategies. Fig. 1. Distribution logistics strategies considered in the GLDP.

The resulting problem does not fit into the description of any classical LRP, as it requires us to select not only which locations have to be used but also, for every location opened, which distribution strategy has to be employed. Since this problem combines location decisions with different distribution methods we call it Generalized Location and Distribution Problem (GLDP). 3. Mathematical formulation The GLDP can be formally described as follows. A complete graph G = (N , A) is given where N is a set of nodes representing the delivery sites and A is a set of arcs representing the connections between them. Two values di and ti are associated with each delivery site i ∈ N: di represents the demand and ti the service time. For each arc (i, j) ∈ A a traveling time cij is given; we assume that all vehicles travel at the same speed. A subset L ⊆ N is given; it represents potential locations for facilities: in each site l ∈ L either a DC or a depot (or neither) can be located. For each DC located in l ∈ L, a service capacity pl and a delivery range rl are given. A cluster for a DC located in l is a set of nodes N¯ ⊆ N such that cli + ti ≤ rl



∀i ∈ N¯

(1)

di ≤ pl .

(2)

i∈N¯



For every l ∈ L, let Wl be the set of all clusters associated with a DC located in l, and let W = l∈L Wl . Due to budget limitations, at most B DCs can be selected. Let H be the set of vehicle types, and for each h ∈ H let Vh and fh be respectively the number of vehicles available for  type h and the capacity of each of them; the total number of vehicles is T = V . A route for a vehicle of type h ∈ H h h∈H

A. Ceselli et al. / Discrete Applied Mathematics 165 (2014) 130–145

133

associated with a depot l ∈ L is a sequence K¯ of sites, starting from l, such that HPK¯ +



ti ≤ M

(3)

i∈K¯



di ≤ fh

(4)

i∈K¯

where HPK¯ is the time required to leave the depot and visit all sites in K¯ in sequence, omitting the time needed to go back to the depot and M is the maximum allowed time to perform the deliveries. Inequalities (3) and (4) model Conditions X and IX, imposing respectively that all sites in the route are reached within the deadline M, and that their total demand does not exceed the vehicle capacity. Hereafter we present a Set Covering formulation of the GLDP, using three sets of binary variables. For each node i ∈ N a variable si represents whether the node is skipped: it takes value 1 if node i is not served and 0 otherwise. A binary variable yw is associated with each cluster w ∈ Wl : it takes value 1 if a DC located in l ∈ L is used to serve cluster w ∈ Wl , and 0 otherwise. Finally a binary variable zk for each k ∈ Klh represents a route operated by a vehicle of type h ∈ H , starting from a depot located at l ∈ L; it takes value 1 if the route is selected and 0 otherwise. Clusters and routes are described by binary coefficients aiw and bik . The former, aiw , takes value 1 if site i ∈ N belongs to cluster w ∈ Wl and 0 otherwise; the latter, bik , takes value 1 if site i ∈ N belongs to route k ∈ Klh and 0 otherwise. The objective (5) of maximizing the demand satisfied is equivalently expressed as the minimization of the unsatisfied demand. The Set Covering formulation of the GLDP is the following: min



d i si

(5)

i∈N

s.t. si +



aiw yw +

 

l∈L w∈Wl



zk ≤ Vh

bik zk ≥ 1 ∀i ∈ N

(6)

l∈L h∈H k∈Klh

∀h ∈ H

(7)

l∈L k∈Klh



zk +

Vh yw ≤ Vh

∀l ∈ L, ∀h ∈ H

(8)

w∈Wl

k∈Klh





yw ≤ B

(9)

l∈L w∈Wl

si ∈ {0, 1}

(10)

yw ∈ {0, 1}

(11)

∀i ∈ N ∀l ∈ L, ∀h ∈ H , ∀k ∈ Klh zk ∈ {0, 1} ∀l ∈ L, ∀w ∈ Wl .

(12)

Constraints (6) state that either a site is served by at least one DC or route, or its associated skip variable must assume value 1. It is worth noting that although these constraints allow for multiple visits of the same site, visiting a site more than once does not improve the value of the objective function. In addition a site cannot be visited more than once in the same cluster or in the same route. These constraints together with the definitions of sets Wl and Klh correspond to Conditions I and II. Constraints (7) ensure Condition VII to be satisfied: an upper bound equal to Vh , for all h ∈ H , is set on the number of vehicles of type h that can be used. Constraints (8) ensure that if a node is used as a depot it cannot be used as a DC and vice-versa, according to Condition III. It is worth noting that no location is forced to be used. The last set of constraints (9), corresponding to Condition IV, put an upper limit on the overall number of DCs. Our approach based on model (5)–(12) is highly flexible: as described in [7], it is possible to include even more features without significantly affecting the corresponding solution algorithms. 4. An exact algorithm We propose an algorithm for the exact solution of the GLDP that is based on the branch-and-cut-and-price paradigm [16]. This method starts from the linear relaxation of the problem, usually referred to as Master Problem (MP), obtained by substituting constrains (10)–(12) with 0 ≤ si ≤ 1 ∀ i ∈ N 0 ≤ yw ≤ 1 ∀l ∈ L, ∀w ∈ Wl . 0 ≤ zk ≤ 1 ∀l ∈ L, ∀h ∈ H , ∀k ∈ Klh . The derived model has an exponential number of columns and so, in order to obtain valid lower bounds, it is solved via column generation (see Desaulniers et al. [11]). When the MP is solved to optimality we look for additional inequalities that

134

A. Ceselli et al. / Discrete Applied Mathematics 165 (2014) 130–145

are violated by the current LP solution. This process of column and row generation is iterated until neither further useful variables nor violated cuts are found. If the solution achieved is fractional, and its value is lower than a known upper bound, then branching is performed. We explore an enumeration tree, recomputing the lower bound at every node. In this section we present the main components of our branch-and-cut-and-price algorithm namely: column generation, cut generation, branching and primal heuristics. 4.1. Column generation The MP contains a number of variables (y and z) whose number is exponential in the cardinality of the sites’ set N , and so a column generation method is applied to optimize it. Such an approach starts from a Restricted MP (RMP) made up of only a small subset of variables. In particular, in our case, the initial RMP is composed of (a) all columns corresponding to skip variables si ; there is a polynomial number of them, that is |N |; (b) a subset of W containing |L| · |N | columns that represents, for every possible DC location, a cluster with only one site; (c) a subset of K containing all |L| · |H | · |N | columns representing, for every possible depot location and every type of vehicle, routes serving only one site. Then, the RMP is solved, and a search for new variables having negative reduced cost is performed. If any such variable is found, it is inserted in the RMP and the process is iterated; otherwise, the RMP solution is optimal for MP as well. Let λ ≥ 0, µ ≤ 0, σ ≤ 0 and ρ ≤ 0 be vectors of dual variables corresponding to constraints (6)–(9). The reduced costs, ζw and ξk , associated with primal variables yw and zk respectively can then be computed, by means of the vectors of dual variables, in the following way:

ζw = −



aiw λi −



Vh σlh − ρ

(13)

h∈H

i∈N

ξk = −



bik λi − µh − σlh .

(14)

i∈N

Two different pricing subproblems have to be solved, respectively for yw and zk variables. It is worth noting that, since these variables represent respectively clusters and routes, each solution to their associated pricing subproblems has to satisfy respectively Conditions VI and V, and Conditions X and IX. Distribution Center pricing subproblem. For each l ∈ L the pricing subproblem associated with variables yw with w ∈ Wl is min ζw = −



ai w λ i −

Vh σlh − ρ

(15)

h∈H

i∈N

s.t. (cli + ti )aiw ≤ rl





∀i ∈ N

d i ai w ≤ p l .

(16) (17)

i∈Rw

aiw ∈ {0, 1}

∀i ∈ N .

(18)

The set of coefficients aiw , encoding a cluster of sites w ∈ Wl served by a DC located in l, has to be determined. Constraints (16) simply require us to fix some aiw value to 0. Since parameters rl , pl and dual variables σlh depend on the facility location where the DC is located, DC pricing can be performed by solving |L| different instances of the 0–1 Knapsack Problem (KP) [17], one for each l ∈ L. Our solution algorithm is based on a classical pseudo-polynomial time algorithm [9]. First, for each l ∈ L, we define an order on the elements of the set N (rl ) = {i ∈ N : cli + ti ≤ rl }. Then, a set of labels (i, u, S , v) is computed, representing the minimum reduced cost v that can be obtained using at most u units of capacity and considering only the first i elements in the set N (rl ), and the set of elements S yielding this reduced cost. The computation of all these labels is done by a dynamic programming algorithm. An optimal solution to the DC pricing problem can then be found in the set {(|N (rl )|, u, S , v)}. It is worth noting that this algorithm allows us to find several solutions in a single run and then it is useful for multiple pricing. Moreover, it allows us to take into account additional features arising after the introduction of strengthening cuts, as explained in the subsequent sections. Route pricing subproblem. Given a vehicle h ∈ H and a potential location l ∈ L for a depot, the pricing subproblem associated with variables zk reads min ξk = −



bik λi − µh − σlh

i∈N

s.t. HPk (b) +

 i∈N

bik ti ≤ M

(19)

A. Ceselli et al. / Discrete Applied Mathematics 165 (2014) 130–145



bik di ≤ Vh .

135

(20)

i∈N

bik ∈ {0, 1}

∀i ∈ N .

(21)

The task is to find the set of values bik defining the set of sites visited by the route k ∈ Klh , together with a minimum cost Hamiltonian Path HP (b) visiting them. Since the values µh and σlh depend on h ∈ H and l ∈ L, in order to optimize this pricing subproblem we need to solve |H | · |L| instances of the Resource Constrained Elementary Shortest Path Problem (RCESPP) [12], one for each h ∈ H and l ∈ L. The literature on solution methods for RCESPP is growing fast: besides k-cycle elimination [13], recent approaches include the ng-routes’ relaxation [2] and the introduction of strong degree constraints [8]. We implemented a procedure similar to that described in [5], that proved to be effective on routing problems with multiple depots and heterogeneous fleets: it includes three different steps: a greedy algorithm, a heuristic dynamic programming algorithm and an exact dynamic programming algorithm. In the remainder we sketch them following the reverse order with respect to their execution, and we refer to [5] for the details of each method. The exact dynamic programming algorithm exploits bidirectional search and decremental state space relaxation [19]. Furthermore, instead of solving an RCESPP for each location l ∈ L of the depot and for each kind h ∈ H of the vehicle employed, we exploit a technique called aggregated pricing [5], that allows us to solve in a single run all the RCESPPs related to the same vehicle type, thus reducing the total number of calls to the dynamic programming algorithm to |H |. Solving the RCESPP by means of the exact dynamic programming algorithm can be too time consuming and thus, following [5], we developed a faster heuristic dynamic programming algorithm. It has the same structure of the exact one, but dominance conditions are relaxed, allowing us to fathom more labels and reducing consequently the computational time, at the price of losing optimality guarantees. Since also the heuristic dynamic programming pricing can be time consuming, at first we search for negative reduced cost routes with a greedy algorithm, that follows a nearest neighbor policy. For each depot l ∈ L and vehicle h ∈ H , a partial path including only l is considered and iteratively extended to a site which has not been visited yet, and yields maximum improvement in the path reduced cost, without violating capacity and distance constraints on vehicle h. When no such site exists, the partial path is closed. 4.2. Additional inequalities In order to strengthen the bound we interleave the column generation process with the dynamic generation of violated inequalities. Many families of valid inequalities have been proposed for the VRP, but not all of them easily extend to the heterogeneous fleet case, that is in turn an essential feature of our application. We modified and extended subset-row [14] and 2-path [15] inequalities. Moreover, we propose a new family of inequalities that are specific to GLDP. 4.2.1. Subset-row inequalities These cuts represent an adaptation of those proposed for the VRP by Jepsen et al. [14], limited to triplets of sites. In every integer feasible solution for every set of three sites there must be at most a single route visiting at least two of them. Since in our problem the sites can be served either by DCs or by routes, in every integer feasible solution for every set of three sites there must be at most either a single DC or a single route serving at least two of them. Subset-row inequalities can then be formalized as follows. Let C = {C ⊆ N : |C | = 3} be the set made of all possible subsets of three sites. For C ∈ C let W (C ) ⊆ W be the subset of all clusters, regardless of the location of their associated DC, containing at least two of the three sites in C ; let K (C ) ⊆ K be the set of all routes, regardless of their starting depot and kind of vehicle used, that visit at least two of the three sites in C . Then the following inequalities are valid for the GLDP:

 w∈W (C )

yw +



zk ≤ 1 ∀C ∈ C .

(22)

k∈K (C )

The separation of these inequalities can be done by complete enumeration since their number is polynomial in the cardinality of the set of sites (O(|N |3 )). However to speed up the separation process we take into account only a subset N¯ ⊆ N of the sites for the generation of the clusters’ set C . In our algorithm a site i ∈ N belongs to N¯ if the value of its associated skip variable si is strictly lower than 1 in the current LP iteration. The introduction of these inequalities changes the pricing subproblems and therefore it requires us to modify the pricing algorithms to take into account the values of the corresponding dual variables. Let ψC be the (non-positive) dual variable associated with the inequality corresponding to the triplet C . The definition of the labels used in both the DC pricer and the dynamic programming procedures for the RCESPP needs to be modified, adding a new resource qC for the cut of every additional triplet C . Handling these additional resources in the route pricing algorithms is done as described in [5]; therefore, in the remainder we only focus on how to handle them in the DC pricing algorithm.

136

A. Ceselli et al. / Discrete Applied Mathematics 165 (2014) 130–145

The labels used in the DC pricing algorithm for each l ∈ L are (i, u, S , q, v), where q is a vector with one component qC for each triplet C whose subset-row inequality is in the RMP. Each qC represents the number of elements of C included in S , and so 0 ≤ qC ≤ 3. Our label pushing algorithm is modified as follows. First, we create an initial label (0, 0, ∅, (0 . . . 0), 0). Then, we iteratively extend labels by increasing values of i, and for each i by increasing values of u. When extending a label (i, u, S , q, v), two labels are created:

• a label (i + 1, u, S ′ , q′ , v ′ ), corresponding to the choice of skipping element i + 1, with S ′ = S , q′ = q, v ′ = v , ′′ ′′ • a label (i + 1, u + di+1 , S ′′ , q′′ , v ′′ ), corresponding  to the choice of including element i + 1, with S = S ∪ {i + 1}, qC = ′′ ′′ |C ∩ S |, for each C ∈ C , and v = v − λi+1 − C ∈C¯ ψC where C¯ = {C ∈ C : q′′C = 2 and qC = 1}. This second label is created only if u + di+1 ≤ pl . Since this method can, in principle, generate for each i and u a different label for each q, when each pair of values i and u is considered, we perform the following dominance checks: a label ℓ′ = (i, u, S ′ , q′ , v ′ ) dominates a label ℓ′′ = (i, u, S ′′ , q′′ , v ′′ ) if

v′ −



ψC ≤ v ′′ ,

(23)

C ∈C˜

where

C˜ = {C ∈ C : C \ {1 . . . i} ̸= ∅ and q′C = 1 and q′′C ∈ {0, 2}}. The rationale behind rule (23) is the following. When C \ {1 . . . i} = ∅, the contribution of ψC cannot be triggered in ℓ′ anymore. When |C ∩ S ′ | ≥ 2, then v ′ already includes the contribution of ψC ; at the same time when |C ∩ S ′ | ≤ |C ∩ S ′′ | ≤ 1 any extension of ℓ′ triggering the contribution of ψC would trigger the same contribution in ℓ′′ . Finally, if q′′C = 2 then at most one element of C follows i in the given order, and therefore if q′C = 0, the contribution of ψC will never be triggered in ℓ′ ; the same applies if q′′C = 3 (that is, no element of C follows i in the given order), and q′C ≤ 1. Only two cases are left open: (q′C = 1, q′′C = 0) and (q′C = 1, q′′C = 2). We remark that only labels having the same i and u values need to be compared during the dominance checks, and that when a particular (i, u) pair is considered, all labels for that pair have already been generated; therefore it is still possible to store labels in a table having one row for each value of i and one column for each value of u, and keep the dynamic programming algorithm very efficient also after the introduction of subset-row cuts. 4.2.2. 2-path inequalities The 2-path inequalities, originally introduced by Kohl et al. [15], were developed for the capacitated VRP with time windows. They impose that at least two routes must be used to visit every set Q ⊆ N of nodes that cannot be served by a single route either due to capacity or time limitations. Given such a set Q , and taking into account that sites can be skipped or served by either routes or DCs, we adapted the 2-path inequalities to the GLDP as follows:



si +

 w∈W (Q )

i∈Q

eQw yw +



Q

gk zk ≥ 2 ∀Q ⊆ N

(24)

k∈K (Q )

where K (Q ) and W (Q ) denote respectively the set of routes and the set of clusters serving at least one site belonging to Q . Q Coefficients eQw and gk are equal to

• two if the variables yw or zk encode respectively a cluster or a route containing all sites in Q , • one if a cluster or route serves Q only partially, • zero otherwise. In this way we obtain a cut that, if consistently updated with the newly generated columns, is always valid for our problem. We observe that, in the traditional VRP, given a set Q that cannot be served by a single route but which is visited by less than two vehicles in a fractional solution it is always possible to find a violated 2-path inequality. This is not always true Q for the GLDP, as all coefficients eQw and gk for the columns in the RMP can be equal to two. In order to avoid useless checks we perform three steps: we identify sets Q potentially leading to a violated cut (identification step), we check necessary conditions for a violated cut based on Q to exist (filtering step) and only for the unfiltered Q sets we proceed to actual separation (coefficient setting). The generation step is carried out as in [15], considering sets Q having

 i∈Q

si +

 l∈L w∈Wl

ai w y w +

  l∈L h∈H k∈Klh

bik zk < 2.

A. Ceselli et al. / Discrete Applied Mathematics 165 (2014) 130–145

137

Then we filter each such set Q by checking the following necessary conditions. Let L¯ ⊆ L be the set of locations hosting DCs serving some customer in Q in the current fractional solution, let L˜ ⊆ L be the set of locations hosting depots, where routes serving customers in Q start and let H˜ ⊆ H be the set of vehicle types used to serve customers in Q . In order for a violated cut to exist, at least one of the following conditions must hold:

• • • •

a DC l ∈ L¯ exists, having not enough capacity to fully serve Q



di > pl ;



i∈Q

a DC l ∈ L¯ exists, having sites in Q outside of its range (maxi∈Q {cli } > rl );   a vehicle type h ∈ H˜ exists, having not enough capacity to fully serve Q i∈Q di > fh ; a depot l ∈ L˜ exists, from which no route can visit all sites in Q within the deadline T , i.e. such that HP (Q ∪ {l}) > T , where HP (Q ∪ {l}) is the time taken to perform a Hamiltonian Path on Q ∪ {l}.

In our algorithm the Hamiltonian Path computation required to check the last condition is performed by means of an exact dynamic programming procedure taken from [5]. We check these conditions in this order. As soon as one of them is found to be true, we move to the coefficient setting step; if all of them are false, no violated cut based on Q can be found, and therefore we skip the coefficient setting step and we consider another set Q . Q Finally, in the coefficient setting step, all columns are checked to fix the values of coefficients eQw and gk according to the definitions given above. As outlined before, when all e and g coefficients take value two, the generated inequality does not help in cutting the current fractional solution; however, since it might become useful when further columns are added to the RMP, it is added to the problem in any case. Let the set Q be the collection of all the subsets Q whose corresponding 2-path inequality is in the RMP. A new set of dual variables γQ ≥ 0 is introduced, one for each Q ∈ Q, that requires the modification of the pricing algorithm. As before, we adapted the exact dynamic programming algorithm for the RCESPP as described in [5], and in the remainder we only focus on the DC pricing algorithm. The labels used in the DC pricing algorithm become (i, u, S , q, o, v), where o is a vector having one component oQ for each set Q ∈ Q whose 2-path inequality is in the RMP, represented the number of elements of Q included in S . Our label pushing algorithm is modified as follows. First, we create an initial label (0, 0, ∅, (0 . . . 0), (0 . . . 0), 0). Then, we iteratively extend labels in order of increasing values of i and u, and when extending a label (i, u, S , q, o, v), two labels are created:

• a label (i + 1, u, S ′ , q′ , o′ , v ′ ) corresponding to the choice of skipping element i + 1, with o′ = o and v ′ = v , • a label (i + 1, u + di+1 , S ′′ , q′′ , o′′ , v ′′ ) corresponding to the choice of including element i + 1, with o′′Q = |Q ∩ S ′′ |, for each set Q ,

v ′′ = v − λi+1 −



ψC −



γQ −

¯′ Q ∈Q

C ∈C¯



γQ

¯ ′′ Q ∈Q

where

¯ ′ = {Q ∈ Q : o′′Q > 1 and oQ = 0} Q and

¯ ′′ = {Q ∈ Q : o′′Q = |Q | and oQ < |Q |}. Q We remark that now the contribution of each γQ can be added to the reduced cost of a label up to two times: when the first element of Q is visited (since the eQ w coefficient in the corresponding column changes from 0 to 1), and when all elements of Q are visited (since eQ changes from 1 to 2). w Resources S and q are updated as before. After this modification to the pricing algorithm, a label ℓ′ = (i, u, S ′ , q′ , o′ , v ′ ) can dominate a label ℓ′′ = (i, u, S ′′ , q′′ , o′′ , v ′′ ) if

v′ −

 C ∈C˜

ψC ≤ v ′′ −

 ˜′ Q ∈Q

γQ −



γQ ,

(25)

˜ ′′ Q ∈Q

where

˜ ′ = {Q ∈ Q : Q \ {1 . . . i} ̸= ∅ and o′′Q = 0 and o′Q > 0} Q and

˜ ′′ = {Q ∈ Q : |Q ∩ {1 . . . i}| = o′′Q and o′′Q > o′Q }. Q In this case, a contribution γQ may be gained by ℓ′′ , but not by ℓ′ in two cases, modeled by the last two terms of expression (25). First, ℓ′′ may include for the first time an element of Q , while ℓ′ cannot; this can happen only if (a) elements of Q can still be included in S ′′ , (b) no element of Q was already included in S ′′ , and (c) at least an element of Q was already included in S ′ . Second, ℓ′′ may include all the elements of Q , while ℓ′ cannot; this can happen only if (a) all those elements 1, . . . , i belonging to Q were selected in S ′′ , while some of them were not selected in S ′ . Even if the dominance criterion is slightly loosened, we experimentally observed that the DC pricing algorithm remains very efficient, also after the introduction of 2-path inequalities.

138

A. Ceselli et al. / Discrete Applied Mathematics 165 (2014) 130–145

(a) A partial fractional solution respecting (8), but violating (26).

(b) A partial fractional solution respecting (26), but violating (8).

Fig. 2. Effect of constraints (8) and (26) on fractional solutions.

4.2.3. Consistency cuts We identified constraints (8) as the most likely cause of a possible poor linear relaxation of the MP since they may promote the use, in a fractional solution, of a same location both as a DC and as a depot. In fact, for every location l ∈ L and every kind of vehicle h ∈ H , it may be convenient to either set several z variables to one assigning only a fractional value to the corresponding yw variable or to assign a fractional value to many y and z variables corresponding to the same location in order to partially cover as many sites as possible. An example of such a fractional solution is shown in Fig. 2(a). For this reason we consider the following set of inequalities: zk +



yw ≤ 1 ∀h ∈ H , ∀l ∈ L, ∀k ∈ Klh .

(26)

w∈Wl

They impose that if any route starting from a location l ∈ L is used, then no DC can be placed in l. We remark that constraints (26) and (8) do not dominate one another. In Fig. 2(b) a fractional solution is depicted, still partially using the same location as both a depot and a DC: it respects (26) but violates (8). Since constraints (26) are exponential in number, a dynamic generation strategy is adopted. In particular every time a variable zk encoding a route starting from the depot l is introduced in the RMP, the corresponding constraint is generated as well. This constraint includes also all yw variables for all w ∈ Wl already in the RMP. In the same way, every time a variable yw encoding a cluster of the DC l is generated, it is added to every constraint of type (26) whose corresponding zk variable encodes a route starting from l. The introduction of these constraints requires a modification of the pricing algorithms. Let P and ν be respectively the set of all inequalities of type (26) that are in the RMP and their associated non-positive dual variables. Let P (l) be the subset of P made up of all the constraints which involve variables yw and zk encoding a route and clusters using location l. Then, a penalty term p∈P (l) νp appears in the DC column reduced cost. These terms are constant in the DC pricing subproblem, as a different subproblem is solved for each l; therefore, no modification to the DC pricing algorithm is needed. On the other hand, each zk variable appears in a single inequality of type (8). In principle, if zk has not already been generated, then its corresponding constraint is not in the RMP and therefore we can consider the corresponding dual variable to have value zero; instead, if zk is already in the RMP, we are not interested in generating it twice. A direct implementation of this principle would require us to provide, by solving the route pricing problem, the most appealing route excluding those already in the RMP. In order to improve the performance of the dynamic programming algorithms we devised the following technique: when checking dominance conditions between each pair of labels ℓ′ and ℓ′′ encoding partial paths of reduced cost v ′ and v ′′ , respectively, we let ℓ′ dominate ℓ′′ only if

v ′ − ν(ℓ′ ) ≤ v ′′ where ν(ℓ) is the variable νp of minimum value whose associated route is in the RMP, and either corresponds to the partial path encoded by ℓ′ or can be obtained by an extension of ℓ′ . If no such route exists, then ν(ℓ) = 0. The joining phase of bidirectional dynamic programming has also to be modified by checking, for each pair of labels, if their join yields a route which is already in the RMP; in this case, the corresponding route is simply discarded. The effectiveness of this method strongly relies in quickly identifying routes in the RMP. To this aim we designed a tree data structure, in which each node corresponds to a partial path, belonging to a route in the RMP, and contains a suitable dual value; one constraint of the set (8) corresponds to each leaf of this tree: before each pricing phase the dual values of these constraints are pushed from the leaves up to the root, in such a way that each node contains the minimum value between those of its children. During the execution of the dynamic programming algorithms, each value ν(ℓ) can then be quickly found by visiting this tree. This data structure is described in detail in [21].

A. Ceselli et al. / Discrete Applied Mathematics 165 (2014) 130–145

139

Fig. 3. The sequence of the pricing algorithms.

4.3. Structure of the price-and-cut loop In order to obtain the dual bound we first perform column generation. Pricing algorithms are called in sequence starting from the dynamic programming procedure for the DC pricing; when no negative reduced cost cluster column is found, then the route pricing algorithms are called. An outline of the invoking sequence of the pricing algorithms is shown in Fig. 3.

140

A. Ceselli et al. / Discrete Applied Mathematics 165 (2014) 130–145

At each column generation iteration we also insert in the RMP each consistency cut (8) corresponding to new routes. When no more columns of negative reduced cost can be found we switch to the cut generation step. A comparison of different cutting strategies is carried out in Section 5. We obtained the best results as follows. First we try to separate subset-row inequalities. As soon as a violated cut is found, the cut generation step is suspended, and we switch back to column generation. If no violated subset-row inequality is found, then 2-path inequalities are considered. Also in this case, the cut generation is stopped as soon as any violated inequality is found. This price-and-cut loop is repeated until neither useful columns nor violated cuts can be found. 4.4. Branching When the RMP solution obtained at the end of the price-and-cut loop is fractional, we recover integrality by branching. We use five different policies, that are considered in the order presented hereafter; branching is executed using the first applicable one. In the remainder we indicate as s¯, y¯ and z¯ the values taken by the variables in the fractional solution. Branching on locations. At first we force the exclusive use of an open location l ∈ L as either a depot or a DC. In detail, let W (l) and  K (l) be the sets ofvariables encoding respectively clusters or routes whose associated DC and depot are located in l. If w∈W (l) y¯ w > 0 and k∈K (l) z¯k > 0 then we branch as follows: two children nodes are generated; we forbid l to be used as a depot in the first, and as a DC in the second. This is done by removing all columns in the RMP belonging to K (l) (resp. W (l)), and avoiding to generate new route (resp. DC) columns pricing. Among all those locations having  for l during this feature, we select for branching the location with maximum w∈W (l) y¯ w + k∈K (l) z¯k value. Branching on skip sites. Second, let s¯i be the variable whose value is fractional and closest to one. Two children nodes are created by branching: in the first site i is left unserviced by fixing s¯i = 1, while in the second site i is forced to be served by fixing s¯i = 0. Branching on clusters. Third, we evaluate the sum of the variables representing clusters serving each site i as

Di =



aiw y¯ w .

l∈L w∈Wl

If this value is not an integer then we identify the set of DCs giving a fractional contribution to Di , that is



 L(i) =

l∈L:



aiw y¯ w > 0 .

w∈Wl

Then we select a site i with maximum |L(i)| value, and we create |L(i)| + 1 children. In the first |L(i)| we impose in turn to serve site i with a particular DC in L(i); in the last child we forbid to serve i with all DCs in L(i). These branching choices have no effect on the route pricing algorithm, and can easily be handled in the DC pricing algorithms by changing the initialization of the first label and forcing us to skip some elements during the extension process.





¯k be the number of routes using a vehicle of type h ∈ H . Branching on the number of vehicles. Fourth, let mh = l∈L k∈Klh z Binary branching is executed on the vehicle type having mh value closest to 0.5. Two children nodes are generated, imposing to use at least mh + 1 and at most mh vehicles of type h, respectively. This is obtained by modifying the left and right hand side of the constraint of type (7) corresponding to the selected h; in particular in the first child node the modified constraint becomes mh ≤



zk ≤ Vh

l∈L k∈Klh

while in the second child node the constraint becomes 0≤



zk ≤ mh .

l∈L k∈Klh

This branching technique leaves the pricing subproblem unchanged. Branching on arcs. Finally, we choose the site i ∈ S that is split among the largest number of routes in the RMP fractional solution and we create two children nodes, forbidding half of the arcs leaving i in each of them. To handle these branching decisions in the route pricing algorithm it is enough to set the travel time of forbidden arcs to +∞ (i.e. a very large constant value). No modification to the DC pricing algorithm is required. 4.5. Primal heuristic algorithm In order to find feasible solutions early during the exploration of the branching tree we have developed a primal heuristic. It looks for feasible solutions by considering only the variables that are in the RMP. The algorithm goes through five different phases: locations’ selection, DCs’ allocation, allocation of vehicles, improvement and feasibility and it is executed once for every node of the branching tree. Let y¯ , z¯ and s¯ be the values taken by the variables in the optimal RMP fractional solution. Hereafter we indicate with I the set of column variables to be fixed to one in the heuristic solution, and with N¯ the set of sites covered by those columns. Both I and N¯ are initialized empty.

A. Ceselli et al. / Discrete Applied Mathematics 165 (2014) 130–145

141

Locations’ selection. In the first phase every location l ∈ L is marked as either potential DC or depot. The algorithm computes    ¯ w and every potential depot udep ¯k ; then every the use of every potential DC location udc = = l l w∈Wl y h∈H k∈Klh z dep

location is marked as potential DC if udc l ≥ ul

and as potential depot otherwise.

DCs’ allocation. In this phase only variables in the RMP that encode clusters are taken into account. Let L(dc ) be the set of dep locations marked as potential DCs, that is L(dc ) = {l ∈ L : udc l ≥ ul }. All locations in L(dc ) are sorted in non-increasing dc order of ul . Then, following this order, each location l ∈ L(dc ) is considered, and the covering potential



cpdc w =

aiw di

i∈N \N¯

of every variable yw ∈ Wl is computed; it corresponds to the sum of the demands of the sites covered by the cluster described by yw that are not already covered by elements in I. A variable yw whose corresponding cpdc w is maximum is selected and added to the solution set I; N¯ is updated accordingly. Then the algorithm moves to the next location in L(dc ) and this process is iterated until either all locations in L(dc ) have been considered or B variables have been added to the set I. If no cluster yw ∈ Wl associated with a potential location l is selected, then l is removed from L(dc ) and marked as potential depot location. Vehicles’ allocation. In this phase only variables in the RMP encoding routes are taken into account. Each vehicle type h ∈ H is considered, in non-increasing value of vehicle capacity fh , and the covering potential



cpro k =

bik · di

i∈N \N¯



of every variable zk ∈ l∈L\L(dc ) Klh is computed; as before, it corresponds to the sum of the demands of the sites covered by the route that are not already covered by elements in I. A variable zk whose cpro k is maximum is selected and added to I, and N¯ is updated accordingly. This process is repeated until either there are no more variables associated with vehicle type h or Vh variables are selected. Then the algorithm moves the next vehicle type. This phase terminates when all types of vehicles have been taken into account. It is worth noting that at this point the value of the solution can be computed as the sum of the demand of the sites i ∈ N \ N¯ . Improvement. We apply a simple local search procedure, aimed at improving the quality of the solution found. In detail, the algorithm evaluates if every variable in I can be replaced by an unselected variable, improving the value of the solution. This is done by removing in turn a single variable yw or zk from I and recomputing the covering potential for all the variables not in I that are compatible with the removed one; by compatible we mean either related to the same DC or to the same pair of depot location and vehicle type. If any compatible variable not in I is found, having better covering potential, we perform the substitution, update set N¯ accordingly and iterate. The whole procedure is iterated twice, first considering the selected variables in their order of insertion in I, and then following the reverse order. Feasibility. Finally, to ensure the feasibility of the solution, the algorithm includes in I all skip variables related to sites in N \ N¯ , that are those uncovered in the current solution. At this point a complete solution for the GLDP is obtained by fixing to one the value of all variables in I, and fixing to zero the value of any other variable. We observe that, while this procedure always generates feasible solutions when executed at the root, it might produce infeasible solutions in the inner nodes of the search tree, due to branching decisions. 5. Experimental analysis We coded our algorithm in C++, using SCIP 2.1 [1] as a branch-and-cut-and-price framework linked to CPLEX 12.2 as an LP solver. SCIP takes care of the management of the columns and rows pool and it is able to remove and re-insert them as needed. All pre-solving algorithms and automatic cut generators embedded in SCIP along with the following general purpose heuristic algorithms: shift-and-propagate, fix-and-infer, int-diving, DINS, RENS, RINS, undercover, mutation, crossover and shifting have been disabled as they were considered either incompatible or useless for our problem. All remaining SCIP parameters were kept at their default values including the branching tree exploration strategies. R The experimental campaign was performed on a single core of a PC Intel⃝ Core 2 DuoTM CPU T7300 (2.00 GHz) with 4 GB RAM and running Ubuntu 10.04 operating system. A time limit of one hour was imposed on each run. Benchmark instances. In order to properly test our approach we considered two datasets drawn from the literature. Dataset 1 is designed to fully test the features of our method, and is used in the experiments detailed in Sections 5.1 and 5.2: we generated 440 instances starting from the approach detailed by Shen et al. [20]. We took into account three different problem sizes: 10, 20 and 50 delivery sites, setting the number of potential facility locations to 2, 3 and 5 respectively. For each problem size we randomly generated five different scenarios, placing sites and location in the Euclidean plane, drawing both their coordinates and the demands of each site from a uniform distribution. We built different instances for every scenario, considering different amounts of available supply D and distribution time deadline M. In particular, being

¯ = D

 i∈N

di

142

A. Ceselli et al. / Discrete Applied Mathematics 165 (2014) 130–145 Table 1 Configurations AC and CC: percentage variation compared with NC.

|N | 10 20 50

AC

CC

B.B.N.

L.B.

Time

B.B.N.

L.B.

Time

−7.59

1.42 0.00 0.00

34.88 33.15 −36.56

−1.90

1.20 0.00 0.00

16.28 21.00 −33.38

5.33 −46.86

16.77 −21.58

¯ , 0.8D¯ , 0.9D¯ and 1.0D. ¯ For the time deadline, being the overall demand, we considered four values for D: 0.7D  ¯ = M

i∈N ,j∈N

cij

|N |2

·

|N | T

¯ , 0.5M ¯ , 0.6M ¯ , 0.8M¯ and 1.0M. ¯ In instances with 10 and 20 the base route length, we considered five values for M: 0.4M ¯ , 1.6D¯ and 2.0D), ¯ and only in instances with 10 delivery sites a delivery sites we considered additional values for D (1.2D ¯ ) was also used. One instance was generated for every problem size and scenario, and for every further value for M (1.2M combination of D and M, obtaining a dataset of 440 instances. In every instance three different types (V1 , V2 and V3 ) of vehicles were considered whose capacity was set, respectively, to 60%, 80%, and 120% of the average delivery capacity computed as DT . In instances with 10 sites T = 4, V1 = 2, V2 = 1, V3 = 1, in instances with 20 sites T = 5, V1 = 2, V2 = 2, V3 = 1 and in instances with 50 sites T = 11, V1 = 5, V2 = 4, V3 = 2. The value of B was always set to 0.5 · |L|. Dataset 2 is taken from [20], and is used in Section 5.3 to perform a computational comparison between our method and that of [20]. The dataset is composed of 45 instances with 50 delivery sites each: the data on delivery site demand and coordinates were kindly provided to us by the authors; euclidean distances are taken between the delivery sites. Each ¯ and D¯ were fixed to instance involves a homogeneous fleet of T = V1 = 10 vehicles. Base demand and deadline values M ¯ , 0.8M ¯ , 1.2M ¯ and 0.7D¯ , 0.9D¯ , 1.2D. ¯ The delivery capacity 500, and three different values for M and D were considered: 0.4M D . In [20] the option of locating DCs is not considered, and therefore we set B = 0. The for every vehicle was then set to 10 dataset includes five instances for each combination of values M and D. 5.1. Cuts’ comparison We have tested our approach on all instances of Dataset 1 using different combinations of the cuts presented in Section 4.2; in particular five different configurations have been tested: NC: CC: SC: 2C: AC:

no cuts are used; only consistency cuts are used; only subset-row inequalities are used; only 2-path inequalities are used; all cuts (consistency, subset-row and 2-path) are used.

The percentage variations, aggregated by instance size, in the number of the branch-and-bound nodes explored (B.B.N.), in the value of the lower bound at the root node (L.B.) and in the computational time required to obtain the lower bound (Time), compared with the NC configuration, are reported in Tables 1 and 2. We first observe that each cut type provides only a small improvement in the lower bound. In fact, the lower bound is affected by the cuts only on instances with 10 sites and by at most 1.42% with AC configuration. At the same time, the introduction of cuts yields a significant change in the number of nodes explored and in the computational time, in particular on the largest instances. In this case all cuts but the subset-row inequalities are able to reduce both the number of the nodes explored and the computational time. Configurations CC and AC turned out to bring the greatest improvement reaching a reduction by more than 30% in both values. Since a consistency cut is added every time a variable zk is generated, many of them are in the final RMP: for example, the average number of consistency cuts added in instances with 50 nodes is about 500. On the other hand, according to SCIP policies, less than 10 2-path and subset-row inequalities turned out to be effective after their introduction in the RMP, considering instances of the same size. Subset-row inequalities are the less efficient cuts in this problem: they introduce almost as much delay as the consistency cuts but are not able to reduce the number of nodes explored on any instance with more than 10 nodes. 2-path inequalities instead do not imply a relevant additional computational effort and they can reduce the number of nodes explored in instances with 10 and 50 nodes. Anyway, all separation procedures are fast: 2-path inequalities are separated in less than one tenth of a second while the exact subset-row separator runs, on average, in less than one second. No significant additional computational effort is required to introduce consistency cuts. In Table 3 we report the number of instances solved to proven optimality with each configuration of cuts. The maximum number of instances is solved using the CC configuration; AC and 2C follow with respectively a total of 434 and 433. Using only the subset-row inequalities does not increment the number of instances solved compared with the NC configuration. In the next section we give a detailed analysis of the results achieved with the CC configuration.

A. Ceselli et al. / Discrete Applied Mathematics 165 (2014) 130–145

143

Table 2 Configurations 2C and SC: percentage variation compared with NC.

|N | 10 20 50

2C

SC

B.B.N.

L.B.

Time

B.B.N.

L.B.s

Time

−1.58

0.07 0.00 0.00

7.75 8.94 −11.29

−7.28

0.19 0.00 0.00

31.01 17.67 19.99

4.08 −7.73

1.34 12.97

Table 3 Number of instances solved to proven optimality for each configuration of cuts.

|N |

NC

AC

CC

2C

SC

10 20 50 Total

210 130 90 430

210 130 94 434

210 130 95 435

210 130 93 433

210 130 90 430

Table 4 Aggregated experimental results: nodes in the branch-and-bound tree, percentage gaps, computing time.

|N |

Gap R%

Gap F%

10

Avg Max

B.B.N. 1.48 19.00

1.27 63.17

0.00 0.00

Time [s] 0.01 0.05

20

Avg Max

28.81 305.00

1.43 28.57

0.00 0.00

0.98 17.27

50

Avg Max

80.27 1046.00

0.21 6.38

0.21 6.38

387.62 3600.19

5.2. General performance Being able to efficiently solve pricing problems is a key issue in column generation methods, and sophisticated heuristic pricing techniques are proposed in the literature [10]. Indeed, in our algorithm we found it useful to proceed to heuristic route generation: as described in Section 4 we first run a very fast greedy algorithm, then we switch to a more accurate dynamic programming heuristic and only afterwards run the exact dynamic programming algorithm. The use of weak pricing routines is particularly controversial, as they are often unable to find near-optimal columns in hard instances. Therefore, we first evaluated the contribution given by the greedy pricing routine to the generation of good routes. In particular, we measured the CPU time needed to complete the column generation process at the root node with and without activating such a procedure, for all instances of Dataset 1. In our tests we observed two distinct behaviors: on 3.95% of the instances the columns produced by greedy pricing actually yielded tailing off effects, producing marginal improvements at each iteration, and slowing down the computation by more than a factor two; on the other hand, on the largest part of the dataset, the greedy pricing produced an average speedup of 24.9%, showing to be in general beneficial in our algorithm. Hence we decided to keep it active in the final version of the code. The final computational results on Dataset 1 are reported in Tables 4 and 5. We report only aggregated values: detailed results for each instance can be found in [21]. For each instance size Table 4 shows the average and maximum number of nodes explored in the branch-and-bound tree (B.B.N.), the gap at the root node and at the end of the computation, taken as the difference between primal and dual bounds divided by the dual bound (Gap R%, Gap F%) and the computational time in seconds (Time [s]). Table 5 reports the percentage of instances solved (Solved), the percentage of instances solved at the root node (R. Solved) and the percentage of instances in which the lower bound at the root node was equal to the value of the optimal solution (Opt LB). Our algorithm was able to optimally solve 435 instances of the 440 included in the dataset. In particular all instances with 10 and 20 sites were closed in less than 20 s. The residual gap considering only the five instances with 50 nodes not solved to optimality ranges from 2.22% to 6.38%. More than half of the instances were solved at the root node; when branching was required, the portion of the search tree explored before an optimal solution was reached did not exceed 100 nodes on average. The last column of Table 5 shows an interesting result: it contains the percentage of instances in which the value of lower bound at the root does not further increase during branching, being equal to the value of the optimal integer solution. This value is about 95% and it is a certification of the strength of the Set Covering formulation. Moreover it also helps in explaining the poor contribution given by the additional cuts: indeed if the possibly fractional value of the RMP is very close (or even equal) to the value of the optimal integer solution very few useful cuts can be added at the root node. Taking that into account, finding good integer feasible solutions seems much harder. Indeed, although the duality gap is almost always null, finding an optimal integer solution may require several branching operations.

144

A. Ceselli et al. / Discrete Applied Mathematics 165 (2014) 130–145 Table 5 Aggregated experimental results: number of instances solved and lower bound improvement.

|N |

Solved %

S. Root %

Opt LB

10 20 50

100.00 100.00 95.00

83.81 36.15 48.00

92.38 93.42 100.00

98.33

55.99

94.92

AVG

Table 6 Results on bioterrorism emergency instances. D

350 450 600

M = 200

M = 400

M = 600

Opt 200 s

Opt 1h

t (s)

Opt 200 s

Opt 1h

t (s)

Opt 200 s

Opt 1h

t (s)

5 5 5

5 5 5

0.18 86.45 11.40

5 2 5

5 5 5

0.46 404.86 0.64

5 1 5

5 4 5

0.75 719.72 0.79

5.3. Results on bioterrorism emergency instances Finally, we performed a computational comparison with the approach presented in [20], that was the starting point for the development of our models. Since [20] does not include the option of DCs, we turned off the DC pricing algorithms, but all other settings were kept as in the full version. In Table 6 we report the computational results on instances of Dataset 2. The table has one row for each value of D, and one vertical block for each value of M. The results reported in [20] were obtained by letting a tabu search algorithm run for 200 s on a machine equipped with a 3.2 GHz Pentium IV processor. For the sake of comparison, in each row and for each block we detail first the number of instances with a particular combination of D and M solved to proven optimality by our method within 200 s of CPU time; we subsequently report the number of instances solved within a time limit of one hour, and the average CPU time needed. When setting a time limit of 200 s, on instances with M = 400 and D = 450, and M = 600 and D = 450, the method described in [20] is competitive, yielding residual gaps of about 0.3%, while our algorithm achieves gaps of about 3.5% and 16.3%, respectively. However, on all the remaining instances, our algorithm yields substantially better results; in particular, instances with D = 600 seem critical for the method described in [20], while our algorithm always reaches proven optimality. In about 85% of the cases it takes less CPU time to solve an instance to proven optimality with our algorithm than to obtain a heuristic solution with tabu search. Furthermore, our algorithm stops on the average well before reaching 200 s of computation, being able to prove optimality early, whereas tabu search keeps trying new solutions until the time limit. Finally, we remark that in all instances but one our algorithm found a proven optimal solution within a time limit of one hour, confirming the results of the previous tests. 6. Conclusions In this paper we have presented a Generalized Location and Distribution Problem, with a possible application to the distribution of drugs in the case of emergency. The GLDP goes beyond the definition of location and routing: not only it represents an extension of the classical multi-depot heterogeneous fleet VRP with profits but it also introduces a new degree of freedom in the design of the logistics system, by allowing for multiple distribution modes to be combined together. We have presented a Set Covering formulation with three types of columns and a branch-and-price-and-cut algorithm to solve it. The algorithm includes multiple pricers, three different types of strengthening cuts, a primal heuristic procedure and several branching strategies. Two of the cut types presented here are a generalization of cuts from the literature, that required us to devise new separation procedures. The third type of cut is new in the context of column generation algorithms. When run on a dataset of 45 instances from the literature, our algorithm proves to outperform the existing methods, offering at the same time more flexibility. Furthermore, our approach has been tested on an additional set of 440 instances showing good results and reaching the optimal solution in almost all the instances within a time limit of one hour; when the optimality was not proved the algorithm provided very tight lower bounds and primal solutions. The structure of the solutions suggests that our algorithm is able to quickly find the optimal value, while the largest part of the computational time is spent in finding the optimal integer solution. This opens the interesting question of devising a fast and reliable specialized heuristic. Acknowledgments This work has been developed in the Operations Research Laboratory ‘‘OptLab’’ supported by ACSU—Associazione Cremasca Studi Universitari and it is part of the project ‘‘Column generation algorithms for routing problems’’ partially funded by the Italian Ministry for Education, University and Research.

A. Ceselli et al. / Discrete Applied Mathematics 165 (2014) 130–145

145

The authors wish to thank Dr. Andrea Bettinelli for the insightful discussions on the topic of the paper, the authors of [20] for readily providing their instance data and the anonymous referees for their useful comments. References [1] T. Achterberg, Scip: solving constraint integer programs, Mathematical Programming Computation 1 (1) (2009) 1–41. http://mpc.zib.de/index.php/ MPC/article/view/4. [2] R. Baldacci, A. Mingozzi, R. Roberti, New route relaxation and pricing strategies for the vehicle routing problem, Operations Research 59 (2011) 1269–1283. [3] R. Baldacci, A. Mingozzi, R. Wolfler-Calvo, An exact method for the capacitated location-routing problem, Operations Research 59 (2011) 1284–1296. [4] J.M. Belenguer, E. Benavent, C. Prins, C. Prodhon, R. Wolfler-Calvo, A branch-and-cut algorithm for the capacitated location routing problem, Computers and Operations Research 38 (2011) 931–941. [5] A. Bettinelli, A. Ceselli, G. Righini, A branch-and-cut-and-price algorithm for the multi-depot heterogeneous vehicle routing problem with time windows, Transportation Research Part C: Emerging Technologies 19 (5) (2011) 723–740. [6] CDCs Division of Strategic National Stockpile. Emergency medkit evaluation study summary - background, key results, and next steps, November 2007. Available online: http://www.bt.cdc.gov/cri/pdf/medkit-evaluation-summary-2007updated.pdf. [7] A. Ceselli, G. Righini, E. Tresoldi, Modeling and solving profitable location and distribution problems, Optimization Letters (2012) in press (available online). [8] C. Contardo, J.-F. Cordeau, B. Gendron, A branch-and-cut-and-price algorithm for the capacitated location-routing problem, Technical Report CIRRELT2011-44, Université de Montrèal, 2011. [9] G.B. Dantzig, Discrete-variable extremum problems, Operations Research (1957). [10] G. Desaulniers, F. Lessard, A. Hadjar, Tabu search, partial elementarity, and generalized k-path inequalities for the vehicle routing problem with time windows, Transportation Science 42 (2008) 387–404. [11] J. Desrosiers, G. Desaulniers, M.M. Solomon, Column Generation, Springer, 2005. [12] D. Feillet, P. Dejax, M. Gendreau, C. Gueguen, An exact algorithm for the elementary shortest path problem with resource constraints: application to some vehicle routing problems, Networks 44 (3) (2004) 216–229. [13] S. Irnich, D. Villeneuve, The shortest path problem with resource constraints and k-cycle elimination for k = 3, INFORMS Journal on Computing 18 (2006) 391–406. [14] M. Jepsen, B. Petersen, S. Spoorendonk, D. Pisinger, Subset-row inequalities applied to the vehicle-routing problem with time windows, Operations Research 56 (2) (2008) 497–511. [15] N. Kohl, J. Desrosiers, O. Madsen, M. Solomon, F. Soumis, 2-path cuts for the vehicle routing problem with time windows, Transportation Science 33 (1999) 101–116. [16] L. Ladányi, T.K. Ralphs, L.E. Trotter, Computational Combinatorial Optimization: Optimal or Provably Near-Optimal Solutions, in: Lecture Notes in Computer Science, vol. 2241, Springer, 2001, pp. 223–260. Chapter branch, cut, and price: sequential and parallel. [17] S. Martello, P. Toth, Knapsack Problems: Algorithms and Computer Implementations, Wiley, New York, 1990. [18] G. Nagy, S. Salhi, Location-routing: issues, models and methods, European Journal of Operational Research 177 (2007) 649–672. [19] G. Righini, M. Salani, New dynamic programming algorithms for the resource constrained elementary shortest path problem, Networks 51 (3) (2008) 155–170. [20] Z. Shen, F. Ordonez, M. Dessouky, A two-stage vehicle routing model for large-scale bioterrorism emergencies, Networks 54 (4) (2009) 255–269. [21] E. Tresoldi, Location and routing problems: an unified approach. Technical report, DTI, Università degli Studi di Milano, Milano, 2011. [22] P. Vansteenwegen, W. Souffriau, D. Van Oudheusden, The orienteering problem: a survey, European Journal of Operational Research 209 (2011) 1–10.