Pickup and delivery problem with hard time windows considering stochastic and time-dependent travel times

Wang, Zheyu; Dessouky, Maged; Van Woensel, Tom; Ioannou, Petros A. Article Pickup and delivery problem with hard time windows considering stochastic and time-dependent travel times EURO Journal on Transportation and Logistics (EJTL) Provided in Cooperation with: Association of European Operational Research Societies (EURO), Fribourg Suggested Citation: Wang, Zheyu; Dessouky, Maged; Van Woensel, Tom; Ioannou, Petros A. (2023) : Pickup and delivery problem with hard time windows considering stochastic and time-dependent travel times, EURO Journal on Transportation and Logistics (EJTL), ISSN 2192-4384, Elsevier, Amsterdam, Vol. 12, Iss. 1, pp. 1-14, https://doi.org/10.1016/j.ejtl.2022.100099 This Version is available at: https://hdl.handle.net/10419/325179 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/ EURO Journal on Transportation and Logistics 12 (2023) 100099 Available online 9 December 2022 2192-4376/© 2022 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/). Contents lists available at ScienceDirect EURO Journal on Transportation and Logistics journal homepage: www.elsevier.com/locate/ejtl Pickup and delivery problem with hard time windows considering stochastic and time-dependent travel times Zheyu Wang a,∗, Maged Dessouky b, Tom Van Woenselc, Petros Ioannou a aMing Hsieh Department of Electrical and Computer Engineering, University of Southern California, Los Angeles, USA bDaniel J. Epstein Department of Industrial and Systems Engineering, University of Southern California, Los Angeles, USA cThe Department of Industrial Engineering and Innovation Sciences, Eindhoven University of Technology, Eindhoven, The Netherlands ARTICLE INFO Keywords: Pickup and delivery problem Stochastic travel times Hard time windows Branch and price ABSTRACT Due to the uncertain nature of the traffic system, it is not trivial for delivery companies to reliably satisfy customers’ time windows. To guarantee the reliability of the pickup and delivery service under stochastic and time-dependent travel times, we consider a pickup and delivery problem with hard time windows considering stochastic and time-dependent travel times. We propose a chance-constrained model where the operational cost and the service’s reliability are considered. To quantify the service reliability, every node is associated with a desired node service level, and there exists a global service level, both measured by success probabilities. We present an estimation method for arrival times and success probabilities under stochastic travel and service times. We propose an exact solution approach based on a branch-price-and-cut framework, where a labeling algorithm generates columns. Computational experiments are conducted to assess the effectiveness of the solution framework, and Monte Carlo simulations are used to show that the proposed method can generate routes that satisfy both node and global service levels. 1. Introduction The pickup and delivery problem (PDP) is an extensively studied routing problem. Customer requests consist of two parts: a pickup at one location and a delivery at another. The pickup and delivery problem has a wide range of applications, including ridesharing services (Wang et al.,2016) and meal delivery services (Aziez et al.,2020). In pickup and delivery services, an important factor that influences customer satisfaction is the punctuality of the service. Nguyen et al. (2019) states that being able to perform the delivery service within a specified time window is one of the most important attributes in customer’s decision-making, given multiple service options. For ridesharing and shuttle services, punctuality is significant to passenger satisfaction. In most used meal delivery apps, like Uber Eats and Meituan, customers can request refunds or cancel the order if the service is later than a specific time. Therefore, companies must perform the service punctually and reliably. However, it is not trivial to satisfy customers’ time windows because of the uncertain nature of the traffic system. Traffic congestion is very common in urban areas where pickup and delivery services are frequently used. Consequently, travel times are usually unpredictable and fluctuate significantly with the time of day (Yazici et al.,2012). It ∗Correspondence to: Ming Hsieh Department of Electrical and Computer Engineering, University of Southern California, 3740 McClintock Ave, Los Angeles, CA, 90089, USA. E-mail address: [email protected] (Z. Wang). is hard for delivery companies to find efficient routing solutions that balance operational costs and service reliability. To provide reliable pickup and delivery solutions with time window constraints in an urban setting, we propose an algorithm that solves pickup and delivery problems while considering the service’s reliability under stochastic and time-dependent travel times. In a significant part of the vehicle routing literature, the travel times between certain customers are assumed to be deterministic. Most of the existing work that studies stochastic travel times adopts the assumption that the time windows are soft, meaning that the time windows can be violated at the expense of some penalty cost in the objective function. At each customer node, a penalty cost is imposed according to how early or late the vehicle arrives. If a vehicle arrives early, it can start its service right away. A limitation of the above method is that it is impractical to measure the cost of late deliveries and determine such penalty cost function in terms of money in real-world applications. The cost can be much more than the refund, for such late deliveries could impact future sales. An alternative method is to assume time windows are hard. Under the setting of hard time windows, a vehicle must wait until the start of a time window to begin its service if it arrives early. This is common for passenger services since they may not be available https://doi.org/10.1016/j.ejtl.2022.100099 Received 18 January 2022; Received in revised form 11 October 2022; Accepted 29 November 2022 EURO Journal on Transportation and Logistics 12 (2023) 100099 2 Z. Wang et al. before their scheduled pickup time. Arriving later than the end of a time window results in a service failure of the whole route. Compared to the soft time window assumption, the hard time window assumption can evaluate the reliability of the service with the probability of satisfying time windows. In this paper, we assume that time windows are hard and introduce the concept of service level, which indicates the desired probability that delivery is on-time. As customers may have different requirements and expectations of punctuality, we assume each node is associated with a node service level. We also assume a global service level measures the reliability of an overall routing solution. The global service level is a direct indicator for the service provider to evaluate its overall service reliability. The precise definitions of these terms are presented in later sections of the paper. We consider a pickup and delivery problem with hard time windows considering stochastic and time-dependent travel times (PDPHTWSTDTT). We assume the travel times between customer locations follow a specific probability distribution, and the distribution may depend on when the travel begins. The distribution is considered known and can be derived from historical data. Each customer is assumed to be associated with a hard time window and the desired node service level. There is a desired global service level. Using a chance-constrained approach, we aim to provide a routing solution that satisfies the desired service levels with the minimum operational cost. We propose an arrival time estimation method based on the works of Jula et al. (2006) and Ehmke et al. (2015). Different from the heuristic solution method proposed by Ehmke et al. (2015), our solution method is an exact solution method based on branch-cut-and-price. We propose a set partitioning formulation of the PDPHTW-STDTT that includes a route success probability constraint. Since it is inefficient and not practical to enumerate all possible routes in the pricing problem, a labeling algorithm is developed to eliminate the less promising routes to accelerate the solution process. However, given the presence of probabilistic information in the routes, deciding if one route is more promising than the other is not trivial. To deal with this challenge, exact dominance rules that handle probabilistic information are proposed. To further accelerate the algorithm, heuristic dominance rules are also proposed. Computational experiments are performed, and Monte Carlo simulations are conducted to show the algorithm’s effectiveness under stochastic travel times. The main contributions of this paper are summarized as follows. 1. We define the pickup and delivery problem with hard time windows considering stochastic and time-dependent travel times and propose a chance-constrained model for the problem. 2. We present an exact algorithm that is based on branch-cut-and- price to solve the proposed problem. New labeling algorithms and dominance rules are proposed in the pricing problem to deal with stochastic travel times and probabilistic information. 3. In the numerical experiments, we evaluate the trade-off between the reliability and the total cost of the routing solution. We also demonstrate the algorithm’s effectiveness in guaranteeing desired service levels. The paper is structured as follows. Section 2presents a review of the related literature. Section 3defines PDPHTW-STDTT and introduces the route service level estimation method. In Section 4, the solution method based on branch-cut-and-price is described. Computational results are reported in Section 5, followed by conclusions in Section 6. 2. Literature review The literature dealing with uncertainty can be categorized into stochastic vehicle routing problems (SVRP) and robust vehicle routing problems. Adulyasak and Jaillet (2016) summarizes that stochastic vehicle routing problems are applied to cases where known probability distributions can describe the uncertainty. In contrast, the robust vehicle routing problem is proposed to deal with scenarios where probability distributions are hard to estimate. Bertsimas et al. (2011) outlined that, compared to Stochastic VRP, the uncertainty model of robust VRP is not stochastic but deterministic and set-based, and its goal is to develop a solution that is feasible for all possible uncertainty realizations in the given set. Agra et al. (2013) proposed two new formulations for the robust vehicle routing problem with time windows and implemented a cutting-plane algorithm for the path inequalities formulation. Braaten et al. (2017) proposed an efficient heuristic for the robust vehicle routing problem with time windows based on adaptive large neighborhood search. Munari et al. (2019) proposed a compact formulation for the robust vehicle routing problem with time windows and first developed a branch-price-and-cut method to solve the set partitioning formulation of the problem. In our work, however, we assume the travel time distribution is known since there is sufficient work that studies the distribution of travel times in urban areas using historical data (Yazici et al.,2012), and information is also available from online databases. Therefore, our work falls into the category of SVRP, and the following literature review will focus on how to deal with time windows in the SVRP context. Recent surveys of SVRP can be found in review papers by Oyola et al. (2017,2018). Some of the work that studies SVRP does not take time windows into account, and the objectives are usually related to the duration of the routes. In the work by Laporte et al. (1992), the VRP with stochastic travel and service time was considered for the first time. Time windows are not introduced, and the intention is to limit the maximal route duration. The authors presented three different formulations based on stochastic programming and solved the problem with a branch-and- cut approach. Kenyon and Morton (2003) proposed two models with different objective functions. The first objective function minimizes the expectation of the route completion time, while the second maximizes the probability of completing all the routes before a certain time. A solution method based on branch-and-cut is proposed. Van Woensel et al. (2008) proposed a vehicle routing problem with dynamic travel times and utilized queueing theory to determine the travel times on the arcs depending on the time. Lecluyse et al. (2009) studied a vehicle routing problem with stochastic time-dependent travel times without time windows. The proposed method adjusts the objective of the classical VRP that minimizes the expected total travel time by adding the standard deviation of the travel times into the objective function. The method’s performance is evaluated on instances of up to 80 customers, which shows that the reliability of the routing solution is improved when the standard deviation of travel times is considered. In most of the literature studying SVRP with time windows, it is assumed that the time windows are soft. Soft time windows imply that vehicles that arrive earlier than the start of time windows do not have to wait, and early and late arrivals will lead to penalty costs. Ando and Taniguchi (2006) studied the VRPTW with uncertain travel times. Late arrivals generate a penalty cost proportional to the delayed time length. The objective is to minimize the total cost, including operational costs and penalties. A genetic algorithm is used to solve the problem. Li et al. (2010) investigated the VRPTW with stochastic travel and service times. One of the two proposed formulations assumes soft time windows, and a stochastic programming model with resource formulation is given. A tabu-search-based heuristic solves the problem. Taş et al. (2013,2014a,b) studied vehicle routing problems with stochastic travel times, including soft time windows under both time-independent and time-dependent settings. The objective function is to minimize a total weighted cost, which includes the penalty for delay and earliness. A method to estimate the expectation and variance of arrival times is proposed. In Taş et al. (2013,2014a), the problems are solved by a tabu-search based heuristic. In Taş et al. (2014b), a solution method based on column generation and branch-and-price solution is proposed. Few papers studied the SVRP with hard time windows, which are the most closely related references to our work. In those papers, the key point is to guarantee a certain service level at each customer node. EURO Journal on Transportation and Logistics 12 (2023) 100099 3 Z. Wang et al. A major challenge is to deal with the truncation of the probability distributions of arrival times caused by the hard time windows, as vehicles have to wait for the start of the time window to begin service. Different methods are proposed to estimate the probability distributions of arrival times after the truncation. The work of Jula et al. (2006) is one of the first efforts that investigate the effect of hard time windows. Both travel times and service times are considered to be stochastic. A Taylor series expansion approximates the mean and variance of arrival time at each customer node. The estimation method deals with nonstationary time-varying travel time distributions. Chebyshev and Chernoff bounds determine if a route meets the required service level. The proposed problem is a traveling salesman problem involving only one vehicle. Li et al. (2010) also proposed a chance-constrained programming formulation for VRP with stochastic travel times, which differs from the above formulation. In the chance-constrained programming model, the service level at each customer node is guaranteed. The probability of route success is derived directly from Monte Carlo simulations without explicit estimation expressions. Ehmke et al. (2015) studied a VRP with stochastic travel times and hard time windows. Using extreme value theory, the authors proposed a different estimation method to approximate the mean and variance of arrival times. When estimating the route success probability after estimating the mean and variance of arrival times, it is approximated under the assumption that the distributions of the arrival times are normal. The authors defined the route feasibility check to solve the problem and integrated it into existing heuristics like tabu-search. To show the effectiveness of the proposed method, the notion of lateness is defined to evaluate the reliability of a route. Miranda and Conceição (2016) proposed a method to approximate the distribution function of arrival times using convolution functions. A metaheuristic is proposed to solve the VRP with stochastic travel times. The authors conducted numerical experiments on instances with up to 100 customers. Gutierrez et al. (2018) extended the arrival times estimation method by Ehmke et al. (2015) by integrating stochastic service times. Several different confidence levels are defined, and a multi-population memetic algorithm is proposed. All of the work above used heuristics to solve the problem. A limitation of the above heuristic solution methods is that even though the service levels at customer nodes can be guaranteed, there is no guarantee on the global service level. Errico et al. (2018) proposed an exact solution based on a branch-cut-and-price method to solve a VRP with hard time windows. The goal is to guarantee a global service level, which is the probability that all the customer nodes are served on time. However, the stochasticity lies in service times, not travel times. Unlike the work mentioned above, the probability distribution functions of service times are assumed to be discrete. Therefore, the probability of a route’s success can be computed directly, and no estimation is needed. The authors presented a set partitioning formulation with a probability constraint and developed a branch-cut-and-price solution framework. Numerical experiments are performed on instances with up to 50 customers. Another related research area is the time-dependent vehicle routing problem (TDVRP), whose survey can be found in Gendreau et al. (2015). Ichoua et al. (2003) first introduced the ‘‘first in, first out’’ (FIFO) property under the assumption of time-dependent travel times and proposed a step-wise function speed model which satisfies the property. Dabia et al. (2013) proposed a branch-and-price algorithm for time-dependent vehicle routing problems with time windows, which can solve instances with up to 100 customers. Sun et al. (2018b) studied the time-dependent pickup and delivery problem with time windows and its variants where some customers can be skipped and the departure time can be flexible. The definition of time-dependent under stochastic travel times assumption is slightly different from the traditional definition used by the above papers and is given in the next section. As for the solution method, our work adopts the branch-cut-and- price framework, which is an exact solution method. The readers Fig. 1. The position of the proposed problem in the literature. interested in the branch-cut-and-price method in VRP applications are referred to a survey by Costa et al. (2019). We study the pickup and delivery problem more restrictive than the vehicle routing problem because the PDP requires the pickup node and the delivery node served by the same vehicle. Therefore, the branch-cut-and-price procedure for solving PDP varies from the VRP version to some degree. The interested readers are referred to the work by Ropke and Cordeau (2009), which provides a detailed branch-cut-and-price solution method for the PDPTW. Fig. 1 summarizes the position of the proposed problem in the literature and its relation with similar problems. The proposed method also belongs to time-dependent vehicle routing problems (TDVRP). Our work contributes to the literature and differs from the other work that studies stochastic travel times VRP with hard time windows in the following ways. First, our work integrates time-dependent travel times. Second, to our knowledge, we are the first to propose an exact solution method by proposing a labeling algorithm with probabilistic information. 3. The PDPHTW-STDTT In this section, we first introduce the definition of the pickup and delivery problem with hard time windows considering stochastic and time-dependent travel times. Then, the method to estimate the means and variances of arrival times is presented. Finally, we show how the route success probability is estimated. 3.1. Problem description Let 𝐺= (𝑁, 𝐴)be a directed graph, where 𝑁= {0,1,…,2𝑛+ 1} is the set of nodes and 𝐴is the set of arcs 𝐴= {(𝑖, 𝑗)|𝑖, 𝑗 ∈𝑁}. Nodes 0 and 2𝑛+ 1 represent the origin and destination depot of the vehicles. The set 𝑁𝑃= {1,…, 𝑛} ∈ 𝑁represents the set of pickup nodes and 𝑁𝐷= {𝑛+ 1,…,2𝑛} ∈ 𝑁represents the set of delivery nodes. There is a set of n requests. Each request 𝑖is associated with a pickup node 𝑖and a delivery node 𝑖+𝑛. For each request, the pickup node must be visited before the delivery node, and both nodes must be visited by the same vehicle. Each node must be visited only once. We consider a capacitated pickup and delivery problem, where each request is associated with load size, 𝑞𝑖. Pickup node 𝑖has a load of 𝑞𝑖and delivery node 𝑛+𝑖has a load of −𝑞𝑖(drop off). As there is no inventory at the depot, 𝑞0=𝑞2𝑛+1 = 0. Each node 𝑖∈𝑁𝑝∪𝑁𝑑has a time window [𝑒𝑖, 𝑙𝑖], where 𝑒𝑖is the start of the time window and 𝑙𝑖is the end of the time window. The vehicle arrives earlier than 𝑒𝑖must wait until 𝑒𝑖 to start its service at node 𝑖.𝑑𝑖𝑗 denotes the distance between node 𝑖 and node 𝑗. A travel cost 𝑐𝑖𝑗 is associated with each arc (𝑖, 𝑗) ∈ 𝐴. For the sake of simplicity, we assume the travel cost is proportional to the travel distance and 𝑐𝑖𝑗 =𝑑𝑖𝑗 in this paper. EURO Journal on Transportation and Logistics 12 (2023) 100099 4 Z. Wang et al. There is a fleet of identical vehicles, each with a maximum capacity of 𝑄. We assume there is no limit on the number of vehicles used. However, a fixed cost 𝐶𝑓𝑖𝑥𝑒𝑑 is paid when a vehicle is used. All the vehicles depart at the depot at time 0 and finish the route at the depot. We assume the travel times between nodes and service times at nodes are stochastic. The travel time on each arc (𝑖, 𝑗) ∈ 𝐴is a Gaussian random process 𝑋𝑖𝑗 (𝑡), with 𝑡representing the time when the vehicle enters arc (𝑖, 𝑗). We assume the mean of 𝑋𝑖𝑗 (𝑡),𝐸[𝑋𝑖𝑗 (𝑡)] is time-dependent, while the variance of 𝑋𝑖𝑗 (𝑡),𝑉 𝑎𝑟[𝑋𝑖𝑗 ]is time-independent, only proportional to the length of the arc, 𝑑𝑖𝑗 . Following Ehmke et al. (2015)’s notation, we define a variation coefficient 𝑐𝑣, so that 𝑉 𝑎𝑟[𝑋𝑖𝑗 ] = 𝑐𝑣𝑑𝑖𝑗 . We assume that stochasticity is only in the travel times but not path choices. Therefore, there is only one path between any two nodes. It is assumed that the travel times on an individual arc at any particular time 𝑡are statistically independent. Most of the literature that studies time-dependent travel times satisfies the ‘‘first in, first out’’ (FIFO) property. The FIFO property under deterministic travel times states that if two vehicles traverse the same arc, the one that leaves first will arrive first. Most literature adopts a speed profile model that consists of step-wise functions and satisfies the FIFO principle proposed by Ichoua et al. (2003). However, when it comes to stochastic and timedependent travel times, the FIFO property cannot be guaranteed at the single-vehicle level because of the uncertainty in travel times. Nevertheless, based on real-life observations and experiences, on average, vehicles that leave early will arrive early. Therefore, the FIFO property should hold at the macro level, meaning that 𝐸[𝑋𝑖𝑗 (𝑡)] satisfies the FIFO principle. We assume that 𝐸[𝑋𝑖𝑗 (𝑡)] is derived from the speed profile model mentioned above. The speed profile model assumes that the scheduling horizon is divided into 𝑘intervals, and the vehicle’s speed is a constant in each of the intervals. Each arc (𝑖, 𝑗) ∈ 𝐴is associated with a speed profile that consists of a constant speed for each interval. We assume all the arcs share the same speed profile, meaning that the speeds on all the arcs are the same at a specific time. The 𝐸[𝑋𝑖𝑗 (𝑡)] derived by the above method is a piece-wise linear function. Details about the FIFO property and the definition of the speed profile can be found in Ichoua et al. (2003) and Sun et al. (2018a). We note that the proposed method can deal with more complicated travel time assumptions where travel speeds differ on different arcs or the travel times follow a more complicated changing pattern. For the sake of simplicity, we adopt the simplified assumption above. The service time at each node is a random variable that is timeindependent. The service time 𝑆𝑖at node 𝑖has mean 𝐸[𝑆𝑖]and variance 𝑉 𝑎𝑟[𝑆𝑖]. It is assumed that 𝑆𝑖follows a normal distribution. Once the service is finished, the vehicle starts to travel to the next node on the route. A route 𝑟is defined by a sequence of nodes 𝑟= (𝑣0, 𝑣1,…, 𝑣𝑚, 𝑣𝑚+1), where 𝑣1,…, 𝑣𝑚∈𝑁𝑝∪𝑁𝑑and 𝑣0, 𝑣𝑚+1 represent the depot. We define a node as successful if the vehicle arrives at the node within its time window. Node success probability is defined by the probability that a node in a route is successful when given the route and all needed distributions. We define that a route is successful if the vehicle arrives at each node on the route within their time windows. Route success probability is the probability that a given route is successful when given all needed distributions. Each node 𝑖∈𝑁𝑝∪𝑁𝑑is associated with a node service level 𝜃𝑗(0≤𝜃𝑗<1), indicating the desired node success probability by the corresponding customer. We also define a global service level 𝛩(0≤𝛩 < 1), indicating the desired probability that all the routes are successful. The objective is to minimize the total costs when the node and global service levels are satisfied. 3.2. Arrival times estimation Several methods for estimating the arrival times under stochastic travel time and hard time windows assumption are mentioned in our literature review. We adopt the method proposed by Jula et al. (2006), because of its capability of dealing with time-dependent travel times. In addition, this estimation method makes it more practical to generate dominance rules for the proposed labeling algorithm, presented later in this paper. A comparison of the adopted method and Ehmke et al. (2015)’s method can be found in Appendix. The main idea of Jula et al. (2006)’s method is to estimate the first and second moments of the arrival times iteratively. Let 𝐴𝑟 𝑖denote the arrival time at node 𝑖on route 𝑟, and 𝐷𝑟 𝑖denote departure time after service at node 𝑖on route 𝑟. The means and variances of 𝑋𝑖𝑗 (𝑡) and 𝑆𝑖are known. When the vehicle travels from node 𝑖to node 𝑗, the expressions below give the means and variances of the arrival and departure times at node 𝑗. The estimations of the means and variances of arrival times and departure times of all the nodes on a route can be derived iteratively by repeating this process. As in Jula et al. (2006), the following equations can be derived: 𝐸[𝐴𝑟 𝑗] ≈ 𝐸[𝐷𝑟 𝑖] + 𝐸[𝑋𝑖𝑗 (𝐸[𝐷𝑟 𝑖])] (1) 𝑉 𝑎𝑟[𝐴𝑟 𝑗] ≈ {1 + 𝐸′[𝑋𝑖𝑗 (𝐸[𝐷𝑟 𝑖])]}2𝑉 𝑎𝑟[𝐷𝑟 𝑖] + 𝑉 𝑎𝑟[𝑋𝑖𝑗 ](2) 𝐸(𝐷𝑟 𝑗) = 𝑔𝑗(𝐸[𝐴𝑟 𝑗]) + 𝐸[𝑆𝑗](3) 𝑉 𝑎𝑟(𝐷𝑟 𝑗) ≈ 1 4(∫𝐸[𝐴𝑟 𝑗]+𝜎[𝐴𝑟 𝑗] 𝐸[𝐴𝑟 𝑗]−𝜎[𝐴𝑟 𝑗] 𝑔′ 𝑗(𝑥)𝑑𝑥)2 +𝑉 𝑎𝑟[𝑆𝑗](4) 𝑔𝑗(𝑡) = ⎧ ⎪ ⎨ ⎪ ⎩ 𝑒𝑗𝑖𝑓 𝑡 ≤𝑒𝑗 𝑡 𝑖𝑓 𝑒𝑗< 𝑡 ≤𝑙𝑗 𝑀(𝑡−𝑙𝑗) + 𝑙𝑗𝑖𝑓 𝑡 > 𝑙𝑗 (5) 𝑔′ 𝑗(𝑡) = ⎧ ⎪ ⎨ ⎪ ⎩ 0𝑖𝑓 𝑡 ≤𝑒𝑗 1𝑖𝑓 𝑒𝑗< 𝑡 ≤𝑙𝑗 𝑀 𝑖𝑓 𝑡 > 𝑙𝑗 (6) where 𝐸′[𝑋𝑖𝑗 (𝑡)] is the derivative of 𝐸[𝑋𝑖𝑗 (𝑡)] with respect to 𝑡.𝜎[𝐴𝑟 𝑗]is the standard deviation of 𝐴𝑟 𝑗.(5) and (6) show that 𝑔𝑗(𝑡)is a piece-wise linear function to model the effect of the hard time window at node 𝑗, and 𝑔′ 𝑗(𝑡)is the derivative of 𝑔𝑗(𝑡)with respect to 𝑡.𝑀represents a very big number. Expression (1) and (2) estimate the expectation and variance of arrival time at a node based on the expectation and variance of departure time at the previous node and travel time between the nodes. Expression (3) and (4) calculate the expectation and variance of departure time at a node based on the expectation and variance of arrival time and service time at the same node. 3.3. Route success probability estimation Given the mean and the variance of the arrival time at a node, the node success probability can be estimated. We adopt the estimation method in Ehmke et al. (2015) that assumes the arrival times conform to a normal distribution. Note that the actual distribution of arrival times may be skewed because of hard time windows, and the actual distribution cannot be computed precisely. Thus, the node success probability at node 𝑖on route 𝑟, denoted by 𝑃𝑟 𝑖, is estimated by 𝑃𝑟 𝑖=𝑃(𝐴𝑟 𝑖≤𝑙𝑖) ≈ 𝛷(𝑙𝑖−𝐸[𝐴𝑟 𝑖] √𝑉 𝑎𝑟[𝐴𝑟 𝑖])(7) where 𝛷(⋅)is the cumulative distribution function of the standard normal distribution. The estimated route success probability 𝑃𝑟is the product of the node success probabilities of the nodes on the route and is given by the following equation: 𝑃𝑟=∏ 𝑖∈𝑟 𝑃𝑟 𝑖(8) A route is considered feasible when all the node service levels of the nodes on the route and all constraints in the pickup and delivery problem are satisfied. Let be the set of all feasible routes. Let EURO Journal on Transportation and Logistics 12 (2023) 100099 5 Z. Wang et al. 𝑐ℎ𝑜𝑠𝑒𝑛 = {𝑟1, 𝑟2,…, 𝑟𝑘}be the set of routes that are chosen in the routing solution. Global success probability 𝑃𝑔𝑙𝑜𝑏𝑎𝑙 is the probability that all chosen routes are successful. To satisfy the desired node service levels and the desired global service level, we have the following constraints: 𝑃𝑟 𝑖≥𝜃𝑖∀𝑟∈𝑐ℎ𝑜𝑠𝑒𝑛,∀𝑖∈𝑟(9) 𝑃𝑔𝑙𝑜𝑏𝑎𝑙 =∏ 𝑟∈𝑐ℎ𝑜𝑠𝑒𝑛 𝑃𝑟≥𝛩(10) where 𝜃𝑖is the node service level of node 𝑖, and 𝛩is the global service level. Expression (9) states that each node’s success probability satisfies its node service level. Expression (10) states that the global success probability satisfies the global service level. 4. Solution approach This section presents a solution method for PDPHTW-STDTT based on branch-cut-and-price. We first provide the set partitioning formulation of the problem. Then, we emphasize how to conduct column generation and solve the pricing problem. Finally, we introduce the branch-cut-and-price framework, the cutting plane method, and the branching strategies. 4.1. Set partitioning formulation We define the cost of a route 𝑟= (𝑣0, 𝑣1,…, 𝑣𝑚, 𝑣𝑚+1)as 𝑐𝑟= 𝐶𝑓𝑖𝑥𝑒𝑑 +∑𝑚 𝑖=0 𝑐𝑣𝑖𝑣𝑖+1 . Let the parameter 𝑎𝑖𝑟 denotes if node 𝑖is visited on route 𝑟.𝑎𝑖𝑟 = 1 if route 𝑟visits node 𝑖, otherwise 𝑎𝑖𝑟 = 0. The binary variable 𝑥𝑟indicates if route 𝑟is chosen in the solution. The PDPHTW-STDTT can be formed in the following set partitioning formulation. 𝑚𝑖𝑛𝑖𝑚𝑖𝑧𝑒 ∑ 𝑟∈𝑐𝑟𝑥𝑟(11) 𝑠.𝑡 ∑ 𝑟∈𝑎𝑖𝑟𝑥𝑟= 1 ∀𝑖∈𝑁𝑝(12) ∑ 𝑟∈𝑥𝑟𝑙𝑛(𝑃𝑟)≥𝑙𝑛(𝛩)(13) 𝑥𝑟∈ {0,1} ∀𝑟∈(14) The objective function (11) minimizes the total cost, including fixed vehicle and arc costs. Constraint (12) indicates that each pickup node is visited exactly once. As we pair a pickup node with its corresponding delivery node in the routes, each delivery node is guaranteed to be visited exactly once. Constraint (13) ensures the global service level is satisfied and is derived from taking the logarithm of both sides of (10). 4.2. Column generation The set of all feasible routes is large and cannot be explicitly enumerated. We utilize a column generation method (Desaulniers et al., 2006) to generate routes iteratively. Only a subset of is maintained in the linear program. In each iteration, a relaxation of the restricted master problem (RMP) is solved, generating the dual multipliers of constraints (12) and (13),𝜋𝑖(∀𝑖∈𝑁𝑝)and 𝜆respectively. Then, the pricing algorithm is called to find routes with negative reduced costs, which are added to the subset of in the RMP. The iteration stops when no such route is found. The subset of is initialized with the set of routes that visit only one pair of pickup and delivery nodes. The reduced cost 𝑐𝑟of a route 𝑟∈is calculated as: 𝑐𝑟=𝑐𝑟−∑ 𝑖∈𝑁𝑝 𝑎𝑖𝑟𝜋𝑖+𝑙𝑛(𝑃𝑟)𝜆(15) =𝐶𝑓𝑖𝑥𝑒𝑑 + 𝑚 ∑ 𝑖=0 𝑐𝑣𝑖𝑣𝑖+1 −∑ 𝑖∈𝑁𝑝 𝑎𝑖𝑟𝜋𝑖+∑ 𝑖∈𝑟 𝑙𝑛(𝑃𝑟 𝑖)𝜆(16) 4.3. Pricing algorithm The pricing algorithm finds routes with negative reduced cost (15). The pricing problem is a variant of the shortest path problems with resource constraints (ESPPRC, Irnich and Desaulniers,2005) and is mainly solved by a labeling algorithm (Righini and Salani,2008). The main idea of the labeling algorithm is to represent every partial route that starts at node 0 and ends at any node 𝑖with a label 𝐿. The label contains information about the route, including the current load on the vehicle, the delivery tasks that need to be finished, the accumulated reduced cost, and the arrival times and probability estimations mentioned in Section 3, etc. Labels are extended from existing labels toward the destination node. During the extension process, we check if the new label satisfies the constraints in the assumption of the problem, like the capacity constraint, node service level constraint, etc. In other words, the routes generated by the labeling algorithm are feasible. To speed up the algorithm, dominance rules are usually proposed to eliminate routes that are not useful. Note that a bidirectional labeling algorithm is computationally more efficient than a unidirectional labeling algorithm. We only consider a forward labeling algorithm. This is because the probability information cannot be extended backward. 4.3.1. Definition of labels A partial route 𝑟𝑘= (𝑣0, 𝑣1,…, 𝑣𝑚)is a route that does not necessarily end at the destination. As the number of partial routes may become large in the labeling algorithm, we use superscript 𝑘to differentiate them. Each partial route 𝑟𝑘is associated with a label 𝐿𝑘. 𝐿𝑘= [(𝐿𝑘),(𝐿𝑘),(𝐿𝑘),(𝐿𝑘), 𝐸𝐴(𝐿𝑘), 𝑉𝐴(𝐿𝑘), 𝐸𝐷(𝐿𝑘), 𝑉𝐷(𝐿𝑘), 𝑃 𝑟𝑜𝑏(𝐿𝑘),(𝐿𝑘),(𝐿𝑘), 𝜋(𝐿𝑘)]. The components of 𝐿𝑘are explained below: •(𝐿𝑘)is the last node visited on partial route 𝑟𝑘. •(𝐿𝑘)is the partial route 𝑟𝑘, which contains all the nodes visited in the order that they are visited. •(𝐿𝑘)is the load in the vehicle after visiting node (𝐿𝑘). •(𝐿𝑘)is the accumulated cost along partial route 𝑟𝑘after visiting node (𝐿𝑘), without considering the fixed cost. •For the necessary information to estimate arrival times and success probabilities, 𝐸𝐴(𝐿𝑘), 𝑉𝐴(𝐿𝑘)are the estimated mean and variance of the arrival time at node (𝐿𝑘). 𝐸𝐷(𝐿𝑘), 𝑉𝐷(𝐿𝑘)are the estimated mean and variance of the departure time at node (𝐿𝑘). 𝑃 𝑟𝑜𝑏(𝐿𝑘)is the estimated route success probability of the partial route 𝑟𝑘. •(𝐿𝑘)⊆ 𝑁𝑝is the set of pickup nodes that have been visited in partial route 𝑟𝑘, whose corresponding delivery nodes have not been visited. •(𝐿𝑘)⊆ 𝑁𝑝is the set of unreachable requests . A node 𝑖is said to be unreachable if: 1. node 𝑖has already been visited on partial route 𝑟𝑘(so that it will not be visited again), or 2. going directly from current node (𝐿𝑘)to pickup node 𝑖 cannot satisfy the service levels, which is 𝑃 𝑟𝑜𝑏(𝐿𝑘)⋅𝑃𝑟𝑘⨁𝑖 𝑖< 𝜃𝑖, or 𝑃 𝑟𝑜𝑏(𝐿𝑘)⋅𝑃𝑟𝑘⨁𝑖 𝑖< 𝛩. Here, 𝑟𝑘⨁𝑖represents the new partial route after extending partial route 𝑟𝑘to node 𝑖. •𝜋(𝐿𝑘)is the sum of the dual values associated with Constraint (12) on the partial route 𝑟𝑘. 4.3.2. Label extension Given a label 𝐿𝑘associated with partial route 𝑟𝑘, its extension to node 𝑖is feasible only if: (𝐿𝑘) + 𝑞𝑖≤𝑄∀𝑖∈𝑁(17) 𝑖∉(𝐿𝑘) ∀𝑖∈𝑁𝑝(18) EURO Journal on Transportation and Logistics 12 (2023) 100099 6 Z. Wang et al. 𝑖−𝑛∈(𝐿𝑘) ∀𝑖∈𝑁𝑑(19) (𝐿𝑘)=∅ 𝑖𝑓 𝑖 = 2𝑛+ 1 (20) 𝑃 𝑟𝑜𝑏(𝐿𝑘)⋅𝑃𝑟𝑘⨁𝑖 𝑖≥𝑚𝑎𝑥{𝜃𝑖, 𝛩}(21) Condition (17) guarantees the vehicle capacity constraint will not be violated. Conditions (18)–(20) make sure the precedence constraint of the pickup and delivery problem is satisfied. Condition (21) ensures that the node and global service levels are satisfied. If it is feasible to extend an existing label 𝐿𝑘to node 𝑖, a new label 𝐿ℎis created, representing the new partial route 𝑟𝑘⨁𝑖. The label updating rules are: (𝐿ℎ) = 𝑖(22) (𝐿ℎ) = {(𝐿𝑘)∪{𝑖}𝑖𝑓 𝑖 ∈𝑁𝑝 (𝐿𝑘)∖{𝑖−𝑛}𝑖𝑓 𝑖 ∈𝑁𝑑 (23) (𝐿ℎ) = {(𝐿𝑘)∪{𝑖}𝑖𝑓 𝑖 ∈𝑁𝑝 (𝐿𝑘)𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒 (24) (𝐿ℎ) = (𝐿𝑘) + 𝑞𝑖(25) (𝐿ℎ) = (𝐿𝑘) + 𝑐(𝐿𝑘)𝑖(26) 𝜋(𝐿ℎ) = {𝜋(𝐿𝑘) + 𝜋𝑖𝑖𝑓 𝑖 ∈𝑁𝑝 𝜋(𝐿𝑘)𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒 (27) 𝐸𝐴(𝐿ℎ) = 𝐸𝐷(𝐿𝑘) + 𝐸[𝑋(𝐿𝑘)𝑖(𝐸𝐷(𝐿𝑘))] (28) 𝑉𝐴(𝐿ℎ) = (1 + 𝐸′[𝑋(𝐿𝑘)𝑖(𝐸𝐷(𝐿𝑘))])2𝑉𝐷(𝐿𝑘) + 𝑉 𝑎𝑟[𝑋(𝐿𝑘)𝑖](29) 𝐸𝐷(𝐿ℎ) = 𝑔𝑖(𝐸𝐴(𝐿ℎ)) + 𝐸[𝑆𝑖](30) 𝑉𝐷(𝐿ℎ) = 1 4(∫𝐸𝐴(𝐿ℎ)+√𝑉𝐴(𝐿ℎ) 𝐸𝐴(𝐿ℎ)−√𝑉𝐴(𝐿ℎ) 𝑔′ 𝑖(𝑥)𝑑𝑥)2+𝑉 𝑎𝑟[𝑆𝑖](31) 𝑃 𝑟𝑜𝑏(𝐿ℎ) = 𝑃 𝑟𝑜𝑏(𝐿𝑘)𝑃𝑟𝑘⨁𝑖 𝑖=𝑃 𝑟𝑜𝑏(𝐿𝑘)𝛷(𝑙𝑖−𝐸𝐴(𝐿ℎ) √𝑉𝐴(𝐿ℎ))(32) Expressions (22)–(27) are derived from the definition of the label. Expressions (28)–(32) are derived from the conclusions in Section 3. 4.3.3. Label dominance rules The concept of label dominance is adopted to restrict the total number of labels and speed up the algorithm. A dominance check is performed within the set of labels with the same ending node. If another dominates a label, the dominated label will be discarded. Definition of dominance: Consider two feasible partial routes 𝑟1 and 𝑟2, with labels 𝐿1and 𝐿2, respectively. If (𝐿1) = (𝐿2), a dominance check with the following principles will be performed. We say 𝐿1dominates 𝐿2so that 𝐿2can be eliminated, if: (1) If extending 𝑟2to any node 𝑖generates a feasible route 𝑟2⨁𝑖, then extending 𝑟1to node 𝑖also generates a feasible route 𝑟1⨁𝑖. (2) For any such extensions 𝑟1⨁𝑖and 𝑟2⨁𝑖, the reduced cost of 𝑟1⨁𝑖is always less than or equal to the reduced cost of 𝑟2⨁𝑖. We first propose and prove Lemma 1, which will be used to prove the exact dominance rules. Lemma 1. For partial route 𝑟1associated with label 𝐿1and 𝑟2associated with label 𝐿2, if 𝐸𝐴(𝐿1)≤𝐸𝐴(𝐿2),𝑃 𝑟𝑜𝑏(𝐿1)≥𝑃 𝑟𝑜𝑏(𝐿2)and (1 + 𝑓′ 𝑚𝑎𝑥)2𝑉𝐷(𝐿1)≤(1 + 𝑓′ 𝑚𝑖𝑛)2𝑉𝐷(𝐿2), then 𝑃𝑟1⨁𝑖 𝑖≥𝑃𝑟2⨁𝑖 𝑖holds for feasible extensions to any node 𝑖, where 𝑓′ 𝑚𝑎𝑥 and 𝑓′ 𝑚𝑖𝑛 are the maximum and minimum value of the derivative of the expected travel time with respect to time, over all arcs and over the scheduling horizon. Proof of Lemma 1.Assume 𝑟1⨁𝑖and 𝑟2⨁𝑖are feasible extensions of 𝑟1and 𝑟2to an arbitrary node 𝑖, and the associated labels are 𝐿1∗ and 𝐿2∗, respectively. We first show that 𝐸𝐴(𝐿1∗)≤𝐸𝐴(𝐿2∗). Given expression (30): 𝐸𝐷(𝐿ℎ) = 𝑔𝑖(𝐸𝐴(𝐿ℎ)) + 𝐸[𝑆𝑖] where 𝑔𝑖(⋅)is non-decreasing, 𝐸𝐴(𝐿1)≤𝐸𝐴(𝐿2), and 𝐸[𝑆𝑗]is constant, we derive that 𝐸𝐷(𝐿1)≤𝐸𝐷(𝐿2). With 𝐸𝐷(𝐿1)≤𝐸𝐷(𝐿2), and expression (28): 𝐸𝐴(𝐿ℎ) = 𝐸𝐷(𝐿𝑘) + 𝐸[𝑋(𝐿𝑘)𝑖(𝐸𝐷(𝐿𝑘))] where we assume the FIFO property of travel time function 𝐸[𝑋𝑖𝑗 (𝑡)], we get 𝐸𝐴(𝐿1∗)≤𝐸𝐴(𝐿2∗). We then show that 𝑉𝐴(𝐿1∗)≤𝑉𝐴(𝐿2∗). In the Lemma statement, we have: (1 + 𝑓′ 𝑚𝑎𝑥)2𝑉𝐷(𝐿1)≤(1 + 𝑓′ 𝑚𝑖𝑛)2𝑉𝐷(𝐿2)(33) Because of the definition of 𝑓′ 𝑚𝑎𝑥 and 𝑓′ 𝑚𝑖𝑛, we have: 𝐸′[𝑋(𝐿1)𝑖(𝐸𝐷(𝐿1))] ≤𝑓′ 𝑚𝑎𝑥 (34) 𝐸′[𝑋(𝐿2)𝑖(𝐸𝐷(𝐿2))] ≥𝑓′ 𝑚𝑖𝑛 (35) From inequalities (33)–(35), it can be derived that (1 + 𝐸′[𝑋(𝐿1)𝑖(𝐸𝐷(𝐿1))])2𝑉𝐷(𝐿1)≤(1 + 𝐸′[𝑋(𝐿2)𝑖(𝐸𝐷(𝐿2))])2𝑉𝐷(𝐿2) (36) Given (36), and expression (29): 𝑉𝐴(𝐿ℎ) = (1 + 𝐸′[𝑋(𝐿𝑘)𝑖(𝐸𝐷(𝐿𝑘))])2𝑉𝐷(𝐿𝑘) + 𝑉 𝑎𝑟[𝑋(𝐿𝑘)𝑖] where 𝑉 𝑎𝑟[𝑋𝑖𝑗 ]is time-independent, we have 𝑉𝐴(𝐿1∗)≤𝑉𝐴(𝐿2∗). Since 𝐸𝐴(𝐿1∗)≤𝐸𝐴(𝐿2∗)and 𝑉𝐴(𝐿1∗)≤𝑉𝐴(𝐿2∗), we have 𝑙𝑖−𝐸𝐴(𝐿1∗) √𝑉𝐴(𝐿1∗) ≥𝑙𝑖−𝐸𝐴(𝐿2∗) √𝑉𝐴(𝐿2∗) . Since 𝛷(⋅)is monotonic increasing, 𝑃𝑟1⨁𝑖 𝑖= 𝛷(𝑙𝑖−𝐸𝐴(𝐿1∗) √𝑉𝐴(𝐿1∗))≥𝛷(𝑙𝑖−𝐸𝐴(𝐿2∗) √𝑉𝐴(𝐿2∗))=𝑃𝑟2⨁𝑖 𝑖.□ Then, the exact pricing rule is proposed. If all eight conditions in Proposition 1 are satisfied, it is guaranteed that one label dominates the other label, according to the definition of dominance. Proposition 1 (Exact Pricing).Label 𝐿2is dominated by 𝐿1, if the following 8 conditions are satisfied: (1) (𝐿1) = (𝐿2). (2) (𝐿1)⊆(𝐿2). (3) (𝐿1) = (𝐿2). Note that under the stochastic and time-dependent assumption, the triangle inequality for travel times does not necessarily hold, so this is different from (𝐿1)⊆(𝐿2)proposed in Ropke and Cordeau (2009). (4) (𝐿1)≤(𝐿2). (5) (𝐿1) − 𝜋(𝐿1) + 𝜆𝑃 𝑟𝑜𝑏(𝐿1)≤(𝐿2) − 𝜋(𝐿2) + 𝜆𝑃 𝑟𝑜𝑏(𝐿2). (6) 𝐸𝐴(𝐿1)≤𝐸𝐴(𝐿2). (7) (1 + 𝑓′ 𝑚𝑎𝑥)2𝑉𝐷(𝐿1)≤(1 + 𝑓′ 𝑚𝑖𝑛)2𝑉𝐷(𝐿2). (8) 𝑃 𝑟𝑜𝑏(𝐿1)≥𝑃 𝑟𝑜𝑏(𝐿2). Proof of Proposition 1.We prove that it is sufficient to derive dominance relationship between 𝐿1and 𝐿2from Proposition 1. In other words, we want to show the definition of dominance is satisfied if given all the conditions in Proposition 1. Assume there exists partial routes 𝑟1 and 𝑟2whose corresponding labels 𝐿1and 𝐿2satisfy conditions (1)–(8) in Proposition 1, and 𝑖is an arbitrary node. Regarding number (1) in the definition of dominance, we want to show if 𝑟2⨁𝑖is feasible based on the definition of feasibility in (17)– (21),𝑟1⨁𝑖is also feasible, for any node 𝑖. For the vehicle capacity constraint, the load associated with partial route 𝑟1⨁𝑖does not violate the capacity constraint (17) if 𝑟2⨁𝑖satisfies the requirement because of condition (4). For the route success probability constraint, condition (8) gives 𝑃 𝑟𝑜𝑏(𝐿1)≥𝑃 𝑟𝑜𝑏(𝐿2), and Lemma 1 gives 𝑃𝑟1⨁𝑖 𝑖≥𝑃𝑟2⨁𝑖 𝑖 for all nodes 𝑖. Therefore, 𝑃 𝑟𝑜𝑏(𝐿1)𝑃𝑟1⨁𝑖 𝑖≥𝑃 𝑟𝑜𝑏(𝐿2)𝑃𝑟2⨁𝑖 𝑖states that 𝑟1⨁𝑖does not violate the probability constraint (21) provided that 𝑟2⨁𝑖satisfies the constraints. As for precedence constraints of the EURO Journal on Transportation and Logistics 12 (2023) 100099 7 Z. Wang et al. pickup and delivery problem, conditions (2) and (3) guarantee that if 𝑟2⨁𝑖satisfies (18)–(20),𝑟1⨁𝑖will also satisfy. Thus, for any node 𝑖, if 𝑟2⨁𝑖is feasible, 𝑟1⨁𝑖is also feasible. Regarding number (2) in the definition, let 𝑐(𝐿𝑘)be the reduced cost of 𝐿𝑘as defined in (15). We assume the labels for new feasible partial routes 𝑟1⨁𝑖and 𝑟2⨁𝑖are 𝐿1∗and 𝐿2∗. We show that 𝑐(𝐿1∗)≤𝑐(𝐿2∗) for all feasible node 𝑖. 𝑐(𝐿1∗) = (𝐿1∗) − 𝜋(𝐿1∗) + 𝜆⋅𝑙𝑛(𝑃 𝑟𝑜𝑏(𝐿1∗)) + 𝐶𝑓 𝑖𝑥𝑒𝑑 (37) = [(𝐿1) + 𝑐(𝐿1)𝑖]−[𝜋(𝐿1) + 𝜋𝑖] +𝜆⋅𝑙𝑛(𝑃 𝑟𝑜𝑏(𝐿1)⋅𝑃𝑟1⨁𝑖 𝑖) + 𝐶𝑓𝑖𝑥𝑒𝑑 (38) = [(𝐿1) − 𝜋(𝐿1) + 𝜆⋅𝑙𝑛(𝑃 𝑟𝑜𝑏(𝐿1))] +𝜆⋅𝑙𝑛(𝑃𝑟1⨁𝑖 𝑖) + 𝑐(𝐿1)𝑖−𝜋𝑖+𝐶𝑓𝑖𝑥𝑒𝑑 (39) ≤[(𝐿2) − 𝜋(𝐿2) + 𝜆⋅𝑙𝑛(𝑃 𝑟𝑜𝑏(𝐿2))] +𝜆⋅𝑙𝑛(𝑃𝑟1⨁𝑖 𝑖) + 𝑐(𝐿1)𝑖−𝜋𝑖+𝐶𝑓𝑖𝑥𝑒𝑑 (40) ≤[(𝐿2) − 𝜋(𝐿2) + 𝜆⋅𝑙𝑛(𝑃 𝑟𝑜𝑏(𝐿2))] +𝜆⋅𝑙𝑛(𝑃𝑟2⨁𝑖 𝑖) + 𝑐(𝐿1)𝑖−𝜋𝑖+𝐶𝑓𝑖𝑥𝑒𝑑 (41) =(𝐿2∗) − 𝜋(𝐿2∗) + 𝜆⋅𝑙𝑛(𝑃 𝑟𝑜𝑏(𝐿2∗)) + 𝐶𝑓𝑖𝑥𝑒𝑑 (42) =𝑐(𝐿2∗)(43) Dual multipliers 𝜆and 𝜋𝑖for 𝑖∈𝑁𝑝are defined in Section 4.2. For 𝑖∉𝑁𝑝, we define 𝜋𝑖= 0. Inequality (40) is derived from condition (5), and inequality (41) is derived from Lemma 1.□ 4.3.4. Heuristic pricing With the exact dominance rules in Proposition 1, we eliminate labels carefully to guarantee that the exact optimal solution can be found, which could be very time-consuming. To speed up the algorithm, we propose heuristic pricing rules which are less restrictive than Proposition 1. As condition (7) in Proposition 1 may be very strict when the travel times change rapidly with time, we loosen it in the heuristic pricing rules. The algorithm may end up with a sub-optimal solution, and the performance of the heuristic pricing is evaluated in Section 5.3. Proposition 2 (Heuristic Pricing).Label 𝐿2is dominated by 𝐿1heuristi- cally, if: (1) (𝐿1) = (𝐿2) (2) (𝐿1)⊆(𝐿2) (3) (𝐿1) = (𝐿2) (4) (𝐿1)≤(𝐿2) (5) (𝐿1) − 𝜋(𝐿1) + 𝜆𝑃 𝑟𝑜𝑏(𝐿1)≤(𝐿2) − 𝜋(𝐿2) + 𝜆𝑃 𝑟𝑜𝑏(𝐿2) (6) 𝐸𝐴(𝐿1)≤𝐸𝐴(𝐿2) (7) 𝑉𝐷(𝐿1)≤𝑉𝐷(𝐿2) (8) 𝑃 𝑟𝑜𝑏(𝐿1)≥𝑃 𝑟𝑜𝑏(𝐿2) 4.4. Branch-cut-and-price framework The branch-cut-and-price algorithm is a branch-and-bound-based method to solve integer programs, which integrates column generation where variables are generated dynamically. Algorithm 1describes the structure of the proposed branch-cut-and- price algorithm. The interested readers are referred to Barnhart et al. (1998) for more details of a branch-and-price algorithm. The column generation and pricing algorithm used that generate new variables are introduced in Sections 4.2 and 4.3 in detail. The following subsections describe the cutting planes and branching strategies used. Algorithm 1 Branch-cut-and-price algorithm while stopping criteria is not satisfied and there exists unprocessed nodes do Select an unprocessed node from the branch-and-bound tree. 𝐹 𝑙𝑎𝑔 ←𝑇 𝑟𝑢𝑒 while True do Solve the linear relaxation of the selected node’s restricted master problem (RMP). Find routes with negative reduced costs by solving the pricing problem repeatedly until no such route can be found. Add the found routes to the RMP. if the current objective value of the RMP is greater than the current upper bound then 𝐹 𝑙𝑎𝑔 ←𝐹 𝑎𝑙𝑠𝑒 break end if Generate cuts using the method introduced in Section 4.4.1. if new cuts are found then Add the constraints to the RMP. else break end if end while if 𝐹 𝑙𝑎𝑔 is 𝐹 𝑎𝑙𝑠𝑒 then continue end if if the solution to the relaxation of the RMP is an integer solution then Update the upper bound. else Generate child nodes using the branching strategies introduced in 4.4.2. Add the generated nodes to the branch-and-bound tree. end if end while 4.4.1. Cutting planes Cuts or valid inequalities strengthen the lower bound after column generation is finished at a given node of the branch-and-bound tree. Jepsen et al. (2008) introduce subset-row cuts (SRC) defined over routing variables 𝑥𝑖, which are: ∑ 𝑟∈⌊𝑝∑ 𝑖∈𝑆 𝑎𝑖𝑟⌋𝑥𝑟≤⌊𝑝|𝑆|⌋ ∀𝑆 ⊆ 𝑁𝑝,0<𝑝<1(44) SRCs are efficient in strengthening the linear relaxation. However, adding the above cut to the RMP makes the pricing problem more challenging, as extra information is needed in the labels to deal with the new constraints. Pecin et al. (2017) propose a generalization of SRCs named limited-memory subset-row cuts (lm-SRC), which have less impact on the pricing problem and make the pricing process faster. We apply the limited-memory subset-row cuts in the following form: ∑ 𝑟∈𝛼(𝐶, 𝑀, 𝑝, 𝑟)𝑥𝑟≤⌊𝑝|𝐶|⌋ ∀𝐶 ⊆ 𝑀 ⊆ 𝑁𝑝,0<𝑝<1(45) where M is an additional memory set, 𝐶 ⊆ 𝑀 ⊆ 𝑁𝑝, and 𝛼is a function of 𝐶, 𝑀, 𝑝, 𝑟. The algorithm to compute coefficient 𝛼and decide memory set 𝑀is explained in Pecin et al. (2017) in detail and is not included here. Even though each label should have an additional dimension for each lm-SRC, the coefficients do not need to be stored in the labels. Therefore, the structure of the pricing problem is not changed, and only minor changes are needed. The interested readers are referred to Pecin et al. (2017) and Sun et al. (2018b). In this paper, we generate cuts with |𝐶|= 3 and 𝑝= 0.5. EURO Journal on Transportation and Logistics 12 (2023) 100099 8 Z. Wang et al. 4.4.2. Branching strategies When the solution to the linear relaxation of the RMP is not integral, we follow the branch-and-bound framework and divide the feasible region of the current node into two regions, generating two child nodes of the current node. The following branching rules are applied. First, we branch on the total number of vehicles used, 𝑛𝑣=∑𝑟∈𝑥𝑟. If 𝑛𝑣is fractional, constraints ∑𝑟∈𝑥𝑟≤⌊𝑛𝑣⌋and ∑𝑟∈𝑥𝑟≥⌈𝑛𝑣⌉are added to the child nodes. The pricing problem will add a dual variable associated with the vehicle number constraint. Then, we branch on the arc-flow variables. The arc (𝑖, 𝑗)whose flow is closest to 0.5 is selected. On the first branch, we set the flow on this arc to 0 by removing arc (𝑖, 𝑗)from arc set 𝐴of this child node. Also, the routes containing arc (𝑖, 𝑗)are removed from the RMP of the child node. On the second branch, we set the flow to 1, and all other arcs starting from 𝑖:(𝑖, 𝑘), 𝑘 ≠𝑗are removed from arc set 𝐴. 5. Computational experiments In this section, we first introduce the design of the computational experiments. Instances for the proposed problem are generated based on the instances introduced by Ropke et al. (2007) for the PDPTW. Then, the performance of the heuristic pricing is evaluated. Finally, we present the results of the experiments and show the proposed method is adequate to satisfy the desired service levels. Besides, the importance of modeling time-dependent travel times is evaluated, and the impact of service levels is analyzed. 5.1. Experiment design Ropke et al. (2007) introduced a set of instances for the PDPTW, where the coordinates of the pickup and delivery nodes are randomly chosen according to a uniform distribution over a [0, 50]×[0, 50] square. The depot is located at (25,25). Our instances adopt the same node coordinates, maximum capacity of vehicles, load quantities, and start times of time windows as those in Ropke et al. (2007). The maximum capacity of a vehicle 𝑄is set to 15. Ropke et al. (2007) assumes the width of time windows is 60 in some instances and 120 in others. We modify Ropke’s assumptions on time window width, making the width of the windows range from 60 to 100. For the fixed cost, we adopt the same parameter as Ropke and Cordeau (2009) that sets a vehicle’s fixed cost 𝐹𝑓𝑖𝑥𝑒𝑑 to 10000. We assume the cost of traversing an arc equals the arc’s length (𝑐𝑖𝑗 =𝑑𝑖𝑗 ). For the parameters regarding stochastic travel and service times, the mean of the service time at each node is randomly chosen from values [8,10,12], and the variance of the service time at each node is randomly chosen from values [1,2,3]. The variance of travel times is decided by variation coefficient 𝑐𝑣, so that 𝑉 𝑎𝑟[𝑋𝑖𝑗 ] = 𝑐𝑣𝑑𝑖𝑗 . In the experiments, 𝑐𝑣 is an adjustable parameter ranging from 0 to 0.4. The means of travel times are decided by the speed profile model (Ichoua et al.,2003), which assumes the travel speed is constant in a time interval on all arcs. We divide the scheduling horizon (0,840) into 4 time intervals: [(0,120),(120,600),(600,720),(720,840)], representing four periods in a day, as suggested in Sun et al. (2018b). Four different speed profiles are used in the experiments, one at a time, presented in Table 1. The mean of the travel times of each arc 𝐸[𝑋𝑖𝑗 (𝑡)] can be computed according to the speed profiles, and 𝐸[𝑋𝑖𝑗 (𝑡)] are piece-wise linear functions. As for service levels, we focus on studying the global service level, rather than node service levels, in our main experiment. This is because some computational experiments studying node service levels can be found in previous work, including Miranda and Conceição (2016) and Gutierrez et al. (2018). The estimation algorithms show that the node service levels can be guaranteed to a large extent. However, no computational study focuses on the global service level. Here, we evaluate how well the global service level can be satisfied using our estimation methods, and we analyze the impact of the different global service levels 𝛩. The desired node service levels are set to 70% in Table 1 Speed profiles. Speed profile number Speed in interval (0, 120) Speed in interval (120,600) Speed in interval (600,720) Speed in interval (720,840) 11111 2 0.9 1.1 0.9 1.1 3 0.75 1.2 0.8 1.25 4 0.67 1.33 0.88 1.33 the experiments. The instance set can be found in the following link: bit.ly/3uX5zF9. We use Monte Carlo simulations to determine the actual success probabilities as benchmarks. We evaluate the algorithm’s effectiveness by checking if the simulated success probabilities satisfy the desired service levels. We also compare the difference between the estimated success probabilities and simulated success probabilities. In each Monte Carlo simulation, travel and service times are sampled from the assumed probability distributions. We compute the simulated success probability after 1000 runs. We compare the proposed algorithm’s routes with those generated by a deterministic routing algorithm that does not consider the stochasticity in travel times and service times (based on Ropke and Cordeau, 2009). The deterministic routing algorithm is introduced in Section 5.2. In this way, the difference in the simulated success rate illustrates the effectiveness of the proposed algorithm. The proposed algorithm is implemented in Java with linear program solver Gurobi 9.0.2. The experiments are performed on a Windows 10 computer with AMD 3.4 GHz CPU and 64 GB RAM. 5.2. An illustrative example In this section, an example is given to show the difference in generated routes between the proposed algorithm and the deterministic routing algorithm. In the example, there are 20 pickup nodes, with the variation coefficient 𝑐𝑣= 0.15, and the width of the time windows is set to 70. The arc travel times follow speed profile number 4 in Table 1. For the deterministic routing algorithm, the travel times are deterministic and computed based on speed profile number 4. Service times are also deterministic, and service time 𝑠𝑖=𝐸[𝑆𝑖]. Service levels are not considered in the deterministic algorithm, and the goal is to minimize the total cost. The proposed algorithm’s desired global service level is set to 70%. The routing solutions generated by the above two algorithms are shown in Fig. 2. Three vehicles are needed in the solution generated by the deterministic routing algorithm, while four vehicles are needed in the solution generated by the proposed routing algorithm. Note that routes #1 in the two generated solutions are identical and are therefore omitted in Fig. 2. Route #2 (blue) and route #3 (red) in the deterministic algorithm are split into three different routes (blue, red, purple) in the proposed algorithm. In the proposed algorithm, one more vehicle is needed to guarantee the desired service rate, leading to a higher total cost. However, compared to the deterministic routing algorithm, the reliability of the routing solution improves significantly, and the simulated global success probability is greater than the desired service level. Some ‘‘compact’’ routes in the deterministic routing algorithm are split into some more ‘‘reliable’’ routes in the proposed algorithm. 5.3. Performance of heuristic pricing As shown in Section 4.3.4, when the travel times are not timedependent, the exact pricing method is identical to the heuristic pricing method, and it is efficient to obtain the exact solution. However, the faster the means of the travel times change throughout the day, the slower the exact pricing algorithm gets. In this section, we evaluate