Simulation-Based Optimization of a Two-Echelon Continuous Review Inventory Model with Lot Size-Dependent Lead Time

Abstract

This study analyzes a stochastic continuous review inventory system (Q,r) using a simulation-based optimization model. The lead time depends on lot size, unit production time, setup time, and a shop floor factor that represents moving, waiting, and lot size inspection times. A simulation-based model is proposed for optimizing order quantity (Q) and reorder point (r) that minimize the total inventory costs (holding, backlogging, and ordering costs) in a two-echelon supply chain, which consists of two identical retailers, a distributor, and a supplier. The simulation model is created with Arena software and validated using an analytical model. The model is interfaced with the OptQuest optimization tool, which is embedded in the Arena software, to search for the least cost lot sizes and reorder points. The proposed model is designed for general demand distributions that are too complex to be solved analytically. Hence, for the first time, the present study considers the stochastic inventory continuous review policy (Q,r) in a two-echelon supply chain system with lot size-dependent lead time L(Q). An experimental study is conducted, and results are provided to assess the developed model. Results show that the optimized Q and r for different distributions of daily demand are not the same even if the associated total inventory costs are close to each other.

1. Introduction

Inventory management plays a vital role in supply chains. It is one of the activities associated with managing customer and vendor relationships and plays a major role in facilitating the balancing of demand and supply. Effective inventory management in supply chains can lead to an increase in profits by cutting inventory-holding costs, and this indicates the need for a general model to manage inventories within a supply chain.

In classical continuous inventory review models ($Q,r$), demand and lead time are dealt with as constant input parameters or as random variables that follow a specific probability distribution. This could lead to an assumption that lead time is not subject to control and given as an input to the model. In practice, this assumption is not realistic since lead time can be controlled through some of its components, such as lot size, unit production cost, production time, and preparation time (setup, moving, and waiting times), as suggested by Kim and Benton [1]. In addition, considering lead time as a controllable factor allows for lead time reductions; this can lead to better customer satisfaction by lowering the risk of shortages and reducing safety stock [2]. Therefore, in many real-world systems, inventory models are stochastic and lead times are dependent and controllable. Several solution methodologies have been proposed to solve stochastic continuous inventory review models ($Q,r$) with a dependent lead time. These methodologies are classified into two main types: analytical models and simulation-based optimization models.

One of the first approximation models developed to study the impact of order quantities in a continuous review inventory policy ($Q,r$) on lead times was proposed by Kim and Benton [1]. The authors assumed that demand is normally distributed, order quantity $Q$ is fixed, and the production lead time varies linearly with the order quantity. In addition, queue-waiting times were assumed to be a fixed fraction of the lead time. Based on an adjusted economic order quantity (EOQ), they developed an iterative algorithm to find the near-optimal order quantity ($Q$) and safety stock ($s$). Results showed that significant savings could be achieved if the interrelationship between order quantity and lead time is taken into account. Hariga [3] suggested a modification to the model by adjusting an expression related to annual backorder cost, in addition to proposing another relationship for the revised lot size. Both models assumed a constant processing time per unit. The model in [4] assumed random processing times and the processing time of a lot was correlated with its size. An exact model was developed, and two cases of lead time dependency were discussed. In the first, the lead time depends on lot size, and in the second, on the reserved capacity, whereby a customer can secure shorter processing times by reserving capacity at the supplier’s manufacturing facility. Ben-Daya and Hariga [5] proposed a model for an integrated production inventory problem for a vendor–buyer inventory system. A continuous review inventory policy was considered. Demand was assumed to be normally distributed and the lead time varied linearly with lot size. The lead time components were lot size-dependent run time and constant delay times, such as moving, waiting, and setup times. Moon and Cha [6] studied lead time reduction by controlling the manufacturer’s regular production rate and taking the risk of paying an additional cost into account. The relationship between the lot size and the production rate of the manufacturer expressed the lead time. An analytical model was proposed to optimize a continuous review inventory policy, based on the assumption that demand is normally distributed. Song et al. [7] proposed an integrated production inventory model for a supplier–buyer inventory system in which the buyer adopted a ($Q,r$) continuous review inventory policy. Demand was assumed to be normally distributed, while the supplier and buyer jointly controlled the lead time. The supplier determines the production rate while the buyer controls the lot size (order quantity). Song et al. [8] presented a Stackelberg game framework to model the interactions between a manufacturer and a retailer. Lead time and demand were distribution-free with only two known parameters: mean and variance. An iterative algorithm was developed to solve the proposed model. Sarkar et al. [9] proposed an integrated vendor–buyer model with shortages. Demand was assumed to be normally distributed and lead time was assumed to be a stochastic variable that depends on the order size and production rate. Sarkar and Giri [10] proposed an algorithm for finding the optimal solution of an integrated inventory model for a single-vendor single-buyer. The model assumed that demand follows a normal distribution and lead time is a linear function of batch size, setup time, and transport time.

Most of the proposed models were approximation models with the assumption that daily demand is normally distributed. However, it is extremely difficult to analytically solve the constrained optimization model that was developed by Hadley and Whitin [11]. Due to these limitations and difficulties, simulation-based optimization models are utilized to model and analyze complex inventory systems in order to obtain more practical results without the need to make restrictive assumptions. Al-harkan and Hariga [2] proposed a simulation-based optimization model for a single-supplier and single-retailer continuous review inventory system assuming a lead time dependent on lot size. The model was implemented for different daily demand distributions (normal, gamma, and log-normal). In addition, Tsai and Chen [12] proposed a simulation-based optimization model for a multi-objective single-supplier and single-retailer inventory system. The proposed model aims to find the values of reorder point and order quantity that simultaneously minimize the total inventory cost, average inventory level, and frequency of inventory shortages. The demand and lead time were assumed to be random and follow a uniform distribution. Avci and Selim [13] proposed a multi-objective simulation-based optimization approach for a three-echelon convergent supply chain system consisting of a customer, a manufacturer, and a number of suppliers. Kleijnen and Wan [14] studied an ($s,S$) inventory system to optimize reorder level $s$ and maximum inventory level $S$. A simulation-based optimization model was developed using Arena and $OptQuest$. The metric to be minimized was the expected inventory costs, which include holding, ordering, and shortage costs while satisfying service (fill rate) constraints. Simulation-based optimization methods in stochastic inventory systems have been used to minimize total inventory costs while satisfying service level constraints as in [15,16,17,18]. A comprehensive survey of simulation-based optimization approaches to solve inventory replenishment problems is presented in [19].

Based on the literature review, we found that previous research solved the single-supplier and single-retailer continuous review inventory system. Most of the proposed models assumed that daily demand is normally distributed. Only studies that utilized simulation-based optimization models considered other distributions for daily demand. The lead time was considered as a variable and depended on several controllable factors in most previous studies. All the proposed models (analytical or simulation) aimed to minimize inventory costs.

In this paper, we solve a two-echelon three-stage supply chain stochastic continuous review inventory system with lot size-dependent lead time. The supply chain system has a supplier in the higher stage, a distributor in the middle stage, and two retailers in the lower stage. Different distributions for daily demand were considered. A simulation-based optimization model was developed to find near-optimal values of reorder point and order quantity that minimize the total inventory cost (holding, backlogging, and ordering costs). Hence, an efficient simulation-based optimization model is proposed for the first time for the problem considered herein.

The following section addresses the methodology. Section 3 presents the results and discussions. Finally, the conclusions are stated in Section 4.

2. Methodology

2.1. Problem Formulation

This study considers a three-stage two-echelon inventory system consisting of a supplier at the higher stage, a distributor at the middle stage, and two retailers at the lower stage, as shown in Figure 1. The main aim of this study is to find the optimal or near-optimal order quantities, as well as reorder points for the distributor and the two retailers that minimize the total inventory costs. The inventory costs include ordering costs, holding costs, and shortage costs in the three-stage supply chain.

A stochastic continuous review inventory control policy ($Q,r$) is employed, which means that retailer 1, retailer 2, or the distributor will place an order whenever the inventory position drops below the reorder point $r$; the order size is $Q$. An Arena model is developed to simulate the policy employed in the three-stage supply chain for ten years (2500 days). The optimization study was conducted using the developed Arena model, as well as the Arena $OptQuest$ optimization tool. Assuming that the demand follows different probability distributions, several scenarios are considered in order to study the performance of the simulation model. Simulated results were validated with results obtained using an analytical method developed in [1]. The following assumptions are considered:

• A customer’s demand is satisfied immediately if the on-hand inventory of the retailer is sufficient. Otherwise, it is considered as a backorder and will be filled by future deliveries.
• A retailer’s order is satisfied immediately if the on-hand inventory of the distributor is sufficient. Otherwise, it is also considered as a backorder and will be filled by future deliveries.
• There is always a sufficient number of items in the supplier’s storage.
• The lead time depends on the order size $Q$; this is explained below.

The replenishment order will be received after a period of time, referred to as the lead time, which depends on the order size ($Q$). Based on the research conducted by Kim and Benton [1], a linear relationship between lead time and lot size can be assumed, and can be expressed mathematically as follows:

$$L\left(Q\right)=\left(\theta +P· Q\right)\delta$$

(1)

where $\theta$ is the setup time (in days); $P$ is the unit production time in days, and $\delta$ is the shop floor queuing factor ($\delta >1)$. We find that the lead time also depends on production time (represented by $P$) and other times related to the order preparation, such as setup (represented by $\theta$), moving, and waiting times (represented by the shop floor factor $\delta$).

In the system described above, the problem is determining the optimal or near-optimal order sizes and reorder points that minimize the total inventory cost in the system shown in Figure 1 under the given assumptions. Our objective function (the total inventory cost) can be represented as the sum of the holding, backorders, and ordering costs. Generally, it is difficult to create an exact formulation of the total inventory cost for a stochastic continuous review inventory system, and several approximate formulations have been proposed in the literature. One of these approximation treatments is the constrained model developed in [11]; the reader can refer to the original publication for more information. The authors reported that it is extremely difficult to analytically solve the constrained optimization model, even with a normally distributed daily demand.

Therefore, we propose a simulation-based-optimization methodology that can be utilized for daily demand distributions of any type. A time-persistent model is developed to calculate the simulated total inventory cost for a period of $n$ days; it is formulated with the following variables:

• $n$ Number of days;
• $n_{Year}$ Number of days per year;
• $t$ Point in time;
• $Ri$ Retailer $i$, where $i=1,2$;
• $DC$ Distributor;
• $I\left(t\right)$ Inventory level at time $t$;
• $I^{+}\left(t\right)$ Physical inventory level at time $t$, $I^{+}\left(t\right)=\max \left\{0, I\left(t\right)\right\}$;
• $I^{-}\left(t\right)$ Backlog amount at time $t$, $I^{-}\left(t\right)=\max \left\{0, -I\left(t\right)\right\}$;
• $N\left(n\right)$ The number of replenishment orders triggered in $\left(0, n\right]$;
• $H$ The holding cost per item per unit time;
• $B$ The backorder cost per item per unit time;
• $O$ The ordering cost for each replenishment order;
• $TC\left(n\right)$ Total inventory cost for a given number of days ($n$);
• $TC{\left(n\right)}_{Ri}$ Total inventory cost for a given number of days ($n$) at retailer $i$;
• $TC{\left(n\right)}_{DC}$ Total inventory cost for a given number of days ($n$) at the distributor ($DC$).

Then, the simulated total cost for $n$ days, based on one simulation replication of the system, can be modeled as follows:

$$TC\left(n\right)=TC{\left(n\right)}_{R1}+TC{\left(n\right)}_{R2}+TC{\left(n\right)}_{DC}$$

(2)

where

$$TC{\left(n\right)}_{R1}=H_{R1}· {{\displaystyle \int }}_0^n{I^{+}}_{R1}\left(t\right)dt+B_{R1}· {{\displaystyle \int }}_0^n{I^{-}}_{R1}\left(t\right)dt+O_{R1}·N{\left(n\right)}_{R1}$$

(3)

$$TC{\left(n\right)}_{R2}=H_{R2}· {{\displaystyle \int }}_0^n{I^{+}}_{R2}\left(t\right)dt+B_{R2}· {{\displaystyle \int }}_0^n{I^{-}}_{R2}\left(t\right)dt+O_{R2}.N{\left(n\right)}_{R2}$$

(4)

$$TC{\left(n\right)}_{DC}=H_{DC}· {{\displaystyle \int }}_0^n{I^{+}}_{DC}\left(t\right)dt+B_{DC}· {{\displaystyle \int }}_0^n{I^{-}}_{DC}\left(t\right)dt+O_{DC}·N{\left(n\right)}_{DC}$$

(5)

When conducting a simulation study, a certain number of replications are needed to estimate a performance measure. These replications are required to find the system’s steady-state and should be determined before the actual run of the simulation model. Several pilot runs, as well as some calculation steps, are conducted using the sequential procedure proposed in [20] to find the number of replications required. Based on the calculated number of replications, $K$, the following model is developed to calculate the expected annual total cost, which is formulated with the following variables:

• β The source of system randomness, which represents the sequence of independent
• $U\left(0 , 1\right)$ random numbers that activate a certain number of simulation replications;
• ${\beta }^1, {\beta }^2, \dots , {\beta }^K$ Independent realizations from the distribution of $\beta$;
• $TC\left(n, {\beta }^k\right)$ The observed total inventory cost ($\mathrm{TC}$) in the interval $\left(0,n\right]$ for the $k^{th}$ simulation replication.

Then, the expected annual total inventory cost (ATC) can be modeled as follows:

$$E\left[ATC\right]=\frac{E\left[TC\left(n\right)\right]}{n}*n_{Year}$$

(6)

where, $n_{Year}$ is the number of days per year and

$$E\left[TC\left(n\right)\right]=\underset{K\to \infty }{\lim }\frac{1}{K}{{\displaystyle \sum }}_{k=1}^{K}TC\left(n, {\beta }^k\right)$$

(7)

In order to optimize the ordering policy ($Q,r$), a mixed integer-programming model is formulated with the following variables:

• $r_{LRi}$ The lower bound of the reorder point for retailer $i$, $i=1, 2$;
• $r_{LDC}$ The lower bound of the reorder point of the distributor;
• $r_{URi}$ The upper bound of the reorder point for retailer $i$;
• $r_{UDC}$ The upper bound of the reorder point for the distributor;
• $Q_{LRi}$ The lower bound of the order quantity for retailer $i$;
• $Q_{LDC}$ The lower bound of the order quantity for the distributor;
• $Q_{URi}$ The upper bound of the order quantity for retailer $i$;
• $Q_{UDC}$ The upper bound of the order quantity for the distributor.

The objective function is given by

$$\mathrm{Min}z=E\left[ATC\right]$$

(8)

subject to

$$r_{LRi}\le r_{Ri}\le r_{URi}\hspace{1em}i=1,2$$

(9)

$$r_{LDC}\le r_{DC}\le r_{UDC}$$

(10)

$$Q_{LRi}\le Q_{Ri}\le Q_{URi}\hspace{1em}i=1,2$$

(11)

$$Q_{LDC}\le Q_{DC}\le Q_{UDC}$$

(12)

Constraints (9)–(12) ensure that the reorder points and order quantities are within the corresponding lower and upper bounds.

2.2. Simulation-Optimization Model

As previously mentioned, in real-world problems, solving inventory systems analytically is challenging or extremely complex. This complexity is mainly due to the variability (or randomness) in demand and lead time. Based on the results reported in [11], most inventory models with a stochastic demand and lead time are a heuristic approximation of real-world systems. In stochastic continuous review systems ($Q,r$), analytical solutions can be obtained under some assumptions, such as normally distributed daily demand, independent lead time, and the out-of-stock time in a cycle being smaller than the cycle length [2]. However, using simulation-based optimization models, these assumptions can be easily met and included. Therefore, simulation is utilized to model and analyze complex inventory systems to obtain more practical results, which is the aim of this study.

2.2.1. Arena Model

As stated above, the system consists of a supplier and a distributor in the higher echelon, and two retailers in the lower echelon (see Figure 1). An Arena network model for the system is shown in Figure 2. The modeled network has five parts: each one represents a complete function of the system. In part 1, several variables related to retailers (1 and 2) are initialized, calculated, and updated. These variables are demand (customers’ orders) creations, inventory levels, inventory positions, and backorders (if any). It should be noted that the simulation model does not accept negative demand values that could be generated in the case of normal distribution. Instead, the negative demand values were regenerated again in a similar way to that of Al-harkan and Hariga [2]. In part 2, decisions related to placing orders (for both retailers) are made. In addition, inventory positions for the distributor and the associated retailer are updated. Then, in parts 3 and 4, replenishments are created for retailers 1 and 2, ordering costs are calculated for retailers 1 and 2, and backorders, if any, are updated. Finally, distributor replenishments are created, distributor ordering costs are calculated (if an order is made), and backorders (if assigned) are updated in part 5. Parameters were defined as variables in the “Variable module” as shown in

Table 1. Parameter definitions in the “Variable module”.

No.	Variable	Definition	No.	Variable	Definition
1	Setup_Cost	Fixed ordering cost	15	STOCK_D	Stock at DC
2	Short_Cost	Stock out cost	16	BACKO_R1	Backorder R1
3	Unit_Holding_Cost_day	Inventory holding cost	17	BACKO_R2	Backorder R2
4	Prod_Time	Unit Preparation time	18	BACKO_D	Backorder DC
5	Theta	Setup time	19	Order_Size_R1	Order size (R1)
6	Zeta	Shop floor queuing factor	20	Order_Size_R2	Order size (R2)
7	n_runs	No. of days (n)	21	Order_Size_D	Order size (DC)
8	n_year	No. of days per year $\left(n_{Year}\right)$	22	Reorder_Point_R1	Reorder point (R1)
9	No_Years	No. of years $\left(n/n_{Year}\right)$	23	Reorder_Point_R2	Reorder point (R2)
10	POS_R1	Inventory position (R1)	24	Reorder_Point_D	Reorder point (DC)
11	POS_R2	Inventory position (R2)	25	Total_Ordering_Cost_R1	Order cost of R1
12	POS_D	Inventory position (DC)	26	Total_Ordering_Cost_R2	Order cost of R2
13	STOCK_R1	Stock at R1	27	Total_Ordering_Cost_D	Order cost of DC
14	STOCK_R2	Stock at R2	28	Total_Ordering_Cost	Total ordering cost

. In addition,

Table 2. Cost calculations using the “Statistic module”.

Name	Type	Expression
Holding Cost	Time-Persistent	n_runs*Unit_Holding_Cost_day*MX(STOCK_R1,0)+n_runs* Unit_Holding_Cost_day*MX(STOCK_R2,0)+n_runs*Unit_Holding_Cost_day*MX(STOCK_D,0)
Annual Holding Cost	Output	DAVG(Holding Cost)/No_Years
Backlog Cost	Time-Persistent	n_runs*Short_Cost * BACKO_R1+n_runs*Short_Cost * BACKO_R2 + n_runs*Short_Cost * BACKO_D
Annual Backlog Cost	Output	DAVG(Backlog Cost)/No_Years
Annual Order Cost	Output	(Total_Ordering_Cost_R1)/No_Years+ (Total_Ordering_Cost_R2)/No_Years + (Total_Ordering_Cost_D)/No_Years
Annual Total Cost	Output	Ovalue(Annual Order Cost) +Ovalue(Annual Holding Cost)+ Ovalue(Annual Backlog Cost)

presents the calculations of the costs reported in Equations (2)–(7) which were calculated using the “Statistic module”.

The model is checked to ensure that it is error-free and running as per requirement. Several pilot runs are also conducted using the required simulation statistics to determine the required number of replications $K$, which is needed to calculate the expected annual total inventory cost in Equations (6) and (7). Results from the pilot runs were also used in the model validation (refer to Section 2.3.2). Based on the selected performance measures of the system, we verify the accuracy of the simulation model in Section 2.3.2.

2.2.2. OptQuest

We use $OptQuest$, which is a useful tool embedded in the Arena software, to find the best (near-optimal) order sizes and reorder points for the system. $OptQuest$ can be used to optimize complex systems that are difficult to formulate or solve analytically using classical optimization tools. It combines several metaheuristic approaches into a single search heuristic. These approaches are based on the tabu search algorithm, neural networks, and scatter search algorithm [14]. This advanced search framework is the key strength of $OptQuest$ [21,22]. $OptQuest$ performs a special (non-monotonic) search, where different evaluations are produced by the progressively generated inputs; these may not always be an improvement, but they provide an efficient path to the best solutions over time. In other words, if a candidate solution is not feasible, $OptQuest$ eliminates it and investigates alternatives that are potentially better [20].

The optimization approach used by $OptQuest$ requires a set of inputs to start the search process for better solutions in the model under investigation. These inputs are related to the model’s decision variables, objective function, and constraints. In our case, decision variables are the order size $Q$ and reorder points $r$ for the two retailers and the distributor (refer to Equations (9)–(12)). $OptQuest$ needs to specify the initial, lower, and upper values of variables to be optimized. The suggested initial values state the search starting point of each optimized variable and have a significant effect on the efficiency and effectiveness of the search process. The related constraints are the lower and upper bound values for each decision variable to be optimized. These lower and upper bound values are based on the analytical solutions obtained in Section 2.3.2 and preliminary experiments. The objective function to be minimized is the annual total cost ($ATC$) (see Equation (8)). For practical cases, $OptQuest$ allows simulation inputs to be defined as integers or continuous inputs, based on their nature. Moreover, a stopping criterion needs to be specified; in our case, the number of simulated candidate solutions is used as the stopping criterion.

2.3. Empirical Evaluation

The performance of the simulation-based optimization model is assessed in this section. The model configurations are first illustrated using an analytically solvable example (the case of normally distributed daily demands) that is used to validate the developed model. An experimental study is then conducted to solve the considered system. The effects of the variability of daily demand and lead time, as well as of the types of fitted distributions for the daily demand, on the inventory system’s performance are investigated.

2.3.1. Configurations

In order to clarify the developed simulation-based optimization model described in the previous section, we fitted the solved distributions into the system. In the existing literature, a similar case (two-stage supply chain with a supplier in the higher stage and a retailer in the lower stage) was solved analytically by [1,3]. Moreover, the same case was solved in [2] using a simulation optimization model for three cases of daily demand distributions: normal, gamma, and log-normal. The authors reported that the approximate problem is analytically solvable only for the normal distribution case, and it is hard to formulate and compute for other distributions. Therefore, we used the approximate analytical model to solve the system under the assumption that daily demands are normally distributed. We also assumed that the two retailers are identical and have the same average demand as the retailer in [1,3]. As a result, the average demand at the distributor was the sum of the demand at the two retailers. The rest of the required data were considered in a way similar to that used in [1,3]. A summary of the data is presented in

Table 3. Data of Retailer 1, Retailer 2, and Distributor.

Retailer 1 and Retailer 2
Average demand per year	$D_{R1,R2}=5000$
Average demand per day (250 days in a year)	${\mu }_{R1,R2}=20$
Inventory holding cost	$H_{R1,R2}=\$0.5$ per unit per day
Stock out cost	$B_{R1,R2}=600$ per unit (out of stock)
Setup time	${\theta }_{R1,R2}=0.125$ days
Ordering cost	$O_{R1,R2} =\$125$ per order
Unit preparation time	$p_{R1,R2}=0.025$ day per unit
Queueing factor	${\delta }_{R1,R2}=10$
Distributor DC
Average demand per year	$D_{DC}=10000$
Average demand per day (250 days in a year)	${\mu }_{DC}=40$
Inventory holding cost	$H_{DC}=\$0.5$ per unit per day
Stock out cost	$B_{DC}=600$ per unit (out of stock)
Setup time	${\theta }_{DC}=0.125$ days
Ordering cost	$O_{DC} =\$125$ per order
Unit production time	$p_{DC}=0.025$ day per unit
Queueing factor	${\delta }_{DC}=10$

. In order to verify the developed simulation-based optimization model, the results of the developed model were compared with the results of the analytical model.

The simulation-based optimization model has several configurations that need to be set before solving the special case under consideration. The lower bounds of the decision variables $Q$ and $r$ for the two retailers ($R1$ and $R2$) and the distributor ($DC$) are (30, 30, 30) and (300, 300, 700), respectively. The upper bounds for the retailers and distributor are (150, 150, 300) and (800, 800, 1700), respectively. These lower and upper bound values of the Q and r are based on the analytical solutions obtained in Section 2.3.2 and preliminary experiments. The developed model is simulated and optimized for a period of 10 years (2500 days). The daily demand is normally distributed with a mean $\mu =20$ and standard deviation (SD)$\sigma =5, 8, 10, 15, \mathrm{and} 22$. We used different values for the SD to study the effect of the variability in demand and lead time on the simulated annual total inventory cost for the specified period (2500 days). The number of replications $K$ is estimated (based on the simulated annual total cost as the performance measure) using the sequential procedure proposed by Law and Kelton [20]. The relative error considered is $\gamma =0.15$ with a confidence level of 95%. We found that 10 replications are enough to obtain the relative error. However, the actual simulation runs are conducted with 20 replications to obtain a better estimate of the annual total cost.

2.3.2. Model Validation

The results of the proposed simulation-based optimization model and the analytical approximation model are reported in

Table 4. Results of the simulation-optimization model and the analytical method.

	$\mathit{\sigma }$	5		8		10		15		22
	Model	Analytical	Simulation	Analytical	Simulation	Analytical	Simulation	Analytical	Simulation	Analytical	Simulation
Retailer 1	Q	77.75	83.00	68.53	67.00	63.56	68.00	53.96	55.00	44.89	45.00
	R	475.67	494.00	462.48	460.00	458.24	480.00	458.43	483.00	474.70	500.00
	AOC	8038.51	7536.25	9119.43	9363.75	9833.48	9446.88	11,582.11	12,829.38	13,923.26	18,874.38
	AHC	12,599.20	12,143.40	16,134.19	16,699.94	18,402.78	17,527.30	23,824.81	21,619.28	30,962.98	22,640.41
	ABC	863.48	766.50	1288.83	224.44	1545.00	1468.31	2120.88	772.10	2818.76	8632.84
	ATC_R1	21,501.19	20,446.15	26,542.45	26,288.14	29,781.26	28,442.48	37,527.80	35,220.75	47,705.00	50,147.62
Retailer 2	Q	77.75	79.00	68.53	71.00	63.56	63.00	53.96	52.00	44.89	45.00
	R	475.67	474.00	462.48	465.00	458.24	446.00	458.43	467.00	474.70	499.00
	AOC	8038.51	7877.50	9119.43	8847.50	9833.48	10,143.75	11,582.11	13,684.38	13,923.26	18,720.63
	AHC	12,599.20	12,140.87	16,134.19	15,022.77	18,402.78	16,409.86	23,824.81	21,192.37	30,962.98	22,867.87
	ABC	863.48	509.22	1288.83	1308.88	1545.00	1738.69	2120.88	1283.83	2818.76	8279.70
	ATC_R2	21,501.19	20,527.59	26,542.45	25,179.15	29,781.26	28,292.30	37,527.80	36,160.58	47,705.00	49,868.19
Distributor	Q	113.40	111.00	103.14	107.00	97.86	107.00	88.17	103.00	79.49	79.00
	R	1292.68	1297.00	1249.32	1310.00	1234.60	1368.00	1229.10	1490.00	1265.00	1476.00
	AOC	11,022.97	11,247.50	12,120.03	11,741.25	12,773.02	11,989.38	14,177.22	13,788.75	15,723.54	21,449.38
	AHC	20,673.45	25,025.53	27,442.41	30,489.11	31,863.16	34,770.58	42,685.87	40,124.62	57,469.79	47,645.67
	ABC	2564.72	762.00	4311.95	681.00	5540.19	1038.00	8778.42	1698.00	13,598.61	1521.00
	ATC_D	34,261.14	37,035.03	43,874.39	42,911.36	50,176.37	47,797.96	65,641.51	55,611.37	86,791.94	70,616.04
ATC		77,263.52	78,008.77	96,959.29	94,378.65	109,738.90	104,532.74	140,697.10	126,992.70	182,201.90	170,631.85

. The terms used in Table 4 are order size (Q), reorder point (R), annual ordering cost (AOC), annual holding cost (AHC), annual backorders cost (ABC), the annual total cost for retailer 1 (ATC_R1), the annual total cost for retailer 2 (ATC_R2), the annual total cost for the distributor (ATC_D), and the total annual cost for the system (ATC).

A comparison between the two models was conducted using the results reported in Table 4. The comparison demonstrated that the proposed model provided accurate results. In fact, the simulation-based optimization model showed (for different values of SD) a slight deviation compared to the analytical model results. The maximum absolute ATC deviation is $9.74\%$ when $\sigma =15$. In almost all scenarios, which are different values of SD, the optimization-based simulation model provided smaller ATC values than the approximate analytical model. Results showed that the optimization-based simulation model outperformed the approximate analytical model, except for $\sigma =5$, when the absolute ATC deviation was $0.96\%$, which is a very small value.

Figure 3 also represents the results of the two models (analytical and simulation). The terms used in Figure 3 are the same as those in Table 4. ATC increases with the increase in the demand variability ($\sigma$), which demonstrates the significance of the effect of demand variability on the ATC. Moreover, a paired t-confidence interval was used to analyze the differences in the means of the ATC generated by the analytical and simulation models. At a confidence level of $\alpha =95\%$, the reported paired t-confidence interval was ($-13994, 1067$) and the p-value was $0.079$. Based on the reported interval and $p$-value, we can perceive that zero is included in the interval and the p-value is larger than 0.05; this means that there is no significant difference between the results generated by the two models. Thus, we can report that our model is accurate and valid.

3. Results and Discussion

As discussed in the previous section, a special case of the problem under consideration was solved analytically to validate the developed model. The solved case was a two-stage supply chain (one supplier and one retailer) with a lot size-dependent lead time. Moreover, it was assumed that the daily demand is normally distributed. However, it has been reported that the analytical solution to this special case is extremely difficult to arrive at [2,6]. Therefore, a simulation-based optimization model was developed to solve more complex real-world problems. In this study, we considered the stochastic inventory continuous review policy ($Q,r$) in a three-stage two-echelon supply chain system with lot size-dependent lead time,$L\left(Q\right)$. The system has a supplier, a distributor, and two identical retailers (see Figure 1).

In this section, an experimental study was conducted to solve the system. The developed model was run for different daily demand distributions to study the effects of the variability of daily demand and lead time on the simulated ATC of the inventory system. We used three distributions for daily demand: normal, gamma, and log-normal. The simulation runs are conducted using different values of the coefficient of variance $cv$ ($cv=0.25, 0.4, 0.5, 0.75, \mathrm{and} 1.1$) for each of the three distributions of daily demand. Based on the relative conversion formulas reported in [23] and the reported parameter values in [2], the used parameters of the gamma and log-normal distributions were chosen so that they provide the same values for the normal distribution parameters ($\mu$ and $\sigma$).

Table 5. The parameters of the three distributions.

$\mathit{c}\mathit{v}$	Normal $\left(\mathit{\mu },\mathit{\sigma }\right)$	Gamma $\left(\mathit{\alpha },\mathit{\beta }\right)$	Log-Normal $\left(\mathit{\alpha },\mathit{\beta }\right)$
0.25	$\left(20, 5\right)$	(16.00, 1.25)	(2.965, 0.246)
0.40	$\left(20, 8\right)$	(6.25, 3.20)	(2.9215, 0.385)
0.50	$\left(20, 10\right)$	(4.00, 5.00)	(2.884, 0.472)
0.75	$\left(20, 15\right)$	(1.78, 11.25)	(2.773, 0.668)
1.10	$\left(20, 22\right)$	(0.83, 24.20)	(2.599, 0.891)

provides a summary of the parameters used for the three distributions. The results of the experimental study are summarized in

Table 6. Results for different types of normal, gamma, and log-normal demand.

Inventory	cv = σ/μ	0.25			0.4			0.5			0.75			1.1
	pdf	Normal	Gamma	Log-Normal	Normal	Gamma	Log-Normal	Normal	Gamma	Log-Normal	Normal	Gamma	Log-Normal	Normal	Gamma	Log-Normal
Retailer 1	Q	83	46	83	67	49	67	68	46	71	55	55	62	45	45	42
	R	494	302	491	460	363	463	480	366	517	483	493	582	500	498	750
	AOC	7536.3	13,516.3	7516.3	9363.8	12,696.9	9313.8	9446.9	13,493.8	8772.5	12,829.4	11,300.0	10,091.3	18,874.4	13,882.5	14,979.4
	AHC	12,143.4	9104.9	11,907.7	16,699.9	14,977.9	17,316.8	17,527.3	17,092.4	21,910.0	21,619.3	27,828.8	34,772.8	22,640.4	33,342.8	65,892.3
	ABC	766.5	822.5	1590.7	224.4	1086.2	1194.6	1468.3	757.6	1750.9	772.1	987.0	6041.0	8632.8	16,519.9	30,001.2
	ATC_R1	20,446.2	23,443.6	21,014.6	26,288.1	28,761.0	27,825.2	28,442.5	31,343.7	32,433.4	35,220.8	40,115.8	50,905.1	50,147.6	63,745.2	110,872.9
Retailer 2	Q	79	55	71	71	51	65	63	51	91	52	50	53	45	45	47
	R	474	352	434	465	358	455	446	397	629	467	474	533	499	500	750
	AOC	7877.5	11,306.3	8782.5	8847.5	12,192.5	9580.0	10,143.8	12,228.1	6886.9	13,684.4	12,496.3	11,750.6	18,720.6	13,912.5	13,273.1
	AHC	12,140.9	10,298.3	11,507.8	15,022.8	13,281.2	17,504.2	16,409.9	18,005.9	24,430.9	21,192.4	28,034.7	33,798.7	22,867.9	33,536.7	63,723.3
	ABC	509.2	910.4	1140.4	1308.9	1521.0	1424.7	1738.7	1096.2	655.0	1283.8	1045.0	2317.1	8279.7	15,886.4	726.7
	ATC_R2	20,527.6	22,515.0	21,430.7	25,179.2	26,994.8	28,508.9	28,292.3	31,330.2	31,972.8	36,160.6	41,575.9	47,866.4	49,868.2	63,335.5	77,723.2
Distributor	Q	111	103	111	107	83	107	107	91	110	103	97	95	79	93	75
	R	1297	1196	1304	1310	1042	1313	1368	1155	1397	1490	1345	1343	1476	1500	1311
	AOC	11,247.5	12,067.5	11,246.3	11,741.3	14,953.8	11,660.0	11,989.4	13,660.0	11,378.1	13,788.8	12,860.0	13,153.1	21,449.4	13,482.5	16,713.1
	AHC	25,025.5	22,417.9	25,986.9	30,489.1	26,783.3	31,858.1	34,770.6	31,345.9	38,374.6	40,124.6	47,547.0	49,432.8	47,645.7	71,014.8	68,653.6
	ABC	762.0	435.0	909.0	681.0	1215.0	1092.0	1038.0	3270.0	4582.5	1698.0	4360.5	7125.0	1521.0	3660.0	4290.0
	ATC_D	37,035.0	34,920.4	38,142.1	42,911.4	42,952.1	44,610.1	47,798.0	48,275.9	54,335.2	55,611.4	64,767.5	69,711.0	70,616.0	88,157.3	89,656.7
ATC		78,008.8	80,879.0	80,587.4	94,378.7	98,707.8	100,944.2	104,532.7	110,949.9	118,741.4	126,992.7	146,459.2	168,482.4	170,631.9	215,238.1	278,252.8

and Figure 4.

As shown in Figure 4, the ATC increases with the increase in the coefficient of variance ($cv$) for all the distributions and increases the most when the daily demand follows log-normal. This indicates a significant effect of the $cv$ and the distribution type on the ATC. The results in Table 6 demonstrate that the solutions generated (order sizes and reorder points) for different distributions of daily demand are not the same even if the associated total costs are close to each other. This explains the difficulty in solving inventory systems with a stochastic daily demand.

Moreover, results have been tested statistically using the analysis of variance (ANOVA) with 0.95 confidence to examine whether the $cv$ or the type of the daily demand distribution has a significant effect on the simulated ATC. p-values less than 0.05 indicate model terms are significant.

From

Table 7. ANOVA table for the ATC.

Source	DF	Seq SS	Contribution	Adj SS	Adj MS	F-Value	p-Value
CV	4	3.76 × 1010	84.60%	3.76 × 1010	9.39 × 109	19.50	0.000
PDF	2	2.98 × 109	6.72%	2.98 × 109	1.49 × 109	3.10	0.101
Error	8	3.85 × 109	8.68%	3.85 × 109	4.82 × 108
Total	14	4.44 × 1010	100.00%

, ANOVA test results revealed a significant effect of $cv$ on ATC since its p-value is 0.00, less than 0.05. Overall, the results showed that $cv$ has the most significant effect on ATC. Figure 5a,b illustrate the interaction between the distribution types and $cv$. From these interaction plots, it is evident that the ATC results difference between the three distributions is high with higher values of $cv$ (Figure 5a). In addition, it is clear that the ATC increases with an increase in $cv$ (Figure 5b).

4. Conclusions

In most inventory models, inventory costs are optimized with the assumption that the lead time is constant. Unfortunately, this is not the case in many real-world inventory systems. In practice, inventory systems are complex and difficult to formulate because of the stochastic nature of demand and lead times. In some cases, the lead time could have a few controllable components, such as lot size, production time, setup time, waiting and moving times, and lot size inspection time, which are related to order preparation and production. This study aimed to solve a complex inventory system that consists of two identical retailers, a distributor, and a supplier. The system is a stochastic inventory continuous review system ($Q,r$) with the lead time depending on lot size, unit production time, setup time, and a shop floor factor that represents moving, waiting, and lot size inspection times. A simulation-based optimization model was developed using Arena and $OptQuest$. A special case of the system was solved analytically to validate the model. Results indicate that the developed model was more accurate. Moreover, the model was tested for a more complex inventory system with different distributions of the system’s daily demand.

The results demonstrated that the optimized solutions (order sizes and reorder points) for different distributions of daily demand are not the same even if the associated total costs are close to each other. This explains the difficulty in solving inventory systems with a stochastic daily demand. Therefore, simulation-based optimization models can be considered as a powerful tool to optimize complex inventory systems without the need to make restrictive assumptions. Extensions of this work can be made by using different inventory control policies, taking into account supplier capacity, multi-products, and finally, cooperation between retailers.