On the comparison of inventory replenishment policies with time-varying stochastic demand for the paper industry

Accepted Manuscript On the comparison of inventory replenishment policies with time-varying stochastic demand for the paper industry David Escu´ın, Lorena Polo, David Cipr´es PII: DOI: Reference:

S0377-0427(16)30154-6 http://dx.doi.org/10.1016/j.cam.2016.03.027 CAM 10581

To appear in:

Journal of Computational and Applied Mathematics

Received date: 25 November 2015 Revised date: 17 March 2016 Please cite this article as: D. Escu´ın, L. Polo, D. Cipr´es, On the comparison of inventory replenishment policies with time-varying stochastic demand for the paper industry, Journal of Computational and Applied Mathematics (2016), http://dx.doi.org/10.1016/j.cam.2016.03.027 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

*Manuscript Click here to view linked References

On the comparison of inventory replenishment policies with time-varying stochastic demand for the paper industry David Escu´ına,∗, Lorena Poloa , David Cipr´esa a Instituto

Tecnol´ogico de Arag´on, Mar´ıa Luna 7-8, 50018 Zaragoza, Spain

Abstract The aim of this paper is the development of a mathematical model to compute the optimal inventory mix to face stochastic demand at minimum cost in a two-level supply chain. The paper addresses a multi-product dynamic lot-sizing problem under stochastic demand subject to capacity and service level constraints. This model is executed to compare a Make To Order (MTO) strategy to a Vendor Managed Inventory (VMI) partnership between the supplier and their customers. Both strategies provide the demand order to be produced. A schedule of production orders is determined over the planning horizon in order to minimize the inventory holding costs of the supply chain, taking into consideration that the supplier is also responsible of initiating the replenishment orders and deliveries of their customers according to the VMI partnership. The simulation model is illustrated empirically using a real case study: a paper manufacturing company that pursues to improve customer service level and supply chain inventory costs through a proper production planning of their paper machines and a suitable VMI order replenishment schedule. Keywords: lot-sizing problem, Vendor Managed Inventory, simulation, optimization, paper industry

1. Introduction Production lot-sizing has a great impact on inventory, particularly under seasonal fluctuations of demand and constrained production capacity. Many companies adopt the MTO (Make To Stock) policy in which products are not built until a confirmed order for products is received by the manufacturer. Other companies maintain high levels of inventory (stock) to face periods of uncertain demand. However, a production schedule which does not adjust accurately the real demand may lead to overstocks for some products and stock-outs for other. Inventory sizing by product is especially important under uncertainty, when the inventory is necessary to guarantee a service level in a stochastic environment. One of the integration practices that can contribute to reduce inventory in the supply chain is Vendor Managed Inventory (VMI). VMI programs allow for consumer demand information to be disseminated up the supply chain, thus mitigating upstream demand fluctuations due to the bullwhip effect [6] and [8] . Due to this demand anticipation, VMI may allow to reduce logistics and manufacturing costs, reduce overall lead-times, improve service level and reduce transportation costs. The aim of this paper is the development of a mathematical model to seek the most effective inventory mix to face stochastic demand at minimum cost in a two-level supply chain. We focus on a multi-product dynamic lot-sizing problem under stochastic demand subject to capacity and service level constraints. Unlike previous studies, this model is executed to compare a MTO strategy to a VMI partnership between the supplier and their customers [4]. Both policies are developed within the model, and their results are compared within the numerical application. The work presents some forecasting nonlinear optimization models that can be brought from Applied Science as done in [7], [9] and [10]. In the problem, a schedule of production orders is determined over the planning horizon in order to minimize the inventory holding costs of the supply chain, taking into consideration that the supplier is also responsible of initiating the replenishment orders and deliveries of their customers according to the VMI partnership. The model also considers features such as service level required, the production capacity at machine level, set up time or product-machine allocation. The integration of stochastic demand in the production/inventory model is performed through the statistical distribution of the forecast accuracy. Historical data are analyzed to select the most suitable forecasting model for each reference (also called SKU). The selected model is triggered to forecast the demand during the rolling horizon. The applicability of the proposed model is illustrated empirically using a real case study: a paper manufacturing company that pursues to improve customer service level and supply chain inventory costs through a proper production ∗ Corresponding

author Email addresses: [email protected] (David Escu´ın ), [email protected] (Lorena Polo ), [email protected] (David Cipr´es)

Preprint submitted to Elsevier

March 17, 2016

planning of their paper machines and a suitable VMI order replenishment schedule. A cost analysis of the supply chain inventory under different service levels and under different adoption rates of VMI, in concordance with local convergence achieved (see [1]]) by the simulation and optimization, demonstrates the potential of this model to improve performance in the supply chain. The paper is organized as follows. In Section 2 we give an overview of previous work in replenishment policies and lot sizing problems. In Section 3, we describe the methodology applied to the paper industry by placing emphasis on the forecasting models which help to estimate the demand of the customers. The mathematical formulation is described in Section 4. Here we will see how the customers and manufactures variables are connected with the replenishment policies. The simulation algorithm is described in Section 5 as a loop where customers and manufacturers trigger orders to be manufactured when the control parameters reach certain values. In Section 6 the results obtained are analyzed and the discussion of the outcomes and conclusions are shown in Section 7. 2. Literature review The work in [18] was one of the first papers that provided insights into the VMI approach by explaining the savings that so often accrue from this strategy. In addition, they describe some underlying technologies required to make the arrangement work. In [6], the authors seek to find the supply chain that minimizes system cost through comparing performance between traditional and VMI systems. A mathematical model is developed, and total supply chain cost is used as the measure of comparison. The models are applied in both traditional and VMI supply chains based on pharmaceutical industry data, and we focus on total cost difference compared through the use of Adjusted Silver Meal (ASM) and Least Unit Cost heuristics. The work in [17] aims at designing a dynamic VMI system in which the entire supply chain performance is optimized in terms of production planning at vendors site, distribution strategy, and inventory management at manufacturers site. The VMI system is modeled as a mixed-integer linear program (MILP) using discrete-time representation and the mathematical representation follows the resource-task net-work formulation. To address the complexity of the problem, the problem is solved directly using an exact detailed model, an iterative procedure combines a novel aggregate model with the detailed model and a novel rolling horizon approach that simultaneously combines the aggregate and the detailed models is designed to solve the problem. In [8] is presented a fractal-based approach for inventory management. A fractal-based echelon does not indicate a functional level or composition of supply chain members but indicates a group of multi- or hetero-functional fractals. The basic fractal unit (BFU) consists of five functional modules including an observer, an analyzer, a resolver, an organizer, and a reporter. They develop mathematical models for the analyzer and resolver to effectively manage supply chain inventories. The paper [15] compares three replenishment strategies in a two-echelon serial supply chain. These strategies include Make-to-Order, Make-to-Stock, and Vendor Managed Inventory. To compare these strategies, they develop several probability models that we use to calculate the expected fill rate, average customer inventory, and average manufacturer inventory. The second study determines the benefit of VMI in a two-echelon arborescent supply chain. [3] compares the expected performance of a VMI supply chain with a traditional serially-linked supply chain. The emphasis of this investigation is the impact these two alternative structures have on the Bullwhip Effect generated in the supply chain. They pay particular attention to the manufacturers production ordering activities via a simulation model based on difference equations. 3. Methodology for the the paper industry problem The case study addressed in this paper is related with the supply chain for the pulp and paper industry. The sector comprises companies that use wood as raw material and produce pulp, paper, board and other cellulose-based products. This work is focused in the two-echelon operations from the manufacturers (facilities that processes the raw material) and the customers that receives the paper and produces strong, lightweight cardboard boxes. The manufacturers must meet the customer’s requirements in terms of paper quality, wide and other paper properties. Several approaches are launched to schedule the production into the paper machines. Board must be dispatched to the customers by different transport modes as in [5]. The study compares the performance of VMI and MTO strategies in the two-echelon serial supply chain with finite production capacity. The performance is investigated using difference equations forming a simulation model. The purpose of using simulation is to provide a simplified environment into which a number of situations and ideas can be tested. From the simulation outputs, the effect of the strategies on manufacturing on-costs, inventory holding costs and transport costs will be quantified and discussed. 2

As said before, with an MTO strategy, the manufacturer produces in response to a customer order, whereas the manufacturer produces in anticipation of orders with an VMI strategy (see [3] and [18]). The paper industry considered in the use case triggers orders in a calendar agreed with the manufacturer (weekly or every two weeks). To calculate the size of the order, the customers estimate the demand during the lead time and dispatched as [13] [14]. It cause large orders to cover the demand under the uncertainty. This is the key aspect that VMI seeks to lower. 3.1. Overall architecture. Firstly, at initial stages, it is necessary to analytically calculate the customers safety stock (how much extra stock that is maintained to mitigate risk of stockouts due to uncertainties in supply and demand), the manufacturer inventory, and fill rate for VMI and MTO strategies. They are based on the lead time, the delay between the time the reorder point is reached and renewed availability. This will be the starting point prior to adjust a iterative approach aimed at performing forecasting and optimization tasks iteratively following similar optimization as in [16]. The overall vision of the paper is reflected in the Figure 1 and is also generalized by [15]. Manufacturers and customers provide qualitative and quantitative information to adopt the forecasting and decision models. Experts work in a collaborative environment and demand forecasts and error estimates (an estimate of how far actual demand may be from forecasted demand) feed the decision model which optimizes the lot sizing (quantity of paper to be produced) of the manufacturers. Historical data stored at the manufactures and warehouses is used as an input to update models and fit forecasts.

information

Forecasting Model

Manufacturers & Consumers

information

parameter update

demand forecast error estimate

Decision Model

Updating Model

sales data historical data

Manufacturer lot sizing

Figure 1: Iterative approach general flow chart

3.2. Demand forecasting The simulation of the supply chain policies requires to previously forecast the demand for each day of the rolling horizon. Capturing the time-varying stochastic pattern within the replenishment policy is a key factor in order to ensure cost minimization. This is done by analyzing univariate and multivariate forecasting methods in a preliminary phase. The methods evaluation analyzes the accuracy of a wide set of methods in order to select the most suitable subset of methods for each data series. The multivariate methods proposed are shown in Table 1. These methods need at least a covariable to perform, contrary to univariate methods, that are only based on history orders. Table 1: Multivariate forecasting methods

Acronym LREG DREG ANN

Name Linear Regression Dynamic Regression Artificial Neural Network

3

Multivariate linear regression (LREG) is a mathematical method that models through a linear equation the relationship between an output variable Yi , dependent variables Xi and a random term  that is supposed to be normally distributed. The general form of the method is as follows: Yt = a + b1 X1t + b2 X2t + · · · + bk Xkt + t ,

t = 1, · · · , T

The multivariate dynamic regression method (DREG) is an advanced linear regression method as it supports w autoregressive terms of the output variable, z lag terms of the independent variables and v Cochrane-Orcutt autoregressive error terms. Yt = a + b1 Yt−1 + · · · + bw Yt−w + c1 Xt−1 + · · · + cz Xt−z + d1 t−1 + · · · + dv t−v

The Cochrane-Orcutt procedure allows to reduce the bias due to an autorregresive structure of the regression error, by transforming the forecasting model with the application of a cuasi-difference. Model fitting in DREG has been carried out by the software [2]. yt − ρyt−1 = α(1 − ρ) + β(Xt − ρXt−1 ) + et

A neural network (ANN) is a multivariate classification method from artificial intelligence that learn complex non linear interrelations between independent and dependent variables. A well known neural network is a Multilayer Perceptron (MLP), that consists of three types of neural layers: the input layer, the hidden layer and the output layer. The idea is that the network filters dependent variables through one or more hidden layers before getting the output. The software WEKA [11] has been used to train the neural network thanks to its powerful in machine learning and classification methods for data mining. In Table 2 we show the set of univariate methods used. The Simple Moving Average method (SMA) is one of the easiest and common methods in time series forecasting, that’s why usually it is considered a naive method. The form of SMA is given by the formulation:

Yt−1 + Yt−2 + · · · + Yt−(n−1) n Whereas in SMA the past observations are weighted equally, the exponential smoothing method family provides exponentially decreasing weights over time. The most common form, the simple exponential smoothing (SES), is recommended when there is not trend nor seasonality and is represented by the following formulation: Yt =

Yt = αYt−1 + (1 − α)Yt−2

Other exponential smoothing methods included in the analysis are the Holt method (HOLT), recommended for series with additive trend and without seasonality, and Holt-Winters (HWIN) applied to series with additive trend and seasonality. An Exponential smoothing state space method (ETS) is also included in the tests. It is a state space model that includes some transition equations that describe how the unobserved components or states (level, trend, seasonal) change over time. The classical decomposition method splits a time series into a trend and a seasonal component and projects them in the forecast horizon. The method used here deals with multiplicative (DES-M) components. It assumes that seasonal components keep constant. The ARIMA (p,d,q) method (ARIMA) includes a combination of autoregressive terms (AR), integration (I) and moving average (MA) in one only model. Since it can represent a correlation structure with the minimum number of parameters, is usually denominated the optimal univariate model. ARIMA can be represented by the following formulation: Yt = −(∆d Yt − Yt ) + φ0 +

p X i=1

φi ∆d Yt−i −

q X

θi εt−i + εt

i=1

Three methods used specifically for intermittent demand (CROSTON, SBA and SBJ) are also included in the tests. Model fitting of all these univariate methods has been carried out through R package [12]. The evaluation is performed by fitting forecasting models on the in-sample data set and testing them on the out-ofsample through some accuracy measures. Accuracy evaluation is usually carried out through the comparison of one or several accuracy indicators. The selection of an accuracy indicator to use is a crucial issue, as different accuracy measures may lead to different conclusions. The most usual error measures are the mean absolute percentage error (MAPE) and the RMSE. However, MAPE gets undefined when demand is zero and provides bad performance in outliers or actual data near to zero. RMSE is scale dependent and is not recommended when comparing the error among series, but it is effective in the identification of the best model among a group of models for single series, that is our case. So, in order to compare the performance of the 4

different forecasting methods, the square root of the mean squared error is selected as accuracy measure. A final remark must be done about the measurement of forecast accuracy. Previous research has shown (rather conclusively) that the best method in terms of forecast accuracy is not necessarily the best when evaluating the impact on resource scheduling. In this case, the 3PL resource planning is based on a direct translation of orders into operators by means of a multiplication factor. So, the implications of forecasting accuracy metrics are directly applicable to resource planning impact kpis. That’s why we still consider RMSE a suitable indicator for this analysis. Table 2: Univariate forecasting methods

Acronym SMA SES DES-M ETS HOLT HWIN ARIMA CROSTON SBA SBJ

Name Simple Moving Average Simple Exponential Smoothing Multiplicative seasonal trend decomposition State Space method Holt method Holt-Winters ARIMA Box-Jenkins Croston method Syntetos-Boylan Approximation Shale-Boylan-Johnston approach

4. Mathematical formulation Once the demand values have been predicted, this section reports on the description of the optimization problem. Let K be the set of customers and M the number of manufacturers. Let S be the set of SKUs and D be the number of days of the simulation. 4.1. Variables for the customers. The Table 3 lists the set of the control variables and parameters for the customers under VMI or MTO replenishment policies. The selected demand forecasting model for each SKU s and customer k returns not only the forecasted demand, dk,s,d , but also the forecasted error as seen in the previous section. mk,s denotes the error ∀k ∈ K, ∀s ∈ S. Once all the demand forecasts are calculated, they are summed up to get the demand average through tdk,s,d and dak,s during the rolling horizon D. The forecast deviation f dk,s refers to the error assumed due to demand fluctuations. This is expressed as the forecasted error multiplied with the demand average and multiplied with the square square root of the sum between the lead time and revision period. The key control parameters for the simulation are the safety stock, ssk,s , and the target stock, tsk,s,d . These variables control the performance of the main loop and respond to the questions of when and how much to order. The safety stock here is calculated based on the inverse of standard deviation of service level multiplied with the forecast deviation. To calculate the target stock, it is necessary previously to get the demand forecast during the lead time. The sum of forecasted demand during the lead time plus the revision period give the demand forecast regardless of the replenishment policy. As denoted above, the lead time and revision period differs from VMI to MTO. As VMI permits to daily check the inventory levels, rpk,s = 1, ∀k ∈ K, ∀s ∈ S whereas rpk,s = 7, ∀k ∈ K, ∀s ∈ S for the MTO policy as the customers are required to make the orders weekly (for the simulations, the Tuesdays). This means that C should be setup with all Tuesdays during the rolling horizon. With regards to lead time, MTO policy makes orders to arrive later than VMI because of manufacturer has to produce the paper and ship the board. The lead time considered to VMI is ltk,s = 2 days as for the MTO case ltk,s = 15 days ∀k ∈ K, ∀s ∈ S. 4.2. Variables for the manufacturers. Similarly, as done with the customers, manufacturers parameters for VMI are calculated as reflected in the Table 4. Note that the manufacturers must control two kinds of inventory levels: stock accumulated for the MTO orders that are awaiting to be shipped (smtom,s,d ) and VMI stock originated by the internal replenishment, svmim,s,d . Binary assignment parameters, amkm,k , represent the network structure of the supply chain. The demand forecast dm,s,d is the sum between the demand forecast for the assigned customers. Likewise, the safety stock and target stock are calculated with the aim of launching orders internally. 5

Table 3: Notations for the customers under VMI or MTO.

stock lead time service level revision period demand forecast forecast error total demand forecast demand forecast average forecast deviation safety stock demand forecast (during lead time) target stock orders calendar

sk,s,d ∀k ∈ K, ∀s ∈ S, ∀d ∈ D ltk,s ∀k ∈ K, ∀s ∈ S slk,s ∀k ∈ K, ∀s ∈ S rpk,s ∀k ∈ K, ∀s ∈ S dk,s,d ∀k ∈ K, ∀s ∈ S, ∀d ∈ D mk,s ∀k P∈ K, ∀s ∈ S tdk,s,d = d∈D dk,s,d ∀k ∈ K, ∀s ∈ S dak,s = tdk,s,d /D ∀k p ∈ K, ∀s ∈ S f dk,s = mk,s · dak,s · ltk,s + rpk,s ∀k ∈ K, ∀s ∈ S ssk,s = φ−1 (slk,s ) · f dk,s ∀k ∈ K, ∀s ∈ S Pd+ltk,s +rpk,s dfk,s,d = i=d+1 dk,s,i ∀k ∈ K, ∀s ∈ S, ∀d ∈ D tsk,s,d = dfk,s,d + ssk,s ∀k ∈ K, ∀s ∈ S, ∀d ∈ D C = the set of days, d ∈ D, in which customers make orders under MTO. Table 4: Notations for the manufacturers.

stock VMI stock MTO assignment manufacturer - customer lead time service level revision period demand average demand forecast forecast deviation safety stock demand forecast (during lead time) target stock

svmim,s,d ∀m ∈ M, ∀s ∈ S, ∀d ∈ D smtom,s,d ∀m ∈ M, ∀s ∈ S, ∀d ∈ D amkm,k ∈ 0, 1 ∀m ∈ M, ∀k ∈ K ltm,s,d ∀m ∈ M, ∀s ∈ S, ∀d ∈ D slm,s ∀m ∈ M, ∀s ∈ S rpm,s P ∀m ∈ M, ∀s ∈ S dam,s = Pk∈K dak,s · amkm,k ∀m ∈ M, ∀s ∈ S dm,s,d = k∈K dfk,s,d · amk p m,k ∀m ∈ M, ∀s ∈ S, ∀d ∈ D f dm,s,0 = mm,s · dam,s · ltm,s + rpm,s ∀m ∈ M, ∀s ∈ S ssm,s,d = φ−1 (slm,s ) · f dm,s,d ∀m ∈ M, ∀s ∈ S, ∀d ∈ D Pd+ltm,s +rpm,s dfm,s,d = i=d+1 dm,s,i ∀m ∈ M, ∀s ∈ S, ∀d ∈ D tsm,s,d = dfm,s,d + ssm,s,d ∀m ∈ M, ∀s ∈ S, ∀d ∈ D

The value of the lead time for the manufacturers depends directly on the machine production cycle. The machines are configured to produce a ordered sequence of each sku. A cycle starts with the manufacturing of a sku and ends when all skus have been manufactured. However, the length of each cycle is variable. As this sequence must be met, the lead time is calculated as the sum of the number of days between two consecutive production date plus 1 (a delay of one day from the order is placed until it is processed). Note that the lead time, safety stock and target stock take values ∀d ∈ D. 4.3. Objective function 2 1 1 , m ∈ M, βs,k ,k ∈ K as the cost of manufacturing a unit of item s in the manufacturer m. Let βs,m Let us define αs,m be the holding cost of unit of item s for the manufacturers and customers respectively. The transportation costs from the 1 manufacturers to the customers are represented by ζm,k . The objective function (1) minimizes the total cost that aggregates production, inventory and transportation costs is:

min(α + β + ζ)

(1)

where α=

XX X

1 αs,m fm,s,d

(2)

d∈D s∈S m∈M

β=

X X X 1 2 ( βs,m (svmim,s,d + smtom,s,d ) + βs,k sk,s,d )

s∈S m∈M

(3)

k∈K

ζ=

X X

1 ζm,k em,k,s,d

(4)

m∈M k∈K

being em,k,s,d the quantity of product s sent from the manufacturers to the customer k and fm,s,d the quantity of product manufactured, ∀m ∈ M, ∀k ∈ K, ∀s ∈ S, ∀d ∈ D. These variables will used later on the section 5.1. 6

5. Simulation algorithm The simulation algorithm (Algorithm 1) executes a loop of D days in which customers and manufacturers must update their internal variables based on the previous periods. Decisions made in each customer individually affect in the decisions taken in the assigned manufacturers in terms of inventory levels and quantity production. Though initial values in d = 0 for the stock are be set properly, the model requires a warm up phase. Some cycles (days) are required to smooth the production and balance the flows of product. In practice, this means that some peaks in the stockouts may be caused or the machine capacities may be reached. It is worthwhile to note here that prior to start running the model, a preliminary stage calculates the error, the safety stock and the target stock for all customers, SKUs and manufacturers. How the forecasting models defined above are initiated (section 3.2), filtered and adjusted is out of the scope of this paper and would certainly deserve other work apart. It is mentioned for the sake of completeness. dk,s,d = demandF orecasting(); mk,s = computeError(); ssk,s = getSaf etyStock(); tsk,s,d = getT argetStock(); for d ∈ D do updateCustomers(); updateM anuf acturers(); end Algorithm 1: Main simulation block for the paper industry.. 5.1. Inventory replenishment cycle for the customers. The customers inventory replenishment is denoted in Algorithm 2. Firstly, depending on the type of replenishment policy, the starting stock, d = 0, is initialized. Under MTO, it is assumed that the lead time multiplied with the demand average may give a adequate value; on the other hand, VMI policy suggests the demand average could be sufficient. During the loop, the stock level, sk,s,d , is updated according to the products received and the demand consumption in the previous day. The virtual stock, vsk,s,d , reflects how much quantity of product is available by considering the pending orders to be received. When the virtual stock is less than target stock, the model produces an order equal to the difference between both variables. Under VMI, it happens regardless the day but MTO policy, as explained earlier, this must be done the days agreed with the manufacturers (weekly for the experiments). sk,s,0 = ltk,s · dak,s = starting stock for MTO; sk,s,0 = dak,s = starting stock for VMI; ok,s,d = quantity of product ordered; rk,s,d = quantity of product received; pk,s,d = quantity of product that is pending to receive; vsk,s,d = virtual stock; C = orders calendar; for d ∈ D do sk,s,d = sk,s,d−1 − dk,s,d−1 + rk,s,d−1 ; pk,s,d = pk,s,d−1 − rk,s,d−1 + ok,s,d−1 ; vsk,s,d = sk,s,d + pk,s,d ; if vsk,s,d < tsk,s,d then if (MTO and d ∈ C) or VMI then ok,s,d = tsk,s,d − vsk,s,d ; end end end Algorithm 2: Inventory replenishment for the customers under MTO or VMI.. 5.2. Lot-sizing approach. From the manufactures side, they have to determine the lot sizes and the sequence of lots while satisfying the demand requirements. The machines must produce the lots for VMI and MTO policies and internal replenishment of the manufacturer. The capacity machines must not be exceed as well as not all SKUs can be manufactured in all machines. Due to the complexity of the lot-sizing and scheduling problem addressed by the manufacturers, the approach proposed assume 7

that setup times are sequence-independent due to the paper manufacturer system and all the new orders are given the same priority (regardless of MTO or VMI policy). Having said this, the manufacturer put into stack all the orders to be manufactured plus the pending orders. The pending orders have priority and they are manufactured in order of size (from the largest to the lowest). It may happen that a order could not be produced completely and must be delayed by adding it into the pending queue. When an order can be manufactured in two or more machines, the machine with the lowest volume of orders is selected. This guides the approach to get balanced load of work which is always well seen by the managers (note that maintenance costs and other derived expenses are not considered either). 5.3. Inventory replenishment cycle for the manufacturers. The Algorithm 3 shows the execution loop for the manufacturers once the parameters controls have been calculated as explained above. The fm,s,d reflects the lots manufactured by the running of the lot-sizing approach described above. The lead time is calculated as ldm,s,d = 1 + nfm,s,d to address the cycle of the machines. Though there has not been any reference to the shipments, this is an important issue that affect to the stock levels on all nodes. The paper company has deals with different transport operators to ship the freight under different circumstances. The explanation of how the transport is managed falls outside of the paper but needs to be cited for the understanding of the problem. The shipments for every day are represented by em,k,s,d as cited previously in 4.3. It is worthwhile to cite that not all days the transport is available: road transport is much more flexible than rail transport. svmim,s,0 = dam,s = starting stock for VMI; smtom,s,0 = 0 = starting stock for MTO; em,k,s,d = quantity of product sent to customer k; fm,s,d = quantity of product manufactured; om,s,d = quantity of product that is ordered internally; pm,s,d = quantity of product that is pending to produce; vsm,s,d = virtual stock; nfm,s,d = number of days left for the machine to produce the SKU s; for d ∈ D do sm,s,d = sm,s,d−1 − em,s,d−1 + fm,s,d−1 ; fm,s,d = lotSizingApproach(); pm,s,d = pm,s,d−1 − fm,s,d−1 + om,s,d−1 ; vsm,s,d = vsm,s,d−1 + pm,s,d−1 ; Pd+1+nf P dfm,s,d = t=d+1 m,s,d k∈K dfk,s,t · amkm,k ; p f dm,s,d = mm,s · dam,s · ltm,s + nfm,s,d ; ssm,s,d = φ−1 (slm,s ) · f dm,s,d ; tsm,s,d = dfm,s,d + ssm,s,d ; if vsm,s,d < tsm,s,d then om,s,d = tsm,s,d − vsm,s,d ; end end Algorithm 3: Inventory replenishment for the manufacturers with lot-sizing approach.. 6. Experimental results The simulation model and lot-sizing approach were applied on a network with 2 manufacturers, 10 customers and 30 skus. The tests seek to compare the behavior and costs by varying the percentage of VMI of the SKUs and customers. 1 1 2 1 The objective function parameters, αs,m βs,m βs,k and ζm,k were instantiated with the data provided by the paper industry under studio. The service level is 95% for all customers, manufacturers and skus. The forecasting horizon (the out-ofsample data set) is 2 months long. The preliminary stage of forecasting methods ran 3-month-data. The experiments were run on a server Intel(R) Xeon(R) X5650 2.67Hhz with Ubuntu Linux and 8GB RAM memory. The simulator and lot-sizing algorithm were implemented in Java. An intensive previous work is done to run the forecasting methods and get the demand values. The Table 5 shows the average RMSE kpi for the above selected methods. LREG proves to be the most suitable on average. Nevertheless, for each sku, the selected method may result in other than LREG. When it comes to varying the percentage of VMI, several approaches were considered. One might think that select randomly a set of sku could be the simplest way to do it. This would affect to all the customers with the skus in stock. Other could be to select a set of representative customers and process all their skus as VMI. However, both approaches 8

Table 5: Average results for the forecasting methods.

Method SMA SES DES-M ETS HOLT HWIN ARIMA CROSTON SBA SBJ LREG DREG ANN

RMSE 2.65 1.65 1.63 1.34 1.76 1.49 1.38 2.06 2.21 2.43 1.23 1.31 2.02

could result in misleading results. VMI policy performs better with low demand variability since it is able to balance the lead time and the other control variables. This suggestion is also supported with the idea that the greater lead time, the greater the probability to suffer stock out. To all this should be added that references with high demand should get more promising results as the scheduling of the machines can be better optimized. Then, and through several meetings with the paper company, the following criteria were considered on the selection of skus for VMI policy. They were weighted with a formulation agreed with the company that is not worth describing. • The references with large values of demand. • The references more shared among the customers. • The references with low variability. • The references with larger lead times. 6.1. Results while varying percentage of VMI

Figure 2: Stock levels while varying percentage of VMI.

9

Figure 3: Costs while varying percentage of VMI.

The Figures 2 and 3 illustrate the stock levels and costs incurred while varying the percentage of VMI policy from 0 to 100%. The stock of the customers decreases when the VMI policy gains volume since the stock of the manufacturer goes up. However, this decrease is higher than the increase so the total inventory levels are reduced accordingly. Quite similar results are illustrated in terms of costs. Note that in the optimal scenario, with 100% VMI, the reduction of the stock levels reaches a value of 24% compared to 0% VMI. The reduction in costs is the 16% between both scenarios. This is motivated by the other incurred costs as transport and manufacturing. 6.2. Lead time performance of MTO versus VMI. The Figure 4 illustrates the performance of the lead time for the 10 customers in a scenario 0% VMI (MTO) and 80% VMI. The reduction in the lead time ranges from 50% to 80%. The highest decrease is produced in the customers with high demands as the manufacturer sends the product almost daily or every other day. This helps to maintain balanced levels of inventory and the shipment can be done frequently.

Figure 4: Lead time results for 10 customers under MTO and 80% VMI.

6.3. Stocks comparison between MTO and VMI. The Figures 5 and 6 depict the classical graphs of the inventory levels, demand and receptions controlled by the safety and target stock. Under MTO, there is high stock peaks when the freight arrives and the inventory is much higher than the VMI scheme. The key factors are the target and safety stock values determined by the lead time. Under a MTO policy, a high value of safety stock is needed to cover the stochastic demand. This is certainly true as lead time increases. However, the safety stock for VMI is lower as the lead time is one or two days in this scenario. In other words, as the forecasted errors are accumulated during the lead time, MTO produces large safety stock to reduce the uncertainty by maintaining the service level.

10

Figure 5: Stock performance under MTO in a customer.

Figure 6: Stock performance under VMI in a customer.

7. Discussion of results and conclusions The paper addresses the comparison between two replenishment policies. The results demonstrate the effectiveness of the model and the algorithms under VMI. The safety stock is reduced and, in turn, the inventory levels. By timing production and shipments, the manufacturer is able to manage the inventory level at the customers. VMI has proved to get a better performance than MTO replenishment strategy. The paper simulates the VMI and traditional supply chain MTO strategy in a scenario which previously has been configured with the forecasted demand of the customer. During the rolling horizon, the designed network of supply chain makes replenishment orders to cover the forecasted demand with the objective of minimizing the total costs. It would be interesting to extend the model to address some aspects more deeply. How the shipments are scheduled based on the availability of the transport modes (rail of road transport) and the influence of the machine cycles derive on a great impact on the supply chain. Besides that, explore other models like SVM can be of interest so that the simulation engine may benefit from other accuracy demand values and analyze their correlations with the lot sizing parameters. Acknowledgment The dissemination of this work has been partly financed by the FSE Operative Programme for Aragon (2007-2013). References [1] Argyros, I. K., Gonza´alez, D. Local Convergence for an Improved Jarratt-type Method in Banach Space. International Journal of Artificial Intelligence and Interactive Multimedia. Volume 3, Number 4 ISSN 1989-1660. [2] Baiocchi, G. and Distaso, W. (2003), GRETL: Econometric software for the GNU generation. Journal of Applied Economics, 11

[3] Disney, S.M. and Towill, D.R., 2003 The effect of VMI dynamics on the bullwhip effect in supply chains, International Journal of Production Economics Volume 85, Number 2, 199-215. [4] Fraile, A. Larrod´e, E., Magre˜na´ n, A. Sicilia, J.A. Decision model for siting transport and logistic facilities in urban environments: A methodological approach. Journal of Computational and Applied Mathematics Volumen 290: 100120, 2015. [5] Gallego, J., Larrod´e, E., Sicilia, J.A., Royo, B. y Fraile, A. Optimization and sizing methodology of necessary infraestructure for the incorporation of hydrogen to the transportation sector. DYNA - Ingenier´ıa e Industria., Volume 89: 405–412, 2014. [6] Govindan, K. The optimal replenishment policy for time-varying stochastic demand under vendor managed inventory. European Journal of Operational Research. Volume 242, Issue 2, 402-423, 2015. [7] Ioannis K. Argyros, Magre˜na´ n, A. On the convergence of an optimal fourth-order family of methods and its dynamics, Applied Mathematics and Computation, Volume 252, 1 February 2015, 336-346. [8] Kwangyeol, R., Ilkyeong M., Seungjin O., Mooyoung J. A fractal echelon approach for inventory management in supply chain networks. International Journal of Production Economics Volume 143, Issue 2, 316-326. 2013. [9] Magre˜na´ n, A. Different anomalies in a Jarratt family of iterative root-finding methods. Applied Mathematics and Computation. Volume 233, 1 May 2014, 29-38. [10] Magre˜na´ n, A., Cordero, A. Gutirrez, J.M., Torregrosa, J.R. Real qualitative behavior of a fourth-order family of iterative methods by using the convergence plane. Mathematics and Computers in Simulation. Volume 105, November 2014, 49-61, [11] Mark Hall, Eibe Frank, Geoffrey Holmes, Bernhard Pfahringer, Peter Reutemann, Ian H. Witten, 2009; The WEKA Data Mining Software: An Update; SIGKDD Explorations, Volume 11, Issue 1. [12] R Core Team. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria, 2015 [13] Royo, B., Fraile, A., Larrod´e, E. y Muerza, M.V. Route planning for a mixed delivery system in long distance transportation and comparison with pure delivery systems. Journal of Computational and Applied Mathematics. Volumen 291: 488-496. 2016. [14] Royo, B. Sicilia, J.A., Oliveros, M.J., Larrod´e, E. Solving a Long-Distance Routing Problem using Ant Colony Optimization. Applied Mathematics & Information Sciences Volumen 9. 415-421. 2015. [15] Peter A. Salzarulo. Vendor managed inventory programs and their effect on supply chain performance. Submitted to the faculty of the University Graduate School in partial fulfillment of the requirements for the degree Doctor of Philosophy in the Kelley School of Business. Indiana University 2006. [16] Sicilia, J.A., Royo, B., Quemada, C. Oliveros, M.J, Larrode, E. An decision support system to long haul freight transportation by means of ant colony optimization. DYNA - Ingenier´ıa e Industria. Volume 90: 105-113. 2015. [17] Tareq A. Al-Ameri, Nilay Shah, Lazaros G. Papageorgiou, Optimization of vendor-managed inventory systems in a rolling horizon framework. Computers & Industrial Engineering, Volume 54, Issue 4, May 2008, Pages 1019-1047, ISSN 0360-8352,. [18] Waller, M., Johnson, M., Davis, T. Vendor-managed inventory in the retail supply chain Journal of Business Logistics 1999, 20, pg. 183

12