The stochastic guaranteed service model with recourse for multi-echelon warehouse management

The guaranteed service model (GSM) computes optimal order-points in multi-echelon inventory control under the assumptions that delivery times can be guaranteed and the demand is bounded. Our new stochastic guaranteed service model (SGSM) with Recourse covers also scenarios that violate these assumptions. Simulation experiments on real-world data of a large German car manufacturer show that policies based on the SGSM dominate GSM-policies.


Introduction
We investigate the multi-echelon inventory control problem of a large German automobile manufacturer. The core of inventory control is to balance cost with service quality. Two main classes of mathematical models have been established in the literature: Stochastic service models (SSM) explicitly take into account the distribution in lead times and demands and account for all actions that have to be made to ful-fill demands. Guaranteed service models (GSM) assume some operational flexibility outside the model that is not accounted for, and the model itself works with bounded demands and bounded lead times instead. The GSM paradigm is motivated by the fact in most companies individuals have the capability to react to unforeseen events successfully in many ways, and no model can possibly capture all these reactions faithfully.
Research on the GSM was initiated by Simpson Simpson (1958). In Graves and Willems (2000) and Magnanti et al. (2006) the investigations were extended to tree structured networks and acyclic networks, respectively. The first application of the GSM to spare parts distribution networks was carried out in Schade (2008) for the spare parts distribution system of a large German car manufacturer. Earlier investigations in other application contexts can be found, e.g., in Inderfurth (1991Inderfurth ( , 1994; Minner (1997).
The GSM is kind of an indirect model whose decision variables are service times that each warehouse in the system (called a node) guarantees to its successor nodes. Thus, in a sense, it computes decisions in the space of event times. These guaranteed service times have to be transformed into decisions on the time at which a replenishment of a certain quantity has to be ordered at each node. Whenever one of the wide-spread (but in general suboptimal) base-stock policies (basically: order up to an stock level S when the stock falls below s) is used, this is easy: Via a bounded-demand assumption, the guaranteed service times can be transformed into minimum stock level requirements at the nodes (see, e.g., Graves and Willems (2000) and Magnanti et al. (2006)).
In Clark and Scarf (1960) base-stock policies are theoretically justified by proving their optimality in elementary special cases. In many environments, (s, S) policies or the like are prescribed despite their suboptimality because the mechanism is easy to understand. Then, only the parameters s and S can be chosen. In such environments, the GSM can be used to determine values for s and S. Finally, in Magnanti et al. (2006) it has been shown, how the GSM including a determination of order points can (approximately) be solved by a mixed integer linear program (MILP). This gives the modeler the opportunity to easily add business constraints in the space of event-times to the GSM without affecting its solvability too much. This model was used in Schade (2008) for the spare parts distribution network that is also considered in this paper.
The main drawback of the GSM is that it cannot keep control of the usage of operational flexibility. The problem is that employees can use operational flexibility (even at no cost) but not beyond a certain capacity. Operational flexibility is used in the GSM to guarantee bounded lead times and bounded demands. That is, deliveries can be expedited and/or outsourced in emergency cases. Assuming a known joint distribution for the lead times and demands, one can keep control of the amount of operational flexibility by prescribing a target service level at the nodes: for example, if we want that during at least 90 % of the time we can deliver on time and enough without using operational flexibility, then we can set up the GSM with the 90 % quantile of the lead time/demand distribution as the bounded lead time and demands. Such target service levels are meaningful at the nodes delivering to end customers (demand nodes). However, prescribed target service levels at internal nodes are hard to justify.
The questions that remain are therefore: How should one decide on the internal target service levels? Should there be at all individual target service levels for the nodes? Do individual quantile-based lead time and demand bounds in a GSM really guarantee that the target service level is achieved in the demand nodes? In a sense, the core task of inventory control (balance cost with service quality) is shifted to the formal choice of the target service levels whose side-effects and correct real-world interpretations are not always straight-forward.
One way to avoid the weaknesses of the GSM is to use a model from the SSM category instead (see, e.g., Sherbrooke (1968) for the METRIC system, Diaz and Fu (2005) for a survey, and Doǧru et al. (2005) for a special version of a stochastic service model). However, in SSM adding further restrictions, e.g., imposed by the business processes of a company, can render a method impractical. This is because many algorithms are based on specialized dynamic programming algorithms that may fail to apply to a system with additional restrictions. In contrast to this, adding restrictions to the ILP model of the GSM to a certain extent does not affect the ILP solution procedure.
To keep the advantages and overcome some of the weaknesses of the GSM, we have introduced the stochastic guaranteed service model with recourse (SGSM) in Rambau and Schade (2010) and applied the first basic version of it to the inventory control problem in a multi-echelon warehouse system of a spare part distributor. The SGSM adds a lead time and demand sampling component and a recourse component to the GSM. (See Birge and Louveaux (1997) for background on stochastic programming with recourse.) This way, each lead time and demand scenario that is covered by the sampling component is accounted for inside the model; in the basic version, generic operational flexibility allowing for smaller safety stock levels leads to additional recourse costs. The SGSM does not need any prescribed service level requirements; it yields service levels as a result of the computation. 1 However, estimates for the recourse costs have to be given for all scenarios where lead times and demands can only be handled with operational flexibility. Since we need not prescribe the service level requirements, the core task of inventory control-balancing cost and service quality-is done inside the model.
In this paper, we go beyond the conference presentation Rambau and Schade (2010) in the following aspects (among others): -We introduce the new SGSM with a non-trivial complete recourse consisting of a transportation option besides the penalty cost for unsatisfied demand, i.e., requested parts that cannot be delivered in time. -We solve the SGSM by a combination of sample average approximation with stateof-the art scenario reduction techniques. This way, a better coverage of unlikely but expensive scenarios is achieved without increasing the computation times in the MILP solver. Our new asymmetric distance function for the asymmetric scenario reduction takes into account the influence of the scenario reduction on the result of the optimization. To the best of our knowledge, this is new.
-We present a more comprehensive documentation of extended computational results, including a new comparison to one representative Doǧru et al. (2005) of the class of stochastic service models that could be implemented to cope with our test data.
In simulations on real data we observe that the new asymmetric scenario reduction technique is able to improve the approximation quality of the SGSM by a large margin. Moreover, the SGSM decreases the inventory holding cost and the recourse costs at the same time compared to the GSM. It would be interesting from a theoretical point of view to also check performances on artificial randomized data. For this work, we focused on the practical impact in realworld applications, for which randomized data is rarely representative. We emphasize that, for this reason, our simulation experiments are completely independent of the assumptions of the tested models -it rather represents our partner's process as closely as possible.
In the following section we introduce the modeling of the GSM and the SGSM before we theoretically compare them in Sect. 3. We develop methods for scenario generation and scenario reduction in Sect. 4. After the description of the simulation method and some computational results in Sect. 5 we end with our conclusions.

Modeling
In this section we first give an introduction to the GSM. We use the ILP modeling approach as in Magnanti et al. (2006). Then we present the SGSM in two different ways. First, in Sect. 2.2 we introduce the SGSM as a two stage stochastic mixedinteger linear program with simple recourse. Second, in Sect. 2.3 we show an extension where the recourse action of the locations supplying the end customers are modeled as a transportation problem. We present all models only for diverging networks, i.e., networks in which all nodes have unique predecessors.

The guaranteed-service-model
In order to make the paper self-contained, we repeat in this section the known MILP formulation for the GSM as can be found in the original work in Magnanti et al. (2006). Since in spare-part systems, there are large expensive parts stored at small stock levels, the order points are required to be integral. Notational conventions are taken from Rambau and Schade (2010). Parameters of the model GSM are: G directed graph describing the warehouse network time period that i needs to bridge with its inventory (i.e., the time between order and delivery of replenishments from the predecessors of i) y i order-point in i The model GSM now reads as follows: With the well-known multiple-choice modeling of piecewise linear functions, this non-linear model can approximately be transformed into an MILP (see Magnanti et al. (2006)).
Note that the GSM is a one-stage model although it deals with a multi-stage decision problem. In a sense, the GSM computes decisions in the space of event times: how long does it take (at most) after an order at node i has been placed by node j until the delivery of node i arrives at node j. These one-stage decisions in the space of event times imply the time x i that has to be bridged with inventory until we can be sure to receive a replenishment. The necessary safety stocks y i at all nodes i can then be computed from x i using the maximal demand quantity i (x i ) during time x i . These decisions can be considered stationary over time. They are transformed by a base-stock policy into the multi-stage sequence of decisions for every node, namely, how much to order given the current stock level.

The stochastic guaranteed-service-model with simple recourse
In this section, we present the simple-recourse version of the SGSM, following Rambau and Schade (2010). The short-comings of the GSM are addressed by turning the deterministic GSM into a stochastic model with recourse. Again, the service times are the decision variables. We first fix our notation for the stochastic data. A lead time/demand distribution consists of the following: -A finite sample set Ω = {1, 2, . . . , |Ω|} of states ω of the world, encoded by natural numbers. -A positive probability for each ω, denoted by p ω > 0, inducing a probability measure P Ω on 2 Ω via P(A) = ω∈A p ω . -For each ω ∈ Ω, a random vector of measurements ξ ω that has the following components: where L ω i ≤ 0 is the lead time in node i in state ω and Ψ ω i is a function that assigns to each time interval x a demand Ψ ω i (x), which denotes the demand presented to node i in state-of-the-world ω.
We call ξ ω the lead time/demand scenario, or scenario for short, of state ω ∈ Ω. The probability that ξ ω = ξ is equal to ω∈Ω:ξ ω =ξ p ω . The induced probability measure on 2 is denoted by P . The set of all scenarios with positive probabilities is denoted by . The complete lead time/demand rate distribution is denoted by ( , p).
We are given a finite lead time/demand distribution ( , p). Whenever a node misses its "guaranteed" service time, then a recourse action has to be taken that delivers the part faster (expediting). Whenever in some state ω of the world a node cannot deliver any of the demanded parts, a recourse action has to be taken that gets the missing parts from somewhere else (outsourcing). In contrast to the GSM, where such actions can be taken tacitly at no cost, in the SGSM operational flexibility is made explicit in the model such that we can assign a cost to it. The remaining approximation of the SGSM is that expediting and outsourcing can always be done by paying a surplus, which leads to a model with simple complete recourse. In the following, we review the formalization of the SGSM from Rambau and Schade (2010).
Scenario and recourse parameters of the SGSM: Ω index set of scenarios p ω probability of scenario ω ∈ Ω t i cost to compensate for one time unit of late delivery, non-negative q ω i recourse variable for missing pieces in scenario ω; "how many pieces should be compensated at a cost of t i per unit?" The SGSM is modeled by the following MILP: Here, constraint (12) makes sure that whenever a guaranteed service time cannot be met, we have to buy r ω i time units of expediting services. Restriction (11) makes sure that no node can expect to receive a piece faster than its predecessor guarantees. Condition (10) indicates that each node delivering to end customers fulfills the service time requirements of the end customer and finally in (13) the reorder point y i plus the outsourced quantities need to be higher than the demand during time x i at node i, as given by Ψ ω i (x i ).
Remark 1 If one wants to prescribe a service level requirement in a node, one has to replace constraints (12) and (13) by the corresponding GSM constraint with the quantile-based demand function i and without recourse variables instead of the scenario-based demand function Ψ ω i with recourse variables.
Again, a linearization of Ψ ω i (x i ) can be carried out by the multiple-choice modeling of piecewise linear functions.
Note that the SGSM is now a two-stage model although we are dealing with a multi-stage decision problem. We are still working in the space of event times and obtain safety stock levels y i for all nodes i. All decisions coming out of the SGSM are considered stationary. However, now we can detect extremal lead time and demand scenarios for which we need to apply a (simple but costly) recourse action.
In the variant of the SGSM with simple recourse, the second stage is equivalent to paying a penalty for missing parts and missed deadlines. The model in the following section contains a non-simple second stage.

Extension with external suppliers and lost sales
The model with simple recourse from the previous section can be extended by modeling an explicit recourse process. We assume that unmet customer demands are lost. However, internal orders are backlogged. The locations that deliver parts to the end customers can order parts from external suppliers to prevent lost sales.
The external suppliers deliver the parts directly to the end customers such that there is no delay in the delivery. The costs of an order from an external supplier depends on the distance between the ordering location and the supplier. Of course, the suppliers do not have unlimited stock such that capacity constraints have to be taken into account. To concentrate on these recourse actions we assume that the lead times in the system are fixed. An extension with lead time uncertainties would be straight-forward. We need some more notation to model the new situation J set of external suppliers C j capacity of the external supplier j q ω ji recourse variable for the number of parts that location i orders at supplier j c ji costs for location i to order one part from supplier j This leads us to the following model: The constraints (20) and (22) are forcing the reorder points y i to be higher than the demand during the time that has to be covered. For locations with end customer demand there is the possibility to procure parts from external suppliers. As the suppliers do not have infinity capacity, constraint (21) must hold. So far, this model does not have complete recourse. Therefore, we introduce another recourse variable. As before, we enable for every location the possibility to pay a penalty for a lost sale if it cannot deliver the ordered parts. For instance one can provide the customer with a replacement vehicle until the spare part can be delivered and the customer's car is fixed. The corresponding penalty recourse variable is denoted by q ω i , as in the first model, and the penalty costs are denoted by c i again. We obtain a two stage stochastic model with complete recourse: Note, that by using the penalty recourse variables we force complete recourse but account for failure by some cost. The computational results in Sect. 5.4 suggest that the SGSM policies with the tested penalty values dominate GSM-policies in terms of both inventory and recourse cost, not only total cost. This means, the resulting SGSM policy, internally using those successful penalty values, will perform better than the corresponding GSM policies also for any other penalty values.

Theoretical comparison of the SGSM with the GSM
One reason for the superior performance of policies based on the SGSM seems to be that it generates structurally different optimal solutions. In this section we give some theoretical evidence for this intuition. We restrict ourselves to linear demand scenarios. Thus, the stochastic data simplifies to the following: For each ω ∈ Ω, the random vector of measurements ξ ω now has the following components: where L ω i ≤ 0 is the lead time in node i in state ω and α ω i > 0 is now the constant demand rate in node i in state ω.
Note that we assume that all demand rates are positive. We do this in order to avoid special cases with limited relevance. We call ξ ω the lead time/demand rate scenario, or scenario for short, of state ω ∈ Ω. As before, the probability that ξ ω = ξ is equal to ω∈Ω:ξ ω =ξ p ω . The set of all scenarios with positive probabilities is again denoted by . The complete lead time/demand rate distribution is again denoted by ( , p).

Lead-times/demand-rates induced by a distribution and a target service level
In order to prepare for a formalization of these considerations, we elaborate in more detail on the chance-constraint interpretation of the GSM restrictions.
Note first that the conclusions of the GSM (optimality of the inventory decisions in the real world) are formally only correct if, in the real world, at each node all lead times and demand rates are bounded by tight choices of L i and α i of the GSM. In that case, a GSM solution will service all customers in time with minimal inventory.
Since this assumption is not realistic over long time-periods, consider the following chance-constraint interpretation: If at some node i ∈ N (G) the lead time of demand rate exceeds the choice of the upper bounds L i and α i chosen for the GSM, then the guaranteed service times computed by the GSM are nevertheless achieved by unspecified actions outside the model (operational flexibility). We require, however, that a prescribed fraction of the demand is handled on time inside the model. If this fraction were set to zero, one would simply shift all demands outside the model, making the GSM pointless. The positive fraction of demands required to be handled inside the model is the target service level mentioned in the introduction. It can in principle be prescribed for each node. In the following we restrict ourselves to the case of a global target service level to reduce notational clutter.
In the chance-constraint interpretation, the GSM arises as a (possibly approximative) deterministic equivalent of the following chance-constraint stochastic program: Since joint chance constraints are difficult to handle in general, one hopes that for a suitable choice of L i and α i the GSM is (close to) a deterministic equivalent of such a chance-constraint stochastic program.
In the presence of a finite lead time/demand rate distribution we obtain: A solution to a GSM with lead times L i and demand rates α i satisfies all demands on time in scenario ξ ω without operational flexibility if the GSM solution is also feasible in a GSM with lead times L ω i and demand rates α ω i . This is the case if L ω i ≤ L i and α ω i ≤ α i . Thus, in order to achieve a target service level of n target ∈ (0, 1], the GSM must choose minimal L i and α i such that we have Even for finite general lead time/demand rate distributions L i and α i may not be uniquely determined by the distribution and the target service level. The GSM may even produce distinct results for distinct choices. Example 1 Consider a single-demand-node network with h = 1 ands out = 0. Let the scenario set be = (L 1 = 1, α 1 = 2), (L 2 = 3, α 2 = 1) with probability 1 2 for both. In order to achieve a target service level of 0.5 tightly, we can either use L = 1, α = 2 or L = 3 and α = 1. The former leads to optimal GSM costs of 2 via x = 1 and y = 2, while the latter achieves optimal GSM costs of 3 via x = 3 and y = 3.
We can, however, find unique lead time and demand rate bounds if we assume the following total-order property of the lead time/demand rate distribution: In this case, we can identify a critical scenario ω * , defined by It has the following property: A feasible solution to the GSM with lead times L ω * i and demand rates α ω * i is feasible for all scenarios ξ ω with ω ≤ ω * . Thus, with these choices any feasible solution to the GSM achieves at least the desired target service level. In an optimal GSM solution, for all i ∈ N (G) constraints (4) and (5) of the GSM are binding. Consequently, the choices of L ω * i and α ω * i are unique. In this situation, we can derive induced lead times and demand rates from a finite lead time/demand rate distribution and a target service level: Definition 1 (Induced Lead-Times and Demand Rates) Let ( , p) be a finite lead time/demand rate distribution with the total-order property (42). Moreover, let n target be a positive target service level that determines a critical scenario ξ * = ξ ω * ∈ . Then the lead time and the demand rate in node i ∈ N (G) induced by ( , p) and n target are L * = L ω * i and α * = α ω * i , respectively. Given an SGSM with input data ( , p), we call the GSM with identical marginal holding costs and induced lead times and demand rates the GSM induced by the SGSM.
Remark 2 In the case of node-individual target service levels, the chance constraints are not joint over all nodes, but lead time and demand rate restrictions still need to hold jointly. (Even prescribing separate target service levels for those is conceivable, which results in purely individual chance constraints.) It is then sufficient that scenarios can be reordered at each node individually such that lead times and demand rates are monotonically increasing simultaneously. We then have to keep track of a critical scenario for each node. The interpretation is that at each node individually at least the target service level fraction of the demand is handled inside the model. The fraction of the total demand over all nodes handled inside the model can be much smaller then. Prescribing more and more target service levels means the need to properly guess more and more parameter values. The SGSM can be seen as an optimization model for this task. Note that a total order property (global or local) was only introduced for a consistent chance-constraint interpretation of the GSM. The SGSM does not need it.

Solution spaces of SGSMs and induced GSMs
Now that we have related the SGSM input data with the GSM input data via the target service level, we can perform a tighter comparison between the two.
In the following, let ( , p) be a finite lead time/demand distribution with the totalorder property. For any target service level n target we denote the critical scenario by ξ * = ξ ω * . Recall that ω * ∈ Ω is minimal such that ω≤ω * p ω = n * ≥ n target . Here, n * is the actual service level of the GSM with induced lead times and demand rates. Moreover, in node i ∈ N (G) let the critical lead time be L * i = L ω * i , and let the critical demand rate be α * i = α ω * i . The following result tells us, given a finite lead time/demand rate distribution, that for any target service level in the GSM, the recourse costs of the SGSM can always be adjusted in such a way that the decisions of the GSM are optimal first-stage decisions in the SGSM.
Theorem 1 Let ( , p) be a finite lead time/demand distribution with the total-order property and positive demands. Let n target be a target service level. Moreover, let (s in ) GSM , (s in ) GSM , x GSM , y GSM be optimal for the GSM with induced lead times L * i ≥ 0 and demand rates α * i > 0. Then there are marginal expediting costs t i and marginal outsourcing costs c i such that for the corresponding SGSM the following solution induced by the GSM is optimal: The proof is an extension of the proof for the equivalence of independent chanceconstraint programs and two-stage stochastic programs with simple recourse (see, e.g., Birge and Louveaux 1997, Section 3.2) via optimality conditions. The classic result cannot be applied directly because our technology matrix depends on the stochastic demand rates. The details of the rather technical proof are given in "Appendix 1".
Next, we show that there are SGSM solutions that can not be found by solving an induced GSM.
Theorem 2 There is a single-node warehouse network G, holding costs h, a lead time/demand-distribution with total-order property, and marginal expediting and outsourcing costs t and c, respectively, with the following property: There is no target service level n target such that the induced GSM has an optimal solution that is optimal for the first stage of the SGSM.
Proof The network consists of a single node with holding cost h = 1. Consider the following data: Depending on an actual service level of 1 3 , 2 3 , or 1, the corresponding induced GSM has the following unique optimal solutions with the indicated GSM costs: The corresponding SGSM solutions with fixed first-stage solutions, the optimal respective recourse values, and their costs read as follows: (s in = 0, s in = 0, x = 3, y = 9, r 1 = 0, q 1 = 0, r 2 = 0, q 1 = 0, r 3 = 0, q 3 = 0) → 9.
The previous example shows that the first-stage part of an optimal SGSM solution need not be extremal in the space of first-stage variables restricted by one scenario only. Since in the example the second scenario equals the average scenario, using the average scenario does not help the GSM to find the optimal first stage either.
In particular, the set of policies that can be generated by SGSM instances is a strict superset of the set of policies that can be generated by induced GSM instances.
In Sect 4.1 we review some basics about Sample-Average-Approximation (SAA) Methods for general finite approximations of probability distributions. The extensive form of the deterministic equivalent problem grows with an increasing number of scenarios. Therefore, we employ scenario reduction as described in Sect. 4.2.

Scenario generation: the SAA-method
To approximate the distributions of the stochastic parameters we generate random numbers according to the assumed distribution. These random numbers build the scenarios in the discrete distribution approximating the real distribution of the stochastic parameters. All samples are assigned probabilities proportional to the number of times they were generated. Sampling techniques like this are quite common in stochastic programming. See for example Birge and Louveaux (1997).
The idea of SAA is to estimate the optimal value function Q(x) of the second stage of a two-stage stochastic program min x∈X c T x + Q(x) with first-stage variables x, second stage variables y, and random parameters ξ . This function is originally defined as Consider a set Ξ of independent, identically distributed samples ξ ∈ Ξ of the random parameters. Then the following is an unbiased estimator for Q(x): and the problem min which can be solved conceptually much easier than the original problem, yields an unbiased estimator for the solution of the original problem. Further information about SAA can be found in Shapiro (2003) for example.

Scenario reduction: the fast forward selection
The goal of scenario reduction is to approximate a discrete distribution with many scenarios by another discrete distribution with significantly fewer scenarios. There are several methods to achieve this goal, usually based on a metric on the space of all possible scenarios (see Heitsch (2007); Heitsch and Römisch (2003); Henrion et al. (2008)).
Let us now sketch the principle of scenario reduction, since we have to make some choices. The approach to reduce the number of scenarios is based on a distance between two scenarios denoted by d(ξ 1 , ξ 2 ), a quantity that we have to define. Let Ξ be any subset of Ξ . Then the distance d(ξ, Ξ ) of a scenario ξ ∈ Ξ to the subset Ξ is given by d(ξ, Ξ ) := min ξ ∈Ξ d(ξ, ξ ).
The discrete distribution with scenario set Ξ that approximates Ξ best with respect to the distance d can be obtained in the following way: Add the probability of all ξ ∈ Ξ \ Ξ to the probability of the scenario ξ (ξ ) ∈ Ξ that is the closest scenario in Ξ to ξ with respect to d. In a sense, the scenarios not in Ξ are replaced by duplicates of the closest scenarios in Ξ . The quality of this approximation can be quantified by the distance d(Ξ, ) between Ξ and Ξ , which is defined as d(Ξ, ) = ξ ∈Ξ p ξ d(ξ, Ξ ). The smaller the distance, the better the approximation. The goal of scenario reduction in general is to find a particular subset Ξ * of a given cardinality such that d(Ξ, * ) is minimal.
The exact scenario reduction problem given any distance function d on scenarios and a cardinality k can be formulated as the following k-median problem: There are several exact and heuristic approaches to this NP-hard combinatorial optimization problem. For scenario reduction, there are special methods available, like the fast forward selection, one of the heuristics introduced in Heitsch (2007); Heitsch and Römisch (2003); Henrion et al. (2008).
The fast forward heuristic works as follows (see Heitsch (2007); Heitsch and Römisch (2003); Henrion et al. (2008) for details). It uses the fact that, by enumeration in time quadratic in the number of scenarios, we can find the scenario that yields the best single-scenario approximation. The idea is now as follows. We sequentially extend the current approximation by one new scenario at a time such that in each step Ξ has minimal distance to the current scenario set union the new one.
Assume we already have constructed a subset Ξ m−1 consisting of m − 1 scenarios, 1 ≤ m ≤ k. Consider the updated distance function By construction, we have for all ξ ∈ Ξ Thus, the best possible one-element extension of a reduced scenario set with respect to d is equal to the best possible one-element scenario approximation with respect tod. Summarized, the fast forward selection works as follows: The approximation of the lead times and demand distributions is split into two parts. First, a finite number of samples ξ ∈ Ξ is generated according to the assumed distribution. These samples build a first discrete approximation where every scenario instance occurs with equal probability p ξ = 1/|Ξ |. Second, the resulting discrete distribution is fed into the fast-forward scenario reduction, i.e., it is approximated by a discrete distribution over a subset of scenarios of prescribed cardinality, which have, in general, non-uniform probabilities.

Symmetric and asymmetric distances of lead-time/demand scenarios
In our computational tests we use two different kinds of distances between two scenarios. The distances are defined by first defining component-wise distances for the lead times and demands of each node, separately. These component-wise distances can be -transformed into a distance between complete scenarios by computing the Euclidean norm of the components' distances or -used for a component-wise scenario reduction whenever the lead time/demand distribution is the product of the components' distributions (this is what we did in the simulations in Sect. 5).
The first distance we will refer to as the symmetric distance. For the lead time component we define Since the demand component of a scenario in a node consists of different demand rates α for every time interval, we have to compare piecewise linear functions. We assume an equidistant discretization of time (e.g., in weeks) and a piecewise linear demand (e.g., constant demand rates in each week). Let the time intervals of linearity be numbered with increasing time by r ∈ R (e.g., week numbers). We will use a discount parameter Δ that depends on the length of the linear pieces of the demand. The motivation is as follows: Since any forecasting error has consequences for the complete remaining time, its influence is smaller if the error occurs later. In our experiments we used Δ = 2 for a time discretization in months. The order is piecewise linear. We separately compare the demand rates α ω,r on the different domains of linearity and assign a weigth to each difference. The higher r , the further in the future the demand rate is realized, and differences for high values of r get a lower weight in the total distance calculation than the ones at the beginning of the forecast period. This way, scenarios that only differ in later time intervals are considered closer than scenarios that differ in earlier time intervals. We formally define the distance between demand Ψ 1 i and Ψ 2 i of two scenarios 1 and 2 as where α ω,r i denotes the demand rate during time interval with index r (e.g., during the r th week) in scenario ξ ω .
There is another option that leads to asymmetric distances. The idea is to anticipate that the approximation is constructed for the use in a stochastic optimization problem. Thus, we would like to find the approximation that yields the least change in the result of the optimization. To decide which scenario is more important for optimization, we need some information about the costs that occur in case of stockholding and in case of stock out. We have this information given as parameter h i , costs for holding one piece in stock, and c i costs for having a stockout of one piece.
This way, we can define the asymmetric distance between the lead time components of two scenarios as The definition of asymmetric distance between the demand components of two scenarios is based on the same idea. We define the following asymmetric distance between the demand components of scenarios: The distance between complete scenarios ξ 1 and ξ 2 , in both the symmetric and asymmetric case, can now be defined as Alternatively, if the components are stochastically independent, we can perform scenario reduction component-wise. This way, a complete scenario in the reduced scenario set is a combination of scenario components from the reduced sets of scenario components. The asymmetric reduction does not necessarily approximate the distribution itself as faithfully as the reduction technique based on symmetric distances. We get a bias in our approximation that depends on the fraction of h i and c i . It will be shown in the next section that this biased reduction is indeed a better approximation to the solution of the optimization problems because it takes into account the cost of a decision in a certain scenario. To the best of our knowledge, this is not a standard method in the Stochastic Programming literature. Although we have not yet any further evidence beyond the problem studied in this paper, we conjecture that objective-aware scenario reduction might be worth a try in other contexts, too.

Simulation experiments on real data
We performed comprehensive computational tests on real-world data from our partner. All computational results report costs incurred by a method in a discrete-time simulation.

Discretization of time
The SGSM can only take finite discrete distributions of lead times and demands. Moreover, all scenarios of the demand distributions must be represented by piecewise linear approximations in order to obtain an MILP formulation for the SGSM.
Our partner forecasts the demand for 1 month. The data includes the expected total demand in the current month, the expected total demand in the coming month and so on. Thus, a straight-forward approach would be to approximate the demand linearly during 1 month. However: If we simply assume linearity of the demand during 1 month, then the rough discretization of time into months leads to demand scenarios with too little variation over time.
We can, of course, choose a finer discretization of time in weeks or days. The finer the discretization is, the more realistic the demand function becomes. In order to get a feeling for this influence, we generated random numbers representing the demand over 1 month or 1 week. Figure 1 shows an example of differences in the scenarios for discretization in months and in weeks.
A problem arises if the discretization of time becomes too small. The shorter the linear pieces in the demand functions, the more variables and constraints in the resulting MILP. This is the reason why the results in Sect. 5.4 are all based on discretization in months or weeks. In our simulations, the discretization in days did not lead to substantial savings compared to the one in weeks. (Our simulation itself does not linearize the demand.) Besides the time discretization, the number of scenarios included in the model is the other quantity that is critical for the mere size and therefore to the computing time of the SGSM. Therefore, we checked the effectivity of the scenario reduction methods in our tests.

Setup of the simulation experiments
We implemented a discrete-time simulation system with independent uniformly distributed random lead times and Poisson distributed demands at all nodes in order to compare the approximated expected long-term dynamic costs incurred by various (s, S)-policies. The policies only differ in the method to compute s for each node.
We investigated two things: 1. How do the parameter settings for scenario generation and scenario reduction influence the long-term cost in the simulation? 2. How does the SGSM with good parameter settings compare with competing algorithms to optimize safety stock?
For the solution of the SGSM we used Sample Average Approximation (SAA) with subsequent scenario reduction according to the old and new techniques Sect. 4.1 with discount parameters Δ = 2 and Δ = 1.25 for a time discretization in months and weeks, respectively. In order to assess the best parameter settings of our scenario reduction, we computed the simulation results for various sizes of sampled and reduced scenario sets (see Sect. 4.2). The GSM was parameterized by the target service level: we investigated GSM with service levels that are frequently required at our partner's: 90 and 96 %. These are denoted by GSM(90 %) and GSM(96 %). The decentralized policies DEZ(90 %) and DEZ(96 %) for the service levels of 90 and 96 %, respectively, served as a benchmark for policies ignoring the network effects: In these models each location tries to reach the given service level target independently at the smallest cost.
All results reported in Sect. 5.4 refer to the average long-term cost returned from the respective ten simulation runs for a method under consideration. Since none of the methods is based on an exact model of the dynamic development of the system, each method produces individual systematic errors in the prediction of the long-term costs as soon as the environment does not satisfy all the respective assumptions. Therefore, we chose to make all cost comparisons in a common simulation environment rather than in one of the safety-stock models. We chose a simulation environment that matches the real situation at our partner's as closely as possible. This way we could assess best the real-world impact of the different model approximations and the different computational approximations made in the various methods to compute safety stocks.

Test data
The data used in the simulation is a real-world data set from our partner. The warehouse network is a star-shaped two-echelon spare parts distribution system. It has one master warehouse (no. 0) and seven warehouses (nos. 1-7) for end customer service. The model SGSM is not restricted to this special structure; it can be applied to any acyclic network structure by straightforward modifications.
For our tests we used the same benchmark data set (inventory costs and demand intensities for 1,127 spare parts) as in Rambau and Schade (2010). We used cost coefficients from cost estimates of our partner for inventory cost and the piece-based so-called "non-sales" cost, that determine the recourse cost coefficients. Since these cost coefficients depend on the part, there are too many of them to be listed here.
We used a simulation horizon of 25 months. The end-customer demands were Poisson distributed. The intensities were estimated from historical data. The lead times were uniformly distributed in an interval from 80 through 120 % around the expected value. The scenario sets were generated as products of independently sampled leadtime and demand scenarios over all nodes. Each policy was confronted with identical sets of lead time and demand samples. Inventory was always controlled by an (s, S)policy. The values for s were chosen by the models under consideration. Demands had to be fulfilled (backorders with higher priority) whenever possible. Depending on the experiment, unmet demand was backlogged or considered lost. The order quantity was taken from our partner data. Our simulation did not allow for faster delivery than the service times computed by the models.
Because of the limited complexity of the network topology, all instances could be solved in less than an hour for the benchmark assortment of 1,127 parts by the MILP solver gurobi 3.0 (set to an optimality gap of five percent). More accuracy made computations slower but did not gain much with respect to the long-term cost of the resulting policy.

Computational results
We first tested the SGSM with many different cardinalities of generated and reduced scenario sets between 1 and 1,000, the range in which computations times were viable. 2 Table 1 shows a small selection of the results that show the effectiveness of the reduction methods. The results presented in this table are the average long-term costs of ten simulation runs. The safety stocks were determined by the SGSM, where "n → k" indicates that at each node n lead time and n demand scenarios were generated and reduced to k lead time and k demand scenarios. The lead times and demands were identical in all the simulations. We can see an enormous reduction in the total costs by applying the reduction techniques introduced in Sect. 4. In the case of generating only three scenarios we observe a very high variability in the costs over the ten simulation runs. During ten simulations, the minimal total costs were 16,379, and the maximal total costs were 33,546. Applying the symmetric/asymmetric reduction technique the minimal total costs were 6,622/3,246 and the maximal total costs were 7,606/3,546, respectively. The costs occurring in the single simulation runs are listed in "Appendix 2".
These results show that applying scenario reduction leads to a much lower variability in the costs because also scenarios with small probability are taken into account. We can see that the results for the asymmetric reduction are quite close to those where all the fifty generated scenarios are included in the model. Table 2 includes the service levels in the different locations during the first of the ten simulation runs.
The service levels in Table 2 show the differences between the symmetric and the (new) asymmetric reduction technique. The asymmetric technique takes into account that for many parts the quotient h i /c i is greater for the leaf warehouses than for the master warehouse. Therefore, for the symmetric technique we get a higher service level at the master warehouse (no. 0), but lower service levels at the warehouses (nos. 1-7).
Simulating the situation modeled in the SGSM with simple recourse leads to the results listed in Table 3. This table includes the average costs over ten simulation runs of the different approaches. Here we calculated the order points s using all the different methods and ran the simulation ten times with different lead times and demands. For all different approaches the lead times and demands in the simulations were identical.
The results for the decentralized method are worse than the results when the order points are calculated by the GSM. Using one of the listed SGSM approaches leads to a cost reduction of 30 % and more. Again, the asymmetric scenario reduction dominates the symmetric one. Another important aspect to notice is that the results using a discretization of time in weeks are remarkably better than results using a discretization in month. Results for each of the ten simulation runs for Model (4) and (10) can be found in "Appendix 2".
In Method (8) a special heuristic is applied (different from the fast forward reduction) that tries to find a critical scenario for the lead time and the demand for every location. This shows that much of the problem's structure can be encoded into a single scenario. This heuristic works properly for the discretization in months and may be extended to finer discretization. This is work in progress. The resulting service levels for the different methods in the first simulations are shown in Table 4.
The differences in the service levels of the symmetric and the asymmetric reduction are no longer substantial. The reason is that now the number of scenarios in the set Ξ * is much higher; thus, both approaches lead to a good approximation of the distribution and its impact on resulting service levels.
As Table 3 shows, the differences in the resulting costs are still remarkable. This is due to more scenarios in the more relevant parts of the distribution in the asymmetric reduction (high lead times and demands if h i /c i is low and vice versa).
The results of simulations of the SGSM with external suppliers from which missing parts can be ordered and lost sales (introduced in Sect. 2.3) are listed in Table 5: The simulation works different to the one applied in Tables 1, 2, 3 and 4. Here the demand that cannot be delivered immediately from the warehouses (nos. 1-7) to the end customers is lost. If the warehouses have not enough stock to deliver the ordered parts, there is the possibility to buy these parts from an external supplier. This recourse action causes costs depending on the distance between the warehouse and the external supplier. The supplier himself has limited stock so that the warehouses are not able to order any amount from them. If a demand at a warehouse can be neither delivered from stock nor ordered from an external supplier, the demand is lost. Internal orders (from a warehouse to the master warehouse) are still backlogged, and the master warehouse delivers the demand as soon as possible to the ordering warehouse. The ordering costs and the capacities of the external suppliers are not included in the data of our partner, so we had to set them artificially.
As we can see in the results of Table 5, the decentralized model and the GSM perform much better in the case with only one kind of uncertainty (demand uncertainty) than in the case of both lead time and demand uncertainty. The SGSM still outperforms the deterministic models achieving 10-20 % of cost savings. Table 6 show the resulting service levels of the different methods.
Here the service levels of the SGSM approaches are very similar to those of the decentralized model and the GSM, both with a prescribed service level of 96 %. The costs in Table 5 tell us that the SGSM treats different parts differently, while the GSM and the decentralized model cover 96 % of the demand for every part, no matter what the costs h i and c ji are. This is the reason why the order points calculated by the SGSM can lead to lower inventory costs and recourse costs at the same time.
Last we want to compare our model to a model that was introduced by Doǧru, de Kok and van Houtum, see Doǧru et al. (2005). In the following we will refer to this model as DoKoHo. In the simulation we need to apply fix lead times as this is one assumption of the DoKoHo model. We simulate a situation that fits to the DoKoHo assumptions where demand is backlogged and there are penalty costs if a location is not able to deliver as demanded. Table 7 shows the results of the simulation for the GSM, the SGSM, and DoKoHo. There are some different parameter settings for DoKoHo where the penalty costs used in the model are multiplied by a factor (γ ). The DoKoHo model outperforms the GSM but causes higher costs than the SGSM. The simulation of the different models lead to the service levels that are shown in Table 8.
The reason for the very low service levels at the master warehouse using the DoKoHo model compared to the ones using GSM or SGSM can be explained easily. In the DoKoHo model there are no explicit service times guaranteed to the successors. But the lower performance of the master warehouse is considered when the successor's safety stock is calculated. In the simulation we use a service time of zero for this case thus the delivery of master warehouse is often late. As all penalty costs at the master warehouse are zero in the simulation, this does not affect the costs that we get applying the DoKoHo models.

Conclusion
We have provided additional evidence that the Stochastic Guaranteed Service Model (SGSM), a stochastic programming version of the Guaranteed Service Model (GSM), first introduced in Rambau and Schade (2010), induces policies that outperform other policies for the computation of safety stock levels in a multi-echelon spare part distribution system of a large German car manufacturer. Moreover, we enhanced the simple-recourse SGSM by a recourse model that covers outsourcing in a more appropriate fashion by a transportation problem. In order to take advantage of the SGSM, we could show in our simulations that the stochasticity needs to be captured by sufficiently large sample sizes: in our example we generated 200 scenarios most of the time and reduced them to 50 applying modified scenario reduction techniques. The resulting MILP models could be solved straightforwardly in our example.
The SGSM makes some assumptions that are only approximations of reality (complete recourse, piecewise linear demand, exogenous lead time and demand distributions in internal nodes). However, our simulation was not restricted by these assumptions; it only checked the resulting policies, no matter what they assumed, and accounted for all the occurring costs. And in this quite realistic simulation experiment, the policies calculated with the SGSM performed extremely well. One reason for this is that the SGSM can have structurally different optimal solutions than the GSM: not all optimal SGSM solutions are extreme in the space of variables of the GSM. Thus the SGSM sometimes finds solutions that the GSM can never provide, no matter which target service level. And such solutions dominated the GSM solutions in our simulations.

Appendix 1: Proof of Theorem 1
We repeat the theorem here for convenience.
Theorem 1 Let (Ξ, p) be a finite lead time/demand distribution with the total-order property and positive demands. Let n target be a target service level. Moreover, let (s in ) GSM , (s in ) GSM , x GSM , y GSM be optimal for the GSM with induced lead times L * i ≥ 0 and demand rates α * i > 0. Then there are marginal expediting costs t i and marginal outsourcing costs c i such that for the corresponding SGSM the following solution induced by the GSM is optimal: Proof We prove the assertion by constructing marginal expediting and outsourcing costs from the complementary slackness condition of a primal-dual optimal solution to the GSM. From this, we construct a primal-dual solution to the corresponding SGSM that satisfies complementary slackness condition of the SGSM. We first list the GSM and its dual DGSM with the assumptions of this section (no integrality constraints, constant demand rates), where all variables appear on the left-hand side and all constants on the right-hand side. With this, the GSM reads as follows: With dual variables π i , ρ i j , σ i , and τ i corresponding in that order to the four sets of restrictions, the dual DGSM of the GSM, with restrictions ordered according to s in i , s out i , x i , y i reads as follows: −π i − j:i j∈A(G) From this we derive the optimality conditions via complementary slackness in primaldual pairs of feasible solutions: Next we do the same for the SGSM. The primal reads as follows: