Stochastic lot-sizing problem with inventory-bounds and constant order-capacities

European Journal of Operational Research 207 (2010) 1398–1409

Contents lists available at ScienceDirect

European Journal of Operational Research journal homepage: www.elsevier.com/locate/ejor

Stochastics and Statistics

Stochastic lot-sizing problem with inventory-bounds and constant order-capacities Yongpei Guan a,*, Tieming Liu b a b

Department of Industrial and Systems Engineering, University of Florida, Gainesville, FL 32611, USA School of Industrial Engineering and Management, Oklahoma State University, Stillwater, OK 74078, USA

a r t i c l e

i n f o

Article history: Received 17 February 2010 Accepted 1 July 2010 Available online 11 July 2010 Keywords: Lot sizing Stochastic programming Inventory bounds

a b s t r a c t In this paper, we study the stochastic version of lot-sizing problems with inventory bounds and order capacities. Customer demands, inventory bounds, and costs are subject to uncertainty and dependent with each other throughout the finite planning horizon. Two models in stochastic programming are developed: the first one has inventory-bound constraints, and the second one has both inventory-bound and order-capacity constraints. We explore structural properties of the two models and develop Oðn2 Þ and Oðn2 T log nÞ dynamic programming algorithms for them, respectively. Our model also generalizes the deterministic lot-sizing problem with inventory bounds. For some cases, when applied to the deterministic versions, our algorithms outperform existing deterministic algorithms. Published by Elsevier B.V.

1. Introduction Lot-sizing (LS) problems with inventory bounds and order capacities (or production capacities) have gained a significant amount of attention recently. A major reason is that production-capacity and inventory-bound constraints are embedded in most of today’s inventory systems. As concepts such as ‘‘lean” and ‘‘just-in-time” are being applied to many stages in supply chains, surplus production capacity and inventory space are trimmed away or given way to other products. As a result, only a limited amount of (production, space, or budget) capacity might be available for a type of product, and this capacity might vary over time due to its own or other products’ seasonalities. In addition, as the structures of the supply chains become more complex than ever, some new issues, such as ‘‘remanufacturing option” (see, e.g., Richter and Sombrutzki, 2000; Richter and Weber, 2001), ‘‘production time-window constraints” (see, e.g., Wolsey, 2006), and ‘‘cumulative capacity constraints” (see, e.g., Sargut and Romeijn, 2007), have come into decision making. In a recent paper, van den Heuvel and Wagelmans (2008) showed that extended LS problems with the above three new conditions can be transformed into an LS problem with inventory bounds. Therefore, efficient algorithms for LS problems with inventory bounds and order capacities could also be used to solve many other practical problems effectively. Compared with their deterministic counterparts, stochastic inventory control problems are even more important in practice. In many real-world problems, customer demands, inventory bounds, and costs are all subject to uncertainty, and to make things worse, they could be dependent on each other in consecutive periods. For example, it is not surprising to see that a high-demand period follows a low-demand one, or a period with high fuel cost is followed by an even higher one. The ability to effectively address inventory control problems with stochastic and dependent parameters will help companies to improve their operational efficiency and thus the bottom line. The stochastic lot-sizing models with inventory bounds and order capacities are thus studied in this paper with this goal in mind. Sampling, an approach commonly used in practice, may work well to estimate the expected total cost for problems with demand uncertainties and correlations. However, it is shown in Ruszczyn´ski and Shapiro (2003) that the sampling approach only provides a lower bound for the expected total cost. Stochastic programming has the ability to address problems with demand uncertainties as well as demand dependencies among different time periods. Therefore, we take this approach to solve the stochastic LS problems with inventory-bound and order-capacity constraints. We would like to point out that the algorithms to solve stochastic LS problems using the stochastic programming approach are by no means easy extensions of their deterministic counterparts. Indeed, the algorithms to solve stochastic LS problems take a dramatically different approach from those to solve the deterministic ones. For a stochastic LS problem, the Wagner–Whitin property, which holds for the deterministic case, does not hold anymore, as described in Ahmed et al. (2003). Therefore, algorithms developed for the deterministic LS

* Corresponding author. Tel.: +1 (352) 392 1464; fax: +1 (352) 392 3537. E-mail addresses: [email protected]fl.edu (Y. Guan), [email protected] (T. Liu). 0377-2217/$ - see front matter Published by Elsevier B.V. doi:10.1016/j.ejor.2010.07.003

Y. Guan, T. Liu / European Journal of Operational Research 207 (2010) 1398–1409

1399

problems based on the Wagner–Whitin property (see, e.g., Wagelmans et al., 1992; Liu, 2008) cannot be applied here. This makes our problem much harder to solve. In this paper, we first analyze the stochastic LS problem with inventory bounds (SB–LS), and then study one with both inventory-bound and order-capacity constraints (SCB–LS). Although the first problem is a special case of the second, they are introduced separately, because the first problem itself has important practical applications (see, e.g., van den Heuvel and Wagelmans, 2008) and special properties that should be emphasized. In both models, except for the order capacities, all other parameters, including demands, inventory, and costs are stochastic and could be dependent on each other in different periods. We derive polynomial time algorithms for the SB–LS and SCB–LS problems in terms of input size. The rest of the paper is organized as follows: In Section 2, we review the literature on lot sizing in a broad scope that includes both deterministic and stochastic modeling approaches. In Section 3, we study the stochastic inventory control problem with inventory bounds. Then, in Section 4, we extend our study to the case for the stochastic inventory control problem with both inventory bounds and order capacities. Finally, we conclude and discuss extensions and special cases of the two models in Section 5.

2. Literature review There is a large amount of literature on lot-sizing problems under various conditions. In the following, we limit our scope to the models for a single product with setup and linear marginal ordering costs. We review the papers that are related to our work in three streams: models with deterministic demands, with stochastic but independent demands, and with stochastic and dependent demands and cost parameters. For each stream, we discuss the unconstrained and constrained models with order capacities and inventory bounds separately. The basic LS model seeks for the optimal order quantity in each period to minimize total ordering costs and holding costs for a single facility facing deterministic demands for a single product over a finite planning horizon. The replenishment and the holding cost functions at each period are deterministic and time-varying. The fundamental model was developed in Wagner and Whitin (1958), with a dynamic programming algorithm that runs in OðT 2 Þ time (where T is the number of time periods) based on the zero-inventory-ordering (ZIO) property, i.e., no replenishment is undertaken if inventory is available. Later on, Wagner (1960) and Zangwill (1968) showed that the ZIO property also holds under the assumption of concave order cost function. Wagelmans et al. (1992) improved the computational complexity for the linear cost case to OðT log TÞ (please also see Federgruen and Tzur, 1991; Aggarwal and Park, 1993), or OðTÞ for the case with non-increasing ordering cost in Wagelmans et al. (1992). The LS problem with order-capacity and inventory-bound constraints was studied by Love (1973). He showed that an optimal order plan has the property that between any two periods with positive but less-than-capacity production, there is at least one period whose ending inventory level is either zero or at the bound. Love also found an OðT 3 Þ dynamic programming algorithm for the problem with only inventory bounds. Gutiérrez et al. (2003) and Gutiérrez et al. (2007) derived a new characterization and OðT 3 Þ algorithms for the problem with bounded inventory, and with and without backlogging, respectively. Toczylowski (1995) developed an OðT 2 Þ algorithm by reducing multigraph edges in a shortest path graph. Recently, the computational complexity was improved to OðT 2 Þ by Liu (2008) and Atamtürk and Küçükyavuz (2008), and OðT log TÞ for the case with non-increasing ordering cost by Liu (2008). Florian and Klein (1971) studied the LS problem with constant production capacities, and they found an OðT 4 Þ algorithm. Later on, the computational complexity was improved to OðT 3 Þ by van Hoesel and Wagelmans (1996). Recently, Wolsey (2006) developed an OðT 4 Þ algorithm for the LS problem with both constant production capacities and production time-windows, which is equivalent to the case with both order capacities and inventory bounds as shown in van den Heuvel and Wagelmans (2008). Other related research works include the cases with demand time windows by Lee et al. (2001), the cases with general order capacities by Kirca (1990) and Chung et al. (1994), and the cases without setup costs by Ahuja and Hochbaum (2008). The stream of literature on inventory control policies with setup cost and stochastic independent demand is led by Scarf (1960). He showed that the expected cost function in each period satisfies the so-called K-convexity property and the existence of a non-stationary (sk, Sk) policy that is optimal for a finite horizon model. Under such a policy, an order is placed to increase the inventory position to Sk (the order-up-to level) if the inventory equals or falls below sk (the re-order point) in period k. Iglehart (1963a,b) proved the existence of stationary (s, S) policies over an infinite planning horizon for the time-discount and the time-average cost cases, respectively. Veinott and Wagner (1965) provided an algorithm to calculate the optimal s and S levels for the infinite-horizon problem. Later on, Zheng (1991) and Zheng and Federgruen (1991) revealed new insights and provided more efficient algorithms for the infinite-horizon problem. Federgruen and Zheng (1992) also studied the continuous-review (r, Q) policy, and they provided an efficient algorithm to calculate the optimal r and Q levels. Stochastic inventory control problems with both setup costs and order capacities and with independent demands are much more complex. Chen and Lambrecht (1996) studied the finite-horizon problem with setup cost and constant capacity constraints. A counter example was provided to show that (sk, Sk) policy is not optimal. They proved that the optimal policy instead satisfies the so-called X  Y band structure. If the inventory level is lower than X, order the full capacity; if the inventory level is above Y, order nothing; however, if the inventory is between X and Y, there is no fixed ordering pattern. Chan and Song (2003) studied the problem with both inventory bounds and order capacities. They further characterized the structure of the X  Y band policy with the so-called (a, b)-convexity, where a is the setup cost and b is the order capacity. They significantly reduced the search space within their dynamic programming framework, and as a result, a more efficient algorithm was found to calculate the optimal ordering policies. More recently, Zhao et al. (2007) studied the continuous-re~ Þ, the optimal view (r, Q) policy with inventory bounds and provided an algorithm to obtain the optimal policy (r*, Q*), starting from ð~r; Q ~ ; r and Q*. Most recent works on policy for the case without inventory bounds. The algorithm runs in polynomial time with respect to ~r; Q stochastic inventory control problems with independent demands also include different approximation algorithms developed by Halman et al. (2009) and Kunnumkal and Topaloglu (2008). Stochastic programming has been applied to solve inventory control problems with stochastic and dependent demands. In this stream, a scenario-tree-based formulation is developed to explore all possible scenarios in the future with the objective to minimize expected total costs. The focus is to develop efficient branch-and-cut algorithms to solve the scenario-tree-based formulation (see, e.g., Guan et al., 2006a,b; Di Summa and Wolsey, 2006), and find polynomial time algorithms in terms of the number of nodes in the stochastic

1400

Y. Guan, T. Liu / European Journal of Operational Research 207 (2010) 1398–1409

programming scenario tree. For the latter one, efficient algorithms have been studied for several unbounded stochastic LS problems using the stochastic scenario-tree-based formulations. Previous studies are the stochastic uncapacitated lot-sizing problem (SU–LS) with setup costs by Guan and Miller (2008), SU–LS without setup costs by Huang and Ahmed (2009), and SU–LS with random lead times by Huang and Küçükyavuz (2008). This paper is an extension of Guan and Miller (2008) with order capacities and inventory bounds. In Guan and Miller (2008), an Oðn2 maxfC; log ngÞ time algorithm was developed for the case with setup costs, where C is the maximum number of children of each node, and an Oðn2 Þ time algorithm was developed for the case without setup costs. This paper contributes to the literature with two models: one stochastic LS model with inventory bounds and the other with both inventory bounds and order capacities. We develop an Oðn2 Þ time algorithm for the first model and an Oðn2 T log nÞ time algorithm for the second model, where n is the number of nodes in the scenario tree. Deterministic LS model with inventory bounds studied by Liu (2008) and Atamtürk and Küçükyavuz (2008) is a special case of a stochastic LS problem with inventory bounds in which there is only one scenario for each period. Liu (2008) and Atamtürk and Küçükyavuz (2008) developed OðT 2 Þ time algorithms for this type of deterministic LS problem. Our algorithm also runs in OðT 2 Þ time, the same computational complexity as they provided. However, our algorithm provides this same computational complexity for the problem under a more general stochastic scenario-tree setting when backlogging is allowed. Deterministic LS model with inventory bounds and constant order capacities studied by Wolsey (2006) is a special case of a stochastic LS problem with both inventory bounds and order capacities, in which there is only one scenario for each period. As described in van den Heuvel and Wagelmans (2008), Wolsey (2006) provided an OðT 4 Þ time algorithm for this type of deterministic LS problem. Our approach improves the computational complexity to OðT 3 log TÞ. 3. Stochastic lot sizing with inventory bounds 3.1. Problem formulation In this section, we consider a single facility facing stochastic and time-varying demands and capacities on the inventory levels for the next T periods. Denote dt as the demand and ut as the capacity on the inventory level in period t, t = 1, . . . , T. Let ct be the unit purchasing cost þ and ft be the fixed ordering cost. At the end of each period, a holding cost, ht , is charged on each unit of inventory carried from period t to  t + 1, and a backlogging penalty ht is charged on each unit of backorders. Demand and other parameters for each time period depend on the history information up to this period. The demand for period t, dt, depends on the demand history from period 1 to period t  1. We adopt a stochastic scenario tree (see, e.g., Ruszczyn´ski and Shapiro, 2003) to represent demand scenarios and parameter realizations from period 1 to period T. Thus, all the problem parameters, including þ  dt ; ut ; ht ; ht ; ct , and ft evolve as a discrete time stochastic process with finite probability space, which is illustrated as a scenario tree of T stages as shown in Fig. 1. In the scenario tree, each node represents a state of the system that can be distinguished by information available up to its stage. For instance, each node in the leaf node set L represents one possible realization of demands from period 1 to period T. Except the root node (indexed as 1), each node i has an unique parent i. Let VðiÞ represent the set of all descendants of node i (including i itself), and Pði; ‘Þ denote the set of nodes on the path from node i to node ‘. For brevity, define V ¼ Vð1Þ and PðiÞ ¼ Pð1; iÞ. Let s(i) denote the time stage  or level of node i in the tree (i.e., sðiÞ ¼ jPðiÞj). Let CðiÞ denote the set of children nodes i has (i.e., CðiÞ ¼ fk 2 V : k ¼ ig) and C be the maximal number of children a node has (i.e., C ¼ maxi2V fjCðiÞjg). Let pi denote the probability associated with the state represented by node i. P Since there are several children explored from node i, the probability pi ¼ j2CðiÞ pj . Denote the cumulative demand along the path from P node i to node j as dij ¼ k2Pði;jÞ dk . We assume the n nodes in the tree are indexed such that d1 = d11 6 d12 6    6 d1n. We would like to point out that the three major issues that cause branchings in the tree, namely, demand forecast, inventory capacity, and production capacity, are all tactical but not operational decisions. They may be adjusted monthly, if not quarterly, but by no means daily. Production quantities could be adjusted weekly or daily in a firm. For example, in Fig. 1, the periods with branching nodes (e.g., time t1 and time t2) are the times when demand forecast, inventory capacity, or production capacity are adjusted. The non-branching nodes (all time indices except t1 and t2) represent the periods when only production decisions are needed. Therefore, for the scenario tree setting shown in Fig. 1, we do not require the scenario tree to be fully explored or each non-leaf node to have at least two children. Under its setting, the number of tree nodes may not go up exponentially as the time period increases. For instance, if there are only two time periods when the problem parameters need to be adjusted up or down for a T period problem, as shown in

Fig. 1. Multi-stage stochastic scenario tree.

Y. Guan, T. Liu / European Journal of Operational Research 207 (2010) 1398–1409

1401

Fig. 1, the number of nodes in the tree is OðTÞ. From an algorithmic point of view, in this paper, we will develop algorithms to obtain optimal solutions for the problem under a general scenario tree structure, i.e., CðiÞ P 1 for each node i 2 V n L. We also refer to an algorithm that can solve the problem under the general scenario tree structure as an algorithm that is regardless of the structure of the scenario tree. We assume zero initial inventory. A general stochastic LS problem whose initial inventory level is unknown and is itself a decision variable can be transformed into another general stochastic LS problem with zero initial inventory by adding a dummy root node as the parent of node 1, with zero order, zero setup and inventory costs, as well as zero demand. Let Yi be the replenishment quantity at time period s(i) corresponding to node i, and Ii be the net inventory left after time period s(i) corresponding to node i. Notice that Ii could be either positive or negative: when Ii P 0, Ii represents the amount of remaining inventory; when Ii < 0,  Ii represents the number of backorders. We can þ  use Hi ðIi Þ ¼ maxfhi Ii ; hi Ii g to represent the inventory/backlogging cost at time period s(i) corresponding to node i. The objective is to determine when to order and how many units to order so as to minimize total expected ordering and holding costs, without violating any inventory capacity constraints. Denote d(Y) as the impulse function, which is equal to 1 if Y > 0, and 0 otherwise (i.e., Y = 0). The problem can be formulated as follows:

ðSB-LSÞ

min

X

pi ðci Y i þ fi dðY i Þ þ Hi ðIi ÞÞ

i2V

s:t: Ii þ Y i ¼ di þ Ii ; Ii 6 ui ; Y i P 0;

i 2 V;

ð1Þ

i 2 V; i 2 V:

The first constraint is for the inventory flow balance, and the second constraint implies that the inventory at the end of each period P cannot exceed the capacity limit. It can be observed from the formulation that i2V pi ¼ T. For notation brevity, we incorporate pi into þ  parameters ci ; f i ; hi and hi , in the further analysis. From this formulation, it can be observed that the deterministic LS with inventory bounds is an SB–LS with only one scenario. 3.2. Optimality condition In this section, we develop a dynamic programming algorithm that runs in polynomial time in terms of the problem size, the number of nodes in the scenario tree. For this goal, we first characterize the structure of the optimal policy. Then we develop dynamic programming recursions within this structure. Finally, we develop an algorithm with polynomial time computational complexity. Before introducing the structure of the optimal policy, we first present an example as described in Fig. 2 to show that Wagner–Whitin property does not hold in general for SB–LS. As shown in the Figure, for each node, we assume that the unit production cost is zero, the unit backlogging cost is infinity, and the inventory bound is infinity. The probabilities, the setup costs, and the unit inventory costs for the six nodes are (1, 1, 0.5, 0.5, 0.5, 0.5), (0, +1+1, 0, 0,+1), and (100, 0, 0, 0, 0, 0), respectively. We assume the demand in each node is a positive constant number. In order to satisfy the demand at each time period and minimize the total cost, the productions are set up in nodes 1, 4, and 5. Accordingly, the production quantities Y1 = 3d, Y4 = d, and Y5 = d, where d is the demand at each node. Then, there is inventory left from parent node and production is set up for node 4 simultaneously. Thus, Y 4 I4 – 0. The Wagner–Whitin property for the deterministic case does not hold here. Note here that from this example, the inventory left from node 2 is the input for both nodes 3 and 4. Therefore, SB–LS is more complicated, as compared to its deterministic counterpart. We first derive the structure of the optimal policy for SB–LS as follows: Proposition 1. For a multi-period stochastic inventory control problem with inventory bounds, there exists an optimal solution (Y*, I*) in which for each node i 2 V, we have 1. Y i ¼ 0, or 2. Ii þ Y i ¼ dik for some node k 2 VðiÞ, or 3. Ii þ Y i ¼ dik þ uk for some node k 2 VðiÞ. Please refer to the Appendix for the proof. Proposition 1 implies that if an order arrives at node i, it is optimal to order an amount so that this amount together with inventory of the previous node will either cover the demands along the path from node i to a descendant of node i, or in addition, fill up the inventory capacity at the end of this descendant node. It is worth noting here that Proposition 1 generalizes the same conclusion derived in Love (1973) and Gutiérrez et al. (2007) into a scenario-tree setting. 3.3. Dynamic programming recursions In the following, we study a backward dynamic programming algorithm to solve the SB–LS problem, based on the results of Proposition 1. Let Ji(I) be the cost-to-go function, which denotes the expected cost of an optimal solution for the SB–LS problem with a planning horizon

Fig. 2. An example to show that Wagner–Whitin property does not hold for SB–LS.

1402

Y. Guan, T. Liu / European Journal of Operational Research 207 (2010) 1398–1409

from period s(i) to T and initial inventory of I, based on current scenario realization at node i (for notation brevity, we will use I rather than Ii for the second argument of this function). In the following, we define Gi(Y) as the replenishment cost function for Y units at period s(i) corresponding to node i

Gi ðYÞ ¼ ci Y þ fi dðYÞ:

ð2Þ

According to Proposition 1, to find an optimal solution, it is sufficient to consider only points in which there is zero or full-capacity inventory at the end of each subsequence. The cost-to-go function at node i with inventory I is calculated as,

( J i ðIÞ ¼ min Hi ðI  di Þ þ

(

X

J ‘ ðI  di Þ; min Gi ðdij  IÞ þ Hi ðdij  di Þ þ j2VðiÞ

‘2CðiÞ

(

min Gi ðdij þ uj  IÞ þ Hi ðdij þ uj  di Þ þ j2VðiÞ

X

) J ‘ ðdij  di Þ ;

‘2CðiÞ

)) J ‘ ðdij þ uj  di Þ

X

ð3Þ

:

‘2CðiÞ

The above three terms correspond to the three situations discussed in Proposition 1, respectively. At period s(i) at node i, the seller may do nothing, and follow an optimal strategy from period s(i) + 1 to T; or replenish the inventory to exactly cover the demand to node j; or replenish the inventory to cover the demand, and in addition, fill up the inventory capacity at node j. Before giving out the details of our algorithms, we assume d10 = 0 and u0 = 0 where node ‘‘0” is a dummy node, let V o ¼ V [ f0g, and introduce another important property. Proposition 2. For a multi-period stochastic inventory control problem with inventory bounds, if the initial inventory is zero, then there exists an optimal solution (Y*, I*) in which for each node i 2 V, we have 1. Ii ¼ d1k  d1i for some k 2 V o , or 2. Ii ¼ d1k þ uk  d1i for some k 2 V o . Proof. This proof can be done by an induction method starting from root node, i.e., node 1. Since the initial inventory is zero, according to Proposition 1, we have I1 ¼ d1k  d1 or I1 ¼ d1k þ uk  d1 for some k 2 V, or I1 ¼ d1 . The conclusion holds. h Assume the conclusion holds for a node i 2 V at time period s(i). That is, Ii ¼ d1k  d1i or Ii ¼ d1k þ uk  d1i for some k 2 V o . We only need to prove that the conclusion holds for each node j 2 CðiÞ. Based on Proposition 1, we verify all three possible scenarios for order at node j as follows: Case 1: Y j ¼ 0. For this case, Ij ¼ Ii  dj ¼ d1k  d1i  dj ¼ d1k  d1j or Ij ¼ Ii  dj ¼ d1k þ uk  d1i  dj ¼ d1k þ uk  d1j for some k 2 V o . 0 0 Case 2: Ii þ Y j ¼ djk0 for some k 2 VðjÞ. For this case, Ij ¼ Ii þ Y j  dj ¼ djk0  dj ¼ d1k0  d1j for some k 2 VðjÞ. 0 0      Case 3: Ii þ Y j ¼ djk0 þ uk0 for some k 2 VðjÞ. For this case, Ij ¼ Ii þ Y j  dj ¼ djk0 þ uk0  dj ¼ d1k0 þ uk0  d1j for some k 2 VðjÞ. In all above three cases, Proposition 2 is preserved in the induction step and therefore our claim is true. h In the following part, we derive a polynomial time algorithm for SB–LS by utilizing Propositions 1 and 2. For each node i 2 V, we first set Ji(I) = + 1 for the interval I > ui . For the case that I 6 ui , we further analyze the cost-to-go function as shown in (3). We have two options: order or do not order. If order occurs, then we define it as the order value function J Pi ðIÞ. As shown in (3), from Proposition 1, we have Yi = dik  I (or Yi = dik + uk  I) for some node k 2 VðiÞ such that dik > I (or dik + uk > I). Then we have

( J Pi ðIÞ

¼ min

þ

X

( min

Gi ðdik  IÞ þ Hi ðdik  di Þ þ

k2VðiÞ:di þui Pdik >I

(

) J ‘ ðdik  di Þ ;

‘2CðiÞ

)) J ‘ ðdik þ uk  di Þ

X

min

k2VðiÞ:di þui Pdik þuk >I

Gi ðdik þ uk  IÞ þ Hi ðdik þ uk  di Þ

:

‘2CðiÞ

To simplify future analysis, we denote þ

/0i ðkÞ ¼ fi þ ci dik þ hi ðdik  di Þ þ

X

J ‘ ðdik  di Þ and

‘2CðiÞ þ

/ui ðkÞ ¼ fi þ ci ðuk þ dik Þ þ hi ðuk þ dik  di Þ þ

X

J ‘ ðuk þ dik  di Þ:

ð4Þ ð5Þ

‘2CðiÞ

Based on (2), (4) and (5), corresponding to each inventory level I, we have JPi ðIÞ ¼ /0i ðkÞ  ci I or JPi ðIÞ ¼ /ui ðkÞ  ci I for some node k 2 VðiÞ. For a given particular k 2 VðiÞ; /0i ðkÞ and /ui ðkÞ are constant numbers. Thus, as shown in the example in Fig. 3, the order value function J Pi ðIÞ is a piecewise linear and right continuous function with the same slope,  ci, for each piece. 0 0 0 For brevity, we further use /i(k ) to represent both /0i ðkÞ and /ui ðkÞ. We have /i ðk Þ ¼ /0i ðkÞ if k0 = k and /i ðk Þ ¼ /ui ðkÞ if k0 = k + n. Accord0 0 ingly, we set dik0 ¼ dik if k = k and dik0 ¼ dik þ uk if k = k + n. Since J Pi ðIÞ is a piecewise linear function with the same slope, ci, for each piece, we assume there are r linear pieces and then each piece of the order value function can be described as follows:

Y. Guan, T. Liu / European Journal of Operational Research 207 (2010) 1398–1409

1403

Fig. 3. The order value function JPi ðIÞ for node i.

0

0

0

0

0

J Pi ðIÞ ¼ /i ðk1 Þ  ci I if I < dik01 ; where k1 2 argminf/i ðk Þ : k 2 VðiÞ or k  n 2 VðiÞg; 0

J Pi ðIÞ ¼ /i ðk2 Þ  ci I if dik01 6 I < dik02 ; 0

0

0

0

where k2 2 argminf/i ðk Þ : k 2 VðiÞ or k  n 2 VðiÞ and dik0 > dik01 g; .. . 0

J Pi ðIÞ ¼ /i ðkr Þ  ci I if dik0r1 6 I < dik0r ; 0

0

0

0

where kr 2 argminf/i ðk Þ : k 2 VðiÞ or k  n 2 VðiÞ and dik0 > dik0r1 g: Otherwise, if no order occurs, then we will call this function the do-not-order value function JNP i ðIÞ. It contains inventory holding cost corresponding to this node and the cost for later periods. This function can be expressed as

J NP i ðIÞ ¼ Hi ðI  di Þ þ

X

J ‘ ðI  di Þ:

ð6Þ

‘2CðiÞ

Therefore, we have (i) If I P maxk2VðiÞ dik , then order is not needed and J i ðIÞ ¼ J NP i ðIÞ. (ii) If I < di, then there is not enough inventory available for the demand at node i. Thus, order is required at this node and J i ðIÞ ¼ J Pi ðIÞ. (iii) If di 6 I < maxk2VðiÞ dik , then both order and do-not-order cases are possible and therefore

n o J i ðIÞ ¼ min J Pi ðIÞ; J NP i ðIÞ :

ð7Þ

3.4. Algorithm and computational complexity analysis For a given zero initial inventory, we calculate and store the value function J1(0). Based on Propositions 1 and 2, for each node i 2 V, it is sufficient to calculate and store Ji(I) for I ¼ d1k  d1i or I ¼ d1k þ uk  d1i for all k 2 V o . In general, as shown in Fig. 4, the inventory level I ¼ d1i for the root node value function J1(I) corresponds to the inventory level I = 0 for Ji(I). Similarly, the inventory level I ¼ d1‘ for the root node value function J1(I) corresponds to the inventory level I = 0 for J‘(I). More specifically, based on (6) and (7), the breakpoints for the value function Ji(I) are obtained by moving the breakpoints for the value function J‘(I) for all ‘ 2 CðiÞ to the right by di. In order to store the breakpoints I ¼ d1k  d1i and I ¼ d1k þ uk  d1i for all k 2 V o for node i, in the calculation of the value functions J‘(I) for all ‘ 2 CðiÞ, we only need to store the breakpoints I ¼ d1k  d1i ¼ d1k  d1‘ and I ¼ d1k þ uk  d1‘ for all k 2 V o . Therefore, the breakpoints I = d1k and I = d1k + uk for all k 2 V o for the root node can be used to calculate and store Ji(I) for each i 2 V, and there are OðnÞ breakpoints. The value functions Ji(I) are calculated and stored backwards starting from the last time period T. Before exploring the details of the algorithm, in the initial step, we (1) set relationship indicator g(i, k) between each pair of nodes (i, k) to indicate if k 2 VðiÞ:gði; kÞ ¼ 1 if node k 2 VðiÞ and g(i, k) = 0, otherwise;

Fig. 4. Shifting of breakpoints.

1404

Y. Guan, T. Liu / European Journal of Operational Research 207 (2010) 1398–1409

P

(2) define h(i, I) to store ‘2CðiÞ J ‘ ðIÞ for which I = d1k  d1i or I = d1k + uk  d1i for all k 2 V o and initialize them to zero; (3) sort the breakpoints I = d1k and I = d1k + uk in a non-decreasing order. From above, it can be observed that g(i, k) can be obtained in Oðn2 Þ time using the tree data structure. It is also easy to observe that (2) and (3) can be completed in Oðn2 Þ time and Oðn log nÞ time respectively. Thus, the initial step can be completed in Oðn2 Þ time.  Now we calculate and store Ji(I) by induction starting from the last time period T. For each node i in time period T, if ui > di and ci < hi , then the value function is

8 > ci ðdi  IÞ þ fi > > > <  hi ðdi  IÞ J i ðIÞ ¼ > > hþi ðI  di Þ > > : þ1

^; if I 6 d i ^ if d < I < d ; i

i

if di 6 I 6 ui ; if I > ui ;

^ ¼ d  f =h  c : If u  6 d or c P h , then we can obtain similar value functions with less number of linear pieces. where d i i i i i i i i i Thus, it is easy to verify that Ji(I) where I ¼ d1k  d1i or I ¼ d1k þ uk  d1i for each node k 2 V o can be calculated and stored in OðnÞ time.   Meanwhile, increase hði ; d1k  d1i Þ and hði ; d1k þ uk  d1i Þ by J i ðd1k  d1i Þ and J i ðd1k þ uk  d1i Þ, respectively for each node k 2 V o . That is,  each individual point in h(i , I) is increased by Ji(I). This step takes OðnÞ time. For each node i in the induction step, the following calculations will be executed:  Step 1: Calculate and store J NP I ¼ d1k þ uk  d1i for each node k 2 V o . After we finish calculating and storing J‘(I) for i ðIÞ for I ¼ d1k  d1i and P each node ‘ 2 CðiÞ, we have obtained hði; IÞ ¼ ‘2CðiÞ J ‘ ðIÞ for I = d1k  d1i and I = d1k + uk  d1i for each node k 2 V o . Thus, for each par0 ticular breakpoint (i.e., corresponding to a particular inventory level I0 ) in J NP i ðIÞ, the corresponding breakpoint (i.e., I  di) is stored  and I ¼ d  for each node k 2 V in h(i, I). Based on (6), J NP ðIÞ for I ¼ d  d þ u  d can be obtained in OðnÞ time. o 1i 1i 1k 1k k i Step 2: Calculate /0i ðkÞ and /ui ðkÞ according to (4) and (5), and store them for each k 2 VðiÞ. As in Step 1, the values of the functions P P hði; dik  di Þ ¼ ‘2CðiÞ J ‘ ðdik  di Þ and hði; dik þ uk  di Þ ¼ ‘2CðiÞ J ‘ ðdik þ uk  di Þ for a given k 2 VðiÞ have been obtained in the previous step when we finish calculating and storing J‘(I) for each node ‘ 2 CðiÞ. Thus this step can be finished in OðnÞ time. Step 3: Calculate and store J Pi ðIÞ. The slope for each piece is fixed (i.e., ci). As mentioned before, we let /i(k0 ) represent /0i ðkÞ and /ui ðkÞ. We initialize / = + 1. Starting from I ¼ d1n þ un  d1i , we calculate and store J Pi ðIÞ for I ¼ d1k  d1i and I ¼ d1k þ uk  d1i for each k 2 V, respectively. (a) Corresponding to each breakpoint I ¼ d1k0  d1i (i.e., I ¼ d1k  d1i or I ¼ d1k þ uk  d1i ) such that either g(i, k0 ) or g(i, k0  n) (if 0 0 0 0 k0  n > 0) is equal to 1, we first set J Pi ðd1k0  d1i Þ ¼ /  ci ðd1k0  d1i Þ. Then let /i ðk Þ ¼ /0i ðk Þ if k0 6 n and /i ðk Þ ¼ /ui ðk  nÞ 0 0 otherwise. If /i(k ) < /, update / = /i(k ). (b) Corresponding to each breakpoint I ¼ d1k0  d1i such that both g(i, k0 ) and g(i, k0  n) are equal to 0, we let J Pi ðd1k0  d1i Þ ¼ /  ci ðd1k0  d1i Þ. This step can be completed in OðnÞ time; Step 4: Calculate and store Ji(I) for I ¼ d1k  d1i and I ¼ d1k þ uk  d1i for each node k 2 V o . This step can be completed in OðnÞ time according to Eq. (7), since the number of breakpoints is bounded by OðnÞ.   Step 5: Update h(i, I) for I ¼ d1k  d1i and I ¼ d1k þ uk  d1i for each node k 2 V o . Increase hði ; d1k  d1i Þ and hði ; d1k þ uk  d1i Þ by J i ðd1k  d1i Þ and J i ðd1k þ uk  d1i Þ, respectively for each node k 2 V o . This step can be completed in OðnÞ time.

It can be observed from above that the total amount of work required for each node is bounded by OðnÞ. Since these operations are required for each node i 2 V, the following conclusion holds. Theorem 1. The general stochastic inventory control problem with inventory bounds can be solved in Oðn2 Þ time, regardless of the structure of the scenario tree. It is also easy to observe that the space complexity for the algorithm is Oðn2 Þ as well. 4. Stochastic lot sizing with inventory bounds and constant order capacities In this section, we investigate the computational complexity to obtain an optimal policy for a case in which both inventory bounds and order capacities exist, regardless of the structure of the scenario tree. We assume there is a constant order capacity q in each period. Compared with the formulation for SB–LS, the formulation for this problem has only one more capacity constraint, Yi 6 q, for any i 2 V, and thus the whole formulation is omitted. Note here, we allow backlogging at each time period including time period T. Therefore, the problem always has a feasible solution. 4.1. Optimality condition We have the following proposition parallel to Proposition 1. Proposition 3. For a multi-period stochastic inventory control problem with inventory bounds and constant order capacity q, there exists an optimal solution (Y*, I*) such that for each node i 2 V, we have 1. Y i ¼ 0 or q or ^ for some node k 2 VðiÞ and an integer m ^ such that 0 6 m ^ 6 sðkÞ  sðiÞ and Y j ¼ 0 or q for each j 2 PðkÞ n PðiÞ, or 2. Y i þ Ii ¼ dik  mq   ^ for some node k 2 VðiÞ and an integer m ^ such that 0 6 m ^ 6 sðkÞ  sðiÞ and Y j ¼ 0 or q for each j 2 PðkÞ n PðiÞ. 3. Y i þ Ii ¼ uk þ dik  mq

1405

Y. Guan, T. Liu / European Journal of Operational Research 207 (2010) 1398–1409

Please refer to the Appendix for the proof. Similar to Propositions 1, 3 generalizes the same conclusion derived in Love (1973) and Gutiérrez et al. (2007) into a scenario-tree setting. Based on this Proposition, set the value function of the problem for node i 2 V; J i ðIÞ ¼ þ1 as I > ui . Otherwise, there are three cases which can be described as follows:  Do not order:

J NP i ðIÞ ¼ Hi ðI  di Þ þ

X

J ‘ ðI  di Þ:

ð8Þ

‘2CðiÞ

 Order at capacity:

J Pq i ðIÞ ¼ fi þ ci q þ H i ðq þ I  di Þ þ

X

J ‘ ðq þ I  di Þ:

ð9Þ

‘2CðiÞ

 Order less than capacity:

J Pi ðIÞ ¼ min ^



min

^

^ sðkÞsðiÞ:dik ðmþ1Þq6I
m;0 J k; ðIÞ; i

min

^ sðkÞsðiÞ:uk þdik ðmþ1Þq6I
^

m;u J k; ðIÞ i

 ð10Þ

^

where Jik;m;0 ðIÞ (or Jik;m;u ðIÞ) is the cost function to order a less-than-capacity amount such that the inventory level after node k is equal to zero (or the inventory bound), before making any other less-than-capacity order. That is, ^ m;0 ^  IÞ þ Hi ðdik  mq ^  di Þ þ J k; ðIÞ ¼ fi þ ci ðdik  mq i

X

^  di Þ; and J ‘ ðdik  mq

‘2CðiÞ ^

m;u ^  IÞ þ Hi ðuk þ dik  mq ^  di Þ þ J k; ðIÞ ¼ fi þ ci ðuk þ dik  mq i

X

^  di Þ: J ‘ ðuk þ dik  mq

ð11Þ ð12Þ

‘2CðiÞ

Pq P Finally, the optimal value function can be described as J i ðIÞ ¼ minfJ NP i ðIÞ; J i ðIÞ; J i ðIÞg:.

As in SB–LS, we assume zero initial inventory. We can obtain the following proposition. Proposition 4. For a multi-period stochastic inventory control problem with inventory bounds and constant order capacity q, if the initial inventory is zero, then there exists an optimal solution (Y*, I*) in which for each node i 2 V, we have 1. Ii ¼ d1k  d1i  mq with  s(i) 6 m 6 s(k)  s(i) for some node k 2 V o , or 2. Ii ¼ uk þ d1k  d1i  mq with  s(i) 6 m 6 s(k)  s(i) for some node k 2 V o .

^ and I1 ¼ d1k  mq ^  d1 , or Y 1 ¼ uk þ d1k  mq ^ and Proof. With zero initial inventory, according to Proposition 3, we have Y 1 ¼ d1k  mq ^  d1 for some node k 2 V o and an integer m ^ such that 0 6 m ^ 6 sðkÞ  sðiÞ, or Y 1 ¼ q and I1 ¼ q  d1 , or Y 1 ¼ 0 and I1 ¼ uk þ d1k  mq I1 ¼ d1 , in an optimal solution. Thus, the conclusion holds for the root node. Assume the conclusion holds for a node i 2 V at time period s(i). That is, Ii ¼ d1k  d1i  mq or Ii ¼ uk þ d1k  d1i  mq with s(i) 6 m 6 s(k)  s(i) for some node k 2 V o in an optimal solution. Then, for each node ‘ 2 CðiÞ, based on Proposition 3, we verify three possible cases for order at node ‘as follows: Case 1: Y ‘ ¼ 0. For this case, since the order amount at node ‘is not equal to q, m 6 s(k)  s(i)  1 = s(k)  s(‘). Then, I‘ ¼ Ii  d‘ ¼ d1k  d1i  mq  d‘ ¼ d1k  d1‘  mq or I‘ ¼ Ii  d‘ ¼ uk þ d1k  d1i  mq  d‘ ¼ uk þ d1k  d1‘  mq for some node k 2 V o with s(‘) <  s(i) 6 m 6 s(k)  s(‘). Case 2: Y ‘ ¼ q. For this case, I‘ ¼ Ii þ q  d‘ ¼ d1k  d1i  ðm  1Þq  d‘ ¼ d1k  d1‘  ðm  1Þq or I‘ ¼ Ii þ q  d‘ ¼ uk þ d1k  d1i  ðm  1Þq  d‘ ¼ uk þ d1k  d1‘  ðm  1Þq for some node k 2 V o with  s(‘) =  s(i)  1 6 m  1 6 s(k)  s(‘). ^ or Ii þ Y ‘ ¼ uk0 þ d‘k0  mq ^ for some node k0 2 Vð‘Þ and an integer m ^ such that 0 6 m ^ 6 sðk0 Þ  sð‘Þ. Then, Case 3: Ii þ Y ‘ ¼ d‘k0  mq 0       ^ ^ 0 0 0 0 0 0 I‘ ¼ Ii þ Y ‘  d‘ ¼ d‘k  mq  d‘ ¼ d1k  d1‘  mq or I‘ ¼ Ii þ Y ‘  d‘ ¼ uk þ d‘k  mq  d‘ ¼ uk þ d1k  d1‘  mq for some k 2 V o and s(‘) < 0 6 m 6 s(k0 )  s(‘). In all above three cases, Proposition 4 is preserved in the induction step, and therefore our claim is true. h Proposition 4 simplifies the algorithm. Based on this proposition, we only need to store the breakpoints for node i that are in the form I ¼ Ii ¼ d1k  d1i  mq or I ¼ uk þ d1k  d1i  mq with  s(i) 6 m 6 s(k)  s(i) for all k 2 V o . Based on (8)–(10), the breakpoints for the value function Ji(I) are obtained by moving the breakpoints for the value function J‘(I) for all ‘ 2 CðiÞ to the right by di or by di  q, plus two new breakpoints I = di and I = di  q. Thus, in order to store the breakpoints I ¼ d1k  d1i  mq and I ¼ uk þ d1k  d1i  mq with  s(i) 6 m 6 s(k)  s(i) for all k 2 V o , in the calculation of the value functions J‘(I) for all ‘ 2 CðiÞ, we only need to store the breakpoints I ¼ d1k  d1i  mq ¼ d1k  d1‘  mq and I ¼ uk þ d1k  d1i  mq ¼ uk þ d1k  d1‘  mq for all k 2 V o and  s(‘) 6 m 6 s(k)  s(‘). Therefore, in the dynamic programming framework to calculate the value function Ji(I) for each node i 2 V backwards starting from leaf nodes, we store the breakpoints I ¼ d1k  d1i  mq and I ¼ uk þ d1k  d1i  mq and their evaluations for all k 2 V o and s(i) 6 m 6 s(k)  s(i). There are at most OðnTÞ breakpoints needed to be stored for all value functions J i ðIÞ; i 2 V.

1406

Y. Guan, T. Liu / European Journal of Operational Research 207 (2010) 1398–1409

4.2. Algorithm and computational complexity analysis Similar to SB–LS, the scenario tree structure can help speed up the algorithm for SCB–LS. For the initial step, we set relationship indicator P

g(i, k) for each pair of nodes i 2 V and k 2 VðiÞ. We set g(i, k) = 1 if k 2 VðiÞ and g(i, k) = 0 otherwise. We use h(i, I) to store ‘2CðiÞ J‘ ðIÞ for which I ¼ d1k  d1‘  mq or I ¼ uk þ d1k  d1‘  mq for all k 2 V o and s(‘) 6 m 6 s(k)  s(‘) and initialize them to zero. Note here this initial

step can be completed in Oðn2 TÞ time. Then, we sort the breakpoints for the root node I0 = d1k  mq and I0 = uk + d1k  mq for each node k 2 V o and T 6 m < T in a nondecreasing sequence. This can be completed in OðnT log nÞ time because T 6 n and Oðlog n2 Þ is equal to Oðlog nÞ. Note here that, these breakpoints are enough to serve as the breakpoints of the value function Ji(I) for each node i 2 V and are defined as the predetermined breakpoint list. Corresponding to each node i 2 V, we set I = 0 at the breakpoint I0 ¼ d1i . Thus, the breakpoints contained in J i ðIÞ; I ¼ d1k  d1i  mq and I ¼ uk þ d1k  d1i  mq for each node k 2 V o and  s(i) 6 m 6 s(k)  s(i), can be stored at the breakpoints I0 = d1k  mq or I0 = uk + d1k  mq for each node k 2 V o and s(i) 6 m 6 s(k). Finally, duplicated breakpoints are considered separately. By taking advantage of both scenario tree structure and the special value functions (8)–(12), the optimal solution for the stochastic inventory control problem with inventory bounds and order capacities can be obtained by induction starting from the last time period T in the following approach.  For each node i in time period T, similar to the case for SB–LS, if ui > di and ci < hi , then the value function is

8  fi þ ci q þ hi ðdi  q  IÞ if I < di  q; > > > > ^; > if di  q 6 I < d > i < fi þ ci ðdi  IÞ  ^ J i ðIÞ ¼ hi ðdi  IÞ if d 6 I < d ; i i > > > > hþi ðI  di Þ if di 6 I 6 ui ; > > : þ1 if I > ui ; ^i ¼ di  fi =ðh  ci Þ. If ui 6 di or ci P h , then we can obtain similar value functions with less number of linear pieces. Then, the where d i i breakpoints I ¼ d1k  d1i  mq and I ¼ uk þ d1k  d1i  mq for all k 2 V o and s(i) 6 m 6 s(k)  s(i) of Ji(I) can be achieved in OðnTÞ time since the number of breakpoints to be stored is bounded by OðnTÞ. The value function h(i, I) is updated in a way that the function value corresponding to each breakpoint I ¼ d1k  d1i  mq and I ¼ uk þ d1k  d1i  mq for all k 2 V o and  s(i) 6 m 6 s(k)  s(i) is increased by Ji(I). This procedure can also be completed in OðnTÞ time. For each induction step for a node i, we execute the following steps: 



P Step 1: Calculate do-not-order value function J NP i ðIÞ. Since ‘2CðiÞ J ‘ ðI  di Þ has been calculated and stored in h(i, I  di), according to formuNP lation (8), J i ðIÞ can be obtained in OðnTÞ time since the number of breakpoints is bounded by OðnTÞ. Note here that the breakpoint I = di for J NP i ðIÞ is in our predetermined breakpoint list and is added into the breakpoint list for Ji(I). Step 2: Calculate order-at-capacity value function J Pq i ðIÞ. As in Step 1, according to (9), this step can be completed in OðnTÞ time since the number of breakpoints is bounded by OðnTÞ. Note here that the breakpoint I = di  q for J Pq i ðIÞ is also in our predetermined breakpoint list and is added into the breakpoint list for Ji(I). Step 3: Calculate and store

^ ¼ fi þ ci ðdik  mqÞ ^ þ Hi ðdik  mq ^  di Þ þ /0i ðk; mÞ

X

^  di Þ and J ‘ ðdik  mq

‘2CðiÞ

^ ¼ fi þ ci ðuk þ dik  mqÞ ^ þ Hi ðuk þ dik  mq ^  di Þ þ /ui ðk; mÞ

X

^  di Þ J ‘ ðuk þ dik  mq

‘2CðiÞ

^ such that 0 6 m ^ 6 sðkÞ  sðiÞ. Note here, the breakpoints I ¼ dik  mq ^  di ¼ for each combination of k 2 VðiÞ and an integer m ^  d1‘ and I ¼ uk þ dik  mq ^  di ¼ uk þ d1k  mq ^  d1‘ are in the set of breakpoints for h(i, I  di). For each combination d1k  mq P ^ such that 0 6 m ^ 6 sðkÞ  sðiÞ, the values of the functions ^ of node k 2 V and an integer m ‘2CðiÞ J ‘ ðdik  mq  di Þ and P ^ ‘2CðiÞ J ‘ ðuk þ dik  mq  di Þ can be obtained by binary search in Oðlog nÞ time, since the number of breakpoints is bounded by ^ Therefore, this entire step can be OðnTÞ and T 6 n. It can also be observed that there are OðnTÞ possible combinations of k and m. completed in OðnT log nÞ time. Step 4: Calculate and store J Pi ðIÞ. Since J Pi ðIÞ is piecewise linear, right continuous, and the slope for each piece is fixed (i.e., ci) for each ^ 0Þ (i.e., this breakpoint is within the interval ½dik  ðm ^ þ 1Þq; dik  mqÞ) ^ breakpoint, the qualified combination ðk; m; or the qualified ^ uÞ that has the smallest /i value (i.e., /0i ðk; mÞ ^ or /ui ðk; mÞ) ^ leads to the corresponding J Pi ðIÞ value for this breakcombination ðk; m; point. To obtain J Pi ðIÞ for each breakpoint, we can start from the breakpoint with the largest inventory value backwards to the breakpoint with the smallest inventory value. Corresponding to each breakpoint, we maintain and update a list of qualified com^ 0Þ and ðk; m; ^ uÞ according to a non-decreasing sequence of their /0i ðk; mÞ ^ and /ui ðk; mÞ ^ values. Denote the list as X binations ðk; m;  ^  ; 0Þ or ðk ; m ^  ; uÞ represents the combination serving as the head of the list. For instance, as shown in Fig. 5, the list of and let ðk ; m qualified combinations for the order value function J Pi ðIÞ at inventory level I ¼ dik1 is X = {(k1, 0, u), (k2, 0, 0), (k3, 0, 0)}. The head of ^  ¼ 0. Since J Pi ðIÞ is right continuous, it can be easily verified that, corresponding to the the list is (k1, 0, u). That is, k* = k1 and m breakpoint with the largest inventory value, there exists no combination in the form as shown in (10) such that this breakpoint  ^  Þ ¼ /ui ðk ; m ^  Þ ¼ þ1: Starting from this breakis within the corresponding interval. Therefore, we let X = ; and initialize /0i ðk ; m point backwards, for each breakpoint I = dik  mq orI = uk + dik  mq for some node k 2 V o and s(i) 6 m 6 s(k)  s(i), we perform the following sub-steps: ^  þ 1Þq (or (1) If X – ;, then check if the current breakpoint value I = dik  mq (or I = uk + dik  mq) is less than dik  ðm ^  þ 1Þq). If so, then the current head is out of range since it reaches the capacity, and the next one in the active uk þ dik  ðm list X will serve as the head. Note here that duplicated breakpoints are considered separately, and width for each piece of

Y. Guan, T. Liu / European Journal of Operational Research 207 (2010) 1398–1409

1407

Fig. 5. The order value function JPi ðIÞ corresponding to I ¼ dik1 .

the order value function line is the same (i.e., q). Thus, the next combination in the active list must be within the range and can serve as the head if it exists. This step takes Oð1Þ time.  ^  Þ  ci I or J Pi ðIÞ ¼ /ui ðk ; m ^  Þ  ci I. This step takes Oð1Þ (2) Calculate and store the corresponding order function value J Pi ðIÞ ¼ /0i ðk ; m time. (3) If g(i, k) = 1 and 0 6 m 6 s(k)  s(i), insert the corresponding combination into the active combination list by binary search to ^ maintain the non-decreasing order of /i(.) values in the list. For instance, if the breakpoint is in the form I ¼ dik  mq, ^ 6 sðkÞ  sðiÞ, then insert the combination ðk; m; ^ 0Þ; otherwise, insert the combination ðk; m; ^ uÞ. We also set this combi06m nation as the tail of the list, because it can be observed that the pairs in the list after this pair are not needed to be considered. If ^ or /ui ðk; mÞ ^ is smaller than the value for the current head of the list, then the combination ðk; m; ^ 0Þ or ðk; m; ^ uÞ serves as /0i ðk; mÞ both the head and the tail of the active list, and it becomes the only combination in the active list. This step takes Oðlog nÞ time. Since there are at most OðnTÞ breakpoints, this procedure can be completed in OðnT log nÞ time. That is, the breakpoints I ¼ d1k  d1i  mqand I ¼ uk þ d1k  d1i  mqfor J Pi ðIÞ with  s(i) 6 m 6 s(k)  s(i) for all k 2 V o can be calculated and stored in OðnTÞ time. Thus, this step takes OðnT log nÞ time. NP Step 5: Calculate and store Ji(I). For each breakpoint, take the minimum of J Pi ðIÞ; J Pq i ðIÞ and J i ðIÞ. This step can be completed in OðnTÞ time since the number of breakpoints is bounded by OðnTÞ. Step 6: Update h(i, I). Corresponding to each inventory breakpoint I ¼ d1k  d1i  mq and I ¼ uk þ d1k  d1i  mq for all k 2 V o and  s(i) 6 m 6 s(k)  s(i), increase the value function h(i, I) by Ji(I). This step can be completed in OðnTÞ time since the number of breakpoints is bounded by OðnTÞ. Therefore, the total amount of work required at each node is bounded by OðnT log nÞ. Since the above operations are required for each node i 2 V, the following conclusion holds. Theorem 2. The general stochastic inventory control problem with inventory bounds and constant order capacities can be solved in Oðn2 T log nÞ time, regardless of the structure of the scenario tree. It can be observed that the space complexity for the algorithm is Oðn2 Þ. 5. Conclusions and future research In this paper, we studied two models to solve general stochastic inventory control problems with inventory bounds and order capacities. These models especially address the dependency among demands and other cost parameters in different time periods. We developed Oðn2 Þ and Oðn2 T log nÞalgorithms for the cases with only inventory bounds and with both inventory bounds and order capacities. When our algorithms are applied to their deterministic counterparts, they have the same or better performance, as compared to the algorithms developed by Liu (2008), Atamtürk and Küçükyavuz (2008), and Wolsey (2006). In our stochastic programming models, replenishment is made after the demand for the same time period is realized. Our model can be extended to address the case in which the replenishment is done before demand realization at each time period. The only change is that we use Y i to represent the order quantity at period s(i) before demand in period s(i) is realized and update Eq. (1) to be Ii þ Y i ¼ di þ Ii . We can obtain a similar dynamic programming framework and the same computational complexity. Our model can also be extended to address the cases that include inventory setup costs. For this case, the inventory value function is not continuous. The value function for the (s, S) policy analysis is not continuous K-convex function, and therefore it is not clear if (s, S) policy is still optimal for the demand-independent cases. However, our approach can handle this case. Properties derived in this paper will still hold, and the overall computational complexity will still be the same. That is, there exists Oðn2 Þ time and Oðn2 T log nÞ algorithms for the problems with inventory setup costs and bounds, for the uncapacitated and constant-capacitated cases respectively. In future research, we will conduct numerical experiments to perform sensitivity analysis for inventory bounds and order capacities. In the experiments, we will test reasonable and extreme values for these parameters. Meanwhile, we will also combine the sensitivity analysis with pricing for capacity expansions to determine the optimal order-capacity and inventory-bound levels. Acknowledgments The authors would like to thank the area editor and the two anonymous referees for their helpful suggestions and comments on improving the quality of this paper. This research was partially supported by the U.S. National Science Foundation under Award CMMI-0700868

1408

Y. Guan, T. Liu / European Journal of Operational Research 207 (2010) 1398–1409

and under CAREER Award CMMI-0748204, and by the Center for Engineering Logistics and Distribution (CELDi), a National Science Foundation sponsored Industry/University Cooperative Research Center. Appendix A Proof of Propositions 1 and 3. First, we observe that Proposition 1 can be considered as a special case of Proposition 3 in which q is a very ^ ¼ 0 according to Proposition 3 if big number. This lies in the fact that, in the optimal solution, we have Ii > d1i , which implies m q > d1n þ maxi2V ui . Thus, Proposition 1 holds as long as Proposition 3 holds. In the following, we show that Proposition 3 holds. Let Z* = d(Y*) and X represent the set of points in the tree such that the replenishment at this node is positive and less than capacity q. That is, X ¼ fi 2 V : 0 < Y i < qg. Thus, we only need to prove the claim that ^ or Ii þ Y i ¼ dik þ uk  mq ^ for some k 2 VðiÞ and an integer m ^ such that 0 6 m ^ 6 sðkÞ  sðiÞ. We prove this for each i 2 X, Ii þ Y i ¼ dik  mq claim by contradiction. For example, suppose (Y*, Z*, I*) that meets neither of conditions in Proposition 3 is indeed optimal. We will show ^ ; Z ^  ; ^I Þ, such that one of them will have a better objective value. that we can generate two alternative solutions ðY  ; Z  ; I Þ and ðY Now assume the claim is not correct. That is, for any optimal solution (Y*, Z*, I*), there exists a subset X0 # X such that for each node ^ nor Ii þ Y i ¼ dik þ uk  mq ^ holds for any pair of k 2 VðiÞ and an integer m ^ such that i 2 X0, we have neither Ii þ Y i ¼ dik  mq ^ 0 6 m 6 sðkÞ  sðiÞ. Let u(i) be those descendants of node i that are the first nodes (except node i) with positive but less than capacity order quantities along a path from node i to a leaf node (note here that u(i) can be empty). Correspondingly, we use notation L1 ðiÞ to represent the set of leaf nodes in VðiÞ such that there is a node (besides node i) with positive and less than capacity order quantity along the path from node i to each leaf node in this set, and notation L2 ðiÞ to represent the set of leaf nodes in VðiÞ such that there is not any node (except node i) with positive and less than capacity order quantity along the path from node i to any leaf node in this set. Then L1 ðiÞ [ L2 ðiÞ ¼ L \ VðiÞ. Thus

n

o

uðiÞ ¼ [‘2L1 ðiÞ argmin sðjÞ : j 2 Pð‘Þ n PðiÞ and 0 < Y j < q : Correspondingly, let U(i) be the set of nodes in the paths from node i to each node in u(i) plus nodes in the paths from node i to each leaf  node in L2 ðiÞ; that is, UðiÞ ¼ [k2uðiÞ[L2 ðiÞ PðkÞ n Pði Þ. Let K(i) be the set of nodes in U(i) whose order amounts are equal to their capacities. That is, K(i) = {j 2 U(i):Yj = q}. Define the cost corresponding to the nodes in V n UðiÞ as

nðY  ; Z  ; I Þ ¼

 X   cj Y j þ fj Z j þ Hj ðIj Þ : j2VnUðiÞ

Then, the objective value for the given optimal solution (Y*, Z*, I*) is

  X   X  Hj Ij þ ðcj q þ fj Þ þ cj Y j þ fj Z j þ Hj Ij :

X

  fðY  ; Z  ; I Þ ¼ nðY  ; Z  ; I Þ þ ci Y i þ fi Z i þ Hi Ii þ

j2UðiÞnðfig[uðiÞÞ

j2KðiÞ

j2uðiÞ

^ ; Z ^  ; ^I Þ such that Now consider two alternative solutions ðY  ; Z  ; I Þ and ðY  Y i ¼ Y i  i , Ij ¼ Ij  i for each j 2 U(i)nu(i) and Y j ¼ Y j þ i for each j 2 u(i), ^  ¼ Y  þ i , ^I ¼ I þ i for each j 2 U(i) nu(i) and Y ^  ¼ Y   i for each j 2 u(i), Y i j j i j j and other components are the same as (Y*, Z*, I*). It is easy to observe that both solutions are feasible for SCB–LS. Note here that Ij – 0 for each j 2 U(i)n u(i) according to our assumption, and the problem parameters are finite. Then, there must exist a very small positive number i such that Ij  i > 0 if Ij > 0, and Ij þ i < 0 if Ij < 0, for each j 2 U(i)nu(i). The following equation holds:

  fðY  ; Z  ; I Þ ¼ nðY  ; Z  ; I Þ þ ci Y i  i þ fi Z i þ Hi ðIi  i Þ þ 

X



¼ fðY ; Z ; IÞ  ci i þ

jj ði Þ þ

j2UðiÞnuðiÞ þ

where jj ði Þ ¼ hj i if Ij > 0 and Similarly, we have

X

j2UðiÞnðfig[uðiÞÞ





  X    X   Hj Ij  i þ ðcj q þ fj Þ þ cj Y j þ i þ fj Z j þ Hj Ij j2KðiÞ

j2uðiÞ

cj i ;

j2uðiÞ

jj ði Þ ¼ hj i if Ij < 0.

  fðY^  ; Z^  ; ^I Þ ¼ nðY  ; Z  ; I Þ þ ci Y i þ i þ fi Z i þ Hi ðIi þ i Þ þ 

X

¼ fðY ; Z ; I Þ þ ci i 

X j2UðiÞnuðiÞ

jj ði Þ 

X

X j2UðiÞnðfig[uðiÞÞ

  X    X   Hj Ij þ i þ ðcj q þ fj Þ þ cj Y j  i þ fj Z j þ Hj Ij j2KðiÞ

j2uðiÞ

c j i :

j2uðiÞ

P P ^ ; Z ^  ; ^I Þg < fðY  ; Z  ; I Þ, contradicting the assumption that (Y*, Z*, I*) is an If ci i  j2UðiÞnuðiÞ jj ði Þ  j2uðiÞ cj i – 0, then minffðY  ; Z  ; I Þ; fðY P P ^Þ ¼ optimal solution. Otherwise, if ci i  j2UðiÞnuðiÞ jj ði Þ  j2uðiÞ cj i ¼ 0, then i can be increased such that Y i ðor Y i     ^ ^ ^ ^ 0; q; dik  mq  Ii ; or dik þ uk  mq  Ii for some node k 2 VðiÞ and 0 6 m 6 sðkÞ  sðiÞ, or Y j ðor Y j Þ ¼ 0 or q for some j 2 u(i). For any case,  0 ^  Ii for some node k 2 U(i) we have the number of elements in X reduced by at least one. Without loss of generality, assume Y i ¼ dik  mq ^ 6 sðkÞ  sðiÞ after increase the value of i. Then ðY  ; Z  ; I Þ is also an optimal solution. and 0 6 m

Y. Guan, T. Liu / European Journal of Operational Research 207 (2010) 1398–1409

1409

A similar argument can be applied to ðY  ; Z  ; I Þ if there exists a node j 2 V such that 0 < Y j < q and Y j – ^  Ij or djk þ uk  mq ^  Ij for any combination of k 2 VðjÞ and 0 6 m ^ 6 sðkÞ  sðiÞ and consecutively the following ones. We will djk  mq ^  Ii or Y i ¼ dik þ uk  mq ^  Ii for some node k 2 VðiÞ and eventually either obtain a solution (Y,Z,I) such that Y i ¼ dik  mq ^ 6 sðkÞ  sðiÞ if 0 < Yi < q corresponding to each node i 2 V, or find a solution with a smaller objective function value, contradicting 06m the original assumption. Therefore, the original conclusion holds. h References Aggarwal, A., Park, J.K., 1993. Improved algorithms for economic lot size problems. Operations Research 41, 549–571. Ahmed, S., King, A.J., Parija, G., 2003. A multi-stage stochastic integer programming approach for capacity expansion under uncertainty. Journal of Global Optimization 26, 3– 24. Ahuja, R.K., Hochbaum, D.S., 2008. Solving linear cost dynamic lot-sizing problems in Oðn log nÞ time. Operations Research 56, 255–261. Atamtürk, A., Küçükyavuz, S., 2008. An Oðn2 Þ algorithm for lot sizing with inventory bounds and fixed costs. Operations Research Letters 36, 297–299. Chan, G.H., Song, Y., 2003. A dynamic analysis of the single-item periodic stochastic inventory system with order capacity. European Journal of Operational Research 146, 529–542. Chen, S., Lambrecht, M., 1996. X–Y band and modified (s, S) policy. Operations Research 44, 1013–1019. Chung, C.S., Flynn, J., Lin, C.H.M., 1994. An effective algorithm for the capacitated single item lot size problem. European Journal of Operational Research 75, 427–440. Di Summa, M., Wolsey, L.A., 2006. Lot-sizing on a tree. Operations Research Letters 105, 55–84. Federgruen, A., Tzur, M., 1991. A simple forward algorithm to solve general dynamic lot sizing models with n periods in Oðn log nÞ or OðnÞ time. Management Science 37, 909– 925. Federgruen, A., Zheng, Y.S., 1992. An efficient algorithm for computing an optimal (r, Q) policy in continuous review stochastic inventory systems. Operations Research 40, 808–813. Florian, M., Klein, M., 1971. Deterministic production planning with concave costs and capacity constraints. Management Science 18, 12–20. Guan, Y., Ahmed, S., Miller, A.J., Nemhauser, G.L., 2006a. On formulations of the stochastic uncapacitated lot-sizing problem. Operations Research Letters 34, 241–250. Guan, Y., Ahmed, S., Nemhauser, G.L., Miller, A.J., 2006b. A branch-and-cut algorithm for the stochastic uncapacitated lot-sizing problem. Mathematical Programming 105, 55– 84. Guan, Y., Miller, A.J., 2008. Polynomial-time algorithms for stochastic uncapacitated lot-sizing problems. Operations Research 56, 1172–1183. Gutiérrez, J., Sedeñ-Noda, A., Colebrook, M., Sicilia, J., 2003. A new characterization for the dynamic lot size problem with bounded inventory. Computers & Operations Research 30, 383–395. Gutiérrez, J., Sedeñ-Noda, A., Colebrook, M., Sicilia, J., 2007. A polynomial algorithm for the production/ordering planning problem with limited storage. Computers & Operations Research 34, 934–937. Halman, N., Klabjan, D., Mostagir, M., Orlin, J., Simchi-Levi, D., 2009. A fully polynomial time approximation scheme for single-item stochastic lot-sizing problems with discrete demand. Mathematics of Operations Research 34, 674–685. Huang, K., Ahmed, S., 2009. The value of multi-stage stochastic programming in capacity planning under uncertainty. Operations Research 57, 893–904. Huang, K., Küçükyavuz, S., 2008. On stochastic lot-sizing problems with random lead times. Operations Research Letters 36, 303–308. Iglehart, D.L., 1963a. Dynamic programming and stationary analysis in inventory problems. In: Scarf, H., Guilford, D., Shelly, M. (Eds.), Multi-Stage Inventory Models and Techniques. Stanford University, Stanford, CA, pp. 1–31. Iglehart, D.L., 1963b. Optimality of (s, S) policies in the infinite horizon dynamic inventory problem. Management Science 9, 259–267. Kirca, Ö., 1990. An efficient algorithm for the capacitated single item dynamic lot size problem. European Journal of Operational Research 45, 15–24. Kunnumkal, S., Topaloglu, H., 2008. Using stochastic approximation methods to compute optimal base-stock levels in inventory control problems. Operations Research 56, 646–664. Lee, C.Y., Cetinkaya, S., Wagelmans, A., 2001. A dynamic lot-sizing model with demand time windows. Management Science 47, 1384–1395. Liu, T., 2008. Economic lot sizing problem with inventory bounds. European Journal of Operational Research 185, 204–215. Love, S.F., 1973. Bounded production and inventory models with piecewise concave costs. Management Science 20, 313–318. Richter, K., Sombrutzki, M., 2000. Remanufacturing planning for the reverse Wagner/Whitin models. European Journal of Operational Research 121, 304–315. Richter, K., Weber, J., 2001. The reverse Wagner/Whitin model with variable manufacturing and remanufacturing cost. International Journal of Production Economics 71, 447– 456. Ruszczyn´ski, A., Shapiro, A. (Eds.), 2003. Stochastic Programming. Handbooks in Operations Research and Management Science. Elsevier Science B.V.. vol. 10. Sargut, F.Z., Romeijn, H.E., 2007. Lot-sizing with non-stationary cumulative capacities. Operations Research Letters 35, 549–557. Scarf, H., 1960. The optimality of (s, S) policies for the dynamic inventory problem. Proceedings of the First Stanford Symposium of Mathematical Methods in the Social Sciences. Stanford University Press, Stanford, CA. Toczylowski, E., 1995. An OðT 2 Þ algorithm for the lot-sizing problem with limited inventory levels. IEEE Symposium on Emerging Technologies & Factory Automation 3, 78– 85. van den Heuvel, W., Wagelmans, A.P.M., 2008. Four equivalent lot-sizing models. Operations Research Letters 36, 465–470. van Hoesel, C.P.M., Wagelmans, A., 1996. An OðT 3 Þ algorithm for the economic lot-sizing problem with constant capacities. Management Science 42, 142–150. Veinott, A.F., Wagner, H.M., 1965. Computing optimal (s, S) inventory policies. Management Science 11, 525–552. Wagelmans, A., van Hoesel, A., Kolen, A., 1992. Economic lot sizing: An Oðn log nÞ algorithm that runs in linear time in the Wagner–Whitin case. Operations Research 40, 145– 156. Wagner, H.M., 1960. A postscript to dynamic problems of the firm. Naval Research Logistics Quarterly 7, 7–12. Wagner, H.M., Whitin, T.M., 1958. Dynamic version of the economic lot size model. Management Science 5, 89–96. Wolsey, L.A., 2006. Lot-sizing with production and delivery time windows. Mathematical Programming 107, 471–489. Zangwill, W.I., 1968. Minimum concave cost flows in certain networks. Management Science 14, 429–450. Zhao, X., Fan, F., Liu, X., Xie, J., 2007. Storage-space capacitated inventory system with (r, Q) policies. Operations Research 55, 854–865. Zheng, Y.S., 1991. A simple proof for optimality of (s, S) policies in infinite-horizon inventory systems. Journal of Applied Probability 28, 802–810. Zheng, Y.S., Federgruen, A., 1991. Finding optimal (s, S) policies is about as simple as evaluating a single policy. Operations Research 39, 654–665.