An Inventory Service-Level Optimization Problem for a Multi-Warehouse Supply Chain Network with Stochastic Demands

Abstract

This research aims to develop a mathematical model and a solution approach for jointly optimizing a global inventory service level and order sizes for a single-commodity supply chain network with multiple warehouses or distribution centers. The latter face stochastic demands, such as most real-world supply chains do nowadays, yielding significant model complexity. The studied problem is of high relevance for inventory management, inventory location, and supply chain network design-related literature, as well as for logistics and supply chain managers. The proposed optimization model minimizes the total costs associated with cycle inventory, safety stock, and stock-out-related events, considering a global inventory service level and differentiated order sizes for a fixed and known set of warehouses. Subsequently, the model is solved by employing the Newton–Raphson algorithm, which is developed and implemented assuming stochastic demands with a normal approximation. The algorithm reached optimality conditions and the convergence criterion in a few iterations, within less than a second, for a variety of real-world sized instances involving up to 200 warehouses. The model solutions are contrasted with those obtained with a previous widely employed approximation, where safety stock costs were further approximated and order sizes were optimized without considering stock-out-related costs. This comparison denotes valuable benefits without significant additional computational efforts. Thus, the proposed approach is suitable for managers of real-world supply chains, since they would be able to attain system performance improvements by simultaneously optimizing the global inventory service level and order sizes, thereby providing a better system cost equilibrium.

1. Introduction

Supply chain management aims to optimize and coordinate decisions and agents to provide products and related services at the right times, quantities, and places, also aiming to minimize system-wide costs and ensuring system service-level requirements (see [1,2,3,4,5]). In this context, inventory control decisions and related costs represent key elements to be addressed, in addition to transport and warehouse operations, among other issues (see [6,7]). One of the most typical supply chain topologies comprises a set of one or more plants or suppliers that distribute products to a set of parallel warehouses or distribution centers that finally serve end customers’ demands, as shown in Figure 1 (see [3,8,9,10,11,12]). Under this topology, inventory control decisions at the warehouse level are of high significance, wherein order sizes, reorder points, and service levels represent the main decisions to be addressed. This has to be achieved while being efficient and following a number of performance indicators, such as inventory availability, stock-out probability, or demand fill rate (see [3,12,13,14]). System costs have been traditionally considered as the main objective in evaluating global performance, where obtaining even marginal savings in relative terms may generate significant improvements in the medium or long term.

From a strategic or long-term perspective, designing the underlying network is an essential problem, which typically consists of locating warehouses and plants for serving the end customers, conforming to the well-known family of facility location problems (FLPs), as shown in [15,16,17,18,19]. Traditionally, plant or warehouse locations and customer assignment decisions are made based on a sequential and hierarchical decision-making structure, where decisions belonging to higher organizational levels are made based on FLP modeling structure, without major concerns regarding the impacts on further lower-level decisions (see [3,11,12]).

Nonetheless, it is worth noticing that, when the supply chain network (SCN) is already optimized and fixed for a specific time frame (e.g., a few weeks, months, or years), the optimization of inventory-related decisions remains a relevant problem to be addressed. In this context, several research works have been developed, where Refs. [20,21] developed approaches to evaluate system performance in a two-level supply chain under different assumptions and considerations. In [22], alternative optimization approaches were developed to deal with inventory decisions, where order sizes and reorder points for a two-level supply chain were optimized. A mathematical model to optimize inventory-related decisions for a complex supply chain network was developed in [23], comprising end customers, retailers, warehouses, productive plants, and raw material suppliers, assuming a fixed service level to be ensured for end customers. Further extensions were presented in [24,25], where the first one addresses a multi-echelon supply chain system with three different demand patterns, while the second one includes product perishability and demands depending on prices and inventory levels. Other recent research works can be found in [26], which presents a deterministic supply chain inventory model involving outsourcing decisions, and [27], which addresses a deterministic single-location inventory model involving product deterioration, a supplier offering price discounts, and payment delays.

Most of the aforementioned works simultaneously optimize decisions at the warehouse level and at the plant or supplier level, consider reorder points and order sizes as decision variables, and yield different service levels for each location. Moreover, some works assume that the inventory service level is a fixed and known parameter, which should be exogenously determined by a supply chain planner or manager. In contrast, this research does not consider plant/supplier decisions and costs, and reorder points are determined through a global service-level optimization (i.e., the same service level for all warehouses). Then, reorder points are individualized for each warehouse mainly depending on lead times, served demand mean and variance, and the optimal global service level. Accordingly, this research aims to address situations in which a global service level along with order sizes are required to be optimized, in contrast to other previous related works.

As can be observed in traditional FLP-related literature, decisions related to the supply chain network design are made without considering specific impacts on inventory management ([9,10,11]). However, inventory management and strategic SCN design must take care of the impacts between each other, where the number of parallel facilities (e.g., warehouses or distribution centers) affects the total system inventory level, resulting in the well-known risk pooling effect, especially in scenarios with stochastic demands. Thus, several works show that the structure of the supply chain network directly relates to the impact that risk pooling has on system costs, conforming to the family of inventory location problems (ILPs) (see [8,9,10,11,28,29,30]), where inventory management issues are jointly addressed with SCN design decisions and costs. Subsequently, Refs. [3,12] developed reviews of ILP-related works. Furthermore, and following the basic inventory location-related aforementioned works, this kind of problem may include different relevant practical issues, such as multi-commodity and multi-period scenarios, inventory and transport capacities, alternative inventory control policies, spare parts supply chains, and supplier selection (see [31,32,33,34,35,36,37]).

In this context, the need for jointly optimizing the inventory-related service level as a problem decision variable appeared to be relevant research that needed to be addressed. Thus, Refs. [8,9] presented approximate strategies to model and solve an ILP facing customer demands with a Poisson distribution and a continuous $(S-1,S)$ inventory control policy. Subsequently, Ref. [38] presented the first research study with an approximate strategy to model and solve an inventory location problem with normal demands and a continuous $(Q,s)$ inventory control policy. Finally, Ref. [39] presented a variation of the model in [38] and considered a heuristic solution approach.

Furthermore, Refs. [37,40] presented a generalized Benders decomposition (GBD)-based solution approach that solved an ILP at optimality with normal demands and continuous $(Q,s)$ inventory control policy in competitive times. A significant contribution of this approach relies on the decomposition process, which yields a sub-problem that fully addresses all inventory-related decisions for a fixed SCN, which was previously determined to be a master problem at each algorithm iteration. These facts represent a significant potentiality, where any inventory model improvement may be integrated into an ILP, and it would be fully captured by the sub-problem following the GBD-based solution approach proposed in [37,40]. Similarly, several widely employed local search heuristics and metaheuristics may be considered to solve a related ILP with inventory service-level optimization, following [32,41,42,43,44]. For these approaches, the design of the SCN is explored by a local search algorithm (i.e., warehouse location and customer assignments), whereas solving a fixed-network inventory optimization problem (as modeled and solved in this research) may be considered a solution evaluation subroutine within the overall algorithm. Consequently, the aim of jointly optimizing order sizes and the service level for a fixed SCN gains attractiveness given the potentiality of integrating it within a broader ILP for optimizing the design of the SCN. Naturally, this and other similar integrative approaches remain as future research to be addressed.

In summary, this paper fills existing research gaps by developing a novel continuous, nonlinear optimization model to jointly optimize the global inventory service level and order sizes for a fixed supply chain network. The addressed network topology comprises multiple parallel warehouses under a normal approximation of end customer demands. A continuous $(Q,s)$ inventory control policy is considered, where a common global service level for all the existing warehouses has to be observed. It worth noting that the proposed formulation involves enhancement of some well-accepted approximations in previous related literature, especially focusing on safety stock and order size optimization. Finally, a multivariate Newton–Raphson-based algorithm is developed to solve the formulated problem at optimality. This well-known method has been successfully employed in supply chain management-related research, such as in [45,46,47].

As discussed in [48], the Newton method, also known as the Newton–Raphson method, is one of the most proper approaches for solving nonlinear equations and is particularly relevant when expressions for the gradient vector and the Jacobian matrix can be computed in an explicit and efficient manner. This approach involves a fine equilibrium between mathematical complexity and tractability while ensuring a quadratic convergence rate. Other approaches widely employed are variations and approximation of the Newton–Raphson approach, such as the Quasi-Newton and the Conjugate Gradient Method. Moreover, given the small- and medium-size instances considered in this research (i.e., at most 201 variables) and the possibility of having a potentially good initial guess for the order size and the global service level (as described in Section 4), an outstanding performance of the employed approach is observed.

The rest of the paper is organized as follows. Section 2 describes the studied problem and the proposed model formulation, including a generic formulation and a specified one considering demands with a normal approximation. Section 3 aims to describe the solution approach developed to solve the proposed formulation, which is derived specifically for demands with a normal distribution. Subsequently, Section 4 describes the computational experimentation carried out, then presents and discusses the obtained results. Section 5 presents a further discussion of the main results obtained in the previous section, along with a brief managerial insight discussion. Finally, conclusions and future research are presented in Section 6.

2. Problem Description and Formulation

2.1. Problem Description

The problem studied in this research consists of optimizing a global inventory control service level and individual order sizes in a system comprising a single plant or supplier that fulfills incoming orders from N parallel warehouses, which finally serve demands from external customers (retailers). This topology is depicted in Figure 1. The global service level consists of a homogeneous maximum stock-out probability that must be observed at each and every warehouse. In this case, it is assumed that each warehouse faces a stochastic demand that is the aggregation of several customer demands. Particularly, following most inventory location-related literature, the aggregated demand at the warehouse level is approximated considering a normal distribution, as in [49].

Note that, although the normal distribution presents a positive probability of observing negative values (which can be considered inadequate for most demand processes), it may still be a suitable approximation in strategic problems, especially when a low variation coefficient of warehouse demands is observed, as indicated in [50]. In these cases, the probability of realizing negative demands may be considered immaterial or insignificant (see Section 5), which provides a simple simulation exercise to denote these facts. In addition, the normal distribution involves clear advantages related to its analytical properties, specially for deriving exact expressions. Moreover, ref. [50] (pp. 156–157) states that, when the variation coefficient is less than $0.5$, the truncated normal distribution does not significantly differ from a normal approximation. Naturally, in cases with significant probabilities of negative demands, the normal approximation may be considered unsuitable, and other distributions or approximations should be considered, as in [49], where a truncated normal distribution was considered.

In the studied system, each warehouse i$(i=1,...,N)$ follows a continuous inventory control policy $(Q_i,R{P}_i)$, where $Q_i$ is the order size to be optimized, and $R{P}_i$ is the reorder point, which is the critical inventory level that generates an order to the plant. The order size must be optimized according to the costs associated with the well-known Wilson model plus penalty costs associated with demands that arise during stock-out events at each warehouse. Finally, the reorder point has to be set in order to observe the global service level or stock-out probability, while the associated safety stock costs are minimized as much as possible. Thus, $R{P}_i$ is set at the minimum value that guarantees the observation of the maximum stock-out probability.

2.2. Model Formulation

The studied formulation relies on the following notation and definitions:

Model parameters:

$O{C}_i$	:	Fixed order costs for each order submitted from warehouse i ($/order).
$H{C}_i$	:	Unit holding cost at warehouse i ($/unit-time).
$P{C}_i$	:	Penalty cost per unit of demand arising during stock-out events
		and served with a delay for each warehouse i ($/unit).
$L{T}_i$	:	Lead time for the incoming orders at warehouse i.
$D_i$	:	Demand mean per time unit (per day) at each warehouse i.
$V_i$	:	Variance of demand per time unit at each warehouse i.

Decision variables:

$Q_i$	:	Order size for each warehouse i.
$R{P}_i$	:	Reorder point at each warehouse i.
$\alpha$	:	Maximum stock-out probability accepted at each and every warehouse.
$\delta$	:	Minimum probability required for not observing stock-outs ($1-\alpha$).
$S{S}_i$	:	Average safety stock observed at each warehouse i.
$SO{D}_i$	:	Average demand arising during stock-out events and served with a
		delay at each warehouse i.

Under the aforementioned definitions and assumptions, each warehouse faces a continuous stochastic demand during lead time $L{T}_i$, defined as $y_i$, with a probability density function ${\varphi }_i\left(y_i\right)$. Then, the problem is formulated as the following nonlinear programming model:

$$\begin{array}{cccc}\hfill & \min \hfill & \hfill & \sum _{i=1}^N\left(O{C}_i{\displaystyle \frac{D_i}{Q_i}}+H{C}_i{\displaystyle \frac{Q_i}{2}}\right)+\sum _{i=1}^{N}H{C}_{i}S{S}_i+\sum _{i=1}^{N}P{C}_{i}SO{D}_i{\displaystyle \frac{D_i}{Q_i}}\hfill \\ \hfill & \mathrm{s}.\mathrm{t}.\hfill & \hfill & \hfill \end{array}$$

(1)

$$\begin{array}{cc}Pr\left(y_i\le R{P}_i\right)=1-\alpha \hfill & \hfill \forall i=1,...,N\end{array}$$

(2)

$$\begin{array}{ccccc}\hfill & \hfill & & S{S}_i={\int }_0^{R{P}_i}(R{P}_i-y_i){\varphi }_i\left(y_i\right)d{y}_i\hfill & \hfill \forall i=1,...,N\end{array}$$

(3)

$$\begin{array}{ccccc}\hfill & \hfill & & SO{D}_i={\int }_{R{P}_i}^{\infty }(y_i-R{P}_i){\varphi }_i\left(y_i\right)d{y}_i\hfill & \hfill \forall i=1,...,N\end{array}$$

(4)

$$\begin{array}{ccccc}\hfill & \hfill & & 0<\alpha <1\hfill & \hfill \end{array}$$

(5)

The first summation of the objective function in Equation (1) is the well-known Wilson model objective function, which considers order and cycle inventory costs associated with the order size at each warehouse (see [51,52]). The second term in the objective function is the expected safety stock costs at all warehouses. Finally, the third summation in (1) is the expected costs associated with demand arising during stock-out events at all warehouses. Note that it is assumed that all demands are eventually served, and these costs represent the additional expenditures of serving demands with some delay in case a stock-out event occurs (also known as backorders).

Equation (2) ensures that $R{P}_i$ is set to guarantee the minimum service level or maximum stock-out probability at each warehouse. Equation (3) determines the expected safety stock just before the orders arrive at each warehouse, which is equivalent to the reorder point minus the expected demand during lead time (integrating over demand values lower than $R{P}_i$; otherwise, it represents a stock-out). Equation (4) allows for computing the expected demand arising since stock-out is observed, until the next order arrives (i.e., it considers demand values over $R{P}_i$; otherwise, a stock-out situation is not observed).

Assuming a normal distribution, Equation (2) allows for setting the minimum reorder point at each warehouse i, as shown in Equation (6). In addition, Equations (3) and (4) can be written as in (7) and (8), respectively (see [38]).

$$\begin{array}{ccc}\hfill R{P}_i& =& D_{i}L{T}_i+Z_{1-\alpha }\sqrt{L{T}_i}\sqrt{V_i}\hfill \end{array}$$

(6)

$$\begin{array}{ccc}\hfill S{S}_i& =& {\int }_0^{R{P}_i}(R{P}_i-y_i){\varphi }_i\left(y_i\right)d{y}_i\approx Z_{1-\alpha }\sqrt{V_{i}L{T}_i}(1-\alpha )+{\displaystyle \frac{\sqrt{V_{i}L{T}_i}}{\sqrt{2\pi }}}{e}^{-{\displaystyle \frac{1}{2}}{Z}_{1-\alpha }^2}\hfill \end{array}$$

(7)

$$\begin{array}{ccc}\hfill SO{D}_i& =& {\int }_{R{P}_i}^{\infty }(y_i-R{P}_i){\varphi }_i\left(y_i\right)d{y}_i={\displaystyle \frac{\sqrt{V_{i}L{T}_i}}{\sqrt{2\pi }}}{e}^{-{\displaystyle \frac{1}{2}}{Z}_{1-\alpha }^2}-Z_{1-\alpha }\sqrt{V_{i}L{T}_i}\alpha \hfill \end{array}$$

(8)

Note that Equation (3) assumes a minimum demand of 0, while the normal distribution is unbounded. However, the right side of the Equation (7) is obtained assuming unbounded demand (i.e., $-\infty$ to $+\infty$), representing an approximation for cases with non-negative demands. This approximation is suitable in cases where the probability of observing negative values is close to 0 and where the approximation error may not be significant.

Replacing (6)–(8) into the model, it can be written as:

$$\begin{array}{cccc}\hfill & \min \hfill & \hfill & \sum _{i=1}^N\left(O{C}_i{\displaystyle \frac{D_i}{Q_i}}+H{C}_i{\displaystyle \frac{Q_i}{2}}\right)+\hfill \\ \hfill & \hfill & & \sum _{i=1}^{N}H{C}_i\left(Z_{1-\alpha }\sqrt{V_{i}L{T}_i}(1-\alpha )+{\displaystyle \frac{\sqrt{V_{i}L{T}_i}}{\sqrt{2\pi }}}{e}^{-{\displaystyle \frac{1}{2}}{Z}_{1-\alpha }^2}\right)+\hfill \\ \hfill & \hfill & & \sum _{i=1}^{N}P{C}_i\left({\displaystyle \frac{\sqrt{V_{i}L{T}_i}}{\sqrt{2\pi }}}{e}^{-{\displaystyle \frac{1}{2}}{Z}_{1-\alpha }^2}-Z_{1-\alpha }\sqrt{V_{i}L{T}_i}\alpha \right){\displaystyle \frac{D_i}{Q_i}}\hfill \\ \hfill & \mathrm{s}.\mathrm{t}.:\hfill & \hfill & \left(5\right)\hfill & \hfill \end{array}$$

(9)

It is worth mentioning that, when the penalty cost is zero $(PC=0)$, then the order size may be determined simply by following the well-known Wilson model, as shown in Expression (10) (see [51,52]). However, given the impact of the order size on the total penalty costs, as observed in the last term of Expression (9), the optimal order size does not coincide with (10). Instead, the order size is computed as Expression (11) by differentiating Equation (9) with respect to $Q_i$ and equalizing it to zero. Moreover, it can be observed in (11) that, when the service level tends to 1, unfilled demand tends to 0, and the optimal order size converges to Expression (10).

$$Q_i^W=\sqrt{{\displaystyle \frac{2O{C}_i{D}_i}{H{C}_i}}}$$

(10)

$$Q_i=\sqrt{{\displaystyle \frac{2D_i\left(O{C}_i+P{C}_{i}SO{D}_i\left(\delta \right)\right)}{H{C}_i}}}=\sqrt{{\displaystyle \frac{2D_{i}O{C}_i}{H{C}_i}}}\sqrt{1+{\displaystyle \frac{P{C}_i}{O{C}_i}}SO{D}_i\left(\delta \right)}$$

(11)

Note that Expression (11) may represent a sort of correction to the traditional Wilson model, where penalty costs between two consecutive orders may be understood as another fixed order cost. Then, the order size computation just corrects the order costs $O{C}_i$ in Expression (10) by $O{C}_i+P{C}_i·SO{D}_i$, as observed in Expression (11).

Although it is natural to require $Q_i>0$, the derivation of optimal order sizes when served demand is non-negative naturally yields a solution that is automatically non-negative, as can be observed in (10) and (11). Otherwise, it is impossible to serve positive demands with negative order sizes. Thus, discarding negative roots in the process is enough to obtain positive order sizes for each existing warehouse with positive demands, and then $Q_i>0$ is not established as a formal constraint of the problem.

3. Solution Approach

The proposed formulation is solved by searching a critical point, i.e., a solution for which the total cost gradient is zero. In other words, the algorithm searches for a solution $\mathbf{x}$ such that $\nabla g\left(\mathbf{x}\right)=\mathbf{0}$, where $g\left(\mathbf{x}\right)$ is the total cost function defined in Equation (9), and $\mathbf{x}$ is the decision variable vector of the problem. Then, the equation $\nabla g\left(\mathbf{x}\right)=\mathbf{0}$ is solved by means of a root-finding problem employing the well-known multivariate Newton–Raphson method [53]. Note that all the involved expressions are explicitly developed for normally distributed demands.

The multivariate Newton–Raphson method allows us to find the approximation of a root for a vector-valued function $\mathbf{F}\left(\mathbf{x}\right)$[53]. This iterative method arises from the Taylor expansion for vector-valued functions around a point ${\mathbf{x}}_0$, as shown in (12), where $J_{\mathbf{F}}\left(\mathbf{x}\right)$ is the Jacobian matrix.

$$\mathbf{F}\left(\mathbf{x}\right)=\mathbf{F}\left({\mathbf{x}}_0\right)+J_{\mathbf{F}}\left({\mathbf{x}}_0\right)(\mathbf{x}-{\mathbf{x}}_0)+\mathcal{O}{(\mathbf{x}-{\mathbf{x}}_0)}^2$$

(12)

The method is based on a linear approximation, ignoring the $\mathcal{O}{\left(\mathbf{x}\right)}^2$ term. If $\mathbf{x}=\mathbf{r}$, where $\mathbf{F}\left(\mathbf{r}\right)=\mathbf{0}$, and ${\mathbf{x}}_0$ is the current guess, then

$$\mathbf{0}=\mathbf{F}\left(\mathbf{r}\right)\approx \mathbf{F}\left({\mathbf{x}}_0\right)+J_{\mathbf{F}}\left({\mathbf{x}}_0\right)(\mathbf{r}-{\mathbf{x}}_0)$$

(13)

or

$$\mathbf{r}-{\mathbf{x}}_0\approx -J_{\mathbf{F}}{\left({\mathbf{x}}_0\right)}^{-1}\mathbf{F}\left({\mathbf{x}}_0\right).$$

(14)

Then, with an initial guess ${\mathbf{x}}_0$, at each iteration k, the new root approximation is

$${\mathbf{x}}_{k+1}\approx {\mathbf{x}}_k-J_{\mathbf{F}}{\left({\mathbf{x}}_k\right)}^{-1}\mathbf{F}\left({\mathbf{x}}_k\right).$$

(15)

Let $EO{Q}_i$, $S{S}_i$, and $SO{D}_i$ be the order size-related costs (i.e., cycle and ordering), safety stock, and expected demand during stock-out events for the i-th warehouse, respectively. Then, for $i\in \{1,\dots ,n\}$:

$$\begin{array}{ccc}\hfill EO{Q}_i& =& {\displaystyle \frac{O{C}_i·D_i}{Q_i}}+{\displaystyle \frac{H{C}_i·Q_i}{2}},\hfill \end{array}$$

(16)

$$\begin{array}{ccc}\hfill S{S}_i& =& \sqrt{V_i·L{T}_i}\left({\displaystyle \frac{e^{-\frac{1}{2}{Z}_{\delta }^2}}{\sqrt{2\pi }}}+Z_{\delta }·\delta \right),\hfill \end{array}$$

(17)

$$\begin{array}{ccc}\hfill SO{D}_i& =& \sqrt{V_i·L{T}_i}\left({\displaystyle \frac{e^{-\frac{1}{2}{Z}_{\delta }^2}}{\sqrt{2\pi }}}-Z_{\delta }(1-\delta )\right),\hfill \end{array}$$

(18)

and, considering $\mathbf{x}={({\mathbf{x}}_1,{\mathbf{x}}_2,...,{\mathbf{x}}_N,{\mathbf{x}}_{N+1})}^{\top }={(Q_1,Q_2,...,Q_N,\delta )}^{\top }\in {\mathbb{R}}^{N+1}$ as the vector of the decision variables to be determined, we define:

$$g\left(\mathbf{x}\right)=\sum _{i=1}^N\left(EO{Q}_i+H{C}_i·S{S}_i+P{C}_i·SO{D}_i·{\displaystyle \frac{D_i}{Q_i}}\right)$$

(19)

where $\delta =1-\alpha$ is the system service level. Function $f\left(\mathbf{x}\right)$ is considered the gradient vector of $g\left(\mathbf{x}\right)$:

$$\begin{array}{cc}\hfill f\left(\mathbf{x}\right)& =\nabla g\left(\mathbf{x}\right)\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill & ={\left[\begin{array}{cccccc}{\displaystyle \frac{\partial }{\partial Q_1}}\left(g\left(\mathbf{x}\right)\right)& \cdots & {\displaystyle \frac{\partial }{\partial Q_i}}\left(g\left(\mathbf{x}\right)\right)& \cdots & {\displaystyle \frac{\partial }{\partial Q_N}}\left(g\left(\mathbf{x}\right)\right)& {\displaystyle \frac{\partial }{\partial \delta }}\left(g\left(\mathbf{x}\right)\right)\end{array}\right]}^{\top }\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill & ={\left[\begin{array}{cccccc}{f}_1\left(\mathbf{x}\right)& \cdots & f_i\left(\mathbf{x}\right)& \cdots & f_N\left(\mathbf{x}\right)& f_{\delta }\left(\mathbf{x}\right)\end{array}\right]}^{\top }\hfill \end{array}$$

(22)

where:

$$\begin{array}{cc}\hfill f_i\left(\mathbf{x}\right)& ={\displaystyle \frac{\partial }{\partial Q_i}}\left(g\left(\mathbf{x}\right)\right)=-{\displaystyle \frac{O{C}_i·D_i}{Q_i^2}}+{\displaystyle \frac{H{C}_i}{2}}-P{C}_i·SO{D}_i·{\displaystyle \frac{D_i}{Q_i^2}}\mathrm{for}i=1,\dots ,N\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill f_{\delta }\left(\mathbf{x}\right)& ={\displaystyle \frac{\partial }{\partial \delta }}\left(g\left(\mathbf{x}\right)\right)=\sqrt{2\pi }{e}^{{\displaystyle \frac{Z_{\delta }^2}{2}}}\left(\sum _{i=1}^N\sqrt{V_i·L{T}_i}\left(H{C}_i·\delta -P{C}_i·{\displaystyle \frac{D_i}{Q_i}}(1-\delta )\right)\right)\hfill \end{array}$$

(24)

We need to compute the optimal value of $\mathbf{x}$, i.e., the fixed point, where the function $g\left(\mathbf{x}\right)$ is minimized, which is equivalent to computing ${\mathbf{x}}^{*}={(Q_1^{*},Q_2^{*},\dots ,Q_N^{*},{\delta }^{*})}^{\top }$, such as $f\left({\mathbf{x}}^{*}\right)=\nabla g\left({\mathbf{x}}^{*}\right)=\mathbf{0}$.

The Newton–Raphson method will be used for addressing this computation. Let ${\mathbf{x}}^{\left(0\right)}$ be an initial guess; then, at the k-th iteration, the current solution ${\mathbf{x}}^{\left(k\right)}$ is given by:

$${\mathbf{x}}^{\left(k\right)}={\mathbf{x}}^{(k-1)}-J_f{\left({\mathbf{x}}^{(k-1)}\right)}^{-1}f\left({\mathbf{x}}^{(k-1)}\right),$$

(25)

where

$$\begin{array}{cc}\hfill J_f\left(\mathbf{x}\right)& =\nabla f\left(\mathbf{x}\right)=\left[\begin{array}{ccccc}{\displaystyle \frac{\partial f\left(\mathbf{x}\right)}{\partial Q_1}}& {\displaystyle \frac{\partial f\left(\mathbf{x}\right)}{\partial Q_2}}& \cdots & {\displaystyle \frac{\partial f\left(\mathbf{x}\right)}{\partial Q_N}}& {\displaystyle \frac{\partial f\left(\mathbf{x}\right)}{\partial \delta }}\end{array}\right]\hfill \\ \hfill & =\left[\begin{array}{ccccc}{\displaystyle \frac{\partial f_1\left(\mathbf{x}\right)}{\partial Q_1}}& {\displaystyle \frac{\partial f_1\left(\mathbf{x}\right)}{\partial Q_2}}& \cdots & {\displaystyle \frac{\partial f_1\left(\mathbf{x}\right)}{\partial Q_N}}& {\displaystyle \frac{\partial f_1\left(\mathbf{x}\right)}{\partial \delta }}\\ {\displaystyle \frac{\partial f_2\left(\mathbf{x}\right)}{\partial Q_1}}& {\displaystyle \frac{\partial f_2\left(\mathbf{x}\right)}{\partial Q_2}}& \cdots & {\displaystyle \frac{\partial f_2\left(\mathbf{x}\right)}{\partial Q_N}}& {\displaystyle \frac{\partial f_2\left(\mathbf{x}\right)}{\partial \delta }}\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ {\displaystyle \frac{\partial f_N\left(\mathbf{x}x\right)}{\partial Q_1}}& {\displaystyle \frac{\partial f_N\left(\mathbf{x}\right)}{\partial Q_2}}& \cdots & {\displaystyle \frac{\partial f_N\left(\mathbf{x}\right)}{\partial Q_N}}& {\displaystyle \frac{\partial f_N\left(\mathbf{x}\right)}{\partial \delta }}\\ {\displaystyle \frac{\partial f_{\delta }\left(\mathbf{x}\right)}{\partial Q_1}}& {\displaystyle \frac{\partial f_{\delta }\left(\mathbf{x}\right)}{\partial Q_2}}& \cdots & {\displaystyle \frac{\partial f_{\delta }\left(\mathbf{x}\right)}{\partial Q_N}}& {\displaystyle \frac{\partial f_{\delta }\left(\mathbf{x}\right)}{\partial \delta }}\end{array}\right]\hfill \\ \hfill & =\left[\begin{array}{ccccccc}{\displaystyle \frac{{\partial }^{2}g(Q,\delta )}{\partial Q_1^2}}& 0& & \cdots & & 0& {\displaystyle \frac{{\partial }^{2}g(Q,\delta )}{\partial \delta \partial Q_1}}\\ 0& {\displaystyle \frac{{\partial }^{2}g(Q,\delta )}{\partial Q_2^2}}& \ddots & & & & {\displaystyle \frac{{\partial }^{2}g(Q,\delta )}{\partial \delta \partial Q_2}}\\ & \ddots & \ddots & \ddots & & ⋮& ⋮\\ ⋮& & \ddots & \ddots & \ddots & & ⋮\\ & & & \ddots & \ddots & 0& ⋮\\ 0& & \cdots & & 0& {\displaystyle \frac{{\partial }^{2}g(Q,\delta )}{\partial Q_N^2}}& {\displaystyle \frac{{\partial }^{2}g(Q,\delta )}{\partial \delta \partial Q_N}}\\ {\displaystyle \frac{{\partial }^{2}g(Q,\delta )}{\partial Q_1\partial \delta }}& {\displaystyle \frac{{\partial }^{2}g(Q,\delta )}{\partial Q_2\partial \delta }}& \cdots & \cdots & \cdots & {\displaystyle \frac{{\partial }^{2}g(Q,\delta )}{\partial Q_N\partial \delta }}& {\displaystyle \frac{{\partial }^{2}g(Q,\delta )}{\partial {\delta }^2}}\end{array}\right]=H_g\left(\mathbf{x}\right),\hfill \end{array}$$

and where $H_g\left(\mathbf{x}\right)$ is the Hessian matrix of $g\left(\mathbf{x}\right)$, and

$$\begin{array}{cc}\hfill {\displaystyle \frac{{\partial }^{2}g(Q,\delta )}{\partial Q_i^2}}& ={\displaystyle \frac{2D_i}{Q_i^3}}\left(O{C}_i+P{C}_i\sqrt{V_{i}L{T}_i}\left({\displaystyle \frac{e^{-\frac{1}{2}{Z}_{\delta }^2}}{\sqrt{2\pi }}}-Z_{\delta }(1-\delta )\right)\right)\mathrm{for}i=1,\dots ,N\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill {\displaystyle \frac{{\partial }^{2}g(Q,\delta )}{\partial Q_i\partial Q_j}}& =0\mathrm{for}i,j=1,\dots ,N\mathrm{such}\mathrm{as}i\ne j\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill {\displaystyle \frac{{\partial }^{2}g(Q,\delta )}{\partial \delta \partial Q_i}}& ={\displaystyle \frac{{\partial }^{2}g(Q,\delta )}{\partial Q_i\partial \delta }}=\sqrt{2\pi }{e}^{{\displaystyle \frac{1}{2}}{Z}_{\delta }^2}\sqrt{V_i·L{T}_i}·{\displaystyle \frac{P{C}_i·D_i}{Q_i^2}}(1-\delta )\mathrm{for}i=1,\dots ,N\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill {\displaystyle \frac{{\partial }^{2}g(Q,\delta )}{\partial {\delta }^2}}& =e_{\delta }\sum _{i=1}^N\left(\sqrt{V_i·L{T}_i}\left(H{C}_i+{\gamma }_i+Z_{\delta }{e}_{\delta }\left(H{C}_i·\delta -{\gamma }_i(1-\delta )\right)\right)\right),\hfill \end{array}$$

(29)

with $e_{\delta }=\sqrt{2\pi }·e^{{\displaystyle \frac{1}{2}}{Z}_{\delta }^2}$ and ${\gamma }_i={\displaystyle \frac{P{C}_i·D_i}{Q_i}}$.

It is important to note that, to fulfill constraint (5), the following check condition is added to the updating step at the k-th iteration:

$${\delta }^{\left(k\right)}=\left\{\begin{array}{cc}{\delta }_{\min }\hfill & \mathrm{if}{\mathbf{x}}_{N+1}^{\left(k\right)}<{\delta }_{\min }\hfill \\ {\delta }_{\max }\hfill & \mathrm{if}{\mathbf{x}}_{N+1}^{\left(k\right)}>{\delta }_{\max }\hfill \\ {\mathbf{x}}_{N+1}^{\left(k\right)}\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(30)

In our implementation, we consider ${\delta }_{\min }=0.5$ and ${\delta }_{\max }=1-{10}^{-m}$, with $m=4$, since it is assumed that $Z_{\delta }\ge 0$ and since $Z_1$ cannot be computed. Naturally, ${\delta }_{\max }$ may consider more significant figures ($m>4$) up to user preferences.

The whole procedure is presented in Algorithm 1. In line 8, Equation (25) is computed in two steps. The first step is computing $\Delta \mathbf{x}=-J_f{\left(\mathbf{x}\right)}^{-1}f\left(\mathbf{x}\right)$ by solving the system $J_f\left(\mathbf{x}\right)\Delta \mathbf{x}=f\left(\mathbf{x}\right)$. The mentioned system is solved by the routine solve of the module numpy.linalg [54]. The second step simply involves updating the current solution by computing $\mathbf{x}=\mathbf{x}+\Delta \mathbf{x}$.

Algorithm 1 Application of Newton–Raphson Method

1: procedure levelService(${\mathbf{Q}}^{\left(0\right)},{\delta }^{\left(0\right)},N,{\delta }_{\min },{\delta }_{\max },\mathtt{tol},{\mathtt{iter}}_{\mathtt{max}}$,params)   2:     $\mathbf{x}\leftarrow \left({\mathbf{Q}}_1^{\left(0\right)},\dots ,{\mathbf{Q}}_N^{\left(0\right)},{\delta }^{\left(0\right)}\right)$   3:     $k\leftarrow 1$   4:     while $k\le {\mathrm{iter}}_{\max }$ do   5:         $\mathrm{error}\leftarrow ∥f\left(\mathbf{x}\right)∥$   6:         if $\mathrm{error}<\mathrm{tol}$ then   7:            break;   8:         $\Delta \mathbf{x}\leftarrow \mathrm{solve}\mathrm{for}-J_f\left(\mathbf{x}\right)\Delta \mathbf{x}=f\left(\mathbf{x}\right)$   9:         $\mathbf{x}\leftarrow \mathbf{x}+\Delta \mathbf{x}$                ▹ Update solution. 10:         $Q^{\left(k\right)}\leftarrow ({\mathbf{x}}_1,{\mathbf{x}}_2,\dots ,{\mathbf{x}}_N)$ 11:         ${\delta }^{\left(k\right)}\leftarrow {\mathbf{x}}_{N+1}$ 12:         if ${\delta }^{\left(k\right)}>{\delta }_{\max }$ then 13:            ${\delta }^{\left(k\right)}\leftarrow {\delta }_{\max }$ 14:            Compute $Q^{\left(k\right)}$ with Equations (18) and (11). 15:         if ${\delta }^{\left(k\right)}<{\delta }_{\min }$ then 16:            ${\delta }^{\left(k\right)}\leftarrow {\delta }_{\min }$ 17:            Compute $Q^{\left(k\right)}$ with Equations (18) and (11). 18:         $k\leftarrow k+1$ 19:     return $\mathbf{x}$

The optimal order size of the Wilson model shown in Equation (10) (see [52]) is considered the initial order size ${\mathbf{Q}}^{\left(0\right)}$ for all warehouses. The initial service level ${\delta }^{\left(0\right)}$ is set to $0.95$, the number of warehouses N depends on the instance, the tolerance of the method tol is set to ${10}^{-6}$ and compared with the norm of the gradient vector $\nabla g\left(\mathbf{x}\right)$, and the number of maximum iterations ${\mathtt{iter}}_{\mathtt{max}}$ is set to 10. The input parameter params contains the demand $D_i$, variance of the demand $V_i$, lead time $L{D}_i$, unit holding cost $H{C}_i$, penalty cost $P{C}_i$, and fixed-order cost $O{C}_i$ for each warehouse i (for more details, see Section 2.2).

4. Computational Experimentation and Numerical Results

In all of these experiments, a service level of ${\delta }_{base}$ is considered as a benchmark solution for the algorithm. In particular, the values considered are computed based on [38], which relies on a previous approximation for optimizing the system service level that assumes an approximated safety stock expression and order sizes computed as in the Wilson model’s solution.

Table 1. Results with $N=5$ warehouses.

CV	PC	${\mathit{\delta }}_{\mathit{base}}$	Base Tot.	${\mathit{\delta }}_{\mathit{op}}$	Op. Tot.	Diff.	%Diff.	t (sec.)
$0.1$	10	$0.637$	$6610.853$	$0.729$	$6593.118$	$17.736$	$0.268$	$0.049$
	50	$0.927$	$6824.682$	$0.931$	$6822.696$	$1.987$	$0.029$	$0.037$
	100	$0.964$	$6914.253$	$0.964$	$6912.805$	$1.447$	$0.021$	$0.039$
$0.2$	10	$0.637$	$7030.687$	$0.724$	$6991.544$	$39.142$	$0.557$	$0.041$
	50	$0.927$	$7458.345$	$0.929$	$7451.254$	$7.091$	$0.095$	$0.039$
	100	$0.964$	$7637.485$	$0.963$	$7631.883$	$5.602$	$0.073$	$0.039$
$0.3$	10	$0.637$	$7450.520$	$0.719$	$7386.327$	$64.194$	$0.862$	$0.047$
	50	$0.927$	$8092.008$	$0.928$	$8076.688$	$15.319$	$0.189$	$0.039$
	100	$0.964$	$8360.718$	$0.963$	$8348.246$	$12.473$	$0.149$	$0.038$

,

Table 2. Results with $N=10$ warehouses.

CV	PC	${\mathit{\delta }}_{\mathit{base}}$	Base Tot.	${\mathit{\delta }}_{\mathit{op}}$	Op. Tot.	Diff.	%Diff.	t (sec.)
$0.1$	10	$0.643$	$13,227.495$	$0.732$	$13,194.073$	$33.423$	$0.253$	$0.051$
	50	$0.929$	$13,646.542$	$0.932$	$13,642.784$	$3.758$	$0.028$	$0.038$
	100	$0.964$	$13,821.493$	$0.965$	$13,818.755$	$2.738$	$0.020$	$0.040$
$0.2$	10	$0.643$	$14,052.249$	$0.727$	$13,978.429$	$73.820$	$0.525$	$0.042$
	50	$0.929$	$14,890.343$	$0.930$	$14,876.928$	$13.414$	$0.090$	$0.040$
	100	$0.964$	$15,240.244$	$0.964$	$15,229.646$	$10.598$	$0.070$	$0.041$
$0.3$	10	$0.643$	$14,877.003$	$0.723$	$14,755.858$	$121.145$	$0.814$	$0.049$
	50	$0.929$	$16,134.143$	$0.929$	$16,105.162$	$28.981$	$0.180$	$0.040$
	100	$0.964$	$16,658.995$	$0.963$	$16,635.402$	$23.593$	$0.142$	$0.040$

,

Table 3. Results with $N=20$ warehouses.

CV	PC	${\mathit{\delta }}_{\mathit{base}}$	Base Tot.	${\mathit{\delta }}_{\mathit{op}}$	Op. Tot.	Diff.	%Diff.	t (sec.)
$0.1$	10	$0.647$	$27,086.481$	$0.734$	$27,017.934$	$68.546$	$0.253$	$0.115$
	50	$0.929$	$27,964.850$	$0.933$	$27,956.983$	$7.867$	$0.028$	$0.042$
	100	$0.965$	$28,330.820$	$0.965$	$28,325.066$	$5.755$	$0.020$	$0.043$
$0.2$	10	$0.647$	$28,818.063$	$0.729$	$28,666.230$	$151.833$	$0.527$	$0.049$
	50	$0.929$	$30,574.802$	$0.931$	$30,546.635$	$28.167$	$0.092$	$0.043$
	100	$0.965$	$31,306.743$	$0.964$	$31,284.438$	$22.304$	$0.071$	$0.047$
$0.3$	10	$0.647$	$30,549.646$	$0.724$	$30,299.868$	$249.777$	$0.818$	$0.049$
	50	$0.929$	$33,184.754$	$0.930$	$33,123.813$	$60.941$	$0.184$	$0.043$
	100	$0.965$	$34,282.665$	$0.964$	$34,232.979$	$49.686$	$0.145$	$0.043$

,

Table 4. Results with $N=50$ warehouses.

CV	PC	${\mathit{\delta }}_{\mathit{base}}$	Base Tot.	${\mathit{\delta }}_{\mathit{op}}$	Op. Tot.	Diff.	%Diff.	t (sec.)
$0.1$	10	$0.662$	$69,521.335$	$0.742$	$69,356.054$	$165.281$	$0.238$	$0.056$
	50	$0.932$	$71,844.193$	$0.935$	$71,823.688$	$20.505$	$0.029$	$0.098$
	100	$0.966$	$72,804.162$	$0.967$	$72,788.943$	$15.219$	$0.021$	$0.045$
$0.2$	10	$0.662$	$74,135.464$	$0.738$	$73,765.225$	$370.239$	$0.499$	$0.049$
	50	$0.932$	$78,781.180$	$0.934$	$78,707.043$	$74.137$	$0.094$	$0.047$
	100	$0.966$	$80,701.118$	$0.966$	$80,641.939$	$59.179$	$0.073$	$0.045$
$0.3$	10	$0.662$	$78,749.593$	$0.733$	$78,134.903$	$614.690$	$0.781$	$0.050$
	50	$0.932$	$85,718.167$	$0.932$	$85,557.140$	$161.027$	$0.188$	$0.045$
	100	$0.966$	$88,598.074$	$0.965$	$88,466.073$	$132.001$	$0.149$	$0.046$

,

Table 5. Results with $N=100$ warehouses.

CV	PC	${\mathit{\delta }}_{\mathit{base}}$	Base Tot.	${\mathit{\delta }}_{\mathit{op}}$	Op. Tot.	Diff.	%Diff.	t (sec.)
$0.1$	10	$0.652$	$139,961.242$	$0.737$	$139,610.759$	$350.483$	$0.250$	$0.062$
	50	$0.930$	$144,609.401$	$0.934$	$144,567.711$	$41.690$	$0.029$	$0.046$
	100	$0.965$	$146,539.884$	$0.966$	$146,509.175$	$30.709$	$0.021$	$0.048$
$0.2$	10	$0.652$	$149,149.699$	$0.732$	$148,369.477$	$780.223$	$0.523$	$0.052$
	50	$0.930$	$158,446.018$	$0.932$	$158,296.041$	$149.977$	$0.095$	$0.049$
	100	$0.965$	$162,306.984$	$0.965$	$162,187.764$	$119.220$	$0.073$	$0.051$
$0.3$	10	$0.652$	$158,338.157$	$0.727$	$157,049.333$	$1288.824$	$0.814$	$0.109$
	50	$0.930$	$172,282.635$	$0.931$	$171,957.523$	$325.112$	$0.189$	$0.048$
	100	$0.965$	$178,074.084$	$0.964$	$177,808.317$	$265.767$	$0.149$	$0.050$

and

Table 6. Results with $N=200$ warehouses.

CV	PC	${\mathit{\delta }}_{\mathit{base}}$	Base Tot.	${\mathit{\delta }}_{\mathit{op}}$	Op. Tot.	Diff.	%Diff.	t (sec.)
$0.1$	10	$0.656$	$270,940.412$	$0.740$	$270,291.532$	$648.880$	$0.239$	$0.142$
	50	$0.931$	$279,788.460$	$0.934$	$279,711.036$	$77.424$	$0.028$	$0.300$
	100	$0.966$	$283,455.235$	$0.966$	$283,398.193$	$57.042$	$0.020$	$0.321$
$0.2$	10	$0.656$	$288,467.066$	$0.735$	$287,021.597$	$1445.469$	$0.501$	$0.354$
	50	$0.931$	$306,163.162$	$0.933$	$305,884.568$	$278.594$	$0.091$	$0.236$
	100	$0.966$	$313,496.712$	$0.965$	$313,275.243$	$221.469$	$0.071$	$0.238$
$0.3$	10	$0.656$	$305,993.720$	$0.730$	$303,604.624$	$2389.096$	$0.781$	$0.259$
	50	$0.931$	$332,537.864$	$0.932$	$331,933.871$	$603.993$	$0.182$	$0.241$
	100	$0.966$	$343,538.189$	$0.965$	$343,044.464$	$493.726$	$0.144$	$0.464$

show the numerical results of the instances with $N=5, 10, 20, 30, 50, 100$, and 200 warehouses, respectively, representing a wide set of reasonable real-world sizes. For each warehouse i, all cases consider $H{C}_i=1$ as a baseline cost parameter, and $P{C}_i$ ranges within $\{10,50,100\}$, representing a significant range of variation with respect to the holding costs, $H{C}_i$. In addition, the coefficient of variation $CV=0.1,0.2$, and $0.3$ for each instance, focused on high-volume, smooth demands. Thus, 9 cases or sub-instances are defined for each instance (i.e., by combining $PC$ and $CV$ variations). Finally, for $i\in \{1,\dots ,N\}$, the order cost $O{C}_i$ is determined in terms of the lead time $L{T}_i$, as $O{C}_i=500+500L{T}_i$, the lead time $L{T}_i$ is randomly set in the range $[2,5]$, and the average demand $D_i$ is considered to be in the range $[250,500]$.

The columns ${\delta }_{base}$ and Base Tot. show the service level and total cost considering the Wilson order size. The columns ${\delta }_{op}$ and Op Tot. show the service level and total cost considering the order size of the presented model. The columns Diff. and %Diff. show the difference between both solutions. The last column t (sec.) shows the computational time in the execution of the presented model.

It is worth noting that the Hessian matrix of the cost function $g\left(\mathbf{x}\right)$ (i.e., the Jacobian matrix of the gradient vector $f\left(\mathbf{x}\right)$) at the end of the algorithm was positive definite in all tested cases. This is shown by the obtained eigenvalues that were strictly positive for all instances, as observed in

Table 7. Minimum (${\lambda }_{\min }$) and maximum (${\lambda }_{\max }$) eigenvalues for the Hessian matrix $H_g\left(\mathbf{x}\right)$ at the approximate solution ${\mathbf{x}}^{*}$, such as $\nabla f\left({\mathbf{x}}^{*}\right)\approx \mathbf{0}$, considering $N=5, 10$, and 20.

$\mathit{CV}$	$\mathit{PC}$	$\mathit{N}=5$		$\mathit{N}=10$		$\mathit{N}=20$
		${\mathbf{\lambda }}_{\min }$	${\mathbf{\lambda }}_{\max }$	${\mathbf{\lambda }}_{\min }$	${\mathbf{\lambda }}_{\max }$	${\mathbf{\lambda }}_{\min }$	${\mathbf{\lambda }}_{\max }$
$0.1$	10	$0.000598$	$3644.7$	$0.000569$	$7263.5$	$0.000570$	$15,378.4$
	50	$0.000599$	$35,620.4$	$0.000570$	$71,653.6$	$0.000571$	$152,555.7$
	100	$0.000600$	$116,503.6$	$0.000571$	$234,758.5$	$0.000572$	$500,358.3$
$0.2$	10	$0.000572$	$7094.3$	$0.000550$	$14,145.9$	$0.000552$	$29,937.9$
	50	$0.000574$	$68,574.7$	$0.000552$	$138,098.8$	$0.000554$	$293,919.9$
	100	$0.000577$	$224,349.8$	$0.000554$	$452,590.9$	$0.000557$	$964,330.5$
$0.3$	10	$0.000547$	$10,361.3$	$0.000532$	$20,671.6$	$0.000534$	$43,731.6$
	50	$0.000549$	$99,025.6$	$0.000535$	$199,643.9$	$0.000537$	$424,741.5$
	100	$0.000554$	$324,036.7$	$0.000538$	$654,444.6$	$0.000541$	$1,393,913.9$

and

Table 8. Minimum (${\lambda }_{\min }$) and maximum (${\lambda }_{\max }$) eigenvalues for the Hessian matrix $H_g\left(\mathbf{x}\right)$ at the approximate solution ${\mathbf{x}}^{*}$, such as $\nabla f\left({\mathbf{x}}^{*}\right)\approx \mathbf{0}$, considering $N=50, 100$, and 200.

$\mathit{CV}$	$\mathit{PC}$	$\mathit{N}=50$		$\mathit{N}=100$		$\mathit{N}=200$
		${\mathbf{\lambda }}_{\min }$	${\mathbf{\lambda }}_{\max }$	${\mathbf{\lambda }}_{\min }$	${\mathbf{\lambda }}_{\max }$	${\mathbf{\lambda }}_{\min }$	${\mathbf{\lambda }}_{\max }$
$0.1$	10	$0.000582$	$42,501.7$	$0.000575$	$82,713.1$	$0.000577$	$159,443.1$
	50	$0.000583$	$431,659.1$	$0.000577$	$827,875.5$	$0.000579$	$1,606,404.9$
	100	$0.000584$	$1,422,147.9$	$0.000578$	$2,719,990.2$	$0.000579$	$5,284,173.3$
$0.2$	10	$0.000566$	$82,647.2$	$0.000559$	$160,914.3$	$0.000562$	$310,311.9$
	50	$0.000569$	$831,004.9$	$0.000562$	$1,593,976.4$	$0.000564$	$3,095,319.7$
	100	$0.000571$	$2,738,989.9$	$0.000564$	$5,238,968.5$	$0.000566$	$10,185,844.2$
$0.3$	10	$0.000551$	$120,589.3$	$0.000543$	$234,894.4$	$0.000547$	$453,144.7$
	50	$0.000554$	$1,199,891.8$	$0.000547$	$2,301,892.7$	$0.000550$	$4,473,392.9$
	100	$0.000557$	$3,956,325.2$	$0.000550$	$7,568,015.3$	$0.000553$	$14,725,563.3$

, which show the minimum and maximum eigenvalues obtained for each instance. Therefore, it can be stated that all obtained solutions at least represent a local minimum for the studied model. In addition, merely as empirical evidence, Figure 2 shows the evolution of the total system cost for $N=5,10,20,50,100$, and 200, with $CV=0.1$ and $PC=10$, for different values of $\delta$, and considering the optimal value of $Q_i$ of expression (11) ($\forall i=1,...,N$) for each fixed value of $\delta$. These figures denote the existence of a single minimum for $\delta$ in the range $(0,1)$. In sum, although it has not been demonstrated that the studied problem has a strictly convex objective function and a single global optimum (given the associated algebraic complexity, especially considering the relationship between $Z_{\delta }$ and $\delta$), it is shown that all the obtained values effectively represent optimal solutions to the problem within the feasible domain of $\delta$.

Nevertheless, convexity may be reasonable and expected, since when the service level $\delta$ approaches 1, the safety stock cost tends to ∞ (i.e., the second term in Expression (1)), and the penalty cost tends to 0 (i.e., the third term in Expression (1)), thus yielding a total cost tending to ∞. In addition, the order size $Q_i$ ($\forall i=1,2,3,\dots ,N$) converges to the Wilson model solution in Expression (10), as observed in Equation (11). This expression at the limit is independent of $\delta$, and the related costs of the first term in (1) do not affect the limit of the total costs (which tends to ∞). On the contrary, when $\delta$ approaches 0, the safety stock cost tends to 0 and the penalty cost to ∞, again yielding a total cost approaching ∞. Moreover, the order size and related costs (the first term in Expression (1)) tends to ∞ when $\delta$ approaches 0.

Furthermore, it may be natural that the penalty cost is strictly and continuously increasing, whereas the safety stock cost is strictly and continuously decreasing, both with $\delta$. Finally, it can be observed that order size-related costs (i.e., the first term in (1)), when order size is determined by Equation (11), are strictly and continuously decreasing from ∞ to the optimal cost of the Wilson model, given by $\sqrt{2·H{C}_i·O{C}_i·D_i}$. In sum, it may be expected that the total cost in expression (1) tends to ∞ when the service level $\delta$ approaches 0 or 1, creating a convex U-shaped function, when also considering that total costs are always positive when $D_i$ and all cost parameters are positive (i.e., 0 represents a lower bound to the total costs), as shown in Figure 2.

As can be observed for all instances in Table 1, Table 2, Table 3, Table 4, Table 5 and Table 6, although the total costs for the proposed approach are slightly lower than those obtained with the benchmark approach, the system service level is significantly different when the penalty cost parameter is the lowest (i.e., $PC=10$). On the contrary, when penalty cost is the highest (i.e., $PC=100$) and the optimal service level tends to 1, the benchmark and the optimal service level tend to be the same. These results are explained by the fact that, when the service level is high, stock-out demands converge to 0 and order sizes tend to Wilson order size; thus, the benchmark solution becomes a very good approximation. In the same way, savings obtained when the proposed approach is employed are more significant when penalty cost parameter $PC$ is the lowest. Accordingly, the proposed approach is more beneficial when $PC$ is not too high and the optimal service level is not near 1. This is especially relevant from a supply chain network design perspective (as discussed in most ILP-related literature described in Section 1), where variations in the global service level may imply significant changes to the SCN configuration, such as different customer assignments or more relevant changes in warehouse locations and the number of located warehouses.

The aforementioned results are contrasted with the very low computing time required by the proposed approach, even for the largest considered instances (i.e., 200 warehouses), always being less than a second. Thus, the fact that the proposed approach has the potential to provide a better system service level, which may even impact the design of the SCN and related costs given the unsubstantial computational efforts required, supports the usage of the proposed approach instead of previous approximations.

As observed in Table 1, Table 2, Table 3, Table 4, Table 5 and Table 6, a significant similitude is observed in the optimal service level when $CV$ varies (i.e., $0.1$, $0.2$, and $0.3$), while in terms of total system costs, a more significant variation is obtained when $CV$ is the highest. For example, for 200 warehouses, if $CV=0.1$ and $PC=10$, then the total savings are about $0.24\%$, whereas when $CV=0.3$ and $PC=10$, the total savings rise to about $0.78\%$ (i.e., more than three times in savings).

In terms of the algorithm behavior, Figure 3, Figure 4 and Figure 5 show the number of iterations executed with $N=5$ warehouses and different values of $CV$ and $PC$. Each graph shows the evolution of the error $\parallel \nabla g\left(\mathbf{x}\right)\parallel$, total cost, i.e., the value of $g\left(\mathbf{x}\right)$, the service level values $\delta$, and the order size for each warehouse. These results denote the good and fast performance of the proposed algorithm, which yield optimal solutions in only 6 iterations. In particular, and in practical terms, it is observed that the solution seems to be a very good approximation of the optimal values after only 3 or 4 iterations. Analogously, these results are very similar for larger instances where, for $N=20, 50, 100$, and 200 warehouses, the algorithm converges in 5 iterations, and the solution after 3 or 4 iterations is practically the same as the final optimal solution, all considering a tolerance of $\mathtt{t}\mathtt{o}\mathtt{l}={10}^{-6}$, as depicted in Figure 6.

All experiments were run on a computer with an 11th Gen Intel Core i5 processor running at $2.6$ GHz with 8 cores using 16 MB of RAM, running Ubuntu version $22.04.4$. The programming language used was Python $3.11.5$ [55]. Although other programming language are suitable for this kind of studies, the selection of Python was based mainly on its free availability and easy use. In addition, the complexity and size of the studied problem and its instances are not so challenging as to require more specialize language programming, such as MATLAB or Mathematica (visited on 15 August 2024).

5. Discussion of Research Results and Managerial Insights

Given the proposed formulation and solution approach, in addition to the obtained results, several insights and implications from a scientific and practical perspective may be obtained. Firstly, this research represents a significant contribution to inventory location and supply chain inventory-related literature through an enhanced formulation along with an exact solution approach to optimizing a global inventory-related service level (i.e., stock-out probability) and order sizes for a multi-warehouse system, along with an exact solution approach based on the well-known Newton–Raphson method. Note that the proposed solution approach is analytically derived considering a normal distribution (i.e., an approximation) for warehouse demands. Nonetheless, the approach may be adapted to other distributions, while analytical derivation may result in a more complex or simpler one depending on the specific distribution considered (e.g., uniform, triangular, log-normal, truncated normal, Poisson, gamma, beta, etc.). In particular, a numerical approximation instead of an analytical derivation of the key algorithm elements (i.e., gradient vector and Jacobian matrix) may be developed considering empirical probabilistic distributions, as described in [48].

It is worth highlighting that, although normal distribution presents a positive probability of observing negative values, it may be insignificant for aggregated high-volume demands with small levels of relative variability (i.e., Coefficient of Variation, $CV$). For example, a normal demand with a positive mean and a $CV$ of $0.5$ implies a negative demand probability of $0.02275$, whereas for all experiments carried out in this research, a maximum $CV$ of $0.3$ is considered, with a negative demand probability of $0.00043$.

Table 9. Summary of normal demand simulation with $\mu =1$.

$\mathit{CV}$	$0.1$	$0.2$	$0.3$	$0.4$	$0.5$
Mean (%)	$0.0000$	$0.0000$	$0.0440$	$0.6710$	$2.2910$
Standard Deviation	$0.0000$	$0.0000$	$0.0641$	$0.2500$	$0.4667$

presents the average percentage of times where negative demand values are obtained for 100 experiments, where each one involves 1000 random demand values following a normal distribution with a $CV=0.1, 0.2, 0.3, 0.4,$ and $0.5$ (note that this experiment is independent of mean demand value, except if it is 0 or negative). It can be observed that, for a $CV$ of $0.1$ and $0.2$, negative demands are not observed, whereas for a $CV$ of $0.3$, only in $0.044\%$ of cases was a negative value observed. In the extreme case of considering a $CV$ of $0.5$, negative demands are observed in $2.29\%$ of cases. Accordingly, potential errors related to the occurrence of negative demands may be considered insignificant, and the proposed approach represents a powerful methodology for practical implementations.

In terms of the modeling structure, it is worth highlighting that previous research within inventory location-related literature fail in modeling in an exact manner both safety stock and unfilled demand costs, where the associated stock-out probability is not considered as a decision variable but as an input parameter. In our research, in contrast, the service level is considered a system decision variable to be jointly optimized, which makes the solution approach difficult. Note that if a different service level for each warehouse is assumed, then the problem may be decoupled into independent problems for each warehouse with two decision variables (i.e., the service level and the order size). On the contrary, as considered in this research, if the required service level is assumed to be the same for all warehouses, then the common system service level and the order sizes for all warehouses must be optimized simultaneously.

Based on the obtained results, it is observed that the proposed model allows for obtaining solutions with a better system performance in contrast to a natural benchmark, where the expected safety stock costs are modeled in the traditional manner observed in inventory management and inventory location literature. This simplified modeling approach yields a different suboptimal service level, slightly different order sizes for each warehouse, and a more expensive system cost compared to the proposed approach. Naturally, if a pre-fixed service level is considered near the optimal value obtained with the proposed model, then the improvements may be immaterial. On the contrary, if a pre-fixed service level is far from the optimal value obtained with our approach, then the improvements may be even more significant.

It is remarkable that the proposed model and solution approach do not imply a significant increase in computational effort, since optimal solutions are obtained after only a few iterations with very low computational times, even for the largest instances considered (i.e., 200 warehouses). Note that observing more than 10 or 20 warehouses for a single company in the real world is very infrequent, except for some global multinational companies.

Moreover, the proposed model and solution approach are particularly relevant as part of a sub-problem within further inventory location problems. This is especially significant when the generalized Benders decomposition is employed, as in [36,37,40], where the full original problem is decomposed into an MIP formulation that deals with warehouse location and customer assignment decisions (i.e., within a master problem), whereas inventory decisions are addressed by an underlying sub-problem (such as the problem addressed in this research). It is remarkable that none of the aforementioned works optimizes the service level, which is considered a fixed parameter. Accordingly, a natural future research direction consists of extending the proposed methodology to address an inventory location problem to design a supply chain network while simultaneously optimizing decisions of warehouse location, customer assignment, order sizes, and the global inventory service level, which enlarge the contributions of this research.

In terms of managerial insights, it worth highlighting that the proposed approach may definitively benefit supply chain managers by providing better solutions in terms of strategic supply chain inventory performance. These benefits may imply significant savings in the long terms. In particular, the proposed methodology (i.e., the model and the solution algorithm) may help in efficiently supporting the decision-making process related to optimally setting significant inventory decision variables in real-world supply chains, such as the global service level and warehouse order sizes. Moreover, very low computational complexity and efforts to implement the proposed approach are required, significantly facilitating its adoption and utilization in real-world settings. Finally, in relation to the existence of third-party logistic providers, they may benefit by reducing emergency shipments, given the optimal setting of the global service level and related safety stock, thus yielding reductions in related costs and improvements in terms of system sustainability and service level.

Additionally, the methodology may be used in small regions where the homogeneity of customers leads to the marketing goal of offering a homogeneous customer service level at the warehouses, since customers may perceive a service level difference or, in small regions, their communication may allow them to share this knowledge. Therefore, it is important to the company to offer the same service at different locations, since the proposed approach allows for optimizing the associated costs.

6. Conclusions and Future Research

In this research, a model to optimize inventory-related decisions for a single-commodity supply chain network comprising a set of parallel warehouses or distribution centers is developed. Inventory decisions in the model are warehouse order sizes and a common global service level consisting of a homogeneous maximum allowed stock-out probability for each warehouse, assuming stochastic demands with a generic formulation. Subsequently, a second specific formulation is derived considering a normal approximation of warehouse demands. Finally, a Newton–Raphson-based algorithm is developed to solve the problem, converging in only a few iterations and seconds for instances with up to 200 parallel warehouses.

The proposed formulation along with the developed solution approach yield significant benefits in comparison with a natural benchmark, without implying substantial additional computational efforts, thus providing significant differences in the obtained service level and lower system costs, which are contributions associated with the proposed approach.

It is worth noting that the studied problem and related variations may arise as a sub-problem within a generalized Benders decomposition for inventory location problems to optimize the supply chain network. Similarly, any two-step local search heuristic or metaheuristics may benefit from employing the proposed approach as a procedure to evaluate every SCN configuration (e.g., fitness evaluation within genetic algorithms). The above arguments enlarge the applicability and relevance of this research.

Relevant future research based on this research may increase its contributions, which involve multi-commodity and multi-period formulations, along with considering other demand distributions. Moreover, empirical distributions may be considered in which expressions involved within the algorithm are derived based on numerical approximations. A simultaneous modeling of decisions and costs at the upstream location level, such as plants or suppliers, denotes a significant extension to be addressed, following most of the literature concerning two-stage supply chain inventory management, as described in the literature review above.