Stochastic Reorder Point-Lot Size (r,Q) Inventory Model under Maximum Entropy Principle

Abstract

This paper takes into account the continuous-review reorder point-lot size (r,Q) inventory model under stochastic demand, with the backorders-lost sales mixture. Moreover, to reflect the practical circumstance in which full information about the demand distribution lacks, we assume that only an estimate of the mean and of the variance is available. Contrarily to the typical approach in which the lead-time demand is supposed Gaussian or is obtained according to the so-called minimax procedure, we take a different perspective. That is, we adopt the maximum entropy principle to model the lead-time demand distribution. In particular, we consider the density that maximizes the entropy over all distributions with given mean and variance. With the aim of minimizing the expected total cost per time unit, we then propose an exact algorithm and a heuristic procedure. The heuristic method exploits an approximated expression of the total cost function achieved by means of an ad hoc first-order Taylor polynomial. We finally carry out numerical experiments with a twofold objective. On the one hand we examine the efficiency of the approximated solution procedure. On the other hand we investigate the performance of the maximum entropy principle in approximating the true lead-time demand distribution.

1. Introduction

When it is necessary to assign probabilities to mutually exclusive events in a sample space and there is not any prior knowledge about them, then we should assume that all these events have equal probability. This is called principle of insufficient reason or principle of indifference. It was originally stated by J. Bernoulli in 1713 and then endorsed by Laplace in 1814 [1].

In the case that we become aware of some information about the probability distribution of the outcomes, then we can adjust the assignment of probabilities accordingly. This is possible by means of the maximum entropy (MaxEnt) principle, which is a variational method of statistical inference originally proposed by Jaynes [2,3,4]. This principle basically works as follows: in the task concerned with determining a probability density (or mass) function subject to constraints, we should use the distribution satisfying those constraints that has the largest entropy. It is possible to prove that the principle of indifference derives from the MaxEnt principle in the case of a finite sample space when no constraint is imposed [5].

The MaxEnt principle has been widely adopted in a variety of fields. Taking into account some of the most recent works, we can cite the following contexts: urban planning [6]; queueing systems [7]; structural dynamics [8]; insurance [9]; computer vision and multimedia [10]; computer algebra [11]; geology [12]; biomechanics [13]; signal processing [14]. The literature gives some reviews concerning applications in ecology, finance, physics, chemistry, and biology [15,16,17,18].

In recent times, the MaxEnt principle has been introduced in the operations management area to approach the problem of evaluating the demand distribution or obtaining a demand forecast when only partial information about the demand are accessible [19,20,21]. Information are partial in the case the exact distribution is not available and only a certain number of observations is given.

With particular regard to the inventory management context, we would observe that, in practice, the decision-maker may actually know only the estimate of the mean and of the variance of the demand, but not the distribution type. In this circumstance, the traditional approach is to consider the demand within a given period as a normally distributed random variable [22]. This also follows from the assumption that individual demands are independent and identically distributed (i.i.d.) random variables, and then, according to the central limit theorem, the gaussianity of their sum can readily be deduced. However, this procedure is hardly valid in reality. In fact, single demands are generally not i.i.d. random variables [20]. In addition, one should also consider that the normal distribution is not recommended for items characterized by demand with large coefficient of variation [23]. A different approach to model the demand distribution is based on adopting the so-called minimax distribution-free procedure [24,25,26]. However, this gives an upper bound of the true cost, which may turn out in a not negligible error with respect to the optimal policy. To overcome the limitations characterizing the previous two methods, researches have recently introduced the MaxEnt principle in the task of optimizing an inventory system under partial information about the distribution of the quantity demanded [20]. We would however observe that, to the authors’ knowledge, the number of applications of the MaxEnt principle in the inventory management field are so far quite limited.

In the inventory management theory, the continuous-review reorder point-lot size (r,Q) model is a well-known policy. The system is continuously reviewed and, whenever the inventory position drops to or below r, an amount of Q units of goods is issued. The two control variables r and Q have different purposes. The replenishment quantity Q affects the trade-off between production or order frequency and inventory: a larger Q leads to few replenishments but higher average inventory levels; a smaller Q turns out to produce low average inventory but more replenishments. The reorder point r affects the stockout probability: a higher r gives larger inventory to assure a smaller stockout probability; a smaller r reduces inventory at the expense of a greater stockout probability.

As real inventory systems are typically subject to demand uncertainties, the (r,Q) policy under stochastic demand is therefore more practical. The literature proposes numerous works concerning stochastic (r,Q) models. In this regard, there exist studies involving Gaussian lead-time demand (e.g., [27,28,29,30]) or exploiting the minimax approach (e.g., [31,32,33]), as well as models that consider single- (e.g., [28,31,32]) or multi-echelon (e.g., [27,29,30]) systems. However, to the best of our knowledge, the MaxEnt principle has never been implemented into the (r,Q) policy.

We would point out that backorders-lost sales mixtures should not be neglected in a stochastic inventory model. In fact, these mixtures are generally adopted to model the different purchasing behaviors of customers when facing stockouts. Actually, some of them may wait until demand is satisfied (such demands are backordered); while others may be impatient (such demands are lost). Numerous studies involve this aspect (see, e.g., [34,35,36]).

Owing to the above observations, we consider the stochastic continuous-review reorder point-lot size (r,Q) model with backorders-lost sales mixtures. We derive the total cost function taking into account the MaxEnt principle. More precisely, we adopt the MaxEnt principle to model the lead-time demand distribution, given certain mean and variance. The purpose is to determine the replenishment policy that minimizes the expected total cost per time unit.

We present an exact algorithm and a heuristic solution procedure. The heuristic algorithm exploits an approximated expression of the total cost function achieved by means of an ad hoc first-order Taylor polynomial (i.e., a first-order truncation of the Taylor series expansion). We finally carry out numerical experiments with a twofold objective. On the one hand we examine the efficiency of the approximated solution procedure. On the other hand we investigate the performance of the MaxEnt principle in approximating the true lead-time demand distribution. In this regard, the MaxEnt principle is firstly compared with the Gaussian approximation and the minimax procedure. Then, a comparison with the approximation provided by the Weibull density is presented, as well. We would, in fact, compare the MaxEnt principle with the modeling capability provided by the Weibull density.

The remainder of the paper is as follows. Section 2 introduces notation, assumptions and optimization problem. In Section 3, we give the exact optimization procedure. In Section 4, we propose an approximated optimization approach. Section 5 presents numerical experiments. Finally, Section 6 deals with conclusions and further remarks.

2. Notation, Assumptions and Problem Definition

The main notation adopted is the following:

Decision variables:
Q	Order quantity (quantity units).
r	Reorder point (quantity units).
Parameters:
h	Unit holding cost rate (monetary unit/quantity unit/time unit).
A	Fixed ordering cost per order (monetary unit/order).
π	Fixed penalty cost per unit short (monetary unit/quantity unit).
π0	Marginal profit per unit (monetary unit/quantity unit).
L	Replenishment lead time (time unit).
$\tilde{\mathsf{\mu }}$	Average demand rate (quantity unit/time unit).
$\tilde{\mathsf{\sigma }}$	Standard deviation of the demand rate (quantity unit/time unit).
β	Fraction of shortage (i.e., demand during the stockout period) that will be lost.
Random variables:
X	Lead-time demand, i.e., quantity demanded during the lead time.
Functions and operators:
$\overline{f}(\cdot )$	Probability density function (p.d.f.) of the lead-time demand.
x+	Maximum between 0 and x, i.e., $x^{+}\equiv \max \left\{0,x\right\}$.

The main assumptions are the following:

• Inventory is continuously reviewed. An order of size Q is issued when the on-hand inventory reaches the reorder point r.
• The reorder point r is positive.
• The distribution of X is unknown/unspecified, but only its mean $\mathsf{\mu }\equiv \tilde{\mathsf{\mu }}L$ and variance ${\mathsf{\sigma }}^2\equiv {\tilde{\mathsf{\sigma }}}^{2}L$ can be evaluated.
• The random variable X is continuous and nonnegative. That is, the lead-time demand can take any value in ${ℝ}_0^{+}$.
• Shortages are allowed and partially backordered with ratio 1 − β. The fraction of shortage with ratio β is lost.
• The time horizon is infinite.

Under our assumptions, the expected total cost per time unit is given by:

$$\overline{C}\left(Q,r\right)=A\frac{\tilde{\mathsf{\mu }}}{Q}+h\left[\frac{Q}{2}+r-\tilde{\mathsf{\mu }}L+\mathsf{\beta }\overline{B}\left(r\right)\right]+\overline{\mathsf{\pi }}\frac{\tilde{\mathsf{\mu }}}{Q}\overline{B}\left(r\right)$$

(1)

where $\overline{\mathsf{\pi }}\equiv \mathsf{\pi }+{\mathsf{\pi }}_0\mathsf{\beta }$, and

$$\overline{B}\left(r\right)=E\left[{\left(X-r\right)}^{+}\right]={\displaystyle \underset{r}{\overset{+\infty }{\mathrm{\int }}}\left(x-r\right)\overline{f}\left(x\right)dx}$$

(2)

is the expected shortage at the end of the cycle. In Equation (2), $E[\cdot ]$ represents the mathematical expectation. Cost function (Equation (1)) consists of the ordering cost, the inventory holding cost, and the shortage cost [37]. Note that Equation (2) can be determined once the density $\overline{f}$ has been specified.

Optimizing replenishments in the inventory system under consideration means solving the following problem:

$$\begin{array}{lll}\left(\text{P1}\right)\hfill & \underset{\left(Q,r\right)}{\min }\hfill & \overline{C}\left(Q,r\right),\hfill \\ \hfill & \text{s.t.}\hfill & Q\in {ℝ}_{+},\hfill \\ \hfill & \hfill & r\in {ℝ}_{+},\hfill \end{array}$$

where $\overline{C}\left(Q,r\right)$ is given by Equation (1). We remind that we are in the case where the actual density $\overline{f}$ of X is unknown/unspecified, but only the mean μ and the variance ${\mathsf{\sigma }}^2$ of X can be assessed. In this circumstance, we cannot solve problem (P1) directly, as the quantity $\overline{B}\left(r\right)$ cannot explicitly be calculated (since $\overline{f}$ is not given). Therefore, the exact solution that optimizes replenishments in the considered system cannot be determined. Consequently, the problem can only be approached adopting a suitable approximation.

In this regard, the procedure that is typically adopted in the literature (see, e.g., [27,28,29,30,31,32,33]) consists in assuming that the density of X is Gaussian or is derived according to the minimax method. We take here a different perspective, which is based on adopting the MaxEnt principle. That is, in Equation (2) we replace $\overline{f}$ with the density f maximizing the entropy over all densities defined on ${ℝ}_0^{+}$ with mean μ and variance ${\mathsf{\sigma }}^2$. This maximization problem can be approached according to the following proposition [5]:

Proposition 1. 

(Maximum entropy distribution) Let us consider the problem of maximizing the entropy $h\left(g\right)=-{\displaystyle {\mathrm{\int }}_{\mathbb{S}}g\left(x\right)\log g\left(x\right)dx}$ over all p.d.f.s g with support S satisfying the following constraints:

• $g\left(x\right)\ge 0$, with equality outside the support S,
• ${\displaystyle {\mathrm{\int }}_{S}g\left(x\right)dx}=1$,
• ${\displaystyle {\mathrm{\int }}_S{p}_i\left(x\right)g\left(x\right)dx}={\alpha }_i$, for $i=1,2,\mathrm{...},m$,

where, for $i=1,2,\mathrm{...},m$, P_i is a (measurable) function and ${\alpha }_i$ is a real number. Let $g^{*}\left(x\right)=\exp \left\{{\lambda }_0+{\displaystyle {\sum }_{i=1}^m{\lambda }_i{p}_i\left(x\right)}\right\}$, $x\in S$, where ${\lambda }_0,\mathrm{...},{\lambda }_m$ are chosen so that $g^{*}$ satisfies the previous constraints. Then $g^{*}$ uniquely maximizes h(g) over all probability density functions g satisfying the previous constraints.

Thanks to Proposition 1, the density f maximizing the entropy over all densities defined on ${ℝ}_0^{+}$ with mean μ and variance ${\mathsf{\sigma }}^2$ is given by

$$f\left(x\right)=\exp \left\{a{x}^2+bx+c\right\},x\in {ℝ}_0^{+}$$

(3)

where the quantities a, b, and c can readily be obtained by solving the following system of equations:

$$\{\begin{array}{l}{\displaystyle \underset{0}{\overset{+\infty }{\mathrm{\int }}}\exp \left\{a{x}^2+bx+c\right\}dx}=1,\hfill \\ {\displaystyle \underset{0}{\overset{+\infty }{\mathrm{\int }}}x\exp \left\{a{x}^2+bx+c\right\}dx}=\mathsf{\mu },\hfill \\ {\displaystyle \underset{0}{\overset{+\infty }{\mathrm{\int }}}{x}^2\exp \left\{a{x}^2+bx+c\right\}dx}={\mathsf{\sigma }}^2+{\mathsf{\mu }}^2.\hfill \end{array}$$

In what follows, the density f will be also referred to as the maximum entropy (or MaxEnt) density.

If we replace $\overline{f}$ with f in Equation (2), the expected shortage at the end of the cycle becomes:

$$B\left(r\right)={\displaystyle \underset{r}{\overset{+\infty }{\mathrm{\int }}}\left(x-r\right)f\left(x\right)dx}$$

(4)

and substituting $\overline{B}\left(r\right)$ with $B\left(r\right)$ in Equation (1) we get:

$$C\left(Q,r\right)=A\frac{\tilde{\mathsf{\mu }}}{Q}+h\left[\frac{Q}{2}+r-\tilde{\mathsf{\mu }}L+\mathsf{\beta }B\left(r\right)\right]+\overline{\mathsf{\pi }}\frac{\tilde{\mathsf{\mu }}}{Q}B\left(r\right)$$

(5)

Ultimately, instead of solving problem (P1) directly (that is not possible, as the quantity $\overline{B}\left(r\right)$ in Equation (1) cannot explicitly be determined) we turn to approach the following:

$$\begin{array}{lll}\left(\text{P2}\right)\hfill & \underset{\left(Q,r\right)}{\min }\hfill & C\left(Q,r\right),\hfill \\ \hfill & \text{s.t.}\hfill & Q\in {ℝ}_{+},\hfill \\ \hfill & \hfill & r\in {ℝ}_{+},\hfill \end{array}$$

where $C\left(Q,r\right)$ is given by Equation (5). Evidently, solving problem (P2) is not equivalent to determining the optimal solution to problem (P1). That is, the solution $\left(Q^{*},r^{*}\right)$ to problem (P2) is not optimal to problem (P1). In fact, the maximum entropy density f is used to “approximate” the true density $\overline{f}$ of the lead-time demand, which is unknown by assumption. Therefore, $\left(Q^{*},r^{*}\right)$ may be considered a “heuristic” (i.e., not optimal) solution to problem (P1). If the decision-maker had full information about $\overline{f}$, then it would clearly be preferable to approach problem (P1) directly.

3. Exact Procedure to Solve Problem (P2)

We can first note that the integrals in the above system of equations can converge if and only if $a<0.$ Under this condition, it is possible to check that the complementary c.d.f. F⁰ associated to the density f given by Equation (3) is expressed as follows:

$$F^0\left(x\right)={\displaystyle {\mathrm{\int }}_x^{+\infty }f\left(t\right)dt}=\frac{1}{2}\sqrt{\frac{\mathsf{\pi }}{-a}}{e}^{-\frac{b^2}{4a}+c}\left[1-\text{erf}\left(-\frac{2ax+b}{2\sqrt{-a}}\right)\right],$$

(6)

where $\text{erf}(\cdot )$ is the Error function [38]. The related c.d.f. F is evidently given by $F\left(x\right)=1-F^0\left(x\right)$. Moreover, the quantity B in Equation (5), i.e., the expected shortage at the end of the cycle, can be determined with the following relation:

$$\begin{array}{c}B\left(r\right)={\displaystyle {\mathrm{\int }}_r^{+\infty }\left(x-r\right)f\left(x\right)dx}={\displaystyle {\mathrm{\int }}_r^{+\infty }{F}^0\left(x\right)dx}\\ =-\frac{1}{2a}\left\{e^{a{x}^2+bx+c}+\frac{1}{2}\sqrt{\frac{\mathsf{\pi }}{-a}}\left(2ar+b\right)e^{-\frac{b^2}{4a}+c}\left[1-\text{erf}\left(-\frac{2ar+b}{2\sqrt{-a}}\right)\right]\right\}\\ =-\frac{1}{2a}\left[f\left(r\right)+\left(2ar+b\right)F^0\left(r\right)\right].\end{array}$$

(7)

The equality ${\displaystyle {\mathrm{\int }}_r^{+\infty }\left(x-r\right)f\left(x\right)dx}={\displaystyle {\mathrm{\int }}_r^{+\infty }{F}^0\left(x\right)dx}$ in Equation (7) is known in literature [39]. Equations (6) and (7) can be obtained exploiting the numerous properties of the Error function given in [40].

From (7), we have that $\frac{d^2}{d{r}^2}B\left(r\right)=f\left(r\right)$. Hence, we can affirm that $B\left(r\right)$ is a convex function of r. Moreover, observing that $E\left[{\left(X-r\right)}^{+}\right]/Q$ is a convex function of $\left(Q,r\right)$, for $Q>0$ and $r>0$, for any p.d.f. of X [41], the convexity of $C\left(Q,r\right)$ in $\left(Q,r\right)$, for $Q>0$ and $r>0$, can readily be deduced. Therefore, the solution to problem (P2) can be obtained by solving the First-Order Conditions of optimality in $\left(Q,r\right)$.

If we impose the First-Order Condition in r, we obtain:

$$F^0\left(r\left(Q\right)\right)=\frac{h}{h\mathsf{\beta }+\overline{\mathsf{\pi }}\frac{\tilde{\mathsf{\mu }}}{Q}}$$

(8)

where $r\left(Q\right)$, which is the optimal r value for any Q, must satisfy Equation (8). In other terms, the optimal r value for any Q can be found from Equation (8) by inverting F (we remind that $F^0\left(x\right)=1-F\left(x\right)$):

$$r\left(Q\right)=F^{-1}\left(1-\frac{h}{h\mathsf{\beta }+\overline{\mathsf{\pi }}\frac{\tilde{\mathsf{\mu }}}{Q}}\right)$$

(9)

where $F^{-1}$ is the inverse function of F. We can further observe that (5) can conveniently be rewritten as follows:

$$C\left(Q,r\right)=A\frac{\tilde{\mathsf{\mu }}}{Q}+h\left(\frac{Q}{2}+r-\tilde{\mathsf{\mu }}L\right)+B\left(r\right)\left(h\mathsf{\beta }+\overline{\mathsf{\pi }}\frac{\tilde{\mathsf{\mu }}}{Q}\right)$$

(10)

Moreover, inserting Equations (7)–(9) into Equation (10) we obtain, after some algebraic manipulations, the following expression:

$$C\left(Q\right)=A\frac{\tilde{\mathsf{\mu }}}{Q}+h\left(\frac{Q}{2}-\tilde{\mathsf{\mu }}L\right)-\frac{1}{2a}\left[\left(h\mathsf{\beta }+\overline{\mathsf{\pi }}\frac{\tilde{\mathsf{\mu }}}{Q}\right)f\left(r\left(Q\right)\right)+hb\right]$$

(11)

In conclusion, the Q-component $Q^{*}$ of the solution to problem (P2) can be found by minimizing Equation (11) in Q. Note that this can only be achieved by means of a numerical technique. In fact, the First-Order Condition in Q imposed to Equation (11) (i.e., the equation $dC\left(Q\right)/dQ=0$) cannot be solved in Q in closed form. In other words, a closed-form formula that solves the equation $dC\left(Q\right)/dQ=0$ in Q does not exist. To minimize Equation (11) in Q, a standard constrained nonlinear minimization algorithm can be used, e.g., the interior-point algorithm. This algorithm can, in fact, be implemented numerically. The r-component $r^{*}$ of the solution to problem (P2) can then be determined inserting $Q^{*}$ into Equation (9).

4. An Approximated Procedure to Solve Problem (P2)

We would observe that optimization approaches that are effective but also difficult to be implemented in practice, e.g., because of large computational time and/or lack of a simple solution procedure (which may consist in a simple formula), may have little relevance under a practical point of view. For example, when a closed-form expression to the First-Order Conditions of optimality lacks, the optimal solution may be obtained with a not negligible effort, which may make unpractical the model itself. This may be the case of a large retailer that typically needs to manage thousands of items, and the control variables are often required to be recalculated frequently. Therefore, the development of efficient and practically applicable approximated optimization procedures is strongly encouraged in these circumstances [42,43,44].

Owing to the above observations, in this section we propose an efficient approximated procedure to approach problem (P2). In fact, we remind that the optimal solution can only be achieved by means of a numerical technique, and this may limit the application of the model in practice. The presented approximation method permits us to obtain a simple formula, which resembles the classic EOQ expression.

The proposed near-optimal procedure is based on an ad hoc approximation of part of cost function (11). In particular, we replace $f\left(r\left(Q\right)\right)$ with a first-order truncation of its Taylor series expansion in Q centered in

$$\overline{Q}\equiv \sqrt{\frac{2\tilde{\mathsf{\mu }}A}{h}}$$

which is the well-known optimal solution in the EOQ model. We would observe that we do not investigate the convergence properties of the Taylor series expansion of $f\left(r\left(Q\right)\right)$. This for two main reasons: (i) this task is particularly hard to be analytically accomplished (in this regard, one should also consider that the model involves several parameters, and the convergence is plausibly strongly dependent on their specific values); and (ii) we rely on the fact that, implementing an approximation around the optimal solution in deterministic conditions, the true cost function is “close” to the approximated expression (in particular around the minimum), with an error that is intuitively even smaller as the variability in the system decreases. Although this last point is actually a heuristic reasoning, it is supported by experimental evidences. The reader is referred to the numerical study section (Section 5), where tests show the efficiency of our approximation method. We would finally remind that a similar approach has successfully been implemented in previous researches [44,45].

Taking the first-order derivative in Q of both side of Equation (8), we have:

$$\begin{array}{ll}\hfill & \frac{d}{dQ}{F}^0\left(r\left(Q\right)\right)=\frac{d}{dQ}\left(\frac{h}{h\mathsf{\beta }+\overline{\mathsf{\pi }}\frac{\tilde{\mathsf{\mu }}}{Q}}\right)\hfill \\ \Rightarrow \hfill & \frac{d{F}^0\left(r\left(Q\right)\right)}{dr\left(Q\right)}\frac{dr\left(Q\right)}{dQ}=\frac{h\overline{\mathsf{\pi }}\tilde{\mathsf{\mu }}}{{\left(\overline{\mathsf{\pi }}\tilde{\mathsf{\mu }}+h\mathsf{\beta }Q\right)}^2}\hfill \\ \Rightarrow \hfill & -f\left(r\left(Q\right)\right)\frac{dr\left(Q\right)}{dQ}=\frac{h\overline{\mathsf{\pi }}\tilde{\mathsf{\mu }}}{{\left(\overline{\mathsf{\pi }}\tilde{\mathsf{\mu }}+h\mathsf{\beta }Q\right)}^2}\hfill \\ \Rightarrow \hfill & \frac{dr\left(Q\right)}{dQ}=-\frac{1}{f\left(r\left(Q\right)\right)}\frac{h\overline{\mathsf{\pi }}\tilde{\mathsf{\mu }}}{{\left(\overline{\mathsf{\pi }}\tilde{\mathsf{\mu }}+h\mathsf{\beta }Q\right)}^2}.\hfill \end{array}$$

We can then note that:

$$\begin{array}{ll}\hfill & \frac{d}{dQ}f\left(r\left(Q\right)\right)=\frac{df\left(r\left(Q\right)\right)}{dr\left(Q\right)}\frac{dr\left(Q\right)}{dQ}\hfill \\ =\hfill & \left(2ar\left(Q\right)+b\right)f\left(r\left(Q\right)\right)\left[-\frac{1}{f\left(r\left(Q\right)\right)}\frac{h\overline{\mathsf{\pi }}\tilde{\mathsf{\mu }}}{{\left(\overline{\mathsf{\pi }}\tilde{\mathsf{\mu }}+h\mathsf{\beta }Q\right)}^2}\right]\hfill \\ =\hfill & -\left(2ar\left(Q\right)+b\right)\frac{h\overline{\mathsf{\pi }}\tilde{\mathsf{\mu }}}{{\left(\overline{\mathsf{\pi }}\tilde{\mathsf{\mu }}+h\mathsf{\beta }Q\right)}^2}.\hfill \end{array}$$

Therefore, with reference to a neighborhood of $\overline{Q}$, we can write:

$$f\left(r\left(Q\right)\right)\approx p_0+p_1\left(Q-\overline{Q}\right)$$

(12)

where:

$$\begin{array}{l}{p}_0\equiv f\left(r\left(\overline{Q}\right)\right),\hfill \\ p_1\equiv -\left(2ar\left(\overline{Q}\right)+b\right)\frac{h\overline{\mathsf{\pi }}\tilde{\mathsf{\mu }}}{{\left(\overline{\mathsf{\pi }}\tilde{\mathsf{\mu }}+h\mathsf{\beta }\overline{Q}\right)}^2}.\hfill \end{array}$$

With some algebraic manipulations, Equation (12) can conveniently be rewritten as follows:

$$f\left(r\left(Q\right)\right)\approx s_0+s_{1}Q,$$

(13)

where:

$$\begin{array}{l}{s}_0\equiv p_0-p_1\overline{Q},\hfill \\ s_1\equiv p_1.\hfill \end{array}$$

Finally, inserting Equation (13) into Equation (11), we can approximate $C\left(Q\right)$ in a neighborhood of $\overline{Q}$ with the following expression:

$$\widehat{C}\left(Q\right)=\frac{u}{Q}+vQ+y,$$

(14)

where:

$$\begin{array}{l}u\equiv A\tilde{\mathsf{\mu }}-\frac{1}{2a}{s}_0\overline{\mathsf{\pi }}\tilde{\mathsf{\mu }},\hfill \\ v\equiv \frac{h}{2}-\frac{1}{2a}{s}_{1}h\mathsf{\beta },\hfill \\ y\equiv -\frac{1}{2a}\left(s_{0}h\mathsf{\beta }+s_1\overline{\mathsf{\pi }}\tilde{\mathsf{\mu }}+hb\right)-h\tilde{\mathsf{\mu }}L.\hfill \end{array}$$

We can observe that $\widehat{C}\left(Q\right)$ resembles the deterministic cost structure plus a constant term in Q. The near-optimal Q-component $\widehat{Q}$ of the solution to problem (P2) can therefore be found solving the equation $d\widehat{C}\left(Q\right)/dQ=0$, which is equivalent to

$$N\left(Q\right)\equiv v{Q}^2-u=0$$

Under the assumption that u and v are positive quantities, $N\left(Q\right)$ admits a unique positive root, which evidently coincides with $\widehat{Q}$. In the numerical section, we will show that this fundamental assumption is licit for a reasonably wide range of parameter values. Therefore, $\widehat{Q}$ is given by

$$\widehat{Q}=\sqrt{\frac{u}{v}}$$

(15)

Once $\widehat{Q}$ has been obtained, the corresponding near-optimal value $\widehat{r}$ of r can be found inserting $\widehat{Q}$ into Equation (9).

We would finally remark that the approximation procedure has allowed the derivation of a near-optimal solution (in Q) in closed form. Moreover, this expression is particularly simple, as it resembles the classic EOQ formula. Under a practical point of view, it is evidently simpler implementing Equation (15) rather than minimizing Equation (11) with a numerical procedure. The approximated solution method may therefore foster the implementation in practice of the model proposed in this paper. In the next section, we will show that the solution found with Equation (15) is efficient for a wide range of parameter values.

5. Numerical Experiments

This section presents numerical experiments whose purpose is that of examining two different questions. First, we numerically evaluate the efficiency of the approximated optimization procedure given in Section 4. The efficiency is assessed in terms of both computational effort required and error achieved. Then, we investigate the performance of the MaxEnt principle in approximating the true lead-time demand distribution under limited information, i.e., when only an estimate of the mean and of the variance of the lead-time demand is available. In the first part of this analysis, we compare the MaxEnt principle with two alternative procedures, i.e., the Gaussian approximation and the minimax approach, taking into account several classes of demand distributions. In the second, and last, part, we present a comparison between the MaxEnt principle and the approximation provided by the Weibull density.

These experiments were performed on a PC with an Intel^® Core^™ i7 processor at 2.4 GHz and with 16 GB of RAM. Moreover, MATLAB^® R2013b was used as computing environment.

5.1. Efficiency of the Approximated Optimization Method

Let us consider the following quantity:

$$\text{APE}\equiv \frac{\left|C\left(Q^{*},r^{*}\right)-C\left(\widehat{Q},\widehat{r}\right)\right|}{C\left(Q^{*},r^{*}\right)}\times 100$$

(16)

which is the Absolute Percentage Error (APE). In Equation (16), cost function C is given by Equation (5), $\left(Q^{*},r^{*}\right)$ is the minimum-cost solution (obtained according to the procedure given in Section 3), and $\left(\widehat{Q},\widehat{r}\right)$ is the near-optimal solution (obtained according to the procedure given in Section 4).

The error is evaluated by means of Equation (16) for different combinations of parameter values. Parameter values are randomly drawn within the intervals shown in Table 1. Although the ranges chosen are purely indicative, values are similar to those typically adopted in the inventory management literature (see, e.g., [24,37,44]). Note that in Table 1 we consider $c_v\equiv \tilde{\mathsf{\sigma }}/\tilde{\mathsf{\mu }}$ instead of $\tilde{\mathsf{\sigma }}$. Once $\tilde{\mathsf{\mu }}$ and c_v are fixed, the corresponding value of $\tilde{\mathsf{\sigma }}$ is therefore given by $\tilde{\mathsf{\sigma }}=c_v\cdot \tilde{\mathsf{\mu }}$.

Table 1. Intervals where parameters take values.

Parameters	Intervals	Units of Measurement
$\tilde{\mathsf{\mu }}$	$\left[100,1000\right]$	(units/year)
cv	$\left[0.05,0.80\right]$	-
L	$\left[0.03,0.17\right]$	(years)
A	$\left[100,250\right]$	($/order)
h	$\left[1,25\right]$	($/unit/year)
π	$\left[20,70\right]$	($/unit)
π0	$\left[80,150\right]$	($/unit)
β	$\left[0.1,0.9\right]$	-

Table 1. Intervals where parameters take values.

Parameters	Intervals	Units of Measurement
$\tilde{\mathsf{\mu }}$	$\left[100,1000\right]$	(units/year)
cv	$\left[0.05,0.80\right]$	-
L	$\left[0.03,0.17\right]$	(years)
A	$\left[100,250\right]$	($/order)
h	$\left[1,25\right]$	($/unit/year)
π	$\left[20,70\right]$	($/unit)
π0	$\left[80,150\right]$	($/unit)
β	$\left[0.1,0.9\right]$	-

Results of the error analysis are shown in Table 2. We can note that the approximated optimization method mainly achieves a very small APE (i.e., APE < 1%). This happens for a reasonable wide range of parameter values. In some experiments, it is less than 0.01%, and it can therefore be considered negligible. In a few cases, the APE reaches values greater than 1%, but smaller than 2%. In general, we can thus affirm that the error achieved is sufficiently small to consider reasonably good the approximation made. It is worth noting that, in the cases where the error is large, i.e., APE > 1%, c_v is high. That is, it seems that a greater variability in the system leads to higher error. This is expected, as the approximation is based on a Taylor series expansion around the optimal solution in deterministic conditions.

Table 2. Results of the error analysis.

Test Number	$\tilde{\mathsf{\mu }}$	cv	L	A	h	π	π0	β	APE
1	142	0.31	0.13	219	14	54	143	0.14	0.04%
2	276	0.59	0.13	232	15	23	145	0.74	<0.01%
3	987	0.59	0.15	165	12	48	99	0.70	0.92%
4	377	0.15	0.10	154	20	59	127	0.21	0.14%
5	371	0.75	0.17	143	21	65	122	0.81	1.93%
6	756	0.48	0.03	167	17	46	106	0.85	0.08%
7	435	0.50	0.15	240	17	30	126	0.16	0.61%
8	941	0.67	0.10	214	11	69	150	0.79	0.32%
9	322	0.69	0.15	237	14	50	90	0.82	0.38%
10	910	0.62	0.15	143	17	53	88	0.43	1.88%
11	355	0.72	0.15	158	13	55	139	0.59	0.54%
12	511	0.59	0.15	208	1	54	111	0.45	0.01%
13	392	0.23	0.08	156	14	48	108	0.42	<0.01%
14	956	0.59	0.09	255	4	23	85	0.23	0.07%
15	110	0.45	0.04	122	16	63	149	0.56	0.01%
16	564	0.30	0.09	174	2	65	84	0.45	<0.01%
17	652	0.66	0.15	240	5	33	143	0.57	0.15%
18	838	0.45	0.06	168	11	69	124	0.66	<0.01%
19	565	0.47	0.05	184	18	41	139	0.69	0.10%
20	448	0.63	0.13	164	18	68	135	0.66	0.67%
21	632	0.39	0.04	134	21	20	141	0.16	0.13%
22	296	0.48	0.05	201	15	22	84	0.22	0.05%
23	849	0.51	0.10	230	3	66	87	0.51	<0.01%
24	104	0.63	0.15	238	25	45	99	0.18	0.34%
25	787	0.11	0.12	178	5	67	121	0.45	<0.01%
26	507	0.68	0.10	183	18	38	96	0.56	0.01%
27	201	0.38	0.07	160	21	40	107	0.39	0.07%
28	178	0.37	0.07	144	11	26	115	0.67	0.02%
29	166	0.35	0.03	133	1	29	90	0.31	<0.01%
30	639	0.73	0.16	133	13	39	117	0.31	1.52%
31	256	0.07	0.16	165	25	58	80	0.64	<0.01%
32	597	0.21	0.14	134	10	65	140	0.42	<0.01%
33	920	0.73	0.11	150	22	42	144	0.13	1.46%
34	261	0.30	0.06	148	11	47	83	0.54	0.01%
35	319	0.17	0.16	241	21	57	92	0.39	<0.01%
36	663	0.46	0.09	143	13	58	134	0.56	0.23%
37	211	0.43	0.08	113	4	30	127	0.45	0.01%
38	108	0.45	0.07	242	23	40	81	0.64	<0.01%
39	151	0.39	0.11	203	18	53	131	0.40	0.06%
40	151	0.78	0.07	189	25	29	93	0.37	0.33%
41	346	0.16	0.09	156	4	42	86	0.59	<0.01%
42	811	0.23	0.09	185	2	45	125	0.28	<0.01%
43	862	0.43	0.07	212	6	68	124	0.58	0.03%
44	329	0.69	0.16	205	19	31	120	0.75	0.84%
45	181	0.29	0.10	109	19	48	117	0.76	<0.01%

Table 2. Results of the error analysis.

Test Number	$\tilde{\mathsf{\mu }}$	cv	L	A	h	π	π0	β	APE
1	142	0.31	0.13	219	14	54	143	0.14	0.04%
2	276	0.59	0.13	232	15	23	145	0.74	<0.01%
3	987	0.59	0.15	165	12	48	99	0.70	0.92%
4	377	0.15	0.10	154	20	59	127	0.21	0.14%
5	371	0.75	0.17	143	21	65	122	0.81	1.93%
6	756	0.48	0.03	167	17	46	106	0.85	0.08%
7	435	0.50	0.15	240	17	30	126	0.16	0.61%
8	941	0.67	0.10	214	11	69	150	0.79	0.32%
9	322	0.69	0.15	237	14	50	90	0.82	0.38%
10	910	0.62	0.15	143	17	53	88	0.43	1.88%
11	355	0.72	0.15	158	13	55	139	0.59	0.54%
12	511	0.59	0.15	208	1	54	111	0.45	0.01%
13	392	0.23	0.08	156	14	48	108	0.42	<0.01%
14	956	0.59	0.09	255	4	23	85	0.23	0.07%
15	110	0.45	0.04	122	16	63	149	0.56	0.01%
16	564	0.30	0.09	174	2	65	84	0.45	<0.01%
17	652	0.66	0.15	240	5	33	143	0.57	0.15%
18	838	0.45	0.06	168	11	69	124	0.66	<0.01%
19	565	0.47	0.05	184	18	41	139	0.69	0.10%
20	448	0.63	0.13	164	18	68	135	0.66	0.67%
21	632	0.39	0.04	134	21	20	141	0.16	0.13%
22	296	0.48	0.05	201	15	22	84	0.22	0.05%
23	849	0.51	0.10	230	3	66	87	0.51	<0.01%
24	104	0.63	0.15	238	25	45	99	0.18	0.34%
25	787	0.11	0.12	178	5	67	121	0.45	<0.01%
26	507	0.68	0.10	183	18	38	96	0.56	0.01%
27	201	0.38	0.07	160	21	40	107	0.39	0.07%
28	178	0.37	0.07	144	11	26	115	0.67	0.02%
29	166	0.35	0.03	133	1	29	90	0.31	<0.01%
30	639	0.73	0.16	133	13	39	117	0.31	1.52%
31	256	0.07	0.16	165	25	58	80	0.64	<0.01%
32	597	0.21	0.14	134	10	65	140	0.42	<0.01%
33	920	0.73	0.11	150	22	42	144	0.13	1.46%
34	261	0.30	0.06	148	11	47	83	0.54	0.01%
35	319	0.17	0.16	241	21	57	92	0.39	<0.01%
36	663	0.46	0.09	143	13	58	134	0.56	0.23%
37	211	0.43	0.08	113	4	30	127	0.45	0.01%
38	108	0.45	0.07	242	23	40	81	0.64	<0.01%
39	151	0.39	0.11	203	18	53	131	0.40	0.06%
40	151	0.78	0.07	189	25	29	93	0.37	0.33%
41	346	0.16	0.09	156	4	42	86	0.59	<0.01%
42	811	0.23	0.09	185	2	45	125	0.28	<0.01%
43	862	0.43	0.07	212	6	68	124	0.58	0.03%
44	329	0.69	0.16	205	19	31	120	0.75	0.84%
45	181	0.29	0.10	109	19	48	117	0.76	<0.01%

With regard to the error analysis, we have carried out an additional investigation. APE has been evaluated over 1000 randomly generated problems, with parameter values drawn within the intervals in Table 1. Results are as follows:

• In about 42.49% of cases, APE < 0.01%;
• In about 51.16% of cases, 0.01% < APE ≤ 1%;
• In about 4.44% of cases, 1% < APE ≤ 2%;
• In about 1.92% of cases, APE > 2%;
• The maximum value achieved is 4.8%.

We can observe that, in more than 93% of cases, APE is smaller than or equal to 1%. While, in more than 98% of cases, we have that APE ≤ 2%.

We now evaluate the computational effort required by the approximated method and by the exact algorithm to solve problem (P2). To this aim, we consider the time required to solve 1000 random problems. In fact, although the time difference on a single problem is practically negligible, the ratio of the computational times may become significant over several problems. In each problem, parameter values are randomly drawn within the intervals shown in Table 1. Both algorithms are tested on the same problems. That is, the comparison is made in terms of time needed to solve the same batch of 1000 random problems, where, in each problem, parameter values are (randomly) drawn within the intervals in Table 1. Results are as follows: the exact algorithm needed 99.42 s, while the approximated solution method spent 6.84 s. That is, over identical 1000 random problems, the percentage of computational time reduction achieved by the approximated solver is more than 93%.

In conclusion, we can assert that the approximated solution method is efficient, in terms of both error achieved and computational effort required. It seems therefore promising for a practical application. We would finally observe that, in every test, the assumption that u and v are positive has been satisfied.

5.2. Comparative Analysis

In this subsection, we investigate the performance of the MaxEnt principle in approximating the true lead-time demand distribution, when only an estimate of the mean and of the variance is given. In the first part of this subsection, this experiment is made taking into account several classes of demand distributions: lognormal, gamma, and Weibull. The MaxEnt principle is compared with two alternative procedures: the Gaussian approximation and the minimax approach. In the second, and last, part of this subsection, the MaxEnt principle is compared with the approximation provided by the Weibull density. In this second analysis, the true density of the lead-time demand is assumed to be a mixture of lognormal distributions.

Let us begin with the first study. Table 3 shows the parameters whose value is kept fixed. These parameter values have been randomly drawn within the ranges in Table 1. The other parameters, i.e., c_v, L and h, take several different values. This is to study the sensitivity of the response with respect to variations in the value of these parameters, which significantly affect the optimal replenishment policy.

Table 3. Parameters whose value is kept fixed.

Parameters	Values	Units of Measurement
$\tilde{\mathsf{\mu }}$	834	(units/year)
A	237	($/order)
π	24	($/unit)
π0	99	($/unit)
β	0.54	-

Table 3. Parameters whose value is kept fixed.

Parameters	Values	Units of Measurement
$\tilde{\mathsf{\mu }}$	834	(units/year)
A	237	($/order)
π	24	($/unit)
π0	99	($/unit)
β	0.54	-

Each single experiment, which is defined for a given set of parameter values, is carried out as follows.

The values of $\tilde{\mathsf{\mu }}$, c_v, and L are used to sample ten observations from the true distribution of the lead-time demand. The parameters of the true lead-time demand distribution have clearly to be determined by imposing that the mean and the standard deviation are $\tilde{\mathsf{\mu }}L$ and $c_v\tilde{\mathsf{\mu }}\sqrt{L}$, respectively. Note that the true distribution of the lead-time demand is unknown to the decision-maker in a real-world application, but only an estimate of its mean and its variance can be obtained. Parameters h, A, π, π₀, L, and β are instead reasonably supposed to be known to the decision-maker. The observations of the lead-time demand are used to calculate a guess of its true mean and variance. These estimates are then exploited to find the (sub-)optimal replenishment policy in each “approximated” model considered (imposing that the mean and the variance in the approximation model are equal to the estimates of the true statistics). We also determine the optimal replenishment policy, and the corresponding minimum cost, we would have adopted under complete information (i.e., if the true lead-time demand distribution were known). Finally, we evaluate the true cost of the ordering rule obtained under Gaussian approximation, minimax approach or MaxEnt principle. These true costs are compared in terms of Absolute Percentage Error with respect to the true minimum cost.

With a certain lead-time demand distribution and for given parameter values, the optimal policy is determined by minimizing Equation (1) in $\left(Q,r\right)$. The quantity $\overline{B}\left(r\right)$ can readily be obtained by means of Equation (2) once the density $\overline{f}$ of the lead-time demand is specified. We remind that the density p of the lead-time demand under the minimax approach is given by [39]:

$$p\left(z\right)=\frac{1}{2}\frac{1}{\sqrt{{\left(1+z^2\right)}^3}}$$

where $z=\left(r-\mathsf{\mu }\right)/\mathsf{\sigma }$.

Concerning the task of assessing the mean and the variance of the lead-time demand, we would observe that ten observations are used as a trade-off between obtaining a good estimate and a too special value. A similar argument was raised in [20].

For each combination of parameter values, under a given true lead-time demand distribution, five independent runs have been made. Results are shown in Table 4 and in Table 5 for the cases $h=5$ $/unit/year and $h=15$ $/unit/year, respectively. In these tables, the smallest error achieved in each run is written in bold.

Table 4. Results of the comparative analysis for the case $h=5$ $/unit/year. L = (days).

		Lognormal			Weibull			Gamma
cv	L	Gauss.	Minimax	MaxEnt	Gauss.	Minimax	MaxEnt	Gauss.	Minimax	MaxEnt
0.05	10	0.08%	1.42%	0.07%	0.10%	1.81%	0.09%	0.65%	4.36%	0.63%
		0.20%	0.74%	0.17%	0.63%	0.82%	0.62%	0.01%	2.08%	0.01%
		0.15%	1.53%	0.13%	0.78%	0.91%	0.74%	0.09%	2.15%	0.09%
		0.23%	0.61%	0.14%	0.20%	2.04%	0.20%	0.43%	0.85%	0.43%
		0.25%	2.91%	0.62%	0.97%	4.76%	0.99%	0.89%	4.99%	0.90%
	30	2.16%	9.74%	2.18%	0.45%	2.21%	0.45%	0.81%	6.66%	0.81%
		0.58%	1.49%	0.58%	0.74%	6.12%	0.74%	1.74%	0.72%	1.74%
		1.94%	0.54%	1.94%	0.28%	2.17%	0.28%	0.32%	1.91%	0.32%
		1.01%	1.29%	1.01%	0.07%	3.66%	0.07%	1.01%	6.88%	1.01%
		0.09%	2.93%	0.09%	0.20%	3.74%	0.20%	0.68%	4.50%	0.68%
0.1	10	2.17%	0.35%	0.88%	2.87%	0.50%	1.88%	8.37%	0.48%	0.34%
		0.48%	1.73%	0.03%	0.08%	3.45%	0.10%	0.83%	1.56%	0.18%
		0.75%	1.41%	0.31%	1.02%	1.59%	0.45%	1.13%	9.12%	2.58%
		0.35%	3.55%	0.32%	1.13%	1.36%	0.20%	0.24%	3.92%	0.20%
		0.67%	1.82%	0.48%	1.03%	1.47%	0.22%	0.35%	2.62%	0.10%
	30	2.17%	13.9%	2.32%	1.41%	3.05%	1.31%	1.83%	2.11%	1.82%
		1.87%	2.07%	1.62%	3.85%	16.29%	4.79%	2.01%	2.31%	1.29%
		4.82%	0.85%	4.75%	4.39%	1.28%	4.39%	0.39%	7.46%	0.39%
		0.21%	4.30%	0.20%	1.71%	2.70%	1.66%	6.47%	21.88%	6.77%
		0.40%	8.38%	0.41%	1.68%	2.62%	1.67%	4.08%	1.18%	3.52%
0.4	10	8.89%	0.14%	6.73%	4.42%	1.00%	2.91%	1.83%	14.07%	0.78%
		1.95%	5.35%	1.94%	2.82%	16.90%	1.25%	12.33%	41.67%	8.31%
		12.65%	1.82%	8.45%	15.54%	44.96%	11.27%	5.72%	35.48%	9.75%
		2.22%	1.46%	1.19%	44.09%	17.38%	38.74%	0.84%	16.47%	2.11%
		17.75%	56.42%	13.32%	3.27%	2.18%	2.12%	33.77%	12.12%	33.13%
	30	30.02%	77.78%	41.28%	3.14%	4.73%	0.47%	1.74%	22.78%	4.97%
		7.30%	1.08%	3.17%	6.71%	2.36%	1.75%	17.78%	2.52%	11.01%
		3.39%	5.79%	2.36%	48.75%	12.09%	35.85%	7.78%	2.73%	1.87%
		4.40%	3.13%	2.07%	3.71%	27.39%	7.84%	2.30%	6.14%	0.07%
		2.41%	5.31%	1.66%	4.54%	3.78%	1.58%	4.03%	4.59%	0.45%
0.8	10	4.28%	4.06%	3.96%	62.38%	29.48	57.97%	59.99%	27.33%	56.03%
		19.38%	44.14%	15.34%	11.06%	9.85%	2.69%	7.36%	9.78%	7.20%
		4.15%	13.00%	5.18%	42.95%	39.25%	16.66%	38.41%	97.31%	32.12%
		18.69%	52.28%	8.56%	5.78%	15.41%	6.80%	15.63%	14.12%	11.44%
		49.02%	21.78%	42.87%	23.90%	4.51%	15.76%	6.47%	24.02%	8.12%
	30	12.76%	14.23%	12.65%	10.02%	19.32%	9.11%	9.48%	26.01%	7.38%
		10.68%	12.15%	10.10%	11.72%	9.36%	5.53%	35.53%	56.80%	24.43%
		22.08%	49.01%	26.87%	13.32%	23.93%	11.69%	16.02%	0.84%	9.35%
		59.61%	95.20%	50.06%	44.93%	79.67%	36.74%	19.49%	13.96%	1.23%
		19.23%	1.42%	15.25%	59.82%	14.67%	47.38%	6.79%	15.96%	6.79%

Table 4. Results of the comparative analysis for the case $h=5$ $/unit/year. L = (days).

		Lognormal			Weibull			Gamma
cv	L	Gauss.	Minimax	MaxEnt	Gauss.	Minimax	MaxEnt	Gauss.	Minimax	MaxEnt
0.05	10	0.08%	1.42%	0.07%	0.10%	1.81%	0.09%	0.65%	4.36%	0.63%
		0.20%	0.74%	0.17%	0.63%	0.82%	0.62%	0.01%	2.08%	0.01%
		0.15%	1.53%	0.13%	0.78%	0.91%	0.74%	0.09%	2.15%	0.09%
		0.23%	0.61%	0.14%	0.20%	2.04%	0.20%	0.43%	0.85%	0.43%
		0.25%	2.91%	0.62%	0.97%	4.76%	0.99%	0.89%	4.99%	0.90%
	30	2.16%	9.74%	2.18%	0.45%	2.21%	0.45%	0.81%	6.66%	0.81%
		0.58%	1.49%	0.58%	0.74%	6.12%	0.74%	1.74%	0.72%	1.74%
		1.94%	0.54%	1.94%	0.28%	2.17%	0.28%	0.32%	1.91%	0.32%
		1.01%	1.29%	1.01%	0.07%	3.66%	0.07%	1.01%	6.88%	1.01%
		0.09%	2.93%	0.09%	0.20%	3.74%	0.20%	0.68%	4.50%	0.68%
0.1	10	2.17%	0.35%	0.88%	2.87%	0.50%	1.88%	8.37%	0.48%	0.34%
		0.48%	1.73%	0.03%	0.08%	3.45%	0.10%	0.83%	1.56%	0.18%
		0.75%	1.41%	0.31%	1.02%	1.59%	0.45%	1.13%	9.12%	2.58%
		0.35%	3.55%	0.32%	1.13%	1.36%	0.20%	0.24%	3.92%	0.20%
		0.67%	1.82%	0.48%	1.03%	1.47%	0.22%	0.35%	2.62%	0.10%
	30	2.17%	13.9%	2.32%	1.41%	3.05%	1.31%	1.83%	2.11%	1.82%
		1.87%	2.07%	1.62%	3.85%	16.29%	4.79%	2.01%	2.31%	1.29%
		4.82%	0.85%	4.75%	4.39%	1.28%	4.39%	0.39%	7.46%	0.39%
		0.21%	4.30%	0.20%	1.71%	2.70%	1.66%	6.47%	21.88%	6.77%
		0.40%	8.38%	0.41%	1.68%	2.62%	1.67%	4.08%	1.18%	3.52%
0.4	10	8.89%	0.14%	6.73%	4.42%	1.00%	2.91%	1.83%	14.07%	0.78%
		1.95%	5.35%	1.94%	2.82%	16.90%	1.25%	12.33%	41.67%	8.31%
		12.65%	1.82%	8.45%	15.54%	44.96%	11.27%	5.72%	35.48%	9.75%
		2.22%	1.46%	1.19%	44.09%	17.38%	38.74%	0.84%	16.47%	2.11%
		17.75%	56.42%	13.32%	3.27%	2.18%	2.12%	33.77%	12.12%	33.13%
	30	30.02%	77.78%	41.28%	3.14%	4.73%	0.47%	1.74%	22.78%	4.97%
		7.30%	1.08%	3.17%	6.71%	2.36%	1.75%	17.78%	2.52%	11.01%
		3.39%	5.79%	2.36%	48.75%	12.09%	35.85%	7.78%	2.73%	1.87%
		4.40%	3.13%	2.07%	3.71%	27.39%	7.84%	2.30%	6.14%	0.07%
		2.41%	5.31%	1.66%	4.54%	3.78%	1.58%	4.03%	4.59%	0.45%
0.8	10	4.28%	4.06%	3.96%	62.38%	29.48	57.97%	59.99%	27.33%	56.03%
		19.38%	44.14%	15.34%	11.06%	9.85%	2.69%	7.36%	9.78%	7.20%
		4.15%	13.00%	5.18%	42.95%	39.25%	16.66%	38.41%	97.31%	32.12%
		18.69%	52.28%	8.56%	5.78%	15.41%	6.80%	15.63%	14.12%	11.44%
		49.02%	21.78%	42.87%	23.90%	4.51%	15.76%	6.47%	24.02%	8.12%
	30	12.76%	14.23%	12.65%	10.02%	19.32%	9.11%	9.48%	26.01%	7.38%
		10.68%	12.15%	10.10%	11.72%	9.36%	5.53%	35.53%	56.80%	24.43%
		22.08%	49.01%	26.87%	13.32%	23.93%	11.69%	16.02%	0.84%	9.35%
		59.61%	95.20%	50.06%	44.93%	79.67%	36.74%	19.49%	13.96%	1.23%
		19.23%	1.42%	15.25%	59.82%	14.67%	47.38%	6.79%	15.96%	6.79%

Table 5. Results of the comparative analysis for the case $h=15$ $/unit/year. L = (days).

		Lognormal			Weibull			Gamma
cv	L	Gauss.	Minimax	MaxEnt	Gauss.	Minimax	MaxEnt	Gauss.	Minimax	MaxEnt
0.05	10	0.61%	0.98%	0.61%	5.81%	1.69%	5.81%	0.10%	1.00%	0.08%
		3.69%	8.90%	4.14%	0.97%	1.15%	0.96%	0.58%	3.22%	0.58%
		5.22%	7.15%	1.54%	0.24%	2.34%	0.23%	0.17%	0.69%	0.14%
		0.68%	0.21%	0.67%	0.22%	2.46%	0.22%	5.97%	1.28%	4.60%
		0.04%	0.95%	0.04%	0.29%	0.55%	0.28%	0.02%	1.49%	0.02%
	30	4.51%	0.87%	4.51%	1.60%	6.43%	1.60%	0.35%	3.58%	0.35%
		0.10%	1.71%	0.10%	0.04%	2.05%	0.04%	0.59%	4.71%	0.59%
		0.05%	1.93%	0.05%	0.39%	4.13%	0.39%	0.05%	1.76%	0.05%
		1.50%	6.22%	0.87%	0.89%	5.08%	0.89%	12.55%	3.55%	12.55%
		0.69%	0.85%	0.54%	1.59%	0.57%	1.59%	0.09%	1.68%	0.09%
0.1	10	1.90%	0.25%	1.01%	7.67%	11.49%	3.53%	1.01%	1.07%	0.90%
		6.06%	11.44%	3.39%	3.57%	11.81%	5.29%	0.55%	3.82%	0.13%
		0.95%	0.47%	0.38%	3.44%	1.05%	3.30%	0.40%	3.92%	0.20%
		0.74%	0.58%	0.18%	1.59%	7.46%	1.79%	4.60%	0.99%	3.79%
		19.46%	9.78%	2.05%	9.64%	3.09%	0.13%	0.52%	0.96%	0.08%
	30	0.08%	3.35%	0.06%	16.25%	5.70%	15.32%	3.22%	1.34%	2.20%
		0.38%	6.00%	0.46%	1.13%	8.15%	1.07%	4.01%	1.46%	3.80%
		1.79%	1.23%	1.77%	3.05%	2.92%	1.67%	0.51%	6.69%	0.08%
		1.11%	1.63%	0.99%	0.19%	3.12%	0.15%	0.32%	5.84%	0.12%
		2.88%	11.44%	2.80%	0.02%	4.12%	0.02%	0.01%	4.30%	0.01%
0.4	10	42.57%	10.27%	5.72%	45.67%	39.91%	25.42%	45.89%	40.65%	20.36%
		7.16%	5.78%	0.50%	3.59%	2.09%	2.99%	24.15%	20.00%	7.81%
		71.38%	63.72%	52.63%	2.13%	6.52%	1.87%	43.84%	38.95%	21.29%
		33.73%	16.26%	29.05%	77.76%	70.58%	53.89%	16.72%	35.84%	12.60%
		7.53%	18.78%	9.79%	78.57%	96.12%	67.02%	42.80%	85.00%	7.28%
	30	62.65%	34.62%	28.38%	60.80%	26.07%	49.36%	6.62%	0.83%	2.57%
		9.05%	0.17%	4.59%	32.01%	48.62%	24.21%	52.75%	40.63%	19.79%
		20.44%	2.18%	12.44%	10.02%	16.37%	7.68%	8.27%	16.60%	5.99%
		7.11%	13.86%	9.19%	10.14%	0.34%	3.81%	84.95%	35.44%	15.45%
		10.49%	15.52%	7.41%	15.16%	25.96%	9.61%	9.56%	0.66%	3.92%
0.8	10	56.30%	52.72%	39.39%	39.65%	60.42%	34.22%	73.52%	19.11%	44.42%
		49.82%	44.38%	30.34%	34.35%	14.10%	9.72%	50.13%	46.24%	18.41%
		32.46%	14.42%	25.41%	50.12%	46.07%	20.51%	56.26%	51.99%	21.75%
		81.37%	76.52%	54.42%	28.73%	8.95%	15.81%	25.50%	6.47%	23.44%
		37.60%	34.12%	15.54%	63.66%	59.47%	34.73%	17.25%	18.67%	16.50%
	30	56.91%	79.41%	63.78%	18.99%	5.04%	15.77%	17.16%	5.73%	14.44%
		17.84%	12.34%	16.83%	72.99%	25.19%	62.53%	57.31%	46.60%	13.48%
		21.47%	15.97%	2.85%	67.05%	54.49%	21.08%	16.78%	13.63%	6.51%
		23.91%	18.87%	2.21%	87.97%	76.52%	30.85%	18.06%	4.48%	12.11%
		16.40%	15.76%	13.40%	25.65%	37.94%	19.86%	20.36%	16.52%	2.61%

Table 5. Results of the comparative analysis for the case $h=15$ $/unit/year. L = (days).

		Lognormal			Weibull			Gamma
cv	L	Gauss.	Minimax	MaxEnt	Gauss.	Minimax	MaxEnt	Gauss.	Minimax	MaxEnt
0.05	10	0.61%	0.98%	0.61%	5.81%	1.69%	5.81%	0.10%	1.00%	0.08%
		3.69%	8.90%	4.14%	0.97%	1.15%	0.96%	0.58%	3.22%	0.58%
		5.22%	7.15%	1.54%	0.24%	2.34%	0.23%	0.17%	0.69%	0.14%
		0.68%	0.21%	0.67%	0.22%	2.46%	0.22%	5.97%	1.28%	4.60%
		0.04%	0.95%	0.04%	0.29%	0.55%	0.28%	0.02%	1.49%	0.02%
	30	4.51%	0.87%	4.51%	1.60%	6.43%	1.60%	0.35%	3.58%	0.35%
		0.10%	1.71%	0.10%	0.04%	2.05%	0.04%	0.59%	4.71%	0.59%
		0.05%	1.93%	0.05%	0.39%	4.13%	0.39%	0.05%	1.76%	0.05%
		1.50%	6.22%	0.87%	0.89%	5.08%	0.89%	12.55%	3.55%	12.55%
		0.69%	0.85%	0.54%	1.59%	0.57%	1.59%	0.09%	1.68%	0.09%
0.1	10	1.90%	0.25%	1.01%	7.67%	11.49%	3.53%	1.01%	1.07%	0.90%
		6.06%	11.44%	3.39%	3.57%	11.81%	5.29%	0.55%	3.82%	0.13%
		0.95%	0.47%	0.38%	3.44%	1.05%	3.30%	0.40%	3.92%	0.20%
		0.74%	0.58%	0.18%	1.59%	7.46%	1.79%	4.60%	0.99%	3.79%
		19.46%	9.78%	2.05%	9.64%	3.09%	0.13%	0.52%	0.96%	0.08%
	30	0.08%	3.35%	0.06%	16.25%	5.70%	15.32%	3.22%	1.34%	2.20%
		0.38%	6.00%	0.46%	1.13%	8.15%	1.07%	4.01%	1.46%	3.80%
		1.79%	1.23%	1.77%	3.05%	2.92%	1.67%	0.51%	6.69%	0.08%
		1.11%	1.63%	0.99%	0.19%	3.12%	0.15%	0.32%	5.84%	0.12%
		2.88%	11.44%	2.80%	0.02%	4.12%	0.02%	0.01%	4.30%	0.01%
0.4	10	42.57%	10.27%	5.72%	45.67%	39.91%	25.42%	45.89%	40.65%	20.36%
		7.16%	5.78%	0.50%	3.59%	2.09%	2.99%	24.15%	20.00%	7.81%
		71.38%	63.72%	52.63%	2.13%	6.52%	1.87%	43.84%	38.95%	21.29%
		33.73%	16.26%	29.05%	77.76%	70.58%	53.89%	16.72%	35.84%	12.60%
		7.53%	18.78%	9.79%	78.57%	96.12%	67.02%	42.80%	85.00%	7.28%
	30	62.65%	34.62%	28.38%	60.80%	26.07%	49.36%	6.62%	0.83%	2.57%
		9.05%	0.17%	4.59%	32.01%	48.62%	24.21%	52.75%	40.63%	19.79%
		20.44%	2.18%	12.44%	10.02%	16.37%	7.68%	8.27%	16.60%	5.99%
		7.11%	13.86%	9.19%	10.14%	0.34%	3.81%	84.95%	35.44%	15.45%
		10.49%	15.52%	7.41%	15.16%	25.96%	9.61%	9.56%	0.66%	3.92%
0.8	10	56.30%	52.72%	39.39%	39.65%	60.42%	34.22%	73.52%	19.11%	44.42%
		49.82%	44.38%	30.34%	34.35%	14.10%	9.72%	50.13%	46.24%	18.41%
		32.46%	14.42%	25.41%	50.12%	46.07%	20.51%	56.26%	51.99%	21.75%
		81.37%	76.52%	54.42%	28.73%	8.95%	15.81%	25.50%	6.47%	23.44%
		37.60%	34.12%	15.54%	63.66%	59.47%	34.73%	17.25%	18.67%	16.50%
	30	56.91%	79.41%	63.78%	18.99%	5.04%	15.77%	17.16%	5.73%	14.44%
		17.84%	12.34%	16.83%	72.99%	25.19%	62.53%	57.31%	46.60%	13.48%
		21.47%	15.97%	2.85%	67.05%	54.49%	21.08%	16.78%	13.63%	6.51%
		23.91%	18.87%	2.21%	87.97%	76.52%	30.85%	18.06%	4.48%	12.11%
		16.40%	15.76%	13.40%	25.65%	37.94%	19.86%	20.36%	16.52%	2.61%

We can first observe that relative performances do not significantly change for different values of h, as well as varying L for fixed c_v. We can also note that the error of all approximation methods is increasing as c_v and h become larger, which confirms results of the error analysis.

For small c_v, the performance of the (r,Q) policy with Gaussian lead-time demand or under MaxEnt principle is substantially similar. In fact, it is known that the Gaussian approximation works well for small c_v, as in such case the normal density is nearly 0 in the negative real semi-axis. In this condition, we can argue that the normal density and the maximum entropy distribution are very close. In contrast, the performance of the Gaussian approximation evidently deteriorates for higher c_v. We can also observe that its performance worsens as h increases for fixed c_v.

With regard to the minimax approach, its performance improves as the coefficient of variation of the lead-time demand increases. That is, when higher variability is involved, it seems a good choice to approximate the true lead-time demand distribution. In contrast, the approximation provided by the normal density appears to be preferable for small c_v.

However, results clearly highlight that using the MaxEnt principle seems the best choice, under the considered conditions. In fact, its performance is very good under all investigated configurations and overcomes that of the other approaches in the majority of cases analyzed. Moreover, it looks to be not much sensitive to variations in the parameter values. This is a significant outcome, as it makes the MaxEnt principle a promising method to model the lead-time demand distribution when the decision-maker is provided with limited information about the true distribution.

We have then carried out additional experiments to compare the performance of the MaxEnt and Weibull distributions in approximating the true lead-time demand distribution. Note that the Weibull density is not typically adopted to represent the lead-time demand (more generally, the demand in a given time interval) [39,41,46]. However, it is known to have great flexibility to model many types of data, in particular thanks to the shape parameter k that allows the density to attain several shapes [47].

These tests are performed similarly to those presented in the first part of this subsection; that is, the procedure is the same. Again, Table 3 shows the parameters whose value is kept fixed. In this session, we have considered only one value for h, i.e., h = 5 $/unit/year. Parameters c_v and L take the same values that have been adopted in the previous experiments. The true density $p\left(x\right)$ of the lead-time demand is assumed to be expressed as follows:

$$p\left(x\right)=0.4\cdot \overline{p}\left(x;{\mathsf{\mu }}_1,{\mathsf{\sigma }}_1^2\right)+0.3\cdot \overline{p}\left(x;{\mathsf{\mu }}_2,{\mathsf{\sigma }}_2^2\right)+0.3\cdot \overline{p}\left(x;{\mathsf{\mu }}_3,{\mathsf{\sigma }}_3^2\right)$$

where $\overline{p}\left(x;{\mathsf{\mu }}_i,{\mathsf{\sigma }}_i^2\right)$, for $i=1,2,3$, is a lognormal density with parameters ${\mathsf{\mu }}_i$ and ${\mathsf{\sigma }}_i^2$. That is, $p\left(x\right)$ is a mixture of lognormal densities. This choice is not based on a specific criterion; basically, we have considered a density that was not “standard”, i.e., belonging to a specific class. Parameters ${\mathsf{\mu }}_i$ and ${\mathsf{\sigma }}_i^2$, for $i=2,3$, are kept fixed and take the following values: ${\mathsf{\mu }}_2=1.07$, ${\mathsf{\mu }}_3=2.69$, ${\mathsf{\sigma }}_2^2=0.23$, and ${\mathsf{\sigma }}_3^2=0.16$. These values are purely indicative. In each experiment, the value of ${\mathsf{\mu }}_1$ and ${\mathsf{\sigma }}_1^2$ is determined to assure that the mean and the standard deviation of $p\left(x\right)$ are equal to $\tilde{\mathsf{\mu }}L$ and $c_v\tilde{\mathsf{\mu }}\sqrt{L}$, respectively (we, in fact, remind that the lead-time demand has mean and standard deviation respectively given by $\tilde{\mathsf{\mu }}L$ and $c_v\tilde{\mathsf{\mu }}\sqrt{L}$, where $c_v\equiv \tilde{\mathsf{\sigma }}/\tilde{\mathsf{\mu }}$).

For each combination of parameter values, three runs have been made. Table 6 shows the results of the comparison between the MaxEnt and Weibull distributions. Performance is measured in terms of Absolute Percentage Error (APE) with respect to the true minimum cost. In addition, note that Table 6 gives the value that the shape parameter k of the Weibull density takes in each run. We can note that the MaxEnt distribution achieves a better performance, as this model has been able to obtain a smaller APE in the majority of tests. In addition, the performance of the Weibull model seems to deteriorate, with respect to the performance of the MaxEnt model, as the variability in the system grows. That is, with increasing c_v, the MaxEnt distribution has turned out to realize the lowest APE more frequently than for small/medium values of c_v. With regard to the APE magnitude, we can observe that it is increasing in c_v, as expected. This result is in accordance to the outcomes in the previous experiments. With regard to the shape parameter k of the Weibull density, this has taken relatively small values. In particular, the greatest value has been observed to be equal to 1.68. A final remark: k appears to be decreasing as the variability in the system increases; that is, in such circumstance, the Weibull density tends to have null mode.

Table 6. Results of the comparison between MaxEnt and Weibull distributions. L = (days).

cv	L	MaxEnt	Weibull	k
0.05	10	3.54%	4.61%	1.68
		2.39%	3.22%	1.41
		0.72%	1.04%	1.24
	30	6.19%	6.84%	1.46
		6.46%	6.32%	1.28
		6.07%	5.98%	1.19
0.1	10	12.15%	15.04%	1.08
		12.22%	11.95%	1.22
		0.64%	0.99%	1.27
	30	4.76%	4.92%	0.58
		0.49%	1.37%	0.75
		1.56%	1.20%	0.67
0.4	10	13.70%	12.45%	1.14
		12.12%	24.11%	0.73
		8.22%	8.56%	1.31
	30	6.86%	7.27%	0.44
		3.05%	5.63%	1.02
		6.84%	13.94%	0.57
0.8	10	80.44%	86.71%	1.54
		59.43%	57.70%	1.14
		7.84%	14.20%	0.72
	30	11.28%	3.31%	0.63
		13.28%	21.06%	0.95
		10.27%	13.27%	0.76

Table 6. Results of the comparison between MaxEnt and Weibull distributions. L = (days).

cv	L	MaxEnt	Weibull	k
0.05	10	3.54%	4.61%	1.68
		2.39%	3.22%	1.41
		0.72%	1.04%	1.24
	30	6.19%	6.84%	1.46
		6.46%	6.32%	1.28
		6.07%	5.98%	1.19
0.1	10	12.15%	15.04%	1.08
		12.22%	11.95%	1.22
		0.64%	0.99%	1.27
	30	4.76%	4.92%	0.58
		0.49%	1.37%	0.75
		1.56%	1.20%	0.67
0.4	10	13.70%	12.45%	1.14
		12.12%	24.11%	0.73
		8.22%	8.56%	1.31
	30	6.86%	7.27%	0.44
		3.05%	5.63%	1.02
		6.84%	13.94%	0.57
0.8	10	80.44%	86.71%	1.54
		59.43%	57.70%	1.14
		7.84%	14.20%	0.72
	30	11.28%	3.31%	0.63
		13.28%	21.06%	0.95
		10.27%	13.27%	0.76

6. Conclusions and Further Remarks

In this paper, we took into account the continuous-review reorder point-lot size (r,Q) policy under stochastic demand, with backorders-lost sales mixture. We modelled the lead-time demand distribution according to the MaxEnt principle. That is, we considered the density that maximizes the entropy over all densities with given mean and variance. This approach is suitable when the decision-maker is provided with limited information about the true distribution of the lead-time demand. This is the case, for example, in which the true distribution is unknown, but only some observations of the demand are available, which allow estimating the mean and the variance of the lead-time demand.

We developed an optimization problem aimed at minimizing the expected total cost per time unit. To approach this problem, we then presented an exact algorithm and a heuristic method. The heuristic algorithm is based on an approximated expression of the cost function achieved by means of an ad hoc first-order truncation of the Taylor series expansion.

Numerical experiments were finally carried out to investigate two aspects. First, the heuristic algorithm was evaluated in terms of both error achieved and computational effort required. Tests proved that the approximated solution procedure is efficient for a wide range of parameter values. Then, we investigated the capability of the MaxEnt principle to approximate the true lead-time demand distribution under partial information. This analysis was firstly made considering the performance of two alternative approaches: the Gaussian approximation and the minimax method. Results highlighted that the MaxEnt principle performed better that the other methods. In fact, it achieved the best outcome in the majority of cases analyzed. A second set of experiments was then carried out to compare the MaxEnt principle with respect to the approximation provided by the Weibull density. Again, the MaxExt principle turned out to be preferable, as it achieved the best result in the majority of tests. In conclusion, the MaxEnt principle seems a promising practical approach to model the lead-time demand distribution under partial information.

The developed model considers the lead time a constant and deterministic quantity, which can therefore be supposed to be known to the decision-maker. However, the case of a random lead time can be taken into account, too. In this regard, it is reasonable to follow a similar argument to what said above concerning the distribution of the lead-time demand. That is, when the lead time is a random variable, full information about its true distribution may not be available to the decision-maker. In this circumstance, only some observations about the lead-time duration can be obtained. To include stochastic lead time into our model, in the relevant case characterized by partial information, we can proceed as follows. A premise is needed: we refer to the condition with sequential deliveries independent of the lead-time demand, which is characterized by orders that cannot cross in time [46]. Note that this is the most common situation in practice, since it is almost always true that orders are received in the same sequence in which they were placed [41,46]. Hence, let us assume that the decision-maker is aware of only some observations about (i) the demand in given time periods; and (ii) the length of lead time. These observations allow estimating the mean and the variance of the lead time (i.e., ${\mathsf{\mu }}_L$ and ${\mathsf{\sigma }}_L^2$, respectively) and of the demand per time unit (i.e., ${\mathsf{\mu }}_D$ and ${\mathsf{\sigma }}_D^2$, respectively). These quantities can then be adopted to assess the mean $\mathsf{\mu }$ and the variance ${\mathsf{\sigma }}^2$ of the lead-time demand according to the following expressions [46]:

$$\begin{array}{l}\mathsf{\mu }={\mathsf{\mu }}_D{\mathsf{\mu }}_L,\hfill \\ {\mathsf{\sigma }}^2={\mathsf{\sigma }}_D^2{\mathsf{\mu }}_L+{\mathsf{\mu }}_D^2{\mathsf{\sigma }}_L^2.\hfill \end{array}$$

Finally. $\mathsf{\mu }$ and ${\mathsf{\sigma }}^2$ can be exploited to estimate parameters a, b, and c of the density maximizing the entropy, which will be used to model the distribution of the lead-time demand, as described in Section 2.

Future researches may be devoted to several studies. For example, the suitability to implement the MaxEnt principle into different inventory systems (e.g., multi-echelon supply chains) may be investigated. In addition, the entropy maximization principle may be applied to generalized entropies, e.g., to Tsallis entropy [48]. Finally, a more general formulation of the expected shortage may also be considered. In particular, the expression known in the risk management literature as “shortfall risk” [49,50] may be taken into account in place of the standard quantity.