New model and heuristics for safety stock placement in general acyclic supply chain networks

Computers & Operations Research 39 (2012) 1333–1344

Contents lists available at ScienceDirect

Computers & Operations Research journal homepage: www.elsevier.com/locate/caor

New model and heuristics for safety stock placement in general acyclic supply chain networks Haitao Li a,n, Dali Jiang b a b

Department of Logistics and Operations Management, College of Business Administration, University of Missouri—St. Louis, One University Blvd, St. Louis, MO 63121, USA Institute of Modern logistics, Logistical Engineering University, Chongqing 401311, P.R.China

a r t i c l e i n f o

abstract

Available online 9 August 2011

We model the safety stock placement problem in general acyclic supply chain networks as a project scheduling problem, for which the constraint programming (CP) techniques are both effective and efficient in finding high quality solutions. We further integrate CP with a genetic algorithm (GA), which improves the CP solution quality significantly. The performance of our hybrid CP–GA algorithm is evaluated on randomly generated test instances. CP–GA is able to find optimal solutions to small problems in fractions of a second, and near optimal solutions of about 5% optimality gap to medium size problems in several minutes on average. & 2011 Elsevier Ltd. All rights reserved.

Keywords: Safety stock placement Inventory positioning General acyclic network Project scheduling Constraint programming (CP) Genetic algorithm (GA)

1. Introduction Supply chains are complex networks integrating suppliers/vendors, manufacturers, distributors and retailers to enable physical entities (raw materials, parts, components, semi-finished and finished products) to be produced and distributed at the right quantities, to the right locations and at the right time. Demand uncertainty requires safety stock to be kept through the supply chain network to meet customers’ order and delivery needs. The problem of determining the location and level of safety stock to keep in the supply chain is known as the inventory positioning or safety stock placement problem (cf. [1,2]). Unlike the traditional strategic supply chain network design problems, the safety stock placement problem addresses the tactical level design for an existing physical supply chain, by optimally balancing system-wide inventory cost and service level. Optimizing the supply chain system-wide safety stock is crucial in today’s complex, competitive and volatile market environment. Since a modern supply chain rarely involves a single isolated firm, the focus here is not limited to coordinating internal materials and inventory management in the traditional material requirement planning (MRP) setting, but also involves external suppliers, vendors, logistics providers and manufacturers ([3], p.1). With the rapid growth of business-to-business (B2B) e-commerce facilitated by the electronic data interchange (EDI) technology and availability of point of sales (POS) data, the ability to redefine inter-firm relationships and coordinate activities among supply chain players becomes a key for supply chain success. Such external coordination calls for

n

Corresponding author. Tel.: þ1 314 516 5890. E-mail address: [email protected] (H. Li).

0305-0548/$ - see front matter & 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.cor.2011.08.001

collaboration among supply chain players, so that decisions can be made in a centralized way. Sometimes the centralized decisionmaking paradigm is natural. For instance, the logistics of the U.S. Army’s spare parts [4] and the U.S. Navy’s depot level repairable (DLR) line items [5] require inventory positioning decisions to be made in a centralized fashion to support military missions. In the business environment, however, external collaboration and coordination need a well-defined incentive and often comes at a price. Examples of real world practices include the continuous replenishment process (CRP), vendor managed inventory (VMI) [6] and supply chain partnership [7]. In the safety stock placement problem addressed in this paper, we assume that supply chain players fully collaborate such that centralized safety stock placement decision can be made. The issues of creating incentive, negotiating and managing such collaborative relationship go beyond the scope of this paper. Safety stock placement is also related to determining the push– pull boundary or implementing the delay differentiation (postponement) strategy in multi-echelon supply chains [8]. An isolated manufacturer often operates either make-to-stock (MTS) or maketo-order (MTO). In an MTS process, components and finished goods are manufactured based on long-term forecasting and planning, and held as inventories to meet the anticipated future demand. MTS is sometimes referred to as the push strategy. In an MTO process, manufacturing is only triggered by a realized order. MTO is often referred to as the pull strategy. In the setting of an integrated supply chain, firms do not have to implement a pure push or pull strategy. By properly coordinating the activities of supply chain players, firms may benefit from storing inventories at certain intermediate stages in the network, giving rise to the hybrid MTS/MTO or push–pull strategy. We refer to [9] for an application of configuring a laptop PC supply chain.

1334

H. Li, D. Jiang / Computers & Operations Research 39 (2012) 1333–1344

We consider a supply chain network consisting of functions, e.g., parts, components or processes, and the demand dependency among them as characterized by a bill-of-materials (BOM). We assume that the network is acyclic, i.e. no directed cycle (‘‘feedback’’) exists. Each sink node faces a Gaussian normal-distributed customer demand to be satisfied. Such external demand will be propagated backward along the supply chain through the demand dependencies in BOM. Each upstream function will quote a guaranteed delivery time, called outbound service time, to satisfy the demand of its immediately downstream function(s). Since the promised delivery time for the finished goods is often smaller than the total supply chain cycle time, certain amount of safety stock needs to be kept at each function to maintain certain service level. We assume that each function operates under a periodic review policy. Accordingly, the level of safety stock for the corresponding function is determined by its net replenishment time. In addition, there might be capacity limit for storage of safety stock that cannot be exceeded at each function. The safety stock placement problem addressed in this paper seeks to find optimal inbound and outbound service times of each function in the supply chain to minimize the system-wide safety stock cost. Complexity of the safety stock placement problem is directly driven by structure of the supply chain. Fig. 1 illustrates three different supply chain structures. Among them, the serial structure in Fig. 1(a) considers only sequential dependencies among supply chain functions, resulting a linear/serial structure. Fig. 1(b) illustrates the spanning tree structure. Clearly, the spanning tree structure is able to model more complex interactions among supply chain functions than the serial structure. For instance, for every unit of component/part assembled at manufacturer C, 2 units from supplier A and 1 unit from supplier B are required. At the same time, manufacturer C cannot start its inbound service until its demand for all its immediate predecessors, A and B, has been satisfied. Through such demand- and timedependencies, the decision at one location will be propagated to other functions in the chain. For instance, by promising a prompt delivery lead time to its customers, a car dealer will not only potentially carry more safety stock himself, but would also like to require an early receipt of the order from his warehouse, which creates pressure for the warehouse to carry more inventory as

A

2 1

C

B

well. Two special cases of the spanning tree structure include the assembly (convergent) network where each node has at most one outgoing arc, and the distribution (divergent) network where each node has at most one incoming arc. The structure of real-life supply chains often goes beyond the complexity of a spanning tree, i.e. the number of arcs in the network will often exceed the number of functions minus one. These supply chains are characterized by a general acyclic network as depicted in Fig. 1(c). More general cases include those where arcs may exist between an arbitrary pair of supply chain functions as shown in Fig. 1(d). It is well-known that the safety stock placement problem in networks with special structure such as the serial and spanning tree structures can be efficiently solved by dynamic programming (cf. [10,11]). The one with the general acyclic structure involves minimizing a concave function, and it is well-known that the general problem of minimizing a concave function is NP-hard [12]. Networks with general structure have been a focal subject in the research field of project scheduling. Since the development of the well-known critical path method (CPM; [13]) and program evaluation and review technique (PERT; [14]) in the early sixties, the project scheduling field has devised numerous efficient solution methods for scheduling activities in large-size general project networks (cf. [15,16]). This has motivated us to resort to the modeling and solution techniques in project scheduling to address the safety stock placement problem in general acyclic supply chain networks. In our modeling framework, an acyclic supply chain network is transformed into a project scheduling network, with each supply chain function being modeled as a project activity. We show that time-related dependencies among inbound and outbound service times can be modeled as temporal relationships among project activities. Then the problem of finding optimal inbound and outbound service times becomes equivalent to the problem of scheduling project activities’ starting times. Such a transformation enables us to borrow a large set of well-developed algorithms in the project scheduling field (cf. [15]) to tackle the safety stock placement problem. In this paper, we devise a constraint programming (CP; [17]) based algorithm to solve the transformed problem. We embed heuristic branching rules exploiting the scheduling characteristics of the problem to enhance CP’s performance. One advantage of the CP algorithm is its ability to generate high-quality feasible solutions fast. When used as initial solutions, these feasible solutions obtained by CP can be improved by a number of local search [18] and metaheuristic methodologies [19]. We design a hybrid framework that integrates CP and genetic algorithm (GA), called CP–GA, to enhance solution quality. Computational results show that our CP–GA algorithm is able to find optimal solutions to small instances of the problems within fractions of a second, and obtain near optimal solutions to medium size problems in reasonable computational time. The remainder of our paper is organized as follows. Section 2 reviews the related literature. We present the project-schedulingbased modeling framework and CP model in Section 3. Section 4 describes the CP solution procedures including constraint propagation and various search strategies devised in this study. The hybrid CP–GA algorithm is presented in Section 5. The computational results are provided in Section 6. Section 7 draws conclusions and suggests future research directions.

2. Related literature

Fig. 1. Supply chain networks with different structures: (a) serial network, (b) spanning tree network, (c) general acyclic network, (d) general acyclic network with arbitrary arcs.

The safety stock problem with simple supply chain structure, such as the serial and spanning tree structures, has been extensively studied. For a serial supply chain, an ‘‘all or nothing’’ strategy is suggested by [20]. The effect of stochastic lead time on inventory positioning is studied in [21]. Inventory positioning in multiple

H. Li, D. Jiang / Computers & Operations Research 39 (2012) 1333–1344

product supply chains with serial structure has been studied in [22]. A serial structure network can be extended to an assembly network where each node has at most one outgoing arc, or a distribution network where each node has at most one incoming arc. Both are special cases of the spanning tree structure [23]. Efficient dynamic programming algorithms have been proposed for serial, divergent (distribution), or convergent (assembly) networks ([10,24]). A pseudo-polynomial dynamic programming algorithm was proposed by [11] for a spanning tree supply chain that consists of a mixture of assembly and distribution sub-networks. Recent research efforts have focused on the more challenging and realistic general acyclic supply chain networks. It is shown that the inventory positioning problem in a general supply chain network is NP-hard [25]. Exact methods such as enumeration and branch-and-bound [25] are often not able to handle real-world problems due to prohibitive amount of computational time for large instances. A piecewise linear approximation algorithm was devised by [1] to iteratively generate linear pieces to approximate the concave objective function. Some redundant constraints are formulated to generate strong flow cover cuts. Two heuristics based on linear approximation and truncated piecewise MIP formulation (with only two linear pieces) were developed in [2], which are shown to be efficient for obtaining near optimal solutions: the iterative LP method can efficiently solve problems with size as large as 8000 nodes and 32,000 arcs. Project scheduling based approaches have been proposed for various manufacturing operations such as order selection, manufacturing planning and scheduling in the make-to-order environment [26]. A project scheduling network has been used to model a generalized supply chain to minimize the supply chain cycle time [27]. By transforming a supply chain network into a project scheduling network, a large set of well-developed project scheduling algorithms can be applied. Effectively and efficiently solving complex scheduling problems can often be achieved by constraint programming (CP) based methods. CP originated in the artificially intelligence (AI) area for solving constraint satisfaction problem (CSP). We refer to [17] for a general introduction to CP. With effective domain reduction (also called constraint propagation) algorithms [28], and efficient search procedures [29], CP based algorithms have been successfully applied for solving a variety of NP-hard scheduling problems (cf. [30]). For large size scheduling problems, however, CP alone is often not efficient due to the burden of searching the entire solutions space (although reduced through constraint propagation). There exist a number of integration schemes to take advantage of the complementary strengths of CP and other solution methods. CP has been used as a preprocessing procedure to reduce the size of project scheduling problems with resource constraints [31]. CP can also be integrated with mixed-integer linear programming (MILP) in the frameworks of Benders decomposition [32], Lagrangian relaxation [33] or column generation [34]. Schemes for integrating CP and local search (LS) have also been proposed [35]. We refer to [36,37] and [38] for surveys on CP-based hybrid algorithms.

3. Problem description We now formally describe the safety stock placement problem in general acyclic supply chain networks, then present its CP formulation. 3.1. Problem description Consider a general acyclic supply chain network G(N,A), with the node set N representing the set of supply chain functions, and

1335

P1 C1 SA1

P2

D1

A1 C2 P3 SA2

D2

P4 C3

P5

A2

SA3

D3

C4

Procurement

Manufacturing

Subassembly

Assembly

Retailers

Note: P – Part, C – Component, SA – Subassembly, A – Assembly, D – Demand.

Fig. 2. An example of general acyclic supply chain network.

the arc set A denoting the set of demand dependencies and temporal relationships among functions. A function can be a part, component or a production/assembly process as characterized by the BOM associated with the product (or family of products). A typical example of general acyclic supply chain network is illustrated in Fig. 2, which consists of 5 stages (echelons) and 17 supply chain functions. It is a ‘‘general’’ network in which it includes the spanning tree network as a special case (i.e. the number of arcs exceeds the number of nodes minus one), and also allows that arbitrary arcs exist among supply chain functions. For instance, Part-1 purchased at the procurement stage is required for manufacturing Component-1 at the manufacturing stage, as well as assembling finished product A1 at the assembly stage. The network in Fig. 2 is ‘‘acyclic’’ as no directed cycle exists. An arc (i,j)AA specifies the demand dependencies between functions i and j. For instance, one unit of Component-1 may require 2 units of Part-1 and 1 unit of Part-2. An arc (i,j) also specifies the temporal (time-related) relationships between function i and j. We assume that the downstream function j cannot start its inbound service time until its upstream (predecessor) function i has been delivered (i’s outbound service has been fulfilled). Following [1] and [2], each function j is assumed to have a constant lead time Tj and a safety stock holding cost rate hj per item per time unit. There is a capacity limit of Cj time periods of safety stock for each function j. The set D of end products at the demand stage are assumed to face stochastic external demands bounded by a non-decreasing concave function, which is needed to establish a meaningful upper bound on demand over varying horizon for each end product. We refer to [11] for detailed justifications for this bounded demand assumption. When the external demand is normally distributed ðm, s2 Þ and i.i.d., the pffiffiffi demand bound can be expressed as tm þ ks t, where t is the net replenishment time, and k is the z-value associated certain pffiffiffi service level. Then the safety stock takes the form ks t and is pffiffiffi proportional to s t, which is commonly used in the literature (cf. [1–3]). In our example, the set D includes three different retailers

1336

H. Li, D. Jiang / Computers & Operations Research 39 (2012) 1333–1344

fD1,D2,D3g. Finished products have to be delivered to them within some customer-specified deadline dj on lead time. The safety stock placement problem is to find optimal inbound and outbound service times of each supply chain function, thus the corresponding optimal safety stock levels at each function, that minimize the system-wide inventory holding cost of the entire supply chain network. Following, we define decision variable Sj as the guaranteed outbound service time that each function j promises to its downstream functions, and SIj as the maximum inbound replenishment lead time for function j, i.e. SIj Zmax{Si9i:(i,j)AA}, which is also the upper bound on the difference between the guaranteed service delivery time Sj and the lead time Tj at j, i.e. SIj ZSj Tj. The mathematical programming formulation of the general capacitated safety stock placement problem can be written as X hj Fj ðSIj þ Tj Sj Þ ð1Þ min jAN

Supply Chain Functions

dummy activity real activity T net replenishment time

Function-1 T

Function-2

SI

S

SI

S

Time

Fig. 3. A Gantt Chart illustrating the project-scheduling-based modeling framework for the safety stock placement problem.

SIj Si Z0,8ði,jÞ A A

ð3Þ

ILOG CP Optimizer [40] to describe the CP formulation. a[j] and b[j] can be declared as below (with bold words being the keywords in OPL and CP Optimizer):

Sj rdj ,8j A D

ð4Þ

Activity a½j in NðTj Þ

ð7Þ

ð5Þ

Activity b½j in Nð0Þ

ð8Þ

s:t:

Sj SIj r Tj ,8j A N

SIj þTj Sj rCj ,

8j A N

Sj ,SIj A Z þ ,8j A N

ð2Þ

ð6Þ

The objective function (1) minimizes the total safety stock cost throughout the entire supply chain, where Fj(U) denotes the level of safety stock at j in terms of days of demand as a function of net replenishment lead time, and is assumed to concave and nondecreasing. Constraint (2) enforces that the net replenishment lead time SIj þTj Sj of each function j must be non-negative. Constraint (3) specifies the precedence relationship between a pair of functions (i,j)AA, i.e. the inbound service time of downstream function j must be no less than the outbound service time of upstream function i. Constraint (4) ensures finished products are delivered to customers within the deadline on lead time. Constraint (5) satisfies the capacity constraint, i.e. the net replenishment time of function j cannot exceed Cj. When the lead time is long and/or deadline is tight, but the available capacity is small, a problem instance might be infeasible. Following [9] and [11], we require that both Sj and SIj be integer, which is in line with the fundamental periodic review policy assumed for each j. We refer to [9] for detailed justifications on the rationale of this integrality requirement. Note that the integrality constraints make our model more restrictive than the formulation in [1] and [2], although integer solutions can often be guaranteed as long as all the input data are integers, given the intrinsic pure network structure of the problem [1]. The safety stock problem described by (1)–(6) is a separable concave minimization problem, which is well-known to be NPhard [12]. Instead of relying on the above nonlinear programming formulation, our approach transforms the addressed safety stock placement problem into a project scheduling problem, and utilizes constraint programming (CP) based methods to solve it. The CP formulation is presented next.

The value in parenthesis specifies the duration of an activity, i.e. the lead time Tj for a real activity a[j], and 0 for a dummy activity b[j]. Such a modeling framework is depicted by the Gantt Chart in Fig. 3. Consider a pair of directly connected functions Function—1 and Function—2. Each function is associated with a real activity representing the actual execution of the function (represented by a rectangle in the chart) and a dummy activity denoting the promised outbound service time (represented by a triangle in the chart). The dummy activity of each function has to start before the end of the corresponding real activity to ensure a non-negative net replenishment time. The length of net replenishment time is represented by a grayed rectangle in the chart. For instance, the net replenishment time of Function—1 equals the difference between the end of the corresponding real activity and the start of the dummy activity, i.e. SI1 þT1  S1. At the same time, the downstream Function—2 cannot start before the start of the dummy activity of the upstream Function—1. Such a timedependency is represented by a solid arrow in the chart. The goal is to find feasible starting times of all real and dummy activities satisfying the temporal constraints, such that the total inventory cost expressed as a nonlinear function of these starting times is minimized. Based on the above modeling framework, we now present the CP formulation. 3.2.1. Non-negative net replenishment time To ensure non-negative net replenishment time for each supply chain function jAN, the associated dummy activity b[j] must start before the end of real activity a[j]: b½j:startsBeforeEndða½jÞ8j A N

ð9Þ

Constraint (9) is equivalent to Constraint (2). 3.2. Constraint programming formulation Our project-scheduling-based framework for the safety stock placement problem models the entire supply chain as a project. For each supply chain function jAN, we define a real activity a[j] representing the execution of function j, and a dummy activity b½j representing the outbound service of j. For illustration purpose, we use the CP constructs available in ILOG OPL Studio [39] and

3.2.2. Precedence relationships For every pair of arc (i,j)AA, the real activity a[j] associated with function j must not start until the start of the dummy activity b[i] associated with function i: a½j startsAfterStart ðb½iÞ8ði,jÞ A A Constraint (10) is equivalent to Constraint (3).

ð10Þ

H. Li, D. Jiang / Computers & Operations Research 39 (2012) 1333–1344

3.2.3. Guaranteed customer delivery lead time For each finished product jAD, the associated dummy activity b[j] must start before the customer-specified deadline dj on delivery lead time: b½j startsBefore ðd½jÞ

8j A D

Definition 1. A constraint is called binary constraint if it involves exactly two decision variables.

3.2.4. Capacity constraints In our CP formulation, the capacity constraint (5) becomes 8j A N

ð12Þ

where Tj  Cj is a minimum time-lag between the start of dummy activity b[j] and the start of real activity a[j] associated with supply chain function j. 3.2.5. Objective function Using CP constructs, the objective function (1) can be written as X min hj Fj ða½j:start þTj b½j:startÞ ð13Þ jAN

where a[j].start and b[j].start denote the starting time of real activity a[j] and dummy activity b[j], respectively.

4. CP algorithm This section describes the CP algorithms to solve the safety stock placement problem. A CP algorithm consists of two main solving techniques: constraint propagation and search. Constraint propagation is a problem reduction technique that transforms problems into equivalent problems, which are easier to solve. It modifies the domains of all variables in a constraint given the change of one of the variables in that constraint. We refer to [41] for theoretical background of constraint propagation and definitions of different consistency concepts. A search procedure is needed to explore the remaining solution space, as problem reduction through constraint propagation alone is often NP-hard. 4.1. Constraint propagation Two types of constraint propagation are implemented in our CP algorithm: I. the initial constraint propagation to reduce the solution domain before the search starts; II. constraint propagation during the search. We use temporal analysis in project scheduling to deduce the time windows of both real and dummy activities, i.e. the type I constraint propagation. Following [16], we solve the following two linear programs (LP): ES_LP ¼ f

min

X

Constraint propagation during the CP search, the type II constraint propagation, is based on the fundamental concept of consistency, and specifically, the arc-consistency (AC) for binary constraints in constraint satisfaction problems (CSP, [17,41]).

ð11Þ

Constraint (11) is equivalent to Constraint (4).

b½j startsAfterStart ða½j,Tj Cj Þ,

1337

ðSIj þ Sj Þ

jAN

s:t: Constraints ð2Þð6Þ g, and LS_LP ¼ X ðSIj þSj Þ f max jAN

s:t: Constraints ð2Þð6Þ g: ES_LP obtains the earliest start (ES) times of real and dummy activities, and LS_LP obtains the latest start (LS) times. The time window ½ESj ,LSj  of real/dummy activity j reduces the domain of decision variables.

It is important to note that our entire constraint system consists of mainly binary constraints, i.e. (9)–(12), for which AC algorithms will be especially effective [37]. We now define AC in the setting of our CP model. Definition 2. Let X’a denotes an assignment of value a to variable X. A binary constraint Cinvolving two start times SIx and SIy of functions x and y, respectively, is arc-consistent, if 8t A ½ESx ,LSx , (t 0 A ½ESy ,LSy  such that fSIx ’t,SIy ’t 0 g satisfies C, and 8t 0 A ½ESy ,LSy , (t A ½ESx ,LSx  such that fSIx ’t,SIy ’t 0 g satisfies C. In other words, constraint C is said to be arc-consistent, if for every value t in the time window of x, there always exists a value t0 in the time window of y satisfying C; at the same time, for every t0 in the time window of y, there exists a value t in the time window of x satisfying C. Define S as the set of all decision variables, i.e. S ¼ Sj A N [SIj A N , and W as the set of time windows (domains) of variables. Let C consist of all binary constraints of (9)–(12). The fundamental AC algorithm, known as AC-1 [41], for achieving arc-consistency of the safety stock placement problem is described below. Procedure: AC-1 ðS,W,CÞ Step 1. Construct constraint list: Q’fCðx,yÞ 9Cðx,yÞ A Cg; Step 2. Examine each constraint in the list: Repeat DomainReduced ’false; For each Cðx,yÞ A Q Do DomainReduced ’ Update_Domain (C(x,y), S,W,C) Until Not DomainReduced Return ðS,W,CÞ. Step 1 constructs a list Q of binary constraints, each of which involves two variables x and y, to be checked for consistency. Step 2 examines each constraint C(x,y) in the list Q, and removes all values, which do not satisfy C(x,y) by calling the sub-procedure Update_Domain. If the domain can be reduced, Step 2 is repeated; otherwise, the AC-1 procedure is terminated. The sub-procedure Update_Domain is described below: Sub-Procedure: Update_Domain (C(x,y), S,W,C) Step 1. Initialization: Deleted ’false; Step 2. Check domain of variables: For each t A W x Do If ) t 0 A W y such that fx’t,y’t 0 g satisfies C(x,y), Then W x ’W x ftg; Deleted ’true; Return Deleted. In Step 2 of the sub-procedure, each value t in the current domain W x of x is scanned. If for a specific t, there does not exist a value t0 in the current domain W y of y, such that the assignment fx’t,y’t 0 g is feasible for C(x,y), the value t is removed from W x . If the maximum number of elements in the domains is M, it is straightforward to show that AC-1 has a time complexity of OðM 3 U9S9U9C9Þ. Several improved AC algorithms are available in the CSP literature [41]. In our implementation, the AC-5 arc-consistency algorithm introduced in [42] is used. AC-5 improves over earlier versions of AC algorithm

1338

H. Li, D. Jiang / Computers & Operations Research 39 (2012) 1333–1344

by working with a queue that iteratively removes a constraint/value pair and performs consistency tests, and is significantly more efficient [43]. 4.2. Branching rules An important component in a CP algorithm is its efficiency to search the solution space. We implement two heuristic branching rules to control the order for selecting variables/values during the CP search. 4.2.1. First-fail rule The first-fail rule selects the most ‘‘difficult’’ variable first, i.e. the one with the smallest cardinality of domain [41]. It is generally believed to result in a smaller search tree, as some variables can be bounded early during the search. As a generic search strategy, however, the first-fail rule does not explore the features of the safety stock placement problem. 4.2.2. Maximal-regret rule We also implement a branching rule based on the idea of maximal-regret principle, which first chooses the variable most promising in improving solution quality. The term regret refers to the difference between what would have been the best decision in a scenario and what was the actual decision. The key of such a branching rule is that it allows one to exploit the structure of the addressed problem, which may lead to finding quality solutions early during the search. It has been shown that finding good feasible solutions early will often reduce the size of CP search tree [41]. The design of maximal-regret rule relies on the following properties due to the scheduling characteristics of the safety stock placement problem. Property 1. Given a fixed schedule for dummy activities, the objective function is regular (non-decreasing) w.r.t. the starting times of real activities. Proof. Let tj represent the net replenishment time of supply chain function jAN, i.e. tj ¼a[j].start þTj b[j].start. Since Fj(U) is a non-decreasing function of tj, the objective function f(U) is also non-decreasing. Thus tj r t0j implies that f ðtj Þ r f ðt0j Þ. Since both Tj and b[j].start are fixed, f is non-decreasing in a[j].start. & Property 2. Given a fixed schedule for real activities, the objective function is non-regular (non-increasing) w.r.t. the starting times of dummy activities. Proof. Similar to the proof of Property 1. & Property 1 and Property 2 imply that an early start time for a real activity, and a late start time for a dummy activity are likely to produce quality solutions. Such insights lead to the following composite maximal-regret branching rule: i) For real activities, we first choose a[j], for which the difference between the earliest possible start time and the next earliest start time is maximal. Then assign the earliest possible start time to a[j].start. ii) For dummy activities, we first choose b[j], for which the difference between the latest possible start time and the next latest start time is maximal. Then assign the latest possible start time to b[j].start. 5. Hybrid CP–GA algorithm Our computational experiences indicate that CP is quite effective for small safety stock placement instances. When the

problem size increases, it becomes expensive for CP alone to search the solution space. This has motivated us to design hybrid algorithms to enhance the quality of CP solutions. Our current research develops a scheme to integrate CP with genetic algorithm (GA). GA is a random search technique that mimics processes observed in natural evolution [44]. It combines survival of the fittest (or best) among solutions with a structured yet randomized exchange (cf. [45,46]). In our hybrid CP–GA algorithm, GA is used as a metaheuristic strategy to improve the CP solutions; on the other hand, CP provides high-quality solutions, which can construct the initial population or serve as elite individuals in the initial population for GA. In addition, CP also provides tighter bounds on the decisions variables, i.e. the inbound and outbound service times, which can significantly reduce the GA search space. The hybrid CP–GA framework is described next, followed by details about the GA components.

5.1. Hybrid CP–GA framework Our proposed hybrid CP–GA framework is depicted in Fig. 4. The initial population is composed of CP solutions and also by a set of randomly generated feasible solutions. CP solutions are obtained by running the CP algorithm described in Section 4 with certain limit on the CP search time. We include the best CP solution together with some intermediate feasible solutions found during the CP search. Randomly generated feasible solutions are also included to introduce diversity to the initial population. Standard penalty function is employed to penalize infeasible solutions. In addition to the initial solutions, CP also provides tight bounds of the decision variables, which can significantly reduce the solution space for GA. Reproduction, crossover and mutation are implemented in a parallel fashion to generate offspring. Crossover enables the algorithm to extract the best genes from different individuals and recombine them into potentially superior children. The crossover operator is applied to each selected pair of parent chromosomes with probability pcross. Mutation adds diversity to the population, and thereby increases the likelihood of generating individuals with better fitness values. The mutation operator is applied to each individual with probability pmut. Thus the entire new generation consists of npcross individuals from crossover, npmut individuals from mutation, and the elite individuals from reproduction. We set Nmax as the maximum number of generations to run. That is, the GA terminates if the number of generations exceeds Nmax .

5.2. Encoding scheme Each chromosome is encoded as a vector of positive integer values, which directly represents the decision variables Sj and SIj. Fig. 5 illustrates such encoding scheme, with each chromosome having a length of 2n, where n is the number of nodes in the supply chain network.

5.3. Genetic operators 5.3.1. Parent selection Parent selection is critical to guide GA toward promising regions in the solution space. In our implementation, we employ the roulette wheel selection method [47] to keep diversity of chromosomes.

H. Li, D. Jiang / Computers & Operations Research 39 (2012) 1333–1344

1339

Constrain Programming

Variable bounds

Best CP solution

Encoding scheme

Feasible CP solutions

Generate random feasible solutions

Initialize population T:=1

Population

Calculate obj. function value

Fitness score Get fitness Select inheritance

Crossover

Reproduction

Mutation

T:=T+1

No T>Nmax Yes Best solution found Fig. 4. Sketch of the hybrid CP–GA framework.

S

......

S

S

SI

SI

……

SI

Fig. 5. Representation of a chromosome.

5.3.2. Crossover operator The crossover operator generates new children by combining information contained in the chromosomes of the parents so that new chromosomes will have the best parts of the parents’ chromosomes. As the solution space is large, it is crucial to avoid spending excessively search effort evaluating infeasible solutions. This has motivated us to consider the arithmetic operator as described in [48] (p. 246), which is a linear (convex) combination of two feasible solutions. When a problem involves only linear constraints, the solution space is convex, so that an arithmetic operator always yields a feasible solution in the convex search space. It has proved to be effective for many linearly constrained problems [49]. Although our safety stock placement problem involves only linear constraints (2)–(5), the solution space is not convex due to the existence of integer constraint (6), thus the elementary arithmetic operator is not immediately applicable. We employ a modified arithmetic crossover operator to create offspring x0 from its parents x1 and x2 in the following way:

feasible to the constraint system (2)–(6), x0 is also feasible. A detailed proof of this property is provided in the Appendix. This modified arithmetic crossover operator enables us to reduce the number of infeasible solutions to be evaluated, enhancing the efficiency and effectiveness of GA significantly.

5.3.3. Mutation operator The role of mutation is to provide a small amount of randomness to explore new solution space and prevent premature convergence. In our GA, the mutation operator first randomly selects one gene of decision variables on a chromosome, and then replaces the selected gene value with a random integer within the solution bound. In order for the child through mutation to be feasible, the following steps are implemented. Step 1. Select a parent and let child ¼parent. Step 2. Set the mutation position: Position¼ l1  2n, where l1 is a real random number in [0.1], and 2n is the length of chromosome. Step 3. Set the value of Positionby following formula: ChildðPositionÞ ¼ LBðPositionÞ þ l2 ðUBðPositionÞLBðPositionÞÞ,

ð14Þ

ð15Þ

where l is a random real number in [0.1]. x0 is rounded down into nearest integer through the floor function bUc. If x1 and x2 are

where l1 is a real random number in [0,1], LB(U) and UB(U) refer to the lower bound and upper bound, respectively, of the value of Position obtained by CP.

0

x ¼ lx1 þ ð1lÞx2 ,

H. Li, D. Jiang / Computers & Operations Research 39 (2012) 1333–1344

Step 4. Check if child is feasible. If so, terminate; otherwise, let child ¼parent and go to step 2. To avoid excessive computational effort to obtain a feasible child, the procedure is terminated when reaching a certain number of iterations.

6. Computational experiments The CP model and algorithms were implemented in Cþþ using constraint programming libraries ILOG CP Optimizer 2.3 [40], which support the implementation of tree search algorithms that apply constraint propagations at the nodes of the tree [50]. The important features of these CP libraries include the following: (1) embedded constraint propagation algorithm, which is a variant of the AC-5 arc consistency algorithm [42]; (2) support for standard backtracking during CP search; (3) Cþþ classes for defining and implementing CP branching schemes. The GA code was implemented in MATLAB using the Genetic Algorithm Toolbox [51]. We generate the problem instances in the following way: The underlying general acyclic network is generated using the project scheduling problem generator called ProGen/Max [52]. We generate networks with 10, 20, 40 and 80 nodes. For each network size, 15 problem instances are randomly generated. pffiffiffiffi We use the standard functional form Fj ðXÞ ¼ X as in the literature (cf. [1,2,9,11]), where X is the net replenishment time. When demand follows the assumed normal distribution Nðm, s2 Þ, the safety stock prequired to achieve certain service level is ffiffiffiffi proportional to s X .

 Following [1] and [2], we generate lead time and capacity from uniform distribution. Specifically, lead time of each supply chain function is generated randomly from U[2,60], and capacity from U[50,70].  There appear to be three ways to generate hj in the literature: (i) it can be independently generated from a uniform distribution as in [2]; (ii) in the main computational experiments of [1], the holding cost rate hk of a downstream function k is obtained by adding a random number U to the largest holding cost rate hj of its immediate upstream function j, i.e. hk ¼hj þU. This is to reflect the fact that the downstream functions are more costly due to the value-adding process; (iii) a more realistic way to model such value-adding process, as in [9], is to calculate hk as the sum of holding cost rate hj of all its upstream functions, i.e. X hj ð16Þ hk ¼

same stopping criteria as in [1] and [2] using the latest version of CPLEX 12.1 [53]. All computations were performed on a Pentium IV PC with 3.2 GHz CPU and 1 G RAM.

6.1. Tuning GA Parameters Our computational experience indicates that, setting the total number of generations Nmax to be 1000 and the population size to be 300 achieves a good balance between solution quality and computational efficiency. The default parallel reproduction operation in the Matlab GA Toolbox is employed for the same pool of population in each generation, such that pcross þ pmut þ #Elite=300 ¼ 1 [51], where #Elite stands for the number of elite solutions directly kept to the next generation. In our implementation, we always keep the current best solution, i.e. we set #Elite ¼ 1. Fig. 6 shows the solution quality of pure GA for five selected instances with different value of pcross ranging from 0.1 to 0.9 (with an increment of 0.1). It appears that a crossover rate of 0.3 leads to better solution quality of pure GA. Thus in our computational experiments, we fix pcross to be 0.3. Fig. 7 illustrates the search process of CP–GA on the 11th 80node instance. The evolution process indicates that it improves the CP initial solutions effectively and quickly. Fig. 8 depicts the

14

Gap %

1340

Instance_1

12

Instance_2

10

Instance_3

8

Instance_4 Instance_5

6 4 2 0 0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Crossover rate Fig. 6. Performance of pure GA with different crossover rate.

j:ðj,kÞ A A

As shown in [1], instances generated through (16) can have a slower convergence rate and lead to a larger LP gap. We use (16) in our experiments to generate hj. As CP is often able to reach high-quality solutions fast, a time limit of 60 s is imposed for the pure CP algorithm: both the Default CP with the first-fail search strategy, and the Max-Reg CP with the maximal-regret search strategy. The CP search in the hybrid CP–GA algorithm is limited to 5 s. Of course, there is a trade-off between the initial CP search effort and solution quality. Our computational experience has shown that 5 s is usually sufficient for CP to generate high-quality initial population for GA. In order to evaluate the solution quality of our algorithms, we solve all problem instances to optimality through the exact piecewise linear MIP method proposed by [1]. We also compare our algorithms with the Iterative LP heuristic by [2]. Both the exact MIP and Iterative LP methods are implemented with the

Fig. 7. Convergence of fitness values for the 11th 80-node instance.

H. Li, D. Jiang / Computers & Operations Research 39 (2012) 1333–1344

average distance between individuals. It shows that our crossover and mutation operators are quite effective in keeping the population diversified during the search.

Fig. 8. Average distance between individuals for the 11th 80-node instance.

1341

6.2. Computational results Deviations of solutions found by our algorithm to optimal solution are reported as the following gap: (Obj  ObjMIP)/ObjMIP, where ObjMIP denotes optimal objective value obtained by the exact MIP method. We also report the CPU time (in seconds) for finding best solutions. Due to the probabilistic nature of GA, the pure GA and hybrid CP–GA were run five times for each instance. Both the best and average solution quality (in the format of [best gap, average gap]) and CPU time are reported. Complete computational results are provided in Tables 1–4. The pure CP algorithm itself performs well for the small 10node instances. Both the Default CP (using the first-fail rule) and the Max-Reg CP (using the maximal-regret rule) find optimal solutions to all fifteen instances. Max-Reg CP is more efficient than Default CP using only fractions of a second on average. The pure GA finds quality solutions within 1% gap, but spends significantly more computational time. The average gap of the Iterative LP method is 1.22%. For the fifteen 20-node instances, Max-Reg CP has an average gap of 1.70%, which clearly outperforms the Default CP in both solution quality and computational time. Such better performance is achieved by the maximal-regret rule for the CP tree search, which often finds good quality solutions early during the search. The hybrid CP–GA

Table 1 Computational results for the 10-node problem instances. Iterative LP

Default CP

Max-Reg CP

Pure GA

Instance #

No. Arcs

Gap %

CPU

# Iter

Gap %

CPU

Gap %

CPU

Gap %

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 Average

17 17 17 19 21 25 21 15 25 15 25 17 16 21 15 19

0.12 1.76 4.05 1.87 0.00 0.00 0.00 1.02 0.00 2.05 0.00 0.00 1.55 0.00 5.83 1.22

o0.01 o0.01 o0.01 o0.01 o0.01 o0.01 o0.01 o0.01 o0.01 o0.01 o0.01 o0.01 o0.01 o0.01 o0.01 o0.01

3 2 2 2 2 2 2 2 2 3 2 2 2 2 2 2

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

40.55 0.14 0.13 0.06 0.02 0.02 0.19 0.20 0.00 0.05 0.02 7.38 0.02 0.00 0.06 3.26

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0.28 0.31 0.13 0.05 o 0.01 o 0.01 0.11 1.77 0.03 0.34 0.00 0.45 0.14 o 0.01 0.06 0.24

[0.12, [1.30, [0.00, [0.00, [0.00, [0.00, [0.00, [1.49, [0.00, [1.07, [0.00, [0.00, [0.66, [0.00, [0.24, [0.33,

CPU 0.62] 2.31] 0.82] 2.67] 0.42] 0.42] 0.16] 1.76] 0.40] 2.95] 0.42] 0.55] 1.34] 0.23] 0.83] 1.06]

22.61 16.33 40.66 10.64 6.203 11.66 18.13 25.95 8.438 50.47 11.64 67.2 25.86 25.02 19.67 24.03

Table 2 Computational results for the 20-node problem instances. Instance #

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 Average

No. Arcs

32 35 36 30 27 50 30 39 24 33 28 35 44 43 40 35

Iterative LP

Default CP

Max-Reg CP

Pure GA

Gap %

CPU

# Iter

Gap %

CPU

Gap %

CPU

Gap %

5.83 2.20 1.41 0.00 3.13 0.21 0.00 3.01 0.03 10.59 1.16 2.24 3.17 3.53 8.40 2.99

0.02 0.05 0.02 0.06 0.11 0.11 0.03 0.06 0.02 0.08 0.14 0.02 0.11 0.03 0.08 0.06

3 3 3 5 4 3 3 3 3 4 3 3 3 2 3 3

0.00 6.17 0.03 3.47 9.20 0.94 8.00 0.03 10.34 4.28 2.90 3.17 0.00 0.00 0.01 3.24

3.08 41.27 0.45 0.75 0.33 43.14 3.16 1.47 0.97 o0.01 0.69 0.67 0.22 36.72 0.27 8.88

0.77 1.41 0.00 1.46 0.31 1.58 2.33 0.03 9.24 1.47 0.12 4.15 1.22 0.33 1.11 1.70

0.56 0.11 0.05 19.09 15.42 0.13 11.06 2.53 58.16 1.42 12.66 o 0.01 0.05 0.94 0.00 8.15

[0.00, [1.12, [0.00, [1.33, [0.11, [0.00, [0.44, [0.03, [5.68, [1.46, [0.00, [3.87, [0.00, [0.00, [0.00, [0.94,

0.10] 1.29] 0.00] 1.41] 0.20] 0.01] 1.22] 0.03] 6.61] 1.46] 0.03] 4.03] 0.31] 0.13] 0.22] 1.14]

CP–GA CPU

Gap %

7.34 47.72 2.00 9.85 53.64 10.85 45.22 2.00 49.39 2.00 7.45 56.16 37.8 22.28 20.03 24.92

[0.00, [1.04, [0.00, [1.33, [0.11, [0.00, [0.44, [0.03, [1.37, [1.46, [0.00, [3.86, [0.00, [0.00, [0.01, [0.64,

CPU 0.00] 1.04] 0.00] 1.46] 0.11] 0.03] 0.94] 0.03] 3.31] 1.46] 0.11] 3.92] 0.82] 0.00] 0.01] 0.88]

2.22 6.00 o 0.01 20.72 o 0.01 1.72 10.56 o 0.01 45.84 0.55 4.80 3.86 3.44 o 0.01 2.02 6.78

1342

H. Li, D. Jiang / Computers & Operations Research 39 (2012) 1333–1344

Table 3 Computational results for the 40-node problem instances. Instance #

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 Average

No. Arcs

75 84 89 70 85 70 77 49 58 64 89 75 69 85 76 74

Iterative LP

Default CP

Max-Reg CP

Pure GA

CP–GA

Gap %

CPU

# Iter

Gap %

CPU

Gap %

CPU

Gap %

CPU

Gap %

1.48 2.24 6.84 1.01 2.46 1.42 12.99 1.94 4.17 0.73 2.99 0.67 3.99 3.41 0.46 3.12

0.16 0.20 0.13 0.11 0.09 0.14 0.17 0.22 0.08 0.06 0.11 0.13 0.20 0.06 0.14 0.13

3 4 3 4 4 3 3 4 5 3 3 3 4 3 3 3

0.53 1.73 18.09 3.67 5.28 1.02 15.74 5.11 5.94 3.41 3.45 1.28 4.55 3.93 1.30 5.00

24.06 12.77 0.02 0.02 7.44 1.72 4.23 0.00 0.02 0.09 0.02 0.02 0.08 0.00 0.03 3.37

0.77 2.71 0.47 2.94 7.01 2.06 5.96 4.86 5.13 3.45 3.24 1.28 3.35 2.09 1.30 3.11

60.25 0.42 1.16 0.22 22.27 0.08 19.75 56.55 49.34 0.03 0.17 0.00 0.08 1.28 0.02 14.11

[2.39, 2.90] [3.63, 5.67] [7.98, 9.56] [2.81, 3.86] [2.48, 3.48] [1.70, 2.17] [8.65, 10.52] [3.81, 4.42] [2.71, 4.86] [3.27, 4.41] [3.79, 5.13] [1.16, 2.06] [5.00, 6.58] [2.18, 3.09] [2.22, 2.98] [3.59, 4.78]

315.86 290.81 160.52 428.58 376.27 489.38 373.25 267.47 453.98 450.33 348.44 553.94 389.41 339.69 531.44 384.62

[0.31, [0.11, [0.35, [0.70, [1.60, [0.39, [0.12, [3.20, [1.36, [1.10, [0.20, [0.66, [0.35, [1.85, [0.40, [0.85,

CPU 0.70] 0.34] 0.43] 1.40] 2.93] 1.43] 1.75] 3.51] 3.68] 2.01] 0.94] 0.88] 1.74] 1.95] 1.30] 1.67]

23.89 44.42 7.91 32.78 38.17 37.70 47.25 80.14 51.89 46.97 45.98 27.77 54.64 32.25 35.00 40.45

Table 4 Computational results for the 80-node problem instances. Instance #

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 Average

No. Arcs

149 134 213 111 142 93 102 139 172 116 141 87 171 108 112 133

Iterative LP

Default CP

Max-Reg CP

Pure GA

CP–GA

Gap %

CPU

# Iter

Gap %

CPU

Gap %

CPU

Gap %

CPU

Gap %

3.46 1.60 3.33 2.33 4.39 5.75 7.70 4.92 1.65 4.63 1.72 1.95 1.67 2.51 4.87 3.50

0.33 0.27 0.27 0.30 0.42 0.34 0.45 0.19 0.36 0.36 0.34 0.75 0.41 0.31 0.27 0.36

4 4 3 5 6 4 6 3 5 5 4 14 4 7 5 5

5.67 7.86 4.20 8.08 11.90 17.75 9.95 12.25 3.38 9.80 18.15 17.76 3.02 12.49 10.98 10.22

0.00 0.00 14.09 0.03 0.00 0.23 0.00 0.02 0.00 0.02 0.63 0.06 0.03 0.02 0.22 1.02

5.67 7.28 5.72 7.89 8.67 14.25 5.78 12.10 2.68 9.52 11.82 17.41 2.69 12.10 9.97 8.90

0.02 0.48 8.22 1.16 4.13 1.89 10.94 55.06 2.09 0.42 31.72 43.83 27.47 0.73 2.97 12.74

[8.98, 10.89] [8.55, 11.24] – – [13.45, 14.43] [18.14, 19.56] [11.38, 12.93] [15.52, 16.33] [5.41, 6.55] [12.52, 13.43] [17.96, 20.98] [18.00, 20.62] [5.13, 6.77] [12.67, 13.72] [15.37, 17.16] [12.54, 14.20]

93.06 82.28 – – 78.38 74.22 78.88 83.77 236.66 81.53 787.64 77.45 87.98 97.89 81.22 149.30

[1.86, [1.13, [2.74, [3.14, [5.74, [7.78, [3.11, [4.63, [1.70, [3.91, [5.39, [9.03, [1.86, [3.90, [4.87, [4.05,

algorithm is able to further improve the CP solution quality to 0.64% on average. Notably, CP–GA finds optimal solutions to four instances for which the pure CP did not find optimal solutions. It also significantly outperforms the pure GA, finding better quality solutions using less than a quarter of GA’s computational time. The average gap of the Iterative LP method is about 3%, which is better than the Default CP, but not as good as the other three. However, the Iterative LP appears to be the most efficient one, spending only 0.06 s (3 iterations) on average. For the fifteen 40-node instances, it takes the exact MIP method 1800 s (8 iterations) on average to reach optimal solutions. MaxReg CP has an average gap of 3.11%, which is significantly better than Default CP. CP–GA further improves Max-Reg CP solution quality to less than 1%. It also clearly outperforms the pure GA. The quality of Iterative LP is better than Default CP and pure GA, comparable to Max-Reg CP, but not as good as CP–GA. The Iterative LP again takes the least amount of time: only 0.13 s on average; the best performed CP–GA spends about 40 s on average. For the fifteen 80-node instances, the exact MIP method spends about 8 h (11 iterations) on average to reach optimality. The pure GA alone could not find feasible solutions to two instances within the running limit of 1000 generations. This suggests that relying purely on randomly generating solutions may not be efficient for large size problems, as much of the searching efforts can be spent

CPU 3.09] 2.91] 4.91] 5.01] 6.53] 9.08] 4.40] 6.88] 2.59] 4.95] 7.30] 11.00] 2.65] 8.71] 7.52] 5.84]

170.44 169.94 186.20 136.81 146.50 153.63 160.59 169.66 63.31 166.11 171.70 154.14 157.77 164.06 152.55 154.89

on exploring the infeasible region. Neither Default CP nor Max-Reg CP alone performs well in terms of solution quality. The hybrid CP– GA is able to significantly improve over both pure GA and CP, achieving an average gap of 4.05%. It is competitive with the Iterative LP method having an average gap of 3.5%. To directly compare with the heuristic methods in [2], we run the Max-Reg CP algorithm, the Iterative LP and the 2-piece iterative MIP (2P Iterative MIP1) algorithms on the 10-node 25-arc dense network instance taken directly from [2]. The comparison results are reported in Table 5. It is observed that the quality of the Iterative LP varies significantly with the delivery deadline of end products: its optimality gap increases from around 1% to 9% as the deadline increases from 0 to 10. The 2P Iterative MIP method performs better, whose gap also tends to increase as the deadline increases. The Max-Reg CP appears to be quite robust with respect to the change in deadline: it finds optimal solutions to all the ten instances.

1 In our implementation of the 2-piece iterative MIP method of Shu and Karimi [2], we corrected a typo in their original paper: the inequality at the bottom of Page 2901 should be Xjk r my instead of y r mXjk . For both methods, the same configuration is used as that in [2]. We further replace the stopping criterion in [2] with a more stringent one: 10  9. The same results were obtained. Additional tuning of algorithm parameters goes beyond the scope of this paper.

H. Li, D. Jiang / Computers & Operations Research 39 (2012) 1333–1344

Table 5 Comparison results on the dense network instance in [2]. Delivery Deadline Iterative LP

0 1 2 3 4 5 6 7 8 9 10

2P Iterative MIP

Max-Reg CP

Gap% # Iter CPU

Gap% # Iter CPU

Gap%

CPU

1.24 2.07 3.52 9.09 9.04 8.98 8.91 8.83 8.75 8.64 8.54

0 0 0 0 1.35 2.02 1.25 1.54 1.24 2.07 3.52

0 0 0 0 0 0 0 0 0 0 0

0.35 0.42 0.27 0.55 0.55 0.39 0.63 0.44 0.44 0.44 0.52

2 2 2 2 2 2 2 2 2 2 2

0.02 0.02 0.02 o 0.01 o 0.01 o 0.01 0.02 0.02 o 0.01 0.02 o 0.01

2 2 2 2 2 2 2 2 2 2 2

0.15 0.13 0.02 0.02 0.24 0.21 0.03 0.02 0.02 0.25 0.22

Despite its efficiency for solving large problems, one potential difficulty of the iterative heuristic methods is the existence of error bounds. We refer to [54] for empirical analysis of the error bounds incurred by linear approximation, and [55] for theoretical analysis on a priori bound for functional approximation in math programming. For the safety stock problem in generalized networks, the impact of such error bounds is evident in the Table 3 (p. 2903) of [2]: it prevents the iterative LP or the 2P Iterative MIP method from finding better solutions for even small size networks (with 10 to 20 nodes) compared to large problems. It is also evident in Table 1 and Table 2 of this paper: the Iterative LP method fails to find optimal solution to any of the fifteen 10-node problems, which have all been optimally solved by Max-Reg CP; for the fifteen 20-node problems, the Iterative LP yields an average gap of about 3%, whereas Max-Reg CP alone has a gap of 1.7% on average. It has been empirically shown in [1], as well as in this paper, that hj generated through (16) can lead to a larger error bounds. The results in Table 5 further suggest that the error bounds of both Iterative LP and 2P Iterative MIP tend to increase as the problem becomes less restrictive with a longer delivery deadline.

1343

algorithm finds optimal solutions to all the 10-node instances within fractions of a second. The hybrid CP–GA significantly improves ether pure CP or pure GA. It achieves an average optimality gap of less than 1% for the 20-node and 40-node instances in less than a minute, and about 5% gap for the 80-node instances in less than 3 min. It takes the exact MIP method about 8 h on average to optimally solve an 80-node instance. Our CP based algorithms are competitive with the state-of-the-art heuristics of [2]. The hybrid CP–GA achieves better solution quality for the 10-node, 20-node and 40-node instances than the Iterative LP method. For the 80-node instances, CP–GA finds competitive solutions using moderately more computational time. While the iterative heuristics of [2] are advantageous in computation efficiency, the hybrid CP–GA offers a novel solution approach that balances solution speed and quality. The CP based methods also appear to be robust with respect to the delivery deadline of end products. Given the fact that the addressed safety stock problem is often strategic or tactical in nature, the planning horizon is usually between twelve months and two years (cf. [9,11]). When it involves sourcing and new product launch decisions, the financial impact of an inventory positioning project can easily reach the scale of millions of dollars (cf. [9,56,57]). Thus algorithms with balanced speed and quality will be attractive to firms. In addition, the hybrid algorithmic framework can easily incorporate the LP based heuristic solutions into the initial population of GA to achieve even better solution quality. The idea of transforming a general acyclic supply chain network to a project scheduling network opens a number of opportunities for modeling and solving the safety stock placement problem. From the modeling perspective, such transformation makes it natural to model various renewable resources ubiquitous in complex manufacturing systems and supply chain networks in general. This provides a feasible and promising way to study safety stock placement in a resource-constrained environment. From the algorithmic perspective, it will be interesting to adapt other scheduling methods and/or local search based metaheuristics for the safety stock placement problem.

Acknowledgments 7. Conclusions and future research We propose a new modeling framework for the safety stock placement problem in general acyclic supply chain networks. Our approach transforms a general acyclic supply chain network into a project scheduling network. With such transformation, it becomes natural to model dependencies among the inbound and outbound service times, and capacity constraints involved, as temporal relationships in project scheduling. More importantly, our approach exploits a new project network perspective to the safety stock placement problem, which makes it possible to employ a large set of methods available in the project scheduling field to deal with the safety stock placement problem. In this paper, we devise constraint programming (CP) based methods for scheduling the inbound and outbound service times. With effective constraint propagation techniques and customized search strategies, our CP algorithm is able to find optimal solutions to small problems, and near optimal solutions for medium size problems fast. In order to further improve solution quality, we develop a scheme to integrate CP with a genetic algorithm (GA). A set of high quality CP solutions are used as the initial population for GA, which significantly reduces the overheads for GA alone to obtain initial solutions. CP also provides reduced time windows (bounds) for the decision variables of inbound and outbound service times, which significantly reduces GA’s search space. Computational study is conducted on randomly generated test instances with size ranging from 10 to 80 nodes. The pure CP

Dr. Jiang’s research was supported by the National Natural Science Foundation of China (Grant # 70871119/G0109), the National Social Science Foundation of China (Grant # 07XJY026), and the Science Innovation Foundation of Logistical Engineering University. The authors also thank three anonymous referees for their comments that helped improving contents and presentation of this paper.

Appendix

Definition 3. A floor function x is defined as the greatest integer that does not exceed a real number x. Definition 4. Let 0 r{x}o1 denote the fractional part of x, such that x¼x þ{x}. Lemma 1. If x1 and x2 satisfy the constraint x rb, and x0 is a convex combination of x1 and x2, then x0 rb holds. Proof. Since x0 is a convex combination of x1 and x2, we have x0 rb. Therefore, x0 r x0 r b. & Lemma 2. If ðx1 ,y1 Þ and ðx2 ,y2 Þ satisfy the constraint x y rb, where b is an integer. Let ðx0 ,y0 Þ be a convex combination of ðx1 ,y1 Þ and ðx2 ,y2 Þ, then x0 y0 rb holds.

1344

H. Li, D. Jiang / Computers & Operations Research 39 (2012) 1333–1344

Proof. Since ðx0 ,y0 Þ is a convex combination of two feasible solutions, it satisfies the constraint, i.e. x0  y0 rb. If suffices to show that x0 y0 rx0 y0 . There are the following two possibilities: (i) If {y0 }r{x0 }, then y0  y0 rx0 x0 , thus x0  y0 rx0  y0 rb. (ii) Consider now {y0 }4{x0 }. Since x0 þ{x0 }  y0 {y0 } rb, then x0  y0 rbþ{y0 }  {x0 }. Both {x0 } and {y0 } are in the range ½0,1Þ, and {x0 } o{y0 }, thus 0 o{y0 }  {x0 }o1. As bUc is a monotonic function, x0  y0 rbþ{y0 }  {x0 }¼ b. That is, x0  y0 rb. & Property 3. Let X1 ¼ ðx11 ,x12 ,. . .,x1,2n Þ and X2 ¼ ðx21 ,x22 ,. . .,x2,2n Þ be two vector of feasible solutions to the constraint system (2)–(6). If X 0 ¼ ðx01 ,x02 ,. . .,x02n Þ is a convex combination of X1 and X2, then X 0 ¼ ðx01 ,x02 ,. . .,x02n Þ is feasible to the constraint system (2)–(6). Proof. By definition, x0i is a convex combination of x1i and x2i for i ¼ 1,2,. . .,2n. Any constraint in the constraint system (2)–(6) takes either the form x rb or the form x  yrb, where b is an integer. Consider any constraint of the form x rb. Without loss of generality, assume it involves the ith decision variable in solution vector. Since x1i and x2i are two feasible solutions to x rb, and x0i is a convex combination of x1i and x2i, it follows from Lemma 1 that x0i rb, i.e. x0i is feasible to the constraint x rb. Consider any constraint of the form x  yrb. Without loss of generality, assume it involves the ith and jth decision variables in the solution vector. Since (x1i ,x1j Þ and (x2i ,x2j Þ are two feasible solutions to x  yrb, and ðx0i ,y0i Þ is convex combination of (x1i ,x1j Þ and (x2i ,x2j Þ, it follows from Lemma 2 that x0i x0j rb, i.e. ðx0i ,x0j Þ is feasible to the constraint x  yrb. Therefore, the solution vector X 0 ¼ ðx01 ,x02 ,. . .,x02n Þ is feasible to the constraint system (2)–(6). & References [1] Magnanti TL, Shen Z-J, Shu J, Simchi-Levi D, Teo C-P. Inventory placement in acyclic supply chain networks. Operations Research Letters 2006;34:228–38. [2] Shu J, Karimi IA. Efficient heuristics for inventory placement in acyclic networks. Computers and Operations Research 2009;36:2899–904. [3] Minner S. Strategic safety stocks in supply chains. Berlin: Springer; 2000. [4] Parlier GH. Transforming u.S. Army supply chains: strategies for management innovation. Business Expert Press; 2011. [5] Burton L. Strategic inventory positioning of navy depot level repairables. Master Thesis. Monterey, CA: Naval Postgraduate School; 2005. [6] Fisher. What is the right supply chain for your product? Harvard Business Review 1997:105–16. [7] Maloni MJ, Benton WC. Supply chain partnerships: opportunities for operations research. European Journal of Operational Research 1997;101:419–29. [8] Simchi-Levi D, Kaminsky P, Simchi-Levi E. Designing and managing the supply chain: concepts, strategies and case studies. New York: McGraw Hill; 2008. [9] Graves SC, Willems SP. Optimizing the supply chain configuration for new products. Management Science 2005;51(8):1165–80. [10] Minner S. Dynamic programming algorithms for multi-stage safety stock optimization. OR Spectrum 1997;19:261–71. [11] Graves SC, Willems SP. Optimizing strategic safety stock placement in supply chains. Manufacturing and Service Operations Management 2000;2(1): 68–83. [12] Sahni S. Computationally related problems. SIAM Journal on Computing 1974;3:262–79. [13] Kelly J. Critical-path planning and scheduling: mathematical basis. Operations Research 1961;9(3):296–320. [14] Malcolm DG, Roseboom JH, Clark CE, Fazar W. Applications of a technique for research and development program evaluation. Operations Research 1959;7(5): 646–69. [15] Demeulemeester EL, Herroelen WS. Project scheduling: a research handbook.Kluwer Academic Publishers; 2002. [16] Neumann K, Schwindt C, Zimmermann J. Project scheduling with time windows and scarce resources: temporal and resource-constrained project scheduling with regular and nonregular objective functions. New York: Springer; 2002. [17] Marriott K, Stuckey PJ. Programming with constraints. Cambridge, Massachusetts: The MIT Press; 1998.

[18] Aarts E, Lenstra JK. Local search in combinatorial optimization. Princeton: Princeton University Press; 2003. [19] Glover F, Kochenberger G. Handbook of metaheuristics. Springer; 2005. [20] Simpson KF. In-process inventories. Operations Research 1958;6(6):863–73. [21] Gallego G, Zipkin P. Stock positioning and performance estimation in serial production-transportation systems. Manufacturing and Service Operations Management 1999;1(1):77–88. [22] Ioannou G, Prastacos G, Skintzi G. Inventory positioning in multiple product supply chains. Annals of Operations Research 2004;126:195–213. [23] Graves SC, Willems SP. Strategic safety stock placement in supply chains. In: Proceedings of the 1996 MSOM Conference; 1996; 299–304. [24] Inderfurth K, Minner S. Safety stocks in multi-stage inventory systems under different service measures. European Journal of Operational Research 1998;106:57–73. [25] Lesnaia E. Optimizing safety stock placement in general network supply chains. PhD thesis. MIT; 2004. [26] Kolisch R. Make-to-order assembly management. Berlin: Springer; 2001. [27] Li H, Womer K. Modeling the supply chain configuration problem under resource constraints. International Journal of Project Management 2008;26(6): 646–54. [28] Baptiste P, Le Pape C, Nuijten W. Constraint-based scheduling: applying constraint programming to scheduling problems. New York: Springer; 2001. [29] Van Hentenryck P. The opl optimization programming language. Cambridge: MIT Press; 1999. [30] Brucker P. Scheduling and constraint propagation. Discrete Applied Mathematics 2002;123(1–3):227–56. [31] Brucker P, Knust S. Lower bounds for resource-constrained project scheduling problems. European Journal of Operational Research 2003;149:302–13. [32] Li H, Womer K. Scheduling projects with multi-skilled personnel by a hybrid milp/cp benders decomposition algorithm. Journal of Scheduling 2009;12(no. 3): 281–98. [33] Sellmann M, Fahle T. Constraint programming based lagrangian relaxation for the automatic recording problem. Annals of Operations Research 2003;118:17–33. [34] Rousseau L-M, Gendreau M, Pesant G, Focacci F. Solving vrptws with constraint programming based column generation. Annals of Operations Research 2004;130:199–216. [35] Pesant G, Gendreau M. A constraint programming framework for local search methods. Journal of Heuristics 1999;5:255–79. [36] Gendreau M. Constraint programming and operations research: comments from an operations researcher. Journal of Heuristics 2002;8:19–24. [37] Hooker J. Logic, optimization and constraint programming. INFORMS Journal on Computing 2002;14(4):285–321. [38] Milano M, Trick M. Constraint and integer programming. In: Milano M, editor. Constraint and integer programming: toward a unified methodology, vol. 27. Springer; 2004. [39] ILOG, Opl studio 6.0 user’s manual. IBM; 2010. [40] ILOG, Ilog cp optimizer 2.3 user manual. IBM; 2010. [41] Tsang E. Foundations of constraint satisfaction.Academic Press; 1993. [42] Van Hentenryck P, Deville Y, Teng C-M. A generic arc-consistency algorithm and its specializations. Artificial Intelligence 1992;57:291–321. [43] Dorndorf U, Pesch E, Phan-Huy T. A time-oriented branch-and-bound algorithm for resource-constrained project scheduling with generalized precedence constraints. Management Science 2000;46(10):1365–84. [44] Holland J. Adaptation in natural and artificial systems. Ann Arbor: University of Michigan Press; 1975. [45] Goldberg DE. Genetic algorithms in search, optimization and machine learning. Addison-Wesley Publishing Company, Inc; 1989. [46] Davis L. Handbook of genetic algorithm. New York: Van Nostrand Reinhold; 1991. [47] Gen M, Cheng R. Genetic algorithms and engineering optimization. New York: Wiley; 2000. [48] Michalewicz Z, Fogel DB. How to solve it: modern heuristics. Berlin: Springer; 2004. [49] Michalewicz Z. Genetic algorithmsþ data structure¼ evolution program. Berlin: Springer; 1996. [50] Le Pape C. Implementation of resource constraints in ilog schedule: a library for the development of constraint-based scheduling systems. Intelligent Systems Engineering 1994;3:55–66. [51] Mathworks—user’s guide for genetic algorithm and direct search toolbox. Mathworks Inc.; 2008. [52] Schwindt C. Generation of resource-constrained project scheduling problems with minimal and maximal time lags. Technical Report WIOR-489. Institute fur Wirtschaftstheorie und Operations Research, University of Karlsruhe; 1996. [53] ILOG, Ilog cplex 12.1 user’s manual. IBM; 2010. [54] Baumol WJ, Bushnell RC. Error produced by linearization in mathematical programming. Econometrica 1967;35(3-4). [55] Geoffrion AM. Objective function approximations in mathematical programming. Mathematical Programming 1977;13:23–37. [56] Amini M, Li H. Supply chain configuration for diffusion of new products: an integrated optimization approach. Omega 2011;39(3):313–22. [57] Li H, Amini M. A hybrid optimization approach to configure supply chains for new products diffusion: a case study of multiple-sourcing strategy. International Journal of Production Research. In press.