Vehicle routing with stochastic demand, service and waiting times—The case of food bank collection problems

European Journal of Operational Research 317 (2024) 111–127 Available online 24 March 2024 0377-2217/© 2024 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Contents lists available at ScienceDirect European Journal of Operational Research journal homepage: www.elsevier.com/locate/eor Production, manufacturing, transportation and logistics Vehicle routing with stochastic demand, service and waiting times — The case of food bank collection problems Meike Reusken a,∗, Gilbert Laporte b,c, Sonja U.K. Rohmerb, Frans Cruijssen a aTilburg University, Department of Econometrics and Operations Research, Zero Hunger Lab, Warandelaan 2, 5037 AB Tilburg, The Netherlands bHEC Montréal, 3000 chemin de la Côte-Sainte-Catherine, Montréal H3T 2A7, Canada cSchool of Management, University of Bath, Bath B2A 2AY, United Kingdom ARTICLE INFO Keywords: Distribution Stochastic vehicle routing Food bank supply chain Matheuristic Probabilistic generalized assignment problem ABSTRACT Food banks play an important role both in combating food waste, and in alleviating hunger. However, due to the many uncertainties that food banks face, they often struggle to effectively collect all food items that donors such as supermarkets are willing to provide. To tackle this problem, we introduce the capacitated vehicle routing problem with travel time restrictions and stochastic demand, service and waiting times, in which the uncertainties are dependent of each other. This problem can be generalized to a large variety of routing applications. The goal of the problem is to determine a minimum number of vehicles, and to plan cost-effective routes for these vehicles so that each route violates the vehicle capacity and the travel time limit with only a very small probability. The resulting problem is highly complex and thus solved by means of a matheuristic, which decomposes the problem into its natural decision components. Thus, it first determines the number of districts into which the service area should be partitioned, before allocating each customer to exactly one district and then plans a route for each district. A set of feedback mechanisms is activated whenever no feasible solution has been found through these steps. Extensive numerical experiments, involving both randomly generated and real-life instances, demonstrate the matheuristic’s effectiveness in solving instances with up to 100 customers. When applying our matheuristic to real-life instances from Dutch and Canadian food banks, we furthermore gain managerial insights to assist in optimizing fleet size and route cost. 1. Introduction Today, 17% of all food that reaches the consumer level (at home, in restaurants, or at grocery stores) is wasted (United Nations,2022). At the same time, a rising number of people, currently estimated at 700 million, are food-insecure, which means that they cannot rely on having an adequate meal every day, and they regularly go to bed hungry (FAO et al.,2021). Achieving food security and reducing food waste are integral parts of the United Nations Sustainable Development Goals, as specified in Targets 2.1 and 12.3, respectively. The link between these two goals is apparent, since reducing food waste ensures, in principle, that more resources become available to meet the needs of vulnerable people who are acutely food-insecure. However, the redistribution of resources is often challenging, and food waste prevention at the retail stage can have unintended negative effects on redistribution initiatives at the end of the chains led by organizations such as food banks and soup kitchens. Being largely dependent on donations from food retailers and restaurants, the European association of food banks FEBA signals, ∗Corresponding author. E-mail address: [email protected] (M. Reusken). in this context, a structural decline in donations due to successful food waste prevention efforts at retailers, leaving fewer products to be donated to food banks (European Food Banks Federation,2022). At the same time, the number of active beneficiaries of food banks is rising due to, among other factors, high inflation levels, putting increasing pressure on food bank operations and on the associated redistribution efforts. Given these pressures, it is crucial for food banks to improve their decision making and streamline their operations in order to optimize redistribution activities and fully utilize the available resources. Transport and logistics remain a major challenge in this context, as highlighted in the recent review of Akkerman et al. (2023). One of the main reasons for this is that efficient route planning, which helps make good use of the available capacities, is often hampered by different types of uncertainty which complicate the decision making. In practice, uncertainty affects several aspects of the problem, such as the quantities donated, the waiting times, the loading times, and the availability of drivers. Focusing on operational real-world challenges in humanitarian logistics, such as the role of uncertainty, Besiou et al. https://doi.org/10.1016/j.ejor.2024.03.031 Received 14 August 2023; Accepted 21 March 2024 European Journal of Operational Research 317 (2024) 111–127 112 M. Reusken et al. (2018) illustrate the potential of operations research models in this context. Inspired by these considerations, this research attempts to capture the uncertain decision environment associated with food bank operations and improve the decision making, giving rise to the Capacitated Vehicle Routing Problem (CVRP) with travel time restrictions, stochastic demand, service and waiting times (CVRP-SDSW). The CVRP-SDSW incorporates two types of decisions: strategic choices for designing delivery districts, which cannot be revised regularly, and operational choices regarding the vehicle route associated with each district, which can be adjusted more easily. The goal of the problem is to determine a minimum number of vehicles (equal to the number of districts in our context), and to plan cost-effective routes for these vehicles, so that each customer location is visited and each route is ‘feasible’, i.e., each route violates the vehicle capacity and the travel time limit with only a very small probability. The stochasticity in the demand and in the service and waiting times at the customer locations makes the problem highly complex. The problem is therefore solved by means of a matheuristic, i.e., a ‘‘heuristic algorithm [ ] made by the interoperation of metaheuristics and mathematical programming techniques’’ (Boschetti et al.,2009), which decomposes the problem into its individual decision components to obtain a feasible solution. Thus, it first determines the number of districts, before allocating each customer to exactly one district and then plans a route for each district. The matheuristic enters a set of feedback mechanisms whenever no feasible solution has been found through these steps. Extensive experiments are carried out to test the computational performance of the matheuristic on a set of benchmark instances derived from the well-known instances of Solomon (1987). In addition, we investigate the impact of different levels of uncertainty and gain valuable managerial insights by applying our methodology to several real-life instances using data gathered from three food banks in the Netherlands and Canada. While this paper is motivated by a food bank setting, our problem can be generalized to other contexts and applied to both collection and distribution problems. Examples of this type of setting, featuring these three different types of uncertainty, outside of the food bank context, can be found in the retail industry (e.g., inventoryrouting problems (Federgruen & Simchi-Levi,1995)) as well as problems encountered in the healthcare sector (e.g., blood bank distribution problems (Kaya & Ozkok,2020)). The remainder of this paper is structured as follows. Section 2 provides a literature review and highlights our scientific contributions, before we formally describe the decision problem in Section 3. Section 4 presents the matheuristic that decomposes the problem to find a solution to the CVRP-SDSW. In Section 5, we present results obtained from numerical experiments on random and real-life instances. Finally, we draw some conclusions in Section 6. 2. Literature review The CVRP-SDSW, introduced in this paper, belongs to the class of stochastic Vehicle Routing Problems (VRPs), which have attracted much attention since their introduction by Tillman (1969). In addition, the CVRP-SDSW relates to studies on food bank logistics, an area of research that has gained prominence in the past 10 years (Mahmoudi et al.,2022). The following subsections review the most relevant studies from both of these fields. 2.1. Vehicle routing under uncertainty Considering different types of uncertainty, the CVRP-SDSW constitutes a new variant of the stochastic VRP, which has been extensively studied for decades (see the reviews of Gendreau et al. (2016) and Oyola et al. (2018)). Over time, this extensive research on the stochastic VRP has led to variants that incorporate several aspects of uncertainty. A variant that has received much attention and that closely relates to the CVRP-SDSW, is the capacitated VRP with stochastic demand, for which customer demands only become known after the creation of the routes. In this case, the vehicle capacity may be exceeded, even if the expected demand along the route does not exceed the vehicle capacity. The scientific literature has proposed various actions to repair the solution when this type of route failure occurs. Examples of such actions are the detour-to-depot recourse policy, where the vehicle returns to the depot when capacity failure occurs in order to restore its capacity before resuming its route (Dror et al.,1989;Laporte et al., 2002;Secomandi & Margot,2009), and the preventive return policy, where the route can be interrupted to travel to the depot to reset the capacity prior to the moment of failure, in the hope of avoiding a failure later on (Dror et al.,1993,1989;Yang et al.,2000). Under consideration of a wide range of application contexts, the VRP with stochastic demand has been extended to incorporate a diverse set of practical requirements. Examples of extensions are the consideration of independent compartments in the vehicles (e.g., Goodson,2015; Mendoza et al.,2010), or, as in our study, a restriction in the duration of the routes (e.g., Goodson,2015;Mendoza et al.,2016). In addition to the literature on stochastic demand, several studies have also considered VRPs with stochasticity in the service times (e.g., Errico et al.,2016;Lei et al.,2012). This type of uncertainty is, moreover, often combined with stochastic travel times and time window constraints, while demand is assumed to be deterministic (Li et al.,2010;Miranda & Conceição,2016;Zhang et al.,2013). These studies generally also assume that the stochastic service and travel times are independent of customer demand. In addition to stochastic service and travel times, some papers also account for uncertainty in the waiting time at customer locations. This is most common in VRPs with time windows, where the vehicle needs to wait if it arrives before the customer is ready to begin the service (Li et al.,2010;Zhang et al.,2013). Additionally, some other ways of modeling waiting time have been proposed, e.g., where waiting times are random variables (Keskin et al.,2021) or time-dependent variables (Keskin et al.,2019). Various methodologies have proven to be effective for exactly or approximately solving stochastic VRPs, including robust optimization (Sungur et al.,2008), integer L-shaped methods (Hoogendoorn & Spliet,2023;Laporte & Louveaux,1993; Laporte et al.,2002) and heuristics such as swarm optimization (Marinakis et al.,2013). See Gendreau et al. (2016) and Oyola et al. (2017) for in-depth surveys. The most important distinction between the CVRP-SDSW and the standard stochastic VRP variants is that the CVRP-SDSW simultaneously tackles the uncertainties in multiple routing components. This combination of uncertainty in demand, service and waiting time has, to the best of our knowledge, not previously been studied. Additionally, the service and waiting times in the CVRP-SDSW depend on the size of the demand. This means that, rather than using random variables for these uncertainties as done in other research, we adopt a linear combination of the random variables for the demand. To provide a solution to the CVRP-SDSW, we develop a matheuristic that bears similarities with the stochastic cluster-first route-second methodologies introduced by Yang et al. (2000) and Haugland et al. (2007). These authors propose algorithms for the CVRP where demand is realized before the routing phase (Haugland et al.,2007) and after the routing phase (Yang et al.,2000). Similarly, we also introduce and solve a stochastic variant of the classical two-phase heuristic of Fisher and Jaikumar (1981). The classical first-phase designs the routing districts by selecting seed points and then clustering the customers to these seeds. For this clustering problem under demand uncertainty, heuristics based on tabu search and multistart principles (Haugland et al.,2007) as well as on cost improvement of repositioning, inner-route exchanges and Or-opt principles (Yang et al.,2000) have been presented. In contrast, we formulate this clustering problem as a tractable programming problem, for which we present a linear and convex quadratic variant, European Journal of Operational Research 317 (2024) 111–127 113 M. Reusken et al. which can be solved to optimality for small instances. This programming problem is a probabilistic variant of the generalized assignment problem, which constitutes a significant methodological contribution as, to the best of our knowledge, the probabilistic variant has not yet been investigated. 2.2. Vehicle routing for food banks Decision support models for food aid supply chains have received growing attention in recent years (Mahmoudi et al.,2022;Rivera et al.,2023). The following focuses specifically on food bank logistics, showing that the majority of studies on VRPs for food banks do not take uncertainty into account. Gunes et al. (2010), for example, consider a deterministic pickup and delivery problem for a food rescue program where surplus food is collected from suppliers and distributed to food agencies. The combination of a deterministic VRP and resource allocation decisions has, moreover, received considerable attention in recent years (Eisenhandler & Tzur,2019a,2019b;Nair et al.,2018,2016,2017;Orgut & Lodree, 2023;Rey et al.,2018). For instance, Eisenhandler and Tzur (2019a) and Eisenhandler and Tzur (2019b) model the effective and equitable allocation of food alongside its delivery to food agencies. Another deterministic problem integrating routing with facility location decisions was studied by Solak et al. (2014) and Davis et al. (2014). Solak et al. (2014) introduced the ‘‘stop-and-drop’’ problem, in which food is delivered by the food bank to ‘‘drop’’ locations, from which agencies can in turn collect it. Their problem particularly applies to food banks in remote areas and determines the routes for the food banks, the optimal placement of these drop locations, and the matching of the agencies to the drop locations. Davis et al. (2014) study a similar problem by using a set of satellite locations as possible transshipment points between a food bank and an agency. Despite this growing interest in VRPs for food banks, there is a notable lack of attention towards the modeling of uncertainty in this context. To the best of our knowledge, Lien et al. (2014) and Balcik et al. (2014) are the only authors who have modeled uncertainty. Lien et al. (2014) include uncertainty in the demand, considering a problem in which the order of visits to the donors and the food agencies, as well as the delivery amounts to the agencies, must be determined for a single vehicle. The demand at each node remains unknown until the arrival of the vehicle. Balcik et al. (2014) extended this problem to a multi-vehicle setting. Table 1 provides an overview of the presented food bank papers, summarizing the types of uncertainty that have been included (column 2) and the country from which the problem or the case study originates (column 3). The table reveals that studies on VRPs for food banks taking into account uncertainty are still rare, while studies focusing on other aspects in the food bank supply chain increasingly highlight the existence and importance of uncertainty (Akkerman et al.,2023). As a result, research on food bank supply chains has started to include uncertainties, for example, in the demand for food assistance (Reusken et al.,2023), the supply of food donations (Davis et al.,2016;Paul & Davis,2022), as well as the available capacity at food agencies (Orgut et al.,2018). Despite this increasing attention, the survey of Mahmoudi et al. (2022) confirms the observed gap, and stresses the significance of considering the multitude of uncertain factors that characterize reallife food bank settings. One goal of this research is therefore to fill this research gap by incorporating the notion of uncertainty into the context of VRPs for food banks. The inspiration for the problem, as well as the real-life data used for this research, originate from food banks in the Netherlands and Canada, two settings that have not yet been studied (see Table 1). Table 1 Table summarizing the literature on VRPs for food banks. Uncertainty Country Gunes et al. (2010) USA Balcik et al. (2014) Demand USA Davis et al. (2014) USA Lien et al. (2014) Demand USA Solak et al. (2014) USA Nair et al. (2016) Australia Nair et al. (2017) Australia Nair et al. (2018) Australia Rey et al. (2018) Australia Eisenhandler and Tzur (2019a) Israel and USA Eisenhandler and Tzur (2019b) Israel and USA Orgut and Lodree (2023) USA 3. Problem description While this paper is inspired by the practical decisions faced by food banks, our problem generalizes the single-period capacitated VRP with travel time restrictions, and stochasticity in the demand, as well as in the service and waiting times at the customer locations. Incorporating both strategic and operational decisions, the overall goal of this problem is to determine a minimum number of districts (equal to the number of vehicles in our context) as well as to identify a costeffective route for each of these districts, serving all customers. Since the strategic design of the districts significantly impacts operational decisions related to the routes, it is important to ensure that the proposed decisions are robust. We therefore prioritize this strategic aspect of the problem, generating solutions that withstand the uncertainties with a controlled probability. Based on these decisions at the strategic level, the problem then aims to generate routes that cope with uncertainty at the operational level through the use of a recourse policy. Considering the respective significance of these two decision levels, the objective of the problem is hierarchical in nature: first minimizing the number of vehicles and then the routing cost, which in our case is expressed in the form of travel time. One of the main challenges of this problem lies in the incorporation of three sources of uncertainty. First, the demand remains unknown until the vehicle arrives at a customer location. Second, the service time is also stochastic since it is considered proportional to the demand and includes both handling activities at the customer locations as well as at the depot. Third, since vehicles may be forced to queue before servicing a customer in case no loading dock is immediately available, this waiting time is uncertain and depends on the number of vehicles in the queue. The length of the queue will only be known upon arrival, and the unit service time for the vehicles waiting in the queue is assumed to be identical for all vehicles at a given location. These time considerations are accounted for in the travel time, which encompasses the time spent on service, waiting, driving and potentially returning to the depot in order to restore capacity in case of route failure. The challenge of incorporating these uncertainty sources is further amplified by the different decision levels involved in our problem, i.e., the combination of strategic and operational decisions. A more detailed overview of these individual decisions is provided in the following. 1. Strategic decisions. The decisions at the strategic level all relate to the design of districts. The first decision at this level concerns the number of districts into which the service area should be divided. For each of the districts, the decision maker requires a vehicle and a driver, and all vehicles are considered to be identical. Since the team of drivers and the fleet of vehicles cannot be regularly revised, the number of districts is a decision made at the strategic level that holds highest priority. Furthermore, we consider the fixed cost associated with each driver as well as each vehicle to be substantially large, making it reasonable to assume that fewer vehicles consistently leads to lower costs. European Journal of Operational Research 317 (2024) 111–127 114 M. Reusken et al. As a result, the ideal number of districts must be as small as possible. Given the number of districts, the second strategic decision concerns the assignment of each customer to precisely one district, with the aim of minimizing the overall route cost. Both of these decisions are dependent on (i) the customer demands to satisfy the vehicle capacity, and (ii) the total travel time which should normally not exceed the maximum duration of the driver’s working day. Hence, the number of districts as well as the customer assignment to districts must be robust to this stochastic setting in which the probability of not satisfying any of these two conditions may not exceed a predetermined threshold. This guarantees that violations of the vehicle capacity and the travel time limit are infrequent occurrences. 2. Operational decisions. Compared to the strategic decisions, the decisions made at the operational level concern the more immediate and shorter-term operations planning, which can be revised more flexibly. The final decision in our problem lies at the operational level and relates to the route associated with each district, i.e., the order in which customers are visited. Each route starts and ends at the depot, which is shared among the districts. Given the demand uncertainty, higher levels than the anticipated demand can result in a route that violates the vehicle capacity, i.e., in a route failure. As in Dror et al. (1993), we assume that this can only occur at the last customer of the route. This is justified by the fact that the probability of exceeding the vehicle capacity twice on the same route is extremely small, and that if at most one failure can occur, the customer location with the highest probability of failure is the last one. One of two actions can be taken in this situation: (i) return to the depot upon failure to reset the capacity, or (ii) plan a preventive return to the depot at a suitable location earlier in the route to avoid a possibly costly failure at the last customer location (Dror et al., 1989;Dror et al.,1993). 4. Methodology Given its stochastic nature and the fact that it combines several NPhard subproblems, the CVRP-SDSW is highly complex and difficult to solve exactly, even for relatively small instances. Therefore, we propose a matheuristic that decomposes the problem into several natural decisions in order to find a solution. In this context, the algorithm first determines the number of districts into which the service area should be divided. In a second step, it defines these districts, using strategically selected seeds and solving a probabilistic generalized assignment problem to cluster the customers. After the creation of the districts, the matheuristic then proceeds with the route planning for each district. Periodically, the matheuristic enters an iterative procedure (IP) to verify or restore the feasibility of the solution. An overview of the structure of the algorithm is depicted in Fig. 1. The following subsections describe each step in detail. The problem is defined over a complete directed graph 𝐺(𝑉 , 𝐴), where the vertices 𝑉= {0} ∪ 𝑁represent all locations to be visited, comprising the depot location {0} and the customer locations 𝑁= {1,…, 𝑛}, while 𝐴= {(𝑖, 𝑗) ∶ 𝑖, 𝑗 ∈𝑉 , 𝑖 ≠𝑗}denotes the set of arcs. Table 2 provides a summary of the notation used. The remaining notations in this table are explained in more detail in the subsections. 4.1. Determining the number of districts The first step of the matheuristic is to select the initial number of districts, denoted by 𝑚. This also determines the number of vehicles and drivers, since there is exactly one vehicle per district. Given the uncertainties involved in the problem, there is a possibility that route failure may occur as a result of exceeding the vehicle capacity or the maximum travel time. As such, the choice of 𝑚should be sufficiently robust against this uncertainty inherent to the problem. We enforce Fig. 1. Flowchart of the structure of our matheuristic. IP1, IP2and IP3are three iterative procedures to be described in Section 4.4. robustness by using chance constraint programming, ensuring that the probability of route failure related to the choice of 𝑚does not exceed given thresholds. For this purpose, we apply two chance constraints in order to respect the vehicle capacity and the travel time limit, determining the corresponding number of districts 𝑚𝑐and 𝑚𝑡, respectively. Considering the values of 𝑚𝑐and 𝑚𝑡, the initial number of districts 𝑚is then equal to max{𝑚𝑐, 𝑚𝑡}. In terms of the vehicle capacity, this means that 𝑚𝑐should be set in such a way that the sum of the stochastic demands exceeds the capacity 𝑄not more than a given proportion 𝛾of the time i.e., 𝑚𝑐is the smallest integer satisfying P(∑ 𝑖∈𝑁 𝜉𝑖> 𝑚𝑐𝑄)≤𝛾, (1) where 𝜉𝑖represents the stochastic demand of customer 𝑖∈𝑁. In addition, uncertainty may also affect the total travel time of a route, which may result in a route failure if the travel time exceeds a limit 𝑇. We consider four individual components which together determine the travel time: 1. Driving. The total duration of driving is denoted by 𝑑𝑣. 2. Service. The service time at a customer 𝑖∈𝑁is defined as 2𝑎𝜉𝑖, where 𝜉𝑖is the demand and 𝑎is the service time per demand unit. The demand is counted twice, as vehicles need to be both loaded and unloaded. Hence, the service time includes handling at the customer location as well as at the depot location. 3. Waiting. The waiting time of a vehicle at a customer’s location 𝑖∈𝑁is considered to be dependent on the number of vehicles already in the queue, and the probability that there are ℎvehicles in the queue at the customer is denoted by 𝑞𝑖ℎ. The service time for all vehicles in the queue is assumed to be identical and proportional to 𝜉𝑖, except for the first vehicle in the queue, which is expected to already be halfway through the service, so that the total service time of the ℎvehicles in the queue at customer 𝑖is 𝑎(1 2+ (ℎ− 1))𝜉𝑖. The expected waiting time at customer 𝑖is then equal to ∑∞ ℎ=1 𝑞𝑖ℎ[𝑎(ℎ−1 2)𝜉𝑖]. The need for 𝜉𝑖arises from the lack of information about the time spent by the other vehicles in the queue, and therefore 𝜉𝑖is used as a proxy for the demand of the other vehicles assuming that this demand is similar in size across all vehicles in the queue. This is justified by the idea that large customers will for example have overall a higher level of activity. 4. Recourse. The final component that influences travel time arises from the recourse policy, which is applied whenever the realized European Journal of Operational Research 317 (2024) 111–127 115 M. Reusken et al. Table 2 Notations. We use lowercase boldface characters to denote vectors. Sets 𝑁Customers, i.e., 𝑁= {1,…, 𝑛} 𝑉Vertices, i.e., 𝑉= {0} ∪ 𝑁 𝐴Arcs, i.e., 𝐴= {(𝑖, 𝑗) ∶ 𝑖, 𝑗 ∈𝑉 , 𝑖 ≠𝑗} 𝐾Districts, i.e., 𝐾= {1,…, 𝑚} 𝑆𝑘Customers in district 𝑘∈𝐾, i.e., 𝑆𝑘= {𝑖∈𝑁∶𝑥𝑖𝑘 = 1} 𝐷𝑘Vertices in district 𝑘∈𝐾used for routing, i.e., 𝐷𝑘= {0} ∪ 𝑆𝑘∪ {|𝑆𝑘|+ 1,|𝑆𝑘|+ 2} 𝐾∗ 𝜆Districts for which the assignment of customers to districts is infeasible at iteration 𝜆, where 𝐾∗ 𝜆⊆ 𝐾 𝑆∗ 𝜆Accumulated family of sets of customer groups that result in infeasible solutions at iteration 𝜆, i.e., 𝑆∗ 𝜆= {𝑆∗ 𝜆−1 ∪𝑆𝑘∶𝑘∈𝐾∗ 𝜆} Parameters 𝛾Safety factor for vehicle capacity to determine 𝑚𝑐 𝛿Safety factor for travel time limit to determine 𝑚𝑡 𝛼Safety factor for vehicle capacity in assignment of customers to districts 𝜂Safety factor for travel time limit in assignment of customers to districts 𝑎Service time per demand unit 𝑠Driving speed 𝑞𝑖ℎ Probability that ℎvehicles are waiting in queue at customer 𝑖∈𝑁 𝑐𝑖𝑗 Length of arc (𝑖, 𝑗) ∈ 𝐴 𝑀Arbitrarily large positive number 𝑄Vehicle capacity 𝑇Travel time limit 𝜆Iteration count 𝜓Termination scalar for the iteration count Unknown parameters 𝑑𝑣 𝑘Duration of driving in district 𝑘∈𝐾 𝑑𝑣Total duration of driving 𝑑𝑟 𝑘Duration of recourse in district 𝑘∈𝐾 𝑑𝑟Total duration of recourse 𝑐𝑖𝑘 Distance associated with assigning customer 𝑖∈𝑁to district 𝑘∈𝐾 Approximations 𝑐𝑖𝑘 Approximated distance of assigning customer 𝑖∈𝑁to district 𝑘∈𝐾  𝑑1Approximation for 𝑑𝑣+𝑑𝑟used to determine the number of districts  𝑑2(𝒙)Approximation for 𝑑𝑣+𝑑𝑟used in the assignment of customers to districts Variables 𝜉𝑖Random variable for stochastic demand of customer 𝑖∈𝑁 𝑥𝑖𝑘 Decision variable for assignment of customers to districts: 1 if customer 𝑖∈𝑁is allocated to district 𝑘∈𝐾, 0 otherwise 𝑦𝑖𝑗𝑘 Decision variable for routing: 1 if arc (𝑖, 𝑗) ∈ 𝐷𝑘is used in district 𝑘∈𝐾, 0 otherwise demand exceeds the vehicle capacity at the last customer. To overcome this type of failure and repair the solution, we consider two possible actions. First, the vehicle may travel from the last customer to the depot to reset its capacity, before revisiting that customer. In this case, the duration of recourse 𝑑𝑟is equal to the expected additional duration of the round trip between the last customer and the depot. As a second option, vehicles can plan a preventive return to the depot at a suitable location (see Section 3), in which case the duration of recourse 𝑑𝑟is equal to zero, and the duration of driving 𝑑𝑣includes the driving time of the extra stop at the depot. These four elements aggregate into the total travel time, which should not exceed the maximum time limit 𝑇with a probability of more than 𝛿. Based on this, we formulate the following chance constraint in order to determine the appropriate number of districts, equal to the smallest integer 𝑚𝑡satisfying P(𝑑𝑣+𝑑𝑟+∑ 𝑖∈𝑁(2𝑎𝜉𝑖+ ∞ ∑ ℎ=1 𝑎(ℎ−1 2)𝑞𝑖ℎ𝜉𝑖)> 𝑚𝑡𝑇)≤𝛿. (2) It should be noted, in this context, that to be able to compute 𝑚𝑡from this constraint, the demands (𝜉𝑖, 𝑖 ∈𝑁) need to be independent and follow a stable distribution (e.g., a normal or a Cauchy distribution, see Fama and Roll (1968)). In this research, we present tractable solutions and experiments for normal independently distributed demand, which is the most frequently considered demand distribution in the stochastic VRP literature (Oyola et al.,2018). From a food bank perspective, the assumption of independent demands is sensible since donors may vary significantly (e.g. supermarkets, farmers, big manufacturers), making it difficult to identify common factors affecting the volume of food donations. Since the actual time spent on driving and recourse are not yet known at this stage of the solution procedure, we introduce  𝑑1as an approximation for 𝑑𝑣+𝑑𝑟. Approximation 1. [ 𝑑1]Due to the limited information that is available regarding the districts at this stage, we use the formula of Beardwood et al. (1959) to approximate the length of an optimal route of a Traveling Salesman Problem (TSP) on the graph 𝐺(𝑉 , 𝐴). Considering that 𝑡𝑖𝑚𝑒 = 𝑑𝑖𝑠𝑡𝑎𝑛𝑐𝑒∕𝑠𝑝𝑒𝑒𝑑, the approximation for the duration is given by  𝑑1=𝛽(𝑛)√𝐸(𝑛+ 1) 𝑠, where 𝑠denotes the speed and 𝐸is the size of the service area, estimating the constant 𝛽(𝑛), as in Tables 1and 2of Franceschetti et al. (2017). 4.2. Districting Now that the number of districts has been determined, we assign each customer to exactly one district. This is achieved in two stages. First, a subset of customers are chosen to serve as seed points, one for each district. In the second stage, a probabilistic generalized assignment problem is solved, where, as in Fisher and Jaikumar (1981), the seeds are used to estimate the distance associated with assigning a customer to a specific district. This problem is stochastic, since demand, service and waiting times are uncertain. Solving these two stages, which are described in more detail below, provides a feasible assignment of customers for each district (or vehicle), considering the limited available information at this stage. European Journal of Operational Research 317 (2024) 111–127 116 M. Reusken et al. 4.2.1. Seed selection While several methods have been proposed for seed selection (e.g., Baker & Sheasby,1999;Koskosidis & Powell,1992;Sultana et al., 2017), this research applies the k-means clustering algorithm, initially proposed by Macqueen (1967). In k-means clustering, several data points are partitioned into clusters so that the within-cluster variances (squared Euclidean distances) are minimized. The method is summarized in Algorithm 1 and adapted to partition the 𝑛customers into 𝑚 clusters. We adopt the common practice of repeating the entire k-means algorithm multiple times with different initial centroids and choose the clustering result that yields the lowest sum of squared distances. For this clustering result, we select for each cluster the furthest customer from the depot as the seed, following the seed selection method of Fisher and Jaikumar (1981). This composite method works well for all types of customer densities, and is also intuitive and computationally efficient. Algorithm 1 k-means clustering 1: Randomly select 𝑚customers as the initial centroids 2: repeat 3: Assign each customer to its closest centroid to create clusters 4: Update centroids by taking means of customer locations within each cluster 5: until 6: The centroids remain unchanged 4.2.2. Assignment of customers to districts Next, we allocate the customers to districts using a probabilistic extension of the classical Generalized Assignment Problem proposed by Fisher and Jaikumar (1981) for the deterministic case. We introduce the binary decision variables 𝑥𝑖𝑘 for the customer assignment to the districts, where 𝑥𝑖𝑘 = 1 if and only if customer 𝑖∈𝑁is allocated to district 𝑘∈𝐾. We formulate the probabilistic generalized assignment problem as follows: min ∑ 𝑘∈𝐾∑ 𝑖∈𝑁 𝑐𝑖𝑘𝑥𝑖𝑘 (3) s.t. ∑ 𝑘∈𝐾 𝑥𝑖𝑘 = 1 𝑖∈𝑁(4) P(∑ 𝑖∈𝑁 𝜉𝑖𝑥𝑖𝑘 > 𝑄)≤𝛼 𝑘 ∈𝐾(5) P(𝑑𝑣 𝑘+𝑑𝑟 𝑘+∑ 𝑖∈𝑁 𝑥𝑖𝑘 (2𝑎𝜉𝑖+ ∞ ∑ ℎ=1 𝑎(ℎ−1 2)𝑞𝑖ℎ𝜉𝑖)> 𝑇 )≤𝜂 𝑘 ∈𝐾(6) 𝑥𝑖𝑘 ∈ {0,1} 𝑖∈𝑁, 𝑘 ∈𝐾. (7) The objective (3) of this model is to minimize the total travel distance, where 𝑐𝑖𝑘 denotes the distance associated with assigning customer 𝑖∈𝑁to district 𝑘∈𝐾. Constraints (4) then ensure that each customer is allocated to exactly one district, while chance constraints (5) guarantee that the probability of exceeding the vehicle capacity in a district is at most 𝛼. Similarly, constraints (6) ensure that the total travel time, consisting of the driving, recourse, service and waiting time in a district, exceeds the travel time limit with a probability of at most 𝜂. Finally, constraints (7) define the domains of variables 𝑥𝑖𝑘. Note that in this formulation the parameters 𝑐𝑖𝑘,𝑑𝑣 𝑘and 𝑑𝑟 𝑘are unknown. In order to estimate them we introduce an updated approximation for 𝑑𝑣+𝑑𝑟. Approximation 2. [ 𝑑2(𝒙)] At this stage of the solution procedure, new information relative to Approximation 1 becomes available, namely the number of districts as well as the seed of each district. Consequently, we propose the updated approximation  𝑑2(𝒙) = ∑ 𝑘∈𝐾(𝑐0𝑗𝑘+𝑐𝑗𝑘0 𝑠+∑ 𝑖∈𝑁 𝑥𝑖𝑘 𝑐𝑖𝑘 𝑠), where 𝒙is the vector of variables 𝑥𝑖𝑘, 𝑖 ∈𝑁, 𝑘 ∈𝐾and 𝑐𝑖𝑘 denotes the approximate distance of assigning customer 𝑖∈𝑁to district 𝑘∈𝐾. This distance corresponds to the increase in distance associated with inserting customer 𝑖∈𝑁in a route that starts and ends at the depot and visits the seed in district 𝑘∈𝐾. Since 𝑐𝑖𝑘 includes only the increase in distance of inserting a customer, the round trip between the seed and the depot is separately added in  𝑑2(𝒙)using 𝑐0𝑗𝑘, which denotes the distance from the depot to the seed customer 𝑗𝑘of district 𝑘∈𝐾. Using this estimate the objective in (3) can be approximated using min  𝑑2(𝒙),(8) which provides a solution that is equivalent to minimizing (i) the approximate total travel distance (i.e., ∑𝑘∈𝐾∑𝑖∈𝑁𝑐𝑖𝑘𝑥𝑖𝑘) and (ii) the approximate total travel time (i.e., adding the expected service and waiting time to  𝑑2(𝒙)). These formulations are identical in terms of the optimal assignment since all components, except 𝑐𝑖𝑘𝑥𝑖𝑘, are constants that are independent of the assignment. In addition to the objective, constraints (6) can be approximated using the expression  𝑑2(𝒙), which gives for 𝑘∈𝐾 P(𝑐0𝑗𝑘+𝑐𝑗𝑘0 𝑠+∑ 𝑖∈𝑁 𝑥𝑖𝑘 (𝑐𝑖𝑘 𝑠+ 2𝑎𝜉𝑖+ ∞ ∑ ℎ=1 𝑎(ℎ−1 2)𝑞𝑖ℎ𝜉𝑖)> 𝑇 )≤𝜂. (9) In summary, the approximated model, consisting of (4),(5),(7), (8) and (9), establishes the assignment of customers to districts. In order to solve this model, we introduce a tractable reformulation of the probabilistic constraints (5) and (9). In the following, a convex quadratic and a linear reformulation are presented. One can implement the assignment problem by choosing one of these reformulations, or alternatively by using another reformulation that is computationally tractable. Convex quadratic reformulation Since the random variables 𝜉𝑖, 𝑖 ∈𝑁are assumed to be mutually independent and to follow normal distributions with means 𝜇𝑖, 𝑖 ∈ 𝑁and variances 𝜎2 𝑖, 𝑖 ∈𝑁, we can reformulate constraints (5) by exploiting the property that normal distributions are invariant under addition (Fama & Roll,1968). Hence, ∑𝑖∈𝑁𝜉𝑖𝑥𝑖𝑘 follows the same distribution as 𝜉𝑖, with mean ∑𝑖∈𝑁𝜇𝑖𝑥𝑖𝑘 and variance ∑𝑖∈𝑁𝜎2 𝑖𝑥2 𝑖𝑘. Constraints (5), therefore, allow the following reformulation ∑ 𝑖∈𝑁 𝜇𝑖𝑥𝑖𝑘 +𝐹−1(1 − 𝛼)√∑ 𝑖∈𝑁 𝜎2 𝑖𝑥2 𝑖𝑘 ≤𝑄 𝑘 ∈𝐾, (10) where 𝐹−1(⋅)denotes the inverse cumulative distribution function of 𝜉𝑖(also frequently named the quantile or percent-point function). Inequalities (10) hold if the following constraints are satisfied for 𝑘∈ 𝐾: (𝐹−1(1 − 𝛼))2∑ 𝑖∈𝑁 𝜎2 𝑖𝑥𝑖𝑘 + 2𝑄∑ 𝑖∈𝑁 𝜇𝑖𝑥𝑖𝑘 −(∑ 𝑖∈𝑁 𝜇𝑖𝑥𝑖𝑘)2 ≤𝑄2(11) 𝑄−∑ 𝑖∈𝑁 𝜇𝑖𝑥𝑖𝑘 ≥0.(12) Note that these inequalities exploit the property that 𝑥𝑖𝑘 is binary, i.e., that 𝑥2 𝑖𝑘 =𝑥𝑖𝑘, 𝑖 ∈𝑁, 𝑘 ∈𝐾. Inequalities (11) then follow directly from squaring both sides of constraints (10). This operation is paired with the requirements of 𝑄−∑𝑖∈𝑁𝜇𝑖𝑥𝑖𝑘 and 𝐹−1(1 − 𝛼)√∑𝑖∈𝑁𝜎2 𝑖𝑥2 𝑖𝑘 being nonnegative. The first requirement is enforced by constraints (12), and the second requirement holds trivially, since the components themselves are nonnegative. A tractable convex quadratic reformulation of constraints (9) can be derived equivalently (see Appendix A). European Journal of Operational Research 317 (2024) 111–127 117 M. Reusken et al. Linear reformulation The nonlinearity arising in the term (∑𝑖∈𝑁𝜇𝑖𝑥𝑖𝑘)2of constraints (11) can, moreover, also be linearized. First, note that ∑ 𝑖∈𝑁 𝜇2 𝑖𝑥2 𝑖𝑘 + 2 ∑ 𝑖∈𝑁∑ 𝑗≠𝑖,𝑗∈𝑁 𝜇𝑖𝜇𝑗𝑥𝑖𝑘𝑥𝑗𝑘 =∑ 𝑖∈𝑁 𝜇2 𝑖𝑥𝑖𝑘 + 2 ∑ 𝑖∈𝑁∑ 𝑗≠𝑖,𝑗∈𝑁 𝜇𝑖𝜇𝑗𝑥𝑖𝑘𝑥𝑗𝑘, since 𝑥𝑖𝑘 is binary. By introducing the new variables 𝑧𝑖𝑗𝑘, 𝑖 ∈𝑁, 𝑗 > 𝑖, 𝑗 ∈𝑁, 𝑘 ∈𝐾, inequalities (11) can be replaced with the following linear constraints in 𝑥𝑖𝑘 and 𝑧𝑖𝑗𝑘: (𝐹−1(1 − 𝛼))2∑ 𝑖∈𝑁 𝜎2 𝑖𝑥𝑖𝑘 + 2𝑄∑ 𝑖∈𝑁 𝜇𝑖𝑥𝑖𝑘 −∑ 𝑖∈𝑁 𝜇2 𝑖𝑥𝑖𝑘 − 2 ∑ 𝑖∈𝑁∑ 𝑗≠𝑖,𝑗∈𝑁 𝜇𝑖𝜇𝑗𝑧𝑖𝑗𝑘 ≤𝑄2𝑘∈𝐾 𝑧𝑖𝑗𝑘 ≥𝑥𝑖𝑘 +𝑥𝑗𝑘 − 1 𝑖∈𝑁, 𝑗 > 𝑖, 𝑗 ∈𝑁, 𝑘 ∈𝐾 𝑧𝑖𝑗𝑘 ≤𝑥𝑖𝑘 𝑖∈𝑁, 𝑗 > 𝑖, 𝑗 ∈𝑁, 𝑘 ∈𝐾 𝑧𝑖𝑗𝑘 ≤𝑥𝑗𝑘 𝑖∈𝑁, 𝑗 > 𝑖, 𝑗 ∈𝑁, 𝑘 ∈𝐾. 4.3. Routing After determining the assignment of customers to districts, we need to solve a stochastic TSP for each district. Since we consider a setting where (i) route failure may only occur at the last customer’s location and (ii) a preventive return to the depot may be scheduled in the hope of avoiding a costlier return at a later stage, the transformation of Dror et al. (1993) can be applied. This transformation reduces our stochastic TSP to a deterministic TSP, by building on the following ideas. First, only one of the two cases (i) and (ii) can occur. Second, the probability of route failure in district 𝑘is equal to 𝑝𝑘=P(∑ 𝑖∈𝑆𝑘 𝜉𝑖> 𝑄)𝑘∈𝐾, where 𝑆𝑘denotes the customers in district 𝑘, i.e., 𝑆𝑘= {𝑖∈𝑁∶ 𝑥𝑖𝑘 = 1}. Note that 𝑝𝑘≤𝛼for all 𝑘. Based on this setting, two artificial vertices |𝑆𝑘|+1 and |𝑆𝑘|+ 2 are introduced for each district 𝑘. Entering |𝑆𝑘|+1 from customer 𝑖indicates a preventive route break, and entering |𝑆𝑘|+2 from customer 𝑖indicates that route failure may occur at a cost of 𝑝𝑘(𝑐𝑖0+𝑐0𝑖). One of these two vertices is visited, which is imposed by setting 𝑐|𝑆𝑘|+1,|𝑆𝑘|+2 and 𝑐|𝑆𝑘|+2,|𝑆𝑘|+1,𝑘∈𝐾, equal to −𝑀, where 𝑀 is an arbitrarily large positive number. For an explanation of the other 𝑐𝑖𝑗 associated with the artificial vertices and a detailed description of the transformation, we refer to Dror et al. (1993), pp. 280–281. Considering the original vertices of a district together with the newly introduced artificial ones, we obtain the sets 𝐷𝑘= {0} ∪ 𝑆𝑘∪ {|𝑆𝑘|+ 1,|𝑆𝑘|+ 2}, 𝑘 ∈𝐾. Using the Dror et al. (1993) transformation of the distance matrix, the routes can now be determined by solving 𝑚 TSPs, i.e., one for every vertex set 𝐷𝑘. By determining these routes, we obtain the expected duration of driving and recourse, i.e., the routing solutions for 𝑘∈𝐾are 𝑑𝑟 𝑘=⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 𝑝𝑘(𝑐𝑖0+𝑐0𝑖) 𝑠if route failure is risked at the last visited customer 𝑖in district 𝑘 0if a preventive return to the depot is scheduled in district 𝑘 (13) 𝑑𝑣 𝑘=∑(𝑖,𝑗)∈𝐷𝑘𝑐𝑖𝑗 𝑦𝑖𝑗𝑘+𝑀 𝑠,(14) where 𝑦𝑖𝑗𝑘 is a binary variable equal to 1 if and only if arc (𝑖, 𝑗) ∈ 𝐷𝑘is used in the TSP solution in district 𝑘. Note that the travel time of the extra trip to the depot, resulting from the recourse action of scheduling a preventive return to the depot, is included in the driving time 𝑑𝑣 𝑘. The expected duration of driving and recourse of the combined districts is then equal to 𝑑𝑟=∑ 𝑘∈𝐾 𝑑𝑟 𝑘(15) 𝑑𝑣=∑ 𝑘∈𝐾 𝑑𝑣 𝑘.(16) A feasible solution to the CVRP-SDSW is obtained when the time chance constraints used to determine the number of districts (constraint (2)) and the assignment of customers to districts (constraints (6)) are satisfied for these routing solutions. 4.4. Iterative procedures to achieve feasibility The last feature of the matheuristic lies in the iterative procedures IP1, IP2and IP3, shown in Fig. 1. These procedures aim to regain feasibility when an infeasible solution is encountered. It may happen, however, that even after several iterations no feasible solution can be identified. This can be for one of two reasons: (i) the approximation  𝑑2(𝒙)used in the assignment of customers to districts is inaccurate for a certain instance and prevents the matheuristic from finding a feasible solution, or (ii) the uncertainty inherent to the problem yields infeasible solutions. In the case of these rare events, two stopping criteria are considered for the algorithm. First, the number of iterations is limited to at most 𝜓|𝐾|, where 𝜓should be large enough to evaluate an adequate number of solutions. This limit depends on |𝐾|, since larger instances are likely to require more iterations in order to reach feasibility. A second stopping criterion is a runtime limit, which is checked before proceeding to another iteration. As long as these limits are not reached, the iterative procedures continue the search for a feasible solution. More precisely, these procedures activate the following actions. IP1:Increase 𝑚when the assignment of customers to districts is infeasible. The number of districts is increased to 𝑚+ 1, when the assignment of customers to districts is infeasible. IP2:Decrease 𝑚when it has been incorrectly initialized. The number of districts is decreased to 𝑚− 1 when the approximation  𝑑1resulted in initializing too many districts. Specifically, we test the feasibility of constraint (2), while adopting the routing solution 𝑑𝑣+𝑑𝑟(obtained from (15) and (16)). If this results in a reduction in the value of 𝑚, the iterative procedure IP2returns to the districting subproblem with 𝑚− 1. Given that IP1yields the opposite output, i.e., iterating while considering 𝑚+ 1, the tested values of 𝑚are stored in order to consider each value at most once. Note that iterating for this reason will rarely happen, since the following scenarios would have to occur simultaneously: (i) the travel time limit is more constraining than the vehicle capacity, i.e., 𝑚𝑐≤𝑚𝑡, and (ii) the adjustment in 𝑚𝑡results in a lower integer value. IP3:Exclusion of infeasible customer combinations when the time chance constraints (6) are violated. Infeasible customer combinations are excluded from forming districts in subsequent iterations of the assignment problem when the time chance constraints (6) are violated while considering the routing solutions 𝑑𝑣 𝑘+ 𝑑𝑟 𝑘, 𝑘 ∈𝐾obtained from (13) and (14). If a district is in violation at iteration 𝜆, it is then added to the set of infeasible districts 𝐾∗ 𝜆, where 𝐾∗ 𝜆⊆ 𝐾. Across iterations, we keep track of the customers belonging to these infeasible districts, using the accumulated family of sets at iteration 𝜆, denoted by 𝑆∗ 𝜆={𝑆∗ 𝜆−1∪𝑆𝑘∶𝑘∈𝐾∗ 𝜆}, where 𝑆𝑘= {𝑖∈𝑁∶𝑥𝑖𝑘 = 1} denotes the customers in district 𝑘∈𝐾. As such, the assignment problem (presented in Section 4.2.2) at iteration 𝜆+ 1 is augmented with the following constraints: ∑ 𝑖∈𝐵 𝑥𝑖𝑘 ≤|𝐵|− 1 𝑘∈𝐾, 𝐵 ∈𝑆∗ 𝜆.(17) European Journal of Operational Research 317 (2024) 111–127 118 M. Reusken et al. Fig. 2. Depot and customer locations in instances c101, c201 and r101. The square depicts the depot, whereas the points correspond to customers, with the lighter colored points corresponding to the first 20 customers. 5. Numerical analysis The proposed matheuristic was coded in Python 3.10.5 using scikit-learn 1.0.2 for seed selection, Gurobi 9.5.2 for the assignment of customers to districts, and python-tsp 0.3.1 for the route planning of each district. For reasons of transparency and generalizability, the data and code used have been made accessible on the Github page of the first author (Reusken,2023). A number of computational experiments were carried out to assess the performance of the matheuristic. All computations were executed on the same machine equipped with two Intel®Xeon®Platinum 8268 processors (2.9 GHz base clock), and a limit of 128 GB of RAM on each processor. In the following, we first describe the instances and inputs that have been used to test the matheuristic, before providing an analysis of the numerical results from the experiments. 5.1. Description of the instances The analysis was performed using a randomly generated data set, inspired by the instances proposed by Solomon (1987), and several real-life data sets, originating from food banks in the Netherlands and Canada. 5.1.1. Random instances The random instances used for our experiments originate from the well-known benchmark instances created by Solomon (1987). These instances consider three types of geographical distributions of the customer locations: clustered (C), random (R) and a mixture of random and clustered (RC). For the main analysis in this research, we selected the first data sets of types C and R with unique customer locations, i.e., data sets c101, c201 and r101. We created subsets to increase the number and variety of instances by extracting the first 20 to 100 customer locations from these data sets, with intervals of 10, resulting in a total of 27 instances. The travel distance (𝑐𝑖𝑗 ) between the locations is computed based on the Euclidean distance. Fig. 2 provides a graphical representation of the three data sets, showing the distribution of customers over the service region. Depicting the first 20 customers in lighter color, this graphical representation also illustrates the impact on the spread of customers within the service region when using only a small subset of the customer locations. In addition to the customer and depot locations, we also used the same vehicle capacity (𝑄) and customer demands as in Solomon (1987). Because the customer demands in these instances are deterministic, we used them as a mean (𝜇𝑖, 𝑖 ∈𝑁) for the stochastic demands in our problem, which we assume to be normally distributed with a standard deviation of 𝜎𝑖= 0.2𝜇𝑖, 𝑖 ∈𝑁. Moreover, since we assume uncertain waiting times, we generated for each customer 𝑖∈𝑁a sensible probability 𝑞𝑖ℎ of having ℎvehicles waiting in the queue. For this purpose, we assume that (i) 𝑞𝑖ℎ varies across customers, and (ii) 𝑞𝑖ℎ Table 3 Parameter settings. Parameter Value 𝑇Travel time limit 8 𝑠Driving speed (km/h) 50 𝛾Safety factor for vehicle capacity to determine 𝑚𝑐0.01 𝛿Safety factor for travel time limit to determine 𝑚𝑡0.01 𝛼Safety factor for vehicle capacity in assignment of customers to districts 0.05 𝜂Safety factor for travel time limit in assignment of customers to districts 0.05 𝑎Service time per demand unit 0.008 is likely to decrease with higher values of ℎ. The precise steps taken to generate the values for 𝑞𝑖ℎ can be found in Appendix B.1. Finally, the remaining input parameters are location-independent and were set to realistic values, summarized in Table 3. In this context, the travel time limit 𝑇, corresponding to a workday of a driver, is set to eight hours, and the average speed of the trucks is assumed to be 50 km/h. Moreover, since in a practical setting the number of districts, and therefore the number of vehicles, cannot be modified on a daily basis, we only allow for a 1% probability of failure to meet the vehicle capacity and travel time limit when determining the number of districts, i.e., 𝛾= 0.01 and 𝛿= 0.01. In contrast, the safety factors for meeting the vehicle capacity (𝛼) and travel time limit (𝜂) in the assignment problem are set more leniently (to values of 0.05), since a route can be changed more easily than the number of vehicles. Furthermore, given that the mean values of the customer demands are between one and 50 for our instances, we consider a service time 𝑎 per demand unit of 0.008 h, which corresponds to approximately five minutes per 10 demand units. 5.1.2. Real-life instances In order to highlight the practical relevance of our matheuristic, we also applied our method to a number of real-life instances, resulting from in-depth field work conducted at a total of six food banks in the Netherlands and Canada (Food Bank Zevenaar,1Food Bank Eemsdelta,2 Food Bank Haaglanden,3Food Bank Almelo,4Moisson Montréal,5and Share the Warmth6). This field work involved several interviews with operational and managerial staff, inspections of equipment and facilities, as well as accompanying a driver during a normal operating day on his collection route. The results of this practical study showed that all of these food banks face a similar collection problem, experiencing high 1See https://www.voedselbankzevenaar.nl/ (in Dutch) 2See https://www.voedselbankeemsdelta.nl/ (in Dutch) 3See https://voedselbankhaaglanden.nl/ (in Dutch) 4See https://www.voedselbankalmelo.nl/ (in Dutch) 5See https://www.moissonmontreal.org/en/ 6See https://sharethewarmth.ca/ European Journal of Operational Research 317 (2024) 111–127 119 M. Reusken et al. Table 4 Summary of the real-life data sets. 1. Organization name Food Bank Zevenaar Food Bank Eemsdelta Moisson Montréal 2. Country The Netherlands The Netherlands Canada 3. Organization type Food agency Food agency Food bank depot 4. Number of instances 1 2 6 5. Number of collection locations 7 17 252 6. Collection area in km2257 773 403 7. Vehicle capacity in kg (𝑄) 2500 1000 4535 8. Travel time limit in hours (𝑇) 4 4 8 9. Range of historic collection demands in kg – [2,64] [5,1241] 10. Service time in hours per kg (𝑎) 0.008 0.008 0.0008 levels of uncertainty in the collection demands, which also affect service and waiting times. Unfortunately only three of the food banks were able to provide sufficient data to generate instances corresponding to the real-life cases, namely Food Bank Zevenaar, Food Bank Eemsdelta and Moisson Montréal. The data provided by these food banks are summarized in Table 4. The first two rows present the organization’s name and the country in which it is based. The third row indicates the role of the organization in the food bank supply chain, where a food agency provides food assistance to the community in need, and a food bank depot supplies food to such food agencies. Consequently, these two organization types typically work with different types of collection locations, i.e., the food agency collects small demands from suppliers in its proximity (e.g, local bakeries) as well as from its supplying food bank depot, while the food bank depot generally collects larger demands from a larger geographical area. Row 4 presents the number of instances that are available for each organization, which corresponds to different sets of collection locations that had to be visited in a given week for the Dutch food agencies, and on each day of the week for Moisson Montréal. Row 5 shows the total number of collection locations. A subset of these locations, along with the depot locations, are mapped in Fig. 3 for illustrative purposes. The travel distance 𝑐𝑖𝑗 between these collection locations is computed based on the geodesic distance, i.e., the length of the shortest path between two points on Earth. The collection area (the circumscribed rectangle encompassing all nodes) is given in square kilometers in row 6 as an indication of the travel distances. The real vehicle capacities in kilograms and the daily travel time limit, corresponding to a workday of a driver, are presented in rows 7 and 8, respectively. Row 9 contains the minimum and maximum values of the historical collection demands (see Appendix B.2 for visualizations of these demands), and we adopt the ranges of Food Bank Eemsdelta for Food Bank Zevenaar as a reference for experimental purposes. The last row of the table displays the service time per demand unit, i.e., the time (in hours) required to collect a kilogram of demand. For the Dutch food agencies, where the collection demands range between two and 64 kg for our instances, we use 𝑎= 0.008, which corresponds to approximately five minutes per 10 kg. This value is significantly lower for the case of Moisson Montréal (i.e., 𝑎= 0.0008), as a food bank depot operates with a higher efficiency, processing larger quantities within a shorter duration (equivalent to approximately 100 kg per five minutes). Considering the historical demands of Food Bank Eemsdelta and Moisson Montréal (see Appendix B.2), we can state that demands do not conform to any established classical distribution, although a vague resemblance to a normal distribution can be observed. For the purpose of computational testing we hence assume demand to be normally distributed using a mean (𝜇𝑖, 𝑖 ∈𝑁) that is randomly sampled from the historic demand ranges indicated in Table 4, and a standard deviation of 𝜎𝑖= 0.2𝜇𝑖, 𝑖 ∈𝑁. Moreover, since we assume uncertain waiting times, yet no historical data were available on this, we generate sensible values for parameter 𝑞𝑖ℎ using the procedure described in Section 5.1.1 (for a detailed description, we refer readers to Appendix B.1). Finally, the values of 𝑠,𝛾,𝛿,𝛼and 𝜂are set to the same realistic values as in the case of the random instances (see Table 3). 5.2. Preliminary experiments on training instances To tune the performance of our matheuristic, a series of preliminary experiments were carried out using random instances that we did not use for our main experiments. Specifically, we used the locations, vehicle capacity and demands of the first RC data set by Solomon (1987), i.e., data set rc101. We generated subsets of sizes {20,30,…,100}, ensuring that each subset consisted of an equal proportion of customers from both the random and clustered groups. The results of these experiments determined which assignment problem reformulation to use, as well as the setup of an efficient procedure reducing the number of variables in the assignment problem, and finally the stopping criteria for our matheuristic. 5.2.1. Reformulation of the assignment problem In Section 4.2.2, we present two computationally tractable reformulations of the assignment problem, a linear and a convex quadratic one. Preliminary experiments showed that the performances of these reformulations are comparable in terms of computational time. However, given the large number of additional variables required for the linear reformulation, we have selected the convex quadratic reformulation for our numerical analysis. 5.2.2. Reduction of the number of variables for the assignment problem Preliminary experiments have shown that most of the computational time is spent on the assignment of customers to districts, as solving this subproblem can take up to several hours depending on the size of the instance. Computational experiments have furthermore indicated that reducing the number of decision variables 𝑥𝑖𝑘, determining the assignment of customer 𝑖to a seed in district 𝑘, can significantly reduce computational times without considerably affecting the solution quality. For this purpose, we set the decision variables 𝑥𝑖𝑘 that correspond to the seeds furthest away from a customer equal to zero, by adding the following constraints to the assignment problem: 𝑥𝑖𝑘 = 0 𝑘∈𝐹𝑖, 𝑖 ∈𝑁, (18) where 𝐹𝑖= {1,…, 𝜅}denotes the set of seeds that are furthest away from customer 𝑖∈𝑁, with 𝐹𝑖⊆ 𝐾. These 𝜅seeds correspond to those that contain the largest insertion costs 𝑐𝑖𝑘. Several values for 𝜅were tested, and the best performance was observed for 𝜅=⌊0.3|𝐾|⌋. It should be noted that we only adopt this procedure of reducing the number of variables for those instances with |𝑁|≥60, since the reduction in computational time for small instances is limited, and the additional constraints may very rarely cause the assignment of customers to districts to be infeasible for small instances. 5.2.3. Stopping criteria The wall clock time limit for the matheuristic is set to 7,200 s, and the termination scalar is set to 𝜓= 10, so that the maximum iteration count is limited to 10|𝐾|. These stopping criteria will only terminate the procedure between iterations, so that an iteration is never interrupted. In addition, since the assignment of customers to districts can be very time consuming for larger instances, several specific stopping criteria European Journal of Operational Research 317 (2024) 111–127 126 M. Reusken et al. Besiou, M., Pedraza-Martinez, A. J., & Van Wassenhove, L. N. (2018). OR applied to humanitarian operations. European Journal of Operational Research,269(2), 397–405. Boschetti, M. A., Maniezzo, V., Roffilli, M., & Röhler, A. B. (2009). Matheuristics: Optimization, simulation and control. In Hybrid metaheuristics: 6th international workshop, HM 2009, udine, Italy, October 16-17, 2009. proceedings 6:vol. 5818, (pp. 171–177). Berlin, Heidelberg: Springer-Verlag. Davis, L. B., Jiang, S. X., Morgan, S. D., Nuamah, I. A., & Terry, J. R. (2016). Analysis and prediction of food donation behavior for a domestic hunger relief organization. International Journal of Production Economics,182, 26–37. Davis, L. B., Sengul, I., Ivy, J. S., Brock, L. G., & Miles, L. (2014). Scheduling food bank collections and deliveries to ensure food safety and improve access. Socio-Econ. Plan. Sci.,48(3), 175–188. Dror, M., Laporte, G., & Louveaux, F. V. (1993). Vehicle routing with stochastic demands and restricted failures. Z. Oper. Res.,37(3), 273–283. Dror, M., Laporte, G., & Trudeau, P. (1989). Vehicle routing with stochastic demands: Properties and solution frameworks. Transportation Science,23(3), 166–176. Eisenhandler, O., & Tzur, M. (2019). A segment-based formulation and a matheuristic for the humanitarian pickup and distribution problem. Transportation Science,53(5), 1389–1408. Eisenhandler, O., & Tzur, M. (2019). The humanitarian pickup and distribution problem. Operations Research,67(1), 10–32. Errico, F., Desaulniers, G., Gendreau, M., Rei, W., & Rousseau, L.-M. (2016). A priori optimization with recourse for the vehicle routing problem with hard time windows and stochastic service times. European Journal of Operational Research,249(1), 55–66. European Food Banks Federation (2022). Assessment of FEBA member’s activities: Technical report. Fama, E. F., & Roll, R. (1968). Some properties of symmetric stable distributions. Journal of the American Statistical Association,63(323), 817–836. FAO, IFAD, UNICEF, WFP, & WHO (2021). The state of food security and nutrition in the world 2021. In Transforming food systems for food security, improved nutrition and affordable healthy diets for all. Rome: FAO. Federgruen, A., & Simchi-Levi, D. (1995). Analysis of vehicle routing and inventoryrouting problems. Handbooks in Operations Research and Management Science,8, 297–373. Fisher, M. L., & Jaikumar, R. (1981). A generalized assignment heuristic for vehicle routing. Networks,11(2), 109–124. Franceschetti, A., Jabali, O., & Laporte, G. (2017). Continuous approximation models in freight distribution management. TOP,25(3), 413–433. Gendreau, M., Jabali, O., & Rei, W. (2016). 50Th anniversary invited article—Future research directions in stochastic vehicle routing. Transportation Science,50(4), 1163–1173. Goodson, J. C. (2015). A priori policy evaluation and cyclic-order-based simulated annealing for the multi-compartment vehicle routing problem with stochastic demands. European Journal of Operational Research,241(2), 361–369. Gunes, C., Van Hoeve, W. J., & Tayur, S. (2010). Vehicle routing for food rescue programs: A comparison of different approaches. In Integration of AI and OR techniques in constraint programming for combinatorial optimization problems: 7th international conference, cPAIOR 2010, bologna, Italy, June 14-18, 2010. proceedings 7(pp. 176–180). Springer, Berlin, Heidelberg. Haugland, D., Ho, S. C., & Laporte, G. (2007). Designing delivery districts for the vehicle routing problem with stochastic demands. European Journal of Operational Research,180(3), 997–1010. Hoogendoorn, Y. N., & Spliet, R. (2023). An improved integer L-shaped method for the vehicle routing problem with stochastic demands. INFORMS Journal on Computing, 35(2), 423–439. Kaya, O., & Ozkok, D. (2020). A blood bank network design problem with integrated facility location, inventory and routing decisions. Networks and Spatial Economics, 20(3), 757–783. Keskin, M., Çatay, B., & Laporte, G. (2021). A simulation-based heuristic for the electric vehicle routing problem with time windows and stochastic waiting times at recharging stations. Computers & Operations Research,125, Article 105060. Keskin, M., Laporte, G., & Çatay, B. (2019). Electric vehicle routing problem with timedependent waiting times at recharging stations. Computers & Operations Research, 107, 77–94. Koskosidis, Y. A., & Powell, W. B. (1992). Clustering algorithms for consolidation of customer orders into vehicle shipments. Transportation Research, Part B (Methodological),26(5), 365–379. Laporte, G., & Louveaux, F. V. (1993). The integer L-shaped method for stochastic integer programs with complete recourse. Operations Research Letters,13(3), 133–142. Laporte, G., Louveaux, F. V., & Van Hamme, L. (2002). An integer L-shaped algorithm for the capacitated vehicle routing problem with stochastic demands. Operations Research,50(3), 415–423. Lei, H., Laporte, G., & Guo, B. (2012). A generalized variable neighborhood search heuristic for the capacitated vehicle routing problem with stochastic service times. TOP,20(1), 99–118. Li, X., Tian, P., & Leung, S. C. (2010). Vehicle routing problems with time windows and stochastic travel and service times: Models and algorithm. International Journal of Production Economics,125(1), 137–145. Lien, R. W., Iravani, S. M. R., & Smilowitz, K. R. (2014). Sequential resource allocation for nonprofit operations. Operations Research,62(2), 301–317. Macqueen, J. B. (1967). Some methods for classification and analysis of multivariate observations. In Proceedings of the 5th berkeley symposium on mathematical statistics and probability (pp. 281–297). University of California Press. Mahmoudi, M., Shirzad, K., & Verter, V. (2022). Decision support models for managing food aid supply chains: A systematic literature review. Socio-Economic Planning Sciences,82. Marinakis, Y., Iordanidou, G. R., & Marinaki, M. (2013). Particle swarm optimization for the vehicle routing problem with stochastic demands. Applied Soft Computing, 13(4), 1693–1704. Mendoza, J. E., Castanier, B., Guéret, C., Medaglia, A. L., & Velasco, N. (2010). A memetic algorithm for the multi-compartment vehicle routing problem with stochastic demands. Computers & Operations Research,37(11), 1886–1898. Mendoza, J. E., Rousseau, L.-M., & Villegas, J. G. (2016). A hybrid metaheuristic for the vehicle routing problem with stochastic demand and duration constraints. Journal of Heuristics,22(4), 539–566. Miranda, D. M., & Conceição, S. V. (2016). The vehicle routing problem with hard time windows and stochastic travel and service time. Expert Systems with Applications, 64, 104–116. Nair, D. J., Grzybowska, H., Fu, Y., & Dixit, V. V. (2018). Scheduling and routing models for food rescue and delivery operations. Socio-Economic Planning Sciences, 63, 18–32. Nair, D. J., Grzybowska, H., Rey, D., & Dixit, V. (2016). Food rescue and delivery: Heuristic algorithm for periodic unpaired pickup and delivery vehicle routing problem. Transportation Research Record,2548, 81–89. Nair, D. J., Rey, D., & Dixit, V. V. (2017). Fair allocation and cost-effective routing models for food rescue and redistribution. IISE Transactions,49(12), 1172–1188. Orgut, I. S., Ivy, J. S., Uzsoy, R., & Hale, C. (2018). Robust optimization approaches for the equitable and effective distribution of donated food. European Journal of Operational Research,269(2), 516–531. Orgut, I. S., & Lodree, E. J. (2023). Equitable distribution of perishable items in a food bank supply chain. Production and Operations Management. Oyola, J., Arntzen, H., & Woodruff, D. L. (2017). The stochastic vehicle routing problem, a literature review, part II: solution methods. EURO Journal on Transportation and Logistics,6(4), 349–388. Oyola, J., Arntzen, H., & Woodruff, D. L. (2018). The stochastic vehicle routing problem, a literature review, part I: models. EURO Journal on Transportation and Logistics,7(3), 193–221. Paul, S., & Davis, L. B. (2022). An ensemble forecasting model for predicting contribution of food donors based on supply behavior. Annals of Operations Research, 319(1), 1–29. Reusken, M. (2023). Matheuristic-for-CVRP-SDSW [source code]. URL https://github. com/MeikeReusken/Matheuristic-for-CVRP-SDSW. Reusken, M., Cruijssen, F., & Fleuren, H. (2023). A food bank supply chain model: Optimizing investments to maximize food assistance. International Journal of Production Economics,261, 108886. Rey, D., Almi’ani, K., & Nair, D. J. (2018). Exact and heuristic algorithms for finding envy-free allocations in food rescue pickup and delivery logistics. Transportation Research Part E: Logistics and Transportation Review,112, 19–46. Rivera, A. F., Smith, N. R., & Ruiz, A. (2023). A systematic literature review of food banks’ supply chain operations with a focus on optimization models. Journal of Humanitarian Logistics and Supply Chain Management,13(1), 10–25. Secomandi, N., & Margot, F. (2009). Reoptimization approaches for the vehicle-routing problem with stochastic demands. Operations Research,57(1), 214–230. Solak, S., Scherrer, C., & Ghoniem, A. (2014). The stop-and-drop problem in nonprofit food distribution networks. Annals of Operations Research,221(1), 407–426. Solomon, M. M. (1987). Algorithms for the vehicle routing and scheduling problems with time window constraints. Operations Research,35(2), 254–265. Sultana, T., Akhand, M. A. H., & Rahman, M. M. H. (2017). A variant of Fisher and jaikumar algorithm to solve capacitated vehicle routing problem. In 2017 - 8th international conference on information technology (ICIT), proceedings (pp. 710–716). IEEE. European Journal of Operational Research 317 (2024) 111–127 127 M. Reusken et al. Sungur, I., Ordóñez, F., & Dessouky, M. (2008). A robust optimization approach for the capacitated vehicle routing problem with demand uncertainty. IIE Transactions, 40(5), 509–523. Tillman, F. A. (1969). The multiple terminal delivery problem with probabilistic demands. Transportation Science,3(3), 192–204. United Nations (2022). The Sustainable Development Goals Report:Technical report. Yang, W. H., Mathur, K., & Ballou, R. H. (2000). Stochastic vehicle routing problem with restocking. Transportation Science,34(1), 99–112. Zhang, J., Lam, W. H., & Chen, B. Y. (2013). A stochastic vehicle routing problem with travel time uncertainty: trade-off between cost and customer service. Networks and Spatial Economics,13(4), 471–496.