Last mile deliveries with lockers: formulations and algorithms

In this paper, we consider the use of lockers in parcel delivery, a recent method used in last mile logistics. Lockers are pickup points made of several cells that are located in several points of a city where customers can collect their parcels as an alternative to home delivery. We study routing problems in which one or multiple vehicles are used to deliver parcels directly to customers or lockers. We also study the influence of the introduction of lockers when these problems include time windows. We propose a set of novel formulations for these problems, some valid inequalities, and a branch-and-cut algorithm. Moreover, we investigate the difference between the routing problems with lockers and the classical routing problems.


Introduction
Retail e-commerce sales worldwide have experienced a continuous growth from 2014 to 2020, the year in which the sales value has more than doubled with respect to 2016, despite the current pandemic, and forecasts expect the sales value to be more than three times as high in 2023 (see, e.g., Statista 2020). One of the main challenges introduced by the e-commerce boom that retailers and logistics companies must overcome is the quick delivery of a large number of parcels at minimum cost, with a particular focus on the last mile delivery. To do so, they are implementing alternative delivery methods with respect to home delivery. One of the strategies put in place by these companies is the use of pickup points (see, e.g., Mitrea et al. 2020) located in places that are convenient to the customers in which couriers deliver the parcels that will be collected by the customers at a later time. The introduction of pickup points for the last mile deliveries as an alternative to home delivery offers two main advantages: the reduction in traveling costs and external costs due to the fact that more parcels are delivered in the same place and thus the couriers can travel shorter routes; moreover, the lack of need of synchronization between the couriers and the customers allows more flexibility in both delivery and collection times that hence decreases the amount of failed deliveries.
Pickup points can be divided into two types: counters and lockers (see, e.g., Amazon 2020). Counters are pickup points located in shops, in which parcels are handled manually by the shopkeeper and thus are subject to opening hours. Lockers are (automated) self-service facilities which are located externally and allow customers to autonomously collect parcels 24 hours a day, seven days a week.
The number of papers considering logistics and optimization problems that include pickup points is limited but increasing. In this paper, we tackle a last mile problem in which one or more identical vehicles can either deliver parcels directly to the customers or to nearby lockers (the adaptation to the counter case is straightforward and can be done by setting more restrictive time windows to the pickup points in order to respect the shopkeepers working hours). We study the cases with and without time windows for which we propose a set of formulations and branch-and-cut algorithms.
The contributions of this paper are the following: (i) We propose enhanced exact methods to solve the traveling salesman problem with lockers, the traveling salesman problem with lockers and time windows, the vehicle routing problem with lockers, and the vehicle routing problem with lockers and time windows. Note that most of the papers in the related literature propose heuristic or metaheuristic algorithms, and that the mathematical formulations are sometimes used just as a detailed description of the problem. (ii) We define a set of valid inequalities for the proposed formulations of the studied problems. (iii) We consider the capacity of the trucks as a constraint that limits the routes, despite this is not considered in the relevant literature. (iv) We provide a comparison among several exact methods for the proposed problems. (v) We improve the results of other methods from the literature that solve similar problems: our methods were able to obtain new optimal solutions for instances of the literature for the first time and improve several best-known solutions. (vi) We provide a discussion on the costs and benefits linked to the introduction of lockers and time windows in the last mile delivery.
In the remainder of the paper, we firstly propose a detailed literature review in Sect. 2. Section 3 provides an explanation on the considered problems, while in Sect. 4 several formulations are proposed to model these problems. The implementation of the branch-and-cut algorithm is defined in Sect. 5, while Sect. 6 presents computational results comparing the proposed formulations with respect to each problem and provides a comment on the impact of the introduction of lockers in the last mile delivery process. Section 7 concludes the paper.

Literature review
The introduction of pickup points in logistics optimization problems is getting more and more attention; moreover, most of the studies that consider the use of lockers have been published in the last five years. For a wide classification on recent studies, we head the reader to the survey by Rohmer and Gendron (2020). Jiang et al. (2020) study the traveling salesman problem with time windows for the last mile delivery in online shopping, a generalization of the traveling salesman problem (TSP) with time windows in which customers can either receive their parcels at home or pick them up from shared delivery facilities (lockers) nearby. Lockers have a maximum radius in which customers can be served and a maximum capacity. No preference between the two types of service is considered; however, the objective function aims at minimizing the truck traveling costs, the traveling cost of the customer to reach the locker, and the cost of opening a locker. The authors propose a Mixed Integer Linear Programming (MILP) formulation and a variable neighborhood descent. The MILP formulation models hard time windows; however, in order to improve the efficiency of the commercial solver, they only solve the case with soft time windows. The largest instance solved exactly counts 60 customers and 6 lockers.
Most of the other works considers multiple vehicles. Oliveira and dos Santos (2020) study a problem in which a limited number of lockers can be opened and used to serve customers. Two different sets of vehicles are employed to carry out the deliveries: one fleet for home delivery and one for servicing lockers. Customer requests are not necessarily unitary. Lockers have a maximum radius and a capacity, but no opening cost is considered; in fact, the problem minimizes the traveling costs. The authors propose a MILP formulation and solve instances with up to 75 customers and 10 lockers; they also introduce a variable neighborhood descent metaheuristic. Osaba et al. (2018) study a multi-objective optimization problem allowing for delivery at third-party drop-off points instead of direct delivery at the customer locations. The problem aims at optimizing bike routes under considerations of total profit, route distance, and the safety of the bikers. They propose six different evolutionary heuristics and solve three instances with 28 customers and 13 drop-off points. Mancini and Gansterer (2021) propose the VRP with private and shared delivery locations where customers can either be served at home within a set time window or at a locker with an unlimited homogeneous fleet of uncapacitated trucks. The truck routes are constrained by a maximum duration. Customers can select one of the two types or let the company choose between the two types of services. The service via locker implies a compensation to the customer and thus has a cost for the company. They define a maximum capacity for the lockers, but the used instances have lockers with a very large capacity. The authors propose a MILP model, a large neighborhood search (LNS)-based matheuristic, and an iterated local search algorithm. They study different time windows preferences, different radii of service, and different compensation schemes. They could solve exactly most of the 50 customer instances, and heuristic solutions are provided for the 75 customer instances. Grabenschweiger et al. (2021) propose the VRP with heterogeneous locker boxes, where customers have multiple requests that can be delivered either at home (respecting a time window) or at a locker by an unlimited homogeneous set of uncapacitated trucks. Trucks must travel for at most a set maximum duration. Customers can be served by the closest locker (if they allow so) for which they receive compensation, that is a cost to the company. If the delivery happens at a locker, then all the parcels of the same customer are stored at the same locker. The authors consider the case with lockers having boxes (or cells) of a single dimension or of different dimensions. They propose two MILP models, each of which considers one of the two features just described and an Adaptive LNS-based metaheuristic. They could solve exactly most of the 25 customer instances, and heuristic solutions are provided for the 75 customer instances, improving some large solutions with respect to Mancini and Gansterer (2021). Orenstein et al. (2019) propose a vehicle routing problem (VRP), called flexible parcel delivery problem, for delivering parcels exclusively to lockers. Truck routes are limited by a maximum shift length. Each parcel is associated with a set of possible destinations (lockers). Parcels have different sizes, and both lockers and trucks are divided into cells of different sizes. No more than one parcel can be accommodated in one cell, and a cell must be sufficiently large to contain a parcel. The objective is to minimize the traveling costs and the penalties incurred in case of failed delivery. The authors propose a MILP formulation and two metaheuristics based on the Clarke and Wright's (C&W) Savings and the Petal heuristics. They consider instances with up to 50 lockers and 1500 parcels, but they solve to optimality only three instances with 20 lockers and 200 parcels. Wen and Li (2016) propose a VRP with soft time windows for last mile delivery and study the impact of the introduction of lockers. The problem aims at minimizing the variable transportation costs, the soft time windows penalty, the C O 2 emission costs, and the fixed cost of labor and vehicles. Once the lockers are introduced, the time windows constraints are not considered anymore, but the locker renting costs and the new small vans acquisition costs are included in the objective function. The authors model the problem with MILP formulations and present a C&W Savings-based heuristic to solve a real-world instance with 12 customers. Dumez et al. (2021) present the vehicle routing problem with delivery options (VRPDO), a VRP with time windows in which several delivery options such as home, workplace, locker, or car trunk are considered. The authors model customer preferences: the customers define their delivery preferences by ordering the options. Both vehicles and lockers are capacitated, and a maximum travel time is imposed on trucks. The objective is to define the truck routes minimizing the traveling costs and the number of used vehicles while guaranteeing the minimum satisfaction level of the customers. The authors propose a formulation, but they solve the problem with a dedicated LNS that heuristically solves instances with up to 200 customers. Tilk et al. (2021) propose a branch-and-price-and-cut algorithm to solve the VRPDO that solves instances with up to 60 customers.
Sitek and Wikarek (2019) consider both pickup and delivery proposing a capacitated VRP (CVRP) called CVRP with pick-up and alternative delivery in which heterogeneous vehicles are used to serve customers at home, at a postal office, or at a capacitated locker. Vehicles must respect a maximum shift duration. The objective is that of minimizing the distance traveled by the trucks and the penalties linked to customer preferences. The authors propose an ILP formulation and a heuristic algorithm, both enhanced by a constraint programming preprocessing. They could solve instances with 20 vehicles, 200 delivery points, and 2000 parcels. Zhou et al. (2018) propose the multi-depot two-echelon VRP with delivery options for the last mile distribution; in the first echelon, the parcels are moved from the main depots to the intermediate ones (satellites), while the parcel delivery is performed in the second echelon. Customers can define a preference for a set of lockers or their home. Two different homogeneous capacitated fleets are used for the two echelons, and a maximum traveling time is imposed on each vehicle. Satellites have a maximum capacity, while lockers do not; however, both can be visited by several vehicles. The objective is to minimize the fixed costs of the depots and the fixed and variable costs of vehicles and satellites. The authors propose a hybrid multi-population genetic metaheuristic. The larger instance used includes two depots, 12 satellites, 30 lockers, and 200 customers. Enthoven et al. (2020) also solve a two-echelon problem: the first echelon accounts for deliveries at satellites or directly to lockers, while in the second echelon they consider deliveries from the satellites to customer houses with zero emission vehicles. The problem minimizes the traveling costs and the cost linked to the lockers. The authors propose a MILP formulation and an adaptive LNS. Their model can solve instances with up to 51 customers, four satellites, and four lockers. Zhou et al. (2019) present the green VRP considering dual services (home delivery and customer pickup) with stochastic travel times. A fleet of capacitated identical vehicles is used to perform the deliveries to the customers that are associated with non-unitary requests and that are already divided into home delivery and self-collection at uncapacitated lockers. Time windows and a maximum traveling time for each vehicle are considered and are subject to stochastic traveling times. The objective is to minimize the transportation costs and the expected recourse cost that includes the violation of the time windows and of the maximum travel time. The problem is formulated as a two-stage stochastic optimization model with recourse strategy and solved by applying a twostage heuristic algorithm. They consider instances with up to 120 customers.
In the following, we focus on the TSP with lockers and time windows and we extend it to the VRP case, as the majority of the papers that consider lockers. Most of the papers described in this section solve the proposed problems with heuristic or metaheuristic algorithms, sometimes only reporting a mathematical formulation to provide a detailed description of the problem. In the following, we propose a set of formulations to solve the tackled problems with exact methods, a set of valid inequalities, and a comparison between our methods and exact methods proposed in the literature to solve similar problems. Moreover, we analyze the obtained results and the impact of lockers and time windows on the solutions.

Problems definition
In this section, we describe the set of problems that we will tackle in the following. We firstly describe the traveling salesman problem with lockers (TSPL), which is a generalization of the TSP in which some customers can be served by lockers. We thus consider the case with multiple capacitated vehicles, the vehicle routing problem with lockers (VRPL). Afterward, we consider the single and the multiple vehicle cases with time windows, the traveling salesman problem with locker and time windows (TSPLTW), and the vehicle routing problem with lockers and time windows (VRPLTW), respectively. The problem description follows this outline.

The single vehicle case: the TSPL
The TSPL is modeled on a digraph G = (V , A). The set of vertices V is partitioned as V = {0}∪ V C ∪ V L , where 0 is the depot, V C = {1, 2, . . . , n} is the set of the n customers and V L = {n + 1, n + 2, . . . , n + m} is the set of the m available lockers. A = {(i, j) : i, j ∈ V } is the given set of arcs, each of which is associated with a distance d i j and a traveling time The problem aims at delivering the parcels to the customers with a truck at minimum cost. The truck must leave from and return to the depot and it can serve the customers either directly at home (home service) or via a set of lockers (self-service). Each customer must be served, but both customers and lockers can be visited at most once by the truck. Each locker b ∈ V L is associated with a maximum capacity defined by the number of customers v b , and with a radius R b beyond which it cannot serve any customer. Note that some customers may lay within the radius of more than one locker, but some may be located beyond the radius of all the lockers, and thus, no self-service is possible for them.
The objective function to be minimized is made of two components: (i) the truck transportation costs and (ii) the cost sustained because the customers must collect their parcels in a locker, each of which is multiplied by a coefficient: α 1 and α 2 , respectively.
A feasible solution for the TSPL is presented in Fig. 1, in which the square is the depot, the circles are the customers, and the triangles are the lockers with their radius. The arrows represent the truck route, and the dashed arrows represent the customer's path to collect their parcels at the lockers.

Extension to the multiple vehicle cases: the VRPL
In order to extend the TSPL to the case that includes multiple capacitated vehicles, the VRPL, we introduce a fleet composed of at most n identical vehicles with capacity Q that, starting and returning from the depot, deliver parcels by visiting lockers or customers. Each customer has a nonnegative request q i . Note that also for the VRPL the capacity of lockers is expressed by the number of customers and is not affected by their requests. All the other characteristics are inherited by the TSPL. Among these, note that each vertex but the depot can be visited at most once. This prevents a

Time windows: the TSPLTW and the VRPLTW
We extend the TSPL and the VRPL to include a time window [e i , l i ] for each customer i ∈ V C within which a customer must be served. Note that a truck can reach a customer i before its earliest time window e i , but it must wait until e i to serve it. The earliest visit time of the depot e 0 is set to the beginning of the time horizon 0 (if e 0 = 0 one can always shift all the time windows by e 0 ), while l 0 defines a maximum time that each route cannot exceed. Pickup points can have time windows [e b , l b ], b ∈ V L ; however, no time window should be imposed to lockers, and for them we set [e b , l b ] = [e 0 , l 0 ]. Time windows constraints can be applied to both the TSPL and the VRPL, defining two new problems: the TSPLTW and the VRPLTW, respectively.
In Table 1, one can find the summary of all the symbols used in this paper.

Formulations
In this section, we first propose a set of different formulations for the VRPL and their adaptation to the VRPLTW. Afterward, we show their adaptation to the TSPL and the TSPLTW. A set of valid inequalities is introduced for all formulations.

Formulations for the VRPL
In the following, we propose three formulations for the VRPL.
Weight of the truck transportation costs α 2 Weight of the self-service costs Variables x i j Binary variable that equals 1 if a vehicle travels the arc (i, j) ∈ A; 0 otherwise

The vehicle flow formulation
The first formulation that we present is based on the wellknown vehicle flow formulation for the VRP (see, e.g., Toth and Vigo 2002), we call it vehicle flow formulation (VFF), and it makes use of two sets of variables. Variables x i j equal 1 if a truck travels arc (i, j) ∈ A, 0 otherwise. Variables y ib equal 1 if the customer i ∈ V C travels to locker b ∈ V L to collect its parcel, 0 otherwise. Note that x ii = 0, i ∈ V and variables y ib are fixed to 0 if d ib > R b for this formulation and those defined in the following.
i∈S j∈S The objective function (1) minimizes two components, each one multiplied by its weight: the distance traveled by the truck and a penalty associated with self-service that depends on the customer and the locker used to serve it. Constraints (2) guarantee the flow conservation in every vertex, while constraints (3) impose that a customer is either served at home by a truck or at a locker. Constraints (4) impose that if a locker is not used, then no vehicle may visit the locker.
Constraints (5) impose that if a locker is used, then one arc entering and one arc leaving the locker must be covered by a truck. Constraints (6) state that if no customer is served via a certain locker, then no truck should visit that locker. Constraints (7) impose a maximum capacity to each locker.
In (8), we adapted the generalized subtour elimination constraints (GSEC) to the VRPL. Variables x and y are defined binary in (9) and (10).
Inequalities A valid inequality imposing a lower bound on the number of vehicles that must leave the depot to serve all customers as in (11) can be added to strengthen the formulation. GSEC (8) are exponentially many and will be included only when violated in a branch-and-cut fashion (for details, see Sect. 5); however, one can statically include GSEC for |S| = 2 as in (12) and (13).
Constraints (14) strengthen the fact that if a locker serves any client, then it must be visited.

The one-commodity formulation
The second formulation that we propose is based on the wellknown one-commodity formulation for the VRP which we call one-commodity formulation (OCF). This formulation uses the previously presented variables x and y and their initial fixing. In addition, it needs a non-negative variable f i j that represents the flow carried by the truck on the arc (i, j) ∈ A. The OCF for the VRPL is given by (1-7), (9), (10) and by the following constraints (15-19) Constraint (15) imposes the amount of parcels returning to the final depot to be zero, which means that all the parcels must be delivered. In constraint (16), the amount of parcels leaving the depot equals the sum of all requests. Constraints (19) guarantee the flow to respect the vehicle capacity on all the arcs and the non-negativity constraint for the flow variables. This formulation differs from the traditional onecommodity formulation for the VRP because we need two separate flow conservation constraints for the customers subject to home delivery or self-service. Constraints (17) guarantee the flow conservation for each customer that is served directly at home, and thus, the request q i must be satisfied by this flow; if the customer is served with a locker, then no request is left by the truck. Constraints (18) guarantee the flow conservation for the lockers that must account for all the requests q i served by the considered locker.

The two-commodity formulation
The third formulation that we present is based on the twocommodity formulation originally proposed for the CVRP by Baldacci et al. (2004). We call this formulation the twocommodity formulation (TCF). TCF inherits the variables x, y, and f from OCF and uses a fourth variable g ji that gives the residual space (the empty space) on the vehicle traveling along arc (i, j) ∈ A, if any. Note that variables f show the load flow along the routes traveled by vehicles from their start to their return to the depot, while variables g show the residual flow in the opposite direction.
The TCF for the VRPL is given by (1-7), (9), (10), (15), (16), (19) and by the following constraints (20-25): Constraints (20) and (21) impose an upper bound on the amount of residual flow leaving and entering the depot, respectively. On each arc used by a truck, the sum of the load flow and the one of the residual flow must equal the capacity of the truck as per constraints (22). Constraints (23) impose that the difference between the load flow and residual capacity flow that enter and leave a vertex is the double of the request if the customer is served directly. Constraints (24) guarantee the same flow conservation for the lockers. Nonnegativity and a maximum capacity are imposed to variables g in (25).  can be included to strengthen the TCF for the VRPL.

Formulations for the VRPLTW
In the following, we define the components needed to include time windows in the previously presented formulations and thus to model the VRPLTW.

The vehicle flow formulation
A non-negative continuous variable t i is used to represent the time at which the truck visiting the vertex i ∈ V \ {0}, that can be either a customer or a locker, can leave. t 0 represents the moment at which the last truck returns to the depot, since each shift begins at the depot at time e 0 . Note that t does not model the arrival time at a customer vertex because the truck is allowed to arrive early and wait the beginning of its time window. The following constraints, added to (1-10), model the VRPLTW and represent its VFF.
Constraints (26) update the time variables when an arc (i, j) ∈ A, i = 0 is traveled and constraints (27) update the same variables when the arcs start from the depot. In all these constraints, Z i j is a sufficiently large number that we compute as l i +τ i j −e j , (i, j) ∈ A. Constraints (28) and (29) strengthen the classical time window constraints e i ≤ t i ≤ l i for the customers and the lockers, respectively. In (30), the time variables are defined non-negative. Inequalities based on time windows In the following, we write some valid inequalities to tighten or easily impose time windows in a VRPLTW.
Constraints (31) fix to zero all the arc variables starting from i arriving in j such that the minimum arrival time in j is larger than l j . Let us define for each couple of vertices max(e k +τ ki , e i )+τ i j > l j } that includes all the vertices k for which a path (k, i, j) is infeasible due to the time windows. Then, inequalities (32) and (33) are valid for the VRPLTW.
If a customer i is visited just after locker b on a route, then the time variable at the locker should be bounded so as to arrive before l i . This is expressed in inequalities (34). On the other hand, if a locker b is visited after vertex i, then the time variable in b cannot be smaller than the time needed to travel from i to b plus the earliest time window of i, e i . This is imposed by inequalities (35).

Adaptation to the TSPL
Previously presented formulations can be adapted to solve the TSPL and the TSPLTW. Let us define constraints (36) that impose the number of vehicle exiting and entering the depot to one and constraints (37) that represent the GSEC specific to the TSPL.

Branch-and-cut implementation
The VFF formulations include an exponential number of constraints, and hence, we solve them with a Branch-and-Cut (B&C) algorithm. Formulations OCF and TCF do not need a dynamic insertion of constraints; even if the introduction of the exponential constraints used in VFF and separated on fractional solutions is still possible, computational evidence showed that their use in a B&C fashion is not profitable. We use the B&C framework of CPLEX 20.10.0, which solves at every node of the enumeration tree the linear relaxation of a MILP model and then, invokes user developed separation procedures to possibly add cuts. Preliminary tests suggested to adopt strong branching as branching rule. At the root node of the B&C procedure, we obtain an initial solution by using a simple greedy heuristic. We first build a TSP or a VRP solution (ignoring the lockers) by using the well-known Nearest Neighborhood (NN) (see, e.g., Toth and Vigo 2002) and the C&W algorithm (see, e.g., Clarke and Wright 1964). We adapted the NN to solve the VRP by starting a new truck route when the maximum capacity of the current truck was reached, or no new customer could be included. The C&W algorithm is applied to solve the TSP without fixing the truck capacity until only one route remains. Moreover, we adapted NN and C&W to provide solutions to problems with time window constraints. For the NN, we check the vertices that can be visited without exceeding their latest time window, at each iteration. In case, the arrival at those vertices precedes their earliest time window, the truck needs to wait the socalled waiting time needed to respect them. The following vertex selected to be visited in the truck sequence minimizes both the distance traveled by the truck and the waiting time. For the C&W, we considered the time window constraints by performing only merging moves that allow respecting them. As in the NN, we have considered both the traveled distance and the waiting time in the evaluating function. Once the solutions of the NN and the C&W (with or without time windows) are obtained, lockers are included in two ways: by performing an assignment based on customers or an assignment based on lockers. In the first case, the algorithm tries to move a customer from the truck route to a locker so that the entire solution remains feasible and the objective function is minimized. A new customer move is checked at each iteration. In the second case, for each unused locker, the algorithm tries to find the smallest set of customers such that their assignment to the locker would improve the objective function. Among these feasible moves, we chose those using the locker that provides the larger improvement in the objective function and among these the one with a smaller number of customers. If possible, we reiterate. Among the possible solutions, we select the one with the smallest value to pass to the corresponding exact algorithm.
Note that it was not always possible to obtain a feasible initial solution for the TSPLTW because of the limitation of having one vehicle and the time window constraints. This is not the case for the TSPL, VRPL, and VRPLTW, for which an initial greedy solution was always at hand. For the TSPL, it is always possible to find a feasible initial solution because any solution that visits all the customers with one truck is feasible, while for the VRPL and the VRPLTW it is always possible to find a feasible initial solution because the fleet is assumed to be as large as the number of customers in such problems.

Separation procedures
The aim of this section is to present the procedures that we use to determine if the valid inequalities that we proposed are violated by a given possibly fractional solution (x, y).

Separation S1
To separate GSEC (37), we first build a supporting graph O(n) max-flows are computed using the depot as a source and any vertex i ∈ V C , in turn, as a sink. If the max-flow value obtained is lower than one, then i is disconnected from the depot, and hence, the cut that corresponds to the set S induced by the min-cut may be added to the model. Moreover, for the VRP case, we obtain a simple algorithmic improvement by checking, for any set S induced by a min-cut, if the minimum number of vehicles required to serve S is lower than the max-flow value obtained, and, in such a case, we add the corresponding (stronger) constraint (8). In this way, the exact separation of the subtour elimination constraints (37) is used as a heuristic for separating (8).
We refer to this procedure as separation S1.

Separation S2
To separate exactly constraints (8), we start by adding an auxiliary vertex a to graph G obtaining a new supporting graph G in which a set of arcs (a, i) links the vertex a to every vertex i ∈ V C , with capacity q i /Q. We thus compute the max-flow on G with vertex a as the source and the depot as the sink. The constraint (8) corresponding to the set S induced by the min-cut is then checked, and, if violated, it is added to the model.
We refer to this procedure as separation S2.
In Table 2, we provide a summary of the presented formulations to solve the proposed problems. Note that, for each cell of the table, the first set of constraints is the one needed to obtain a feasible solution for the problem, while the set of constraints after a + includes the non-necessary inequalities that improve formulations performances. These additional inequalities are used to speed up each formulation, except for constraints (14), which are not included in VFF.
The described separation procedures are invoked at every node of the enumeration tree in case of integer solutions and at every node with depth not larger than 8 for fractional solutions, basing our decision on preliminary tests. S1 is used to separate constraints (37) in formulations VFF for the TSPL and TSPLTW. S1 and S2 are used to separate constraints (37) and (8) in formulations VFF for the VRPL and VRPLTW. On the basis of computational evidence on the instances that we tested, we stop the separation process as soon as we find and add a violated cut, if any. The order of the separation procedures, when both are used, follows the one in which they are presented, because S1 can separate exactly constraints (37) and heuristically constraints (8), and thus, separation S2 is invoked only in case no violated cut is found by S1, to separate constraints (8) exactly.

Computational results
In this section, we show the computational results obtained by running the proposed methods on a single thread of CPU Intel i7-6850K @ 3.60 GHz with 64 GB RAM. All our methods have been implemented in C++ and compiled with gcc 7.5.0. As already mentioned, the used solver is CPLEX 20.10.0.

Benchmark instances
Routing problems with lockers are very recent, and no standard set of benchmark instances has been established yet. We make use of the instances used by Jiang et al. (2020), who adapted 40 benchmark instances for the TSP with time windows to the TSPLTW by including a set of locker vertices to the graph. In particular, the original instances are divided into two sets: the first set with time windows of amplitude of 20 was initially proposed by Dumas et al. (1995); the second set with time windows of 100 was proposed by Gendreau et al. (1998) and Dumas et al. (1995). Each of these sets can be divided into four groups of five instances each, with the same number of customers n = 20, 40, 60, 100. The number of added lockers is m = n/10 . Each instance is modified to include m lockers by adding the first m vertices of the following set of coordinates:  (0, 0). Note that the instances include n + m + 1 vertices. The distance d i j is computed as an Euclidean distance between each couple of vertices (i, j) ∈ A, rounded to the second decimal. The distance matrix is corrected so as to respect the triangular inequality, as in Jiang et al. (2020). The following characteristics are also derived from the instances proposed by Jiang et al. (2020). The traveling time is set to equal the distance: τ i j = d i j , (i, j) ∈ A. The cost of serving customers via lockers also depends on the distance: (1-10), S1, S2+ (11-13) ( 1-7), (9), (10), (15-19) + (11-14) ( 1-7), (9), (10), (15) (1-7), (9), (10), (36), (37), S1 + (13) ( 1-7), (9), (10), (36), (38-41) + (13), (1-7), (9), (10) For each instance and problem, we show the needed constraints and separation procedures. After the + sign, we report the valid inequalities For all instances, the parameters are set as follows: α 1 = 1, α 2 = 0.5, and R b = 20, v b = 5, b ∈ V L . We extended time windows to the lockers impos- Note that by imposing tighter time windows one can model counters instead of lockers. For the multiple vehicle case, we considered q i = 1, i ∈ V C and a capacity Q = n/2 .

TSPL and VRPL
In Tables 3 and 4, we show the results obtained by running the proposed formulations to solve the TSPL and the VRPL, respectively. We grouped the results per instance size (number of customers) and for each formulation, we report the percentage gap (computed as 100 · (U B − L B)/U B, being U B the best-known upper bound of each instance and L B the lower bound obtained at the end of the iterations), the seconds needed to reach the optimal solution or the time limit of 2 hours, and the number of optima. The formulation providing the best results when solving the TSPL is OCF: it could solve almost all the instances to the proven optimality in less than 5 minutes, on average. TCF and VFF could solve, respectively, 34 and 32 out of 40 instances with very small gaps but with larger computing times with respect to OCF, on average. However, for the smaller instances, VFF is the fastest method among all.
The proposed formulations provide good results when solving the VRPL, showing small optimality gaps also for the largest instances; however, only OCF could solve one instance with n = 100 to optimality within the time limit. OCF is the best performing of our formulations, being able to solve 19 out of 40 instances and showing an average %gap of 1.64. Note that TCF shows better results than VFF, even if VFF is, once more, quicker in solving the smaller instances.

TSPLTW and VRPLTW
In Tables 5 and 6, we show the results for the TSPLTW and the VRPLTW, respectively. The results are displayed by instance size (n) and time window size (tw). For each formulation, we show the percentage gap (computed as explained above), the seconds needed to solve the instance to optimality or to reach the time limit of 2 hours, and the number of obtained optima. Note that the results are averaged on those instances for which an upper bound was available, that is not always the case for the TSPLTW, because small time windows and the number of vehicles limited to one makes it hard to find feasible solutions, if any; indeed, we suspect that few instances do not have one.
Considering the results for the TSPLTW, one can see that OCF provides the largest number of optima (22 out of 40), the smallest average computing times, and the smallest average optimality gap. The VRPLTW results show a similar behavior, with OCF that can solve 23 out of 40 instances and that provides the best gaps and solving times, on average. Almost all the instances with up to 40 customers are solved to optimality when considering both the TSPLTW and the VRPLTW, but the instances with larger time windows are harder to solve because of the wider solution space. VFF is once more the fastest method when solving the smaller instances for both the single vehicle and the multi-vehicle case. Jiang et al. (2020) propose a Miller-Tucker-Zemlin formulation for a problem that they call the traveling salesman problem with time windows for the last mile delivery in online shopping, a slightly different version of the TSPLTW in which the time windows are imposed also to the service at the lockers. In particular, if a customer is served with a pickup point, then the corresponding parcel must be delivered early enough to allow the customer to travel and collect it within his/her time window; for instance, if customer i is assigned to locker b, b must be visited within l i − τ ib . This  (48), which state that the time in which the customer i served at a locker b is at least the time t b in which the truck leaves the locker plus the time needed by the customer to reach the locker τ ib .

Comparison with Jiang et al.
Moreover, their objective function computes locker setup costs that are proportional to the number of customer served by a locker: for each customer served via locker b ∈ V L a penalty π b must be paid. The authors solve the relaxation of the problem which they propose by including soft time windows constraints, i.e., allowing the latest time window to be exceeded at a cost of a penalty. In order to keep track of the violation of time windows, we introduced a new non-negative variable s i , i ∈ V . The new objective function is thus (49), and constraints (50) and (51) must be imposed to guarantee that s i represents the exceeded time in each vertex, if any. Note that the two new components of the objective function are multiplied by two coefficients α 3 = 5 and α 4 = 1. The In addition, some time window constraints must be modified. Constraints (26)-(48) can be inherited by changing their 'big-M' value to a large enough number Z . One such large number is the value of the initial heuristic upper bound provided to the solver. Inequalities (35) can be inherited as they are while constraints (28) and (29) must be changed as follows: The following inequalities, aimed at bounding from above the value of variables s, can also be added: Constraints (53) provide an upper bound on each variable s i , i ∈ V , while constraints (54) bound variables s b to a large value for each locker used, while it sets the same variable to 0 if the locker is not used. Note that a valid large number Z is, once more, the value of the initial heuristic upper bound on z . Despite the described modification, it was not possible to obtain the same results showed by Jiang et al. (2020), so we implemented their formulation, which we call MT Z in the following, to provide a comparison with our methods. The results are shown in Table 7 in which we display the instance size n, the time window size (tw), and, for each formulation, the percentage gap, the seconds needed to get to the end of iterations (proven optimality or the time limit of 2 hours), and the number of optima.
One can see that all the proposed formulations provide better results than MTZ, being able to solve to the proven optimality 18 instances compared to the 14 instances solved by MTZ. In particular, all instances of size 20 and 8 out of 10 of size 40 could be solved to optimality with our methods. VFF is the formulation with shorter times on average and the one that presents the best %gap for the instances with up to n = 60 and time windows of value 20. From n = 60, tw = 100, the best method is OCF. The percentage gaps shown for instances with n = 100 are large also because no very good quality solutions could be found. Mancini and Gansterer (2021) propose the vehicle routing with private and shared delivery locations and solve it both exactly and heuristically. We are interested in their exact results. The problem that they solve is similar to the VRPLTW with some modifications. First of all, the objective function includes one more term that aims at minimizing the number of trucks used (multiplied by coefficient α 5 = 1). The new objective function can thus be rewritten as in (55).

Comparison with Mancini and Gansterer
The values of the used parameters are: α 1 = 1, α 2 = 1, Another modification is needed because Mancini and Gansterer consider a fixed service time σ i = 5 for serving a customer i ∈ V C and an una tantum service time σ b = 10 when visiting a locker b ∈ V L . To include that, we can update the service times as follows: modification is about the truck capacity: the authors use an unlimited set of uncapacitated trucks, and thus, we can set Q = n and q i = 1, i ∈ V C . Mancini and Gansterer propose three sets of 10 instances available online (Mancini and Gansterer 2020) with 5 lockers and 25, 50, and 75 customers, respectively, randomly generated in a 10*10 km area. The maximum travel time for each truck is 12 h (thus l 0 = 720), and the time window of a customer can be one of the 12 hours. The travel time is equal to the travel distance that is computed as 3 times the Euclidean distance between two points. In our adaptation, the travel times will also include the service times as explained previously. All customers can be served either with a truck or via a nearby locker; indeed, only customers within a travel time radius of 15 can be served by that locker (R b = 15, b ∈ V L ). Note that we must fix variables y ib = 0 if τ ib > R b , service times excluded. Table 8 reports the results of the MILP proposed by Mancini and Gansterer (MG) and the results of our methods adapted to solve their version of the VRPLTW on their set of instances. We show the number of customers, the number of lockers, and the name of the instances. Moreover, we report the best-known solution (BKS) for each instance. The BKS is bold when it is the proven optimum and in italic when the optimum is shown for the first time in this paper. We also show, for each method, the percentage gap between the BKS and the lower bound at the end of the iterations, the seconds needed to prove optimality if the time limit of 1 hour is not reached, and the number of optima for each set of instances. Note that the %gap of MG is updated (improved) with respect to what is shown in Mancini and Gansterer (2021), because we used the improved BKS values.
Our methods could provide a larger number of proven optima with respect to MG and most of the times with shorter times. In particular, OCF can solve 28 out of 30 instances, obtaining two more 50 customer instances optima and eight new 75 customer instances optima. This is reflected by the average gap of our best formulation, OCF, which is 0.15, with respect to 3.74 of MG. Note that for smaller instances (25 customers) VFF is the fastest method.

Results analysis
In this section, we compare the solutions of the problems previously presented to inspect the impact of the lockers and time windows on the solutions and their costs.
In Figure 2, we compare the value of the solutions of different TSP-based problems on the same instances. We display all values in boxplots, in which the lower and the upper bounds are the first and the third quartile, the middle line shows the median, and the diamond shows the average.
The first comparison we perform is between the classical TSP and the TSPL; the %di f f erence of the costs is computed as 100 · (c T S P − c T S P L )/c T S P . One can see in Fig. 2a that the introduction of lockers can improve the cost of the solutions, because it can simplify the path of the truck, by at most about 10% on the used instances, but on average the improvements on the total cost are between 3.5 and 7.2%, approximately. Moreover, this value seems to decrease with the size of the instance. On the other hand, the routing costs in the TSP are much higher than in the TSPL, with up to 30% difference, and this is due to the fact that more customers are served with lockers and that the routing cost is not the only one considered in the TSPL.
In Figure 2b, we show the percentage difference of the classical TSP with time windows (TSPTW) and of the costs of the TSPLTW computed as 100 · (c T S PT W − c T S P LT W )/c T S PT W , which means that we show the cost improvement that occurs when using lockers in a TSP with time windows. This improvement can reach more than 40%, which shows how the introduction of lockers can help in decoupling the times of customers and couriers and improve the solution cost overall. In fact, the improvement is more evident when the time windows are more restrictive; note, for instance, the very large difference between n20tw20 and n20tw100. The routing costs show an even higher improve-ment that is due to the fact that more customers are served with lockers and thus do not affect the routing component.
In Figure 3, we compare the value of the solutions of different VRP-based problems on the same instances. In Fig. 3a, we show the %di f f erence in costs between the VRPL and the classical VRP, while in Figure 3b, we show the %di f f erence on costs between the VRPTW (the classical VRP with time windows) and the VRPLTW. They present a trend very similar to what was shown in 2a, 2b; however, one can see that in the multiple vehicle case the cost improvement is lower Fig. 2 Increase or decrease in the total costs and routing costs with respect to the several versions of the TSP based problems solved. In black the instances with n = 20 and with the time windows of 20, in gray the instances with n20tw100, in blue n40tw20, in red n40tw100, in orange n60tw20, in green n60tw100, in light blue n100tw20, and in yellow n100tw100

(b)
than that of the single vehicle case, even if still very high. On average, the %di f f erence is between 0 and 11% when considering the VRP and the VRPL and between 3% and 31% when considering the VRPTW and the VRPLTW. This is due to the fact that multiple vehicles can operate in parallel, and this makes it easier to serve customers within their time windows with less costly routing solutions.
In Figure 4, we show the solutions of the TSPs and VRPs on one instance with n = 20, m = 2, and tw = 20. In the graphical representation, the circles are the customers, the triangles are the lockers, and the square is the depot. Note the effect of the time windows on the solutions of the TSPTW and VRPTW: the trucks must travel back and forth, running on costly routes, in order to respect the time windows of the customers. One can see that the solutions of the TSPLTW and VRPLTW take advantage of the lockers and the truck routes are visibly less costly. Moreover, one can see that the TSPL and the VRPL solutions provide an improvement of the TSP and VRP solutions by using lockers to serve customers that are located far from the routes and making the routes effectively shorter.
In Table 9, we show, for all the studied problems, the percentage of lockers used (%lock) and the percentage of customers served via locker (%self) averaged on the number of customers and time window size. One can see that problems with time window constraints use more lockers, and this is also the case of instances with more restrictive time windows because lockers help in desynchronizing by allowing the courier to leave the parcel at the locker before the time window of a customer is opened. For the same reason, the number of lockers that are used is higher in TSPs solutions than in VRPs ones; indeed, the fewer the vehicles, the harder it is to match the time windows. Nevertheless, the number of used lockers increases with the size of the instance.
In Table 10, we show the number of trucks used in the optimal solutions of the several VRPs studied, averaged on the number of customers and the time window size. In the solutions of the problems including time windows, more than twice the number of customers is served via lockers than in the other problems, and this number increases with the increase in the instance size. More restrictive time windows also require more customers to be served via lockers .
The classical VRPTW, having no lockers to be used to serve customers, needs a larger number of vehicles in all cases with respect to the other problems, and this underlines once more the relevance of the use of lockers. On the other hand, the number of vehicles used is minimal when no time window is imposed .

Analysis on the number of lockers
In this section, we analyze the difference in the total costs of the optimal solutions when allowing only a subset of lockers to be open. In particular, we study the instances with n = 40 for all the problems (TSPL, TSPLTW, VRPL, VRPLTW) by allowing from 0 lockers (thus solving the TSP, TSPTW, VRP, VRPTW, respectively) up to 4 available lockers (as in the solutions presented in previous sections) and we analyze the difference in the optimal solutions. In Fig. 5, we analyze the difference in costs: we propose a boxplot that shows, for each problem, the percentage difference between the value of the optimal solution with 0 lockers and the value of the optimal solutions with 1, 2, 3, and 4 lockers computed as follows: %di f f erence = 100 · (z 0 − z i )/z i , being z 0 the solution with 0 lockers and z i the value of the solution with several lockers, i = 1, . . . , 4. Note that for a few instances, the optimal solution was not available and we did not include the %di f f erence for them. One can clearly see that the problems with time windows benefit greatly from the increase in available lockers, more evidently in the single  vehicle case, improving the solutions costs up to more than 40% when using 4 available lockers. On the other hand, the problems without time windows show a limited improvement in the solution costs with up to 10% when using up to 4 lockers. Moreover, the instances with more restrictive time windows benefit more than the others when the number of available lockers is increased. Note that we divided the instances between tw=20 and tw=100 also for the TSPL and the VRPL because they represent two different sets and one can easily compare the results with those of the TSPLTW and the VRPLTW.
In Table 11, we show the percentage number of lockers used (%lock) and the percentage number of customers served via lockers (%sel f ) in the optimal solutions of the studied problems averaged on the number of available lockers (#lock, from 1 to 4) and on the time window size. The problems with time windows make use of a larger number of lockers and serve more customers via lockers, being the solutions of the TSPLTW those showing the highest numbers. The smaller the time windows, the more lockers are used and the more the customers are served via lockers. In the TSPL and the VRPL solutions, one can see that for the used instances the solutions with 3 and 4 lockers serve the same number of customers via lockers, while the more constrained problems keep benefiting from the increase in available lockers.
In Table 12, we show the number of trucks used in the VRPL and the VRPLTW solutions averaged on the instances with n = 40 for number of available lockers from 0 to 4. The number of trucks used in the optimal solution of the VRPL is always two, while it decreases with the increase in available lockers in the optimal solutions of the VRPLTW, because lockers help in decoupling the truck and the time windows, thus allowing to use fewer trucks to perform deliveries while respecting the time windows.

Conclusion
In this paper, we have considered single and multiple vehicle routing problems with lockers and with time windows. We proposed three mixed integer linear programming formulations, one of which needs a branch-and-cut implementation, that are applied to all the problems, and a set of valid inequalities. The proposed one-commodity formulation is the one that performed better when solving problems without time windows, being able to solve all the instances but one with up to 100 customers and 10 lockers for the single vehicle case, and many with up to 60 customers and 6 lockers for the multiple vehicles case. The same formulation provide the smallest gaps on problems with time windows, but the proposed vehicle flow formulation is faster when solving smaller instances. We also improve the literature results on two related problems obtaining new optima and improved best-known solutions. By inspecting the solutions, we can conclude that the introduction of lockers diminishes the total costs and provides simpler routes for the vehicles; this is particularly true when time windows are considered. In such cases, the lockers help in decoupling the couriers and the customers' times allowing less costly routes to be taken. The use of multiple vehicles has a similar, even if less evident, effect .
Possible future research is focused on the improvements of the newly proposed formulations by defining new valid inequalities to include in a branch-and-cut fashion. Moreover, we intend to foster the ideas of this paper to similar applications and problems, including the following features: delivery preference for customers; lockers modeled with cells of different capacities; pickup and delivery; etc. Eventually, we intend to develop dedicated metaheuristic and branchand-price algorithms .