1 INTRODUCTION
Inventory management pervades the decisionmaking in many logistics networks (LNs), being thus a topic of great interest in academia. The key questions that inventory management aims to answer are how often the inventory position should be verified, when to place an order, how large an order should be, and what is the amount to keep as safety stocks in case of probabilistic demands (^{Namit & Chen, 1999}). Maintaining safety stocks in multiechelon LNs still lead to other basic questions, such as, for example, what should be their amount in the entire system and how much stock should be allocated at different levels of the echelons (^{Axsäter, 2006}).
A replenishment control policy establishes rules and courses of action to answer key questions relating to inventory management. In particular, it means that safety stocks can be managed in distinct ways. For example, the decision on replenishment at each facility of a LN can be based directly on stock positions or on the echelon stock level of each facility (the sum of stock positions of the facility and of all the downstream facilities in the LN). Regardless of the policy adopted, the aim is to determine the best level of inventory investment to achieve the desired service level, i.e., to guarantee a minimal percentage of demand fulfillment.
In the literature, there are several inventory policy proposals for singleechelon LNs with probabilistic demands. Among them, one can mention the continuous review systems (s, Q) and (s, S), and the periodic review systems (R, S) and (R, s, S), where s is the order point, Q is the fixed order quantity, R represents the time interval between orders, and S is the maximum stock position or target level (^{Hadley & Whitin, 1963}; ^{Silver et al., 1998}; ^{Zipkin, 2000}).
Inventory control systems with periodic review are widely used both in retail and in manufacturing, as they require less transactional effort; involve easier planning for calculating workload needs; facilitate customer services and receiving from suppliers; allow better replenishment coordination, especially when they involve multiple items; and generate more stability for LNs. Furthermore, when dealing with a singleechelon LN with stationary demand, the periodic review yields the best results, and, in the case of a multiechelon network, it has the advantage of being easily implementable. In the latter case, although this system is not necessarily optimal, it is capable to provide nearly optimal solutions (^{Federgruen & Zipkin, 1984}).
In general, inventory management models available in the literature deal with twoechelon LNs with centralized information, which means that all the information on demand and inventory levels at retailers is shared with the distribution center (DC). (^{Clark & Scarf 1960}) presented the first study in this area. They proposed a model for a twoechelon distribution series system, a DC and a single retailer, with a preset review interval, and developed a solution method to obtain the optimal replenishment control policy. Using the same concept, (^{Axsäter 2006}) described a model considering a twoechelon distribution arborescent system, a DC and multiple retailers, and developed an approximate solution method for obtaining an optimal policy with preset review intervals and equal target levels for all retailers. In those models, the costs of carrying stocks and shortages of the item are considered, excluding the ordering cost, as the review intervals at the DC and at retailers are equal and known. (^{Van der Heijden et al. 1997}) proposed a central study in this area, addressing inventory allocation policies in a distribution system with n echelons, where it is possible to keep stocks at all levels of the structure to meet preset service levels. (^{Chu & Shen 2010}) provided an approximate solution for safety stocks policy with periodic review for all facilities in a twoechelon distribution system. In their model, the ratio between the review intervals at DC and at retailers is restricted to a power of two, and it is necessary to preset the service level in all facilities, including the DC.
In practice, besides the need for an inventory control policy, LNs with more than one echelon require the definition of a rationing rule for quantities of items in shortage. When a DC does not have sufficient stock of an item to completely and simultaneously meet all orders from retailers in a period, the rationing rule defines the shortage distribution that the DC perform to all retailers in that period.
Fair Share (FS) is the most wellknown rationing rule. According to (^{Jonsson et al. 1987}), its central idea is to minimize the quantity of the item on backlog of orders (when postponements of demand fulfillment in later periods are possible) by imposing equal shortage probabilities to all retailers. To overcome this limitation, (^{De Kok 1990}) proposed the Consistent Appropriate Share (CAS), in which rationing fractions are effectively fixed based on the demands during the replenishment time at retailers, generalizing FS. However, CAS may cause imbalances or allocations of negative shortages, which is the case when the allocated volume of a retail shortage is greater than its order placed at the DC. Later, (^{Van der Heijden 1997}) gave an important contribution to the development of rationing rules, determining the rationing fractions by minimizing an imbalance average measure through the introduction of the Balanced Stock (BS) rule.
(^{Lagodimos et al. 2008}) emphasize that existing inventory modeling assumptions in the literature vary according to the rationing rule considered. Studies related to FS generally assume that the demand follows a normal distribution, while the CAS and BS rules assume either Erlang or Gamma distribution, directly affecting the developed model. While studies related to FS usually have detailed analytical models, studies related to CAS and BS are more general, requiring both numerical integration as a special approach technique, as well as optimization techniques.
Moreover, to model inventory management in a more realistic viewpoint, one can use stochastic programming as an alternative. In practice, LNs may have stochastic phenomena that depend on timedependent market conditions, such as customer demands and prices of raw materials or freight. In fact, stochastic programming models are used when solutions of the corresponding problems are sensitive to changes in their uncertain parameters (^{Birge & Louveaux, 1997}), and especially when the premises with respect to the probabilistic phenomena are restrictive.
According to (^{Higle 2005}), the most applied stochastic programming model is the twostage with resource. In this technique, the firststage variables, often referred to as design variables, correspond to those decisions that must be made before the actual realization of the uncertain parameters is observed, also known as hereandnow decisions. Then, based on these decisions and realization of random events, resource variables are considered in the second stage, which in turn are linked to control decisions, also known as waitandsee decisions.
In these models, uncertainty is represented by a discrete and finite set of scenarios that approximates the original (continuous or discrete) stochastic phenomena. To represent uncertain demand levels of a singleitem LN with a single echelon, (^{Cunha et al. 2017}) used the approach known as Sample Average Approximation (SAA) based on the Monte Carlo simulation method, which allows stochastic phenomena to be considered in a more adherent manner, without relying on any specific scenario generation method. The size of the set of scenarios is closely linked to the quality of the representation of the stochastic phenomenon, but it is important to note that the higher the cardinality of this set, the more challenging is the problem in terms of computational resources. Therefore, it is important to use appropriate techniques that allow obtaining good solutions in computational times that are acceptable in practice.
Research on LN design integrated to inventory management with uncertainty in one or more parameters is relatively new. Most of the literature focuses on singleechelon LNs. Using stochastic programming, some studies already consider multiechelon systems, as in the case of (^{Gupta & Maranas 2000}), (^{Santoso et al. 2005}), (^{Oliveira & Hamacher 2012}) and (^{Oliveira et al. 2013}). Nevertheless, despite considering inventory management and LN design jointly, these works did not address directly replenishment and inventory control policies. On the other hand, (^{Daskin et al. 2002}), (^{Shen et al. 2003}) and (^{You & Grossmann 2008}) addressed LN design and inventory policy optimization without the use of stochastic programming technique. Using stochastic programming, (^{Fattahi et al. 2015}) proposed a replenishment and inventory control methodology for a twoechelon LN in series, based on a continuous review policy (s, S), considering a single item with uncertain demand, and (^{Cunha et al. 2017}) proposed a replenishment policy for singleitem singleechelon LNs with uncertain demands, with periodic review and variable order quantities in regard to ordering, holding and shortage costs.
In this paper, we present a new methodology based on twostage stochastic programming with the use of SAA to support decisionmaking regarding the inventory management policy for a single item in an arborescent twoechelon LN, whose levels of demand are uncertain, over a finite time horizon. To this end, the proposed models seek to define the optimal parameters of a (R, S) control system, with periodic review and variable order quantity, to achieve minimum costs, introducing rationing fractions for quantities of the item in shortage, as well. Moreover, we performed computational experiments with generated instances to illustrate the potential of the proposed methodology.
The main contributions of this paper are detailed as follows:
Here we propose twostage stochastic programming models to determine the optimal parameters (R, S) of the replenishment policy for singleitem twoechelon LNs with a single DC and multiple retailers with uncertain demands. The deterministic equivalent models here proposed are mixedinteger nonlinear programming models, which are linearized in a novel technique.
Regarding the cited works of the literature that deals with the problem here addressed, we point out that in the proposed models: (a) the service levels are not considered as preset parameters; (b) the periodic review intervals can differ among the retailers and CD, since they are considered variables in the model; and (c) hypothesis to the stochastic phenome can be relaxed.
Moreover, we introduce a variable rationing rule to determine the fractions for shortage quantities to meet the retailers’ demands as possible, which indirectly improve the service levels. One remarkable feature of the proposed rationing rule is that it is capable of addressing negative allocations of shortage, an issue often observed when applying currently available rules.
Computational experiments with the proposed methodology are presented for randomly generated instances whose results are compared with the results obtained with a simulation technique.
It is worth mentioning that, in the proposed models, the relevant costs refers to ordering, carrying and shortage, which are considered deterministic, but possibly changing along the planning horizon, and the demand fulfilment being postponed if necessary (referred hereinafter as backlog case). In addition, the proposed approach does not require restrictive assumptions, such as time independence, normal distribution, nor fixed costs throughout the time horizon, to determine the optimal parameters (R, S) of the inventory policy, as well as any other assumption concerning the determination of the rationing fractions. It can thus be applied to a wider range of problems arising in this context.
In what follows, Section 2 shows the description of the problem under study. Section 3 introduces the proposed stochastic programming models to determine the optimal parameters of the control system (R, S). Section 4 briefly presents the technique used for discrete representation of random phenomenon. Section 5 presents the results of numerical experiments with instances randomly generated. Conclusion and future developments are presented in Section 6.
2 PROBLEM DESCRIPTION
Considering a singleitem twoechelon LN, the problem is to define an inventory control policy based on periodic review and variable order quantity under demand uncertainty, aimed at minimizing relevant costs and obtaining a desired service level. A DC and a set of retailers compose the arborescent distribution system. The DC places its orders to an external supplier, stores the item, and fulfills the retailers. Each retailer places its orders with the DC, stores the item and serves customers. The problem at hand does not consider costs, delays or capacities relative to transportation between the external supplier and the DC, between the DC and retailers, nor between retailers and customers.
The DC and each retailer i, i = 1, ..., _{ NI } adopt respectively inventory control policies (R _{0}, S _{0}) and (_{ Ri } , _{ Si } ), where _{ NI } is the total of retailers in the distribution system, R _{0} and _{ Ri } represent the review interval at the DC and at each retailer i, and S _{0} and _{ Si } represent the inventory target level at the DC and at each retailer i. Notice that the echelon stock referring to a retailer is equal to the stock position, as there are no facilities downstream. Thus, the problem consists in determining the optimal target levels S _{0} and _{ Si } , and the optimal review intervals R _{0} and _{ Ri } , i = 1, ..., _{ NI } , of the inventory control policies of the DC and the retailers related to a single item, whose demands follow known probability density functions, along a time horizon composed of a finite and discrete number of uniformsized time periods.
Let _{ NP } be the total of periods of a given planning horizon. The review intervals, or the times between orders, R _{0} and _{ Ri } , i = 1, ..., _{ NI } , are modeled as multiples of the period p. The lead times L _{0} and _{ Li } , i = 1, ..., _{ NI } , are fixed and known a priori, and are also multiples of p. The amount ordered by the DC to the external supplier is given by the difference between the target level S _{0} and the stock position at the moment of the ordering. Likewise, each order placed by retailer i with the DC has a required amount given by the difference between the target level _{ Si } and the stock position at the moment of the ordering. In this control system, the first order of the DC, as well as the first order of each retailer, is placed respectively at the beginning of the first period of the time horizon, which will be delivered respectively at the beginning of periods 1 + L _{0} and 1 + _{ Li } , i = 1, ..., _{ NI } . The quantity of the item received at the beginning of a period can be consumed in the same period. Moreover, the external supplier has always sufficient stock tofulfill the DC. It is possible to store the item at all facilities, without inventory capacity limitations. The orders placed by the retailers with the DC, as well as the demands from the customers for the item, can be partially fulfilled, and the quantities in short are assumed to be fulfilled as soon as possible.
In the twoechelon arborescent distribution system being considered, the quantity of the item at the DC in the upstream echelon can be not sufficient to completely and simultaneously meet the demand of all retailers in the downstream echelon in a period. In this case, the decision maker must define in advance a rationing rule to be applied for distributing the shortage among the DC and retailers. To overcome this difficult, two distinct rules are considered alternatively. In the first rule, the rationing results from the application of a fixed and equal percentage to the quantities of the item in short throughout the time horizon. The second rule applies dynamically variable percentages, depending on the need of each retailer.
The relevant costs determining the optimal control parameters (R
_{0,}
S
_{0}) and (_{
Ri
}
_{,}
_{
Si
} ) at DC and at each retailer i = 1, ..., _{
NI
} , are: the carrying costs per item per period,
Figure 1 schematically illustrates a periodic review inventory control policy of a single item in a network with a DC and a single retailer. The lead times and the review intervals are set to L
_{0 =}
L
_{1} = 1 period, R
_{0} = 2 periods and R
_{1} = 1 period. In this control system, the DC places its orders with the external supplier at the beginning of periods p and p + 2, denoted by P
It is worth noticing that no assumption is made on the random phenomenon driving the demand levels along the time horizon. In particular, all existing methods capable of solving the problem at hand require that the demand levels should be independent random variables and that the stochastic process should be stationary. We point out, however, that these assumptions do not need to be enforced for the applicability of the methodology proposed in this paper.
3 STOCHASTIC PROGRAMMING MODELS
In this section, we propose twostage stochastic programming models to find the optimal parameters of the inventory policies (R _{0,} S _{0}) and (_{ Ri } _{,} _{ Si } ), i = 1, ..., _{ NI } , related to an arborescent singleitem twoechelon LN with uncertain demands, whose deterministic equivalent models are firstly formulated as mixed integer nonlinear programming (MINLP) models and then approximately linearized, resulting in mixed integer linear programming (MILP) models. Due to the consideration of rationing fractions for shortages in the inventory policies, a new linearization scheme is introduced based on the binary representation of decimal numbers. In the proposed models, the firststage decisions consist of the optimal parameters (R _{0,} S _{0}) and (_{ Ri } _{,} _{ Si, fi } ), i = 1, ..., _{ NI } , where _{ fi } represents the fraction between the accumulated orders of retailer i that were not fulfilled by the DC and the accumulated orders of all retailers that were not fulfilled by the DC. The secondstage decisions are related to the optimal inventory levels and orders quantities along the time horizon, which are directly affected by the firststage decisions and by the realization of the demand uncertainty.
Therefore, for determining the optimal parameters in the context of uncertain demands, we consider two distinct approaches. The first takes into account the minimization of relevant costs (ordering cost, carrying cost and shortage cost) and it is denoted by B _{3}. In the second approach, which is appropriated for the cases in which it is difficult to quantify the shortage cost, a condition related to the demands directly fulfilled from the stocks is imposed, by defining a minimal service level by means of setting a minimal demand fraction that must be promptly fulfilled (fill rate). This approach, denoted as P _{2}, consider the condition as a constraint in the model that minimizes the ordering and carrying costs. We remark that the notation B _{3} and P _{2} are consistent with that used in the literature (see, for example, ^{Silver et al., 1998}).
As follows, in Subsections 3.1 and 3.2, we present in detail the models for the problem at hand with the B _{3} approach considering fixed and variable rationing rules. Then, in Subsection 3.3, the model with the P _{2} approach is presented with fixed and variable rationing rules, as well.
3.1 Minimization of relevant costs with fixed rationing fractions (f _{i} ) : Model B _{3}  F
We propose a twostage stochastic programming model to determine the optimal values R _{0,} S _{0}, _{ Ri } , _{ Si } and _{ fi } , i = 1, ..., _{ NI } , for the periodic review and variable order quantity inventory control policies, aiming at minimizing the relevant costs (ordering cost, carrying cost and shortage cost) and satisfy the demands while considering the balance of stocks along the time horizon. When it is not possible to completely fulfill the orders placed by the retailers in a period, it is possible to postpone the fulfillment to future periods (backlog) with the definition of a fixedfraction rationing rule. Later, we will compare the results of this police with fixed rationing percentages with the one with variable percentages.
The model uncertainty is relative to the demand levels of the customers for the single item along the planning horizon, which are represented as random variables that follow a known continuous probability distribution. To represent the uncertain demands as discrete and finite phenomena, that is, as a number of finite scenarios with discrete values, we use a sampling technique based on the Sample Average Approximation, which will be briefly described in the following section. Hereinafter, we assume that this uncertainty representation is made available.
3.1.1 Notation
The notation used from now on for parameters and variables intends to represent the operation of the inventory control system being modelled, and so some of it could not match the usual notation that is found in the literature.
In addition to the notation previously used, we define the following:
Sets and indexes
B sizes of binary representations; tb ∈ B = {1, ..., _{ NB } }; where _{ NB } is the total of digits in the binary expansion (e.g., to represent 7 into binary base we need 3 digits, 111, and in this case _{ NB } =3);
I retailers, i ∈ I = {1, …, _{ NI } };
P time periods, p ∈ P = {1, …, _{ NP } };
Ω scenarios, ξ ∈ Ω;
T
_{0} possible review intervals at DC, r
_{0} ∈ T
_{0} = {1, …,
_{
Ti
} possible review intervals at retailer i, _{
ri
} ∈ _{
Ti
} = {1, …,
Parameters
D
_{
Vtb
} auxiliary parameter, _{
Vtb
} ∈
W auxiliary _{
Np
} ×
Variables
A
A
A
F
F
F
_{ fi } shortage fraction for retailer i, that is, orders placed by retailer i not fulfilled by the DC (_{ fi } = Σ_{ B ji,tbVtb) } ;
I
I
_{
Ie
}
_{
Ie
}
JF
_{ Ji,tb } auxiliary binary variable in approximating the binary representation of _{ fi } ;
P
P
X
3.1.2 Firststage problem
The firststage problem is related to the decisions with respect to the review intervals R _{0,} _{ Ri } , the target levels S _{0,} _{ Si } , and the fractions _{ fi, i = } 1_{ , ...,NI } which must be made before the realization of uncertainty and aiming at minimizing ordering costs and the expected holding and shortage costs. The firststage problem is formulated as follows:
Expression (1) models the total costs to be minimized. The first two terms refer to the sum of the costs of ordering along the planning horizon at the DC and retailers, while the third term represents the expected value of the total cost relative to secondstage problem.
Constraints (2) and (3) enforce that a single value for the review intervals R
_{0} and _{
Ri
} must be determined (R
_{0} = r
_{0} ∈ T
_{0} = {1, …,
3.1.3 Secondstage problem
The secondstage problem consists of minimizing the carrying and shortage costs along the time horizon facing the choices for R _{0,} _{ Ri, } S _{0,} _{ Si } and _{ fi, i = } 1_{ , ..., NI } , for a given realization of uncertain demands in each scenario ξ, to satisfy as possible the customers’ demands in all time periods. For each scenario ξ ∈ Ω, the secondstage problem is formulated as:
In the objective function (10), the first two terms model the carrying costs at the DC and all retailers, which considers existing stocks on hand at the end of each period p, while the last term represents the costs related to unmet demands by all retailers, that is, the shortage costs along the time horizon. The sum of the carrying and shortage costs is minimized over all time periods.
Constraints (11), (12), (13) and (14) represent the balance of the stocks on hand between two successive periods, in each scenario ξ, at the DC and all retailers, respectively. Similarly, constraints (15) and (16) represent the echelon stock balance between two successive periods, in each scenario ξ, at DC and all retailers, respectively. The order quantity is defined based on the echelon stock, enforcing the (R, S) control system.
Constraints (17) and (18) impose the DC fulfillment to be the sum of the retailers’ fulfillment. Constraints (19), (20) and (21) model the demand fulfillment in each period for each scenario. Constraint (22) imposes that the accumulated orders of retailer i not fulfilled by the DC must be equal to the fraction _{ fi } of the accumulated orders of all retailers not fulfilled by the DC.
Constraint (23) models the order quantity of retailer i at the beginning of each period p as the target level _{
Si
} minus the echelon stock of retailer i at the beginning of period p (that coincides with the echelon stock at the end of period p  1), at the beginning of each cycle (indicated when
Constraints (25) and (26) indicate if the stock on hand at the DC is sufficient to fulfill the orders of the retailers. When X(ξ) = 1, it means that all retailers’ orders are fulfilled, while, when X(ξ) = 0, at least one retailer is not fulfilled by the DC. Constraint (27) imposes nonnegativity for the decision variables, while constraint (28) set their initial value to zero.
Concerning F
The deterministic equivalent model related to the twostage stochastic programming problem is given by (1)(9) and by  Ω  replications of (10)(28). Observe that constraints from (22) to (24) turn the model to be characterized as a MINLP, which is a class of problems known as being computational challenging. Regarding constraints (23) and (24), we introduce the linearized version of the problem model as follows. First, we introduce the variables
The exact linearization of constraints (23) and (24), and the approximate linearization of constraint (22) result in the substitution of expressions (23) and (24) by expressions from (29) to (44), and in the substitution of constraint (22) by the expressions from (45) to (48), apart the introduction of nonnegativity constraints for the auxiliary variables JF
minimize (10)
subject to (11)(21); (25)(30)
As we consider stock positions at the end of periods, constraints (15) and (16) are divided into constraints (50)(53), such that the echelon stock of retailer i in the first period is equal to the target level _{ Si } minus the average demand per period, and that the echelon stock of the DC in the first period is equal to the target level S _{0} minus the sum of the average demand of the retailers, as follows:
Moreover, the order of retailer i not fulfilled by DC, denoted by F
3.2 Minimization of relevant costs with variable rationing fractions (
An alternative for rationing the quantity in short among the retailers is to set a fraction based on the retailers’ needs (accumulated unmet orders) for each period and each scenario, f
where
and Nec(ξ)^{
P
} represents the sum of all the retailers’ demands until period p and Nec
3.2.1 Notation
Consider also the following notation:
Variables
dif (ξ)^{
P
} auxiliary variable to compute F
f
LF
LNec
I
Nec(ξ)^{ P } total of all retailers’ demands until period p;
Nec
Using the binary expansion and a predefined precision of 1/10^{
y
} , where y is a fixed known value, we observe that f
As the approximate value given by the binary expansion ∑_{
B
}
_{
Vtbl
}
is rationed equally among all retailers.
From (56), we note that _{
Vtbl
}
As
Thus, the approximate linearization of constraint (55) results in its substitution by the expressions (61) to (68):
3.3 Minimization of costs with level service condition: Model P _{2}
In this approach, the objective is to minimize the ordering and carrying costs along the time horizon with the additional condition on the minimum service level, that is, the requirement that the average demand fulfillment is greater than or equal to a preset value given by the decision maker.
The firststage problem is equal to the corresponding model of approach B _{3}, except for the shortage cost that is removed from the objective function. Thus, the secondstage problem model has the following additional constraints for a given scenario ξ:
where
For this approach, the fixed and variable rationing rules are considered, resulting into the models named as P _{2}  F and P _{2}  V, respectively.
4 SAMPLING TECHNIQUE
Probabilistic parameters that follow continuous distributions impose some difficult in solving stochastic optimization problems. In particular, in the addressed problem the difficult is associated to the evaluation of the firststage objective function (1) given in general terms by
where Ω is the finite set of scenarios of the demands at retailers along the time horizon, with ξ ∈ Ω, and, for simplification, the rationing fractions are considered fixed for all scenarios.
To overcome this difficult, we use a sampling technique based on the Monte Carlo simulation method, known as Sample Average Approximation (SAA), (^{Santoso et al., 2005}). Its main idea is to seek an approximation of the objective function value, considering the average of solutions for a problem instance composed by M subsets of N scenarios, which are successively and independently sampled.
Thus, for each subset M, the firststage objective function can be approximated by the following problem:
where Q(R, S, f, ξ^{
n
} ) is the objective function of the secondstage problem to be evaluated in each subset M and scenario ξ^{
n
} . Given a collection of subsets of scenarios independently generated by sampling
where the value of the firststage objective function is approximated by:
According to (^{Santoso et al. 2005}), the expected value _{ ĝN } is less than or equal to the minimal optimal value and, since _{ ĝN, M } is a biased estimator of the expected value _{ ĝN } , the expected value _{ ĝN, M } is also less than or equal to the minimal optimal value of the problem. Hence, with this technique we can consider the minimal optimal value obtained as a lower bound (LB) of the optimal value of the original objective function.
Choosing good feasible firststage solutions (R', S', f'), the objective function value of the firststage model can be approximated by:
Given
and thus,
where N' is the size of the sampling that does not depend on the sampling of size N used to compute (R', S', f'). Indeed, a total of M'N' deterministic and independent secondstage subproblems are being solved with lower computational effort. The advantage of using larger values for N' and M' refers to the attempt to reduce errors in the estimates.
^{Santoso et al. (2005}) also state that
(^{Linderoth et al. 2006}) show that with this technique we can obtain lower and upper bounds of the optimal value and that such bounds converge to the optimal value as N increases.
From the Central Limit Theorem, the confident interval for LB and UB, with level α, where P(z ≤ z _{α}) = 1  α, could be expressed respectively as
where
There are many ways to estimate the number of scenarios in order to obtain an approximation of the optimal value with a certain margin of error (^{Kleywegt et al., 2002}; ^{Shapiro and HomemdeMelo, 1998}). One of them is to use some of the ideas of the sampling technique, which gives statistical basis to get the number of scenarios.
To this end, we have from (74) that _{ ĝN } is the minimal expected value of the objective function, which is a random variable. Moreover, _{ ĝN } is also an estimator for the minimal value of the objective function. Thus, for each scenario of {ξ^{1}, ξ^{2}, ..., ξ^{ N } }, we have, for n = 1, …, N, that the expected values of the deterministic objective function values are:
with variance given by:
From the Central Limit Theorem, we define a confident interval with level α/2 for the estimator _{ ĝN }:
Now, using the reverse form of this interval, an estimate for the lower bound of the number of scenarios needed in the approximation of the objective function is
where β ∈ [0, 1]. Besides the theoretical result stated in (81), in practice, the determination of the number of scenarios should consider the tradeoff between the computational effort and the quality desired for the solution.
5 NUMERICAL EXPERIMENTS
In this section, we show two computational experiments conducted with the proposed stochastic programing models for randomly generated instances. The twostage stochastic programming models and the sampling technique were implemented using AIMMS version 3.13, and the corresponding MILP models were solved using CPLEX version 12.5. We performed the experiments in an AMD DuoCore 1.9 GHz processor with 4 GB RAM.
The first computational experiment is conducted with B
_{3}  F model, which minimizes the relevant costs, including the shortage cost and considering the fixedfractions rationing rule. To this purpose, we generated instances for the problem with 1 DC and 3 retailers, with parameters values set to:
In this experiment, the N value was defined according to (81). Therefore, for a 5% confidence interval, we set α =0.05 and β =0.1. Then, to get an approximation of the objective function value we considered N = 50σ_{50}, with _{
ĝN
} = 298.89 and
Assuming that a stationary stochastic process of second order represents the demands for the item along the periods, the N scenarios were generated as follows:
where a is the (constant) demand level and εp is the error associated to the model at each period, which follows a normal distribution with zero average and variance σ^{2}. For each scenario ξ and for all the periods of the time horizon, the demands for the item of the retailers follow a normal distribution with mean 27, 81 and 54, respectively, and variance 23, 39 and 31, respectively, having 54 as reference for the mean (54/2, (3/2)54 and 54) and 31 as reference for the variance (31  8, 31+8 and 31). The demands along the periods of each retailer follow identical probability distributions, without correlations among them.
To obtain the lower bound (LB), as in (75), for the approximation of the objective function we performed 10 runs (M = 10) considering 10 scenarios (N = 10) and 20 periods (_{ Np } = 20), which showed to be sufficient to observe in all runs the same result: r _{0} = 2. To obtain the upper bound (UB), as in (78), we considered a candidate solution given by the mean values for S _{0}, _{ Si , fi } and r _{0} of the 10 solutions already obtained. Then, we performed 100 runs (M = 100), with 50 scenarios (N = 50) and 50 periods (_{ Np } = 50) to estimate the upper bound. The carrying and the shortage costs were set to zero in the first 3 periods to allow the initialization of stocks.
Table 1 describes the experiment with the B _{3}  F model in terms of the instance size of its deterministic equivalent MILP model, and the computational time taken to solve it, given in terms of CPU time.
Table 2 shows the estimates for the optimal values for r _{0} , S _{0}, _{ Si ,fi, LB, UB } , with precision of 0.1 (y = 1) for the binary expansion, where column lbe shows the percentage error for the estimate of LB ((σ_{ LB } /LB)*100), and column ube shows the percentage error for the estimate of UB ((σ _{ UB } /UB)*100).
r _{ 0 }  S _{0}  _{ Si }  S _{2}  S _{3}  f _{1}  f _{2}  f _{3}  LB  σ_{ LB }  UB  σ_{ UB }  gap  lbe  ube 
2  661.7  59.7  168.6  114.0  0.26  0.38  0.36  299.4  5.6  301.7  1.6  2.3  1.9  0.5 
The obtained results suggest that the configuration of the experiment with 10 runs (M = 10) for LB is reasonable, since the percentage error is 1.9%, while for the UB the percentage error is 0.5%, resulting in 0.7% for the percentage error of the gap.
The second computational experiment aims to analyze the values S _{0}, _{ Si , fi, LB, UB } , and the expected value of the fraction relative to the unmet demand (shortage fraction) obtained by models P _{2}  F and P _{2}  V. To this end, we generated 4 instances, denoted by I1, I2, I3 and I4, where I4 differs from the others with respect to the nature of the stochastic process of the demands levels.
Considering a LN with 1 DC and 3 retailers, we have the following parameters setting: the costs of carrying per unit of the item in period p are
For all instances, N was set according to (81), obtaining _{ ĝN } with 50 scenarios. To generate N scenarios, for I1, I2 and I3 we supposed a stationary stochastic process of second order for the demands as (83), and for I4 we supposed a nonstationary stochastic process given by the random walk:
where εp is the corresponding error for each period, which follows a normal distribution with zero mean and variance σ^{2}. In this case, the presence of D(ξ)^{p  1}in (83) indicates dependence of the demand on previous period, which makes the process nonstationary.
For the instances I1, I2 and I3, the demands for the item in each retailer follow a normal distribution with mean 27, 81, and 54, respectively, and variance 23, 39 and 31, respectively, in each period. For I4, the initial demand in each retailer is set to 81, 54 and 67, and the corresponding error in each period follows a normal distribution with zero mean and variance equal to 1 for all retailers.
To obtain the lower bounds, as in (75), for the approximation of the objective function, we performed 10 runs (M = 10), with 10 scenarios (N = 10) and 30 periods _{ NP } = 30. To get the upper bound, as in (78), a candidate solution is set with the mean values for S _{0}, _{ Si } and _{ fi } of the 10 solutions obtained. Then, we performed 100 runs (M = 100), considering 30 scenarios (N = 30) and 30 periods (_{ NP } = 30 to get the upper bound.
Table 3 describes the experiment with P _{2}  F and P _{2}  V models in terms of the instance size of their equivalent deterministic MILP model, and the computational time in terms of CPU time.
Model  N  _{ NP }  Total of variables  Total of constraints  Time (s) 

P_{2}  F  10  30  11724 (1209 integer)  17721  248.45 
P_{2}  V  10  30  18605 (3900 integer)  28788  312.97 
Tables 4, 5, 6 and 7 show the comparative results of the models P
_{2}  F and P
_{2}  V with respect to S
_{0}, _{
Si, fi
} , LB, UB, with precision of 0.1 (y = 1) for the binary expansion, and to the expect value of the fractions of the unmet demand (shortage fractions) for instances I1, I2 I3 and I4. The first column (column 1 
for P _{2}  V, and S _{0}, _{ Si } and_{ , fi } for P _{2}  F, after 100 runs (M = 100), with 30 scenarios (N = 30) and 30 periods (_{ NP } = 30).
P _{ 2 }  F with precision 0.1  
1

S _{0}  S _{1}  S _{2}  S _{3}  f _{1}  f _{2}  f _{3}  LB  σ_{ LB }  UB  σ_{ UB }  gap  lbe  lbe  1

1

1

15%  754  56  162  108  0.17  0.51  0.32  153  2.2  151  1.7  1.4  1.4%  1.1%  15.6  15.3  15.5 
10%  778  57  163  109  0.17  0.52  0.31  174  1.5  173  2.3  1.5  0.8%  1.3%  10.0  10.1  10.3 
5%  805  59  166  112  0.18  0.49  0.33  208  3.0  208  3.0  0.  1.4%  1.4%  5.2  5.0  4.9 
1%  839  64  173  119  0.24  0.43  0.33  283  8.6  283  4.3  0.5  3.0%  1.5%  1.2  1.1  1.0 
P _{ 2 }  V with precision 0.1  
1

S _{0}  S _{1}  S _{2}  S _{3} 



LB  σ_{ LB }  UB  σ_{ UB }  gap  lbe  lbe  1

1

1

15%  754  55  162  108  0.14  0.53  0.33  152  1.9  151  1.6  1.09  1.3%  1.1%  15.8  15.3  15.4 
10%  778  56  163  109  0.14  0.53  0.33  174  2.9  172  2.1  2.38  1.6%  1.2%  10.9  10.2  10.2 
5%  806  58  167  112  0.14  0.53  0.33  209  3.0  210  2.7  0.70  1.4%  1.3%  5.3  4.6  4.8 
1%  840  64  173  118  0.15  0.53  0.33  286  7.6  285  2.9  1.0  2.7%  1.0%  1.0  1.1  1.1 
P _{ 2 }  F with precision 0.1  
1

S _{0}  S _{1}  S _{2}  S _{3}  f _{1}  f _{2}  f _{3}  LB  σ_{ LB }  UB  σ_{ UB }  gap  lbe  ube  1

1

1

15%  738  65  181  120  0.2  0.49  0.31  415  3.0  414  3.8  1.3  0.7%  0.9%  15.2  15.1  15.1 
10%  762  65  181  120  0.2  0.49  0.31  474  2.3  472  3.9  2.1  0.5%  0.8%  10.1  10.3  10.3 
5%  793  65  179  124  0.2  0.45  0.35  549  4.8  550  3.5  1.7  0.9%  0.6%  4.9  4.8  4.9 
1%  829  71  181  126  0.26  0.4  0.34  666  10.8  660  6.4  5.2  1.6%  1.0%  1.1  1.1  1.1 
P _{ 2 }  V with precision 0.1  
1

S _{0}  S _{1}  S _{2}  S _{3} 



LB  σ_{ LB }  UB  σ_{ UB }  gap  lbe  ube  1

1

1

15%  739  60  183  120  0.14  0.53  0.33  418  3.1  417  3.8  0.91  0.7%  0.9%  15.1  15.2  15.1 
10%  763  61  182  120  0.14  0.53  0.33  477  2.5  475  3.8  2.44  0.5%  0.8%  9.8  10.3  10.3 
5%  794  62  184  122  0.14  0.53  0.33  552  4.9  553  4.4  1.17  0.9%  0.8%  5.0  4.8  4.9 
1%  830  67  186  125  0.14  0.53  0.33  671  11.7  666  5.7  5.8  1.7%  0.9%  1.0  1.1  1.1 
P _{ 2 }  F with precision 0.1  
S _{ 0 }  S _{1}  S _{2}  S _{3}  f _{1}  f _{2}  f _{3}  LB  σ_{ LB }  UB  σ_{ UB }  gap  lbe  ube  1

1

1

785.7  55.2  164.3  112.0  0.22  0.58  0.19  179.9  2.8  183.2  2.2  3.3  1.6%  1.2%  13.7  9.1  4.7 
P _{ 2 }  V with precision 0.1  
S _{ 0 }  S _{1}  S _{2}  S _{3} 



LB  σ_{ LB }  UB  σ_{ UB }  gap  lbe  ube  1

1

1

784.9  53.3  161.2  114.6  0.14  0.53  0.33  180.7  2.8  181.8  2.2  1.2  1.5%  1.2%  15.2  10.4  4.8 
P _{ 2 }  F with precision 0.1  
S _{ 0 }  S _{1}  S _{2}  S _{3}  f _{1}  f _{2}  f _{3}  LB  σ_{ LB }  UB  σ_{ UB }  gap  lbe  ube  1

1

1

967.2  161.7  108.7  135.1  0.65  0.24  0.11  196.3  5.3  198.5  5.2  2.3  2.7%  2.6%  14.9  9.4  4.6 
P _{ 2 }  V with precision 0.1  
S _{ 0 }  S _{1}  S _{2}  S _{3} 



LB  σ_{ LB }  UB  σ_{ UB }  gap  lbe  ube  1

1

1

969.7  154.8  109.0  138.1  0.40  0.29  0.31  202.5  4.0  202.8  5.0  0.3  2.0%  2.5%  16.2  9.84  5.6 
Upon analyzing the numerical experiments regarding the instances I1, I2 I3 and I4, we observe that the service levels obtained with the minimization of the carrying costs are similar for the both rationing rules, when compared with the desired levels of service for the retailers. Also, we observe that the obtained results are associated with small percentage errors for all instances, even in the case of considering a small number of periods and scenarios. Still, from Tables 4 and 5, we observe that the values obtained for the fixed rationings fractions _{
fi
} change according to the preset values for the service levels, whereas the values obtained for the mean fractions of unmet demands
The results suggest in general that the configuration of the experiment, taking 10 runs for the estimation of the lower bound, is reasonable, since the gap and the standard deviation are small. We observe that the percentage errors lbe and ube are close to 1% for the instances I1 (except, when 1 
6 CONCLUSION
Based on stochastic programming, this research proposed a new comprehensive approach regarding the representation of demand uncertainty in the problem of determining the optimal parameters of an inventory control system of a single item with periodic review and assessment of shortage in an arborescent twoechelon logistics network over a finite time horizon. Indeed, this approach allows the relaxation of assumptions on the behavior of uncertain parameters, and in particular on the stochasticity behavior of the demand for the item. In relation to shortages of the item that may occur, we introduce a variable rationing rule at the DC that deals with imbalances and indirectly improves the service levels. Specifically, we proposed mixedinteger linear programming models that are deterministic equivalent formulations of the proposed twostage stochastic programming models to obtain optimal review interval and target level of the control system (R,S) of a logistics network with a DC and many retailers. To obtain approximate solutions for the deterministic models, we used a sampling technique to generate finite scenarios with discrete values for demand levels. Both stationary and nonstationary processes were considered to represent the uncertain demand.
Computational experiments were conducted with the proposed models for randomly generated instances considering different values for the carrying and shortage costs, and service levels. The results of the experiments suggest that the configuration taking 10 repetitions to estimate the lower bound was reasonable, since both gap and standard deviation are substantially reduced in all of the tests. Also, the results obtained for the service level targets are very similar to the ones preset for the retailers against the corresponding results obtained with the optimization of the costs for the fixed and the variable rationing shortage policies, even considering a small number of periods and scenarios.
This study demonstrated that it is possible to determine approximately the optimal inventory system parameters (R, S) through twostage stochastic programming approach. In most of the computational tests, the proposed approach gives optimal solutions with percentage errors close to 1% for the lower and upper bounds. This fact confirms that the proposed methodology has the advantage of being potentially applicable to a wide range of situations, since its application is independent of assumptions imposed with respect to stochastic phenomena.
It is worth to emphasize that this methodology has a limitation with respect to the large computational effort required, when the number of scenarios grows. One option would be to use decomposition methods to reduce this effort. In addition, we suggest for further studies to consider systems in other configurations, such as the extension to the multiechelon case, and the extension of computational assessments for additional instances and alternative configurations.