Ind. Eng. Chem. Res. 2004, 43, 1485-1498
1485
PROCESS DESIGN AND CONTROL Global Optimization for the Cyclic Scheduling and Operation of Multistage Continuous Plants A. Alle† and J. M. Pinto*,†,‡ Department of Chemical Engineering, University of Sao Paulo, Avenida Prof. Luciano Gualberto, t.3, 380, Sa˜ o Paulo SP 05508-900, Brazil, and Department of Chemical and Biological Sciences and Engineering, Polytechnic University, Six Metrotech Center, Brooklyn, New York 11021
This work addresses the global optimization of the simultaneous problem of the cyclic scheduling and operation of multistage continuous plants. In this problem, production rates and yields are additional optimization variables for plant scheduling. The representation proposed for this problem is a mixed-integer nonlinear programming (MINLP) model that has a nonconvex feasible region and a nonconvex objective function. To address nonconvexity, a spatial branch-and-bound global optimization algorithm is developed to solve the model. An illustrative example shows that the global approach is effectively able to yield a more profitable solution than a local optimization algorithm. Moreover, it is shown that modifications in the steps of the global optimization algorithm, such as preferential branching at a variable, can significantly improve its performance. Results also show that local optimization can provide very good estimates for the global solution when processing conditions have narrow variability ranges and plants operate at nearly full capacity. 1. Introduction As previously reported,1 the cyclic scheduling of multistage continuous plants involves tradeoffs between the cycle time, inventory, and changeover costs. More complex tradeoffs arise when process models of the stages are introduced into the scheduling problem.2 In this case, processing rates and yields are additional variables that must be optimized to set the optimal operating point for the plant. The TSPFLOP model,2 which incorporates variable processing rates and yields, describes the simultaneous problem of cyclic scheduling and operational optimization of mulstistage continuous plants. TSPFLOP is a mixed-integer nonlinear programming (MINLP) model with a nonconvex objective function and a nonconvex feasible region. These features prevent solvers based on local optimization techniques from obtaining the global optimum. In fact, a locally optimal schedule can differ to a great extent from a globally optimal one because the two solutions might lie in completely different regions of the solution space. As a consequence, solutions from TSPFLOP can be subject to significant improvement. To avoid suboptimal schedules, global optimization methods are required to solve the TSPFLOP model. The goals of this work are (1) to introduce a generic formulation for the simultaneous problem of the cyclic scheduling and operation of multistage continuous plants with finite intermediate storage and (2) to develop a global optimization algorithm based on the general formulation for spatial branch-and-bound methods3 to globally solve this problem. * To whom correspondence should be addressed. E-mail:
[email protected]. Tel.: 1-718-260-3569. Fax: 1-718-260-3125. † University of Sao Paulo. ‡ Polytechnic University.
The remainder of this paper is structured in the following manner. Section 2 extends the TSPFLOP model, which only accounted for variations in the firststage processing rates and yields, to a generalized formulation in which the rates and yields of every stage are variables for optimizing the operation of multiproduct continuous plants. Section 3 presents a global optimization algorithm, based on the general formulation of spatial branch-and-bound methods,3 to be used for solving TSPFLOP. In section 4, an example with a three-product, two-stage plant illustrates the difference between local and global optimal solutions for the TSPFLOP model. This section also shows how parameters of the global optimization algorithm can be tuned to improve its performance. Finally, section 5 presents the main conclusions of the paper. 2. Extended Model for Simultaneous Scheduling and Operational Optimization The TSPFLOP model was previously proposed2 for the case in which only one stage of the plant has variable yields. However, the model can be extended to the most general case in which all stages present variable yields as functions of the processing rates. To extend the operational optimization to all stages, the production rates are assumed to vary within a range up γlo im e γim e γim
∀ i, m
(1)
Moreover, the mass balance coefficients of product i at every stage m, Rim, are assumed to be exponential functions of the processing rate2
Rim ) exp(γim/bim)
10.1021/ie034118a CCC: $27.50 © 2004 American Chemical Society Published on Web 02/14/2004
∀ i, m
(2)
1486 Ind. Eng. Chem. Res., Vol. 43, No. 6, 2004
Because the relations between the processing times and product amounts were previously considered in an implicit manner by Alle and Pinto2, some constraints of the original TSPFLOP model must be rewritten, such as timing constraints, and others must be included, such as the stagewise mass balance equations. The constraints of the generic TSPFLOP model are as follows:
Transition constraints Zij ) 1 ∑ j*i
∀i
(3)
Zij ) 1 ∑ j*j
∀j
(4)
Cost transition constraints CT )
∑i ∑j ZijCtrij
(5)
Mass balance constraints ∀ i, m
Wim ) γimTPim
(6)
∀ i, m ) 1, ..., M - 1
Wim ) Rim+1Wim+1 ∀i
Fi ) Ri1Wi1
(7) (8)
Timing constraints TS11 )
∑i Zi1τi11 up
-(1 - Zij)Tc
(9)
TSim e TSim+1
∀ i, j * i, j ) 2, ..., NP, ∀ m (10) ∀ i, m ) 1, ..., M - 1
(11)
TSim + TPim e TSim+1 + TPim+1 ∀ i, m ) 1, ..., M - 1 (12) Tc g
∑i TPim + ∑i ∑j τijmZij
∀m
(13)
Demand satisfaction constraints WiM g diTc
∀i
(14)
Intermediate inventory constraints ∀ i, m ) 1, ..., M - 1 (15)
Imaxim ) Wim - Invim
0 e Invim - δim(TSim + TPim - TSim + 1) e ∀ i, m ) 1, ..., M - 1 (16)
(1 - Xim)Wup im
∀ i, m ) 1, ..., M - 1
0 e Invim e Wup imXim Imaxim e Imaxup im
∀ i, m ) 1, ..., M - 1
(17) (18)
Operating constraints ∀ i, m
Rim ) exp(γim/bim) up γlo im e γim e γim
δim e γim
∀ i, m
∀ i, m
δim eRim+1γim+1
(19) (20) (21)
∀ i, m ) 1, ..., M - 1
type
(22)
The objective function of TSPFLOP is given by the gross profit, which is defined as the difference between the revenues of the final products and the overall
terms
fractional bilinear
WiM/Tc, CT/Tc, Fi/Tc, Imaxim/Tc γimTPim, Rim+1Wim+1, δim(TSim + TPim - TSim+1), Rimγim, Ri1Wi1 TPiMWiM/Tc γimRimWim/Tc Rim ) exp(γim/bim)
fractional bilinear fractional trilinear exponential
production costs (raw material, operating, intermediate inventory, changeover, and final inventory costs). Intermediate inventory costs are calculated from the maximum level reached by every intermediate product during the plant cycle, whereas final inventory costs are calculated from the average inventory levels of the final products. To account for variable processing rates at all stages, operating costs for every product i at every stage m, OPim, must be included in the objective function. As previously stated,2 OPim is proportional to the production rate, γim, and to the total amount processed of product i at stage m, RimWim
OPim ) COimγimRimWim
∀ i, m
(23)
where COim is an operating cost coefficient. Thus, the objective function of the extended TSPFLOP model is as follows
max Profit )
eTSjm - (TSim + TPim + Zijτijm) e
(1 - Zij)Tcup
Table 1. Nonlinear Terms Present in the TSPFLOP Model
1
∑i (PiWiM - CfiFi -
[
Tc
M-1
COimγimRimWim - ∑ CinvimImaxim) - CT] ∑ m m 1
(
∑CinvfiWiM 1 2 i
)
TPiM Tc
(24)
The extended TSPFLOP model can be stated as the maximization of the objective function in eq 24 under constraints 1-22. Note that this is a mixed-integer model, because of the binary variables Zij and Xim, as well as a nonlinear model, because of the nonlinear objective function and several nonlinear constraints. Table 1 shows all nonlinear terms included in the TSPFLOP model. The objective function and the feasible region in TSPFLOP are not only nonlinear but also nonconvex, even in the relaxed form in which all binary variables can assume continuous values for the [0, 1] interval. Consequently, global optimality is not guaranteed when TSPFLOP is solved by methods based on the KarushKuhn-Tucker conditions, such as outer approximation13 or generalized Benders decomposition.14 Global optimization methods are required to solve TSPFLOP to proven optimality. 3. Global Optimization Algorithm The global optimization algorithm to be used in the solution of TSPFLOP is based on the general formulation of spatial branch-and-bound methods,3 which was further extended to bilinear and fractional programs4 as well as to mixed-integer nonlinear programs.5,6 The main operations of the algorithm are (1) convex relaxation, (2) the generation of upper and lower bounds for the objective function, and (3) spatial branching.7 The key idea of the spatial branch-and-bound method is the construction of a convex approximation of the
Ind. Eng. Chem. Res., Vol. 43, No. 6, 2004 1487
model. This operation is known as convex relaxation and is obtained from the upper and lower bounds of every variable. Therefore, each variable must be defined within an interval [xlo, xup], even when these bounds are only estimates. Thus, the space defined by the constraints is contained in a hyperrectangle whose edges are the domain intervals of the optimization variables. There exists a straightforward method for obtaining the convex relaxation7 of any nonconvex MINLP model that contains well-behaved functions (bilinear, fractional, exponential, and differentiable univariate functions). Initially, the original problem is reformulated so that it contains only linear constraints and special definitions of the nonconvex terms. In this way, any mixed-integer problem (where the vector x represents continuous variables and y represents binary variables)
max f(x,y) x,y subject to
g(x,y) ) 0 h(x,y) e 0 xl e x e xu y ∈ [0, 1]
(25)
can be transformed into the following standard form
max wobj w subject to
A‚w ) b wlo e w e wup y ∈ [0, 1] wk ≡ wiwj wk ≡
wi wj
wk ≡ wni wk ≡ fn(wi)
∀ (i, j, k) ∈ Tbt
(26)
∀ (i, j, k) ∈ Tlft ∀ (i, k) ∈ Tet ∀ (i, k) ∈ Tuft
where the vector w comprises the original continuous variables x and binary variables y, as well as slack and auxiliary variables that have been introduced to convert functions f and g into the linear equality system and nonconvex terms definitions in problem 26. A and b are the matrix and coefficient vector that define the linear set of constraints in problem 26, respectively. The index obj denotes the position of a single variable that corresponds to the objective function variable within the vector w. This reformulation procedure generates a model that is identical to the original one, i.e., problem 26 is equivalent to problem 25. Therefore, any feasible point from one model can be transformed into a feasible point of the other. The main difference between problems 25 and 26 is that all nonlinear terms in problem 25 are replaced by auxiliary variables that are described by the sets Tbt, Tlft, Tet, and Tuft, which correspond to the bilinear, fractional, exponential, and univariate function
terms, respectively. It should be noted that every element of these sets is either a triplet or a pair of indices that relates the corresponding elements of the vector w in the appropriate manner. For instance
(12, 15, 21) ∈ Tlft S w21 )
w12 w15
The definitions of these new variables generate nonconvex constraints that cannot be present at the convex relaxation. However, they can be replaced by a set of inequalities that are under- and overestimators and that already exist for the bilinear, fractional, and exponential terms (see section 3.2.3). The under- and overestimators make the feasible region of the relaxed model an overestimate of the feasible region of the original model. In the case of a maximization problem such as the TSPFLOP model, the objective function of the relaxed model is an overestimate of the nonconvex objective function of the original model. Next, the spatial branchand-bound method is described for the case of global maximization. After the convex relaxation is obtained, the spatial branch-and-bound algorithm solves the relaxed model to obtain an upper bound (UB) for the global optimum. A lower bound (LB) can be obtained from any feasible point of the original problem. If the difference between the bounds is smaller than a given tolerance , then the search stops. Otherwise, the feasible region is branched through the convenient partition of the domain interval of variables that are present in nonconvex terms. In each region created during the branching, the convex relaxation is again generated and solved, as well the local optimization of the original problem. If an improved LB is found, this bound is updated. If the solution of the convex relaxation is greater than the current LB, then the region is branched, and the algorithm repeats the previous steps of finding the solution of the convex relaxation and a feasible solution of the nonconvex problem in the created regions. Otherwise, the region is discarded. The UB value decreases because of the improved approximation as the search regions become smaller. This iterative procedure goes on until the difference between the UB and LB satisfies the tolerance . In fact, the algorithm works similarly to the branchand-bound (B&B) method for mixed-integer programming (MIP) problems. The main difference is in the branching procedure, which consists of the partitioning of a domain of continuous variables and not of the assignment of a discrete value to a variable, as is done in the solution of mixed-integer programs. The bounding procedures are identical, as both algorithms generate upper and lower bounds at every node, with the upper bounds being obtained from a relaxation of the problem, which in the case of MIPs consists of allowing discrete variables to assume continuous values and of solving the problem as a continuous optimization problem. The following section outlines the spatial branch-andbound algorithm developed for the global optimization of the TSPFLOP model. 3.1. Algorithm Outline. For a maximization problem given by general formulation 25, the outline of the algorithm is as follows:7 Step 1. Search Initialization. Set f lo r -∞ and the list L of regions to the single domain R ≡ [wlo wup].
1488 Ind. Eng. Chem. Res., Vol. 43, No. 6, 2004
Figure 1. Scheme for the global optimization algorithm.
Step 2. Initial Bound Tightening. Tighten the bounds in the region R ≡ [wlo wup]. If any lower bound exceeds the respective upper bound, go to step 9. Step 3. Subregion Selection. If L ) Ø, go to step 9; otherwise, choose a subregion R from the list L such that
R ) arg max fR,up R∈L Step 4. Bound Tightening in R. Tighten the bounds for subregion R.
If the lower bounds for one or more variables exceed the corresponding upper bounds, go to step 8. Step 5. Determination of Objective Function Upper Bound in R. Generate the convex relaxation of problem 25 in R. Solve this convex relaxation to determine the objective function upper bound, f R,up. If the convex relaxation is infeasible or if f R,up - e f lo, go to step 8. Step 6. Determination of Objective Function Lower Bound in R. Solve (locally) problem 25 in R to generate a lower bound, f R,lo.
Ind. Eng. Chem. Res., Vol. 43, No. 6, 2004 1489
If a local maximum cannot be determined or if f R,lo < f lo, go to step 7. Otherwise, set f lo r f R,lo and delete from L all subregions S such that f S,up - e f lo. If f R,up - f R,lo e , go to step 8. Step 7. Partition of Region R. Apply a branching rule to subregion R to choose a branch variable and its corresponding value on which to branch. Create two new subregions, R lo and Rup, by partitioning R at the branch point. Add these to the list L. Step 8. Deletion of Region R. Delete subregion R from the list of regions L. Go to step 3. Step 9. Search Termination. If f lo ) -∞, the problem is globally infeasible. Otherwise, the global optimum objective function value is f lo. The algorithm keeps a list L of subregions (hyperrectangles) that must be examined for the global optimum. It also keeps a record of the best lower bound, f lo, for the global optimum of the objective function. The main iteration loop (steps 3-8) processes a single region R from the list L and calculates the corresponding lower and upper bounds of the objective function (steps 5 and 6, respectively). From these values, the algorithm can either discard (step 8) or partition the region into two subregions (step 7), which are then placed for further examination. Figure 1 shows the general form of the spatial branchand-bound method implemented to solve the TSPFLOP model. The following section contains a more detailed explanation on the operations performed in the main steps of the algorithm. 3.2. Details of the Algorithm Implementation. 3.2.1. Identification of a Good Initial Solution. A first modification introduced in the spatial branch-andbound algorithm developed here is the solution of the original problem by the solver DICOPT,8 to generate a local solution that serves as a good first estimate of the global optimum, instead of initializing Profitlo ) -∞. This estimate of Profitlo is used in the second step (section 3.2.2) within the contraction procedure,9 which allows for a reduction of the variable intervals. 3.2.2. Bounds Tightening (Steps 2 And 4). The spatial branch-and-bound algorithm has the property of convergence to the global optimum in a finite number of iterations given a tolerance even without the bounds-tightening steps. However, the latter are applied to accelerate convergence. There is a tradeoff between the time spent in these steps and the time saved in the remaining steps. Basically, there are two different approaches of such tightening:7 one based on optimization and the other on feasibility. 3.2.2.1. Optimization-Based Bounds Tightening. The maximum and minimum values that a variable wi can assume can be calculated by solving two convex optimization problems
This technique can be applied to each variable wi separately, updating the convex relaxation constraints whenever possible. A single pass of the n variables requires the solution of 2n optimization problems. As the maximum tightening would be achieved only after several passes and significant computational effort, it is recommended3 that this type of tightening be performed only in step 2. The TSPFLOP model has binary variables. Previous experience has shown that it is not beneficial to fix these variables in step 2 because of the time spent in solving the optimization problems and the little progress made to eliminate some of the 0-1 values for the next steps of the algorithm. Therefore, only continuous variables are subject to bound tightening in step 2 for the TSPFLOP model. The constraint Profit g Profitlo is included in the convex relaxation constraint set. Therefore, the boundtightening procedure is similar to the branch-andcontract algorithm.9 3.2.2.2. Feasibility-Based Bounds Tightening. A computationally cheaper, although less effective, approach is used in step 4 to tighten the bounds of the variables in the nonconvex terms of the TSPFLOP model. The nonconvex terms in TSPFLOP are transformed such that they belong to only one of the two sets: bilinear terms and fractional terms. The exponential term can be directly included in the convex relaxation as is shown further on. For the triplet (i, j, k) of the Tbt set of bilinear terms, the following expressions to tighten bounds of wi are suggested10
[ [
wup i
) min
subject to convex relaxation constraints wlo ew ewup and wlo i ) minwi w
subject to convex relaxation constraints wlo ew ewup (27)
wup i ,
)] )]
lo up up wlo k wk w k w k max lo, up, lo , up wj wj wj wj
(28)
(29)
Expressions for wj are similar. For fractional terms (i, j, k) ∈ Tlft, the following expressions are suggested1
wlo i ) lo lo lo up up lo up up max[wlo i , min(wk wj , wk wj , wk wj , wk wj )] (30)
wup i ) lo lo lo up up lo up up min[wup i , max(wk wj , wk wj , wk wj , wk wj )] (31)
wup i ) maxwi w
( (
lo up up wlo k wk w k w k , , , up wlo wlo wup j wj j j
lo wlo i ) max wi , min
[ [
( (
lo wlo j ) max wj , min
wup j
) min
wup j ,
)] )]
lo up wup wlo i wi w i i , , , lo up lo up wk wk wk wk
lo up wup wlo i wi w i i max lo, up, lo , up wk wk wk wk
(32)
(33)
In eqs 28, 29, 32, and 33, any term in which the denominator is zero must be omitted from the list of
1490 Ind. Eng. Chem. Res., Vol. 43, No. 6, 2004 Table 2. Estimators of McCormick11 underestimators bilinear terms (i, j, k) ∈ Tbt wk ≡ w i wj linear fractional terms (i, j, k) ∈ Tlft wk ≡ wi/wj
up lo up wk e wlo i wj + w i wj - w i w j lo up lo wk e wup w + w w w w j i j i i j
lo lo lo wi g wlo k wj + wkwj - wk wj up up up wi g wup w + w w w w j k k j k j
up lo up wi e wlo k wj + wkwj - wk wj lo up lo wi e wup w + w w w w j k k j k j
wk g wk g
min(‚) and max(‚) operators. The algorithm developed in this work makes use of constraints 28-33 for the bound tightening performed within step 4. 3.2.3. Convex Relaxation: Upper Bound for the Objective Function (Step 5). Step 5 of the algorithm determines an upper bound for the objective function of the exact problem in formulation 25 for the subregion R through the solution of its convex relaxation. In this work, two types of convex relaxation are compared. The first, linear relaxation (LR), contains only linear constraints, whereas the second, nonlinear relaxation (NLR), also contains nonlinear albeit convex constraints. The estimators from McCormick11 shown in Table 2 are used in the LR. To build the LR, the TSPFLOP model must be initially reformulated so as to contain only linear constraints and nonconvex term definitions. This can be done in the following manner
wim1 ) Wim
∀ i, m
(34)
wim2 ) Rim
∀ i, m
(35)
wim3 ) γim
∀ i, m
(36)
wim4 ) TPim wim5 ) δim
∀ i, m
(37)
∀ i, m ) 1, ..., M - 1
(38)
wim6 ) TSim + TPim - TSim+1 ∀ i, m ) 1, ..., M - 1 (39) wim7 ) δim(TSim + TPim - TSim+1) ) wim6wim5 ∀ i, m ) 1, ..., M - 1 (40) wim8 ) Rimγim ) wim2wim3 ∀ i, m ) 2, ..., M (41) wim9 ) RimWim ) wim2wim1 wim10 ) γimRimWim ) wim3wim9 wi11 ) WiMTPiM ) wiM1wiM4 w12 )
∀ i, m ∀ i, m ∀i
(42) (43) (44)
∑i (PiwiM1 - ∑j ZijCtrij - CfiFi M-1
COimwim10 - ∑ CinvimImaxim - 0.5Cinvfiwi11) ∑ m m w13 ) Tc w14 )
w12 w13
(45) (46) (47)
From the mass balance constraint in eq 6, it follows that
Wim ) γimTPim w wim1 ) wim3wim4
overestimators
lo lo lo wlo i wj + w i wj - w i wj up up up wup w + w w w w j i j i i j
∀ i, m (48)
Constraints 34-39 and 45-46 only substitute variables of the original model by new variables w. Most of
Table 3. Bilinear and Fractional Terms in the Reformulated TSPFLOP Model equation
classification
triplet
40 41 42 43 44 47 48
bilinear bilinear bilinear bilinear bilinear linear fractional bilinear
(im6, im5, im7), ∀ i, m ) 1, ..., M - 1 (im2, im3, im8), ∀ i, m ) 2, ..., M (im2, im1, im9), ∀ i, m (im3, im9, im10), ∀ i, m (iM1, iM4, iM11), ∀ i (12, 13, 14) (im3, im4, im1), ∀ i, m
these w variables have three indices because they are defined in the domain of product i vs stage m, for instance, Wim, which is defined in eq 34. Thus, the triplets that compose the sets Tbt and Tlft are given by expressions of the type (imi′, imj′, imk′). In the case of eq 44, the new variable does not need to be defined in the domain of the stages (index m) because it is formed from variables only of the last stage (m ) M). Table 3 lists the bilinear and linear fractional terms in the reformulated TSPFLOP model. The sets Tbt and Tlft for the TSPFLOP model are as follows
Tbt ) {{(im6, im5, im7)|i ) 1, .... NP; m ) 1, ..., M - 1} ∪ {(im2, im3, im8)|i ) 1, ..., NP; m ) 2, ..., M} ∪ {(im2, im1, im9), (im3, im9, im10), (im3, im4, im1)|i ) 1, ..., NP; m ) 1, ..., M}} ∪ {(iM1, iM4, i11)|i ) 1, ..., NP} (49) Tlft ) {(12,13,14)}
(50)
The linear relaxation (LR) can be easily obtained from the reformulated model by the straightforward application of the estimators of McCormick11 listed in Table 2 to the sets Tbt and Tlft defined in eqs 49 and 50. The tightest estimators for a given function are called convex envelopes. McCormick’s estimators11 are convex envelopes for bilinear functions but not for linear fractional functions.4 For the latter, the convex envelopes are given by the overestimators of Quesada and Grossmann,4 listed in Table 4. The nonlinear convex relaxation (NLR) uses the estimators of Quesada and Grossmann4 for the fractional terms and McCormick’s estimators11 for the bilinear and fractional terms to achieve a tighter relaxation of the original problem. To fulfill the convexity conditions given in Table 4, an assumption is made that the variable w12 is always positive. This assumption is reasonable, as the plant should be profitable at least in the vicinity of the operating point. In this case, the underestimators are valid, but the overestimators are not. Thus, only the underestimators of Quesada and Grossmann4 can be used in the NLR. Nonlinear eq 2 defines a nonconvex surface in the continuous space. Therefore, it is a nonconvex constraint that must be relaxed. The relaxation of eq 2 can be done
Ind. Eng. Chem. Res., Vol. 43, No. 6, 2004 1491
in two different ways. In the case of LR, the equation is replaced by a set of underestimators that are tangent lo up up , wim3 to the surface in the points wim3 , and 0.5(wim3 + lo wim3) and by a single overestimator that is a secant to lo up and wim3 ) the surface in the points (wim3
( ) ( )
lo wim3 1 lo lo exp (wim3 - wim3 ) + wim2 bim bim
(51)
up wim3 1 up up exp (wim3 - wim3 ) + wim2 wim2 g bim bim
(52)
wim2 g
wim2 g
1 exp bim
wim2 e
[
up 0.5(wim3
(
+ bim
]
lo wim3 )
[
up lo wim3 - wim3
)
]
lo lo (wim3 - wim3 ) + wim2 (54)
However, in the case of the nonlinear relaxation (NLR), eq 2 can be transformed into the following inequality
∀ i, m
(55)
Reformulated TSPFLOP
∑i CinvfiwiM1
subject to
Zij ) 1 ∑ j*i
∀i
Zij ) 1 ∑ i*j
∀j
mass balance constraints wim1 ) wim+1,9 Fi ) wi,1,9
∑i wim4 + ∑i ∑j τijmZij
∀i
∀ i, m ) 1, ..., M - 1
∀ i, m ) 1, ..., M - 1
∀m
demand constraints ∀i
wiM1 g diw13
maximum intermediate inventory constraints ∀ i, m ) 1, ..., M - 1
Imaxim ) wim1 - Invim
up 0 e Invim - wim7 e (1 - Xim)wim1 ∀ i, m ) 1, ..., M - 1
∀ i, m ) 1, ..., M - 1
up 0 e Invim e wim1 Xim
Imaxim e Imaxup im
∀ i, m ) 1, ..., M - 1
(56)
∀ i, m
wim2 g exp(wim3/bim) wim5 e wim3
(55)
∀ i, m ∀ i, m ) 1, ..., M - 1
wim5 e wim+1,8 w12 definition w12 )
∑i (PiwiM1 - ∑j ZijCtrij - CfiFi M-1
COimwim10 - ∑ CinvimImaxim - 0.5Cinvfiwi,11) ∑ m m (45) nonconvex term definitions wim7 ) wim6wim5
∀ i, m ) 1, ..., M - 1 (40)
wim8 ) wim2wim3
∀ i, m ) 2, ..., M
(41)
wim9 ) wim2wim1
∀ i, m
(42)
wim10 ) wim3wim9 wi11 ) wiM1wiM4 w14 )
transition constraints
∀ i, m ) 1, ..., M - 1
TSim eTSim+1
operating constraints
This substitution does not affect the global optimum because inequality 55 must be active in the global optimum given that wim2 is inversely proportional to the processing yield. Therefore, wim2 is maximized for a given processing rate to obtain the minimum raw material consumption. Constraint 55 defines a convex region as the epigraph of a function exp(x), which is convex.12 Therefore, it need not be substituted by either a linear constraint or a nonlinear definition of a new variable. In summary, the LR and NLR differ only in two points: (1) the NLR contains the Quesada and Grossmann4 estimators for fractional terms and (2) the NLR makes use of a nonlinear convex constraint that substitutes for eq 2. The reformulated model containing only convex constraints and nonconvex term definitions is given as follows
max Profit ) w14 - 0.5
∑i Zi1τi11
-(1 - Zij)Tcup e TSjm - wim6 e (1 - Zij)Tcup ∀ i, j * i, j ) 2, ..., NP; ∀ m
w13 g
up [wim3 - 0.5(wim3 +
wim2 g exp(wim3/bim)
TS11 )
TSim + wim4 e TSim+1 + wim +1,4
up lo 0.5(wim3 + wim3 ) lo )] + exp (53) wim3 bim
up lo wim2 - wim2
timing constraints
∀ i, m
(43)
∀i
(44)
w12 w13
wim1 ) wim3wim4
(47) ∀ i, m
(48)
upper and lower bounds wlo e w e wup
(57) ∀ i, m
0 e Invim e Invup im 0 e Imaxim e Imaxup im up Flo i e Fi e Fi
∀i
∀ i, m
(58) (58a) (58b)
1492 Ind. Eng. Chem. Res., Vol. 43, No. 6, 2004 Table 4. Nonlinear Estimators of Quesada and Grossmann4 for Linear Fractional Terms wk ≡ wi/wj, (i, j, k) ∈ Tlft estimators
convexity conditions
underestimators lo up lo wk g wup i /wj + wi/wj - wi /wj up lo up wk g wlo i /wj + wi/wj - wi /wj
overestimators up up up wk e wup i /wj + wi/wj - wi /wj lo lo lo wk e wlo i /wj + wi/wj - wi /wj
lo wup i e 0, wj > 0 or g 0, wup wup i j < 0 lo lo wi e 0, wj > 0 or up wlo i e 0, wj < 0
bilinear terms ∀ (i, j, k) ∈ Tbt
(59)
linear fractional terms AElft ijk
)
|wR k
-
() wR i
wR j
|
R R n AEet ik ) |wk - (wi ) |
∀ (i, k) ∈ Tet
(61)
∀ (i, k) ∈ Tuft
(62)
univariate function terms lo wup i g 0, wj > 0 or wup e 0, wlo i j < 0 lo wlo g 0, w i j > 0 or up wlo i e 0, wj < 0
Note that the mass balance in eq 6 in the original model is equivalent to the nonconvex term definition in eq 48 in the reformulated model. The bounds defined in constraints 57 and 58 are obtained from the plant capacity limits (inventory tanks and production rates) and from process parameters (transition costs and times). It is assumed that the storage capacity for the final products implicitly provides the upper bound for the cycle time, even though the model does not explicitly consider final inventory. 3.2.4. Lower Limit for the Objective Function (Step 6). The goal of step 6 is to provide a lower bound of the original problem in the current subregion. Every feasible point of the original problem provides this bound. Two alternatives exist. (1) The first alternative is local optimization of the MINLP problem. The reformulated model is solved by an MINLP method, as the outer approximation8,13 or the generalized Benders decomposition.14 In the case of the outer approximation method, every binary variable must take part only in linear constraints, as is the case for the TSPFLOP model. (2) The second alternative is local optimization with fixed binary variables. The goal of this alternative is to reduce the computational effort with respect to this step. A simpler solution procedure rather than the local optimum is sought as the lower bound for the objective function. In this approach, the most direct method is to fix the values of the binary variables obtained by the solution of the convex relaxation in step 5. In the spatial branch-and-bound algorithm implemented in this work, the first alternative was chosen because the MINLP solver DICOPT8 has previously been shown to be very robust in solving TSPFLOP locally.2 3.2.5. Branching (Step 7). Step 7 of the algorithm chooses a variable from the nonconvex sets at which to branch. This choice results in the partitioning of the subregion R that has just been investigated. The objective of this choice is to minimize the feasibility gap between the convex relaxation and the exact problem.7 To do that, the approximation error, AE, for every nonlinear term of problem 26 must be sought. The AE is defined as follows
R R R AEbt ijk ) |wk - wi wj |
exponentiation terms
∀ (i, j, k) ∈ Tlft
(60)
R R AEuft ik ) |wk - fn(wi )|
In the case of the NLR, only eqs 59 and 60 are used. In the case of the LR, eq 62 is also used for the exponential function. There is no need to branch at variables that take part only in linear terms because doing so cannot achieve any improvement in the quality of the convex relaxation.15 Smith and Pantelides7 suggested that branching be performed at variables of the nonconvex term with the largest AE
AE* ) max[AEbt ijk
∀ (i, j, k) ∈ Tbt, AElft ijk
∀ (i, j, k) ∈ Tlft, AEet ik
∀ (i, k) ∈ Tet, AEuft ik ∀ (i, k) ∈ Tuft] (63)
This work adopted this rule with a small modification: it is possible to branch preferentially at the cycle time Tc (or w13). In fact, Tc participates in a single nonconvex term, namely, the w14 definition, which is the only fractional term in TSPFLOP. Nevertheless, Tc directly bounds most variables in the model, e.g., the processing times and the amounts to be produced. Therefore, the impact of branching at Tc cannot be measured as a function of the decrease in the error of w14 because, as a result of the loop in step 4, the partition of the Tc interval influences the contraction of the bounds for almost all remaining variables. Hence, the algorithm sets priority on branching in Tc as long as the relative difference between the relaxation solution in the current node, Profitc,up, and that in its parent node, Profitp,up, is larger than the predetermined parameter con
|Profitc,up - Profitp,up| g con Profitc,up
(64)
In the case of nonconvex terms defined by two variables, the following expression is used to determine the branching variable
[| (
wp ) arg min 0.5 i,j
R,lo wR i - wi
wR,up - wR,lo i i
|
0.5 -
(
)|
,
R,lo wR j - wj
wR,up - wR,lo j j
)|]
(65)
i.e., the variable whose value wR p in the solution of step 5 (convex relaxation) is closer to the middle point of its R,up ]. In this case, the branchdomain interval [wR,lo p , wp R ing point is the wp value in the convex relaxation solution. However, the wR p value might not represent a sigR,up nificant reduction of the interval [wR,lo ], as it can p , wp be very close to one of its bounds. Thus, Zamora and Grossmann9 suggested a small modification of this criterion: in cases when the partition in wR p does not R,up reduce the domain [wR,lo , w ] to a certain fraction, p p dim, then the branching is done at the middle point of the domain of the wp variable. This criterion is adopted
Ind. Eng. Chem. Res., Vol. 43, No. 6, 2004 1493 Table 5. Data for the Illustrative Example product
Pi ($/ton)
di (ton/h)
Cfi ($/ton)
COi1 ($ h/ton2)
COi2 ($ h/ton2)
Cinvi1 ($/ton)
Cinvfi [$/(ton h)]
Imaxi1 (ton)
A B C
134 269 210
0.0160 0.0158 0.0145
23.3 25.5 16.1
3.66 3.96 4.59
6.00 8.22 4.71
6.84 3.96 4.50
0.089 0.076 0.038
4.814 3.604 7.812
stage 1
stage 2
product
bi1
lo γi1 (ton/h)
up γi1 (ton/h)
bi2
lo γi2 (ton/h)
up γi2 (ton/h)
A B C
5.10 8.13 3.39
0.028 0.021 0.045
0.560 0.420 0.907
6.00 9.00 6.48
0.025 0.019 0.038
0.491 0.383 0.767
τijm (h) stage 1
stage 2
Ctrij ($)
product
A
B
C
A
B
C
A
B
C
A B C Tclo ) 0 Tcup ) 2000 h
0 5.0 4.6
6.8 0 6.2
3.7 3 0
0 2.6 5.0
5.6 0 6.9
1.7 2.0 0
0 3326 3069
1943 0 1245
1694 1137 0
in the spatial branch-and-bound algorithm developed here. Therefore
If 1 - dim e
(
R,lo wR p - wp
wR,up - wR,lo p p
)
e dim, then branch at wR p
otherwise, branch at
wR,up + wR,lo p p 2
Several other branching strategies are available, as shown in Adjiman et al.16 4. Results The spatial branch-and-bound method was implemented in GAMS.17 The DICOPT solver8 was used to solve the original MINLP model and the NLR model. The CONOPT218 and XPRESS 12.519 solvers were used to solve the NLP subproblem and MILP master problem, respectively. The XPRESS solver was also used to solve the LR problem, which is an MILP (mixed-integer linear programming) problem. All problems were solved on a Pentium III platform (500 MHz, 512 MB of RAM) under Windows NT 4.0. Bounds for processing rates were not tightened within the optimization method in step 2 because they can be obtained in a straightforward way through the direct substitution of the bounds for Wim and TPim into the mass balance equations. Bounds for the mass balance coefficients Rim are also easily obtained from the application of the bounds on γim to eq 2. To save time, only Tc, Wim, TPim, w12, and the auxiliary variables w14 and wi11 are chosen to have their bounds tightened during step 2. The variable wi11 is chosen because of its potential impact on the objective function given that it is the product between two variables of large magnitude. This section is divided as follows: Initially, an illustrative example is presented to show how different the global optimal solution can be from a locally optimal one. Then, a sensitivity analysis is done to show how the different types of branching strategies influence the performance of the algorithm. The third subsection shows how the algorithm solves problems of various levels of difficulty and how the DICOPT solver can
provide the global optimum when problem parameters reflect real-world conditions of operation of continuous plants. 4.1. Illustrative Example. The goal of this example is to show that the difference between a locally optimal schedule and the globally optimal schedule can be significant. Table 5 contains data for a two-stage, threeproduct continuous plant. Figure 2 shows the difference between the locally and globally optimal schedules of this plant. Note that the inventory level profiles in this figure correspond to the intermediate products between stages 1 and 2. Note in Figure 2 that the global profitability is approximately 30% greater than the local value. Note also the significant difference between the processing times in the two schedules, as well as the different inventory profiles. Table 6 reports the production rates obtained in the local and global optima. In Table 6, note that the production rates of the local optimum are at their upper bounds, except for the production rate of C in the first stage. In the globally optimum case, however, the production rates are at their upper-bound values to obtain the maximum productivity, whereas the rates in the first stage are set at values slightly lower than the upper bound to reduce the raw material consumption and the operating and intermediate inventory costs. 4.2. Strategies for the Global Optimization Algorithm. The performance of the global optimization algorithm can be significantly improved through the appropriate tuning of its internal operations. These variants were mentioned along the description of the algorithm (section 3.2). Here, the impact of the proposed variants in the algorithm performance is studied for the solution of the illustrative example. The CPU times and numbers of iterations correspond to solutions obtained for a tolerance of ) 0.10 × Profitlo. A maximum number of 200 iterations is set for attaining the global optimum within the tolerance. The gap between the local optimum and the solution of the convex relaxation in the first node of the branchand-bound tree is given by
Gap(1) )
(1) Profit(1) conv - Profitlocal
Profit(1) conv
(66)
1494 Ind. Eng. Chem. Res., Vol. 43, No. 6, 2004
Figure 3. Effect of repetitions in the loop of step 2. Table 7. Effect of Repetitions in the Loop of Step 2 repetitions
CPU time in step 2 (s)
Gap(1) (%)
overall CPU time (s)
iterations
2 3 4 5
23.1 35.1 48.1 59.9
82.5 77.0 72.3 68.2
409.0 355.6 369.0 372.3
>200 184 180 174
Table 8. Influence of the Minimum Relative Interval Contraction Parameter (Edim)
Figure 2. Local and global optima for the multistage continuous plant of the illustrative example. Table 6. Optimal Production Rates local optimum
global optimum
product
γi1 (ton/h)
γi2 ($/ton)
γi1 (ton/h)
γi2 ($/ton)
A B C
0.560 0.420 0.877
0.491 0.383 0.767
0.547 0.411 0.863
0.491 0.383 0.767
where Profit(1) conv is the objective function value of the convex relaxation solution and Profit(1) local is the value of the local solution of the original problem in the first
dim
overall CPU time (s)
iterations
0.0 0.8 0.9 1.0
465.1 377.4 277.7 284.3
188 196 144 148
node. The greater the value of Gap(1), the greater the initial distance between the original problem and its convex relaxation. Initially, the influence of the tightening of the initial bounds in step 2 is analyzed. The LR is used, with branching priority at Tc with con ) 2% and dim ) 0.70. From Table 7 and Figure 3, it can be seen that there is a tradeoff between the number of repetitions of the loop in step 2 and the overall time for solving the problem. As the number of repetitions in step 2 increases, the Gap(1) value decreases, and the algorithm requires fewer iterations to achieve the global optimum. On the other hand, the time spent in step 2 increases. In this case, three repetitions require the smallest overall solution time. The influence of the dim parameter, which measures the minimum fraction for contracting the branching interval, is shown in Table 8. The LR with preferential branching at Tc with con ) 2% and with three repetitions in step 2 is used. In Table 8, dim ) 0.0 means bipartition of the domain interval, whereas dim ) 1.0 means branching always at the solution of the convex relaxation. Note that there is a tradeoff between dim and the CPU time. The optimum performance occurs at dim ) 0.90, i.e., branching at the convex relaxation value unless this means a contraction of less than 10% of the domain interval. The con value for priority branching at Tc also influences the algorithm performance as Table 9 shows. In this table, con ) 0.0 means branching only in Tc, whereas con ) 1.0 means that there is no preference for branching at Tc. The tests used the LR, dim ) 0.90, and three repetitions in step 2. Table 9 shows that con is a crucial parameter for the solution of the illustrative example. In fact, the values
Ind. Eng. Chem. Res., Vol. 43, No. 6, 2004 1495 Table 9. Influence of the Minimum Relative Upper Bound Improvement for Preferential Branching at Tc (Econ) con
CPU time (s)
iterations
0.0 0.01 0.02 0.03 0.05 1.0
414.3 414.3 277.7 257.8 423.5 550.0
>200 >200 144 138 >200 >200
Table 10. Algorithm Performance Dependence on Tc Bounds Tcup
Gap(1)
CPU time (s)
iterations
500 1000 1500 2000
0.2 27.9 56.7 77.0
40.8 51.7 164.8 277.7
2 10 74 144
in the top and bottom lines of the table did not present convergence within the maximum number of iterations allowed. The best performance was achieved at around 0.03 (or 3%). This means that preferential branching at Tc is a good strategy even at deeper levels of the spatial branch-and-bound tree. The performance of the global optimization algorithm also depends on the specific instance of the problem to be solved: its size (in terms of the numbers of products and stages), the width of the domain intervals for the continuous variables, and the nature of the data (prices, demands, exponential function coefficients) all influence the computational effort required to obtain the global optimum. The model parameters influence the computational time for the convergence of the spatial branch-andbound algorithm. The larger the bi coefficient, the more nonlinear the problem becomes, and as a consequence, the worse the linear relaxation (LR) becomes. Small demands imply greater slacks in the production capacity. This affects the performance of the algorithm because the generation of feasible schedules becomes easier. The larger the domain intervals, the longer the solution times. Table 10 shows the influence of the cycle time bounds on the Gap(1) value and on the solution time of the illustrative example. All data from Table 5 are kept, except for the upper bound on the cycle time. The algorithm is performed with the LR, three repetitions in step 2, dim ) 0.90, and con ) 0.02. Note in Table 10 that Gap(1) linearly increases whereas the CPU time and the number of iterations exponentially increase with Tcup. Figure 4 illustrates these trends. The exponential dependence of the gap on the CPU time indicates the importance of step 2 until Gap(1) assumes a reasonable value. Figure 5 shows the search tree for the illustrative example when Tcup ) 1000. Note that the algorithm builds a binary tree to store the partial solutions. Therefore, the parent node k generates two child nodes, 2k and 2k + 1. Because of its priority branching, Figure 5 shows that only the bounds for Tc were partitioned during the global search. Initially, the Tclo value is tightened and goes from 1 to 171.2 h. Nodes 4, 5, 12, and 13 are pruned because the difference between the convex solution and the lower bound for the problem (120.7 $/h) is smaller than the relative tolerance of 10% (or 12.07 $/h in
Figure 4. Influence of Tcup on the initial convex relaxation and on the CPU time for the solution of the illustrative example.
absolute figures). The global optimum is thus found at node 7, whose value is 120.7 $/h. 4.3. Comparison between the Types of Relaxation for the TSPFLOP Model. In this section, the linear relaxation and nonlinear convex relaxation approaches are compared for a set of problems with different degrees of difficulty. The data for the problems were generated randomly using a technique similar to that of Dobson,20 who studied a reduced version of the problem described in this work containing a single stage that operates at fixed yields and rates for every product. Dobson used three basic parameters to generate problem data: the utilization level u, the product diversity factor r, and the sequence dependence factor d. The utilization level u is the fraction of the total cycle time Tc that is effectively used for the production. The diversity factor r represents the extent of the variation of parameters such as prices, inventory costs, and demands. The sequence dependence diversity factor d describes the extent of the variation in the sequence-dependent transition times and costs. To precisely describe the data generation procedure, let U[a, b] represent a random variable uniformly distributed in the [a, b] interval. Minimum values for the product prices, demands, yield decay coefficients, final and intermediate inventory costs, operating and raw material costs, and transition times and costs are denoted as Pmin, dmin, bmin, Cinvfmin, Cinvmin, COmin, Cfmin, Ctrmin, and τmin, respectively. The parameters that are independent of the sequence are obtained as follows: Pi ) PminU[1, r], di ) dminU[1, r], Cinvfi ) CinvfminU[1, r], Cinvim ) CinvfminU[1, r], Cfi ) CfminU[1, r], COim ) COminU[1, r], and bmin ) bminU[1, r]. A basic transition cost, Ctri ) CtrminU[1, r], and a basic transition time, τi ) τminU[1, r], are chosen for every product i. Using these values, the transition times and costs that are dependent on the sequence are generated from the d parameter as follows: Ctrij ) CtriU[1, d] and τijm ) τiU[1, d], respectively. Here, the values chosen are r ) 3 and d ) 3. To simplify the procedure, it is assumed that only the last stage can be the slowest or bottleneck stage. Therefore, the production rates at stage M are no larger than the rates in the previous stages. The production rates are obtained from the utilization level u and the difference between the productivities of consecutive stages. The minimum utilization level of the plant is defined as the fraction of the cycle time that is spent by the bottleneck stage to satisfy the minimum
1496 Ind. Eng. Chem. Res., Vol. 43, No. 6, 2004
Figure 5. Evolution of the spatial branch-and-bound search for the global optimum of the illustrative example with Tcup ) 1000.
product demand
umin )
∑i
di
γup iM
(67)
The parameter umin is a function of the demands di and the rates γup iM, which are randomly generated. As a consequence, it is not straightforward to estimate the value of umin a priori. To circumvent this difficulty, a small range of u, i.e., [ulo, uup], is used instead of an exact value. The procedure to generate γup iM values for a given number of products, NP, compatible with the minimum degree of utilization, umin, is as follows: 1. Set the utilization range estimate, [ulo, uup]. 2. Obtain u j min from a uniform probabilistic distribution function in this interval. 3. For every i, calculate
γup iM ) dmin
U[1, r] NP u j min
where NP is the number of products processed in the plant. Observe that the numerator has the same uniform distribution function as is used to obtain the demand di. 4. If
ulo e
∑i
di up
e uup
γiM
then stop; otherwise, go to step 2. This procedure has proved very efficient for generating γup iM values for given NP and umin. In most cases, step 2 is done only once.
The utilization range used was umin ) [0.79, 0.80], because the plants are commonly designed to have little idle capacity. up From the γup iM values, the other γim values are generated as follows up up γup im ) fbottleRim+1 γim+1
∀ i, m ) 1, ..., M - 1 (68)
where fbottle is a given factor, denoted the bottleneck factor, that defines the quotient between the production rates and the consumption rates of consecutive stages, given that the last stage is the slowest. This work used fbottle ) 1.05 because flowshop plants usually have stages designed to run at compatible rates. The γlo im values are obtained as follows up γlo im ) frangeγim
∀ i, m
(69)
where frange is the quotient between the minimum and maximum rates of a given product in a given stage. The value frange ) 0.80 was chosen because real-world continuous plants have little flexibility to operate outside their nominal processing capacity. A real-world lube oil plant, on which this work has been based, cannot experiences changes in rate of more than 10% from their nominal values. It is assumed that intermediate storage tanks were designed to keep the subsequent unit running at its minimum consumption rate even when the previous unit is operating at its maximum production rate during the maximum cycle time. Therefore, the maximum storage capacity for product i at the intermediate stage m is given by up up lo lo Imaxup im ) (γim - Rim+1 γim+1)Tc ∀ i, m ) 1, ..., M - 1 (70)
Ind. Eng. Chem. Res., Vol. 43, No. 6, 2004 1497 Table 11. Parameters for Generation of Problem Data parameter
value
parameter
value
Pmin ($/ton) Ctrmin ($) COmin ($ h/ton2) dmin (ton/h) Cinvmin ($/ton)
100 800 3.00 0.010 3.00
Cfmin ($/ton) τmin (h) Cinvfmin [$/(h ton)] bmin (ton/h)
10.00 1.0 0.03 0.5
Table 13. Gaps of the Linear and Nonlinear Convex Relaxations in the First Node Gap(1) (%) problem
LR
NLR
1 2 3 4 5 6 7
43.065 53.752 33.361 54.313 55.942 53.881 88.201
43.064 53.752 33.361 54.312 55.942 53.872 88.201
Table 12. Comparison between the DICOPT Algorithm with the Global Optimization Algorithm Using the Linear (LR) and Nonlinear (NLR) Relaxations DICOPT
global/ LR
CPU profit time profit prob NP M ($/h) (s) ($/h) iter 1 2 3 4 5 6 7 a
3 4 5 6 5 7 6
2 2 2 2 3 2 3
3.58 3.53 6.37 8.47 3.23 8.24 5.64
1.34 1.67 2.50 3.45 4.62 6.68 8.26
3.58 3.53 6.37 8.47 3.23 8.24 5.64
30 90 22 64 70 56 82
CPU time (s) 60.4 157 153 436 576 732 862
global/ NLR profit ($/h) 3.58
8.47 3.23 8.24 5.64
iter
CPU time (s)
30 111.6 NCa a NC 64 556 70 1194 56 1474 82 1820
No convergence within 200 iterations
The minimum costs are chosen such that all components of the objective function have a significant share of the overall cost. The parameters for data generation are listed in Table 11. The value of Tcup was chosen as 4320 h (6 months). Table 12 shows the performances of the local optimization algorithm DICOPT and of the global optimization algorithm using LR and NLR in solving problems of different sizes, with data randomly generated by the previously described procedure. The global optimization algorithm uses ) 0.10 × Profitlo, con ) 0.03, dim ) 0.90, and three repetitions of step 2. Table 12 shows that DICOPT is very efficient in obtaining the global solution for the problems created through the procedure and parameters previously described. The reasons for this good performance are the relatively narrow ranges of variation of the processing rates and the minimal idle capacity used. These factors cause the bilinear terms to be almost constant and the problem to be nearly convex, which is the main condition for global optimality with DICOPT. In fact, previous experience in random problem generation shows that DICOPT tends to fail only for problems that violate the assumed real-world conditions, that is, for problems with large idle capacities, very long cycle times, and ranges of processing rates. This situation can be found in the illustrative example. Despite the accuracy of DICOPT, the global optimization algorithm plays an important role in validating the results from the local solver. Table 12 also shows that the algorithm using the LR is able to solve all of the proposed problems, whereas the algorithm with the NLR is not able to obtain solutions for problems 2 and 3. In these problems, the global search stopped as a result of a failure during the solution of an instance of the convex relaxation. Therefore, the LR is more robust for use in the spatial branchand-bound algorithm, as previously suggested.7 Another advantage of the LR is the shorter solution time, around 50% lower than that of the NLR. The two relaxations require the same number of iterations to solve the problems. Furthermore, their relaxation gaps for the first node are almost identical, as shown in Table 13. The ranges of the production rates are narrow, so that the linear approximation of the
exponential term is very accurate. Moreover, the underestimators of Quesada and Grossmann4 for w13/w14 are redundant, as this value is maximized in the objective function. 5. Conclusions The proposed global optimization algorithm has been shown to be effective in the solution of the simultaneous problem of the scheduling and operation of a continuous multistage plant. Results obtained using the DICOPT algorithm, that is, a local solver, could be checked and verified to be global solutions for all problems generated in which the idle capacities and ranges of production rates were small. An illustrative example was used to show that a local solution can be very different from the global one in case in which the plant allows large ranges of production rates and large idle capacities. The performance of the global optimization algorithm was tested by analyzing the impacts of several parameters, namely, the number of repetitions in the bound-tightening step, the minimum relative domain contraction per iteration, and the minimum relative improvement in the upper bound for preferential branching at a chosen variable. These results reveal that there is a critical tradeoff between the number of iterations of the global solver and the CPU time, which is mostly affected by preprocessing of the bounds and the branching strategy. Moreover, two convex relaxation strategies were tested, and the linear convex relaxation exhibited shorter computational times and greater robustness when compared to the nonlinear relaxation for the solution of problems with variable sizes. In summary, the local solver DICOPT proved to be very effective in finding the global solution in plants with small idle capacities and narrow ranges of variation of processing rates, which are assumed to be realworld conditions for operating multistage continuous plants. Despite the accuracy of the local solver in the assumed real-world conditions, the global optimization algorithm plays an important role in validating the global solution results. Acknowledgment The authors acknowledge support received from FAPESP (Fundac¸ a˜o de Amparo a` Pesquisa do Estado de Sa˜o Paulo) under Grants 99/02657-8 and 98/14384-3. Nomenclature Indices i, j ) 1, ..., NP ) products k ) 1, ..., K ) auxiliary variables m ) 1, ..., M ) stages
1498 Ind. Eng. Chem. Res., Vol. 43, No. 6, 2004 Continuous Variables CT ) total transition cost Fi ) amount of raw material consumed during the cycle for the production of product i Imaxim ) maximum intermediate inventory level for product i at stage m Invim ) difference between amount produced and maximum inventory level for product i between stages m and m+1 OPim ) operating cost to process product i at stage m Tc ) cycle time TPim ) processing time of product i at stage m TSim ) start time of processing product i at stage m Wim ) amount of product i produced at stage m wimk ) auxiliary variable for problem reformulation Rim ) mass balance coefficient (inverse of yield) of product i at stage m γim ) production rate of product i at stage m δim ) coefficient for the calculation of the maximum inventory of product i at stage m Binary Variables Xim ) indication of whether productions of i at stages m and m + 1 overlap Zij ) indication of whether product i is immediately preceded by product j Model Parameters bim ) parameter for the mass balance coefficient of product i at stage m Cfi ) raw material cost coefficient for product i Cinvim ) cost coefficient for inventory of product i at stage m Cinvfi ) cost coefficient for inventory of final product i COim ) operating cost coefficient to process product i at stage m Ctrij ) transition cost between product i and product j di ) minimum demand rate of product i Pi ) price of product i τijm ) transition time between product i and product j at stage m Global Optimization Algorithm Parameters con ) minimum relative improvement in the upper bound for preferential branching at a chosen variable dim ) minimum relative domain contraction per iteration ) relative tolerance for global optimality
Literature Cited (1) Pinto, J. M.; Grossmann, I. E. Optimal cyclic scheduling of multistage multiproduct plants. Comput. Chem. Eng. 1994, 18, 797.
(2) 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. (3) Horst, R.; Tuy, H. Global Optimization: Deterministic Approaches 2, 3rd ed.; Springer: Berlin, 1994. (4) Quesada, I.; Grossmann, I. E. A global optimization algorithm for linear fractional and bilinear programs. J. Global Optim. 1995, 6, 39. (5) Ryoo, H. S.; Sahinidis, N. V. Global optimization of nonconvex NLP’s and MINLP’s with applications in process design Comput. Chem. Eng. 1995, 19, 551. (6) Ryoo, H. S.; Sahinidis, N. V. A branch-and-reduce approach for global optimization. J. Global Optim. 1996, 8, 107. (7) Smith, E. M. B.; Pantelides, C. C. A symbolic reformulation/ spatial branch and bound algorithm for the global optimisation of nonconvex MINLPs. Comput. Chem. Eng. 1999, 23, 457. (8) Viswanathan, J.; Grossmann, I. E. A combined penalty function and outer approximation method for MINLP optimization. Comput. Chem. Eng. 1990, 14, 769. (9) Zamora, J. M.; Grossmann, I. E. Continuous global optimization of structured process systems models. Comput. Chem. Eng. 1998, 22, 1749. (10) Smith, E. M. B.; Pantelides, C. C. Global optimization of general process models. In Global Optimization in Engineering Design; Grossmann, I. E., Ed.; Kluwer Academic Publishers: Dordrecht, The Netherlands, 1996; pp 355-386. (11) McCormick, G. P. Computability of global solutions to factorable nonconvex programs: Part IsConvex underestimating problems. Math. Prog. 1976, 10, 147. (12) Floudas, C. A. Nonlinear and Mixed-Integer Optimization. Fundamentals and Applications; Oxford University Press: Oxford, U.K., 1995. (13) Duran, M. A.; Grossmann, I. E. An outer approximation algorithm for a class of mixed integer nonlinear programs. Math. Prog. 1986, 36, 307. (14) Geoffrion, A. M. Generalized Benders decomposition. J. Optim. Theory Appl. 1972, 10, 237. (15) Adjiman, C. S.; Dallwig, S.; Floudas, C. A.; Neumaier, A. Global optimization method, RBB, for general twice-differentiable constrained NLPssII. Implementation and computational results. Comput. Chem. Eng. 1998, 22, 1159. (16) Adjiman, C. S.; Dallwig, S.; Floudas, C. A.; Neumaier, A. Global optimization method, RBB, for general twice-differentiable constrained NLPssI. Theoretical Advances. Comput. Chem. Eng. 1998, 22, 1137. (17) Brooke, A.; Kendrick, D.; Meeraus, A. GAMSsA Users’ Guide; The Scientific Press: Redwood City, CA, 1998. (18) Drud, A. S. CONOPTsA large-scale GRG Code. ORSA J. Comput. 1992, 6, 207. (19) XPRESS-MP Optimisation Subroutine Library. Reference Manual, release 12; Dash Associates: Blisworth, U.K., 1999. (20) Dobson, G. The cyclic lot scheduling problem with sequence dependent setups. Ann. Oper. Res. 1992, 40, 736.
Received for review September 8, 2003 Revised manuscript received January 7, 2004 Accepted January 13, 2004 IE034118A
