# Robust Inventory System Optimization Based on Simulation and Multiple Criteria Decision Making

Authors: Ahmad Mortazavi; Alireza Arshadi khamseh
Journal: Modelling and Simulation in Engineering
 (2014)

---

## Abstract

Inventory management in retailers is difficult and complex decision making process which is related to the conflict criteria, also existence of cyclic changes and trend in demand is inevitable in many industries. In this paper, simulation modeling is considered as efficient tool for modeling of retailer multiproduct inventory system. For simulation model optimization, a novel multicriteria and robust surrogate model is designed based on multiple attribute decision making (MADM) method, design of experiments (DOE), and principal component analysis (PCA). This approach as a main contribution of this paper, provides a framework for robust multiple criteria decision making under uncertainty.


---

## Body

## 1. Introduction

Supply chain is complicated system that includes many components such as suppliers, factories, distribution centers, and retailers. These components are connected to each other by three streams of financial, information, and material flow. In many references, there is a fundamental hypothesis that financial flow is upstream from customers to suppliers and material flow is downstream from suppliers to customers while information flow is mutual. By the way, sometimes, this hypothesis is neglected about material flow (e.g., reveres logistic).Material flow plays dominant role in supply chain and is defined as inventory problem. Inventory actually is the bridge that connects material handling and production planning to the supply chain [1]. On the other hand, retailers in highly competitive market face dynamic change of demand including seasonal cyclic change and long-term trend. These situations caused multiproduct inventory systems to change into complex multicriteria systems.Complex systems are challenging in the case of modelling and computation. In most of the problems, modelling of complex systems is very time consuming and also is not error free. Furthermore, modelling of complex systems needs too much computational effort to solve and sometimes they are not solvable in reasonable time. In this situation, importance of data-driven methods emerges [2–4].Although model-based approach (i.e., simulation model) has many advantages, such as what if analysis and ability of scenario generation, data-driven approach is less time consuming and also easy to implement. So, these approaches can be combined to achieve tailored approach which inherits advantages of model-based and data-driven approaches. In this research we used discrete event simulation model and data-driven methods like principal component analysis and multiple attribute decision making to design an efficient framework for robust and multiobjective optimization of retailer inventory system with multiproduct as complex system.The rest of the paper is structured as follows. Literature review is presented in Section2. Problem statement is presented in Section 3, proposed framework which contains simulation modeling and surrogate design is described in Section 4, Section 5 provides numerical result, and finally Section 6 is dedicated to conclusion.


## 2. Literature Review

As the investigated problem encompasses two issues of simulation optimization and inventory problem, related literature is categorized in two separated parts for providing more supportive literature review.


### 2.1. Inventory Problem

Basic problems of inventory system are studied thoroughly by [5–7]. Economical order quantity (EOQ) is the simplest model of inventory problem; in this model, demand of each period is constant and time independent. EOQ does not consider lost sale, back order, and other cost for reason of simplification and emphasizes only on holding and ordering cost. Cheng [8] developed a model of inventory system with cost-dependent demand and included production cost in proposed model. Chen et al. [9] studied back order by fuzzy technique and Zhao et al. [10] suggested analytic model with demand according to time series. They concluded that, in time-dependent demand, efficiency of EOQ model increases with shorter lead time and weaker autocorrelation.Basic models that were offered for inventory problem only have one objective function, including different inventory costs, while advancement in technology and intensification of business competition caused necessity of other criteria. To the best of our knowledge, in recent years, service level has been observed in variety of supply chain and inventory problem as performance criterion. Adding service level criterion, classic definition of inventory problem is changed into multiobjective optimization. In such problem, inventory costs should be minimized while service level should be maximized. Available models for this problem are divided in two groups of deterministic and stochastic. Both of these models can be solved by three approaches of analytic methods, for example, mathematical programming, metaheuristic methods, and simulation. Lau et al. [11] benefited from simulation to compare four inventory management policies with two criteria of cost and service levels. They also surveyed preorder and information sharing impacts on their models. Xu and Zhao [12] used fuzzy rough simulation to optimize multiobjective problem to minimize wasted cost and maximize expected value of revenue. Hnaien et al. [13] Surveyed two-level assembly system with two objectives of service level and maintenance cost. They considered stochastic lead time and applied genetic algorithm to solve this problem.Although cost and service levels are important criteria for inventory system performance, there are other criteria that should be considered such as amount of systems inventory, which is important factor with significant impact on inventory system behavior. Because of lead time uncertainty, that origins from natural disasters and transportation problems, organizations face delay in delivery, so most of them hold safety stocks. This phenomenon is the main cause that leads to the increase of inventory in hand. As inventory in hand increases, the inventory system is faced with holding cost and other problems like decrease in quality, lack of flexibility, and so on. So, amount of systems inventory can be considered as performance factor and inventory in hand should be minimized as mentioned in the just in time (JIT) philosophy. Purnomo et al. [14] researched about influence of periodic replenishment and continue replenishment inventory policies on supply chain and considered both inventory in hand and work in progress as performance factors.Now in recent years modeling of inventory system as multiobjective and stochastic problem is an interesting area for research.


### 2.2. Simulation Optimization

One of the well-known simulation methods is discrete event simulation that is based on stochastic processes and could be efficient tool to capture stochastic behaviour of different systems. While discrete event simulation has several advantages, it is not optimization tool individually [15]. By the way, because of its flexibility, simulation can be coupled with other techniques such as metaheuristic algorithms or stochastic methods. This synthesizes makes the powerful and advantageous approach of simulation optimization with vast area of research. Simulation optimization is powerful arsenal for optimization of complex systems such as military, aerospace, and supply chain [16]. To the best of our knowledge, there are three main optimization techniques that were reported as suitable techniques for simulation optimization. These considered methods are metaheuristic optimization, stochastic approximation (SA) methods, and surrogate models.Fu [17] extensively described role of applied methods in simulation optimization and also surveyed techniques employed in optimization package of simulation software. Wang [18] used hybrid approach including genetic algorithm and artificial neural network for simulation optimization and Keskin et al. [19] applied discrete event simulation and scatter search algorithm for optimization of inventory system and vendor selection. Mazhari et al. [20] developed a simulation optimization framework based on hybrid simulation model (system dynamic and agent based model) and metaheuristic algorithm. Also, Duan and Liao [21] applied metaheuristic approach for developing simulation optimization framework in order to optimize replenishment policy of inventory system in capacitated supply chain.Although using metaheuristic optimization algorithms is straight forward approach for simulation optimization, it is time consuming and needs high level of computational effort. So, it is inefficient in case of simulation optimization with more than one objective function.While metaheuristic algorithms use stochastic searching methods, SA is based on gradient search. Because of noisy situation of observations, SA algorithms consider expected value of objective function. SA family includes attractive methods because their convergence is guaranteed theoretically. Simultaneous perturbation stochastic approximation (SPSA) is noteworthy algorithm of SA family. The theoretical aspects of SPSA are deeply described by Spall’s [22–25] proposed simulation optimization framework for inventory control in supply chain based on SPSA.In contrary to two former methods, surrogate modelling is postprocessing method, so it is less time consuming. In surrogate modelling, the main idea is to fit single surface to the decision space and use this surface instead of simulation model for optimization. In this area, response surface methodology (RSM) [26] and supervised learning methods (e.g., artificial neural network or support vector machine) are considerable. For instance, Can and Heavey [27] applied artificial neural network to develop surrogate model for discrete event simulation. Azadeh et al. [28] used artificial neural network for the designing of simulation optimization framework and they applied proposed framework for optimization of waiting time in tandem queue systems. Bornatico et al. [29] proposed a surrogate model based on redial basis function for simulation optimization of energy systems. Wan et al. [30] designed simulation optimization framework using least square support vector machine (LSSVM) for optimization of inventory level in three-stage supply chain. They also showed that proposed framework leads to better solution with less number of simulation runs in comparison with SPSA algorithm. Surrogate modelling is less time consuming in comparison to metaheuristic or SPSA approaches, but this approach loses accuracy in multiple objectives problem solving.Although a bunch of papers published in simulation optimization area, to the best of our knowledge a tiny number of them are dedicated to the multiobjective optimization [31] and robustness [32, 33]. In this case, using metaheuristic and SPSA approaches is very time consuming and it is not economic for optimization of simulation model with accurate details. On the other hand, all of the reviewed approaches lose their accuracy when there are multiple objectives. With these considerations, this paper purposes a framework for optimization of detailed simulation model of inventory system with multiple objectives. Proposed framework is less time consuming in comparison with metaheuristic or SPSA approaches while it provides robust and accurate solutions. So, the proposed framework is relatively new and contribution of this research entails threefold as follows.(i)
We modelled cyclic and long-term demand based on nonparametric time series modelling for more realistic consideration.
(ii)
We proposed surrogate model for robust and multiobjective optimization of multiproduct inventory system based on discrete event simulation, full factorial design of experiments (DOE), and multiple attribute decision making (MADM) technique.
(iii)
Due to the stochastic nature of objective function, we employed principal component analysis (PCA) as statistical method to improve MADM performance.


## 2.1. Inventory Problem

Basic problems of inventory system are studied thoroughly by [5–7]. Economical order quantity (EOQ) is the simplest model of inventory problem; in this model, demand of each period is constant and time independent. EOQ does not consider lost sale, back order, and other cost for reason of simplification and emphasizes only on holding and ordering cost. Cheng [8] developed a model of inventory system with cost-dependent demand and included production cost in proposed model. Chen et al. [9] studied back order by fuzzy technique and Zhao et al. [10] suggested analytic model with demand according to time series. They concluded that, in time-dependent demand, efficiency of EOQ model increases with shorter lead time and weaker autocorrelation.Basic models that were offered for inventory problem only have one objective function, including different inventory costs, while advancement in technology and intensification of business competition caused necessity of other criteria. To the best of our knowledge, in recent years, service level has been observed in variety of supply chain and inventory problem as performance criterion. Adding service level criterion, classic definition of inventory problem is changed into multiobjective optimization. In such problem, inventory costs should be minimized while service level should be maximized. Available models for this problem are divided in two groups of deterministic and stochastic. Both of these models can be solved by three approaches of analytic methods, for example, mathematical programming, metaheuristic methods, and simulation. Lau et al. [11] benefited from simulation to compare four inventory management policies with two criteria of cost and service levels. They also surveyed preorder and information sharing impacts on their models. Xu and Zhao [12] used fuzzy rough simulation to optimize multiobjective problem to minimize wasted cost and maximize expected value of revenue. Hnaien et al. [13] Surveyed two-level assembly system with two objectives of service level and maintenance cost. They considered stochastic lead time and applied genetic algorithm to solve this problem.Although cost and service levels are important criteria for inventory system performance, there are other criteria that should be considered such as amount of systems inventory, which is important factor with significant impact on inventory system behavior. Because of lead time uncertainty, that origins from natural disasters and transportation problems, organizations face delay in delivery, so most of them hold safety stocks. This phenomenon is the main cause that leads to the increase of inventory in hand. As inventory in hand increases, the inventory system is faced with holding cost and other problems like decrease in quality, lack of flexibility, and so on. So, amount of systems inventory can be considered as performance factor and inventory in hand should be minimized as mentioned in the just in time (JIT) philosophy. Purnomo et al. [14] researched about influence of periodic replenishment and continue replenishment inventory policies on supply chain and considered both inventory in hand and work in progress as performance factors.Now in recent years modeling of inventory system as multiobjective and stochastic problem is an interesting area for research.


## 2.2. Simulation Optimization

One of the well-known simulation methods is discrete event simulation that is based on stochastic processes and could be efficient tool to capture stochastic behaviour of different systems. While discrete event simulation has several advantages, it is not optimization tool individually [15]. By the way, because of its flexibility, simulation can be coupled with other techniques such as metaheuristic algorithms or stochastic methods. This synthesizes makes the powerful and advantageous approach of simulation optimization with vast area of research. Simulation optimization is powerful arsenal for optimization of complex systems such as military, aerospace, and supply chain [16]. To the best of our knowledge, there are three main optimization techniques that were reported as suitable techniques for simulation optimization. These considered methods are metaheuristic optimization, stochastic approximation (SA) methods, and surrogate models.Fu [17] extensively described role of applied methods in simulation optimization and also surveyed techniques employed in optimization package of simulation software. Wang [18] used hybrid approach including genetic algorithm and artificial neural network for simulation optimization and Keskin et al. [19] applied discrete event simulation and scatter search algorithm for optimization of inventory system and vendor selection. Mazhari et al. [20] developed a simulation optimization framework based on hybrid simulation model (system dynamic and agent based model) and metaheuristic algorithm. Also, Duan and Liao [21] applied metaheuristic approach for developing simulation optimization framework in order to optimize replenishment policy of inventory system in capacitated supply chain.Although using metaheuristic optimization algorithms is straight forward approach for simulation optimization, it is time consuming and needs high level of computational effort. So, it is inefficient in case of simulation optimization with more than one objective function.While metaheuristic algorithms use stochastic searching methods, SA is based on gradient search. Because of noisy situation of observations, SA algorithms consider expected value of objective function. SA family includes attractive methods because their convergence is guaranteed theoretically. Simultaneous perturbation stochastic approximation (SPSA) is noteworthy algorithm of SA family. The theoretical aspects of SPSA are deeply described by Spall’s [22–25] proposed simulation optimization framework for inventory control in supply chain based on SPSA.In contrary to two former methods, surrogate modelling is postprocessing method, so it is less time consuming. In surrogate modelling, the main idea is to fit single surface to the decision space and use this surface instead of simulation model for optimization. In this area, response surface methodology (RSM) [26] and supervised learning methods (e.g., artificial neural network or support vector machine) are considerable. For instance, Can and Heavey [27] applied artificial neural network to develop surrogate model for discrete event simulation. Azadeh et al. [28] used artificial neural network for the designing of simulation optimization framework and they applied proposed framework for optimization of waiting time in tandem queue systems. Bornatico et al. [29] proposed a surrogate model based on redial basis function for simulation optimization of energy systems. Wan et al. [30] designed simulation optimization framework using least square support vector machine (LSSVM) for optimization of inventory level in three-stage supply chain. They also showed that proposed framework leads to better solution with less number of simulation runs in comparison with SPSA algorithm. Surrogate modelling is less time consuming in comparison to metaheuristic or SPSA approaches, but this approach loses accuracy in multiple objectives problem solving.Although a bunch of papers published in simulation optimization area, to the best of our knowledge a tiny number of them are dedicated to the multiobjective optimization [31] and robustness [32, 33]. In this case, using metaheuristic and SPSA approaches is very time consuming and it is not economic for optimization of simulation model with accurate details. On the other hand, all of the reviewed approaches lose their accuracy when there are multiple objectives. With these considerations, this paper purposes a framework for optimization of detailed simulation model of inventory system with multiple objectives. Proposed framework is less time consuming in comparison with metaheuristic or SPSA approaches while it provides robust and accurate solutions. So, the proposed framework is relatively new and contribution of this research entails threefold as follows.(i)
We modelled cyclic and long-term demand based on nonparametric time series modelling for more realistic consideration.
(ii)
We proposed surrogate model for robust and multiobjective optimization of multiproduct inventory system based on discrete event simulation, full factorial design of experiments (DOE), and multiple attribute decision making (MADM) technique.
(iii)
Due to the stochastic nature of objective function, we employed principal component analysis (PCA) as statistical method to improve MADM performance.


## 3. Problem Statement

The problem is concerning retailer who sales office furniture and facility. The retailer sales four products, respectively,A, B, C, and D; the aim is the optimization of inventory system according to information which is adapted from local business. Key features of retailer products from inventory view point are as shown in Table 1.Table 1
Product type
Average demand
Holding cost in planning period
Reorder cost
Lost sale cost
A
987
30000
50000
20000
B
1520
30000
50000
15000
C
1598
30000
50000
25000
D
1569
30000
50000
18000In this table, second column gives average demand of each product type in a year, third column provides holding cost of each product in a planning period, fourth column is dedicated to ordering cost of each type of products, and finally numbers of fifth column are cost of lost sales which are incurred to retailers when they cannot satisfy the demand of customers for each type of products.The fundamental assumption that should be considered in this problem is as follows.(1)
Order cost for each type of products includes transportation and order registration cost.
(2)
There is no backlog inventory so inventory level is nonnegative all the time.
(3)
Profit of each product is considered as lost sale cost because unavailable products incur lost profits that are interpreted as cost of lost sale.
(4)
According to the adapted information, these products have five years life cycle and then will be substituted with new products.
(5)
Planning periods for system under study are as long as 20 days.The notations that will be used to describe the problems are as follows.Indices. Consider the following:t: index of planning periods t
=
1,2
,
3
,
…
,
T;
i: index of demands in planning period i
=
1,2
,
3
,
…
,
N;
j: index of orders in planning period j
=
1,2
,
3
,
…
,
M.Parameters. Consider the following:C
r: reorder cost;
C
h: holding cost of each product in planning period;
C
l: cost of lost sale for each product.Variables. Consider the following:I
t
l: inventory level in tth planning period;
I
t
p: inventory position in tth planning period;
d
i
t: quantity of ith demand in tth planning period;
n
t: number of reorder in tth planning period;
x
i
t: 1 if d
i
t is less than I
t
l, 0 otherwise;
Q
j
t: quantity of jth order in tth planning period;
L: lead time for organized orders;
R
p: reorder point.Demand of each product in the planning period (d
i
t) is stochastic variable and is generated by nonhomogeneous Poisson distribution. So, total number of arrived demand in planning period (N) is probabilistic. C
r, C
h, and C
l are different costs of inventory system according to Table 1. I
t
l is inventory level and refers to physical quantity of inventory which is available in retailer while I
t
p is the position of inventory and includes quantity of on-order inventory in addition to inventory level in planning period. x
i
t is binary variable which is one if inventory level is greater than the arrived demand. So, if x
i
t is one, d
i
t can be satisfied and otherwise it is lost sale. R
p is the reorder point for organizing of new order. In other words, if inventory level reaches R
p or less, a new order would be organized with quantity of Q
j
t. Each organized order reaches the retailer and increases inventory level after passing of lead time (L). Considering assumptions and described notations, the following equations are the main objective functions of defined problem:
(1)
Min
⁡
∑
t
=
1
T
C
r
n
t
+
∑
t
=
1
T
C
h
I
t
l
+
∑
t
=
1
T
∑
i
=
1
N
C
B
(
1
-
x
i
t
)
d
i
t
(2)
Max
⁡
∑
t
=
1
T
∑
i
=
1
N
x
i
t
d
i
t
∑
t
=
1
t
∑
i
=
1
N
d
i
t
(3)
Min
⁡
∑
t
=
1
T
I
t
p
T
.
In (1), the objective is minimizing the total cost including costs that depend on reordering, handling, and lost sales. Lost sales not only incur excess cost but also decrease retailer credit. So, (2) is considered to maximize service level independently. In (2), the objective function increases, where x
i
t is 1 for ith demand in tth period and such situation is possible if I
t
l is greater than d
i
t. In fact, (2) causes increase in inventory level while (1)–(3) causes decrease in inventory level. State of inventory level depends on number of orders in each planning period (n
t), quantity of orders (
Q
j
t
), and quantity of demands (d
i
t) while state of inventory position depends on inventory level and lead time (L), so (3) is responsible for minimizing average of inventory position including inventory in hand (
I
t
l
) and on-order inventory. Average of inventory position should be minimized in order to improve flexibility of retailer and approach to the just in time (JIT).Demand of each product follows different pattern with both long-term and cyclic trend. So, mentioned objective functions are considered individually for each product type.


## 4. Proposed Framework

In the defined problem, demand of products exposes highly dynamic pattern and as time passes, demand and its variation increase; hence, multiresolution method is employed for demand modelling. Also, three different policies for inventory control are considered, which are reordered based on fixed quantity (FQ), fixed interval (FI), and demand forecasting (DF). Simulation of developed model is implemented in Arena 13.5. Optimization of simulation model is performed by surrogate model that is based on full factorial design of experiment (DOE). For construction of decision space, DOE factors include inventory policy, reorder point, and lead time with three levels for each of them. So, there are3
3
=
27 different combinations of decision variables to form feasible scenarios. Ranking of produced scenarios is accomplished by MADM technique. For ranking of scenarios, three objective function values are considered (i.e., cost, service level, and average of inventory position) and, in addition, robustness of service level against demand fluctuation is considered. Also PCA is applied for more realistic weighting of objective function values based on their statistical influence on improvement of other objectives. Finally, interacting plot is employed for sensitivity analysis and investigation of solutions in detail.


### 4.1. Simulation Modelling

Simulation modelling of problem consists of two parts which are modelling of demand and modelling of inventory policies. In this paper, demand is nonhomogeneous Poisson process and three different inventory policies based on continues reviewing, periodic reviewing, and periodic reviewing with forecasting of future demand are considered.


#### 4.1.1. Modelling of Demand

For customers demand modeling, multiresolution method is applied. Kuhl and Wilson [34] developed this method for simulation of nonhomogeneous Poisson process with trend and cyclic changes. This method estimates mean intensity function and the nonparametric nature of this method is one of the most important advantages in comparison with other methods. So, it is independent of statistical parameters and applicable in variety of problems. Furthermore, multiresolution method can support combination of multiple cyclic changes simultaneously. Another advantage of multiresolution is its capability in the modelling of nonsymmetric cyclic pattern. As in our case, demand has nonsymmetric pattern with cyclic changes and long-term trend, and multiresolution approach is a reasonable choice. More theoretical and application of used method are provided in [34].


#### 4.1.2. Modelling of Inventory Policies

As the main effort of this paper is inventory system optimization, modelling of inventory policies plays an important role in this problem. In this problem, optimization is manipulated by selection of appropriate inventory policy and configuration of its parameters to the way that leads to the optimal state of inventory system. In the inventory management, three approaches are common strategies which are fixed order quantity, fixed time interval, and forecasting methods [35]. In the first strategy, inventory level should be reviewed continuously until it reaches below predetermined quantity (reorder point) then order would be organized with fixed quantity of inventory. In the second strategy, reviewing period is a fixed time interval, but quantity order is variable for each order that is based on consumption rate. While first strategy needs more effort for continuously reviewing of inventory level, the second strategy is easier to handle, but the risk of shortage in fixed interval strategy is more in comparison to fixed order quantity. So, due to the mitigation of shortage risk in fixed interval strategy, order quantity is slightly more than fixed order quantity [35].In the third strategy, reviewing period is fixed as second policy, but reorder quantity is based on forecasting of future demands. As the mentioned strategies are fundamental in inventory management literature and are also common among retailers of office furniture, in this study three policies based on fixed order quantity, fixed time interval, and demand forecasting are developed as follows.(1)
Continuous reviewing with economic quantity order.
(2)
Periodic reviewing with order quantity based on demand confidence interval during the consumption period.
(3)
Periodic reviewing with order quantity based on forecasting of future demands.Policy 1 (Fixed Quantity). Based on this policy, each demand will be satisfied if there is sufficient inventory in hand. After satisfaction of each demand, inventory in hand will be checked to see if the inventory level reaches to reorder point (
R
p
). If inventory level has reached to reorder point, economic order quantity (
Q
j
t
) would be organized, otherwise system waits for next demand. Economic order quantity is derived by Wilson formula [36]. If there is not sufficient inventory to satisfy arrived demand, quantity of demand is considered as lost sale and lost profit treated as cost. If no order has been organized, system reorder inventory otherwise waits for arriving of organized order according to adjusted lead time (L). Logic of this policy is visible in Figure 1(a) and is labelled as FQ policy.(a)
FQ policy
(b)
FI policy
(c)
DF policyPolicy 2 (Fixed Interval). In this policy, criterion for reordering is fixed time interval that is known as planning period. After this period, inventory level would be examined and reorder will be organized on condition that inventory level has reached to reorder point (R
p). For more realistic consideration, order quantity (Q
j
t) is calculated based on demand cumulative distribution function in planning period. For example, in the inventory system, that is planned for 10% lost sale with demand which is distributed based on exponential distribution function, it should be ordered as much as cumulative probability of exponential distribution equals to 0.9. In this policy, inventory level would be updated after lead time (L). Logic of this policy is shown in Figure 1(b) and is labelled as FI policy.Policy 3 (Demand Forecasting). This policy is similar to policy 2 with some differences. In policy 2, probability distribution function of demand in planning period is estimated by historical data and then reorder is organized based on service level (cumulative probability of demand satisfaction). But, in this policy after fixed interval, forecasting of future demand will be performed. If inventory level has reached to reorder point (R
p), order would be organized and inventory level would be updated after lead time (L), otherwise only information would be updated. If forecasted demand is less than inventory capacity, order quantity (Q
j
t) would be as much as forecasted quantity otherwise inventory order will be as much as inventory capacity for each product. As demands follow time series with autocorrelation, forecasting is implemented by exponential smoothing method [37].In this procedure,A
t is defined as demand estimation of tth planning period. This term is calculated according to real demand of period t indicated by D
t and estimation of demand in previous period which is indicated by A
t
-
1 plus trend adjustment value of previous period which is indicated by T
t
-
1. Impact of current real demand is considered by parameter α that is defined between zero and one. Trend adjustment value is derived based on A
t and A
t
-
1 plus trend adjustment value of previous period with consideration of β as weight parameter which is defined between zero and one. Equations related to A
t and T
t are expressed in (4) and (5), respectively. Finally, demand estimation plus trend adjustment form of future demand forecasting is indicated by F
t
+
1 which is described in (6). Logic of this policy is shown in Figure 1(c) and is labelled as DF policy. Consider
(4)
A
t
=
α
×
D
t
+
(
1
-
α
)
×
(
A
t
-
1
+
T
t
-
1
)
,
(5)
T
t
=
β
×
(
A
t
-
A
t
-
1
)
+
(
1
-
β
)
×
T
t
-
1
,
(6)
F
t
+
1
=
A
t
+
T
t
.


### 4.2. Surrogate Modelling

Surrogate model is designed for optimization of simulated inventory system. In this framework, design of experiments is responsible for producing different scenarios. Each scenario has four criteria including cost, service level, average of inventory position, and robustness against demand fluctuation. Importance of each criterion is determined by PCA and finally these scenarios are ranked by MADM technique as multiple criteria decision making tool.


#### 4.2.1. Design of Experimentation

Although DOE roots back to the statistical quality control, nowadays it is a powerful tool for analysing complex systems. DOE is statistical method and organizes structured experiments with several factors [38]. This method not only determines effect of each factor on response variable, but also considers multiple effects simultaneously. DOE is extended technique that includes several designs such as full factorial, fractional factorial, and nested design. Specific design varies for each problem and should be selected based on problem condition. Our problem consists of tree effective factors, namely, inventory policy, reorder point, and lead time, these factors have nonlinear effect on objective function values. As the aim of DOE in this research is producing decision space with few factors (in this research, there are 3 factors), using full factorial design is preferable. On the other hand, because nonlinear behaviour of response variables (i.e., objective functions) centre point is considered. Each factor has high level “+,” low level “−,” and centre point “0.” Factors and levels configurations are shown in Table 2 while Table 3 shows experimental design. In both of these tables, minus sign is representative of low level of related factors while plus sign is indicator of high level of factors and zero determines centre point. For example, in the case of lead time (L), 3 is low level, 7 is high level, and 5 is centre point. Lead time is considered with days as measuring unit and begins from time order is organized until the arrival of inventory. Reorder point is labelled as R
p and is intended to satisfy product demand during lead time (L). Eventually, inventory policy is labelled as IP and refers to inventory policies that were described in Section 4.1.2. The mentioned design includes 27 runs and there are four response variables: cost, service level, average of inventory position, and robustness. In fact, first three response variables are expected values of defined objective functions (i.e., resp., (1), (2), and (3)), but robustness is a defined criterion which guarantees that optimal solution remains valid in situation where demand variation and intensification increase. For evaluation of robustness, demands are multiplied by random variable with normal distribution (μ
=
1, σ
=
1.5). This random variable is also restricted to positive values to prevent negative demands. Expected value of 4000 generated numbers with mentioned constraint is equal to 1.48; this result implies that multiplied random number increases both demand variation and intensification, simultaneously. Then, the percent of decrease in service level is considered as robustness for all 27 scenarios produced by DOE. It is obvious that a lower decrease in service level is desirable.Table 2
DOE configuration.
Products
Factors
IP
R
t
L
−
0
+
−
0
+
−
0
+
A
FI
FQ
DF
10
20
30
3
5
7
B
FI
FQ
DF
15
30
45
3
5
7
C
FI
FQ
DF
15
30
45
3
5
7
D
FI
FQ
DF
15
30
45
3
5
7Table 3
Experimental design.
Number of trial
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
R
p
−
−
−
−
−
−
−
−
−
0
0
0
0
0
0
0
0
0
+
+
+
+
+
+
+
+
+
L
−
−
−
0
0
0
+
+
+
−
−
−
0
0
0
+
+
+
−
−
−
0
0
0
+
+
+
IP
−
0
+
−
0
+
−
0
+
−
0
+
−
0
+
−
0
+
−
0
+
−
0
+
−
0
+


#### 4.2.2. Weight Assignment Based on PCA Method

PCA is the abbreviation of “principal component analysis” and is known as powerful data-driven method successfully applied in several areas [2]. PCA analysis can determine key variables that govern most of the variability. This technique helps to remove less important variables to reduce dimension of data sets. PCA considers two criteria for each variable in data sets: the first is variability and the second is correlation with other variables. So, input for this method consists of variance covariance matrix [39, 40]. Although other methods like entropy is applicable in multiple attribute decision making problem, PCA advantage is considerable because variables correlations are not ignored. As our problem is constructed of highly correlated variables, PCA is the preferable choice. For correlation justification, remember that improvement of inventory position criterion causes negative influence on shortage cost or improvement of service level criterion causes negative influence on inventory in hand and holding cost. In these circumstances, there is no reason for applying variation based methods like entropy, so in this paper PCA is applied for weight assignment. These weights are used by MADM technique (described in Section 4.2.3) to improve quality of ranking. Because of relative utility of scenarios, linear scaleless method employed to scaleless data obtained from simulation.(1)
Calculate scaleless matrix ofD
=
(
d
1
,
d
2
,
…
,
d
p
)
n
×
p which d
i
j is scaleless value of ith scenario in the view point of jth criterion.
(2)
Calculate sample mean vector;d
-
=
(
d
-
1
,
d
-
2
,
…
,
d
-
p
)
p and covariance matrix based on (7) and (8), respectively,
(7)
d
-
k
=
1
n
∑
i
=
1
n
d
i
j
,
(8)
S
=
(
s
k
q
)
p
×
p
=
1
n
-
1
(
D
-
d
-
)
T
(
D
-
d
-
)
.
(3)
C
1
/
s
k
k is diagonal p
×
p matrix, kth diagonal element is 1
/
s
k
k, and sample correlation matrix (
R
) is calculated as mentioned in the following equation:
(9)
R
=
C
1
s
k
k
·
S
·
C
1
s
k
k
.
(4)
Solve the following equation, whereI
p is a p
×
p identity matrix:
(10)
|
R
-
λ
I
p
|
=
0
.
Solving (10) results in P individual eigenvalues regarding λ
1
≥
λ
2
≥
⋯
≥
λ
p and ∑
k
=
1
p
λ
k
=
P. These eigenvalues have P-specific eigenvectors, (
l
1
k
,
l
2
k
,
…
,
l
p
k
) and k
=
1
,
…
,
p. These vectors create principal components as expressed in the following equation:
(11)
PC
k
=
∑
q
=
1
p
l
q
k
d
q
j
.
(5)
Select sufficient number of principal components. In this problem, first,M components whose cumulative percentage of their eigenvalues (C
M) is greater than 95% are selected. Cumulative percentage is calculated based on the following equation:
(12)
C
m
=
∑
k
=
1
m
λ
k
∑
k
=
1
p
λ
k
=
∑
k
=
1
m
λ
k
P
.
(6)
Finally, weights vector that is indicated byZ is derived from the mathematical relation in (13) while W
k is equal to λ
k
/
P, if the entire elements of PC
k are positive, otherwise -
λ
k
/
P if entire component are negative
(13)
Z
=
∑
k
=
1
M
W
k
PC
k
.In other cases, positive or negative W
k should be defined in the order that final weight vector is positive. For more detailed discussion about W
k, refer to [41]. Finally, for each product, there are four Z parameters that determine relative importance for the following criteria: cost, service level, average of inventory position, and robustness; each element of Z vector is divided by the sum of vector elements for normalization.


#### 4.2.3. VIKOR Method

This method was developed by Opricovic and is the abbreviation of expression that in Serbian means “multiple criteria optimization and compromising solutions” [42]. VIKOR method is based on compromising solutions which is the result of Yu [43] and Zeleny [44] studies. This approach considers closeness to ideal solution. In comparison with other MADM methods, VIKOR is able to present compromising solutions and substitutes, these solutions with best one in the case of necessity [45]. VIKOR method is applicable for problems with multiple discrete criteria. Generally, this method can be considered as ranking tool for scenarios which are generated by DOE. Ranking is based on L-P metric function according to (14) with definition of f
i
j as the value of ith criterion for jth scenario. In construction of L-P metric function, f
i
* indicates best value of all scenarios in the view point of ith criterion and f
i
- represents the worst value of all scenarios with consideration of ith criterion. w
i is the weight parameter for ith criterion and is derived by PCA method. Consider
(14)
L
p
=
{
∑
i
=
1
n
[
w
i
(
f
i
*
-
f
i
j
)
(
f
i
*
-
f
i
-
)
]
p
}
1
/
p
1
≪
P
≪
∞
;
j
=
1,2
,
3
,
…
,
J
,
L
∞
=
max
⁡
i
[
w
i
(
f
i
*
-
f
i
j
)
(
f
i
*
-
f
i
-
)
]
VIKOR method only involves L
1 and L
∞ which are known as S
j and R
j, respectively. S and R index refer to scenarios so both of these functions are regarded for each scenario. VIKOR steps are as follows.(1)
Calculatef
i
- and f
i
* for entire of criteria.
(2)
CalculateS
j and R
j for entire of scenarios.
(3)
With the assumption thatS
* and R
* are the smallest values among all of S
j and R
j, respectively, and also S
- and R
- are the greatest values in comparison with all of S
j and R
j, respectively, Q
j is calculated for jth scenario based on (15). In this equation, v is strategic weight that determines the importance of individual criterion (S
j) against group importance of criteria (R
j) for each scenario. Consider
(15)
Q
j
=
v
(
S
j
-
S
*
)
S
-
-
S
*
+
(
1
-
v
)
(
R
j
-
R
*
)
R
-
-
R
*
.
(4)
Scenarios should be sorted in descending order based onS, R, and Q specifically that result in three different ranking lists.
(5)
According toQ, the scenario which is indicated by a
′ is the best scenario, if two conditions, namely, (C1) and (C2) are satisfied.(C1). Acceptable Advantage. This condition implies that the second best solution (a
′
′) should be far enough from the best solution (a
′) to be accepted as unique solution. The mathematical equation is defined in (16), considering J as the total number of scenarios. Consider
(16)
Q
(
a
′′
)
-
Q
(
a
′
)
≫
1
(
J
-
1
)
.(C2). Acceptable Stability. This condition validates a
′ as the best stable solution if a
′ is also the best solution based on S or R.In case that one or both of the mentioned conditions are not satisfied, compromised solutions are involved. VIKOR method suggests that in such case compromising solutions have equal value to the best solution and could be considered interchangeably by decision makers.If (C1) is not satisfied, scenariosa
′
,
a
′
′
,
…
,
a
M are considered. a
′ is the best solution and a
′
′
,
…
,
a
M are compromising solutions with priority of Q. on the other hand, a
′, a
′
′ are considered if (C2) is not satisfied. M superscript is determined by (17)
(17)
Q
(
a
M
)
-
Q
(
a
′
)
≫
1
(
J
-
1
)
.


## 4.1. Simulation Modelling

Simulation modelling of problem consists of two parts which are modelling of demand and modelling of inventory policies. In this paper, demand is nonhomogeneous Poisson process and three different inventory policies based on continues reviewing, periodic reviewing, and periodic reviewing with forecasting of future demand are considered.


### 4.1.1. Modelling of Demand

For customers demand modeling, multiresolution method is applied. Kuhl and Wilson [34] developed this method for simulation of nonhomogeneous Poisson process with trend and cyclic changes. This method estimates mean intensity function and the nonparametric nature of this method is one of the most important advantages in comparison with other methods. So, it is independent of statistical parameters and applicable in variety of problems. Furthermore, multiresolution method can support combination of multiple cyclic changes simultaneously. Another advantage of multiresolution is its capability in the modelling of nonsymmetric cyclic pattern. As in our case, demand has nonsymmetric pattern with cyclic changes and long-term trend, and multiresolution approach is a reasonable choice. More theoretical and application of used method are provided in [34].


### 4.1.2. Modelling of Inventory Policies

As the main effort of this paper is inventory system optimization, modelling of inventory policies plays an important role in this problem. In this problem, optimization is manipulated by selection of appropriate inventory policy and configuration of its parameters to the way that leads to the optimal state of inventory system. In the inventory management, three approaches are common strategies which are fixed order quantity, fixed time interval, and forecasting methods [35]. In the first strategy, inventory level should be reviewed continuously until it reaches below predetermined quantity (reorder point) then order would be organized with fixed quantity of inventory. In the second strategy, reviewing period is a fixed time interval, but quantity order is variable for each order that is based on consumption rate. While first strategy needs more effort for continuously reviewing of inventory level, the second strategy is easier to handle, but the risk of shortage in fixed interval strategy is more in comparison to fixed order quantity. So, due to the mitigation of shortage risk in fixed interval strategy, order quantity is slightly more than fixed order quantity [35].In the third strategy, reviewing period is fixed as second policy, but reorder quantity is based on forecasting of future demands. As the mentioned strategies are fundamental in inventory management literature and are also common among retailers of office furniture, in this study three policies based on fixed order quantity, fixed time interval, and demand forecasting are developed as follows.(1)
Continuous reviewing with economic quantity order.
(2)
Periodic reviewing with order quantity based on demand confidence interval during the consumption period.
(3)
Periodic reviewing with order quantity based on forecasting of future demands.Policy 1 (Fixed Quantity). Based on this policy, each demand will be satisfied if there is sufficient inventory in hand. After satisfaction of each demand, inventory in hand will be checked to see if the inventory level reaches to reorder point (
R
p
). If inventory level has reached to reorder point, economic order quantity (
Q
j
t
) would be organized, otherwise system waits for next demand. Economic order quantity is derived by Wilson formula [36]. If there is not sufficient inventory to satisfy arrived demand, quantity of demand is considered as lost sale and lost profit treated as cost. If no order has been organized, system reorder inventory otherwise waits for arriving of organized order according to adjusted lead time (L). Logic of this policy is visible in Figure 1(a) and is labelled as FQ policy.(a)
FQ policy
(b)
FI policy
(c)
DF policyPolicy 2 (Fixed Interval). In this policy, criterion for reordering is fixed time interval that is known as planning period. After this period, inventory level would be examined and reorder will be organized on condition that inventory level has reached to reorder point (R
p). For more realistic consideration, order quantity (Q
j
t) is calculated based on demand cumulative distribution function in planning period. For example, in the inventory system, that is planned for 10% lost sale with demand which is distributed based on exponential distribution function, it should be ordered as much as cumulative probability of exponential distribution equals to 0.9. In this policy, inventory level would be updated after lead time (L). Logic of this policy is shown in Figure 1(b) and is labelled as FI policy.Policy 3 (Demand Forecasting). This policy is similar to policy 2 with some differences. In policy 2, probability distribution function of demand in planning period is estimated by historical data and then reorder is organized based on service level (cumulative probability of demand satisfaction). But, in this policy after fixed interval, forecasting of future demand will be performed. If inventory level has reached to reorder point (R
p), order would be organized and inventory level would be updated after lead time (L), otherwise only information would be updated. If forecasted demand is less than inventory capacity, order quantity (Q
j
t) would be as much as forecasted quantity otherwise inventory order will be as much as inventory capacity for each product. As demands follow time series with autocorrelation, forecasting is implemented by exponential smoothing method [37].In this procedure,A
t is defined as demand estimation of tth planning period. This term is calculated according to real demand of period t indicated by D
t and estimation of demand in previous period which is indicated by A
t
-
1 plus trend adjustment value of previous period which is indicated by T
t
-
1. Impact of current real demand is considered by parameter α that is defined between zero and one. Trend adjustment value is derived based on A
t and A
t
-
1 plus trend adjustment value of previous period with consideration of β as weight parameter which is defined between zero and one. Equations related to A
t and T
t are expressed in (4) and (5), respectively. Finally, demand estimation plus trend adjustment form of future demand forecasting is indicated by F
t
+
1 which is described in (6). Logic of this policy is shown in Figure 1(c) and is labelled as DF policy. Consider
(4)
A
t
=
α
×
D
t
+
(
1
-
α
)
×
(
A
t
-
1
+
T
t
-
1
)
,
(5)
T
t
=
β
×
(
A
t
-
A
t
-
1
)
+
(
1
-
β
)
×
T
t
-
1
,
(6)
F
t
+
1
=
A
t
+
T
t
.


## 4.1.1. Modelling of Demand

For customers demand modeling, multiresolution method is applied. Kuhl and Wilson [34] developed this method for simulation of nonhomogeneous Poisson process with trend and cyclic changes. This method estimates mean intensity function and the nonparametric nature of this method is one of the most important advantages in comparison with other methods. So, it is independent of statistical parameters and applicable in variety of problems. Furthermore, multiresolution method can support combination of multiple cyclic changes simultaneously. Another advantage of multiresolution is its capability in the modelling of nonsymmetric cyclic pattern. As in our case, demand has nonsymmetric pattern with cyclic changes and long-term trend, and multiresolution approach is a reasonable choice. More theoretical and application of used method are provided in [34].


## 4.1.2. Modelling of Inventory Policies

As the main effort of this paper is inventory system optimization, modelling of inventory policies plays an important role in this problem. In this problem, optimization is manipulated by selection of appropriate inventory policy and configuration of its parameters to the way that leads to the optimal state of inventory system. In the inventory management, three approaches are common strategies which are fixed order quantity, fixed time interval, and forecasting methods [35]. In the first strategy, inventory level should be reviewed continuously until it reaches below predetermined quantity (reorder point) then order would be organized with fixed quantity of inventory. In the second strategy, reviewing period is a fixed time interval, but quantity order is variable for each order that is based on consumption rate. While first strategy needs more effort for continuously reviewing of inventory level, the second strategy is easier to handle, but the risk of shortage in fixed interval strategy is more in comparison to fixed order quantity. So, due to the mitigation of shortage risk in fixed interval strategy, order quantity is slightly more than fixed order quantity [35].In the third strategy, reviewing period is fixed as second policy, but reorder quantity is based on forecasting of future demands. As the mentioned strategies are fundamental in inventory management literature and are also common among retailers of office furniture, in this study three policies based on fixed order quantity, fixed time interval, and demand forecasting are developed as follows.(1)
Continuous reviewing with economic quantity order.
(2)
Periodic reviewing with order quantity based on demand confidence interval during the consumption period.
(3)
Periodic reviewing with order quantity based on forecasting of future demands.Policy 1 (Fixed Quantity). Based on this policy, each demand will be satisfied if there is sufficient inventory in hand. After satisfaction of each demand, inventory in hand will be checked to see if the inventory level reaches to reorder point (
R
p
). If inventory level has reached to reorder point, economic order quantity (
Q
j
t
) would be organized, otherwise system waits for next demand. Economic order quantity is derived by Wilson formula [36]. If there is not sufficient inventory to satisfy arrived demand, quantity of demand is considered as lost sale and lost profit treated as cost. If no order has been organized, system reorder inventory otherwise waits for arriving of organized order according to adjusted lead time (L). Logic of this policy is visible in Figure 1(a) and is labelled as FQ policy.(a)
FQ policy
(b)
FI policy
(c)
DF policyPolicy 2 (Fixed Interval). In this policy, criterion for reordering is fixed time interval that is known as planning period. After this period, inventory level would be examined and reorder will be organized on condition that inventory level has reached to reorder point (R
p). For more realistic consideration, order quantity (Q
j
t) is calculated based on demand cumulative distribution function in planning period. For example, in the inventory system, that is planned for 10% lost sale with demand which is distributed based on exponential distribution function, it should be ordered as much as cumulative probability of exponential distribution equals to 0.9. In this policy, inventory level would be updated after lead time (L). Logic of this policy is shown in Figure 1(b) and is labelled as FI policy.Policy 3 (Demand Forecasting). This policy is similar to policy 2 with some differences. In policy 2, probability distribution function of demand in planning period is estimated by historical data and then reorder is organized based on service level (cumulative probability of demand satisfaction). But, in this policy after fixed interval, forecasting of future demand will be performed. If inventory level has reached to reorder point (R
p), order would be organized and inventory level would be updated after lead time (L), otherwise only information would be updated. If forecasted demand is less than inventory capacity, order quantity (Q
j
t) would be as much as forecasted quantity otherwise inventory order will be as much as inventory capacity for each product. As demands follow time series with autocorrelation, forecasting is implemented by exponential smoothing method [37].In this procedure,A
t is defined as demand estimation of tth planning period. This term is calculated according to real demand of period t indicated by D
t and estimation of demand in previous period which is indicated by A
t
-
1 plus trend adjustment value of previous period which is indicated by T
t
-
1. Impact of current real demand is considered by parameter α that is defined between zero and one. Trend adjustment value is derived based on A
t and A
t
-
1 plus trend adjustment value of previous period with consideration of β as weight parameter which is defined between zero and one. Equations related to A
t and T
t are expressed in (4) and (5), respectively. Finally, demand estimation plus trend adjustment form of future demand forecasting is indicated by F
t
+
1 which is described in (6). Logic of this policy is shown in Figure 1(c) and is labelled as DF policy. Consider
(4)
A
t
=
α
×
D
t
+
(
1
-
α
)
×
(
A
t
-
1
+
T
t
-
1
)
,
(5)
T
t
=
β
×
(
A
t
-
A
t
-
1
)
+
(
1
-
β
)
×
T
t
-
1
,
(6)
F
t
+
1
=
A
t
+
T
t
.


## 4.2. Surrogate Modelling

Surrogate model is designed for optimization of simulated inventory system. In this framework, design of experiments is responsible for producing different scenarios. Each scenario has four criteria including cost, service level, average of inventory position, and robustness against demand fluctuation. Importance of each criterion is determined by PCA and finally these scenarios are ranked by MADM technique as multiple criteria decision making tool.


### 4.2.1. Design of Experimentation

Although DOE roots back to the statistical quality control, nowadays it is a powerful tool for analysing complex systems. DOE is statistical method and organizes structured experiments with several factors [38]. This method not only determines effect of each factor on response variable, but also considers multiple effects simultaneously. DOE is extended technique that includes several designs such as full factorial, fractional factorial, and nested design. Specific design varies for each problem and should be selected based on problem condition. Our problem consists of tree effective factors, namely, inventory policy, reorder point, and lead time, these factors have nonlinear effect on objective function values. As the aim of DOE in this research is producing decision space with few factors (in this research, there are 3 factors), using full factorial design is preferable. On the other hand, because nonlinear behaviour of response variables (i.e., objective functions) centre point is considered. Each factor has high level “+,” low level “−,” and centre point “0.” Factors and levels configurations are shown in Table 2 while Table 3 shows experimental design. In both of these tables, minus sign is representative of low level of related factors while plus sign is indicator of high level of factors and zero determines centre point. For example, in the case of lead time (L), 3 is low level, 7 is high level, and 5 is centre point. Lead time is considered with days as measuring unit and begins from time order is organized until the arrival of inventory. Reorder point is labelled as R
p and is intended to satisfy product demand during lead time (L). Eventually, inventory policy is labelled as IP and refers to inventory policies that were described in Section 4.1.2. The mentioned design includes 27 runs and there are four response variables: cost, service level, average of inventory position, and robustness. In fact, first three response variables are expected values of defined objective functions (i.e., resp., (1), (2), and (3)), but robustness is a defined criterion which guarantees that optimal solution remains valid in situation where demand variation and intensification increase. For evaluation of robustness, demands are multiplied by random variable with normal distribution (μ
=
1, σ
=
1.5). This random variable is also restricted to positive values to prevent negative demands. Expected value of 4000 generated numbers with mentioned constraint is equal to 1.48; this result implies that multiplied random number increases both demand variation and intensification, simultaneously. Then, the percent of decrease in service level is considered as robustness for all 27 scenarios produced by DOE. It is obvious that a lower decrease in service level is desirable.Table 2
DOE configuration.
Products
Factors
IP
R
t
L
−
0
+
−
0
+
−
0
+
A
FI
FQ
DF
10
20
30
3
5
7
B
FI
FQ
DF
15
30
45
3
5
7
C
FI
FQ
DF
15
30
45
3
5
7
D
FI
FQ
DF
15
30
45
3
5
7Table 3
Experimental design.
Number of trial
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
R
p
−
−
−
−
−
−
−
−
−
0
0
0
0
0
0
0
0
0
+
+
+
+
+
+
+
+
+
L
−
−
−
0
0
0
+
+
+
−
−
−
0
0
0
+
+
+
−
−
−
0
0
0
+
+
+
IP
−
0
+
−
0
+
−
0
+
−
0
+
−
0
+
−
0
+
−
0
+
−
0
+
−
0
+


### 4.2.2. Weight Assignment Based on PCA Method

PCA is the abbreviation of “principal component analysis” and is known as powerful data-driven method successfully applied in several areas [2]. PCA analysis can determine key variables that govern most of the variability. This technique helps to remove less important variables to reduce dimension of data sets. PCA considers two criteria for each variable in data sets: the first is variability and the second is correlation with other variables. So, input for this method consists of variance covariance matrix [39, 40]. Although other methods like entropy is applicable in multiple attribute decision making problem, PCA advantage is considerable because variables correlations are not ignored. As our problem is constructed of highly correlated variables, PCA is the preferable choice. For correlation justification, remember that improvement of inventory position criterion causes negative influence on shortage cost or improvement of service level criterion causes negative influence on inventory in hand and holding cost. In these circumstances, there is no reason for applying variation based methods like entropy, so in this paper PCA is applied for weight assignment. These weights are used by MADM technique (described in Section 4.2.3) to improve quality of ranking. Because of relative utility of scenarios, linear scaleless method employed to scaleless data obtained from simulation.(1)
Calculate scaleless matrix ofD
=
(
d
1
,
d
2
,
…
,
d
p
)
n
×
p which d
i
j is scaleless value of ith scenario in the view point of jth criterion.
(2)
Calculate sample mean vector;d
-
=
(
d
-
1
,
d
-
2
,
…
,
d
-
p
)
p and covariance matrix based on (7) and (8), respectively,
(7)
d
-
k
=
1
n
∑
i
=
1
n
d
i
j
,
(8)
S
=
(
s
k
q
)
p
×
p
=
1
n
-
1
(
D
-
d
-
)
T
(
D
-
d
-
)
.
(3)
C
1
/
s
k
k is diagonal p
×
p matrix, kth diagonal element is 1
/
s
k
k, and sample correlation matrix (
R
) is calculated as mentioned in the following equation:
(9)
R
=
C
1
s
k
k
·
S
·
C
1
s
k
k
.
(4)
Solve the following equation, whereI
p is a p
×
p identity matrix:
(10)
|
R
-
λ
I
p
|
=
0
.
Solving (10) results in P individual eigenvalues regarding λ
1
≥
λ
2
≥
⋯
≥
λ
p and ∑
k
=
1
p
λ
k
=
P. These eigenvalues have P-specific eigenvectors, (
l
1
k
,
l
2
k
,
…
,
l
p
k
) and k
=
1
,
…
,
p. These vectors create principal components as expressed in the following equation:
(11)
PC
k
=
∑
q
=
1
p
l
q
k
d
q
j
.
(5)
Select sufficient number of principal components. In this problem, first,M components whose cumulative percentage of their eigenvalues (C
M) is greater than 95% are selected. Cumulative percentage is calculated based on the following equation:
(12)
C
m
=
∑
k
=
1
m
λ
k
∑
k
=
1
p
λ
k
=
∑
k
=
1
m
λ
k
P
.
(6)
Finally, weights vector that is indicated byZ is derived from the mathematical relation in (13) while W
k is equal to λ
k
/
P, if the entire elements of PC
k are positive, otherwise -
λ
k
/
P if entire component are negative
(13)
Z
=
∑
k
=
1
M
W
k
PC
k
.In other cases, positive or negative W
k should be defined in the order that final weight vector is positive. For more detailed discussion about W
k, refer to [41]. Finally, for each product, there are four Z parameters that determine relative importance for the following criteria: cost, service level, average of inventory position, and robustness; each element of Z vector is divided by the sum of vector elements for normalization.


### 4.2.3. VIKOR Method

This method was developed by Opricovic and is the abbreviation of expression that in Serbian means “multiple criteria optimization and compromising solutions” [42]. VIKOR method is based on compromising solutions which is the result of Yu [43] and Zeleny [44] studies. This approach considers closeness to ideal solution. In comparison with other MADM methods, VIKOR is able to present compromising solutions and substitutes, these solutions with best one in the case of necessity [45]. VIKOR method is applicable for problems with multiple discrete criteria. Generally, this method can be considered as ranking tool for scenarios which are generated by DOE. Ranking is based on L-P metric function according to (14) with definition of f
i
j as the value of ith criterion for jth scenario. In construction of L-P metric function, f
i
* indicates best value of all scenarios in the view point of ith criterion and f
i
- represents the worst value of all scenarios with consideration of ith criterion. w
i is the weight parameter for ith criterion and is derived by PCA method. Consider
(14)
L
p
=
{
∑
i
=
1
n
[
w
i
(
f
i
*
-
f
i
j
)
(
f
i
*
-
f
i
-
)
]
p
}
1
/
p
1
≪
P
≪
∞
;
j
=
1,2
,
3
,
…
,
J
,
L
∞
=
max
⁡
i
[
w
i
(
f
i
*
-
f
i
j
)
(
f
i
*
-
f
i
-
)
]
VIKOR method only involves L
1 and L
∞ which are known as S
j and R
j, respectively. S and R index refer to scenarios so both of these functions are regarded for each scenario. VIKOR steps are as follows.(1)
Calculatef
i
- and f
i
* for entire of criteria.
(2)
CalculateS
j and R
j for entire of scenarios.
(3)
With the assumption thatS
* and R
* are the smallest values among all of S
j and R
j, respectively, and also S
- and R
- are the greatest values in comparison with all of S
j and R
j, respectively, Q
j is calculated for jth scenario based on (15). In this equation, v is strategic weight that determines the importance of individual criterion (S
j) against group importance of criteria (R
j) for each scenario. Consider
(15)
Q
j
=
v
(
S
j
-
S
*
)
S
-
-
S
*
+
(
1
-
v
)
(
R
j
-
R
*
)
R
-
-
R
*
.
(4)
Scenarios should be sorted in descending order based onS, R, and Q specifically that result in three different ranking lists.
(5)
According toQ, the scenario which is indicated by a
′ is the best scenario, if two conditions, namely, (C1) and (C2) are satisfied.(C1). Acceptable Advantage. This condition implies that the second best solution (a
′
′) should be far enough from the best solution (a
′) to be accepted as unique solution. The mathematical equation is defined in (16), considering J as the total number of scenarios. Consider
(16)
Q
(
a
′′
)
-
Q
(
a
′
)
≫
1
(
J
-
1
)
.(C2). Acceptable Stability. This condition validates a
′ as the best stable solution if a
′ is also the best solution based on S or R.In case that one or both of the mentioned conditions are not satisfied, compromised solutions are involved. VIKOR method suggests that in such case compromising solutions have equal value to the best solution and could be considered interchangeably by decision makers.If (C1) is not satisfied, scenariosa
′
,
a
′
′
,
…
,
a
M are considered. a
′ is the best solution and a
′
′
,
…
,
a
M are compromising solutions with priority of Q. on the other hand, a
′, a
′
′ are considered if (C2) is not satisfied. M superscript is determined by (17)
(17)
Q
(
a
M
)
-
Q
(
a
′
)
≫
1
(
J
-
1
)
.


## 4.2.1. Design of Experimentation

Although DOE roots back to the statistical quality control, nowadays it is a powerful tool for analysing complex systems. DOE is statistical method and organizes structured experiments with several factors [38]. This method not only determines effect of each factor on response variable, but also considers multiple effects simultaneously. DOE is extended technique that includes several designs such as full factorial, fractional factorial, and nested design. Specific design varies for each problem and should be selected based on problem condition. Our problem consists of tree effective factors, namely, inventory policy, reorder point, and lead time, these factors have nonlinear effect on objective function values. As the aim of DOE in this research is producing decision space with few factors (in this research, there are 3 factors), using full factorial design is preferable. On the other hand, because nonlinear behaviour of response variables (i.e., objective functions) centre point is considered. Each factor has high level “+,” low level “−,” and centre point “0.” Factors and levels configurations are shown in Table 2 while Table 3 shows experimental design. In both of these tables, minus sign is representative of low level of related factors while plus sign is indicator of high level of factors and zero determines centre point. For example, in the case of lead time (L), 3 is low level, 7 is high level, and 5 is centre point. Lead time is considered with days as measuring unit and begins from time order is organized until the arrival of inventory. Reorder point is labelled as R
p and is intended to satisfy product demand during lead time (L). Eventually, inventory policy is labelled as IP and refers to inventory policies that were described in Section 4.1.2. The mentioned design includes 27 runs and there are four response variables: cost, service level, average of inventory position, and robustness. In fact, first three response variables are expected values of defined objective functions (i.e., resp., (1), (2), and (3)), but robustness is a defined criterion which guarantees that optimal solution remains valid in situation where demand variation and intensification increase. For evaluation of robustness, demands are multiplied by random variable with normal distribution (μ
=
1, σ
=
1.5). This random variable is also restricted to positive values to prevent negative demands. Expected value of 4000 generated numbers with mentioned constraint is equal to 1.48; this result implies that multiplied random number increases both demand variation and intensification, simultaneously. Then, the percent of decrease in service level is considered as robustness for all 27 scenarios produced by DOE. It is obvious that a lower decrease in service level is desirable.Table 2
DOE configuration.
Products
Factors
IP
R
t
L
−
0
+
−
0
+
−
0
+
A
FI
FQ
DF
10
20
30
3
5
7
B
FI
FQ
DF
15
30
45
3
5
7
C
FI
FQ
DF
15
30
45
3
5
7
D
FI
FQ
DF
15
30
45
3
5
7Table 3
Experimental design.
Number of trial
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
R
p
−
−
−
−
−
−
−
−
−
0
0
0
0
0
0
0
0
0
+
+
+
+
+
+
+
+
+
L
−
−
−
0
0
0
+
+
+
−
−
−
0
0
0
+
+
+
−
−
−
0
0
0
+
+
+
IP
−
0
+
−
0
+
−
0
+
−
0
+
−
0
+
−
0
+
−
0
+
−
0
+
−
0
+


## 4.2.2. Weight Assignment Based on PCA Method

PCA is the abbreviation of “principal component analysis” and is known as powerful data-driven method successfully applied in several areas [2]. PCA analysis can determine key variables that govern most of the variability. This technique helps to remove less important variables to reduce dimension of data sets. PCA considers two criteria for each variable in data sets: the first is variability and the second is correlation with other variables. So, input for this method consists of variance covariance matrix [39, 40]. Although other methods like entropy is applicable in multiple attribute decision making problem, PCA advantage is considerable because variables correlations are not ignored. As our problem is constructed of highly correlated variables, PCA is the preferable choice. For correlation justification, remember that improvement of inventory position criterion causes negative influence on shortage cost or improvement of service level criterion causes negative influence on inventory in hand and holding cost. In these circumstances, there is no reason for applying variation based methods like entropy, so in this paper PCA is applied for weight assignment. These weights are used by MADM technique (described in Section 4.2.3) to improve quality of ranking. Because of relative utility of scenarios, linear scaleless method employed to scaleless data obtained from simulation.(1)
Calculate scaleless matrix ofD
=
(
d
1
,
d
2
,
…
,
d
p
)
n
×
p which d
i
j is scaleless value of ith scenario in the view point of jth criterion.
(2)
Calculate sample mean vector;d
-
=
(
d
-
1
,
d
-
2
,
…
,
d
-
p
)
p and covariance matrix based on (7) and (8), respectively,
(7)
d
-
k
=
1
n
∑
i
=
1
n
d
i
j
,
(8)
S
=
(
s
k
q
)
p
×
p
=
1
n
-
1
(
D
-
d
-
)
T
(
D
-
d
-
)
.
(3)
C
1
/
s
k
k is diagonal p
×
p matrix, kth diagonal element is 1
/
s
k
k, and sample correlation matrix (
R
) is calculated as mentioned in the following equation:
(9)
R
=
C
1
s
k
k
·
S
·
C
1
s
k
k
.
(4)
Solve the following equation, whereI
p is a p
×
p identity matrix:
(10)
|
R
-
λ
I
p
|
=
0
.
Solving (10) results in P individual eigenvalues regarding λ
1
≥
λ
2
≥
⋯
≥
λ
p and ∑
k
=
1
p
λ
k
=
P. These eigenvalues have P-specific eigenvectors, (
l
1
k
,
l
2
k
,
…
,
l
p
k
) and k
=
1
,
…
,
p. These vectors create principal components as expressed in the following equation:
(11)
PC
k
=
∑
q
=
1
p
l
q
k
d
q
j
.
(5)
Select sufficient number of principal components. In this problem, first,M components whose cumulative percentage of their eigenvalues (C
M) is greater than 95% are selected. Cumulative percentage is calculated based on the following equation:
(12)
C
m
=
∑
k
=
1
m
λ
k
∑
k
=
1
p
λ
k
=
∑
k
=
1
m
λ
k
P
.
(6)
Finally, weights vector that is indicated byZ is derived from the mathematical relation in (13) while W
k is equal to λ
k
/
P, if the entire elements of PC
k are positive, otherwise -
λ
k
/
P if entire component are negative
(13)
Z
=
∑
k
=
1
M
W
k
PC
k
.In other cases, positive or negative W
k should be defined in the order that final weight vector is positive. For more detailed discussion about W
k, refer to [41]. Finally, for each product, there are four Z parameters that determine relative importance for the following criteria: cost, service level, average of inventory position, and robustness; each element of Z vector is divided by the sum of vector elements for normalization.


## 4.2.3. VIKOR Method

This method was developed by Opricovic and is the abbreviation of expression that in Serbian means “multiple criteria optimization and compromising solutions” [42]. VIKOR method is based on compromising solutions which is the result of Yu [43] and Zeleny [44] studies. This approach considers closeness to ideal solution. In comparison with other MADM methods, VIKOR is able to present compromising solutions and substitutes, these solutions with best one in the case of necessity [45]. VIKOR method is applicable for problems with multiple discrete criteria. Generally, this method can be considered as ranking tool for scenarios which are generated by DOE. Ranking is based on L-P metric function according to (14) with definition of f
i
j as the value of ith criterion for jth scenario. In construction of L-P metric function, f
i
* indicates best value of all scenarios in the view point of ith criterion and f
i
- represents the worst value of all scenarios with consideration of ith criterion. w
i is the weight parameter for ith criterion and is derived by PCA method. Consider
(14)
L
p
=
{
∑
i
=
1
n
[
w
i
(
f
i
*
-
f
i
j
)
(
f
i
*
-
f
i
-
)
]
p
}
1
/
p
1
≪
P
≪
∞
;
j
=
1,2
,
3
,
…
,
J
,
L
∞
=
max
⁡
i
[
w
i
(
f
i
*
-
f
i
j
)
(
f
i
*
-
f
i
-
)
]
VIKOR method only involves L
1 and L
∞ which are known as S
j and R
j, respectively. S and R index refer to scenarios so both of these functions are regarded for each scenario. VIKOR steps are as follows.(1)
Calculatef
i
- and f
i
* for entire of criteria.
(2)
CalculateS
j and R
j for entire of scenarios.
(3)
With the assumption thatS
* and R
* are the smallest values among all of S
j and R
j, respectively, and also S
- and R
- are the greatest values in comparison with all of S
j and R
j, respectively, Q
j is calculated for jth scenario based on (15). In this equation, v is strategic weight that determines the importance of individual criterion (S
j) against group importance of criteria (R
j) for each scenario. Consider
(15)
Q
j
=
v
(
S
j
-
S
*
)
S
-
-
S
*
+
(
1
-
v
)
(
R
j
-
R
*
)
R
-
-
R
*
.
(4)
Scenarios should be sorted in descending order based onS, R, and Q specifically that result in three different ranking lists.
(5)
According toQ, the scenario which is indicated by a
′ is the best scenario, if two conditions, namely, (C1) and (C2) are satisfied.(C1). Acceptable Advantage. This condition implies that the second best solution (a
′
′) should be far enough from the best solution (a
′) to be accepted as unique solution. The mathematical equation is defined in (16), considering J as the total number of scenarios. Consider
(16)
Q
(
a
′′
)
-
Q
(
a
′
)
≫
1
(
J
-
1
)
.(C2). Acceptable Stability. This condition validates a
′ as the best stable solution if a
′ is also the best solution based on S or R.In case that one or both of the mentioned conditions are not satisfied, compromised solutions are involved. VIKOR method suggests that in such case compromising solutions have equal value to the best solution and could be considered interchangeably by decision makers.If (C1) is not satisfied, scenariosa
′
,
a
′
′
,
…
,
a
M are considered. a
′ is the best solution and a
′
′
,
…
,
a
M are compromising solutions with priority of Q. on the other hand, a
′, a
′
′ are considered if (C2) is not satisfied. M superscript is determined by (17)
(17)
Q
(
a
M
)
-
Q
(
a
′
)
≫
1
(
J
-
1
)
.


## 5. Numerical Result

Simulation process was executed on laptop with 1.8 GH CPU and 4 GB of RAM. Running of each replication with animation and maximum speed takes 7 minutes and 35 seconds. As simulation of each scenario contains 10 replications and there are 27 scenarios, the total simulation time will be 2046 minutes and 36 seconds.


### 5.1. Multiresolution Result for Demand Modelling

Result of demand modelling is presented for productB because this analysis for other products is similar. For this problem, polynomial functions are applied for yearly and monthly demand rates. Equations (18) and (19) are estimated by nonlinear regression model. In these functions t
y and t
m are time variables
(18)
f
(
t
y
)
=
0.008942
×
t
y
3
-
0.03884
×
t
y
2
+
0.1867
×
t
y
-
0.08018
and 95% confidence interval is presented in Table 4 for each coefficient of (18).Table 4
95% confidence interval for estimated coefficient off
(
t
y
).
‍Coefficients
Lower bund
Upper bund
0.008942
0.004835
0.01305
−0.03884
−0.07604
−0.001652
0.1867
0.08654
0.287
−0.08018
−0.1568
−0.003523In this table, coefficients of (18) are visible in the first column while lower bund and upper bund of each coefficient are presented in the second and third columns, respectively. There is no interval including zero, so it is concluded that coefficient estimation is valid.Estimation error forf
(
t
y
) is reported based on the sum of square errors (SSE) which is 1.504
×
10
-
6. Consider
(19)
f
(
t
m
)
=
4.768
×
10
-
5
×
t
m
5
-
0.001692
×
t
m
4
+
0.0214
×
t
m
3
-
0.1159
×
t
m
2
+
0.3474
×
t
m
-
0.2286
.
Monthly demand rate function that is indicated by f
(
t
m
) and its coefficients confidence interval are presented in Table 5. SEE is equal to 0.0005825. Based on local business information, total demand of product B is estimated to be about 6822 units in five years. According to the information about product B, Figure 2 shows final multiresolution result for this product.Table 5
95% confidence interval for estimated coefficient off
(
t
m
).
‍Coefficients
Lower bund
Upper bund
4.768 × 10−5
4.451 × 10−5
5.086 × 10−5
-
0.001692
−0.001795
−0.00159
0.0214
0.0202
0.0226
-
0.1159
−0.1221
−0.1097
0.3474
0.334
0.3608
-
0.2286
−0.2371
−0.22Figure 2
Comparison between real and estimated demand.


### 5.2. Simulation Parameters and Configuration

For simulation of problem, it is necessary to configure inventory policy parameters and variables. We have three different policies, namely, FQ, FI, and DF. Variables of FQ policy are reorder point (R
p), lead time (L), and economic quantity orders (Q
j
t). First two variables are configured based on Table 2 and third variable is derived from Wilson formula, according to average of demand, holding cost (C
h), and ordering cost (C
r). The calculated values are 57, 71, 72, and 72 for product A, B, C, and D, respectively.Orders quantity is unknown for FI policy and is calculated based on probability distribution function which is obtained by demand simulation of each product for 20 days’ time interval. Distribution fitting is performed by input analyser of Arena software and graphical result for productB is shown in Figure 3. Results for other products are reported in Table 6.Table 6
Estimation of distribution functions of demand and reorders quantity.
Product type
Distribution function
P value
Order quantity
A
278* BETA(1.13,3.28)
>0.15
100
B
GAMMA(67.8,1.18)
>0.15
92
C
EXPO(81.8)
>0.15
100
D
446* BETA(0.957,3.21)
>0.15
115Figure 3
Distribution function of customer demand for reorder interval of productB.In this table, fitted probability distribution function for each product is presented in second column whileP value of fitting which is greater than 0.15 for fitted distribution grantees goodness of fitting in 95% confidence level. Finally, according to fitted distribution functions, for FI policy, proper ordering amount of each product is presented in the fourth column.DF policy has two unknown parameters, namely,α and β. The first parameter is weight of current real demand in (4) and second one is the weight of current trend in (5). These parameters are adjusted by simulation of five years demand. Demand for the entire four products simultaneously is considered and an average of lost sale percentage is considered as response variable. Different values of response variables are reported in Table 7. According to the different values of α and β, the best values for α and β are 0.9 and 0.1, respectively, which results in minimum lost percentage in average.Table 7
Average of lost sale percentage for different values ofα and β.
α and β values
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
Response toα
0.380
0.356
0.334
0.320
0.306
0.292
0.284
0.276
0.268
Response toβ
0.296
0.302
0.306
0.306
0.306
0.304
0.305
0.304
0.302


### 5.3. Weight Assignment

As mentioned before, PCA is applied for weight assignment. For this purpose, DOE result is needed in the form of decision matrix that is shown in Table8. In this table, four criteria: cost, average of inventory position, service level, and robustness are evaluated for each scenario of product B. Decision matrix should be scaleless that is performed by linear method. Then, the result of eigenvalues is calculated and shown in Table 9. Also result of principal components is presented in Table 10.Table 8
Decision matrix for productB.
Scenario number
Cost
Average of inventory position
Service level %
Robustness %
1
32984000
43.52
76
6
2
29214710
39.7514
79
3
3
37114636
33.9623
69
2
4
34840720
41.0948
71
5
5
36820570
34.4738
72
6
6
39191837
33.1091
67
2
7
37037180
39.1012
68
4
8
40994680
29.8712
64
6
9
41407045
31.9678
63
1
10
28282930
52.4462
83
7
11
26422145
48.9243
85
5
12
31023829
41.6049
79
3
13
30395190
48.7146
79
6
14
35393680
39.5812
74
5
15
32838055
39.5646
77
3
16
33285615
46.3341
76
6
17
39652020
33.5668
66
7
18
35269609
38.3944
73
2
19
26227680
62.5612
87
7
20
25565820
58.6688
88
6
21
28711069
49.449
85
4
22
27301580
58.5472
85
7
23
34016415
47.0961
78
7
24
30601546
46.3975
82
3
25
29752815
56.0421
82
7
26
41024410
38.2994
68
6
27
31604741
44.3053
78
3Table 9
Result of eigenvalues.
Eigenvalues
3.1204
0.7428
0.1111
0.0257
Proportion
0.780
0.186
0.028
0.006
Cumulative
0.780
0.966
0.994
1.000Table 10
Result of principle components.
Variables
PC1
PC2
PC3
PC4
Inventory
0.543
0.089
−0.818
0.167
Cost
−0.538
−0.277
−0.511
−0.610
Service level
−0.551
−0.210
−0.230
0.774
Robustness
0.335
−0.933
0.126
0.023As it is visible in Table9, cumulative proportion of the first and second components is 0.966 which is greater than 0.95, so the first two components are considered as principal components. Their relative weights are 0.78 and 0.186, respectively. Final weights are normalized and results are reported in Table 11.Table 11
Weights of criteria for different products.
Criteria
Products
A
B
C
D
Average of inventory position
0.1012
0.0590
0.1942
0.0462
Cost
0.2852
0.1836
0.1742
0.3081
Service level %
0.1195
0.1392
0.0966
0.1206
Robustness %
0.4941
0.6182
0.5350
0.5251


### 5.4. Ranking and Selection of Scenarios

Weights of criteria obtained by PCA method in addition to decision matrix are input of VIKOR method.S, Q, and R values are calculated based on (14) and (15) while strategic weight indicated by v is equal to 0.5, so individual and group utility have same importance. Calculated values are reported in Table 12 for entire scenario as VIKOR output.Table 12
VIKOR output for productB.
Run number
S
R
Q
1
0.634311
0.796269
0.71529
2
0.003184
0.185077
0.09413
3
0.059804
0.042283
0.051044
4
0.536946
0.592539
0.564742
5
0.718862
0.796269
0.757565
6
0.116278
0.089887
0.103082
7
0.428814
0.388808
0.408811
8
0.86078
0.796269
0.828525
9
0.02043
0.140653
0.080541
10
0.677386
1
0.838693
11
0.266254
0.592539
0.429396
12
0.043962
0.185077
0.11452
13
0.571695
0.796269
0.683982
14
0.515097
0.592539
0.553818
15
0.09173
0.185077
0.138404
16
0.648693
0.796269
0.722481
17
1
1
1
18
0
0
0
19
0.630698
1
0.815349
20
0.4239
0.796269
0.610084
21
0.139527
0.388808
0.264167
22
0.6581
1
0.82905
23
0.819338
1
0.909669
24
0.022245
0.185077
0.103661
25
0.726183
1
0.863091
26
0.849516
0.796269
0.822893
27
0.072768
0.185077
0.128922Product
B. According to Table 12, values of R, S, and Q are zero for scenario 18. Regarding (C1) (Acceptable Advantage), the difference between the best scenario and second ranked scenario should be more than 0.0038. The second ranked scenario regarding Q value is scenario 3 and its value is 0.051044 that satisfies both (C1) and (C2). Final result for product B is unique scenario that is number 18 and consists of the following configuration. Inventory policy is DF with R
p being equal to 30 products and 7 days lead time. According to decision matrix which is shown in Table 8, the performance of this scenario for product B results in 38.3944 products as average of inventory position and expected cost is equal to 35269609 for product life cycle. 73% of customer demand will be satisfied and, in the case of increasing demand intensity and variation, 2% decrease is expected in service level, so it will change into 71%. Analyses for other products are similar, so only final results are presented in this paper.Product
A. In this case, VIKOR method ranks scenario number 15 as the best one. This scenario consists of DF inventory policy, 20 products for R
p and 5 days for lead time. Performance of this scenario results in 37.4778 products as average of inventory position. Expected cost is 31603294 and service level will be 66%. This scenario is robust against demand intensification and variation because changes in demand affect service level as small as 1%. But scenario number 15 does not satisfy (C1), so compromising solutions are considered. Although the best scenario is number 15, VIKOR method implies that compromising solutions have the same values. So, in real situation, compromising solution can be substituted. This solution regarding their priority is as follows.Scenario number 24: inventory policy is DF,R
p is equal to 30 products, and lead time is 5 days.
Scenario number 18: inventory policy is DF,R
p is equal to 20 products, and lead time is 7 days.
Scenario number 12: inventory policy is DF,R
pis equal to 20 products, and lead time is 3 days.Product
C. For this product, scenario number 1 is ranked as the best one. This scenario consists of FI inventory policy, R
p is 15 products, and lead time is 3 days. Result of this scenario is 43.87 products as the average of inventory position. Expected cost is 66143000 and service level will be 61% with robustness of 1%. Compromising solutions are as follows.Scenario number 6: inventory policy is DF,R
p is 15 products, and lead time is 5 days.
Scenario number 3: inventory policy is DF,R
p is 15 products, and lead time is 3 days.
Scenario number 7: inventory policy is FI,R
p is 15 products, and lead time is 7 days.
Scenario number 4: inventory policy is FI,R
p is 15 products, and lead time is 5 days.Product
D. Scenario number 3 has the best rank for this product. This scenario consists of DF inventory policy, R
p is 15 products, and lead time is 3 days. This scenario leads to 36.3896 products as average of inventory position with 62648337 expected cost and 61% service level with 1% robustness. Compromising solutions are as follows.Scenario number 6: this scenario is mentioned before as a compromising solution for product C.
Scenario number 17: inventory policy is FI,R
p is equal to 15 products, and lead time is 7 days.


### 5.5. Sensitivity Analysis

As mentioned in the previous section, all scenarios related to productA prefer DF inventory policy, so it is concluded that DF policy is unique optimal policy while lead time includes values of 3, 5, and 7 in selected scenarios. This situation implies that product A is insensitive to the delivery time. It can be wise to decide lead time about 5 days, but there is no need for strict control on lead time while R
p should be controlled strictly to prevent inventory in hand from reaching lower than 20 products.In the case of productB, (C1) is satisfied so there is a unique scenario while there are four compromising solution for product C. In this situation, based on VIKOR method, compromising solutions have same value while interaction plot can be applied for investigation of existing difference between compromising solutions. As mentioned before, the best scenario suggests FI inventory policy for this product but in compromising solutions there are both FI and DF policies. Lead time varies from 3 to 7 days. So, this confusing situation is resolved by deeper analysis and applying interaction plot. Figure 4 shows interaction plot for DOE factors and average of inventory position as response variables.Interaction plot for productC: average of inventory position.
(a)
(b)
(c)Figure4 shows that with either DF or FI policy, there is no significant difference among different lead times in the view point of inventory position, but DF leads to less average of inventory position. With R
p equaling 15, there is less sensitivity to different lead times. So, configuring R
p which is equal to 15 products, the average of inventory position criterion will be robust against lead time variations.Figure5 regards cost criterion and shows that FI and DF policies are equal and make least sensitivity related to lead time. On the other hand, three plots of R
p against lead time are parallel. The parallel plots are interpreted as lack of relation between R
p and lead time. With R
p equaling 15, increase in inventory cost could be mitigated when lead time sets to 3 days.Interaction plot for productC: cost.
(a)
(b)
(c)In Figure6, with 3 days for lead time, there is no significant difference between inventory policies in the view point of service level. R
p and lead time are independent and the best service level is related to FI policy, which is also very sensitive to R
p. As R
p increases, service level will grow. According to the interaction plot analyses, optimal decision for product C is FI inventory policy, 3 days lead time, and 15 or more products for R
p.Interaction plot for productC: service level %.
(a)
(b)
(c)Decision about productD is simple. Based on the information obtained from solutions, it is recommended to employ DF inventory policy with R
p being equal to 15 products and 3 days for lead time.


## 5.1. Multiresolution Result for Demand Modelling

Result of demand modelling is presented for productB because this analysis for other products is similar. For this problem, polynomial functions are applied for yearly and monthly demand rates. Equations (18) and (19) are estimated by nonlinear regression model. In these functions t
y and t
m are time variables
(18)
f
(
t
y
)
=
0.008942
×
t
y
3
-
0.03884
×
t
y
2
+
0.1867
×
t
y
-
0.08018
and 95% confidence interval is presented in Table 4 for each coefficient of (18).Table 4
95% confidence interval for estimated coefficient off
(
t
y
).
‍Coefficients
Lower bund
Upper bund
0.008942
0.004835
0.01305
−0.03884
−0.07604
−0.001652
0.1867
0.08654
0.287
−0.08018
−0.1568
−0.003523In this table, coefficients of (18) are visible in the first column while lower bund and upper bund of each coefficient are presented in the second and third columns, respectively. There is no interval including zero, so it is concluded that coefficient estimation is valid.Estimation error forf
(
t
y
) is reported based on the sum of square errors (SSE) which is 1.504
×
10
-
6. Consider
(19)
f
(
t
m
)
=
4.768
×
10
-
5
×
t
m
5
-
0.001692
×
t
m
4
+
0.0214
×
t
m
3
-
0.1159
×
t
m
2
+
0.3474
×
t
m
-
0.2286
.
Monthly demand rate function that is indicated by f
(
t
m
) and its coefficients confidence interval are presented in Table 5. SEE is equal to 0.0005825. Based on local business information, total demand of product B is estimated to be about 6822 units in five years. According to the information about product B, Figure 2 shows final multiresolution result for this product.Table 5
95% confidence interval for estimated coefficient off
(
t
m
).
‍Coefficients
Lower bund
Upper bund
4.768 × 10−5
4.451 × 10−5
5.086 × 10−5
-
0.001692
−0.001795
−0.00159
0.0214
0.0202
0.0226
-
0.1159
−0.1221
−0.1097
0.3474
0.334
0.3608
-
0.2286
−0.2371
−0.22Figure 2
Comparison between real and estimated demand.


## 5.2. Simulation Parameters and Configuration

For simulation of problem, it is necessary to configure inventory policy parameters and variables. We have three different policies, namely, FQ, FI, and DF. Variables of FQ policy are reorder point (R
p), lead time (L), and economic quantity orders (Q
j
t). First two variables are configured based on Table 2 and third variable is derived from Wilson formula, according to average of demand, holding cost (C
h), and ordering cost (C
r). The calculated values are 57, 71, 72, and 72 for product A, B, C, and D, respectively.Orders quantity is unknown for FI policy and is calculated based on probability distribution function which is obtained by demand simulation of each product for 20 days’ time interval. Distribution fitting is performed by input analyser of Arena software and graphical result for productB is shown in Figure 3. Results for other products are reported in Table 6.Table 6
Estimation of distribution functions of demand and reorders quantity.
Product type
Distribution function
P value
Order quantity
A
278* BETA(1.13,3.28)
>0.15
100
B
GAMMA(67.8,1.18)
>0.15
92
C
EXPO(81.8)
>0.15
100
D
446* BETA(0.957,3.21)
>0.15
115Figure 3
Distribution function of customer demand for reorder interval of productB.In this table, fitted probability distribution function for each product is presented in second column whileP value of fitting which is greater than 0.15 for fitted distribution grantees goodness of fitting in 95% confidence level. Finally, according to fitted distribution functions, for FI policy, proper ordering amount of each product is presented in the fourth column.DF policy has two unknown parameters, namely,α and β. The first parameter is weight of current real demand in (4) and second one is the weight of current trend in (5). These parameters are adjusted by simulation of five years demand. Demand for the entire four products simultaneously is considered and an average of lost sale percentage is considered as response variable. Different values of response variables are reported in Table 7. According to the different values of α and β, the best values for α and β are 0.9 and 0.1, respectively, which results in minimum lost percentage in average.Table 7
Average of lost sale percentage for different values ofα and β.
α and β values
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
Response toα
0.380
0.356
0.334
0.320
0.306
0.292
0.284
0.276
0.268
Response toβ
0.296
0.302
0.306
0.306
0.306
0.304
0.305
0.304
0.302


## 5.3. Weight Assignment

As mentioned before, PCA is applied for weight assignment. For this purpose, DOE result is needed in the form of decision matrix that is shown in Table8. In this table, four criteria: cost, average of inventory position, service level, and robustness are evaluated for each scenario of product B. Decision matrix should be scaleless that is performed by linear method. Then, the result of eigenvalues is calculated and shown in Table 9. Also result of principal components is presented in Table 10.Table 8
Decision matrix for productB.
Scenario number
Cost
Average of inventory position
Service level %
Robustness %
1
32984000
43.52
76
6
2
29214710
39.7514
79
3
3
37114636
33.9623
69
2
4
34840720
41.0948
71
5
5
36820570
34.4738
72
6
6
39191837
33.1091
67
2
7
37037180
39.1012
68
4
8
40994680
29.8712
64
6
9
41407045
31.9678
63
1
10
28282930
52.4462
83
7
11
26422145
48.9243
85
5
12
31023829
41.6049
79
3
13
30395190
48.7146
79
6
14
35393680
39.5812
74
5
15
32838055
39.5646
77
3
16
33285615
46.3341
76
6
17
39652020
33.5668
66
7
18
35269609
38.3944
73
2
19
26227680
62.5612
87
7
20
25565820
58.6688
88
6
21
28711069
49.449
85
4
22
27301580
58.5472
85
7
23
34016415
47.0961
78
7
24
30601546
46.3975
82
3
25
29752815
56.0421
82
7
26
41024410
38.2994
68
6
27
31604741
44.3053
78
3Table 9
Result of eigenvalues.
Eigenvalues
3.1204
0.7428
0.1111
0.0257
Proportion
0.780
0.186
0.028
0.006
Cumulative
0.780
0.966
0.994
1.000Table 10
Result of principle components.
Variables
PC1
PC2
PC3
PC4
Inventory
0.543
0.089
−0.818
0.167
Cost
−0.538
−0.277
−0.511
−0.610
Service level
−0.551
−0.210
−0.230
0.774
Robustness
0.335
−0.933
0.126
0.023As it is visible in Table9, cumulative proportion of the first and second components is 0.966 which is greater than 0.95, so the first two components are considered as principal components. Their relative weights are 0.78 and 0.186, respectively. Final weights are normalized and results are reported in Table 11.Table 11
Weights of criteria for different products.
Criteria
Products
A
B
C
D
Average of inventory position
0.1012
0.0590
0.1942
0.0462
Cost
0.2852
0.1836
0.1742
0.3081
Service level %
0.1195
0.1392
0.0966
0.1206
Robustness %
0.4941
0.6182
0.5350
0.5251


## 5.4. Ranking and Selection of Scenarios

Weights of criteria obtained by PCA method in addition to decision matrix are input of VIKOR method.S, Q, and R values are calculated based on (14) and (15) while strategic weight indicated by v is equal to 0.5, so individual and group utility have same importance. Calculated values are reported in Table 12 for entire scenario as VIKOR output.Table 12
VIKOR output for productB.
Run number
S
R
Q
1
0.634311
0.796269
0.71529
2
0.003184
0.185077
0.09413
3
0.059804
0.042283
0.051044
4
0.536946
0.592539
0.564742
5
0.718862
0.796269
0.757565
6
0.116278
0.089887
0.103082
7
0.428814
0.388808
0.408811
8
0.86078
0.796269
0.828525
9
0.02043
0.140653
0.080541
10
0.677386
1
0.838693
11
0.266254
0.592539
0.429396
12
0.043962
0.185077
0.11452
13
0.571695
0.796269
0.683982
14
0.515097
0.592539
0.553818
15
0.09173
0.185077
0.138404
16
0.648693
0.796269
0.722481
17
1
1
1
18
0
0
0
19
0.630698
1
0.815349
20
0.4239
0.796269
0.610084
21
0.139527
0.388808
0.264167
22
0.6581
1
0.82905
23
0.819338
1
0.909669
24
0.022245
0.185077
0.103661
25
0.726183
1
0.863091
26
0.849516
0.796269
0.822893
27
0.072768
0.185077
0.128922Product
B. According to Table 12, values of R, S, and Q are zero for scenario 18. Regarding (C1) (Acceptable Advantage), the difference between the best scenario and second ranked scenario should be more than 0.0038. The second ranked scenario regarding Q value is scenario 3 and its value is 0.051044 that satisfies both (C1) and (C2). Final result for product B is unique scenario that is number 18 and consists of the following configuration. Inventory policy is DF with R
p being equal to 30 products and 7 days lead time. According to decision matrix which is shown in Table 8, the performance of this scenario for product B results in 38.3944 products as average of inventory position and expected cost is equal to 35269609 for product life cycle. 73% of customer demand will be satisfied and, in the case of increasing demand intensity and variation, 2% decrease is expected in service level, so it will change into 71%. Analyses for other products are similar, so only final results are presented in this paper.Product
A. In this case, VIKOR method ranks scenario number 15 as the best one. This scenario consists of DF inventory policy, 20 products for R
p and 5 days for lead time. Performance of this scenario results in 37.4778 products as average of inventory position. Expected cost is 31603294 and service level will be 66%. This scenario is robust against demand intensification and variation because changes in demand affect service level as small as 1%. But scenario number 15 does not satisfy (C1), so compromising solutions are considered. Although the best scenario is number 15, VIKOR method implies that compromising solutions have the same values. So, in real situation, compromising solution can be substituted. This solution regarding their priority is as follows.Scenario number 24: inventory policy is DF,R
p is equal to 30 products, and lead time is 5 days.
Scenario number 18: inventory policy is DF,R
p is equal to 20 products, and lead time is 7 days.
Scenario number 12: inventory policy is DF,R
pis equal to 20 products, and lead time is 3 days.Product
C. For this product, scenario number 1 is ranked as the best one. This scenario consists of FI inventory policy, R
p is 15 products, and lead time is 3 days. Result of this scenario is 43.87 products as the average of inventory position. Expected cost is 66143000 and service level will be 61% with robustness of 1%. Compromising solutions are as follows.Scenario number 6: inventory policy is DF,R
p is 15 products, and lead time is 5 days.
Scenario number 3: inventory policy is DF,R
p is 15 products, and lead time is 3 days.
Scenario number 7: inventory policy is FI,R
p is 15 products, and lead time is 7 days.
Scenario number 4: inventory policy is FI,R
p is 15 products, and lead time is 5 days.Product
D. Scenario number 3 has the best rank for this product. This scenario consists of DF inventory policy, R
p is 15 products, and lead time is 3 days. This scenario leads to 36.3896 products as average of inventory position with 62648337 expected cost and 61% service level with 1% robustness. Compromising solutions are as follows.Scenario number 6: this scenario is mentioned before as a compromising solution for product C.
Scenario number 17: inventory policy is FI,R
p is equal to 15 products, and lead time is 7 days.


## 5.5. Sensitivity Analysis

As mentioned in the previous section, all scenarios related to productA prefer DF inventory policy, so it is concluded that DF policy is unique optimal policy while lead time includes values of 3, 5, and 7 in selected scenarios. This situation implies that product A is insensitive to the delivery time. It can be wise to decide lead time about 5 days, but there is no need for strict control on lead time while R
p should be controlled strictly to prevent inventory in hand from reaching lower than 20 products.In the case of productB, (C1) is satisfied so there is a unique scenario while there are four compromising solution for product C. In this situation, based on VIKOR method, compromising solutions have same value while interaction plot can be applied for investigation of existing difference between compromising solutions. As mentioned before, the best scenario suggests FI inventory policy for this product but in compromising solutions there are both FI and DF policies. Lead time varies from 3 to 7 days. So, this confusing situation is resolved by deeper analysis and applying interaction plot. Figure 4 shows interaction plot for DOE factors and average of inventory position as response variables.Interaction plot for productC: average of inventory position.
(a)
(b)
(c)Figure4 shows that with either DF or FI policy, there is no significant difference among different lead times in the view point of inventory position, but DF leads to less average of inventory position. With R
p equaling 15, there is less sensitivity to different lead times. So, configuring R
p which is equal to 15 products, the average of inventory position criterion will be robust against lead time variations.Figure5 regards cost criterion and shows that FI and DF policies are equal and make least sensitivity related to lead time. On the other hand, three plots of R
p against lead time are parallel. The parallel plots are interpreted as lack of relation between R
p and lead time. With R
p equaling 15, increase in inventory cost could be mitigated when lead time sets to 3 days.Interaction plot for productC: cost.
(a)
(b)
(c)In Figure6, with 3 days for lead time, there is no significant difference between inventory policies in the view point of service level. R
p and lead time are independent and the best service level is related to FI policy, which is also very sensitive to R
p. As R
p increases, service level will grow. According to the interaction plot analyses, optimal decision for product C is FI inventory policy, 3 days lead time, and 15 or more products for R
p.Interaction plot for productC: service level %.
(a)
(b)
(c)Decision about productD is simple. Based on the information obtained from solutions, it is recommended to employ DF inventory policy with R
p being equal to 15 products and 3 days for lead time.


## 6. Conclusion

In this paper, a simulation optimization framework is proposed for robust optimization of retailer inventory system. This framework is less time consuming in comparison with metaheuristic or SPSA approach, it also provides robust solution for multiple objective functions. In this research also trend and cyclic change of customers demand are considered for more realistic view of problem. Proposed framework consists of discrete event simulation and surrogate model. In this framework, full factorial design of experiment is employed for producing decision space. Based on the three decision factors, 27 different scenarios are produced and simulated. In simulation modelling, multiresolution method is applied for simulation of demand with highly dynamic pattern. Because there are multiple criteria for inventory system performance, VIKOR method is employed as multicriteria decision making technique. In addition, robustness is considered as performance criterion and PCA is used for improving performance of VIKOR method. Finally, developed framework gave the best ranked and compromising solutions, so interaction plot applied for more investigation and sensitivity analysis of obtained solutions.Due to the novelty of the proposed framework, it is recommended that future studies encompass optimization of inventory system of integrated supply chain. Also improvement of developed framework can be considered as future study.



---

*Source: 305120-2014-09-02.xml*