A Robust Optimization Method for Location Selection of Parcel Lockers under Uncertain Demands

: Parcel lockers have continuously growing in popularity as an alternative mode for last-mile delivery services due to their capability of effectively alleviating the risk of a delivery failure, increas-ing the possibility of delivery consolidation, and reducing the number of drop-off sites. However, poorly located of parcel lockers be less efﬁcient. When determining the parcel locker location, inade-quate consideration of uncertain demands can potentially increase the risk of unsatisﬁed demands. To remedy this issue, a robust optimization model is proposed in this paper with consideration of the demand uncertainties, including the large and small parcels to be received and sent. Not only can the collection locations be optimally determined, but so can the number of large and small parcel lockers for each location at the same time under various robust levels. Meanwhile, the sites whose demands are covered by one of the collection locations are also determined by the constraints of acceptable walking distance. A series of numerical experiments has been performed to evaluate the proposed model, with the main focus being on the solution robustness. Since the set of non-linear constraints are transformed into the linear counterparts, the robust solution can be obtained by the existing solvers within a reasonable time with moderate computing power. The experimental results also provide useful guidance for the practical application of the method, as slightly more conservative decision making can secure the solution robustness with only a marginal increase in costs.


Introduction
With the development of information technology and e-commerce, the market share of the B2C shopping mode continues to grow, which brings huge demand for express delivery. Statistics from China Post show that China's express delivery service companies handled 108.3 billion pieces of business in 2021, up 29.92% from 2020 [1]. However, the costliest stage in the supply chain is the last-mile delivery, which accounts for approximately 28% of the total cost [2]. In the traditional door-to-door delivery service, the failed delivery attempts is one of the most expensive resources for the last-mile delivery. In addition, customer security concerns over door-to-door delivery potentially degrade their satisfaction.
As compared with the traditional door-to-door delivery mode, the parcel locker mode is convenient, secure, and accessible 24/7 for customers collecting and sending parcels, while high degree of consolidation and a wide service time window can effectively delivery services extricate from the cost predicament of the last-mile delivery. At present, more than 20 countries around the world have deployed parcel lockers [3]. From 2011 to 2019, DHL added more than 10,000 parcel locker service points in Germany [4,5]. In recent years, the InPost parcel locker service has gained rapid popularity over Europe [4]. In China, this paper, we employed robust optimization to tackle the uncertainty issue related to the delivery demand, as no assumption on the probability distribution of uncertain demand should be made.
The robust optimization method proposed in this paper contributes to knowledge by extending the location optimization method for parcel lockers under uncertain demands. More specifically, we consider two sources of uncertainty for both the receiving and sending parcel demand. In addition, the optimal locations and number of lockers can be optimized simultaneously with various sizes of lockers. After formulating the robust counterpart for the bounded demands, a set of non-linear constraints are transformed into the corresponding linear constraints to accelerate the solving speed.
The remaining sections of this paper are structured as follows. The next section is devoted to formulating the location problem of parcel lockers under a number of assumptions. In Section 3, the models to determine the optimal location of parcel lockers are successively proposed under deterministic and stochastic demands, respectively. Following that, a series of experiments and the corresponding results are described to evaluate the robust optimization model proposed. This paper is concluded with a number of findings, which are summarized in Section 5.

Problem Formulation and Assumptions Problem Description and Model Assumptions
In this study, we present the distribution network as G(I, E), comprising a set of demand sites I and edges E, with the edge weight being the shortest distance between two connected sites. Although all demand sites could be equipped with a number of parcel lockers, the operating cost would be unacceptably high. Therefore, we need to select a set of demand sites to accommodate parcel lockers and determine the number of lockers for each collection location, at the minimal cost and within an acceptable walking distance for those customers whose communities are not equipped with parcel lockers. In other words, the customers can go to the nearest parcel locker site to collect their parcels within the acceptable walking distance. A diagram of this problem is shown in Figure 1. During the planning phase, it is unrealistic to accurately predict the exact demand for each site, but the demand, varying in a range, can be normally estimated with some knowledge of the future trend. Obviously, it rarely happens that all demand sites have the maximum demand, except for in boom seasons. Due to the fact that most parcel lockers are operated independently by a third party, we assume that the delivery task is accomplished by the courier companies. The determination of the number of parcel lockers also needs to consider the operation process of couriers in collecting and uploading parcels. Based on our survey on the courier operation, it is assumed that a courier empties the parcel lockers by collecting the parcels to be sent out and not be taken within a certain period (i.e., 48 h in this study) at first, and then stores the newly arrived parcels when reaching a parcel locker site. Further assumptions are made as follows: (1) The delivery demand, including the number of parcels to be received and sent by customers, of each site fluctuates within a known range.
(2) Customers are always willing to go to the closest collection location of parcel lockers to collect and send parcels, and the distance to the closest parcel lockers must be within an acceptable walking distance.
(3) Each collection location is equipped with two sizes of lockers for large and small parcels, respectively, and the small parcel can also be stored in a large locker if no small lockers are free.
(4) For simplicity, we assume that the cost (including manufacturing, construction and maintenance cost, etc.) and the size of a large parcel locker are δ (δ > 1, δ ∈ Z) times those of a small locker on average. The cost and the size of a small locker are known values (the parameters used for this study are provided in Section 4.1).
(5) A parcel collection location can be equipped with ϑ large parcel lockers or, equivalently, δϑ small ones. (4) For simplicity, we assume that the cost (including manufacturing, construction and maintenance cost, etc.) and the size of a large parcel locker are ( 1, ∈ ) times those of a small locker on average. The cost and the size of a small locker are known values (the parameters used for this study are provided in Section 4.1).
(5) A parcel collection location can be equipped with large parcel lockers or, equivalently, small ones. Before presenting the optimal location models, the variables and parameters involved in the models are summarized in Table 1.

Sets
Descriptions I The set of demand points E The edge set, and each edge associated with a weight being the shortest distance between two connected sites Parameters Descriptions The purchase and maintenance cost of a small locker The land rent cost per unit at the collection location i A positive integer greater than 1 A positive integer The maximum number of lockers allowed at each site A sufficiently large positive number

M
A sufficiently large positive number The distance from demand site j to collection site i The maximum acceptable walking distance The demand for large parcel services at site j The demand for small parcel services at site j The number of large parcels uncollected during service time at site j The number of small parcels uncollected during service time at site j ̅ The average demand for large parcel services at site j The maximum variation in demand variation for large parcel services at site j Before presenting the optimal location models, the variables and parameters involved in the models are summarized in Table 1.  The number of large parcels to be set at site i x i s The number of small parcels to be set at site i y ij A binary variable. If demand site j is assigned to demand site i, y ij = 1, or 0 otherwise

Optimal Location Models
With the problem description and assumptions provided in the previous section, this section presents the optimization model to determine the collection locations and the

Optimal Location Model under Deterministic Demands
The optimization goal is to minimize the total operating costs, including the cost of purchase, maintenance and land rent of the parcel locker network. Obviously, the cost of purchase and maintenance of a parcel locker is only related to the type and number of parcel lockers, which will be the same for different sites. However, this is not the case for the rental cost of the land used to accommodate a set of parcel lockers, since different rental rate may be applied at different areas. Based on the assumptions summarized in the previous section, the objective function can be expressed as (1).
where the first part of the objective function represents the total costs of purchase and maintenance for all parcel lockers, and the second part the total land rental costs. x i b and x i s are decision variables, representing the number of large and small parcel lockers to be provided at site i. δ is a positive integer greater than 1, indicating that the size, purchase and maintenance cost of a large parcel locker is δ times that of a small one. c g i is the land rent of a parcel locker at site i. Since it is normal that several parcel lockers are stacked one above another, the land area occupied by a set of stacked parcel lockers (thereafter, called a unit of parcel lockers) is equal to the land needed to accommodate one single locker. It is assumed that a unit of parcel lockers consisting of ϑ (a positive integer) large parcel lockers or, equivalently, δϑ small parcel lockers and, therefore, x i s + δx i b /(δϑ), represents the number of parcel locker units to be deployed at site i. Note that we round up to the nearest integer when estimating the area to rent for a set of parcel lockers. However, for simplicity, we relax * and then obtain (2).
The constraints associated with the optimization model are shown in (3)-(11).
Constraints (3) and (4) are non-negative constraints on the decision variables x i b and x i s . When x i s + x i b > 0, demand point i is selected as a collection site. Due to the limited area available to accommodate parcel lockers at a site, the total number of large and small parcel lockers cannot be equipped without an up-bound, which is indicated by Constraint (5). In Constraint (6), y ij = 1 indicates that those customers at demand site j need to go to site i to collect or send their parcels, or 0 otherwise. Furthermore, each demand site can only be associated to one collection location, which is expressed by Constraint (7). Constraint (8) states the fact that, if site i is selected as a collection location, the provision of parcel lockers should at least satisfy the demand from site i. Constraint (9) ensures that collection site i, associated with demand point j, has been equipped with parcel lockers. We assume that customers always tend to collect and send parcels at the nearest collection site, and this restriction is realized by Constraint (10). As shown in (11), we impose a constraint that the distance from demand site j to collection location i cannot exceed the maximum acceptable walking distance. Constraint (12a) ensures that the unoccupied large parcel lockers at collection site i can satisfy all demand on the large parcel delivery service from all demand sites associated to this collection site. As small parcels can also be accommodated by the unoccupied space typically reserved for large parcels, the requirements on storing small parcels can be satisfied jointly by the unoccupied small lockers and the large lockers. This is ensured by Constraint (12b).
Note that Constraint (10) contains min function, which can raise difficulty when solving this optimization problem. Therefore, we transform it into the corresponding linear constraint, as shown in (13).
In (13), the sufficiently large positive number M should satisfy max When y ij j(1) = 1, demand site j is allocated to collection site i and "(1)" indicates that site j is the nearest site for the customers from site j to collect and send parcels. Correspondingly, "j(k)" and "j(m)" denote that site i represents the kth and mth nearest sites from site j, respectively.
To prove the equivalence of (13) to (10), Lemma 1 and Lemma 2 are given first.
. . , |I|}, it can be easily con- Lemma 2. If ∑ i∈I y ij = 1 and y ij j(k) = 1, y ij j(s) = 0, s = k, ∀ i, j ∈ I, ∀ k, s ∈ K must hold. This is to say, demand site j can only be allocated to one collection site.
In the following proof, three cases are presented based on Lemma 1 and 2. Case 1: When site I, which is the nearest site to site j, is selected as a collection location and all parcels of the customers at site j are stored in the parcel lockers at site i, i.e., x i b,j(1) + x i s,j(1) ≥ 1, and y ij ≤ 0 holds based on Lemma 1. Therefore, Case 2: When site i, which is kth nearest site to site j, is chosen as an collection location and the customers at site j need to go to site i to receive the delivery services, i.e., , and y ij j(k) = 1. (1) When k = 2 (i.e., s = 1) y ij j(2) = 1, according to Lemma 2, y ij j(1) = 0. Additionally, since is true.
(3) When |I| > h > k Similarly, we can obtain y ij j(h) = 0 based on Lemma 2, with y ij j(k) = 1. Since Case 3: When h = k and h, k = 1, . . . , |I| Suppose x i b,j(k) + x i s,j(k) ≥ 1, x i b,j(s) + x i s,j(s) = 0 (1 ≤ s ≤ k − 1) (i.e., site i, which is kth nearest site to site j, is selected as a site to provide parcel lockers) and y ij j(h) = 1 (i.e., site i is associated to site j, which is hth nearest site to site i), y ij j(k) = 0 holds, as y ij j(h) = 1 and ∑ i∈I y ij = 1. Furthermore, it is apparent that

Demand Uncertainties Considered
This subsection provides a summary of the demand uncertainties considered in this study, and each of them is depicted as a random variable for the delivery demand from each site. Since it is not easy to accurately determine the probability distribution of these random variables during the planning phase, this study assumes that the delivery demands of large and small parcels to be received and sent by customers at each site vary within a known interval. More specifically, it is assumed that the number of large parcels d j b and small parcels d j s to be received by the customers at site j are in the range d j b −d j b , d j b +d j b and d j s −d j s , d j s +d j s respectively, while the large and small parcels to be sent by the customers are estimated to vary in the ranges of r j b −r j b , r j b +r j b and r j s −r j s , r j s +r j s .

Optimization Location Model under Uncertain Demands
Based on the robust optimization approach proposed by Bertsimas and Sim (2004), this study extends the robust optimization model to the parcel locker location optimization problem containing multiple random variables representing the uncertainties in the number of small and large parcels demanded and uncollected. To this end, the constraints for the uncertain demands are formulated as (14a) and (14b), which correspond to (12a) and (12b) for the deterministic demands.
Note that the uncollected parcels at collection site i are split for all sites associated with collection site i for simplicity. That is, ∑ j∈J y ij r j b = x i b and ∑ j∈J y ij r j s = x i s . Constraint (14a) ensures that the number of large parcel lockers provided at site i should be adequate to accommodate the maximum potential demand on the large parcels (including the large parcels to be collected and uncollected) from all demand sites associated to site i (J i denotes the set of demand sites from which customers have to go to site i to collect or send parcels) for given Γ i . On the other hand, the supply of parcel lockers for all demands, including small and larger parcels, is ensured by (14b) in the same way as that in (14a). Among them, Γ i , which can vary in the range of [0, |J i |], is an important parameter to adjust the robustness. More specifically, the larger Γ i is, the smaller the chance the demand constraint will be violated. When Γ i = |J i |, the most robustness can be obtained for the problem. However, it is too conservative, as it hardly happens when the demands from all sites approach the upper bounds at the same time. On the other extreme, the optimization model will deteriorate to the deterministic one when Γ i = 0. If Γ i is not an integer, the integer part (i.e., Γ i = |S i |, S i denotes that the demand at site i is assumed to be the maximum and {S i } is the set containing all sites where the demand reaches the maximum) represents the number of sites which have the maximum demand, while the demand variation at site t i (t i ∈ J i \S i ) is restricted to the fractional times its up-bound i unoccupied (i.e., t i represents that the demand at site i is below the up-bound, and {t i } is the set containing all the sites at which the demand is below the maximum).
, we can re-write it as, with an auxiliary variable z 1 ij (which varies in the range of [0, 1] when y * ij is one and 0 otherwise) (15)- (17) indicate that the variation in the total demand (including the large parcels to be received and sent by customers) at the collection location i is maximized for the demand sites associated with y * i under parameter Γ 1 i . To convert the stochastic constraint into the corresponding deterministic form, we write the dual form for (15)-(17) as follows ) and the association indicator y * i for site i, with an auxiliary variable z 2 ij (which varies in the range of [0, 1] when y * ij is one and 0 otherwise). For convenience, we rearrange it as follows: The dual form of (22)- (24) is Therefore, the location problem of parcel lockers under uncertain demands can be expressed as an integer linear programming (ILP) model, presented as follows.
The above models can be solved by the existing commercial solvers (e.g., CPLEX or Gurobi, etc.). As we coded this model in the MATLAB environment and, therefore, the solver provided for ILP in MATLAB is used in this study.

Numerical Experiments
Before presenting the experiments designed to evaluate the proposed location method for parcel lockers, a set of parameters involved in the simulated experiments are introduced first.

Parameter Settings
In this study, the primary parameters, shown in Table 1, for a unit of 100 parcel lockers provided by a parcel locker management company are adopted as the basis for the numerical experiments. Based on the information shown in Table 2, it can be directly derived that δ is 2, ϑ is 60, the land rent cost c g i of a parcel locker is about 16.44 yuan per day, the purchase and maintenance cost c f of a small locker, and communication cost and electricity charging cost is 0.22 yuan/day. Since the land rent rate at each demand point may not necessarily be the same, we assume that the land rent rates in the experiments varies randomly with the mean value of 16.44 yuan/parcel locker/day. In the experiments, the developed robust optimization model has been evaluated based on the location problems consisting of 50 demand sites, with 150 m as the maximum acceptable walking distance. The demand generation process for each experiment is as follows. The large parcel demand for each site is randomly generated in the range [30, 50], with the demand variation restricted to the range [2,12]. On the other hand, it is assumed that the number of large parcels uncollected at each site varies in the range [4,7], with minimum and maximum variation bounded to 2 and 5, respectively. The same generation process is adopted for the small parcel demand for each site, with the mean demand and variation restricted to the range [50, 150] and [10,40], respectively. Furthermore, the ranges [15,25] and [10,20] are used to generate the mean and variation values in the small parcels which go uncollected.
In each experiment, the time for searching a location solution is restricted to 10 min and the parameters remain unaltered apart from the one whose influence on the solution is examined. All experiments were performed on a x64-PC with Intel Core i7-8550U 1.80 GHz CPU and 8.0 GB RAM.

Experimental Results
It should be noted that, due to the need to analyze the total number of demand points from which the demand fluctuation value is obtained, this paper uses Γ (that is, the sum of all pick-up points Γ i and Γ i varies in the range [0, |J i |]) as an indicator in the experimental analysis.
For ease of explaining the experimental results, we first present an example for the experiment conducted with Γ = 20 (i.e., the summation of Γ i for all sites) for a set of delivery requirements, randomly generated based on the parameters summarized in the previous sub-section. The results obtained are shown in Table 3. Not only the solution of the collection locations, the demand sites associated with each collection site, and the number of large and small parcel lockers, but also unsatisfied demands are listed in Table 3. For example, site 8 has been selected as a collection location which covers only the demand from its own site, and the number of large and small parcel lockers determined as 36 and 75, respectively. Since the actual small parcel demand at site 8 is 76, the unsatisfied demand is 1. Although Γ is set as 20 which indicates only 40% of uncertain demands are considered during the location optimization, most demands are satisfied by the solution generated from our approach in general.
To gain an insight into the converging process, the major operation of the solving process for the above example is briefly summarized here. The problem solving starts by relaxing the ILP problem into a linear programming (LP) problem, which is solved by the dual simplex algorithm. In this example, the optimal objective value obtained for the relaxed problem is 2279.38, which is an initial lower bound and can never be reached for the ILP problem. Subsequently, the cut generation is performed to tighten the LP relaxation of the ILP problem so that the solutions get closer to integers. After obtaining a feasible solution or one which exceeds expectations, the branch-and-bound procedure is invoked with the iterative construction of a set of subproblems which attempt to converge to a solution of the ILP. The number of nodes explored, running time, number of integer solutions found, objective value and relative gap are briefly summarized in Table 4. It is apparent that the objective value decreases as more nodes are explored. Although a quick convergence can be observed at the initial stage after applying the branch-and-bound algorithm, no further reduction in the objective value can be found at the later stage.

Robustness Analysis
The role of parameter Γ is to regulate the robustness degree for the optimal solution, since the amount of uncertain demand is controlled implicitly by Γ, thus resulting in different violating degree. To evaluate the adaptability and robustness of the proposed method under various uncertain demand levels, a series of numerical experiments have been conducted by varying parameter Γ from 0 to 45 (corresponding to 0 to 0.9|J i |) with an interval of 5, and the unsatisfied parcel requirements for the solution under each Γ value have been obtained. Note that the influence of each Γ value on solution has been examined by performing 20 independent runs under different demands randomly generated and the averaged results are presented in Table 5. Relative cost ratio listed in Table 5 denotes the gap between the cost of the optimal solution obtained for a certain Γ value (i.e., Γ > 0) and that for the deterministic demand (i.e., Γ = 0). Similarly, a relatively large parcel locker ratio and relatively small parcel locket ratio represent the difference between the number of large and small parcel locker under a certain Γ value (i.e., Γ > 0) and those under deterministic demand (i.e., Γ = 0).
As shown in Table 4, the unsatisfied numbers of large and small parcel demands under the deterministic situation (i.e., Γ = 0) are 209.5 and 697.0, respectively. When Γ increases to 5, the cost increases by 24.21% (614.81 RMB), but only 11 large parcel demands and 74.5 small parcel demands are not satisfied on average, as compared to those under the deterministic demand. When Γ continuously increases from 5 to 45, the total cost rises slightly, but the satisfied demand diminishes rapidly. Furthermore, it can be seen that all demands can be satisfied when Γ ≥ 25. This implies that it would be too conservative if the worst case (i.e., demands at all sites reach to the maximum at the same time) is taken into account when determining the parcel locker location, though the increase in the total cost is relatively moderate. Although it seems that the value of Γ can be set as half of the number of sites (i.e., Γ = 25) used in this experiment, a slightly over-conservative value is recommended in practice to secure that all demands can be satisfied with high probability.
Therefore, based on the experimental results, the following conclusions can be made. First of all, the robust optimization method proposed in this study can deliver a flexible solution to the location problem of parcel lockers in terms of robustness to uncertain demands. Furthermore, high robustness can be achieved for the location solution without considerably sacrificing the total cost. However, it would be recommended that higher value should be selected for parameter Γ for the cost-insensitive cases. Due to the transform of the nonlinear constraints into the corresponding linear counterparts, the existing commercial solver can be employed to generate a high-quality solution for the parcel locker location problem within a reasonable time with moderate computing power.

Conclusions
The parcel locker delivery mode is an effective solution to the problem of last-mile delivery, which has been promoted and deployed in many countries around the world. A reasonable location is critical to promoting the further development of parcel lockers and better to meet the customer demand. However, no reports on the robust optimization model for the parcel locker location problem with consideration of multiple uncertain demands have been found in the literature. Therefore, the study presented in this paper attempts to extend the robust optimal method to the parcel locker location problem by formulating an optimization model to determine both the collection locations and the number of parcel lockers for each collection, with consideration of demand uncertainties. The optimization model has been evaluated through a series of experiments, with the main focus being on the analysis of the solution robustness. The following findings can be obtained from this study.
(1) The optimization model developed is capable of dealing with a set of uncertainties in the numbers of small and large parcels to be received and sent, respectively.
(2) The collection locations and the number of parcel lockers can be determined simultaneously by the robust optimization model proposed, with all demand sites to be covered by the nearest collection site within acceptable walking distance.
(3) The solution robustness and unsatisfied demand risk can be reconciled by adjusting parameter Γ. Although the experimental findings imply that a moderately high level of Γ is able to provide high robustness for the solution, it is recommended that Γ could take a slightly large value to avoid any risk of unmet demands, since the cost resulted from the improved robustness is not significantly exacerbated.
(4) The robust solution can be obtained by the existing solver within a reasonable timescale with moderate computing power thanks to the linearization of several nonlinear constraints, which considerably simplify the programming model.
However, this study also has several limitations, which should be addressed in future studies. First of all, a sharp increase in the delivery demand during special periods (e.g., COVID-19 pandemic) is not considered as a source of demand uncertainty in this study, as it is not easy to accurately estimate the potential bounds of demand. On the other hand, the predication bias on the long-term trend of the parcel locker demand is not explicitly taken into account when constructing the robust optimization model. In practice, several delivery modes, such as attended self-pick-up points and parcel lockers, may be jointly employed in the last-mile delivery service, a practice which also deserves a further investigation.