ABSTRACT

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 simulation-based 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%.

Keywords:
replenishment and control systems; two-echelon logistics networks; stochastic programming; shortage rationing rules

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}19 NAMIT K & CHEN J. 1999. Solutions to the ??Q, r?? inventory model for gamma lead-time demand. International Journal of Physical Distribution & Logistics Management, 29(2): 138-154.). 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}1 AXSÄTER S. 2006. Inventory Control. 2nd ed., Springer.).

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}12 HADLEY G & WHITIN TM. 1963. Analysis of Inventory Systems. Prentice-Hall.; ^{Silver et al., 1998}26 SILVER EA, PYKE DF & PETERSON R. 1998. Inventory Management and Production Planning and Scheduling. 3 ed. John Wiley & Sons.; ^{Zipkin, 2000}30 ZIPKIN PH. 2000. Foundations of Inventory Management. McGraw-Hill.).

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}10 FEDERGRUEN A & ZIPKIN PH. 1984. Computational issues in an infinite-horizon, Multiechelon Inventory Model. Operations Research, 32(4): 818-836.).

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}4 CLARK AJ & SCARF HE. 1960. Optimal policies for a multi-echelon inventory problem. Management Science, 6(4): 475-490.) 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}1 AXSÄTER S. 2006. Inventory Control. 2nd ed., Springer.) 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}28 VAN DER HEIJDEN MC, DIKS EB & DE KOK AG. 1997. Stock allocation in general multi-echelon distributionsystems with (R,S) order-up-to-policies. International Journal of Production Economics, 49(2): 157-174.) 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}3 CHU LY & SHEN ZM. 2010. A Power-of-Two ordering policy for one-warehouse multiretailer systems with stochastic demand. Operations Research, 58(2): 492-502.) 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}14 JONSSON H & SILVER EA. 1987. Analysis of a two-echelon inventory control system with complete redistribution. Management Science, 33(2): 215-227.), 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}7 DE KOK AG. 1990. Hierarchical production planning for consumer goods. European Journal of Operational Research, 45(1): 55-69.) 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}27 VAN DER HEIJDEN MC. 1997. Supply rationing in multi-echelon divergent systems. European Journal of Operational Research, 101(3): 532-549.) 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}16 LAGODIMOS AG & KOUKOUMIALOS S. 2008. Service performance of two-echelon supply chains under linear rationing. International Journal of Production Economics, 112(2): 869-884.) 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}2 BIRGE JR & LOUVEAUX F. 1997. Introduction to Stochastic Programming. Springer.), and especially when the premises with respect to the probabilistic phenomena are restrictive.

According to (^{Higle 2005}13 HIGLE JL. 2005. Stochastic Programming: Optimization When Uncertainty Matters. In: INFORMS Tutorials in Operations Research.), 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}5 CUNHA PSA, OLIVEIRA F & RAUPP FMP. 2017. A two-stage stochastic programing model for periodic replenishment control system under demand uncertainty. Computers & Industrial Engineering, 107: 313-326.) 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}11 GUPTA A & MARANAS CD. 2000. A Two-Stage Modeling and Solution Framework for Multisite Midterm Planning under Demand Uncertainty. Industrial & Engineering Chemestry Research, 39(10): 3799-3813.), (^{Santoso et al. 2005}22 SANTOSO T, AHMED S, GOETSCHALCKX M & SHAPIRO A. 2005. A stochastic programming approach for supply chain network design under uncertainty. European Journal of Operational Research, 167(1): 96-115.), (^{Oliveira & Hamacher 2012}21 OLIVEIRA F & HAMACHER S. 2012. Optimization of the petroleum product supply chain under uncertainty: A case study in northern Brazil. Industrial & Engineering Chemistry Research, 51(11): 4279-4287.) and (^{Oliveira et al. 2013}20 OLIVEIRA F, GUPTA V, HAMACHER S & GROSSMANN IE. 2013. A Lagrangean decomposition approach for oil supply chain investment planning under uncertainty with risk considerations. Computers & Chemical Engineering, 50(3): 184-195.). 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}6 DASKIN MS, COULLARD CR & SHEN Z-JM. 2002. An inventory-location model: formulation, solution algorithm and computational results. Annals of Operations Research, 110(1): 83-106.), (^{Shen et al. 2003}24 SHEN Z-JM, COULLARD C & DASKIN MS. 2003. A joint location-inventory model. Transportation Science, 37(1): 40-55.) and (^{You & Grossmann 2008}29 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.) addressed LN design and inventory policy optimization without the use of stochastic programming technique. Using stochastic programming, (^{Fattahi et al. 2015}9 FATTAHI M, MAHOOTCHI M, MOATTAR HUSSEINI SM, KEYVANSHOKOOH E & ALBORZI F. 2015. Investigating replenishment policies for centralised and decentralised supply chains using stochastic programming approach. International Journal of Production Research, 53(1): 41-69.) 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}5 CUNHA PSA, OLIVEIRA F & RAUPP FMP. 2017. A two-stage stochastic programing model for periodic replenishment control system under demand uncertainty. Computers & Industrial Engineering, 107: 313-326.) 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.

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:

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 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, ..., _^NI adopt respectively inventory control policies (R ₀, S ₀) and (_^Ri , _^Si ), where _^NI is the total of retailers in the distribution system, R ₀ and _^Ri represent the review interval at the DC and at each retailer i, and S ₀ 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 ₀ and _^Si , and the optimal review intervals R ₀ 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 uniform-sized time periods.

Let _^NP be the total of periods of a given planning horizon. The review intervals, or the times between orders, R ₀ and _^Ri , i = 1, ..., _^NI , are modeled as multiples of the period p. The lead times L ₀ 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 ₀ 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 ₀ 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 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 ₀) and (_{^Ri , ^Si} ) at DC and at each retailer i = 1, ..., ^_NI , are: the carrying costs per item per period, $h_0^p$ and $h_i^p$, and the ordering costs $C_{F_0}^p$ and $C_{F_i}^p$ 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_i^p$ 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 period, R ₀ = 2 periods and R ₁ = 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${(\xi )}_0^p$ and P${(\xi )}_0^{p+2}$, 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${(\xi )}_i^p$, P${(\xi )}_i^{p+1}$ and P${(\xi )}_i^{p+2}$. The first two orders are completely fulfilled by the DC at the beginning of periods p + 1 and p + 2, denoted by A${(\xi )}_0^p$ and A${(\xi )}_0^{p+1}$, 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 A${(\xi )}_0^{p+2}$, with the amount of the item in short given by F${(\xi )}_0^{p+2}$. This backorder is then fulfilled by the DC at the beginning of period p + 4. The demands from the customers denoted by D${(\xi )}_i^p$, D${(\xi )}_i^{p+1}$ and D${(\xi )}_i^{p+2}$ arrive in periods p, p + 1, and p + 2, being completely met in periods p and p + 2, and partially met in period p + 1. Hence, the retailer has a backorder in the amount F${(\xi )}_i^{p+1}$ 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.

3 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 ₀) and (_{^Ri , ^Si} ), i = 1, ..., ^_NI , 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 ₀) 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 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 ₃. 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 ₂, consider the condition as a constraint in the model that minimizes the ordering and carrying costs. We remark that the notation B ₃ and P ₂ are consistent with that used in the literature (see, for example, ^{Silver et al., 1998}26 SILVER EA, PYKE DF & PETERSON R. 1998. Inventory Management and Production Planning and Scheduling. 3 ed. John Wiley & Sons.).

As follows, in Subsections 3.1 and 3.2, we present in detail the models for the problem at hand with the B ₃ approach considering fixed and variable rationing rules. Then, in Subsection 3.3, the model with the P ₂ 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 ₃ - F

We propose a two-stage stochastic programming model to determine the optimal values R _0, S ₀, ^_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 fixed-fraction 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₀ possible review intervals at DC, r ₀ ∈ T ₀ = {1, …, $N_{R_0}$};

_^Ti possible review intervals at retailer i, _^ri ∈ _^Ti = {1, …, $N_{R_i}$};

Parameters

$b_i^p$ unit cost of the item in short in retailer i in period p;

$C_{F_0}^p$ ordering cost at DC in period p;

$C_{F_i}^p$ ordering cost at retailer i in period p;

D${(\xi )}_i^p$ demand at retailer i in scenario ξ in period p;

$h_i^p$ cost of carrying one unit in period p at retailer i;

$h_0^p$ cost of carrying one unit in period p at DC;

$\overline{ITI}$ auxiliary parameter to compute the quantities of the item ordered given in terms of the upper bound for the stock position;

$\overline{S}$ auxiliary parameter to compute the quantities of the item ordered given in terms of the upper bound for the inventory target level;

_^Vtb auxiliary parameter, _^Vtb ∈ $\left\{\frac{2^0}{{10}^y},\frac{2^1}{{10}^y},\frac{2^2}{{10}^y},\dots ,\frac{2^{N_B}}{{10}^y}\right\}$ and 2 > Σ_B$\frac{2^{tb}}{{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_{\mathrm{0,}p}^{r_0}$ auxiliary parameter that indicates the period that occurs an order at DC depending on the value r ₀; $w_{\mathrm{0,}p}^{r_0}$ ∈{0, 1}; r ₀ = 1, ..., $N_{R_0}$; p = 1, ..., ^_NP ;

$w_{i,p}^{r_i}$ auxiliary parameter that indicates the period that occurs an order at retailer i depending on the value _^ri ; $w_{i,p}^{r_i}$ ∈ {0, 1}; _^ri = 1, ..., $N_{R_i}$; p = 1, ..., _^NP;

W auxiliary _^Np × $N_{R_k}$ matrix of $w_{k,p}^{r_k}$ parameters, k = 0, ..., _^NI , defined as

Variables

A${(\xi )}_0^p$ 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${(\xi )}_i^p$ 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${(\xi )}_{\mathrm{0,}i}^p$ accumulated orders from retailers i fulfilled by DC in scenario ξ in period p;

F${(\xi )}_0^p$ accumulated orders from all retailers not fulfilled by DC in scenario ξ in period p;

F${(\xi )}_i^p$ accumulated demands not met by retailer i in scenario ξ in period p;

F${(\xi )}_{\mathrm{0,}i}^p$ accumulated orders of retailer i unmet by DC in scenario ξ in period p;

_^fi shortage fraction for retailer i, that is, orders placed by retailer i not fulfilled by the DC (_^fi = Σ_{^{B ji,tbVtb)}} ;

I${(\xi )}_i^p$ stock on hand of retailer i in scenario ξ at the end of period p;

I${(\xi )}_0^p$ stock on hand of the DC in scenario ξ at the end of period p;

_^Ie${(\xi )}_i^p$ echelon stock of retailer i in scenario ξ at the end of period p;

_^Ie${(\xi )}_0^p$ echelon stock of the DC in scenario ξ at the end of period p;

$I_I^e{\left(\xi \right)}_0^p$ echelon stock of the DC in scenario ξ at the beginning of period p;

$I_I^e{\left(\xi \right)}_i^p$ echelon stock of retailer i in scenario ξ at the beginning of period p;

$I{V}_I^e{\left(\xi \right)}_0^p$ auxiliary variable for the echelon stock of the DC in scenario ξ at the beginning of period p;

$I{V}_I^e{\left(\xi \right)}_i^p$ auxiliary variable for the echelon stock of retailer i in scenario ξ at the beginning of period p;

JF${(\xi )}_{i,tb}^p$ auxiliary variable representing the amount of unmet orders placed by retailers in scenario ξ in period p

_^Ji,tb auxiliary binary variable in approximating the binary representation of _^fi ;

P${(\xi )}_i^p$ quantity of the item ordered by retailer i in scenario ξ at the beginning of period p;

P${(\xi )}_0^p$ quantity of the item ordered by the DC in scenario ξ at the beginning of period p;

$S{V}_0^p$ auxiliary variable for the maximal inventory level of the item at DC in period p;

$S{V}_i^p$ auxiliary variable for the maximal inventory level of the item at retailer i in period p;

$u_0^{r_0}$ auxiliary binary variable for determining R ₀;

$u_i^{r_i}$ auxiliary binary variable in the determination of _^Ri ;

$v_0^p$ indicates if there exists an order for the item at the DC in period p; $v_0^p$ ∈ {0, 1};

$v_i^p$ indicates if there exists or not an order for the item at retailer i in period p; $v_i^p$ ∈ {0, 1};

X${(\xi )}_0^p$ indicates if there exists a shortage or stock on hand at DC in scenario ξ at the end of period p; X${(\xi )}_0^p$ ∈ {0, 1};

3.1.2 First-stage problem

The first-stage 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 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 ₀ and ^_Ri must be determined (R ₀ = r ₀ ∈ T ₀ = {1, …, $N_{R_0}$} and ^{_{Ri = ri}} ∈ Ti = {1, …, $N_{R_i}$}, when $u_0^{r_o}$ = 1 and $u_i^{r_i}$ = 1). Constraints (4) and (5) indicate that orders occur every interval R ₀ at the DC and every interval ^_Ri 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_{\mathrm{0,}p}^{r_0}$ and $w_{i,p}^{r_i}$). 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_k^{r_k},v_k^p$, k = 0, ..., ^_NI , are defined as binary.

3.1.3 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, ^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 second-stage 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 $v_i^p$ = 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 ₀ 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_0^p$ = 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${(\xi )}_{\mathrm{0,}i}^p$ in (21), we make two remarks. By definition of F${(\xi )}_{\mathrm{0,}i}^p$ in (22), depending on the value of _^fi , F${(\xi )}_{\mathrm{0,}i}^p$ can be lower than or equal to P${(\xi )}_i^p$ + F${(\xi )}_{\mathrm{0,}i}^{p-1}$, and thus A${(\xi )}_{\mathrm{0,}i}^p$ is positive and there is no imbalance (allocation of negative shortage to retailers). Otherwise, if F${(\xi )}_{\mathrm{0,}i}^p$ is greater than P${(\xi )}_i^p$ + F${(\xi )}_{\mathrm{0,}i}^{p-1}$, then A${(\xi )}_{\mathrm{0,}i}^p$ 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${(\xi )}_{\mathrm{0,}i}^p$. Moreover, in case A${(\xi )}_{\mathrm{0,}i}^p$ 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), we introduce the linearized version of the problem model as follows. First, we introduce the variables $I_I^e{(\xi )}_0^p$ and $I_I^e{(\xi )}_i^p$, ∀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_I^e{(\xi )}_0^p=I^e{(\xi )}_0^{p-1}$ and $I_I^e{(\xi )}_i^p=I^e{(\xi )}_i^{p-1}$, 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) by the expressions from (45) to (48), apart the introduction of non-negativity constraints for the auxiliary variables JF${(\xi )}_{i,bin}^p$, _^SVP , $I_I^e{(\xi )}_0^p$ and $I{V}_I^e{(\xi )}_0^p$. Finally, the second-stage problem corresponds to the following MILP model:

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 ₀ 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${(\xi )}_{\mathrm{0,}i}^p$, must be equal to the corresponding order P${(\xi )}_i^p$ plus the quantity in short in the earlier period F${(\xi )}_{\mathrm{0,}i}^{p-1}$ during L ₀:

3.2 Minimization of relevant costs with variable rationing fractions ( $f{(\xi )}_i^p$ ):Model B₃ - 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${(\xi )}_i^p$. Thus, model B ₃ - V differs from model B ₂ - V with respect to first-stage decisions. While in model B ₃ - F, the fraction ^_fi is a first-stage decision, in model B ₃ - V the fraction f${(\xi )}_i^p$ can 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${(\xi )}_i^p$ represents the demand of retailer i until period p.

3.2.1 Notation

Consider also the following notation:

Variables

dif (ξ)^P auxiliary variable to compute F${(\xi )}_i^p$

f${(\xi )}_i^p$ fraction between the need of retailer i and the needs of all retailers until period p of scenario ξ;

LF${(\xi )}_{i,tb}^p$ auxiliary variable for the amount of unmeet orders placed by retailers in scenario ξ in period p;

LNec${(\xi )}_{i,tb}^p$ auxiliary variable to compute Nec(ξ)^P ;

I${(\xi )}_{i,tb}^p$ auxiliary binary variable for the approximation of the binary representation of f${(\xi )}_i^p$;

Nec(ξ)^P total of all retailers’ demands until period p;

Nec${(\xi )}_i^p$ demand of retailer i until period p;

Using the binary expansion and a predefined precision of 1/10^y , where y is a fixed known value, we observe that f${(\xi )}_i^p$ is in the following interval:

As the approximate value given by the binary expansion ∑_{B ^Vtbl}${(\xi )}_{i,tb}^p$ of the fraction f${(\xi )}_i^p$ is always less than or equal to the corresponding true value plus the precision term, on the right-hand side of the inequality (57), thus ∑_{B ^Vtbl}${(\xi )}_{i,tb}^p$ 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 _^Vtbl${(\xi )}_{i,tb}^p$ considers the need of each retailer and thus a new binary representation is proposed:

As $Nec{(\xi )}^{p}l{(\xi )}_{i,tb}^p=LNec{(\xi )}_{i,tb}^p$ and $F{(\xi )}_0^{p}l{(\xi )}_{i,tb}^p=LF{(\xi )}_{i,tb}^p$, it follows that:

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 ₂

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 ₃, 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 ${\overline{f}}_i$ is the expected value of the fractions of the demands of retailer i promptly fulfilled and $F\text{'}{(\xi )}_i^p={\displaystyle \sum }_{p}F{(\xi )}_i^p-{\displaystyle \sum }_{p}F{(\xi )}_i^{p-1}$, and, as already defined, X${(\xi )}_i^p$ 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 ₂ - F and P ₂ - 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 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}22 SANTOSO T, AHMED S, GOETSCHALCKX M & SHAPIRO A. 2005. A stochastic programming approach for supply chain network design under uncertainty. European Journal of Operational Research, 167(1): 96-115.). 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, ξⁿ ) is the objective function of the second-stage problem to be evaluated in each subset M and scenario ξⁿ . Given a collection of subsets of scenarios independently generated by sampling $\left({\xi }_j^1,\dots ,{\xi }_j^N\right)$, j = 1, ..., M, we have:

where the value of the first-stage objective function is approximated by:

According to (^{Santoso et al. 2005}22 SANTOSO T, AHMED S, GOETSCHALCKX M & SHAPIRO A. 2005. A stochastic programming approach for supply chain network design under uncertainty. European Journal of Operational Research, 167(1): 96-115.), 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 first-stage model can be approximated by:

Given ${\xi }_{j\text{'}}^1,\dots ,{\xi }_{j\text{'}}^{N\text{'}}$, 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}22 SANTOSO T, AHMED S, GOETSCHALCKX M & SHAPIRO A. 2005. A stochastic programming approach for supply chain network design under uncertainty. European Journal of Operational Research, 167(1): 96-115.) also state that ${\widehat{\phi }}_{N\text{'}}$ is an unbiased estimator for $\widehat{\phi }$(R', S', f'). Since (R', S', f') is a feasible solution of the problem, $\widehat{\phi }$(R', S', f') is greater than or equal to the minimal value of the problem. And, since $\widehat{\phi }$_{N ', M '} is an unbiased estimator for $\widehat{\phi }$_{N '}, $\widehat{\phi }$_{N ', M '} is also greater than or equal to the optimal minimum value of the problem. Thus, $\widehat{\phi }$_{N ', M '} can be considered as an upper bound (UB) for the optimal minimum value of the original objective function.

(^{Linderoth et al. 2006}18 LINDEROTH J, SHAPIRO A & WRIGHT S. 2006. The empirical behavior of sampling methods for stochastic programming. Annals of Operations Research, 142(1): 215-241.) 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 ${\sigma }_{LB}^2$ and ${\sigma }_{UB}^2$ 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}15 KLEYWEGT AJ, SHAPIRO A & HOMEM-DE-MELO T. 2002. The sample average approximation method for stochastic discrete optimization. SIAM Journal on Optimization, 12(2): 479-502.; ^{Shapiro and Homem-de-Melo, 1998}23 SHAPIRO A & HOMEM-DE-MELLO T. 1998. A simulation-based approach to two-stage stochastic programming with recourse. Mathematical Programming, 81(3): 301-325.). 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 {ξ¹, ξ², ..., ξ^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 trade-off 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 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 ₃ - 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: $N_{R_0}$ = 3, r ₀ ∈ {1, 2, 3}, $N_{R_i}$ = 1, ^_ri =1, ∀i, $C_{F_0}^p$ = 200 and L ₀ = L ₁ = 1, ∀i. Also, the costs of carrying per unit of the item in period p are $h_0^p$ = 1, ∀p, and $h_i^p$ = 4, ∀i, p; and the shortage cost per unit of the item per period is $b_i^p$ = 10, ∀i, p.

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σ₅₀, with ^_ĝN = 298.89 and ${\widehat{\sigma }}_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 σ². 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 ₀ = 2. To obtain the upper bound (UB), as in (78), we considered a candidate solution given by the mean values for S ₀, _{^{Si , fi}} and r ₀ 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 ₃ - 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 ₀ , S ₀, ^{_{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).

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 ₀, ^{_{Si , fi, LB, UB}} , and the expected value of the fraction relative to the unmet demand (shortage fraction) obtained by models P ₂ - F and P ₂ - 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_0^p$ = 1 and $h_i^p$ = 4, ∀i, for I1, I3 and I4, and $h_0^p$ = 3.5 for I2; and the expected value of the fraction for the demand fully met at the retailers are ${\overline{f}}_i$ = 85%, 90%, 95% and 99%, ∀i, for I1 and I2, and for I3 and I4 we set ${\overline{f}}_1$ = 85%, ${\overline{f}}_2$ = 90% and ${\overline{f}}_3$ = 95%. In this experiment, the review intervals at DC and retailers were fixed respectively at R ₀ = 3 and R ₀ = 1, ∀i, with lead times L _{0 =} L ₁ = 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 σ². 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 _^NP = 30. To get the upper bound, as in (78), a candidate solution is set with the mean values for S ₀, _^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 ₂ - F and P ₂ - V models in terms of the instance size of their equivalent deterministic MILP model, and the computational time in terms of CPU time.

Tables 4, 5, 6 and 7 show the comparative results of the models P ₂ - F and P ₂ - V with respect to S ₀, ^{_{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 - ${\overline{f}}_i$) indicates the maximum limit for the expect value of the fraction of unmet demand. The last 3 columns, 1 - ${\overline{f}}_1$, 1 - ${\overline{f}}_2$ and 1 - ${\overline{f}}_3$, show the service level achieved for each instance with the experiment, considering a candidate solution with mean values of S ₀ and ^_Si , and the mean fraction of unmet demand given by

for P ₂ - V, and S ₀, ^_Si and^{_{, fi}} for P ₂ - F, after 100 runs (M = 100), with 30 scenarios (N = 30) and 30 periods (^_NP = 30).

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 pre-set values for the service levels, whereas the values obtained for the mean fractions of unmet demands $f_i^m$ do not change for the majority of the values of the service levels. This is due to the fact that _^fi is a first-stage decision variable whose value is determined through the optimization process, while the value of $f_i^m$ is imposed by the variable-fractions rationing rule. Also, from Tables 4 and 5, we verify that S ₀ values decrease and _^Si values increase as $h_0^p$ increases for P ₂ - F, as well as for P ₂ - V. Table 6 shows that the results for P ₂ - 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 - ${\overline{f}}_i$ = 1%) I2 and I3, and between 2% and 3% in the other cases.

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 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.