2nd. Eng. Chem. Res. 1995,34, 4378-4394
4378
PROCESS DESIGN AND CONTROL Optimal Process Synthesis in a Modular Simulator Environment: New Formulation of the Mixed-Integer Nonlinear Programming Problem Jean-Michel F. Reneaume, Bernard M. Koehret, and Xavier L. Joulia’ INP-ENSIGC, IGPIAFP, 18 Chemin de La Loge, 31078 Toulouse Cedex, France
This paper deals with t h e formulation and the solution of a class of mixed-integer nonlinear programming (MINLP) problems, applied to optimal process synthesis, in a modular simulator environment. What is the first set out is the general solution strategy in a modular simulator environment. The proposed strategy leads to a reduced optimization problem and permits rigorous modeling of both logical and unit nodes of the superstructure. However, it involves implicit relations between the interest variables calculated by the simulator and both continuous decision and binary topological variables. I n order to handle these implicit relations, a new formulation of the MINLP optimization problem is proposed. This formulation is based on the introduction of a new set of optimization variables and constraints: the “pseudovariables” and the “pseudotorn streams”. The proposed formulation leads to a reduced set of additional variables. The solution method is based on the outer approximatiodequation relaxation (ON ER)algorithm. The master problem (MILP) is constructed by partial application of the modeling decomposition (M/D) strategy. Three examples illustrate the capabilities of the process synthesizer.
I. Introduction Optimal process synthesis is still one of the most important challenges in chemical engineering (Stephanopoulos, 1981; Biegler and Grossmann, 1985). Numerous works are dedicated t o this field and can be classified into one of these three following general approaches (Diwekar et al., 1992): the hierarchical decomposition approach, (Douglas, 1988) involving the use of heuristics; the thermodynamic target (Linnhoff, 1981), based on fundamental physical principles; and the algorithmic approach, resorting to mathematical programming theories. This paper focuses on the third approach: process structure and operating conditions can be simultaneously optimized and, under some hypotheses, optimality of the solution can be guaranteed. In the algorithmic approach, a superstructure embedding a finite number of structural alternatives is first constructed. Several processes are then defined, and the best one must be chosen. The definition of the superstructure is still one of the critical points in this approach (Grossmann, 1990). In order to reduce the size of the superstructure, Daichendt and Grossmann (Daichendt et al., 1992, and Grossman, 1994) propose a preliminary screening procedure. Friedler and coworkers (Friedler et al., 1992) propose a procedure t o generate the superstructure. Once the superstructure is defined, the optimization problem is formulated. In the most general case, this problem is modeled as a mixed-integer nonlinear program (MINLP). The process topology is described by binary variables. The operating conditions, such as reflux ratio, and the design parameters, such as exchange area, involve continuous
* Author t o whom all correspondence should be addressed.
variables. Objective function and constraints are nonlinear with respect to continuous optimization variables. Solving this problem with one of the appropriate algorithms is the last step. Recently, considerable research has been devoted to the development of such algorithms. Three of the most important MINLP optimization techniques are random search (Salcedo, 1992; Floquet et al., 1994), branch and bound (Beale, 19771, and techniques based on the projection into the space of the complicating variables (i.e., binary variables). The third one generally seems to be the most efficient. The initial MINLP problem is broken down into two subproblems: First is a primal problem, where binary variables are fixed, which results in a nonlinear program (NLP). A particular structure is extracted from the superstructure, and the optimization of its operating conditions is performed. Such a problem is solved using codes such as successive quadratic programming (SQP)(Biegler and Cuthrell, 1985)or reduced projected gradient (RPG) (Pibouleau, 1985). This problem leads t o the determination of an upper bound on the optimal solution. Second is a mixed-integer linear program (MILP) called master problem. When solved by a method such as branch and bound, it leads to a new value of the binary variables and a n increasing lower bound on the optimal solution. These two subproblems are solved alternatively up to convergence (for example, lower bound greater than or equal to upper bound). The solution is guaranteed to be a global optimum only if the initial MINLP problem is convex. There are two methods based on this principle which only differ with respect to the master problem formulation: generalized benders decomposition (GBD) (Geoffrion, 1972) which involves a dual representation of the
Q888-5885l95/2634-4378$09.0QlQ 0 1995 American Chemical Society
Ind. Eng. Chem. Res., Vol. 34, No. 12, 1995 4379 master problem; outer approximation (OA) (Duran and Grossmann, 1986) with a primal representation. Even if it leads to the solution of a larger master problem, the OA algorithm is generally considered more efficient: its lower bound is greater than the lower bound performed by the GBD master problem; thus, convergence is achieved in fewer iterations. This is an important feature since NLP subproblem solution is the stage that requires the greater amount of computational effort. This is why numerous works have been dedicated to the development of this algorithm: Kocis and Grossmann (1987) introduced equality constraints (outer approximation/equation relaxation (OA/ER)). Yuan (1988) proposed a formulation where binary variables appear nonlinearly. Viswanathan and Grossmann (1990) developed the “augmented penalty/outer approximation/ equation relaxation” (AP/OA/ER) in order to limit the effects of nonconvexity. Raman and Grossmann (1992, 1993) proposed the integration of logic and heuristic knowledge. Finally, in order to reduce the size of the master problem without affecting convergence properties of the OA algorithm, Diwekar et al. (1992) suggested a mixed formulation of the master problem (GBD/OA/ EFUAP). All this research gave rise to the development of numerous codes such as APROS (Paules and Floudas, 19891, DICOPTSS (Kocis and Grossmann, 1989a; Viswanathan and Grossmann, 1990), and PROSY” (Kravanja and Grossmann, 1990). These packages are generally conceived in an equation-oriented approach. Process synthesis generates a large scale mixed-integer nonlinear program: all constraints and variables are treated at the optimization level. Then simplified models are used t o describe the different phenomena involved in the process. At the same time, modular simulators (ASPEN, HYSIM, PROII, ProSim, etc.) have been developed. Over the years, they have accumulated rigorous models. Our aim is to integrate MINLP techniques into a modular simulator environment, in order to conceive a process synthesizer capable of puting these models to good use. For this work, the modular simulator ProSim (1994)has been chosen since analytical information about gradients is easily accessible. Little work has been done on this approach except Diwekar and co-workers (Diwekar et al., 1992; Diwekar and Rubin, 1993) in the ASPEN simulator. In a modular simulator environment the solution strategy includes three hierarchical levels. The first level, called “superstructure level”, is essentially dedicated to the construction and solution of the master problem. The primal problem is solved a t the second level, called “structure level”. The third level, called “module level”, consists of open loop simulation of the process. The major problem with this solution strategy is that output variables, or interest variables, that appear in the MINLP problem are implicit functions of the continuous optimization variables, or decision variables. To circumvent this problem, Diwekar et al. (1992) propose making the objective function and the constraints explicit by adding new variables called pseudovariables. These pseudovariables are equated to the interest variables. Since the number of variables increases and gradients are calculated by numerical perturbations, the authors propose a decomposition scheme of the set of optimization variables in order t o reduce the computational load on the NLP subproblem: gradients of the objective function and constraints with respect to all the pseudovariables are analytically computed, and gradients with respect to the other
variables are obtained by numerical perturbations. Concerning the IGCC example (integrated gasification combined cycle), the authors show that this strategy entails a reduction of computational time of about 70% in comparison with the time neded if perturbations are systematically used. The main differences between this work and the work of Diwekar and co-workers are as follows. The MINLP optimization problem is formulated so that the logical nodes, allowing the superstructure construction, are treated at the module level. This involves an important reduction of the size of the initial MINLP problem and especially of the primal NLP problem, which is the most expensive stage in terms of computational time. In addition, rigorous mixer models can be used easily when necessary. Gradients are calculated from analytical sensitivities (Wolbert, 1992; Wolbert et al., 1994); therefore, the decomposition proposed by the Diwekar and co-workers scheme is not used. In the proposed formulation, interest variables are implicit functions not only of continuous decision variables but also of binary topological variables. Then, a new formulation of the master problem t o handle these implicit relations is used. This formulation is based on the introduction of new pseudovariables and “pseudotorn streams”. Pseudotorn streams are used in order to break down the superstructure into independent structural alternatives. As in the work of Diwekar and coworkers, the pseudovariables are used in order to make the objective function and the constraints explicit. A n analysis of the optimization problem is proposed in order to reduce the number of additional variables. The modeling/decomposition (M/D) strategy proposed by Kocis and Grossmann (1989b) is partially applied to the master problem construction: constraint linearizations are modified with topological variables. But, in our formulation, some constraints are modified even if they are still active. This formulation is illustrated though three examples, one of which is the toluene hydrodeallrylation (HDA) process. 11. Initial MINLP Problem: Variables and
Constraints Definition In this section, the general MINLP optimization problem for process synthesis is described. We are particularly interested in the mathematical formulation. Our purpose is to define the different functions appearing in the objective function as well as the constraints and the set of variables. A partition of the set of variables is proposed in order to introduce later the formulation in a modular simulator environment. The mathematical formulation of the initial MINLP problem can be stated as
+
min at*y fld,s) Y d.X
s.t.
B*y+ h(d,s)I0 Cay
+ g(d,s)= 0
m(d,x>= 0 D.y
+ E*[:]
Ie
(PI)
4380 Ind. Eng. Chem. Res., Vol. 34, No. 12, 1995
x=
xE { d x E Rm2,x1 Ix IXU} $’is the set of binary topological variables which define the activity of the different units in the superstructure . @ is the set of the continuous decision variables of the superstructure. It generally concerns independent variables that define the process feed streams and chosen design or operating parameters. X’is the set of continuous independent variables of the mathematical models of the superstructure units: mass balances, energy balances, equilibrium relations, etc. These models are described by the equality constraints m(d,x)= 0 (1) Several elements of the vector x arise in the objective function cf, and in the equality (g) and inequality (h) constraints. These elements are merged in the vector s, which is a subvector of x. s is the interest variables vector. sE
s={s/s E P3,s1Is Is”}
Figure 1. Superstructure for 4REAC example.
rate; xi,,, B partial flow rate; xi,3, C partial flow rate; The mathematical formulation of the MINLP optimization problem is given by (Pz):
xi,4,temperature.
+ 2000y2 + 500y3+ 800y4+ 50d1 + d,x 100d2+ s1 + s2+ s3+ s4 - 4ooS5 (PJ The operating constraints that form the set B-y+ h(d,s) min lOOOy, Y
I0
are
0.9 - s6 I0
100 - S6 I0 In this example, there is no operating equality constraints: Cay + g(d,s)= 0. The constraints that form the set m(d,x)= 0 are feed
The initial MINLP problem is supposed to be linear with respect to binary variables. Functions f (objective function), h (inequality Constraints), g (equality constraints), and m (models equations) are assumed t o be convex. The next constraints
D-y
+ E-[:]
x ~-,d2~= 0 x1.2 = 0 ‘1,3
x ~ -,d, ~= 0
Ie
are of two different types: “integer cuts” which describe topological constraints in the superstructure (for pure integer constraints, E = 0); “bound constraints” (these constraints ensure that inactive variables are equated to zero; this is particularly important during MILP master problem solution; during NLP primal problem solution, inactive variables are set equal to zero by construction). Illustrative Example. The following example (4REAC) illustrates the approach we have been adopting throughout this paper. Component C has to be synthesized. The two reactions (A B) and (B C) take place in two different continuous stirred tank reactors. For each reaction, two different reactors are available. Kinetics are supposed to be first order in reactives. We assume that reactors are isothermal and that their volumes are known. Reaction heats and pressure drops are neglected. For any given reaction, activation energy and preexponential factor vary from one reactor to another. The physical relevance of this example is not significant t o the purpose here. The problem addressed is the selection, a t each step of the synthesis, of the reactor and the determination of process operating conditions to maximize annual profit. Operating conditions can be optimized by modifying raw material A temperature and flow. Annual cost includes fixed charges, raw material costs, and utilities costs (feed preheating and isothermal reactors). Revenue is based on product C sales. A minimal production and purity are specified. The flow sheet superstructure for the 4REAC example is shown in Figure 1. Each stream of the superstructure is characterized by four variables: xi,^, A partial flow
-
=
-
logical splitters
logical mixers
+
x4 x5 - X6 = 0
reactors
d2
+ c2 exp(-(cl/dl)) = o i = 9 , 1 0 - xi,l + xipl,l = 0 i = 4, 5
- x i - 2 3 2 d, XQ Xi,J
- xi,2 ’43
+ xiw1,2= 0
i=4,5
- ’i-2,3 -0
i=9,10
xi,,- xi-2,1 --0 ’i,4
- xi-2,4
-0
i = 9,10
i = 4 , 5 , 9, 10
si - M / 1 0 0 0 = 0
i = 1,2
Ind. Eng. Chem. Res., Vol. 34,No. 12,1995 4381 Table 1. 4REAC Kinetic Data
M1
I M3
M2
c1
500
800
700
c2
4000
5400
200 000
Superstructure level
M4 900 750 000
I'
I L
Y Slruclure level
si - A?I12000 = 0
i = 3,4
outlet s5
s6
Module level (Opm Imp p r m s s simulatim)
- '11,3 = Objective lunRm
- xl1,3Id2 = 0
The following constraints form the set D-y e:
+ E.[d,sltI
integer cuts Y1 +Y2
=1
Y3 + Y 4
=1
bound constraints
inf < YlSl -
s1 < - y1sS;lp
inf < s Y3S3 - 3
y4.Xf"f
xE sE
< y3sy -
I X.1 - y4.xqUp
i = 8, 10
x={ d x E R44, x1I x IXU}
Sc X= (91s E R6,s1 I 8 I s"}
The operating cost of each reactor is computed from the reaction heat AH. The numerical values of the constants c1 and c2 are given in Table 1. c1 squares to an energy of activation and c2 to a product (preexponential factor x reactor volume).
111. General Solution Strategy in a Modular Simulator Environment 111.1. Introduction. The MINLP initial problem has been formulated. Different solving strategies are available. The two limit ones are an equation-oriented strategy, where all constraints and variables are, theoretically, treated at the same hierarchical level, and a modular strategy, which is of particular interest in this paper. In the modular strategy, the superstructure unit model equations m(d,x)= 0 are solved a t each optimization step. For fixed values of y and d, the vector x
Figure 2. Hierarchical levels of the solution strategy.
and the subvector of interest variables s are computed by the simulator. Vector x (and s) is no longer part of the independent variable set of the optimization problem. Subsequently they will be denoted by g and E. We distinguish three hierarchical levels (see Figure 2): superstructure level (MILP master problem construction and solution, NLP primal problem construction); structure level (NLP primal problem solution); module level (solution of the different unit model equations, objective function, and constraints of primal problem computation). 111.2. Superstructure Level. This level is the highest in the general solution strategy. It includes three distinct steps: construction of the master problem according to the O U R algorithm (Kocis and Grossmann, 1987); solution of the master problem using the branch and bound method; construction of the primal problem, which will be solved a t the structure level. At each iteration, the master problem is constructed by adding the constraints linearized at the optimal point of the NLP primal problem. Linearization of the objective function is also added. The fourth section will be dedicated to the master problem formulation since it is the critical point in this strategy: problems induced by implicit relations between interest variables and continuous decision and binary topological variables need to be solved. In the second step, the master problem is solved using a branch and boundhimplex method. This code (MOMILP) has been developed in our laboratory. In the third step, the optimization problem of the particular structure, described by the fixed values of the topological variables, is constructed. This step is entirely automated in our software and has required major modifications in the simulator architecture. Several tasks, usually assigned to the simulator preprocessor, are now brought back a t the executive level. Among these tasks we quote the following: the generation of the sequential calling list as existing modules only are included in the current list, and the generation of the NLP optimization problem as only existing variables and constraints are considered. Generally speaking, a continuous optimization variable is associated with either a module (operating or design parameter) or a stream (partial flow, temperature, etc.). A variable is then active if and only if the associated module or stream is present in the current structure. A constraint is active if and only if at least one of the decision or interest variables occurring in its expression is active. Now, we discuss the problem of the first iteration. The initial value of the binary variables is either assigned by the user or automatically generated by the code. In both cases, the structural consistency of this
4382 Ind. Eng. Chem. Res., Vol. 34, No. 12, 1995
initial value is tested. Otherwise, a new value is computed. Suppose that the initial value of y is generated, the corresponding NLP primal problem is solved, and the MILP master problem is constructed. Three possibilities are available: direct solution of the master problem (in this case, inactive constraints in the initial structure are not linearized, and the master problem results in a poor approximationof the nonlinear search space; a large number of NLPMILP iterations may be necessary t o find the optimal structure); NLP optimization of a minimal number of structures in order to study all the unit models of the MINLP problem; suboptimization of the non-existing units in the initial structure (Kocis and Grossmann, 1989b). In the two last cases, the first master problem results in a good approximation of the MINLP problem. For the moment, only the first two methods are directly available in the software, since the third one is specially adapted to the equation-oriented approach. The principle of the second method is to focus attention on the first iteration, optimizing complete structures, in order to avoid having to carry out a large number of iterations subsequently. Of course, the global objective is to minimize the number of NLP optimization problems. Consider a superstructure with N serial logical splitters. Each logical splitter has only two outlet streams. Suppose that there are no topological constraints linking binary variables associated with different mixers. Example 4REAC verifies all these assumptions with N = 2. Whatever the value of N , the study of all the unit models only needs NLP optimization of two structures: an initial one and its complement. The larger N is, the more efficient the method is. Generally speaking, the number of necessary structures at the first iteration is linked with the number of mixer outlet streams. Work is under way t o a more general superstructure: logical splitters are not assumed to be serial and have more than two outlet streams. 111.3. Structure Level. In the general solution strategy, the structure level is dedicated to the optimization of current structure operating conditions: NLP primal problem. At this level, the topological variables are fixed and the resulting NLP optimization problem has been constructed at the superstructure level according to the infeasible path approach (Biegler and Cuthrell, 1985). The optimization variables are the continuous decision (d)and torn stream variables (zt,). We call them “natural” torn streams in opposition to the “pseudo”torn streams that we will introduce during the implicit relation treatment. In order t o ensure the continuity of the torn streams, equality constraints are also added: optimization variables zt, are equated to output variables zt, calculated by a sequential pass through the open loop process. It is important to remember that only active constraints and variables are taken into account. At iteration k, the mathematical formulation of the NLP primal problem is given by
s.t.
Bey
+ hk(dk ,Ek ) I0
Cay
+ $(dk,sk)= 0
y
E
q={yly fixed in {0,1}”}
dkE @ = { a d E Rml,d‘ Id zt; E
z,,= {Zts/Zts
Id”}
E Rm2,zt; IZts IZt;}
sk E S= {E computed
at module level}
gtt E &. = {gtscomputed at module level} The exponent k means that only the active elements of the different vectors at the iteration k occur in the expression of the optimization problem. As is the set of independent variables (partial molar flows, temperature, pressure) which define the natural torn streams of the superstructure. A,is the set of outlet variables associated with the torn streams, which are implicit functions of the optimization variables. J is the set of the interest variables, implicit functions of the optimization variables. For the 4REAC example and y equated to [0101It,the NLP primal problem formulation is min 2800 d
+ 50d, + 100d2+ z2(d)+ S4(d)- 400g,(d) (P,)
s.t. 0.9 - S&d)I0
100 - ;&d) I0
d E Gi7 = { d d E R2, d’ Id Id”} s E S= {;
computed a t module level}
As this step is now standard in all process simulators, it is not necessary to detail it in this paper. The NLP primal problem is solved using an SQP algorithm (Biegler and Cuthrell, 1985). Here, the important point is that exact gradient information is available in the ProSim simulator (Wolbert, 1992; Wolbert et al., 1994): in every module not only process unit model equations are solved, but outlet variable derivatives with respect to inlet variables are also calculated. These inputoutput sensitivities are calculated by analytical differentiation. A chain-ruling procedure is used in order to generate the optimization problem gradients from module sensitivities. The computational time is obviously reduced in comparison to numerical perturbations: 70430% when the SQP algorithm is use& As the primal problem is the limiting step of the MINLP procedure, this feature is particularly important. For a MINLP optimization like the HDA one, more than 95% of the total computational time is used during the NLP optimization of the different structures. Therefore we can expect about the same computational time gains for MINLP problems as for the NLP ones. Note that the more the number of natural optimization variables exceeds that of pseudovariables, the more efficient the sensitivity strategy will be, compared to the decomposition scheme (Diwekar et al.,1992; Diwekar and Rubin, 1993). The general structure of the software allows an easy implementation of new codes at the different stages of the procedure: implantation of different commercial improved codes at the superstructure level (SCICONIC,
Ind. Eng. Chem. Res., Vol. 34, No. 12, 1995 4383
L Y2
Figure 3. Logical splitter. C!
Yl
Figure 4. Logical mixer.
ZOOM, OSL, ...) or a t the structure level (reduced projected gradient (Pibouleau, 1983, SRQPD (Chen and Macchietto, 1989),QPKWIK (Schmid and Biegler, 1993)). The master problem construction stage can also be easily modified in order to implant new MINLP algorithms, such as GBD, AF’/OA/ER, etc. 111.4. Module Level. In a modular simulator environment, the different model equations are partitioned into modules. Generally speaking, there is equivalence between a module and a process unit. The sequential solution of these equation groups, for fxed values of y, d, and zts, leads to the computed values of the interest variables s and the torn stream outlet variables Zt,. The objective function and the constraints of the current primal problem can then be evaluated. We distinguish three module types: process unit modules; objective function and constraints evaluation module (Muti on Figure 2) (the user has to supply information about the MINLP objective function and constraints; torn stream constraints are directly treated by the structure level; and integer cuts and bound constraints are treated a t the superstructure level); logical nodes. We are now going to detail the logical node modeling problem. These particular nodes are of two types: logical splitters and logical mixers. These nodes are used to describe the structural alternatives embedded in the superstructure. There are two cases of operating conditions. First Case. Only one outlet (respectively inlet) stream of the splitter (respectively mixer) is active. In this case the splitter (Figure 3) directly assigns inlet stream information to the active outlet stream (partial flows, temperature, pression, enthalpy, etc.).
Cl=C,
ify,= 1
C, = C,
ify, = 1
In the same way, the mixer (Figure 4) assigns active inlet stream information to the outlet stream.
C, = C,
ify, = 1andy, = 0
C, = C,
ify, = 0 and y, = 1
Second Case. At least two outlet (respectively inlet) streams of the splitter (respectively mixer) are active. It is the case of the logical mixer LM3 in the HDA example when the y4 value is 1 (see Figure 8). In this case, the logical module automatically calls on the appropriate process unit module (mixing or splitting). In this way, rigorous models are used without increasing the optimization problem size. In previous works,
constant heat capacities were assumed. In the equation-oriented formulation, this assumption was made (Kocis and Grossmann, 1989b) in order to handlenonconvexities and to limit the size of the MINLP problem. In our approach, since we do not simplify process unit models, there is no objective reason to simplify logical node models when they act as a process unit. For the logical splitter, an additional problem arises from the management of the continuous optimization variables: split fractions of the different active outlet streams. Note that the equations of the logical unit models are treated at the module level. This is a crucial point of the formulation proposed. It allows easy management of the rigorous models, as seen previously, as well as a large reduction in the MINLP (and consequently NLP) problem size. In opposition to the formulation proposed by Kocis and Grossmann (198913) and Diwekar, et al. (1992), the variables of the logical units models are not independent variables of the optimization variables. In the 4REAC example, this leads t o 44 less variables and 16 less constraints. We will see, in fact, that the reduction of the problem size is not that large because the interest and torn streams variables are, due to this formulation, implicit functions of the binary topological variables. This problem will be tackled in the next section. Finally, note that, as Figures 3 and 4 suggest, logical units are not limited to two streams. The only limitation arises from the combinatorial aspect of the MINLP optimization problem.
IV. Formulation of the Optimization Problem IV.l. Introduction. This section is dedicated to the two problems that arise from MINLP optimization in a modular simulator environment, when the above described strategy is used. The first problem is related to constraints where occurring interest variable activities are described by the same binary variable. Then, when the constraint is active, all interest variables occurring in its expression are active too. We will call them “type I constraints”. The second one is related to constraints and an objective function where occurring interest variable activities are described by different binary variables. Then the constraint can be active, even if at least one of its occuring variables is not. We will call them“type I1 constraints”. IV.2. Type I Constraints. It is important to keep in mind that the MILP master problem is constructed by adding constraint linearizations at the solution point of the successive NLP primal problems. It is solved using a branch and bound method. The study of a terminal vertex of the tree (linear programming (LP)), corresponds to the optimization of a particular structure of the superstructure. This structure is partially represented by the linearizations of its objective function and part of its constraints. These linearizations can occur at different points. Constraints associated with other structures corresponding to other terminal vertices of the tree must be removed from the problem. If not, the master problem can result in a suboptimal solution. Consider the 4REAC example. The purity constraint and the production constraint are always active and are implicit functions of the decision variables dl and dz, which are always active too. The operating conditions’ optimization of the structure described by y1 = [0101It
4384 Ind. Eng. Chem. Res., Vol. 34, No. 12, 1995
(see problem (P4)in section 111.3) leads to the following point: objective at the solution:
-13 215.1
(dl’)*= 941.3 (dZ1)*= 254.2 The master problem is constructed by adding the following linearizations:
+
-8.1 x 10-~d, 3.5 x 10-~d,I1.3 x (purity) (3)
-2.1 x 10-’d1 - 0.81d2 I-97
(production)
(4)
The NLP optimization of the structure described by y2 = [1001It leads to the following point: objective a t the solution:
-17 925.4
(d12)*= 714.1 (d;)* = 218.4 This point is the solution of the MINLP problem. But it could not be reached a t the superstructure level, since it involved violation of the linearized purity constraint obtained for the structure y’ = [0101It: eqs 3 and 4. During the master problem solution, the criterion of the structure [1001It is overestimated and the algorithm fails to find the optimum. Note that, in this example, the violation is not very significant since the reactor models are nearly the same. This result can seem to be in opposition with the ON ER algorithm principle: two linearizations of a same constraint, at two different points, do not lead to a linear search space that a overestimates nonlinear search space. We are now going to demonstrate that this problem does not arise from the convexity of the constraint but from the modular simulator environment itself. Definition of Pseudotorn Streams. Consider an inequality constraint h(d,g) I 0. d is the decision variable vector and the interest variable vector. g is an implicit function of d, but-this is the crucial point-this implicit function varies from one structure to another. In other words, s is an implicit function of the topology which describes the path leading from the decision variables to the interest variables. During the master problem construction stage, the linearization of one Constraint can actually lead to the linearization of different functions. The purity constraint of the 4REAC example is
0.9 - &(dl,d2) 5 0
(5)
If y = [1010It, Sg(dl,d2) results from the composition of the models of units M I and M3. If y = [0101It,g6(dl,d2) results from the composition of models M2 and M4. The implicit function 96(dl,d2) depends on the topology described by y. In other words :6(dl,d2) depends on the path leading from the decision variables (dl and dz) to the interest variable (56). Even if the purity constraint is always active, the functions that are linearized vary from one MINLP iteration to another. Now note the particular function of the logical mixers in the superstructure. Streams from different struc-
tural alternatives converge into these nodes. As seen previously, we have chosen to treat these nodes a t the module level: the size of the MINLP problem is reduced, and it is a simple matter to use rigorous models. But, with this formulation, interest variables become implicit functions of the binary variables. This results in the failure of the algorithm. If the outlet streams of the logical mixers are treated as torn streams, a new set of independent variables Zlm and constraints (Zlm - Zlm(d,Zk,Zlm)= 0) are introduced into the problem formilation. These constraints ensure stream continuity at the convergence point. The important point is that interest variables situated after the torn stream, with the flow, are no longer direct functions of upstream binary or continuous variables. These interest variables are simply functions of pseudotorn stream variables (Zlm) and downstream decision variables. We have already introduced the idea of a path connecting the decision variables t o the interest ones. This formulation, introduction of the pseudotorn stream, is equivalent to splitting the path into structurally independent subpaths. Figure 5 shows the calculation graph of a structure k included in a superstructure with n serial structural alternatives. The structural alternatives are figured with numbered rectangles (from 1 to n). The type I constraint, hlk, is a function of the interest variable gk, which is an implicit function of the decision variable d?, the pseudotorn stream variables Zlm,? and possibly the torn stream variables zts,$. Thus, -s?, and hlk are independent of the other structural alternatives. With this new formulation, each interest variable is associated with a unique function as in the equationoriented approach. During the master problem construction stage, linearizations are treated as in the equation-oriented approach: if the constraint is not always active, its linearization is modified by the binary variable that describes its activity. On the other hand, if the constraint is always active, its linearization is not modified; as we shall see later, this is not true for pseudotorn stream constraints. Linearization Modification. It is probably necessary t o recall what “linearization modification” means. In the equation-oriented approach, Kocis and Grossmann (1989b) propose a modeling/decomposition WD) strategy. As part of this strategy, the constraint linearizations are modified as follows. Consider the nonlinear constraint h(d) I 0. This constraint is active in the current structure, but its activity is described by the binary variable y. The master problem is constructed adding the modified linearized constraint:
Uy
+ (h(d*)+ V,ht(d - d*))IU
(6)
U is a sufficiently large scalar. Ify = 0, the linearized constraint is satisfied whatever the value of d . The linearized constraint is inactive. If y = 1, it is active. The reason why Kocis and Grossmann introduce this strategy is as follows. During the master problem solution, particularly a t each terminal vertex, inactive optimization variables are set to zero. That can lead to a constraint violation since we are out of the valid linearization area. In the equation-oriented approach, this problem only arises for constraints that can be inactive and not for the ever active constraints if they are convex. Since this strategy allows constraint deactivation, it could be used to handle the problem of implicit relations between interest and binary topological variables. The
Figure 5. Calculation graph.
.i"""'
'4
'\awn)
C4
Figure 6. General logical mixer.
Figure 7. 4REAC. Introduction of the pseudotorn streams.
4REAC purity constraint, linearized a t the optimum of the structure described by y = [0101It, is modified as follows:
of Blmk.U is a real large enough and L represents the linearized constraint. In this way, the linear constraints are active, if and only if the terminal vertex corresponds t o a structure that includes the logical mixer with the same operating conditions, that means the same active inlet streams. Illustrative Example. The strategy of the pseudotorn stream introduction is now illustrated with the 4REAC example (Figure 7). The mathematical formulation of the MINLP problem is
UY,
+ UY, - 8.1 x 10-~d,+ 3.5 x 10-~d,I 1.3 x + 2U (7)
U is a sufficiently large scalar. During the master problem solving, this constraint is active if and only if y = [0101It. At the terminal vertex [1001It, this linear constraint is then inactive and the solution point d2* can be reached. But this strategy leads to a poor linear approximation of the search space. The optimization of all structures included in the superstructure can be necessary to find the optimal one. The linearization modification is not a sufficient a strategy for handling implicit relations in a modular simulator environment. The nonlinear function itself has to be modified by introducing the pseudotorn stream. But this new formulation involves a set of equality constraints associated with the pseudotorn stream. These constraints are particularly important since they ensure process continuity between alternatives. In fact, in a general way, inlet streams of logical mixers can come from different structural alternatives or, in some cases, can remain active. Figure 6 shows a logical mixer with four inlet streams: C1 and C2 come from the same alternative and are associated with the binary variables y1 and y2 submitted to the constraint y1+ y2 I1. C3 comes from an other alternative and is associated with y3. C4 is always active. Ylm denotes the subvector of y whose elements are the binary variables that describe the operating conditions of the mixer (here ylm = b0v31t). Consider the constraints associated with the pseudotorn stream leaving the logical mixer. The current structure is described by the topological vector yk that includes ylmk.The pseudotorn stream constraints are linearized at the optimal point of the current structure. The MILP master problem is constructed by adding the linearizations modified as follows:
C JeBlrnk
UYlmJk -
C
UYlmJk + L 5 CB1mkU
(8)
J a h k
Blmkis the set of ylmkelements equal to 1. Nlmkis the set of the ylmk elements equal to 0. CBlmkis the cardinal
min lOOy, YAzlrn
+ 2000y2 + 500y3 + 800y4 + 50d1 +
1 0 0 4 + s,(d) + s2(d>+ S3(Zlm1) + S4(Zlml) 4ooS5(Zlm1)
(p5)
set.
Y3 + Y 4
=1
s E J = (2 computed at module level} qmE gm = {%m computed at module level} Note that it is not necessary to introduce a pseudotorn stream associated with the logical mixer (LMd. The functions associated with the interest variables ~5 and -S6 can be characterized by a unique binary variable (y3 or y4). At the master problem construction stage, purity and production constraints are modified by the variables y3 or y4 and the pseudotorn stream constraints (LM1) are modified by y1 or y2.
4386 Ind. Eng. Chem. Res., Vol. 34, No. 12, 1995
Some Remarks. To conclude this section, two remarks are called for. First, if the strategy proposed by Diwekar et al. (1992) is applied (constraints made explicit by adding new pseudo decision variables), the algorithm fails. We introduce a pseudovariable d3 in the 4REAC purity constraint. This constraint is now explicit with respect to the decision variables: 0.9 - d3 I0
(9)
We also introduce a new constraint:
In fact, the problem previously set out (implicit relation between interest variables and the binary topological variables) has been transferred to this constraint. The reason why this strategy fails is that logical nodes are treated at the module level. This leads to a shorter MINLP formulation, even if we introduce the pseudotorn streams. The second remark concerns the calculation of the real U. As discussed by Raman and Grossmann (19941, this modification strategy generally leads to an ill-conditioned problem: U must be a large real. Instead of taking an arbitrary large value, the lowest value of U can be easily computed. We know the linear constraint coefficients (a,)and the optimization variable bounds:
IV.3. Type I1 Constraints. In this section, we consider the functions of the optimization problem involving interest variables whose activities are described by independent binary variables. It is the case of the constraint ha in Figure 5. It is often the case of the MLNLP objective function. For instance, the 4REAC objective function involves interest variables SI,22, t 3 , and 5 4 which come from different structural alternatives. The problem comes from the fact that it is necessary, but not possible, to compute the gradient of such a function: this would involve the differentiation of an inactive interest variable. It is false to say that the gradient of the function with respect to the inactive interest variable is zero. Definition of Pseudovariables. The adopted solution is the one proposed by Diwekar and co-workers (Diwekaret al., 1992; Diwekar and Raman, 1993). Type I1 constraints are made explicit by adding new pseudovariables s. These continuous optimization variables are equated to the corresponding interest variables. Note that this strategy is not applied t o the whole set of interest variables but only to those arising in type I1 constraints. The number of additional variables is then greatly limited. The interest variables vector is divided into two subvectors:
s': elements of g involving pseudovariables -s": such t h a t = [,"'It
Vector v denotes here all optimization variables (decision, torn streams, pseudotorn streams, etc.). The lowest value of U is given by
U = ZaiZi + K
(13)
with
Illustrative Example. Using this strategy, we introduce four pseudovariables SI,s2, s3, and s4 respectively associated with the interest variables SI,92, 93, and s4 in the 4REAC example. The objective function in (Pi) becomes min lOOOy, y&,zlrn
+ 200y2 + 5ooy3 + 800y4+ 50d, + 100d2 + s1 + s2+ s3+ s4 - 400S5(Zlml)
We also introduce four new equality constraints:
For each constraint, this value is automatically computed. Note that the bounds determination is a crucial point in this method. But another problem arises. We have seen previously that if a constraint is not always active (the activity is determined by the binary variable y), the master problem is obtained by addition of the modified linear constraint: U y + X q v i + K ~U
(17)
The linear constraint is active if and only if y = 1. During the branch and bound analysis, binary variables are able to take intermediate values (strictly included between 0 and 1). For such values, the constraint activity makes no sense. This seems to be a theoretical limitation of the modification technique. In practice, for the examples that we have studied, we have never found any failure case due to the master problem formulation. That is the reason why we have not introduced another technique such as those proposed by Balas (1985) or Raman and Grossmann (1994). Note that many numerical approaches have intermediate steps that do not make physical sense.
Whatever the current structure, the objective function gradients can be computed. If y4 = 0, the interest variable 9 is inactive and the pseudovariable is equal to 0. But the function derivative with respect to s4 is known and is equal to 1. On the other hand, the ) involves no difficulties constraint (s4 - ~ 4 ( ~ 1 ~=1 0) since it is inactive: it has not to be linearized. Generally speaking, the new constraints are type I constraints and are consequently treated the way presented in the previous section: they are linearized if they are active, and their linearizations are modified, if necessary, by the binary variable describing their activities. Some Remarks. We have already seen that it is not necessary to introduce a pseudotorn stream for the second logical mixer LM2: the functions associated with the interest variables 55 and 96 can be characterized by a unique binary variable (y3 or y4). In order to avoid four additional variables and constraints, we did not do this. In the same way, it is not necessary to introduce
Ind. Eng. Chem. Res., Vol. 34, No. 12, 1995 4387 the pseudovariable s5 since the interest variable 5 5 is always active: at each MINLP iteration, objective function linearization can be modified by y 3 or y 4 , which are the binary variables that characterize the function associated with the interest variable 55. In these case, we introduce the peudo-variable s5 because the objective function is made linear with respect to the optimization variables. This particular form is very interesting: the master problem objective function is the same as the MINLP one. The master problem formulation is then greatly simplified. Addition of this pseudovariable calls for a new equality constraint.
extended formulation of the primal problem [ O l O l I . min 2800 d,zlrn
+ sod, + 100d2.+g(d)+
(pfj)
s.t.
100 - sJ(Zlrn) 5 0
- %m(d) =0
=0
For the 4REAC example, the interest variable vector is partitioned as follows:
-s’ = [51&y4S51t
-
400S5(Zlrn)
zlrn s 5 - s&lrn,)
54(zlrn)
and
g” = [g,I
IV.4. Torn Streams Treatment. As seen earlier, there are two classes of torn streams: natural torn streams (zd involved by the infeasible path approach for the solution of the NLP primal problem; pseudotorn streams ( ~ 1 involved ~ ) by type I constraints treatment. By analogy with the general solution strategy, the former can be considered as structure torn streams (operating condition optimization of a particular structure) when the latter are considered as superstructure torn streams (MILP master problem construction). During the NLP primal problem resolution, the structure torn streams are necessary since we have chosen the infeasible path approach. On the other hand, superstructure torn streams are unnecessary from the primal problem convergence point of view. Two strategies are consequently available. The primal problem is solved without superstructure torn streams. The optimal operating point of the current structure can be reached using a classical NLP optimization algorithm. This formulation has the important advantage of resulting in a reduced formulation of the NLP primal problem: superstructure torn streams are not present. This formulation will be referred to as the “reduced” formulation. In section 111.3 we give the reduced formulation of the primal optimization problem for the structure [ O l O l I . This problem only includes two variables and two constraints. On the other hand, the gradients of the objective function and the constraints, obtained at the solution, are useless for the next stage (master problem construction). As shown earlier, the different structural alternatives have t o be independent in order to construct a master problem that efficiently handles the implicit function problem. As can be seen in Figure 5 , a decision variable associated with a structural alternative has no direct influence on an interest variable associated with another alternative. This sensitivity, if it exists, is transmitted indirectly by superstructure torn streams. Therefore, if the reduced formulation is chosen for the solving of the NLP primal problem, an intermediate stage has to be introduced before the master problem construction stage: the formulation of the primal problem with the superstructure torn streams and gradient calculation at the previously reached NLP optimum. This stage will be referred to as the “structure intermediate stage” (SIS). In the second strategy, the NLP primal problem is directly constructed with superstructure torn streams. The optimization problem is then increased in size. This strategy will be referred to as the “extended” formulation. For the sake of comparison, we give the 4REAC
d E Q={d/dER2,dL5d
