PERIODIC REVIEW SYSTEM FOR INVENTORY REPLENISHMENT CONTROL FOR A TWO-ECHELON LOGISTICS NETWORK UNDER DEMAND UNCERTAINTY: A TWO-STAGE STOCHASTIC PROGRAMING APPROACH

Here, we propose a novel methodology for replenishment and control systems for inventories of two-echelon logistics networks using a two-stage stochastic programming, considering periodic review and uncertain demands. In addition, to achieve better customer services, we introduce a variable rationing rule to address quantities of the item in short. The devised models are reformulated into their deterministic equivalent, resulting in nonlinear mixed-integer programming models, which are then approximately linearized. To deal with the uncertain nature of the item demand levels, we apply a Monte Carlo simulationbased method to generate finite and discrete sets of scenarios. Moreover, the proposed approach does not require restricted assumptions to the behavior of the probabilistic phenomena, as does several existing methods in the literature. Numerical experiments with the proposed approach for randomly generated instances of the problem show results with errors around 1%.


INTRODUCTION
Inventory management pervades the decision-making 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 multi-echelon 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 single-echelon 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 single-echelon LN with stationary demand, the periodic review yields the best results, and, in the case of a multi-echelon 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 two-echelon 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 two-echelon distribution series system, a DC and a single retailer, with a pre-set 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 two-echelon distribution arborescent system, a DC and multiple retailers, and developed an approximate solution method for obtaining an optimal policy with pre-set 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 pre-set service levels.Chu & Shen (2010) provided an approximate solution for safety stocks policy with periodic review for all facilities in a two-echelon 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 pre-set 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 well-known 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 time-dependent 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 two-stage with resource.In this technique, the first-stage 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 here-and-now 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 wait-and-see 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 single-item 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 single-echelon LNs.Using stochastic programming, some studies already consider multi-echelon systems, as in the case of Gupta & Maranas (2000), Santoso et  In this paper, we present a new methodology based on two-stage stochastic programming with the use of SAA to support decision-making regarding the inventory management policy for a single item in an arborescent two-echelon 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: 1.Here we propose two-stage stochastic programming models to determine the optimal parameters (R, S) of the replenishment policy for single-item two-echelon LNs with a single DC and multiple retailers with uncertain demands.The deterministic equivalent models here proposed are mixed-integer nonlinear programming models, which are linearized in a novel technique.
2. 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 pre-set 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.
3.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.
4. 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.

PROBLEM DESCRIPTION
Considering a single-item two-echelon 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, . . ., N I , adopt respectively inventory control policies (R 0 , S 0 ) and (R i , S i ), where N I is the total of retailers in the distribution system, R 0 and R i represent the review interval at the DC and at each retailer i, and S 0 and S i 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 S i , and the optimal review intervals R 0 and R i , i = 1, . . ., N I , 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 uniform-sized time periods.
Let N P be the total of periods of a given planning horizon.The review intervals, or the times between orders, R 0 and R i , i = 1, . . ., N I , are modeled as multiples of the period p.The lead times L 0 and L i , i = 1, . . ., N I , 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 S i 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 + L i , i = 1, . . ., N I .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 to fulfill 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 two-echelon 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 (R i , S i ) at DC and at each retailer i = 1, . . ., N I , are: the carrying costs per item per period, h p 0 and h p i , and the ordering costs C p F 0 and C p F i per period, which is independent of the corresponding order quantities.Any demand that is not fully met by retailer i will be penalized by a shortage cost b p i per period, proportionally to the amount of items in short.There is no shortage cost for the case when the DC could not fully meet the amounts ordered by the retailers.The shortage cost is only applied to retailers that could not meet the customers' demands.
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(ξ ) p 0 and P(ξ ) p+2 0 , which are completely fulfilled at the beginning of periods p +1 and p +3, since the external supplier has no capacity fulfill limitations.In turn, the retailer places its orders with the DC at the beginning of periods p, p + 1 and p + 2, denoted by P(ξ ) , respectively.The order placed by the retailer with the DC in period p + 2 is partially fulfilled by the DC at the beginning of period p + 3, in the amount of that will be met in period p + 2.
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.

STOCHASTIC PROGRAMMING MODELS
In this section, we propose two-stage stochastic programming models to find the optimal parameters of the inventory policies (R 0 , S 0 ) and (R i , S i ), i = 1, . . ., N I , related to an arborescent single-item two-echelon 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 first-stage decisions consist of the optimal parameters (R 0 , S 0 ) and (R i , S i , f i ), i = 1, . . ., N I , where f i 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 second-stage decisions are related to the optimal inventory levels and orders quantities along the time horizon, which are directly affected by the first-stage 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.

Minimization of relevant costs with fixed rationing fractions
We propose a two-stage stochastic programming model to determine the optimal values R 0 , S 0 , R i , S i and f i , i = 1, . . ., N I , 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.

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: 10 y > 1, y ∈ N * , where 1/10 y represents the desired precision (e.g., if y = 1, the precision is decimal; if y = 2, it is centesimal; and so forth); w r 0 0, p auxiliary parameter that indicates the period that occurs an order at DC depending on the value r 0 ; w

Sets and indexes
auxiliary parameter that indicates the period that occurs an order at retailer i depending on the value r i ; w accumulated orders from all retailers fulfilled by the DC in scenario ξ in period p, where accumulated orders refer to the orders in period p plus all unmet orders of earlier periods; indicates if there exists a shortage or stock on hand at DC in scenario ξ at the end of period p; X (ξ ) p 0 ∈ {0, 1};

First-stage problem
The first-stage problem is related to the decisions with respect to the review intervals R 0 , R i , the target levels S 0 , S i , and the fractions f i , i = 1, . . ., N I , which must be made before the realization of uncertainty and aiming at minimizing ordering costs and the expected holding and shortage costs.The first-stage 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 second-stage problem.
Constraints ( 2) and ( 3) enforce that a single value for the review intervals R 0 and R i must be determined (R 0 = r 0 ∈ T 0 = {1, . . ., N R 0 } and R i = r i ∈ T i = {1, . . ., N R i }, when u r o 0 = 1 and u r i i = 1).Constraints (4) and (5) indicate that orders occur every interval R 0 at the DC and every interval R i at retailer i, and that the first orders occur always at the beginning of the first period of the time horizon (according to the parameters values w r 0 0, p and w r i i, p ). Constraints ( 6) and ( 7) define lower and upper bounds for the variables that represent the maximal inventory levels at the DC and retailers, respectively.Note that they can be set as storage capacities available at the CD and retailers.Finally, in (8) and ( 9), the first-stage variables u r k k , v p k , k = 0, . . ., N I , are defined as binary.

Second-stage problem
The second-stage problem consists of minimizing the carrying and shortage costs along the time horizon facing the choices for R 0 , R i , S 0 , S i and f i , i = 1, . . ., N I , 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 second-stage problem is formulated as: (10) subject to I (ξ ) Pesquisa Operacional, Vol.37(2), 2017 I e (ξ ) A(ξ A(ξ

P(ξ
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 f i 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 S i 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 v p i = 1); it is equal to zero, otherwise.Constraint (24) models the order quantity of the DC, at the beginning of period p, as the target level S 0 minus the echelon stock of DC 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 v p 0 = 1); it is equal to zero, otherwise.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 non-negativity for the decision variables, while constraint (28) set their initial value to zero.
Concerning F(ξ ) p 0,i in (21), we make two remarks.By definition of F(ξ ) p 0,i in (22), depending on the value of f i , F(ξ ) p 0,i can be lower than or equal to P(ξ is positive and there is no imbalance (allocation of negative shortage to retailers).Otherwise, if F(ξ ) p 0,i is greater than P(ξ p 0,i is negative and imbalance occurs.In this case, for the proposed model to be correct it is necessary to consider the relaxation of the constraint that imposes the non-negativity for variable A(ξ ) p 0,i .Moreover, in case A(ξ ) p 0,i is negative, this means that a shortage occurs and consequently a cost will be charged.If this consideration is not true, a solution with higher cost can be obtained in a scenario with lower probability of occurrence.
The deterministic equivalent model related to the two-stage 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 and I e I (ξ ) , constraints ( 23) and ( 24) are rewritten as 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)   ∀ p (43) 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 S i 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: I e (ξ ) Moreover, the order of retailer i not fulfilled by DC, denoted by F(ξ ) p 0,i , must be equal to the corresponding order P(ξ ) p i plus the quantity in short in the earlier period F(ξ ) p−1 0,i during L 0 :

Minimization of relevant costs with variable rationing fractions ( f (ξ )
p i ): Model B 3 − V 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 (ξ ) be revised in the second stage.Therefore, we rewrite constraint (22) as (55) to model the variable-fractions rationing rule: where and Nec(ξ ) p represents the sum of all the retailers' demands until period p and Nec(ξ ) p i represents the demand of retailer i until period p.

Notation
Consider also the following notation: fraction between the need of retailer i and the needs of all retailers until period p of scenario ξ ; Using the binary expansion and a predefined precision of 1/10 y , where y is a fixed known value, we observe that f (ξ ) p i is in the following interval: As the approximate value given by the binary expansion B V t bl(ξ ) p i,tb of the fraction f (ξ ) p i is always less than or equal to the corresponding true value plus the precision term, on the righthand side of the inequality (57), thus B V t bl(ξ ) p i,tb related to all retailers will often be less than 1.To overcome this drawback, the difference is rationed equally among all retailers.
From (56), we note that V tb l(ξ ) p i,tb considers the need of each retailer and thus a new binary representation is proposed: As Nec(ξ ) p l(ξ Thus, the approximate linearization of constraint (55) results in its substitution by the expressions (61) to (68):

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 pre-set value given by the decision maker.
The first-stage 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 second-stage problem model has the following additional constraints for a given scenario ξ : where fi is the expected value of the fractions of the demands of retailer i promptly fulfilled and F (ξ ) , and, as already defined, X (ξ ) p i indicates if there is or not shortage or stock on hand at retailer i at the end of period p.
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.

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 first-stage 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 first-stage objective function can be approximated by the following problem: where Q(R, S, f, ξ n ) is the objective function of the second-stage problem to be evaluated in each subset M and scenario ξ n .Given a collection of subsets of scenarios independently generated by sampling (ξ 1 j , . . ., ξ N j ), j = 1, . . ., M, we have: where the value of the first-stage 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 first-stage solutions (R , S , f ), the objective function value of the firststage model can be approximated by: Given ξ 1 j , . . ., ξ N j , j = 1, . . ., M , we have: 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 second-stage 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 φN is an unbiased estimator for φ(R , S , f ).Since (R , S , f ) is a feasible solution of the problem, φ(R , S , f ) is greater than or equal to the minimal value of the problem.And, since φN ,M is an unbiased estimator for φN , φN ,M is also greater than or equal to the optimal minimum value of the problem.Thus, φN ,M can be considered as an upper bound (UB) for the optimal minimum value of the original objective function.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 σ 2 L B and σ 2 U B are respectively the variance estimators of LB and UB.Besides the obtained bounds, the estimates of the optimality gap and its variance are: 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 Homemde-Melo, 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 Pesquisa Operacional, Vol.37(2), 2017 where β ∈ [0, 1].Besides the theoretical result stated in (81), in practice, the determination of the number of scenarios should consider the trade-off between the computational effort and the quality desired for the solution.

NUMERICAL EXPERIMENTS
In this section, we show two computational experiments conducted with the proposed stochastic programing models for randomly generated instances.The two-stage 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 Duo-Core 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 fixed-fractions 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.89and σN = 17.52, resulting in N > 5.28.
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 (N P = 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 , S i , f i and r 0 of the 10 solutions already obtained.Then, we performed 100 runs (M = 100), with 50 scenarios (N = 50) and 50 periods (N P = 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 , S i , f i , L B, U B, with precision of 0.1 (y = 1) for the binary expansion, where column lbe shows the percentage error for the estimate of LB ((σ L B /L B) * 100), and column ube shows the percentage error for the estimate of UB ((σ U B /U B) * 100).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 , S i , f i , L B, U B, 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 h p 0 = 1 and h p i = 4, ∀i, for I1, I3 and I4, and h p 0 = 3.5 for I2; and the expected value of the fraction for the demand fully met at the retailers are fi = 85%, 90%, 95% and 99%, ∀i, for I1 and I2, and for I3 and I4 we set f1 = 85%, f2 = 90% and f3 = 95%.In this experiment, the review intervals at DC and retailers were fixed respectively at R 0 = 3 and R i = 1, ∀i, with lead times L 0 = L i = 1, ∀i.The carrying and the shortage costs were set to zero in the first 3 periods to allow the initialization of stocks.
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 non-stationary 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 non-stationary.
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 (N P = 30).To get the upper bound, as in (78), a candidate solution is set with the mean values for S 0 , S i and f i of the 10 solutions obtained.Then, we performed 100 runs (M = 100), considering 30 scenarios (N = 30) and 30 periods (N P = 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.
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 f i change according to the pre-set values for the service levels, whereas the values obtained for the mean fractions of unmet demands f m i do not change for the majority of the values of the service levels.This is due to the fact that f i is a first-stage decision variable whose value is determined through the optimization process, while the value of f m i is imposed by the variable-fractions rationing       rule.Also, from Tables 4 and 5, we verify that S 0 values decrease and S i values increase as h p 0 increases for P 2 − F, as well as for P 2 − V .Table 6 shows that the results for P 2 − V model are associated with a smaller gap.From Table 7 we verify that both rationing rules are similar even considering a non-stationary stochastic process to generate the 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 − fi = 1%) I2 and I3, and between 2% and 3% in the other cases.

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 two-echelon 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 mixed-integer linear programming models that are deterministic equivalent formulations of the proposed two-stage 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 non-stationary 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 pre-set 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 two-stage 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 multi-echelon case, and the extension of computational assessments for additional instances and alternative configurations.10 y > 1, y ∈ N * , where 1/(10) y represents the desired precision (e.g., if y = 1, the precision is decimal; if y = 2, it is centesimal; and so forth); Pesquisa Operacional, Vol.37(2), 2017 P.S.A. CUNHA, F. OLIVEIRA and F.M.P. RAUPP 275 w r 0 0, p auxiliary parameter that indicates the period that occurs an order at DC depending on the value r 0 ; w r 0 0, p ∈ {0, 1}; r 0 = 1, . . ., N R 0 ; p = 1, . . ., N P ; w r i i, p auxiliary parameter that indicates the period that occurs an order at retailer i depending on the value r i ; w auxiliary variable representing the amount of unmet orders placed by retailers in scenario ξ in period p j i,tb auxiliary binary variable in approximating the binary representation of f i ; 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 two-echelon 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 single-item single-echelon LNs with uncertain demands, with periodic review and variable order quantities in regard to ordering, holding and shortage costs.

.
The first two orders are completely fulfilled by the DC at the beginning of periods p + 1 and p + 2, denoted by A(ξ ) p 0 and A(ξ ) p+1 0

Figure 1 -
Figure 1 -Scheme of the dynamics of a fixed control system.
), we introduce the linearized version of the problem model as follows.First, we introduce the variables I e I (ξ ) p 0 and I e I (ξ ) p i , ∀ p, i, that represent, for the DC and each retailer i, the stock position at the beginning of period p for given scenario ξ .As I e I Also, the costs of carrying per unit of the item in period p are h p 0 = 1, ∀ p, and h p i = 4, ∀i, p; and the shortage cost per unit of the item per period is b p i = 10, ∀i, p.

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 , S i , f i , L B, U B, 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 − fi ) indicates the maximum limit for the expect value of the fraction of unmet demand.The last 3 columns, 1 − f1 , 1 − f2 and 1 − f3 , show the service level achieved for each instance with the experiment, considering a candidate solution with mean values of S 0 and S i , and the mean fraction of unmet demand given by

[ 29 ]B
YOU F & GROSSMANN IE. 2008.Mixed-integer nonlinear programming models and algorithms for large-scale supply chain design with stochastic inventory management.Industrial & Engineering Chemistry Research, 47(20): 7802-7817.[30] ZIPKIN PH. 2000.Foundations of Inventory Management.McGraw-Hill.size of binary representations; tb ∈ B = {1, . . ., N B }; where N B 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 N B = 3); I retailers, i ∈ I = {1, . . ., N I }; P time periods, p ∈ P = {1, . . ., N P }; scenarios, ξ ∈ ; T 0 possible review intervals at DC, r 0 ∈ T 0 = {1, . . ., N R 0 }; T i possible review intervals at retailer i, r i ∈ T i = {1, . . ., N R i }; the item in short in retailer i in period p; retailer i in period p; D(ξ ) p i demand at retailer i in scenario ξ in period p; h p 0 cost of carrying one unit in period p at DC; h p i cost of carrying one unit in period p at retailer i; I T I auxiliary parameter to compute the quantities of the item ordered given in terms of the upper bound for the stock position; L 0 lead time at DC L i lead time at retailer i M total of independent subsets of N demand scenarios N total of demand scenarios at retailers along the time horizon S auxiliary parameter that gives an upper bound for the inventory target level V tb auxiliary parameter, V tb ∈ 2 ∈ {0, 1}; r i = 1, . . ., N R i ; p = 1, . . ., N P ; Variables A(ξ ) p 0 accumulated orders from all retailers fulfilled by the DC in scenario ξ in period p, where accumulated orders refer to the orders in period p plus all unmet orders of earlier periods;A(ξ )p i accumulated demands met by retailer i in scenario ξ in period p, where accumulated demand stands for the demand in period p plus all unmet demands of earlier periods;A(ξ ) p 0,i accumulated orders from retailers i fulfilled by DC in scenario ξ in period p; di f (ξ ) p auxiliary variable to compute F(ξ ) all retailers not fulfilled by DC in scenario ξ in period p; F(ξ ) p i accumulated demands not met by retailer i in scenario ξ in period p; F(ξ ) p 0,i accumulated orders of retailer i unmet by DC in scenario ξ in period p; f i shortage fraction for retailer i, that is, orders placed by retailer i not fulfilled by the DC ( f i = B j i,tb V tb ); f (ξ ) p i fraction between the need of retailer i and the needs of all retailers until period p of scenario ξ ; I (ξ ) p 0 stock on hand of the DC in scenario ξ at the end of period p; I (ξ ) p i stock on hand of retailer i in scenario ξ at the end of period p; I e (ξ ) p 0 echelon stock of the DC in scenario ξ at the end of period p; I e (ξ ) p i echelon stock of retailer i in scenario ξ at the end of period p; the DC in scenario ξ at the beginning of period p; I e I (ξ ) p iechelon stock of retailer i in scenario ξ at the beginning of period p; the echelon stock of the DC in scenario ξ at the beginning of period p; the echelon stock of retailer i in scenario ξ at the beginning of period p; Sauxiliary parameter to compute the quantities of the item ordered given in terms of the upper bound for the inventory target level; sizes of binary representations; tb ∈ B = {1, . . ., N B }; where N B 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 N B = 3); I retailers, i ∈ I = {1, . . ., N I }; P time periods, p ∈ P = {1, . . ., N P }; scenarios, ξ ∈ ; T 0 possible review intervals at DC, r 0 ∈ T 0 = {1, . . ., N R 0 }; T i possible review intervals at retailer i, r i ∈ T i = {1, . . ., N R i }; I T I auxiliary parameter to compute the quantities of the item ordered given in terms of the upper bound for the stock position;

Table 1 -
B 3 − F equivalent deterministic model data.

Table 2 -
Numerical results with B 3 − F model instance with precision 0.1.

Table 3 -
P 2 − F and P 2 − V equivalent deterministic data model.

Table 4 -
Comparative numerical results for I1 instance.

Table 5 -
Comparative numerical results for I2 instance.

Table 6 -
Comparative numerical results for I3 instance.