A simulation optimization approach to apply value at risk analysis on the inventory routing problem with backlogged demand

Article history: Received January 26 2014 Received in Revised Format June 2


Introduction
To enhance the supply chain performance, coordination and integration of various components of the supply chain have been attracting attention.In fact, an effective method to improve supply chain performance is necessary for the vendor for determining the quantity which is ordered by its downstream customers.This approach, which is known as vendor managed inventory (VMI), can reduce inventories and stock-outs by using advanced online messages and data retrieval systems (Aviv & Federgruen, 1998).In VMI the supplier monitors the inventory of each customer and decides the replenishment policy of each customer.Vendor managed re-supply is an example of value creating logistics and creates value for both suppliers and customers, i.e., a win-win situation (Campbell & Savelsbergh, 2004).Hence, the supplier is responsible for the inventory level of each customer and acts as a central decision-maker, therefore, the supplier has to solve an integrated inventory-routing or production-inventory-routing problem.
VMI is seen as an effective mean for managing inventory through the strategic and effective use of technologies, which enable the flow of information throughout the entire supply chain.Notwithstanding its potential benefits, only a relatively small number of articles have analytically approached the issue of integrating decisions (Cáceres-Cruz et al., 2012).This phenomenon is probably due to the complexity of the process.In the body of literature this issue is known as Inventory Routing Problem (IRP) (Campbell et al., 2002).
The combination of vehicle routing and inventory management is the common definition for IRP.The formulation of this problem often takes the form of a complex mixed integer model, hence it is a NP-Hard problem (Chen & Lin, 2009).In IRP a supplier or distributor has to deliver products to a number of geographically dispersed customers, subject to some constraints (Coelho et al., 2012).To tackle this problem, some studies have been proposed where some of them proposed deterministic models (Bard & Nananukul, 2009;Bertazzi et al., 2002) and some others have dealt with stochastic models (Schedl & Strauss, 2011;Bertazzi et al., 2011;Chen & Lin, 2009).Deterministic approaches commonly assume that the knowledge supporting the problem are perfectly accessible.Since this assumption is not always true, the deterministic approaches are not always appropriate.In other hand, routing problem stochastic inventory (SIRP) can cope with the intrinsic characteristic of demands swhich is stochastic.
Our aim in this paper is to propose a new simulation optimization procedure to solve IRP problems.This study considers one or a combination of some risk factors to design the inventory routing problem model.Some variations of risks in inventory management and routing problem are as follows:  Disruption risks where the supply of products may be disrupted based on some disasters. Stock-out risks, which means the supplier or the customers ran out of stock in some periods. Cost risk, which implies the maximum risk we can take on cost expenditure.
The simulation optimization procedure, handles a group of scenarios.
The remainder of the paper is organized as follows: In Section 2, we introduce and study the inventory routing problem.In Section 3, we take a closer look at simulation optimization method that we used to tackle the problem.Finally, in Section 4, a numerical example is presented.

Literature Review
This section presents a review on some of the researches in the field of IRP.In this section we devide IRP into five major parts: (1) exact algorothms that try to solve IRP problems, (2) heuristic algorithm, (3) stochastic inventory routing problem (SIRP), (4) stochastic and dynamic inventory routing problem and (5) risk in inventory literature.Archetti et al. (2007) considered a distribution problem in which a product has to be shipped from a supplier to several retailers over a given time horizon.They presented a mixed-integer linear programming model and derived new additional valid inequalities used to strengthen the linear relaxation of the model.They implement a branch-and-cut algorithm to solve the model into optimality.Another branch-and-cut algorithm was presented by Coelho and Laporte (2012) in which they extended the model of Archetti et al. (2007) and formulated the problem as multi-product, multi-vehicle.

Heuristic algorithm
Most of the early papers on the IRP have applied simple heuristics to simplify various versions of the problem.Current heuristic approaches are more powerful and usually use Meta-Heruristic methods, which apply a local search procedure and policies to avoid local optimum.Over the last years, interest in hybrid metaheuristics has risen considerably among researchers in combinatorial optimization.The best results found for many practical or academic optimization problems are obtained by hybrid algorithms (Talbi, 2002).Archetti et al. (2012) presented a heuristic, which combines a Tabu Search scheme with ad hoc designed mixed integer programming models.They considered an inventory routing problem in discrete time slots where a supplier has to serve a set of customers over a time horizon.The objective is the minimization of the summation of the inventory and transportation costs.Popović et al. (2012) developed a Variable Neighborhood Search (VNS) heuristic for solving a multi-product multi-period IRP in fuel delivery with multi-compartment homogeneous vehicles, and deterministic consumption that varies with each petrol station and each fuel type.

Stochastic Inventory Routing Problem (SIRP)
In the SIRP, the supplier considers customer demand only in a probabilistic sense.Demand stochasticity means that shortages may occur.In order to discourage them, a penalty is imposed whenever a customer runs out of stock, and this penalty is usually modeled as a proportion of the unsatisfied demand.Unsatisfied demand is typically considered to be lost, that is, there is no backlogging.The objective of the SIRP remains the same as in the deterministic case, but is written so as to accommodate the stochastic and unknown future parameters: the supplier must determine a distribution policy that minimizes its expected discounted value (revenue minus costs) over the planning horizon, which can be finite or infinite (Coelho et al., 2012).

Heuristic algorithm
Geiger and Sevaux (2011) studied a problem with unknown demand varying within 10% of a mean value.They proposed several policies based on delivery frequencies for each customer.They provided the Pareto front approximation of such policies when moving from a total routing-optimized solution to an inventory-optimized one.In order to solve the problem for several periods, they apply the record-torecord travel heuristic of (Li et al., 2007).et al. (2012) introduced a robust inventory routing problem where a supplier distributes a single product to multiple customers facing dynamic uncertain demands over a finite discrete time horizon and backlogging of the demand was allowed.Their problem was to determine the delivery quantities as well as the times and routes to the customers while ensuring feasibility regardless of the realized demands and minimizing the total cost composed of transportation, inventory holding, and shortage costs.

Stochastic and dynamic Inventory Routing Problem
In the dynamic IRP, customer demand is gradually revealed over time (e.g., at the end of each period) and one must solve the problem repeatedly with the available information.In Dynamic and Stochastic Inventory-Routing Problems, the demand of the customers is known in a probabilistic sense and would be revealed over time, thus yielding a dynamic and stochastic problem.Tirado et al. (2012) introduced a dynamic and stochastic maritime transportation problem originating in industrial shipping.Three heuristic methods adapted to this problem were considered and their performance in minimizing transportation costs was evaluated.In the first heuristic, called the myopic dynamic heuristic (MDH), a deterministic sub problem was solved using an adaptation of Tabu Search presented by Korsvik et al. (2010).The second heuristic was based on multiple scenario approach with consensus (MSAC) of Bent and Van Hentenryck (2004).It relied on generating a set of scenarios that consist of the currently known vehicle along with a sampled realization of future cargoes.Each of these scenarios was solved using Tabu Search.The third heuristic that they proposed was an adaptation of the branch-and-regret heuristic (BRH) of Hvattum et al. (2007).Again, scenarios were generated and solved using Tabu Search.They tested the three heuristics on a number of realistic, but randomly generated test instances.

2.5
Risk in Inventory Literature Wu et al. (2009) studied risk-averse newsvendor model with mean-variance objective function.They showed that stock out had a significant impact on newsvendor ordering decisions and the risk-averse newsvendor did not necessarily order less than the risk-neutral newsvendor did.Yang et al. (2007) extended the classical newsvendor problem by presenting a downside risk constraint from the viewpoint of inventory control.They analyzed the form of the optimal order quantity when they restricted that the probability that the cost level was larger than or equal to a fixed cost constant was less than a fixed value of probability.Luciano et al. (2003) built a decision model, where the choice concerns the quantity to be ordered to face a random demand, in order to optimize the expected result (an expected cost to be minimized or an expected profit to be maximized) the model explored the probability distribution of the result, both via analytical methods and via simulation methods.Chen et al. (2007) proposed a framework for incorporating risk aversion in multi period inventory models.They considered two related problems, one with exogenous demand and the other one with demand depending on price.In both cases, they differentiated between models with fixed-ordering cost and without fixed-ordering cost.Zhang et al. (2009) presented some convex stochastic programming models for single and multi-period inventory control problems where the market demand is random and order quantities need to be decided before demand is realized.

Problem Description and Model Formulation
Consider an inventory routing problem that contains one distributor and one type of product.The distributor 0 delivers products to a set of N retailers over a given time horizon T. The product could be sent from the distributor to the retailers at each period t T  by a vehicle or a group of vehicles that have the capacity of Q.In each period a quantity it q of products is sent to the retailer i that have the demand it  .In some cases, we might face shortage of product it it q   .
Our contribution to the literature is considering a risk measure in the model to control the cost that depot can afford.The cost consists of holding costs, carrying costs and backlogging cost.We consider a depot center or distributor who has financial limits.In other words, due to the capital limitations, we cannot satisfy every demand sent to our center.The risk measure changes the problem in three possible ways:  The route structure can change.
 The amount of product sent to retailers can be determined.
 Controllable factors such as number of vehicles to satisfy the constraint related to risk measure must be tunes.
The above-mentioned assumptions will be clarified after running the model.The assumptions are summarized as follows:  The time horizon is finite and the period unit is day.
 The distributor has the capacity C 0 and the retailers have the capacity i C is assumed to be discrete.
 A fleet of homogeneous vehicles K is considered and each with capacity Q.  Transportation cost of each vehicle ( ij c ) on each period t is constant.
 The initial inventory for each distributor is considered. The model structure is one to many. The routing is multiple. Unmet demand is backlogged as back order. The model assumes holding cost for both distributors and retailers. Inventory policy is the maximum level (ML). Vehicles routes are cyclic, meaning that each the vehicle that travels its route returns back to the distributor's node.
 The basic IRP is defined on a graph = ( , ) where V is the vertex set of nodes ( ) is the set of arcs.Vertex indexed zero presents the supplier and

Mathematical Model
subject to ( ) , 1 1 Eq. ( 1) gives the optimum total cost including the expected inventory cost for all customers, the variable transportation cost, the expected holding cost, the expected shortage cost for customers and the fixed vehicle cost, respectively.Constraint (2) is a risk measure applied on cost.Constraints (3) are the inventory balance constraints for individual customers.Constraints (4) ensure that the quantity sent to each customer in each period should be less than the depot capacity.Constraints (5) ensure that every customer's warehouse inventory capacity should be no less than its maximum inventory level in the period 1. Constraints (6) ensure that every customer's warehouse inventory capacity should be no less than its maximum inventory level in the period (t).Constraints (7) ensure that loading of each vehicle in each period does not exceed respective capacities.Constraints (8) ensure that the number of vehicles used for delivery in each period does not exceed the number of vehicles.Constraints (9) ensure that the number of vehicles leaving from a customer or the depot is equal to that of arriving vehicles.Constraints (10) ensure that each vehicle can travel maximum once in each period.Constraints (11) ensure that each vehicle from current node can only travel to only another node.Constraint ( 12) and ( 13) ensure that at end of each period, each retailer cannot have both the inventory and back-order.Constraint ( 14) is logical relationships between ( ikt q ), ( jikt X ).Constraints (15) avoid sub-tours for each vehicle at each period.Finally, constraints ( 16) to (20) show the type of variables.

Simulation Optimization
The solution approach framework of this research is shown in Fig. 2. First, according to this figure, the probability distribution of the customers must be identified, and based on this, the customers' demands is generated.The next step is to determine each scenario and number of runs that simulation must be implemented.Based on the scenarios, total cost for each run is calculated and then the probability distribution must be fitted to the total costs.In the next step, the VaR measure is identified based on significance level.For validating the solution procedure, some runs are randomly selected and supposing that their demands are deterministic, the exact solutions are obtained using the GAMS software.In this research, if the gap between the proposed and exact solution is less than 10% then the next scenario is simulated.Based on the VaR measure of all scenarios, scenarios are ranked and the best is chosen as the best solution of the problem.Since the Stochastic IRP is an NP-Hard problem, and with a VaR constraint becomes more complex, it is not possible to solve the large scale (and usually real-world cases) with an exact method.Simulation is a technique that measures and describes various characteristics of the bottom-line performance measure of a model when one or more values for the independent variables are uncertain."Stochastic simulation optimization" (often shortened as simulation optimization) refers to stochastic optimization using simulation Specifically, the underlying problem is stochastic and the goal is to find the values of controllable parameters (decision variables) to optimize some performance measures of interest, which are evaluated via stochastic simulation, such as discrete-event simulation or Monte Carlo simulation (Fu, 2002;Fu et al., 2008).The act of obtaining the best decision under given circumstances is known as optimization.

Fig. 3. Simulation Framework
Simulation allows one to specify a system accurately by using the logically complex, and often nonalgebraic, variables and constraints.Complex and stochastic systems can therefore be modeled through simulation.While the development of the new technology has dramatically increased computational power, efficiency is still a big concern when using simulation for stochastic optimization problems.
There are two important issues in the efficiency: (1) at the level of stochastic simulator, a large number of simulation replications must be performed in order to capture the randomness and obtain a sound statistical estimate at a specified level of confidence.
(2) at the level of optimization, many design alternatives must be well evaluated via stochastic simulation in order to determine the best design or to iteratively converge to the optimal design.
The objective in simulation is to describe the distribution and characteristics of the possible values of the bottom-line performance measure (Y ) and the independent variables ( 1 2 , ,..., k X X X ).The idea behind simulation is similar to the notion of playing out many what-if scenarios.The difference is that the process of assigning values to the cells in the spreadsheet that represent random variables is automated so that: (1) the values are assigned in a non-biased way, and (2) the spreadsheet software user is relieved of the burden of determining these values.With simulation, repeatedly and randomly, sample values are generated for each uncertain input ( 1 2 , ,..., k X X X ) in the model, and then the result value of our bottom-line performance measure   Y is computed.Then, the sample values of (Y ) can be used to estimate the true distribution and other characteristics of the performance measure ( Y ).For example, sample observations are used to construct a frequency distribution of the performance measure, to estimate the range of values over which the performance measure might vary, to estimate its mean and variance, and to estimate the probability of the actual value of the performance measure greater (or less) than a particular value.All these measures provide greater insight into the risk associated with a given decision than a single value calculated based on the expected values for the uncertain independent variables.In particular, we now have some idea of the best-case and worst-case total cost outcomes for the company.We have a better idea of the distribution and variability of the possible outcomes, and a more precise idea about where the mean of the distribution is located.It is also known a way of determining how likely it is for the actual outcome to fall above or below some value.

Distribution fitting
The critical value ( , N D  ) is pre-determined.After the data are collected, data can be plotted using a histogram, and from the histogram or the problematic nature of the data, it can be determined which distribution is the most appropriate.Then, the parameters for that distribution can be estimated from the data using some principles such as Maximum Likelihood Estimation (MLE) or Minimum Squared Errors (MSE).

Goodness of fit test
After identifying the distribution, it is necessary to check how representative the fitted distribution is.It can be done by using either a heuristic approach or formal statistical tests.In the heuristic approach, a probability plot can be used, and from the plot, it can visually be determined if the distribution really fits well with the data.As for statistical tests, either chi-square tests or Kolmogorov-Smirnov tests can be used to check whether the distribution fits the data well.The idea of the Chi-Square Test is to compare the histogram with the probability mass function of the fitted distribution.The Kolmogorov-Smirnov Test focuses on the comparison of the distribution functions between the data and the fitted one.

Chi-square tests
This is an old but popular test.A chi-square test is considered as a formal comparison of a histogram with the probability density or mass function of the fitted distribution.The steps of doing the chi-square test are (Rao & Scott, 1981): 1. Draw the frequency histogram of the data.2. Compute the probability content of each interval, i p based on the fitted distribution.

Compute the expected frequency in each interval ˆi
where n is the total number of observations.4. Let i f be that of observations in the interval.Compute the test statistic: 5. Determine the critical value at the level of significant  and 2 ,    with d degrees of freedom, where d = (number of intervals) − (number of estimated parameters) 6. Conclude that the data does not fit the distribution if the test statistics is greater than the threshold.

Kolmogorov-Smirnov tests
Instead of grouping data in each interval as chi-square tests, Kolmogorov-Smirnov tests intend to compare the empirical distribution function with the fitted distribution ( F ). Kolmogorov Smirnov tests do not require users to group data in any way and so no information is lost.The steps of doing the Kolmogorov-Smirnov test are (Massey Jr, 1951): 1. Rank the data from the smallest to the largest, i.e.,      

Random Number and Variables Generation
When running a simulation, it is necessary to generate random numbers that follow some certain distributions for capturing the system's stochastic behaviors, such as service times or customer interarrival times.In computer simulation, a procedure to generate a sequence of numbers is used in order to behave similarly to the random number, which are called pseudo random numbers.The process that generates these numbers is called random number generator.
In order to generate numbers fitting to a certain distribution, it is needed to first generate a number following a uniform distribution between zero and one, and then some random variable generation methods are used to transform the uniform random number into a random number that follows the required distribution.

4.4.Output Analysis
Output analysis aims at analyzing the data generated by simulation, and estimates the performance of the system.Since the simulation is stochastic, multiple simulation replications must be performed in order to have a good estimate.The required number of simulation replications depends on the uncertainties of the simulation output.In output analysis, it is useful to look into ways to analyze the simulation output and then decide how long each simulation should be run and how many times is needed to replicate the simulation.

4.5.Verification and Validation
After the simulation model is constructed, it must be verified and validated to ensure it represents the actual system.Verification refers to check the model to see whether the conceptual model is accurately coded into a computer program while validation involves with checking the model to see whether it sufficiently represents the actual system so that the objectives of the simulation study can be achieved.The purpose of validation is to compare the behavior of the simulation model to the real system.However, for validation, it should be noticed that although no model is completely representative of the system under study, some models could be useful.Hence, it is always important to know the objective of the simulation study so that the simulation model is able to help to achieve the proposed objective.There are two different types of tests, which can be used for validation.
Subjective test: It is also known as a qualitative test.People who are knowledgeable on one or more aspects of the system should be involved, and let them to check if the model is reasonable.
Objective test: Is also known as a quantitative test.Statistical tests are used to compare some aspects of the system data set with the model data set.The model assumptions must be validated, and then compare the model input-output relationships to corresponding input-output relationships in the real system.Test 2 is utilized in this study.

Experimental results
In this section, some computational examples are presented to illustrate possible applications of the proposed model for considering the risk of lack of capital in an IRP.In this section some computational examples are presented to illustrate possible applications of the proposed model for considering the risk of lack of capital in an IRP.It is assumed that all customer demands are following the uniform distribution with different bounds as shown in Table 2.These distributions remain steady as the system goes forward.The Distributions of the demand at each customer is as in Table 2.We run each scenario for over 100 times.Each run contains seven working days ( 7 T  ).Considering this, all surplus and shortage demand until the seventh day are retained and then for next run the history of the system is erased (i.e., surplus and shortage demand in 0 t  (beginning of the period) is equal to zero).Performing these runs, the cost of each run is calculated.It should be noticed that in a scenario cost of transportation and vehicle fixed cost don't change, because the transportation cost is dependent on the route structure and qi and since neither of them change it remains fixed in every run.The cost matrix for transportation is as Table 3.The fixed cost of vehicles is considered 1000 unit.So the routing cost that contains a fixed cost of vehicles and the cost of transportation is calculated.The Table 4 Total Transportation Cost for each Scenario (first 15) shows these costs for all scenarios.For each period (day), the back order cost or holding cost is also calculated and added to transportation cost.The back order and holding cost per unit is as Table 6.By calculating all of these costs, a number is dedicated to every run.So for each run total cost is considered.These total costs follow a distribution that should be identified.Thorough a chi-square test explained in previous section, the best distribution of total costs is fitted.Different distributions can be obtained from the chi-square test.The distributions for the example scenarios are shown in Table 7.Note that the colored cells in P-Value column are the ones that the distribution related to them is chosen through a Kolmogorov Smirnov test.In the state, where chi square test is not reliable enough, hence, for choosing the best distribution Kolmogorov Smirnov test is utilized.After the distribution of each scenario is identified.VaR measure is calculated considering four different significance levels.
The significance levels considered are 0.005, 0.025, 0.05, and 0.10.The significance levels may change by user preference.
The scenarios are designed by best possible routes and different set of i q .Each i q is calculated based on the holding cost and back order cost.The logic behind determining each i q is that if the holding cost of a customer is bigger than its backorder cost ( i i h p  ) a certain ratio of standard deviation is subtracted to the mean of that customer demand.Otherwise if the holding cost is lower than the back order cost ( i i h p  ) a certain ratio of standard deviation is added to the customer demand.This makes sure that the system bear lower cost.Table 8 and Table 9 show the results of 30 scenarios.For each significance level the rank of each scenario is calculated and regard to manager preference for significance level the best scenario is chosen.Fig. 4 shows the box plot for all scenarios that can be used to recognize best scenarios especially when the VaR is not important.As it is clear from the, in 0.9 significance level, the scenario number 30 is the best scenario.For 0.95 and 0.975 significance level, again scenario 30 with 5-2.3-1.4 route shape shown in figure 5 is the best among all scenarios.However, for 0.995 significance level, the scenario number 13 comes to top.According to Fig. 5, the orange lines are return routs that are the routs each vehicle travels for coming back to the depot.To validate the simulation random runs are chosen from each scenario and considering the demand related to that run is deterministic it is solved through an exact model represented in problem definition section.If the result of two answers is not deviated more than 10 percent then the scenario is accepted.The table below shows some example of validation for the problem we solved.
To validate the solution some runs have been picked up from each scenario.The runs are supposed to be deterministic and are solved by an exact model and as we discussed in the model if the gap between two answers is more than 10% the solution is rejected.Table 10 shows the results for scenario 10.The exact solution is done through Cplex solver.The platform in which the software is run is an Intel 2.5 GHz Core i5 CPU with 4 GB RAM.As it is clear in the table below that the time for solving the problem by an exact method even in deterministic form is timely infeasible.

Conclusion and Future Research
Inventory Routing Problem addresses the issue of coordinating inventory replenishment policies and distribution plans in a cost effective manner.More precisely, the IRP integrates inventory and distribution aspects in the same planning process.Considering the risk of allocating capital to this matter has been quite an issue to the managers.In this study a model has been proposed that deals with the risk of allotting too much capital in an IRP with stochastic demand.A VaR measure has been considered to be able to model the risk.This measure acts like a constraint in the model and does not let the cost of an inventory routing exceeds from VaR in a known significance level.To solve the model, a simulation optimization approach has been applied.The approach helps dealing with uncertainty of demand.To validate the simulation optimization some samples randomly are chosen and have been solved thorough exact method.The results of validation in experimental results showed the simulation optimization approach is quite adequate.
For further works, one may consider the problem by heterogeneous vehicles.Also, modeling a dynamic stochastic inventory routing problem (DSIRP) while considering the risk can be an appropriate field of study.For more efficient simulation optimization, one can make a pre-selection for scenarios.

Fig. 1
Fig. 1 shows a schematic structure example of IRP.

Fig. 2 .
Fig. 2. The proposed simulation optimization structure Fig. 4. Total Cost Box Plots for each Scenario

Fig. 5 .
Fig. 5. Route Structure for Optimized Scenarios (Scenario 30) Popović et al. (2012)the expected losses subordinate to risk aversion constraints expressed through Value at Risk (VaR) and Conditional Value at Risk (CVaR) as risk measures.Popović et al. (2012)considered the single period stochastic inventory (newsvendor) problem with downside risk constraints.They used Value at Risk as a risk measure in newsvendor framework and extended it to multi product newsvendor.
Cumulative probability distribution function of demand

Table 5
Total Transportation Cost for each Scenario (second 15)

Table 6
Holding cost and Backorder Cost for Customer i

Table 7
Distribution for scenarios total costs

Table 8
Results of Simulation Optimization for First 15 Scenarios

Table 9
Results of Simulation Optimization for Second 15 Scenarios