Integrated location and inventory planning in service parts logistics with customer-based service levels

Integrated Location and Inventory Planning in Service Parts Logistics with Customer-based Service Levels

Journal Pre-proof

Integrated Location and Inventory Planning in Service Parts Logistics with Customer-based Service Levels Mehmet Ferhat Candas, Erhan Kutanoglu PII: DOI: Reference:

S0377-2217(20)30096-5 https://doi.org/10.1016/j.ejor.2020.01.058 EOR 16312

To appear in:

European Journal of Operational Research

Received date: Accepted date:

27 September 2018 29 January 2020

Please cite this article as: Mehmet Ferhat Candas, Erhan Kutanoglu, Integrated Location and Inventory Planning in Service Parts Logistics with Customer-based Service Levels, European Journal of Operational Research (2020), doi: https://doi.org/10.1016/j.ejor.2020.01.058

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Elsevier B.V. All rights reserved.

Highlights • A location & inventory problem with customer-based service levels is analyzed • Network design and stock levels for service parts systems are jointly optimized • A strengthened formulation of the associated integer program is developed • A scalable methodology based on Lagrangian relaxation is developed and tested • Results show the proposed method solves large problems near optimally

1

Integrated Location and Inventory Planning in Service Parts Logistics with Customer-based Service Levels Mehmet Ferhat Candas, Erhan Kutanoglu∗ Operations Research and Industrial Engineering Graduate Program, Cockrell School of Engineering The University of Texas at Austin, Austin, TX 78712

Abstract We analyze the integrated logistics network design and inventory stocking problem in service parts logistics, where each customer requires a certain response time-based service level. We introduce a non-linear mixed integer optimization model that explicitly captures the interdependency between network design (locating facilities, and allocating customers’ demands to facilities) and inventory stocking decisions (stock levels and their corresponding stochastic fill rates). We then provide two different linearized mixed integer formulations for this problem that can solve small and medium size problems. We show that it is still a challenging problem in general, for which we have a Lagrangian-relaxation based approach that provides extremely tight lower and upper bounds. Our extensive computational study shows the effectiveness of the proposed approach for a variety of problem instances, producing provably near-optimal integrated solutions. Keywords: network design, transportation, inventory, service parts

1. Introduction An important part of post-sales service is providing necessary replacement parts to existing, geographically dispersed customers when they experience problems with their product in use. For the capital intensive, high technology or heavy machinery product manufacturers, supplying fast, high-quality post-sales service parts to customers has become increasingly important. Service is triggered to respond to a customer, potentially experiencing downtime at a critical product due to a failed part, hence the importance of time-responsiveness of the part replacement service. Therefore, as part of effective post-sales service, designing a responsive network of part stocking facilities, allocating demands from customers to these facilities and deciding part stock levels in a cost-effective and integrated manner (collectively called service parts logistics or SPL) are very important and relevant for both practice and research. ∗

Corresponding author Email addresses: [email protected] (Mehmet Ferhat Candas), [email protected] (Erhan Kutanoglu )

Preprint submitted to Elsevier

February 10, 2020

According to a survey by [22], 40 to 50% of profits of manufacturing companies come from postsales services such as replacement parts, maintenance and servicing. In certain industries, customers spend 5 to 20 times the initial sale price on subsequent service and consumables (for more information on growing importance of SPL, see [4], [13], [3]). As a result of this trend, there is a significant potential benefit in improving decision making in SPL with efficient network designs and inventory systems. Candas and Kutanoglu [4] show that integrating the SPL decision sets (locating part stocking facilities, allocating customer demands to these facilities, and selecting stock levels maintained at these facilities simultaneously) leads to significant cost savings as opposed to the traditional approach of making location and allocation decisions (collectively called logistic network design or LND) and then stocking decisions sequentially. Time-based responsiveness in SPL networks is defined through service levels. In this paper, we consider customer-based service levels: Each customer’s demand is served by a facility within a certain time window (typically measured in several hours in SPL networks) and the facility supplies the needed replacement parts with at least a certain percentage of the time given the unexpected nature of the failure of the parts. This, in a way, requires the SPL service provider to satisfy a customer’s demand within the time window with a probability that is at least equal to the customer-based service level. This is in contrast to the earlier integrated SPL studies where an overall network-wide service level was targeted [4, 13]. Note that actual SPL systems may have a mix of different service types depending on the service contracts with customers (including network-wide and customer-based). Here, we study the customer-based service level problem and develop solution methodologies specific to this problem. Customer based service levels are especially prevalent in capital intensive heavy equipment and business machinery industry where users of the equipment seeking very high uptime spends significant funds for service contracts that speculate customized attention to response time based service and the manufacturers or service providers cannot ignore the differences in service levels while designing the network and deciding the stock levels. Motivated by these real challenges in today’s SPL systems, we model the integrated network design and inventory stocking problem with customer-based service levels focusing on one critical part, where inventory decisions are explicitly considered. We provide two different formulations, each applicable to small and medium size problems. For large problems, we provide a relaxation-based algorithm that provides near optimal solutions within reasonable computation times. We note that an actual SPL system may involve many service parts and service levels that may be aggregated over multiple parts (say a product) or over multiple customers (say a region). In this paper, we focus on the single-part problem with customer-specific (disaggregate) service levels with the goal of understanding the structure of the problem and developing a powerful optimization-based heuristic. More specifically, with our decidedly simplistic problem definition focusing on customerbased service levels for one part, we are able to achieve a significant methodological development that will guarantee near-optimal solutions and structural insights. This algorithm will help in designing more scalable solution methods for more realistic problems that involve multiple parts and more aggregated service levels. 3

2. Literature Review and Contributions 2.1. Literature Review We focus our literature review on integrated network design (facility location and allocation) and inventory problems, with a particular attention on low-demand SPL systems. Given the significant issues specific to the slow moving spare/service parts systems, there has been a focused effort to model and optimize such systems explicitly. Literature on spare/service parts logistics systems includes automotive industry [8], computer and other electronic equipment service [5, 6, 7], and military [24]. Focused inventory texts on service parts logistics include [28] and [18]. A more recent and focused review [3] on service parts logistics focusing mainly on inventory issues covers research progress until 2014 and cites 125 papers. Given the vast number of studies in this domain even since 2014, we refer, as examples, to two more recent multi-facility but inventoryfocused studies in [30] and [29] and references therein. One of the early papers modifies the uncapacitated facility location problem to implicitly consider limited inventory levels [2]. Nozick and Turnquist [20] approximate inventory costs as part of the fixed facility costs and propose a model that maximizes service coverage. This is extended to minimize costs subject to service coverage constraints in [19]. The paper by Daskin et al. [10] is probably the first study that explicitly includes inventory holding cost as part of an uncapacitated facility location model. Their model assumes economic order quantity (EOQ)-based ordering and constant fill rate-based safety stocks across all facilities. The total cost function including the inventory related terms makes the overall model a nonlinear integer program which is then solved using Lagrangian relaxation. Shen et al. [27] develop a column generation-based method to solve a similar problem. A parallel and very similar model is analyzed in [16]. Yet another solution method based on conic quadratic integer programming for the same joint location and inventory problem is developed in [1]. Extensions to these earlier models analyze the tradeoffs between service and cost in joint location and inventory problems [26] and consider two echelons instead of typically assumed one-echelon networks [15]. Given the demonstrated interaction between location and inventory decisions, and the particular issues specific to service logistics systems such as time-based service levels, we see an increasing interest in modeling and solving joint network design and inventory problems for service systems. To the best of our knowledge, the first studies in this area are by Candas and Kutanoglu [4] and Jeet et al. [13], which we believe are the most relevant to our work in this paper. Candas and Kutanoglu [4] investigate the integrated SPL system and show the benefits of network design and inventory integration explicitly by solving the integrated problem using a fill rate function approximated by a step function. Jeet et al. [13] model, analyze and develop solution techniques for the multi-product version of the integrated problem using a different fill rate approximation scheme in the process. Other studies focusing on similar integration issues include [23] that studies an integrated problem focusing on repairable parts (instead of replaceable) and modeling facilities’ repair capacities explicitly. A similar capacitated and integrated SPL model with wait time constraints 4

was proposed and analyzed in [15] that proposed a Lagrangian relaxation approach to solve the underlying integer programming model. A more recent and relevant paper by Gzara et al. [11] compare relatively similar integrated SPL models with earlier models in [4, 13] and show the benefits of integration as well as the effectiveness of the focused attention on solving the integrated problem. Especially their warehouse specific service level model is relevant to our work, as our customer-specific service targets lead to facility-based service levels. Since Gzara et al. [11] solve their models directly using off-the-shelf solvers, our proposed Lagrangian relaxation-based methodology can alternatively be used to solve their models with slight modifications if needed. Finally, Sen et al. [25] reports an actual implementation of a relatively integrated model while designing Applied Materials’ SPL system. Finally, although it is not the focus of this paper, the important issue of the interaction between location and inventory decisions and aggregate multi-part service levels has not been studied extensively. We are aware of one such study, [12], that investigates the problem in the context of part commonality. Here, we define and analyze a realistic extension of an integrated problem (i.e., explicit customerbased service levels incorporated into the formulation), which lead us to develop effective algorithms based on the relaxation and decomposition of this new problem class producing probably near optimal solutions. 2.2. Contributions This paper has several contributions to the literature in terms of the problem considered, modeling, and methodology as follows: We analyze the integrated logistic network design and inventory stocking problem in service parts logistics, where service requirements are defined for each customer. Modeling and analyzing problems with specific service requirements for critical customers seen in demanding industries is very important. We introduce a non-linear mixed integer optimization model explicitly capturing the interdependency between network design and inventory stocking problems. We then provide a linearized integer formulation that uses exact fill rates (without any approximation). We also propose a strengthened integer formulation with a tighter solution space. Finally, we develop a Lagrangian-relaxation based methodology providing extremely tight lower and upper bounds to the problem. We introduce a new relaxation method called “Hybrid-relaxation”, benefiting from the strength of the Lagrangian-relaxation and the speed of the LP relaxation. We formalize the method and show the worst case performance of the Hybrid-relaxation compared to the Lagrangian-relaxation. We finally develop an efficient subgradient-based algorithm that leads to provably near-optimal solutions.

3. Problem Definition and Modeling 3.1. Notation and Modeling Assumptions We are given a set of candidate facility locations I (indexed by i), and a set of customers J (indexed by j). When we open facility i (or more precisely, locate a facility at candidate location i), we 5

incur a one-time fixed cost of fi and we incur holding cost of hi per unit time for each unit kept in stock. The unit transportation cost between facility i and demand point j is cij . Let τij be the transportation time from facility i to customer j. Comparing τij to the service time window w, we obtain δij , the identifier which takes value 1 if facility i can supply a part requested at customer j within the specified service time window (τij ≤ w), 0 otherwise (τij > w). Customer j has a service level αj as the percentage of the demand to be satisfied within the time window. The mean demand rate (for a given time unit) is dj at customer j. We make the following assumptions to facilitate model development: • We assume that network design involves the stocking facilities that are all in one echelon facing the direct demand from geographically dispersed customers. We assume that these facilities to be located are replenished from a central warehouse with infinite capacity (that is, the central warehouse can replenish the stocking facilities anytime without any delay). The lead times from the central warehouse to all facilities are known and constant (denoted by ti for facility i in the models). • Due to the low-demand nature of the motivating SPL problem, we assume that the facilities use the continuous review, one-for-one base-stock (or (S − 1, S)) replenishment policy for one critical part. This is typical as demand is low, inventory holding cost is high compared to the shipping cost, and lead time is relatively short in SPL systems. Moinzadeh and Lee [17] develop an analytical model and show that base-stock replenishment policy is optimal for these settings. Since parts are replenished one at a time, the total replenishment setup cost is constant, hence is not considered explicitly for simplicity. Main references in this area (Muckstadt [18] and Sherbrooke [28]) also rely on the wide applicability and optimality of the base stock policies in SPL systems. • We assume that demand (replacement part requests) from each customer occurs one at a time according to an independent Poisson process, which is typical in low-demand settings. We further assume that we know the mean demand rates obtained from the part failure distributions and the number of parts in use at each customer. Any unsatisfied demand due to a stockout at a facility is backordered. The part availabilities, hence the fill rates, are calculated based on the steady state analysis using the stochastic inventory theory (see [32]). See also [18] for standard assumptions in SPL inventory systems. • We assume that we have one service time window defined. Although typical SPL systems have multiple “tiered” service time windows (with increasing service level requirements for longer time windows), usually one of the time windows is most restrictive, hence assuming one window should not hinder the model’s value. Besides, modifying the model (as will be seen) for multiple windows is straightforward. In the experiments, however, we vary the time window as a control factor to see its effect on the results. • We assume that we know which customers a facility can serve within the service time window. As this is usually a function of distance and the mode of transportation available to the 6

facility and customer, we assume that this processing of transportation times is performed for each customer and facility pair a priori. We further assume that each customer’s part request is satisfied by a single direct shipment from its facility on a First Come First Served basis, without shipment consolidation or bundling. Not only this is an actual practice in SPL systems, but also it is very unlikely to have time or opportunity to consolidate multiple shipments due to low demands and strict time windows. We also assume that an open facility serving multiple customers does not reserve stock or prioritize customers, emphasizing the FCFS nature of serving customers from the facility. Given that customers have different service levels, reserving stock for high-service-level customers or rationing inventory in general may be appealing. However, calculating fill rates even in most stylized models with one facility, with already allocated Poisson demand from just two customer classes, is not possible in closed form. Given the additional issues considered in our model (design a network with multiple facilities, with demand to be still allocated from virtually many customer classes), considering inventory rationing as part of our integrated model is not viable. Our no-prioritization FCFS assumption can be viewed as a special case of an inventory rationing policy where the reservation threshold is set to 0. That is, the inventory levels and the total costs from our model are upper bounds for the ultimate problem with inventory rationing. A post processing step can potentially adjust the individual stock levels (typically lower) to take advantage of inventory rationing. We highlight this opportunity as a future research direction later. • We finally assume that it is not possible to split a customer’s demand while allocating it to facilities, i.e., each customer has to be served by a single stocking facility (single sourcing). 3.2. On the Fill Rate Function For a facility facing Poisson-distributed demand, the fill rate function is defined as β(s, λ) = G(s − 1, λ) for (base) stock level s and mean lead time demand λ, where G is the cumulative (Poisson) Ps−1 k −λ distribution function, G(s − 1, λ) = k=0 λ e /k! [32]. The models formulated later make use of the following results (See Figure 1 (later in the paper) for an illustration of fill rate as a decreasing function of the mean lead time demand λ for stock levels 1 and 2): Proposition 1. For a given stock level s, β(s, λ) is strictly decreasing function of λ. Proposition 2. The fill rate is a concave function of λ for given stock level s, when λ < s − 1, and it is convex when λ ≥ s − 1. As assigning more demand to a facility with a given stock level causes its fill rate to decrease, there is a limit on the amount of demand the facility can serve if a target fill rate (a lower bound) is to be achieved for the assigned demand. Suppose now that we need to guarantee a fill rate level of b at a facility with stock level s. We define parameter µ(s) as the maximum mean lead time demand that can be assigned to this facility to guarantee a fill rate level of b or higher. That is, µ(s) is solved from the equation 7

b=

s−1 X µ(s)k e−µ(s)

k!

k=0

.

(1)

With the results above, we note that µ(s) is an increasing function of s for a given b. Although µ(0) and µ(1) can be solved explicitly from equation (1) as a closed form function of other parameters, µ(s) for s ≥ 2 does not have a closed-form equation. However, they can be calculated numerically using mathematical software. More specifically, the fill rate as a function of stock level s and mean lead time demand λ is β(s, λ) = Q(s, λ), where Q(s, λ) is the incomplete regularized Ps−1 r −λ is the incomplete gamma function, gamma function: Q(s, λ) = Γ(s,λ) r=0 λ e Γ(s) . Here, Γ(s, λ) = and Γ(s) is the complete gamma function with Γ(s) = (s − 1)!. (We calculate µ(s) parameters for given fill rates and stock levels numerically, using the inverse of the regularized gamma function in Mathematica). 3.3. Mathematical Formulation In this section we introduce two formulations for the integrated location and inventory problem with customer-based service levels. To facilitate both formulations, we assume that there is a finite number of alternatives for stock levels Si . We denote its largest element by Li , which can be calculated by considering maximum possible demand that can be assigned to the facility. Hence, we consider Si ∈ {0, 1, 2, . . . , Li } for all i (for brevity, we use Li for both largest element of this set, and also for the set itself). We make this set large enough to guarantee practically 100¯ α% fill rate, where α ¯ = maxj∈J {αj }, at the largest stock level Li for the maximum total demand that can be assigned to a facility. Hence, constraints (11) allow the stock levels to be selected from the initially developed set of integer stock levels. The formulation below (denoted as the Customer Based Problem, or CBP ) uses stock levels (Si ), fill rates (βi ) and mean lead time demands (λi ) as decision variables which are simultaneously optimized along with the facility locations (Yi ) and customer assignments (Xij ): X XX X min fi Yi + cij dj Xij + hi Si (2) i∈I

i∈I j∈J

X

i∈I

Xij = 1, ∀j ∈ J

(3)

Xij ≤ Yi ∀i ∈ I, j ∈ J

(4)

Si ≤ Li Yi , ∀i ∈ I X λi = ti dj Xij , ∀i ∈ I

(5)

subject to

i∈I

(6)

j∈J

SX i −1

e−λi , ∀i ∈ I r!

(7)

βi ≥ αj Xij , ∀i ∈ I, j ∈ J

(8)

Xij = 0 or 1, ∀i ∈ I, j ∈ J

(9)

Yi = 0 or 1, ∀i ∈ I 8

(10)

βi =

r=0

λri

Si ∈ {0, 1, 2, . . . , Li } , ∀i ∈ I.

(11)

The first term in the objective is the total fixed facility costs, calculated with fi associated with facility opening decisions Yi , the second term is the transportation/shipping costs for assignments Xij , and the third term is the inventory costs, which applies unit holding cost hi to stock level Si . With this in mind, constraints (3) guarantee all demand at each customer to be fully assigned to a facility. Constraints (4) mean that any facility serving a customer should be open. Constraints (5) allow stock levels to be greater than 0 only for open facilities. Constraints (8) are the time-based service level constraints, each for a customer, stating that each customer has to be served by a facility having the fill rate higher than the customer’s target service level. Here βi (calculated in constraints (7)) is the fill rate at facility i at its stock level Si and at its mean lead time demand λi (calculated in constraints (6)), where ti is the lead time for facility i. Note that, in the customer based service problem (CBP), customer j has to be served with the promised “positive” service level αj > 0. Therefore, we do not consider assignments of the customers to the facilities outside the time window. For simplicity, we define set of facilities Ij for each customer j, as facilities that can serve customer j within the time window, Ij = {i ∈ I : δij = 1}, and similarly we define set of customers Ji for each facility i, as customers that can be served by facility i within the time window, Ji = {j ∈ J : δij = 1}. Here assignment variables Xij are defined only for j ∈ J and i ∈ Ij . With this simplification, we can rewrite the service level constraints (8) as βi ≥ αj Xij for all j ∈ Ji . Briefly, constraints (8) state that we should choose the minimum inventory level Si∗ producing a fill rate higher than the maximum target service level of the served customers Jˆi , (minSi ∈Li βi ≥ maxj∈Jˆi {αj }). The model (2)-(11) is a nonlinear mixed integer programming problem. The source of nonlinearity is the service level constraints (8), since the fill rate variables βi are nonlinear functions of two other decision variables, stock levels Si and mean lead time demands λi . Note that constraints (8) have two functions: (1) deciding the required fill rate at facility i, and (2) deciding the inventory level at the facility to satisfy this fill rate (along with constraints (7). To simplify these constraints further, we handle decisions of fill rate and inventory levels in two sets of constraints. First, we predefine a set of potential fill rates that a facility should achieve, depending on the customers potentially served by that facility. To find the minimum potential required fill rates, called “base fill rates”, at the facilities, we define set Bi of potential fill rates for each facility i. If there are Ni (for brevity we use Ni for both largest element of this set, and also the set itself) distinct service level requirements in set Ji (usually Ni is much smaller than |Ji |), then Bi = {bi1 , bi2 , . . . , bi,Ni } is the set of distinct service level requirements of the customers in Ji (without loss of generality we assume that bin ’s are in increasing order in Bi ). We now define binary variable Win , which takes value 1 if the base fill rate at the facility i is equal to bin for every i and n. To guarantee that the facilities have fill rates greater than or equal to the service level of their served customers, we define the constraints: 9

X

n∈Ni : bin ≥αj

Win ≥ Xij , ∀j ∈ J, i ∈ Ij

P (Note that, this constraint can be alternatively written as n∈Ni bin Win ≥ αj Xij , Below, we use Ni (αj ) = {n ∈ Ni : bin ≥ αj } to simplify the model.

∀j ∈ J, i ∈ Ij ).

To facilitate the second function of the service level constraints (deciding the stock levels), we calculate the parameters µins for all i ∈ I, n ∈ Ni and s ∈ Li from the fill rate function such
P k −µins /k!. Hence, µ that bin = s−1 ins denotes the maximum mean lead time demand k=0 (µins ) (e) that can be assigned to facility i with stock level s to achieve a fill rate of at least bin . For an illustration of the calculation of µins , see Figure 1, which shows fill rate βi for changing mean lead time demand λi and for stock levels Li = {1, 2}, assuming two distinct service levels (Ni = 2) where Bi = {bi1 = 0.5, bi2 = 0.7}. As an example from the figure, to guarantee at least 70% fill rate (bin = 0.7 with n = 2), with 1 unit of stock, facility i can accept at most µi,n=2,s=1 = 0.36 P units of mean lead time demand (ti j∈Ji dj Xij ), and with 2 units of stock, facility i can accept up to µi,s=2,n=2 = 1.1 units of mean lead time demand. We now define binary variable Qis , ∀i, s, which takes value 1 if the stock level at facility i is equal to s. Then, the following constraints guarantee the required stock level s to be selected: X X µins Qis + Mi (1 − Win ), ∀i ∈ I, n ∈ 1, ..., Ni dj Xij ≤ ti s∈Li

j∈Ji

P where Mi = n (hence Wiˆ n = 1), the j∈Ji dj . Note that, if the base fill rate at facility i is biˆ P P constraint will remain as ti j∈Ji dj Xij ≤ s∈Li µiˆns Qis , which as a forcing constraint ensures that the assigned demand for facility i is such that the mean lead time demand for the corresponding stock level and fill rate is not exceeded. We now write the complete, linearized and simplified model:   X XX X X min fi Yi + cij dj Xij + hi  sQis  (12) i∈I

i∈I j∈Ji

i∈I

X

Xij = 1, ∀j ∈ J

(13)

Xij ≤ Yi , ∀i ∈ Ij , j ∈ J X Win , ∀j ∈ J, i ∈ Ij

(14)

µins Qis + Mi (1 − Win ), ∀i ∈ I, n ∈ Ni

(16)

subject to

i∈Ij

Xij ≤ ti

X

j∈Ji

dj Xij ≤

X

s∈Lin

s∈Li

(15)

n∈Ni (αj )

X

Qis ≤ 1, ∀i ∈ I

(17)

Win ≤ 1, ∀i ∈ I

(18)

Xij = 0 or 1, ∀i ∈ I, j ∈ Ji

(19)

s∈Li

X

n∈Ni

10

Win = 0 or 1, ∀i ∈ I, n ∈ Ni

(20)

Qis = 0 or 1, ∀i ∈ I, s ∈ Li .

(21)

βi 0.9 0.8

bi,n=2=

0.7 0.6

bi,n=1=

0.5

s=2

0.4 0.3 0.2 0.1

s=1 0

0.5

µi,n=2,s=1

1

µi,n=1,s=1

µi,n=2,s=2

1.5

2

µi,n=1,s=2

λi

Figure 1: Calculation of µins assuming two distinct service levels (Ni = 2) in customer set Ji , where Bi = {bi1 = 0.5, bi2 = 0.7} for stock levels Li = {1, 2}

Here, the first term in the objective (12) is the fixed cost associated with facility opening decisions Yi , the second term is the transportation cost for assignments Xij , and the third term is the inventory stocking costs captured with stock choice (binary) variables Qis . Constraints (13) assign all the demands to a facility within the time window while constraints (15)-(16) guarantee the service level requirements. Note that in the formulation above, we are using different maximum potential inventory levels Lin for each facility i and the nth base fill rate level bin . Since the problem size increases with Li , we use the tightest value, depending on the facility (since each facility may cover different amount of total demand within its time window) and on the fill rate level (since the required inventory level increases with fill rate requirements). Therefore, we calculate Lin by solving for it P in −1 k −λi
P in L λi e /k! ≥ bin , where λi = ti j∈Ji dj . Using Figure 1 as an illustration, if the total k=0 mean lead time demand within facility i’s time window coverage is 0.5 units (λi = 0.5), to guarantee a bin = 50% fill rate, facility i may need 1 unit of stock. However, to guarantee a bin = 70% fill rate, facility i may need up to 2 units of stock. Hence, for the example in Figure 1, Li1 = 1 and Li2 = 2. We use the notation Li for the maximum of Lin ’s for facility i (for this example Li = 2). With this formulation (which we call formulation F 1 as we introduce a strengthened version called F 2 later), we can solve small and medium size problems (i.e., up to 30 facilities and 100 customers), but scalability to larger problems is an issue. From our experiments, we see that F 1 has weak LP relaxation bounds, resulting in long solution times to close the optimality gap. To improve this formulation, we seek to tighten the formulation, by adding cuts and/or having stronger constraints with higher dimensions if possible. In the following section, we propose an alternative formulation, 11

denoted by F 2 with a tighter feasible space, and then compare it with F 1 in terms of solution times and initial LP relaxation bounds. 3.4. Improved Formulation Instead of using separate binary variables Yi , Win and Qis for choosing facilities, base fill rate levels and stock levels, we now define composite variables Zins which is 1 when facility i is open with fill rate bin and stock level s, and 0 otherwise, for each i, n, and s. We also expand assignment variables Xij to Xinsj . The idea behind this notation is to view each (i, n, s) combination as a “network node” with predefined fill rate and stock level that a customer may be assigned. We also introduce fis = fi + hi s as the fixed cost of opening facility i with stock level s, and µ ¯ins = µins /ti as the maximum annual mean demand that can be served by facility i with fill rate of bin and stock level of s. The new formulation is as follows: min

XX X

XX XX

cij dj Xinsj

(22)

Xinsj = 1 ∀j ∈ J

(23)

Xinsj ≤ Zins , ∀j ∈ J, i ∈ Ij , n ∈ Ni (αj ), s ∈ Li X dj Xinsj ≤ µ ¯ins ∀i ∈ I, n ∈ Ni (αj ), s ∈ Li

(24)

Zins ≤ 1 ∀i ∈ I

(26)

Xinsj = 0 or 1, ∀i ∈ I, n ∈ Ni (αj ), s ∈ Li , j ∈ Ji

(27)

Zins = 0 or 1, ∀i ∈ I, n ∈ Ni , s ∈ Li .

(28)

i∈I n∈Ni s∈Li

subject to

fis Zins + X

i∈I n∈Ni s∈Li j∈Ji

X

X

i∈Ij n∈Ni (αj ) s∈Li

j∈Ji

X X

n∈Ni s∈Li

(25)

Similar to the objective in previous model (12), the first term in the objective (22) combines the first and last terms of the original objective, representing the fixed cost associated with both the P P P facility opening decisions (Yi = n∈Ni s∈Li Zins ), and stock level decisions (Qis = n∈Ni Zins ). The second term in the objective is the transportation cost for assignments Xinsj to network nodes (i, n, s). The constraints are analogous to the previous model (13) - (18). Constraints (23) assign all demand to a network node within the time window guaranteeing the service level requirements. Constraints (24) mean that any node serving a customer should be open. Constraints (25) guarantee that total assignment to a network node does not exceed the capacity of the node. Constraints (26) are the limiting constraints, preventing any facility to be open for more than one fill rate or stock level. 3.5. Comparing Formulations F1 and F2 We now discuss advantages and disadvantages of formulations F 1 and F 2 and compare them in terms of their problem sizes (numbers of variables and constraints), solution times, and quality of 12

the LP relaxations with a set of computational experiments. Then, we select one of the formulations as the basis for the methodology development later in the paper. The formulation with composite and expanded variables are known to have stronger relaxation bounds (see [9] for an example of such use of composite variables). Although we do not provide a formal proof for F 2 having stronger relaxation than F 1, we can find feasible solutions for the LP relaxation of F 1 which are not feasible for the LP relaxation of F 2. As an example, if customer j 0 is fully assigned to facility i0 with fill rate bi0 n0 and stock level s0 , assuming µi0 n0 s0 = 2ti dj 0 , solution (Yi0 , Xi0 j 0 , Wi0 s0 , Qi0 s0 ) = (1, 1, 1, 0.5) (keeping half of stock level s0 ) is feasible in the LP relaxation of F 1. However, the same solution (opening facility i0 with fill rate bi0 n0 and half of stock level s0 ) is not feasible in the LP relaxation of F 2, due to forcing constraints Xinsj ≤ Zins . Formulation (12)-(21) for CBP , F 1, has |I|(1 + |J| + |N | + |S|) variables and |I|(2 + 2|J| + |N |) + |J| constraints, whereas the composite variable formulation, F 2, has |I||J||N ||S| binary and |I||N ||S| continuous variables, with |I| + |J| + |I||N ||S|(1 + |J|) constraints. For the medium size instance we use in our experiments (with 50 facilities, 200 customers, low mean demand, and 2 hours time window), F 1 has about 1,500 variables and 3,000 constraints, whereas F 2 has about 11,000 variables and 11,000 constraints. The main reason for F 2’s size to be much larger than F 1 is the assignment variables for all (i, n, s, j) combinations. Table 1 summarizes the experimental data set we use to compare F 1 and F 2. We generate problem instances using our industrial partner’s real data set as basis, especially for network sizes, time windows, demand rates, service levels, and cost tradeoffs. For the sake of generating multiple spatial spreads of customers and facilities, we randomize “maps” of candidate facilities and customers generating them in a fixed size stylized grid. To this end, we randomly assign facility and customer locations to a 100 by 100 grid and generate two different data sets, denoted by a and b, with different size instances: small |I| = 25, |J| = 100, medium |I| = 50, |J| = 200, and large |I| = 100, |J| = 400. Using a year as the common time unit, mean demand rates are generated uniformly, with a mean of 2.5 (continuous uniform between 1 and 4) for low demand instances, and a mean of 6 (continuous uniform between 4 and 8) for high demand instances. (These rates are a bit larger than those in the actual data, but to see the effect of inventory on network decisions, the magnification of these rates is necessary, given that we are analyzing single part problems. When a large number of parts is explicitly modeled, the actual demand rates can be used as the inventory costs aggregated across multiple parts will naturally lead to the tradeoff between network design and inventory). Inspired from the real data, we have two time windows, 2 hours and 4 hours, which are assumed to be 20 units and 40 units of Euclidean distance in the 100x100 grid. Each instance is run with low ($250), medium ($500) and high ($1000) holding costs. Fixed facility costs are $1000 and transportation costs are distance-based (5$ per unit distance) and average transportation cost is around $150. We randomly split customers into 3 groups, requiring high (90%), medium (70%) and low (50%) target service levels, in a way simulating the so called gold, silver and bronze customer segments. We limit the maximum solution time to 3600 seconds. Tables 2 and 3 summarize the results of the experiments for data sets a and b respectively. We 13

Table 1: Experimental data set

Data Size (|I|, |J|) Mean demand, d Time window, w Holding cost, h Fixed cost, f Transportation cost, c

a, b small (25, 100), medium (50, 200), large (100, 400) low [1, 4], high [4, 8] 2 hours, 4 hours $250, $500, $1000 $1000 $5 per unit distance

Total

72 instances (named as a1-a36 and b1-b36)

provide initial LP relaxation solution values, best feasible solution values (which is optimal if completed within 3600 seconds), the solution times and the status of the final solutions obtained (where N/F means “not finished,” and Opt means optimal is found). We also provide final gaps between the best lower and upper bounds obtained from the optimization solver. Note that this is zero when an optimal solution is found. We solve these problems with Xpress-MP by FICO (www.fico.com). Note that, Xpress-MP uses presolve and cut generation techniques together with branch-and-bound, which have significant effect on the solution time (e.g., with formulation 2, the instance a13 is solved to optimality within 36 seconds, whereas it has more than 2% gap at the end of 3600 seconds when cut generation is disabled). As shown in Table 2, LP relaxation of F 1 is solved very quickly (less than 2 seconds in all of the instances), however they provide weak bounds. Meanwhile, the average LP solution time of F 2 is much longer, and F 2 can not even find the LP relaxation solution within 3600 seconds for some large instances (a35, a36 and b36). On the other hand, the average improvement of F 2 LP bounds (when it has one) over F 1 bounds is about 120%, ranging from 46% (for instance a7) to 210% (for instance a30). We see that F 1 cannot be solved to optimality in any instances. Actually, F 1 solutions for four instances are optimal (a1, a2, and b1, b2), but they are not proven optimal. On the other hand, F 2 finishes almost all of the small and some of the medium size instances to optimality within 3600 seconds (average improvement over F 1 upper bound is 17% for these instances). However, for the large instances F 2 looses its superiority over F 1, and for 16 of the large instances (a28-a30, a33-a36, and b26-b30, b33-b36) F 1 upper bounds are better. We conclude that formulation F 2 has faster convergence and better quality lower bounds in general, but it struggles as problem size increases. Note that, for the largest problem (a36), F 1 has about 5000 variables and constraints, but F 2 has more than 70000 variables and constraints. In addition, F 2 also has structural advantages for decomposition. Therefore, we use F 2 for a basis of our methodology development to solve larger instances.

14

15

ins a1 a2 a3 a4 a5 a6 a7 a8 a9 a10 a11 a12 a13 a14 a15 a16 a17 a18 a19 a20 a21 a22 a23 a24 a25 a26 a27 a28 a29 a30 a31 a32 a33 a34 a35 a36

large

medium

small

size

high

low

high

low

high

low

d

4

2

4

2

4

2

4

2

4

2

4

2

w (hrs)

h ($) 250 500 1000 250 500 1000 250 500 1000 250 500 1000 250 500 1000 250 500 1000 250 500 1000 250 500 1000 250 500 1000 250 500 1000 250 500 1000 250 500 1000

LP soln 15790 16766 18507 12446 12446 12446 34392 35774 38124 30214 30214 30214 28371 30213 33454 22004 22004 22004 59058 62157 66603 51892 51892 51892 31491 31671 31769 29860 29860 28985 73157 73664 74439 69902 69902 69902

Time 0.1 0.0 0.0 0.1 0.1 0.1 0.0 0.0 0.0 0.1 0.1 0.1 0.1 0.1 0.1 0.3 0.3 0.3 0.1 0.1 0.1 0.3 0.3 0.3 0.5 0.5 0.6 1.4 1.4 1.4 0.6 0.6 0.6 1.4 1.4 1.4

Formulation 1 Best UB Time 30210 3600 35460 3600 46107 3600 29992 3600 33801 3600 43242 3600 52962 3600 61611 3600 78867 3600 55735 3600 63343 3600 81619 3600 56782 3600 70984 3600 90660 3600 56655 3600 68430 3600 84328 3600 102804 3600 119603 3600 151092 3600 109947 3600 126925 3600 158591 3600 104797 3600 122430 3600 167634 3600 120452 3600 149273 3600 203801 3600 183793 3600 223200 3600 283862 3600 195209 3600 234224 3600 346458 3600 Status N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F

Gap 2.6 3.1 6.0 25.2 30.8 43.9 4.8 6.8 10.8 16.7 19.0 32.4 21.2 30.7 35.0 40.2 50.4 59.0 19.3 23.0 31.0 32.8 39.4 51.2 54.2 59.5 69.6 63.0 69.9 79.6 42.9 52.2 61.2 46.5 55.6 69.8

LP soln 28848 32793 40233 27217 29551 33455 50318 56246 67781 50205 55592 64499 50650 56515 67558 47707 51535 57504 91249 100756 117951 90008 98459 111211 75499 82276 93788 75127 81173 89241 138954 153504 175875 138954

Table 2: Comparison of F 1 and F 2 with data a Time 0.4 0.3 0.3 2.7 3.6 5.9 0.4 0.4 0.4 2.8 5.0 6.9 2.7 3.0 3.8 45.2 71.2 119.5 2.5 3.2 3.8 70.5 98.8 185.9 105.4 164.5 242.2 1587.8 2392.6 3425.2 82.7 185.0 331.3 1801.4

Formulation 2 Best UB Time 30210 1 35460 2 45960 4 28919 18 32889 53 39906 57 52775 1 61275 6 78016 163 52775 11 61079 604 75628 2769 53499 36 62749 1004 80723 3644 50602 1252 57939 3620 74209 3615 96196 45 111332 3600 142984 3600 95116 968 111791 3600 137158 3600 82899 3600 118778 3600 154298 3600 199372 3600 245070 3600 3600 150592 3600 187268 3600 332097 3600 3600 3600 3600 Status Opt Opt Opt Opt Opt Opt Opt Opt Opt Opt Opt Opt Opt Opt N/F Opt N/F N/F Opt N/F N/F Opt N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F

0.6 6.0 34.4

Gap 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.2 0.0 0.3 4.4 0.0 0.5 3.2 0.0 2.4 1.5 1.5 20.2 22.4 60.3 64.9

16

ins b1 b2 b3 b4 b5 b6 b7 b8 b9 b10 b11 b12 b13 b14 b15 b16 b17 b18 b19 b20 b21 b22 b23 b24 b25 b26 b27 b28 b29 b30 b31 b32 b33 b34 b35 b36

large

medium

small

size

high

low

high

low

high

low

d

4

2

4

2

4

2

4

2

4

2

4

2

w (hrs)

h ($) 250 500 1000 250 500 1000 250 500 1000 250 500 1000 250 500 1000 250 500 1000 250 500 1000 250 500 1000 250 500 1000 250 500 1000 250 500 1000 250 500 1000

LP soln 13931 14294 14747 11874 11874 11874 31997 32710 33689 28832 28832 28832 21953 21992 22063 20989 20989 20989 49385 49448 49506 47299 47299 47491 32311 32311 32311 31599 31599 31599 75791 75791 75791 74116 74116 74116

Time 0.0 0.0 0.0 0.1 0.1 0.0 0.0 0.0 0.0 0.1 0.0 0.1 0.1 0.1 0.1 0.3 0.3 0.3 0.1 0.1 0.1 0.3 0.3 0.3 0.6 0.5 0.5 1.4 1.4 1.4 0.5 0.6 0.5 1.3 1.3 1.3

Formulation 1 Best UB Time 31983 3600 37521 3600 49400 3600 29582 3600 35402 3600 42476 3600 54905 3600 64159 3600 82019 3600 55646 3600 64794 3600 82011 3600 55881 3600 69385 3600 94409 3600 57533 3600 69132 3600 86882 3600 100573 3600 120554 3600 152238 3600 104764 3600 128708 3600 155684 3600 103519 3600 130165 3600 166456 3600 118857 3600 154471 3600 206096 3600 184609 3600 221796 3600 282311 3600 208510 3600 251581 3600 354921 3600 Status N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F

Gap 5.3 6.6 12.7 20.6 29.8 33.1 7.3 9.7 14.2 14.6 21.7 32.7 23.9 33.6 43.7 42.6 51.6 60.9 22.3 28.9 35.4 30.4 41.3 53.0 54.4 63.1 69.5 62.1 70.9 78.1 43.5 52.2 60.6 49.9 58.9 70.7

LP soln 30073 33853 41268 27405 29846 33775 52149 58337 69938 51104 56419 64260 48232 53733 63913 46220 50496 57102 85871 95181 111608 85769 94272 105628 75281 81615 92554 75095 81268 91327 139622 152793 173579 139622 152792

Table 3: Comparison of F 1 and F 2 with data b Time 0.4 0.3 0.3 1.1 1.6 2.2 0.5 0.4 0.4 1.9 3.4 3.6 2.5 3.1 3.4 45.6 57.5 118.2 2.4 3.4 5.3 51.0 97.5 209.2 90.2 125.8 190.7 1248.8 1885.7 2935.0 92.1 174.4 271.4 1968.2 3246.6

Formulation 2 Best UB Time 31983 3 37521 2 48521 4 28881 8 32656 9 39686 25 54629 2 63415 15 80467 238 53411 27 61369 281 76269 3600 51519 1894 59807 1367 76830 3600 48994 360 56520 3600 73142 3600 90890 190 105503 2394 136857 3600 90884 2409 105705 3600 134057 3600 81872 3600 134395 3600 179903 3600 151397 3600 244420 3600 315842 3600 154180 3600 181316 3600 305612 3600 3600 3600 3600 Status Opt Opt Opt Opt Opt Opt Opt Opt Opt Opt Opt N/F Opt Opt N/F Opt N/F N/F Opt Opt N/F Opt N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F N/F

2.4 3.0 28.4

Gap 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.4 0.0 0.0 1.4 0.0 0.4 4.6 0.0 0.0 2.7 0.0 0.6 3.1 0.5 29.7 33.8 47.0 64.3

4. Methodology We propose two lower bounding techniques based on Lagrangian relaxation in Section 4.1. Then, we develop an upper bound algorithm in Section 4.2. Finally, we combine these techniques in a subgradient algorithm to obtain solutions close to optimum for large problems. 4.1. Lower Bound In this section, we first provide a lower bounding algorithm based on Lagrangian relaxation. Next, we provide a “hybrid Lagrangian-LP relaxation” technique which uses the advantages of Lagrangian relaxation and LP relaxation together. 4.1.1. Lagrangian Relaxation In formulation (22 - 28), only the assignment constraints (23) tie the facilities together. The other constraints can be decomposed across facilities, defining a subproblem for each facility. We relax constraints (23) with Lagrangian multipliers denoted by πj . The objective function of the relaxed problem is the following:   XX X XX XX X X X X fis Zins + cij dj Xinsj + πj 1 − Xinsj  . i∈I n∈Ni s∈Li

i∈I n∈Ni s∈Li j∈Ji

i∈Ij n∈Ni s∈Li

j∈J

Combining Lagrangian penalties and transportation costs as c¯ij = cij dj −πj , we obtain the following as the Lagrangian relaxation problem, LR(π): X XX XX XX X c¯ij Xinsj + πj (29) fis Zins + min j∈J

i∈I n∈Ni s∈Li j∈Ji

i∈I n∈Ni s∈Li

subject to Xinsj ≤ Zins , ∀j ∈ J, i ∈ Ij , n ∈ Ni (αj ), s ∈ Li X dj Xinsj ≤ µ ¯ins ∀i ∈ I, n ∈ Ni (αj ), s ∈ Li

(30)

Zins ≤ 1 ∀i ∈ I

(32)

Xinsj = 0 or 1, ∀i ∈ I, n ∈ Ni (αj ), s ∈ Li , j ∈ Ji

(33)

Zins = 0 or 1, ∀i ∈ I, n ∈ Ni , s ∈ Li .

(34)

j∈J

X X

n∈Ni s∈Li

(31)

A solution method for LR is described with the next proposition. First, note that the relaxed problem is decomposed into facility subproblems. Subproblem Pi for each facility i is stated as: X X X XX Pi : min fis Zins + c¯ij Xinsj (35) n∈Ni s∈Li

n∈Ni s∈Li j∈Ji

subject to Xinsj ≤ Zins , ∀j ∈ Ji , n ∈ Ni (αj ), s ∈ Li X dj Xinsj ≤ µ ¯ins ∀n ∈ Ni (αj ), s ∈ Li j∈Ji

17

(36)

(37)

X X

Zins ≤ 1

(38)

Xinsj = 0 or 1, ∀n ∈ Ni (αj ), s ∈ Li , j ∈ Ji

(39)

Zins = 0 or 1, ∀n ∈ Ni , s ∈ Li .

(40)

n∈Ni s∈Li

With the optimal solution values v(Pi ) to problems Pi , we can compute the optimal value of the P P Lagrangian relaxation problem as v(LR(π)) = i v(Pi ) + j∈J πj . The following proposition shows that each facility problem Pi can be further decomposed into base fill rate and stock level problems, Pins : Proposition 3. For each i, v(Pi ) = minn∈Ni ,s∈Li {0, fis + v(Pins )} where Pins :

min

X

c¯ij Xinsj

(41)

dj Xinsj ≤ µ ¯ins

(42)

Xinsj = 0 or 1, ∀j ∈ Jin

(43)

X

subject to

j∈Jin

j∈Jin

where Jin = {j ∈ Ji : αj ≤ bin }. Proof. Proof of Proposition 3. For each fixed i, Pi can be further decomposed into sub-problems Pins for each n and s, with the additional requirement that (via constraints (38)) at most one of those sub-problems can be “active.” Hence we have: v(Pi ) = min

XX n

Zins (fis + v(Pins ))

s

XX

subject to

n

s

Zins ≤ 1

(44) (45)

Since setting Zins to 0 is feasible (corresponding to not opening facility (i, n, s)), problem Pi can be solved by solving problems Pins , for all n ∈ Ni , and s ∈ Li , and setting v(Pi ) = minn∈Ni ,s∈Li {0, fis + v(Pins )}. Corollary 4.1. LR(π) is solved by finding solutions to a series of knapsack problems for each facility i, base fill rate level n, and inventory level s as Pins , (41)-(43). Using Corollary 4.1, we develop Algorithm 1 to solve Pi . The complexity of Algorithm 1 is determined by the number of knapsack problems to be solved and also by the complexity of the knapsack problem itself, which is known to be NP-hard. However, there are practically efficient algorithms such as the dynamic programming algorithm in [14] and [21] for the knapsack problems. One well known algorithm for the knapsack problems is constructed on the idea of dynamically increasing the knapsack capacity, which has a complexity of O(|J|Di logDi ), P where Di = j∈Ji dj is the capacity of the knapsacks (i.e., the total demand in our case). Note P that, overall, we solve n∈Ni Lin knapsack problems; thus, keeping Lin as tight as possible has 18

significant effect on the solution time. Then, the complexity of solving problem Pi with Algorithm P 1 is O(( n∈Ni Lin )|J|Di logDi ).

We need to solve many knapsack problems to compute the LR lower bound, e.g., for a problem instance of 100 facilities and 400 customers, on average we solve 15 knapsack problems per facility. Therefore, to control the computation time, we focus on decreasing the number and size of knapsack problems to be solved. We employ the following ideas: i. Note that customers with nonnegative c¯ij can not be in the solution. Hence, we construct set − Ji− = {j ∈ J : c¯ij < 0} (and similarly Jin = {j ∈ Ji− : αj ≤ bin }) in Algorithm 1 for facility i. This results in smaller knapsack problems (each with smaller number of customers). ii. Recall that we calculate Lin by considering the total lead time demand λi of the customers in set Ji (see discussion in Section 3.3). Now, with the simplification above, we recalculate P − the total lead time demand λ− i = ti j∈Ji− dj where λi ≤ λi , which leads to a new set of − stock levels L− in with |Lin | ≤ |Lin |. Hence, the number of knapsack problems to be solved is P P reduced from n∈Ni Lin to n∈Ni L− in .

iii. For the largest stock level for base fill rate n, L− ) = in , we can simply calculate v(Pi,n,L− in P − P ¯ij , since, for s = Lin , j∈J − dj ≤ µ ¯ins (hence constraints (42) is satisfied). Then, we j∈Ji− c i solve knapsack problems Pins up to stock level L− in − 1. iv. Similarly, we construct the set of base fill rates Ni− by considering the customers in Ji− (there may not be any customer in set Ji− requiring a specific fill rate bin ). In addition to the ideas above, the key idea that will further reduce the number of knapsack problems is “fathoming” the knapsack problems by using lower and upper bounds (on the knapsack problems) (The term fathoming refers to skipping solving a knapsack problem altogether). Since we seek the best (n, s) pair for each facility i, by using a good incumbent solution for i (i.e., upper bound of a good (ˆ n, sˆ) pair with upper bound U B(Piˆnsˆ), U Bi∗ = U B(Piˆnsˆ)), we can fathom the knapsack problems for all (n, s) with LB(Pins ) ≥ U Bi∗ , where LB() and U B() are lower and upper bound values for the corresponding knapsack problems. Note that we can find very tight lower and upper bounds for knapsack problems efficiently. A well known lower bound technique for the knapsack problem is filling up the capacity starting with the customers having lowest c¯ij /dj ratio, and ignoring the integrality requirement for the last customer in the knapsack (so we can include some portion of the last customer to fully utilize the capacity). This algorithm gives the same bound as the LP relaxation. Note that if we remove the last customer from the knapsack, we obtain a feasible solution, hence an upper bound. We can further improve this solution by checking the remaining customers to fill the remaining capacity. We apply all these ideas to Algorithm 1 in an efficient manner as shown in the new version, called Algorithm 2. Note that, the time spent for actually solving the knapsack problems to optimality dominates the overall computational burden. With the ideas used in Algorithm 2, the number of knapsacks solved 19

is reduced drastically (see results in Section 4.3 for details). Still, we need to solve many knapsack problems, particularly for solving large instances. In the next section, we provide a hybrid scheme, as an alternative to conventional Lagrangian relaxation, to further control the computational cost due to knapsack problems. In Section 4.3 we compare these alternative lower bounding techniques to finalize our methodological development. 4.1.2. Hybrid Relaxation. The purpose in this section is to take advantage of the easily computable LP relaxation bounds of the knapsack subproblems and avoid solving an integer knapsack problem whenever its LP relaxation bound “signals” that the improvement (in the facility-level subproblem’s value) due to this knapsack’s integer solution is very small. Ultimately, this results in a small loss in the tightness of the lower bounds obtained through Lagrangian relaxation, but the loss (and the gain in computational time due to unsolved integer knapsack problems) can be controlled through a tolerance parameter which is used to decide whether it is sufficient to use the LP relaxation bound or if solving the integer knapsack problem is needed. As the decisions about which knapsack problems should be solved to optimality depend on the values of the Lagrangian multipliers and “dynamically” determined in each subgradient iteration, we can view the underlying model as a further relaxation of the Lagrangian relaxed problem, LR(π). In this new relaxation, some of the binary variables are treated as continuous (as in the conventional LP relaxation), hence one can view this as a hybrid between LP relaxation and Lagrangian relaxation, in which the subset of variables to be relaxed is decided in each iteration according to the value of the tolerance parameter and the quality of the solutions. The overall idea is to benefit from both the simplification arising from the LP relaxation (and computational savings) and the tightness of the Lagrangian relaxation. We can formulate the so called hybrid relaxation (HR) by replacing the constraints (33) in the LR(π) (29)-(34), with constraints (46) and (47), allowing some Xinsj variables to be non-integer. With this modification, we conclude that v(HR(π)) ≤ v(LR(π)) ≤ v(CBP ). 0 ≤ Xinsj ≤ 1, ∀(i, n, s, j) ∈ (I, Ni (αj ), Li , Ji )LP Xinsj = 0 or 1, ∀(i, n, s, j) ∈ (I, Ni (αj ), Li , Ji ) \ (I, Ni , Li , Ji )

LP

(46) (47)

Note that the subset of the variables that are relaxed to be continuous, Xinsj (i, n, s, j) ∈ (I, Ni , Li , Ji )LP is determined in every iteration (hence dependent on π) and potentially change over the iterations, so we cannot compare the π-specific LR and HR bounds with the LP relaxation bound of the original problem, v(LP ). The challenging part of the HR formulation is the definition of set (I, Ni , Li , Ji )LP in constraints (46). If we define this set to be empty, then we have v(LR(π)) = v(HR(π)) for all π and hence even their Lagrangian relaxation duals would be the same. On the other hand, if this set has all possible elements (all Xinsj variables are non-integer), the HR bound will be closer to the LP bound (if the same variables are kept continuous). We address the challenge of determining this 20

subset by deciding it dynamically (changing it over the iterations), depending on the tolerance parameter and tightness of the bounds. However, even with changing subsets of variable relaxations, we can show that the bound obtained through HR is weak as compared to the LR bound, but, as mentioned above, we control the degradation of the lower bound. Recall that, Algorithm 2 solves knapsack problem Pins to optimality if it has potential to improve the incumbent for facility i, U Bi∗ (that is, fis + LB(Pins ) < U Bi∗ ). We apply the HR scheme at this point stating that if the potential improvement by a knapsack is less than a tolerance, denoted by η, (that is, U Bi∗ − (fis + LB(Pins )) < η), then we accept the lower bound solution (which is coming from continuous relaxation of selected Xinsj variables in this knapsack problem) and use it as the value of problem Pins . Through η, we control the amount by which the lower bound value deteriorates: Proposition 4. The worst case performance of HR with respect to LR is v(LR(π)) − v(HR(π)) ≤ η|I| Proof. Proof of Proposition 4. If facility i is open with the lower bound solution instead of its optimal solution, then the loss is v(Pins ) − LB(Pins ). We know that U Bi∗ − (fis + LB(Pins )) < η (from the definition of HR), and fis + v(Pins ) ≤ U Bi∗ (from optimality). Hence v(Pins ) − LB(Pins ) ≤ U Bi∗ − (fis + LB(Pins )) < η. Then, if all facilities are open, and if all of their solutions come from the LP relaxations of the knapsack problems, the total loss in the overall lower bound across all facilities is bounded by η|I|. Note that the worst case loss given in Proposition 4 is very pessimistic, since (1) only a small subset of the facilities are open over all candidates, and (2) the difference LB(Pins ) − v(Pins ) is much smaller than η in general. Hence, HR performs much better in practice. We compare HR with LR in Section 4.3. 4.2. Upper Bound In the previous section, we propose methods to find lower bounds for CBP . The next step in our solution methodology is to find an upper bound for the overall problem. Later in Section 4.3, we show how we combine lower and upper bound methods in an iterative subgradient algorithm. At the end of the lower bounding algorithm, we obtain a solution for the relaxed problem, which is infeasible since some customers may be assigned to multiple facilities whereas others may not be assigned to any facility at all (note that if the lower bound solution is feasible, then it is optimal). Our purpose in the upper bounding methodology is to use the information obtained from this infeasible lower bound solution and develop a feasible, low-cost solution and the method is based on two steps: (1) feasibility restoration, and (2) improvement search. 21

4.2.1. Step 1. Feasibility Restoration In this step we develop a feasible solution, by modifying the lower bound solution. We first categorize the customers across the following three groups: i. ii. iii.

P

i∈I

P

i∈I

P

i∈I

P

n∈Ni (αj )

P

n∈Ni (αj )

P

n∈Ni (αj )

P

P P

s∈Li

Xinsj = 1, customers assigned to a single facility only.

s∈Li

Xinsj > 1, customers assigned to multiple facilities.

s∈Li

Xinsj = 0, customers not assigned to any facility.

The customers in the first group are already feasible for the assignment constraints, hence we keep their assignments as they are in the lower bound solution. The customers in the second group cause infeasibility because of the multiple assignments. To overcome the infeasibility, for each customer j in this group, we choose the best facility among the facilities that the customer j is assigned, Iˆj . We first calculate the cost of assignment cˆij for the facilities in Iˆj . We calculate cˆij by considering the benefit of cancelling the assignment (i, j). Note that cˆij may have different components: • If assignment (i, j) is cancelled, we save the transportation cost; hence, cˆij + = c¯ij • According to the other customers that are assigned to i, Jˆi , and the capacity usage at facility P i, j∈Jˆi dj , cancelling the assignment (i, j) may lead to reducing the stock level and, hence, benefiting from the holding cost. Suppose that the capacity of facility i is µ ¯iˆnsˆ, where biˆn ≥ P αj , ∀j ∈ Jˆi , and µ ¯iˆnsˆ ≥ j∈Jˆi dj . When the (i, j) assignment is cancelled, we recalculate the required base fill rate level bi˜n as minn∈Ni bin ≥ αj ∀j ∈ Jˆi \ j (if j is the only customer such that biˆn = αj , then required base fill rate level reduces), and required stock level s˜ as P mins∈Li˜n µ ¯i˜ns ≥ j∈Jˆi \j dj . Then, we update cˆij as cˆij + = hi (ˆ s − s˜). • If j is the only customer that facility i serves, when the assignment (i, j) is cancelled, we also close facility i, hence, cˆij includes fixed facility opening cost recovered by closing facility i, cˆij + = fi .

Then, we keep the assignment (i∗ , j) for customer j, where cˆi∗ ,j = mini∈Iˆj cˆij , and cancel the assignments (i, j) for all i ∈ Iˆj \ i∗ . The customers in the third group are handled similarly. Since they are not assigned to any facility, we consider all facilities in Ij as candidates for these customers, and calculate the cost of assignment cˆij as follows: • For candidate facility i, we have the transportation cost, hence cˆij + = c¯ij • According to the other customers that are assigned to i, Jˆi , and the capacity usage at facility P i, j∈Jˆi dj , assigning customer j to facility i may require increasing the base fill rate and 22

stock level at the facility. Note that the capacity of facility i is µ ¯iˆnsˆ, where biˆn ≥ αj ∀j ∈ Jˆi P and µ ¯iˆnsˆ ≥ j∈Jˆi dj . When assignment (i, j) is included, we recalculate the required base fill rate level bi˜n as minn∈Ni bin ≥ αj ∀j ∈ Jˆi ∪ j, and required stock level s˜ as mins∈Li˜n µ ¯i˜ns ≥ P ˆij as cˆij + = hi (˜ s − sˆ). j∈Jˆi ∪j dj . Then, we update c

• If facility i is not open, hence j is the only customer assigned to it, cˆij includes fixed facility opening cost too, cˆij + = fi .

Then, we assign each of the customers in this group to their lowest cost facility i∗ , where cˆi∗ ,j = mini∈Ij cˆij . 4.2.2. Step 2. Improvement Search Note that in the feasibility restoration step, we do not check whether the customers in the first group have lowest cost allocations or not. In addition, we evaluate each customer sequentially, based on the calculations performed on previous customers. Here, we search all possible reallocations of the customers to improve the upper bound solution. We calculate the cost cˆij of assigning each customer j to facilities in Ij as described in the feasibility restoration step, and then choose the lowest cost facility among them. Note that when we reallocate customer j from one facility to another, say from i1 to i2 , the total cost will improve as much as cˆi1 j − cˆi2 j . We repeat this search until no more improvement is possible through customer reassignment. Note that cˆij reflects the change in the cost if assignment (i, j) is canceled (or added), but it does not reflect the cost of customer j to the overall system. For example if only two customers j1 and j2 are assigned to facility i with one unit of stock, canceling assignment (i, j1 ) will save cˆij1 = c¯ij1 and canceling assignment (i, j2 ) will save cˆij2 = c¯ij2 . However, canceling both customers j1 and j2 will save fi + hi + c¯ij1 + c¯ij2 . Hence, we name this type of cost calculation as “partial cost” (considering customer j partially, assuming all other customers are fixed in their assignments). To have a more “system-oriented” view, we define another cost term, called “marginal cost,” which is calculated considering the share of each customer in the total system cost. To do this, we P calculate the cost of a unit demand at each facility as pi = fiˆs / j∈Jˆi dj (note that the facilities with lower “utilization” have higher cost per unit demand). Then, the marginal cost of customer j is calculated as cˇij = (cij + pi )dj . Marginal and partial costs have opposite influences: We obtain P P P P P the total system cost by summing the marginal costs i∈Iˆ j∈Jˆi cˇij = i∈Iˆ fiˆs + i∈Iˆ j∈Jˆi cij dj , but when we move customer j, say from facility i1 to facility i2 , the difference between marginal costs of customer j at i1 and at i2 (ˇ ci1 j − cˇi2 j ), does not give the change in the total cost. Therefore, while decreasing the marginal cost of a customer, we may increase the total system cost. By using marginal costs in our improvement step, we benefit by moving customers from facilities with low utilization (high marginal cost) to facilities with high utilization (low marginal cost). Hence, this approach has tendency to decrease the stock levels and number of facilities. Note that neither cost calculation method is superior to the other, but both have benefits. Therefore, we use both methods, depending on the situation (see Algorithm 3 for more details). 23

Note that our upper bound methodology handles customers sequentially, one at a time; hence, the evaluation for a customer is based on the calculations and assignments made up to this customer. Therefore, the order in which we evaluate customers affects the solution. Our experiments showed that considering the customers with lower service level requirements (αj ) and high demands (dj ) first has good overall performance. 4.3. Subgradient Optimization Algorithm In this section, we present the overall subgradient-based algorithm (Algorithm 3) which uses the lower and upper bound techniques. Recall that we have two lower bound methods based on different relaxation techniques (LR and HR), each decomposed into facilities having lower bound cost of LBi . Later in Section 4.4, we compare the performances of the subgradient algorithm with each lower bound technique. We go into feasibility restoration of the lower bound solution (and calculate upper bounds) in a subgradient iteration only if the actual cost of the current lower bound solution, denoted as P P P ˆ V , is less than 2U B ∗ (see Step 14 in Algorithm 3). Here, V = s∈Li fis Zins + i∈I n∈Ni P P P P ˆ insj . The idea here is not to spend time for upper bound calcucij dj X i∈I

n∈Ni (αj )

s∈Li

j∈Ji

lations when V > 2U B ∗ , expecting that the upper bound algorithms will not be able to reduce the cost of the relaxed solution by more than its half. In an iteration in which the feasibility restoration condition is satisfied, we first obtain a feasible solution from the lower bound solution by feasibility restoration. If the feasible upper bound U B is less than 1.5 times best upper bound, then we spend more effort on this upper bound to improve its value further through an improvement search. We stop the subgradient search when we reach any of the following stopping conditions (listed with how they are represented in Table 4): • When the number of iterations reaches the maximum allowed (i.e., 500 iterations; shown as “iter ” in Table 4). • When σ is too low (i.e., if σ < 0.001; shown as “low σ” in the table). • When the gap between best upper and lower bound is small (i.e., if (U B ∗ −LB ∗ )/LB ∗ < 0.003; shown as “gap” in the table). • When the run time reaches the limit (i.e., 3600 seconds; shown as “time” in the table). 4.4. Comparison of Hybrid and Lagrangian Relaxations In this section, we conduct experiments to compare the lower bound techniques proposed in Sections 4.1.1 and 4.1.2. We test Algorithm 3 with Lagrangian relaxation and Hybrid relaxation using the same data sets described for comparing formulations F 1 and F 2 in Section 3.5. The tolerance value η is fixed to be $25 in the HR-based lower bound calculation. 24

We provide lower and upper bound values, the number of subgradient iterations performed, the number of knapsack problems solved, the reason for terminating the algorithm, and the total solution time for both lower bounding techniques in Tables 4 and 5. In the last two columns of the tables, we calculate the percentage difference between lower bounds 100(LBLR − LBHR )/LBLR and upper bounds 100(U BLR − U BHR )/U BLR . With the improvements made in Algorithm 2 (instead of using Algorithm 1), the number of knapsack problems solved in the lower bound calculations is reduced dramatically, e.g., if we were to solve instance a36 using Algorithm 1 as subroutine, more than 5000 knapsack problems need to be solved at a single iteration (hence, the instance needs about 2,500,000 knapsack problems to be solved within 500 iterations). With Algorithm 2, we need to solve about 55 knapsack problems per iteration with LR, and 30 knapsack problems per iteration with HR, which means we save about 99% of the effort. Note that the LR based lower bound algorithm still requires too many knapsack problems to be solved, resulting 28 instances not completing 500 iterations before reaching the run time limit (hence stopping criteria is time for these instances, e.g., a24-a36). However, HR solves even fewer knapsack problems, so that only 12 instances terminate because of the time limit (e.g., a31-a36). Resulting from the reductions in the number of knapsack problems solved, the HR-based lower bound algorithm has an average time per iteration of about 3 seconds, as compared with 7.5 seconds for the LR-based algorithm. Hence we can conclude that the HR-based algorithm is about 2.5 times faster than the LR-based algorithm, which especially shows the benefits with large instances. Note that, in all of the large instances (a25-a36, b25-b36), HR has both better lower bounds (up to 3.2%) and also better upper bounds (up to 7.1%). Even for the small and medium instances HR performs as good as LR; and, in the worst case (see instance a15), the HR upper bound is 0.9% worse than the LR upper bound. Note that LR spends 3072 seconds, while HR spends only 1011 seconds on this instance. Overall, HR gives either better bounds in the same or shorter time, or it gives bounds very close to those of LR within much shorter times. Hence, we conclude that HR is superior to LR. 4.4.1. Comparison of F1 and F2 vs HR In Tables 6 and 7, we pick the best solution between F 1 and F 2 (so called best of formulations as BF ), and compare it with the HR solution for each instance. In these tables, we provide best ˆ the total stock levels upper bound solutions, their solution times, the number of open facilities |I|, ˆ and the upper bound gaps as (100(U BHR − U BBF )/U BBF ). |S| For small instances, the upper bounds given by BF and HR are very close and, in almost all of these instances (other than b3 and b12), the HR upper bound is worse than the BF upper bound with a gap always less than 1%. For medium size instances, there is no winner in terms of upper bounds; HR has better upper bounds for some instances (up to 2.5%, for a18) and BF has better 25

upper bounds for some others (up to 1.3%, for a15). However, in terms of solution times, HR performs much better, and, on average, HR spends 777 seconds and BF spends 2450 seconds per medium size instance. For large instances, HR gives better upper bounds in almost all instances (other than a31 with a gap of 0.1%, and b25 with a gap of 0.6%), and improves the upper bound of BF up to 41.2%. This shows that HR excels in large scale problems with additional savings in computation times.

5. Insights and Conclusions In this paper, we study the integrated location and inventory problem with customer-based service levels. We consider inventory decisions (stock levels) and their costs explicitly in a network design model, thus making service levels (part availabilities and fill rates) vary across facilities to achieve each customer’s service level requirement. We first provide different MIP formulations, applicable to small and medium size problems. We compare these formulations and highlight the advantages of each. For large problems, we develop a Lagrangian relaxation based algorithm, enhanced with improvement techniques. We further improve the Lagrangian relaxation based algorithm, by applying the so-called “hybrid relaxation,” which uses the LP lower bounds of the knapsack sub-problems when possible to attack truly large scale problems. Our results show the effectiveness of the overall approach producing provably near optimal solutions with relatively short computation times. Although our focus has been on the methodological development for the customer-based SPL problem, we highlight several insights from the computational study here. Analyzing mainly the HR upper bound solutions and their values in Tables 4 - 7, we conclude the following: 1. As expected, larger networks with more candidate facilities and more customers to serve, higher demand rates (less reliable parts), more expensive parts, and shorter service time windows lead to higher-cost solutions for given customer-based service levels, regardless of the formulation or the relaxation used. As the results show, this is mainly due to the number of facilities that need to open and the stock levels kept at these facilities, both being on the high side with such changes in network, demand, cost and time window levels. 2. In general, for a given problem size, service time window, and holding cost, proposed solutions’ total costs increase less than the rate of increase in demand rates, observed while comparing the low and high demand levels. That is, total system costs do not even double even though demand rates more than double between the low and high demand levels. Hence, proposed solutions for high demand levels somehow lead to economies of scale in dealing with higher demand levels. This is also due to the behavior of the open locations and stock levels, as the proposed solution does not have to double the number of facilities or stock levels to deal with the increased demand rates. 3. A similar observation is made in increasing holding costs: Unit holding costs quadruple between the low and high levels, but the costs increase at a lower rate. The main interaction between location and inventory decisions is observed with changing holding costs. Increasing 26

the holding costs causes the model to more carefully open the facilities, sometimes leading to the same or fewer open facilities and definitely lower stock levels to save inventory holding costs. The overall costs go up due to higher unit holding costs, but we end up with a smaller network and lower total stock levels. 4. Service time window’s effect on costs is less predictable, in the sense that having a longer time window to satisfy a given set of service levels typically produces lower cost solutions, but does not universally. Even when the total cost is lower, the savings is a lot smaller than the change in the time window suggests. Due to the tradeoffs between facility and transportation costs, the instances with longer time windows for given service levels may lead to solutions with very similar costs or with very small relative cost savings. For example, the long time window solutions suggest closing off very few of the open facilities and reducing a handful of stock levels relative to the short time window solutions. This means that to reasonably serve the customers within 4 hours instead of, say, 2 hours, we still have to establish a sizable network with a sizable stock of parts. There are many future research endeavors one can take after this study. One avenue for future research involves considering different service level definitions such as service requirements for customer-groups, regions, or the whole network. One can also consider these service definitions simultaneously, modeling critical customers/items (requiring specific service levels) and non-critical (requiring service levels for a group of customers/items). Another future research plan may address interesting extension of the problem here to include multiple parts, multiple products, productbased service levels and their allocation to individual parts, along with part commonality and potentially allowing inventory sharing among facilities. Maybe one of the most natural and potentially most challenging extensions of the current study is to explicitly consider customer prioritization and inventory rationing (and stock thresholds or reserves) for high service-level customers. Even when a single facility faces with only two types of Poisson-demand customers (say high and low service levels), expressing the service level for the high-service level customer type in closed form even for a given stock level seems challenging [31]. Hence simultaneous optimization of facility locations, demand allocations, and stock levels along with stock thresholds for customer classes would be an extremely challenging, if not impossible, problem, just from a modeling point of view. Solving the modeled problem in a reasonable time will be another challenge. Maybe, one practical way to tackle this extension is to consider it as a post-processing step after a method (such as ours) which does not allow stock thresholds. Given that with stock thresholds allowed, we can carefully lower the stock levels and increase the stock thresholds from 0, to lower the inventory costs, reduce service to low service-level customers and maintain service to the high service level customers. Of course, this would not jointly optimize the network, stock levels and stock thresholds. We would not be able to observe, for example, the effects of stock thresholds on the facility locations. Towards that end, a potential first step could be to consider a more simplified setup, say two customer classes and a simple stock threshold policy (say a threshold of a constant level of 1 or 2 units of stock for all parts and stock levels), and approximate the service level for a handful of stock levels and stock thresholds. 27

28

ins a1 a2 a3 a4 a5 a6 a7 a8 a9 a10 a11 a12 a13 a14 a15 a16 a17 a18 a19 a20 a21 a22 a23 a24 a25 a26 a27 a28 a29 a30 a31 a32 a33 a34 a35 a36

large

medium

small

size

high

low

high

low

high

low

d

4

2

4

2

4

2

4

2

4

2

4

2

w (hrs)

h ($) 250 500 1000 250 500 1000 250 500 1000 250 500 1000 250 500 1000 250 500 1000 250 500 1000 250 500 1000 250 500 1000 250 500 1000 250 500 1000 250 500 1000

LB 30160 35392 45657 28867 32856 39817 52740 61098 77198 52705 60901 75436 53362 62264 79076 50524 57737 71021 96085 110757 138652 94992 109273 135675 81490 94505 117599 80748 93243 113132 148733 175666 223672 148517 171301 216875

UB 30210 35491 45980 28919 32951 39932 52775 61275 78275 52775 61179 75734 53520 62769 81085 50693 58274 72389 96226 111206 140778 95354 110051 137562 83085 98880 126452 85722 100067 123096 151975 181803 235713 152729 187401 236969

Iter 72 135 500 183 311 123 61 88 500 174 500 500 150 500 500 474 500 500 193 500 429 500 500 423 187 160 127 139 116 116 113 103 101 96 83 80

LR Time 30 65 407 80 172 90 33 68 549 109 546 835 347 2203 3072 1165 2094 2772 908 3065 3600 2482 3267 3600 3600 3600 3600 3600 3600 3600 3600 3600 3600 3600 3600 3600 Reason gap gap iter gap gap gap gap gap iter gap iter iter gap iter iter low σ iter iter gap iter time iter iter time time time time time time time time time time time time time

# Knp 437 993 6344 1144 2465 1348 520 1099 8838 1576 7868 12341 1735 11338 15764 5412 9939 12807 4710 15818 18743 11412 15055 16754 4883 4921 5084 4715 4762 4695 5007 4998 4944 4274 4411 4422

LB 30141 35377 45640 28841 32801 39839 52747 61103 77167 52678 60874 75329 53339 62146 78892 50483 57687 70885 95997 110648 138371 94881 109211 135406 81530 94721 119602 81316 93987 115011 149540 175890 224752 149484 175897 223793

UB 30221 35480 45982 29002 32951 39906 52775 61275 78275 52775 61179 76102 53499 62778 81795 50798 58312 72389 96345 111711 140673 95171 109903 137562 82660 97567 124568 82944 97482 119915 150713 180917 234464 151046 181081 236319

Table 4: Comparison of LR and HR with data set a Iter 80 184 490 500 374 312 45 144 482 83 500 500 368 500 500 500 500 500 500 500 500 500 500 500 500 500 500 500 500 500 411 213 167 321 195 134

HR Time 14 40 179 30 29 41 12 73 326 28 314 342 277 558 1011 220 253 364 493 1506 1762 539 1477 1635 1130 1370 2541 1010 1515 2002 3600 3600 3600 3600 3600 3600 Reason gap gap low σ iter low σ gap gap gap low σ gap iter low σ gap iter iter iter iter iter iter iter iter iter iter iter iter iter iter iter iter iter time time time time time time

# Knp 179 514 2516 106 217 393 165 1131 4996 369 4153 4461 1120 2483 4851 539 668 1179 2119 7437 8802 1769 6241 7056 579 885 2480 317 812 1465 4201 4707 4754 3289 3844 4168

LB Gap 0.1 0.0 0.0 0.1 0.2 -0.1 0.0 0.0 0.0 0.1 0.0 0.1 0.0 0.2 0.2 0.1 0.1 0.2 0.1 0.1 0.2 0.1 0.1 0.2 0.0 -0.2 -1.7 -0.7 -0.8 -1.7 -0.5 -0.1 -0.5 -0.7 -2.7 -3.2

UB Gap 0.0 0.0 0.0 -0.3 0.0 0.1 0.0 0.0 0.0 0.0 0.0 -0.5 0.0 0.0 -0.9 -0.2 -0.1 0.0 -0.1 -0.5 0.1 0.2 0.1 0.0 0.5 1.3 1.5 3.2 2.6 2.6 0.8 0.5 0.5 1.1 3.4 0.3

29

ins b1 b2 b3 b4 b5 b6 b7 b8 b9 b10 b11 b12 b13 b14 b15 b16 b17 b18 b19 b20 b21 b22 b23 b24 b25 b26 b27 b28 b29 b30 b31 b32 b33 b34 b35 b36

large

medium

small

size

high

low

high

low

high

low

d

4

2

4

2

4

2

4

2

4

2

4

2

w (hrs)

h ($) 250 500 1000 250 500 1000 250 500 1000 250 500 1000 250 500 1000 250 500 1000 250 500 1000 250 500 1000 250 500 1000 250 500 1000 250 500 1000 250 500 1000

LB 31914 37371 47964 28835 32619 39598 54467 63226 79750 53341 61324 75751 51336 59618 75576 48953 56173 69867 90764 105396 133336 90769 105294 130250 81303 94270 117819 78626 93363 116141 147318 172082 222093 150241 168030 218334

UB 32006 37760 49016 28908 32658 39686 54629 63415 81275 53520 61484 77347 51670 60433 77374 49068 56400 71614 91032 105577 135277 91038 105565 133406 82894 97852 127970 90799 97831 127582 156129 185253 230676 152534 192676 236130

Iter 236 500 500 130 112 194 61 185 457 500 242 500 500 500 500 228 500 500 128 474 457 217 464 444 190 157 122 133 126 112 112 97 101 105 92 86

LR Time 131 288 431 48 57 90 36 216 757 373 229 711 1505 2384 2837 517 1661 2307 526 2837 3600 1054 3128 3600 3600 3600 3600 3600 3600 3600 3600 3600 3600 3600 3600 3600 Reason gap iter iter gap gap gap gap gap low σ iter gap iter iter iter iter gap iter iter gap gap time gap gap time time time time time time time time time time time time time

# Knp 2012 4987 6671 785 839 1485 570 3935 12190 5476 3366 10120 7590 12038 15205 2417 8204 10738 2704 14735 18635 4854 14596 17258 4885 5042 5078 4874 4779 4751 5118 5144 4986 4316 4523 4438

LB 31863 37360 47894 28824 32572 39568 54476 63215 79597 53313 61316 75629 51265 59470 75405 48922 56131 69806 90673 105240 133170 90670 105151 130081 81370 94439 119069 81179 93876 117558 150464 175865 222439 150431 175641 222206

UB 32006 37871 49016 28908 32656 39686 54629 63415 80889 53520 61467 77161 51650 60424 77073 49079 56476 71808 90944 105961 135032 90926 105584 132658 82385 96619 124797 82236 96443 123148 151501 177889 230194 152062 179023 231705

Table 5: Comparison of LR and HR with data set b Iter 500 500 394 139 149 99 215 500 500 500 331 500 500 500 500 500 500 500 240 500 500 226 500 500 500 500 500 500 500 500 500 405 179 371 220 140

HR Time 135 198 185 13 23 26 61 288 414 198 153 232 128 191 446 286 385 781 260 1015 1991 362 1248 1458 848 1207 3431 998 1296 2360 3600 3600 3600 3600 3600 3600 Reason iter iter low σ gap gap gap gap iter iter iter gap iter iter iter iter iter iter iter gap iter iter gap iter iter iter iter iter iter iter iter time time time time time time

# Knp 1906 2904 2790 123 278 336 868 4494 6481 2676 2073 3148 193 526 1840 969 1430 3197 1158 4888 10009 1392 5290 6098 150 750 3953 150 616 2007 3854 4125 4753 3378 3899 4209

LB Gap 0.2 0.0 0.1 0.0 0.1 0.1 0.0 0.0 0.2 0.1 0.0 0.2 0.1 0.2 0.2 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 -0.1 -0.2 -1.1 -3.2 -0.5 -1.2 -2.1 -2.2 -0.2 -0.1 -4.5 -1.8

UB Gap 0.0 -0.3 0.0 0.0 0.0 0.0 0.0 0.0 0.5 0.0 0.0 0.2 0.0 0.0 0.4 0.0 -0.1 -0.3 0.1 -0.4 0.2 0.1 0.0 0.6 0.6 1.3 2.5 9.4 1.4 3.5 3.0 4.0 0.2 0.3 7.1 1.9

30

ins a1 a2 a3 a4 a5 a6 a7 a8 a9 a10 a11 a12 a13 a14 a15 a16 a17 a18 a19 a20 a21 a22 a23 a24 a25 a26 a27 a281 a291 a301 a31 a32 a331 a341 a351 a361

large

medium

small

size

high

low

high

low

high

low

d

4

2

4

2

4

2

4

2

4

2

4

2

w (hrs)

h ($) 250 500 1000 250 500 1000 250 500 1000 250 500 1000 250 500 1000 250 500 1000 250 500 1000 250 500 1000 250 500 1000 250 500 1000 250 500 1000 250 500 1000

Best UB 30210 35460 45960 28919 32889 39906 52775 61275 78016 52775 61079 75628 53499 62749 80723 50602 57939 74209 96196 111332 142984 95116 111791 137158 82899 118778 154298 120452 149273 203801 150592 187268 283862 195209 234224 346458

F1 and F2 Time Status 1 Opt 2 Opt 4 Opt 18 Opt 53 Opt 57 Opt 1 Opt 6 Opt 163 Opt 11 Opt 604 Opt 2769 Opt 36 Opt 1004 Opt 3644 N/F 1252 Opt 3620 N/F 3615 N/F 45 Opt 3600 N/F 3600 N/F 968 Opt 3600 N/F 3600 N/F 3600 N/F 3600 N/F 3600 N/F 3600 N/F 3600 N/F 3600 N/F 3600 N/F 3600 N/F 3600 N/F 3600 N/F 3600 N/F 3600 N/F ˆ |I| 11 11 11 6 6 5 11 11 11 11 11 9 16 16 16 11 10 10 19 18 19 18 17 13 21 29 26 35 42 52 31 33 43 68 64 62 |ˆ s| 21 21 21 16 15 14 34 34 33 34 32 28 37 37 34 31 28 27 61 58 55 61 59 52 59 78 78 63 77 86 115 114 105 162 158 135 UB 30221 35480 45982 29002 32951 39906 52775 61275 78275 52775 61179 76102 53499 62778 81795 50798 58312 72389 96345 111711 140673 95171 109903 137562 82660 97567 124568 82944 97482 119915 150713 180917 234464 151046 181081 236319

HR Time 14 40 179 30 29 41 12 73 326 28 314 342 277 558 1011 220 253 364 493 1506 1762 539 1477 1635 1130 1370 2541 1010 1515 2002 3600 3600 3600 3600 3600 3600

Table 6: Comparison of best of F 1 and F 2 vs HR with data a ˆ |I| 11 11 11 6 6 5 11 11 11 11 10 9 16 16 16 11 10 9 19 18 17 17 14 13 19 18 16 19 18 12 32 29 22 30 28 24 |ˆ s| 21 21 21 16 15 14 34 34 34 34 33 27 37 37 37 31 30 27 61 59 57 60 56 54 58 55 54 58 57 46 118 113 102 116 113 100 UB Gap 0.0 0.1 0.0 0.3 0.2 0.0 0.0 0.0 0.3 0.0 0.2 0.6 0.0 0.0 1.3 0.4 0.6 -2.5 0.2 0.3 -1.6 0.1 -1.7 0.3 -0.3 -17.9 -19.3 -31.1 -34.7 -41.2 0.1 -3.4 -17.4 -22.6 -22.7 -31.8

31

ins b1 b2 b3 b4 b5 b6 b7 b8 b9 b10 b11 b12 b13 b14 b15 b16 b17 b18 b19 b20 b21 b22 b23 b24 b25 b261 b271 b281 b291 b301 b31 b32 b331 b341 b351 b361

large

medium

small

size

high

low

high

low

high

low

d

4

2

4

2

4

2

4

2

4

2

4

2

w (hrs)

h ($) 250 500 1000 250 500 1000 250 500 1000 250 500 1000 250 500 1000 250 500 1000 250 500 1000 250 500 1000 250 500 1000 250 500 1000 250 500 1000 250 500 1000

Best UB 31983 37521 48521 28881 32656 39686 54629 63415 80467 53411 61369 76269 51519 59807 76830 48994 56520 73142 90890 105503 136857 90884 105705 134057 81872 130165 166456 118857 154471 206096 154180 181316 282311 208510 251581 354921

F1 and F2 Time Status 3 Opt 2 Opt 4 Opt 8 Opt 9 Opt 25 Opt 2 Opt 15 Opt 238 Opt 27 Opt 281 Opt 3600 N/F 1894 Opt 1367 Opt 3600 N/F 360 Opt 3600 N/F 3600 N/F 190 Opt 2394 Opt 3600 N/F 2409 Opt 3600 N/F 3600 N/F 3600 N/F 3600 N/F 3600 N/F 3600 N/F 3600 N/F 3600 N/F 3600 N/F 3600 N/F 3600 N/F 3600 N/F 3600 N/F 3600 N/F ˆ |I| 11 11 11 7 6 5 13 12 13 11 11 9 14 14 14 11 11 11 18 17 19 16 16 16 18 39 39 35 46 47 33 29 46 80 59 55 |ˆ s| 23 22 22 17 15 14 37 35 33 33 31 27 34 33 32 30 30 28 60 57 56 59 57 52 55 70 67 76 69 67 113 102 105 180 117 108 UB 32006 37871 49016 28908 32656 39686 54629 63415 80889 53520 61467 77161 51650 60424 77073 49079 56476 71808 90944 105961 135032 90926 105584 132658 82385 96619 124797 82236 96443 123148 151501 177889 230194 152062 179023 231705

HR Time 135 198 185 13 23 26 61 288 414 198 153 232 128 191 446 286 385 781 260 1015 1991 362 1248 1458 848 1207 3431 998 1296 2360 3600 3600 3600 3600 3600 3600

Table 7: Comparison of best of F 1 and F 2 vs HR with data b ˆ |I| 11 11 12 6 6 5 13 12 12 11 11 9 14 14 14 11 10 8 18 17 17 17 16 12 18 18 17 17 17 13 29 24 22 27 24 19 |ˆ s| 23 23 22 15 15 14 37 35 34 34 31 29 35 35 33 30 29 26 61 58 57 60 58 50 58 57 55 57 55 51 111 102 99 109 102 97 UB Gap 0.1 0.9 1.0 0.1 0.0 0.0 0.0 0.0 0.5 0.2 0.2 1.2 0.3 1.0 0.3 0.2 -0.1 -1.8 0.1 0.4 -1.3 0.0 -0.1 -1.0 0.6 -25.8 -25.0 -30.8 -37.6 -40.2 -1.7 -1.9 -18.5 -27.1 -28.8 -34.7

Appendix A. Algorithms Algorithm 1 Solving Pi 1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12:

for (n, s) ∈ (Ni , Lin ) do SOLVE Pins for j ∈ Jin Construct the solution set as Jˆins = {j ∈ Jin : Xinsj = 1} end for Find (n∗ , s∗ ) such that V (Pi,n∗ ,s∗ ) = minn∈Ni ,s∈Lin V (Pins ) if fis + V (Pi,n∗ ,s∗ ) ≤ 0 then V (Pi ) = fi,s∗ + V (Pi,n∗ ,s∗ ) Zi,n∗ ,s∗ = 1 . Open facility i with fill rate bi,n∗ and stock level s∗ Xi,n∗ ,s∗ ,j = 1 for all j ∈ Jˆin∗ s∗ . Assign customers to facility i else Set V (Pi ) = 0 . Do not open facility i end if

In lines 1-4 of Algorithm 1 we solve a knapsack problem for each base fill rate and stock level, and we perform the required updates in lines 5-10 if the facility is open. If the facility is not open (by setting all Xinsj and Zins variables for this facility to 0) we set its value to zero in line 11. Algorithm 2 Improved algorithm to solve problem Pi 1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12: 13: 14: 15:

U Bi∗ = 0 , (s∗ , n∗ ) = (∅, ∅) Ji− = {j ∈ Ji : c¯ij < 0} Calculate Ni− and L− in − − Ji = {j ∈ Ji : c¯i[k] ≤ c¯i[l] ∀k ≤ l} − Jin = {j ∈ Ji− : αj ≤ bin } A = {(n, s) : n ∈ Ni− , s ∈ L− in } − − for n ∈ Ni , s = Lin do P V (Pins ) = j∈J − c¯ij in if fis + V (Pins ) < U Bi∗ then U Bi∗ = fis + V (Pins ) s∗ = s and n∗ = n − Jˆins = Jin end if A = A \ (n, s) end for

. Initialize incumbent solution . Set of customers having negative cost . Order the customers with respect to their costs . Set of active knapsack problems to be solved . recall that µins ≥

P

− j∈Jin

dj for s = L− in

. update incumbent solution . record the customers in the solution of (i, n, s) . FATHOM

In line 1 of Algorithm 2, we initialize the incumbent knapsack solution to be 0, since v(Pi ) can not be positive (see discussion in Proposition 3), and we also initialize stock and base fill rate levels to be empty. In lines 2-3 we calculate Ji− , Ni− , and L− in by only considering customers having negative costs. We order the customers with respect to their c¯ij /dj ratios in line 4, and construct ordered 32

Algorithm 2 Improved algorithm to solve the problem Pi (continued) 16: 17: 18: 19: 20: 21: 22: 23: 24: 25: 26: 27: 28: 29: 30: 31: 32: 33: 34: 35: 36: 37: 38: 39: 40: 41: 42: 43: 44: 45: 46: 47: 48: 49: 50: 51: 52: 53: 54:

for (n, s) ∈ A do Calculate LB(Pins ) if fis + LB(Pins ) ≥ U Bi∗ then A = A \ (n, s) . FATHOM else Calculate U B(Pins ) if (U B(Pins ) − LB(Pins ))/|LB(Pins )| <  then U Bi∗ = fis + U B(Pins ) . update incumbent solution ∗ ∗ s = s and n = n Jˆins = j ∈ U B(Pins ) . record the customers in the solution of (i, n, s) A = A \ (n, s) . FATHOM else if fis + U B(Pins ) < U Bi∗ then U Bi∗ = fis + U B(Pins ) . update incumbent solution end if end if end if end for while A 6= ∅ do for (n, s) ∈ A do if fis + LB(Pins ) ≥ U Bi∗ then A = A \ (n, s) . FATHOM end if end for − SOLVE Piˆnsˆ for (ˆ n, sˆ) : LB(Piˆnsˆ) ≤ LB(Pins ) ∀(n, s) ∈ A j ∈ Jin A = A \ (ˆ n, sˆ) . FATHOM ∗ if fis + V (Piˆnsˆ) < U Bi then U Bi∗ = fis + V (Piˆnsˆ) . update incumbent solution ∗ ∗ s = sˆ and n = n ˆ ˆ Jiˆnsˆ = j ∈ V (Piˆnsˆ) . record the customers in the solution of (i, n ˆ , sˆ) end if end while if (s∗ , n∗ ) 6= (∅, ∅) then V (Pi ) = Ui∗ Zi,n∗ ,s∗ = 1 . Open facility i with fill rate bin∗ and stock level s∗ Xi,n∗ ,s∗ ,j = 1 for all j ∈ Jˆi,n∗ ,s∗ . Assign customers to facility i else Set V (Pi ) = 0 . Do not open facility i end if 33

− sets Jin for each n ∈ Ni− in line 5. Note that while calculating lower bound (LB(Pins ) in line 17), upper bound (U B(Pins ) in line 21), and solving for optimality (in line 40) we only consider the − customers in set Jin . In line 6, we initialize the active knapsack problems to be solved for facility i.

As discussed before (in item iii above), in lines 7-15, we calculate v(Pins ) for the stock level s = L− in − (which can cover demands of all customers in Jin ) for all base fill rates. Then, we update the incumbent by picking up the best (most negative) among those. The problems that are solved in this for loop are then removed from the active set of knapsacks. In the second for loop (lines 16-33), we calculate a lower bound for each active knapsack, and, if they do not have potential to improve the current incumbent (see the condition in line 18), we fathom these knapsacks, else we calculate their upper bounds, and update the incumbent if needed. If the percentage difference between upper and lower bounds is very small (less than ), then we conclude that the upper bound solution is in fact optimal. For the same fill rate level, while the stock level is increasing, we use the smaller stock level lower bound solution as a base and add new customers for the remaining capacity. In other words, we calculate the lower bound for the knapsack problem of: min subject to

X

X

c¯ij Xinsj

(A.1)

dj Xinsj ≤ µ ¯ins − µ ¯in,s−1

(A.2)

− Xinsj = 0 or 1, ∀j ∈ J˜in

(A.3)

− j∈J˜in

− j∈J˜in

− where J˜in = {j ∈ Jin \ Jˆin,s−1 } and Jˆin,s−1 are customers in the solution of the knapsack problem with the smaller stock level (note that the last customer in the knapsack problem with stock level s − 1 can be partially assigned, and hence, J˜in will have remaining part of this customer).

In lines 35-39, we fathom the knapsacks by using their lower bounds and the current incumbent. If there are still active knapsacks, we pick a candidate problem with the lowest lower bound, and solve it to optimality by using Xpress-MP solver. Then, we update the incumbent and repeat this procedure until no active knapsack remains. In the final part of the algorithm (lines 48-54), we open facility if it is profitable, i.e., the best knapsack solution value is negative.

34

Algorithm 3 Subgradient algorithm 1: LB ∗ = 0 2: U B = f ind f easible(LB ∗ ) 3: U B1 = improve(U B, partial) 4: U B2 = improve(U B1 , marginal, partial) 5: U B3 = improve(U B, marginal, partial) 6: U B ∗ = min{U B1 , U B2 , U B3 } 7: Initialize: π, σ 8: while stopping conditions not satisfied do 9: Run Algorithm 2 ∀i ∈ I P 10: LB = i∈I V (Pi ) 11: if LB > LB ∗ then 12: LB ∗ = LB 13: end if 14: if V ≤ 2U B ∗ then 15: U B = f ind f easible(LB) 16: if U B ≤ 1.5U B ∗ then 17: if U B < U B ∗ then 18: U B1 = improve(U B, partial) 19: U B2 = improve(U B1 , marginal, partial) 20: U B3 = improve(U B, marginal, partial) 21: else 22: U B1 = improve(U B, marginal, partial) 23: end if 24: end if 25: U B ∗ = min{U B1 , U B2 , U B3 , U B ∗ } 26: end if 27: if LB ∗ not improved 20 consecutive iterations then 28: σ = 0.5σ 29: end if P P P P 30: γ := σ(U B ∗ − LB)/ j (1 − i n s Xinsj )2 31: for j ∈ J do P P P 32: πj := πj + γ(1 − i n s Xinsj ) 33: end for 34: end while

35

. returns V (Pi )

In the first part of the Algorithm 3 we have initializations: we set the initial lower bound to 0 in line 1, we calculate the initial upper bound by using several combinations of improvement methods in lines 2-6. Note that we do not have a starting lower bound solution while finding the initial upper bound in line 2; hence, no customer is considered to have an assignment in the upper bound algorithm (see Section 4.2, item iii). We refer to the feasibility restoration step as “f ind f easible(LB)” where LB stands for the lower bound solution; it returns the feasible upper bound solution. And we refer to the improvement step as “improve(U B, method),” where U B stands for the upper bound solution, and method stands for the cost calculation method (note that we use improve(U B, method1, method2) in several steps, which implies we first apply the improvement search with cost calculation method1 and then with cost calculation method2). We initialize step size reduction coefficient σ and Lagrangian multipliers πj in line 7. Initialization of πj ’s has an important role in determining a good start in the subgradient search. For the dual P problem, πj represents the profit arising from assigning customer j, and j∈J πj (dual objective) represents the total profit where the right hand sides of the assignment constraints are unit vector. P Ignoring the duality gap, we calculate πj ’s from the primal-dual optimality condition j∈J πj = U B ∗ , where U B ∗ is the initial upper bound value obtained in line 6. Hence, πj corresponds to the cost savings we can have from not assigning customer j. We use the marginal costs of customers calculated in the upper bound solution, and set dual variables (Lagrangian multipliers) as πj = cˇˆij where ˆi is the facility to which customer j is assigned in the upper bound solution. Starting from line 8 until the last step of the Algorithm, we iteratively calculate a lower bound (lines 9-13) and an upper bound (lines 14-26) for a given π, and update parameters σ, γ and π (lines 27-33).

References [1] A. Atamturk, G. Berenguer, and Z.M. Shen. A conic integer programming approach to stochastic joint location-inventory problems. Operations Research, 60(2):366–381, 2012. [2] F. Barahona and D. Jensen. Plant location with minimum inventory. Mathematical Programming, 83:101–111, 1998. [3] R.J.I. Basten and G.J. van Houtum. System-oriented inventory models for spare parts. Surveys in Operations Research and Management Science, 19:34–55, 2014. [4] M.F. Candas and E. Kutanoglu. Benefits of considering inventory in service parts logistics network design problems with time-based service constraints. IIE Transactions, 39:159–176, February 2007. [5] M.A. Cohen, P.R. Kleindorfer, and H.L. Lee. Service constrained (s,S) inventory systems with priority demand classes and lost sales. Management Science, 34(4):482–499, 1988.

36

[6] M.A. Cohen, P. Kamesam, P.R. Kleindorfer, H.L. Lee, and A. Tekerian. OPTIMIZER: IBM’s multi-echelon inventory system for managing service logistics. Interfaces, 20(1):65–82, 1990. [7] M.A. Cohen, Y-S. Zheng, and Y. Wang. Identifying opportunities for improving teradyne’s service parts logistics system. Interfaces, 29(4):1–18, 1999. [8] M.A. Cohen, C. Cull, H.L. Lee, and D. Willen. Saturn’s supply chain innovation: High value in after sales service. MIT Sloan Management Review, 41(4, Summer):93–101, 2000. [9] A.M. Cohn and C. Barnhart. Composite-variable modeling for service parts logistics. Annals of Operations Research, 145(1):383–383, Jul 2006. ISSN 1572-9338. doi: 10.1007/ s10479-006-0055-2. URL https://doi.org/10.1007/s10479-006-0055-2. [10] M.S. Daskin, C.R. Coullard, and Z.-J.M. Shen. An inventory-location model: Formulation, solution algorithm and computational results. Annals of Operations Research, 110:13–106, 2002. [11] F. Gzara, E. Nematollahi, and A. Dasci. Linear location-inventory models for service parts logistics network design. Computers and Industrial Engineering, 69:53–63, 2014. [12] V. Jeet and E. Kutanoglu. Part commonality effects on integrated network design and inventory models for low-demand service parts logistics systems. International Journal of Production Economics, 206:46 – 58, 2018. ISSN 0925-5273. [13] V. Jeet, E. Kutanoglu, and A. Partani. Logistics network design with inventory stocking for low-demand parts: Modeling and optimization. IIE Transactions, 41:389–407, 2009. [14] H. Kellerer, U. Pferschy, and D. Pisinger. Knapsack Problems. Springer, August 2005. [15] H. Mak and Z.M. Shen. A two-echelon inventory-location problem with service considerations. Naval Research Logistics, 56(8):730–744, 2009. [16] P.A. Miranda and R.A. Garrido. Incorporating inventory control decisions into a strategic distribution network design model with stochastic demand. Transportation Research Part E, 40:183–207, 2004. [17] K. Moinzadeh and H.L. Lee. Batch size and stocking levels in multi-echelon repairable systems. Management Science, 32(12):1567–1581, December 1986. [18] J.A. Muckstadt. Analysis and Algorithms for Service Parts Supply Chains. Springer, New York, NY, 2005. [19] L.K. Nozick. The fixed charge facility location problem with coverage restrictions. Transportation Research Part E, 37:281–296, 2001. [20] L.K. Nozick and M.A. Turnquist. Integrating inventory impacts into a fixed-charge model for locating distribution centers. Transportation Research Part E, 34(3):173–186, 1998. 37

[21] U. Pferschy. Dynamic programming revisited: Improving knapsack algorithms. Computing, 63(4):419–430, 1999. [22] K. Poole. Seizing the potential of the service supply chain. Supply Chain Management Review, pages 54–61, 2003. [23] J.A. Rappold and B.D. Van Roo. Designing multi-echelon service parts networks with finite repair capacity. European Journal of Operational Research, 199(3):781–792, 2009. [24] W.D. Rustenburg, G.J. van Houtum, and W.H.M. Zijm. Spare parts management at complex technology-based organizations: An agenda for research. International Journal of Production Economics, 71:177–193, 2001. [25] A. Sen, D. Bhatia, and K. Dogan. Applied Materials uses operations research to design its service and parts network. Interfaces, 40(4):253–266, 2010. [26] Z.-J.M. Shen and M.S. Daskin. Trade-offs between customer service and cost in integrated supply chain design. Manufacturing & Service Operations Management, 7(3):169 – 271, 2005. [27] Z.-J.M. Shen, C. Coullard, and M.S. Daskin. A joint location-inventory model. Transportation Science, 37(1):40–55, 2003. [28] C.C. Sherbrooke. Optimal Inventory Modeling of Systems. Wiley, New York, NY, 1992. [29] E. Topan, Z.P. Bayindir, and T. Tan. Heuristics for multi-item two-echelon spare parts inventory control subject to aggregate and individual service measures. European Journal of Operational Research, 256(1):126 – 138, 2017. ISSN 0377-2217. [30] G.J. van Houtum and B. Kranenburg. Spare parts inventory control under system availability constraints. Springer, New York, NY, 2015. [31] O. Vicil and P. Jackson. Computationally efficient optimization of stock pooling and allocation levels for two-demand-classes under general lead time distributions. IIE Transactions, 48(10): 955 – 974, 2016. [32] P. Zipkin. Fundemantals of Inventory Management. McGraw Hill - Irwin, Boston, MA, 2000.

38