Stochastic optimization for transshipment problems

Int. J. Production Economics 135 (2012) 61–72

Contents lists available at ScienceDirect

Int. J. Production Economics journal homepage: www.elsevier.com/locate/ijpe

Stochastic optimization for transshipment problems with positive replenishment lead times b, ¨ Yeming (Yale) Gong a,1, Enver Yucesan a b

EMLYON Business School, 23 avenue Guy de Collongue, 69134 Ecully Cedex, France Technology and Operations Management Area, INSEAD, Boulevard de Constance, 77305 Fontainebleau Cedex, France

a r t i c l e in fo

abstract

Article history: Received 5 June 2009 Accepted 16 September 2010 Available online 25 September 2010

Transshipments, monitored movements of material at the same echelon of a supply chain, represent an effective pooling mechanism. Earlier papers dealing with transshipments either do not incorporate replenishment lead times into their analysis, or only provide a heuristic algorithm where optimality cannot be guaranteed beyond settings with two locations. This paper uses infinitesimal perturbation analysis by combining with a stochastic approximation method to examine the multi-location transshipment problem with positive replenishment lead times. It demonstrates the computation of optimal base stock quantities through sample path optimization. From a methodological perspective, this paper deploys a duality-based gradient computation method to improve computational efficiency. From an application perspective, it solves transshipment problems with non-negligible replenishment lead times. A numerical study illustrates the performance of the proposed approach. & 2010 Elsevier B.V. All rights reserved.

Keywords: Supply chain management Stochastic optimization Transshipment Simulation

1. Introduction Physical pooling of inventories (Eppen, 1979) has been widely used in practice to reduce cost and improve service. For example, CIBA Vision has consolidated all of its country-based warehouses in Europe into a single European Logistics Center near Frankfurt to supply contact lenses throughout the continent. On the other hand, the practice of transshipments, the possibility of sharing stock between pairs of locations at the same echelon through enhanced visibility, but without the need to put the stock physically in the same location, has been less popular. A largescale example of transshipping spare parts among the 29 power generating plants of a utility company is provided by Kukreja et al. (2001). To emphasize the requirement for supply chain transparency at the same echelon, we will refer to this practice as information pooling. Transshipments provide an effective mechanism for correcting discrepancies between the locations’ observed demand and their available inventory. As a result, transshipments may lead to cost reductions and improved service (through higher flexibility and responsiveness) in the supply chain without increasing system-wide inventories (see, e.g., Lee et al., 2007; Guinet, 2001). Kukreja et al. (2001) report a 70% reduction in total system cost over an electric utility’s

Corresponding author. Tel.: +33 160724017; fax: + 33 160745500.

E-mail addresses: [email protected] (Y. (Yale). Gong) ¨ [email protected] (E. Yucesan). 1 Tel.: +33 478682796; fax: +33 478337928. 0925-5273/$ - see front matter & 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.ijpe.2010.09.020

decentralized policy. Herer et al. (2002) show that transshipment can be used to enhance supply chain agility and leanness. Allen (1958), Krishnan and Rao (1965), and Lee (1987) are among the earliest works on transshipment problems. Since then the literature on transshipments has been growing rapidly where restrictive assumptions are being relaxed. The reader is referred to Herer et al. (2006) for a comprehensive review of this rich literature. Although they are often overlooked in the literature, replenishment lead times constitute one of the critical factors impacting the effectiveness of a transshipment system. Consider, for example, the Normandy landing where the military logistics system can be viewed as a two-echelon supply chain with the main base as a ‘‘supplier’’ in England and five bases on Normandy beaches in France. On D-Day, Allied Forces landing on Utah Beach met with much less Nazi resistance than those landing on Omaha Beach, which enabled them to move troops and material from Utah Beach to Omaha Beach. This transfer can be viewed as a transshipment, whereby ignoring replenishment lead times, i.e., the time to move new troops and material across the English Channel, would have had disastrous consequences. Similarly, ASML, a Dutch manufacturer of photolithography equipment, reports that its customers in Japan, which manufacture microprocessors, regularly transship spare parts among themselves in order to avoid downtime – hence, lost throughput – due to long replenishment lead times for spare parts shipped from Holland. Finally, the recent financial crisis highlighted the importance of replenishment lead times and transshipments in the banking system. Banks, which need to cover their open positions, have two

62

Y. (Yale) Gong, E. Y¨ ucesan / Int. J. Production Economics 135 (2012) 61–72

sources for additional cash: either they can participate at the auctions organized by the Federal Reserve on the last Thursday of each month or they can borrow cash from each other overnight based on the current prime rate. In all three examples, replenishment lead times play an important role in the operational efficiency and the service level of transshipment systems. Yet, there exist only a handful of papers that explore the transshipment problem with replenishment lead times. Most of these papers are heuristic based. For example, under the assumption that no further transshipments ¨ will take place, Axsater (2003) introduces a heuristic algorithm for a lateral transshipment system with constant positive lead times, and derives a new decision rule especially valuable in more complex decision situations. There are also studies that use enumeration to compute an optimal stock level. For example, under the assumption of Poisson demand and exponentially distributed lead times, Alfredsson and Verrijdt (1999) develop a birth and death model for a three-location transshipment system. Tagaras and Cohen (1992) use grid search over a two-dimensional discrete state space to provide a optimal value. They subsequently use the result of this complete enumeration process as a benchmark for their heuristic method, which has satisfactory performance. Our contribution to this rich literature is the introduction of a practical algorithm that is able to compute optimal parameter values of a policy that minimizes the expected average total cost. More specifically, given a base stock policy, this paper will provide a method to determine the optimal base stock quantities that minimize the expected average total cost in a multi-location transshipment system with positive replenishment lead times. In terms of positive replenishment lead times, this paper extends Herer et al. (2006), who study the multi-location transshipment without replenishment lead times. In terms of a multi-location setting, this paper generalizes Tagaras and Cohen (1992), who consider non-negligible replenishment lead times in two-location inventory systems. To compute the optimal base stock values for a multi-location system with positive replenishment lead times, we consider a sample path optimization (SPO) technique, with the advantage of high efficiency of optimization and wide applicability of simulation. Banerjee et al. (2003) conducted a simulation study on the transshipment problem. SPO has also been applied to optimize a manufacturing system with transportation delay (Mourania et al., 2008). However, SPO requires an estimate of the gradient. There exist various gradient estimation techniques such as infinitesimal perturbation analysis (IPA), likelihood ratios, finite differences, and simultaneous perturbation (Fu, 2002). In this paper, we adopt IPA as an efficient gradient estimation technique (Ho et al., 1979). Successful applications of IPA have been reported in inventory modeling, manufacturing systems analysis (Glasserman, 1994), and supply chain problems (Glasserman and Tayur, 1995). IPA has also been successfully applied in fluid-queue settings. For example, Suri and Fu (1994) is the first paper to apply IPA in the continuous tandem production line. Cassandras et al. (2002) have applied Perturbation Analysis for online control and optimization of stochastic fluid models. Wardi et al. (2002) applies IPA to loss-related and workload-related metrics in a class of stochastic continuous flow models. In this paper, we deploy a linear programming (LP) and network flow model (LP/network), use sample path optimization with IPA, and demonstrate the computation of the optimal base stock quantities. A duality-based gradient computation method is deployed to improve algorithm efficiency and to simplify the proofs of key results. The remainder of the paper is organized as follows: In the following section, we introduce the multi-location transshipment model with positive replenishment lead times, including the

network flow representation, the stochastic programming formulation and its deterministic counterpart. Section 3 is devoted to the details of the algorithm and its implementation. In Section 4, we present partial results from our extensive numerical experimentation, comparing our algorithm to alternative approaches that have already appeared in the literature. We conclude with final comments in Section 5.

2. Model 2.1. Problem description We consider a system with one supplier and N distinct stocking locations (e.g., retailers), which face customer demand. The retailers may differ in their cost and demand parameters. The demand distribution at each retailer in a period is assumed to be known and stationary over time. The system inventory is reviewed periodically and replenishment orders are placed with the supplier. The replenishment order will be delivered after a positive replenishment lead time L, which may be different for different retailers. In each period, the replenishment and transshipment quantities must be determined to minimize the expected average total cost. The total cost is the sum of the replenishment, transshipment, holding, backlog penalty, and expediting costs. Based on the classical result of Veinott (1965), Herer et al. (2006) prove that, in the absence of fixed costs per replenishment order, if transshipments are made to compensate for an actual shortage (instead of building up inventory at another stocking location), there exists an optimal base stock S ¼(S1, S2,y,SN) policy for all possible myopic transshipment policies. In our case, we will continue to adhere to the base stock replenishment policy. In period t, events occur in the following order: First, retailers observe demands. Demand realizations represent the resolution of the only uncertain event in a period. Once demand is observed, decisions about transshipment quantities are made. The transshipment transfers are then made immediately. Subsequently, demand is satisfied. Any unsatisfied demand is backlogged or expedited through an expensive secondary source. At this point, backlogs and inventories are assessed, and penalty and holding costs, respectively, are incurred. Second, replenishment orders that had been placed at the supplier in period t L are received. These orders are used to satisfy the backlog in period t L and, if possible, to increase the inventory level in period t. The decision on the replenishment quantity is then made and a replenishment order is placed with the supplier to be delivered in period t +L. Any remaining inventory is carried to the next period, t + 1. To describe the operation of the system, we use the following notation: L T T N N Dti dti f(D) hi pi

positive and deterministic replenishment lead time planning time horizon f1,2, . . . ,Tg, the time index set number of stocking locations (retailers) f1,2, . . . ,Ng, the location index set random variable associated with demand at retailer i in period t, 8i A N ,t A T actual demand realization at retailer i in an arbitrary period t, 8i A N ,t A T joint probability density function for the demand D holding cost incurred at retailer i per unit held per period, 8i A N penalty cost incurred at retailer i per unit backlogged in the first T L periods, 8iA N


li

ci cîj

cij

expediting cost for the shortfall at retailer i satisfied from a secondary source to avoid lost sales in the last L periods, 8iA N replenishment cost per unit at retailer i, 8i A N direct transshipment cost per unit transshipped from retailer i to retailer j, 8i,j A N , reflecting the logistics cost of moving units from retailer i to retailer j effective transshipment cost per unit transshipped from retailer i to retailer j, cij ¼ cîj þ ci cj , 8i,j A N . Usually, we consider the setting where cij 4 0,8i aj, and cii ¼0

The quantity of interest is Si, base stock quantity at the location i,8i A N . 2.2. Modeling assumptions We will make the following assumptions, which are necessary to avoid pathological cases. Assumption 1 (Lead times). We assume that replenishment lead times are positive and deterministic. Lateral transshipment lead times are negligible between any pair of stocking locations. Although Alfredsson and Verrijdt (1999) consider exponentially distributed lead times, they indicate that, in practice, ‘‘these shipment times are close to deterministic’’. Kukreja et al. (2001) further show that service performance is identical for all lead time distributions with the same finite mean. The negligible lateral transshipment lead time is a regular assumption in transshipment literature since it is usually uneconomical to employ transshipments when lateral transshipment lead time is too large (Kukreja et al., 2001). More specifically, in the presence of long transshipment lead times, the practice may backfire as the items in transit cannot help either the destination (that has a shortage) or the source (that has a surplus), eroding the pooling advantage of transshipments. Assumption 2 (Demand). Customer demand at each retailer is generated by a stochastic process. Demand is backlogged when a retailer is out of stock in t¼1,2,y,T L. In the last L periods, however, since the replenishment orders cannot arrive on time within the finite horizon, any shortfall is expedited through a fast but expensive secondary (local) source to avoid lost sales. Demand has a continuous CDF, but is not necessarily independent across retailers. Here, we assume that, during the last L periods, any shortfall is expedited through a fast but expensive secondary source to avoid lost sales. Such a practice is getting increasingly popular; for example, catalog mail-order firms, which source in large quantities from Southeast Asia before or at the beginning of the selling season, activate secondary sources in Mexico late in the selling season due to long replenishment lead times from Asia. Assumption 3 (Replenishment policy). The base stock quantity is non-negative, which also implies a non-shortage inducing replenishment policy (Herer et al., 2006). A replenishment quantity ordered in period t L arrives in period t and satisfies the backorders observed and generated at period t L; any remaining units go to the next period, t +1. This is a reasonable assumption when backlogged customers are quoted a future delivery date and there is no perceived additional benefit to beat the promised delivery date. The assumption is also valid for customized products where the customer would like to receive exactly what he or she had ordered at the promised time. A non-shortage-inducing replenishment policy eliminates the pathological cases of negative base stocks. No-buildup property is

63

guaranteed by hi rcij þ hj ,8i,j A N , which means that it is not economical to transfer a unit from location i to location j so that it would be held in location j rather than i. Assumption 4 (Transshipment policy). The transshipment policy is myopic, that is, the transshipment quantities are independent of the period in which they are made; they depend only on the pre-transshipment inventory and the observed demand. We assume complete pooling, which means that a location offers its entire available inventory when another location is experiencing a stockout. On the other hand, as stated earlier, we also assume that transshipments are never made to build up inventory at the receiving location, but only made to satisfy current actual shortage. 2.3. Model formulation We present the network flow formulation first, and then introduce the stochastic programming (SP) and its deterministic counterpart. 2.3.1. Network flow representation Given a base stock policy for the replenishment quantities, the optimal transshipment quantities need to be determined each period between every pair of retailers. We develop a linear cost network flow model as follows. In the presence of positive lead times, the inventory position will not always be equal to inventory on hand since there exists inventory on order in the pipeline. Let Hti be inventory on h and at location i at the beginning of period t, and fHt ,t A T g be the stochastic inventory-on-hand process with Ht A R þ . In the presence of positive replenishment lead times, this transshipment system cannot be reduced to a single-period problem. We therefore model the system in a finitehorizon setting. During the planning horizon, we model the movement of stock as a network flow problem. Fig. 1 presents a network with three stocking locations in a setting with a 5-period planning horizon and a 2-period replenishment lead time. At the beginning of each period t, the excess inventory from the previous period is available. This b eginning inventory at retailer i in period t is depicted by a source node, Bti . This stock can be used in one of three different ways: satisfy demand at retailer i, satisfy demand at retailer j (i.e., transshipment from retailer i to retailer j), and hold in inventory at retailer i. The satisfaction of the dem and at retailer i in period t is reflected by the middle sink node denoted by Mti , while the transshipment to retailer j is reflected by another sink node Mtj . Similarly, we will denote by Eti the e nding inventory at retailer i at the end of period t. At the end of the period, the material will be used in two ways: to satisfy backorders or to build up inventory at the retailer. Finally, we have a node Rt to represent the r eplenishment requested from the supplier in period t to be delivered in period t +L. Note that the stock at the beginning of the horizon, and the replenishments made during the first T L periods are the only two sources of material. This setup captures the material flow from the demand side (i.e., the sinks) as well. In period t, the demand at retailer i, can be satisfied in one of two different ways: from the inventory at retailer i, or from the inventory at another retailer j (i.e., through a transshipment from retailer j to retailer i). Another sink for material is the requirement that each retailer i starts the subsequent period with inventory position equal to Si. These units can come from one of two sources: the inventory at retailer i or replenishment arrivals during the period. The arcs in the network flow problem correspond exactly to the activities described above and are summarized (with the associated unit costs) in Table 1. In our LP formulation, we therefore denote the

64


Fig. 1. Network flow representation for a 5-period horizon and 2-period lead time system.

Table 1 Definition of the arcs in the network flow problem. Arc

Variable

Unit cost

(Bti , Eti )

FBt Et

hi

Inventory held at retailer i in period t

(Bti , Mti )

FBt Mt

0

Stock at retailer i used to satisfy demand at retailer i in period t

(Bti , Mtj )

FBt Mt

cij

Transshipment from retailer i to retailer j in period t

(Eti + L, Mti )

FEt þ L Mt

pi

In the first T L periods, shortages backlogged at retailer i in period t

(Eti , Bti + 1)

FEt Bt þ 1

li 0

In the last L periods, shortfall at retailer i in period t Inventory on hand at the end of period t carried to the next period t+ 1

(Rt, Eti + L)

FRt Et þ L

0

Inventory at retailer i increased through replenishment in period t

i

i

i

i

i

j

i

i

i

Meaning

i

i

starting and ending nodes of the f low in the network through such variables as FBt Mt . In other words, FBt Mt is the flow in the network

of the demands in the period. This relationship can be expressed P PN t ¨ as follows: N 2009, i ¼ 1 FRt Et þ L ¼ i ¼ 1 di (see Gong and Yucesan,

from node Bti to Mtj (that is, the quantity transshipped from stocking location i to location j in period t). Unit cost here is the cost per unit flow (that is, the transshipment cost, cij). Given dti , the actual demand realization at retailer i in period t, the system behavior is captured through the auxiliary variable Iti , inventory level at retailer i immediately after the execution of transshipments and the satisfaction of demand. We denote: ðIit Þ þ ¼ maxfIit ,0g, ðIit Þ ¼ maxfIit ,0g. Thus, the realized average cost per period, C, of the system over a horizon T is equal to 2 0 1 T N N N N X X X X 1 4X þ t t @ c F t tþ h ðI Þ þ ci di A C¼ T t ¼ 1 i ¼ 1 j ¼ 1,j a i ij Bi Mj i ¼ 1 i i i¼1 # T L X N T N X X X t t pi ðIi Þ þ li ðIi Þ : þ

Proposition A-1). The flow balance equations at nodes Rt will therefore simplify our network flow representation. The state spaces in the first L periods, in the last L periods, and in the middle [L+ 1,T L] periods are different. Let us consider the following four cases. (i) The first period: t ¼1. There are no replenishment order arrivals. For Period 1, inventory on hand at the beginning of period, H1i , is equal to the base stocks Si. (ii) The first L 1 periods: t ¼2,y,L. There are no replenishment order arrivals. Inventory on hand at the beginning of period, Hti , is the ending inventory of the previous period, namely FEt1 Bt . (iii) The

i

j

i

t ¼1i¼1

j

t ¼ TL þ 1 i ¼ 1

The terms in this cost function reflect the total transshipment cost, holding cost of excess inventory, total replenishment cost, total backlog cost, and total shortfall cost, respectively, over the entire planning horizon. For a base stock policy, replenishment quantities at location i in period t can be computed as follows: 8 Pt1 < ðSi ðFEt Bt þ 1 þ m ¼ tL FRm Em þ L FEt þ L Mt ÞÞ þ , t 4L, i i i i i FRt Et þ L ¼ ð1Þ P þ i : ðSi ðFEt Bt þ 1 þ t1 F , 1 r t r L: m þ L F t þ L t ÞÞ m m¼1 R E E M i

i

i

i

i

With a base stock policy and the assumptions above, the sum of the replenishment orders at all the locations is equal to the sum

i

i

i

middle periods: t ¼L+1,y,T L. This is the ‘‘typical’’ period; inventory on hand at the beginning of period, Hti , is the inventory from the last period, namely FEt1 Bt . Unmet demand is backlogged, i

i

and replenishment orders placed in period t will arrive in period t+ L. (iv) The last L periods: t ¼T L+ 1,y,T. Different from the middle periods, the shortfall is satisfied from a secondary (and more expensive) source as the replenishment orders, FRt Et þ L , in the i

last L periods cannot arrive in time within the finite horizon.

2.3.2. SP and LP representations We now introduce an SP model and an LP model, where the LP model is the deterministic counterpart of the SP model once demand is generated. The reason for building two models here is that we will use a stochastic approximation algorithm to compute the optimal base stock values. Since demand is a random variable, our problem is essentially a stochastic programming problem. We formulate this stochastic programming model in Problem (1) whose objective is to minimize


the expected average total cost per period E½CðS,DÞ given a status variable, the random demand D and for a decision variable, the base stock S¼(S1,S2,y,SN) in the system. Note FBt Mt Z 0,8i,j A N ,t A T . i

3. Algorithms and implementation 2

1 T

E½CðS,DÞ ¼ min E4

S ¼ ðS1 ,S2 ,...,SN Þ

S

þ

TL X

T X N X

N X

cij FBt Mt þ

t ¼ 1 i ¼ 1 j ¼ 1,j a i N X

T X

pi ðIit Þ þ

t ¼1i¼1

i

j

N X

T X N X

hi ðIit Þ þ

t ¼1i¼1

# li ðIit Þ :

t ¼ TL þ 1 i ¼ 1

Based on the stochastic programming Problem (1), and the network flow model presented in Fig. 1, we construct an LP formulation in Problem (2). Once demand D is observed as d in every period over the planning horizon T, Problem (2) becomes the deterministic counterpart of Problem (1). Through this LP formulation, we compute the derivative using duality – hence, avoiding cumbersome derivative recursions – simplifying the computation and improving algorithm efficiency. In addition, highly efficient LP packages exist to solve large-scale LP problems to implement our sample-path-based algorithm. Problem 2. 2 T X N N T X N X X 1 4X min CðS,dÞ ¼ min cij FBt Mt þ hi FBt Et i j i i S ¼ ðS1 ,S2 ,...,SN Þ S T t ¼ 1 i ¼ 1 j ¼ 1,j a i t ¼1i¼1 þ

T L X N X

N X

FBt Mt þ i

i

N X

FBt Mt þ i

i

i

i

FBt Mt þ FBt Et ¼ Si , i

j ¼ 1,j a i

j

i

i

FBt Mt þ FBt Et FEt1 Bt ¼ 0, i

j ¼ 1,j a i

j

i

i

T X

pi FEt þ L Mt þ

t ¼1i¼1

s:t:

i

i

N X

# li FEt þ L Mt

t ¼ TL þ 1 i ¼ 1

i

i

i ¼ 1, . . . ,N, t ¼ 1;

ð2Þ

i ¼ 1, . . . ,N, t ¼ 2, . . . ,T; ð3Þ

N X

FBt Mt þ i

i

N X

FBt Mt þ FEt þ L Mt ¼ dti , j

j ¼ 1,j a i N X

FRt Et þ L ¼ i

i¼1

i

i

i

dti ,

i

i ¼ 1, . . . ,N, t ¼ 1, . . . ,TL:

i ¼ 1, . . . ,N, t ¼ 1, . . . ,L;

i

i

FBt Et þFRtL Et FEt MtL FEt Bt þ 1 ¼ 0, i

i

i

i

i

i

i

FBt Et ,FBt Mt ,FBt Mt ,FEt Mt ,FRt Et Z0, i

i

i

i ¼ 1, . . . ,N, t ¼ 1, . . . ,T;

ð4Þ

ð5Þ

i¼1

FBt Et FEt Bt þ 1 ¼ 0, i

i

Proof. The proof is presented in the Appendix.

j

Problem 1. min

65

i

j

i

i

i

i ¼ 1, . . . ,N, t ¼ L þ 1, . . . ,T;

i ¼ 1, . . . ,N, t ¼ 1, . . . ,T:

ð6Þ ð7Þ ð8Þ

Eqs. (3)–(5), and (7) represent the physical inventory balance constraints at the nodes Bti , Mti , Rt, and Eti , respectively. This LP formulation will be at the heart of our algorithm; its feasibility is therefore a necessary condition for successful implementation. If all cost parameters, demand, and base stock levels are finite, then Problem (2) is feasible and has a finite optimum. This is established by Proposition 1. Proposition 1. If demand Di ,8i A N has a density on ð0,1Þ with E½Di o 1,8i A N , unit cost hi ,cij ,pi ,li A R þ and hi ,cij ,pi ,li o 1,8i,j A N , base stock Si A R þ and Si o1,8iA N then Problem (2) is feasible, and has a finite optimum with probability 1.

3.1. Algorithms To compute the optimal base stock values, we adopt a samplepath-based optimization algorithm, where we use IPA to compute the gradient value. In particular, we start with an arbitrary base stock level, Si, for stocking location i. After randomly generating an instance of the demand over T periods for each location, we construct and solve Problem (2) in a deterministic fashion. We then estimate the unbiased gradients of the average total cost. The LP is used not only to compute the optimal transshipment quantities, but also to help accumulate gradients ð@C=@Si ,i ¼ 1, . . . ,NÞ. With the unbiased gradients of the average total cost, we use a stochastic approximation algorithm to determine the optimal base stock levels. The procedure is summarized in a pseudo-code format, where K denotes the total number of steps taken in a path search, U represents the total number of iterations of the inner loop, dCu is the desired gradients of the average total cost at the uth step of the inner loop, ak represents the step size, which satisfies P1 P1 k 2 Condition 1: k ¼ 1 ak ¼ 1, and k ¼ 1 ak o 1, Si represents the base stock level for retailer i at the kth step. Algorithm 1. The algorithm. (I) (I) Initialization. Initialize K and U. Set initial base stock levels S0i ,8i A N . (II) Set k’1. Repeat. Set u’1. Repeat. (1) Generate the demand of T periods at each location from f(D). (2) Solve Problem (2) to determine optimal transshipment quantities. (3) Compute and accumulate the desired empirical gradients of the average cost dCu by invoking duality. (4) u’u þ 1. Until u ¼U. (5) Estimating unbiased gradients of the average cost by calculating the desired gradients of the average cost as PU 1 u ¼ 1 dC u . U (6) Update the base stock levels P ak U1 U Ski ’Sk1 u ¼ 1 dC u . i (7) k’k þ 1. Until k¼ K. (III) Return fSki ,8i A N ,k ¼ 1 : Kg and objective function value.

In Step (I), the algorithm starts with an arbitrary value for the base stock levels, S0i . K and U are also set by the experimenter. In general, the values of these parameters are problem dependent and can be determined, for instance, by a pilot study to mitigate the following trade-off: with a small number of steps, K, the experiment cannot provide sufficient data leading to a large variance in the output. On the other hand, an excessively large K can be wasteful in improving the optimal value. The outer loop in Step (II) includes the inner loop computations for {u ¼1:U}, the desired gradient calculations, and the

66


updating of base stock values. In Step (1), demand is generated at each retailer. Note that any covariance structure across the retailers is allowed in f(D). In Step (2), once the demand is observed, Problem (2) is solved in a deterministic fashion to compute the optimal transshipment quantities and the minimum-cost flow. In Step (3), the gradients of the average total cost with respect to the base stock levels are computed. Our LP formulation greatly simplifies these computations. Furthermore, the implementation of the derivative computation in this step is very efficient, as established in Theorem 1. Theorem 1. Based on the LP structure in our problem, the derivative of the average total cost with respect to base stock, @C=@Si , is equal to the corresponding dual optimal solution u*w, where w is determined by the position of Si in the LP formulation. For an N-location problem, w¼N + i. Proof. The proof is presented in the Appendix. The fact that, in a linear program, the dual value of a constraint is the derivative of the objective function with respect to the right-hand side of that constraint is well known (see, for instance, Swaminathan and Tayur, 1999). From Proposition 1 and Theorem 1, we have Corollary 1. This corollary will subsequently support Proposition 2. Corollary 1. If demand Di ,8i A N has a density on ð0,1Þ with unit cost hi ,cij ,pi ,li A R þ and hi ,cij ,pi ,li o 1,8i,j A N , base stock Si A R þ and Si o 1,8i A N , then the gradient of average cost with respect to base stock exists and is bounded.

E½Di o 1,8i A N ,

Proof. The proof is presented in the Appendix. In Step (5), we estimate the gradients of the average total cost, P which are actually random variables, by ð1=UÞ U u ¼ 1 dC u , where the gradient of interest is rS E½CðSÞ whereas our numerical procedure estimates E½rS CðSÞ. In other words, we use the expected value of the sample path derivative obtained via simulation instead of using the derivative of the expected cost. We therefore need to establish the unbiasedness of the gradient estimator, namely that the average of the changes represents the change in expectations. For our setting, in order to establish the unbiasedness of the gradient estimator, i.e., E½rS CðSÞ ¼ rS E½CðSÞ, we need to first prove several basic properties of the expected average total cost function. Lemmas 1 and 2 provide the basic conditions directly required by Proposition 2. Lemma 1. If demand Di ,8i A N has a density on ð0,1Þ with E½Di o 1,8i A N , finite unit cost hi ,cij ,pi ,li A R þ , base stock Si A R þ and Si o 1,8i A N , then E½CðSÞ is a proper convex function. Proof. The proof is presented in the Appendix. We also need to show that E½CðSÞ is continuously differentiable since, if E½CðSÞ is a proper convex function and the CDF of demand F(D) is continuous, then, from Rockafellar (1970, Theorem 25.2 and Corollary 25.5.2), we can obtain the following lemma. Lemma 2. If CDF of demand F(D) is continuous, E½CðSÞ is continuously differentiable everywhere. Proof. The proof is presented in the Appendix. As shown by Glasserman (1991), provided that the objective function E½CðSÞ is convex and continuously differentiable with respect to the base stock levels, IPA estimators will be unbiased. We can now establish Proposition 2. Proposition 2. If demand Di ,8i A N has a density on ð0,1Þ with PU u ¼ 1 dC u is

E½Di o 1,8i A N , then the gradient estimator ð1=UÞ

unbiased in the transshipment system with positive replenishment lead times, i.e., E½rS CðSÞ ¼ rS E½CðSÞ. Proof. The proof is presented in the Appendix. The base stock level Si is updated through Ski ’Sk1 i P ak ð1=UÞ Uu ¼ 1 dC u in Step (6). Also note that since the algorithm stops at k ¼K, we do not need an additional stopping rule. A key issue in this step is the selection of a suitable step size, ak . The rule of thumb for choosing ak is to let the step size go to zero fast enough so that the algorithm actually converges to a value of S, but not so fast that it will induce a wrong value. A widely accepted condition to meet that criterion is Condition 1. For instance, ak ¼ a=k for some fixed a 40 satisfies Condition 1. While we can theoretically use any step size satisfying Condition 1, the actual implementation requires further care; this is discussed in Section 3.2. Step (III) returns fSki ,8iA N ,k ¼ 1 : Kg and the objective function value at each step. This enables us to conduct the output analysis, which will also be discussed in Section 3.2. 3.2. Implementation We implemented the algorithm in Matlab, together with an LP solver. We wrote a main program to implement the sample path based algorithm, and input LP characteristic information to the LP solver. The LP solver returns the optimal base stock values and optimal average cost at each step to the main program. In implementing the algorithm, we need to select the following parameters: the initial value for base stock levels, Si0, the number of steps for the algorithm, K (the outer loop), the step size, ak , and the number of iterations, U (the inner loop). To examine the impact of the initial base stock values, we have experimented with different initial values along with different lead time settings. It turns out that initial values do not affect the final result, which shows the robustness of our algorithm. However, initial values do affect the convergence rate of the algorithm. There are several approaches proposed in the literature to determine the total number of steps for the algorithm to strike a balance between convergence and efficiency. For instance, Yin (1988) determines K through an ellipsoidal confidence region F A RN for S, with volume less than or equal to a c, so that, for each 0 o o o 1, limc-0 PfSi A F,i ¼ 1, . . . ,Ng ¼ 1o. His approach assumes a step size 1/k. Alternatively, one can stop as soon as the length of the confidence interval CIC for the sample means of average costs mC is less than or equal to 2d for some (small) d 4 0, i.e. CIC =mC r 2d, or equivalently, as soon as the estimated standard deviation of the sample mean is sufficiently small. This approach (e.g., see Garthwaite and Buckland, 1992; Fabian, 1968; Yin, 1990) assumes asymptotic normality of the sample means. Unfortunately, the effectiveness of any of these approaches is problem dependent. We have therefore conducted a comprehensive search, experimenting with different strategies using the illustrative examples from Krishnan and Rao (1965), Tagaras (1989), and Tagaras and Cohen (1992), where optimal and heuristic solutions are available. Based on this experimentation, we set the total number of steps for the path search K¼2000 even though the base stock values converged around k¼500 in most our experiments (in the slowest cases, it converged around k¼950). In this paper we define the convergence as jðSk þ 1 Sk Þ=Sk j r e, for a sufficiently small e which can satisfy the condition 0 o e o1. Even though all step sizes satisfying Condition 1 will theoretically lead to a correct result, different step sizes impact the convergence rate of the implemented algorithm (Kesten, 1952; L’Ecuyer et al., 1994). With the objective of optimally allocating a computational budget to successive SA iterations,


L’Ecuyer and Yin (1998) provide an asymptotically optimal choice of a matrix for step sizes. Their method can be used in situations with a tighter budget or with a larger scale. For our problem, the base stock decision is not a real-time decision. Since it is a planning decision, our computational budget is relatively large. After extensive experiments with step sizes from 0.003/k to 1000/ k, we have adopted ak ¼ 3=k as our step size in subsequent analysis. When the number of stocking locations and/or demand variance are very large, the algorithm will converge relatively slowly. In this case, we need to either increase K, which is timeconsuming, or increase step sizes. After extensive experiments, we retained K ¼2000 and adopted a step size of ak ¼ 3y=k, where y is used to capture the impact of the number of locations, N. Finally, we need to determine the number of iterations, U. Law and Kelton (2000, Chapter 9), offers an approximation for the minimum number of replications, U, that achieves a relative error smaller than g with probability 1b. However, in an SA algorithm, there is little need to control the relative error on the gradient estimator as the errors are averaged out across the SA iterations. Tang et al. (1999, p. 1829) therefore provide a very efficient method to set U as 1 in a computing environment same with ours, but they also state that setting U¼1 ‘‘appears natural, but not necessary for convergence’’. One can then adjust the step size to promote rapid convergence. Based on our pilot study, we mainly use U¼20 as the number of iterations and ak ¼ 3=k as the step size in our experiments. 4. Results In Section 4.1, we first validate our algorithm by comparing our results with those from Tagaras and Cohen (1992). Section 4.2 presents our solution and main experiments in more extensive settings with varying lead times and number of locations. In Section 4.3, we further explore the impact of demand variance. Section 4.4 studies the impact of correlated demand on transshipments with positive replenishment lead times. 4.1. Comparison studies We compare our optimal base stock values with the heuristic values reported in Tagaras and Cohen (1992). Their heuristic

160

S1, Tagaras & Cohen S2, Tagaras & Cohen

67

pffiffiffiffiffiffiffiffiffiffiffiffi algorithm is based on the model, Si ¼ mi ðLi þ 1Þ þ xi si Li þ1, to compute the base stock values, where mi ¼ EðDi Þ is the mean and s2i is the variance of the single-period demand at location i, while xi is computed by xi ¼ ðSLi ¼ 0 mi Þ=si , where SLi ¼ 0 is the optimal value for the problem with zero lead time. SLi ¼ 0 can be efficiently obtained from the existing literature (see Krishnan and Rao, 1965; Tagaras, 1989). In this experiment, we compare our optimal base stock values with the heuristic values from Tagaras and Cohen (1992). Since the Tagaras and Cohen algorithm is conceived for a two-location setting, we also apply our algorithm to the two-location setting. Base stock values in two locations determined by the two algorithms are presented in Fig. 2 (left panel), where we observe that: (i) Tagaras and Cohen (1992) heuristic tends to overestimate the optimal base stock values. All of their base stock quantities are consistently higher than ours. (ii) With increasing lead times, the overestimation by the Tagaras and Cohen algorithm becomes even more significant. Under a complete pooling transshipment policy, we further compare the objective function values from both algorithms. To ensure fairness, we use common random numbers across the two algorithms, and input optimal base stock values obtained by the Tagaras and Cohen algorithm and by our algorithm, and run 1000 independent replications to estimate the average total costs for both algorithms. It turns out that, for each lead time, the performance of our algorithm is better, consistently achieving ¨ lower average cost (Gong and Yucesan, 2009). We further compare our results with the results from a system without transshipments. The results are presented in Fig. 2 (right panel), where we observe that (i) for each lead time, the performance of the system with transshipments is always better than that of the system without transshipments even in the presence of replenishment lead times. The average total costs of systems without transshipments are always higher, provided that transshipment costs are not excessively large. (ii) When lead times are shorter, the results from the Tagaras and Cohen algorithm are closer to our results. But when lead times get longer, the results from the Tagaras and Cohen algorithm are closer to those from a system without transshipments, failing to reflect the added value of information pooling. From the comparison experiments, we find that the increased variability of inventory levels caused by replenishment lead times could be partially offset by the

120

S1, IPA & SPO S2, IPA & SPO

System without transshipment Tagaras & Cohen Algorithm

140 100

SPO & IPA Algorithm

120

Average cost

Base Stock

80 100 80 60

60

40 40 20 20 0

2

4

6

8 10 Lead Time

12

14

0

2

4

6

Fig. 2. Comparison of base stock values (left) and three settings (right).

8 10 Lead Time

12

14

68


opportunity to transship among the retailers, and would increase the attractiveness of transshipments. For a more robust comparison, we contrast our results with grid enumeration and with Tagaras and Cohen’s heuristic algorithms with the same cost parameters in Tagaras and Cohen (1992). The results are presented in Table 2, where data by ‘‘Grid’’ and ‘‘Heuristic’’ algorithms are reproduced from Tables 2 and 5, respectively, in Tagaras and Cohen (1992). We can observe that (i) all of the optimal values computed by Tagaras and Cohen’s grid enumeration (although the results of grid enumeration are still local optimal, they are relatively reliable) are covered by our 95% confidence intervals, implying that our results are relatively reliable in these experiments; (ii) we confirm that Tagaras and Cohen’s heuristic method performs well in cases where lead times are relatively short; and (iii) Tagaras and Cohen’s heuristic method tends to overestimate the values of base stocks. We have repeated the comparison experiments with different cost parameters (hi,cij,pi,li), experiment parameters (K,U), lead times L. The observations remain consistent.

We also present our optimal base stock and optimal average cost in Table 3, including the mean and the half width of a 95% confidence interval. From Table 3, we observe that both the average total cost and the optimal base stock levels will increase with longer replenishment lead times. In closing, we should note that our algorithm is equally effective in environments with different lead times in different locations. For detailed numerical results, the reader is referred to ¨ Gong and Yucesan (2009). 4.2.2. Varying the number of locations For a setting with a 2-period replenishment lead time and a 50-period horizon, we varied the number of locations from 2 to 12. The convergence of the optimal base stock values with the 12location scenario is depicted in Fig. 3, which illustrates the search process for the optimal base stock levels in those 12 locations. Table 3 Optimal base stock and average cost for each lead time.a L

4.2. Broader experiments 4.2.1. Varying lead times We consider a three-location problem with D1 Normalð10,5Þ, D2 Normalð5,2:5Þ, and D3 Normalð5,2:5Þ over a horizon of T¼50, and implement the simulation experiments with different replenishment lead times varying from 2 periods to 14 periods. For each simulation run (the complete optimization process for each lead time case), the elapsed computation time ranges from 4.56 to 7.91 min of wall clock time with a 2.6 GHz Pentium processor with 2 GB RAM memory.

2 4 6 8 10 12 14

S*1b

S*2c

S*3d

Average cost

34.42( 70.01) 53.73( 70.01) 71.72( 70.02) 88.51( 70.03) 104.80(7 0.02) 119.81( 70.03) 134.11( 70.04)

18.06( 70.01) 27.93( 70.01) 37.61( 70.01) 47.44( 70.02) 56.95( 70.02) 66.21( 70.02) 75.71( 70.02)

18.11( 7 0.01) 27.84( 7 0.01) 37.78( 7 0.02) 47.40( 70.02) 56.87( 7 0.03) 65.95( 7 0.02) 75.57( 7 0.02)

36.22( 7 0.21) 46.63( 7 0.27) 57.71( 7 0.32) 70.66( 70.35) 84.90( 70.38) 101.52( 7 0.40) 120.87( 7 0.42)

h ¼ 2,p ¼ 12,cij ¼ 5,l ¼ 20,K ¼ 2000, ak ¼ 3=k,S0 ¼ ð0,0,0Þ. Optimal base stocks in Location 1 with demand D1 Normalð10,5Þ. c D2 Normalð5,2:5Þ. d D3 Normalð5,2:5Þ. a

b

Table 2 Comparison with Tagaras and Cohen’s algorithms. Example No.

a

S1 (L1, L2)

b

SPO

S2 c

Grid

d

Heur.

e

Cost

SPO

Grid

Heur.

SPO

Grid

Heur.

Diff.f

1 2 3 4

(1,1) (1,1) (3,3) (3,3)

466.817 0.42 474.327 0.35 892.737 0.54 900.347 0.62

467 474 893 900

468 478 897 910

467.377 0.42 474.227 0.38 893.177 0.67 901.09 70.70

467 474 893 900

468 478 897 910

228.127 0.18 258.427 0.19 316.147 0.21 348.617 0.21

228.13 258.38 316.15 348.67

228.17 258.63 316.36 349.74

0.022% 0.081% 0.070% 0.320%

5 6 7 8

(1,1) (1,1) (3,3) (3,3)

418.457 0.42 421.777 0.38 827.117 0.58 827.567 0.57

419 422 827 828

420 422 829 831

477.987 0.43 480.77 7 0.41 904.78 7 0.64 911.237 0.74

477 480 905 911

476 482 907 916

160.02 70.17 171.857 0.18 223.827 0.21 236.327 0.22

160.03 171.84 223.88 236.53

160.03 171.88 223.92 236.79

0.006% 0.017% 0.045% 0.200%

9 10 11 12

(1,1) (1,1) (3,3) (3,3)

222.967 0.24 225.327 0.25 431.03 7 0.35 435.05 7 0.45

223 225 431 434

223 226 432 437

417.347 0.40 418.237 0.41 824.897 0.61 825.897 0.70

417 419 825 826

418 421 826 829

68.41 7 0.09 77.34 7 0.10 95.04 70.11 104.69 7 0.14

68.42 77.35 95.09 104.71

68.43 77.43 95.13 104.99

0.015% 0.120% 0.095% 0.290%

13 14 15 16

(1,1) (1,1) (3,3) (3,3)

224.127 0.21 225.897 0.28 434.597 0.37 437.01 7 0.41

224 226 434 436

225 227 435 438

472.777 0.41 479.347 0.39 902.89 7 0.74 907.56 7 0.76

473 479 903 908

474 481 904 915

164.287 0.16 178.527 0.17 229.427 0.18 244.687 0.21

164.31 178.57 229.44 244.73

164.33 178.64 229.48 245.11

0.030% 0.067% 0.026% 0.180%

17 18 19 20

(1,1) (3,1) (1,1) (3,1)

215.897 0.22 409.11 7 0.42 210.77 7 0.23 399.01 7 0.41

215 408 210 398

217 424 214 419

427.347 0.39 434.897 0.41 491.797 0.42 499.347 0.43

428 435 492 500

427 427 492 492

78.70 70.07 92.81 7 0.08 178.897 0.17 190.69 7 0.19

78.69 92.79 179.88 190.63

78.76 93.98 180.06 192.19

0.076% 1.260% 0.650% 0.790%

21 22 23 24

(1,1) (3,1) (1,1) (3,1)

187.567 0.25 348.777 0.47 171.997 0.19 339.01 7 0.42

187 348 171 338

196 394 181 373

452.237 0.40 469.897 0.44 524.217 0.48 532.897 0.49

453 470 525 533

444 444 516 516

73.65 7 0.06 84.52 7 0.08 170.30 70.16 178.687 0.17

73.62 84.43 170.34 178.54

74.04 87.94 170.75 180.84

0.530% 4.050% 0.260% 1.210%

a

‘‘No.’’, the index of examples, corresponding to Tagaras and Cohen (1992, Table 1). (L1, L2) denotes the lead times at two locations. c Data in ‘‘SPO’’ are the values computed by our algorithm. An estimator includes a mean and a HW for 95% CI. d Data in ‘‘Grid’’ are the values computed by Tagaras and Cohen’s grid enumeration. e Data in ‘‘Heur.’’ are the values computed by Tagaras and Cohen’s heuristic method. f ‘‘Diff.’’ is defined as jc1 c1 j=c2 100%, where c1 is the value from Tagaras and Cohen heuristic algorithm, and c2 is the mean value from SPO algorithm. b


We present the optimal base stock and the average total cost estimates for scenarios with different number of locations in Table 4, where we observe that the participation of a larger number of locations in transshipments enables each location to lower its own inventory and the total cost. We conclude that the algorithm works well with different lead times and different number of locations. In addition to the cases reported here, we have experimented with different demand distributions, different horizons, and different cost structures. We also conducted

69

experiments in a setting with different replenishment lead times for different retailers. In all these cases, the algorithm converges fast, and provides estimated values with a small variance. 4.3. The effect of demand variances In this section, we examine the performance of our algorithm in an environment with high demand variance and explore the effect of demand variance on optimal base stock values. As an

550 500

Base Stocks Si in 12 locations

450 400 350 300 250 S S S S S S S S S S S S

200 150 100 50 0

0

200

400 600 800 1000 1200 1400 1600 Steps in Sample−path−based Optimization

1800

2000

Fig. 3. Search paths of base stock values in a 12-location case.

Table 4 The optimal base stock and the average total cost for the scenarios with different number of locations. N

No. 1 N(5,2.5)a

2b 19.10cd 7 0.02 3 17.57 7 0.03 4 17.97 7 0.02 5 17.84 7 0.03 6 16.94 7 0.03 7 16.96 7 0.03 8 16.61 7 0.03 9 16.83 7 0.03 10 16.78 7 0.03 11 16.48 7 0.04 12 16.60 7 0.04

No. 2 (10,5) 35.46 7 0.02 35.67 7 0.01 35.71 7 0.02 35.40 7 0.02 34.13 7 0.03 34.05 7 0.03 34.02 7 0.03 33.49 7 0.03 33.70 7 0.03 33.24 7 0.04 32.61 7 0.05

No. 3 (15,7.5)

52.80 7 0.04 52.40 7 0.02 52.92 7 0.02 51.12 7 0.03 51.62 7 0.03 50.60 7 0.04 50.52 7 0.04 50.01 7 0.03 49.57 7 0.04 49.43 7 0.05

No. 4 (20,10)

67.68 70.03 68.56 70.02 68.23 70.03 67.11 70.03 67.85 70.03 67.58 70.04 66.83 70.04 66.40 70.05 66.45 70.05

No. 5 (25,12.5)

82.19 7 0.02 85.60 7 0.03 84.42 7 0.03 84.59 7 0.03 84.57 7 0.04 84.33 7 0.03 83.01 7 0.04 83.36 7 0.04

No. 6 (30,15)

101.07 7 0.04 100.81 7 0.03 100.77 7 0.04 101.79 7 0.04 100.37 7 0.03 99.66 7 0.05 99.85 7 0.05

No. 7 (35,17.5)

116.42 7 0.03 116.58 7 0.04 116.98 7 0.03 117.02 7 0.04 116.51 7 0.05 116.29 7 0.05

No. 8 (40,20)

131.41 7 0.05 132.03 7 0.04 132.96 7 0.03 133.60 7 0.04 132.68 7 0.05

No. 9 (45,22.5)

145.95 70.05 147.49 70.04 147.66 70.05 147.56 70.05

No. 10 (50,25)

163.53 7 0.05 163.53 7 0.06 164.62 7 0.04

No. 11 (55,27.5)

179.78 7 0.05 177.64 7 0.05

No. 12 (60,30)

Ave. cost

195.17 7 0.08

32.59 7 0.20 55.23 7 0.33 81.72 7 0.47 111.34 7 0.66 144.79 7 0.83 180.43 7 1.03 220.00 7 1.25 261.99 7 1.44 306.71 7 1.65 352.71 7 1.96 402.54 7 2.24

a It shows the demand distribution in one of 12 locations, for example, Nð60,30Þ denotes a N ormal distribution with a mean of 60 and a standard deviation of 30 or a variance of 900. b This scenario with two locations considers the locations ‘‘No. 1’’ and ‘‘No. 2’’. c h ¼2, p ¼12, cij ¼ 5, l ¼20, K¼ 2000. When both the number of locations and demand variance are large, we need to increase either K or ak . After extensive experiments, we recommend a heuristic algorithm of step sizes ak ¼ 3y=k, y ¼ ðN þ 1Þ=2 when both the number of locations and demand variances are large (N Z 6,CV ¼ s=m Z 0:5), enabling faster convergence of the Ski (within K¼2000). d Each estimated value includes a mean and a half width of a 95% confidence interval.

70


illustration, consider a three-location problem with three demand profiles, including a demand profile of high variance with D1 Normalð5,5Þ, D2 Normalð10,10Þ, and D3 Normalð15,15Þ, a demand of medium variance with standard deviations of 2.5, 5, 7.5 for locations 1, 2, and 3, respectively, and a demand of low variance with a standard deviation of 1 at all locations, over a horizon of 50 periods. We implemented the simulation experiments with the same cost parameters and different replenishment lead times of 2, 6, and 12 periods. We note that our algorithm’s performance is satisfactory in a high-variance environment; {Ski , k¼1:K} converge to the correct values within very short running time (the first 1000 steps) in all cases. We further note that demand variances significantly influence the optimal base stock values. For example, for the case with a 12-period lead time, the optimal base stock values (means of base stocks and their HWs of a 95% CI) are, respectively, 130.07 ( 70.01), 136.25 ( 70.02) and 151.23 ( 70.02) for a low, medium and high-variance demand profile, respectively. Note that in the Tagaras and Cohen algorithm, Si ¼ mi ðLi þ 1Þ þ pffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffi xi si Li þ 1 ¼ mi ðLi þ 1Þ þ ðSLi ¼ 0 mi Þ Li þ 1, since xi ¼ ðSLi ¼ 0 mi Þ=si . So SLi Z 1 is not a direct function of si , although SLi ¼ 0 could be a function of si . Therefore, it may not capture the change in demand variance when the lead time increases. This experiment shows another advantage of our algorithm, which can easily capture the impact of a change in demand variance. 4.4. The effect of demand correlation To assess the impact of correlated demand across retailers, we considered scenarios with different demand correlations and different lead times. A case with zero correlation is also added as a benchmark. Unlike the previous section, the demand faced by the retailers is modeled as a multivariate normal random variable with mean m and standard deviation si . The (i,j)th entry of the variancecovariance matrix is given by si sj rij , where rij denotes the level of demand correlation being investigated when i aj and one when i¼j. We adopt the system configuration 5 in Herer et al. (2006), where transshipments are possible between any pair of locations. We examine the impact of demand correlation on the average total cost per period. Table 5 shows the impact of demand correlation on the average total cost and the base stock values for a 3-retailer configuration for different lead times. We observe that, for each lead time, when demand correlation gets smaller (or negative), the average cost always decreases, which implies that the effectiveness of transshipments in matching demand and supply is enhanced. In general, we can conclude that positive

Table 5 The optimal average cost (and S*1) with different lead times L and correlation coefficients r.a L

r ¼ 0:5

r ¼ 0:25

r¼0

r ¼ 0:25

r ¼ 0:5

2

23.94b (300.77) 93.70 (599.78) 217.59 (899.70) 397.98 (1198.90) 631.45 (1498.30)

68.31 (305.88) 146.87 (596.04) 273.92 (885.13) 450.58 (1171.01) 677.10 (1452.40)

88.89 (307.81) 171.63 (595.51) 298.16 (879.38) 472.64 (1158.10) 697.74 (1433.81)

104.37 (308.72) 188.92 (594.39) 315.90 (873.98) 490.16 (1147.40) 712.88 (1419.61)

117.27 (310.05) 205.25 (594.32) 331.67 (869.63) 503.95 (1140.21) 719.93 (1405.60)

5 8 11 14

a

h ¼ 1, p ¼ 4, cij ¼ 0:5, l ¼ 10, m ¼ 100, si ¼ 20. The data outside of parentheses are average costs and inside of parentheses are the optimal base stocks S*1.

correlation reduces the effectiveness of transshipments while negative correlation enhances it.

5. Concluding remarks In this paper, we consider the multi-location transshipment problem with positive replenishment lead times. The main contributions of this paper include the following: (1) using simulation optimization combined with an LP/network flow formulation and IPA, we provide a flexible and efficient algorithm to compute the optimal base stock quantities for the multi-location transshipment problem with positive replenishment lead times; (2) experimenting with scenarios of high and low levels of demand correlation along with different lead times, we show the inverse relationship between the benefits of transshipments and demand correlation at different lead time settings; (3) we are also able to provide better objective function value than an existing heuristic algorithm; and (4) we introduce a duality-based gradient computation method to significantly improve computational efficiency. Our analysis of transshipment problems with positive replenishment lead times lends itself to potential extensions in several directions. For example, the research could be extended to transshipment problems with stochastic replenishment lead times. It may be possible to analyze this problem with a setting of discrete inventory quantities. It is also interesting to compare the IPA-based stochastic optimization algorithm with other stochastic programming algorithms.

Acknowledgements This research is supported by NSFC (no. 70901028). The authors would like to thank the editor and two referees for their constructive comments that have led to many improvements in the exposition of this paper.

Appendix

Proof of Proposition 1. (i) Problem (2) is feasible. We can always find a feasible solution. Let FRt Et þ L ¼ FEt þ L Mt ¼ dti for the first T L i i i periods and FEt þ L Mt ¼ dti for the last L periods. Si ¼ FB1 E1 ¼ i i i i FE1 B2 ¼ FB2 E2 ¼ , , ¼ FBt Et ¼ FEt Bt þ 1 ,8i A N . All transshipment quani i i i i i i i tities FBt Mt ¼ 0 with i,j A N . This set of values can always satisfy all i j constraints in Problem (2), so Problem (2) is always feasible. (ii) The optimum of the Problem (2) is finite w.p.1. Since cost vector hi ,cij ,pi ,li A R þ are finite and all decision variables are nonnegative, the objective value satisfies cT x Z0. Since it is a minimization problem, the Problem (2) has a finite optimum. & Proof of Theorem 1. Proposition 1 has shown that Problem (2) always has a finite optimum, and the optimal objective total cost value C* ¼cTBB 1b*, where B is a basis matrix, and b* is the righthand side column associated with the basis. Let b be the righthand side column, which consists of (3N + 1)T components. From the structure of the LP formulation, we note that the parameter matrix A is full rank with rank(A)¼(3N + 1)T. Since every component in right-hand side column is also in the b* associated with the basis B, that is b* ¼b. We have

b

C ¼ cTB B1 b:

ð9Þ


Also note that cTBB 1 ¼u, and u is the dual optimal solution. So we obtain C ¼ ðu ÞT b ¼

ð3N þ 1ÞT X

ðui ÞT bi :

ð10Þ

i¼1

By checking the LP formulation, we observe that bN + 1 ¼S1, bN + 2 ¼S2,y,b2N ¼SN, and S1,S2,y,SN only appear at these positions. Perturbing Si by an infinitesimal amount, b is perturbed to ~ because the difference bb~ is sufficiently small, B 1b remains b, positive w.p.1. and we still have a basic feasible solution. The reduced costs c are not affected and remain non-negative. Thus, the optimal basis B will not change, and Eq. (9) still holds. Then we take the derivatives on both sides of Eq. (10), and obtain @C @C @C @C @C @C ¼ ¼ uN þ 1 , ¼ ¼ uN þ 2 , . . . , ¼ ¼ u2N : @S1 @bN þ 1 @S2 @bN þ 2 @SN @b2N ð11Þ Therefore, the derivative @C=@Si equals the corresponding dual optimal solution u*w, where w is determined by the position of Si in the LP formulation. For N-location problem, w¼ N +i. & Proof of Corollary 1. Proposition 1 shows that both Problem (2) and its dual are feasible, our demand is both continuous and bounded, and @C=@Si is differentiable up to period T with probability 1. So the derivative of average cost with respect to base stock @C=@Si exists, and equals the corresponding dual optimal solution u*w. The boundedness of @C=@Si follows from the boundedness of its dual solution u*w. & Proof of Lemma 1. We firstly show E½CðSÞ is a convex function, then show it is a proper one. Let GðbÞ ¼ fxjAx ¼ b,x Z 0g be the feasible set, Q ¼ fbjGðbÞis non emptyg, and for any b A Q , we define CðbÞ ¼ minx A GðbÞ cT x, which is the optimal cost as a function b. Bertsimas and Tsitsiklis (1997) have shown that the objective functions C(b) of linear programs are convex functions of their right-hand sides b. In our LP formulation, all Si variables appear on the right-hand side of the linear programming, Si A b,8iA N . The convexity of the average cost in base stock levels S follows this property. The proper convexity of E½CðSÞ is followed by two points: (i) E½CðSÞ 41 for all S A R þ . If demand is finite and non-negative, unit costs are finite and non-negative, we have CðSÞ Z0 for all S A R þ . (ii) E½CðSÞo þ1 for some S A R þ (By contradiction). If it does not hold, Problem (2) will have no finite optimum, which contradicts Corollary 1. & Proof of Lemma 2. We use Rockafellar (1970, Theorem 25.2 and Corollary 25.5.1) to prove this lemma. We need two conditions: (i) E½CðSÞ is a proper convex function, which has already been proven by Lemma 1. (ii) Partial derivatives @E½C=@Si exist and are finite everywhere, which will be proven in the following. We have 2 0 1 Z T N N N N X X X X 1 4X þ t t @ E½CðSÞ ¼ cij FBt Mt þ hi ðIi Þ þ c i di A i j RNþ T t ¼ 1 i ¼ 1 j ¼ 1,j a i i¼1 i¼1 # T L X N T N X X X þ pi ðIit Þ þ li ðIit Þ dFðDÞ: t ¼1i¼1

t ¼ TL þ 1 i ¼ 1

Let W(S,D) be the term inside the integral; we have E½CðSÞ ¼ R RNþ WðS,DÞ dFðDÞ. Its partial derivative is Z @E½CðSÞ WðSi þ eei ,Di ÞWðSi ,Di Þ ¼ lim ð12Þ dFðDi Þ: e-0 RNþ @Si e

71

In order to prove the integral in Eq. (12) is absolutely convergent, now we will show that the term inside the brackets is bounded (Rockafellar, 1970). From Corollary 1, we know that the gradient of the average cost with respect to base stock exists and is equal to a finite u*w for any Si. We have j½WðSi þ eei ,Di Þ WðSi ,Di Þ=ejo 1. Since the term inside the brackets is bounded, the integral in Eq. (12) is absolutely convergent. We can put the limit inside the integral in Eq. (12). Also note that W(S,D) is convex and its partials must exist everywhere except at a countable number of points (Rockafellar, 1970, Theorem 25.3). Since F(D) is continuous, these points will have measure zero. Then we have Z @E½CðSÞ @WðSi ,Di Þ ¼ dFðDi Þ, ð13Þ @Si @Si RNþ Y where Y is the countable points set at which the partials @WðSi ,Di Þ=@Si do not exist; the probability measure of this set is zero. Then the partial derivatives exist everywhere outside of this zero-measure set Y, i.e., partial derivatives exist everywhere. For a proper convex function E½CðSÞ, if the partial derivatives @E½CðSÞ=@Si exist and are finite everywhere, then E½CðSÞ will be continuously differentiable everywhere. & Proof of Proposition 2. Lemma 1 has shown that E½CðSÞ is convex. Lemma 2 has shown that the derivatives of the average total cost with respect to the base stock levels, @E½C=@Si ,8iA N are continuous. Corollary 1 has shown that @E½C=@Si ,8iA N are finite. As shown by Glasserman (1991), IPA estimators will be unbiased provided that the objective function E½CðSÞ is convex and smooth (gradient is both continuous and finite) with respect to the base stock levels. & References Alfredsson, P., Verrijdt, J., 1999. Modelling emergency supply flexibility in a twoechelon inventory system. Management Science 45, 1416–1431. Allen, S.G., 1958. Redistribution of total stock over several user locations. Naval Research Logistics Quarterly 5, 337–345. ¨ Axsater, S., 2003. A new decision rule for lateral transshipments in inventory systems. Management Science 49, 1168–1179. Banerjee, A., Burton, J., Banerjee, S., 2003. A simulation study of lateral shipments in single supplier, multiple buyers supply chain networks. International Journal of Production Economics 81–82, 103–114. Bertsimas, D., Tsitsiklis, J.N., 1997. Introduction to Linear Optimization. Athena Scientific, Belmont. Cassandras, C.G., Wardi, Y., Melamed, B., Sun, G., Panayiotou, C.G., 2002. Perturbation analysis for online control and optimization of stochastic fluid models. IEEE Transactions on Automatic Control 47, 1234–1248. Eppen, G., 1979. Effects of centralization on expected costs in a multi-location newsboy problem. Management Science 25, 498–501. Fabian, V., 1968. On asymptotic normality in stochastic approximation. Annals of Mathematical Statistics 39, 1327–1332. Fu, M.C., 2002. Optimization for simulation: theory vs practice. INFORMS Journal on Computing 14, 192–215. Garthwaite, P.H., Buckland, S.T., 1992. Generating Monte Carlo confidence intervals by the Robbins–Monro process. Applied Statistics 41, 159–171. ¨ Gong, Y., Yucesan, E., 2009. The transshipment problem with positive replenishment lead times. Working Paper, INSEAD, France. Glasserman, P., 1991. Gradient Estimation via Perturbation Analysis. Kluwer Academic Publishers. Glasserman, P., 1994. Perturbation analysis of production networks. In: Yao, D.D. (Ed.), Stochastic Modeling and Analysis of Manufacturing Systems. SpringerVerlag, New York, pp. 233–280. Glasserman, P., Tayur, S., 1995. Sensitivity analysis for base stock levels in multiechelon production-inventory systems. Management Science 41, 263–281. Guinet, A., 2001. Multi-site planning: a transshipment problem. International Journal of Production Economics 74, 21–32. ¨ Herer, Y.T., Tzur, M., Yucesan, E., 2006. The multi-location transshipment problem. IIE Transactions 38, 185–200. ¨ Herer, Y.T., Tzur, M., Yucesan, E., 2002. Transshipments: an emerging inventory recourse to achieve supply chain leagility. International Journal of Production Economics 80, 201–212. Ho, Y.C., Eyler, M.A., Chien, T.T., 1979. A gradient technique for general buffer storage design in a serial production line. International Journal of Production Research 17, 557–580.

72


Kesten, H., 1952. Accelerated stochastic approximation. Annals of Mathematical Statistics 23, 463–467. Krishnan, K.S., Rao, V.R.K., 1965. Inventory control in N warehouses. Journal of Industrial Engineering 16, 212–215. Kukreja, A., Schmidt, C.P., Miller, D.M., 2001. Stocking decisions for low usage items in a multilocation inventory system. Management Science 47, 1371–1383. Law, A.M., Kelton, W.D., 2000. Simulation Modeling and Analysis, third ed. McGraw-Hill. Lee, H.L., 1987. A multi-echelon inventory model for repairable items with emergency lateral transshipments. Management Science 33, 1302–1316. Lee, Y.H., Jung, J.W., Jeon, Y.S., 2007. An effective lateral transshipment policy to improve service level in the supply chain. International Journal of Production Economics 106, 115–126. L’Ecuyer, P., Giroux, N., Glynn, P.W., 1994. Stochastic optimization by simulation: numerical experiments with the M/M/1 queue in steady-state. Management Science 40, 1245–1261. L’Ecuyer, P., Yin, G., 1998. Budget-dependent convergence rate of stochastic approximation. SIAM Journal on Optimization 8, 217–247. Mourania, I., Hennequina, S., Xie, X., 2008. Simulation-based optimization of a single-stage failure-prone manufacturing system with transportation delay. International Journal of Production Economics 112, 26–36.

Rockafellar, R.T., 1970. Convex Analysis. Princeton University Pre0ss, Princeton, NJ. Suri, R., Fu, B.-R., 1994. On using continuous flow lines to model discrete production lines. Discrete Event Dynamic Systems 4, 129–169. Swaminathan, J., Tayur, S., 1999. Stochastic programming models for managing product variety. In: Tayur, S., Ganeshan, R., Magazine, M. (Eds.), Quantitative Models for Supply Chain Management. Kluwer, pp. 585–622. Tagaras, G., 1989. Effects of pooling on the optimization and service levels of twolocation inventory systems. IIE Transactions 21, 250–257. Tagaras, G., Cohen, M., 1992. Pooling in two-location inventory systems with nonnegligible replenishment lead times. Management Science 38, 1067–1083. Tang, Q.-Y., L’ Ecuyer, P., Chen, H.-F., 1999. Asymptotic efficiency of perturbationanalysis-based stochastic approximation with averaging. SIAM Journal on Control and Optimization 37, 1822–1847. Veinott, A.F., 1965. Optimal policy for a multi-period dynamic non-stationary inventory problem. Management Science 12, 206–222. Wardi, Y., Melamed, B., Cassandras, C.G., Panayiotou, C.G., 2002. Online IPA gradient estimators in stochastic continuous fluid models. Journal of Optimization Theory and Applications 115, 369–405. Yin, G., 1988. A stopped stochastic approximation algorithm. Systems & Control Letters 11, 107–115. Yin, G., 1990. A stopping rule for the Robbins–Monro method. Journal of Optimization Theory and Applications 11, 107–115.