1466
Ind. Eng. Chem. Res. 1992,31, 1466-1477
PROCESS ENGINEERING AND DESIGN An Outer-Approximation Method for Multiperiod Design Optimization Dimitrios K. Varvarezos, Ignacio E. Grossmann,* and Lorenz T. Biegler Department of Chemical Engineering, Carnegie Mellon University, Pittsburgh, Pennsylvania 15213
This paper addresses the development of an efficient optimization method for convex nonlinear and mixed-integer nonlinear multiperiod design optimization problems. An outer-approximation-based decomposition method for solving these problems is proposed. The method is applied to multiperiod multiproduct batch plant problems operating with single product campaigns. Multiperiod models are presented for the design and future capacity expansions of such plants. Numerical results are compared with existing general solution methods such as MINOS and SQP for the nonlinear programming case and DICOPT++ for the mixed-integer nonlinear programming case. The proposed method is advantageous in both robustness and time efficiency, with savings up to 90%.
Introduction Over the past decade there has been an increased interest in the development of systematic methods for the design of flexible chemical plants (Grossmann and Straub, 1991). The motivation for this comes from the fact that in practice flexibility is usually introduced by applying empirical overdesign, a practice that does not guarantee optimality or even feasibility over the desirable range of conditions. A major class of flexibility problems, the multiperiod design problem, involves designing plants which are capable of operating under various specified conditions in a sequence of different time periods (Groasmann and Sargent, 1979;Groasmann and Halemane, 1982). Another class deals with the uncertainty involved in some of the design parameters, which is the problem of optimal design under uncertainty (see Grossmann et al., (1983)for a review). As has been shown by Grossmann and Sargent (1978)and Halemane and Groasmann (1983), this problem requires the iterative solution of multiperiod design problems which involve nonlinear programming (NLP) or more generally mixed-integer nonlinear programming problems (MINLP), where the number of decision variables and constraints can become rather large. Typical examples of multiperiod plants are refineries that process different types of crudes and batch plants that produce a variety of different specialty chemicals. Other situations in which multiperiod operation schemes appear are planta that have seasonal demands. One characteristic of the multiperiod design problem is that it involves two distinct classes of variables, the design and the state variables. The design variables, representing equipment elements such as reactor and vessel volumes or exchanger areas, are the same for all periods of operation. The state variables, representing operating conditions such as temperatures, flow rates, or concentrations, can be different for different periods of operation. The goal in multiperiod design is to decide on the values of the design and state variables so that the plant will be feasible to operate at all specified conditions while being optimal over the specified time horizon.
* Author to whom correspondence should be addressed.
The biggest problem in solving a multiperiod model lies in the fact that as the number of periods increases there is a disproportionate increase in the computational complexity in terms of both solution time requirementa and robustness. This behavior prohibits the solution of industrially relevant multiperiod problems using standard general-purpose optimization methods. It is the major objective of this paper to develop efficient techniques for solving convex versions of these problems which have applications in the design and planning of multiperiod batch plants. The paper will begin with the statement of the problem and a brief review of decomposition methods. We then proceed to the development of an outer-approximationbased decomposition method for the NLP case in which the design and state variables are assumed to be continuous. We then present an MINLP extension in which some of the design variables are restricted to 0-1 values. The application of these methods is illustrated in the multiperiod models for the design and capacity expansions of multiproduct batch plants operating with single product campaigns. Finally, numerical results are presented, drawing comparisons with solutions when no decomposition is performed.
Problem Statement For the multiperiod design problem to be addressed in this paper, it will be assumed that the plant is subjected to constant operating conditions in N successive time periods. The dynamics of the process will be neglected, since the length of the transients is considered to be much smaller than the time periods of the successive steady states. The multiperiod design problem can be mathematically formulated as a nonlinear programming problem (NLP) when the topology of the process is fixed. More generally, however, its formulation will correspond to a mixed-integer nonlinear programming problem (MINLP) in which the topology of the plant is also subject to optimization for a given superstructure of alternatives. Consider first the case in which the topology is fixed. The variables are partitioned into two categories: the vector d , of design variables, is associated with the sizing of the units and remains fixed once the design is imple-
0888-5885/92/ 2631-1466$O3.O0/0 0 1992 American Chemical Society
Ind. Eng. Chem. Res., Vol. 31,No. 6, 1992 1467 d
xi
x2
x3
x5
x4
m -
hi
CPU time (s)
81
h2
g* h3
loo0
R3 h4 g4
0
Figure 1. Block diagonal structure in the constraints of model 1for a five-period problem.
mented for all the different periods of operation; the second class are the state and control variables, vectors xi for each different period i, that can be manipulated in each period in order to meet the production specifications. Thus,the general NLP mathematical formulation becomes N
(1)
subject to
xi
E X i = ( x i E Rn(xFI d E
x i I xiu), i = 1, ..., N
D = (d E RqldL 5 d
5 du)
Note that in the above formulation the order of the periods can be arbitrary since the operation of each period is independent of its relative position in the sequence. An important characteristic of this problem can be seen in the matrix representation of Figure 1,where the variables and the constraints form a block diagonal structure. The design variables, d , are the complicating variables and interconnect all different periods of operation. More generally, if we allow for structural changes to be a subject of optimization, the problem will become a multiperiod MINLP. Here a different type of variables, the binary variables y, are introduced, and they are primarily associated with the existence or nonexistence of a unit, or in general with a decision concerning the design of the plant. Structurally, they are similar to the design variables, since they also remain fixed once the design is implemented for the different periods. In this formulation all the binary variables are treated as design variables while the state variables remain strictly continuous. Assuming that the 0-1 binary variables can be represented linearly in the model, the multiperiod MINLP problem can be formulated as N
min z = fo(d) + Cfi(d,xi)+ cTy i=l
subject to
20
40 60 80 Number of periods
100
120
Figure 2. Solution times for two small example problems solved with MINOS without decomposition (MicroVAX 11).
r
i=l
E 0
hs gs
min z = f o ( d )+ Cfi(d,xi)
i
(2)
For large industrial multiperiod design problems there are two kinds of problems that arise as the size of the problem increases. Firstly, the computational requirements for solving the NLP in (1)or the MINLP in (2)can be very expensive, so that ultimately the solution of the problem may not be achieved in a reasonable amount of time. Secondly, the complexity and the size of the problems can be such that standard algorithms for NLP and MINLP may fail. The major reason for this is that the number of variables and constraints increases with the number of periods. Despite the fact that this increase is linear in the number of periods, the computational demand increases in most cases at least quadratically. This behavior is illustrated in the solution of two small mathematical problems described in the Appendix. In both cases, by using general-purpose optimization methods to solve the undecomposed multiperiod problem, there is a quadraticlcubic dependence of the solution time to the number of periods, as seen in Figure 2.
Decomposition Strategies From the above discussion, the need for an alternative solution strategy becomes apparent, especially as the number of time periods increases. The most promising direction is one that exploits the block diagonal structure of the problem through a decomposition strategy. Note that, once the design variables are fixed, the subproblems for each period become decoupled, and hence they can be optimized independently for xi, i = 1, N . The basic idea behind the NLP decomposition strategy by Grossmann and Halemane (1982) is based on a projection restriction procedure (Grigoriadis, 1971) by which the problem at the level of design variables is reduced into one in which variables are eliminated by using the active constraints at each time period. A common characteristic in related methods is that they exploit the problem structure through an algebraic-matrix reformulation scheme. The idea in decomposition techniques that eliminate variables from active sets is currently embedded implicitly in some of the existing reduced gradient based NLP optimization packages, such as MINOS (Murtagh and Saunders, 1985). Therefore, in this work, an alternative decomposition scheme will be developed in order to address the problem at a higher algorithmic level. At this level, an existing method is the generalized Benders decomposition by Goeffrion (1972),which can be applied to NLP and MINLP problems. Here, the original problem is partitioned into two subproblems, the primal and the master, and the solution is found by iterating between them. The primal problem is the original problem in which the complicating variables, the design variables d (and the binary variables in the MINLP case) are fixed, and therefore the optimization is done with respect to the x i variables. The master problem is a dual problem in
1468 Ind. Eng. Chem. Res., Vol. 31, No. 6, 1992
which the optimization is done with respect to the complicating variables, subject to approximations obtained from the Lagrangian of the primal problem. In each iteration the master problem, due to its relaxed form, predicta a lower bound to the objective function that increases monotonically with the number of iterations, while the primal problem provides an upper bound to the objective. Generalized Benders decomposition was initially explored in this work, but the results were not very encouraging. Firstly, the master problem is an NLP subproblem of increasing size, and secondly a large number of iterations between subproblems and the master problem may be required.
Proposed Method The new solution method is motivated from the ideas of the cutting plane method for convex NLP's (Kelley, 1960) and the outer-approximation method for MINLP's ( D u m and Grossmann, 1986). The method also involves an alternating sequence of primal and master problems and is as follows for the multiperiod NLP in (1). The primal problem is the NLP with fixed values of the complicating variables d, identical to the one in Benders decomposition. The main difference lies in the formulation of the master problem which is s i m h in nature to the OA algorithm for MINLP. Here, the information from the primal problem is used to approximate the feasible region by accumulating linear approximations of the constraints of the original problem at the optimal points, x?, in each iteration k of the primal problem. The linearization of the constraints and the objective is being done with respect to both x i and d, leading to an LP master representation for the NLP on the fullproblem space. This then provides a polyhedral approximation of the continuous feasible space of problem 1where the objective is to determine new values of the complicating variables d. As for the multiperiod MINLP problem (2), the optimization for the discrete variables y will be performed as in the OA algorithm by Duran and Grossmann (1986). However, since the NLP subproblem will be solved by the outer-approximation scheme described above, additional cutting planes will be included in the MILP master problem for the MINLP. All the analyais that follows is based on convexity assumptions for all the functions involved; i.e., f i and gi should be convex as well as the relaxed inequality form of hi at the optimum, for i = 1, N. NLP Decomposition Strategy. In the case of the NLP multiperiod problem in (l),the primal and the master subproblem have the following form: primal NLP (optimize x i for fixed dk) N
min z = Cfi(dk,xi) i=l
(3)
subject to
xi = { X i E R"IXiL
I x i I XIUJ, i = 1, ..., N Since this subproblem might be infeasible for some values d k ,a more general formulation that allows for a feasible subproblem follows by relaxing the above formulation using a penalty function approach for minimizing the violation of the constraints: primal feasibility NLP (optimize x i for fixed dk) xi
E
N
N
i= 1
i=l
min z = Cfi(dk,xi)+ pCui
(3')
subject to
1
hi(&&) - u i 5 0 - h i ( d ' ~ j )- ui 5 0
- uj 5 0
xi E Xi = {xi E R"IxiL I
xi
i = 1, ..., N
i = 1, ..., N
I xiu),
ULO, uER'
where u is a scalar variable and p is a positive number that should be greater than the absolute value of the largest Lagrangian multiplier of the constraints, according to the exact penalty scheme (Han and Mangasarian, 1979). master problem OA/LP (optimize x i and d for K approximation points) min
(4)
ci
subject to N
c
i=l
+ V&ik,dk)T(Xi - z t ) + Vdfi(Xik,dk)T(ddk)] + fo(dk) + Vdfo(dk)T(d- dk) I a k = 1, ..., K
[fi(.ik,dk)
subject to T;{hi(x;,,dk)+ V&hi(xi k ,dk )T (xi - x ; )
+
- d')) 5 0
Vdh,+&"T(d
gj(.r;,,d')
+ v,gi(x:,d)Tdri
- $1 +
V&j(X:,d)T(d - d')
5
0
1
i = 1 , ...,N; k = l ) ...)K
r ( d k )+ Vdr(dk)T(d- dk) I0 xi
E x = (Xi E
R"IXiL I xi
d E D = (d E R*ldL I d I d'),
I
XiU) ci
E R'
In the above master problem the equations of the original problem (1) are relaxed as inequalities following the equality relaxation scheme by Kocis and Grossman (1987), where Tiis a diagonal matrix with elements = sign ( A I ) and Ap is the Lagrangian multiplier of the Ith equation, at iteration k. In the proposed decomposition scheme the primal problem provides an upper bound for the original objective function (provided that it is feasible), while the master problem predicts a monotonically increasing lower bound to the original objective if convexity conditions hold. The convergence criterion in this case will be the difference between the upper and the lower bound in each iteration; the two values will be within a small tolerance t after a finite number of interations. The most important feature of the LP master formulation in (4) is the fact that it forms a linearized representation of the original NLP (11,in the space of both the d and xi variables. This would suggest that successive solution of the LP master problem in which new linearizations are generated and added to the master problem can help to obtain a better representation of the original NLP problem. In this way the solution of NLP subproblems can be avoided, a t least at the initial stages of the algorithm. In particular, it is proposed that successive master problem solutions should be within a convergence tolerance {, before resorting to the solution of a primal NLP. Clearly, if a sufficiently small tolerance is specified, the algorithm would reduce to a pure cutting plane without requiring any primal NLP solution. Another step toward the development of a more efficient algorithm is the elimination of redundant constraints. This idea exploits the fact that convergence can be guaranteed
tjt
Ind. Eng. Chem. Res., Vol. 31, No. 6, 1992 1469 by dropping, under some conditions, the inactive constraints after some iterations, which results in a reduction of the computational requirements. In particular, Eaves and Zangwill (1971) have proved that the constraint rejection scheme, described below, does not affect the convergence properties of a general cutting plane method, provided that certain conditions hold for dropping the constraints. These conditions relate to a specific improvement in the objective function with respect to a separator function. The proof for the constraint dropping scheme also holds without the above provisions when the objective of the master problem is nonlinear and the problem is solved as a linearly constrained NLP (Topkis, 1982). The proposed outer-approximation/repetitivelinear programming algorithm (OA/RLP) can now be formally stated, assuming that convexity conditions hold in the multiperiod NLP problem (1). The main steps in the proposed OA/RLP algorithm are as follows: Step 1. Select vectors xi1and d’. Set K = 1; set the upper bound zu = and the lower bound z’ = -=. Select tolerances t, {. Step 2. Set zIold = 2.Set up and solve the LP master problem (4). Set new lower bound zl = a. If lzu - z’l Ie or lz’ - zIoldl Ic convergence is achieved, STOP. If lz‘ 2’Old( I{then go to step 3; else set K = K + 1and repeat step 2. Step 3. Solve primal NLP (3’). If the solution is feasible for the original problem (u = 0) set zK = z; else zK = m. Set upper bound zu = min (zK,zU).If lzu - z’l IE convergence is achieved, STOP; else set K = K + 1 and return to step 2. An additional overall convergence criterion can been introduced in step 2 of the above implementation. This involves an e satisfaction of the constraints. Convergence is achieved if Ihi(xf,dk)lI4 and gi(x;,dk) I e’’, where d and t” are vectors whose elements are small positive numbers. The overall convergence parameters c, e’, and d’, and the inner convergence parameter {, can be chosen as a small fraction of the current bounds and the value of the individual constraints. The actual number of the NLP calls could be as low as zero, if a small value is used for the inner tolerance {. In some cases it is advisable not to use a very small value. The reason is that the convergence of the L P s might be slow close to the optimum. In this case the primal NLP solution can accelerate convergence in the final steps. Typical values for { and t used in this work were in the range of 10-6-10-6 and 104-10-’ as a fraction of the absolute value of the objective function, respectively (i.e., c = iO-’((ZI 1)). From an algorithmicstandpoint, the constraint rejection scheme is applied after the solution of each master problem in step 2. If the difference between the value of the objective in the current and in the previous iteration is greater or equal to the value of a nonnegative separator function 6(x:,dk) at the current point, that is, zk 1 zk-’ + 6(xt,dk),the inactive constraints from all previous iterations are dropped. Possible choices for separator functions are gi(xi,d)or g?(xi,d) for scalar functions g, as suggested by Eaves and Zangwill (1971). In this algorithm, the separator function was chosen to be the maximum of the violated inequality constraints, 6(xi,d) maxi (gi(xi,d)J* The above method converges to the optimum if the original NLP problem (1)has a finite optimal solution and its objective function and the constraints are differentiable and convex. We can always assume in (1)that the objective function is linear by introducing an additional variable a. Also,based on the premise that all the equa-
+
tions can be relaxed as convex inequalities, (1) can be equivalently written as min a
(1’)
subject to zi(d,q) - a s 0 hi(d$i) 5 0 gidzli) 5 0
I
i = 1, ...,N
r(d) I O and the inequality constraints of (1’)can be collectively denoted as g. Since we accumulate constraints in the master problem (4), the objective function of the master problem forms a nondecreasing sequence ah that approaches the optimum from below, as shown in the following lemma: Lemma 1. The value of the objective function a of the master problem (4) forms a nondecreasing sequence ah,k = 1, ..., K, which is always less than or equal to the objective a* of the transformed original problem (1’). Proof. (i) We define the feasible region of problem (1’) and (4)(at the Kth iteration) as F = (xi,dk(xi,d)I0, x i E Xi, i = 1, ..., N , d E DJand SK = ( ~ i , d k ( ~ ? , d ~ ) V g T ( x t , d k )A(x:,dk) I 0 for k = 1, ...,K, x i E X i , i = 1, ...,N , d E DJ,respectively. In both cases g collectively denotes all the constraints. From the convexity assumptions for all the involved functions we have g(Xt,dk) + VgT(Q,dk) A(x:,dk) Ig(xi,d)
+
Therefore, from the definition of the two sets we have F SK and Sk 1Sk+’, k = 1, ..., K - 1. (ii) Since the constraints in the master problem (4) are accumulated, at iteration K we have for all k IK: z(x:,dk) + V z T ( x t , d k )A(xF,dk) IaK for k = 1,..., K where A(x:,dk) is the vector of the differences x i - x k k and d - d k ,at the optimal point of (4). Since we minimize a: aK = min(z(xt,dk)+ V T z ( x t , d k )A(x:,dk)) k
Therefore &, k = 1,...,K is a nondecreasing sequence and ax Ia*.
The convergence of this method derives from the following theorem, similar to the one for the standard cutting plane method of Kelly (Minoux, 1986): Theorem 1. If the functions g are convex and continuously differentiable, and if (1’) has a finite optimal solution, then any cluster point of the sequence (zP,dk), generated by the above method, is an optimal solution of problem 1’. Proof. If problem 1’has a finite optimal solution, then after some step K the sequence of points (x:,dk), k = 1, ..., K, is contained in a bounded and hence compact set N. Let (xil,dl)lEL(N 3 L)be a subsequence which converges to ( x i , d ?. Consider the subsequence (xil,dlJIET(L 3 2‘) of points for which a cutting plane, with respect to the active constraints, is generated. If at each step we add cutting planea for all the active constraints, then we notice that either g(xJ,dl)I 0 after some step 1 1 lo, or the subsequence (xi,djllETis infinite.
1470 Ind. Eng. Chem. Res., Vol. 31, No. 6, 1992
r(dj) + V d r ( d j ) T ( d- d j ) - u I0
In the case when (xil,dL)lET is infinite, we have for all t' E T (t'> t ) g(x,t,dt)
+ VgT(x,t,dt)A(x:,dt) I
-
-
xi E X i = (xi E R"IxiL Ixi Ixi"),
O
d E
g(x,t,d') IIlvg(xi',d')llllA(x~,~t)ll
Since 11 A(x:,dt) 11 follows that
0 and IIVg(x:,dt))ll
g(x,t,d')
-
and hence ( x / , d') is a solution of (1'). Now, if (xi*, d*) is an optimal solution of (l'),then from lemma 1 we have at every step at I a*,from which we deduce that d Ia*,which shows that ( x / , d? is an optimal solution of the minimization problem (1'). MINLP Decomposition Strategy. In the MINLP case, the decomposition strategy presented in the previous section can be directly applied to the NLP phase of the OA algorithm by Duran and Grossmann (1986). However, since additional linearizations are generated through the proposed OA/RLP scheme, these cutting planes will be also used in the master problem of OA. More specifically, the original MINLP problem (2) is initially decomposed into the NLP subproblem where the binary variables y are fixed, and the mixed-integer linear programming (MILP) maeter problem which supplies the new binary vector. The NLP subproblem is solved with the algorithm stated in the previous section. For the MILP master problem all the linearizationsgenerated from the NLP subproblem are included. The advantage in this modified method comes from the fact that with the primal NLP that is solved with the OA/RLP approach, the additional linearizationsin the MILP master provide a tighter lower bound, and therefore a better representation of the original problem at each iteration. It should be noted, however, that this will not always translate into computational savings due to the increased number of rows in the MILP. The subproblems and the master problem are as follows: primal problem-feasible NLP NLP phase (optimize x i for fixed dj and y k ) N
min z = fo(dj) + Cfi(d',xi) + cTyk + pu i=l
(5)
subject to
Xi
hX8,xi) + Ayk = 0 i = 1, ...,N gi(d'&i) + Byh - u 5 0 dB)-uSO X E i = (xi E RnlxiLIx i Ixiu), i = 1, ..., N
d E D = (d E RqldL Id Idu) u 1 0 ,u E R1
LP phase (optimize x i and d for fixed y k ) min z = a + pu subject to N
(6)
+ VJi(Xi)T(Xi- x,) + Vdi(d')T(d - d')) + fo(dj)+ VJo(dJ)T(d - dj) + cTyk I cy j = 1, ..., Kink
C(fi(X,,dj) i=l
subject to
T:(hlCx;',d') + V & ~ ( X / , ~- ) ! z ) ~ +( X ~ VdhiF/,&T(d - d')]+ Ayk - u I 0 g,(x;',d) + V&(zj,d')T(xi - x i ) + V & i ( ~ / , d ' ) ~ -( dd')+ Byk - u 5 0 i = 1,
...)
i
1, ...,N
D = ( d E RqldL I d I d') uL0, uER1
llVg(x[,d?ll, it
g(x[,d? I0
j = 1, ...,Kink
In the above implementation both phases of the NLP subproblem are relaxed by using a penalty function scheme in order to prevent infeasibilities that may occur for some of the binary vectors yk. master problem-MILP (optimize xi, d , and y ) min z = a
(7)
subject to N
+
C(fi(~;l,d') Vj i ( X i ) T ( X i - x i )
i=l
+ Vdi(d')T(d - d')) +
+
fo(d') + Vd0Jo(dj)T(d- d') cTy 5 cy j = 1, ..., Kink; k = 1, ..., KO,,
subject to
Ti(hl(xi,d') + V & i ( ~ / , d ) ~ (-XI) xi + vdhi(Z/,d)'(d - d')]+ AY S O g,Cx;i,d) + vzg&4x/,d')T(xi- x j ) + Vdgi(~:,&~(d- 8)+ By < 0
I
k i = l , ..., N, j = 1 , ..., Kin; k = l , ..., KO,,
r(dj) + Vdr(dj)T(d- d') I O
j = 1, ..., Kink;k = 1, ..., KO,,
k = 1, ...,KO,,
C yi - C yi ,< lBkl - 1 iEB'
xi E Xi
i€M
(xi E R"IX? I X ; I xiu),
i = 1, ...,N
d E D = (d E RQldLI d IduI Y E y = (o,l)m
Bk = (ily; = 1) and Nk = lily; = 0) In the above scheme the primal problem provides an upper bound for the original objective function (provided that it is feasible) since it is a restricted form of the original MI" (2). The w t e r problem provides a monotonically increasing lower bound to the original objective provided that the convexity assumptions hold. Since the original multiperiod problem is a convex MINLP, the convergence criterion in this case will be the crossing of the bounds, provided that after each master iteration an integer cut of the previous point yk is added to the master problem, as indicated in problem 7. Another option in initializing the solution procedure, which was actually used in this algorithm, is to solve a relaxed version of the original MINLP (problem 8) through the OA/RLP scheme, in the first iteration, by relaxing the integrality requirements on y . This exploits the efficiency of OA/RLP for big problems and avoids the problem of having to specify initial values of the binary variables (see also Viswanathan and Grossmann (1990)). This relaxed MINLP subproblem has the form N
1
N; j = 1, ..., K,"
k
min z = fo(d) + Cfi(d,xi)+ cTy i=l
subject to h,(d,~J+ A y = 0
g i ( d j i )+ By r(d)I 0
5
0
]
= 1p ...)
(8)
Ind. Eng. Chem. Res., Vol. 31, No.6, 1992 1471 xi
E X i = {xi E R"IxiL 5 x i 5 xiu), i = 1, ...,N d E D = (d E @IdL5 d 5 dul y E Y,= ly E R"I0 5 y 5 1)
The a h outer-approximation/mixed-integer repetitive linear programming algorithm (OA/MIRLP) can now be formally stated, assuming that all convexity conditions hold. The main steps in the algorithm are as follows: Step 1. Select vectors x:, d', y' and solve the relaxed MINLP (8)by using the OA/RLP approach. If y is integer, the solution is found, STOP. Otherwise, set Kout= 1, zy =
-.
Step 2. Set up the MJLP master problem (7)and solve to fmd the integer vector yK with objective value zKa. Set zI = zkwt. If 2' 2 z" optimum is found, STOP. Step 3. Solve the NLP subproblem for yKoUt via the sequence OA/RLP of (5) and (6)to find a point (xiKwt, dKout, yKm) with objective value zKout. If the problem is feasible (u = 0) set z" = zKwt. If z1 2 z' the optimum, z", is found, STOP. Otherwise go to step 2. The convergence criterion is the croeaing of the bounds, since the problem is convex and at each iteration an integer cut is added to the master problem. One interesting feature of the above algorithm is that the direct solution of NLP problems in step 3 is not needed. Instead, the NLP subproblems can be solved by a series of LP's, provided that the inner convergence tolerance (is small enough, as indicated previously in the paper. In this scheme, phase one of the NLP subproblem (5) is eliminated. This particular scheme is not limited to multiperiod problems but can be generally applied to any MINLP problem. The constraint rejection step is also utilized in this implementation, in the form earlier discussed, resulting in smaller size MILP's. It has been noted that this step does not affect the number of major iterations, while it significantly reduces the solution time requirements. From a theoretical standpoint the above algorithm will converge, provided that the continuous functions are differentiable and convex and there is a finite optimum to the original problem (2). Convergence properties become identical to the ones stated in Duran and Groeamann (1986). since by theorem 1the individual NLP subproblems will converge to the optimum, and the MILP representation has at least all the linearizations of the original OA form.
Remarks On the basis of above OA/RLP decompoeition method, there are some alternative schemes that mbe proposed. Firstly, the master problem of OA/RLP can be formulated in a slightly different way, namely, as a linwly constrained NLP; linearization could be done only in the space of the constraints while the objective remains nonlinear. In that ca8e linear convergence can be guaranteed and the inactive constraints can be dropped after each iteration (see Topkis (1982)), without the conditions of a separator function stated in previous section. Secondly, the new point (xp, dk) at which linearization is performed in the repetitive scheme could be found in a slightly different way. Instead of using the solution of the previous master problem, a linear combination can be used based on the current solution point of the master problem and a strict interior point provided by the most recent feasible primal subproblem. This scheme might - further improve the speed bf convergence. At the imdementation level. on the other hand. the successive L P s can be solved tdrough the dual LP problem, in which case the addition of linearizations at each
0
2
4
6
8
10
Expansion timer (yrs)
Figure 3. Demand prediction over a 10-year horizon.
iteration would be more efficient, since the simplex tableau would be dual feasible and hence fewer pivots would be required. Application t o Design of Batch Processes ks an application of the algorithm to multiperiod design, we will consider multiperiod batch plants which produce different products with varying demands over different time periods. The key problem in the design of these plants, which stems from the nature of their products (pharmaceuticals, food products, etc.), is the changing pattern of the demand over a finite horizon time. Often, demand forecasts for a product are given over the next 1, 2, or 5 years. Planning for the expansion of production capacity in the initial design stage is then required to satisfy the predicted demand changes so as to minimize the discounted investment costs. Here the issue is to consider the tradeoff between the economies-of-scale savings of large initial capacities versus the cost of installing the capacity before it is needed. The major decisions in capacity expansion problems are the expansion sizes and the expansion times. The discounted cost of all expansions depends on the expansion cost function, the discount rate, and the demand growth over time. The expansion cost function is usually concave, exhibiting economies of scale. As pointed out by Luss (1982), the discount rate has a significant impact on the optimal policy, not only because of the different opportunity cost of money and inflation, but also because it must reflect reductions in expansion cost due to technological innovations. The MINLP model, by Gmsmann and Sargent (1979), for the optimal design of multiproduct batch plants with fixed demands, involves several processing stages with possible parallel units in each of them. In this work, this model will be extended to include different periods of operations to account for the demand changes during each period. The problem of expansions has been previously addressed through a singlestep expansion model (Wellons and Reklaitis, 1989). and where the time of expansions is undetermined. In this paper, a model will be developed to allow for several expansions to occur (or not) at different prespecified times, so that the optimal expansion strategy over a finite horizon can be determined. As an example, Figure 3 shows the case when fluctuation of the demands is considered with 20 time periods over a 10-year horizon. The initial installation and expansions are considered in 2-year intervals. Another extension of thii model toward a more realistic approach involves simultaneous optimization with resped to the amount of production. Traditionally, in batch plant models of this type, the objective function involves only the investment cost which is minimized for fixed amounts of production. In this paper we will consider for the ob-
1472 Ind. Eng. Chem. Res., Vol. 31, No. 6, 1992
jective function that additional revenues can be obtained if there are time periods in which there is an excess of capacity in which the production can be increased beyond the specified demands. The inclusion of such a revenue function with variable production amounts in the multiperiod model, however, can be shown to introduce a concave term in the objective function which cannot be convexified. In this work, a special convex revenue function will be considered for the objective in terms of the slack time for production, and which assumes a uniform production increase for all products at each time period. Two multiperiod multiproduct batch plant models will be developed with single product campaigns. In the first model the goal is to determine optimal vessel sizes Vj,for a multiperiod and multiproduct plant with fixed topology consisting of M one-unit stages, producing N products all in T different periods of operation. This leads to the following multiperiod nonconvex NLP which, however, can be convexified with an exponential transformation and effectively solved with the proposed OA/RLP algorithm. model (M.1): (a) net cost (investment cost - revenue for additional production) M
min z = CajVF - Cptj=l
t=l
Ht
(c) production horizon time for each period N QitCit
+O,IHt
(b) volume for each stage Vj h SijtBit i = 1,N, j = 1, M; t = 1, T (c) cycle time for each product at each period Nj'Cit h tijt 7 = 1, II; i = 1, N, = 1, M, t = 1, T; 2'2 II and t mod(7) = 0 (d) production horizon time
H,- 0,
(b) volume for each stage i = 1, N j = 1, M; t = 1, T Vj h SijtBit
c-Bit
are denoted by 7 = 1,..., II. The planning problem is then to determine the optimal number of parallel units Nj' at each potential expansion period 7 ( 7 = 1, ..., II), as well as the volumes Vjfor each unit which remain constant with the time periods due to the fact that they operate out of phase. One assumption concerning the parallel units for each stage is that they are of the same volume. model (M.2): (a) net cost (investment cost + expansion cost - revenue for additional production) M n~ min z = CajlNjlV,aiC Caj'(N; - Nj'-')V;8j + ;=l r=2 j = l
t=l,T
(e) number of parallel units K
N;' = CkYjk' k=l
j = 1,M;
7
= 1, n
(f) single choice of Nj'
i=l
K
CY$' = 1
(d) bounds
j = 1,M,
7
= 1, II
k=l
Vju IV j IV t
j = 1, M
i = 1, N t = 1, T Bitu IBit 5 B i t t = 1, T 0 I0, 5 6,Ht where the variables of the problem are the volume of each unit, V , ,the batch size for the product i in each period, Bit,and the slack production time in each period, 0,. The parameters of this problem are the cost coefficients for the units aj and Pj, the units size factors Si.t,the cycle time for each product at each period Cit,the demand for every product at each period Qit,the sales marginal profit at each period p t , and the horizon time for each period Ht. As noted before, it is assumed that the production of the different products can be increased uniformly in proportion to the demands Qit during the extra production time 0,. Therefore,p t is given by the marginal profit per period, based on the demands Qit,during each time period t. For realistic purposes the extra time 0, is bounded between zero and a given fraction 6, of the total horizon time for each period. The meaning of the lower bound is that the given demand Qit must always be satisfied, while the allowable increase in production cannot exceed a certain percentage of this demand. In the second model, the above batch design problem is further extended in order to include the design of the optimal topology and to address the expansion problem of the batch plant. The resulting problem is a multiperiod MINLP which can be convexified and solved by the proposed OA/MIRLP algorithm. In the following multiperiod formulation the periods of operation are denoted by t = 1,..., T, while the periods at which initial installation and expansions can take place
(g) nondecreasing number of expansion units Nj'+l h Nj' 7 = 1, I I- 1 (h) logical conditions for each expansion NjT+lIN.'I KZ'+l 7 = 1, II - 1
+
(i) bounds 1 5 Nj' 5 Njru
j = 1, M;
7
= 1, II
i = 1,N t = 1, T C i t 5 Cit 5 Citu Yjk' = 0, 1 j = 1, M; 7 = 1, II; k = 1, K Z'=O,l 7=2,II Here the additional variables of the problem are the number of parallel units for each stage at each potential expansion period 7 , N;, the cycle time for each time period Cit which is a variable in this model, the binary variable indicating the number of parallel units at each stage Yikr (k = 1, ...,K ) , and an additional binary variable 2' associated with the decision of an expansion at each potential expansion period. The parameter tijtrepresents the processing times for each product i at stage j at time period t. The parameter y' represents the fixed cost associated with every expansion. The values of aj' and y' implicitly involve the discount rate. They are based on given present values a; and y 1 and an annual discount rate of 10%. Example Problems Example problems will be presented to illustrate the application of the proposed decomposition techniques. First the convexifed NLP multiperiod multiproduct batch
Ind. Eng. Chem. Res., Vol. 31, No. 6,1992 1473 Table I. Data for Example 1 Case I: Products, N = 5; Stages, M = 3 a.= 250 ($) Bj 0.6 = 250 Vj” = 25000 Ht = BOOO/T Q a t E (210000- 290000) Qct E (160000- 200000) QEt E (100000 - 150000)
cost coefficients bounds on volumes (L) horizon time (h) production ranges (kg/year)
rf,.
j=l,M j=l,M t=l,T Q B ~E (110000 - 195000) QDt E (130000 - 185000)
A
1.9 0.7 0.7 4.7 1.2
B C D E cycle time cit for product i (all periods t ) (h)
addnl production ranges (kg/year)
t=l
(b) volume for each stage vi L In (Sijt) bit i = 1,N
+
H,- e,
j = 1,M, t = 1, T
(c) production horizon time N
CQitCi,exp[-bit] + 8, IH, i=l (d) bounds In [Vju] Ivi IIn [VjL]
t = 1, T
j = 1, M
i= 1, N, t = 1, T b , ~ , t = 1, T
In [Bitu]Ibit I In [Bi>]
o Ie, I
B 6.5
Case 11: Products, N = 10; Stages, M = 6 Q U t E (140000 - 320000) QCCt E (110000 - 210000) Q m i E [115000-225000)
plant model will be solved for a small and a medium size plant, over several time periods in each case. Then the convexified MINLP model will be also solved for the two cases and for different operation and expansion periods. The models were developed and solved within the modeling system GAMS (Kendrick and Meeraus, 1985). The package used for LP’s and NLP’s was MINOS, for MILP’s it was SCICONIC, for the undecomposed NLP’s (additionally) it was the academic code SQP (Biegler and Cuthrell, 1985),and f d y for the MINLP it was the academic code DICOPT++ (Viswanathan and Grossmann, 1990). Example 1. In this example the NLP case will be considered. First we convexify the model through an exponential transformation as suggested by Kocis and Grossmann (1988). The idea is to express the nonconvex product terms as the sum of exponential functions which are convex. This requires the definition of the transformed variables vi = In [Vi], bit = In [Bit]. Using the above and performing the necessary transformations yields the following convex NLP model: model (M/C.l): (a) net cost (investment cost - revenue for additional production) M Ht min z = Caj exp[fiivj] - C p , J-1
A 8.3
Note that the second summation in the objective function is only a function of 8, and that each term is convex in this variable. Case I. It is assumed that the plant consists of three processing stages (M= 3) with a single unit each and that it produces five products (N= 5). Its operation will be designed for a series of different time periods, namely, from 1(T= 1) up to 20 (T= 20). The total horizon time will be 1year. Table I contains the data for this example. The
2 2.0 0.8 2.6 2.3 3.6
1
stage size factor Sijt for product i at stage j (same for all periods t ) (L/kg)
3 5.2 0.9 1.6 1.6 2.4 E 4.2
D 3.5
C 5.4
E (120000- 230000) Q D D ~E
(125000 - 270000)
Table 11. Comparison of the Optimal Multiperiod Solution and the Worst Case Solution for a Five-Period Problem of Case I (a) Optimal Multiperiod Design stage 1 2 3 volume (L) 3658 2140 2408 period 1 2 3 4 5 batch size for A 463.0 463.2 463.2 463.2 463.2 each B 2675.5 2675.5 2675.5 2675.5 2675.5 product C 823.2 823.2 823.2 823.2 823.2 (kg) D 778.4 778.4 778.4 178.4 778.4 E 594.6 594.6 594.6 594.6 594.6 5.3 8.6 14.9 0.6 0.0 addnl production (% ) net cost ($) 59 820 (b) Worst Case Design 1 2 4010 2346 1 2 3 A 507.6 507.6 507.6 B 2933.0 2933.0 2933.0 C 902.5 902.5 902.5 D 853.3 853.3 853.3 E 651.8 651.8 651.8 addnl production (%) 15.5 19.0 26.0 62 180 net cost ($) stage volume (L) period batch size for each product (kg)
3 2640 4 5 507.6 507.6 2933.0 2933.0 902.5 902.5 853.3 853.3 651.8 651.8 10.3 9.0
Table 111. Problem Sizes for Cases I and I1 in Example 1 no. of total no. of total no. of periods variables constraints Case I 5 33 81 10 63 161 20 123 321 30 183 481 5 10 20 30
Case I1 61 116 226 336
306 611 1221 1831
optimal solution for a typical five-period problem is shown in Table IIa. The number of constraints in this problem is 16T + 1and the number of variables 6T + 3, shown for typical cases in Table 111. The computational results and a comparison with other methods are listed in Table IV. Case 11. In this case the plant consists of 6 processing stages (M = 6) with a single unit each and produces 10 different products ( N = 10);its operation will be designed
1474 Ind. Eng. Chem. Res., Vol. 31, No. 6, 1992 Table IV. Computational Results for Example 1 MINOS CPU OA/RLP CPU SQP CPU no. of time" ( 8 ) time" (s) time" (s) periods Case I 5 7.1 5.6 5.3 10 34.6 14.6 11.1 20 1140.5 38.9 14.6 30 3566.8 69.2 36.4 5 10 20 30
43.6 279.5 1579.1 C
Case I1 23.0 63.2 b d
11.2
30.2 84.8 222.6
"CPU time results on a 6320 VAX mainframe. *MINOS failed from the given initial point. c~~~ failed due to working set size limitations. MINOS failed to solve.
for a series of time periods ranging from 1 ( T = 1)up to 20 (T = 20), over a 1-year horizon. The number of constraints in this problem is 615" + 1 and the number of variables 115"+ 6. Input data and size information are listed in Tables I and 111,while the computational results are shown in Table IV. The use of the proposed multiperiod NLP model (M/ C.1) for demand variation yields an optimal design and production plan. The advantage of a multiperiod design versus a worst case approach is shown in Table 11. The worst case design approach results from considering a single period design consisting of the maximum demand for each product over the different time periods. Although this scheme guarantees feasibility for all potential periods, it is generally more expensive. As shown in Table IIb, the worst case design has an investment of $62 180 vs $59 820 of the optimal multiperiod design. Depending on the actual value of the marginal profit of overproduction p,, the worst case design gives a higher net cost ranging from 5% to 20%. Also, depending on the value of p,, there are two different trends in the optimal multiperiod solution. If the value of the marginal profit p , is low, there is one bottleneck period in which the additional production is zero, 8, = 0, which means that in this period it is not worth producing more that it is suggested. In this case the batch sizes are the same for all different periods due to the fact that the horizon constraint is active in all periods (see Table 11). If the value of the marginal profit p , is high, then the additional production 8, is driven to its upper bound; hence there is not bottleneck period and the batch sizes are different for each time period. Another characteristic of the multiperiod solution is the fact that, apart from the additional production in the above case, no variables lie at their bounds at the optimal solution. The sizes of the different problems with respect to the number of variables and constraints, and the respective trends as the number of periods increases, are shown in Table 111. From the computational standpoint the performance of the proposed algorithm is compared to the case when a reduced gradient method (MINOS, Murtagh and Saunders (1985)) and a successive quadratic programming method (SQP, Biegler and Cuthrell (1985)) were applied directly to these methods without performing decomposition. As can be seen from the results in Table IV, OA/RLP outperforms MINOS and especially SQP. CPU time savings were achieved in all cases, particularly as the number of periods increases. Moreover, as it is seen in the case I1 problem, which is larger, MINOS failed to solve the 20-period problem from the same starting point, succeeding from a "better" one, and it failed on all problems with more than 20 periods. On the other hand the proposed OA/RLP method was successful in all different sizes
of the problem in both cases, showing a remarkable robustness. Another important characteristic of this method is the fact that the number of the inner LP calls is constant (five and six for case I and case 11, respectively) and independent of the number of periods and hence the size of a particular problem. For the first case no NLP subproblem solution was required, while in the second case only one NLP subproblem was solved. Example 2. In this example the MINLP batch plant design will be considered. First we convexify model (M.2) through an exponential transformation. This requires the additional definition of the transformed variables uj = In [ V;], bit = In [Bit],n ,r = In [ N j r ] cit , = In [Cit]. Using the above yields the foliowing convex MINLP model: model (M/C.2): (a) net cost (investment cost + expansion cost - revenue for additional production) l l M
min z = C C(ajr-l- ajr)exp(nj'l r=2 j=1
+ Pjuj) +
Cajl M exp(njl + Ojuj) + Cn y r Z r - C p , - Ht 7x2 t=l H,- e,
j=1
(b) volume for each stage i = 1, N; j = 1,M; t = 1, T u; I In (Sijt)+ bit (c) cycle time for each product at each period njr + cit 1 In (ti;,) 7 = 1, II; i = 1, N 3 II and t mod(7) = 0 j = 1, M; t = 1, T; T (d) production horizon time N
CQ, exp(ci, - bit) + 8,IH,
t = 1, T
i=l
(e) number of parallel units K
C l n (k)Yjkr = njr k=l
7
= 1, T; j = 1,M
(f) single choice of NjT K
CYjkr=l
j=l,M; 7 = f , n
k=l
(g) nondecreasing number of expansion units n,r+l I njr 7 = 1, I I - 1; j = 1, M I (h) logical conditions for each expansion njr+l I n,' + In ( K ) Z ' + l 7 = 1, II - 1; j = 1, M (i) bounds 0 Injr IIn [Njru]
i = 1,N; 7 = 1, n In [ C i p ] Icit I In [C,,"] Case I. The plant consists of three processing stages (M= 3), and each of them may consist of up to five parallel units ( K = 5). It produces five products (N = 5), and ita operation will be designed for given forecasts of demands over a series of different time periods, namely, from 2 (T = 2), up to 20 (T = 20). The total time horizon will be 10 years, and expansions will be allowed to occur over a variety of different expansion periods (311 = 2,3, 5). Table V contains the data for this example, while Table VI11 contains information on the sizes of several problems for cases I and II. The number of constraints in this problem is 16T + 27311 - 5, the total number of variables is 6T + 24II + 3, and the number of binary variables is 16II - 1, where II is the total number of possible expansions. The
Ind. Eng. Chem. Res., Vol. 31, No. 6, 1992 1475 Table V. Additional Data for Example 2 Case I Products, N = 5; Stages, M = 3 aj' = 250 ($) y' = $20000 r = 0.1 Nj'u = 5 j=l,M Ht= 10(8000/r) t-l,T 10 years Q AE ~ 1100000 - 580 000) Qct E 180000 - 6400001 QEt E (120000 - 575000)
cost coefficients annual discount rate max no. of parallel units horizon time (h) total horizon production ranges (kg/year)
Case I1 Products, N = 10; Stages, M = 6 VjL = 250 Vj" = 45000 QM~ E { 140000 - 450 000) Qcct E (110000 - 600000) Q E E~ { 140 000 - 510 0001
bounds on volumes (L) addnl production ranges (kg/year)
aj = 0.6
j = 1, M
QBt E 150000-6800001 QDt E (160000 - 635000)
j=l,M QBet E (150000 - 595000) Q D D ~ E (150000 - 575000)
Table VI. Optimal Solution for a Five-Expansion Twentu-Period Batch Plant stage 1 2 3 2 2 parallel units for each 1 (OP 2 2 3 3 expansion period" 2 (2) 3 3 4 3 (4) 3 3 4 4 (6) 3 3 4 5 (8) volume per unit (L) 24288 18381 15987 period 1 5 10 15 20 batch size for A 3074 3074 3074 3074 3074 each product B 17763 17763 17763 17763 17763 7070 6978 7070 7070 (kg) C 620 D 4222 5168 2929 5168 5168 E 5106 5106 5106 5106 5106 addnl production (%) 33.3 14.8 33.3 18.5 0.0 OIn parentheses the actual year of operation is shown. 655500 934 800
Table VIII. Problem Sizes for Cases I and I1 in Example 2 no. of no. of total no. of total no. no. of binary expansions periods constraints variables variables Case I 2 2 81 63 31 111 31 2 10 209 2 20 369 171 31 3 3 115 93 47 5 5 200 153 79 5 10 280 183 79 5 20 440 243 79 10 10 425 303 159
Table VII. Comparison between a Period-by-Period Expansion Design and the Proposed Optimal Expansion Design (a) Optimal Expansion Design staae 1 2 3 parallel units for each 1 1 1 1 expansion period 2 2 2 3 3 3 3 4 volume per unit (L) 22471 16950 14794 net coat ($) 655 500
The optimal solution of a 20-period batch plant with five potential expansions over a 10-year horizon is listed in Table VI. Both the expansion model (M/C.2) and a period-by-period capacity expansion model were investigated in a problem with three possible expansion periods. In the period-by-period design approach the optimization is being done for each individual expansion period separately, while in the proposed model the expansion decision is made simultaneously with the optimization of the other design variables. The results from a comparative study for a three-period plant with three potential expansion periods over a 10-year horizon are presented in Table VII. As seen in this table, the period-by-period design gives a suboptimal solution that results in a substantially higher net present value cost ($934800 vs 655 500), which represents an increase of 42% when compared to the design through the proposed expansion model. Note here that even though the two approaches predict the same number of parallel units for the initial installation, the corresponding volumes for each stage are different. In the myopic period-by-period approach, the predicted volumes are smaller, since only the current demands are considered and future conditions are not taken into account for the initial design. Hence, these smaller volume units result into larger, and therefore more expensive, capacity expansions (in terms of number of units), as seen in Table VII. The difference in the objective function for the two design approaches, in this example, comes from the two factors that form the net cost according to the scheme: net cost = investment coat - revenue from additional production. As seen in the example above, the proposed model is superior to a period-by-period design not only because of the economies of scale, but also due to the additional production that can be generated. As can be seen in Table
(b) Period-bv-Period Expansion Design stage 1 2 parallel units for each 1 1 1 expansion period 2 5 6 3 9 9 volume per unit (L) 7911 5603 net coat ($) 934 800
3 1 7 12 5208
optimal solution for a typical problem is listed in Table VI and a comparison case is listed in Table VII, while computational results are listed in Table IX. Case 11. In this case the size of the problem is larger consisting of six processing stages (A4 = 6), and each of them may consist of up to five parallel units. It produces 10 products (N= lo), and its operation will be designed for a series of different time periods, namely, from 2 (T = 2) up to 20 (T= 20). The total time horizon will be 10 years, and expansions will be allowed to occur every year. The number of constraints in this problem is 612' + 84n - 11,the number of variables is 112' 77n+ 6, and the number of binary variables is 61n- 1,where lI is the total number of possible expansions. Computational results are listed in Table IX.
+
2 2 2 5 5 5 10
2 10 20 5 10 20 10
Case I1 279 784 1377 714 1019 1629 1439
182 270 380 446 501 611 886
121
121 121 304 304 304 609
1476 Ind. Eng. Chem. Res., Vol. 31, No. 6, 1992 Table IX. Comwtational Results for Examde 2 parallel units stage 1 2 3 Case I 2 2 3 2 2 3 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 2 2 2 2 3 2 2 3 1 1 1 2 2 3 2 2 3 2 2 3 2 2 3 2 2 2 2 3 3 3 3 4 3 3 4 3 3 4 1 1 2 1 1 2 1 1 2 2 2 3
no. of no. of exDansions Deriods 2
10
3
3
5
5
5
10
5
20
10
10
DICOPT++
CPU time" (s) (maior iter)
OA/MIRLP CPU time" ( 8 ) (maior iter)
163.9 (2)
87.7 (2)
105.9 (2)
54.4 (2)
631.5 (3)
594.4 (3)
404.7 (2)
189.2 (2)
1982.9 (3)
1526.4 (3)
b
1584.2 (2)
.............. *
2 2 3
2
10
5
5
2 3 1 1 2 2 2
2 3 1 1 2 2 2
2 3 2 2 2 2 2
Case I1 2 2 4 3 2 2 2 2 2 2 2 2 2 2
2 3 1 1 2 2 2
3722.2 (3) C
substantial portion of the total solution time. However, overall the proposed method is more robust and efficient in solving large multiperiod NLP's and MINLP's. This performance can be partially attributed to the fact that it mainly relies on the solution of LP's and MILPs, and also because whenever an NLP subproblem is solvent, it is done on the reduced space of state and control variables and for each period independently.
Conclusions A decomposition algorithm based on outer approximation has been proposed for efficiently solving convex NLP and MINLP multiperiod design problems. The proposed methods have been successfully applied to the design of a multiperiod multiproduct batch plant with deterministic variations in the demands. For this purpose multiperiod multiproduct batch plant models were developed, whose major components include optimal productions policy and potential capacity expansions over several time periods. With these formulations, expansion policies that involve expensive overdesigns and/or large numbers of expansions are avoided. The computational performance was illustrated for the NLP and MINLP models on examples with up to 20 periods and 10 potential time expansions. The results show that substantial savings in computation time, and increased robustness can be achieved with the proposed method.
2631.6 (2)
Acknowledgment
3818.2 (2)
We acknowledge financial support from the Department of Chemical Engineering at Carnegie Mellon University.
a CPU time results on a 6320 VAX mainframe. Internal MILP solver failed. cLimit of 5000 s exceeded in the first iteration.
VII, in the optimal design there is more capacity installed in the beginning, which can generate additional production. In fact, the plant is operating under capacity in all periods of operation, but one (one bottleneck period), whereas in the other approach most of the periods are bottleneck periods. It is important to note that when solving the above type of multiperiod models the basic idea is to make an expansion decision for the current time period while accounting for the effect of future changes in the demands that are likely to occur. Thus, the goal is not to implement the specific future plant expansions, but rather to provide better support for the current decision. Such formulations are also valuable in assessing the sensitivity of current decisions to different long-term scenarios. The computational performance of the proposed OA/ MIRLP scheme is better than DICOFT++ in all of the examined cases as seen in Table IX. In the above examples the lower bound predicted at each iteration by OA/MIRLP was always higher than the one predicted by DICOFT++. The number of major iterations for the proposed method were fewer or equal to the ones of DICOPT++. This is attributed to the more accurate representation of the master problem in this method. However, since the major part of the computations are devoted to the solution of the master MILP, the results are not always showing increasing time savings as the number of time periods increases. This can be attributed to the different degree of difficultyin solving the primal NLP and the master MILP. Since the advantage of this method lies mainly in the NLP part of the problem, the difference only became more pronounced in problems where the NLP phase shares a
Appendix The two analytic examples solved to identify the computational nature of multiperiod design problems are presented here. example 1 N
min C(aix: i=l
+ bidl + cid2)+ dS2+ d4 - d5
+ d I 2+ d,2 + d? - d4 + d52I 0 + 2di2 - d2 + 3d: + d42 - d5 I 0 i = 1, ..., N
subject to -cpi -aixi
example 2 N
min C(aix? + bidl + cid2)+ d,2 i=l
d, subject to -cixi
+ d4 -
+ 4de + 4d7 + 4dg + dg - dlo
+ d12+ dz2+ d,2 - d4 + d,2 +
+ dT2+ dg2+ dg45 0 - a i ~ i+ 2di2 - d2 + 3d,2 + d42 - d5 + d7 + de + dg2+ dlo IO i = 1, ...,N d:
where the parameters are ai E { 2 , 6 ] ,bi E (-25, -101, and ci E {lo, 2'7) for different periods i. As stated under Problem Statement, x i are the variables for each period i and dl, ..., dlo are the complicating variables.
Literature Cited Biegler, L. T.;Cuthrell, J. E. Improved Infeasible Path Optimization for Sequential Modular Simulators. Comp. Chem. Eng. 1985,9, 257. Duran, M. A.; Grossmann, I. E. An Outer Approximation Algorithm for a Class of Mixed-Integer Nonlinear Programs. Math. Program. 1986,36,307. Eaves, C. B.; Zangwill, W. I. Generalized Cutting Plane Algorithms. SIAM J. Control 1971, 9, 529.
Znd. Eng. Chem. Res. 1992,31,1477-1490 Geoffrion, A. M. Generalized Benders Decomposition. J. Optim. Theory Appl. 1962,10, 237. Grossmann, I. E.; Sargent, R. W. H. Optimum Design of Chemical Plants with Uncertain Parameters. AZChE J. 1978, 24, 1021. Grossmann, I. E.; Sargent, R. W. H. Optimum Design of Multipurpose Chemical Plants. Ind. Eng. Chem. Process Des. Dev. 1979, 18, 343.
Grossmann, I. E.; Halemane, K. P. Decomposition Strategy for Designing Flexible Chemical Plants. AZChE J. 1982, 28, 686. Groesmann, I. E.; Straub, D. A. Recent Developments in the Evaluation and Optimization of Flexible Chemical Processes. Paper presented at COPE/91, Barcelona, 1991. Grossmann, I. E.; Halemane, K. P.; Swaney, R. E. Optimization Strategies for Flexible Chemical Processes. Comput. Chem. Eng. 1983, 7, 439.
Grigoriadis, M. D. A Projective Method for Structural Nonlinear Programs. Math. Prog. 1971, I, 321. Halemane, K. P.; Grossmann, I. E. Optimal Process Design under Uncertainty. AZChE J. 1983,29,425. Han, S. P.; Mangasarian, 0. L. Exact Penalty Functions in Nonlinear Programming. Math. Program. 1979,17,73. Kelley, J. E. The Cutting Plane Method for Solving Convex Problems. J. SOC. Znd. Appl. Math. 1960,8, 703. Kendrick, D.; Meeraus, A. GAMS, an Introduction; Technical Report Development and Research Department at the World Bank Washington, DC, 1985.
1477
Kocis, G. R.; Grossmann, I. E. Relaxation Strategy for the Structural Optimization of Process Flowsheets. Znd. Eng. Chem. Res. 1987, 26, 1869.
Kocis, G. R.; Grossmann, I. E. Global Optimization of Nonconvex Mixed-Integer Nonlinear Programming (MINLP) Problems in Process Synthesis. Znd. Eng. Chem. Res. 1988,27, 1407. Luss, H. Operations Research and Capacity Expansion Problems: A Survey. Oper. Res. 1982, 30, 907. Minoux, M. Mathematical Programming: Theory and Algorithms; Wiley: New York, 1986. Murtagh, B. A.; Saunders, M. A. MINOS User’s Guide; Technical Report SOL 83-20, Systems Optimization Laboratory, Department of Operations Research, Stanford University: Stanford, CA, 1985.
Topkis, M. D. A Cutting Plane Algorithm with Linear and Geometric Rates of Convergence. J. Optim. Theory Appl. 1982,36, 1. Viswanathan, J.; Grossmann, I. E. A Combined Penalty Function and Outer Approximation Method for MINLP Optimization. Comput. Chem. Eng. 1990,14, 769. Wellons, H. S.; Reklaitis, G. V. The Design of Multiproduct Batch Plants under Uncertainty with Staged Expansion. Comput. Chem. Eng. 1989,13, 115.
Received f o r review October 28, 1991 Revised manuscript received February 24, 1992 Accepted March 23, 1992
Reburning Chemistry: A Kinetic Modeling Study Pia Kilpinen,*YtPeter Glarborg,*and Mikko Hupat Department of Chemical Engineering, Abo Akademi University, SF-20520 Turku, Finland, and Department of Chemical Engineering, Technical University of Denmark, DK-2800 Lyngby, Denmark
NO reduction chemistry in natural gas (methane) reburning was studied using detailed kinetic modeling. A reaction set including 225 reversible elementary gas-phase reactions and 48 chemical species was applied to an ideal plug flow reactor, and the most important reactions leading to NO reduction were identified and quantified for a number of conditions relevant for natural gas reburning. In addition, the influence of different process parameters on the NO reduction was investigated in the reburn zone and burn-out zone, respectively. Further, comparison of the calculations to available laboratory-scale data on reburning is made. The impact of various fluid dynamic, mixing, and chemical effects-not accounted for in the calculations-on the NO reduction and the optimum reburning conditions predicted is discussed. 1. Introduction Reburning (Wendt et al., 1973)is a combustion modification method which uses fuel to reduce NO emission. About 10-20% (energy) of the fuel-the reburn fuel-is added above the main combustion zone, to form a fuel-rich zone where NO is partially reduced to N2 and to other nitrogenous species (HCN, NHJ. After the reburn zone, additional air is supplied to complete the combustion. In this burn-out stage, the remaining fixed-nitrogen species (NO, HCN, NH,) are converted back to NO and/or to N2 (Figure 1). Reburning has been studied in a number of bench and pilot scale experiments (Chen et al., 1986;Lanier et al., 1986;Bortz and Offen, 1987;McCarthy et al., 1987;Mechenbier and Kremer, 1987;Chen et al., 1988;Mereb and Wendt, 1990),and it has been shown to be an attractive method for NO reduction, particularly for combustion systems with high initial NO emission, such as pulverized coal combustion. Furthermore, it has been found that under most conditions the most effective reburn fuels are
* T o whom correspondence should be addressed. + Abo
Akademi University. *Technical University of Denmark.
the gaseous hydrocarbons, like methane. This is partly due to the ability of these fuels to readily produce hydrocarbon radicals, which can rapidly react with NO (Myerson et al., 1957;Myerson, 1974),and partly due to the fact that these fuels do not contain any organic fuel-nitrogen. The maximum NO reductions obtained have typically been in the range 50-70% (Chen et al., 1986;Lanier et al., 1986; Bortz and Offen, 1987;Mechenbier and Kremer, 1987)or even higher-about 85%-if reburning is used together with ammonia or other additive injection (Seeker, 1990). Recently, reburning has also been demonstrated in large scale in Japan, USA, and Europe. The purpose of the present work has been to study the complex chemistry occurring during the NO reduction in natural gas (methane) reburning in order to give guidance for further optimization of the process. The work was performed using detailed kinetic modeling with a reaction set of 225 reversible elementary gas-phase reactions and 48 chemical species. The primary fuel was assumed to be coal. No heterogeneous gas-solid reactions, Le., NO-char, were however considered. It is suggested in the literature that even in the case when coal is used as a reburn fuel, the NO reduction chemistry in reburning is governed mainly by the gas-phase reactions; see, e.g., Bose and Wendt (1988).
0888-5885/92/2631-1477$03.00/00 1992 American Chemical Society
