Integrated facility location and capacity planning under uncertainty

We address a multi-period facility location problem with two customer segments having distinct service requirements. While customers in one segment receive preferred service, customers in the other segment accept delayed deliveries as long as lateness does not exceed a pre-specified threshold. The objective is to define a schedule for facility deployment and capacity scalability that satisfies all customer demands at minimum cost. Facilities can have their capacities adjusted over the planning horizon through incrementally increasing or reducing the number of modular units they hold. These two features, capacity expansion and capacity contraction, can help substantially improve the flexibility in responding to demand changes. Future customer demands are assumed to be unknown. We propose two different frameworks for planning capacity decisions and present a two-stage stochastic model for each one of them. While in the first model decisions related to capacity scalability are modeled as first-stage decisions, in the second model, capacity adjustments are deferred to the second stage. We develop the extensive forms of the associated stochastic programs for the case of demand uncertainty being captured by a finite set of scenarios. Additional inequalities are proposed to enhance the original formulations. An extensive computational study with randomly generated instances shows that the proposed enhancements are very useful. Specifically, 97.5% of the instances can be solved to optimality in much shorter computing times. Important insights are also provided into the impact of the two different frameworks for planning capacity adjustments on the facility network configuration and its total cost.


Introduction
Location analysis provides a framework for simultaneously finding sites for facilities and assigning spatially distributed demand points to these facilities to optimize some measurable criterion (e.g., cost, profit). Since market conditions and customer service requirements evolve over time, it is also relevant to determine the optimal timing for opening new facilities, in addition to selecting their best locations (Nickel and Saldanha-da-Gama 2019). Furthermore, to be better prepared to market changes (e.g., demand fluctuations), the decision process should also comprise capacity sizing decisions. The latter should specify a schedule for capacity acquisition, capacity expansion, and capacity contraction at the selected locations, in the course of the planning horizon.
Generally, facility location research relies on the assumption that future conditions can be predicted with a reasonable degree of accuracy. This has given rise to the development of deterministic models, which are by far the most common. Even though such models have the advantage of often being amenable to mathematical analysis, the mathematical clarity comes at the cost of decreased realism. In practice, many factors are uncertain (e.g., demands, costs, resource availability), thereby motivating the development of stochastic models (Correia and Saldanha-da-Gama 2019). Ignoring uncertainty entails the risk that once the true nature of the uncertain parameters reveals itself, adjustments in the facility network configuration may be necessary that are time-consuming and require a prohibitively high expense.
In this paper, we consider a variant of the multi-period capacitated facility location problem in which some data are uncertain. We assume that there is a finite set of customers with demand for a single product family over a multi-period finite planning horizon, and a finite set of potential locations for the facilities that will offer the commodity. The goal is to choose facility locations, determine a schedule for facility opening and capacity adjustment (i.e., expansion and contraction), and to decide which customers to serve and from which facilities. A distinctive feature of our problem is that customers are sensitive to delivery lead times. This means that some customers require their demands to be met on time (i.e., these customers impose a zero delivery lead time), whereas the remaining customers accept delayed deliveries as long as lateness does not exceed a pre-specified threshold. This form of customer segmentation can be encountered in a variety of industries, such as retailing (Duran et al. 2008;Li et al. 2015) and aftermarket services (Alvarez et al. 2015;Hung et al. 2012;Wang et al. 2002). Naturally, firms that follow this type of customer segmentation also adopt a price differentiation policy as an incentive for customers to accept longer delivery lead times. Future demands are not known a priori. Instead, they are assumed to be described by a finite discrete random vector, following a probability distribution that is known in advance (e.g., estimated using historical data). Under the assumption that the decision-maker follows a risk neutral strategy, the objective is to minimize the total expected cost. The latter includes fixed costs for opening and operating facilities, fixed costs for capacity scaling, and variable processing, distribution, and tardiness penalty costs. We note that our problem arises in the context of facility sizing decisions being reversible in the medium term. This is the case, for example, when on-demand warehousing is used, which is an emerging trend that allows relatively short commitment horizons for lease contracts (Unnu 2020;Pazour and Unnu 2018). Technological innovations in the chemical process industry also provide flexible production facilities whose capacities are easily redeployed through adding or removing modular units (Becker et al. 2019). Therefore, the three planning levels are integrated through considering strategic decisions (involving facility location), tactical decisions (regarding capacity scalability), and operational decisions (pertaining to demand allocation). Accordingly, facility location and capacity sizing decisions can only be made at selected time periods, whereas demand allocation decisions can be made at any time period of the planning horizon.
A variant of our problem was studied by Correia and Melo (2017), who have shown that this is already a challenging NP-hard problem when all parameters are assumed to be known with certainty. To the best of the authors' knowledge, our work is the first attempt to understand the additional complexity that arises from incorporating stochasticity into the underlying difficult-to-solve deterministic facility location and capacity scalability problem in a multi-period setting. For this reason, we focus on demand uncertainty and defer the consideration of additional uncertain parameters.
Given that demand uncertainty is captured by a finite set of scenarios with known probabilities, we follow a scenario-based approach and formulate two-stage stochastic programs for two different strategies regarding capacity sizing. In both models, the first-stage problem determines the opening schedule for facilities and their initial capacities over the entire planning horizon before uncertainty about demand is revealed. Additionally, operational decisions are captured in the second-stage problem by specifying the commodity flow between the selected locations and the customers. Hence, distribution decisions are made after the demand uncertainty has been resolved. The proposed stochastic programs differ in the way capacity scalability decisions are treated. In the first model, we follow the approach usually taken in the literature (Correia and Saldanha-da-Gama 2019;Govindan et al. 2017;Snyder 2006), and assume that these decisions have a strategic nature. Hence, we include them in the set of ex ante or first-stage decisions. By contrast, in the second model, decisions on the timing and sizing of capacity adjustments are assumed to have a tactical nature. In this case, they become ex post or second-stage decisions. This setting is suitable when, e.g. on-demand warehousing services are available, thus enabling the capacities of operating facilities to be upgraded and downgraded on a medium-term basis (Unnu 2020;Pazour and Unnu 2018).
The contributions of the present work are summarized as follows.
(1) We propose two alternative modeling strategies for a new facility location and capacity planning problem with uncertain demand. (2) For each strategy, we formulate the Deterministic Equivalent Program and develop valid inequalities that significantly improve the polyhedral description of the feasible region. (3) We evaluate the potential benefit obtained from solving the stochastic programs instead of their deterministic counterparts. (4) We assess the validity of the proposed two-stage stochastic mixed-integer linear programming (MILP) models by reporting and analyzing the results of an extensive computational study using a generalpurpose optimization solver. In particular, we provide insights into the impact of considering capacity planning decisions as second-stage (tactical) decisions as opposed to the typical approach of treating them as first-stage (strategic) decisions. Without the support of the two models proposed in this paper it would otherwise be difficult for stakeholders (typically from higher management) to perceive the trade-offs associated with different strategies for facility configuration.
The remainder of this paper is organized as follows. Section 2 provides an overview of the relevant literature. In Sect. 3, we formally describe our problem and present two different stochastic programs along with the corresponding deterministic equivalent forms. Several sets of valid inequalities are developed in Sect. 4 to enhance the original formulations and, therefore, make them more amenable to standard off-the-shelf optimization software. Computational results are reported in Sect. 5 and the proposed formulations are compared using various metrics. Finally, Sect. 6 includes a summary of our findings and gives directions for future research.

Literature review
In this section, we review selected research articles on facility location. Our aim is not to provide an exhaustive review of this field but to analyze the extent to which the features captured by our models have been addressed in the literature. Accordingly, special emphasis will be given to stochastic and multi-period location problems.
Different sources of uncertainty have been considered, either individually or simultaneously, in the literature dedicated to stochastic facility location, with demand uncertainty being the most often considered factor. Early contributions addressing uncertain demand include the works by Louveaux (1986) and Laporte et al. (1994). Over the past decades, the stochastic dimension in facility location modeling has evolved towards different areas, e.g., hub location (Alumur et al. 2012;Contreras et al. 2011;Correia et al. 2018;Rahmati and Neghabi 2021). From an application perspective, recent advances in stochastic modeling have been reported, e.g., in humanitarian logistics (Döyen et al. 2012;Grass et al. 2020;Kınay et al. 2018), and supply chain network design (Ghafarimoghadam et al. 2019;Govindan et al. 2017;Mohamed et al. 2020).
While static (or single-period) stochastic location problems have attracted increasing attention in the literature, there has been limited research on the multi-period stochastic counterpart. Albareda-Sambola et al. (2013) study a multi-period facility location problem with various sources of uncertainty. Marín et al. (2018) propose a stochastic programming modeling framework for a general class of covering location problems and develop a Lagrangian relaxation-based heuristic to tackle large-scale instances. In the context of humanitarian logistics, Kim et al. (2019) investigate the problem of finding optimal locations for drone facilities taking into account that flight distance is uncertain due to battery consumption being affected by weather conditions. Khodaparasti et al. (2018) address the problem of designing a network of nursing homes under uncertain demand and a budget constraint on each time period. Marković et al. (2017) develop a multi-stage stochastic program for the problem of locating law enforcement facilities on a road network to intercept stochastic vehicle flows that try to avoid these facilities. In Mohamed et al. (2020), a multi-stage stochastic program is proposed to locate intermediate facilities in a two-echelon distribution network facing uncertain demand.
Capacity planning decisions are often incorporated into facility location problems (Owen and Daskin 1998;Verter and Dincer 1995). In particular, the gradual adjustment of capacities over a planning horizon is an important driver for improving flexibility within a facility network, a feature that is captured in this paper. Ideally, the provision of capacity should be in line with actual demand needs. Depending on the application context, different strategies can be adopted along the planning horizon, such as capacity expansion (Correia et al. 2013;Cortinhal et al. 2015), capacity expansion combined with capacity reduction (Antunes and Peeters 2001;Wilhelm et al. 2013), relocation of capacities to different locations (Melo et al. 2006), and temporary facility closing and reopening (Dias et al. 2007;Jena et al. 2016). Among these strategies, capacity expansion is the predominant policy in a stochastic setting. Early contributions in this area focus on determining a schedule for capacity acquisition and expansion for a set of facilities already operating at fixed locations (Ahmed and Garcia 2004;Ahmed et al. 2003). Any continuous amount of capacity can be installed. By contrast, Alonso-Ayuso et al. (2003) assume a finite set of capacity expansion levels in a multi-commodity production system. For a multi-echelon logistics network, Aghezzaf (2005) propose a decomposition-based algorithm to determine a facility location and capacity expansion schedule that is robust to the uncertain realization of customer demand. In the problem studied by Hernández et al. (2012), the location and size of prisons are to be determined along with a schedule for capacity upgrading at both existing and new facilities under stochastic demand. Recently, a stochastic modeling framework was proposed for a multi-period hub location problem with uncertain demand ). In addition to finding the location of hub facilities and setting their initial capacities, an expansion plan is also determined. This entails deciding on the time period and the number of additional modular units to be installed in a hub facility. Even though modular capacities are relevant in many contexts, e.g. to represent modular equipment (Delmelle et al. 2014;Gourdin and Klopfenstein 2008), they have been hardly considered in the literature. This feature is captured in our work in conjunction with capacity expansion and capacity reduction. To the best of our knowledge, there is only one other study that has addressed both capacity upgrading and downgrading in a stochastic facility location setting, namely Zhuge et al. (2016). In their two-stage stochastic program, first-stage decisions include locating new distribution centers (DCs) in the first time period and selecting their initial capacities from a set of three options. The second-stage decisions define the timing for capacity expansions and contractions as well as the flow of multiple products. Uncertainty in future demand and the budget available for replacing capacity levels are considered. Small-sized problems using aggregated data from an automotive logistics firm are solved by a Lagrangian-based heuristic. The only other works that also model capacity scalability decisions as second-stage decisions, the former being limited to capacity expansion decisions, are Correia et al. (2018) for a hub location problem and Heckmann (2016) for a facility location problem with disruptive events.
Despite considerable progress on stochastic location modeling, there are still a number of relevant features which have not been addressed so far. The present paper contributes to filling this gap by proposing two stochastic models for location and capacity planning in a multi-period context with service-differentiated customer segments.

Problem statement and stochastic formulations
The stochastic multi-period facility location and capacity planning problem addressed in this paper builds upon a deterministic variant studied in Correia and Melo (2017). For convenience, we restate the assumptions that are common to the deterministic and stochastic settings, and introduce additional information pertaining to the stochastic problem. Subsequently, two stochastic formulations are proposed that differ in the framework under which capacity planning decisions can be made.
A planning horizon comprising a finite number of time periods is considered along with a set of customers with demands for a single commodity in each period. Customers are differentiated according to their service requirements. Specifically, some customers require their demands to be met on time, whereas the remaining customers accept delayed deliveries as long as lateness does not exceed a pre-specified threshold. In this case, an order can be split over multiple time periods for the same customer. Typically, customers belonging to the first segment are willing to pay higher prices to have their demands satisfied in the same period in which they occur. To customers in the second segment, discount prices are often offered in exchange for longer response times. In our modeling framework, this scheme is translated into a tardiness penalty cost that is incurred per unit of backorder and per period. Customer service-differentiation offers greater flexibility to a firm in designing and managing its facility network.
A set of potential locations where facilities can be established is assumed to be available at the beginning of the planning horizon. When a facility is opened, its initial capacity level must also be specified. The latter is expressed by the number of modular units that are installed, all having the same size. The adoption of a base modular unit is common in different application areas, such as in the public sector (Delmelle et al. 2014) and in telecommunications (Gourdin and Klopfenstein 2008). Over the planning horizon, the capacity of an operating facility can be adjusted to increase its responsiveness to changes in the level of customer demand. Two forms of capacity scalability are considered: capacity expansion and capacity contraction. This means that one or several modular units can be added to or removed from a facility. Naturally, these options are mutually exclusive for the same facility at the same time period. A limit on the total number of modular units that a facility can hold at any time is prespecified. Due to the sizeable investment associated with opening facilities, locations cannot be temporarily closed and reopened.
The objective is to find a schedule for the deployment of facilities and the expansion and contraction of their capacities to satisfy customer demands at minimal total cost. Different time scales are considered for decision-making. The capital-intensive decisions involving facility location and capacity planning decisions can be made at the beginning of selected time periods, hereafter termed design periods. All other decisions concerning the processing of the commodity at the open facilities and the allocation of customer demands can be made at any time period. This setting is also adopted by other authors, e.g., Mohamed et al. (2020) and Bashiri et al. (2012). Uncertainty about future demand is described by a finite random vector whose joint probability distribution is assumed to be foreknown. The finitely many possible realizations of this random vector are called scenarios. Specifically, each scenario describes all customer demands over all periods of the planning horizon.
Before detailing the two stochastic programming models, we first introduce the notation that is common to both formulations.

Notation
Sets T Set of time periods along the planning horizon. T L Subset of design time periods at which location and capacity planning decisions can be made (T L ⊂ T ). I Set of potential facility sites. J 0 Set of customers requiring timely deliveries (i.e., their demands must be satisfied in the same time periods in which they occur). J 1 Set of customers that tolerate delayed demand fulfillment. J Set of all customers (J = J 0 ∪ J 1 ; J 0 ∩ J 1 = ∅). Let = 1, resp. max = max{ ∈ T L }, be the first, resp. last, design period in which decisions can be made on opening facilities at potential locations and adjusting their capacities through adding or removing modular units. The size of a facility at time period ∈ T L is the outcome of the design and scalability decisions taken until that period. For < max , this size remains unchanged over all intermediate periods between two consecutive design periods, and ( < ). For = max , the capacity of the facility is maintained until the last period of the planning horizon. Let φ( ) denote the last time period between two consecutive design periods and . This parameter is defined as follows: Deterministic parameters Q Capacity of one modular unit. n i Maximum number of modular units that can be available at location i at any time period (i ∈ I ). ρ j Maximum allowed delay (expressed by the total number of time periods) to satisfy the demand of customer j ( j ∈ J 1 ). F O ik Fixed cost of opening facility i with k modular units at design period (i ∈ I ; k = 1, . . . , n i ; ∈ T L ). F E ik Fixed cost of adding k modular units to the capacity of facility i at design period (i ∈ I ; k = 1, . . . , n i − 1; ∈ T L \{1}). F R ik Fixed cost of removing k modular units from the capacity of facility i at design ik Unit processing cost charged by facility i operating with k modular units at time The customer-specified time lag for demand satisfaction ρ j ( j ∈ J 1 ) imposes that an order placed by customer j for time period t (t ∈ T ) must be filled over periods t, t +1, . . . , t +ρ j . In case t + ρ j > |T |, then the last delivery must occur at period |T |. This condition ensures that demand cannot be carried over to periods beyond the planning horizon. If ρ j = 0 for every customer j ∈ J 1 then our problem reduces to the classical setting in multi-period facility location.
As shown in Appendix B, some cost parameters capture economies of scale, namely all fixed costs (F O ik , F E ik , F R ik , M t ik ) and the variable processing costs (o t ik ). Stochastic parameters: It is assumed that the joint probability distribution of the random vector is foreknown.

A stochastic formulation
Given the problem assumptions, a natural approach to the decision-making process is to develop a two-stage stochastic programming model, where the facility location and capacity planning decisions are made in the first stage, and the remaining processing and distribution decisions are deferred to the second stage. In other words, the first-stage (or ex ante) decisions are associated with the definition of a schedule for location, capacity acquisition, and capacity adjustment to be implemented over the entire planning horizon, before a realization of the demand becomes known. This results from the strategic nature assumed for these decisions as they have a long-term impact. The following binary variables represent the first-stage decisions.
z ik = 1 if a new facility is established with k modular units at potential location i at time period , 0 otherwise (i ∈ I ; k = 1, . . . , n i ; ∈ T L ).
(2) e ik = 1 if the capacity of a facility in location i is expanded with k modular units at r ik = 1 if the capacity of a facility in location i is reduced by removing k modular After customer demand is observed, additional decisions are made by specifying the quantities of product to be processed at the operating facilities and distributed to customers over the planning horizon. Therefore, second-stage (or ex post) decisions are represented by the following continuous variables: y tt i j (ξ ) : Amount of product distributed from facility i to customer j at time period t to (partially) satisfy demand of period t (i ∈ I ; j ∈ J 1 ; t ∈ T ; t = t, t + 1, . . . , min{t + ρ j , |T |}).
(7) w t ik (ξ ) : Total quantity of product handled by facility i operating with k modular units at time period t (i ∈ I ; k = 1, . . . , n i ; t ∈ T ).
The implicit representation of the two-stage stochastic program is as follows: (z, e, r , v, ξ)] denoting the recourse function, which represents the expected value of the second-stage problem. The optimal value of this problem, Q (z, e, r , v, ξ), is given by In the first-stage problem (9)-(16), the objective is to minimize the sum of fixed costs and expected processing and demand allocation costs. The fixed costs account for facility location, capacity expansion, capacity contraction, and facility operating costs. By including the expected total cost in (9), a neutral approach to risk is assumed, which is a common criterion to address risk in decision-making problems, often leading to computationally tractable models (Birge and Louveaux 2011). Constraints (10) ensure that at most one facility can be opened at a potential location over the planning horizon. Equalities (11) impose that facilities can only be operated provided they have previously been established. Constraints (12) guarantee that an open facility can be subjected to at most one type of capacity adjustment (either expansion or contraction) at each design period. The number of available modular units at open facilities is defined by equalities (13) for each design period. This number is the outcome of the sizing choice made when the facility was opened, and the number of modular units that were added or removed afterward. Finally, constraints (14)-(16) set the binary conditions on the first-stage (strategic) variables.
Every particular realization of customer demand ξ of yields a second-stage LP-model (17)-(24). Observe that in this model, the design variables z ik , v ik , e ik , and r ik take on fixed values. The objective function (17) minimizes the sum of processing costs at operating facilities, distribution costs to customers, and tardiness penalty costs for orders delivered with delay to customer segment J 1 . Demand satisfaction is imposed by constraints (18) and (19). Capacity constraints are enforced by inequalities (20). Equalities (21) state that the product flow leaving a facility at a given time period is split into deliveries to priority customers and to customers that tolerate late shipments. Non-negativity conditions on the second-stage (operational) variables are set by constraints (22)-(24). Observe that we are facing a stochastic problem with fixed recourse, since the coefficients of the second-stage variables (i.e., the unit cost parameters in (17) and the capacity of a modular unit in (20)) are all known in advance.
Recall that we assume that the random variable ξ has a finite support, i.e., ξ is defined by a finite probability distribution. Accordingly, let s ∈ S be the index of the possible realizations of ξ (called scenarios), where S is a finite set, and let π s represent the associated (positive) probabilities such that s∈S π s = 1. In this case, E [Q(z, e, r , v, ξ)] = s∈S π s Q(z, e, r , v, ξ s ). It follows that we can rewrite formulation (9)-(16) as a MILP model by defining scenario-indexed demands and recourse variables. The former are given by d t js (t ∈ T ; j ∈ J ; s ∈ S), while the latter replace the second-stage variables (6)-(8) and are represented by y tt i js : Amount of product distributed from facility i to customer j at time period t to (partially) satisfy demand of period t under scenario Using the above redefinition of the recourse variables, we obtain the so-called deterministic equivalent MILP model, or extensive form, hereafter denoted (P DE−1 ).
We note that the extensive form of the stochastic program becomes very large even when the random variable ξ has a moderate number of possible realizations. Since all the strategic decisions made in the first stage do not depend on the realization of the random demand process, and so the location and capacity scalability schedule determined for the entire planning horizon is the same for all scenarios, formulation (P DE−1 ) satisfies the non-anticipativity principle. For this reason, (P DE−1 ) will be called the scenario-independent location and capacity scalability model.

A stochastic formulation for an alternative strategy
Discrepancies between the capacity of a firm and the demands of its customers result in inefficiencies, either in underutilized resources or dissatisfied customers. The ability of a firm to adjust the capacity of its resources in response to varying customer demand is one of the key factors in decreasing these discrepancies and achieving competitiveness. Therefore, an alternative strategy to the settings considered in the previous section is to assume that capacity adjustment decisions have a tactical nature, and thus can be deferred to the secondstage problem. In this strategy, the first-stage problem consists of defining a schedule for locating facilities and setting their initial capacities over the planning horizon. Accordingly, the strategic facility location decisions are the here-and-now decisions. The second-stage problem addresses the wait-and-see (tactical and operational) decisions by prescribing a scheme for adjusting the capacities of operating facilities (through expansion or contraction) at the design periods, and specifying the product flows from facilities to customers over all time periods. In contrast to the approach presented in Sect. 3.2, the capacity adjustment measures depend now on the realization of the uncertain customer demand. As mentioned earlier, the resulting two-stage stochastic program is suitable for those cases in which sizing decisions can be reverted in the medium-term. This arises, for example, in the context of leasing or renting space and equipment. On-demand warehousing is an emerging trend that enables firms to adapt quickly to variable demand and cost conditions through deploying storage space at different locations, for different volumes, in a dynamic fashion (Unnu 2020;Pazour and Unnu 2018). In particular, e-commerce retailers profit from this pay-per-use model. Another setting that also allows flexible facility reconfiguration arises in the chemical process industry through modular plants whose capacities can be expanded, contracted or even relocated over a short planning horizon (Becker et al. 2019). Malladi et al. (2020) and Faugère et al. (2020) also report on different logistics systems that rely on modular production capacity that can be redeployable in the short term.
To formulate a two-stage stochastic program for this alternative strategy, the binary location variables (2) and the continuous flow variables (6), (7), and (8) are used. Moreover, the binary variables associated with the sizing decisions, i.e., (3), (4), and (5), are replaced by e ik (ξ ), r ik (ξ ), and v ik (ξ ), respectively. Each one of these new variables depends on the random vector ξ ruling future customer demands. The total expected cost to be minimized is given by (10) and (14). The optimal value Q (z, ξ) of the second-stage problem is determined by Due to capacity adjustment decisions being now recourse decisions, the new constraints (37)-(39) are the counterpart of constraints (11)-(13). Likewise, the capacity constraints (40) replace inequalities (20).
By imposing the same conditions on the random demand vector ξ as in Sect. 3.2, we can rewrite the recourse function as a function of all scenarios, i.e., E [Q (z, ξ)] = s∈S π s Q (z, ξ s ). Accordingly, the set of scenario-dependent flow variables (25)-(27) is extended with the new scenario-indexed capacity scalability variables, v iks , e iks , and r iks . The definition of these variables and the resulting deterministic equivalent MILP formulation, hereafter denoted (P DE−2 ), are detailed in Appendix A. Since all capacity adjustment decisions depend on the realizations of the random parameters, (P DE−2 ) will be called the scenario-dependent capacity scalability model. (P DE−2 ) defines a large-scale MILP model with a significant number of binary variables and constraints, even for a small-sized scenario set.

Additional inequalities
In an attempt to improve the polyhedral description of the set of feasible solutions of models (P DE−1 ) and (P DE−2 ), we develop in this section several sets of additional inequalities. The proposed enhancements may contribute to more problem instances and with larger sizes being solved to (near-)optimality with general-purpose optimization software within reasonable computing time.

Inequalities for model (P DE−1 )
Recall from (28)-(35) that in the extensive form of the deterministic equivalent model (P DE−1 ), all facility location and capacity-related decisions are made in the first stage, before uncertainty about future demand is disclosed. Next, we successively show how the information conveyed by the various scenarios can be used in this case to derive a lower bound on the total number of facilities that must be open at a particular design period.
Let us assume, without loss of generality, that the length of the planning horizon is a multiple of the number of design periods, and let us define σ = |T |/|T L |. Hence, parameter σ is integer-valued and gives the total number of periods between two consecutive design periods. Since all demand requirements of the preferred customer segment J 0 must be serviced without delay in every scenario, they are used to set minimum capacity requirements between two consecutive design periods. The latter are given by for every ∈ T L and φ( ) defined in (1). Hence, the following inequalities must hold for every design period: Dividing the above inequalities by the size of a modular unit, Q, and taking into account that the left-hand side must be integer-valued, it follows that Furthermore, for a particular design period ( ∈ T L ), let D denote the largest minimum demand quantity over all scenarios that has to be covered from period through period φ( ).
The above expression includes all orders from customer segment J 0 as well as the minimum quantity of demand that must be delivered to customers J 1 accepting late shipments. Since demand cannot be lost, in the last time interval (which comprises all periods from max through |T |), the orders of both customer segments must be completely satisfied. Accordingly, the following inequalities ensure that facilities open at time period (and therefore at all subsequent periods up to period φ( )) must have sufficient capacity to serve the demand D : Dividing these inequalities, first by Q and then by the time lag σ , we obtain Combining inequalities (43) and (44), and defining By applying the Chvátal-Gomory rounding procedure a finite number of times p ( p ≥ 1) to the above inequalities (Wolsey 1998), we obtain i∈I The reasoning followed to derive (45) is itself a proof of the following result:

Proposition 1 The set of inequalities (45) is valid for formulation (P DE−1 ).
The enhanced formulation is called hereafter (P + DE−1 ) for the scenario-independent location and capacity scalability model.

Inequalities for model (P DE−2 )
For the extensive form of the deterministic model given by (52)-(58) (see Appendix A), we develop additional inequalities that serve a purpose similar to (45). In this case, the variables v iks (cf. (51)) representing the number of modular units available in a location at time period ( ∈ T L ) are scenario-indexed. Hence, the minimum quantity of demand that must be satisfied in each period between two consecutive design periods is scenario dependent and defined by Given a scenario s and a design period , it follows that the minimum number of open facilities must satisfy the following conditions: The above inequalities are the counterpart of (43) for the scenario-dependent capacity scalability model.
Additionally, let s denote the minimum total demand requirements to be satisfied in scenario s for all periods from through φ( ), with ∈ T L .
Therefore, the following inequalities must also hold and using inequalities (46) and (47) yields the following relations: where the lower bound M s is defined as the maximum of the two parameters on the righthand sides of inequalities (46) and (47). Finally, we apply once again the Chvátal-Gomory rounding method to the above inequalities and obtain the following result: is valid for formulation (P DE−2 ).
The steps presented above serve as proof for this proposition. If for a particular design period and scenario s there exist several inequalities (48), all having the same right-hand side for different values of p, then we only need to consider the inequality with the strongest left-hand side. A similar remark also applies to (45). The enhanced formulation for the scenario-dependent capacity scalability case is called hereafter (P + DE−2 ).

Computational study
In this section, we present the results of an extensive computational study. The computational experiments were guided by five objectives, namely: (i) to evaluate the usefulness of using a general-purpose MILP solver to identify optimal or near-optimal solutions to the original deterministic equivalent models within reasonable computing time; (ii) to analyze the effectiveness of the proposed additional inequalities; (iii) to compare the two strategies for deciding on the timing and sizing of capacity adjustments; (iv) to discuss relevant insights into the characteristics of the (near-)optimal solutions identified; and (v) to assess the benefits of using a stochastic programming approach for the problem under study.

Characteristics of test instances
Since benchmark data sets are not available for the problem at hand, a set of instances was randomly generated. The size of each instance is mainly dictated by the length of the planning horizon and the total number of customers according to the choices specified in Table 1. Planning horizons with 12 and 24 time periods are considered. In both cases, there are three design periods for making location and capacity scalability decisions. Specifically, for |T | = 12 (resp. |T | = 24), these opportunities occur at the beginning of periods 1, 5, and 9 (resp. 1, 9, and 17). Two different sizes for the customer set are assumed and three different partitions of J are considered. For example, a set of 40 customers can be partitioned into |J 0 | = 10 and |J 1 | = 30, |J 0 | = |J 1 | = 20, and |J 0 | = 30 and |J 1 | = 10. In all cases, customers in segment J 1 specify the same maximum number of periods for late deliveries, i.e., ρ j = ρ for all j ∈ J 1 . Instances with ρ = 1, ρ = 2, and ρ = 3 were generated.
In total, five scenarios are considered that reflect different variations in demand. This choice is mainly driven by the very large size of our models, as displayed in Table 2, which is impacted by the number of periods and the number of customers selected. This difficulty is reflected in other similar research in stochastic facility location as well, thus limiting the total number of scenarios (Aghezzaf (2005) 2014) for a stochastic fixed-charge transportation problem, we have combined these stochastic distributions in different ways to create five scenarios. In scenario 1, all customer demands are generated from interval (i) in the first time period; in scenario 2, all demand levels are taken from interval (ii) in the first period; in scenario 3, interval (iii) is used to generate all demands in period 1; in scenarios 4 and 5, a demand level is initially chosen at random for each customer, among the three intervals. This sets the choice for period t = 1. Each interval has an equal probability of being selected. In every scenario s, the demands of customer j in all periods t > 1 are calculated as β t js d t−1 js with β t js ∈ U[0.8, 1.2]. This scheme creates demand fluctuations over the planning horizon that will call for adjustments in the capacity of the facilities. We notice that these scenarios differ significantly from each other, originating a high demand variability, and thus making them a suitable choice in a setting with uncertain demand.
The capacity of a modular unit is defined as a function of the demand scenarios, namely , with n i = 4 for every i ∈ I . Appendix B provides further details on the procedure for randomly generating all cost parameters. Five instances were randomly generated for each combination of the parameters shown in Table 1, giving rise to a total of 180 different instances. Table 2 summarizes the sizes of these instances according to the formulation chosen. In particular, the number of binary variables increases significantly in the scenario-dependent capacity scalability model, namely more than 360%. Together with the considerable number of continuous variables, especially for the instances with 40 customers, we notice that the models' dimensions are challenging in a stochastic MILP context. With the valid inequalities proposed in Sect. 4, the total number of constraints slightly increases in both enhanced formulations.

Numerical results
Formulations (P DE−1 ) and (P DE−2 ) along with their enhancements were coded in C++ using IBM ILOG Concert Technology and run with IBM ILOG CPLEX 12.6. All experiments were performed on a workstation with a multi-core Intel Xeon E5-2650V3 processor (2.3 GHz, 10 cores), 32 GB RAM, and running the Ubuntu operating system (64-bit). Due to the strategic nature of the problem that we study, fast solution times are not of paramount importance. Therefore, a limit of 10 hours of CPU time was set for each solver run. Furthermore, CPLEX was used with default settings under a deterministic parallel mode.
Regarding goal (i) listed at the beginning of Sect. 5, we evaluate the effectiveness of using CPLEX to identify (near-)optimal solutions to the test instances for the proposed original formulations, (P DE−1 ) and (P DE−2 ). To do so, we compare the optimality gaps and the computing times reported by CPLEX. Table 3 summarizes the relevant information according to the size of the customer set, while Table 4 provides information from a different perspective, namely for varying values for the maximum delivery delay tolerated by customer Table 2 Average size of the test instances for formulations segment J 1 . Column 1 gives the parameter selected (|J | in Table 3 and ρ in Table 4). Columns 2-5 report the number of instances solved to optimality (# opt sol.) and the number of instances not solved to proven optimality within the specified time limit (# non-opt sol.) for all models (i.e., the original formulations and their enhancements). For those instances not solved optimally, the minimum, average, and maximum optimality gaps achieved by CPLEX after 10 h of computing time are given for each formulation in columns 7-10 (MIP gap (%) = (z U B − z L B )/z U B × 100%, with z U B denoting the objective value of the best feasible solution available and z L B representing the best lower bound identified during the branchand-cut procedure). Columns 11-14 display the minimum, average, and maximum computing times, in seconds (MIP CPU (s)). The last row refers to all 180 instances. The best (average) values with respect to the evaluation criteria are given in bold type. Table 3 shows that the capability of CPLEX to identify an optimal solution within the given time limit is affected by the formulation that is used. Regarding formulations (P DE−1 ) and (P DE−2 ), optimality is achieved in all test instances with 20 customers, but deteriorates as the total number of customers grows to 40. Not surprisingly, the larger formulation (P DE−2 ) poses an additional challenge to CPLEX. In particular, for this formulation, optimal solutions are only available for less than half of the instances with 40 customers (44.4% or 40/90), whereas 80% (72/90) of the instances with the same number of customers could be solved to optimality with model (P DE−1 ). These findings are also supported by the growth observed in the computing times. Interestingly, for those instances that could not be solved to proven optimality, the associated gaps are, on average, relatively small (cf. columns 7 and 9 in Table 3). The challenges posed by formulation (P DE−2 ) result in slightly larger gaps compared to formulation (P DE−1 ). Table 3 also reveals the positive effect of the proposed valid inequalities (45) and (48), an aspect that was mentioned at the beginning of Sect. 5 as goal (ii) to be investigated. Specifically, the additional constraints prove to be very useful since they allow the identification of optimal solutions for a significantly larger number of problem instances in shorter computing times. Adding inequalities (45) to formulation (P DE−1 ) results in obtaining optimal solutions to all the instances. Regarding the enhanced model (P + DE−2 ), and despite the greater size of this formulation, 95% (171/180) of all instances are solved to optimality. In this case, the valid inequalities (48) prove to be very effective for the more challenging instances with 40 customers (90% or 81/90 optimal solutions become available). Table 4 shows that increasing the maximum lead time for customers accepting late shipments results in slightly more challenging problems for CPLEX. This fact is evidenced by the larger number of instances not solved to proven optimality as the value of ρ grows. We remark that increasing ρ brings more flexibility to designing the facility network and adjusting the capacities of the operating facilities over the planning horizon. Interestingly, the quality of the best solutions returned by CPLEX upon reaching the time limit does not seem to be considerably affected by parameter ρ. The enhanced formulations clearly succeed in solving more instances to optimality and in providing lower optimality gaps under a substantial reduction of the computational burden.
As evidenced by the significantly shorter computing times of the tighter formulations (P + DE−1 ) and (P + DE−2 ) (cf. last column in Tables 3, 4), the exploration of the branch-and-cut tree generated by CPLEX is less time consuming, despite the greater number of constraints. However, additional computational effort may be required to solve the associated linear relaxations. This feature is present in solving the LP-relaxation of model (P + DE−2 ) but not for the LP-relaxation of model (P + DE−1 ), as displayed in Tables 5 and 6. These tables report in the last four columns the minimum, average, and maximum solution times of the     (Table 5) and parameter ρ (Table 6). In addition, columns 3-6 show the minimum, average, and maximum integrality gaps (LP gap (%)). The latter give the relative percent deviation between the objective value of the best feasible solution available (z U B ) and the lower bound provided by the linear relaxation. The best average values over all 180 instances are highlighted in bold type in the last row of these tables. Tables 5 and 6 indicate that the linear relaxations of the enhanced formulations provide substantially better lower bounds, and as a result lower integrality gaps, as opposed to the original models. Even though the number of available optimal solutions decreases with a growing size of J , and accordingly the quality of z U B declines, the average integrality gap becomes significantly smaller when inequalities (45) and (48) are added to the associated formulations. By contrast, the magnitude of the maximum allowed tardiness for orders in customer segment J 1 does not seem to have a major effect on the average integrality gap in any of the four formulations, as shown in Table 6.
Finally, we note that good integrality gaps such as the ones we have obtained with the enhanced formulations can be very useful in evaluating the quality of a feasible solution produced, for example, by a tailored heuristic method. In particular, this feature is relevant for problems with large customer sets since in this case it becomes increasingly difficult for a state-of-the-art optimization solver to achieve optimality within acceptable time. From a computational viewpoint, the optimal solution to the linear relaxation is inexpensive (< 200 seconds for the enhanced models, see Tables 5 and 6).

Deployment of facility location and capacity scalability strategies
In this section, we focus on goals (iii) and (iv) stated at the beginning of Sect. 5, and compare the characteristics of the two strategies proposed for locating facilities and adjusting their capacities over the planning horizon. For this purpose, Table 7 gives the minimum, average, and maximum relative deviation of the objective values of the best solutions reported by CPLEX for the enhanced formulations. Specifically, columns 2-4 (resp. 6-8) display the value of (z + DE−1 − z + DE−2 )/z + DE−1 × 100% according to the size of the customer set (resp. the maximum delivery delay). We denote by z + DE−1 and z + DE−2 the best objective values identified by CPLEX to formulation (P + DE−1 ) and (P + DE−2 ), respectively, within the given time limit. Each row associated with the size of a customer set aggregates the results for 90 instances, while each row associated with parameter ρ aggregates the results for 60 instances. Since the relative deviations presented in this table are all positive, they indicate that cost savings can be achieved with the scenario-dependent capacity scalability model (strategy 2) over the scenario-independent location and scalability model (strategy 1). This finding suggests that modeling capacity adjustment decisions as second-stage variables in formulation (P + DE−2 ) allows a more agile response to temporal changes in demand. In other words, these decisions can be made taking into account the realization of demand. Therefore, in general, the resulting schedule for facility opening and capacity adjustment is cheaper compared to the case in which all design and capacity planning decisions must be made in the first stage. In the latter case, such decisions do not depend on the particular realization of the random demand process.
Over all instances, strategy 2 allows for up to 10.92% reduction in total cost, with an average reduction of 4.5%. Even when an instance is not solved to optimality with formulation (P + DE−2 ) and an optimal solution is available to formulation (P + DE−1 ), strategy 2 always Table 5 Integrality gaps and solution times of the LP relaxation of all formulations according to the total number of customers  has lower total cost than strategy 1. This particular case occurs in 9 instances (5% of the instances).
To further analyze the differences between the two strategies, Table 8 provides additional information about the selected facilities for each combination of potential facilities (|I |), number of customers (|J |), and value of ρ (columns 1-3). Average results for strategy 1, resp. strategy 2, are summarized in columns 4-7, resp. columns 8-11. Since adjustments in the number of modular units at open facilities are second-stage decisions in strategy 2, the average values shown in columns 8-10 are taken over all scenarios. The average delayed demand is also obtained over all scenarios for both strategies (columns 7 and 11).
Not surprisingly, the total number of open facilities is, on average, similar in both strategies. This is explained by facility opening decisions being modeled as first-stage variables in both approaches. Furthermore, decisions on upgrading and downgrading the capacities of operating facilities are driven by the type of strategy used. When such choices are included in the first-stage (strategy 1), there is a tendency to provide more capacity through installing additional modular units. Capacity contraction hardly ever occurs. By contrast, capacity adjustments are more noticeable in strategy 2, especially with regard to removing modular units. The additional cost incurred by reducing the capacity of a facility is offset by a lower maintenance cost for that facility. As expected, strategy 2 offers greater flexibility to respond to demand changes over the planning horizon, thereby yielding network configurations with lower total cost (cf. Table 7).
In general, and independently of the strategy followed to determine a schedule for facility opening and capacity adjustment, the percentage of late deliveries increases with the maximum allowed delay. This feature is more evident in the instances with 20 customers.

The relevance of the stochastic setting
This section focuses on goal (v) of the computational study, namely on the evaluation of the benefits of the stochastic programming approach. This is a relevant issue because of the additional complexity posed by the stochastic setting compared to its deterministic counterpart. Even though there are no standard 'robust measures' for assessing the relevance of using a stochastic model instead of a (simpler) deterministic model (Birge and Louveaux 2011), two metrics are often used. The first is the expected value of perfect information (EVPI), while the second is the value of the stochastic solution (VSS). Many studies in stochastic facility location have used these metrics Diglio et al. 2020;Heckmann 2016;Hinojosa et al. 2014;Marín et al. 2018;Moreno et al. 2018). We have calculated the EVPI and VSS to measure the benefit of working with the stochastic programs (P + DE−1 ) and (P + DE−2 ). As discussed by Birge and Louveaux (2011), EVPI and VSS are distinct indicators that assess different types of uncertainty. The EVPI represents the maximum amount that the decision-maker would be willing to pay in return for complete information about the future. It is defined as the difference between the optimal value of the stochastic program (in our case, z + DE−1 or z + DE−2 ) and the optimal value of the so-called 'wait-and-see' (WS) problem. For each demand scenario s, formulation (P + DE−1 ) (or (P + DE−2 )) is solved under the assumption that the scenario occurs with probability 1. In this way, we obtain the optimal decisions that should be made if it were possible to fully predict future demands. Let z + DE−1,s (resp. z + DE−2,s ) denote the optimal value of the deterministic problem associated with a particular scenario s ∈ S and strategy 1 (resp. strategy 2). It follows that the optimal value of the WS problem is given by s∈S π s z + DE−1,s for strategy 1. A similar calculation is performed for strategy 2.
The VSS measures the difference between the optimal objective value of the so-called 'expected value problem' (EEV) and the optimal objective value of the stochastic program. The former is obtained by creating a single scenario in which all random variables are replaced by their expected values, and then determining the optimal values of the first-stage variables for that single scenario. Next, the expected cost over all scenarios using those first-stage variable values is calculated. In the case of our problem, the average scenario would be constructed by taking the average demand of every customer in each time period, d t j = s∈S π s d t js ( j ∈ J ; t ∈ T ). Unfortunately, the associated values of the first-stage variables may have no feasible completion in the second-stage, thus resulting in an infeasible stochastic problem, either for model (P + DE−1 ) or model (P + DE−2 ), or even both. To overcome this difficulty, a 'reference scenario' is often used (Birge and Louveaux 2011). Typically, the worst-case scenario, i.e. the scenario with the highest total demand level for all periods is adopted. Due to the features of our data, it is not possible to identify such scenario (recall that the demand of every customer may decline or grow between two consecutive periods, and in scenarios 4 and 5, a random choice of the interval from which demand is generated takes place for each customer and t = 1). Hence, a so-called 'modified reference scenario' was determined by taking d t j = max s∈S d t js for every j ∈ J and t ∈ T . Summarizing, the EVPI measures the value of knowing the future with certainty, while the VSS indicates how good the optimal solution to the EEV is as an approximation to the optimal solution to the stochastic program.
Tables 9 and 10 report the average relative values of EVPI (column 5) and modified VSS (column 7) for strategy 1 and strategy 2, respectively. The relative values are calculated by dividing each metric by the optimal objective value of the associated stochastic problem. Column 4 indicates the number of instances for which optimal solutions are available and thus were used for calculating these average relative gaps. In addition, columns 6 and 8 give the average total computing time required to solve the WS and EEV problems to optimality, respectively. The EVPI and modified VSS gaps reported in the tables indicate that it is appropriate to consider a stochastic approach rather than its deterministic counterpart. For the scenario-independent location and capacity scalability model (strategy 1), the relative EVPI gap ranges from 12.59% to 28.83%, with an average of 21.12%. Concerning the scenario-dependent capacity scalability model (strategy 2), slightly smaller EVPI gaps are obtained (average: 17.35%) that vary between 9.83% and 25.83%. In both cases, these relatively large values reveal that having access to perfect information about future customer demands would be quite valuable. This is an indication of a high degree of uncertainty, which, in turn, demonstrates the usefulness of considering a stochastic approach. Interestingly, the average EVPI gap does not seem to be greatly affected neither by the size of the customer set nor by the maximum delivery delay tolerated by customer segment J 1 . For the modified VSS gap, we observe that it gradually increases with growing value of ρ. Recall that this metric measures how far away the value of the optimal solution of a deterministic simplified version of the problem is from the optimal value of a stochastic problem. By increasing ρ, more opportunities arise to satisfy the demand of some customers with delay, which, in turn, enable a lower investment in designing the facility network (e.g., through reducing the need to establish more facilities). This feature becomes even more relevant in a stochastic setting, whereas the quality of the optimal solution to the deterministic EEV problem deteriorates. In other words, the approximation of the stochastic model by one where expected values are used is poor. Finally, and as expected, the computational burden required to optimally solve the WS and EEV problems is considerably lower than solving the stochastic programming models.

Conclusions
We have studied a novel multi-period facility location and capacity planning problem under uncertainty taking into account the sensitivity of each individual customer to delivery lead times. In addition to determining a schedule for opening facilities and setting their initial number of modules, this new problem also allows for two different ways to adjust capacity over the planning horizon. Specifically, capacities can be increased incrementally by installing additional modules or decreased by removing existing modular units. For the case of uncertain future customer demand, we have proposed two stochastic models that differ in the framework in which capacity adjustments can be planned. The first two-stage stochastic model follows the natural approach of defining all decisions related to the design of the facility network as first-stage decisions. Demand allocation decisions are deferred to the second stage. This model is suitable to real-world settings where opening facilities and sizing their capacities are time-consuming and require substantial capital investment. The second model presents an alternative strategy which is particularly relevant when capacity adjustment decisions have a tactical nature, e.g., when short commitment horizons for lease or rental contracts are offered. In this case, first-stage decisions include the selection of locations to open facilities and setting their initial capacities. The phased capacity expansion and contraction of facilities along with demand allocation represent the recourse decisions. Naturally, it is up to the decision-maker to decide which of the approaches should be followed depending on the problem at hand. Under the assumption that demand uncertainty can be captured by a finite set of scenarios, we have presented the extensive form of the deterministic equivalent for the two stochastic programs. Additionally, valid inequalities were developed to enhance the resulting large-scale formulations.
Numerical experiments with randomly generated instances indicate that the proposed additional inequalities substantially increase the capability of an off-the-shelf MILP solver to identify optimal solutions in much shorter computing times. Moreover, improved solution quality is also achieved for instances for which the specified time limit is reached without guaranteed optimality. Deferring the decisions on capacity expansion and capacity contraction to the second-stage problem results in network configurations with noticeable lower total cost compared to the approach where such decisions belong to the first-stage problem. This cost advantage arises due to the possibility of adapting the network configuration to particular demand realizations. However, this approach leads to considerably larger models, which, in turn, require higher computational burden. Our comparative analysis enables the decision-maker to understand the opportunities that arise from different strategies and helps raise awareness of the potential benefits of conducting this type of study. Despite the average number of open facilities being similar in the two approaches, decisions on upgrading or downgrading the capacities of operating facilities differ. In particular, greater flexibility to respond to demand fluctuations is displayed by the second strategy. Thus, a firm having access to e.g. on-demand warehousing services will be aware of the benefits to be gained. Finally, we have also shown that it is worth adopting a modeling framework which captures the uncertainty in future customer demand despite the additional complexity.
From a practical viewpoint, the NP-hard problems that we study pose significant challenges as they involve a very large number of variables and constraints. Hence, obtaining an optimal solution within reasonable computing time with a general-purpose optimization solver becomes extremely difficult for large-sized problems. This characteristic is also encountered in many other multi-period facility location problems under uncertainty (e.g., Correia et al. (2018)), and calls for the development of specially tailored decomposition approaches and heuristic methods. In this case, the enhancements proposed in Sect. 4 would also be helpful as they sharpen the linear relaxation. The LP-bounds would also be useful to evaluate the quality of the solutions returned by a heuristic procedure. Moreover, it would be interesting to investigate the application of a sampling-based approach to handle problems with a very large number of scenarios. In further work, the proposed models could also be extended by considering the uncertainty present in other parameters (e.g., costs) along with the unknown demand. Finally, another line of research could be directed toward the development of a multi-stage stochastic framework to enable the revision of the planning decisions as more information regarding the uncertain data is revealed. Even though a multi-stage model yields a better characterization of the dynamic planning process, solution techniques remain computationally challenging (Bakker et al. 2020). However, recent developments in stochastic nested decomposition techniques seem to be a promising venue for future research (Zou et al. 2019). Recall that σ = |T |/|T L | gives the total number of periods between two consecutive design periods. Since |T L | = 3, the planning horizon is divided into three sections. For |T | = 12 (resp. |T | = 24), each section spans 4 (resp. 8) time periods. The unit distribution costs are assumed to be constant in all time periods within a section. However, the costs increase between 1% and 3% from one section to the next. This type of pattern is also present in the unit variable processing costs. The latter reflect economies of scale by considering the number of modular units that are available at a particular location. For i ∈ I , we set Recall that at most four modular units may be available at a location, i.e. n i = 4 (i ∈ I ).
Regarding the generation of the fixed costs, the following procedure was employed.
In the remaining strategic periods ∈ T L \{1}, we set F O ik = U[1.01, 1.03] F O −1 ik for every k = 1, . . . , n i . Once again, the capacity size of a facility at the time it is established also affects the above costs.
• All remaining fixed costs at location i (i ∈ I ) represent a given percentage of the fixed opening costs: F E ik = 0.25 F O ik k = 1, . . . , n i − 1; ∈ T L \{1}, F R ik = 0.10 F O ik k = 1, . . . , n i − 1; ∈ T L \{1}, M t ik = 0.05 F O ik k = 1, . . . , n i ; ∈ T L ; t = , . . . , t ; t ∈ T : t < + 1. Finally, we adapt the scheme proposed in Correia and Melo (2017) to generate the unit tardiness penalty costs for late deliveries to customers j ∈ J 1 as follows: p tt j = 0.01 θ t j (t − t) 2 for t ∈ T and t = t, . . . , min{t + ρ j , |T |}, with Observe that the unit tardiness penalty cost is a function of the average fixed facility operating costs, and the average variable distribution and processing costs.