A memetic algorithm for the capacitated location-routing problem with mixed backhauls

Author's Accepted Manuscript

A memetic algorithm for the capacitated location-routing problem with mixed backhauls Ismail Karaoglan, Fulya Altiparmak

www.elsevier.com/locate/caor

PII: DOI: Reference:

S0305-0548(14)00167-1 http://dx.doi.org/10.1016/j.cor.2014.06.009 CAOR3591

To appear in:

Computers & Operations Research

Received date: 30 June 2013 Revised date: 25 March 2014 Accepted date: 7 June 2014 Cite this article as: Ismail Karaoglan, Fulya Altiparmak, A memetic algorithm for the capacitated location-routing problem with mixed backhauls, Computers & Operations Research, http://dx.doi.org/10.1016/j.cor.2014.06.009 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

A Memetic Algorithm for the Capacitated Location-Routing Problem with Mixed Backhauls

Ismail Karaoglana1, FulyaAltiparmakb

a

Department of Industrial Engineering,

Faculty of Engineering, Selcuk University, 42075 Selcuklu/Konya Turkey e-mail: [email protected]

b

Department of Industrial Engineering,

Faculty of Engineering, Gazi University, 06570 Maltepe/Ankara Turkey e-mail: [email protected]

A Memetic Algorithm for the Capacitated Location-Routing Problem with Mixed Backhauls Abstract The design of distribution networks is one of the most important problems in supply chain and logistics management. The main elements in designing a distribution network are location and routing decisions. As these elements are interdependent in many distribution networks, the overall system cost can decrease if location and routing decisions are simultaneously tackled. In this paper, we consider a Capacitated Location-Routing Problem with Mixed Backhauls (CLRPMB) which is a general case of the capacitated location-routing problem. CLRPMB is defined as finding locations of the depots and designing vehicle routes in such a way that pickup and delivery demands of each customer must be performed with the same vehicle and the overall cost is minimized. Since CLRPMB is an NP-hard problem, we propose a memetic algorithmto solve the problem. To evaluate the performance of the 1

Corresponding Author

1

proposed approach, we conduct an experimental study and compare its results with the lower bounds obtained by the branch-and-cut algorithm on a set of instances derived from the literature. Computational results indicate that the proposed approach is able to find optimal or very good quality solutions in a reasonable computation time. Keywords: Capacitated location-routing problem, mixed backhauls, memetic algorithm. 1. Introduction The distribution network design (DND) is one of the most important problems in supply chain and logistics management. Distribution costs often represent an important portion of overall system cost and substantial savings can be achieved by improving distribution system. Among the others, location of facilities (i.e. factory, depot, warehouse etc.) can be thought of a crucial component in DND, because of the large setup and operating costs. In general, the classical facility location problem requires that the customers must be served directly. This situation is true if the demand of more than one customer cannot be loaded to the same vehicle because of the vehicle capacity. However, in many applications in practice, demand of customers is less than the vehicle capacity and deliveries are made on a route in which two or more customers are visited sequentially. Therefore, locating the facilities without considering vehicle routes may lead to suboptimal solutions [1]. Location Routing Problem (LRP) overcomes this drawback by considering vehicle routes while locating the facilities. In the general form, LRP deals with determining the location of facilities and the routes of the vehicles for serving the customers under some constraints such as facility and vehicle capacity, route lengths, etc. to satisfy demands of all customers and to minimize total cost including routing costs, vehicle fixed costs, facility fixed costs and facility operating costs. Some of the application areas of LRP in practice can be given as food and drink distribution, military equipment location, parcel delivery and telecommunication network design. The facility location problem (FLP) and vehicle routing problem (VRP) are two main components of LRP. Since both problems belong to the class of NP-hard problems, LRP is also NP-hard. Several mathematical models and exact solution procedures were developed for small- and medium-size LRPs in the literature [2-4]. Because of its complexity, classical heuristic approaches and meta-heuristic approaches are successfully implemented for the problem. Comprehensive reviews of location-routing models and their applications are provided in [5-7] and more recently in [8, 9].

2

In the LRP literature, it is seen that the researchers have considered classical VRP (i.e. each vehicle starts from a depot, traverses through a number of customers, delivers goods to each customer and returns the same depot). However, in practice, customers can have pickup and delivery demands and they request that both demands should be met at the same time. This kind of problem is known in the literature as vehicle routing problem with mixed backhauls (VRPMB). This problem has been studied extensively in the literature. We refer the interested readers to the review papers of Berbeglia et al. [10] and Parragh et al. [11] for the VRPMB and its variants. Because of the increasing importance of LRP and VRPMB in practice and scientific researches, in this paper, we deal with a variant of LRP called capacitated LRP with mixed backhauls (CLRPMB). CLRPMB is an extension of LRP in terms of types of the customers’ demand and the depot capacities. CLRPMB can be considered as an extension of the travelling salesman location problem with pickup and delivery (TSPPD) introduced by Moshiev [12] in terms of the number of depots to be located and the capacity of vehicles. CLRPMB can also be considered as a special case of many to many LRP introduced by Nagy and Salhi [13] in which several customers wish to send goods to others. Some of the application areas of CLRPMB in practice can be found in the reverse logistics context. For example, in food and beverage industry, firms are responsible for both distributing the foods (beverages) to customers and collecting the packing materials (empty bottles) for reusing or disposing. Thus, the decision maker should take into account the decisions of where to locate the depots and how to serve customers within routes. In this paper, we propose a memetic algorithm (abbreviated as MA) based on genetic algorithm (GA), simulated annealing (SA) and integer programming formulation (IPF) to solve the CLRPMB instances. The MA implements SA and IPF as a local search algorithm recursively, after each recombination procedure of GA. We investigate the performance of MA on a set of instances derived from the literature. This paper is organized as follows: Section 2 considers related literature for the problem. Section 3 presents problem definition and mathematical formulation. Section 4 gives details of memetic algorithm. Section 5 reports computational results and finally Section 6 includes conclusion and future directions on this subject. 2. Related Literature After the importance of the considering location and routing decisions simultaneously was stressed by Webb [14] and Christofides and Eilon [15], different models and solution 3

approaches for the LRP have been proposed in the literature to formulate and solve distribution network design problems. Because of the problem complexity, very limited number of exact solution procedure has been implemented for LRP. One of the earliest exact algorithm for LRP have been introduced by Laporte and Nobert [16] where one uncapacitated depot is to be located and fixed number of capacitated vehicles is to be routed. Slightly different version of this problem has been studied by Laporte et al. [17]. More complicated and realistic version of LRP has been considered by Ambrosino and Scutella [18]. This problem consists of four layered network (i.e. plants, distribution centers, large customers and customers) and complex routing plan between different layers. Authors have been proposed an MIP formulation for the problem and investigated the complexity of the problem on several test problems. More recently, several exact solution procedures have been found in literature for LRP. Berger et al. [3] have developed a set partitioning formulation, new family of valid inequalities and a branch-and-price algorithm for LRP with route length constraints. Belenguer et al. [4] have studied CLRP and presented two new MIP formulations with exponential number of constraints and branch-and-cut algorithm based on these formulations, separately. Baldacci et al. [19] have proposed an efficient lower bounding procedure and an exact solution method for CLRP. The effectiveness of the method has been investigated on several capacitated/uncapacitated test instances. Computational results have showed that their method outperforms the current best-known exact methods for both upper and lower bounds. More recently, Contardo et al. [20] have developed an exact solution procedure for CLRP based on column and cut generation. They have formulated the problem using set partitioning type formulation, developed several new valid inequalities.The results have showed that their bounding procedures are stronger than those of Baldacci et al. [19] for some instances. Due to the increasing complexity of LRP, there are several papers in literature which have been referred to heuristic procedures. Jacobsen and Madsen [21] have consider the two level LRP where the first (second) level consists of factory and transfer points (transfer point and customers), and routing operation may occurs in both levels. The authors have proposed a greedy heuristic based on steiner tree for the problem. Nambiar et al. [22] developed a hierarchical method for the problem. This algorithm begins by opening a single depot. Then the number of opened depot is increased one by one and saving algorithm is called after each increment. These steps are repeated until no improvement is obtained. Similar approaches have been used by Srivastava and Benton [23] and Srivastava [24]. Perl and Daskin [25, 26] have used iterative approach, firstly. This approach consists of two phases (i.e. location and

4

routing phases) and these phases communicate with each other. This approach have been improved by Hansen et al. [27]. More recently, meta-heuristic approaches have been used successfully to solve LRPs. Tuzun and Burke [28] have developed a tabu search algorithm for LRP. This algorithm has also two phases and works based on iterative approach. In the first phase, one depot is opened randomly, the customers are assigned to this opened depot and Clarke and Wright’s saving algorithm [29] is called to construct the vehicle routes. Then the solution is improved by using tabu search algorithm which is used in both phases, independently. Similar approach is also used by Bouhauf et al. [30]. In that paper, simulated annealing and ant colony optimization approach are used for location and routing phases, separately. This problem have also been solved by using hybrid meta-heuristic which is based on simulated annealing and tabu search [31] and lagrange relaxation [32]. Different types of LRP have also been studied in the literature. Chan and Baker [33] have considered the distance-constrained LRP where customers can be served from multiple route/depot. The authors have proposed a mathematical model and a heuristic approach for the problem. Albareda-Sambola et al.[34] have consider the LRP where one route for each opened depot. The author has developed a mathematical model, lower bound evaluation approach and a tabu search approach for the problem. More general case of LRP has also been considered in the literature. Prins et al. [35] have studied the both depot and route capacitated LRP (CLRP), proposed a mathematical model and a hybrid heuristic approach which includes randomized saving algorithm, learning process, local search and path relinking operations. The same problem has been solved by using memetic algorithm which consists of genetic algorithm and local search procedure [36]. Moreover, a distance function which evaluates the similarity between individual pairs has been developed to diversify the solutions in the population. Prins et al. [37] have studied CLRP by using a two-phase algorithm which exchanges information between the location and routing phases. In the location phase, the routes including their customers are combined into super customers, and the corresponding problem is solved by using a lagrangean relaxation. In the routing phase, a granular tabu search procedure was used to solve the resulting multi depot vehicle routing problem. Barreto et al. [38] have developed a heuristic procedure based on cluster first route second approach for the same problem. In this heuristic several proximity measures have been developed. Proximity measures are used for grouping customer and evaluate the similarity of the customers in each group. Firstly, customers are grouped using proximity measures and feasible routes are constructed, then these routes are improved using 5

3-opt search mechanism, and finally these routes considered as a single super-node and single source facility location problem is solved which leads to a feasible CLRP solution. Marinakis and Marinaki [39] have solved the problem by using hybrid meta-heuristic solution procedure which is a particle swarm optimization hybridized with greedy randomized adaptive search procedure, expanding neighborhood search and path relinking. In a similar work, Marinakis et al. [40] have used honey bees mating optimization algorithm instead of particle swarm optimization. In another interesting paper, Duhamel et al. [41] have used evolutionary local search hybridized with greedy randomized adaptive search procedure. In that paper, permutation encoding has been used to represent the solution, special splitting algorithm developed to obtain feasible CLRP solutions and greedy procedure used to obtain initial population of evolutionary algorithm. The authors have reported that the approach outperforms all previously published methods and provides numerous new best solutions. This problem has also been solved by using an adaptive large neighborhood search heuristic including new neighborhood operators and adapted operators from the literature [42], multiple ant colony optimization [43], a modified granular tabu search meta-heuristic which considers five granular neighborhoods, three different diversification strategies and a perturbation procedure [44], and a heuristic procedure including greedy randomized local search procedure followed by local search procedures, an integer-linear programming formulation and column generation procedure [45]. Contardo et al. [46] have considered an exciting version of CLRP called as two echelon CLPR. The authors have developed an exponential size mathematical formulation, several families of valid inequalities and a branch-and-cut algorithm which can solve the test instances to optimality including 50 customers.The authors have also introduced an adaptive large neighborhood search meta-heuristic algorithm to obtain good quality upper bounds for large size problems containing up to 200 customers. A very similar problem has been solved by Nguyen et al. [47] using five different heuristic methods in which four constructive heuristic and a hybrid meta-heuristic including greedy randomized adaptive search procedure, learning process and path relinking. The authors have reported that the meta-heuristic outperforms the constructive heuristics. Toyoglu et al. [48] have introduced more complex of CLRP including several levels (i.e. depots, fixed/mobile transfer points and customers), multiproduct, heterogeneous vehicle fleet and time window restriction for customers. The authors have developed a mathematical formulation and several valid inequalities for the problem and investigate the performance on the test instances. The same problem has also been considered by Toyoglu et al. [49]. In that 6

paper, a new mathematical formulation (including fewer constraints and variables), valid inequalities and a clustering based heuristic solution procedure have been developed. The authors have showed that new formulation gives better bounds than the previous one and the heuristic solution procedure gives better upper bounds than the formulation for all test problems. Despite the fact that CLRP and VRPMB are two important problems for scientific literature and practice, CLRPMB has not received appreciable attention from researchers so far. In our previous study [50], we developed two mixed integer programming (MIP) formulations, which are node-based and flow-based formulations, for the problem and presented several polynomial-size valid inequalities adapted from literature to strengthen the formulations. Additionally, we proposed a heuristic approach consisting of two-phases (i.e. location and routing) to solve CLRPMB instances. Our experimental studies revealed that the flow-based formulation with valid inequalities, i.e. strong formulation, solves small-size CLRPMB instances to optimality within two hours of computation time by CPLEX, and the proposed heuristic is a computationally efficient approach to solve CLRPMB. Moreover, motivated by the success of flow-based formulation, we proposed a branch-and-cut algorithm using this formulation [51]. Computational results showed that optimal solutions obtained for some instances with up to 88 customers and 8 potential depots, in a reasonable computation time. 3. Problem Definition and Mathematical Formulation The distribution network design problem under investigation can be formally defined as follows; let G(N,A) be a complete directed network where N=N0∪NC is a set of nodes in which N0 and NC represent the potential depot nodes and customers, respectively.

A = {(i, j ) : i, j ∈ N} is the set of arcs, and to each arc (i,j) is associated a nonnegative cost (distance) cij, with cij = cji for each i,j∈N and triangular inequality holds (i.e., cij+ cjk≥ cik). A capacity CDk and a fixed opening cost FDk are associated with each possible depot node k∈N0. Each customer has delivery demands (di) coming from one selected depot and pickup demands (pi) going back to the same depot. There are unlimited number of homogeneous vehicles with known capacities (CV) and fixed operating cost (FV) at depots. The problem is to simultaneously determine the locations of depots, assignment of customers to opened depots and the corresponding delivery routes with minimum total cost subject to the following constraints: 7

•

Each vehicle can perform at most one route.

•

Each customer must be served by exactly one vehicle.

•

Each route must begin and end at the same depot,

•

The total load on vehicles must not exceed the vehicle capacity

•

The total pickup and total delivery load of the customers assigned to an opened depot must not exceed the capacity of the depot.

We utilize following decision variables in order to formulate CLRPMB. ⎧1 if a vehicle travels directly from node i to node j (∀i,j ∈ N ) xij = ⎨ ⎩0 otherwise

⎧1 if depot k is opened (∀k ∈ N 0 ) yk = ⎨ ⎩0 otherwise ⎧1 if customer i is assigned to depot k (∀i ∈ N C ,∀k ∈ N 0 ) zik = ⎨ ⎩0 otherwise

U ij ; demand to be delivered to customers routed after node i and transported in arc (i, j) if a vehicle travels directly from node i to node j (∀i, j ∈ N ) , otherwise 0 Vij ; demand to be picked-up from customers routed up to node i (including node i) and

transported in arc (i, j) if a vehicle travels directly from node i to node j (∀i, j ∈ N ) , otherwise 0 The flow-based formulation of CLRPMB is given as follows:

min

∑ ∑ c x + ∑ FD y + ∑ ∑ FVx i∈N j∈N

s.t.

∑x j∈N

ij

∑x j∈N

ji

∑U j∈N

ji

ij ij

ij

k

k

k∈N 0 i∈NC

(1)

ki

=1

∀i ∈ NC

(2)

= ∑ xij

∀i ∈ N

(3)

− ∑ U ij = di

∀i ∈ NC

(4)

∀i ∈ NC

(5)

∀i, j ∈ N ,i ≠ j

(6)

j∈N

j∈N

∑V − ∑V j∈N

k∈N 0

j∈N

ji

= pi

U ij + Vij ≤ CVxij

8

∑U

kj

=

∑z

∑U

jk

=0

∀k ∈ N 0

(7)

∀k ∈ N 0

(8)

∀k ∈ N 0

(9)

∀k ∈ N 0

(10)

U ij ≤ (CV − di ) xij

∀i ∈ NC , ∀j ∈ N

(11)

Vij ≤ (CV − p j ) xij

∀i ∈ N , ∀j ∈ NC

(12)

U ij ≥ d j xij

∀i ∈ N , ∀j ∈ NC

(13)

Vij ≥ pi xij

∀i ∈ NC , ∀j ∈ N

(14)

∑z

∀i ∈ NC

(15)

j∈NC

j∈NC

∑V

j∈NC

=

kj

=0

j∈NC

∑V

j∈NC

k∈N 0

∑z

jk

j∈NC

ik

jk

jk

dj

pj

=1

∑dz

≤ CDk yk

∀k ∈ N 0

(16)

∑ pz

≤ CDk yk

∀k ∈ N 0

(17)

xik ≤ zik

∀i ∈ NC , ∀k ∈ N 0

(18)

xki ≤ zik

∀i ∈ NC , ∀k ∈ N 0

(19)

∀i, j ∈ NC ,i ≠ j , ∀k ∈ N 0

(20)

xij ∈ {0 ,1}

∀i, j ∈ N

(21)

zik ∈ {0,1}

∀i ∈ NC , ∀k ∈ N 0

(22)

yk ∈ {0,1}

∀k ∈ N 0

(23)

U ij ,Vij ≥ 0

∀i, j ∈ N

(24)

i∈NC

i∈NC

i ik

i ik

xij + zik +

∑

m∈N0 ,m ≠ k

z jm ≤ 2

where xij is set to zero when max {d i + d j ; pi + p j ;d j + pi} > CV , ∀i, j ∈ NC ,i ≠ j . This restriction guarantees that any incompatible customer pair whose total demands are greater than the vehicle capacity is not appeared in the same route. In this formulation, objective function (1) minimizes the total system cost including transportation, depot opening and vehicle operating costs. Constraint sets (2) and (3) are known as degree constraints. While constraint set (2) ensures that each customer must be visited exactly once, constraint set (3)implies that entering and leaving arcs to each node are equal. Constraint sets (4) and (5) are

9

flow conservation constraints for delivery and pickup demands, respectively.

These

constraints eliminate subtours and satisfy the pickup and delivery demands of each customer. The role of constraint set (6) is to prevent the total load on any arc exceeding the vehicle capacity. While constraint set (7) ensures that the total delivery load dispatching from each depot equals to the total delivery demand of customers which are assigned to the corresponding depot, constraint set (8) provides that the total amount of delivery load returning to the depots must be equal to zero. Similarly, while constraint set (9) ensures that the total pickup load entering to each depot equals to the total pickup demand of customers which are assigned to the corresponding depot, constraint set (10) provides that the total amount of pickup load dispatching from the depots must be equal to zero. Constraint sets (11)-(14) are bounding constraints for additional variables. Constraint set (15) assigns each customer to only one depot. Constraint sets (16) and (17) avoid the total delivery and pickup loads on any depot exceeding the depot capacity, respectively. Constraint sets (18)-(20) forbid the illegal routes, i.e. routes which do not start and end at the same depot. Finally, constraint sets (21)-(24) are known as integrality constraints which define the nature of the decision variables. This formulation includes O(|NC|+|N0|2) binary variables, O(|NC|+|N0|2) additional variables and O(|NC|2|N0|) constraints. 4. Memetic Algorithm

Since CLRPMB is an NP-hard problem, MIP formulations are not directly applicable to find optimal solutions for large-size problems. Hence, heuristics are necessary to quickly obtain solutions for these problems. In this paper, we propose a memetic algorithm (MA) based on genetic algorithm (GA), simulated annealing (SA) and integer programming formulation (IPF) to solve the CLRPMB instances. GA refers to a class of adaptive search procedures based on the principles derived from natural evolution and genetics. GA offer significant advantages over conventional methods by using simultaneously several search principles and heuristics. The most important ones include a population-wide search, a continuous balance between convergence and diversity, and the principles of building-block combination. GA can be implemented in several different ways to solve any problem [52]. In our implementation, we use steady-state replacement strategy in which only one new solution is obtained by genetic operators and inserted to the current population, and then the worst chromosome is removed from population. For further improve the solution after each recombination procedure of GA, SA and IPF are used as a 10

local search algorithm recursively. This section presents the detailed description of components of the proposed algorithm, MA, for CLRPMB.

4.1. Representation Representation is one of the important issues that affect the performance of GA. Usually different problems have different data structures or genetic representations. In this paper, we implement encoding scheme originally developed by Prins et al. [36] for CLRP. In this encoding scheme, a chromosome consists of two parts: i) depot status and ii) customer sequence. The depot status is a vector of |N0| potential depots. Each gene in the vector represents the status of corresponding depot, i.e. closed or opened. If the depot i is closed, its gene takes the value of zero, otherwise it has an integer number that represents the index in customer sequence of the first customer assigned to that depot. The customer sequence is a permutation of |NC| customers in which lists of customers assigned to opened depots are connected in series. The decoding of a chromosome includes: i) determining the opened depots using depot status part and ii) assigment of customers to the opened depots considering depot and customers status and also depot capacities iii) obtaining optimal routes for each opened depot and its customers by using partitioning procedure. The detailed steps of decoding is given in Figure 1. In Step 1, assignment lists and auxiliary graph are initialized. By using both parts of chromosome, depots are opened and customers are assigned to opened depots without considering the depot capacities (Step 2 and 3). It should be noted that all customers must be assigned to any opened depot. Therefore, depot status part must include index “1”. Step 2 checks this necessity and repairs if it does not include. If assignment of customers to the opened depots is failed because of the depot capacity restrictions, last customer is removed from the depot and assigned to a depot which has available capacity (Step 4). After feasible assignment of customers to the opened depots, CLRPMB is reduced to independent VRPMBs. Finally optimal routes for each VRPMB are obtained using modified version of route partitioning procedure proposed by Prins[53] for VRP (Step 5). The evaluation of a given partial solution Ωk = {c1 ,c2 ,… ,cs} where s is the number of customers assigned to opened depot k, is concerned with partitioning customer orders along Ωk into routes so that the load on any arc of each route does not exceed vehicle capacity, while the total cost is minimized. An auxiliary graph is constructed as follows: let G ( Ν k ,Αk ) be a directed

acyclic

graph

where

N k = {k ∪ Ωk }

is

the

set

of

vertices

and

11

Ak = {(ci ,c j ) : ci ,c j ∈ N k ;i < j} is the set of directed arcs. Each arc (ci ,c j ) ∈ Ak represents a

feasible route where the vehicle departs from depot (node k) and visits nodes ci +1 ,ci + 2 ,...,c j −1 , and c j , consecutively, and turns back to the depot, such that load on any arc does not exceed vehicle capacity. After the construction of auxiliary graph, a minimum costpath on G from node k to last node (cs) defines the optimal partition of the sequence which leads to a feasible VRPMB solution originating from depot k. Figure 1 here An illustrative example for a chromosome is given in Figure 2. In this example, CLRPMB consist of 3 potential depots and 11 customers. When the depot status part of the chromosome is considered, it is seen that depots 1 and 3 are opened while depot 2 is closed (Figure 2a). The indexes of first customers assigned the depots 1 and 3 are 7 and 1, respectively. It means that the customers indexed between 1 and 6 in the customer sequence, i.e. (7-6-10-11-4-8), are assigned to depot 3 and the customers indexed between 7 and 11, i.e. (9-2-1-3-5), are assigned to depot 1 (Figure 2b). After these assignments, an acyclic graph is constructed (Figure 2c) and optimal route for each opened depot is obtained by using adapted route partitioning procedure briefly explained above (Figure 2d). In this example, vehicle capacity is set to 50. Figure 2 here

4.2. Initialization of Population In classical GA, the initial population is usually generated randomly. However, selection a suitable initial population accelerates the convergence of GA. To reach optimal/near optimal solution for the problem, we generate an initial population including heuristic solutions as well as random solutions to explore the different regions of solution space. While the first half of the initial population consists of heuristic solutions obtained by a greedy heuristic approach called Extended Clarke and Wright Algorithm (ECWA), the second half includes randomly generated solutions. ECWA is an extension of the well-known Clarke and Wright heuristic [29] and firstly introduced by Prins et al. [36] for CLRP. The steps of the algorithm for generating an initial feasible solution are given in Figure 3. Firstly, assignment lists are initialized (Step 1). Then, a penalty cost ( pci ) arising when the customer is assigned to the second closest depot instead of the closest one is calculated for each customer, and the 12

customers are listed in non-increasing order of their penalty costs (Step 2). In Step 3, each customer is assigned to the closest depot starting from top of the list. If some of the customers are not assigned to their closest depots because of the depots’ capacity, their penalty costs are recalculated considering closest ones with enough capacity and the assignment of these customers to appropriate depots is tried again. These steps, i.e. calculating penalty costs and assignment of customers to closest depots with enough capacity, are repeated until all customers are assigned. Step 5 includes the construction of vehicle routes for each opened depot. Depots without customers are closed and a simple route for each customer from whose depot is constructed. Then, each pair of routes (R1, R2) in the feasible solution are examined in terms of cost saving obtained by their combination and a pair of nodes providing largest cost saving is combined (see merge operator explained in moving strategies). This strategy is repeated until no capacity-feasible combination is found. In order to generate the different good solutions for the initial population, randomized version of ECWA (RECWA) is applied. In this procedure route combination alternatives are listed in non-increasing order according to their cost savings and one of them is selected randomly from the first N0 alternatives in this list. Figure 3 here

4.3. Crossover and Mutation Operators The crossover and mutation operators are used to explore new solutions in the search space. In general, crossover operator is applied to obtain new solutions by exchanging some information between selected parents. The chromosome in MA composes of two different encoding types, i.e. integer number encoding in depot status and permutation encoding in customers’ sequence. Therefore, different crossover operators should be used for each part of chromosome. After preliminary experiments, we utilize one-point crossover for the depot status and order crossover (OX) for the customer sequence. One-point crossover is a wellknown classical genetic operator in which a cut point is selected randomly between 1 and

N 0 − 1 , and then the depot status of the offspring is obtained by taking left and right parts of the cut point from the first and second parents, respectively. In OX, two different crossing points are randomly selected. These two points enclose a section of one or more genes of a chromosome. This section of the first parent is directly transferred into the offspring without changing the index of each gene. Then, the offspring is filled up by copying the elements of the second parent, without considering the genes copied from first parent. It is worthy to note 13

that the parents undergoing crossover operation in each generation of MA are selected by binary tournament selection strategy. Unlike the crossover, mutation is usually done by modifying gene within a chromosome. After preliminary experiments, we implement classical swap mutator for the second part of chromosome. In this operator, randomly selected two customers are exchanged. An illustrative example for crossover and mutation operators is given in Figure 4. Figure 4 here Crossover and mutation operators may lead to infeasible solution. Therefore, after each recombination operation, the feasibility check should be performed. Firstly, it must be checked that all the customers are assigned to any opened depot (i.e. depot status includes the index 1). If it fails, then the first closed depot is opened and the depot status of this depot is changed from 0 to 1. Secondly, capacity violation of each opened depot must be checked (i.e. total delivery or pickup demand of customers assigned to an opened depot is greater than the capacity of the depot). If a capacity violation exists for any opened depot, then the last assigned customer is removed from the corresponding depot and assigned to the first opened depot having enough capacity. If none of the depots have enough capacity to assign the corresponding customer, one of the closest depots is opened randomly.

4.4. Local Search In each generation of MA, we implement a heuristic procedure as a Local Search (LS) algorithm to improve the offspring obtained by genetic operators. LS algorithm includes two main phase called “Improving the Solution” and “Reconfiguration of Current Assignment of Routes”. While the first phase uses Simulated Annealing (SA) procedure to improve the solution, the second phase converts the current feasible solution to the equivalent Facility Location Problem (FLP) and solves it using Integer Programming Formulation (IPF). The detailed descriptions of these phases are given below. It should be noted that these steps are also implemented to the feasible CLRPMB solution obtained by RECWA or decoding of the chromosome.

4.4.1. Improving the Solution We implement SA algorithm to improve the solution in the first phase. SA, which stems from the simulation of the annealing of solids, is a stochastic search technique that is able to 14

escape local optima using a probability function [54]. It starts with an initial feasible solution and improves it until stopping criterion is met. In our SA algorithm neighborhood operators of CLRPs (i.e. swap, insert, move etc.) are performed. For the CLRP the feasibility of movement can easily be evaluated before examining by considering the statistics of the current solution, such as the total load of routes and depots, demand of moving customer etc. However, because of the pickup and delivery structure of the problem, we need to additional data structure and a special evaluation scheme for CLRPMB to speed up SA algorithm. For this purpose, we used following notations and definitions in the remainder. A feasible CLRPMB solution includes a set of routes (R) and each route r (∀r ∈ R) starts from an opened depot (i.e. depot k), visits customers and turns back to the depot k. Each route is represented by an ordered list of NCr customers and two copies of depot k. Additionally, each customer is depicted with ri: ith customer in route r, and ti: position of customer i. Pickup and delivery demand of each route are calculated as the sum of pickup and delivery demands of customers (i.e. Dr = ∑ i =1 r d ri and Pr = ∑ i =1 r pri , respectively). Moreover, the pickup and NC

NC

delivery demand of each depot k are calculated as the sum of pickup and delivery demand of routes assigned to depot k, respectively (i.e. DDk = ∑ r∈R ark Dr and DPk = ∑ r∈R ark Pr , where

ark is an indicator variable and ark = 1 if route r is assigned to depot k, otherwise ark = 0 ). Furthermore, to model the load fluctuations on a vehicle, we use following metrices: DL+rt (

PL+rt ) is the total delivery (pickup) load on vehicle of route r just after node positioned at t, DL−rt ( PL−rt ) is the total delivery (pickup) load on vehicle of route r just before node positioned at t, ML+rt ( ML−rt ) is the maximum load on vehicle of route r after (before) node positioned at t. In each iteration of SA procedure, some of these metrics which are affected from movement are updated in a constant time. SA algorithm uses following neighborhood operators which use the same data structure as explained above to check the feasibility of movement as fast as possible. RouteMove: One randomly selected route (r) is removed from its current depot and assigned

to the best opened depot (k) which having enough remaining capacity. The moving of route r to depot k is feasible if the following condition is satisfied. - max {DDk + Dr ,DPk + Pr } ≤ CDk

(25)

15

In this operator, substantial saving can be obtained (i.e. saving may includes depot opening cost if there is only one route of depot k). An illustrative example for this operator can be seen in Figure 5. In this figure route (2-7-6-2) is removed from depot 2 and assigned to depot 3 (Figure 5a). After this assignment, depot 2 is closed, since there is no other route assigned to this depot (Figure 5b). Figure 5 here CustMove: One randomly selected customer (i) is removed from its current position and

inserted to the best position (tB) of the same route (r1) or different route (r2) belonging to the same depot(k1) or different depot (k2). This operation is feasible under following conditions. - if k1 ≠ k2 which means that routes r1 and r2 belongs to different depots, then it must be checked the feasibility of the depot load and the load fluctuation of inserted route r2. max {DDk2 + di ,DPk2 + pi} ≤ CDk2

(26)

max {ML−r2tB + d i ,ML+r2 ,tB −1 + pi} ≤ CV

(27)

- if k1 = k2 and r1 ≠ r2 which means that routes r1 and r2 belongs to the same depot but different routes, then it is enough to check the feasibility of load fluctuation of the inserted route r2 (inequality (27)). - if k1 = k2 and r1 = r2 = r which means that the selected customer is moved into the same route, thenthe load fluctuation of each position between the current and the new position must be checked as follows.

if ti t B

DL+rt + PL+rt + di − pi ≤ CV DL−rt + PL−rt − di + pi ≤ CV

t = ti + 1,ti + 2 ,...,t B t = t B ,t B + 1,...,ti − 1

(28)

In this operator, substantial saving can be obtained (i.e. saving may include vehicle fixed cost if there is only one customer in route r1, and depot opening cost if there is only one route in depot k1). An illustrative example for this operator can be seen in Figure 6. In this figure customer 11 is removed from its current position (Figure 6a) and inserted to the third position of route (3-10-4-8-3) (Figure 6b). Figure 6 here

16

Swap: Two randomly selected customers (i and j), which are in the same route (r1) or in two

different routes (r2) belonging to the same depot (k1) or two different depots (k1 and k2), are exchanged. The swap operation is feasible under following conditions. - if k1 ≠ k2 which means that routes r1 and r2 belongs to different depots, then it must be checked the feasibility of the load variation of both k1 and k2, and the load fluctuation of

r1 and r2. max {DDk1 + d j − di ,DPk1 + p j − pi} ≤ CDk1

(29)

max {DDk2 + d i − d j ,DPk2 + pi − p j} ≤ CDk2

(30)

max {ML−r1ti − d i + d j ,ML+r1ti − pi + p j} ≤ CV

(31)

{

}

max ML−r2t j − d j + d i ,ML+r2t j − p j + pi ≤ CV

(32)

- if k1 = k2 and r1 ≠ r2 which means that routes r1 and r2 belongs to the same depot but different routes, then it is enough to check the feasibility of the load fluctuation of r1 and

r2 (inequality (31) and (32)). - if k1 = k2 and r1 = r2 = r which means that selected customers areswapped in the same route, then the load fluctuation of each position between the current and the new position must be checked as follows.

if ti t j

DL+rt + PL+rt + di − pi − d j + p j ≤ CV DL+rt + PL+rt + d j − p j − di + pi ≤ CV

t = ti ,ti + 1,...,t j − 1 t = t j ,t j + 1,...,ti − 1

(33)

An illustrative example for swap operator can be seen in Figure 7. In this figure, customers 11 and 2 are removed from their current positions (Figure 7a) and swapped (Figure 7b). Figure 7 here Merge: Two routes (r1 and r2) belonging to the same or different depots are selected and

merged into the new one (r=r1∪r2) considering following alternatives: new route is assigned to the depot of first (k1) or second route (k2), or to the another depot (k) opened or not. Hence, after these alternatives are evaluated in each merge operation, the one with largest cost saving is applied. Merge operation is feasible under following simple depot and route conditions. max {DDk + Dr1 + Dr2 ,DPk + Pr1 + Pr2 } ≤ CDk max {DDk + Dr2 ,DPk + Pr2 } ≤ CDk max {DDk + Dr1 ,DPk + Pr1 } ≤ CDk

if k ≠ k1 and k ≠ k2 if k = k1 and k ≠ k2

(34)

if k ≠ k1 and k = k2

17

{

}

max ML−r1NCr + Dr2 ,ML+r2 1 + Pr1 ≤ CV 1

(35)

In this operator, considerable saving can be obtained (i.e. saving exactly includes vehicle fixed cost, and may include depot opening cost if there is only one route of depots k1 and/or

k2). An illustrative example for this operator can be seen in Figure 8. In this figure routes (310-11-3) and (2-4-8-2) are merged into new one and connected to opened depot 3 as (3-10-114-8-3). Figure 8 here Opt: Two non-consecutive arcs (i-j and m-n), which are in the same route (r1) or in two

different routes (r1 and r2) belonging to the same depot (k1) or different depots (k1 and k2), are deleted. If the deleted arcs are in the same route ( r1 = r2 = r ), then two new arcs are created (i-

m and j-n) and the path lying between the created arc pair is reversed. If deleted arcs are in two different routes of the same depot ( r1 ≠ r2 and k1 = k2 = k ), each route is divided into two parts as starting and terminating part. Then two new arcs are created in such way that starting and terminating parts of two different routes are connected ( k ,...,i,n,...,k and k ,...,m, j,...,k ). Otherwise, i.e. the deleted arcs belong to two different routes of different depots ( r1 ≠ r2 and

k1 ≠ k2 ), after two different parts coming from different routes are connected two new arcs, they are revised such that each route starts and finishes at the same depot ( k1 ,...,i,n,...,k1 and

k2 ,...,m, j,...,k2 ). Opt operation is feasible under following conditions. - if k1 ≠ k2 which means that routes r1 and r2 belongs to different depots, then it must be checked that the feasibility of the load variation of depots k1 and k2, and the load fluctuation of routes r1 and r2 which are as given in (36)-(39). max {DDk1 + DL−r2tn − DL+r1ti ,DPk1 + ( Pr2 − PL−r2tn ) − ( Pr1 − PL+r1ti )} ≤ CDk1

{ max {ML max {ML

(

}

)

max DDk2 + DL−r1t j − DL+r2tm ,DPk2 + Pr1 − PL−r1t j − ( Pr2 − PL+r2tm ) ≤ CDk1

}

− r1t j

− DL−r1t j + DL−r2tn ,ML+r1ti − PL+r1ti + PL+r2tm ≤ CV

− r2 tn

− DL−r2tn + DL−r1t j ,ML+r2tm − PL+r2tm + PL+r1ti ≤ CV

}

(36) (37) (38) (39)

- if k1 = k2 and r1 ≠ r2 which means that routes r1 and r2 belong to the same depot, thenit is enough to check the feasibility of the load fluctuation of r1 and r2 (inequality (38) and (39)). 18

- if k1 = k2 and r1 = r2 = r which means that selected customer is moved in the same route, then the load fluctuation of the vehicle between selected arcs must be checked as follows. DL+rti + PL+rti + ∑ sm=t t

d rs − ∑ sm=t t

m −t

m −t

prs ≤ CV

t = 0 ,1,...,tm − t j

(40)

An example for opt operator can be seen in Figure 9. In this figure arcs (9-4) and (11-2) are removed (Figure 9a) and new arcs (9-2) and (11-4) are added (Figure 9b). Please note that the depot connection of the last customers is also changed (i.e. arcs (8-1) and (1-3) are replaced with (8-3) and (1-1), respectively). Figure 9 here In our SA algorithm, best improving strategy is used. According to this strategy, all moving operators are called and the best one among the neighbors is chosen as a new solution

(Snew ) for the problem. Moreover, SA algorithm implements a candidate list strategy in generating the neighbors since searching whole neighborhood of the current solution by a moving strategy is a very time consuming process. According to this strategy, each moving strategy generates randomly a subset of neighbors (NNx) that satisfying vehicle and depot capacity constraints and they are gathered in a pool to select a new solution. If the new solution is better than current solution ( Scur ) then it is accepted as a current solution, otherwise it is accepted with probability of exp (−Δs / Titr ) as a current solution. Δs is a relative percent deviation of quality of the new solution from the current solution and calculated by

[( f ( S new ) − f ( Scur )) / f ( Scur )]* 100 .

In each iteration of SA, temperature (Titr ) is reduced

using a geometric cooling schedule, i.e. Titr = αTitr −1 where α is a cooling rate. SA stops when the temperature reaches to final temperature.

4.4.2. Reconfiguration of Current Assignment of Routes In this phase, we convert the current CLRPMB solution to the equivalent Facility Location Problem (FLP) without changing the structure of routes. Then we solved the converted problem using general purpose mixed integer programming solver (i.e. CPLEX). Although the FLP is an NP-hard problem, optimal solution can be obtained in a very short computation time for the small-size instances. The detailed description of conversion procedure is as follows: 19

First and last arcs of each route, which connect the first and the last customers to depot

k, are removed and customers are aggregated to a single supernode (r). The assignment cost of each supernode r to any depot k ( Crk ) includes the vehicle fixed cost (FV), transportation cost between customers in r, the transportation cost from depot k to the first customer of r and the transportation

cost

from

the

last

customer

of

r

to

depot

k

(i.e.

Crk = FV + ckr1 + ∑ i =1 r cri ri+1 + crNC k ). After this transformation, it is determined the locations NC −1

r

of depots and assignment of supernodes to opened depots, such that each supernode can be assigned only one depot, and the total pickup and delivery load of the supernodes assigned to a depot must not exceed the capacity of the depot. Under these assumptions, the decision variables and the mathematical model of FLP is as follows:

⎧1 if supernode r is assigned to depot k (∀r ∈ R,∀k ∈ N 0 ) ⎩0 otherwise

ξ rk = ⎨

⎧1 if depot k is opened (∀k ∈ N 0 ) βk = ⎨ ⎩0 otherwise min ∑ ∑ Crk ξ rk + r∈R k∈N 0

s.t.

∑ξ

k∈N0

rk

∑ FD β

k∈N 0

k

(41)

k

=1

∀r ∈ R

(42)

∑D ξ

≤ CDk β k

∀k ∈ N 0

(43)

∑ Pξ

≤ CDk β k

∀k ∈ N 0

(44)

ξ rk ∈ {0,1}

∀r ∈ R , ∀k ∈ N 0

(45)

β k ∈ {0,1}

∀k ∈ N 0

(46)

r∈R

r∈R

r rk

r rk

In this formulation objective function (41) minimizes the total cost including supernode assignment and depot fixed costs. Constraints (42) ensure that each supernode must be assigned to only one depot. Constraints (43) and (44) guarantee that total delivery and pickup loads on any depot must not exceed the depot capacity, respectively. Finally, constraints (45) and (46) are integrality constraints. The contributions of this phase are twofold: i) improve the solution and ii) give us a new chance to improve the solution using the first phase. Please note that optimal solution of the given model is always same or better solution than corresponding CLRPMB solution. Therefore, we implement these two phases recursively, if better solution is obtained in the second phase. Steps of local search procedure is given in Figure 10. 20

Figure 10 here The overall pseudo-code procedure for solving CLRPMB is outlined in Figure 11. Figure 11 here 5. Computational Results

This section presents the results from our computational experiments which were organized into two stages. While the first stage investigates the performance of MA on CLRP test instances, the second stage evaluates its performance on CLRPMB test instances derived from the literature. MA was coded using C++ programming language and all experiments were performed on Intel Xeon 3.16 GHz equipped with 1 GB RAM computer. The state-of-the-art LP/MIP solver CPLEX (version10.2) was used to solve the formulations in local search. After some preliminary experiments, following parameter values were implemented in MA for solving both CLRP and CLRPMB test instances: population size was taken as |NC|+|N0|, crossover and mutation rates were set to 1 and 0.05, respectively, finally the number of generations was set to |NC|*|N0|. The initial temperature (T0) in SA was taken as 665 in which an inferior solution (inferior by 70% relative to current solution) is accepted with a probability of 0.90, the final temperature (Tf) was taken as 0.15 such that a solution which is inferior by 1% relative to current solution is accepted with a probability of 0.001. The number of neighbors generated by each moving strategy was set as NN rm = N 0 and NN cm = NN sw = NN me = NN opt = N C , and finally cooling rate was considered as 0.95. Because of the stochastic nature of MA, it was run 5 times for each instance with different random seeds.

5.1. Performance of MA on CLRP Instances Since CLRPMB is the general case of the CLRP, the developed algorithm can be used to solve CLRP instances by considering the original demand of each customer i as delivery demand ( di = qi ) and the pickup demand as zero ( pi = 0) . To investigate the performance of MA, we solved the original CLRP test instances of Prodhon [55] and compared the results with the upper bounds obtained by Contardo et al. [20] which are the best bounds found so far. Table 1 reports the results for MA on CLRP instance. In this table, the first column represents the name of CLRP instances and the second column gives the best upper bound

(UBCon )

reported in [20]. The text in bold characters indicates that the optimal solution is

obtained for the corresponding problem. Next three columns give the minimum, average and maximum percentage gap value between the best upper bound and the solution value of MA,

21

and calculated as 100* [( MA − UBCon ) / UBCon ] . Finally, the last column gives the average computation time per run over the five runs. Table 1 here The results for CLRP test instances are quite promising. Optimal solutions for ten problems are obtained in a very short computation time. Moreover, the best upper bound for one instance namely 100-10-1b is improved. Finally, average gap value between MA and the best upper bounds are 0.43%, and the average computation time is about 87 seconds. These results indicate that MA is a preferable approach for solving CLRP instances successfully.

5.2. Performance of MA on CLRPMB Instances Before the computation results, firstly we give a brief information about the test instances. In order to investigate the performance of MA, CLRPMB test instances were generated from the test set for CLRP presented by Barreto [56] and Prodhon [55]. Barreto’s set includes 20 test instances which were taken from CLRP literature or obtained by adding depots to classical VRP instances. While the number of customers varies between 8 and 318, the number of depots is between 2 and 15. We consider the names of original CLRP and VRP instances to denote an instance from this set in sequel. The name of each instance includes the information about the name of the author, who generated the instance, the publication year of the instance, the number of customers, |NC|, and the number of potential depots, |N0| (i.e.Author-Year-|NC|x|N0|). Prodhon’s set consists of 28 CLRP instances with capacitated routes and depots. The main characteristics of the set are given as follows: the number of customers |Nc| in {20, 50, 100, 200}, the number of potential depots |N0| in {5, 10}, uniform demands in [11, 20], the number of clusters clu in {1, 2, 3} (1 means that all nodes scatter on Euclidean plane), vehicle capacity CV in {a, b} where a = 70 and b = 150. In the set, depot capacities have been determined in such way that at least two or three depots are opened. In rest of the paper, we utilize following coding structure to denote an instance from this set: |NC|-|N0|-cluCV. In order to generate the delivery and pickup demands of the customers in each test instance, we utilized demand separation approaches of Salhi and Nagy [57] and Angelelli and Mansini[58] (denoted as SN and AM, respectively). These approaches are briefly defined as follows. In SN approach, a ratio ri = min ( xi / yi ; yi / xi ) , where xi and yi are the coordinates of customer i, is calculated for each customer i, and then the delivery and pickup demands are 22

generated as di = ri * qi and pi = qi − di , where qi is the original demand of customer i. We refer this type of problems as X. Similarly, another type of problem, called Y, is generated by exchanging delivery and pickup demands of each customer. In AM approach, the original demand of each customer i is considered as delivery demand ( di = qi ) and the pickup demand of the corresponding customer is generated as pi = ⎣⎢(1 − γ ) qi ⎦⎥ if i is even and pi = ⎣⎢(1 + γ ) qi ⎦⎥ if i is odd. In this paper, we consider two γ values as 0.2 and 0.8 to generate two different types of problems called Z and W, respectively. We take into account the test problems including up to 100 customers as in [51] to compare with our results. As a result, 148 test instances were generated by using 37 original CLRP test problem and 4 different separation strategies (X, Y, Z and W type). In order to investigate the performance of MA, the results are compared with the lower bounds obtained by Karaoglan et al. [51] on the test instances of Barreto [56] and Prodhon [55]. In the comparison, following performance measures are used: i) Average solution time (CPU), ii) Percentage gap (gap): it is calculated as 100* [(UB − LB) / LB] where LB is the lower bound obtained by Karaoglan et al. [51] and upper bound (UB) is the solution obtained by MA, iii)Robustness Index (RI): this index is calculated as 100* [(UBwost − UBbest ) / UBwost ] where UBwost and UBbest are the worst and best solution obtained from five runs of MA, respectively. Tables 2 – Table 5 present computational results for MA on each CLRPMB test set, respectively. In these tables, while the first column represents the name of CLRPMB instances, the second column gives the demand separation strategy (DSS). Next column shows the percentage gap value of branch-and-cut (B&C) algorithm given in [51]. Subsequent five columns report the minimum, average and maximum values of percentage gap, average computation time and robustness index (RI) value over the five runs obtained by MA. To investigate the effect of the integer programming formulation in local search procedure, we compare our algorithm with MA without IPF which is obtained by removing reconfiguration steps from local search. The results are reported in last five columns of these tables. Table 2 and Table 5 here As seen from tables, good quality solutions are obtained by MA in a short computation time. On small-size problems (test problems up to 36 customer in Barreto’s and 20 customers 23

in Prodhon’s test set), two problems with 88 customers in Barreto’s test set, and two problems with 50 customers in Prodhon’s test set, optimal solutions are obtained in a very short computation time (maximum solution time for these problems is 36.27 seconds). Moreover, satisfactory results are obtained for other test problems. While average gaps of MA are 1.42% and 6.51% for the SN and AM demand separation procedures for Barreto’s test set, these values reduce to 0.15% and 1.32% for Prodhon’s test set. According to these results, good quality solutions for Barreto’s and Prodhon’s test sets are obtained in a short computation time. Moreover, instances obtained by using AM approach are more difficult to solve than those derived from same test sets by SN approach. While average gaps are 1.42% and 0.15% for test problems with SN approach, respectively, this value increases to 6.51% and 1.32% for test problems with AM approach, respectively. Similarly, CPU times also support these results. While average CPU times 45.27 and 95.03 seconds for test problems with SN approach, respectively, this value increases to 77.72 and 141.23 seconds for test problems with AM approach, respectively. Results show also that the removing the reconfiguration step of local search procedure leads to worse average gaps (on average 3.43% to 8.16% for Barreto’s, 0.75% to 2.64% for Prodhon’s test set), reduction on the number of instances solved to optimality (49 instead of 56) over 148 instances, and reduction on solution times from 95.17 seconds to 56.57 seconds. Robustness Index (RI) is another important criterion to evaluate the behavior of heuristic algorithms in stochastic nature. It is well known that most of the heuristic algorithms work with stochastic numbers and different solutions can be obtained with different seeds. In these algorithms, it is expected to obtain optimal/near optimal solutions which are close to each other (robust solutions). Therefore we can evaluate the robustness of the algorithm by using RI, where the lower this is, the better the robustness. In our experiments, therobust solutions are obtained by using MA (average robustness indexes are 0.16 for Barreto’s and 0.06 for Prodhon’s test set). These results are an indicator that MA is robust to initial solutions. Many memetic algorithms in the literature do not use mutation operator where the local search procedure play the role of a mutation.We also analyze this special structure of MA and present the results in Table 6. In this table, while the first column represents the name of CLRP instance sets, the second column gives the demand separation strategy (DSS). Subsequent five columns report the minimum, average and maximum values of percentage

24

gap, average computation time and robustness index (RI) value over the five runs obtained by MA and MA without mutation, respectively. Table 6 here These results show that the mutation operator for the algorithm slightly improves the performance of the algorithm. The average gap of MA without mutation is 2.57% this value reduces to 2.35% for the MA with the cost of tolerable increase on solution time (44.16 seconds versus 89.81 seconds). 6. Conclusion

In this paper, we considered CLRPMB and developed a memetic algorithm to solve the problem. We investigated the performance of the algorithm on both CLRP and CLRPMB test instances derived from the literature, and evaluate its performance using lower and upper bounds found in the literature. Computational results indicated that MA is a viable approach for solving CLRP instances. We obtained feasible solutions with the average percentage gap of 0.43, and improved the best known feasible solution for an instance 100-10-1b. Additionally, for CLRPMB instances, good quality solutions (average gaps are 3.97% and 0.74%) are obtained by the algorithm in a reasonable computation time (average CPU times are 61.49 and 118.13 seconds). Further research can be performed on the new crossover and mutation operator for MA and new moving strategies for SA can be introduced. Finally, proposed algorithm can be used to obtain an initial solution for any exact algorithm (i.e. branch-andbound, branch-and-cut) to shorten the optimization process of exact algorithm. Acknowledgement

We thank to the Associate Editor and anonymous referee for their helpful comments and suggestions, which greatly contributed to the improvement of the paper. This research is fully supported by The Scientific and Technological Research Council of Turkey as Scientific Research Project (No: 108E069).

25

References

[1] Salhi S, Rand GK. The effect of ignoring routes when locating depots. European Journal of Operational Research. 1989;39(2):150-6. [2] Laporte G, Nobert Y, Arpin D. An exact algorithm for solving a capacitated locationrouting problem. Annals of Operations Research. 1986;6:293-310. [3] Berger RT, Coullard CR, Daskin MS. Location-Routing Problems with Distance Constraints. Transportation Science. 2007;41(1):29-43. [4] Belenguer JM, Benavent E, Prins C, Prodhon C, Calvo RW. A Branch-and-Cut method for the Capacitated Location-Routing Problem. Computers & Operations Research. 2011;38(6):931-41. [5] Min H, Jayaraman V, Srivastava R. Combined location-routing problems:A synthesis and future research directions. European Journal of Operational Research. 1998;108:1-15. [6] Laporte G. Location-routing problems. In: Golden BL, Assad AA, (Eds.) Vehicle Routing: Methods and Studies. Amsterdam: North-Holland; 1988. p. 163-98. [7] Nagy G, Salhi S. Location-routing: Issues, models and methods. European Journal of Operational Research. 2007;177(2):649-72. [8] Lopes RB, Ferreira C, Santos BS, Barreto S. A taxonomical analysis, current methods and objectives on location-routing problems. International Transactions in Operational Research. 2013;20(6):795-822. [9] Prodhon C, Prins C. A survey of recent research on location-routing problems. European Journal of Operational Research. 2014. [10] Berbeglia G, Cordeau J-F, Gribkovskaia I, Laporte G. Static pickup and delivery problems: a classification scheme and survey. TOP. 2007;15(1):1-31. [11] Parragh SN, Doerner KF, Hartl RF. A survey on pickup and delivery problems. Journal für Betriebswirtschaft. 2008;58(2):81-117. [12] Mosheiov G. The Travelling Salesman Problem with pick-up and delivery. European Journal of Operational Research. 1994;79(2):299-310. [13] Nagy G, Salhi S. The many-to-many location-routing problem. TOP. 1998;6(2):261-75. [14] Webb MHJ. Cost Functions in the Location of Depots for Multiple-Delivery Journeys. Operational Research Quarterly. 1968;19(3):311-20. [15] Christofides N, Eilon S. An Algorithm for the Vehicle-Dispatching Problem. Operational Research Quarterly. 1969;20(3):309-18. [16] Laporte G, Nobert Y. An exact algorithm for minimizing routing and operating costs in depot location. European Journal of Operational Research. 1981;6(2):224-6. [17] Laporte G, Nobert Y, Pelletier P. Hamiltonian location problems. European Journal of Operational Research. 1983;12(1):82-9. [18] Ambrosino D, Scutella MG. Distribution network design: New problems and related models. European Journal of Operational Research. 2005;165(3):610-24. [19] Baldacci R, Mingozzi A, Calvo RW. An Exact Method for the Capacitated LocationRouting Problem. Oper Res. 2011;59(5):1284-96. [20] Contardo C, Cordeau JF, Gendron B. An exact algorithm based on cut-and-column generation for the capacitated location-routing problem. INFORMS Journal on Computing. 2014;26(1):88-102. [21] Jacobsen SK, Madsen OBG. A comparative study of heuristics for a two-level routinglocation problem. European Journal of Operational Research. 1980;5:378-87. [22] Nambiar JM, Gelders LF, Van Wassenhove LN. A large scale location-allocation problem in the natural rubber industry. European Journal of Operational Research. 1981;6(2):183-9. 26

[23] Srivastava R, Benton WC. The location-routing problem: Considerations in physical distribution system design. Computers & Operations Research. 1990;17(5):427-35. [24] Srivastava R. Alternate solution procedures for the location-routing problem. Omega. 1993;21(4):497-506. [25] Perl J, Daskin MS. A unified warehouse location-routing methodology. Journal of Business Logistics. 1984;5(1):92-111. [26] Perl J, Daskin MS. A warehouse location-routing problem. Transportation Research Part B. 1985;19(5):381-96. [27] Hansen PH, Hegedahl B, Hjortkjær S, Obel B. A heuristic solution to the warehouse location-routing problem. European Journal of Operational Research. 1994;76(1):111-27. [28] Tuzun D, Burke LI. Two-phase tabu search approach to the location routing problem. European Journal of Operational Research. 1999;116(1):87-99. [29] Clarke C, Wright JQ. Scheduling of vehicle from a central depot to a number of delivery points. Oper Res. 1964;12(4):568-81. [30] Bouhafs L, Hajjam A, Koukam A. A combination of simulated annealing and Ant Colony System for the capacitated location-routing problem. 2006. p. 409-16. [31] Wu TH, Low C, Bai JW. Heuristic solutions to multi-depot location-routing problems. Computers & Operations Research. 2002;29(10):1393-415. [32] Özyurt Z, Aksen D. Solving the Multi-Depot Location-Routing Problem with Lagrangian Relaxation. In: Baker EK, Joseph A, Mehrotra A, Trick MA, (Eds.) Extending the Horizons: Advances in Computing, Optimization, and Decision Technologies: Springer US; 2007. p. 125-44. [33] Chan Y, Baker SF. The multiple depot, multiple traveling salesmen facility-location problem: Vehicle range, service frequency, and heuristic implementations. Mathematical and Computer Modelling. 2005;41(8-9):1035-53. [34] Albareda-Sambola M. A compact model and tight bounds for a combined locationrouting problem. Computers & Operations Research. 2005;32(3):407-28. [35] Prins C, Prodhon C, Calvo RW. Solving the capacitated location-routing problem by a GRASP complemented by a learning process and a path relinking. 4or. 2006;4(3):221-38. [36] Prins C, Prodhon C, Calvo R. A Memetic Algorithm with Population Management (MA|PM) for the Capacitated Location-Routing Problem. In: Gottlieb J, Raidl G, (Eds.) Evolutionary Computation in Combinatorial Optimization: Springer Berlin / Heidelberg; 2006. p. 183-94. [37] Prins C, Prodhon C, Ruiz A, Soriano P, Calvo RW. Solving the capacitated locationrouting problem by a cooperative Lagrangean relaxation-granular tabu search heuristic. Transportation Science. 2007;41(4):470-83. [38] Barreto S, Ferreira C, Paixao J, Santos B. Using clustering analysis in a capacitated location-routing problem. European Journal of Operational Research. 2007;179(3):968-77. [39] Marinakis Y, Marinaki M. A Particle Swarm Optimization Algorithm with Path Relinking for the Location Routing Problem. Journal of Mathematical Modelling and Algorithms. 2008;7(1):59-78. [40] Marinakis Y, Marinaki M, Matsatsinis N. Honey Bees Mating Optimization for the Location Routing Problem. 2008. [41] Duhamel C, Lacomme P, Prins C, Prodhon C. A GRASP×ELS approach for the capacitated location-routing problem. Computers & Operations Research. 2010;37(11):191223. [42] Hemmelmayr VC, Cordeau JF, Crainic TG. An adaptive large neighborhood search heuristic for Two-Echelon Vehicle Routing Problems arising in city logistics. Computers and Operations Research. 2012;39(12):3215-28.

27

[43] Ting CJ, Chen CH. A multiple ant colony optimization algorithm for the capacitated location routing problem. Int J Prod Econ. 2013;141(1):34-44. [44] Escobar JW, Linfati R, Toth P. A two-phase hybrid heuristic algorithm for the capacitated location-routing problem. Computers & Operations Research. 2013;40(1):70-9. [45] Contardo C, Cordeau JF, Gendron B. A GRASP + ILP-based metaheuristic for the capacitated location-routing problem. Journal of Heuristics. 2014;20(1):1-38. [46] Contardo C, Hemmelmayr V, Crainic TG. Lower and upper bounds for the two-echelon capacitated location-routing problem. Computers & Operations Research. 2012;39(12):318599. [47] Nguyen VP, Prins C, Prodhon C. Solving the two-echelon location routing problem by a GRASP reinforced by a learning process and path relinking. European Journal of Operational Research. 2012;216(1):113-26. [48] Toyoglu H, Karasan OE, Kara BY. Distribution network design on the battlefield. Naval Research Logistics. 2011;58(3):188-209. [49] Toyoglu H, Karasan OE, Kara BY. A New Formulation Approach for Location-Routing Problems. Networks and Spatial Economics. 2012;12(4):635-59. [50] Karaoglan I, Altiparmak F, Kara I, Dengiz B. The location-routing problem with simultaneous pickup and delivery: Formulations and a heuristic approach. Omega-Int J Manage S. 2012;40(4):465-77. [51] Karaoglan I, Altiparmak F, Kara I, Dengiz B. A branch and cut algorithm for the location-routing problem with simultaneous pickup and delivery. European Journal of Operational Research. 2011;211(2):318-32. [52] Gen M, Cheng R. Genetic algorithms & engineering design. NewYork: Wiley; 1997. [53] Prins C. A simple and effective evolutionary algorithm for the vehicle routing problem. Comput Oper Res. 2004;31(12):1985-2002. [54] Kirkpatrick S, Gelatt CD, Vecchi MP. Optimization by Simulated Annealing. Science. 1983;220(4598):671-80. [55] Prodhon C. Le probl`eme de localisation-routage [PhD thesis]. France: Universit´e de Technologie deTroyes; 2006. [56] Barreto S. Analise e modelização de problemas de localização-distribuição [Analysis and modelling of location-routing problems] [PhD thesis]. Portugal: University of Aveiro; 2004. [57] Salhi S, Nagy G. A cluster insertion heuristic for single and multiple depot vehicle routing problems with backhauling. Journal of the Operational Research Society. 1999;50(10):1034-42. [58] Angelelli E, Mansini R. The vehicle routing problem with time windows and simultaneous pick-up and delivery. In: Klose A, Speranza MG, Van Wassenhove LN, (Eds.) Quantitative Approaches to Distribution Logistics and Supply Chain Management. Berlin: Springer-Verlag; 2002. p. 249-67.

28

Figures Procedure: Input:

Decoding (decoding of a given choromosome) d(t): depot status part of chromosome c(t): customer sequence part of chromosome problem data CLRPMB solution

Output: //Initialization Step 1. Ωk ← ∅ , N k ← ∅ , Ak ← ∅ //Determine the opened depots and assignment of the customers to the opened depots Step 2. Find “1” in d(t). If it fails find first closed depot t* (i.e. d(t*)=0) and d(t*)←1 Step 3. Apply following steps for assigning customers to the opened depots Step 3.1. i ← arg min {d (t ) | d (t ) ≠ 0} Step 3.3.

j ← arg min {d (t ) | d (t ) ≠ 0 and d (t ) > i} if j = ∅ then j = Nc

Step 3.4.

k* ← k | d (k ) = i

Step 3.2.

Ωk ← Ωk ∪ c (m) ,i ≤ m < j Step 3.6. if j = Nc then goto Step 4 Step 3.7. i ← j and goto Step 3.2 //Check the customer assignments Step 4. for each k ∈ N0 apply following steps Step 3.5.

(∑

)

di ; ∑ i∈Ω pi > CDk then

Step 4.1.

if max

Step 4.2.

i ← last element of Ωk

Step 4.3.

k * ← first depot satisfying max di + ∑ j∈Ω d j ; pi + ∑ j∈Ω p j > CDk*

Step 4.4.

if k = ∅ then k * ← k | d (k ) = 0

Step 4.5.

Remove i from Ωk , add to the end of Ωk* and goto Step 4

i∈Ωk

k

(

k*

k*

)

*

//Obtain optimal routes for each opened depot Step 5. for each k ∈ N0 apply following steps to obtain a feasible solution (S) Step 5.1.

if Ωk ≠ ∅ then construct an auxiliary graph G ( Ν k ,Αk ) using following steps

Step 5.1.1. Add node k and Ωk to N k

Step 5.1.2. for each (i, j ) ∈ N k , add (i, j ) arc to Ak , if it satisfy load constraints for vehicles. Step 5.1.3. Find least cost path from the first node to the last node of G ( Ν k ,Αk ) which leads to a optimal partition of corresponding customer sequence. Step 6. Calculate objective function value f ( S ) by using equation (1) and return the initial feasible solution, S Figure 1. Decoding procedure of a chromosome

29

(a)

(b)

(c)

(d)

Figure 2. A sample chromosome and corresponding solution

30

Procedure: InitFeasSol (generating an initial feasible solution) Input:CLRPMB problem data ′ : initial feasible CLRPMB solution Output: Sinit //Initialization Step 1. Ωc ← Nc , DelCDk ← CDk , PicCDk ← CDk //Determine penalty costs Step 2. Perform following steps for each customer i, i ∈ Ωc Step 2.1.

Find nearest depot k * (i) and second nearest depot j* (i) to customer i

Step 2.1.1.

k * (i ) ← arg mini∈NC ,k∈N0 {cik | d i ≤ DelCDk and pi ≤ PicCDk }

c d ≤ DelCD j and pi ≤ PicCD j } ( ) { ij i Step 2.2. Calculate the penalty costs for customer i Step 2.2.1. pci ← ci , j* (i) − ci ,k* (i)

Step 2.1.2.

j* (i ) ← arg mini∈N

C

, j∈ N 0 \ k*

//Assign customers to the nearest available depot Step 3. Perform following steps for each customer i, i ∈ Ωc Step 3.1.

Select a customer i* with largest penalty cost , i* ← arg maxi∈Ωc { pci}

Step 3.2.

if di* ≤ DelCDk* (i) and pi* ≤ PicCDk* (i) then

Step 3.2.1. assign customer i* to depot k * (i)

zi* ,k* (i) ← 1 and yk* (i) ← 1 , Step 3.2.2. update remaining capacity of depot k and remove customer i from Ωc

DelCDk* (i) ← DelCDk* (i) − di* , PicCDk* (i) ← PicCDk* (i) − pi* , Ωc ← Ω c \ i* , Step 4. if Ωc ≠ ∅ then goto Step 2 else goto Step 5 //Obtain vehicle routes for each opened depot Step 5. Perform following steps to obtain a feasible solution (S) Step 5.1. Construct a simple route for each customer i assigned depot k , xki ← 1 and xik ← 1 .

Calculate objective function value f ( S ) by using equation (1), Apply merge operator to the feasible solution (S) to obtain an initial feasible ′ solution Sinit ′ Step 6. Return the initial feasible solution, Sinit Step 5.2. Step 5.3.

Figure 3. The algorithm to generate an initial feasible solutions

31

Figure 4. A sample for crossover operator

(a)

(b) Figure 5. A sample for RouteMove operator

32

(a)

(b) Figure 6. A sample for CustMove operator

(a)

(b) Figure 7. A sample for Swap operator

33

(a)

(b) Figure 8. A sample for Merge operator

(a)

(b) Figure 9. A sample for Opt operator

34

Procedure: LocalSearch Input:An initial CLRPMB solution ( Sinit ) Output: An improved CLRPMB solution ( Sbest ) //Initialization Step 1. To ← Tinit , Sbest ← Sinit , Simp ← Sinit , Scur ← Sinit ,

f ( Sbest ) ← f ( Sinit ) , f ( Simp ) ← f ( Sinit ) , f ( Scur ) ← f ( Sinit ) //FirstPhase (Improving the Solution) Step 2. Perform following steps to obtain Snew from Scur Step 2.1. Crm ← ∅ , Ccm ← ∅ , Csw ← ∅ , Cme ← ∅ , Copt ← ∅ Step 2.2.

Generate NN rm neighbors of Scur by RouteMoveOperator and put them to Crm

Step 2.3.

Generate NN cm neighbors of Scur by CustMoveOperator and put them to Ccm

Step 2.4.

Generate NN sw neighbors of Scur by Swap Operator and put them to Csw

Step 2.5.

Generate NN me neighbors of Scur by Merge Operator and put them to Cme

Step 2.6.

Generate NN opt neighbors of Scur by Opt Operator and put them to Copt

Step 2.7.

Select a new solution Snew among neighbors

S new ← arg min { f ( S ′ )| S ′ ∈ (Crm ∪ Ccm ∪ Csw ∪ Cme ∪ Copt )}

Step 3. if f ( S new ) < f ( Simp ) then,

Simp ← S new , f ( Simp ) ← f ( S new ) ,

Scur ← Snew , f ( Scur ) ← f ( Snew ) and gotoStep 2 Step 4. if f ( Snew ) < f ( Scur ) then Scur ← Snew , f ( Scur ) ← f ( Snew ) else begin Δs = [( f ( S new ) − f ( S cur )) / f ( Scur )]* 100 if U (0,1) ≤ exp (−Δs / Titr ) then Scur ← Snew and f ( Scur ) ← f ( Snew ) end Step 5. Titr ← α Titr −1 Step 6. if Titr < T f then goto Step 7 else goto Step 2

//SecondPhase (Reconfiguration of Current Assignment of Routes) Step 7. Convert Simp to equivalent FLP Step 8. Solve FLP to obtain Sbest

Step 9. if f ( Sbest ) < f ( Simp ) then

Sinit ← Sbest , f ( Sinit ) ← f ( Sbest ) and goto Step 1 else Report Sbest and f ( Sbest ) Figure 10. The steps of Local Search procedure

35

Procedure: Memetic Algorithm for CLRPMB Input:CLRPMB data, algorithm parameters Output: Best solution ( Sbest ) Step 1. k ← 0 Step 2. Initialize P (k ) using InitFeasSol procedure Step 3. Evaluate P (k ) using Decodingprocedure Step 4. Improve P (k ) using LocalSearch procedure Step 5. Whilenot (termination condition) do Step 5.1. Select P1 and P2 by binary tournament from P (k ) Step 5.2. Apply crossover to P1 and P2 for obtaining C Step 5.3. Apply mutation to C Step 5.4. Improve C using LocalSearch procedure Step 5.5. Update P (k ) by deleting the worst solution and adding C Step 5.6. k ← k + 1 Step 6. Output best solution Figure 11. The overall pseudo-code procedure of the MA

36

Tables Table 1. Computational results for the CLRP instances of Prodhon’s test set Instances

UBCon

20-5-1a 54793 20-5-1b 39104 20-5-2a 48908 20-5-2b 37542 50-5-1a 90111 50-5-1b 63242 50-5-2a 88298 50-5-2b 67308 50-5-3a 86203 50-5-3b 61830 100-5-1a 274814 Avg ( gapmin ) = 0.27 ,

gap gap CPU Instances UBCon min avg max (sec.) min avg max 0.00 0.00 0.00 3.85 100-5-1b 213568 0.96 1.31 1.55 0.00 0.00 0.00 3.88 100-5-2a 193671 0.75 1.03 1.24 0.00 0.00 0.00 3.54 100-5-2b 157095 0.07 0.10 0.15 0.00 0.00 0.00 3.24 100-5-3a 200079 0.63 0.95 1.13 0.00 0.00 0.00 10.28 100-5-3b 152441 0.65 0.87 1.12 0.00 0.00 0.00 7.41 100-10-1a 289017 0.46 0.87 1.17 0.00 0.18 0.39 14.18 100-10-1b 234641 -0.22 0.06 0.51 0.05 0.31 0.68 9.29 100-10-2a 243590 0.18 0.58 1.04 0.00 0.03 0.14 12.20 100-10-2b 203988 0.00 0.20 0.52 0.00 0.00 0.00 9.05 100-10-3a 252421 0.92 1.18 1.50 1.33 1.60 1.84 59.95 100-10-3b 204597 0.10 0.16 0.27 Avg ( gapavg ) = 0.43 , Avg ( gapmax ) = 0.60 , Avg ( CPU ) = 86.96 seconds

CPU (sec.) 53.74 28.40 19.34 37.12 22.41 860.94 472.99 86.60 35.52 118.83 40.42

37

Table 2. Computational results for the instances derived from Barreto’s test set by SN demand separation procedure B&C Instances Srivastava86-8x2 Perl83-12x2 Gaskell67-21x5 Gaskell67-22x5 Min92-27x5 Gaskell67-29x5 Gaskell67-32x5(1) Gaskell67-32x5(2) Gaskell67-36x5 Ch69-50x5 Perl83-55x15 Ch69-75x10 Perl83-85x7 Daskin95-88x8 Ch69-100x10 Average

DSS X Y X Y X Y X Y X Y X Y X Y X Y X Y X Y X Y X Y X Y X Y X Y

gap

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 2.46 2.76 3.03 3.24 11.52 9.77 6.25 5.22 12.61 13.40 8.05 6.58 2.83

MA min 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 2.46 2.91 2.72 2.66 7.77 6.57 3.69 3.67 0.00 0.00 4.88 3.38 1.36

gap avg 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 2.94 2.93 2.72 2.66 7.93 6.62 3.78 3.68 0.42 0.69 4.88 3.38 1.42

max 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 3.20 2.95 2.72 2.66 8.15 6.79 4.12 3.74 2.60 1.64 4.88 3.38 1.56

CPU (sec.) 7.62 9.10 1.97 2.57 13.87 15.55 7.48 7.23 12.78 13.28 13.68 13.77 25.15 25.21 26.27 27.26 17.85 18.64 41.64 41.96 143.68 143.46 76.79 77.74 147.21 143.57 35.21 30.62 113.53 103.49 45.27

RI

min 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.72 3.01 0.04 3.47 0.00 3.22 0.00 2.99 0.35 8.51 0.20 7.02 0.41 5.28 0.07 5.70 2.53 13.75 1.61 19.95 0.00 6.45 0.00 5.02 0.20 2.81

MA without IPF gap CPU RI avg max (sec.) 0.00 0.00 7.62 0.00 0.00 0.00 9.10 0.00 0.00 0.00 2.85 0.00 0.00 0.00 2.89 0.00 0.00 0.00 11.56 0.00 0.00 0.00 11.48 0.00 0.00 0.00 7.28 0.00 0.00 0.00 7.52 0.00 0.00 0.00 13.23 0.00 0.00 0.00 13.12 0.00 0.00 0.00 16.45 0.00 0.00 0.00 17.23 0.00 0.00 0.00 23.07 0.00 0.00 0.00 23.16 0.00 0.13 0.64 23.67 0.64 0.13 0.64 24.75 0.64 0.00 0.00 19.54 0.00 0.00 0.00 19.44 0.00 4.00 5.28 45.88 2.15 4.47 5.17 46.39 1.61 3.45 3.68 53.89 0.44 3.40 3.72 54.60 0.70 9.98 10.98 83.24 2.23 8.49 9.29 85.27 2.07 5.87 6.94 99.72 1.56 5.99 6.72 103.31 0.96 21.45 28.37 31.65 11.39 23.00 26.16 29.26 4.92 7.29 8.28 135.85 1.69 5.33 6.21 136.72 1.12 3.43 4.07 38.66 1.07

38

Table 3. Computational results for the instances derived from Barreto’s test set by AM demand separation procedure B&C Instances Srivastava86-8x2

DSS

W Z Perl83-12x2 W Z Gaskell67-21x5 W Z Gaskell67-22x5 W Z Min92-27x5 W Z Gaskell67-29x5 W Z Gaskell67-32x5(1) W Z Gaskell67-32x5(2) W Z Gaskell67-36x5 W Z Ch69-50x5 W Z Perl83-55x15 W Z Ch69-75x10 W Z Perl83-85x7 W Z Daskin95-88x8 W Z Ch69-100x10 W Z Average

gap

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 4.73 3.81 26.46 26.27 21.29 14.69 30.81 30.56 28.14 18.45 14.79 10.20 7.67

MA min 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 4.73 3.63 26.10 24.79 16.67 11.38 27.61 27.88 20.08 16.55 8.42 5.98 6.46

gap avg 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 4.74 3.67 26.11 24.79 16.67 11.49 27.61 27.96 20.56 16.55 8.72 6.52 6.51

max 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 4.78 3.78 26.12 24.79 16.67 11.64 27.61 28.15 21.50 16.55 8.89 7.28 6.59

CPU (sec.) 6.49 10.16 2.37 2.88 13.59 14.71 11.36 11.00 17.07 17.92 15.53 13.36 22.86 24.87 22.39 24.67 18.24 18.80 44.27 42.37 429.14 429.08 86.07 84.37 299.91 307.91 33.80 29.73 136.98 139.72 77.72

RI

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.05 0.15 0.01 0.00 0.00 0.23 0.00 0.20 1.17 0.00 0.43 1.21 0.12

min 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.08 0.80 0.00 0.00 0.00 0.00 4.73 3.86 26.91 25.65 19.02 16.20 29.95 29.71 29.11 23.23 11.31 8.90 7.65

MA without IPF gap CPU avg max (sec.) 0.00 0.00 6.49 0.00 0.00 10.16 0.00 0.00 3.13 0.00 0.00 3.08 0.23 1.15 9.90 0.57 1.42 10.78 0.00 0.00 9.98 0.00 0.00 10.20 0.00 0.00 13.26 0.00 0.00 14.30 0.00 0.00 14.16 0.00 0.00 15.25 0.08 0.08 21.08 1.35 2.26 21.95 0.00 0.00 20.99 0.00 0.00 23.86 0.00 0.00 19.43 0.00 0.00 18.69 5.56 6.68 45.16 4.12 4.45 45.26 27.06 27.19 52.77 25.99 26.21 54.15 20.06 20.66 80.44 16.20 16.20 86.04 30.13 30.49 100.29 29.80 29.97 102.06 33.43 38.65 30.07 27.76 33.46 26.80 12.86 14.10 140.20 9.71 10.36 140.46 8.16 8.78 38.35

RI

0.00 0.00 0.00 0.00 1.13 1.40 0.00 0.00 0.00 0.00 0.00 0.00 0.00 1.43 0.00 0.00 0.00 0.00 1.83 0.57 0.22 0.44 1.36 0.00 0.41 0.20 6.88 7.67 2.44 1.32 0.91

39

Table 4. Computational results for the instances derived from Prodhon’s test set by SN demand separation procedure B&C Instances 20-5-1a 20-5-1b 20-5-2a 20-5-2b 50-5-1a 50-5-1b 50-5-2a 50-5-2b 50-5-3a 50-5-3b 100-5-1a 100-5-1b 100-5-2a 100-5-2b 100-5-3a 100-5-3b 100-10-1a 100-10-1b 100-10-2a 100-10-2b 100-10-3a 100-10-3b Average

DSS X Y X Y X Y X Y X Y X Y X Y X Y X Y X Y X Y X Y X Y X Y X Y X Y X Y X Y X Y X Y X Y X Y

gap

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.18 0.22 0.04 0.04 0.10 0.09 0.00 0.00 0.24 0.25 0.12 0.03 0.18 0.17 0.12 0.12 0.11 0.11 0.08 0.07 0.17 0.22 0.11 0.10 1.08 1.08 0.08 0.07 47.23 2.01 1.09 1.09 1.12 0.33 0.07 0.10 1.32

min 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.16 0.16 0.04 0.04 0.08 0.09 0.00 0.00 0.23 0.20 0.09 0.03 0.11 0.11 0.08 0.08 0.08 0.08 0.05 0.05 0.12 0.11 0.07 0.06 0.15 0.15 0.05 0.05 0.09 0.09 1.08 1.08 0.12 0.13 0.04 0.04 0.12

MA gap avg max 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.17 0.18 0.17 0.17 0.04 0.04 0.04 0.04 0.08 0.08 0.09 0.09 0.00 0.00 0.00 0.00 0.23 0.24 0.22 0.24 0.09 0.09 0.03 0.06 0.12 0.13 0.12 0.13 0.08 0.08 0.08 0.08 0.08 0.09 0.08 0.08 0.05 0.05 0.05 0.06 0.12 0.13 0.12 0.12 0.07 0.08 0.07 0.08 0.84 1.02 0.84 1.02 0.05 0.05 0.05 0.05 0.10 0.13 0.10 0.11 1.08 1.08 1.09 1.09 0.14 0.15 0.14 0.16 0.05 0.06 0.04 0.05 0.15 0.17

CPU (sec.) 19.25 17.38 12.78 12.56 14.85 15.46 13.01 11.46 55.09 51.49 51.94 47.70 49.35 46.69 36.27 34.36 44.75 40.64 43.03 43.76 153.49 142.50 124.82 117.79 244.52 250.51 180.52 211.16 72.51 88.07 53.63 57.37 208.71 210.60 189.04 197.01 173.37 169.04 136.75 115.82 123.72 105.08 101.00 92.58 95.03

RI

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.03 0.01 0.00 0.00 0.00 0.01 0.00 0.00 0.01 0.04 0.00 0.03 0.02 0.02 0.00 0.00 0.01 0.01 0.00 0.00 0.01 0.01 0.01 0.01 0.86 0.86 0.00 0.00 0.04 0.02 0.00 0.00 0.03 0.04 0.02 0.00 0.05

min 0.00 0.03 0.00 0.00 0.00 0.00 0.00 0.00 0.33 0.34 0.04 0.08 0.11 0.12 0.01 0.00 0.35 0.29 0.09 0.06 0.21 0.21 0.13 0.14 0.16 0.13 0.08 0.08 1.96 2.71 0.13 0.11 1.09 1.05 0.08 0.08 1.45 1.41 1.13 1.13 0.50 1.23 0.25 0.18 0.40

MA without IPF gap CPU avg max (sec.) 0.01 0.02 8.14 0.03 0.03 8.27 0.00 0.00 7.35 0.00 0.00 7.57 0.02 0.06 8.33 0.01 0.03 8.48 0.00 0.00 8.64 0.00 0.00 8.73 3.89 7.49 44.97 1.97 6.59 44.12 0.07 0.10 51.15 0.10 0.15 48.22 0.11 0.12 35.97 0.14 0.15 35.62 0.03 0.07 30.31 0.05 0.10 31.95 0.50 0.68 43.73 0.37 0.41 43.53 0.13 0.20 43.35 0.08 0.10 42.31 0.22 0.23 153.64 0.26 0.35 111.61 0.19 0.33 89.42 0.18 0.31 90.89 0.62 1.47 150.01 0.45 1.50 143.56 0.09 0.10 140.97 0.18 0.53 141.88 2.92 4.38 75.22 3.03 3.19 77.39 0.15 0.18 54.06 0.14 0.16 54.23 1.63 2.79 123.44 2.13 2.95 124.36 0.37 1.08 114.39 0.88 1.10 114.76 2.41 2.95 108.34 2.33 2.71 110.52 1.45 1.61 89.32 1.44 1.55 90.92 1.50 2.32 62.94 1.77 2.79 65.92 0.48 0.96 42.58 0.75 1.18 43.87 0.75 1.21 66.70

RI

0.02 0.00 0.00 0.00 0.06 0.03 0.00 0.00 6.66 5.86 0.07 0.07 0.02 0.03 0.06 0.10 0.33 0.12 0.11 0.05 0.02 0.14 0.20 0.17 1.29 1.35 0.03 0.44 2.32 0.47 0.05 0.05 1.65 1.84 0.99 1.01 1.46 1.26 0.47 0.41 1.78 1.53 0.71 0.99 0.78

40

Table 5. Computational results for the instances derived from Prodhon’s test set by AM demand separation procedure B&C Instances DSS 20-5-1a 20-5-1b 20-5-2a 20-5-2b 50-5-1a 50-5-1b 50-5-2a 50-5-2b 50-5-3a 50-5-3b 100-5-1a 100-5-1b 100-5-2a 100-5-2b 100-5-3a 100-5-3b 100-10-1a 100-10-1b 100-10-2a 100-10-2b 100-10-3a 100-10-3b Average

W Z W Z W Z W Z W Z W Z W Z W Z W Z W Z W Z W Z W Z W Z W Z W Z W Z W Z W Z W Z W Z W Z

gap

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 14.00 0.30 0.21 0.18 2.61 0.16 0.13 0.11 4.88 0.26 0.21 0.13 1.38 0.74 0.79 0.77 39.28 43.64 41.97 41.95 3.79 1.93 2.18 1.12 1.43 0.54 0.10 0.04 28.70 28.90 1.84 29.50 26.38 2.50 27.52 1.39 7.99

min 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 10.76 0.19 0.16 0.14 2.55 0.12 0.12 0.11 4.74 0.19 0.18 0.12 0.73 0.09 0.75 0.13 1.77 2.44 1.01 0.98 1.94 1.88 1.10 1.09 0.96 0.49 0.04 0.03 1.24 0.85 0.57 1.38 7.41 1.30 7.57 1.31 1.28

MA gap avg max 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 10.80 10.83 0.19 0.20 0.16 0.16 0.14 0.14 2.56 2.56 0.12 0.12 0.12 0.12 0.11 0.11 4.74 4.75 0.20 0.22 0.19 0.21 0.12 0.12 0.74 0.75 0.10 0.12 0.75 0.76 0.50 0.75 1.77 1.78 2.44 2.45 1.01 1.01 0.99 0.99 1.94 1.95 1.89 1.90 1.11 1.12 1.09 1.09 0.96 0.96 0.49 0.50 0.04 0.04 0.03 0.03 1.25 1.25 1.09 1.44 1.05 1.18 1.38 1.38 7.54 8.00 1.66 1.91 7.57 7.57 1.31 1.31 1.32 1.36

CPU (sec.) 13.17 15.17 14.85 14.44 12.83 14.08 14.20 14.32 81.60 57.18 74.62 44.39 86.88 47.91 42.70 40.93 61.30 73.53 58.82 52.82 139.97 150.44 143.08 127.69 101.19 126.83 139.81 151.74 90.08 105.71 81.45 95.55 674.83 727.38 603.49 406.96 262.82 177.54 207.67 156.47 275.00 171.60 164.89 96.38 141.23

RI

min 0.00 0.04 0.00 0.02 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.06 14.02 0.01 0.30 0.00 0.19 0.00 0.17 0.00 2.62 0.00 0.18 0.00 0.13 0.00 0.13 0.01 9.47 0.03 0.93 0.02 0.77 0.00 0.20 0.02 0.83 0.03 0.19 0.01 0.85 0.61 0.78 0.01 1.87 0.00 2.58 0.01 1.10 0.01 1.05 0.01 3.94 0.02 3.29 0.02 1.33 0.00 1.18 0.00 3.22 0.01 1.08 0.01 0.53 0.00 0.05 0.01 1.89 0.57 1.51 0.60 1.25 0.00 1.42 0.55 10.57 0.59 2.87 0.00 7.66 0.00 1.69 0.07 1.86

MAwithoutIPF gap CPU RI avg max (sec.) 0.11 0.31 6.38 0.26 0.07 0.12 7.34 0.10 0.00 0.01 7.83 0.01 0.00 0.00 6.78 0.00 0.00 0.00 7.51 0.00 0.00 0.00 7.60 0.00 0.03 0.05 9.07 0.05 0.00 0.00 7.79 0.00 14.75 17.46 35.92 2.93 0.34 0.41 41.35 0.11 0.23 0.32 45.37 0.13 0.19 0.21 42.34 0.04 2.65 2.70 31.33 0.08 0.77 2.60 35.77 2.36 0.21 0.26 37.11 0.12 0.14 0.16 35.11 0.02 13.65 14.93 39.63 4.74 3.92 5.62 44.42 4.44 1.30 1.95 45.26 1.15 0.23 0.27 45.49 0.06 0.84 0.85 140.11 0.02 0.20 0.22 150.59 0.03 0.85 0.86 143.22 0.01 0.78 0.78 103.97 0.00 1.88 1.88 101.30 0.01 3.38 3.93 136.16 1.29 1.64 2.18 140.29 1.06 1.14 1.25 138.64 0.20 5.38 5.99 86.72 1.94 3.47 3.66 95.80 0.36 2.28 2.78 78.42 1.41 1.48 1.60 73.85 0.42 3.53 3.90 104.87 0.65 1.54 1.96 120.45 0.86 0.84 1.16 126.00 0.62 0.13 0.39 120.08 0.35 3.01 3.83 104.50 1.86 2.33 2.68 112.54 1.14 1.86 2.31 104.65 1.05 1.60 1.97 100.65 0.54 22.58 26.33 88.42 12.48 3.07 3.18 89.37 0.31 11.87 26.31 66.50 14.76 1.70 1.71 60.61 0.02 2.64 3.39 71.07 1.32

41

Table 6. Computational results for MA without mutation operator MA Test Set Barreto Prodhon Average

DSS SN AM SN AM

min 1.36 6.46 0.12 1.28 2.30

gap avg 1.42 6.51 0.15 1.32 2.35

max 1.56 6.59 0.17 1.36 2.42

CPU (sec.) 45.27 77.72 95.03 141.23 89.81

RI

0.20 0.12 0.05 0.07 0.11

min 1.80 6.47 0.16 1.34 2.44

MA withoutmutation gap CPU (sec.) avg max 1.98 2.28 9.16 6.68 6.94 60.73 0.20 0.26 24.42 1.40 1.53 82.32 2.57 2.75 44.16

RI

0.45 0.42 0.10 0.18 0.29

Highlights

- We consider the capacitated location-routing problem with mixed backhauls. - We propose a memetic algorithm based on genetic algorithm, simulated annealing and integer programming formulation. - Memetic algorithm obtains the optimal or good quality solutions in a short computation time.

42