The heterogeneous vehicle routing problem with time windows and a limited number of resources

1 This is an Accepted Manuscript of an article published by Elsevier on Engineering Applications of Artificial Intelligence, Volume 94, September 2020, available at: https://doi.org/10.1016/j.engappai.2020.103745. Copyright Elsevier 2020. En idUS: Licencia Creative Commons CC BY-NC-ND To cite this article: Molina, J. C., Salmeron, J. L., Eguia, I., & Racero, J. (2020). The heterogeneous vehicle routing problem with time windows and a limited number of resources. Engineering Applications of Artificial Intelligence, 94, 103745.. 2 The heterogeneous vehicle routing problem with time windows and a limited number of resources Abstract This paper introduces the heterogeneous vehicle routing problem with time windows and a limited number of resources (HVRPTW-LR), a practical extension of the classical vehicle routing problem in which routes to be designed share common scarce resources. The HVRPTW-LR arises when a limited number of resources, such as vehicles, drivers, instruments, and so on, are available but are insufficient to serve all customers in a route planning. Therefore, the route design involves the selection of customers to be visited at each route and resources to be used. Applications to this problem are found in real services companies with high seasonal demand which attend to different types of works and have to decide on how to effectively manage their resources. For designing optimal routes, a hierarchical objective function is considered, maximizing the total number of served customers as the primary objective, and minimizing the travel costs as secondary. A mathematical model of linear programming is introduced to describe and understand all constraints clearly. The problem is first heuristically solved by a semi-parallel insertion heuristic. Then, solutions are improved by a hybrid variable neighborhood descent metaheuristic based on a Tabu Search algorithm for the exploration of the neighborhood and a holding list. Experiments are conducted on numerous sets of benchmark instances from the literature to evaluate the performance of the proposed algorithm. Results show that the algorithm proposed in this paper has a good performance and can be easily applied for solving numerous vehicle routing problem variants from the literature. A new set of benchmark cases for the HVRPTW-LR are also presented and solved. 3 Keywords Vehicle Routing Problem, VRP with time windows, Fixed Heterogeneous Fleet, Variable Neighborhood Descent, Analytical and Numerical Modelling, Performance Analysis. 1 Introduction The world of logistics is changing at breakneck speed. Excess of inventory or overstock is a supply chain problem costing companies large amounts of money every year. Often such inventory is reduced by requiring smaller and more frequent order quantities. Moreover, some logistics sectors are also exposed to high demand seasonality during the year. Thus, in specific months companies have a higher demand than the rest of the year requiring an increase of their logistical resources to complete all customers’ orders. Most distribution companies need resources for the provision of their services. Resources represent sources of supply, support, or aid which are needed to complete the work. In this work, the resources are classified into two types: renewable resources (RRs) and non-renewable or consumable resources (CRs). RRs are needed during the execution of the service action but are not consumed. They are available on a period by period basis, i.e. the available amount is renewed from period to period. Typical examples are the instruments needed to perform the service, vehicles, skilled drivers and others. On the other hand, CRs are consumed during the service action and they are also available on a period by period basis but rather have limited consumption availability for the route planning. Examples are raw materials, wires, electric counters, etc. Obviously, both types of resources are limited and the decision on how to efficiently use them to distribute products or to perform services in a route planning is difficult to take. 4 Fleet managers play an important role in logistics and transportation companies. They are responsible for routing drivers and vehicles in order to keep distributions on schedule and within established budgets. The growing demand for efficiency in all areas of companies induces fleet managers to make their fleets as effective and efficient as possible, using the least amount of resources. Based on this observation, if on certain occasions the total customer demand in a specific time period exceeds the capabilities of the company, some of the received orders from customers may be postponed to the next planning horizon as the cost of increasing the logistical resources is too high. On the other hand, the fleet of vehicles in a company is usually heterogeneous as the company incorporates vehicles of different characteristics over time (Hoff et al., 2010) and also it provides a better adaptation to the customer demand (Yepes and Medina, 2002). The heterogeneous vehicle routing problem with time windows (TW) and a limited number of resources (HVRPTW-LR) is an extension of the classical heterogeneous vehicle routing problem with time windows (HVRPTW) where some customers of a specific planning may not be served due to an excess in the required resources (RRs and/or CRs). In this situation, the HVRPTW-LR requires maximizing the use of those resources, postponing as few orders as possible to the next planning horizon. This paper makes two main scientific contributions. Firstly, the HVRPTW-LR is formally defined by a mixed-integer linear programming model that considers the resource limitations in the route planning. To the best of our knowledge, the HVRPTWLR has not yet been introduced in the literature and better defines real-world vehicle routing problems (VRPs). The second contribution is the development of a methodology for solving the HVRPTW-LR. For this purpose, a hybrid variable neighborhood descent (VND) metaheuristic based on a Tabu Search (TS) algorithm for the exploration of the neighborhood is proposed in this paper. The algorithm introduces new neighbor- 5 hood structure definitions and a holding list (HL) that contains the list of unserved customers. The HL also prevents the algorithm from being trapped at local optimal and achieves to explore a larger search space. To the best of our knowledge, hybrids of VND and TS with HL (VNDTS-HL) are not proposed in the literature for solving VRPs. The rest of the paper is structured as follows: Section 2 provides a review of the literature. The problem description and the mathematical formulation are defined in Section 3. Section 4 presents a detailed description of the main components of the solution approach. Computational experiments are presented in Section 5 and finally, the conclusions are given in Section 6. 2 Literature survey The Heterogeneous Fleet VRP (HF-VRP) is a variant of the classical vehicle routing problem (VRP) that appears when a fleet of vehicles (limited or unlimited), characterized by different capacities and costs is available for the distribution activities (Koç et al., 2016). There are two main divisions of such problems depending on the available fleet. The HF-VRP with unlimited fleet, known as the Fleet Size and Mix VRP (FSMVRP), consists of determining the best fleet composition and routing when there is no limitation on the number of available vehicles of each type. On the other hand, the variant with a limited number of vehicles, called the Heterogeneous VRP (HVRP) consists of optimizing the routes with the available fixed fleet. Two different objective functions have been considered to compute route costs to be minimized. Both consist of minimizing the sum of fixed vehicle costs and dependent routing costs. The latter can be related to the total trip duration without service times or based on the total distance trav- 6 eled. The HF-VRP with TW (HF-VRPTW) appears when additional restrictions are introduced to force customers to be served by a vehicle in a predefined time interval [EWi , LWi], where EWi and LWi are the earliest and latest times to start the service respectively. The HF-VRPTW has been widely studied and a large number of papers have been published in the literature, as can be observed in the surveys of Baldacci et al. (2008), Irnich et al. (2014) and Koç et al. (2016). The majority of these works address the FSMVRP with TW (FSMVRPTW), maybe because it is easier to find feasible solutions than in HVRP with TW (HVRPTW). In the area of HVRPTW, Paraskevopoulos et al. (2008) proposed a two-phase solution framework based on a semi-parallel construction heuristic and a hybridized TS integrated within a new reactive variable neighborhood search (VNS) metaheuristic algorithm. The algorithm solved both HVRPTW and FSMVRPTW. Later, Koç et al. (2015) proposed a hybrid evolutionary algorithm to also solve both problems. The algorithm combines advanced procedures and several metaheuristics principles from Vidal et al. (2012, 2014), such as adaptive large-scale neighborhood search and population search. They obtained some new solutions on benchmark instances outperforming all previous algorithms. Both the FSMVRPTW and HVRPTW have given rise to a multitude of variants which have received particular attention in the last recent years. We refer to the survey of Koç et al. (2016) for further information on HF-VRPTW. All these variants assume that a feasible solution occurs when every customer of a route plan is served by a vehicle on a route satisfying all problem constraints, whereas a solution of the HVRPTW-LR may contain unserved customers. Following this consideration, Lau et al. (2003) proposed the VRP with TW and a limited number of vehicles 7 (m-VRPTW). This variant considers a homogeneous limited fleet that is available at the depot. In contrast to the HVRPTW, a solution to m-VRPTW may contain some unserved customers due to the limited amount of vehicles. The objective in this problem is a hierarchical function where maximizing the total number of served customers is considered as the primary objective and minimizing the total traveled distance as secondary. Our research not only considers the limitation in the number of vehicles but also incorporates limitations in RRs and/or CRs, which are used and/or consumed at customer locations. Thus, the different routes designed in a specific day have to share these scarce resources. In order to solve the aforementioned m-VRPTW, Lau et al. (2003) proposed a TS approach characterized by a HL and a mechanism to force dense packing within a route. The algorithm was tested on a set of nine benchmark instances derived from the clustered set of Solomon (1987) for the VRPTW; they also showed good results for the standard VRPTW. Extending this research, many authors have presented different approaches for the m-VRPTW (Lim and Wang 2004; Lim and Zhang 2007; Wang et al 2008). Jiang et al. (2014) defined the VRP with heterogeneous fleet and time windows to generalize the three variants existing in the literature (m-VRPTW, FSMVRPTW and HVRPTW). They considered a hierarchical objective function where the total number of served customers is considered as the primary objective and minimizing the total number of vehicles and total traveled distance as secondary and tertiary respectively. The authors developed a two-phase TS algorithm using the TS developed by Lau et al. (2003). Relating to the use of limited resources, Molina et al. (2018) present the problem of designing routes in service companies that are responsible for the metrological control 8 of measuring equipments. They solve a rich vehicle routing problem that includes fixed heterogeneous fleet of vehicles, time windows for customers and depot, resource synchronization between tours, driver-customer and vehicle-customer constraints, customer priorities and unserved customers. Hempsch and Irnich (2008) introduced the VRPs with inter-tour resource constraints in which the feasibility of a solution depends on properties of several tours and cannot be decided by considering the individual tours separately. Examples of inter-tour resource constraints are the presence of a limited number of docking stations at depot, a limited number of long tours, the number of stops, the arrival time at the depot, etc. For further information, an extended survey on VRPs with multiple synchronization constraints is due to Drexl (2012). Our problem differs from the previous one in the types of resources and interdependencies considered. In the HVRPTW-LR, the utilization or consumption of a specified resource is always performed at customer sites and there exists a global limitation for all vehicles. In order to design a methodology in accordance with the characteristics of the tackled problem, a hybrid VND metaheuristic based on a TS algorithm for the exploration of the neighborhood is proposed in this paper. Although the algorithm is specifically designed to consider the specific assumptions, objectives and constraints of the HVRPTW-LR, it could also be very useful in other research areas, such as in ShuttleBased Storage and Retrieval Systems (SBS/RS) (see e.g. Carlo et al., 2012; Borovinšek et al., 2017; Lerher et al., 2017, 2015; Dukic et al., 2015; Lerher, 2016; Rosi et al., 2016). 9 3. Problem Description and Mathematical Model The HVRPTW-LR is defined on a graph G= {N, A} with N ={0,1,…,N} as a set of nodes for a route planning, where node 0 represents the depot, and A is a set of arcs defined between each pair of nodes. A set of K heterogeneous vehicles is represented by K={1,…,K} and is available from the depot. With the purpose of containing those nodes that cannot be served in the specific routes planning, a HL is available. Each customer i ∈ {1,…,N} has a fixed demand and a time window for starting the service The sets of RRs, except vehicles, and CRs are represented by R={1,…,R} and C={1,…,C} respectively. Moreover, there exists a latest returning time to the depot for each vehicle k ∈ {1,…,K} and a limited number of each resource. The aim of the HVRPTW-LR is to design a set of vehicle routes making the following assumptions:  Each vehicle serves a single route during the route planning.  All vehicle routes start and end at the depot.  The vehicle fleet is heterogeneous, i.e. vehicles have different capacities, variable costs and latest returning times to the depot.  Customers of a specific planning may not be served.  Each customer has a positive demand which has to be fully satisfied once by exactly one vehicle, if the customer is served in the route planning.  Customers may require RR and/or CR to be served.  The total demand of a route does not exceed the vehicle capacity.  The service of every visited customer starts within its time windows; if the vehicle arrives before the earliest time, it must wait.  A vehicle cannot exceed its latest returning time to the depot. 16         ujuujijiijuij TTSTYLWTTSTYLWC 3, (4d) u uNi Nj ijkikuij DXDQC               4, (4e)               Ni i Nji ij Ni i Nji ij DN TT N D NN TT 1 )1( , ,  (4f) A point of primary importance, for the HVRPTW-LR, is the effective utilization and/or consumption of the available resources. For this reason, two different penalty costs are introduced in the “greedy” function. First, the renewable penalty cost (PCR) gives priority to insert customers which do not modify the RR utilization on the route. Thus, the value of the PCR will be equal to one if the candidate customer for insertion requires the same types of RRs as customers belonging to the partially constructed route. On the other hand, the consumable penalty cost (PCC) only penalizes those customers requiring a great number of CRs. In the PCC, the evaluation criterion for penalizing a customer insertion is based on the third quartile. Although the formulae for percentiles are uniquely defined for continuous random variables, in the discrete case, quartiles divide the customers, previously sorted from lowest to highest consumptions, into four quarters having the same number of customers in each quarter. The third quartile represents the value for which the resource consumption of the 75 per cent of the customers is less than that value. Therefore, those customers with a greater value in the consumable resource (CR) consumption than the quartile will be penalized. As observed, the main idea of this method is to design routes that require the least number of RRs and CRs. Recall that penalty costs only will be considered according to the problem specification. 17 Finally, as the primary objective of the HVRPTW-LR consists of maximizing the total number of served customers, a new metric (Mk) is introduced in Eq. 5 for the routes selection to prioritize the objective of the tackled problem. The main idea is to measure the number of customers that is served in a route using a vehicle of type k. In the case of having multiple routes serving the same number of customers, the metric will choose the route with smaller travel costs. At the end of the procedure, if there are any unassigned customers left, a HL is generated to locate them. 𝑀𝑘=∑𝑋𝑖𝑗𝑘 (𝑖,𝑗)𝜖𝑟𝑜𝑢𝑡𝑒 𝑘 1+𝜌∙∑𝑇𝐷𝑖𝑗∙𝑋𝑖𝑗𝑘 (𝑖,𝑗)𝜖𝑟𝑜𝑢𝑡𝑒 𝑘 (5) 4.2 The optimization process The metaheuristic method we propose, the VNDTS-HL, is originally inspired in VNS. VNS is a metaheuristic proposed by Mladenovic and Hansen (1997) that is composed of three phases: shaking, local search and move. Thus, VNS explores a set of neighborhoods of a current solution, makes a local search and moves to another solution only if there has been an improvement. Since a local optimum for a given type of move (neighborhood structure) is not necessarily so for another, VNS changes the neighborhood structure during the search in order to escape from local optima. The search continues until a local minimum with respect to all neighborhood structures is reached (Hansen et al., 2010). The shaking procedure is a diversification mechanism that consists in perturbing a solution by applying a random move, to provide a new starting point for the local search. Since the shaking mechanism is used to allow a more efficient and effective intensification local search, random moves are not implemented in our approach because they may 18 lead to different route solutions for different executions of the algorithm without changing the parameters, which is not desirable for real distribution companies. For this reason, the deterministic version of VNS, called VND, is considered. Other variants of the VNS are presented in Hansen and Mladenovic (2003). To avoid termination at a local minimum and in order to modulate the intensification and diversification of the search, a TS mechanism, that allows non-improving moves, is incorporated to the local search phase. Moreover, a new mechanism, called RRM, is applied to HVRPTW-LR to escape from local minimum. Finally, a HL is also introduced to tackle unserved customers in the final solution. Unlike evolutionary algorithms, VNDTS-HL is not population-based, and successively moves from one solution to another. This is an advantage to efficiently update the objective function and to check the feasibility in the resource limitations restrictions, without recalculation of the overall solution. Moreover, the choice of VNDTS-HL is motivated by the high complexity of the HVRPTW-LR, which requires algorithms with substantial diversification possibilities as the VND scheme. Furthermore, the utilization of TS results in an intensification of the search, which is of vital importance to find promising solutions. On the other hand, implementation of TS requires specific definition of a neighborhood, which is given by the VND. In the literature, VNS scheme has proved to be very adaptable to VRP variants and it has been successfully applied combined with TS for solving the HVRPTW. Some authors (Molina et al., 2019; Paraskevopoulos et al., 2008) performed the local search of VNS by TS while VNS controls the neighborhood changes. In contrast, Brandão (2006, 2009, 2011) presented a hybrid algorithm where different neighborhood structures are considered in a TS mechanism. Recent research on applications of VNS, or of hybrids of VNS combined with other metaheuristics is diverse and numerous (Bräysy, 2001, 19 2003; Chen et al., 2010; De Armas et al., 2015; Hansen and Perez., 2010; Paraskevopoulos et al., 2008). For further information, Hansen and Perez (2010) present an extensive review on VNS applications. More precisely, the VNDTS-HL starts by defining a set of neighborhood structures Nλ (λ = 1,…, λmax). The iterative process starts from an initial solution s. Then, a local search based on TS is performed to determine a new solution s´ in Nλ. If f(s´) is better than the best solution f(s), then s is replaced by s´, and the search returns to N1, otherwise the search explores the next neighborhood Nλ+1. This is repeated until all neighborhood structures are examined (λ = λmax). At this point, the RRM phase is executed once to escape from a local optimum and the VNDTS-HL restarts from the modified solution until a new local optimum solution s´ is reached. If the obtained solution f(s´) improves f(s), s is replaced by s´ and the algorithm starts from a new initial solution from the set F. The pseudocode of the proposed algorithm is presented in Algorithm 1. Variable Neighborhood Descent Tabu Search algorithm 1 F ← Semi-parallel_Insertion_Heuristic; 2 Define a set of neighborhood structures Nλ, λ=1, 2,…,λmax; 3 For all solutions s of set F do: 4 While (CPU time consumed ≤ γ) do: 5 Set k ← 1 ; RRM ← False; 6 While (k ≤ kmax) do: 7 s’ ← TabuSearch (s, Nλ); 8 If f (s’) improves f (s) then 9 s ← s’; k ← 1; 10 Else 11 If k< kmax then 12 k ← k+1; 13 Else 14 If RRM=False then 15 ŝ ← s; s ← RRMPhase( s) ; k ← 1; RRM ← True; 16 Else 17 If f (s) improves f (ŝ) then 18 ŝ ← s 19 EndIf 20 EndIf 21 EndIf 22 EndIf 23 EndWhile 24 UpdateBestSolution (ŝ); 25 EndWhile 26 EndFor 20 4.2.1 Tabu Search The TS is a widely used metaheuristic due to Glover (1989) that carries out the exploration of the solution space moving successively, in each iteration from one solution s to the best or first improving solution of its neighborhood Nλ(s) even if it may cause a deterioration in the objective function. The central mechanism in the TS is a short-term memory known as the Tabu List (TL) which stores the solutions explored throughout the search or, more commonly, some relevant attributes of these solutions. In this problem type, due to the computational effort, the proposed TL consists of storing the attributes of the performed movements such as the involve nodes and their initial positions in the routes before moving to other solution. To prevent the search from returning to recently visited solutions and to drive the search toward regions of the solution space not yet explored, these selected attributes are declared as tabu and remain in the TL for a specific number of iterations (tabu tenure) unless the aspiration criterion is satisfied. The latter is applied when a move declared as tabu builds a solution that overcomes the best solution found so far. In this case, the mechanism dismisses the TL and the move is accepted. Literature shows several schemes to determine and control the size of the TL during the search (Glover, 1989; Paraskevopoulos et al. 2018). Obviously, the use of a small tabu tenure results in a more effective intensification search since it allows the TS to explore some areas of the solution space allowing the cycling of small periods. In contrast, a large tabu tenure will drive the search toward a part of the solution space that has not been explored yet, usually escaping from a current local optimum solution. This work proposes a scheme for the tabu tenure that provides a balance between diversifica- 21 tion and intensification search strategies using a particular mechanism that has been successfully applied in the works of Molina et al. (2019) and Paraskevopoulos et al. (2008). Initially, the tabu tenure tsize is set equal to a lower value tmin. A diversification mechanism is provided by incrementing at each iteration the tabu tenure in one unit up to an upper bound tmax while no improvement is observed. In contrast, the intensification mechanism is performed when an improvement in the objective function is achieved. For this purpose, the tabu tenure is reinitialized to tmin. The pseudocode of the proposed TS is presented in Algorithm 2. Tabu Search algorithm 1 Given a solution s and a neighborhood structure k; 2 Initialize tabu list TLz of tmin size; 3 elite ← s, counter ← 0, tsize ← tmin, AspirationCondition(s); 4 While (counter ≤ MaxIters) do: 5 Find s’ ∈ Nλ(s) | s subject to tabu & aspiration conditions 6 AllowedSet(s) ← s’; 7 s ← ChooseFirstImproving(AllowedSet(s)); 8 UpdateTabulist(); 9 If f (s) improves f (elite) then 10 elite ← s; counter ← 0, tsize ← tmin, AspirationCondition(elite); 11 Else 12 counter ← counter + 1; 13 If (tsize < tmax) then 14 tsize ← tsize + 1; 15 EndIf 16 EndIf 17 EndWhile 18 s ← elite 4.2.2 Neighborhood structures The VND scheme implemented in this work oscillates between seven neighborhood structures (λmax=7) designed for transforming solutions with the purpose of finding improved configurations. As opposed to other metaheuristics based on VNS, where the neighborhood structures are defined by a single type of move, in our implementation, most of the neighborhood structures are defined by a set of types of moves. 22 Thus, the TS explores at each iteration the neighboring solutions of a current solution s. These are obtained applying the different types of moves defined on a neighborhood Nλ. Then, the first improving solution or the best solution of Nλ(s) is selected as the new current solution. As observed, it is obtained at each iteration by performing only one type of move of the neighborhood Nλ. Moreover, the selected move for finding the new solution in the neighborhood Nλ(s) may be different at each iteration of the TS. The order of the neighborhoods was selected after some experiments considering the impact in the final solution and their cardinality. They are briefly described in the order adopted as follows:  MIXTURE: This neighborhood structure is defined by Relocate (Savelsbergh, 1992) and Exchange (Kindervater and Savelsbergh, 1998) operators which are only applied on pairs of routes (inter-route). The Relocate operator aims to generate a solution by removing a customer from a route and inserting it into another route while the Exchange operator consists of swapping a pair of customers from two different routes.  CROSSINGS: This neighborhood structure is also applied only on pairs of routes and is defined by the CROSS-exchange (Taillard et al., 1997) operator. The CROSSexchange swaps segments of customers between two routes. The different segments may contain an arbitrary number of customers but due to the typically vast number of neighbors that would result, the segment length is limited to three customers. Thus, sets of 1-2, 2-2, 1-3, 2-3 and 3-3 swaps are defined and executed in the listed order.  λ-OPT: This neighborhood structure is applied only on single routes (intra-route) and it aims to generate feasible solutions by examining all possible moves defined in 2OPT (Croes, 1958), and 3-OPT (Lin, 1965) operators. 23  INTERCHANGES: This neighborhood structure is similar to the Mixture type but is applied only on single routes (intra-route).  INSERTIONS: This neighborhood structure is applied only on pairs of routes and is composed by double and triple insertion operators (Brandão, 2011). In a double or triple insertion move, the operation is similar to a single insertion except for removing a segment length of two or three consecutive customers respectively in a route.  GENI INSERTION: This neighborhood structure is applied only on pairs of routes and is only composed of the Generalized Insertion (GENI) operator for the time window variant (Gendreau et al., 1998). It basically consists of removing a customer from a route and inserting it into any two customers of another route. If these customers are not consecutive, different moves in the adjacent customers of the route are performed to make the insertion possible. To reduce the computational effort required for the evaluation of the insertions a p-neighborhood is defined. Thus, the allowable insertions are restricted to between any two customers from the p-neighborhood of the customer to be inserted.  GENI SWAP: This neighborhood structure is applied only on pairs of routes. It consists of removing a pair of customers belonging to two different routes and swapping them but performing the insertions through a GENI operator. It is important to note that the proposed solution scheme includes a HL containing the list of the customers that are not served. The HL is similar to a “phantom” route which participates in the regular local search inducing an extended neighborhood search space for every inter-route operator (Lau et al., 2003). Consequently, there exist some additional moves: 24  Relocate from holding list: Transferring customers from the holding list to an existing route.  Relocate to holding list: Transferring customers from an existing route to the holding list.  Exchange with holding list: Exchanging customers from an existing route with another group of customers in the holding list. Thus, customers of a selected route will be searched completely for possible transfer to/from or for exchange with customers in the HL. The hierarchical cost structure of this problem, favors the transfer of customers from the HL to the routes, increasing the chances of finding high quality solutions to the problem. In addition, the HL favors the procedure to search for better solutions by going through the infeasible solution space (Lim and Zhang, 2007). 4.2.3 Resources Readjustment Mechanism The Resources Readjustment Mechanism (RRM) introduced in this section is only applied to VRPs with limitations on RRs and/or CRs. After the VNDTS-HL terminates at a local optimum, the RRM is executed to modify the solution in order to escape from the local optimum and to diversify the search. The purpose of the RRM is to preserve some favorable features of the current solution in a similar way to the Reformation Phase in Paraskevopoulos et al. (2008) or the Adaptive Diversification Mechanisms in Wei et al. (2014). The RRM is a Ruin and Recreate approach (Schrimpf et al., 2000) which partially destructs a solution and then tries to rebuild it to obtain a new admissible solution. Thus, 25 in the first phase, a part of the solution is destroyed by deleting several chosen customers. It is reasonable to expect that the sequences of nodes of the routes requiring the same types of RRs are relatively good. In the same manner, as the feasibility of the solution is also determined by the CRs consumption, customers with higher demands in CRs are not desired. For these reasons, those customers belonging to the current solution with a greater value in the CR consumption than the third quartile are first penalized and moved to the HL. Then, the utilization of each RR r ∈ {1,…,R} is counted for each vehicle k ∈ {1,…,K} in each route identifying the resource less used on each route (r_mink). In order to change the configuration in the use of RRs on the routes, all customers of a route k using r_mink are extracted from the solution and also moved to the HL. We can reasonably hope that it is possible to find again an admissible solution with a new distribution in the use of RRs and CRs. The second phase consists in the reconstruction of the ruined solution. For this purpose, the semi-parallel construction heuristic, introduced in Section 4.1 is applied but considering all current constructed routes at the same time. More specifically, the list of unserved customers is formed by the customers placed in the HL. Next, all customers in the HL are evaluated in all possible positions between two adjacent customers in all partially constructed routes from the ruined solution using the “greedy” function. The customer with the lowest value in the “greedy” function is placed in the specified position of the selected route and the procedure is repeated until no customer can be inserted in any route. Finally, the VNDTS-HL restarts the optimization procedure from the new solution until a new local optimum solution is reached. Despite the RRM apparently degrading the quality of the solution, in terms of number of served customers, re-starting the optimization procedure from this state can guide the search toward better solutions from the solution space. 32 achieving to serve one more customer. According to figures 2(a)-2(b) it can be seen that the search progress of the TS_imp tends to generate better quality solutions when compared to the classical TS. A possible explanation to this fact is that the intensification mechanism introduced in the TS achieves a better exploitation of the local search in the neighborhood of a new obtained solution, avoiding getting trapped in local optima. (a) C101-LR(D) (b) RC101-LR(D) Figure.-2 Comparison between classical TS and TS implemented. 5.3.2 VNDTS-HL performance comparison in HVRPTW-LR instances 17 19 21 23 25 27 29 31 0 0,5 1 1,5 2 2,5 3 3,5 4 Objective function Time (minutes) Classical TS TS implemented 13 15 17 19 21 23 25 27 0 0,5 1 1,5 2 2,5 3 3,5 4 Objective Function Time (minutes) Classical TS TS implemented 33 This section compares the VNDTS-HL with the VND, proposed by Mladenović and Hansen (1997), and the TS algorithm. The basic VND considers a steepest descent procedure (known also as best improvement local search) for the local search, while the TS algorithm is implemented following the scheme described in section 4.2.1 with the described neighborhood structures. (a) C101-LR(D) (b) RC101-LR(D) Figure.-3 Comparative performance of VND, TS and VNDTS-HL. 17 19 21 23 25 27 29 31 0246810 Objective function Time (minutes) VNDTS-HL TS VND 13 15 17 19 21 23 25 27 0 1 2 3 4 5 Objective Function Time (minutes) VNDTS-HL TS VND 34 Figures 3(a)-3(b) illustrate the search progress of VND, TS and VNDTS-HL. They show that the VNDTS efficiently improves the performance of the VND and TS algorithm in both instances. Although the VND algorithm is the most effective in terms of computational time consumption, the obtained solutions are of lower quality. Contrary, TS requires larger computational times, since it performs a more thorough exploration of the neighborhood at each iteration. It is observed that the VNDTS-HL obtains solutions of higher quality in similar computational times. Moreover, the RRM allows to escape from the local minimum and to diversify the search in reasonable computing times. 5.4 Comparative analysis As there are not benchmark instances for the HVRPTW-LR, the performance of the VNDTS-HL is assessed in test cases from two different VRPTW variants, comparing the results obtained with those presented in the literature. Thus, experimental results for m-VRPTW instances are shown in Section 5.4.1. In Section 5.4.2, a comparative analysis with the best results from the literature corresponding to the HVRPTW benchmark instances is performed. Finally, the test cases for the HVRPTW-LR are solved in Section 5.4.3. 5.4.1 Performance on m-VRPTW test cases Table 2 compares the results obtained by the proposed metaheuristic, denoted as VNDTS-HL with the results obtained by Lim and Zhang (2007) (LZ), Lim and Wang (2004) (LW) and Lau et al. (2003) (LST) in C1 instances. Since LST did not propose results for either R1 or RC1 sets, Tables 3 and 4 only compares the results with those obtained by LZ and LW. These tables illustrate the maximum number of served cus- 35 tomers for each instance as it is the primary objective function of the m-VRPTW. The first line of the table indicates the number of available vehicles while the first column of the table shows the instance category. The numbers in the brackets correspond in the following order to LZ, LW and LST published results. As mentioned above, results for R1 and RC1 instances are only provided for the first two works. For each problem a bold face refers to match with current best-known solution (BKS), whereas a bold face with a ‘*’ indicates new BKS. For C1 instances, the performance of the VNDTS-HL is superior to LW and LST algorithms, with 2 and 7 more customers served respectively. Only LZ approach produces a large number of total served customers than VNDTS-HL. Moreover, the VNDTS-HL achieved 1 new solution with an improvement of 1 customer. For set R1 and RC1 the performance of VNDTS-HL is similar to C1 instances. The average performance of VNDTS-HL is better than LW approach, serving 16 more customers, except for the algorithm of LZ. Nevertheless, our algorithm found 4 new best solutions in RC instances. Resuming, the VNDTS-HL found 5 new best solutions and matched 122 out of 246 instances. Problems of large dimension impose several additional difficulties for the VNDTSHL, especially when the size of the fleet is reduced when it is compared with the total demand of the customers. The adaptation of random mechanisms to the local search phase of the VNDTS-HL, such as those used in the works published by LZ, LW or LST, could also have a rather positive influence in the final solution. They can help the VNDTS-HL to explore the solution space in a more effective way, allowing finding better solutions in m-VRPTW instances. Random mechanisms would allow insertions and/or exchanges of unserved customers from the HL, even if they may cause a decrease in the total number of served customers. However, they would not be as effective 36 when solving other HVRPTW variants. Recall that the VNDTS-HL is executed in a deterministic way, and the main contribution of this work is to provide an effective method to solve the HVRPTW-LR, being flexible to provide reasonably good results for other variants. 5.4.2 Performance on HVRPTW test cases Table 5 summarizes the average results obtained by VNDTS-HL compared to the current state-of-the-art solution approaches for the HVRPTW in Paraskevopoulos et al. (2008) instances. The first line of the table lists the authors using the following abbreviations: LS for Liu and Shen (1999), ReVNTS for Paraskevopoulos et al. (2008) and HEA for Koç et al. (2015). The results show that the VNDTS-HL produced reasonability good results generating 9 BKSs (problems R104A, C102A, C103A, C104A, RC101A, RC102A, RC103A, RC104A, RC204A) and obtaining the same BKS on 4 problems (C101A, C201A, C202A, C204A). Although VNDTS-HL was capable of finding 13 BKSs in a reasonable computational time, the average percentage deviation (0.35%) shows that our algorithm is not as effective for the HVRPTW as HEA (0.12%). However, it is observed that the three algorithms are very accurate with an average deviation that is less than 0.50% for the 24 test problems. The VNDTS-HL presents a worst case performance of 1.63%.. Moreover, the VNDTS-HL even slightly outperforms HEA on C1A and RC1A instances (Table 5), where it produces competitive quality solutions. These results can also be extended to C2A instances, where the difference with HEA is less than 0.01% on average. Furthermore, table 7 shows the results obtained by VNDTS-HL compared to the solutions obtained by Jiang et al. (2014). To the best of our knowledge, we are the first au- 37 thors to apply the algorithm in their instances. The first line of the table lists the authors using the abbreviation of TSJ for Jiang et al. (2014). The first column of the table shows the instance category. Then, the fixed costs (FC), the variable costs (VC), the total costs (TC) and the percentage deviation (% Dev) of the total costs with respect to VNDTSHL are showed. The last rows indicate the minimum, average and maximum percentage deviation over all problem instances and the average values for fixed and variable costs for each method. The VNDTS-HL yields quality solutions with an average deviation of -0.51% and a worst-case performance of 1.59%. Moreover, the VNDTS-HL found 44 new best solutions and outperforms TSJ in almost all RC instances. Looking at the results obtained, on average the VNDTS-HL presents -1.35% of lower vehicle fixed costs than TSJ, which means that our algorithm obtains better fleet compositions to reduce the total costs. On the other hand, this fact adversely affects to the average variable costs, which only yields an increase of 0.40%. 5.4.3 Performance on HVRPTW-LR test cases The details of the results obtained by VNDTS-HL for all HVRPTW-LR instances can be found in Appendix C. 6. Conclusions Real-world VRPs present a variety of constraints and attributes that are not considered in traditional problem model formulations. In this study, a new variant of the HVRPTW is introduced to consider the limitation of resources in the design of the routes (HVRPTW-LR). It considers a fixed fleet of heterogeneous vehicles and the need to use 38 RRs and/or CRs for performing the services at customer sites. Due to the resource limitation (vehicles, RRs and CRs), the HVRPTW-LR includes the possibility of being unable to serve all customers in a route-planning and defines a hierarchical objective function, where maximizing the total number of served customers is considered as the primary objective, while minimizing the total travel cost as secondary. The complexity of the HVRPTW-LR is tackled by a hybrid variable neighborhood descent metaheuristic based on a TS algorithm for the exploration of the neighborhood and a HL (VNDTS-HL). In the first phase, several solutions are obtained with different sequences of parameter values, using a semi-parallel construction heuristic. In the second phase, a hybrid VNDTS-HL is proposed for processing a subset of initial solutions. A HL which holds temporarily unserved customers is integrated within the VNDTS-HL to introduce flexibility for solving instances in which it is impossible to serve all customers due to a limitation in the available resources. To improve the performance of the algorithm, a scheme that provides a balance between diversification and intensification search strategies, is incorporated to the TS. Moreover, a new mechanism, called RRM, is proposed when considering RRs and/or CRs. The positive impact of both mechanisms has been observed through computational experiments. Due to a lack of benchmark data sets for the problem under investigation, computational experiments were carried out on benchmark instances with a limited number of vehicles. Experimental results show that the algorithm is very robust, showing a good performance on a wide range of very different problems. To fully evaluate the effectiveness of the algorithm, new test cases for the HVRPTW-LR were generated and solved. The HVRPTW-LR provides a research gap into some real-world practical applications; i.e. service companies that are responsible for the metrological control of measur- 39 ing equipment at customer sites. They design routes accordingly to the daily limitation on the number of RRs (measure instruments or vehicles) and CRs in the central laboratory. Thus, the different routes designed in a specific day have to share these scarce resources. Moreover, some customers must be served in a predefined TW, early in the workday, in order to not interfere with the opening hours. The home care providers industry or the industrial maintenance services providers sector are other real-life examples that includes the consideration of the attributes presented in the HVRPTW-LR. Managers from these types of companies can use the proposed algorithm to increase the companies´ capabilities (reducing the number of unserved orders in a route planning) and also to reduce the total costs. Performance improvement in operation management in environments with limited resources needs much attention and study. The further research will be the development of more effective algorithms to consider skill requirements in customer orders as well to incorporate a weekly route planning. Table 2. Results for m-VRPTW for C1 instances Nº Veh. 10 9 8 7 6 5 4 C101 100 (100-100-100) 92 (92-92-92) 84 (84-84-84) 75 (75-75-75) 66 (66-66-66) 57 (57-57-57) 47 (47-47-47) C102 100 (100-100-100) 99 (99-99-99) 92 (93-92-92) 83 (84-83-84) 72 (72-72-72) 61* (60-60-60) 49 (49-49-49) C103 100 (100-100-100) 99 (99-99-99) 93 (94-93-93) 84 (84-84-84) 72 (73-72-72) 61 (61-60-60) 49 (49-48-48) C104 100 (100-100-100) 99 (99-99-99) 93 (94-94-94) 84 (84-84-84) 73(73-73-72) 61 (61-61-60) 49 (49-49-48) C105 100 (100-100-100) 92 (93-92-92) 84 (84-84-84) 75 (75-75-75) 66 (66-66-66) 57 (57-57-57) 47 (47-47-47) C106 100 (100-100-100) 93 (94-93-93) 86 (86-86-86) 77 (77-77-77) 68 (68-68-68) 58 (58-58-58) 47 (47-47-47) C107 100 (100-100-100) 93 (94-93-93) 85 (86-85-85) 76 (77-76-76) 67 (67-67-66) 57 (57-57-57) 47 (47-47-47) C108 100 (100-100-100) 96 (98-96-96) 87 (89-87-87) 78 (80-78-78) 69 (70-69-69) 59 (59-59-59) 48 (48-48-48) C109 100 (100-100-100) 99 (99-99-99) 92 (92-92-92) 82 (83-82-82) 72 (72-72-72) 61 (61-61-60) 49 (49-49-48) TOTAL LZ=4878; LW=4859 ; LST=4854; VNDTS-HL=4861 Table 3. Results for m-VRPTW for R1 and RC1 instances Nº Veh. 19 18 17 16 15 14 13 12 11 10 9 8 7 6 5 4 R101 100(100-100) 98(99-98) 97 (98-97) 95 (96-94) 91(93-91) 88(91-88) 86(87-86) 81(82-81) 76(77-75) 70(71-69) 65(65-65) 59(59-59) 53(53-53) 46(46-46) 39(39-39) 32(32-32) R102 100(100-100) 99(99-99) 97(98-98) 95(96-95) 92(94-92) 90(91-89) 88(89-85) 84(85-81) 79(81-78) 74(75-72) 67(68-67) 60(61-59) 53(53-53) 45(45-45) R103 100(100-100) 99(100-99) 98(99-98) 97(98-96) 93(96-93) 89(91-88) 82(85-83) 74(78-74) 67(69-67) 58(60-58) 48(49-48) R104 100(100-100) 98(100-98) 90(92-90) 83(84-83) 73(74-73) 62(64-62) 52(53-52) R105 100(100-100) 98(99-98) 95(97-96) 92(92-92) 87(87-86) 80(81-80) 73(74-73) 66(67-66) 59(59-59) 50(51-51) 41(42-41) R106 100(100-100) 99(100-99) 97(99-96) 93(96-92) 87(90-86) 80(83-80) 74(75-73) 66(67-66) 57(58-57) 48(48-48) R107 100(100-100) 99(100-99) 95(96-95) 88(89-88) 80(81-81) 71(72-71) 62(62-62) 50(50-50) R108 100(100-100) 92(94-93) 85(85-85) 74(75-74) 63(64-63) 53(53-53) R109 100(100-100) 94(97-95) 89(90-89) 82(83-82) 75(75-74) 67(66-67) 58(57-58) 47(47-47) R110 100(100-100) 94(95-94) 86(87-86) 76(79-76) 68(69-69) 59(60-59) 49(48-49) R111 100(100-100) 99(100-99) 94(96-95) 87(89-88) 80(80-80) 70(71-70) 61(61-61) 50(50-50) R112 100(100-100) 98(100-99) 91(92-92) 83(83-83) 73(73-73) 62(63-62) 51(51-51) RC101 100(100-100) 97(98-97) 93(94-93) 89(90-89) 84(84-84) 78(78-78) 71(71-71) 63(63-64) 56(56-56) 49(49-49) 40(40-40) RC102 100(100-100) 95(97-95) 89(92-89) 85(86-85) 79(79-77) 71(72-71) 63(63-63) 54(54-54) 44(44-44) RC103 100(100-100) 97(99-98) 91(94-91) 85(86-85) 77(77-77) 67(68-67) 58(58-57) 47(47-47) RC104 100(100-100) 96(98-95) 89(91-89) 81(82-81) 71(72-72) 61(61-61) 50(50-50) RC105 100(100-100) 99(100-99) 97(98-98) 94(96-94) 89(91-89) 83(85-82) 77(77-77) 70(70-70) 62(61-62) 53(53-53) 42(44-42) RC106 100(100-100) 98(100-99) 94(95-95) 87(89-87) 80(80-80) 71(72-71) 62(63-62) 53(53-53) 43(43-43) RC107 100(100-100) 97(99-97) 91(93-90) 84(84-84) 77*(76-76) 67*(65-66) 56(56-56) 46(46-46) RC108 100(100-100) 94(96-93) 87(88-86) 78(78-78) 70*(68-69) 59*(58-58) 47(47-47) TOTAL LZ=14380; LW=14220 ; VNDTS-HL=14236 12 Dukic, G., Opetuk, T., & Lerher, T. (2015). A throughput model for a dual-tray Vertical Lift Module with a human order-picker. International journal of production economics, 170, 874-881. Gendreau, M., Hertz, A., Laporte, G., & Stan, M. (1998). A generalized insertion heuristic for the traveling salesman problem with time windows. Operations Research, 46(3), 330-335. Glover, F. (1989). Tabu search-part I. ORSA Journal on Computing, 1(3), 190-206. Hansen, P., & Mladenovic, N. (2003). A tutorial on variable neighborhood search. Groupe d'études et de recherche en analyse des décisions, HEC Montréal. Hansen, P., Mladenović, N., & Pérez, J. A. M. (2010). Variable neighborhood search: methods and applications. Annals of Operations Research, 175(1), 367-407. Hoff, A., Andersson, H., Christiansen, M., Hasle, G., & Løkketangen, A. (2010). Industrial aspects and literature survey: Fleet composition and routing. Computers & Operations Research, 37(12), 2041-2061. Ioannou, G., Kritikos, M., & Prastacos, G., (2001). A greedy look-ahead heuristic for the vehicle routing problem with time windows. Journal of the Operational Research Society, 52(5), 523-537. Irnich, S., Schneider, M., & Vigo, D. (2014). Four variants of the vehicle routing problem. Vehicle Routing: Problems, Methods, and Applications, 18, 241-271. Jiang, J., Ng, K. M., Poh, K. L., & Teo, K. M. (2014). Vehicle routing problem with a heterogeneous fleet and time windows. Expert Systems with Applications, 41(8), 3748-3760. Kindervater, G.A.P.,& Savelsbergh, M.W.P., (1998). Vehicle routing: handling edge exchanges. In: Aarts, E.,Lenstra, J.K. (eds.) Local Search in Combinatorial Optimization. pp. 337–360. Wiley, UK. 13 Koç, Ç., Bektaş, T., Jabali, O., & Laporte, G. (2015). A hybrid evolutionary algorithm for heterogeneous fleet vehicle routing problems with time windows. Computers & Operations Research, 64, 11-27. Koç, Ç., Bektaş, T., Jabali, O., & Laporte, G. (2016). Thirty years of heterogeneous vehicle routing. European Journal of Operational Research, 249(1), 1-21. Kritikos, M. N., & Ioannou, G. (2013). The heterogeneous fleet vehicle routing problem with overloads and time windows. International Journal of Production Economics, 144(1), 68-75. Lau, H. C., Sim, M., & Teo, K. M. (2003). Vehicle routing problem with time windows and a limited number of vehicles. European Journal of Operational Research, 148, 559-569. Lerher, T., Ekren, B. Y., Sari, Z., & Rosi, B. (2015). Simulation analysis of shuttle based storage and retrieval systems. International Journal of Simulation Modelling, 14(1), 48-59. Lerher, T. (2016). Travel time model for double-deep shuttle-based storage and retrieval systems. International Journal of Production Research, 54(9), 2519-2540. Lerher, T., Borovinsek, M., Ficko, M., & Palcic, I. (2017). Parametric study of throughput performance in SBS/RS based on simulation. International Journal of Simulation Modelling, 16(1), 96-107. Lim, A., & Wang, F. (2004). A smoothed dynamic tabu search embedded GRASP for m-VRPTW. In Tools with Artificial Intelligence, 2004. ICTAI 2004. 16th IEEE International Conference on (pp. 704-708). IEEE. Lim, A., & Zhang, X. (2007). A two-stage heuristic with ejection pools and generalized ejection chains for the vehicle routing problem with time windows. INFORMS Journal on Computing, 19(3), 443-457. 14 Lin, S. (1965). Computer solutions of the traveling salesman problem. Bell System Technical Journal, The, 44(10), 2245-2269. Liu F.H. & Shen S.Y. (1999). The fleet size and mix vehicle routing problem with time windows. Journal of the Operational Research Society, 50(7),721-732, 1999. Mladenović, N., & Hansen, P. (1997). Variable neighborhood search. Computers & Operations Research, 24(11), 1097-1100. Molina, J. C., Eguia, I., & Racero, J. (2018). An optimization approach for designing routes in metrological control services: a case study. Flexible Services and Manufacturing Journal, 30(4), 924-952. Molina, J. C., Eguia, I., & Racero, J. (2019). Reducing pollutant emissions in a waste collection vehicle routing problem using a variable neighborhood tabu search algorithm: a case study. TOP, 27(2), 253-287. Paraskevopoulos, D. C., Repoussis, P. P., Tarantilis, C. D., Ioannou, G., & Prastacos, G. P. (2008). A reactive variable neighborhood tabu search for the heterogeneous fleet vehicle routing problem with time windows. Journal of Heuristics, 14(5), 425-455. Rosi, B., Grasic, L., Dukic, G., Opetuk, T., & Lerher, T. (2016). Simulation-based performance analysis of automated single-tray vertical lift module. International Journal of Simulation Modelling, 15(1), 97-108. Savelsbergh, M. W., (1992). The vehicle routing problem with time windows: Minimizing route duration. ORSA Journal on Computing, 4(2), 146-154. Schrimpf, G., Schneider, J., Stamm-Wilbrandt, H., & Dueck, G. (2000). Record breaking optimization results using the ruin and recreate principle. Journal of Computational Physics, 159(2), 139-171. Solomon MM (1987) Algorithms for the vehicle routing and scheduling problems with time window constraints. Operational Research 35(2), 254-265. 15 Steuer, R. E. (1986). Multiple criteria optimization: theory, computation, and applications. Wiley. Taillard, É., Badeau, P., Gendreau, M., Guertin, F., & Potvin, J. Y., (1997). A tabu search heuristic for the vehicle routing problem with soft time windows. Transportation Science, 31(2), 170-186. Vidal, T., Crainic, T. G., Gendreau, M., Lahrichi, N., & Rei, W. (2012). A hybrid genetic algorithm for multidepot and periodic vehicle routing problems. Operations Research, 60(3), 611-624. Vidal, T., Crainic, T. G., Gendreau, M., & Prins, C. (2014). A unified solution framework for multi-attribute vehicle routing problems. European Journal of Operational Research, 234(3), 658-673. Wang, X. P., Xu, C. L., & Hu, X. P. (2008). Genetic algorithm for vehicle routing problem with time windows and a limited number of vehicles. In Management Science and Engineering. ICMSE 2008. 15th Annual Conference Proceedings., International Conference on (pp. 128-133). IEEE. Wei, L., Zhang, Z., & Lim, A. (2014). An adaptive variable neighborhood search for a heterogeneous fleet vehicle routing problem with three-dimensional loading constraints. IEEE Computational Intelligence Magazine, 9(4), 18-30. Yepes, V., & Medina, J. R. (2002). Optimización de rutas con flotas heterogéneas y múltiples usos de vehículos VRPHEMSTW. In Actas del V Congreso de Métodos Numéricos en Ingeniería. Madrid (pp. 3-6).