Heuristics for the multi-item capacitated lot-sizing problem with lost sales

Computers & Operations Research 40 (2013) 264–272

Contents lists available at SciVerse ScienceDirect

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

Heuristics for the multi-item capacitated lot-sizing problem with lost sales  a Nabil Absi a,n, Boris Detienne b, Ste´phane Dauze re-Pe´res a b

Ecole des Mines de Saint-Etienne, CMP-Site Georges Charpak-880 route de Mimet, F-13541 Gardanne, France Universite´ d’Avignon et des Pays de Vaucluse, 339, chemin des Meinajaries, 84911 Avignon Cedex 9, France

a r t i c l e i n f o

a b s t r a c t

Available online 6 July 2012

This paper deals with the multi-item capacitated lot-sizing problem with setup times and lost sales. Because of lost sales, demands can be partially or totally lost. To find a good lower bound, we use a Lagrangian relaxation of the capacity constraints, when single-item uncapacitated lot-sizing problems with lost sales have to be solved. Each subproblem is solved using an adaptation of the OðT 2 Þ dynamic programming algorithm of Aksen et al. [5]. To find feasible solutions, we propose a non-myopic heuristic based on a probing strategy and a refining procedure. We also propose a metaheuristic based on the adaptive large neighborhood search principle to improve solutions. Some computational experiments showing the effectiveness and limitation of each approach are presented. & 2012 Elsevier Ltd. All rights reserved.

Keywords: Lot-sizing Lost sales Lagrangian relaxation Probing heuristics Adaptive local search

1. Introduction Production planning consists in deciding how to transform raw materials or semi-finished products into final products in order to satisfy demands on time and at minimal cost. Determining lot sizes is a crucial decision in production planning; which consists in calculating the quantity to produce for each item at each time period. In industrial contexts, several constraints may complicate the problem. In particular, the fact that items need a resource makes the problem more complex. Indeed, this can lead to the impossibility to entirely satisfy demands when there is not enough capacity. Such an amount of unsatisfied demands is referred to as shortage on demand or lost sales. In this paper, we address the single-level, single-resource, Multi-item Capacitated Lot-Sizing problem with setup times and Lost Sales called MCLS-LS. MCLS-LS consists in planning the production of N items over a horizon of T periods in order to satisfy time-varying demands. Demands are given for each item at each period. The different parameters of the problem are as follows. Producing a unit of item i in period t incurs a production cost ait and requires vit units of the total capacity Ct. The holding cost of a unit of product i at the end of period t is git . Product dependent setup costs bit and setup times fit are incurred at each period where production takes place. The problem has the distinctive feature of allowing demand shortages (also called lost sales). Lost sales are particularly relevant in problems with tight capacities. Indeed, when the available resources are not sufficient to produce the total demand, the capacity is spread among the

n

Corresponding author. Tel.: þ33 442616656; fax: þ33 442616591. E-mail address: [email protected] (N. Absi).

0305-0548/$ - see front matter & 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.cor.2012.06.010

items by minimizing the total lost sales. Thus, we introduce in the model a unitary lost sales cost jit for item i at period t. These costs should be viewed as penalty costs and their values are high compared to other unitary cost components. The objective is to minimize total production, setup, holding, and lost sale costs. The classical Multi-item Capacitated Lot-Sizing problem with setup times (MCLS) is strongly NP-Hard [10]. Even the single item capacitated lot-sizing problem is NP-Hard in the ordinary sense [8]. The MCLS-LS problem is then NP-Hard. In fact, if we set lost sale costs to sufficiently high values, the problem to solve becomes the classical MCLS problem. Lot-sizing problems have been studied for five decades, with numerous references dealing with capacitated lot-sizing problems. For a recent survey on lot-sizing problems, the reader can refer to [17,18] for Multi-item Capacitated Lot-Sizing problems, and [9] for single-item capacitated lot-sizing problems. Different approaches were addressed in the literature in order to find near optimal heuristic solutions for the MCLS problem. Trigeiro et al. [26] were among the first to solve the MCLS problem. They propose a smoothing heuristic based on Lagrangian relaxation of the capacity constraints. At each step of the subgradient method, a heuristic is called to obtain a feasible solution from the current Lagrangian lower bound. Belvaux and Wolsey [6], Pochet and Wolsey [24] and Miller et al. [20] propose exact methods to solve MCLS problems by strengthening the Linear Programming (LP) relaxation using valid inequalities. Recently, Degraeve and Jans [12] propose a new Dantzig–Wolfe reformulation and branch-and-price algorithm for the MCLS problem. There has been little research on lot-sizing problems with demand shortages or lost sales. Sandbothe and Thompson [25] introduced the concept of shortages to the classical uncapacitated

N. Absi et al. / Computers & Operations Research 40 (2013) 264–272

single-item lot-sizing problem, they called it stockouts. The authors propose an OðT 3 Þ dynamic programming algorithm to solve the problem optimally. Aksen et al. [5] address the same problem but called this concept lost sales. They improved the previous complexity by proposing an OðT 2 Þ dynamic programming algorithm. Hwang and van den Heuvel [14] propose an OðT 4 Þ dynamic programming algorithm to solve optimally the classical uncapacitated single-item lot-sizing problem with lost sales, upper bounds on stocks and concave costs. They also propose an algorithm in OðT log TÞ and O(T) to solve respectively the uncapacitated lot-sizing problem with lost sales and non-speculative cost structure, and the same problem with nonincreasing selling prices. Absi et al. [4] address the uncapacitated single-item lot-sizing with time windows, early productions, lost sales and backlogging. The authors develop OðT 2 Þ solving algorithms for different variants of this problem. Berk et al. [7] study the single-item lot-sizing problem for a warm/cold process with immediate lost sales by defining some properties of the optimal solutions. Recently, Absi and KedadSidhoum [2] propose a branch-and-cut algorithm to solve the MCLS-LS problem with production time windows. They use a generalized version of Miller et al. [20] valid inequalities to strengthen the LP relaxations in the branch-and-bound tree. Absi and Kedad-Sidhoum [1,3] develop respectively MIP-Based heuristics to solve the MCLS-LS problem with additional industrial constraints and a Lagrangian relaxation approach to solve the MCLS-LS problem with safety stocks. The main contributions of this paper are twofold. First, we adapt the dynamic programming algorithm of Aksen et al. [5] to take any cost structure into account and use this algorithm in a Lagrangian relaxation approach to find good lower bounds. Second, we develop new Lagrangian heuristics based on a smoothing algorithm and a probing strategy to find near-optimal solutions. Generally, the smoothing heuristics that are presented in the literature are myopic. We show through our computational experiments (Section 5.1) that these heuristics are no longer competitive with recent versions of commercial mathematical programming solvers. The algorithm proposed in this paper is non-myopic, i.e. a probing heuristic is used at each step to evaluate promising moves. It provides better results than classical smoothing heuristics. We also show the efficiency of the adaptive local search principle in improving lot-sizing solutions. This principle already shows its performance in solving several vehicle routing problems [22,23]. Section 2 describes MIP formulations of the MCLS-LS problem. In Section 3, we present a Lagrangian relaxation approach based on the relaxation of capacity constraints, an adaptation of the dynamic program proposed by Aksen et al. [5] to solve the single item uncapacitated version of MCLS-LS, and a subgradient method. Section 4 describes the principle of the Lagrangian heuristics based on a probing strategy as well as a refining procedure. We propose 14 neighbor operators and integrate them in a metaheuristic (Adaptive Local Search) based on the selection principle of the adaptive large neighborhood search [22] to improve solutions. Computational experiments are shown in Section 5 to evaluate the effectiveness and limit of our procedures. Finally, Section 6 provides a short conclusion and future research directions.

2. Mathematical formulation of the MCLS-LS problem Different mathematical formulations were presented in the literature for the MCLS problem. Two mathematical formulations of the MCLS-LS problem are recalled in this section. The first one is a generalization of the classical formulation of the MCLS

265

problem usually called aggregate formulation. The second formulation is based on the facility location formulation initially proposed by Krarup and Bilde [19] for the uncapacitated singleitem problem, which is often called disaggregate formulation. 2.1. An aggregate formulation In the following, we present an aggregate formulation of the MCLS-LS problem. This formulation is addressed in several papers such as Trigeiro et al. [26]. The notations are given below: Sets and indices T: Number of periods. T ¼ f1, . . . ,Tg. N: Number of items. I ¼ f1, . . . ,Ng. i: index of an item, i ¼ 1, . . . ,N. t: index of a period, t ¼ 1, . . . ,T. Parameters: dit: demand for item i at period t. Ct: available capacity at period t. fit: setup time for item i at period t. vit: unitary resource consumption for item i at period t. ait : production unit cost for item i at period t. bit : setup cost for item i at period t. git : inventory unit cost for item i at period t. jit : lost sale unit cost for item i at period t. Variables: xit: the quantity of item i produced in period t. yit: binary setup variable, equal to 1 if item i is produced at period t (i.e. xit 40Þ, and 0 otherwise. sit: the inventory level of item i at the end of period t. rit: the lost sales of item i at period t. The aggregate formulation (noted AGG) of the MCLS-LS model is stated as follows: X min ðait xit þ bit yit þ jit r it þ git sit Þ ð1Þ i A I ,t A T

si,t1 þ r it þxit ¼ dit þ sit X ðvit xit þf it yit Þ rct

8i,t 8t

ð2Þ ð3Þ

iAI

xit rM it yit r it rdit

8i,t

8i,t

xit ,r it ,sit Z 0 yit A f0; 1g

ð4Þ ð5Þ

8i,t 8i,t

ð6Þ ð7Þ

The objective function (1) minimizes the total cost that aggregates production, setup, inventory and shortage costs. Constraints (2) are the inventory balance equations. Constraints (3) are the capacity constraints; the overall consumption must be lower than or equal to the available capacity. Constraints (4) relate the binary setup variables to the continuous production P variables. Mit can be set to minf Tt0 ¼ t dit0 ,ðct f it Þ=vit g. Constraints (5) define upper bounds on the lost sale variables. The domain definitions of the variables are defined in Constraints (6) and (7).

266

N. Absi et al. / Computers & Operations Research 40 (2013) 264–272

2.2. A disaggregate formulation

3.1. Lagrangian relaxation of capacity constraints

In the following, the MCLS-LS problem is formulated as a Facility Location problem. Generally, this formulation provides a strong lower bound for capacitated lot-sizing problems. Contrary to the previous formulation, the production variable is defined by the production period and the delivery period. We define a new production variable zitt0 that corresponds to the quantity of item i produced at period t to satisfy demand at period t 0 . Using these new variables, the inventory variables may be eliminated. In fact, the inventory sit is equal to the sum of all production quantities of item i at period q ðq rtÞ to satisfy a demand at period t 0 ðt 0 4tÞ (i.e. P P sit ¼ tq ¼ 1 Tt0 ¼ t þ 1 ziqt0 ). This formulation has more production variables than the previous one (NT2 vs. NT production variables), but generally provides better lower bounds. Let us define a0itt0 , the cost of producing a unit of item i at period t to satisfy a unit of demand at period t 0 : a0itt0 ¼ ait þ Pt01 q ¼ t giq . The disaggregate formulation (noted FAL) of MCLS-LS is stated as follows: ! T X X 0 min aitt0 zitt0 þ bit yit þ jit r it ð8Þ

Since the seminal paper of Trigeiro et al. [26] on using Lagrangian relaxation for the classical MCLS problem, different authors used or improved this approach to address Multi-item Capacitated Lot-Sizing problems. The main idea behind the method is to decompose a Multi-item Capacitated Lot-Sizing problem into N single-item uncapacitated lot-sizing problems by relaxing the linking constraints (the capacity constraints in this case). The capacity violation is penalized in the objective function using Lagrange multipliers pt Z0. When using Lagrangian relaxation on capacity constraints to derive a lower bound for the MCLS-LS problem, we obtain N ULS-LS problems that must be solved optimally. Since each item is considered separately, it is convenient to remove the item index i to facilitate the reading. The aggregate formulation of the ULS-LS problem is: X X min ðat xt þ bt yt þ jt r t þ gt st Þ þ pt ðvt xt þ f t yt ct Þ ð15Þ

i A I ,t A T

r it0 þ

t0 X

t0 ¼ t

zitt0 ¼ dit0

8i,t 0

st1 þ r t þ xt ¼ dt þ st

ð9Þ

r t rdt

t0 ¼ t

8t

ð10Þ

8i,t 0 ,t rt 0

8i,t 0

ð11Þ ð12Þ

r it0 ,zitt0 Z0

8i,t ,t rt

yit A f0; 1g

8i,t

0

0

ð16Þ

8t

ð17Þ

yt A f0; 1g

ð18Þ 8t

ð19Þ

8t

ð20Þ 2

iAI

zitt0 rdit0 yit

8t

8t

xt ,r t ,st Z 0

T X X X vit zitt0 þ f it yit r ct

r it0 r dit0

tAT

xt rM t yt

t¼1

iAI

tAT

ð13Þ ð14Þ

The objective function (8) minimizes the total cost that aggregates production, setup, inventory and shortage costs. Constraints (9) are the inventory balance equations. Constraints (10) are the capacity constraints. Constraints (11) relate the binary setup variables to the new continuous production variables. It means that, if production takes place at period t for item i then the quantity of item i produced at period t to satisfy demand at period t 0 cannot be higher than the demand of period t 0 . Constraints (12) define upper bounds on the lost sale variables. The domain definitions of the variables are defined in Constraints (13) and (14).

Aksen et al. [5] propose an OðT Þ dynamic programming algorithm called ACC to solve the ULS-LS problem. The ACC algorithm works only if the assumption at 4 jt is true for all t. This assumption is used in [5] to derive the following property. Property 1 (Aksen et al., 2003 [5]). There is an optimal solution such that the demand at any period t will be fully satisfied if procurement is made during that period (i.e. xt nr t ¼ 0). This property does not hold if the assumption at 4 jt is not satisfied. When performing a Lagrangian relaxation of the capacity constraints, the variables xt are penalized in the objective function and this penalty increases the unitary production cost, which may become larger than the shortage cost, even if it is not the case in the original instance. The following counterexample shows that Property 1 is not valid in this case. Consider the following ULS-LS problem with two periods. The parameters of the example are presented in Table 1. The details of the only optimal solution for this counterexample are given in Table 2. In this solution, r 1 nx1 ¼ 100, which contradicts Property 1. In the next section we propose an adaptation of the recursive dynamic programming algorithm ACC.

3. Lagrangian based lower bounds 3.2. Solving the ULS-LS problem When considering the MCLS-LS problem, capacity constraints (3) are the only constraints linking the items. By relaxing the capacity constraints and penalizing them in the objective function, the resulting problem can be decomposed into N single-item uncapacitated lot-sizing problems with lost sales (ULS-LS). Under the assumption that the unitary lost sale costs are larger than the unitary production costs, ULS-LS can be solved using an OðT 2 Þ dynamic programming algorithm [5]. In this section, we first present the resulting problem using the Lagrangian relaxation of capacity constraints, based on the aggregate formulation. Second, an adaptation of the dynamic program proposed by Aksen et al. [5] is developed. Finally, a subgradient optimization algorithm is proposed to find a lower bound for MCLS-LS.

Let us present the ACC algorithm and its adaptation to take into account solutions that are no longer dominated when Table 1 Counterexample for the ULS-LS problem. Periods (t)

1

2

dt

10 10 10 0 9

10 100 100 0 12

at bt

gt jt

N. Absi et al. / Computers & Operations Research 40 (2013) 264–272

Table 2 Optimal solution for the counterexample of Table 1. Periods (t)

1

2

xt yt rt st

10 1 10 10

0 0 0 0

267

(d) Update the step length of the subgradient method so that it satisfies the convergence conditions of the subgradient algorithm. The algorithm ends when a time limit is exceeded or when the step length calculated at Step 2 becomes smaller than a threshold E.

4. Lagrangian based upper bounds Property 1 is not valid. However, the other properties of dominant solutions are still valid (st1 nxt ¼ 0 8t, and r t nðdt r t Þ ¼ 0 8t). Definition 1. Let us recall the definitions used in the ACC algorithm:

 MCjt denotes the minimum between the marginal cost of  

losing the demand in period t and meeting this production P by a production in j before t, MC jt ¼ minðjt , aj þ t1 q ¼ j gq Þ Cj(t) is the minimum cost of handling periods /1,tS with st ¼0 when period j is the last production period ðj r tÞ. C(t) denotes the minimum possible cost of handling periods /1,tS with st ¼0. C(T) is the optimal objective function, and is computed by a forward recursion starting from Cð0Þ ¼ 0. The forward recursion of the ACC algorithm is recalled below.

Algorithm 1. The forward recursion of ACC.

2: 3: 4: 5: 6: 7:

P C 0 ðtÞ’ tq ¼ 0 jq dq ,8t A /0,TS for t ¼1 to T do C t ðtÞ’Cðt1Þ þ bt þ at dt for j¼ 1 to t1 do C j ðtÞ’C j ðt1Þ þ MC jt dt end for CðtÞ ¼ min fC j ðtÞg and jnt ¼ arg min fC j ðtÞg

8:

end for

1:

j A /0,tS

j A /0,tS

In order to adapt the ACC algorithm to solve the general version of ULS-LS, we need to define the parameter MC tt ¼ min ðjt , at Þ and to replace Line 3 of the previous forward recursion by C t ðtÞ’Cðt1Þ þMC tt dt þ bt . To obtain the optimal policy of procurements and losses over the time horizon, we adapt the backtracking procedure of the ACC algorithm. 3.3. Lagrangian relaxation algorithm The Lagrangian relaxation algorithm mainly consists of the following steps (see [13] for a general survey on Lagrangian relaxation). The Lagrangian relaxation parameters are fixed according to the scheme proposed by Held et al. [15]. 1. Initialization: All Lagrange multipliers pt are set to 0. 2. For a given iteration k: (a) Solve the relaxed problem using the dynamic programming algorithm described in Section 3.2 for each ULS-LS sub-problem. A current lower bound is found. If the lower bound value improves the current best one, then it is saved. (b) Update the Lagrange multipliers using the subgradient optimization method. (c) Stopping criteria: If any stopping condition is met, then save the best solution.

Different methods were proposed in the literature to obtain good feasible solutions. The most interesting one was proposed by Trigeiro et al. [26] to solve the classical MCLS problem. To obtain a feasible solution, the authors developed a smoothing heuristic starting from a usually infeasible Lagrangian solution. In this section, we first present the general idea behind a general smoothing heuristic and the limits of this method. Then, we present the general scheme of a heuristic approach based on an improved smoothing procedure. The different components of our heuristic are finally detailed.

4.1. Smoothing heuristics principle The general idea behind a smoothing heuristic is to start from the Lagrangian solution, which is most often infeasible, to construct a feasible solution. Production quantities are moved from one period to another in order to avoid capacity violations. Trigeiro et al. [26] propose a smoothing heuristic called TTM for the classical MCLS problem. The TTM smoothing heuristic is myopic, it performs slight modifications of the primal solution so as to fit the available capacity. TTM has a maximum of four passes and stops either when a feasible solution is obtained for the original MCLS problem or when the maximum allowed number of passes is reached (in this case, the procedure may not manage to remove overcapacity). TTM starts with a backward pass that moves production quantities from the end to the beginning of the horizon to eliminate capacity violations. If the solution is still infeasible, a forward pass is used to move production from the beginning to the end of the horizon. The third and fourth passes are identical to the first and second passes. When a feasible solution is found, a fix-up pass is used to eliminate unnecessary inventory. The advantage of this heuristic is its computational time which represents an average of less than one percent of the total computing time of the Lagrangian relaxation heuristic. Nevertheless, the TTM heuristic is myopic and does not consider the impact of moving a production quantity from one period to another. This decision is performed a priori, based on the Lagrangian cost divided by the quantity of eliminated overcapacity. Since the computational time of the heuristic is very small, we propose in Section 4.3 a non-myopic heuristic to obtain a good upper bound for the MCLS-LS problem.

4.2. General scheme of our heuristic In order to get good feasible solutions, we extend the approach of Trigeiro et al. [26] to the following three-phase procedure: Smoothing phase (Improved_Smooth). Starting from an optimal solution of the current Lagrangian subproblem (which usually does not respect the resource constraints), this phase aims at reducing the overconsumption of the resource. See Section 4.3. Cutting phase (Cut). This phase provides a feasible solution, starting from the potentially infeasible solution obtained at the smoothing phase. See Section 4.4.

268

N. Absi et al. / Computers & Operations Research 40 (2013) 264–272

Improving phase (Improve). This phase tries to improve the feasible solution by applying a hill-climbing approach. See Section 4.5. The general organization of the heuristic is summarized in Algorithm 2. Algorithm 2. Lagrangian heuristic. Let sL be a Lagrangian solution sn ’CutðsL Þ s’sL ; iter’1 while s is not feasible and iter rMaxIter do s’Improved_SmoothðsÞ s’CutðsÞ if costðsÞ o costðsn Þ then sn ’s end if end while if costðsn Þ is smaller than the best known upper bound then sn ’Improveðsn Þ end if return sn

4.3. Non-myopic smoothing heuristic This part of our method is an improvement of the algorithm proposed by Trigeiro et al. [26]. As already stated, in their paper, the authors propose to remove the overconsumption of the resource at a given period by moving production quantities to an earlier period (backward phase) or to a later period (forward phase). Let us summarize the backward phase used by Trigeiro et al. [26] in Algorithm 3. Greedy_Choice returns an item in and a time period t n o t according to a greedy criterion which estimates the Lagrangian cost increase per unit of resource saved at period t. The forward phase is similar, except that the greedy choice is made among time periods after t. Algorithm 3. Backward phase of Trigeiro et al. Smoothing procedure TrigeiroBackwardðÞ. for t from T down to 1 do while resource constraint at t is violated do Greedy_Choiceðt,in ,t n ,backwardÞ Move the correct amount of production quantity of in from t to tn end while end for

Although it is modulated by the Lagrangian multipliers, this criterion is very local since it considers only the current state of the plan and the state when the modification is made. The idea of our non-myopic heuristic is to use an adaptation of the algorithm of Trigeiro et al. [26] to estimate the quality of a move (probing strategy), instead of the initial criterion. Algorithm 4. Evaluate_Moveðt,i,t 0 ,directionÞ. Move the correct production quantity of i from t to t 0 if direction¼Backward then Perform TrigeiroBackward end if Perform TrigeiroForward Perform TrigeiroBackward Perform the Cut heuristic return the value of the solution obtained

This adapted version is summarized in Algorithm 4, and consists of a truncated run of the classical algorithm, followed by the cutting heuristic. It returns the cost of the solution that would be obtained if the classical criterion was applied until the end of the procedure after the considered move, thus an upper bound of the cost of the final solution after the move. To design an improved backward or forward phase of the smoothing algorithm, we replace the use of Greedy_Choice by a procedure that finds the best item to move according to Evaluate_Move. Finally, an improved smoothing pass Improved_Smooth() consists in an improved forward pass followed by an improved backward pass. One can notice that, although Trigeiro et al. [26] aim at obtaining a feasible solution at this step, this phase is, in our approach, devoted to finding solution that is a good tradeoff between reducing infeasibility and minimizing costs. 4.4. Cutting phase Unlike the classical MCLS problem, finding a feasible solution to MCLS-LS is easy: any solution that does not satisfy the resource constraints can be extended to a feasible solution by losing all unsatisfied demands. Moreover, if we assume that the production periods for each item are known (thus, the cost and resource consumption of the corresponding setups are fixed), then the resulting sub-problem consists in solving a linear program that finds an optimal production plan, i.e. quantities and lost sales, that is consistent with the fixed setups. Hence, the cutting phase of our algorithm is based on a simple variable-fixing method: for each item, the production periods in the starting solution (e.g. derived from the smoothing procedure) ~ r~ , s~ Þ be are fixed while the others are forbidden. Formally, let ðx, the starting solution. The cutting heuristic sets the values y it ¼ 1 if x~ 4 0, and y it ¼ 0 if x~ ¼ 0, 8i,t, and then solves the linear program (1)–(7) where variables y are fixed to the values of y. The resulting solution is then iteratively improved using the dominance Property 2. In the sequel, this step is called Refining procedure. This property allows setups that were fixed but are not used at a profitable level in the current solution to be removed. It is used in a loop that solves the linear program according to the current set of fixed variables, removes all setups satisfying the dominance property in the resulting solution, and iterates until no setup variable has changed. Property 2. Let ðx,y,r,sÞ be a feasible solution of (1)–(7) whose cost is z, and let iA I ,t A T . If y it ¼ 1 and bit þ ait x it 4 jit x it , then there exists at least one feasible solution ðx 0 ,y 0 ,r 0 ,s 0 Þ such that y 0it ¼ 0 and y 0i0 t0 ¼ y i0 t0 , 8ði0 ,t 0 Þ A I  T fði,tÞg, whose cost is smaller than z. Proof. Straightforward: the solution where all variables are identical except y 0it ¼ 0, x 0it ¼ 0 and r 0it ¼ x it is trivially feasible and has a cost smaller than z. & 4.5. Improvement phase This phase consists in a metaheuristic (called Adaptive Local Search ‘‘ALS’’) based on 14 neighborhoods that are selected according to probabilities that are updated throughout the algorithm. The idea of each neighborhood is similar to the one used in the cutting phase: The set of setup variables is fixed, and the resulting LP is solved. As it is easy to compute an optimal solution corresponding to a vector of fixed variables, a solution is coded as its vector of setup variables. A neighbor of a solution is obtained by a slight change in this vector. As the number of neighborhoods used in this method is relatively large, we used the same principle of the adaptive large

N. Absi et al. / Computers & Operations Research 40 (2013) 264–272

neighborhood search (ALNS), described in [23] and adapted by [21] for lot-sizing problems, to choose the most promising neighborhood. The choice of neighborhood to use is done dynamically using probabilities that are updated according to the performance of the neighborhoods. This principle is called roulette wheel. The course of the algorithm is divided into time segments. In our implementation, a segment corresponds to 1000 iterations. During each segment t, the algorithm calculates the probability Pit of choosing the ith neighborhood as follows. For each segment t and for each neighborhood i, we compute an observed performance wit ¼

269

demands. This distinction allows to consider different performance and probabilities for their two alternatives, which is beneficial when solving instances with different characteristics. 1. Remove one setup: choose ði,tÞ A I  T such that y it ¼ 1, and set y0 it ¼ 0. The two next neighborhoods add a setup: choose ði,tÞ A I  T such that y it ¼ 0 and set y 0it ¼ 1. They differ from each other by the way (i,t) is chosen: 2. Chose an item i with lost demand at t: r it 4 0. 3. Chose an item i with no lost demand at t: r it ¼ 0.

w it ait

w it corresponds to the amount of objective function improvement induced by the moves chosen among neighborhood i. ait is the computing time spent during segment t evaluating neighborhood i. Notice that we use this criterion instead of the number of calls to neighborhood i because our neighborhoods may be very different from each other from a computational complexity point-of-view. A relative performance is then calculated as

The two next neighborhoods add a setup for an item at one period, after releasing some production resource at this period by removing the setup of another item: choose ði1 ,tÞ A I  T such that y i1 t ¼ 0, ði2 ,tÞ A I  T such that y i2 t ¼ 1 and set y 0i1 t ¼ 1 and y 0i2 t ¼ 0. Once again, they differ from each other by the way ði1 ,tÞ is chosen:

w P it ¼ P it j wjt

4. Chose an item i1 with lost demand at t: r i1 t 4 0. 5. Chose an item i1 with no lost demand at t: r i1 t ¼ 0. 6. Try to remove all the production quantity of a given item at one period, and to balance it by creating a new setup for this item at another period: choose ði,t 1 Þ and ði,t 2 Þ such that y it1 ¼ 1 and y it2 ¼ 0, and set y 0it1 ¼ 0 and y 0it2 ¼ 1. 7. Try to group the production of one item split over two periods into a single third period: choose ði,t 1 Þ, ði,t 2 Þ and ði,t 3 Þ such that t 1 at 2 , y it1 ¼ 1, y it2 ¼ 1 and y it3 ¼ 0 and set y 0it1 ¼ 0, y 0it2 ¼ 0 and y 0it3 ¼ 1.

At the end of each segment, we update the probability Pit based on P it and previous values of the probability. A reaction factor r controls how much the roulette weight depends on the score in the most recent time segment t: P i,t þ 1 ¼ rP it þ ð1rÞPit

ð21Þ

Moreover, we slightly modified the classical principle of the ALNS procedure by centering the probabilities when no solution improvements have been made during segment t: 0

P i,t þ 1 ¼ r 0 Pit þ

1r #neighborhoods

The idea behind this is that, to improve the current solution, the algorithm may need to use a neighborhood which has been discarded because it was not efficient at the beginning of the search. So, it could be beneficial to increase the weights of lowweighted neighborhoods in this case. Another argument is that, if we only use formula (21) to update smoothened weights when no improvements are made, weights are not modified although all neighborhoods were as good as the others in the last segment, and should tend to have an equal chance to be selected by the roulette wheel. Preliminary experiments led us to choose r¼ 0.1 and r 0 ¼ r=20. Inside our main algorithm, ALS search can be called several times, to improve different basic solutions. The first time it is called, probabilities are uniformly initialized: P i1 ¼

1 #neighborhoods

To initialize subsequent ALS runs, we base on the global performance of the neighborhoods defined as Pt w 0 0 wGit ¼ Ptt ¼ 1 it 0 t0 ¼ 1 ait The probabilities are finally initialized using the relation: wG 1r 00 P i,t þ 1 ¼ r 00 P it G þ #neighborhoods j wjt r 00 is empirically fixed to r 00 ¼ 34. At each step of the ALS procedure, a neighbor is chosen among the following neighborhoods of the current solution. Some of them are similar but deal with items either with or without lost

The two next neighborhoods swap the production periods of two items: choose ði1 ,t 1 Þ and ði2 ,t 2 Þ such that i1 a i2 , y i1 t1 ¼ 1, y i2 t2 ¼ 1, y i1 t2 ¼ 0 and y i2 t1 ¼ 0, and set y 0i1 t1 ¼ 0, y 0i2 t2 ¼ 0, y 0i1 t2 ¼ 1 and y 0i2 t1 ¼ 1. They differ from each other by the way ði1 ,t 1 Þ is chosen: 8. Chose an item i1 with lost demand at t1: r i1 t1 40. 9. Chose an item i with no lost demand at t: r i1 t1 ¼ 0. The last neighborhoods share the following principle: add a setup for an item at one period, after releasing some production resource at this period by shifting the setup of another item to another period. Formally, choose ði1 ,t 1 Þ and ði2 ,t 2 Þ such that y i1 t1 ¼ 0 and y i2 t1 ¼ 1 and y i2 t2 ¼ 0. Set y 0i1 t1 ¼ 1, y 0i2 t1 ¼ 0 and y 0i2 t2 ¼ 1. They differ by the way ði1 ,t 1 Þ and ði2 ,t 2 Þ are chosen:

10. Add a setup for an item with lost demand at the period considered, and randomly shift the setup of another item to a nearby time period: r i1 t1 40, t 2 ¼ t 1 1 or t 2 ¼ t 1 þ 1. 11. Add a setup for an item with lost demand at the period considered ðr i1 t1 4 0Þ, and shift the setup of another item to the best possible time period t2 between t 1 2 and t 1 þ2, t 2 at 1 . More precisely, each neighbor with t 2 A ft 1 2,t 1 1, t 1 þ1,t 1 þ 2g is evaluated using the evaluation function of the ALS procedure. The value of t2 corresponding to the best neighbor is kept. 12. Similar to neighborhood 11, except that the setup is added for an item with lost demand at the next time period ðr i1 t1 þ 1 4 0Þ. The idea is the same, but we use storage for one period. 13. Similar to neighborhood 11, except that we add a setup for an item without lost demand at the period considered ðr i1 t1 ¼ 0Þ. 14. Add a setup for an item with or without lost demand at the period considered, and shift the setup of another period to a faraway period: t 22 = ½t 1 2,t 1 þ2.

270

N. Absi et al. / Computers & Operations Research 40 (2013) 264–272

5. Computational experiments In this section, we present computational experiments resulting from the application of our Lagrangian relaxation based heuristics to the MCLS-LS problem. Our algorithms are coded in the Cþþ programming language using Microsoft Visual Studio 2008. For the sake of comparison, we use the callable CPLEX 12.1 library [16] to solve the MILP problems. CPLEX 12.1 is used with the default settings. In order to solve the linear programs in our heuristic methods, we use the Open Source Clp 1.11 solver [11] from the COIN-OR library. The tests were done on an Intel Xeon CPU 2.67 GHz PC with 3.48 GB RAM under the Window XP operating system. CPU times are given in seconds. We report computational tests on extended instances from the 540 instances of the lot-sizing library LOTSIZELIB Trigeiro et al. [26]. These instances are characterized by 10, 20 or 30 items and 20 periods. These instances are also characterized by enough capacity to satisfy all demands over the planning horizon. Since these instances have enough capacity to satisfy all demands over the planning horizon, we made some modifications to induce shortages. We derived 1080 new instances from the classical ones by decreasing the resource capacity. The available resource capacity ct at period t is multiplied by a factor (capacity factor) equal to 0.8 or 0.6. Lost sales costs are considered as penalty costs and their values must be larger than other cost components. Therefore, jit is fixed such that jit b maxi0 ,t0 fai0 t0 ; gi0 t0 g. We carried out a comparison between the following methods:

 The CPLEX12 method implements the standard branch-and-cut



of CPLEX 12.1 solver. We add the parameter AGG (respectively FAL) if the aggregate formulation (respectively the disaggregate formulation) of the MCLS-LS problem is used. The LR algorithm implements the Lagrangian relaxation of capacity constraints which gives a lower bound for the MCLS-LS problem. During the Lagrangian relaxation, an upper bound is computed using the non-myopic smoothing heuristic described in the last section. This smoothing heuristic has three parameters. The first parameter (Prob) corresponds to the number of improved smoothing passes performed by the method (probing depth). The remaining passes (up to a total of four passes) are classical passes from the TTM algorithm. The second parameter (Refine) corresponds to activating or deactivating the refining procedure during the Lagrangian relaxation. The last parameter (ALS) corresponds to activating or deactivating the ALS method described in the previous section.

In our implementation, the threshold for the stopping criterion of the subgradient procedure is set to E ¼ 105 . This relatively small value allows generating a large number of different sets of near-optimal Lagrangian multipliers and Lagrangian solutions. For every 1000 iterations, the relative improvement of the best lower bound for the last 1000 iterations is calculated. The heuristic procedure is called when the sub-gradient procedure has nearly converged, i.e. when this value is less than 1%. This prevents wasting time computing a solution from bad quality Lagrangian multipliers. To keep the presentation of the computational results concise, average gaps for all instances sharing the same value of given parameters are reported on a single line. 5.1. Impact of probing depth using the probing procedure In this section, we study the impact of the depth parameter when using the probing procedure. Table 3 summarizes the computational behavior of the LR method by varying the depth

Table 3 Impact of probing depth on the probing procedure. Probing depth CF 0.6

0.8

0

1

2

3

N

Gap (%) Time Gap (%) Time Gap (%) Time Gap (%) Time

10 20 30 10 20 30

10.29 6.64 3.79 15.75 8.16 6.49

0.6 1.2 1.8 0.6 1.0 1.5

7.75 4.25 2.54 10.48 5.41 4.42

1.1 3.3 7.0 0.9 2.1 4.1

6.09 3.44 2.18 9.00 4.71 3.93

2.4 8.4 18.6 1.6 4.4 9.7

5.48 3.28 2.10 8.45 4.49 3.89

3.6 13.3 29.4 2.2 6.4 15.1

8.52

1.1

5.81

3.1

4.89

7.5

4.62

11.6

Mean

Table 4 Impact of the refining and local search procedures. Improving procedures

LR

CF

N

gap (%)

0.6

10 20 30

0.8

10 20 30

Mean

LR þRefine

LRþ Refineþ ALS

time

gap (%)

time

gap (%)

time

6.47 3.72 2.31

2.7 9.7 21.6

5.48 3.28 2.10

3.6 13.3 29.4

2.97 1.96 1.21

10.1 31.5 63.1

9.84 4.70 4.01

1.7 4.7 11.1

8.45 4.49 3.89

2.2 6.4 15.1

4.84 2.32 2.10

8.8 23.2 44.3

5.17

8.6

4.62

11.6

2.57

30.2

of the probing heuristic (0, 1, 2 or 3). Average gaps (gap) and average CPU times (time) are computed for each set of instances characterized by the following parameters: A number of items (N) and capacity factor (CF). The LR method stops if a time limit of 100 s is achieved or the stopping criteria of the subgradient algorithm are satisfied. In these computations, the refining procedure is activated and the local search is not activated. In this section, gaps are calculated as the percentage difference between the best upper bound and the best lower bound compared to the best upper bound. From Table 3, we note that the differences between the gaps of LR with probing depths 1 and 2 and LR without probing are rather large (gaps are improved respectively by 32% and 43%), while the differences between the gaps of LR with probing depth of 2 and probing depth of 3 are relatively small (5%). The computational time of the LR method increase almost linearly between probing depths 1 and 3. Note also that LR with probing depth of 1 and 2 give a good trade-off between gaps and solution times, and probing depth 3 provides the best results. From these computations we conclude that the probing improves considerably the classical TTM smoothing heuristic. TTM gaps are obtained by setting the probing depth to 0.

5.2. Impact of refining and local search procedures In this section, we discuss the advantage of using the refining and local search procedures. Recall that the refining procedure uses the dominance property (Property 2) to improve the obtained solution. Tables 4 summarizes the computational behavior of the LR method without improvement procedures, when only the refining procedure is activated and when both the refining and local search procedures are activated. In these computations, the depth of the LR method is fixed to 3 (since this parameter gives the best results). The stopping criterion is the same as in the previous section. Average gaps and average CPU times are computed for

N. Absi et al. / Computers & Operations Research 40 (2013) 264–272

Table 5 Relative performance of the different neighborhoods with respect to Capacity Factor values. Neighborhood

Avg. relative performance

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

CF ¼ 0:6

CF ¼ 0:8

0.17 0.06 0.05 0.10 0.09 0.05 0.05 0.08 0.07 0.04 0.12 0.06 0.03 0.03

0.14 0.06 0.11 0.06 0.12 0.06 0.05 0.07 0.07 0.03 0.08 0.04 0.06 0.05

Max. relative performance

0.62 0.75 0.69 0.64 0.58 0.43 0.69 0.44 0.70 0.69 0.54 0.66 0.31 0.38

CF

N

LR/AGG

LR/FAL

LR time

AGG time

FAL time

0.6

10 20 30 10 20 30

1.58 1.51 1.52 2.96 2.96 2.98

0.99 0.99 1.00 0.97 0.99 0.98

0.2 0.4 0.6 0.2 0.4 0.5

0.0 0.1 0.1 0.0 0.1 0.1

0.2 0.5 0.6 0.2 0.5 0.6

2.25

0.99

0.4

0.1

0.4

Mean

criteria of the subgradient algorithm are satisfied. In these computations, the improvement procedures are not activated. Table 6 shows that the LR method provides better lower bounds than LBCPLEX12AGG with an average CPU time that is less than one second. Note also that the LR method obtains almost equivalent lower bounds with an average CPU time that is less than one second compared to CPLEX12FAL . 5.4. Repeated subgradient procedure

Table 6 Quality of LR lower bounds.

0.8

271

each set of instances characterized by the following parameters: The number of items (N) and the capacity factor (CF). From Table 4, it can be noted that, when activating the refine procedure, gaps are considerably improved whereas CPU times do not increase very much. We can also observe that gaps are considerably reduced when activating the adaptive local search procedure. Here again, CPU times do not increase so much. Table 5 shows the interest of using adaptive probabilities for selecting the neighborhood used at each step of the ALS procedure. The relative performance reported for neighborhood i is calculated as the relative global performance at the end of the P LR_MS algorithm (described in Section 5.4): wGit = j wGjt . For some neighborhoods, one can observe significantly different average performance according to the different values of the Capacity Factor. Moreover, each of them proves to be very useful for at least one instance, even if its average performance is modest. Likewise, the minimum relative performance is null for every neighborhoods, showing that each of them is useless for at least one instance of our test bed.

5.3. Quality of Lagrangian relaxation lower bounds Table 6 shows a comparison between the Lagrangian relaxation lower bounds (LBLR), and the lower bounds obtained by solving the linear relaxation of the aggregate formulation ðLBCPLEX12AGG Þ and linear relaxation of the disaggregate formulation ðLBCPLEX12FAL Þ using CPLEX 12.1. The comparison is shown through the average ratio between (LBLR) and ðLBCPLEX12AGG Þ denoted (LR/AGG), the average ratio between (LBLR) and ðLBCPLEX12FAL Þ denoted (LR/FAL), and the average CPU time of LR (LR time), LBCPLEX12AGG (AGG time) and LBCPLEX12AGG (FAL time). The LR algorithm stops if a time limit of 1 seconds is achieved or the stopping

Since the LR Method obtains good lower bounds within relatively small CPU times, we developed a new variant of the LR approach that we call multi-start Lagrangian relaxation approach (denoted LR_MS). Contrary to the LR method that stops when the subgradient conditions are satisfied, the LR_MS method restarts with new Lagrangian multipliers. These new multipliers are calculated from the last obtained multipliers increased by randomly generated values. This method stops when a given CPU time is reached. Table 7 summarizes the computational behavior of LR_MS with and without the metaheuristic procedure (ALS) for different probing depths (0, 1, 2, or 3) and CPLEX12AGG and CPLEX12FAL methods subject to a 100-s time limit. Table 7 reports average gaps that are computed for each set of instances characterized by a number of items (N) and by the same capacity factor (CF). Gaps are calculated as the percentage difference between the best upper bound and the best known lower bound compared to the best known lower bound. In these computations, the refining procedure and the ALS are always activated for the LR_MS method. From Table 7, it can be noted that, for the LR_MS method, the best gaps are obtained with the probing depth 1, although best results are obtained with the probing depth 3 for the LR method. A simple explanation is that the probing depth 1 is the best tradeoff between the quality of the non-improved heuristic solution and the time spent to find it: Finding and improving different solutions obtained with probing depth 1 is preferable to finding and improving fewer better quality solutions obtained with probing depths 2 and 3. We can also notice that LR_MS outperforms CPLEX12AGG and CPLEX12FAL . CPLEX12FAL gives better results than LR_MS on one type of instances characterized by CF¼0.8 and N ¼30 but the difference is only 0:01%. The standard deviation of the gaps obtained by LR_MS, CPLEX12FAL and CPLEX12AGG are respectively 4:13%, 5:23% and 6:39%. In addition, the number of optimal solutions obtained by the three methods are respectively 9, 30 and 30. A solution is considered optimal if its gap is lower than 0:01%. These last observations show that the LR_MS method is more stable than CPLEX12FAL and CPLEX12AGG but provides less optimal solutions. One of the main advantage of the probing heuristic is that it can be adapted for lot-sizing Table 7 Computational results: LR_MS vs. CPLEX (100 seconds). Method

CPLEX12AGG (%)

CPLEX12FAL (%)

Probing depth CF

LR_MS (%) 0

1

2

3

N 0.6

10 20 30

2.74 2.32 1.64

2.54 1.81 1.19

2.22 1.75 1.13

2.16 1.62 1.03

2.21 1.68 1.15

2.27 1.79 1.21

0.8

10 20 30

4.04 2.26 2.58

3.87 2.02 1.85

3.65 1.97 2.04

3.48 1.84 1.86

3.95 1.89 1.90

3.85 1.97 1.97

2.60

2.21

2.13

2.00

2.13

2.18

Mean

272

N. Absi et al. / Computers & Operations Research 40 (2013) 264–272

problems with additional constraints such as backlogging, time windows and safety stocks.

6. Conclusion In this paper, we addressed the Multi-item Capacitated Lot-Sizing problem with setup times and lost sales (MCLS-LS). To obtain good lower bounds, we used a Lagrangian relaxation of the capacity constraints. We also developed an algorithm that solves the induced single-item sub-problems optimally using an adaptation of the OðT 2 Þ dynamic programming algorithm of Aksen et al. [5]. To find feasible solutions, we proposed an original approach based on a non-myopic heuristic that uses a probing strategy and a refining procedure. We also proposed an adaptive local search principle to improve the upper bound during the Lagrangian relaxation. We obtain bounds of the same quality as a tight formulation solved using a standard solver (CPLEX 12.1). We also proposed a variant of the classical Lagrangian relaxation method that is based on a multi-start strategy. Computational experiments show that the probing strategy improves the solution of the classical smoothing heuristic. The method that combines the probing, the multi-start and adaptive local search strategies provides very good and stable results. The main advantage of this heuristic is its flexibility to integrate additional constraints. As a perspective to this work, it will be interesting to use Lagrangian lower bounds and upper bounds in a branch and bound approach. It will also be interesting to use the same approach to solve Multiitem Capacitated Lot-Sizing problems with production time windows, early deliveries and lost sales by using the dynamic programming algorithm proposed in [4] and adapting again the smoothing heuristic proposed by Trigeiro et al. [26]. Another perspective would be to design our algorithm as a multi-thread program to take advantage of multi-core platforms.

Acknowledgments The authors would like to thank anonymous referees for their careful reading, thorough reviews and valuable comments that helped us improving the quality of the paper. References [1] Absi N, Kedad-Sidhoum S. MIP-based heuristics for multi-item capacitated lot-sizing problem with setup times and shortage costs. RAIRO-Operations Research 2007;41(2):171–92. [2] Absi N, Kedad-Sidhoum S. The multi-item capacitated lot-sizing problem with setup times and shortage costs. European Journal of Operational Research 2008;185(3):1351–74.

[3] Absi N, Kedad-Sidhoum S. The multi-item capacitated lot-sizing problem with safety stocks and demand shortage costs. Computers & Operations Research 2009;36(11):2926–36. [4] Absi N, Kedad-Sidhoum S, Dauze re-Pe´re s S. Uncapacitated lot-sizing problem with production time windows, early production, backlog and lost sale. International Journal of Production Research 2010;49(9):2551–66 2011. [5] Aksen D, Altinkemer K, Chand S. The single-item lot-sizing problem with immediate lost sales. European Journal of Operational Research 2003;147: 558–66. [6] Belvaux G, Wolsey LA. bc-prod: a specialized branch-and-cut system for lot-sizing problems. Management Science 2000;46:724–38. [7] Berk E, Toy AO, Hazir O. Single item lot-sizing problem for a warm/cold process with immediate lost sales. European Journal of Operational Research 2008;187(3):1251–67. [8] Bitran G, Yanasse HH. Computational complexity of the capacitated lot size problem. Management Science 1982;28:1174–86. [9] Brahimi N, Dauze re-Pe´re s S, Najid NM, Nordli A. Single item lot sizing problems. European Journal of Operational Research 2006;168(1):1–16. [10] Chen WH, Thizy JM. Analysis of relaxations for the multi-item capacitated lot-sizing problem. Annals of Operations Research 1990;26:29–72. [11] COIN-OR Clp 1.11 solver callable library. /http://www.coin-or.org/ S; 2011. [12] Degraeve Z, Jans R. A new Dantzig-Wolfe reformulation and branch-andprice algorithm for the capacitated lot-sizing problem with setup times. Operations Research 2007;55(5):909–20. [13] Fisher ML. The Lagrangian relaxation method for solving integer programming problems. Management Science 1981;27:1–18. [14] Hwang H-C, van den Heuvel W. Economic lot-sizing problem with bounded inventory and lost-sales. Econometric Institute Report EI 2009-01, Erasmus University Rotterdam, The Netherlands; 2009. [15] Held M, Wolfe P, Crowder HD. Validation of subgradient optimization. Mathematical Programming 1974;6:62–88. [16] IMB ILOG CPLEX 12.1. Callable Library. IBM ILOG S.A. /http://www.ilog.com S; 2011. [17] Jans R, Degraeve Z. Meta-heuristics for dynamic lot sizing: a review and comparison of solution approaches. European Journal of Operational Research 2007;177:1855–75. [18] Jans R, Degraeve Z. Modeling industrial lot sizing problems: a review. International Journal of Production Research 2008;46(6):1619–43. [19] Krarup J, Bilde O. Plant location, set covering and economic lot sizes: an O(mn)-algorithm for structured problems. In: Collatz L, et al., editors. Optimierung bei graphentheoretischen and ganzzhligen Problemen. Basel: Birkhauser Verlag; 1977. p. 155–80. [20] Miller AJ, Nemhauser GL, Savelsbergh MWP. On the polyhedral structure of a multi-item production planning model with setup times. Mathematical Programming 2003;94:375–405. [21] Muller LF, Spoorendonk S, Pisinger D. A hybrid adaptive large neighborhood search heuristic for lot-sizing with setup times. European Journal of Operational Research 2012;218(3):614–23. [22] Røpke S, Pisinger D. An adaptive large neighborhood search heuristic for the pickup and delivery problem with time windows. Transportation Science 2006;40:455–72. [23] Pisinger D, Røpke S. Large neighborhood search. Handbook of Metaheuristics— International Series in Operations Research & Management Science 2010;146:399–419. [24] Pochet Y, Wolsey LA. Solving multi-item lot-sizing problems using strong cutting planes. Management Science 1991;37:53–67. [25] Sandbothe RA, Thompson GL. A forward algorithm for the capacitated lot size model with stockouts. Operations Research 1990;38(3):474–86. [26] Trigeiro W, Thomas LJ, McLain JO. Capacitated lot-sizing with setup times. Management Science 1989;35:353–66.