A memory-based iterated local search algorithm for the multi-depot open vehicle routing problem

Journal Pre-proof

A memory-based iterated local search algorithm for the multi-depot open vehicle routing problem ˜ Jose´ Brandao PII: DOI: Reference:

S0377-2217(20)30027-8 https://doi.org/10.1016/j.ejor.2020.01.008 EOR 16260

To appear in:

European Journal of Operational Research

Received date: Accepted date:

24 April 2019 7 January 2020

˜ , A memory-based iterated local search algorithm for the Please cite this article as: Jose´ Brandao multi-depot open vehicle routing problem, European Journal of Operational Research (2020), doi: https://doi.org/10.1016/j.ejor.2020.01.008

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. 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. © 2020 Published by Elsevier B.V.

1 Highlights



The iterated local search algorithm (MBILSA) which is presented to solve the multidepot open vehicle routing problem (MDOVRP) is based on the history of the search.

 

The MBILSA is effective and efficient, outperforming all existing algorithms for the MDOVRP. The perturbation procedures are based on the memory of the search and are simple



and effective. We propose a fast and effective heuristic method to improve the open routes, which is based on their structure.

2

A memory-based iterated local search algorithm for the multi-depot open vehicle routing problem José Brandão Departamento de Gestão, Escola de Economia e Gestão, Universidade do Minho, Largo do Paço, 4704 – 553 Braga, Portugal CEMAPRE, ISEG, University of Lisbon, Portugal. E-mail: [email protected] Abstract The problem studied in this paper is the multi-depot open vehicle routing problem, which has the following two differences in relation to the classical vehicle routing problem: there are several depots; the vehicles do not return to the depot after delivering the goods to the customers, i.e., the end of the route is not the starting point. There are many practical applications for this problem, however the great majority of the studies have only addressed the open vehicle routing problem with a single depot. In this paper, we present an iterated local search algorithm, in which the moves performed during the local search are recalled and this historical search information is then used to define the moves executed inside the perturbation procedures. Therefore, it is recorded the number of times that each customer is moved during the local search. Since this information is continuously updated and changes in each iteration, the search is driven to potentially better regions of the solution space, and increases the chance of avoiding cycling, even when using deterministic perturbations. The performance of this algorithm was tested using a large set of benchmark problems and was compared with other algorithms, and the results show that it is very competitive.

Key words: transportation; iterated local search; memory; multi-depot; vehicle routing.

1. Introduction Given a set of customers with known demands and geographical locations and a set of depots that store the goods for servicing customers’ demands, the multi-depot open vehicle routing problem (MDOVRP) consists of determining a set of routes where each route is assigned to exactly one vehicle, in order that all customers’ requirements are satisfied and a given objective is reached. Each route consists of a sequence of customers, which starts at one depot and finishes at one of the customers. A more formal definition of this problem is given in the paragraph below. The MDOVRP can be represented by a directed graph G = (V, A), where V is the set of vertices and A is the set of arcs that connect the vertices of V. V = D  N, where D = {0, ..., d-1} is the set of d depots and N = {d, d + 1, ..., d + n-1} is the set of n customers. The travel distance between i and j, (i, j)  A is dij. We assume that the matrix (dij) is symmetric

3 and satisfies the triangle inequality, i.e., dij = dji, and dij ≤ dik + djk, for all i, j, k  V. However, as the vehicles are not allowed to return to the depot, we set dij = + for i  N, j  D. Additionally, dij = + for i, j  D, i j, because a vehicle is not allowed to travel from one depot to another. Q is the capacity of each vehicle and qi is demand of vertex i, i  V; qi = 0, i  D and 0 < qi ≤ Q, i  N. The minimum number of vehicles required, m, is given by (1) and the maximum is n, assuming that each vehicle travels exactly one route. m=

 qi / Q .

(1)

iN

Therefore, the constraints considered in this problem are the following: all the vehicles have the same capacity; the total demand of all the customers on a route must not exceed the capacity of the vehicle; each customer is visited just once by one of the vehicles that fulfils his demand; the number of vehicles available at each depot is unlimited. As will be explained in more detail below, the objective is to minimise the total distance (or time) travelled by the vehicles. More generally, if s is a given solution, then the value of the objective function is called f(s). The MDOVRP is an extension of the classical vehicle routing problem (VRP) which considers more than one depot and routes that are open because the vehicles do not return to the starting depot after visiting the last customer of the route or, if this occurs, the vehicles make exactly the same trip, but in the opposite order for the collection of goods from the customers. The MDOVRP or its variants, like the OVRP or the OVRP with time windows are found quite frequently in the real world. This occurs because many companies do not have their own fleet of vehicles for the distribution of goods or, in some cases, the existing fleet is not enough and, therefore, some external vehicles also have to be hired. This problem also appears in other kinds of transportation problems, such as, for example, school bus routing, and in the planning of train or airplane routes, as long as the routes start and end at different points. It is commonly acknowledged that Schrage (1981) was the first author to bring attention to the VRP with open routes, as well as to its applications. A description of many practical applications of this type of problem can be found in the literature, which contain a high variety of different constraints, although most of them have only one depot. A few examples of these real applications are the following: Bodin et al. (1983) solved a problem of express airmail distribution in the USA; Fu and Wright (1994) modelled the plan of the freight services through the Channel Tunnel as an OVRP; Li and Fu (2002) studied a problem of school bus routing in Hong Kong; Tarantilis and Kiranoudis (2002) studied a case of the

4 delivery of fresh meat from two depots to a set of customers; Russell et al. (2008) addressed the newspaper distribution problem of a large newspaper producer in the USA; Yu et al. (2011) studied the transportation of materials for coal mines; López-Sánchez et al. (2014) investigated the case of the transportation of the employees of a bus company; Bauer and Lysgaard (2015) solved the cable layout problem of an offshore wind farm; Şevkli and Güler (2017) studied the delivery of newspapers; Atefi et al. (2018) solved the distribution problem of a Canadian major producer of cookies and snacks, and; Shen et al. (2018) solved a MDOVRP with time windows, with the objective to minimise the total cost of travelling, including the cost of carbon emissions. From the theoretical point of view, Sariklis and Powell (2000) and Brandão (2004) were the first to develop an algorithm specifically designed for the OVRP. Sariklis and Powell (2000) proposed a classical heuristic that takes into account the properties of the symmetric OVRP, without maximum route length constraints, while Brandão (2004) devised a tabu search algorithm which also deals with these constraints. In both cases, two hierarchical objectives were used: primary - the minimisation of the total number of vehicles required, and; secondary - the minimisation of the total travel distance (or time). Since then, many other papers have been published for solving the OVRP, all of which consider the same objective function as, for example, the following: Li et al. (2007) used a record to record algorithm; Pisinger and Ropke (2007) presented an adaptive large neighbourhood search algorithm that was also used to solve other types of VRP; Letchford et al. (2007) proposed an exact algorithm; Fleszar et al. (2009) developed a variable neighbourhood search algorithm; Vidal et al. (2014) apply a hybrid genetic search algorithm, which also solves several VRP variants; Atefi et al. (2018) used an iterated local search algorithm, although the main purpose of their study was to solve a practical problem, as stated before. Amongst those authors who have stated as objective function the minimisation of the total distance travelled, we cite the following examples: Tarantilis et al. (2005) used a list based threshold accepting algorithm; Zachariadis and Kiranoudis (2010) proposed a tabu search algorithm that explores large neighbourhoods, although in this paper the test problems are also solved by taking into account those two hierarchical objectives; Soto et al. (2017) developed an algorithm that combines multiple variable neighbourhood search with tabu search. As far as we know, the only papers that have been published for solving the MDOVRP, not considering the real applications, are the following: Liu et al. (2014), who proposed an exact method and a hybrid genetic algorithm; Lalla-Ruiz et al. (2016), who present an exact method based on a new integer programming formulation, and; Soto et al.

5 (2017). All of these authors used as objective function the minimisation of the total travel distance and do not take into account travelling time constraints. Therefore, as our algorithm is designed for the MDOVRP and as we wish to compare its performance with these ones, the same objective function is used. However, since we aim to strengthen the conclusions regarding the performance of our algorithm, we will also solve the OVRP and, in this case, both types of objective functions will be used. It should be noted that the main reason for assuming the objective of minimising the number of vehicles is their fixed cost (capital expenditure, maintenance, taxes, etc.), as the cost of using one more vehicle exceeds, in general, the savings obtained by the consequent reduction of the distance travelled. However, in the MDOVRP, the vehicles are hired, and therefore this cost does not exist directly for the company that contracts the service. Consequently, minimising the total distance travelled may be more appropriate in practice. Furthermore, we should also note that all these papers assume that the number of vehicles is identical to the number of routes, i.e., that each vehicle performs exactly one route. However, in practice this may not be the case, especially when the vehicles are owned by the company that sells the goods and the vehicles return to the depot, i.e., when the problem is what is usually called a ‘multi-trip VRP’ (see, for example, Brandão and Mercer, 1997). In this case, the allocation of the fixed cost of a vehicle to a route is rather different. Considering the similarities between the MDOVRP and the MDVRP, as well as the VRP, we suggest to interested reader the excellent literature reviews of Montoya-Torres et al. (2015), Laporte (2009), and Braekers et al. (2016) for the MDVRP, for the VRP and for many VRP variants, respectively. Representing a few examples of papers for the MDVRP we cite the following: Baldacci and Mingozzi (2009), Contardo and Martinelli (2014), who proposed exact algorithms; and the more recent papers of Alinaghian and Shokouhi (2018), Lahyani et al. (2018), Li et al. (2018), and Zhou et al. (2018), who studied several variants of the MDVRP. We can say that solving the MDOVRP comprises the following three interrelated decisions: assignment of the customers to the depots; definition of the clusters of customers assigned to each depot, and; definition of the sequence of routing the customers inside each cluster. This last sub-problem corresponds to the problem of finding the optimal Hamiltonian path (HP), which is NP-hard and, consequently, the whole problem is also NP-hard. Therefore, approximate algorithms are still required for the vast majority of practical applications.

6 In this paper, we present a memory-based iterated local search algorithm (MBILSA) for solving the MDOVRP. Therefore, the most relevant and innovative contribution of this paper is the use of the history of the search for the definition of the perturbation procedures, in a manner that has shown to be very effective. Thus, the perturbation is driven by the previous search, in the way that is explained in Section 2.5. As far as we know, this is the first iterated local search algorithm that incorporates memory of the search. Furthermore, the algorithm comprises two features that also play an important role in the quality of the final solutions: elite solutions for intensifying the search around good local optima, and; a constraint on the local search moves that can be performed, based on the vicinity among the vertices, in order to reduce the size of the neighbourhoods explored, whilst directing the search towards promising regions. Additionally, we used a simple algorithm for improving the Hamiltonian path, which takes into account the structure of this problem and has proven to be rather effective. Given that the MBILSA uses the memory of the search as well as some other features that appear in the Tabu Search (TS) algorithms, some readers may wonder whether the MBILSA is not simply a TS algorithm. The use of memory is the most distinctive feature of the TS methodology and is the main reason for its wide and successful application to many different types of optimisation problems. In TS, the memory is implemented in the algorithms through a tabu list, which gives rise to its name. In a given iteration, a solution, or a move, cannot be accepted for the search in the next iteration if it is in the tabu list, unless it complies with what is called an ‘aspiration criterion’, which, for example, means that it is better than any solution found so far. In most cases, the tabu list does not contain the solutions already used in the search, but rather some of their attributes, and the tabu status of a solution is removed after a given number of iterations. The type of memory implemented through the tabu list is known as ‘short term memory’, and it is compulsory for any TS algorithm, otherwise the algorithm cannot be classified as a TS algorithm. Many TS algorithms also use another kind of memory, which is called ‘long term memory’, or ‘frequency-based memory’, which is the kind of memory used in the MBILSA, although it is implemented differently from the usual in the TS algorithms. In conclusion, the MBILSA has this characteristic in common with some TS algorithms and includes some other features – such as elite solutions, which are also present in some TS algorithms, although the metaheuristic itself is rather different. One of the characteristics that contributes most to the distinction between metaheuristics is the mechanism used for escaping from the local optima. In the TS, this is the tabu list, while in the iterated local search (ILS), it is the perturbation, as is explained below.

7 In this sense, the ILS is closer to the variable neighbourhood search than to the TS. In conclusion, our algorithm is an ILS algorithm, as it follows the ILS methodology which is explained in the next section, and it is not a TS algorithm, because it does not use tabu restrictions for preventing cycling. However, it uses memory of the previous search similar to TS, but this memory is used in a very different way in both cases: in the TS the memory is used for preventing the acceptance of the solutions which have already been used a few iterations before (short term memory) or for promoting the acceptance of solutions which are distant from the region currently under exploration (long term memory); whereas in the our ILS algorithm, the memory is used for orienting the perturbation mechanism, although it certainly could be used in many other different ways without transforming it into a TS algorithm. Other distinct methodologies also use memory, such as scatter search and neural networks. The structure of the remainder of this paper is the following: Section 2 describes the MBILSA, starting with the description of its main features and finishing with the outline of the whole algorithm; in Section 3, we present the computational experiments and compare the quality of our algorithm with the best algorithms published, and; final conclusions are presented in Section 4. 2. The memory-based iterated local search algorithm 2.1. General iterated local search framework The iterated local search is classified as metaheuristic, because it allows the exploration of different regions of the solution space, without being immediately trapped in a local optimum. As the name suggests, the search consists of moving iteratively from a given local optimum to a new one. If the algorithm is well designed and the cycling is avoided, a new local optimum is found at each iteration, through a set of local search procedures, which are also called operators. This means that given a solution s, a local optimum, s*, is obtained by applying these improving operators. In order to escape from s*, this solution is perturbed, i.e., it is transformed by the execution of a set of moves, generating sp. Then, the local search operators are applied to sp originating another local optimum sp*. Note that the solution perturbed at each iteration can be the last local optimum that has just been generated, or any other of the set of local optima, S*, found during the previous search. The rules used for making this choice are called acceptance criteria. Therefore, each iteration comprises three steps: local search; choice of the solution to be perturbed, and; the perturbation. The whole process starts with the generation of an initial solution, and ends when the stopping criterion

8 (or criteria) is reached, such as, for example, the maximum number of iterations. The performance of an algorithm depends on the characteristics of these three steps. A detailed explanation of the iterated local search is provided by Lourenço et al. (2003). The abovedescribed algorithm is presented below, where sb is the best solution found so far. Algorithm 1: General structure of the iterated local search Generate the initial solution s, set sb = s. while the stopping criterion is not met do: s* = Local_search(s); if f(s*) < f(sb), set sb = s*. s = Acceptance_criteria (s*, S*); sp = Perturbation(s), set s = sp; end while In the next four sections, we describe how the general steps of the above-mentioned algorithm were defined in the MBILSA, following the same order: Section 2.2 – methods used for obtaining the initial solution; Section 2.3 – set of operators for performing the local search; Section 2.4 – the acceptance criterion; Section 2.5 – the methods used for executing the perturbation. Sections 2.6 and 2.7 describe other MBILSA features, which are mainly related to the acceptance criteria and perturbation, respectively.

2.2. Initial solution In this paper we use two methods for obtaining the initial solution: a nearest neighbour heuristic (NNH) and an insertion heuristic (IH). The NNH is always used, except if the first objective is to minimise the number of routes and the number of routes of the solution found by the NNH is greater than the value m given by (1). If this happens, another solution is determined by the IH and this one is used to initiate the iterated local search, if the solution generated by this method contains fewer routes than the one yielded by the NNH. Both methods are constructive and sequential, which means that the solution is obtained by generating one route after another until all the customers are routed. The nearest neighbour heuristic starts with the assignment of each customer to the depot nearest to him. This means that disjoint sets are created containing the customers whose distance to a given depot is lowest, i.e., the number of sets may be equal to the number of depots or may be just one. Next, each route is started with the un-routed customer nearest to one of the depots. After this, the admissible customer belonging to the same set and nearest to one of the vertices already in the route, is added to this route under construction just after the vertex to which he is closest. This process is repeated until all customers are routed or no other customer is admissible. A customer is said to be admissible if his insertion in the route

9 does not violate any of the constraints. When no customer is admissible in the current route, a new route is started, or the process stops if all the customers are already routed. After finishing each route, an improving nearest procedure is applied to it in order to try to reduce its distance, as will be explained later. The insertion heuristic is applied as follows. First, the customers are placed in a list by decreasing order of their demand. Then the un-routed customer with the largest demand is chosen for starting each route, which is then assigned to the nearest depot. This route is then filled by inserting the un-routed and admissible customers one-by-one, following the order of that list. Each customer is placed in the route in the position where the distance travelled is increased least. After ending each route, the route is submitted to the same improving procedure used in the NNH. 2.3. Local search operators Given a solution s, the neighbourhood of s, N(s), is the set of all solutions that can be reached from s, through a given kind of modification or move. Some of these solutions may be better than s and, if this happens, they are said to be local optima. The procedure that defines the rules of the moves is called local search operator or neighbourhood structure. The operators that were used are the following: i) cross over; ii) insertion; iii) interchange of chains; iv) swap; v) intra swap; vi) shift; vii) 3-opt, and; viii) nenbhp. The first four operators execute moves between the routes of s, while the other four perform moves inside each route of s. A move can only be executed if the resulting solution is feasible. If a move improves the current solution, then it is immediately executed and the search goes on, following the same initial order of the trial moves. When designing the MBILSA, we naturally took advantage of the knowledge that we acquired from our previous iterated local search algorithm for another type of vehicle routing problem (Brandão, 2018). Therefore, the MBILSA has some characteristics in common with this algorithm, in particular regarding the local search operators. However, 3-opt was not applied in it as well as the nenbhp, which plays an important role in the MBILSA. Furthermore, the most innovative feature of the MBILSA – the perturbations based on the memory of the search – is used in this algorithm for the first time. The cross over between two routes of s, ri and rj, consists of deleting one edge (i, j)  ri and another edge (m, n)  rj and then replacing these two edges by the edges (i, n) and (m, j), as shown in Figure 1. All the possible combinations are tried, i.e., every pair of routes and every edge in each route. It should be noted that in this procedure no neighbourhood

10 restriction is imposed because, if it was required that i and n were δ-nearest neighbours (δneighbours, for short) as well as m and j, then this would be rather restrictive. On the contrary, if the neighbourhood was defined between any customer of the two segments of the routes that are joined, then that would be almost as loose so as not to impose any constraint, but would rather require computing time for verifying the neighbourhood compliance. Furthermore, this may be a way of correcting an inappropriate choice of δ.

ri

i

j

rj

m

n

Figure 1. Representation of the cross over. In the insertion move, a customer i  ri is removed from its current route and a trial insertion is made in any of the other routes of the solution, containing at least one of the δneighbours of i, between any two customers, or near the depot. If the objective is to minimise the number of routes and the minimum calculated by (1) has not been reached yet, then no neighbourhood restriction is imposed if ri contains only one customer, as this move would generate a better solution, as the number of routes is reduced by one unit. When the objective is to minimise the total distance travelled, i can also be placed in a new route, if ri contains at least two customers and if this move improves the current solution. This operator is the only one that allows the creation of new routes and only if the objective is to minimise the total distance. The interchange of chains between two routes consists of exchanging chains of edges between them. This operator has a high flexibility and can comprise an enormous amount of combinations. However, we decided to apply only the following three different cases: 1) an edge (i, j)  ri is removed from ri and placed in any other route between two customers or the depot (exchange (2, 0), for short); 2) an edge (i, j)  ri exchanges its position with a customer m  rj (exchange (2, 1)), and; 3) an edge (i, j)  ri exchanges its position with another edge (m, n)  rj (exchange (2, 2)). In any of these cases, the edges moved are placed in the reverse order. The exchange (2, 2) is represented in Figure 2, which shows that the edges (i1, i), (j, j1), (m1, m) and (n, n1) were deleted and replaced by the edges (i1, n), (m1, j), (i, n1) and (m, j1). This example also shows that after the move the edges (i, j) and (m, n) are travelled in the reverse order. The neighbourhood restriction requires that at least one customer of the edge moved is a δ-neighbour of the customers adjacent to its new location. For example, the edges (i, j) and (m, n) of Figure 2 can only be moved if m1 is a δ-neighbour of j (or vice versa), or i

11 is a δ-neighbour of n1, and i1 is a δ-neighbour of n, or m is a δ-neighbour of j1. In case 1, if ri contains only two customers, then this constraint is not imposed.

ri

i1

rj

m1

i

j

j1

ri'

i1

m

n

n1

r j'

m1

i

j

j1

m

n

n1

Figure 2. Representation of the interchange of chains. A swap move consists of exchanging two customers belonging to two different routes. A trial swap of a customer i  ri with a customer j  rj requires the following two conditions: i) at least one customer of rj is a -neighbour of i, and at least one customer of ri is a neighbour of j, and; ii) at least one of the routes ri or rj contains more than one customer. The trial move works as follows. First, the customers i and j are removed from their current route, then they are inserted in the routes rj and ri, respectively. The search of the position for inserting each customer in the corresponding route is performed sequentially, starting just after the depot. If a better solution is found, the move is immediately executed. The intra swap, shift, 3-opt and nenbhp operators perform route improvement inside each route of the solution. The intra swap is applied to every route of s and consists of exchanging the position of two customers i and j of the same route. This is applied to every pair (i, j) in a systematic way, which follows the sequence of the customers inside the route, starting with the customer i next to the depot. If there is an improvement, the move is then immediately performed and the process goes on with the new route, but without restarting from the beginning of the route. For example, if i and j are, respectively, in the positions 3 and 5 of the route, and this contains 10 customers, then the swap implies that i goes to position 5 and j to position 3, and that the next trial swap is between the customers in the positions 3 and 6. The shift consists of moving a customer i, forward or backward to another position in the route. The trial moves are performed in a systematic way just as in the intra swap, such that for each customer of the route all the trial moves in both directions are evaluated without any repetition. The 3-opt is the well-known 3-optimisation for the travelling salesman problem (TSP), but adapted to the HP. This consists of removing 3 edges from the current route and replacing them by three other edges. It is programmed in such a way that the 2-opt is also included.

12 The nenbhp is similar to the classical nearest neighbour algorithm for the TSP. The nenbhp is applied to a given route and changes the sequence of its customers, in order to try to reduce the distance travelled. Taking into account that a route is a set of customers, C, plus the depot, this algorithm works as follows: starting at the depot, the next customer of the sequence is the customer of C which is not yet in the sequence and is nearest to the last customer that entered in the sequence. Extensive computational tests reported in the literature (see, for example, Syslo et al., 1983) have shown that this algorithm is not one of the best performers amongst the simple constructive heuristics for the TSP, in spite of being intuitively appealing. This may be explained by the fact that as long as fewer vertices remain to be visited, the freedom of choice is reduced and no freedom exists regarding the choice of the edge linking the last vertex of the tour to the initial one. However, in the HP, this last edge does not exist, which may contribute to better performance. The nenbhp is much faster than the 3-opt, although the experiments proved that, on average, better MDOVRP solutions are found if the nenbhp is used instead of the 3-opt. It should be noted that this is not a direct comparison of performance between both methods for the HP, as these operators are integrated in the whole algorithm. Due to these experiments, the nenbhp is applied in every iteration of the MBILSA, while the 3-opt is only applied when a new best solution is discovered. 2.4. Acceptance criterion As was explained in Section 2.1, the acceptance criterion is the set of rules that defines which solution is chosen for applying the following perturbation. The proximate optimality principle – POP – (Glover and Laguna, 1997) states that, in general, the good solutions are close to each other. Therefore, it seems that a good criterion would be to always choose the best known solution. However, depending on the type and the strength of the perturbation procedure, this creates the risk of going back to previous local optima, and could even start cycling. This risk is much higher in MBILSA, because it is deterministic. On account of this, we decided to use what is called the random walk criterion, i.e., the perturbation is applied in every iteration to the last solution found by the local search procedure, s*. This criterion, by itself, is not a guarantee that no cycling will occur, but reduces very significantly that risk, if it is combined with a good perturbation. Complementary to this, in order to take advantage of the POP and the history of the search, we stored a given number of best solutions found during the previous search, and, after a

13 given number of iterations, these solutions are those used for applying the perturbation, as will be explained in more detail in Section 2.6. 2.5. Perturbation As explained above, the role of the perturbation is to drive the search towards a different region of the solution space and, in this way, escape from a given local optimum, s*. Without the perturbation it would be impossible to go on with the search. The goal of any algorithm is to find the optimal solution, and, by the POP, it should be expected that it is not far apart from a given s*, especially if this is already good. In other words, this means that it is more likely to find final good solutions performing small steps, i.e., going progressively to different regions of the solution space, than walking in large steps, making big changes at each iteration. Therefore, based on our experiments, we are convinced that strong perturbations are more unlikely to produce good results, although they need to have enough strength to allow the exploration of different regions. This means that the perturbed solution, sp, should not be easily reversed by the local search. However, this does not depend only on sp, as it also depends on the order of applying the local search operators, as will be explained in Section 2.8. In the MBILSA, the perturbation is performed by applying the following three different operators to a given solution s: i) shift; ii) swap, and; iii) filling. All these operators have one characteristic in common: the moves performed by them are based on the memory of the previous search. This memory feature is based on the frequency of move of each customer and is implemented in two different ways: i) Every time a customer is moved from one route to another in the local search, this move is counted, but in the cross over only the moves of j and n of Figure 1 are considered. The moves inside each route are not taken into account. For a better understanding of the explanation below, we call FREQ to the array where the number of these moves of each customer is registered. ii) Every time that a customer is moved inside one variant of the filling perturbation (called filling 2), the move is registered in another frequency array which is named FREQ_FILL. Therefore, we consider two types of frequency of move, although the first is much more important for the MBILSA than the second one, as this is only used with one variant of the operator filling, as will be explained below.

14 The use of memory has two important consequences: i) the search is directed to the most promising regions of the solution space, and; ii) since the contents of this memory changes from one iteration to another, hardly one perturbed solution is identical to those obtained in the previous iterations. Considering the role of the perturbation, it can easily be concluded that all the perturbation operators execute non-improving moves. However, albeit rarely, it can come to pass that a new best solution is found, especially in the case of filling. This operator may eliminate one route and the resulting solution may have one route less than the best known one, which is a new best if the objective is to minimise the total number of routes. We describe in detail below how each perturbation procedure works. As it happens in the local search, a move can only be executed if the resulting solution is feasible. The operator shift works in the same way as in the local search, but follows three different rules: i) the moves are non-improving; ii) the order of moving the customers is defined by the decreasing order of FREQ, while in the local search the order of the trial moves is defined by their sequence in the route; iii) there is a maximum number of moves allowed. For example, the perturbation of the route r = [0, 12, 37, 15, 45, 33, 10, 39], taken from a test problem, yielded rp = [0, 10, 33, 45, 12, 37, 15, 39]. The perturbation swap consists of performing a given number of non-improving moves, which involve the exchange of two vertices belonging to two different routes. The swap (i, j), if i is a depot and j a customer, is only allowed if j is in a route with a depot different from i; otherwise, this would be an insertion move and not a swap. The other rules followed in the execution of this move are the following: i) the vertices are ordered by the decreasing order of FREQ and the trial moves are performed in this order; ii) each vertex can be moved only once; iii) there is a threshold on the number of moves that can be performed; iv) when both vertices are customers each one takes the position of the other in the respective route, but if one vertex is a depot and the other is a customer, then the depot is placed at the first position of the route and the customer is placed just after the depot. The move involving one or two depots causes a rather strong perturbation. It should be noted that rule i) implies that the vertices with higher frequency have priority for being moved. However, the vertices with lower frequency may also be moved, although this is more unlikely. There is another kind of perturbation swap, which is called swap complete. This move is different from the swap just described above, as follows: i) the trial moves follow the lexicographic order of the customers and the depots cannot be moved; ii) all the feasible moves are performed, but each customer can be moved only once. This move, which creates a

15 strong perturbation, is only applied at the beginning of Phases 4 and 5, which are described below, when the objective is to minimise the number of routes. Given a solution s, the perturbation filling consists of transferring customers from the other routes to one route of s. The route that receives the customers is the route with the least number of customers (least route, for short), or alternately, the second least number. This operator comprises three variants: filling 1 – the order of choosing the customers to be moved is defined by FREQ or, alternately, by the sum of the distance of the edges of s that link that customer to the other vertices (in both cases, it is used the decreasing order); filling 2 – the order of choosing the customers to be moved is defined by the increasing order of FREQ_FILL; filling 3 – the order of choosing the customers to be moved is defined by the decreasing order of FREQ, although in this case the route filled is the one with least demand or, alternately, the second least demand. The perturbation filling is only applied in the restart, whose meaning is explained in Section 2.6. Assuming that k is the counter of the number of iterations without improving the best known solution, filling 1 is applied in the restart if k mod 3 is 0 or 1, otherwise, filling 2 is applied. Filling 3 is only applied in Phases 4 and 5, when there is no move after applying filling 1 or 2, if the objective is to minimise the number of routes. 2.6. Elite solutions and restart During the search, a set of solutions containing the best solutions found up until then is saved. These solutions are selected at each iteration after finishing the local search. Given two solutions, s1 and s2, if the objective is to minimise the number of routes, s1 is better than s2 if it contains fewer routes. If they have the same number of routes, s1 is better than s2 if its travelling distance is smaller than that of s2. When the objective is to minimise the travelling distance, only this criterion is taken into account. These best solutions are called elite solutions, and they are ordered in the set by decreasing quality. It should be noted that the composition of this set is dynamic, as it may change at each iteration. The purpose of the elite solutions is to intensify the search around good solutions. At each restart, one of the elite solutions is used for applying the perturbation filling. The frequency of the restart is determined by K, whose value is defined later in Section 3.2, i.e., the restart happens every time that k is a multiple of K. At each restart, a different elite solution is used, following their order in the set. After using all of them, the process begins with the first one of the set.

16 We should note that the restart comprises the following actions: i) the current solution is replaced by an elite solution; ii) the perturbation filling is applied, and; iii) the value of  is changed, according to the rule defined in Section 3.2. 2.7. Phases of the MBILSA The basic phases of the MBILSA are Phase 1 and Phase 2. The only difference between them is that the perturbation shift is applied in Phase 1, while the perturbation swap is applied in Phase 2. These perturbations are applied in each iteration, during the respective phase, except in the iteration of the restart, in which case the perturbation filling is applied. These two phases may then be repeated alternately as many times as wanted, to try to improve the current best solution. However, as in each new application each phase may be different from the other in a few details, as is explained below, they were given different names, according to their sequence of application: Phases 3, 4, 5, 6, and 7 (odd numbers correspond to Phase 1, while even numbers to Phase 2, regarding the type of perturbation applied). The initial solution of Phase 1 is generated by the method defined in Section 2.2, and the initial solution of any other phase is the best solution found in the previous phase. The stopping criterion of any phase is a given number of iterations, without improving the best known solution. The specific value that was used is defined in Section 3.2. At the beginning of each phase, the value of the parameters are reset to their original value, the counters are reset at zero, the elite solutions are deleted, and the array FREQ is emptied. As was explained above, the main purpose of the MBILSA is to solve the MDOVRP, taking as objective function the minimisation of the total distance travelled. This means that for evaluating the performance of the MBILSA, the most important thing is to compare it with the other published algorithms conceived for the same purpose. Nevertheless, in order to strength this evaluation, it was also used for solving the OVRP and the results were compared with other algorithms created for solving this problem. As, in this case, some authors used as objective function the minimisation of the number of routes (more rigorously, the two hierarchical objectives defined in the Introduction), while other authors used the minimisation of the total distance, the MBILSA took into account both objectives, separately, using only the first five phases. However, since the large OVRP created by Golden et al. (1998) are more difficult, they were solved by the MBILSA in the following two stages: Stage 1 – where the five phases defined above were first applied, considering the hierarchical objectives; Stage 2 – then, starting with the best solution found in Stage 1, Phases 6 and 7 were applied with the objective of minimising the total distance, however, if a better solution for the hierarchical

17 objectives is found it is also recorded. Consequently, the solutions reported for each of the objectives are those obtained at the end of Stage 2 and the computing time includes both stages, i.e., it is identical for both objectives. In Stage 1, the number of routes is equal or smaller than in Stage 2, and no new routes can be created, which implies that the moves are more constrained. The consequences of this are that in Stage 1, the average computing time per iteration is slightly higher and the diversification is more difficult than in Stage 2. In order to increase the diversification in Stage 1, more elite solutions were used, and, at the beginning of Phases 4 and 5, the perturbation swap complete was applied. It should be noted that, as was explained in Section 2.3, when the objective is to minimise the distance, the insertion operator can create a new route if that move improves the current solution. Therefore, the creation of new routes in Stage 2 is allowed. 2.8. Programming decisions In this section, we present a set of programming features that contribute significantly to the computational efficiency and to the quality of the solutions yielded by the MBILSA. These features include data structures, the storage of information from one iteration to another, and also inside each local search operator during each iteration, the order of applying the local search operators, and optimising the application of the swap operator. Most of the information required during the execution of the programme is stored in rectangular arrays of different sizes and data types, according to the needs, which guarantees that the time of storage and retrieval of any item of information is O(1). For example, the δneighbourhood among the vertices is represented by a square array of size |V| of Boolean data type. This means that if vertex j is a δ-neighbour of vertex i, then the value in the entry (i, j) is true, otherwise it is false. Therefore, verifying whether a trial move complies with the neighbourhood restriction is computationally very efficient. There is some information which is valid from one iteration to the next, or at each iteration within each local search procedure, which can be easily kept in the central computer memory during the execution of the programme. As a consequence, a substantial amount of computing time can be saved as many calculations and operations are avoided. This information consists of registering the global and local status of each route of the current solution, s in two separate Boolean arrays. The global status indicates if each route of s has been modified by the current procedure, or by the previous ones, whereas the local status only records whether each route of s has been modified by the current procedure. This is explained in the next paragraph, through an example.

18 Let us assume, for example, that in one iteration only the routes ri and rj have been transformed by the local search and by the perturbation. This means that in the next iteration, any potential move of each local search operator should only be tried if this involves ri, or rj, or both, otherwise it is a waste of time. However, if a trial move involving ri and rk is executed, from now on, in this iteration, any potential move involving rk should also be considered. A similar reasoning is valid for the transformations that happen within each local search operator, as this is applied as long as there is an improvement, i.e., the complete cycle of all possible trial moves may be repeated several times. For example, let us assume that in the insertion move, in the beginning of a given iteration, all the routes of s can be used. Additionally, let us assume that the first cycle only transformed ri, and rj. If this happens, then the potential moves that should be evaluated in the next cycle are only those that involve ri, or rj, or both. Another feature that influences computing time is the order of generating and evaluating the trial moves inside each operator. Each potential move must pass through several filters successively: the global state, the local state, the neighbourhood constraint and the capacity constraint. The last step is the evaluation of the cost of the move, as this is the one that requires more computing time. The order in which the local search operators are applied is presented in the outline of the algorithm in Section 2.9. The first aim of the sequence used is to avoid reversing the moves performed by the perturbation, in order to guarantee the exploration of different regions of the solution space. It should be recalled that the moves executed by the perturbation operators are of the type insertion (applied in the restart), shift, and swap. Therefore, if it is convenient that these types of moves are not applied in the local search to the solution that has been perturbed immediately before, in order to reduce the chance of reversing the moves executed in the perturbation. Consequently, for each iteration in the local search, the operators interchange of chains and cross over are applied before insertion, shift is the fourth and swap is the sixth. Another aim of the sequence used is to avoid consecutively applying operators that execute moves with more similarities. For example, interchange of chains has more similarities with insertion, than with cross over. This is why cross over is placed between both. For the same reason, insertion and swap are separated by two operators. On the other hand, when a new best solution is found, the local search operators are applied in a different order, and 3-opt is applied instead of nenbhp, just to try to drive the search towards regions not yet explored.

19 Since the swap is the MBILSA procedure that requires the most computing time, its execution has been analysed in detail, in order to try to make it more efficient and effective. In the first set of computational experiments all the local search procedures were excluded, except swap. The main conclusions drawn from these experiments are the following: i) each time that the swap is applied, the total number of cycles can be rather large (remember that the cycle is repeated as long as there is an improvement), even if the problems are of small size; ii) in general, most of the improvement is obtained in the first four cycles, although the computing time increases a lot if there is no limit on the number of cycles. Another set of experiments analysed what happens when all the procedures of the algorithm are active. In this case, on average, the final solutions are better if the number of cycles is bound by four, but the computing time is a bit lower if there is no limit on the number of cycles. On the other hand, if the swap operator is not applied, the computing time is reduced rather significantly, however, on average, the solutions are considerably worse. Consequently, taking all these experiments into account, the best alternative is to apply the swap operator with a limit of four cycles, as the MBILSA is rather fast and it is most important to find high quality solutions. 2.9. Outline of the MBILSA The operation of the MBILSA can now be summarised in the following two algorithms. When the objective is to minimise the number of routes, f(s) < f(sb) if the number of routes of s is lower than that of sb or, if both solutions have the same number of routes, but the distance of s is lower. When the objective is to minimise the distance travelled, f(s) < f(sb) if the distance of s is lower than that of sb. The stopping conditions in Algorithm 3 are based on f(s), taking into account the objective function, as has just been explained above, i.e., the execution of the Local Search is repeated while f(s’) < f(s), where s’ is the solution at the end of the cycle, and s is the solution at its beginning. It should be noted that both in Algorithm 2 and Algorithm 3, the current solution, s, may be transformed by each of the operators, but for simplifying the presentation, we call it s throughout all the steps. The parameters used by the MBILSA are the following:

and  – maximum number of moves allowed in each perturbation shift, swap, and filling, respectively; δ – number of nearest neighbours; E – maximum number of different elite solutions explored; K – frequency of the restart; T – limit of the number of iterations without improvement allowed in each phase; k – counter of the number of iterations without

20 improving the best known solution in each phase; P – Maximum number of phases allowed; phase – number of the phase that is currently being executed. Algorithm 2: The MBILSA 1: Data: initial solution s, parameters , δ, E, K, T, P. 2: Set sb = s, k = 1, phase = 1. 3: while phase ≤ P do 4: while k ≤ T do 5: Apply Local Search to s // Algorithm 3 6: if(k mod K = 0) go to step 12 // Apply restart 7: Save s as a new elite solution if it is better than any of the existing ones 8: if s was improved in the previous step, set k = 1 9: if(phase mod 2 = 1) apply perturbation Shift to s 10: else apply perturbation Swap to s 11: Set k = k + 1 and go to step 4 12: do{ Update δ; 13: Set s = best elite solution not yet explored, or restart with the first of the list; 14: if(k mod 3 ≤ 1) apply perturbation Filling 1 to s; 15: else apply perturbation Filling 2 to s; 16: if s was not transformed by the previous procedures and phase = 4 or 5, apply Filling 3} 17: end while // Cycle of k ≤ T 18: Set s = sb, k = 1, phase = phase + 1, reset the original value of δ, delete the elite solutions, empty the array FREQ. 19: end while // Cycle of phase ≤ P

Algorithm 3: Local Search 1: Input: s, sb, and stopping conditions. 2: do{ 3: Apply Interchange of chains to s 4: Apply Cross over to s 5: do{ 6: Apply Insertion to s 7: Apply Shift to s 8: Apply Intra swap to s 9: } while there is an improvement of f(s) 10: Apply Swap to s 11: Apply Nenbhp to s 12: if(f(s) < f(sb)) do{ 13: Set 1 = ,  = n/2 14: Apply Cross over to s 15: Apply Swap to s 16: Apply Interchange of chains to s 17: do{ 18: Apply Insertion to s 19: Apply Shift to s

21 20: 21: 22: 23: 24:

Apply Intra swap to s } while there is an improvement of f(s) Apply 3-opt to s Set  = 1} } while stopping conditions are not met 3. Results of computational experiments 3.1. Benchmark problems In order to evaluate the performance of the MBILSA the following three sets of

benchmark problems from the literature were used: Set 1) problems taken from Christofides et al. (1979) and from Fisher (1994), which are identified by their original number, and are prefixed, respectively, with letters C and F; Set 2) problems generated by Golden et al. (1999), which are identified by the prefix O (notation used by Li et al., 2007); Set 3) problems proposed by Cordeau et al. (1997). The first two sets were originally proposed for the VRP, and the third set for the MDVRP, although they have been used by many authors as test problems for the OVRP and MDOVRP, respectively. The first set contains examples of small and medium size, the second set contains large size examples, while the third set contains a mix of the three sizes. The examples of this last set define a maximum number of vehicles available at each depot, but this constraint is not considered in the MDOVRP, as this was the assumption taken by the other authors who solved this problem. All these test problems are Euclidean, and the distances are calculated in double precision.

3.2. Parameters setting The MBILSA contains a small number of parameters because in its design we tried to keep this number as low as possible. Nevertheless, the total number of combinations of possible values is almost unlimited. Furthermore, a major difficulty is that the parameters are interdependent, i.e., the influence of the value of one parameter in the results depends on the value of each of the other parameters. Therefore, in our experiments we did not try to find the “best” combination of values, but simply good values that could be used with all test problems. This means that the value of each parameter was defined after some experiments, and that all results presented in the tables were obtained with this same set of parameters. The only exception to this is the “our best” column, which contains the best solutions found during all computational experiments. It should be noted that different sets of parameter values, or even changing the value of only one parameter, may yield different solutions.

22 However, the MBILSA is not sensitive to small changes in most of the parameters, and, therefore, the average results are, in general, similar, when these small changes are introduced. The definition of the final value of the parameters included three steps: i) theoretical analyses; ii) preliminary experiments; iii) systematic experiments. The first step takes into account the role of the parameters in the algorithm, the range of values that could be taken by them, and their possible influence on the final results. The conclusions of this step are also based on the knowledge acquired with other algorithms of the same type. The second step is a trial and error process, in which several values, based on the conclusions of the first step, were tried. In this step only MDOVRP examples were solved, and, through it, a short list of values of each parameter was defined, which was used in the systematic experiments. In the following paragraphs these points are explained with more detail. The -neighbourhood constraint was introduced to reduce the size of the search space, whilst trying to avoid the exclusion of high quality solutions. Therefore,  can be any integer between zero and n-1. If it is zero, then no other solutions can be explored, apart from the initial solution, and; if it is n-1, it does not impose any reduction of the search space. If R is the number of routes of the optimal solution, then each customer will be in a route with approximately n/R customers. Therefore, this ratio was taken as a guide for choosing the value of . However, if a route contains many customers, then each of them has its own neighbours, which increases the potential number of customers that can be inserted in that route. Due to this, we decided to use the same  for all the test problems. That ratio is, on average, 7.4 for the multi-depot test problems, and 38.6 for the large problem with one depot. However, since the first purpose of the MBILSA is to solve the MDOVRP, we decided to experiment the values 5 and 10. This kind of analysis was performed in relation to the other parameters. After the preliminary experiences, the following values of the parameters, which seem to demonstrate a good performance, were defined for testing the influence of the changes in one parameter at a time, while maintaining constant all the other parameters (the non integer values are rounded to the nearest integer): δ = 10,  = 3,  = 10,  = n/20, E = 5, K = 100 and T = 2000. The systematic experiments described below were performed with the single depot problems, C5, C11 and F12, and with the multi-depot problems, p06, p07, p10, p11, pr03 and pr05, described in the previous section. The objective was the minimisation of the distance and only Phases 1 and 2 were used. The results are summarised in Table 1, which presents the total distance and the total computing time, in seconds. Besides, some other experiments

23 included all the benchmark problems, and a few assumed that the objective was the minimisation of the number of routes. Table 1 – Results for different parameters combinations.



E K    5 10 2 4 15 20 n/10 10 200 Value Distance 11694.7 11665.0 11665.3 11681.3 11720.2 11696.2 11665.9 11699.0 11712.7 238 235 189 268 269 245 266 269 217 CPU (s)

We found that better results could be obtained by varying the value of δ during the search. Therefore, the initial value was set at 10, and was increased by 5 units at each restart if it is less than 50, otherwise, it is set at 13, and after that, it is increased again by 5 at every restart. After reaching 48, it is reset to 10. Furthermore, when a new best solution is found at the end of the local search (see Algorithm 3, step 13), δ is set at n/2, and then all the interroute local search operators are applied. After this, δ is reset to its previous value. As a consequence of this dynamic variation, the quality of the solutions has improved, and the computation time increased slightly. It should be noted that this variation also contributes to avoid cycling. Table 1, shows that the initial value of 5 gives slightly worse results than 10 (in both cases with dynamic variation). Accordingly, this was the value that was used for testing all the other parameters. The value defines the number of moves allowed in the perturbation shift. One move can, for example, originate the reversal of the order of two adjacent customers in the current route, although it can also happen that, for example, a customer from the end of the route goes to the beginning (as explained before, this depends on FREQ). This means that just one move can produce a rather strong perturbation of the route. Furthermore, this perturbation is applied to every route of the current solution. This explains why the perturbation shift is rather effective and only small values of should be used. The values tried for were 2, 3 and 4. Table 1 shows that, on average, 3 gives a very slightly better solution cost than 2, and that it is better than 4. Note that 3 is in the set of parameters defined above for performing the systematic experiments. Therefore, the corresponding solution cost is found in the second cell of Table 1, under δ = 10. As expected, in the case of 2, the computing time is lower, because the perturbation is weaker and, consequently, on average, the local optimum (algorithm 2, step 5) is reached quicker in each iteration. Since the results for 2 and 3 are so close, further tests were performed with the other benchmark problems. These tests enabled us to conclude that in the case of 3, the average results were slightly better. With regards to  the values 10, 15 and 20 were tried. From Table 1 we can conclude that, on average,  equal to 10 yields better solutions than 15 and 20, and this one gives worse

24 solutions than 15. The computing time is also better for 10. This means that when this parameter is larger, the diversification increases, which implies that the computing time also increases, as would be expected. Nevertheless, this is not true when comparing the results obtained with 15 and 20. Therefore, in order to draw safer conclusions, all the problems of Sets 1 and 3 were solved by using these three values of . In this case, the trend observed was the same: equal to 10 produced significantly better results, but the difference of using 15 or 20 is much smaller, i.e., the total distance is only a little bit higher with15, and the total computing time is almost the same. The limit  is only applied in the perturbations filling 1 and filling 3, while in filling 2 there is no limit. The aim of using a larger limit dependent on the size of the problem in the restart is just to allow a stronger diversification. The values of  tried were n/20 and n/10. It should be noted that when the objective is to minimise the number of routes, and the free capacity in the m vehicles is small, the number of moves hardly reaches . From Table 1, we can conclude that, on average, these two values yield almost the same results, but the execution is slightly faster with n/20. By performing additional tests, it was confirmed that, on average, n/20 gives slightly better results. The higher the value of E, the greater the diversification, and the lower the intensification originated by the restart. We experimented two values: 5 and 10. We observed that, on average, the solutions are better with E equal to 5 when the objective is to minimise the total distance (see Table 1), and they are better with 10 when the objective is to minimise the number of routes. These results are explained by the fact that, in this last case, the moves are more constrained and, consequently, it is more difficult to diversify the search, and a higher E helps for that. Therefore, both values of E have been used, according to the objective function. Restarting very frequently reduces the benefits of a random walk. On the contrary, a restart allows an intensification of the search around a good solution. Considering this, two values were experimented for K – 100, 200. Good results were obtained with both, but since with 100 they were slightly better, this is the value used. T is the stopping criterion of each phase. Therefore, the total number of iterations executed and, consequently, the computing depends on its value. However, the total number of iterations is not known in advance, as reaching the threshold T depends on the evolution of the search. This is just the benefit of using such stopping criterion since it is worth carrying on with the search while there are improvements, however, if this does not happen during many iterations, it is unlikely that another improvement will occur again. In the definition of

25 T, the main purpose is to find a good compromise between solution quality and computing time. In general, if there is enough diversification, then the higher the number of iterations executed, the better will be the solutions found. The experiments have demonstrated that an improvement can occur even after a very great number of iterations, especially if n is large. Nevertheless, in general, these improvements are small and, consequently, they do not justify spending a long computing time. This means that, as long as more iterations are executed, the likelihood of an improvement decreases and, if this happens, then there is a trend of this being smaller, on average. After some experiments, we found that T equal to 2000 gives a good compromise between computation time and solution quality. Therefore, this has been the value used in Phases 1, 2 and 3. Phases 4 and 5 are applied when the solutions are already quite good. Therefore, T needs to be higher than before, as otherwise hardly any new improvements would occur. Therefore, we decided to use T equal to 5000 in these two phases. Phases 6 and 7 are performed with an objective function that is different from the previous five, and, consequently, this is a source of the diversification. Therefore, considering this fact and the need of not spending too much computing time, these two phases were executed with T equal to 3000.

3.3. Comparison of the results Our algorithm was programmed in the C language, and it was executed on a desktop computer with an Intel Core i7-3820 processor at 3.6 GHz, and 32 GB of RAM. The results yielded by the MBILSA are compared, in terms of solution cost and computing time, with the those given by the following algorithms: Liu et al. (2014) (LJG14), Soto et al. (2017) (SSRR17), Li et al. (2007) (LGW07), Zachariadis and Kiranoudis (2010) (ZK10), Vidal et al. (2014) (VCGP14), and Atefi et al. (2018) (ASCR18). The first two are the most relevant for this paper, as they solved the MDOVRP, although Soto et al. (2017) also solved the OVRP, considering as objective to minimise the distance travelled. The other four solved the OVRP only, taking as objective to minimise the distance (ZK10), or to minimise the number of routes (LGW07, ZK10, VCGP14 and ASCR18). These algorithms seem to be among the best ones published for solving these problems and LGW07 was chosen because of its performance in solving the second set. The algorithms were executed on the following computers: LJG14 – 3.2 GHz Dual Core, 4 GB of RAM (we assume that it is a Pentium IV, but the complete reference was not given); SSRR17 – Intel Core i7-2600 at 3.40 GHz with 16 GB RAM; LGW07 – Athlon at 1 GHz; ZK10 – T5500 at 1.66 GHz; VCGP14 – Opteron 250 at 2.4 GHz; ASCR18 – Dual core 2.6 GHz, 2 GB of RAM (we assume that it is a Pentium IV,

26 since no other information is provided). A very rough idea of the speeds of the different computers can be obtained by using Dongarra’s (2014) table, which presents the number (in millions) of floating-point operations per second (Mflop/s) processed by each computer. These estimated speeds are given in Table 2.

Table 2 - Relative speeds of the computers used by the algorithms. Algorithm MBILSA LJG14 SSRR17 LGW07 ZK10 VCGP14 ASCR18

Processor Intel i7-3820, 3.6 GHz Pentium IV, 3.2 GHz Dual Core Intel i7-2600, 3.40 GHz Athlon, 1 GHz T5500, 1.66 GHz Opteron 250, 2.4 GHz Pentium IV, 2.6 GHz Dual Core

Mflop/s 1960 1450 1785 450 795 1385 1400

Speed scaled 1 0.7398 0.9107 0.2296 0.4056 0.7066 0.7143

The MBILSA is deterministic, and therefore, the results presented in the tables were obtained with one execution of the algorithm with the set of parameters defined in Section 3.2. The SSRR17 and the LGW07 are also deterministic. The other algorithms contain several stochastic parameters and, therefore, each test problem was solved several times and the authors presented the best solutions found and the average solutions, as well as, the average computing time. The LJG14 was executed 20 times, the ZK10 and the VCGP14 were executed 10 times, and the ASCR18 was executed 5 times. The average results were used to compare with our algorithm, although the best results are also used for defining the best solutions published. In the case of ZK10, the computing time is not directly comparable with the time taken by the other algorithms, as they do not present the total time, but just the time until finding the best solution. In the tables, the headers of the columns or rows not yet defined mean the following: BKS – best known solution, taken from LJG14, SSRR17, ZK10 and VCGP14, although it should be noted that these authors collected their data from several sources. Some solutions for the MDOVRP and also the OVRP are proven to be optimal. These solutions are marked with an asterisk. The MBILSA solutions equal to or better than the best known are written in bold, and also in italic if they are better. CPU – total computing time, in seconds. The original time taken by the other algorithms is converted using Table 2, in order to simplify the comparison. CPU deviation – the ratio between the total CPU time taken by each algorithm and the total CPU time spent by the fastest. Dist – the travel distance of a solution or, for some algorithms, average of several solutions generated by different executions, as explained above.

27 Distance deviation (%) – the difference between the total distance, for each set of problems, yielded by each algorithm and the total best known distance and then divided by this and multiplied by 100. MBILSA-Tph – a solution found by the MBILSA, when only the first two phases are executed. Our best – the best solution found by the MBILSA in the course of all experiments. R – the number of routes of a solution, or, for some algorithms, the average of several solutions generated by different executions. This number is irrelevant when the objective is to minimise the total distance. However, for the sake of completeness, it is presented for the OVRP, since ZK10 also does so. When the objective is to minimise the number of routes (Tables 6 and 7), this number is not written if it is identical to m in the solution given by the algorithm. Table 3 – Results for the MDOVRP test problems.

BKS Name n d Dist p01 50 4 386.18* p02 50 4 375.93* p03 75 5 474.57* p04 100 2 662.22 p05 100 2 608.73 p06 100 3 611.99 p07 100 4 613.96 p08 249 2 2852.56 p09 249 3 2625.04 p10 249 4 2539.60 p11 249 5 2508.45 p12 80 2 953.26* p15 160 4 1883.53* p18 240 6 2818.36* pr01 48 4 647.03* pr02 96 4 978.82* pr03 144 4 1423.48 pr04 192 4 1517.80 pr05 240 4 1746.57 pr06 288 4 2009.97 pr07 72 6 821.25* pr08 144 6 1254.45 pr09 216 6 1615.42 pr10 288 6 2013.86 Total 33943.03 Distance deviation (%) CPU deviation # best solutions # new best solutions

Our best Dist 386.18 375.93 474.91 662.22 607.53 611.99 608.90 2788.96 2578.49 2491.44 2468.45 953.26 1885.81 2818.36 647.03 979.82 1424.90 1517.60 1700.19 1985.59 821.25 1254.45 1592.77 1974.56 33610.59 -0.98  20 11

LJG14 SSRR17 Dist CPU Dist CPU 388.82 10.2 386.18 1.3 376.90 10.7 376.78 1.1 477.69 30.2 474.57 5.5 669.49 63.6 670.13 5.5 616.35 67.3 609.81 13.5 621.13 67.6 616.90 23.9 623.91 64.9 615.98 19.6 2913.63 805.9 2852.56 50.4 2755.28 815.5 2625.04 54.5 2616.02 827.7 2539.60 44.3 2598.56 824.0 2508.45 77.4 955.80 36.5 953.25 1.1 1899.98 249.3 1885.80 8.4 2854.81 768.1 2818.34 6.1 647.03 0.5 647.03 9.7 985.69 62.8 979.82 8.3 1458.75 223.6 1426.95 25.8 1586.67 272.2 1517.80 12.8 1755.26 757.9 1758.51 87.4 2085.62 1307.6 2009.97 62.0 824.09 30.1 821.25 1.1 1277.79 186.2 1266.17 5.1 1666.02 575.9 1615.42 15.8 2076.08 1287.9 2013.86 159.8 34731.37 9355.4 33990.17 691.2 2.32   0.14 15.6 1.2   1 14      

MBILSA-Tph Dist CPU 386.18 3.0 376.44 3.7 476.86 7.5 662.22 18.0 609.98 8.3 611.99 10.2 610.98 7.7 2797.49 73.2 2597.28 49.1 2503.96 44.5 2487.80 48.0 953.26 2.9 1885.81 12.9 2818.36 26.5 647.03 1.5 979.82 9.8 1431.44 21.3 1524.18 22.6 1703.53 44.7 2003.83 74.1 821.25 4.7 1258.83 11.0 1592.77 48.6 1992.93 44.1 33734.22 597.9 -0.62    16  9 

MBILSA Dist CPU 7.6 386.18 8.7 376.44 16.2 475.79 662.22 30.6 609.04 22.2 611.99 22.0 609.60 30.6 2794.10 148.1 2582.65 151.6 2500.43 115.3 2478.31 133.3 953.26 10.4 1885.81 36.6 2818.36 76.0 5.0 647.03 979.82 19.1 1429.38 42.2 1524.18 60.5 1700.93 129.9 1992.26 179.3 821.25 11.8 1258.64 32.4 1592.77 85.7 1975.36 174.5 33665.80 1549.6 -0.82  2.6  16  9 

28

The results of Table 3 show that the MBILSA performs better, in terms of distance, than the SSRR17 and substantially better than the LJG14, either if two phases or five phases are executed. Furthermore, even if only two phases are executed, our algorithm improves the best known total distance substantially. Similar conclusions can be taken by observing the number of best solutions (16 out of 24) and the number of new best solutions (9). Additionally, the “our best” column shows that our algorithm was able to match, or improve, 20 of the best known solutions, finding 11 new best ones. Additionally, the number of solutions that are proven to be optimal is 9, and our algorithm found 6 of them. In terms of computing time, executing two phases of the MBILSA is a bit faster than the SSRR17, and is 2.6 times faster than executing five phases of the MBILSA, and is much faster than the LJG14. The Tables 4 and 5 present the results for the OVRP, when the objective is to minimise the total distance travelled, for Set 1 and Set 2, respectively. Table 4 – Results for the OVRP, Set 1 (distance minimisation).

BKS Name n R Dist C1 50 6 412.96 C2 75 11 564.06 C3 100 9 639.26 C4 150 12 733.13 C5 199 17 868.11 C11 120 10 678.54 C12 100 10 534.24 F11 71 4 177.00 F12 134 8 761.68 Total 87 5368.98 Distance deviation (%) CPU deviation

Our best R Dist 6 412.96 11 564.06 9 639.26 12 733.13 17 867.89 10 679.98 10 534.24 4 177.00 8 761.45 87 5369.97 0.02   

R 6.0 11.0 8.8 12.0 17.0 9.3 10.0 4.0 8.0 86.1  

ZK10 SSRR17 Dist CPU Dist CPU 412.96 12.2 412.96 0.2 564.06 25.1 564.06 0.7 641.20 39.7 639.26 6.9 734.92 72.6 733.13 34.0 871.67 108.7 868.11 887.7 679.27 41.8 685.08 317.4 534.24 13.0 534.24 41.8 177.00 41.8 177.00 22.4 763.10 92.9 769.55 43.2 5378.42 447.8 5383.39 1354.3 0.18 0.27   1.5 4.5  

R 6 11 9 12 17 11 10 4 8 88  

MBILSA Dist CPU 412.96 6.0 564.06 13.6 639.88 22.6 733.13 43.1 868.93 92.0 680.76 50.4 534.24 16.6 177.70 9.4 761.63 44.9 5373.29 298.6 0.08   

From Table 4 it can be concluded that all the algorithms have a good performance, as their total distance is very close to the best known. However, the MBILSA is a bit better than the SSRR17, and is also slightly better than the ZK10. Additionally, the MBILSA found two new best solutions. In terms of computing time, the MBILSA is substantially faster than the other two. Table 5 – Results for the OVRP, Set 2 (distance minimisation).

BKS Name n Dist O1 200 5988.35 O2 240 4549.46 O3 280 7705.52

Our best ZK10 SSRR17 R Dist R Dist CPU Dist CPU 6 5988.35 6 5996.71 236.5 6018.5 31.8 9 4552.55 9 4560.32 312.3 4574.12 530.9 7 7707.36 7 7733.65 341.9 7705.52 256.8

R 6 9 7

MBILSA Dist CPU 5999.75 118.9 4552.55 241.9 7707.36 387.4

29 O4 320 7251.30 10 O5 360 9138.82 9 O6 400 9793.72 9 O7 440 10338.42 10 O8 480 12390.15 10 Total 67155.74 70 Distance deviation (%) CPU deviation

7242.14 9152.47 9794.83 10362.67 12398.89 67199.26 0.06 

10 9 9 10 10 70  

7258.07 404.4 7269.38 131.6 10 7269.74 366.7 9176.85 449.0 9138.82 318.9 8 9152.68 696.9 9819.56 494.8 9801.86 288.1 9 9799.80 1029.8 10362.44 489.2 10338.42 851.0 10 10381.88 1279.9 12428.91 573.5 12390.15 760.3 10 12398.89 1211.5 67336.51 3301.6 67236.77 3169.4 69 67262.65 5333.0 0.27 0.12 0.16     1.0 1.7     

The results of Table 5 show that all the algorithms generate solutions very close to the best known. Furthermore, the solutions given by the SSRR17 are slightly better than those given by the MBILSA and these ones are a bit better than the solutions of the ZK10. In terms of computing time, the MBILSA is slower than the other two, but it should be noticed that the CPU reported by ZK10 has a different meaning, as not all computing time is included. In turn, Tables 6 and 7 present the solutions obtained by the MBILSA for the OVRP for the Sets 1 and 2, respectively, when the objective is to minimise the number of routes. Table 6 – Results for the OVRP, Set 1 (minimisation of the number of routes). BKS Our best LGW07 ZK10 VCGP14 ASCR18 MBILSA Name n m Dist Dist Dist CPU Dist CPU Dist CPU Dist CPU Dist CPU C1 50 5 416.06* 416.06 416.06 1.4 416.06 11.4 416.06 22.4 416.58 11.7 416.06 6.8 C2 75 10 567.14 567.60 567.14 7.2 568.38 29.2 568.15 27.9 567.14 9.7 567.77 14.1 C3 100 8 639.74* 639.74 639.74 9.1 639.98 39.3 639.74 46.4 639.88 138.3 639.74 20.5 C4 150 12 733.13 733.13 733.13 29.5 733.93 82.7 733.13 94.5 733.13 166.7 733.13 36.3 C5 199 16 878.94 886.75 924.96 87.4 895.62 134.7 890.15 225.7 880.46 127.1 893.61 108.4 C11 120 7 682.12 683.63 682.54 27.9 682.34 30.8 682.12 82.5 682.12 117.9 689.49 43.7 C12 100 10 534.24* 534.24 534.24 7.6 534.24 19.1 534.24 33.3 534.24 109.7 534.24 16.2 F11 71 4 177.00* 177.00 177.00 4.5 177.00 53.5 177.00 35.5 177.00 96.6 177.00 9.5 F12 134 7 769.55 769.55 769.66 36.3 770.57 112.8 769.68 93.4 769.55 152.1 770.51 51.9 Total 79 5397.92 5407.70 5444.47 210.9 5418.12 513.5 5410.27 661.6 5400.10 929.9 5421.55 307.4 0.78 0.04 Distance deviation (%) 0.10 0.29 0.14 0.35      CPU deviation 4.4 1.5 2.4 3.1       

The conclusions that can be drawn from Table 6, bearing in mind that the best solutions were obtained by several state-of-the-art algorithms (4 solutions are optimal), are as follows. The five algorithms found a solution with the minimum number of routes for every problem. In terms of distance, all the algorithms yield good solutions, as all of them are close to the best known and most of them match the best known. However, on average, the ASCR18 performs a little better than VCGP14 and this performs a little better than the ZK10, which is very slightly better than the MBILSA. The LGW07 is worse than the other four. In terms of computation time, the LGW07 is faster than the other three and it is a bit faster than the MBILSA. The ASCR18 is substantially slower than the others. Table 7 – Results for the OVRP, Set 2 (minimisation of the number of routes).

Name n m O1 200 5 O2 240 9

BKS Our best Dist Dist 6018.52 6018.52 4549.46 4552.55

LGW07 Dist CPU 6018.52 83.9 4584.55 100.9

ZK10 Dist CPU 6018.52 257.6 4562.88 337.5

MBILSA Dist CPU 6024.23 118.9 4552.55 241.9

30 O3 280 7 7705.52 O4 320 10 7251.30 O5 360 8 9193.15 O6 400 9 9793.72 O7 440 10 10338.42 O8 480 10 12390.15 Total 68 67240.24 Distance deviation (%) CPU deviation

7707.36 7242.14 9152.47 9794.83 10362.67 12398.89 67229.43 -0.02 

7732.85 113.1 7735.10 373.6 7291.89 131.7 7264.32 409.3 9197.61 176.0 9243.69 644.9 9803.80 224.4 9824.44 449.4 10374.97 214.8 10363.28 443.7 12429.56 258.7 12430.06 516.3 67433.75 1303.5 67442.29 3432.3 0.29 0.30   2.6   

7707.36 7269.74 9152.68 9799.80 10381.88 12398.89 67287.13 0.07 

387.4 366.7 696.9 1029.8 1279.9 1211.5 5333.0  4.1

The first conclusion that should be taken from Table 7 is that the three algorithms found a solution with the minimum number of routes for every problem. The total distance of the solutions of the three algorithms is rather close to the best known, which is a good indication of the quality of them. However, the MBILSA performs a little better than the other two, even improving the best known total distance a little bit and it discovered two new best solutions. The ZK10 seems to be a bit faster than the MBILSA, although we should be aware that the CPU reported by the authors does not include all the computing time spent by the ZK10. The LGW07 is faster than the other two. In the next paragraphs we carry out a statistical analysis of the main conclusions made about the performance of our algorithm. The primary purpose of our algorithm is to solve the MDOVRP. Accordingly, the most important statistical to test to be applied is to verify whether there is a significant difference between the quality of the solutions for the MDOVRP yield from using the MBILSA and those produced by the LJG14 and the SSRR17. In order to do this, we started by calculating the deviation to the best known solution (BKS), i.e., the percentage difference between the travel distance of BKS for each problem of Table 2 and the solution distance given by each of these algorithms, as follows: 100(z – z*) / z*, where z* is the distance of BKS and z is the corresponding solution distance produced by each algorithm (simply called deviation, for short). This deviation metric has the advantage of not being affected by the units in which the distance is measured. The mean deviation is 1.82, 0.20, and -0.50 for LJG14, SSRR17, and MBILSA, respectively. Our objective is to test whether the mean deviation given by the MBILSA is statistically significantly different from that given by both the LJG14 and SSSRR17. We therefore have to compare separately two differences of the means: i) the mean of the LJG14, less the mean of the MBILSA, and; ii) the mean of the SSRR17, less the mean of the MBILSA. To do this, we define as the null hypothesis (H0) that the difference is less or equal to zero, and as the alternative hypothesis (Ha) that the difference is greater than zero. Since the three algorithms solved the same problems, we obtained pairs

31 of solution values, and the purpose is precisely to perform a pairwise comparison of the MBILSA with the other two algorithms. Furthermore, having 24 problems solved is a reasonably large number, and thus we can apply the matched-samples t-test. The p-values obtained, assuming H0, were 3.5 x 10-6 and 5.6 x 10-4, for the comparisons i) and ii), respectively. This means that there is a strong evidence for rejecting H0 and for accepting Ha, i.e., that the difference of the means is statistically significant. We can therefore conclude with a high degree of confidence that the MBILSA performs better than the other two algorithms. The previous statistical test requires the assumption that the sample deviations come from a Normal population, but this assumption can be avoided by performing the Wilcoxon signed rank test. As in the previous test, we used matched-samples of the deviations, however, since this is a rank test, the distances could be used, instead of deviations, which gives approximately the same results. This test includes two cases: i) the median difference between the deviations of the LJG14 and the deviations of the MBILSA; ii) the median difference between the deviations of the SSRR17 and the deviations of the MBILSA. The hypotheses are also similar to those used before: whereby for H0, the median difference is less or equal to zero, and for Ha, the median difference is greater than zero. Given that the size of the samples after removing the zero differences is 23 in i), and 20 in ii), the Normal approximation is valid for calculating the p-values, which are 1.4 x 10-5 and 1.8 x 10-3, respectively. Accordingly, there is a strong evidence against H0, and, consequently, there is evidence for accepting Ha. This means that the results of both types of statistical tests are consistent. This, therefore, strengthens the level of confidence that the MBILSA has a better performance than either of the other two algorithms. To conclude, we present the summary of the same statistical tests, but considering the MBILSA-Tph, instead of the MBILSA. The mean deviation for the MBILSA-Tph is -0.36. The t-test gives the p-values 3.0 x 10-6 and 2.0 x 10-3 when comparing the MBILSA-Tph with the LJG14 (Case i) and with the SSRR17 (Case ii), respectively. The Wilcoxon signed rank test gives p-values of 1.4 x 10-5 and 4.0 x 10-3 for the cases i) and ii), respectively, using the Normal approximation. These results provide strong evidence that even the MBILSA-Tph performs substantially better than both the LJG14 and the SSRR17. 4. Conclusions In this research we developed an iterated local search algorithm that uses the memory of the search for the definition of the perturbation procedures. This is the basis of the creation

32 of a rather simple, effective, and efficient algorithm. In fact, the MBILSA was able to find very good solutions for the MDOVRP in a short computing time. In fact, it yielded 16 solutions out of 24 that match or improve the best known solutions, 9 of which are new best solutions. In addition, the best known total travelling distance was improved by 0.82%, and, considering the eight problems for which the travelling distance is worse than the best known, the maximum deviation is 0.42%. Furthermore, considering all the computational experiments, the MBILSA matched or improved 20 solutions of the same set of 24, generating 11 new best solutions and improving the best known total travelling distance by 0.98%. Additionally, the MBILSA was also tested in the solution of the OVRP and the results are also quite competitive with the best existing algorithms. The good behaviour of the MBILSA may also be explained by the following features: the restriction of the search to the most promising regions of the solution space through the application of a neighbourhood constraint; the use of a simple heuristic for improving each route which takes into account the structure of the Hamiltonian path, and; the use of elite solutions for intensifying the search around good local optima. References Alinaghian, M., & Shokouhi, N. (2018). Multi-depot multi-compartment vehicle routing problem, solved by a hybrid adaptive large neighborhood search. Omega, 76, 85–99. Atefi, R., Salari, M., L, Coelho, L., & Renaud, J. (2018). The open vehicle routing problem with decoupling points. European Journal of Operational Research, 265, 316-327. Baldacci, R., & Mingozzi, A. (2009). A unified exact method for solving different classes of vehicle routing problems. Mathematical Programming, 120, 347–380. Bauer, J., & Lysgaard, J. (2015). The offshore wind farm array cable layout problem: a planar open vehicle routing problem. Journal of the Operational Research Society, 66 (3), 360-368. Bodin, L., Golden, B., Assad, A., & Ball, M. (1983). Routing and scheduling of vehicles and crews: the state of the art. Computers and Operations Research, 10, (2), 63-211. Braekers, K., Ramaekers, K., & Nieuwenhuyse, I. (2016). The vehicle routing problem: State of the art classification and review. Computers and Industrial Engineering, 99, 300-313. Brandão, J., & Mercer, A. (1997). A tabu search algorithm for the multi-trip vehicle routing and scheduling problem. European Journal of Operational Research, 100 (1), 180-191. Brandão, J. (2004). A tabu search algorithm for the open vehicle routing problem. European Journal of Operational Research, 157, 552-564. Brandão, J. (2018). Iterated local search algorithm with ejection chains for the open vehicle routing problem with time windows. Computers and Industrial Engineering, 120, 146–159. Christofides N., Mingozzi, A., & Toth, P. (1979). The vehicle routing problem. In: Christofides, N., Mingozzi, A., Toth, P. and Sandi, C. (Eds), Combinatorial Optimization, Willey, Chichester, pp. 313-338. Contardo, C., & Martinelli, R. (2014). A new exact algorithm for the multi-depot vehicle routing problem under capacity and route length constraints. Discrete Optimization, 12, 129–146. Cordeau, J.F., Gendreau, M., & Laporte, G. (1997). A Tabu search heuristic for periodic and multi-depot vehicle routing problems. Networks, 30 (2), 105–119. Dongarra, J. (2014). Performance of various computers using standard linear equations software. Report CS-8985, University of Tennessee. Fleszar, K., Osman, I., Hindi, K. (2009). A variable neighbourhood search algorithm for the open vehicle routing problem. European Journal of Operational Research, 195(3), 803–809. Fisher, M. (1994). A polynomial algorithm for the degree-constrained minimum k-tree problem. Operations Research, 42 (4), 775-779.

33 Fu, Z., & Wright, M. (1994). Train plan model for British Rail freight services through the Channel Tunnel. Journal of the Operational Research Society, 45, 384–391. Golden, B., Wasil, E., Kelly, J., Chao, I-M. (1998). The impact of metaheuristics on solving the vehicle routing problem: algorithms, problem sets, and computational results. In: Crainic, T., Laporte, G. (Eds.). Fleet Management and Logistics. Kluwer, Boston, 33–56. Glover, F. and Laguna, M. (1997). Tabu Search. Kluwer Academic Publishers. Lahyani, R., Coelho, L., & Jacques R. (2018). Alternative formulations and improved bounds for the multidepot fleet size and mix vehicle routing problem. OR Spectrum, 40, 125–157. Lalla-Ruiz, E., Expósito-Izquierdo, C., Taheripour, S., & Voß, S. (2016). An improved formulation for the multi-depot open vehicle routing problem. OR Spectrum, 38, 175-187. Laporte, G. (2009). Fifty years of vehicle routing. Transportation Science, 43, 408–416. Letchford, A., Lysgaard, J., & Eglese, R. (2007). A branch-and-cut algorithm for the capacitated open vehicle routing problem. Journal of the Operational Research Society, 58, 1642-1651. Li, L., & Fu, Z. (2002). The school bus routing problem: a case study. Journal of the Operational Research Society, 53, 552–558. López-Sánchez, A., Hernández-Díaz, A., Vigo, D., Caballero, R., & Molina, J. (2014). A multi-start algorithm for a balanced real-world open vehicle routing problem. European Journal of Operational Research, 238, 104113. Li, F., Golden, B., & Wasil, E. (2007). The open vehicle routing problem: Algorithms, large-scale test problems, and computational results. Computers and Operations Research, 34, 2918-2930. Li, J., Wang, R., Li, T., Lu, Z., & Pardalos, P. (2018). Benefit analysis of shared depot resources for multidepot vehicle routing problem with fuel consumption. Transportation Research Part D: Transport and Environment, 59, 417-432. Liu, R., Jiang, Z., & Geng, N. (2014). A hybrid genetic algorithm for the multi-depot open vehicle routing problem. OR Spectrum, 36, 401-421. Lourenço, H., Martin, O., & Stützle, T. (2003). Iterated local search. In: Glover, F., Kochenberger, G. (Eds.), Handbook of Metaheuristics. Kluwer, Boston, 321–353. Montoya-Torres, J., Franco, J., Isaza, S, Jiménez, H., & Herazo-Padilla, N. (2015). A literature review on the vehicle routing problem with multiple depots. Computers and Industrial Engineering, 79, 115–129. Pisinger, D., & Ropke, S. (2007). A general heuristic for vehicle routing problems. Computers and Operations Research, 34, 2403-2435. Russell, R., Chianga, W-C., & Zepeda, D. (2008). Integrating multi-product production and distribution in newspaper logistics. Computers and Operations Research, 35, 1576 – 1588. Schrage, L. (1981). Formulation and structure of more complex/realistic routing and scheduling problems. Networks 11, 229-232. Sariklis, D., & Powell, S. (2000). A heuristic method for the open vehicle routing problem. Journal of the Operational Research Society, 51, 564-573. Şevkli, A., & Güler, B. (2017). A multi-phase oscillated variable neighbourhood search algorithm for a realworld open vehicle routing problem. Applied Soft Computing, 58, 128–144. Shen, L., Tao, F., & Wang, S. (2018). Multi-depot open vehicle routing problem with time windows based on carbon trading. International Journal of Environmental Research and Public Health, 15, 2025. Soto, M., Sevaux, M., Rossi, A., & Reinholz, A. (2017). Multiple neighborhood search, tabu search and ejection chains for the multi-depot open vehicle routing problem. Computers and Industrial Engineering, 107, 211–222. Syslo, M., Deo, N., & Kowaklik, J. (1983). Discrete Optimization Algorithms with Pascal Programs. Prentice Hall, Inc. Tarantilis, C., & Kiranoudis, C. (2002). Distribution of fresh meat. Journal of Food Engineering, 51, 85–91. Tarantilis, C., Ioannou, G., & Kiranoudis, C. (2005). Solving the open vehicle routing problem via a single parameter metaheuristic algorithm. Journal of the Operational Research Society, 56, 588-596. Vidal, T., Crainic, T., Gendreau, M., Prins, C. (2014). A unified solution framework for multi-attribute vehicle routing problems. European Journal of Operational Research, 234, 658–673. Yu, S., Ding, C., & Zhu, K. (2011). A hybrid GA–TS algorithm for open vehicle routing optimization of coal mines material. Expert Systems with Applications, 38, 10568–10573. Zachariadis, E., & Kiranoudis, C. (2010). An open vehicle routing problem metaheuristic for examining wide solution neighborhoods. Computers and Operations Research, 37, 712-723. Zhou, L., Baldacci, R., Vigo, D., & Wang, X. (2018). A multi-depot two-echelon vehicle routing problem with delivery options arising in the last mile distribution. European Journal of Operational Research, 265, 765–778.