The Economic Lot Scheduling Problem under Performance Decay

Department of Chemical Engineering, University of Sao Paulo, Sao Paulo SP 05508-900, Brazil, Othmer Department of Chemical and Biological Sciences and...
0 downloads 0 Views 170KB Size
Ind. Eng. Chem. Res. 2004, 43, 6463-6475

6463

The Economic Lot Scheduling Problem under Performance Decay Alessandro Alle,† Jose M. Pinto,*,†,‡ and Lazaros G. Papageorgiou§ Department of Chemical Engineering, University of Sao Paulo, Sao Paulo SP 05508-900, Brazil, Othmer Department of Chemical and Biological Sciences and Engineering, Polytechnic University, Six Metrotech Center, Brooklyn, New York 11021, and Centre for Process Systems Engineering, Department of Chemical Engineering, University College London, Torrington Place, London WC1E 7JE, U.K.

The aim of this work is to propose a mathematical programming model for the economic lot scheduling problem (ELSP) with performance decay. First, the problem is formulated as a mixed integer nonlinear program (MINLP), which is found to be nonconvex. The model is then transformed into a mixed integer linear programming (MILP) model through the discretization of the cycle time. This model is tested in a wide set of randomly generated problems with various degrees of difficulty. The results show that the MILP model can achieve the optimum within satisfactory CPU times. An illustrative example demonstrates the applicability of the MILP model and its potential benefits in comparison with a hierarchical approach. 1. Introduction The economic lot scheduling problem (ELSP) is a wellknown problem in the operations research literature. It can be described in the following way: Given a set of products that have continuous fixed demands over an infinite planning horizon, schedule them in a cyclic way in the plant such that the overall production cost, comprising setup and holding costs (see Figure 1), is a minimum. Kuik et al.1 and Drexl and Kimms2 presented extensive literature reviews of the ELSP and other similar problems. Wolsey3 surveyed work that can be used to improve lot scheduling formulations either by a priori reformulation or by the addition of valid inequalities. Previous works consider basically two approaches for determining the product cycle time:4 (1) The common cycle (CC) (Hanssman5) approach assumes that all products have the same cycle time (see Figure 2). This approach has the advantage of always finding a feasible schedule, and it consists of a simplification of the more comprehensive problem where every product has its own cycle time. (2) The basic period (BP) approach (Bomberger,6 Elmaghraby7) allows different cycle times for different products, but restricts each product’s cycle time to be an integer multiple, k, of a basic period of time w (see Figure 3). All lots of an item are of the same size. In general, this approach generates better solutions than the CC approach. However, its main drawback is in ensuring feasibility. Jones and Inman8 have shown that very simple CC scheduling approaches produce optimal or near optimal schedules in many realistic situations. Sahinidis and Grossmann4 proposed a mixed integer nonlinear program (MINLP) model to solve an ELSP in multiple lines with sequence-dependent setups using the CC approach. Pinto and Grossmann9 studied the case of the common cycle schedule for multistage plants with intermediate inventory. * To whom correspondence should be addressed. E-mail: [email protected]. † University of Sao Paulo. ‡ Polytechnic University. § University College London.

Figure 1. Tradeoffs in the economic lot scheduling problem (ELSP).

Figure 2. Common cycle approach.

Figure 3. Basic period approach.

With respect to basic period approaches, Grznar and Riggle10 developed an algorithm that finds global optimal solutions, systematically reducing the (k, w) integer space until an exhaustive search can be done. Gonc¸ alves and Leachman11 presented a hybrid heuristic and a linear programming approach to solve lot scheduling

10.1021/ie049956z CCC: $27.50 © 2004 American Chemical Society Published on Web 09/01/2004

6464

Ind. Eng. Chem. Res., Vol. 43, No. 20, 2004

problems in a continuous time domain. Meyr12 introduced an algorithm that combines a local search with continuous dual optimization to solve the ELSP with sequence-dependent setups in continuous time. Ouenniche and Boctor13 adapted the powers-of-two heuristic14,15 to solve an ELSP that assumes that cycle times of all products are multiples of powers of 2 of a basic period. Yao and Elmaghraby16 presented an analysis of the ELSP without capacity constraints using the powers-of-two heuristic. Recently, Oh and Karimi17 proposed a mathematical programming model as well as a heuristic solution procedure for solving the ELSP with a fixed planning horizon in a single machine with sequence-dependent setups. Ben-Daya18 dealt with an integrated model for the joint determination of the economic lot size and preventive maintenance level for an imperfect process for a single item having a general deterioration distribution with increasing hazard rate. Many processes in the industrial environment are subject to performance decay. For instance, in the chemical processing industry, typical cases are catalyst deactivation in FCC reactors, coking of ethylene furnaces, fouling of heat-exchangers, etc. These units require periodic maintenance to have their performance restored. The problem of stopping these units can be posed as an optimization problem, because there is a tradeoff between allowing units either to operate at lower yield rates or stopping them for maintenance to restore performance. There is an optimal campaign time for which the production losses due to yield decay are balanced by the production gains for maintaining unit operation. To our knowledge, there is no work in the literature that considers ELSP in units subject to performance decay. Jain and Grossmann19 developed an MINLP model for scheduling the cyclic operation of multiproduct parallel units with performance decay. However, this model did not include inventory costs. This same gap exists in the MINLP model proposed by Alle et al.,20 which is aimed at scheduling production and cleaning operations in multiproduct serial plants. These works cannot be strictly classified as studies of the ELSP because they lack one of its main features, which is the consideration of holding costs. The aim of this work is to propose a mathematical programming model for the ELSP with performance decay and to study efficient methods for its solution. The main difficulties related to this new model are the nonlinear and nonconvex terms associated with the dynamic behavior of the system with cycle time. This paper is structured as follows: The problem of economic lot scheduling under performance decay is formally stated in section 2. Section 3 presents an MINLP model for the ELSP which is then linearized and simplified in several aspects, giving rise to a mixed integer linear programming (MILP) model. A study on the computational performance of the MILP model is described in section 4. Section 5 presents two examples of the ELSP with performance decay that are solved with the MILP model. The first example shows the impact of the variations of inventory costs on the optimal schedule, whereas the second shows the benefits of the proposed model when compared with a hierarchical scheduling procedure. Finally, section 6 summarizes the conclusions of the work.

Figure 4. Yield of product i and total amount produced in a stage under decaying performance.

2. Problem Statement The plant processes NP products in the same order in a single stage (permutation flowshop plant). For each product i, a fixed demand rate, di, should be satisfied. Only one run of each product is allowed per cycle (CC approach). The process yield ηi is assumed to decay linearly with time as follows

ηi ) ai - bit

(1)

where ai represents the yield at the beginning of the processing of product i and bi is the linear decay coefficient for product i. The choice of a linear decay function is due to its simplicity and frequent occurrence. Figure 4 shows the yield and amount of product i produced under linearly decaying performance given by ηi ) 1 - 0.02t. As performance decreases, the processing unit requires maintenance to restore the productivity to its initial state. The decision to stop for maintenance interferes with the production schedule. Therefore, it is important to consider both production and cleaning scheduling decisions simultaneously. The overall objective of the planning model is to minimize the overall cost over a given cycle time. The overall cost is given by the sum of the raw material, cleaning, and holding costs. As this work adopts the CC approach, cleaning operations take place during the setup of the unit for processing a different product campaign. Therefore, henceforth, cleaning, setup, and changeover are taken as synonyms. Setup times and costs are incurred whenever the line is changed from one product to another. As transition times and costs are sequencedependent, there are tradeoffs concerning the product sequence. A sequence that minimizes the total transition costs might not be the one that maximizes the availability of the plant because transition times and costs are not necessarily proportional. Simultaneous cleaning and production scheduling involves tradeoffs between raw material consumption and cleaning and inventory costs. A policy of frequent stops might be adopted to minimize raw material consumption because the unit operates at higher yields. An added advantage is that shorter runs mean lower inventory costs. On the other hand, if the objective is to minimize cleaning costs, longer cycles are necessary, which implies larger inventory costs. In addition, raw

Ind. Eng. Chem. Res., Vol. 43, No. 20, 2004 6465

Invfi(t) ) Wi - dit if TPi < t e Tc

(7)

Substituting eq 1 into eq 6 and eq 4 into eq 7, it follows that

1 Invfi(t) ) (aiGi - di)t - biGit2 if 0 e t e TPi (8) 2 Invfi(t) ) di(Tc - t) if TPi < t e Tc

Figure 5. Final inventory during a cycle under performance decay.

(9)

Thus, the overall inventory cost, CH, is given by material consumption increases because of operation at lower yields. Therefore, it is evident that cleaning decisions should be considered in the ELSP with performance decay. Overall, the performance-decay ELSP can be stated as follows: (i) Given: demand rates (di); decay functions [ηi(t)]; raw material (Cfi), setup (Ctrij), and inventory (Cinvfi,) costs; setup times (τij); and feeding rates (Gi) (ii) Determine: the product sequence (Zij), start times (TSi), processing times (TPi), cycle time (Tc), and product amounts (Wi) (iii) So as to minimize: the overall cost (OC) over the cycle time (Tc).

CH )

[

1

1

∑i Cinvfi 2(aiGi - di)TPi2 - 6biGiTPi3 +

As performance decays, the required consumption of raw material for producing a given amount of product is larger. Therefore, the raw material cost should be incorporated into the objective function. Let Cfi be the unit cost of raw material. Then, the overall cost of raw material consumed, CF, is written as

CF ) 3. Modeling of the ELSP with Linear Performance Decay This section presents a mathematical programming model for the problem described in the previous section. In section 3.1, the problem is formulated as an MINLP that is found to be nonconvex. Then, section 3.2 presents a reformulated MILP through the discretization of cycle time. 3.1. MINLP Formulation. The amount of product i produced in the cycle, Wi, is given by the integral of the product between the yield, ηi, and the feed rate of raw material for product i, Gi, along the processing time, TPi

Wi )

∫0TP ηi(t)Gi dt i

(3)

Therefore, the total amount of product i, Wi, is a quadratic function of the processing time, TPi. The amount of product i produced must satisfy the continuous demand, di, over the cycle time, Tc

Wi ) diTc

(4)

The inventory profile of product i during the cycle is given as in Figure 5. The area below the curve in Figure 5 is proportional to the cost of holding product i during the cycle, chi, given by

chi ) Cinvfi

∫0

Tc

Invfi(t) dt

(5)

where

Invfi(t) )

∫0tηi(θ)Gi dθ - dit

if 0 e t e TPi

∑i CfiGiTPi

(11)

Setup costs are sequence-dependent; therefore, the overall setup cost, CT, is given by

CT )

∑i ∑j CtrijZij

(12)

where the variable Zij defines the product ordering in the following way (TSP-like formulation, see Alle and Pinto21)

Zij )

{

1 if product i immediately precedes product j 0 otherwise

(2)

Upon substitution of the yield function given in eq 1 into eq 2, the amount of product i produced in this cycle is given as follows

1 Wi ) aiGiTPi - biGiTPi2 2

]

1 di(Tc - TPi)2 (10) 2

(6)

Thus, in the CC approach, the production sequence can be determined by

∑i Zij ) 1

∀j

(13)

∑j Zij ) 1

∀i

(14)

Equation 13 states that exactly one product i precedes a given product j, whereas eq 14 states that only one product j follows a given product i. Storage capacity might be finite. In this case, it is necessary to determine the maximum amount of each product to be stored during the cycle. The maximum inventory level, Imaxi, reached during the production run of product i is given by

1 Imaxi ) (aiGi - di)TLi - biGiTLi2 2

∀i (15)

where TLi denotes the time instant at which the inventory profile reaches the maximum quantity stored of product i. This maximum inventory level must not exceed the maximum storage limit, Imaxup i . Thus

6466

Ind. Eng. Chem. Res., Vol. 43, No. 20, 2004

leads to

∀i

Imaxi ) IMi

Case B:

(23)

To determine the maximum inventory level, Imaxi, a new set of binary variables, Yi, is introduced to indicate whether case A holds, together with the following constraints up -TPup i Yi +  e TPi - tni e TPi (1 - Yi)

Imaxi g IMi(1 - Yi)

∀i (24)

∀i

Imaxi g di(Tc - TPi) - Imaxup i (1 - Yi)

(25) ∀i (26)

where TPup i is an upper bound for the processing time and  is a small number. It is straightforward to verify that, if Yi ) 1, then constraint 24 becomes Figure 6. Maximum inventory levels.

Imaxi e Imaxup i

-TPup i +  e TPi - tni e 0 w tni g TPi ∀i

(16)

Constraint 15 introduces one more set of nonlinear equations and a new set of variables, TLi, into the model. However, it is possible to determine the maximum amount stored using only linear inequalities. As shown in Figure 6, there are two cases for the final inventory profile: In case A, the maximum inventory level occurs at the end of the production of i, TPi. After that, the inventory is depleted at the continuous demand rate di. In case B, the maximum inventory level Imaxi occurs at tni, which denotes the point of null derivative of inventory as a function of time

tni )

ai di bi biGi

∀i

(17)

Note that eq 17 holds only if bi * 0. Otherwise, performance decay would not occur, and the maximum inventory, Imaxi, would be placed at TPi. Moreover, tni is a constant parameter because it depends only on the decay parameters and the processing rate of each product. From Figure 6, it follows that

Case A:

TPi e tni w TLi ) TPi

Case B:

TPi > tni w TLi ) tni

∀i ∀i

Case A:

Imaxi ) di(Tc - TPi)

∀i

1 Imaxi ) (aiGi - di)tni - biGitni2 2

(20)

∀i (21)

The value of Imaxi given by eq 21 is constant because tni is a constant parameter. Defining parameter IMi as follows

1 IMi ) (aiGi - di)tni - biGitni2 2

∀i

0 +  e TPi - tni e TPup i w tni < TPi

(22)

∀i

which defines case B. Together, constraints 25 and 26 yield Imaxi g IMi and Imaxi g di(Tc - TPi) - Imaxup i , respectively. The latter becomes redundant because its right side is negative; thus Imaxi ) IMi. Therefore, constraints 25 and 26 sufficiently represent the maximum inventory levels, so nonlinear eq 15 can be removed. The cycle time must be larger than the sum of the processing and transition times

Tc g

∑i TPi + ∑i ∑j τijZij

(27)

The model should provide optimal starting times of each process at the unit. Because scheduling is cyclic, product 1 is arbitrarily chosen as the first in the sequence

TS1 )

(19)

In case B, the maximum inventory Imaxi occurs at tni

Case B:

which defines case A. Together, constraints 25 and 26 give Imaxi g 0 and Imaxi g di(Tc - TPi), respectively. Otherwise, if Yi ) 0, then

(18)

The maximum inventory in case A is placed at TPi. Then, the maximum inventory is depleted at the continuous demand rate di. Thus

∀i

∑j τj1Zj1

(28)

Each product starts to be processed only after the previous processing and setup are completed

TSj )

∑i Zij(TSi + TPi + τij)

j ) 2, ..., NP (29)

Equation 29 can be linearized through the use of the following constraint21

-Tcup(1 - Zij) e TSj - (TSi + TPi + τij) e Tcup(1 - Zij)

∀i * j, j ) 2, ..., NP (30)

It must be pointed that constraint 30 prevents solutions that present subcycles (Alle and Pinto21). The objective function to be minimized consists of the sum of the cleaning, holding, and raw material costs, divided by the cycle time. Thus, the ELSP model with performance decay can be stated as follows:

Ind. Eng. Chem. Res., Vol. 43, No. 20, 2004 6467

Model M1

and CT

min OC )

CT + CH + CF Tc

(31) OC )

subject to

CT + B2

-

Tc

∀j

∑j Zij ) 1

∀i

2 3

B3i ) Fi + Bi,

∑i ∑j CtrijZij ∑i

CH )

CF )

∑i Dix1 - CiTc - D + ETc

[

1 1 Cinvfi (aiGi - di)TPi2 - biGiTPi3 + 2 6 1 di(Tc - TPi)2 2

∑i CfiGiTPi

D)

]

Di )

∑i Di,

B2 )

E)

aidiCinvfi , bi

1 Wi ) aiGiTPi - biGiTPi2 2

ai Fi ) Cfi Gi bi

(33)

Xl ) 1 ∑ l)1

(34)

where T ˆ cl represents alternative discrete values for Tc such that

∀i

up -TPup i Yi +  e TPi - tni e TPi (1 - Yi)

∀i

ˆ cl < T ˆ cl+1 e Tcup Tclo e T

∀i

Imaxi g di(Tc - TPi) - Imaxup i (1 - Yi)

∀i

Cycle time constraint

∑i TPi + ∑i ∑j τijZij

∀l

CTloXl e CTXl e CTupXl

∑j τj1Zj1

CT )

-(1 - Zij)Tcup eTSj - (TSi + TPi + τij) e (1 - Zij)Tcup

l ) 1, ..., L - 1 (35)

The derivation of tight lower and upper bounds for the cycle time is given in Appendix C. The first term on the right-hand side of eq 32 presents the ratio of Tc and CT. It can be linearized by defining the variable CTXl as the product of the continuous variable CT and the binary variable Xl. Hence

LHS of the constraint TS1 )

diCinvfi , 2

L

Maximum inventory constraints

Tc g

3bi2

T ˆ clXl ∑ l)1

Tc )

∀i

Imaxi g IMi(1 - Yi)

Bi )

ai3GiCinvfi

L

∀i

Demand satisfaction constraint

Imaxi e Imaxup i

∑i Ei,

Ei )

bidi Ci ) 2 ai2Gi

∑i B3i,

Note that eq 32 is nonconvex. One step toward making it convex is to discretize the cycle time Tc into L values as follows

Mass balance constraint

Wi ) diTc

(32)

where

Cost definition constraints CT )

+

Tc

Sequence constraints

∑i Zij ) 1

∑i B3ix1 - CiTc

∀i * j, j ) 2, ..., NP

(37)

The objective function can then be rewritten as

Integrality constraints Yi, Zij ∈ {0,1} Model M1 corresponds to a MINLP model with a nonconvex (pseudoconvex) objective function and nonconvex constraints 3 and 10. 3.2. Reformulated MILP. First, objective function 31 can be rewritten after proper substitutions (see Appendix A). Moreover, mass balance and inventory satisfaction constraints as well as inventory and raw material cost definitions can be eliminated (see Appendices A and B). Consequently, the objective function can be expressed as a function of only the variables Tc

∑l CTXl

(36)

min OC )

∑l

(

CTXl T ˆ cl

)

+ RlXl - D

(38)

where

Rl )

B2 T ˆ cl

-

∑i

B3ix1 - CiT ˆ cl T ˆ cl

+ 2 3

∑i Dix1 - CiTˆ cl + ETˆ cl

6468

Ind. Eng. Chem. Res., Vol. 43, No. 20, 2004

The variable TPi can also be expressed in discrete values as

Integrality constraints Zij, Xl ∈ {0, 1}

L

TPi )

TP ˆilXl ∑ l)1

∀i

(39)

ˆil is obtained from T ˆ cl (see eq A.3) in Apwhere TP pendix A) as

ˆil ) TP

ai (1 - x1 - CiT ˆ cl) bi

∀i,l

The model can be further simplified from the analysis of inventory constraints. These constraints can be eliminated by the use of bounds for the processing times that consider storage capacity limitations (see Appendix B). After the simplifications above, model M1 for ELSP with performance decay has been transformed into a MILP model, named M2:

Model M2 min OC )

(

∑l

CTXl T ˆ cl

)

+ RlXl - D

subject to

Discrete cycle time constraints L

Xl ) 1 ∑ l)1 L

Tc )

T ˆ clXl ∑ l)1 L

TP ˆilXl ∑ l)1

TPi )

∀i

Cost definition constraints CT )

∑i ∑j CtrijZij

CTloXl e CTXl e CTupXl CT )

∀l

∑l CTXl

Sequence constraints

∑i Zij ) 1

∀j

∑j Zij ) 1

∀i

Cycle time constraint Tc g

∑i TPi + ∑i ∑j τijZij

Timing constraints TS1 )

∑j τj1Zj1

-Tcup(1 - Zij) e TSj - (TSi + TPi + τij) e Tcup(1 - Zij)

∀i * j, j ) 2, ..., NP

4. Results In this section, model M2 is tested in a set of randomly generated problems. To obtain the reference optimal solution of the problem, a benchmark is set by solving each instance for L ) 100. The deviation from the benchmark optimum, OD, is given by

OD )

f - fL)100 fL)100

(40)

where fL)100 is the objective function value found in the benchmark solution and f is the value obtained by solving model M2 with L ) 10 or L ) 5 discrete values for Tc. The computational effort is measured as the CPU time required to solve each instance of the problem. A PC Pentium III 500 MHz computer with 512 Mb of RAM was used to perform the computational tasks. Model M2 was implemented with GAMS (Brooke et al.22). The MILP solver was CPLEX 7.0 (Ilog23). 4.1. Data Generation. This paper relies on the approach of Dobson15 that is used by Oh and Karimi17 to generate problems with random data. Dobson15 used three parameters to represent the degree of difficulty associated with the randomly simulated problems. These are the machine utilization level u, the product diversity factor r, and the sequence dependence diversity factor d. The utilization level u denotes the fraction of the total cycle time Tc that is actually used for production. Because the overall machine time is spent in either production or setup or idleness, (1 - u) indicates the extent of the setup and idle times. The product diversity factor r describes the extent of variation in setup times; decay parameters; maximum inventory capacities; and setup, raw material, and inventory costs across products, whereas the sequence dependence diversity factor d describes the range of variation in sequence-dependent setup costs and times. The values of r and d are assigned a priori for each problem and then used to generate the problem data. The values used in our experiments are from Oh and Karimi:17 r ∈ {3, 5} and d ∈ {3, 5}. To describe the problem of data generation precisely, let U[a, b] represent a uniformly distributed random variable in the range [a, b]. The first value selected is di ) U[15, 25] (ton/day). Then, minimum values for the transition times and costs, raw material and inventory costs, and linear decay coefficientsτmin, Ctrmin, Cinvfmin, Cfmin, and bmin, respectivelysare chosen. Sequence-independent parameters are chosen as follows: Cinvfi ) CinvfminU[1, r], Cfi ) CfminU[1, r], and bmin ) CfminU[1, r]. A base changeover cost, Ctri ) CtrminU[1, r], and a base changeover time, τi ) τminU[1, r], are chosen for each product i. Using these base setup times and costs, the actual sequence-dependent setup costs and times are generated: Ctrij ) CtriU[1, r] and τij ) τiU[1, r]. To facilitate the generation procedure, ai ) 1 was stated for each product i. The bi parameter is generated in the same way as the cost parameters with bmin ) 0.01. This value was selected such that the performance decay had a sensible influence on the process yield, considering around 7 days for product campaigns.

Ind. Eng. Chem. Res., Vol. 43, No. 20, 2004 6469

Maximum inventory capacities were randomly generated as follows

Imaxup i ) IMiU[0.5, 1.5]

(41)

where IMi is the null derivative value of the inventory function Invfi(t), given in eq 22. Dobson,15 working with sequence-dependent ELSPs, asserts that it is difficult to generate a problem with an exact value of machine utilization u because this parameter is a function of randomly generated demands and production rates. In the case of constant yields, utilization is given in the following way

di

∑i η G

u)

i

(42) i

To overcome this problem, the author used small ranges of u, [umin, umax], rather than exact values. For problems with performance decay, there is an additional difficulty of time-varying yields that prevents the utilization from being estimated precisely prior to problem solution. Therefore, it is satisfactory to make use of a rough estimate of u0 in data generation, which can be obtained from eq 42

u0 )

di

∑i a G i

(43) i

where u0 is a lower bound for the utilization in the case of performance decay. The actual utilization is given by

u)

∑i TPi

(44)

Tc

and might be significantly different from the estimation obtained with eq43. The idleness s is defined as the time fraction during which the unit remains rigorously idle, i.e., without cleaning or production

(

s) 1-

)

∑i TPi + ∑i ∑j Zijτij Tc

× 100%

(45)

The idleness measures the slack that exists in the production capacity of the plant. The procedure used to generate Gi values for a given number NP of products, compatible with a given estimated utilization, u0, is as follows: (1) State the machine utilization u0 in a range [u0min, u0max]. (2) Take a random value u j 0 that lies in this range. (3) For every i, calculate

Gi )

U[15, 25] NP aiu j0

where NP is the number of products. Note that the numerator has the same uniform probability distribution as was used to generate di. e ∑idi/(aiGi) e umax (4) Verify that umin 0 0 . If not, go to step 2.

This procedure has proven to be very efficient for the generation of Gi values for a given problem size of NP and a range of u0. In most cases, step 2 is performed just once. The second parameter, the actual u, is calculated from eq 44 with the optimization results and compared with u0 in the tests presented in the next sections. The CC approach always yields a feasible schedule for the common ELSP when the utilization given in eq 43 is lower than 1 (Jones and Inman8). However, the ELSP with performance decay does not guarantee a feasible solution to the problem even when u0 is less than 1 because this parameter represents an underestimate of the true u value required to satisfy the demand. The closer u0 gets to 1.0, the less likely the problem is to be feasible. Checking whether the bounds (see Appendix C) are coherent (i.e., each lower bound is smaller than the respective upper bound) is an easy and straightforward procedure for detecting infeasible problems. Model M2 was tested on a number of randomly generated problems obtained with variable values of r and d, as well as with the variation of the number of products and ranges of u0. For all problems, the minimum setup time (τmin) is 0.4 days. This value was chosen so that changeovers are limited to less than 20% of the total processing time, given that the average product campaign duration was a priori set to last around 7 days. Minimum cost parameters are Cinvfmin ) 0.003 $‚day-1‚ton-1, Cfmin ) 0.3 $‚day-1, and Ctrmin ) 800 $. These values were chosen so that all cost parcels make significant contributions to the overall production costs. The computational experiments were divided into two sections. The goal of section 4.2 is to evaluate the impacts of parameters r, d, and u0 on the performance of the scheduling models. The expectation is that the more diversified the data are, the more difficult the problem is to solve. Section 4.3 aims at showing how an increase in the number of products influences the computational effort required to solve the problem. 4.2. Influence of Data Variability and Machine Utilization. This section analyzes the influence of data variability and estimated utilization, u0, on the computational performance of the models. Overall, 10 problems were randomly generated for NP ) 20 at every range of u0. The results in Table 1 refer to the average CPU times and solution values for these problems. The actual utilization u and idleness s result from the application of constraints 44 and 45, respectively, to the optimal times obtained from the model benchmark solution. 4.3. Influence of Problem Size. The results in Table 2 show the effect of the number of products to be scheduled on the performance of the proposed MILP model in terms of CPU solution times and deviations from the benchmark optimal solution (OD). The data generation parameters were selected as follows: u0 ) [0.74, 0.75], r ) 3, and d ) 5. Ten problems were generated for each row of the table, which shows the average CPU times and ODs for these problems. As can be seen in Table 2, the CPU times required for model M2 increase exponentially with problem size. The relative measure OD, although small, seems to increase with NP as well. 5. Examples In this section, the benefits of the MILP model are demonstrated through the analysis of two case studies.

6470

Ind. Eng. Chem. Res., Vol. 43, No. 20, 2004

Table 1. Influence of Data Variability and Utilization for Model M2a L)5

L ) 10

L ) 100

r

d

u0

u

s (%)

CPU (s)

σCPU (s)

OD (%)

σOD (%)

CPU (s)

σCPU (s)

OD (%)

σOD (%)

CPU(s)

σCPU (s)

3

3

0.79-0.80 0.74-0.75 0.69-0.70 0.64-0.65

0.867 0.802 0.748 0.689

0.03 1.21 6.09 12.5

5.6 2.4 2.8 1.5

3.8 0.6 2.4 1.0

0.27 0.37 0.47 0.33

0.305 0.430 0.445 0.380

8.6 5.5 2.5 1.9

5.5 3.5 1.5 0.8

0.15 0.07 0.02 0.11

0.155 0.074 0.021 0.154

52.9 32.9 17.6 17.1

37.2 22.3 10.7 7.8

3

5

0.74-0.75 0.69-0.70 0.64-0.65 0.59-0.60

0.820 0.751 0.695 0.635

0.10 1.01 3.98 9.90

10.3 3.2 3.0 2.8

4.7 2.8 2.3 1.9

0.49 0.37 0.33 0.67

0.292 0.444 0.280 0.538

8.1 4.8 3.1 2.9

4.6 3.0 1.5 1.3

0.17 0.15 0.05 0.15

0.149 0.122 0.047 0.132

59.4 25.4 27.4 22.1

44.0 12.6 14.6 12.9

5

5

0.69-0.70 0.64-0.65 0.59-0.60 0.54-0.55

0.795 0.722 0.659 0.599

0.14 0.94 1.39 3.15

7.1 5.1 5.3 2.8

5.6 4.5 3.2 1.6

0.14 0.04 0.34 0.16

0.289 0.107 0.301 0.249

9.5 9.3 5.8 3.2

8.4 7.4 3.7 2.1

0.06 0.06 0.07 0.03

0.132 0.121 0.090 0.056

67.4 35.4 31.4 25.8

113.5 46.2 22.3 19.9

a Note in Table 1 that the more diversified the problem, the higher the actual utilization u for a given range of u . As parameters r and 0 d increase, u0 must be reduced to generate feasible problems, i.e., problems where u is less than 1. This is the reason for different u0 ranges in Table 1. Note also that the values of OD are significantly small. These results imply that 10 discrete values (L ) 10) are enough to obtain very accurate optimal schedules within satisfactory CPU times. As expected, the CPU time increases with data variability and decreases with the slack in plant capacity, s.

Table 2. Effect of Problem Size in the Computational Performance of the Models L)5

L ) 10

NP

u

s (%)

CPU (s)

OD (%)

CPU (s)

OD (%)

15 20 25 30

0.83 0.82 0.79 0.79

0.13 0.10 0.00 0.00

2.0 8.1 20.8 61.5

0.41 0.49 0.83 0.82

1.7 10.3 37.2 99.4

0.37 0.17 0.54 0.56

In the first one, it is shown how model M2 schedules the optimal inventory policies for a unit with performance decay with different inventory costs. In the second case study, the effectiveness of the MILP model is shown through a comparison of the optimal schedule with another obtained from a procedure based on the hierarchical decomposition of the problem into transition cost minimization and lot-sizing subproblems. 5.1. Relation between Inventory Costs and Optimal Inventory Policies. Inventory cost is the sum of several costs related to the maintenance of storage, the capital depreciation cost calculated from product prices and a given interest rate, the storage rental costs, and the energy costs (if products must be kept under specific conditions, for example, frozen goods). Therefore, the uncertainty in this parameter can be high. Next, the impact of the inventory cost on the optimal production scheduling and inventory policy is illustrated. Consider the problem of scheduling economic lots of three products (A-C) in a plant subject to performance decay. Table 3 provides the basic problem data. To study the impact of inventory costs under the optimal schedule, two scenarios are proposed as follows: For scenario I, Cinvfi ) 0.04 $/(ton‚day). For scenario II, Cinvfi ) 0.01 $/(ton‚day). Table 4 shows the costs and the objective function values for the optimal schedules obtained by model M2 (L ) 10) applied to both scenarios. Because of its higher holding costs, scenario I has larger inventory costs per unit cycle time than scenario II, which has a lower (37%)

objective function value. Figure 7 shows the inventory optimal profiles obtained by model M2 for scenarios I and II. Note that the maximum storage capacity is respected for every product. In scenario I, the inventory costs (Cinvfi) are larger. Note that product B follows case A, shown in Figure 6. In scenario II, the inventory costs are smaller. Thus, the cycle time is longer (90 vs 60 days in scenario I). Consequently, the areas under the curves are larger. Note that product B reaches its maximum level before the end of its processing, as in case B in Figure 6. 5.2. Mathematical Programming vs Hierarchical Procedure. This illustrative example shows a comparison between the proposed MILP and a hierarchical procedure for scheduling a unit that processes 20 products. The hierarchical approach comprises two steps: In the first step, the product sequence is obtained by minimizing the overall setup cost. The second step determines the cycle and processing times of every product for this sequence. The hierarchical procedure simulates the practical approach by schedulers of using heuristic rules and problem decomposition to solve schedules more rapidly. The idea is to first determine the sequence using the changeover times that depend on the sequence and then solve the timing problem for a fixed sequence. Problem data were randomly generated as shown in Table 5, where U[a, b] represents a uniformly distributed random variable in the range [a, b]. Model M2 was implemented in GAMS and solved in 9.4 s on a PC Pentium III 500 MHz computer with 512 Mb RAM with the MILP solver CPLEX 7.0. The model contains 868 constraints and 502 variables, of which 420 are binary. Table 6 and Figure 8 compare the schedules obtained from the two approaches. Figure 8 shows that each approach yields a different production sequence. The hierarchical approach has the minimum setup cost, whereas the schedule obtained from model M2 has a shorter cycle time (90 vs 225 days).

Table 3. Basic Problem Data Ctrij ($) i 1 2 3

bi

(day-1)

0.0100 0.0080 0.0095

di (ton/day)

Gi (ton/day)

1.20 12.50 1.00

20 22 18

Imaxup i

(ton)

150 300 150

τij (days)

Cfi ($/ton)

j)1

2

3

j)1

2

3

1.0 1.0 1.0

0 400 400

600 0 600

320 320 0

0 2 2

3.5 0 3.5

1.8 1.8 0

Ind. Eng. Chem. Res., Vol. 43, No. 20, 2004 6471

150% lower inventory costs, which yields a 20.0% lower overall cost per unit of cycle time. The same example is run for the case in which performance decay is negligible. In other words, product yield values do not decrease over time. The implementation was done by setting the values of parameter bi, i ) 1, ..., 20, to 5 × 10-7 day-1. The objective function value for this case is OC ) 101.5 × 103 $/h, which is 3.6% lower than the optimal solution reported in Table 6 and, as expected, underestimates the operating costs. Interestingly, the cycle time is the same as that obtained for the case of performance decay and shown in Figure 8. 6. Conclusions

Figure 7. Optimal inventory policies for scenarios I and II in the case study. Table 4. Costs in the Schedules Obtained with Model M2 for Scenarios I and II costs (103 $/h)

scenario I

scenario II

CF/Tc CT/Tc CH/Tc OC

17.2 22.1 19.8 59.1

19.7 14.7 2.7 37.1

Table 5. Parameters for the Illustrative Example parameter

expression

parameter

expression

di Cinvfi ai Ctrij τij Imaxi

U[15, 25] ton/day 0.01U[1, 3] ton/day 1 CtriU[1, 5] $ τiU[1, 5] $ IMiU[0.5, 1.5] ton

Gi Cfi bi Ctri τi Tcl

U[410, 680] ton/day 100U[1, 3] ton/day 0.02U[1, 3] day-1 3000U[1, 3] $ 0.4U[1, 3] $ 75 + (360 - 75)(l - 1)/19, l ) 1, ..., 20

Table 6. Costs in the Schedules Obtained with the Hierarchical Procedure and with Model M2 costs (103 $/h)

hierarchical

M2

benefı´t (%)

CF/Tc CT/Tc CH/Tc OC

84.4 2.1 45.1 131.7

79.0 8.2 18.1 105.3

6.8 -74.2 149.5 20.0

As a consequence, the MILP approach generates higher process yields as well as lower inventory levels. In fact, Table 6 shows that the MILP has a 74.2% larger setup cost. However, it presents 6.8% lower raw material and

This work introduced an original problem into the operational research area, based on a classical problem (ELSP), in which a recurrent chemical process feature was incorporated: performance decay with processing time. ELSP with performance decay was first formulated as a nonconvex MINLP that was simplified and transformed into an MILP model. This model was tested for a large set of randomly generated ELSPs with various degrees of difficulty measured by data diversification, machine utilization level, and number of products. Our results show that model M2 provides very accurate results with low values of discrete cycle times and relatively short CPU times even for problems with wide data variability and high machine utilization levels. However, the computational effort increases exponentially with the problem size, as expected because of the consideration of sequence-dependent setups, which causes the ELSP to be an NP-hard problem, similarly to the TSP (traveling salesman problem). Finally, case studies show the model’s ability to predict optimal schedules under different scenarios and its advantages in comparison with a hierarchical approach. Future work concerns the global solution of the nonconvex MINLP model, whose objective is to compare its solution to that obtained from the reformulated MILP. Moreover, production rates can be considered as decision variables provided that there are tradeoffs between the rates and yield decay. Finally, the stochastic version of the problem that considers uncertainty in both the yield performance and the product demands must be addressed. Acknowledgment The authors acknowledge support received from FAPESP (Fundacao de Amparo a Pesquisa do Estado de Sao Paulo) under Grants 99/02657-8 and 98/ 14384-3. Appendix A: Elimination of Mass Balance Equations To eliminate the mass balance equations, first, eq 4 is substituted into eq 3

1 aiGiTPi - biGiTPi2 ) diTc 2

(A.1)

Equation A.1 can be solved for TPi as a function of Tc

TPi )

x

ai ai ( bi bi

bidiTc 1-2 2 ai Gi

(A.2)

6472

Ind. Eng. Chem. Res., Vol. 43, No. 20, 2004

Figure 8. Schedules for the (a) hierarchical and (b) MILP approaches.

However, from eq C.1 in Appendix C below, it follows that

TP(1) i )

x

ai ai + bi bi

1-2

bidiTc 2

ai Gi

> TPup i )

ai bi

Therefore, the only meaningful solution for eq A.1 is

x

ai ai TPi ) bi bi TPi )

1-2

bidiTc ai2Gi

diCinvfi , 2

ai Fi ) Cfi Gi bi

Therefore, the holding and raw material cost definition constraints can also be eliminated from the model through the substitution of eqs A.5 and A.6 into eq 31, which yields eq 32. Appendix B: Elimination of Inventory Constraints

w

A sufficient condition that ensures that the inventory bound is never violated is as follows

ai (1 - x1 - CiTc) bi

(A.3)

0 e t e TPi w Imaxup i - Invfi(t) g 0

(B.1)

Substituting eq 8 into inequality B.1 and rearranging, it follows that

where

bidi Ci ) 2 2 ai Gi

0 e t eTPi w Imaxup i +

From A.3, the choice of Tc defines all processing times TPi. Note that eq A.3 holds only if the square root argument is nonnegative. Thus

0 e 1 - CiTc e 1 w 0 < Tc e

1 Ci

∀i

∑i

(

2 Bi - Bix1 - CiTc + DiTcx1 - CiTc 3

)

DiTc + EiTc2 (A.5) CF )

∑i Fi(1 - x1 - CiTc)

where

ai3GiCinvfi , Bi ) 3bi2

Di )

aidiCinvfi bi

(A.6)

bi 2 G t - (aiGi - di)t g0 2 i (B.2)

To obtain the values of t, 0 e t e TPi, that satisfy inequality B.2, a root analysis of this inequality is required. The roots are

(A.4)

Note that Tc cannot be null in eq A.4, because it is the denominator of the objective function in eq 31. Equation A.3 can be directly substituted into eqs 10 and 11 to give

CH )

Ei )

t(1) i )

aiGi - di - x(aiGi - di)2 - 2biGiImaxup i b i Gi

t(2) i )

aiGi - di + x(aiGi - di)2 - 2biGiImaxup i b i Gi

and

(2) It is evident that t(1) i e ti and that the left-hand side of inequality B.2 is a parabola with positive convexity (a positive quadratic term). Thus, Figure 9 shows the behavior of the quadratic function in inequality B.2. From Figure 9, it can be concluded that inequality and t g t(2) and it is B.2 is satisfied when t e t(1) i i (1) (2) violated when ti < t < ti . A sufficient condition for inequality B.2 to never be violated is that the parabola does not intercept the t axis or that intercepts it at a single point, i.e., that the quadratic function does not have real roots or has one single real root. This occurs

Ind. Eng. Chem. Res., Vol. 43, No. 20, 2004 6473

guarantee that the inventory capacity is never violated, TPil values must respect bounds given in inequality B.5. Appendix C: Tight Bounds for the ELSP with Performance Decay

Figure 9. Positive and negative values assumed by the quadratic function in the maximum inventory constraint B.2.

Bound determination is crucial for solving the problem because the tighter the bounds on variables, the better the relaxation of the MILP model. This appendix focuses on the determination of the tightest bounds for the optimization variables. Obviously, the production yield is nonnegative. Therefore, an upper bound for TPi can be obtained from eq 1 as follows

if the argument of the square root in the expressions (2) for t(1) i and ti is nonpositive, which means that

Imaxup i

(aiGi - di)2 g 2biGi

(B.3)

Therefore, for every product i such that the inequality B.3 is satisfied, the maximum inventory constraints (inequality 16) are unnecessary. These products can be identified before the model is solved, and then their corresponding capacity constraints (in inequality16) can be eliminated. The other products can violate constraint 16. These products define the set Ilim, given by

{|

Ilim ) i Imaxup i e

}

(aiGi - di)2 2biGi

(B.4)

In this way, constraints 24-26 can be written only for ∀i ∈ Ilim instead of φοF ∀i. However, even for these products, it is possible to determine the values to be assumed by binary variable Yi before solving the scheduling model, as will be shown next. When i ∈ Ilim, the maximum storage capacity can be violated if the processing time is long enough that TPi > t(1) i , as can be seen in Figure 9. Therefore, a sufficient condition to guarantee the satisfaction of constraint 16 is that TPi e t(1) i for i ∈ Ilim, i.e.

aiGi - di - x(aiGi - di) b i Gi 2

0 e TPi e

2biGiImaxup i ∀i ∈ Ilim (B.5)

From inequality B.5, it follows that

TPi e

aiGi - di , b i Gi

∀i ∈ Ilim

(B.6)

Therefore, it follows from eq 17 that case A occurs for i ∈ Ilim

TPi e tni,

∀i ∈ Ilim

(B.7)

Therefore, case B shown in Figure 6 will never occur for i ∈ Ilim. Consequently, the corresponding binary variables Yi can be fixed at Yi ) 1. Thus, inventory constraints are unnecessary because, for i ∉ Ilim, there is no reason for determining the maximum inventory, given that the storage capacity cannot be violated. For i ∈ Ilim, then necessarily TPi e tni and Yi ) 1, which allow these constraints to be eliminated. However, to

ai - biTPi g 0 w TPup i )

ai bi

∀i

(C.1)

When there is a limitation on the storage capacity, an upper bound of TPi is obtained from inequality B.5

TPup i )

aiGi - di - x(aiGi - di)2 - 2biGiImaxup i b i Gi ∀i ∈ Ilim (C.2)

The cycle time Tc depends on the processing times TPi and vice versa. The following expression for Tc as a function of a given TPi is derived from eq A.3

Tc )

[ (

bi 1 1 - 1 - TPi Ci ai

)] 2

∀i

(C.3)

Equation C.3 expresses Tc as an increasing function of TPi under the conditions imposed by eq C.1, i.e., TPi(bi/ai) e 1. Therefore, it follows that

Tc e

[ (

bi 1 1 - 1 - TPup Ci ai i

)] 2

∀i

(C.4)

Thus, the upper bound for Tc is obtained as

Tcup ) min i

{[ (

bi 1 1 - 1 - TPup Ci ai i

) ]} 2

(C.5)

This bound is as tight as possible because it takes into account the processing time limit due to the performance decay as well as that due to the storage capacity. However, the upper bounds for TPi given by eqs C.1 and C.2 are not the tightest ones because it follows from eq A.3 that

TPup i )

ai (1 - x1 - CiTcup) bi

(C.6)

Therefore, the substitution of Tcup from eq C.5 into eq C.6 can result in TPup i values that are lower than the ones given in eqs C.1 and C.2. Thus eq C.6 must also be used to impose tighter bounds for TPi. Obviously, that is obtained when substituted again every TPup i into eq C.5 yields the same Tcup value previously found because the use of eq C.6 is a direct substitution procedure. Consequently, TPup i does not admit further adjustment, thus being the tighest values for this variable.

6474

Ind. Eng. Chem. Res., Vol. 43, No. 20, 2004

In addition to the tightest upper bounds for Tc and TPi, lower bounds must be determined. The processing time of a given product in the case of performance decay must be higher than the processing time in the case that the initial yield is constant. Consequently

diTclo TPi g a i Gi

On the other hand

Tc g

∀i

{

(C.7)

{

ai (1 bi

τlo

Tclo ) max

1-

∑i x1 - CiTclo)

∀i

(C.8)

Combining expressions C.7 and C.8, it follows that

TPlo i ) max

∑i x

1 - CiTclo),

}

diTclo a i Gi

∀i (C.9)

The lower bound for Tc is obtained from the substitution of inequality C.7 into eq 27

diTclo

(C.14)

Comparing expressions C.13 and C.14, a tighter bound for Tc is found as follows

From eq A.3, it follows that

ai TPi g (1 bi

∑i TPloi + τlo

,

di

∑i a G i

}

∑i TPlo + τlo

i

(C.15)

Note that the bound for Tc depends on the bounds for TPi and vice versa. Therefore, they must be determined iteratively in a loop that must proceed until no further improvement can be achieved. If the number of products to be processed is very large, the determination of valid lower bound for τ, given by eq C.11 subject to constraints 13, 14, and C.12, requires significant computational effort because this problem is a TSP (traveling salesperson problem), well-known to be NP-hard (Lawler24). Nevertheless, there are other ways to determine bounds for τ, such as those given in eqs C.16 and C.17

(C.10)

τlo )

{τij} ∑i min j

(C.16)

where τlo is a lower bound for the total transistion time, given by

τup )

{τij} ∑i max j

(C.17)

Tc g

∑i a G i

min τ )

+ τlo

i

∑i ∑j τijZij

(C.11)

subject to

∑i Zij ) 1

∀j

(13)

∑j Zij ) 1

∀i

(14)

and constraints of subcycle elimination, such as an adaptation of expression 30 with the assignment of positive values to TPi or TPi ) 0, in case all τij are strictly positive up

-(1 - Zij)Tc

e TSj - (TSi + τij) e ∀i * j, j ) 2, ..., NP (C.12)

(1 - Zij)Tcup

Tc )

diTclo

∑i a G i

Tc ) 1-

CTup )

{Ctrij} ∑i max j

(C.19)

Because product 1 is arbitrarily chosen as the first product to be processed in the cyclic schedule, a lower bound for variable TS1 can be obtained as follows

TSlo 1 ) minτi1 i

(C.20)

Thus, lower bounds for the remaining variables TSi, i ) 2, ..., NP, are obtained as lo lo TSlo i ) TS1 + TP1 + τi1

i ) 2, ..., NP

i

+τ w

(C.21)

(C.22)

However, upper bounds for the remaining variables TSi, i ) 2, ..., NP, are obtained as follows

) τlo w

up - TPlo TSup i ) Tc i

i

i ) 2, ..., NP

(C.23)

Literature Cited

τlo

lo

(C.18)

TSup 1 ) maxτi1

di i

{Ctrij} ∑i min j

lo

i

∑i a G

Tclo - Tclo

CTlo )

An upper bound for TS1 is given by

From inequality C.10, it follows that lo

Bounds for the variable CT can be sought in a similar way

(C.13)

di

∑i a G i

i

(1) Kuik, R.; Salomon, M.; Wassenhove, L. N. Batching Decisions: Structure and Models. Eur. J. Oper. Res. 1994, 75, 243. (2) Drexl, A.; Kimms, A. Lot Sizing and SchedulingsSurvey and Extensions. Eur. J. Oper. Res. 1997, 99, 221.

Ind. Eng. Chem. Res., Vol. 43, No. 20, 2004 6475 (3) Wolsey, L. A. MIP Modelling of Changeovers in Production Planning and Scheduling Problems. Eur. J. Oper. Res. 1997, 99, 154. (4) Sahinidis, N. V.; Grossmann, I. E. MINLP Model for Cyclic Multiproduct Scheduling on Continuous Parallel Lines. Comput. Chem. Eng. 1991, 15, 85. (5) Hansmann, F. Operations Research in Production and Inventory; Wiley: New York, 1962. (6) Bomberger, E. A Dynamic Programming Approach to a Lot Size Scheduling Problem. Manage. Sci. 1966, 12, 778. (7) Elmaghraby, S. E. The Economic Lot Scheduling Problem (ELSP): Review and Extensions. Manage. Sci. 1978, 24, 587. (8) Jones, P. C.; Inman, R. P. When is the Economic Lot Scheduling Problem Easy? IIE Trans. 1989, 21, 11. (9) Pinto, J. M.; Grossmann, I. E. Optimal Cyclic Scheduling of Multistage Multiproduct Plants. Comput. Chem. Eng. 1994, 18, 797. (10) Grznar, J.; Riggle, C. An Optimal Algorithm for the Basic Period Approach to the Economic Lot Scheduling Problem. Omega Int. J. Manage. Sci. 1997, 25, 335. (11) Gonc¸ alves, J. F.; Leachman, R. C. A Hybrid Heuristic and Linear Programming Approach to Multi-Product Machine Scheduling. Eur. J. Oper. Res. 1998, 110, 548. (12) Meyr, H. Simultaneous Lotsizing and Scheduling by Combining Local Search with Dual Reoptimization. Eur. J. Oper. Res. 2000, 120, 311. (13) Ouenniche, J.; Boctor, F. F. The Multi-Product, Economic Lot-Sizing Problem in Flow Shops: The Powers-of-Two Heuristic. Comput. Oper. Res. 2001, 28, 1165. (14) Singh, H.; Foster, J. Production Scheduling with Sequence Dependent Setup Costs. IIE Trans. 1987, 19, 43. (15) Dobson, G. The Cyclic Lot Scheduling Problem with Sequence Dependent Setups. Oper. Res. 1992, 40, 736.

(16) Yao, M. J.; Elmaghraby, S. E. The Economic Lot Scheduling Problem under Power-of-two Policy. Comput. Math. Appl. 2001, 41, 1379. (17) Oh, H. C.; Karimi, I. A. Planning Production on a Single Processor with Sequence Dependent Setups. Part I: Determination of Campaigns. Comput. Chem. Eng. 2001, 25, 1021. (18) Ben-Daya, M. The Economic Production Lot-Sizing Problem with Imperfect Production Processes and Imperfect Maintenance. Int. J. Prod. Econ. 2002, 76, 257. (19) Jain, V.; Grossmann, I. E. Cyclic Scheduling of Continuous Parallel Process Units with Decaying Performance. AIChE J. 1998, 44, 1623. (20) Alle, A.; Pinto, J. M.; Papageorgiou, L. G. A Mathematical Programming Approach for Cyclic Production and Cleaning Scheduling of Multistage Continuous Plants. Comput. Chem. Eng. 2004, 28, 3. (21) Alle, A.; Pinto, J. M. Mixed-Integer Programming Models for the Scheduling and Operational Optimization of Multiproduct Continuous Plants. Ind. Eng. Chem. Res. 2002, 41, 2689. (22) Brooke, A.; Kendrick, D.; Meeraus, A. GAMSsA Users’ Guide; The Scientific Press: Redwood City, CA, 1998. (23) CPLEX 7.0 User’s Manual; ILOG Co: Gentilly, France, 2000. (24) Lawler, E. L. The Traveling Salesman Problem: a Guided Tour of Combinatorial Optimization; Wiley-Interscience: New York, 1985.

Received for review January 12, 2004 Revised manuscript received June 18, 2004 Accepted June 30, 2004 IE049956Z