A metaheuristic for a time-dependent vehicle routing problem with time windows, two vehicle fleets and synchronization on a road network

Guillen Reyes, Fernando O.; Gendreau, Michel; Potvin, Jean-Yves Article A metaheuristic for a time-dependent vehicle routing problem with time windows, two vehicle fleets and synchronization on a road network EURO Journal on Transportation and Logistics (EJTL) Provided in Cooperation with: Association of European Operational Research Societies (EURO), Fribourg Suggested Citation: Guillen Reyes, Fernando O.; Gendreau, Michel; Potvin, Jean-Yves (2024) : A metaheuristic for a time-dependent vehicle routing problem with time windows, two vehicle fleets and synchronization on a road network, EURO Journal on Transportation and Logistics (EJTL), ISSN 2192-4384, Elsevier, Amsterdam, Vol. 13, Iss. 1, pp. 1-20, https://doi.org/10.1016/j.ejtl.2024.100143 This Version is available at: https://hdl.handle.net/10419/325215 Standard-Nutzungsbedingungen: Die Dokumente auf EconStor dürfen zu eigenen wissenschaftlichen Zwecken und zum Privatgebrauch gespeichert und kopiert werden. Sie dürfen die Dokumente nicht für öffentliche oder kommerzielle Zwecke vervielfältigen, öffentlich ausstellen, öffentlich zugänglich machen, vertreiben oder anderweitig nutzen. Sofern die Verfasser die Dokumente unter Open-Content-Lizenzen (insbesondere CC-Lizenzen) zur Verfügung gestellt haben sollten, gelten abweichend von diesen Nutzungsbedingungen die in der dort genannten Lizenz gewährten Nutzungsrechte. Terms of use: Documents in EconStor may be saved and copied for your personal and scholarly purposes. You are not to copy documents for public or commercial purposes, to exhibit the documents publicly, to make them publicly available on the internet, or to distribute or otherwise use the documents in public. If the documents have been made available under an Open Content Licence (especially Creative Commons Licences), you may exercise further usage rights as specified in the indicated licence. https://creativecommons.org/licenses/by-nc-nd/4.0/ Contents lists available at ScienceDirect EURO Journal on Transportation and Logistics journal homepage: www.elsevier.com/locate/ejtl A metaheuristic for a time-dependent vehicle routing problem with time windows, two vehicle fleets and synchronization on a road network Fernando O. Guillen Reyesa,c, Michel Gendreau b,c, Jean-Yves Potvin a,c,∗ aDépartement d’informatique et de recherche opérationnelle, Université de Montréal, Montréal, Canada bDépartement de mathématiques et de génie industriel, Polytechnique Montréal, Montréal, Canada cCentre Interuniversitaire de Recherche sur les Réseaux d’Entreprise, la Logistique et le Transport (CIRRELT), Montréal, Canada ARTICLE INFO Keywords: Vehicle routing problem Road network Time-dependent Time windows Transfer points Synchronization Metaheuristic Slack induction by string removals ABSTRACT In this work, we extend the time-dependent vehicle routing problem with time windows on a road network by considering two types of vehicles, large and small, to serve customers. Motivated from city logistics applications, large vehicles are forbidden from the downtown area. Accordingly, goods must be transferred from large to small vehicles to serve downtown customers. This leads to synchronization issues at transfer points, which are special locations without storage capacity. The problem is not a pure two-echelon vehicle routing problem, since customers outside of the downtown area can be served directly by large vehicles. The problem is further compounded by the presence of time-dependent travel times that are defined on the arcs of the road network and are used to model congestion periods. To solve this difficult problem, we propose an adaptation of the Slack Induction by String Removals metaheuristic, which is state-of-the-art for the classical capacitated vehicle routing problem. Computational results on a set of test instances with different characteristics empirically demonstrate the optimization capabilities of this new metaheuristic on a problem which is much more complicated than the capacitated vehicle routing problem. 1. Introduction Although the vehicle routing problem (VRP) has been widely studied for a long time, time-dependent variants have spurred the interest of researchers only recently. Time-dependency is an important issue, since the time to travel from one point to another in a network often depends on the departure time (c.f., rush hours). Furthermore, not only does the time to travel along a path between two customers may change depending on the departure time, but even the best path to use may also change. Thus, recent studies have exploited the additional information available in a road network to account for multiple possible paths between two customers, which is often referred to as the timedependent vehicle routing problem with time windows on a road network (TDVRPTWRN). In this paper, we consider an extension of this problem where both large (black) and small (green) delivery vehicles are involved and where some parts of the road network are forbidden to one type of vehicles or the other. For example, the downtown area is not accessible to large vehicles, whereas areas far from downtown are not accessible to small vehicles (e.g., bicycles). Since the goods to be delivered are initially loaded in large vehicles, a customer located in an area not accessible to them can only be served through a transfer of its demand from a large to a small vehicle. This transfer takes place ∗Corresponding author. E-mail addresses: [email protected] (F.O.G. Reyes), [email protected] (M. Gendreau), [email protected] (J.-Y. Potvin). at special locations with no storage capacity, known as transfer points (TPs). This also leads to synchronization issues between the two types of vehicles at transfer points. In the following, this difficult delivery problem will be referred to as the TDVRPTWRN with transfer points or TDVRPTWTPRN. Our problem needs to be distinguished from problems with intermediate facilities, with or without storage capacity, since there is no facility as such to transfer goods. It must also be distinguished from twoor multi-echelon VRPs where vehicles are organized into a strict hierarchical structure to deliver goods to customers. In our problem, black vehicles can very well serve customers directly, as long as they do not belong to forbidden areas. Our contribution lies in the adaptation of a state-of-the-art metaheuristic for the capacitated VRP (CVRP) for a much more difficult problem that involves two types of vehicles, three types of customers, time-dependent travel times and synchronization between the two types of vehicles to transfer goods at transfer points. As far as we know, this problem has never been addressed in the literature. In the following, Section 2first reviews problems related to ours, namely VRPs on road networks, time-dependent VRPs and VRPs with intermediate facilities. Section 3then precisely describes our problem. https://doi.org/10.1016/j.ejtl.2024.100143 Received 11 March 2024; Received in revised form 30 August 2024; Accepted 13 September 2024 EURO Journal on Transportation and Logistics 13 (2024) 100143 Available online 18 September 2024 2192-4376/© 2024 The Authors. Published by Elsevier B.V. on behalf of Association of European Operational Research Societies (EURO). This is an open access article under the CC BY-NC-ND license ( http://creativecommons.org/licenses/by-nc-nd/4.0/ ). F.O.G. Reyes et al. The original implementation of the metaheuristic Slack Induction by String Removals (SISR) for solving the CVRP is described in Section 4. Then, Section 5introduces time issues that arise in the TDVRPTWTPRN, in particular calculation of time bounds to check insertion feasibility in constant time and synchronization between black and green vehicles at transfer points. Specific modifications to the original SISR implementation that are required to address the TDVRPTWTPRN are presented in Section 6. Then, computational results obtained on test instances derived from a benchmark for the TDVRPTWRN are reported. Finally, a conclusion follows. 2. Literature review The main features of our problem are (1) the consideration of a full road network, (2) time-dependent travel times and (3) intermediate points to transfer goods from one type of vehicle to another. Problems with these characteristics are briefly reviewed in the following. 2.1. VRPs on road networks In many VRPs, it is implicitly assumed that the best path between two customers (or customer and depot) can be uniquely identified in the underlying road network. Then, a so-called customer-based graph is constructed, with nodes that correspond to the customers plus the depot and with an arc between each pair of nodes that stands for the corresponding best path. However, it is often the case that a single best path cannot be identified a priori, for example when multiple attributes (distance, time, cost) are associated with each road segment in the road network. That is, the best path is not necessarily the same depending on the attribute considered. Furthermore, trade-offs between different attributes are discarded if a single path is considered. Accordingly, working with a customer-based graph reduces the solution space and may lead to an overestimation of the optimum cost. To address this issue, two approaches are reported in the literature: representing the road network as a multigraph or working directly with the full road network. In multigraphs, customer-based graphs are extended by introducing parallel arcs between each pair of nodes, where each arc stands for a different path that is worth considering in the underlying road network. To the best of our knowledge, the first use of a multigraph for a multi-attribute vehicle routing problem, namely a dial-a-ride problem, is reported in Garaix et al. (2010). The authors in Ben Ticha et al. (2018), who surveyed a number of papers based on multigraphs (up to 2018), indicate that the latter provide average benefits between 5% and 15% when compared to solutions obtained on customer-based graphs. One difficulty with multigraphs comes from their construction, since the set of parallel arcs between two nodes can be large and may be difficult to obtain (multicriteria shortest path problems must be solved when more than one attribute is associated with each road segment). One may settle for only a subset of all possible parallel arcs, but at the expense of a reduced solution space. On the other hand, working with the full road network preserves the entire solution space. A few studies compare the use of multigraphs versus road-network graphs with somewhat different observations. In Ben Ticha et al. (2017,2019), the authors empirically demonstrate the benefits of using a multigraph representation versus a customer-based graph for a branch-and-price algorithm and an adaptive large neighborhood search (ALNS) applied to a bi-objective VRP with time windows (VRPTW) that accounts for both travel time and travel cost. However, they note the considerable computation times needed to compute the multigraph. In the same context, the authors in Letchford et al. (2014) propose to work directly with the road network since the pricing problem in the branch-andprice algorithm can be solved more quickly, while the construction of a multigraph is avoided. Also, working with a road network appears to be more natural and straightforward. These observations must be contrasted with those in Ben Ticha et al. (2019), where a comparison between a branch-and-price algorithm applied to a road network and a multigraph representation for a bi-objective VRPTW (travel time, travel cost) shows that both approaches are competitive, with a slight advantage for multigraphs on the more realistic instances. The reader interested in those issues, as well as some additional ones, concerning VRPs defined on multigraphs and road networks, is referred to the survey in Ben Ticha et al. (2018). In the latter survey, the authors also note that identifying the best path between two customers in a road network is not that simple, even if a single attribute is considered. This is the case in particular for time-dependent travel times or travel costs, as discussed in the next section. 2.2. Time-dependent VRPs In TDVRPTWs, the travel time along an arc depends on the departure time from the origin node. When the objective function is not related to time, like the classical minimization of total traveled distance, the shortest path in the road network between each pair of customers (or customer and depot) can be computed to construct a customer-based graph, given that the shortest path does not change over time. The variations in travel times are then accounted for along these shortest paths. This is the approach used in the first studies about TDVRPTWs. To the best of our knowledge, the first work that addressed timedependency (although without time windows at customers) using a customer-based graph is found in Beasley (1981). In this work, the time horizon is divided into periods with a different travel time matrix for each period. This is equivalent to defining a step function to model the travel time at different periods between any given pair of nodes. A similar approach is proposed in Malandraki and Daskin (1992) for the TDVRPTW. Since the travel time is constant within a period, but may abruptly change from one period to the next, the two previous models do not satisfy the First-In-First-Out (FIFO) property, where it is required that a vehicle traveling earlier on an arc must arrive at the destination node earlier than any other vehicle traveling later on the same arc. A different model is proposed in Hill and Benton (1992), where each node is assigned a speed at a given period of time, which can be interpreted as the average speed around the node. Then, the travel time on a given arc between two nodes is based on the average speed around these two nodes. Once again, since the travel time on a given arc is constant within a period, the FIFO property is not satisfied. A model that satisfies the FIFO property was finally proposed in Ichoua et al. (2003). Here, the authors use a step function to model speed (rather than travel time) at different time periods. That is, the speed is constant within a given time period, but may change from one period to the next. The travel time along an arc is then computed by taking into account the new speed when the time boundary between two periods is crossed. Thus, every vehicle that travels along the same arc within the same period has the same speed and the speed of every vehicle changes similarly when the boundary between two given periods is crossed. This way to model time-dependency has been largely adopted in the following years to solve TDVRPTWs using exact methods and metaheuristics (Donati et al., 2008;Dabia et al.,2013;Pan et al.,2021a,b). In a few cases, continuous functions with special characteristics to satisfy the FIFO property have also been used to model time-dependent travel times (Haghani and Jung,2005;Balseiro et al.,2011). For a detailed literature review on time-dependent VRPs using customer-based graphs (up to 2015), the reader is referred to Gendreau et al. (2015). When more realistic objective functions based on travel times or travel costs are considered, a new difficulty occurs that prevents the use of customer-based graphs. That is, not only does the travel time change along a given path in the road network depending on the departure time from the origin node, but even the best (fastest or least-cost) path may change. This is accounted for by using either a multi-graph representation or the full road network. In the case of multi-graphs, EURO Journal on Transportation and Logistics 13 (2024) 100143 2 F.O.G. Reyes et al. parallel arcs between any given pair of customers stand for different best paths in the underlying road network depending on the departure time from the origin customer node (Ben Ticha et al.,2017,2019). In more recent works, the authors tend to use the full road network to handle these multiple alternative paths (Ben Ticha et al.,2019,2021; Gmira et al.,2021), for reasons already mentioned in Section 2.1. 2.3. VRPs with intermediate facilities Due to the presence of transfer points in our problem, we provide here an overview of the literature on VRPs with intermediate facilities, which are referred to as satellites, hubs, transshipment points or crossdocks. They all represent intermediate points where goods can be transferred while they move from their origin to their destination. In Speranza et al. (2016), the authors provide a survey about intermediate facilities in freight transportation, while a survey dedicated to cross-docking is found in Van Belle et al. (2012). In Speranza et al. (2016), the problems are divided into two classes, that is, two-echelon VRPs (2E-VRPs) and pickup and delivery problems with cross-docks (PDPCDs). With regard to 2E-VRPs, the intermediate facilities are called satellites and have some storage capacity. At the first-level or echelon, vehicles carry goods from a depot to satellites, while at the second-level goods are transported by other vehicles from satellites to customers. Typically, a strict hierarchy is observed, that is, direct deliveries from the depot to customers is forbidden. In this survey, no work requires synchronization between vehicles at satellites. In the case of the surveyed PDPCDs, however, cross-docks have no or little capacity and synchronization is required. Two works are worth mentioning, since there is no real intermediate facility, only transfer points (like in our work). In Bouros et al. (2011), transfers can take place at arbitrary locations and the vehicle that arrives first at the transfer point waits as long as necessary to transfer goods to the other vehicle, although a waiting penalty is incurred. In Minic and Laporte (2006), transshipment points are predetermined and synchronization is achieved by setting a time window at these transfer points. In Crainic et al. (2009), a 2E-VRP is proposed in the context of city logistics, where it is called a two-tier city logistics system. City Distribution Centers (CDCs), located at the outskirts of the city, form the first tier of the system where freight is sorted and consolidated. The second tier of the system is made of satellites located close to or within the city-center area. Different vehicle fleets are used to transport freight from CDCs to satellites and from satellites to customers. In particular, vehicles of the second tier must be adapted for utilization in dense city zones. Since it is assumed that satellites operate according to a cross-dock transshipment operational model, vehicle synchronization is required. That is, vehicles of the first and second tier must meet at satellites at a given time, with very short waiting time permitted. This work proposes only a modeling framework and no algorithmic solution is developed. Different exact and heuristic algorithms were later proposed in Crainic et al. (2011) and Perboli et al. (2010,2009) to solve variants of the initial model (e.g., no synchronization, storage capacity at satellites). A recent work in Jia et al. (2022) addresses a two-commodity 2E-VRP with synchronization at satellites. Two types of vehicles are considered at the first level, one for each commodity, as opposed to the second level where only one type of vehicles is considered. Synchronization is only established between the two types of first-level vehicles, which have to meet at satellites to favor efficiency at the second level. The problem is solved with ALNS where, at each iteration, destroy and a repair operators are applied to the second-level tours, and a reconstruction procedure is then applied to the first-level tours in case of infeasibility, followed by an improvement procedure. The 2E-VRP reported in Anderluh et al. (2017) is particularly interesting because is shares similarities with our problem. In this work, the first-level vehicles are called vans and the second-level vehicles are called bicycles. Similarly, there are two classes of customers depending on their location: customers located at the city center are called bike-customers, while customers outside of the center are called van-customers. Transfers take place at satellites, with no storage capacity. These satellites are located at the boundary of the city center (a van can cross the city center, but a penalty is incurred). In the proposed heuristic methodology based on GRASP and path-relinking, the second-level tours are constructed before the first level tours. In this way, information about the arrival times of bikes at satellites can be used to construct the first-level tours and account for synchronization. The authors in Grangier et al. (2016) address a similar problem called the two-echelon multi-trip VRP with satellite synchronization. The proposed methodology also constructs second-level tours before first-level tours to produce the initial solution. Then, an ALNS is applied with features aimed at improving synchronization, like a destroy operator that removes trips with the worst synchronization. In Kafle et al. (2017), the authors describe a crowdsourced system for urban parcel deliveries, where truck-carriers visit intermediate facilities called relay points and where parcels are transferred to pedestrians or cyclists that are close to the end customers. However, if no pedestrian or cyclist is available, the truck can perform the deliveries itself. In the proposed system, the delivery tasks, as well as candidate relay points, are broadcast on-line. Then, pedestrians and cyclists bid for these delivery tasks. Thus, a bid selection problem must be solved in addition to the routing problem. This is done in both cases with a tabu search. Another similar crowdsourced system is described In Sampaio et al. (2020), where goods can be dropped at transfer points to be picked up later by other vehicles (i.e., transfer points have storage capacity). 3. Problem definition This paper addresses the time-dependent vehicle routing problem with time windows and transfer points on a road network or TDVRPTWTPRN. As previously mentioned, two types of vehicles with different capacities are considered: black (large) and green (small) vehicles. The two sets of vehicles are denoted 𝐾𝐵and 𝐾𝐺, respectively. There are also three types of customers: black customers that can be served by black vehicles only; green customers that can be served by green vehicles only; and neutral customers that can be served by both types of vehicles. A road network in this context is a directed graph 𝐺= (𝑉 , 𝐴), where 𝑉is the set of nodes of cardinality 𝑛and 𝐴the set of arcs or road segments. The set of nodes is then partitioned as follow: -𝐷= {𝑑𝑏, 𝑑𝑔}is the set of depots with 𝑑𝑏the depot for black vehicles and 𝑑𝑔the depot for green vehicles; -𝐶𝐵is the set of black customers of cardinality 𝑛𝐵; -𝐶𝐺is the set of green customers of cardinality 𝑛𝐺; -𝐶𝐸is the set of neutral customers of cardinality 𝑛𝐸; -TP is the set of transfer points of cardinality 𝑛TP; -RJ is the set of road junctions (i.e., any node that is not a depot, a customer or a transfer point). The TDVRPTWTPRN can be characterized as follow: •Each customer 𝑖has a demand (load) 𝑑𝑖and a service (dwell) time 𝑠𝑡𝑖; •Each customer 𝑖has a time window [𝛼𝑖, 𝛽𝑖]to constrain the service start time. If a vehicle arrives at customer 𝑖before 𝛼𝑖, then it must wait until 𝛼𝑖to start the service. On the other hand, a vehicle cannot arrive after 𝛽𝑖; •The demand of all customers is assumed to be loaded into black vehicles at the start; •Each black vehicle performs a single route that starts and ends at the black depot; each green vehicle performs a single route that starts and ends at the green depot; EURO Journal on Transportation and Logistics 13 (2024) 100143 3 F.O.G. Reyes et al. •The black and green depots have a time window [0, 𝑇 ], where 𝑇 is the end of the time horizon; all vehicles must be back at their depot before or at time 𝑇; •A green vehicle can serve only one customer at a time. Thus, after departing from the green depot, a green vehicle visits repeatedly a transfer point, to get a load, followed by the corresponding customer to deliver that load. At the end, the vehicle returns to the green depot. •A black vehicle has a capacity 𝑄𝑏that allows it to serve many customers and to carry loads that will be transferred to green vehicles; •Black and green vehicles have different speeds. Thus, each arc (𝑖, 𝑗) ∈ 𝐴is associated with two time-dependent travel speed functions 𝑣𝑏 𝑖,𝑗 (𝑡)and 𝑣𝑔 𝑖,𝑗 (𝑡)for black and green vehicles, respectively. •Transfer points are fixed locations without storage capacity, where loads are transferred from a black vehicle to one or more green vehicles (we assume, without loss of generality, that the transfer time is null). Thus, synchronization between the two types of vehicles is required at transfer points, which may lead to waiting time for a black vehicle (if one or more green vehicles that must receive loads from the black vehicle have not yet arrived at the transfer point) or green vehicles (if the black vehicle has not yet arrived at the transfer point). A black vehicle can visit the same transfer point multiple times along its route; the same is true of green vehicles. Each visit to transfer point tp ∈TP in a route is represented by a copy which is unambiguously denoted 𝑡𝑝𝑘 𝑗, where 𝑘is a vehicle and 𝑗is the copy (or visit) index. That is, copy 𝑡𝑝𝑘 𝑗corresponds to the 𝑗th visit of transfer point tp in the route of vehicle 𝑘; •Each black customer is served directly and exactly once by a black vehicle; each green customer is served exactly once by a green vehicle after its demand has been transferred from a black vehicle at a transfer point; each neutral customer is served exactly once, either directly by a black vehicle or by a green vehicle after its demand has been transferred from a black vehicle at a transfer point; •The objective is to determine routes of minimal total duration such that all customers are served and all constraints are satisfied. Fig. 1 shows a typical solution, with one black route starting from the black depot (square). This route is identified by arcs (1) to (10). At the first transfer point 𝑡𝑝1(triangle), there is a connection with the green route with arcs identified with broken lines. This route starts at the green depot (square), gets a load from the black vehicle at 𝑡𝑝1, delivers the load to neutral customer 𝑛𝑐1, gets another load from the same black vehicle at the second transfer point 𝑡𝑝2, delivers the load to green customer 𝑔𝑐1and returns to the green depot. The arcs of the second green route are identified with dotted lines. This small route starts at the green depot, gets a load from the black vehicle at 𝑡𝑝2, delivers the load to green customer 𝑔𝑐2and returns to the green depot. It should be noted that the black vehicle transfers two loads, one for each green vehicle, at transfer point 𝑡𝑝2. 4. SISR for the CVRP The methodology for solving our problem is the Slack Induction by String Removals (SISR) metaheuristic (Christiaens and Vanden Berghe, 2020), which has proven to be state-of-the-art for the CVRP. It is based on the ALNS framework, initially proposed in Ropke and Pisinger (2006). Accordingly, SISR also exploits the ruin-and-recreate principle where, at each iteration, a number of nodes are first removed from the routes of the current solution (ruin) and reinserted (recreate) to produce a new solution. A simulated annealing-based criterion is then applied to decide if the new solution should be accepted or not as the current solution. This framework has been used with success to address many different vehicle routing problems, see Pisinger and Ropke (2019). In the following, we first describe the original SISR metaheuristic for the CVRP. This description is quite detailed to allow the reader to fully understand later the modifications that we have performed to this algorithm to address our problem. The basic idea of SISR is to remove strings of consecutive customers from a solution, with at most one string removed from any given route. Algorithm 1shows the pseudo-code of SISR for the CVRP. First, two parameter values are set : 𝐿𝑚𝑎𝑥, which is used to determine the maximum length of a string, and 𝑐, which corresponds to the average number of customers to be removed from a solution. By appropriately setting these values, many strings of small length or only a few strings of large length can be removed. Given that simulated annealing principles guide the search through an exponential cooling schedule, the starting temperature 𝜏0, final temperature 𝜏𝑓,𝜏0> 𝜏𝑓>0, and number of iterations 𝑓are defined, with the current temperature 𝜏initially set to 𝜏0, see statements 2 and 3. Then, the cooling factor 𝜌is defined in statement 4 in such a way that 𝑓ruin-and-recreate iterations are performed. An adjacency list 𝑎𝑑𝑗(𝑖)is then created for each customer 𝑖in statement 5. This list contains all customers ordered from closest to farthest in distance from 𝑖, with 𝑖as its first element. This adjacency list is used to favor the removal of strings that are relatively close to each other, even if they come from different routes. Before proceeding with the main loop, an initial solution is created in step 6 in a straightforward way, by creating an individual route for each customer. This initial solution becomes the current solution 𝑠as well as the best solution known to date 𝑠𝑏𝑒𝑠𝑡. The main loop corresponds to statements 8 to 18. At each iteration, a ruin operator and a recreate operator are applied to a copy 𝑠of current solution 𝑠, see statements 9 and 10. Note that the set 𝐴−is used to store the removed customers. The resulting solution of the ruin-and-recreate process 𝑠is accepted as the new current solution 𝑠 if it satisfies the simulated annealing-based criterion in statement 11 (see Christiaens and Vanden Berghe (2020)). It also replaces 𝑠𝑏𝑒𝑠𝑡 if it is the best solution found thus far. The current temperature 𝜏is then updated before the next iteration starts. After 𝑓iterations of the main loop, the whole procedure stops and returns the best solution found. Algorithm 1 SISR for CVRP 1: Set 𝐿𝑚𝑎𝑥 and 𝑐 2: Set 𝜏0,𝜏𝑓and 𝑓 3: 𝜏←𝜏0 4: 𝜌←(𝜏𝑓 𝜏0)1∕𝑓 5: Generate adjacency list 𝑎𝑑𝑗(𝑖)for each customer 𝑖 6: Generate initial solution 𝑠(with set of routes 𝑅𝑠) 7: 𝑠𝑏𝑒𝑠𝑡 ←𝑠 8: for 𝑓iterations do 9: 𝑠, 𝐴−←𝑅𝑢𝑖𝑛(𝑠) 10: 𝑠←𝑅𝑒𝑐𝑟𝑒𝑎𝑡𝑒(𝑠, 𝐴−) 11: if 𝐶𝑜𝑠𝑡(𝑠)<(𝐶𝑜𝑠𝑡(𝑠) − 𝜏ln(𝑈(0,1))) then 12: 𝑠←𝑠; 13: end if 14: if 𝐶𝑜𝑠𝑡(𝑠)< 𝐶𝑜𝑠𝑡(𝑠𝑏𝑒𝑠𝑡)then 15: 𝑠𝑏𝑒𝑠𝑡 ←𝑠 16: end if 17: 𝜏←𝜌𝜏 18: end for 19: Return 𝑠𝑏𝑒𝑠𝑡 4.1. Ruin In the ruin procedure described in Algorithm 2, the maximum length of a string to be removed 𝑙𝑚𝑎𝑥 𝑠is first set to the minimum of 𝐿𝑚𝑎𝑥 and the average number of nodes in a route of the current solution EURO Journal on Transportation and Logistics 13 (2024) 100143 4 F.O.G. Reyes et al. Fig. 1. An example of a solution to the TDVRPTWTPRN. Algorithm 2 Ruin(𝑠) 1: 𝑙𝑚𝑎𝑥 𝑠←min{𝐿𝑚𝑎𝑥, 𝐴𝑣𝑔𝑁𝑜𝑑𝑒𝑠𝐼𝑛𝑅𝑜𝑢𝑡𝑒𝑠(𝑠)} 2: Calculate 𝑛𝑚𝑎𝑥 𝑠with 𝑙𝑚𝑎𝑥 𝑠and 𝑐 3: 𝑛𝑠←⌊𝑈(1, 𝑛max 𝑠+ 1)⌋ 4: 𝑅−←∅ 5: 𝐴−←∅ 6: 𝑠←𝑠 𝑅𝑠←𝑅𝑠 7: Select randomly a seed customer 𝑖𝑠𝑒𝑒𝑑 in 𝑠 8: for 𝑖∈𝑎𝑑𝑗(𝑖𝑠𝑒𝑒𝑑 )and |𝑅−|< 𝑛𝑠do 9: 𝑟←route of customer 𝑖 10: if 𝑖∉𝐴−and 𝑟∉𝑅−then 11: 𝑙max 𝑟←min{𝑙max 𝑠,|𝑟|} 12: 𝑙𝑟←⌊𝑈(1, 𝑙max 𝑟+ 1)⌋ 13: 𝑅𝑢𝑖𝑛𝑂𝑝 ←𝑅𝑎𝑛𝑑𝑜𝑚(String,Split-String) 14: 𝐴−←𝐴−∪𝑅𝑢𝑖𝑛𝑂𝑝(𝑠, 𝑟, 𝑙𝑟, 𝑖) 15: 𝑅−←𝑅−∪ {𝑟} 16: if 𝑟is empty then 17: 𝑅𝑠←𝑅𝑠∖{𝑟} 18: end if 19: end if 20: end for 21: Return 𝑠,𝐴− Algorithm 3 Recreate(𝑠,𝐴−) 1: Sort(𝐴−)⊳Recreate 2: for 𝑖∈𝐴−do 3: 𝑝𝑏𝑒𝑠𝑡 ←𝑁𝑈𝐿𝐿;𝐶𝑜𝑠𝑡𝐼𝑛𝑠𝑒𝑟𝑡𝑏𝑒𝑠𝑡 ←∞ 4: for 𝑟∈𝑅𝑠and 𝑟feasible with insertion of 𝑖do 5: for 𝑝𝑟in 𝑟do 6: if 𝑈(0,1) <1 − 𝛾then 7: if 𝑝𝑏𝑒𝑠𝑡 =𝑁𝑈𝐿𝐿 or 𝐶𝑜𝑠𝑡𝐼𝑛𝑠𝑒𝑟𝑡(𝑖, 𝑝𝑟)< 𝐶𝑜𝑠𝑡𝐼𝑛𝑠𝑒𝑟𝑡𝑏𝑒𝑠𝑡 then 8: 𝑝𝑏𝑒𝑠𝑡 ←𝑝𝑟 9: 𝐶𝑜𝑠𝑡𝐼𝑛𝑠𝑒𝑟𝑡𝑏𝑒𝑠𝑡 ←𝐶𝑜𝑠𝑡𝐼𝑛𝑠𝑒𝑟𝑡(𝑖, 𝑝𝑟) 10: end if 11: end if 12: end for 13: end for 14: if 𝑝𝑏𝑒𝑠𝑡 =𝑁𝑈𝐿𝐿 then 15: 𝑅𝑠←𝑅𝑠∪ {new empty route 𝑟} 16: 𝑝𝑏𝑒𝑠𝑡 ←first position in 𝑟 17: end if 18: Insert 𝑖in position 𝑝𝑏𝑒𝑠𝑡 19: end for 20: Return 𝑠 𝐴𝑣𝑔𝑅𝑜𝑢𝑡𝑒𝑁𝑜𝑑𝑒𝑠(𝑠), see statement 1. Then, in statement 2, the maximum number of removed strings 𝑛𝑚𝑎𝑥 𝑠is calculated using 𝑙𝑚𝑎𝑥 𝑠and 𝑐, see the exact formula in Christiaens and Vanden Berghe (2020). The actual number of removed strings 𝑛𝑠is chosen from a continuous uniform distribution defined between 1 and 𝑛𝑚𝑎𝑥 𝑠+ 1, as indicated in statement 3. The set of ruined routes 𝑅−and the set of removed customers 𝐴−are then initialized with the empty set. After creating a copy 𝑠of the current solution 𝑠, a random seed customer 𝑖𝑠𝑒𝑒𝑑 is chosen and its adjacency list is processed (from closest to farthest customers) until all customers have been considered or the number of ruined routes is reached, see the main loop in statements 8–20. Note that the number of ruined routes is the same as the number of removed strings 𝑛𝑠, since each string is removed from a different route. If the current customer 𝑖in the adjacency list of 𝑖𝑠𝑒𝑒𝑑 has not been previously removed and if the route 𝑟that serves 𝑖has not been previously ruined (statement 10) then the ruin operator is applied to route 𝑟. In statements 11 and 12, the actual length of the removed string 𝑙𝑟is chosen from a continuous uniform distribution defined between 1 and 𝑙max 𝑟+ 1, where 𝑙max 𝑟is the minimum of 𝑙max 𝑠and cardinality of 𝑟(since the length of the removed string cannot exceed the number of nodes in 𝑟). Then, a random choice between two ruin operators takes place in statement 13. These operators are: •String: A random string of length 𝑙𝑟that contains the current customer 𝑖in the adjacency list of 𝑖𝑠𝑒𝑒𝑑 is removed from route 𝑟. This is illustrated in Fig. 2(a) for a string of length four with the gray node 𝑖3as current customer 𝑖; •Split-String: A random string of length 𝑙𝑟+𝑚that contains the current customer 𝑖in the adjacency list of 𝑖𝑠𝑒𝑒𝑑 is chosen in route 𝑟(where the procedure to select a value for 𝑚is precisely described in Christiaens and Vanden Berghe (2020)). Then, a random substring of 𝑚consecutive customers within the chosen string is kept in the route, so that only 𝑙𝑟customers are removed. The substring of length 𝑚cuts the string of length 𝑙𝑟+𝑚in two parts, unless the substring is at the very beginning or very end of the string of length 𝑙𝑟+𝑚. An example is provided in Fig. 2(b) for a string of length five with the gray node 𝑖3as current customer 𝑖. In this example, 𝑚= 2, so that only three customers are removed from the route. The removed customers are then added to 𝐴−and the ruined route to 𝑅−in statements 14 and 15. If route 𝑟becomes empty, then it is deleted from the set of routes in the solution, as indicated in statements 16 and 17. At the end, the ruined solution 𝑠and the set of removed customers are returned. 4.2. Recreate A new complete solution is then produced with the recreate operator by reinserting the removed customers. This operator is described EURO Journal on Transportation and Logistics 13 (2024) 100143 5 F.O.G. Reyes et al. Fig. 2. Examples of the two SISR ruin operators : (a) String (b) Split-String. in Algorithm 3. In statement 1, the removed customers in 𝐴−are first sorted using a sorting criterion chosen with a particular distribution probability among: random, decreasing demand, increasing distance from the depot, decreasing distance from the depot. Based on the chosen order, the customers in 𝐴−are considered one by one and reinserted in the set of routes 𝑅𝑠of solution 𝑠, see the main loop in statements 2–19. Each insertion place in each route that can accommodate the demand of customer 𝑖is considered and the best encountered insertion place 𝑝𝑏𝑒𝑠𝑡 is identified. It should be noted, however, that the chosen insertion place is not necessarily the best one among all feasible insertion places, due to blinks that correspond to a small probability 𝛾 of skipping a position, see statement 6. In particular, if the best position is skipped then only the second best position can be chosen (as long as this position is not skipped too). Statements 14–17 cover the situation when no feasible insertion place is found for customer 𝑖. In this case, a new route is created for that customer. At the end, the recreated solution 𝑠is returned. 5. Time dependency The SISR metaheuristic for the CVRP, as described in the previous section, needs to be considerably modified to address the TDVRPTWTPRN. In particular, the time dimension must now be taken into account; furthermore, routes are not independent anymore since they interact through transfer points. In this section, we introduce the basics of our time-dependent travel time model and explain how time bounds can be derived at each node along the routes of black and green vehicles. 5.1. Time-dependent travel times The IGP model proposed in Ichoua et al. (2003) is used to model time dependency. In this model, the time horizon [0, 𝑇 ]is partitioned into a number 𝑙of time periods [0, 𝑡1),[𝑡1, 𝑡2), ...., [𝑡𝑙−2, 𝑡𝑙−1),[𝑡𝑙−1, 𝑇 ], where 𝑡1,𝑡2, ..., 𝑡𝑙−1 are time boundaries between two periods. For any given arc, a travel speed is associated with each period and a speed change occurs when a vehicle crosses a time boundary. The algorithmic procedure to compute the travel time along an arc for a given departure time based on this model is provided in Ichoua et al. (2003). Although speed is modeled as a step function of time, the corresponding travel time function is a piecewise linear function. Fig. 3 shows an example of a travel speed function on a given arc (𝑖, 𝑗)and the corresponding travel time function, assuming that the arc is of length 4. 5.2. Dominant shortest-path structure The dominant shortest-path structure (DSPS), as described in Gmira et al. (2021), is useful to quickly identify the fastest path between any given pair of nodes (either customers, depots or transfer points) in the road network for any given departure time. First, a number of good paths between two given nodes 𝑖and 𝑗are identified by applying a time-dependent Dijkstra’s algorithm (Gmira et al.,2021) using different departure times from 𝑖, like time boundaries between two periods. The travel time function of each one of those paths is obtained by combining the travel time functions of all arcs along that path (which also produces a piecewise linear function). Fig. 4 shows an example of a DSPS based on three different fastest paths between two nodes. In this figure, the arrival time is represented as a function of the departure time, so that the corresponding travel time is simply the difference between arrival and departure times. Since the IGP model satisfies the FIFO property, this piecewise linear function is non decreasing. By overlapping the three paths, it is possible to identify the fastest among the three paths for any given departure time. It is worth noting that the DSPS is exact only if the time-dependent Dijkstra’s algorithm is applied with a sufficiently large number of departure times to cover all fastest paths between two nodes, which is rarely the case in practice. But better accuracy is obtained with more departure times. Two different travel time functions are associated with each arc, depending if a black or a green vehicle follows that arc, because they do not have the same speed. This leads to two different DSPSs between each pair of nodes made of either customers, depots or transfer points. Accordingly, the following notation will be used: •𝐴𝑇 𝑏(𝑖, 𝑗, 𝑑𝑡)is the arrival time at 𝑗when a black vehicle departs from 𝑖at time dt and follows the fastest path to reach 𝑗, as determined by the DSPS of nodes 𝑖and 𝑗for a black vehicle; •𝐷𝑇 𝑏(𝑖, 𝑗, 𝑎𝑡)is the inverse of 𝐴𝑇𝑏(𝑖, 𝑗, 𝑑𝑡)and is the departure time at 𝑖that allows a black vehicle to arrive at 𝑗at time at; •𝐴𝑇 𝑔(𝑖, 𝑗, 𝑑𝑡)is the arrival time at 𝑗when a green vehicle departs from 𝑖at time dt and follows the fastest path to reach 𝑗, as determined by the DSPS of the pair of nodes 𝑖and 𝑗for a green vehicle; •𝐷𝑇 𝑔(𝑖, 𝑗, 𝑎𝑡)is the inverse of 𝐴𝑇 𝑔(𝑖, 𝑗, 𝑑𝑡)and is the departure time at 𝑖that allows a green vehicle to arrive at 𝑗at time at. 5.3. Synchronization at a transfer point Let us consider 𝑡𝑝𝑘𝑏 ⋅a copy of transfer point tp ∈TP in the route of black vehicle 𝑘𝑏∈𝐾𝐵, where the subscript ⋅corresponds to a particular EURO Journal on Transportation and Logistics 13 (2024) 100143 6 F.O.G. Reyes et al. Fig. 3. (a) Travel speed function of arc (𝑖,𝑗) (b) Corresponding travel time function assuming that arc (𝑖,𝑗) is of length 4. Fig. 4. (a) Three different fastest paths between two nodes obtained at different time points using a time-dependent Dijkstra’s algorithm (b) Corresponding dominant shortest path structure. copy (visit) index of transfer point tp in the route of vehicle 𝑘𝑏. Let us also consider 𝑡𝑝𝑘𝑔 1 ⋅,𝑡𝑝𝑘𝑔 2 ⋅, ..., 𝑡𝑝𝑘𝑔 ℎ ⋅,ℎcopies of transfer point tp in the routes of green vehicles 𝑘𝑔 1,𝑘𝑔 2, ..., 𝑘𝑔 ℎ∈𝐾𝐺. We assume that black vehicle 𝑘𝑏needs to be synchronized with green vehicles 𝑘𝑔 1,𝑘𝑔 2, ..., 𝑘𝑔 ℎ at these copies of transfer point tp. We first account for the arrival time of the last green vehicle: 𝑎𝑡𝑚𝑎𝑥 = max 𝑙=1,…,ℎ{𝑎𝑡 𝑡𝑝𝑘𝑔 𝑙 ⋅ }(1) Then, the departure time of black vehicle 𝑘𝑏from 𝑡𝑝𝑘𝑏 ⋅is : 𝑑𝑡𝑡𝑝𝑘𝑏 ⋅ = max{𝑎𝑡𝑡𝑝𝑘𝑏 ⋅ , 𝑎𝑡𝑚𝑎𝑥}(2) That is, if all green vehicles arrive at the transfer point before black vehicle 𝑘𝑏, then the latter can depart immediately (given that the time to transfer loads from the black vehicle to green vehicles is null) with 𝑑𝑡𝑡𝑝𝑘𝑏 ⋅ =𝑎𝑡𝑡𝑝𝑘𝑏 ⋅ . Otherwise, vehicle 𝑘𝑏will depart at the arrival time 𝑎𝑡𝑚𝑎𝑥 of the last green vehicle. The departure of each green vehicle 𝑘𝑔 𝑙from 𝑡𝑝𝑘𝑔 𝑙 ⋅,𝑙=1, ..., ℎ, is: 𝑑𝑡 𝑡𝑝𝑘𝑔 𝑙 ⋅ = max{𝑎𝑡𝑡𝑝𝑘𝑏 ⋅ , 𝑎𝑡 𝑡𝑝𝑘𝑔 𝑙 ⋅ }, 𝑙 = 1,…, ℎ (3) That is, if green vehicle 𝑘𝑔 𝑙arrives at the transfer point before black vehicle 𝑘𝑏, it must wait for the arrival of vehicle 𝑘𝑏before it can depart from 𝑡𝑝𝑘𝑔 𝑙 ⋅. Otherwise, it can depart immediately with 𝑑𝑡 𝑡𝑝𝑘𝑔 𝑙 ⋅ =𝑎𝑡 𝑡𝑝𝑘𝑔 𝑙 ⋅ . 5.4. Time bounds In the following, we define earliest and latest time bounds for the arrival at and departure from each node in the route of a black or green vehicle, where a node can be a customer, a transfer point or a depot. That is, a vehicle must arrive at (depart from) a node before its latest arrival (departure) time to guarantee that the rest of the route satisfies the time constraints. For simplifications purposes, the forward and backward propagation procedures described below focus on a single black or green route and do not account for possible complex interactions among routes (see Section 5.5) . 5.4.1. Green route Here, we explain how to propagate the earliest and latest arrival and departure times in a green route. For this purpose, let us consider the route of green vehicle 𝑘𝑔∈𝐾𝐺which is made of (1) a copy 𝑑𝑔 0of the green depot to start the route, (2) a sequence of copies of one or more transfer points 𝑡𝑝𝑘𝑔 𝑙⋅,𝑙=1, ..., 𝑝, each followed by a green (or neutral) customer 𝑖𝑙∈𝐶𝐺∪𝐶𝐸,𝑙=1, ..., 𝑝and (3) a copy 𝑑𝑔 𝑝+1 of the green depot to end the route. That is, the route of green vehicle 𝑘𝑔is 𝑑𝑔 0,𝑡𝑝𝑘𝑔 1⋅,𝑖1,𝑡𝑝𝑘𝑔 2⋅,𝑖2, ..., 𝑡𝑝𝑘𝑔 𝑝⋅,𝑖𝑝,𝑑𝑔 𝑝+1. Earliest arrival and departure times First, the earliest departure time from 𝑑𝑔 0is set equal to 0. Then, we go forward by first computing the earliest arrival time eat at transfer point 𝑡𝑝𝑘𝑔 1⋅, using function 𝐴𝑇 𝑔: 𝑒𝑎𝑡𝑡𝑝𝑘𝑔 1⋅ =𝐴𝑇 𝑔(𝑑𝑔 0, 𝑡𝑝𝑘𝑔 1⋅,0) (4) Now, to determine the earliest departure time edt, we need to account for the corresponding black vehicle 𝑘𝑏from which green vehicle 𝑘𝑔should receive a load. Accordingly, if 𝑒𝑎𝑡𝑡𝑝𝑘𝑔 1⋅ < 𝑒𝑎𝑡𝑡𝑝𝑘𝑏 1⋅ , then 𝑒𝑑𝑡𝑡𝑝𝑘𝑔 1⋅ =𝑒𝑎𝑡𝑡𝑝𝑘𝑏 1⋅ , since green vehicle 𝑘𝑔cannot depart from the transfer point before the earliest arrival time of black vehicle 𝑘𝑏. Otherwise, 𝑒𝑑𝑡𝑡𝑝𝑘𝑔 1⋅ =𝑒𝑎𝑡𝑡𝑝𝑘𝑔 1⋅ , given that the time to transfer a load is null. EURO Journal on Transportation and Logistics 13 (2024) 100143 7 F.O.G. Reyes et al. Still going forward, we now consider customer 𝑖𝑖and compute its earliest departure time as : 𝑒𝑎𝑡𝑖1=𝐴𝑇 𝑔(𝑡𝑝𝑘𝑔 1⋅, 𝑖1, 𝑒𝑑𝑡𝑡𝑝𝑘𝑔 1⋅ )(5) Now, if 𝑒𝑎𝑡𝑖1< 𝛼𝑖1, then the earliest departure time 𝑒𝑑𝑡𝑖1=𝛼𝑖1+𝑠𝑡𝑖1, otherwise 𝑒𝑑𝑡𝑖1=𝑒𝑎𝑡𝑖1+𝑠𝑡𝑖1. This forward procedure is repeated until the end depot 𝑑𝑔 𝑝+1 is reached and its earliest arrival time is determined. Latest arrival and departure times We start by setting the latest arrival time lat at the end depot 𝑑𝑔 𝑝+1 to be the end of time horizon 𝑇, that is 𝑙𝑎𝑡𝑑𝑔 𝑝+1 =𝑇. Then, we go backward by first computing the latest departure time ldt at customer 𝑖𝑝that allows vehicle 𝑘𝑔to arrive at 𝑑𝑔 𝑝+1 at time 𝑙𝑎𝑡𝑑𝑔 𝑝+1 , using function 𝐷𝑇 𝑔: 𝑙𝑑𝑡𝑖𝑝=𝐷𝑇 𝑔(𝑖𝑝, 𝑑𝑔 𝑝+1, 𝑙𝑎𝑡𝑑𝑔 𝑝+1 )(6) Now, if 𝑙𝑑𝑡𝑖𝑝> 𝛽𝑖𝑝+𝑠𝑡𝑖𝑝, then 𝑙𝑑𝑡𝑖𝑝is reset to 𝛽𝑖𝑝+𝑠𝑡𝑖𝑝, because vehicle 𝑘𝑔cannot depart from 𝑖𝑝later than 𝛽𝑖𝑝+𝑠𝑡𝑖𝑝without violating the time window constraint (i.e., the arrival time cannot exceed 𝛽𝑖𝑝). Then, the latest arrival time at customer 𝑖𝑝is simply computed as 𝑙𝑎𝑡𝑖𝑝=𝑙𝑑𝑡𝑖𝑝−𝑠𝑡𝑖𝑝. Still going backward, we now consider the transfer point 𝑡𝑝𝑘𝑔 𝑝⋅and compute its latest departure time as : 𝑙𝑑𝑡𝑡𝑝𝑘𝑔 𝑝⋅ =𝐷𝑇 𝑔(𝑡𝑝𝑘𝑔 𝑝⋅, 𝑖𝑝, 𝑙𝑎𝑡𝑖𝑝)(7) To determine the latest arrival time of the green vehicle 𝑘𝑔, we must account for the corresponding black vehicle 𝑘𝑏that transfers a load to vehicle 𝑘𝑔. That is, the green vehicle cannot arrive after the latest departure time of black vehicle 𝑘𝑏through the following formula : 𝑙𝑎𝑡𝑡𝑝𝑘𝑔 𝑝⋅ = min{𝑙𝑑𝑡𝑡𝑝𝑘𝑔 𝑝⋅ , 𝑙𝑑𝑡𝑡𝑝𝑘𝑏 𝑝⋅ }(8) This backward procedure is applied until the starting depot 𝑑𝑔 0is reached and its latest departure time is determined. It should be noted that a forward propagation starting from the latest departure time at 𝑑𝑔 0, until 𝑑𝑔 𝑝+1 is reached, would produce the latest feasible schedule (i.e., latest possible arrival and departure times at each node along the green route). 5.4.2. Black route Here, we explain how to propagate the earliest and latest arrival and departure times in a black route. For this purpose, let us consider the route of black vehicle 𝑘𝑏∈𝐾𝐵which is made of (1) a copy 𝑑𝑏 0of the black depot to start the route, (2) an arbitrary sequence of length 𝑝of black (or neutral) customers and copies of one or more transfer points and (3) a copy 𝑑𝑏 𝑝+1 of the black depot to end the route. Earliest arrival and departure times The procedure to compute the earliest arrival and departure times in a black route is similar to the one described for the green route, but two differences are noteworthy: (1) the function 𝐴𝑇 𝑏is used to compute the arrival time at a given node from the earliest departure time of the previous node and (2) the earliest departure time at a copy of a transfer point is computed differently, because green vehicles that visit the same transfer point to get a load from the black vehicle must be accounted for. Considering case (2), let us suppose that black vehicle 𝑘𝑏visits copy 𝑡𝑝𝑘𝑏 ⋅of transfer point tp ∈TP and that ℎgreen vehicles 𝑘𝑔 1, 𝑘𝑔 2,…, 𝑘𝑔 ℎ visit copies 𝑡𝑝𝑘𝑔 1 ⋅, 𝑡𝑝𝑘𝑔 2 ⋅,…, 𝑡𝑝𝑘𝑔 ℎ ⋅of the same transfer point and that synchronization is required (i.e., black vehicle 𝑘𝑏must transfer a load to each green vehicle). To compute the earliest departure time of vehicle 𝑘𝑏at the transfer point, we first consider the maximum earliest arrival time over all green vehicles, that is: 𝑒𝑎𝑡𝑚𝑎𝑥 = max 𝑙=1,…,ℎ{𝑒𝑎𝑡 𝑡𝑝𝑘𝑔 𝑙 ⋅ }(9) Then, the earliest departure time of vehicle 𝑘𝑏at 𝑡𝑝𝑘𝑏 ⋅can be computed from its earliest arrival time as follow: 𝑒𝑑𝑡𝑡𝑝𝑘𝑏 ⋅ = max{𝑒𝑎𝑡𝑚𝑎𝑥, 𝑒𝑎𝑡𝑡𝑝𝑘𝑏 ⋅ }(10) That is, black vehicle 𝑘𝑏cannot depart earlier than the earliest arrival time of the last green vehicle, otherwise one or more green vehicles will not get their load. Latest arrival and departure times The procedure to compute the latest arrival and departure times of a black vehicle is similar to the one described for a green vehicle, although two differences are noteworthy: (1) the function 𝐷𝑇 𝑏is used to compute the departure time from a given node to reach the next node at its latest arrival time and (2) the latest arrival time at a copy of a transfer point is computed differently, because green vehicles that visit the same transfer point to get a load from the black vehicle must be accounted for. Considering case (2), let us suppose that black vehicle 𝑘𝑏visits copy 𝑡𝑝𝑘𝑏 ⋅of transfer point tp ∈TP and that ℎgreen vehicles 𝑘𝑔 1, 𝑘𝑔 2,…, 𝑘𝑔 ℎ visit copies 𝑡𝑝𝑘𝑔 1 ⋅, 𝑡𝑝𝑘𝑔 2 ⋅,…, 𝑡𝑝𝑘𝑔 ℎ ⋅of the same transfer point and that synchronization is required (i.e., black vehicle 𝑘𝑏must transfer a load to each green vehicle). To compute the latest arrival time of vehicle 𝑘𝑏at the transfer point, we first consider the minimum latest departure time over all green vehicles, that is: 𝑙𝑑𝑡𝑚𝑖𝑛 = min 𝑙=1,…,ℎ{𝑙𝑑𝑡 𝑡𝑝𝑘𝑔 𝑙 ⋅ }(11) Then, the latest arrival time of vehicle 𝑘𝑏at 𝑡𝑝𝑘𝑏 ⋅can be computed from its latest departure time as follow: 𝑙𝑎𝑡𝑡𝑝𝑘𝑏 ⋅ = min{𝑙𝑑𝑡𝑚𝑖𝑛, 𝑙𝑑𝑡𝑡𝑝𝑘𝑏 ⋅ }(12) That is, black vehicle 𝑘𝑏cannot arrive at the transfer point later than the minimum latest departure time over all green vehicles that require synchronization, otherwise one or more green vehicles will not get their load. 5.5. Interaction among routes In the previous section, our description of forward and backward propagation procedures to derive time bounds has focused on a single black or green route. However, complex interactions may occur when multiple black and green routes are involved. Fig. 5 shows an example where customer 𝑖is inserted between nodes prev and next in black route 𝑘𝑏 2(as it occurs during the recreate procedure of SISR). Forward propagation is illustrated in Fig. 5(a). First, the earliest arrival and departure times of the newly inserted customer 𝑖are calculated from the earliest departure time at transfer point prev. Then, forward propagation is triggered along the black route. However, a green route that connects the three black routes is encountered at the transfer point just after customer next. Thus, another forward propagation is triggered at this transfer point along the green route which, in turn, leads to a transfer point that connects the green route to black route 𝑘𝑏 3, thus triggering another forward propagation along that black route. It should be noted that this illustration is a worst case, because forward propagation along a route terminates as soon as the earliest departure time from a node does not change. Backward propagation is illustrated in Fig. 5(b). First, latest arrival and departure times at the newly inserted customer 𝑖are calculated from the latest arrival time at customer next. Then, backward propagation is triggered along black route 𝑘𝑏 2. Since node prev is a transfer point, it triggers another backward propagation along the corresponding green route. Still going backward along the black route, another transfer point is met that triggers backward propagation along another green route. Finally, both green routes connect to black route 𝑘𝑏 1at the same transfer point, thus triggering a backward propagation along that black route also. Once again, this illustration is a worst case, because EURO Journal on Transportation and Logistics 13 (2024) 100143 8 F.O.G. Reyes et al. Fig. 13. Example of an instance generated with road network 𝑅𝑁3. (For colors of the nodes, the reader is referred to the web version of this article). only. For this reason, black customers are located in the outside and boundary regions, while green customers are located in the downtown and boundary regions. Neutral customers and transfer points are only found in the boundary region, since they can be visited by both black and green vehicles. •The black depot is randomly located in the outside region, while the green depot is randomly located in the downtown region. There are also 10 transfer points that are randomly located in the boundary region. •Arcs are characterized by the regions where they are found. Let (𝑖, 𝑗)be an arc in the network. If both 𝑖and 𝑗are in the boundary (or frontier) region, then the arc is of type F (and is accessible to both black and green vehicles). If both 𝑖and 𝑗are in the downtown region, or one is in the downtown region and the other in the boundary region, then the arc is of type D (and is accessible only to green vehicles). If both 𝑖and 𝑗are in the outside region, or if one is in the outside region and the other is in the boundary region, then the arc is of type O (and is accessible only to black vehicles). It should be noted that no arcs connect the downtown and outside regions. Given that the regions are the same for a given network, then the arc type also stays the same for all instances generated from a given road network. •The test instances are duplicated by considering two different sets of travel speed multipliers (scenarios), where a speed multiplier depends on the time period, vehicle (black or green) and arc type (D, F, O), see Tables 1 and 2. In the two scenarios 𝑆𝐼and 𝑆𝐼𝐼 , the second and fourth periods correspond to rush hours. The speed multipliers of green vehicles are fixed at 1 everywhere, which means that they are not affected by congestion since they are small (e.g., bicycles). Black vehicles are faster than green vehicles when there is no congestion. However, they are slower than green vehicles in the boundary region during rush hours in scenario 𝑆𝐼 while they have the same speed than green vehicles in scenario 𝑆𝐼𝐼 . The second scenario is aimed at evaluating the impact of increasing the speed of black vehicles when compared to green vehicles. •The capacity of black vehicles is set to 40. •The demand of each customer is randomly selected from {1,…, 5}. Overall, there are 4 road networks ×3 numbers of customers (𝑛𝑐)× 4 customer distributions ×2 types of time windows ×2 scenarios for a total of 192 instances, that is, 96 instances for each scenario. Table 1 Speed multipliers for (a) black vehicles and (b) green vehicles under scenario 𝑆𝐼. (a) Black vehicles Arc type Time period 𝜏1= [0,20) 𝜏2= [20,30) 𝜏3= [30,70) 𝜏4= [70,80) 𝜏5= [80,100) 𝑇 𝑦𝑝𝑒 F 1.2 0.8 1.2 0.8 1.2 𝑇 𝑦𝑝𝑒 O 1.5 1.0 1.5 1.0 1.5 (b) Green vehicles Arc type Time period 𝜏1= [0,20) 𝜏2= [20,30) 𝜏3= [30,70) 𝜏4= [70,80) 𝜏5= [80,100) 𝑇 𝑦𝑝𝑒 D 1.0 1.0 1.0 1.0 1.0 𝑇 𝑦𝑝𝑒 F 1.0 1.0 1.0 1.0 1.0 Table 2 Speed multipliers for (a) black vehicles and (b) green vehicles under scenario 𝑆𝐼𝐼 . (a) Black vehicles Arc type Time period 𝜏1= [0,20) 𝜏2= [20,30) 𝜏3= [30,70) 𝜏4= [70,80) 𝜏5= [80,100) 𝑇 𝑦𝑝𝑒 F 1.5 1.0 1.5 1.0 1.5 𝑇 𝑦𝑝𝑒 O 2.0 1.5 2.0 1.5 2.0 (b) Green vehicles Arc type Time period 𝜏1= [0,20) 𝜏2= [20,30) 𝜏3= [30,70) 𝜏4= [70,80) 𝜏5= [80,100) 𝑇 𝑦𝑝𝑒 D 1.0 1.0 1.0 1.0 1.0 𝑇 𝑦𝑝𝑒 F 1.0 1.0 1.0 1.0 1.0 7.2. Parameter tuning Four parameters have a significant impact on the performance of our SISR, namely, 𝑐(average number of removed customers), 𝑛𝑝𝑜𝑠 (number of best insertion positions, based on the approximation), 𝑛𝑓𝑡𝑝 (number of nearest feasible transfer points) and 𝐿𝑚𝑎𝑥 (maximum length of removed strings). To adjust their values, we selected a subset of 16 tuning instances with 100 customers by randomly selecting only one of the four networks, for each possible configuration of customer distribution (𝐷1,𝐷2,𝐷3,𝐷4), time window (NTW, WTW) and scenario (𝑆1,𝑆2). Given that solution quality tends to improve with increasing values of 𝑛𝑝𝑜𝑠 and 𝑛𝑓𝑡𝑝, at the expense of computation time, these two parameters were first set to high values, that is, 𝑛𝑝𝑜𝑠 =7 and 𝑛𝑓𝑡𝑝 = 10 (the latter value cannot be larger, since there are only 10 transfer points in each instance). In other words, we did not care at this point about computation time. Then, we focused on parameters 𝑐and 𝐿𝑚𝑎𝑥 EURO Journal on Transportation and Logistics 13 (2024) 100143 15 F.O.G. Reyes et al. Fig. 14. Convergence curves observed for instances with 200 customers generated with road network 𝑅𝑁4under scenario 𝑆𝐼. (For colors of the curves, the reader is referred to the web version of this article). Table 3 Impact of parameter values on solution quality and computation times. Parameter values 𝑐=9 11 13 15 17 19 Avg. cost 2552.1 2532.7 2523.3 2520.3 2520.1 2520.0 Avg. time 1.29 1.54 1.81 2.10 2.36 2.56 𝑛𝑝𝑜𝑠 =12345678910 Avg. cost 2537.9 2524.6 2521.4 2521.9 2520.3 2521.0 2520.3 2520.2 2519.6 2520.6 Avg. time 1.76 1.92 1.93 1.99 2.00 2.05 2.1 2.16 2.16 2.21 𝑛𝑓𝑡𝑝 =12345678910 Avg. cost 2531.2 2525.5 2522.0 2520.2 2520.6 2518.9 2520.1 2520.3 2520.4 2520.3 Avg. time 0.83 1.15 1.41 1.61 1.75 1.88 1.97 2.02 2.06 2.10 𝐿𝑚𝑎𝑥 =12345678910 Avg. cost 2542.9 2522.1 2520.2 2519.7 2520.4 2520.0 2520.1 2520.3 2519.8 2520.5 Avg. time 2.02 2.14 2.09 2.09 2.08 2.10 2.11 2.10 2.08 2.08 and tuned them with the IRACE software (López-Ibánez et al.,2016), using 𝑐= {5,…,17} and 𝐿𝑚𝑎𝑥 = {3,…,13}. The best values returned by IRACE were 𝑐= 15 and 𝐿𝑚𝑎𝑥 = 8. Based on the default configuration 𝑐= 15,𝑛𝑝𝑜𝑠 =7, 𝑛𝑓𝑡𝑝 =10 and 𝐿𝑚𝑎𝑥 =8, we then modified the value of one parameter at a time, keeping the other parameters at their default value. The values considered for each parameter were: 𝑐= {9,11,13,15,17,19},𝑛𝑝𝑜𝑠 = {1,…,10},𝑛𝑓𝑡𝑝 = {1,…,10} and 𝐿𝑚𝑎𝑥 = {1,…,10}. Since our algorithm is non deterministic, we show the average results (solution quality, computation time in hours) obtained over 10 runs on each tuning instance in Table 3. As expected, increasing the values of parameters 𝑛𝑝𝑜𝑠 and 𝑛𝑓𝑡𝑝 leads to an increase in computation time, although the impact is more significant in the case of 𝑛𝑓𝑡𝑝, since it increases the number of possible insertions of green and neutral customers in Phase III and Phase IV of the Recreate method (these insertions are quite complicated). The computation times increase even more with increasing values of parameter 𝑐because more removed customers simply mean more customers to be reinserted. On the other hand, parameter 𝐿𝑚𝑎𝑥 has no impact on computation time. With regard to solution quality, we observe an improvement in solution quality for the first values of each parameter, but then some kind of stagnation is observed. Accordingly, the parameter setting 𝑐=15, 𝑛𝑝𝑜𝑠 = 3,𝑛𝑓𝑡𝑝 = 4 and 𝐿𝑚𝑎𝑥 = 4 was chosen for the experiments reported in the following sections. We also checked that this particular combination of parameter values led to good solutions on the tuning instances, which turned to be true with an average solution cost of 2522.2 and average computation time of 1.47 h. Some experiments were also performed with regard to the number of iterations. We observed that convergence is obtained, even on the largest instances with 200 customers, after a maximum of 300,000 iterations. That is, a plateau is reached and no further significant improvement in solution quality is observed. Fig. 14 shows an example of convergence curves for the best solutions found on instances with 200 customers generated with road network 𝑅𝑁4, using the four customer distributions, and both narrow and wide time windows, under scenario 𝑆𝐼. In this figure, black, green, gray and blue curves are associated with customer distributions 𝐷1,𝐷2,𝐷3and 𝐷4, respectively, while full lines and broken lines are associated with instances with narrow and wide time windows, respectively. Based on the results obtained, the number of iterations was set to 300,000 for all instances (which is admittedly too much for instances with 50 and 100 customers). 7.3. Results on test instances We report in this section the results produced by our algorithm on the whole set of test instances, based on 10 different runs on each instance. Table 8 in the Appendix reports the best and average costs, as well as the average computation times in hours for each instance. Each line of this table corresponds to a particular type of instance using the notation 𝑅𝑁𝑥_𝑛𝑦_𝐷𝑧, where 𝑥is the road network index, 𝑦is the number of customers and 𝑧is the customer distribution index. For each type, we show the results obtained on instances associated with narrow time windows (NTW) and wide time windows (WTW) under scenarios 𝑆𝐼and 𝑆𝐼𝐼 .Table 4 in this section is a reduced version, where averages are taken over the four road networks. That is, 𝑛𝑦_𝐷𝑧 in Table 4 encompasses 𝑅𝑁1_𝑛𝑦_𝐷𝑧,𝑅𝑁2_𝑛𝑦_𝐷𝑧,𝑅𝑁3_𝑛𝑦_𝐷𝑧 and 𝑅𝑁4_𝑛𝑦_𝐷𝑧, so that the numbers in Table 4 correspond to the bold Avg. lines in the full table in the Appendix. Since there are no similar instances in the literature that we could refer to for comparison purposes, Table 5 reports both the average best EURO Journal on Transportation and Logistics 13 (2024) 100143 16 F.O.G. Reyes et al. Table 4 Solution cost and computation time in hours for each subset of instances. Instances NTW WTW Avg. best cost Avg. cost Avg. time Avg. best cost Avg. cost Avg. time 𝑆𝐼𝑆𝐼𝐼 𝑆𝐼𝑆𝐼𝐼 𝑆𝐼𝑆𝐼𝐼 𝑆𝐼𝑆𝐼𝐼 𝑆𝐼𝑆𝐼𝐼 𝑆𝐼𝑆𝐼𝐼 n50_D1 1397.8 1125.3 1398.9 1126.6 0.27 0.28 1219.0 929.5 1220.5 930.0 0.26 0.26 n50_D2 2002.5 1841.0 2004.5 1844.5 0.77 0.87 1801.4 1608.4 1805.2 1615.5 0.76 0.81 n50_D3 1368.4 1112.9 1374.8 1118.0 0.66 0.69 1117.0 899.2 1117.0 901.8 0.59 0.54 n50_D4 1748.9 1517.1 1752.0 1520.4 0.53 0.56 1542.4 1279.1 1543.8 1280.3 0.51 0.51 n100_D1 2498.1 2030.9 2505.7 2034.1 0.64 0.70 2128.9 1702.6 2140.4 1710.0 0.64 0.65 n100_D2 3684.7 3298.5 3696.0 3308.6 2.18 2.37 3314.9 2909.0 3326.5 2922.0 2.06 2.13 n100_D3 2368.0 1976.9 2372.9 1979.9 1.86 1.92 1916.3 1564.2 1921.6 1572.9 1.42 1.39 n100_D4 3138.3 2734.4 3148.2 2756.3 1.34 1.55 2743.1 2303.3 2754.8 2316.5 1.35 1.34 n200_D1 4186.9 3461.1 4221.5 3489.7 1.71 1.88 3582.6 2822.7 3610.5 2855.8 1.42 1.57 n200_D2 6980.9 6356.7 7028.4 6408.2 5.76 6.47 6321.5 5691.5 6380.1 5758.7 5.47 5.69 n200_D3 4399.9 3638.8 4439.4 3672.7 4.54 5.03 3636.4 2918.6 3663.4 2941.5 3.44 3.52 n200_D4 5613.9 4894.5 5642.2 4930.7 3.70 4.21 4927.8 4255.0 4967.9 4294.4 3.41 3.57 Overall Avg. 3282.4 2832.3 3298.7 2849.1 2.00 2.21 2854.3 2407.0 2871.0 2425.0 1.78 1.83 Table 5 Average best improvements and average improvements for each subset of instances. Instances NTW WTW Avg. best Impr. Avg. Impr. Avg. best Impr. Avg. Impr. 𝑆𝐼𝑆𝐼𝐼 𝑆𝐼𝑆𝐼𝐼 𝑆𝐼𝑆𝐼𝐼 𝑆𝐼𝑆𝐼𝐼 n50_D1 13.02 14.14 10.48 11.19 16.03 18.38 13.48 15.03 n50_D2 13.28 15.33 11.21 13.15 16.75 19.36 14.28 16.51 n50_D3 33.53 41.27 30.95 38.51 35.42 43.14 32.02 39.05 n50_D4 11.20 14.87 9.49 11.64 15.06 18.21 12.62 15.52 n100_D1 17.59 20.28 15.96 18.46 22.32 24.56 20.33 22.70 n100_D2 17.55 19.68 15.34 17.67 19.48 22.70 17.63 20.68 n100_D3 41.76 46.49 39.68 44.93 44.31 48.55 42.14 45.87 n100_D4 14.67 17.62 13.27 15.28 20.16 22.18 17.98 20.56 n200_D1 25.57 27.73 24.00 25.46 29.57 31.24 27.73 29.54 n200_D2 18.76 22.39 17.32 21.05 21.01 23.97 19.50 22.20 n200_D3 46.53 52.75 45.27 51.05 48.17 54.24 46.51 52.07 n200_D4 19.75 22.61 18.46 20.89 23.95 25.71 21.95 23.81 Overall Avg. 22.77 26.26 20.95 24.11 26.02 29.35 23.85 26.96 improvements and average improvements provided by our algorithm over the initial solutions, to measure its optimization power. Denoting 𝑠𝑖and 𝑠𝑓the initial and final solutions produced by our algorithm on a given instance, with 𝑐𝑜𝑠𝑡(𝑠𝑖)and 𝑐𝑜𝑠𝑡(𝑠𝑓)their respective cost, the percentage of improvement of the final solution over the initial one is calculated as follow: 𝐼𝑚𝑝𝑟 = 100 (𝑐𝑜𝑠𝑡(𝑠𝑖) − 𝑐𝑜𝑠𝑡(𝑠𝑓) 𝑐𝑜𝑠𝑡(𝑠𝑖))(14) Tables 4 and 5will be referred to in the following sections when we analyze in more detail the behavior of our algorithm for different characteristics of the test instances. For now, we can observe the obvious increase in solution cost and computation time with instance size in Table 4. We also note the overall average improvements over the initial solutions in Table 5 that stand between 20% and 30%, which is substantial. It is clear that the greedy insertion heuristic is limited with regard to solution quality, but still, these percentages of improvement show that our algorithm can take advantage of optimization opportunities. 7.4. Synchronization In this section, we examine if synchronization between black and green vehicles at transfer points is efficient. For this purpose, we define the percentage of solution cost (duration) that corresponds to the total time that black and green vehicles spend at transfer points. For a given solution 𝑠, this percentage is denoted as 𝜌𝑠. In Eq. (15), this percentage Table 6 Minimum, maximum and average values of 𝜌𝑠for each subset of instances. Instances Min Max Avg. n50_D1 0.12 2.18 0.77 n50_D2 1.05 3.71 2.08 n50_D3 0.03 2.55 1.02 n50_D4 0.18 3.23 1.60 Avg. 0.34 2.92 1.37 n100_D1 0.01 3.03 0.90 n100_D2 0.82 3.11 1.91 n100_D3 0.19 1.97 0.95 n100_D4 0.45 3.13 1.48 Avg. 0.37 2.81 1.31 n200_D1 0.45 2.45 1.25 n200_D2 0.84 3.37 1.91 n200_D3 0.25 2.20 0.86 n200_D4 0.53 2.85 1.63 Avg. 0.52 2.71 1.41 Overall Avg. 0.01 3.71 1.36 is defined using 𝛥TP𝑏 𝑠and 𝛥TP𝑔 𝑠, which are the total time spent at transfer points by black and green vehicles, respectively. 𝜌𝑠= 100 (𝛥TP𝑏 𝑠+𝛥TP𝑔 𝑠 𝑐𝑜𝑠𝑡(𝑠))(15) Table 6 reports the minimum, maximum and average values of 𝜌𝑠 on different subsets of instances. The format of this table is reduced when compared to Tables 4 and 5by also averaging over the two types of time windows and the two scenarios. We observe that the 𝜌𝑠 values are smaller for customer distributions 𝐷1and 𝐷3. In fact, if we compute the average 𝜌𝑠values for 𝐷1,𝐷2,𝐷3and 𝐷4, we obtain 0.97%, 1.97%, 0.94% and 1.57%, respectively. The fact that 𝐷1and 𝐷3lead to better synchronization than 𝐷2and 𝐷4can be explained by their small percentage of green customers (20%), since they are the only ones for which synchronization at a transfer point is mandatory. In any case, the differences observed are small in absolute terms. The fact is that only 1.36% (overall average) of solution cost is due to synchronization, which indicates that synchronization is well achieved. 7.5. Customer distributions Table 4 shows that customer distributions 𝐷1and 𝐷3are the best for solution cost, while 𝐷1outperforms the three other distributions for computation time. Both 𝐷1and 𝐷3have only 20% of green customers, which is beneficial for solution cost because these customers require a detour at a transfer point for both a black and a green vehicle. Furthermore, 𝐷1has a large percentage of 60% of black EURO Journal on Transportation and Logistics 13 (2024) 100143 17 F.O.G. Reyes et al. customers (versus only 20% of neutral customers), which is helpful for computation time because only simple insertions in black routes need to be considered for black customers. Although 𝐷3is competitive with 𝐷1for solution cost, this is not the case for computation time. Distribution 𝐷3has 60% of neutral customers (versus only 20% of black customers) and their potential insertion in both black and green routes need to be considered. As opposed to black routes, insertion in green routes is more complicated and more computationally expensive. With the largest percentage of 60% of green customers, distribution 𝐷2is consequently the worst for solution cost and computation time. When considering improvements over initial solutions in Table 5, the largest improvements are associated with distribution 𝐷3. This distribution has the largest percentage of neutral customers (60%) and these customers offer more flexibility for optimization because they can be inserted either in black or green routes. In this case, the average percentage of neutral customers belonging to black routes in the initial solutions ranges between 51.1% and 54.4%. However, in the final solutions, these percentages drastically increase between 94.6% and 98.6%. That is, the optimization algorithm finds ways to move a large proportion of neutral customers in black routes, which decreases solution cost. Accordingly, more neutral customers means more opportunities for improvement. 7.6. Time windows Tables 4 and 5show that better solution costs and larger improvements over initial solutions are associated with instances with wide time windows. Clearly, these time windows offer more feasible insertion places for customers and, consequently, greater flexibility for the optimization procedure to move them around. For example, the percentage of neutral customers in black routes for the instances with wide time windows is 98.3%, as compared with 92.8% for instances with narrow time windows. We also observed a reduced number of routes in the solutions obtained on instances with wide time windows, when compared with narrow time windows, namely 25.2% less black routes and 15.1% less green routes. Finally, black vehicles visit 2% more transfer points on average in the presence of wide time windows. 7.7. Scenarios Better solution costs are observed in Table 4 under scenario 𝑆𝐼𝐼 , when compared with scenario 𝑆𝐼. Since black vehicles now travel faster, the time windows at black and neutral customers become easier to satisfy. We also observed that, on average, a larger percentage of neutral customers is served by black vehicles in the final solutions under scenario 𝑆𝐼𝐼 (97.4%) when compared to scenario 𝑆𝐼(93.7%). It is generally less costly to serve neutral customers with black vehicles and the latter can get a larger share of neutral customers when their speed increases. In other words, they win more often the ‘‘battle’’ for neutral customers against green vehicles in the boundary region. We also observed a reduction of 14.7% and 6.2% in the number of black and green routes, respectively, under scenario 𝑆𝐼𝐼 , with corresponding increases of 20.5% and 3.5% in the average number of customers in black and green routes, respectively. Finally, black vehicles visit 4.4% more transfer points on average, when compared to scenario 𝑆𝐼. 7.8. Number of transfer points We propose here an experiment where we gradually reduce the number of transfer points to see the impact on the solutions obtained. To this end, some transfer points are randomly selected and transformed into simple nodes (road junctions). In some cases, infeasible instances may be created (e.g., it may not be possible for any green vehicle to visit a transfer point and serve a green customer before the upper bound of its time window). Table 7 Average gaps between average solution costs with 10 transfer points and 𝑘= 8, 6, 4 transfer points. Instances Avg. Gap 𝑘= 8 𝑘= 6 𝑘= 4 RN1_n100_D1 −0.04 6.02 – RN2_n100_D1 −0.10 −0.11 0.47 RN3_n100_D1 −0.02 −0.04 0.51 RN4_n100_D1 0.51 1.74 1.76 RN1_n100_D2 2.06 10.95 27.45 RN2_n100_D2 1.56 1.58 5.86 RN3_n100_D2 0.62 1.56 9.08 RN4_n100_D2 1.81 7.09 7.30 RN1_n100_D3 0.48 7.32 16.08 RN2_n100_D3 −0.05 −0.09 0.48 RN3_n100_D3 1.12 2.34 6.55 RN4_n100_D3 1.24 1.79 1.82 RN1_n100_D4 1.00 9.72 – RN2_n100_D4 −0.86 −0.75 −0.02 RN3_n100_D4 0.19 1.21 4.45 RN4_n100_D4 0.86 1.47 1.68 Here, we focus on the test instances with 100 customers. First, two transfer points are randomly chosen and removed from the 10 original ones to obtain instances with 8 transfer points. From these 8 transfer points, the procedure is repeated to get instances with 6 transfer points. Finally, two additional transfer points are removed to get instances with 4 transfer points. Ten runs were performed on each feasible instance with a reduced number of transfer points. Table 7 shows the average gaps between the average cost of solutions obtained with 10 transfer points and with 𝑘= 8,6,4transfer points, based on Eq. (16). Each line in this table is the average over four instances, considering that there are two types of time windows and two scenarios. 𝐺𝑎𝑝 = 100 (𝑎𝑣𝑒|𝑇 𝑃 |=𝑘 𝑐𝑜𝑠𝑡 −𝑎𝑣𝑒|𝑇 𝑃 |=10 𝑐𝑜𝑠𝑡 𝑎𝑣𝑒|𝑇 𝑃 |=10 𝑐𝑜𝑠𝑡 )(16) We observed no infeasible instance with 𝑘=8 transfer points, one infeasible instance with 𝑘=6 transfer points and 19 infeasible instances (out of 64 instances) with 𝑘=4. Thus, some averages are computed with less than four gaps. When there is no value, the four instances are infeasible. A few small negative values appear in the table, which means that slightly better solutions are obtained with fewer transfer points. This situation can occur, due to the randomized nature of our algorithm, when the removed transfer points are not or seldom used in the original instances, thus producing no or little impact on solution quality. Obviously, removing transfer points generally lead to worse solutions and this trend is more pronounced when more transfer points are removed. Customer distribution 𝐷2is more affected than the other distributions due to its large percentage of green customers (which must use transfer points). We also see that instances associated with road network 𝑅𝑁1are greatly affected when 4 or 6 transfer points are removed. In the solutions obtained with 10 transfer points, we observed that some transfer points are much more exploited than others. Clearly, such critical transfer points are more likely to disappear when more transfer points are removed, which in turn greatly impact solution quality. 8. Conclusion In this work, we have proposed a new SISR metaheuristic for the TDVRPTWTPRN. To the best of our knowledge, this challenging problem where customers can be served either directly by black vehicles or indirectly by green vehicles through transfer points has not been EURO Journal on Transportation and Logistics 13 (2024) 100143 18 F.O.G. Reyes et al. Table 8 Solution cost and computation time in hours for each instance. Instances NTW WTW Best cost Avg. cost Avg. time Best cost Avg. cost Avg. time (10 runs) (10 runs) (h) (10 runs) (10 runs) (h) 𝑆𝐼𝑆𝐼𝐼 𝑆𝐼𝑆𝐼𝐼 𝑆𝐼𝑆𝐼𝐼 𝑆𝐼𝑆𝐼𝐼 𝑆𝐼𝑆𝐼𝐼 𝑆𝐼𝑆𝐼𝐼 RN1_n50_D1 1563.47 1217.90 1563.91 1217.90 0.30 0.29 1331.90 1030.23 1332.03 1030.23 0.26 0.29 RN2_n50_D1 1365.48 1096.91 1365.48 1096.91 0.25 0.23 1159.71 898.09 1165.04 898.09 0.26 0.24 RN3_n50_D1 1430.09 1159.05 1432.17 1159.65 0.25 0.27 1292.61 950.84 1292.81 950.85 0.25 0.24 RN4_n50_D1 1232.26 1027.51 1233.85 1031.95 0.28 0.34 1091.71 839.01 1092.13 840.98 0.27 0.28 Avg. 1397.83 1125.34 1398.85 1126.60 0.27 0.28 1218.98 929.54 1220.50 930.04 0.26 0.26 RN1_n50_D2 2295.75 2067.62 2295.83 2069.44 0.81 0.90 2018.20 1737.63 2018.72 1742.08 0.82 0.89 RN2_n50_D2 1929.28 1752.89 1932.68 1753.13 0.74 0.81 1664.44 1476.99 1664.84 1478.37 0.71 0.77 RN3_n50_D2 1996.08 1864.24 1996.14 1875.20 0.69 0.71 1843.62 1696.06 1855.69 1713.47 0.63 0.71 RN4_n50_D2 1788.88 1679.38 1793.19 1680.13 0.83 1.06 1679.18 1523.08 1681.66 1527.97 0.88 0.88 Avg. 2002.50 1841.03 2004.46 1844.48 0.77 0.87 1801.36 1608.44 1805.23 1615.47 0.76 0.81 RN1_n50_D3 1440.17 1183.37 1449.07 1184.69 0.65 0.68 1234.84 983.12 1234.86 985.14 0.61 0.53 RN2_n50_D3 1365.98 1087.79 1378.58 1105.24 0.72 0.66 1059.19 845.14 1059.19 853.41 0.57 0.53 RN3_n50_D3 1353.17 1120.44 1356.42 1121.83 0.57 0.61 1094.61 931.62 1094.61 931.62 0.59 0.53 RN4_n50_D3 1314.31 1059.94 1314.97 1060.22 0.70 0.78 1079.52 836.87 1079.52 837.14 0.60 0.58 Avg. 1368.41 1112.88 1374.76 1118.00 0.66 0.69 1117.04 899.19 1117.05 901.83 0.59 0.54 RN1_n50_D4 1786.66 1587.00 1786.66 1587.00 0.49 0.53 1679.56 1365.57 1679.56 1365.71 0.56 0.47 RN2_n50_D4 1944.47 1629.15 1947.70 1630.98 0.54 0.51 1622.76 1328.29 1622.76 1330.34 0.57 0.61 RN3_n50_D4 1535.58 1376.85 1537.07 1386.52 0.48 0.47 1355.06 1195.54 1355.06 1197.90 0.41 0.45 RN4_n50_D4 1728.69 1475.37 1736.48 1477.26 0.61 0.72 1512.20 1227.06 1517.89 1227.06 0.49 0.51 Avg. 1748.85 1517.09 1751.98 1520.44 0.53 0.56 1542.39 1279.12 1543.82 1280.25 0.51 0.51 RN1_n100_D1 2526.91 2085.89 2528.87 2088.53 0.64 0.70 2189.86 1662.22 2198.75 1675.03 0.65 0.56 RN2_n100_D1 2559.89 2067.83 2570.32 2071.26 0.53 0.65 2161.96 1726.10 2177.80 1727.87 0.56 0.59 RN3_n100_D1 2540.37 2113.93 2540.48 2115.17 0.64 0.65 2200.53 1889.07 2202.66 1891.56 0.70 0.78 RN4_n100_D1 2365.21 1856.08 2383.25 1861.46 0.73 0.81 1963.41 1533.17 1982.42 1545.58 0.63 0.69 Avg. 2498.09 2030.93 2505.73 2034.11 0.64 0.70 2128.94 1702.64 2140.41 1710.01 0.64 0.65 RN1_n100_D2 3804.11 3452.62 3812.33 3461.66 2.07 2.22 3426.28 2988.77 3434.23 3006.73 2.03 1.99 RN2_n100_D2 3537.95 3237.98 3543.50 3246.02 1.84 2.37 3256.69 2866.68 3265.31 2873.14 2.05 2.03 RN3_n100_D2 3919.34 3520.71 3936.62 3534.99 2.12 2.18 3538.40 3158.21 3551.07 3168.52 1.79 1.98 RN4_n100_D2 3477.32 2982.62 3491.33 2991.65 2.69 2.72 3038.18 2622.38 3055.51 2639.41 2.36 2.53 Avg. 3684.68 3298.48 3695.95 3308.58 2.18 2.37 3314.89 2909.01 3326.53 2921.95 2.06 2.13 RN1_n100_D3 2480.34 1996.41 2488.51 1997.98 1.83 1.7 2124.82 1635.38 2127.03 1645.98 1.58 1.29 RN2_n100_D3 2315.24 1990.81 2319.78 1996.14 1.70 1.71 1923.52 1612.65 1930.14 1618.16 1.37 1.44 RN3_n100_D3 2435.62 2051.60 2440.28 2054.78 2.09 2.03 1851.77 1569.09 1856.39 1573.41 1.32 1.36 RN4_n100_D3 2240.82 1868.86 2243.17 1870.65 1.82 2.23 1764.91 1439.53 1772.65 1454.09 1.40 1.47 Avg. 2368.00 1976.92 2372.94 1979.89 1.86 1.92 1916.26 1564.16 1921.55 1572.91 1.42 1.39 RN1_n100_D4 3534.04 2929.31 3538.32 2940.50 1.56 1.65 3044.64 2398.28 3048.74 2415.50 1.43 1.31 RN2_n100_D4 2807.13 2505.56 2827.05 2573.06 1.09 1.27 2444.41 2141.93 2456.91 2151.48 1.07 1.16 RN3_n100_D4 3161.54 2876.90 3168.41 2882.75 1.28 1.64 2742.29 2428.80 2765.12 2448.75 1.48 1.50 RN4_n100_D4 3050.61 2625.61 3058.81 2628.78 1.43 1.65 2741.18 2244.28 2748.61 2250.31 1.41 1.38 Avg. 3138.33 2734.35 3148.15 2756.27 1.34 1.55 2743.13 2303.32 2754.84 2316.51 1.35 1.34 RN1_n200_D1 4586.99 3673.59 4607.01 3698.37 1.97 1.84 3890.26 2964.06 3911.65 3002.76 1.59 1.63 RN2_n200_D1 4006.30 3338.27 4037.61 3380.21 1.41 1.76 3420.59 2774.14 3444.88 2805.14 1.32 1.62 RN3_n200_D1 4235.06 3605.87 4287.61 3632.18 1.83 1.97 3632.22 3011.05 3649.54 3032.49 1.27 1.43 RN4_n200_D1 3919.20 3226.55 3953.88 3247.98 1.62 1.95 3387.47 2541.72 3436.05 2582.96 1.52 1.59 Avg. 4186.89 3461.07 4221.53 3489.68 1.71 1.88 3582.63 2822.74 3610.53 2855.84 1.42 1.57 RN1_n200_D2 7428.48 6656.70 7468.67 6667.52 5.64 6.09 6650.45 5893.97 6692.08 5936.45 5.23 5.60 RN2_n200_D2 6883.99 6320.66 6938.59 6428.07 5.29 6.62 6220.31 5650.36 6267.01 5712.54 5.40 5.85 RN3_n200_D2 7143.03 6595.11 7201.80 6636.35 6.35 6.69 6615.83 6020.32 6681.27 6079.96 5.38 5.24 RN4_n200_D2 6468.03 5854.25 6504.63 5900.71 5.76 6.48 5799.34 5201.53 5879.89 5306.00 5.86 6.07 Avg. 6980.88 6356.68 7028.42 6408.16 5.76 6.47 6321.48 5691.54 6380.06 5758.74 5.47 5.69 RN1_n200_D3 4333.66 3631.80 4385.00 3671.93 4.12 4.57 3647.47 2941.04 3684.66 2960.16 2.97 2.99 RN2_n200_D3 4336.39 3638.83 4385.59 3681.60 4.46 5.27 3622.43 2916.56 3644.45 2943.40 3.56 4.12 RN3_n200_D3 4600.69 3735.66 4628.11 3767.51 4.49 4.84 3818.68 3096.20 3835.78 3112.43 3.53 3.35 RN4_n200_D3 4328.94 3548.88 4358.82 3569.77 5.06 5.45 3457.12 2720.75 3488.66 2749.87 3.69 3.63 Avg. 4399.92 3638.79 4439.38 3672.70 4.54 5.03 3636.43 2918.64 3663.39 2941.46 3.44 3.52 RN1_n200_D4 5863.31 5023.25 5872.71 5048.40 4.19 4.35 5056.61 4284.61 5081.22 4323.29 3.37 3.26 RN2_n200_D4 5302.13 4717.90 5314.75 4770.30 2.82 3.42 4661.26 4074.33 4685.65 4125.05 2.93 3.15 RN3_n200_D4 5842.90 5220.16 5902.48 5258.46 4.38 4.93 5141.24 4588.17 5213.82 4629.61 3.86 3.95 RN4_n200_D4 5447.38 4616.55 5478.77 4645.50 3.42 4.17 4851.89 4073.05 4890.97 4099.56 3.50 3.93 Avg. 5613.93 4894.47 5642.18 4930.67 3.70 4.21 4927.75 4255.04 4967.92 4294.38 3.41 3.57 Overall avg. 3282.36 2832.34 3298.69 2849.13 2.00 2.21 2854.27 2406.95 2870.99 2424.95 1.78 1.83 previously addressed in the literature. Computational results on test instances with different characteristics show that our algorithm performs well, in particular by finding ways to transfer more neutral customers into black routes, which lead to solutions of better quality. Furthermore, we observed that the time spent by vehicles at transfer points is very low, thus indicating that good synchronization is achieved. For the future, new ruin operators could be developed to enhance the performance of our algorithm and solve other difficult variants of vehicle routing problems. It would also be interesting to consider the integration of learning into our algorithm, for example to identify the most promising transfer points, based on the topology of the network and distribution of customers. EURO Journal on Transportation and Logistics 13 (2024) 100143 19 F.O.G. Reyes et al. CRediT authorship contribution statement Fernando O. Guillen Reyes: Data curation, Investigation, Methodology, Software, Validation, Writing – original draft. Michel Gendreau: Conceptualization, Funding acquisition, Methodology, Project administration, Resources, Supervision, Validation, Writing – review & editing. Jean-Yves Potvin: Conceptualization, Funding acquisition, Methodology, Supervision, Writing – review & editing, Project administration, Resources. Declaration of competing interest The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Jean-Yves Potvin reports financial support was provided by Natural Sciences and Engineering Research Council of Canada. If there are other authors, they declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgments The work reported in this paper was financially supported by the Natural Sciences and Engineering Research Council (NSERC) of Canada. This support is gratefully acknowledged. Appendix See Table 8. References Anderluh, A., Hemmelmayr, V., Nolz, P., 2017. Synchronizing vans and cargo bikes in a city distribution network. CEJOR Cent. Eur. J. Oper. Res. 25, 345–376. Balseiro, S., Loiseau, I., Ramonet, J., 2011. An ant colony algorithm hybridized with insertion heuristics for the time dependent vehicle routing problem with time windows. Comput. Oper. Res. 38, 954–966. Beasley, J., 1981. Adapting the savings algorithm for varying inter-customer travel times. Omega 9, 658–659. Ben Ticha, H., 2017. Vehicle Routing Problems with Road-Network Information (Ph.D. thesis). Université Clermont Auvergne. Ben Ticha, H., Absi, N., Feillet, D., Quilliot, A., 2017. Empirical analysis for the VRPTW with a multigraph representation for the road network. Comput. Oper. Res. 88, 103–116. Ben Ticha, H., Absi, N., Feillet, D., Quilliot, A., 2018. Vehicle routing problems with road-network information: State of the art. Networks 72, 393–406. Ben Ticha, H., Absi, N., Feillet, D., Quilliot, A., 2019. Multigraph modeling and adaptive large neighborhood search for the vehicle routing problem with time windows. Comput. Oper. Res. 104, 113–126. Ben Ticha, H., Absi, N., Feillet, D., Quilliot, A., Van Woensel, T., 2019. A branchand-price algorithm for the vehicle routing problem with time windows on a road network. Networks 73, 401–417. Ben Ticha, H., Absi, N., Feillet, D., Quilliot, A., Van Woensel, T., 2021. The time-dependent vehicle routing problem with time windows and road-network information. SN Oper. Res. Forum 2, 4. Bouros, P., Sacharidis, D., Dalamagas, T., Sellis, T., 2011. Dynamic pickup and delivery with transfers. In: Gertz, M., Renz, M., Zhou, X., Hoel, E., Ku, W.-S., Voisard, A., Zhang, C., Chen, H., Tang, L., Huang, Y., Lu, C.-T., Ravada, S. (Eds.), Advances in Spatial and Temporal Databases. In: Lecture Notes in Computer Science 10411, pp. 112–129. Christiaens, J., Vanden Berghe, G., 2020. Slack induction by string removals for vehicle routing problems. Transp. Sci. 54, 417–433. Crainic, T.G., Mancini, S., Perboli, G., Tadei, R., 2011. Multi-start heuristics for the two-echelon vehicle routing problem. In: Merz, P., Hao, J.-K. (Eds.), Evolutionary Computation in Combinatorial Optimization. In: Lecture Notes in Computer Science 6622, Springer, pp. 179–190. Crainic, T.G., Ricciardi, N., Storchi, G., 2009. Models for evaluating and planning city logistics systems. Transp. Sci. 43, 432–454. Dabia, S., Röpke, S., Van Woensel, T., de Kok, T., 2013. Branch and price for the time-dependent vehicle routing problem with time windows. Transp. Sci. 47, 380–391. Donati, A.V., Montemanni, R., Casagrande, N., Rizzoli, A.E., Gambardella, L.M., 2008. Time dependent vehicle routing problem with a multi ant colony system. European J. Oper. Res. 185, 1174–1191. Garaix, T., Artigues, C., Feillet, D., Josselin, D., 2010. Vehicle routing problems with alternative paths: An application to on-demand transportation. European J. Oper. Res. 204, 62–75. Gendreau, M., Ghiani, G., Guerriero, E., 2015. Time-dependent routing problems: A review. Comput. Oper. Res. 64, 189–197. Gmira, M., Gendreau, M., Lodi, A., Potvin, J.-Y., 2021. Tabu search for the timedependent vehicle routing problem with time windows on a road network. European J. Oper. Res. 288, 129–140. Grangier, P., Gendreau, M., Lehuédé, F., Rousseau, L.-M., 2016. An adaptive large neighborhood search for the two-echelon multiple-trip vehicle routing problem with satellite synchronization. European J. Oper. Res. 254, 80–91. Haghani, A., Jung, S., 2005. A dynamic vehicle routing problem with time-dependent travel times. Comput. Oper. Res. 32, 2959–2986. Hill, A.V., Benton, W.C., 1992. Modelling intra-city time-dependent travel speeds for vehicle scheduling problems. J. Oper. Res. Soc. 43, 343–351. Ichoua, S., Gendreau, M., Potvin, J.-Y., 2003. Vehicle dispatching with time-dependent travel times. European J. Oper. Res. 144, 379–396. Jia, S., Deng, L., Zhao, Q., Chen, Y., 2022. An adaptive large neighborhood search heuristic for multi-commodity two-echelon vehicle routing problem with satellite synchronization. J. Ind. Manage. Optim. 19, 1187–1210. Kafle, N., Zou, B., Lin, J., 2017. Design and modeling of a crowdsource-enabled system for urban parcel relay and delivery. Transp. Res. B 99, 62–82. Letchford, A.N., Nasiri, S.D., Oukil, A., 2014. Pricing routines for vehicle routing with time windows on road networks. Comput. Oper. Res. 51, 331–337. López-Ibánez, M., Dubois-Lacoste, J., Pérez Cáceres, L., Stützle, T., Birattari, M., 2016. The IRACE package: Iterated racing for automatic algorithm configuration. Oper. Res. Perspect. 3, 43–58. Malandraki, C., Daskin, M., 1992. Time dependent vehicle routing problems: Formulations, properties and heuristic algorithms. Transp. Sci. 26, 185–200. Minic, S., Laporte, G., 2006. The pickup and delivery problem with time windows and transshipment. INFOR 44, 217–227. Pan, B., Zhang, Z., Lim, A., 2021a. A hybrid algorithm for time-dependent vehicle routing problem with time windows. Comput. Oper. Res. 128, 105193. Pan, B., Zhang, Z., Lim, A., 2021b. Multi-trip time-dependent vehicle routing problem with time windows. European J. Oper. Res. 291, 218–231. Perboli, G., Tadei, R., Tadei, R., 2010. New families of valid inequalities for the two-echelon vehicle routing problem. Electron. Notes Discrete Math. 36, 639–646. Perboli, G., Tadei, R., Vigo, D., 2009. The two-echelon capacitated vehicle routing problem: Models and math-based heuristics. Transp. Sci. 45, 364–380. Pisinger, D., Ropke, S., 2019. Large neighborhood search. In: Gendreau, M., Potvin, J.- Y. (Eds.), Handbook of Metaheuristics, Third Edition. In: International Series in Operations Research & Management Science 272, Springer, pp. 99–127. Ropke, S., Pisinger, D., 2006. An adaptive large neighborhood search heuristic for the pickup and delivery problem with time windows. Transp. Sci. 40, 455–472. Sampaio, A., Savelsbergh, M., Veelenturf, L., Van Woensel, T., 2020. Delivery systems with crowd-sourced drivers: A pickup and delivery problem with transfers. Networks 76, 232–255. Speranza, M., Guastaroba, G., Vigo, D., 2016. Intermediate facilities in freight transportation planning: A survey. Transp. Sci. 50, 763–789. Van Belle, J., Valckenaers, P., Cattrysse, D., 2012. Cross-docking: State of the art. Omega 40, 827–846. EURO Journal on Transportation and Logistics 13 (2024) 100143 20