# Exploring the Potential Use of the Birnbaum-Saunders Distribution in Inventory Management

Authors: Peter Wanke; Víctor Leiva
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

Choosing the suitable demand distribution during lead-time is an important issue in inventory models. Much research has explored the advantage of following a distributional assumption different from the normality. The Birnbaum-Saunders (BS) distribution is a probabilistic model that has its genesis in engineering but is also being widely applied to other fields including business, industry, and management. We conduct numeric experiments using the R statistical software to assess the adequacy of the BS distribution against the normal and gamma distributions in light of the traditional lot size-reorder point inventory model, known as (Q, r). The BS distribution is well-known to be robust to extreme values; indeed, results indicate that it is a more adequate assumption under higher values of the lead-time demand coefficient of variation, thus outperforming the gamma and the normal assumptions.


---

## Body

## 1. Introduction and Bibliographical Review

Inventory management permeates decision-making in countless firms. The topic has been extensively studied in academic and corporate spheres, for example, Braglia et al. [1] and Cai et al. [2]. The key questions which the inventory management seeks to answer, usually influenced by a variety of circumstances, are as follows: when to order, determining an economic order quantity (EOQ) or lot size, and how much safety stock (SS) to keep, establishing a reorder point (ROP); see Namit and Chen [3] and Porras and Dekker [4].According to Wanke [5], inventory management involves a set of decisions whose objective is to match existing demand with the supply of products and materials over space and time. This objective allows us to achieve specified costs and service levels, considering product, operation, and demand characteristics. It is known that the inventory total cost (TC) is a function of ordering, holding, and shortage costs; see Hillier and Lieberman [6].The importance attached by firms to inventory management can be attributed to the following: first and foremost, it is the need to ensure that products, given the competitive pressure exercised by markets, are always supplied to customers at the least possible cost; see Eaves [7]. Second, some other factors contribute to a high concern with inventory management, such as product diversity or behavior; see Huiskonen [8]. High opportunity costs also contribute to this concern, thus affecting the financial indicators on which assessments of firm performance are based; see Wanke [5].The inventory management models are frequently classified in two types: pull and push. On the one hand, according to Ballou and Burnetas [9], pull-type planning models range from those that set inventory levels based on the EOQ to those fixed in proportion to forecasted demand. The EOQ model is the simplest and most fundamental of all inventory models because it describes important trade-offs between fixed ordering and holding costs; see Nahmias [10]. Despite its shortcomings, the basic EOQ model is the cornerstone of several software packages for inventory control; see Lee and Nahmias [11]. Interested readers can refer to Yan and Wang [12] and Min et al. [13] for more details about the EOQ model. Today EOQ is used in conjunction with ROP in inventory control models to determine cycle and SS under demand per unit of time (DPUT) and lead-time (LT) uncertainty. These models are well described in most logistics and operations textbooks, as are their underlying assumptions; see Nahmias [10]. Models that set inventory levels in proportion to forecasted demand constitute a particular form of the periodic review model, except that replenishment quantities are not based on the EOQ model; see Ballou [14]. On the other hand, push-type planning models take place when inventory decisions are based on the demand or its forecast at multiple downstream stocking locations, similar to the resource planning system of logistics.According to Silver et al. [15], Eaves and Kingsman [16], Syntetos et al. [17], and Boylan et al. [18], demand is one of the main factors in inventory management models. In general, demand may be classified from two settings. First, it can be deterministic or random; if random, the demand follows a statistical distribution (also known as probabilistic model); otherwise, the demand is constant, which implies a degenerate statistical distribution; that is, its variance is equal to zero. Second, demand might be independent or dependent. For more details, see Disney et al. [19], Porras and Dekker [4], Wanke [5], and Rojas et al. [20].Demand uncertainties directly affect the operation of the physical system of logistics. Moreover, to be closer to reality, single or multiple period inventory models must take into account that demand is occurring in a random fashion, which is explained by several factors. Thus, DPUT is taken to be a random variable (RV). Furthermore, during LT, due to the mentioned randomness, the corresponding demand (LTD) is also a RV; therefore, the behavior of DPUT and LTD must be described by statistical distributions; see Johnson et al. [21, 22]. The Gaussian (or normal) distribution is often used for describing the data of these two RVs (DPUT, LTD) involved in inventory models. However, it is well-known that the normal distribution is validly used for RVs that take negative and positive values with a symmetrical behavior. Hence, first, quantities less than zero could be admitted when the modeling is carried out under the normal distribution, which is not possible in real-world situations for DPUT and LTD, because they only admit values greater than zero; see Nahmias [10]. Second, another drawback using the normal model is that DPUT and LTD data often follow asymmetric distributions; see Moors and Strijbosch [23]. Mentzer and Krishnan [24] studied the nonnormality effect on inventory models and found that the normal distribution is appropriate in few practical cases; see also Eppen and Martin [25]. A recent case study with DPUT data of 89 food products supports such nonnormality; see Leiva et al. [26] and Rojas et al. [20]. In any case, the normality assumption must be checked by goodness-of-fit methods; see Barros et al. [27]. Thus, the use of the normal distribution to model DPUT and LTD and then to determine the ROP and SS can lead to wrong results, resulting in shortages or excess inventories. Nonnormal distributions with positive support that have been used for describing DPUT in inventory management include models such as gamma or Erlang, inverse Gaussian, log-normal, Pearson, Poisson, uniform, and Weibull; see Burgin [28], Tadikamalla [29], Lau [30], Wanke [31], Cobb et al. [32], and Pan et al. [33].A unimodal, two-parameter probability model with positive support and asymmetry to the right that is receiving considerable attention is the Birnbaum-Saunders (BS) distribution; see Birnbaum and Saunders [34] and Johnson et al. [22, pages 651–663]. The BS distribution has good properties and is related to the normal distribution and implemented in the R statistical software (http://www.r-project.org) via a package calledgbs; see R Team [35]. Although the BS distribution has its genesis from engineering, its applications range across diverse fields as business, industry, and management, which have been conducted by an international, transdisciplinary group of researchers; see, for example, Jin and Kawczak [36], Podlaski [37], Bhatti [38], Lio et al. [39], Paula et al. [40], Marchant et al. [41], and Leiva et al. [42, 43]. In addition, although originally conceived as a count model, the BS distribution includes the duration of the counting period (daily or weekly), which obviates having to collect additional data, among other properties; see Fox et al. [44]. In sum, the BS distribution is a good candidate for describing demand data in inventory models; see Leiva et al. [26] and Rojas et al. [20].Our main objective is to explore the use of the BS distribution in inventory management. Differently from previous studies that exclusively considered the effects of one given distribution on inventory decision-making, we also analyze its adequacy in light of different operating characteristics and costs. Specifically, we assess how the BS, gamma, and normal LTD distributions interact with relevant product characteristics and affect the optimal EOQ and SS inventory indicators in terms of the optimization of the TC function. We minimize this function using stochastic programming, a technique where constraints and/or objective function of the problem to be optimized contain RVs that can follow any distribution; see Shapiro et al. [45] and Thangaraj et al. [46]. We solve the problem of stochastic programming with a search heuristic called differential evolution (DE), which is a global numerical optimization approach based on genetic algorithm concepts; see Storn and Price [47] and Price et al. [48]. We implement our results inR code, which is available upon request from the authors.Section2 reviews general aspects of inventory management models and the statistical distributions used in this study. Section 3 explores the effect of different LTD distributions in inventory management, introducing the simulation scenario, formulating the stochastic programming model, discussing the DE algorithm, and providing a numerical study. Section 4 concludes the study making some considerations on the management of our findings and on future research.


## 2. Background

In this section, we discuss general aspects of inventory management models and demand statistical distributions used for the methodology presented in Section3.


### 2.1. Inventory Management Models

The(
Q
,
r
) model is based on the ROP (r) and the EOQ model given by(1)
Q
=
2
λ
P
H
1
/
2
,where λ is the DPUT rate in units of the product and P, H are the ordering and holding costs, respectively; see Yan and Wang [12] and Min et al. [13]. However, as mentioned, DPUT is a RV. Then, λ given in (1) must be calculated as the mean (expected value) of the DPUT distribution that adequately fits the data. Specifically, let D
t be a RV corresponding to the DPUT at time t, forming a sequence of independent and identically distributed RVs with mean E
(
D
t
)
=
λ and variance V
a
r
(
D
t
)
=
σ
D
2. In addition, let L be a RV corresponding to the LT between the ordering of a product and its delivery (expressed in time units) with mean E
(
L
)
=
μ
L and variance V
a
r
(
L
)
=
σ
L
2. Then, the LTD is given by(2)
X
=
∑
t
=
1
L
D
t
,with probability density function (PDF) f
X
(
·
) and whose expectation and variance are, respectively, defined as(3)
μ
=
E
X
=
μ
L
λ
,
σ
2
=
V
a
r
X
=
σ
L
2
λ
2
+
μ
L
σ
D
2
.The ROP can be computed from μ expressed in (3). However, to be protected from randomness of the LTD, it is necessary to include a SS, which allows the ROP to become(4)
r
=
μ
+
k
σ
,where μ, σ are defined in (3) and k is the safety factor (SF) or number of standard deviations (SDs) σ of the LTD. Note that although the (
Q
,
r
) model is based on (1) and (4), it is possible to see that r is obtained by k. Thus, we refer to the (
Q
,
r
) model as (
Q
,
k
) thereafter.The expected TC of the inventory is given by(5)
C
Q
,
k
=
H
Q
2
+
k
σ
+
λ
Q
P
+
S
∫
μ
+
k
σ
∞
f
X
u
d
u
,where Q (in units of the product) is given in (1) and μ, k, and σ are given in (4); H is the holding cost (in $ per $ per unit of time); P is the ordering cost (in $ per each replenishment order placed); S is the shortage cost (in $ incurred whenever a stock-out occurs); and f
X
(
·
) is the PDF of the LTD given in (2). In order to minimize the expected TC defined in (5), we optimize the indicator Q given in (1) altogether with k given in (4).


### 2.2. Demand Distributions

Notice that it is necessary to specify the LTD distribution to determine the SS given in (4), which allows the SF to be established; see Porras and Dekker [4]. In order to facilitate the calculation of the SF, the LTD has been traditionally modeled with the normal distribution; see Silver and Peterson [49]. Thus, the SF k for a specific service level can be obtained from a percentile of the standard normal distribution, denoted by N
(
0,1
). However, as mentioned, various studies criticize the normality assumption. Therefore, the use of the normal distribution to determine ROP and SS given in (4) is questionable, leading to possible stock shortage or excess.Silver [50] pointed out that in most models leading to inventory management decisions some assumptions are made, often in an implicit way. The effects of these assumptions on costs and service levels should be taken into account. The most common ones are (i) to assume a demand distribution (e.g., normal) and (ii) to suppose that the distribution parameters are known (e.g., the mean and SD) or estimated from the demand data. Lau [30] presented a model for computing EOQs and SSs given in (1) and (4), respectively, using the first four moments, that is, mean, variance, third moment reflecting skewness, and fourth moment reflecting kurtosis of any given LTD distribution. Lau [30] also pointed out the risk of misleading decisions regarding ROP and customer service level when one considers a normally distributed LTD. The 95th percentile of the distribution is often used to set service levels. Next, we present some mathematical features for the three LDT distributions to be considered in this study, that is, the BS, gamma, and normal models.The Normal Distribution. A RV X following a normal distribution with mean E
(
X
)
=
μ
∈
R and variance V
a
r
(
X
)
=
σ
2
>
0 is denoted by X
~
N
(
μ
,
σ
2
), where “~” means “distributed as”. In this case, PDF, cumulative distribution function (CDF), and (QF) quantile function of X are, respectively, (6)
f
X
x
=
1
2
π
σ
2
1
/
2
exp
⁡
-
1
2
x
-
μ
2
σ
2
,
x
∈
R
,
F
X
x
=
Φ
x
-
μ
σ
,
x
∈
R
,
x
q
=
z
q
σ
+
μ
,
0
<
q
<
1
,where Φ
(
z
)
=
∫
-
∞
z
‍
ϕ
(
u
)
d
u and z
(
q
)
=
Φ
-
1
(
q
), for 0
<
q
<
1, with (7)
ϕ
z
=
1
2
π
1
/
2
exp
⁡
-
1
2
z
2
,
z
∈
R
,and Φ
-
1
(
·
) being the N
(
0,1
) inverse CDF or QF. In addition, the coefficients of variation (CV), skewness or asymmetry (CS), and kurtosis (CK) of X
~
N
(
μ
,
σ
2
) are, respectively, (8)
CV
X
=
σ
μ
,
CS
X
=
0
,
CK
X
=
3
.The BS Distribution. A RV X following a BS distribution with shape α
>
0 and scale β
>
0 parameters is denoted by X
~
B
S
(
α
,
β
). In this case, the PDF, CDF, and QF of X are, respectively,(9)
f
X
x
=
1
2
π
1
/
2
exp
⁡
-
1
2
α
2
ξ
2
x
β
x
+
β
2
α
β
1
/
2
x
3
/
2
,
x
>
0
,
F
X
x
=
Φ
1
α
ξ
x
β
,
x
>
0
,
x
q
=
F
-
1
q
=
β
α
z
q
2
+
α
z
q
2
2
+
1
1
/
2
2
,
0
<
q
<
1
,where ξ
(
y
)
=
y
1
/
2
-
y
-
1
/
2
=
2
sinh
⁡
(
log
⁡
(
y
1
/
2
)
), for y
>
0, Φ
(
·
) is the N
(
0,1
) CDF, z
(
q
) is the N
(
0,1
) QF, and F
-
1
(
·
) is the inverse CDF of X. Note that x
(
0.5
)
=
β; that is, β is also the median or 50th percentile of the distribution. The mean, variance, CV, CS, and CK of X
~
B
S
(
α
,
β
) are, respectively, (10)
E
X
=
β
1
+
α
2
2
,
Var
X
=
β
2
α
2
1
+
5
α
2
4
,
CV
X
=
α
4
+
5
α
2
1
/
2
2
+
α
2
,
CS
X
=
4
α
6
+
11
α
2
4
+
5
α
2
3
/
2
,
CK
X
=
3
+
6
α
2
40
+
93
α
2
4
+
5
α
2
2
.In addition, the RVs X
~
B
S
(
α
,
β
) and Z
~
N
(
0,1
) are related by (11)
X
=
β
α
Z
2
+
α
Z
2
2
+
1
1
/
2
2
,
Z
=
1
α
ξ
X
β
.Also, note that W
=
Z
2 follows a chi-squared distribution with one degree of freedom. The BS distribution holds the scale and reciprocation properties; that is, (i) c
X
~
B
S
(
α
,
c
β
), with c
>
0, and (ii) 1
/
X
~
B
S
α
,
1
/
β, respectively.The Gamma Distribution. A RV X following a gamma distribution with shape α
>
0 and scale β
>
0 parameters is denoted by X
~
G
a
m
m
a
(
α
,
β
). In this case, the PDF and CDF of X are, respectively, (12)
f
X
x
=
β
α
Γ
α
x
α
-
1
exp
⁡
-
β
x
,
x
>
0
,
F
X
x
=
γ
α
,
β
x
Γ
α
,
x
>
0
,where Γ
(
·
) and γ
(
·
) stand for the usual and incomplete gamma functions, respectively. The corresponding QF given by x
(
q
)
=
F
-
1
(
q
), for 0
<
q
<
1, must be obtained by solving this equation with an iterative numerical method. The mean, variance, CV, CS, and CK of X are, respectively, (13)
E
X
=
α
β
,
Var
X
=
α
β
2
,
CV
X
=
1
α
1
/
2
,
CS
X
=
2
α
1
/
2
,
CK
X
=
3
+
6
α
.The gamma distribution also shares the scale property; that is, c
X
~
G
a
m
m
a
(
α
,
c
β
), with c
>
0.


## 2.1. Inventory Management Models

The(
Q
,
r
) model is based on the ROP (r) and the EOQ model given by(1)
Q
=
2
λ
P
H
1
/
2
,where λ is the DPUT rate in units of the product and P, H are the ordering and holding costs, respectively; see Yan and Wang [12] and Min et al. [13]. However, as mentioned, DPUT is a RV. Then, λ given in (1) must be calculated as the mean (expected value) of the DPUT distribution that adequately fits the data. Specifically, let D
t be a RV corresponding to the DPUT at time t, forming a sequence of independent and identically distributed RVs with mean E
(
D
t
)
=
λ and variance V
a
r
(
D
t
)
=
σ
D
2. In addition, let L be a RV corresponding to the LT between the ordering of a product and its delivery (expressed in time units) with mean E
(
L
)
=
μ
L and variance V
a
r
(
L
)
=
σ
L
2. Then, the LTD is given by(2)
X
=
∑
t
=
1
L
D
t
,with probability density function (PDF) f
X
(
·
) and whose expectation and variance are, respectively, defined as(3)
μ
=
E
X
=
μ
L
λ
,
σ
2
=
V
a
r
X
=
σ
L
2
λ
2
+
μ
L
σ
D
2
.The ROP can be computed from μ expressed in (3). However, to be protected from randomness of the LTD, it is necessary to include a SS, which allows the ROP to become(4)
r
=
μ
+
k
σ
,where μ, σ are defined in (3) and k is the safety factor (SF) or number of standard deviations (SDs) σ of the LTD. Note that although the (
Q
,
r
) model is based on (1) and (4), it is possible to see that r is obtained by k. Thus, we refer to the (
Q
,
r
) model as (
Q
,
k
) thereafter.The expected TC of the inventory is given by(5)
C
Q
,
k
=
H
Q
2
+
k
σ
+
λ
Q
P
+
S
∫
μ
+
k
σ
∞
f
X
u
d
u
,where Q (in units of the product) is given in (1) and μ, k, and σ are given in (4); H is the holding cost (in $ per $ per unit of time); P is the ordering cost (in $ per each replenishment order placed); S is the shortage cost (in $ incurred whenever a stock-out occurs); and f
X
(
·
) is the PDF of the LTD given in (2). In order to minimize the expected TC defined in (5), we optimize the indicator Q given in (1) altogether with k given in (4).


## 2.2. Demand Distributions

Notice that it is necessary to specify the LTD distribution to determine the SS given in (4), which allows the SF to be established; see Porras and Dekker [4]. In order to facilitate the calculation of the SF, the LTD has been traditionally modeled with the normal distribution; see Silver and Peterson [49]. Thus, the SF k for a specific service level can be obtained from a percentile of the standard normal distribution, denoted by N
(
0,1
). However, as mentioned, various studies criticize the normality assumption. Therefore, the use of the normal distribution to determine ROP and SS given in (4) is questionable, leading to possible stock shortage or excess.Silver [50] pointed out that in most models leading to inventory management decisions some assumptions are made, often in an implicit way. The effects of these assumptions on costs and service levels should be taken into account. The most common ones are (i) to assume a demand distribution (e.g., normal) and (ii) to suppose that the distribution parameters are known (e.g., the mean and SD) or estimated from the demand data. Lau [30] presented a model for computing EOQs and SSs given in (1) and (4), respectively, using the first four moments, that is, mean, variance, third moment reflecting skewness, and fourth moment reflecting kurtosis of any given LTD distribution. Lau [30] also pointed out the risk of misleading decisions regarding ROP and customer service level when one considers a normally distributed LTD. The 95th percentile of the distribution is often used to set service levels. Next, we present some mathematical features for the three LDT distributions to be considered in this study, that is, the BS, gamma, and normal models.The Normal Distribution. A RV X following a normal distribution with mean E
(
X
)
=
μ
∈
R and variance V
a
r
(
X
)
=
σ
2
>
0 is denoted by X
~
N
(
μ
,
σ
2
), where “~” means “distributed as”. In this case, PDF, cumulative distribution function (CDF), and (QF) quantile function of X are, respectively, (6)
f
X
x
=
1
2
π
σ
2
1
/
2
exp
⁡
-
1
2
x
-
μ
2
σ
2
,
x
∈
R
,
F
X
x
=
Φ
x
-
μ
σ
,
x
∈
R
,
x
q
=
z
q
σ
+
μ
,
0
<
q
<
1
,where Φ
(
z
)
=
∫
-
∞
z
‍
ϕ
(
u
)
d
u and z
(
q
)
=
Φ
-
1
(
q
), for 0
<
q
<
1, with (7)
ϕ
z
=
1
2
π
1
/
2
exp
⁡
-
1
2
z
2
,
z
∈
R
,and Φ
-
1
(
·
) being the N
(
0,1
) inverse CDF or QF. In addition, the coefficients of variation (CV), skewness or asymmetry (CS), and kurtosis (CK) of X
~
N
(
μ
,
σ
2
) are, respectively, (8)
CV
X
=
σ
μ
,
CS
X
=
0
,
CK
X
=
3
.The BS Distribution. A RV X following a BS distribution with shape α
>
0 and scale β
>
0 parameters is denoted by X
~
B
S
(
α
,
β
). In this case, the PDF, CDF, and QF of X are, respectively,(9)
f
X
x
=
1
2
π
1
/
2
exp
⁡
-
1
2
α
2
ξ
2
x
β
x
+
β
2
α
β
1
/
2
x
3
/
2
,
x
>
0
,
F
X
x
=
Φ
1
α
ξ
x
β
,
x
>
0
,
x
q
=
F
-
1
q
=
β
α
z
q
2
+
α
z
q
2
2
+
1
1
/
2
2
,
0
<
q
<
1
,where ξ
(
y
)
=
y
1
/
2
-
y
-
1
/
2
=
2
sinh
⁡
(
log
⁡
(
y
1
/
2
)
), for y
>
0, Φ
(
·
) is the N
(
0,1
) CDF, z
(
q
) is the N
(
0,1
) QF, and F
-
1
(
·
) is the inverse CDF of X. Note that x
(
0.5
)
=
β; that is, β is also the median or 50th percentile of the distribution. The mean, variance, CV, CS, and CK of X
~
B
S
(
α
,
β
) are, respectively, (10)
E
X
=
β
1
+
α
2
2
,
Var
X
=
β
2
α
2
1
+
5
α
2
4
,
CV
X
=
α
4
+
5
α
2
1
/
2
2
+
α
2
,
CS
X
=
4
α
6
+
11
α
2
4
+
5
α
2
3
/
2
,
CK
X
=
3
+
6
α
2
40
+
93
α
2
4
+
5
α
2
2
.In addition, the RVs X
~
B
S
(
α
,
β
) and Z
~
N
(
0,1
) are related by (11)
X
=
β
α
Z
2
+
α
Z
2
2
+
1
1
/
2
2
,
Z
=
1
α
ξ
X
β
.Also, note that W
=
Z
2 follows a chi-squared distribution with one degree of freedom. The BS distribution holds the scale and reciprocation properties; that is, (i) c
X
~
B
S
(
α
,
c
β
), with c
>
0, and (ii) 1
/
X
~
B
S
α
,
1
/
β, respectively.The Gamma Distribution. A RV X following a gamma distribution with shape α
>
0 and scale β
>
0 parameters is denoted by X
~
G
a
m
m
a
(
α
,
β
). In this case, the PDF and CDF of X are, respectively, (12)
f
X
x
=
β
α
Γ
α
x
α
-
1
exp
⁡
-
β
x
,
x
>
0
,
F
X
x
=
γ
α
,
β
x
Γ
α
,
x
>
0
,where Γ
(
·
) and γ
(
·
) stand for the usual and incomplete gamma functions, respectively. The corresponding QF given by x
(
q
)
=
F
-
1
(
q
), for 0
<
q
<
1, must be obtained by solving this equation with an iterative numerical method. The mean, variance, CV, CS, and CK of X are, respectively, (13)
E
X
=
α
β
,
Var
X
=
α
β
2
,
CV
X
=
1
α
1
/
2
,
CS
X
=
2
α
1
/
2
,
CK
X
=
3
+
6
α
.The gamma distribution also shares the scale property; that is, c
X
~
G
a
m
m
a
(
α
,
c
β
), with c
>
0.


## 3. Assessing the Impact of Different Distributions

In this section, we introduce our simulation scenario and formulate the stochastic programming used to optimize the expected TC associated with the(
Q
,
k
) model. Then, we discuss the DE algorithm, which allows us to solve the problem of stochastic programming, and provide a numerical study performed with theR software. We evaluate how different LTD distributions (BS, gamma, and normal) interact with different inventory indicators (demand, cost, and LT) and their underlying EOQs and SSs. We assess under what circumstances a distributional assumption is preferable to the other one in terms of the expected TC.


### 3.1. Scenario of the Simulation Study

Assume BS, gamma, and normal distributions for the LTD. Then, fix values for the parameters of these distributions by considering values for means and SDs of DPUT and LT generated from uniform distributions. Now, generate holding, ordering, and shortage costs also from uniform distributions. This allows us to establish the expected TC to be minimized. Ten thousand (10000) different simulated scenarios of means and SDs for DPUT and LT as well as holding, ordering, and shortage costs are generated using anR package calledstats. The values of these uniformly distributed inventory indicators used to build the scenarios are presented in Table 1. They were chosen based on values proposed in selected papers focused on managerial and industrial applications compiled by Wanke [5]; see Table  2 in this reference for more details about values that are frequently used to generate simulation scenarios in inventory management problems. A discussion about this can be found in Wanke [51].Table 1
Range of the mentioned uniformly distributed indicator for simulations.
Indicator
Minimum
Maximum
DPUT mean (units/day)
80
120
DPUT SD (units/day)
3
30
LT mean (in days)
1
5
LT SD (in days)
0.50
2.00
Holding cost ($/unit/day)
0.00
0.68
Ordering cost ($/order)
17
60
Shortage cost ($/shortage)
0
100


### 3.2. Stochastic Programming

Once the values for the inventory indicators are defined and the distributional assumptions (BS, gamma, and normal) for the LTD established, stochastic programming is performed on the expected TC function given in (5); see Namit and Chen [3]. According to the values provided in Table 1, we assume values for means and SDs of the DPUT and LT, as well as for the holding, ordering, and shortage costs, to be uniformly distributed in the objective function corresponding to the expected TC given in (5). The decision variables of the programming are Q and k given in (1) and (4), respectively, whereas λ, σ
D, μ
L, σ
L, H, P, and S are given from (5). Therefore, our optimization problem of the expected TC can be visualized like a model of stochastic programming formulated as(14)
Minimize
Z
=
C
Q
,
k
Subject  to:
Q
>
0
,
k
>
0
,
λ
~
U
80,120
,
σ
D
~
U
3,30
,
μ
L
~
U
1,5
,
σ
L
~
U
0.50,20
,
H
~
U
0.00,0.68
,
P
~
U
17,60
,
S
~
U
0,100
.The problem of stochastic programming formulated in (14) is aimed (i) at identifying the most adequate inventory policy for each of the distributional assumptions and (ii) at minimizing the expected TC. The optimized problem can provide useful information for academics and practitioners on how these assumptions interact with product, operation, and demand characteristics. To solve the problem of stochastic programming formulated in (14), we use the DE algorithm detailed in the next section.


### 3.3. Differential Evolution

DE is a member of the family of genetic algorithms, which mimic the process of natural selection in an evolutionary manner; see Holland [52]. A genetic algorithm solves optimization problems with biology-inspired operators of crossover, mutation, and selection, generating successive populations of individuals (solutions or generations). Then, the DE algorithm optimizes problems by evolving a population of candidate solutions employing the mentioned operators. The DE algorithm uses floating-point techniques for obtaining the solutions and arithmetic operations in their mutation, in contrast to classic genetic algorithms. In addition, the DE algorithm finds the global optimum of the objective function, which is not required to be either continuous or differentiable; see Thangaraj et al. [46] and Mullen et al. [53].The DE algorithm has also been used to optimize problems that arise in inventory management, such as joint replenishment, replenishment coordination, and inventory location-allocation; see Qu et al. [54, 55] and Wang et al. [56]. In the present study, the problem of stochastic programming formulated in (14) is solved with the DE algorithm, which allows us to find the optimal values of Q and k that minimize the expected TC function C
(
·
) given in (5), for 10000 simulated scenarios.In what follows, we first discuss the DE algorithm used in our research work; see Storn and Price [47] and Price et al. [48] for more details about the algorithm. Then, we sketch the content of anR package calledDEoptim, which implements the DE algorithm and was first published on CRAN in 2005. Since becoming publicly available, it has been used by several authors to solve optimization problems arising in diverse domains. We refer interested readers to Ardia et al. [57] and Mullen et al. [53] for a detailed description of the package.LetN be the number of members y
∈
R
d (tuning parameter vector) in the population, where d denotes the dimension of the vector y, in our case given by y
=
(
Q
,
k
)
⊤
∈
R
2. The DE algorithm needs a starting population, which is obtained by sampling the objective function at multiple randomly chosen initial points (generation or solution zero (0)) for Q and k. Before the population is initialized, both lower and upper bounds for each tuning parameter must be specified. Parameter bounds establish the domain from which the N vectors at the generation 0 are chosen. To establish the generation 0, N guesses for the optimal value of y must be provided, using either random values within a range defined by the practitioner or values fixed by him/her. Each generation creates a new population from the current population members {
y
j
,
g
,
j
=
1
,
…
,
N
,
g
=
1,2
,
…
}, where j indexes the tuning parameter vector that make up the population and g indexes the generation. The new generation is obtained using differential mutation of the population members. In our research, N
=
10000 denotes the number of simulated scenarios and y denotes the tuning parameters for the lot size and the number of demand SDs. An initial mutant parameter vector {
v
j
,
g
,
j
=
1
,
…
,
N
,
g
=
1,2
,
…
} is created by choosing three members of the population y
j
1
,
g, y
j
2
,
g, and y
j
3
,
g, at random. Then, the elements of the initial mutant parameter vector v
j
,
g
=
(
v
j
i
,
g
) are generated by(15)
v
j
i
,
g
=
y
j
1
,
g
+
ζ
y
j
2
,
g
-
y
j
3
,
g
,
i
=
1
,
…
,
d
,where d
=
2 in our case and ζ is a positive scale factor whose effective value is usually less than one (usual default: ζ
=
0.8). After the first mutation operation, it continues until k mutations have been made or until certain crossover probability (CP) in [
0
,
1
] is less than u, where u is a random number from the uniform distribution in [
0
,
1
]. CP controls the fraction of the tuning parameter values that are copied from the mutant and approximates the probability that a parameter value is inherited from the mutant, since at least one mutation always occurs. Mutation is applied in this way to each member of the population.Calculations in our simulation study were performed with the aid of theDEoptim package, which consists of the core functionDEoptim() whose arguments are as follows:(i)
fn is the function to be minimized, which must have as its first argument the vector of real-valued parameters to optimize and return a scalar real result.
(ii)
lower,upper correspond to two vectors establishing scalar real lower and upper bounds on each tuning parameter to be optimized; the ith element of thelower andupper vectors corresponds to the ith parameter; the implementation searches the global optimum offn betweenlower andupper.
(iii)
… allows the user to pass additional arguments to the functionfn.
(iv)
control is a list whose default value is the return value ofDEoptim.control() but whose main elements are interpreted as follows:
(a)
VTR specifies the global minimum offn if it is known or if you wish to cease optimization after having reached a certain value (default =-Inf).
(b)
strategy defines the differential evolution strategy used in the optimization procedure, described in detail by Mullen et al. [53].
(c)
NPis the number of population members (d
e
f
a
u
l
t
=
10
×
d or 50), in our case denoted by N.
(d)
bs: ifbs isFALSE, then every mutant is tested against a member in the previous generation, and the best value survives into the next generation; ifbs isTRUE, then the old generation andNP mutants are sorted by their associated objective function values, and the bestNP vector proceeds into the next generation (default =FALSE).
(e)
itermax is the maximum number of iterations (i.e., population generations) to be allowed (default = 200).
(f)
CR is the CP from the interval [
0
,
1
] (default = 0.9).
(g)
F is the stepsize from interval [
0
,
2
], in our case denoted by ζ (default = 0.8).
(h)
trace,initialpop,storepopfrom,storepopfreq,checkWinner, andavWinner are other elements of the listcontrol, which are described in detail by Mullen et al. [53]; see also Price et al. [48].The return value of the functionDEoptim() is a member of the S3 classDEoptim. Members of this class have aplot and asummary that allow us to analyze the optimizer’s output. In our application of theDEoptim package, decision variables Q and k were lower-bounded at zero. As regards the list of tuning parameters used in our research, the default values for theDEoptim package were considered.


### 3.4. Numerical Results and Discussion

The underlying idea of performing a sensitivity analysis on the testing variables related to product, demand, and operational characteristics is to discriminate between groups where the three distributional assumptions led to minimal TCs. Table2 summarizes the numerical experiments conducted with theR software for 10000 simulated scenarios using the DE algorithm. Note that the BS distribution yields, on average, smaller TCs, EOQs, SSs, and, therefore, smaller inventory levels, in comparison to the normal and gamma assumptions.Table 2
Summary of the simulations for the mentioned distribution and management indicator.
Indicator
Distribution
Gamma
Normal
BS
Sum of TCs (in $)
205097.63
209175.49
202134.31
Sum of EOQs
767449.00
764816.00
762843.00
Sum of SSs
41166.66
53425.51
38103.58
Sum of average inventory level
424891.16
435833.51
419525.08Table3 depicts the adequacy of the distributional assumption built upon the empirical evidence presented in Silver et al. [15]. We conclude that the BS assumption for the LTD was more adequate in 81.60% of cases in terms of TCs and inventory levels; moreover, it prevails in circumstances of a higher LTD CV in comparison to gamma and normal assumptions, that is, a CV of 0.56 for the BS distribution against CVs of 0.44 and 0.34 for the gamma and normal distributions, respectively. It is worth mentioning that, due to the numeric integral for computing the BS CDF, which diverged in less than 1% of cases, the gamma distribution was a preferable assumption in terms of TCs and inventory levels in such cases.Table 3
Number and percentage of times that the indicated distribution yielded minimal TC for the simulations.
Distribution
Number
%
Average CV of the LTD
BS
8160
81.60
0.56
Normal
1780
17.80
0.34
Gamma∗
60
0.60
0.44
∗Number of cases where numeric integral for computing the BS CDF diverged.


## 3.1. Scenario of the Simulation Study

Assume BS, gamma, and normal distributions for the LTD. Then, fix values for the parameters of these distributions by considering values for means and SDs of DPUT and LT generated from uniform distributions. Now, generate holding, ordering, and shortage costs also from uniform distributions. This allows us to establish the expected TC to be minimized. Ten thousand (10000) different simulated scenarios of means and SDs for DPUT and LT as well as holding, ordering, and shortage costs are generated using anR package calledstats. The values of these uniformly distributed inventory indicators used to build the scenarios are presented in Table 1. They were chosen based on values proposed in selected papers focused on managerial and industrial applications compiled by Wanke [5]; see Table  2 in this reference for more details about values that are frequently used to generate simulation scenarios in inventory management problems. A discussion about this can be found in Wanke [51].Table 1
Range of the mentioned uniformly distributed indicator for simulations.
Indicator
Minimum
Maximum
DPUT mean (units/day)
80
120
DPUT SD (units/day)
3
30
LT mean (in days)
1
5
LT SD (in days)
0.50
2.00
Holding cost ($/unit/day)
0.00
0.68
Ordering cost ($/order)
17
60
Shortage cost ($/shortage)
0
100


## 3.2. Stochastic Programming

Once the values for the inventory indicators are defined and the distributional assumptions (BS, gamma, and normal) for the LTD established, stochastic programming is performed on the expected TC function given in (5); see Namit and Chen [3]. According to the values provided in Table 1, we assume values for means and SDs of the DPUT and LT, as well as for the holding, ordering, and shortage costs, to be uniformly distributed in the objective function corresponding to the expected TC given in (5). The decision variables of the programming are Q and k given in (1) and (4), respectively, whereas λ, σ
D, μ
L, σ
L, H, P, and S are given from (5). Therefore, our optimization problem of the expected TC can be visualized like a model of stochastic programming formulated as(14)
Minimize
Z
=
C
Q
,
k
Subject  to:
Q
>
0
,
k
>
0
,
λ
~
U
80,120
,
σ
D
~
U
3,30
,
μ
L
~
U
1,5
,
σ
L
~
U
0.50,20
,
H
~
U
0.00,0.68
,
P
~
U
17,60
,
S
~
U
0,100
.The problem of stochastic programming formulated in (14) is aimed (i) at identifying the most adequate inventory policy for each of the distributional assumptions and (ii) at minimizing the expected TC. The optimized problem can provide useful information for academics and practitioners on how these assumptions interact with product, operation, and demand characteristics. To solve the problem of stochastic programming formulated in (14), we use the DE algorithm detailed in the next section.


## 3.3. Differential Evolution

DE is a member of the family of genetic algorithms, which mimic the process of natural selection in an evolutionary manner; see Holland [52]. A genetic algorithm solves optimization problems with biology-inspired operators of crossover, mutation, and selection, generating successive populations of individuals (solutions or generations). Then, the DE algorithm optimizes problems by evolving a population of candidate solutions employing the mentioned operators. The DE algorithm uses floating-point techniques for obtaining the solutions and arithmetic operations in their mutation, in contrast to classic genetic algorithms. In addition, the DE algorithm finds the global optimum of the objective function, which is not required to be either continuous or differentiable; see Thangaraj et al. [46] and Mullen et al. [53].The DE algorithm has also been used to optimize problems that arise in inventory management, such as joint replenishment, replenishment coordination, and inventory location-allocation; see Qu et al. [54, 55] and Wang et al. [56]. In the present study, the problem of stochastic programming formulated in (14) is solved with the DE algorithm, which allows us to find the optimal values of Q and k that minimize the expected TC function C
(
·
) given in (5), for 10000 simulated scenarios.In what follows, we first discuss the DE algorithm used in our research work; see Storn and Price [47] and Price et al. [48] for more details about the algorithm. Then, we sketch the content of anR package calledDEoptim, which implements the DE algorithm and was first published on CRAN in 2005. Since becoming publicly available, it has been used by several authors to solve optimization problems arising in diverse domains. We refer interested readers to Ardia et al. [57] and Mullen et al. [53] for a detailed description of the package.LetN be the number of members y
∈
R
d (tuning parameter vector) in the population, where d denotes the dimension of the vector y, in our case given by y
=
(
Q
,
k
)
⊤
∈
R
2. The DE algorithm needs a starting population, which is obtained by sampling the objective function at multiple randomly chosen initial points (generation or solution zero (0)) for Q and k. Before the population is initialized, both lower and upper bounds for each tuning parameter must be specified. Parameter bounds establish the domain from which the N vectors at the generation 0 are chosen. To establish the generation 0, N guesses for the optimal value of y must be provided, using either random values within a range defined by the practitioner or values fixed by him/her. Each generation creates a new population from the current population members {
y
j
,
g
,
j
=
1
,
…
,
N
,
g
=
1,2
,
…
}, where j indexes the tuning parameter vector that make up the population and g indexes the generation. The new generation is obtained using differential mutation of the population members. In our research, N
=
10000 denotes the number of simulated scenarios and y denotes the tuning parameters for the lot size and the number of demand SDs. An initial mutant parameter vector {
v
j
,
g
,
j
=
1
,
…
,
N
,
g
=
1,2
,
…
} is created by choosing three members of the population y
j
1
,
g, y
j
2
,
g, and y
j
3
,
g, at random. Then, the elements of the initial mutant parameter vector v
j
,
g
=
(
v
j
i
,
g
) are generated by(15)
v
j
i
,
g
=
y
j
1
,
g
+
ζ
y
j
2
,
g
-
y
j
3
,
g
,
i
=
1
,
…
,
d
,where d
=
2 in our case and ζ is a positive scale factor whose effective value is usually less than one (usual default: ζ
=
0.8). After the first mutation operation, it continues until k mutations have been made or until certain crossover probability (CP) in [
0
,
1
] is less than u, where u is a random number from the uniform distribution in [
0
,
1
]. CP controls the fraction of the tuning parameter values that are copied from the mutant and approximates the probability that a parameter value is inherited from the mutant, since at least one mutation always occurs. Mutation is applied in this way to each member of the population.Calculations in our simulation study were performed with the aid of theDEoptim package, which consists of the core functionDEoptim() whose arguments are as follows:(i)
fn is the function to be minimized, which must have as its first argument the vector of real-valued parameters to optimize and return a scalar real result.
(ii)
lower,upper correspond to two vectors establishing scalar real lower and upper bounds on each tuning parameter to be optimized; the ith element of thelower andupper vectors corresponds to the ith parameter; the implementation searches the global optimum offn betweenlower andupper.
(iii)
… allows the user to pass additional arguments to the functionfn.
(iv)
control is a list whose default value is the return value ofDEoptim.control() but whose main elements are interpreted as follows:
(a)
VTR specifies the global minimum offn if it is known or if you wish to cease optimization after having reached a certain value (default =-Inf).
(b)
strategy defines the differential evolution strategy used in the optimization procedure, described in detail by Mullen et al. [53].
(c)
NPis the number of population members (d
e
f
a
u
l
t
=
10
×
d or 50), in our case denoted by N.
(d)
bs: ifbs isFALSE, then every mutant is tested against a member in the previous generation, and the best value survives into the next generation; ifbs isTRUE, then the old generation andNP mutants are sorted by their associated objective function values, and the bestNP vector proceeds into the next generation (default =FALSE).
(e)
itermax is the maximum number of iterations (i.e., population generations) to be allowed (default = 200).
(f)
CR is the CP from the interval [
0
,
1
] (default = 0.9).
(g)
F is the stepsize from interval [
0
,
2
], in our case denoted by ζ (default = 0.8).
(h)
trace,initialpop,storepopfrom,storepopfreq,checkWinner, andavWinner are other elements of the listcontrol, which are described in detail by Mullen et al. [53]; see also Price et al. [48].The return value of the functionDEoptim() is a member of the S3 classDEoptim. Members of this class have aplot and asummary that allow us to analyze the optimizer’s output. In our application of theDEoptim package, decision variables Q and k were lower-bounded at zero. As regards the list of tuning parameters used in our research, the default values for theDEoptim package were considered.


## 3.4. Numerical Results and Discussion

The underlying idea of performing a sensitivity analysis on the testing variables related to product, demand, and operational characteristics is to discriminate between groups where the three distributional assumptions led to minimal TCs. Table2 summarizes the numerical experiments conducted with theR software for 10000 simulated scenarios using the DE algorithm. Note that the BS distribution yields, on average, smaller TCs, EOQs, SSs, and, therefore, smaller inventory levels, in comparison to the normal and gamma assumptions.Table 2
Summary of the simulations for the mentioned distribution and management indicator.
Indicator
Distribution
Gamma
Normal
BS
Sum of TCs (in $)
205097.63
209175.49
202134.31
Sum of EOQs
767449.00
764816.00
762843.00
Sum of SSs
41166.66
53425.51
38103.58
Sum of average inventory level
424891.16
435833.51
419525.08Table3 depicts the adequacy of the distributional assumption built upon the empirical evidence presented in Silver et al. [15]. We conclude that the BS assumption for the LTD was more adequate in 81.60% of cases in terms of TCs and inventory levels; moreover, it prevails in circumstances of a higher LTD CV in comparison to gamma and normal assumptions, that is, a CV of 0.56 for the BS distribution against CVs of 0.44 and 0.34 for the gamma and normal distributions, respectively. It is worth mentioning that, due to the numeric integral for computing the BS CDF, which diverged in less than 1% of cases, the gamma distribution was a preferable assumption in terms of TCs and inventory levels in such cases.Table 3
Number and percentage of times that the indicated distribution yielded minimal TC for the simulations.
Distribution
Number
%
Average CV of the LTD
BS
8160
81.60
0.56
Normal
1780
17.80
0.34
Gamma∗
60
0.60
0.44
∗Number of cases where numeric integral for computing the BS CDF diverged.


## 4. Conclusions

This paper has assessed how different demand distributions during lead-time interact with relevant characteristics of the product. We have considered the choice of the optimal inventory policy in terms of total costs of the inventory management. Differently from most previous studies on the topic that exclusively explored the effects of one single distributional assumption, this paper also analyzed its adequacy in light of some inventory management key elements for decision-making, such as cost, demand, and lead-time. It was shown that the Birnbaum-Saunders distribution outperformed the normal and gamma assumptions with respect to demand uncertainty during the lead-time. The contributions of this paper are on both the academic and the practical sides. Departing from what was found in previous studies, the obtained results provide a guidance on the selection of the most appropriate distributional assumption for the demand during the lead-time. Specifically, while the Birnbaum-Saunders distribution is more adequate to handle demand during lead-time with high coefficients of variation, the gamma and the normal assumptions should be restricted to well-behaved patterns. This paper also can present a practical contribution by means of numerical analyses conducted with the aid of a computational code developed for such purposes in the R statistical software.With respect to future research endeavors, one possible extension of this work could be to analyze the results presented here in terms of inventory management in decentralized systems. The basic idea would be to assess, via numerical studies, the extent to which the decision to pool inventories is affected by different distributional assumptions and by their levels of skewness and kurtosis. In addition, in practical case studies, the usage of diagnostic statistical methods to detect influential data of demand or LT can be studied for this type of inventory models. In that case, demand distributions that provide parameter estimators to be robust to atypical data of demand must be considered; see, for example, Paula et al. [40].



---

*Source: 827246-2015-10-25.xml*