On the Optimum Choice of Decision Variables for Equation-Oriented

4742

Ind. Eng. Chem. Res. 1999, 38, 4742-4758

On the Optimum Choice of Decision Variables for Equation-Oriented Global Optimization Romualdo L. Salcedo* and Ricardo M. Lima Departamento de Engenharia Quı´mica, Faculdade de Engenharia da Universidade do Porto, Rua dos Bragas, 4099 Porto, Portugal

In design problems, where the set of variables is larger than the set of equations, the difference corresponds to the degrees of freedom available to the designer. The use of equation-oriented simulators is particularly useful for the global optimization of nonconvex problems, such as those that usually describe chemical processes. This paper shows that, by coupling a combinatorial optimizer with a tearing/partitioning algorithm, the simulation step of an optimization problem can be posed as a combinatorial optimization problem, with the objective of minimizing the cost of the simulation step. The functional form in which the variables appear in the equations can easily be taken into account as constraints to the optimization problem. The concept is described in detail for several examples found in the chemical engineering literature, showing that the proposed method may be a useful preprocessing tool for the global optimization of nonconvex NLP or MINLP problems, where SQP-based methods may not be adequate. Introduction For an optimization problem, where there are always more independent variables than equations, the decision variables may be all continuous, viz., with nonlinear programming (NLP) problems, all discrete, such as in combinatorial optimization, or discrete and continuous, such as in mixed-integer nonlinear programming (MINLP) problems. Chemical processes are usually described by both linear and nonlinear functions, and the resulting optimization problems are usually nonconvex. For the global optimization of these problems, the simulation step has to be solved prior to the optimization step. This, in turn, requires that the set of equations describing the process, i.e., mass and energy balances, stoichiometry, rate and transport equations, and equilibrium and design relations, must be determinate.1,2 The simulation step is usually performed in one of two distinct modes: sequential modular or equationoriented.2-4 With sequential modular simulators, the flow of material in the actual process is paralleled by the flow of unit calculations. Here, unit and physicochemical (or thermodynamic) models remain selfcontained as distinct procedures, which are called at a higher level in order to force convergence of the stream connectivity equations that reflect the flowsheet topology.3 Although the modular nature allows for an easy expansion of the simulator’s capability, convergence problems may occur because of difficult iteration loops in the calculation sequence.5 The iteration loops are usually converged using successive substitution with Wegstein or dominant eigenvalue acceleration.1-3 Nowadays, the modular approach incorporates successive (or sequential) quadratic programming (SQP) to solve simultaneously the recycle and design optimization problems.4,6 These algorithms employ linearization techniques for the constraints and may follow infeasible paths, such that the simulation step does not have to * To whom all correspondence should be addressed. Email: [email protected]. Phone: +351-2-2041644. Fax: +3512-2000808.

be solved. However, SQP-based algorithms are local search methods,7 and linearizations may cut parts of the feasible region of the original problem.7,8 Thus, multiple solutions for the subproblems and nonconvex search regions may fail to produce the global optimum. Recently, simultaneous solution of stream connections, at the flowsheet topology level, using quasiNewton or Broyden updates,3,5 has been added to sequential modular simulators. This feature mostly distinguishes simultaneous modular strategies from the traditional sequential simulation approach.4,6 However, adding design constraints and specifications requires additional iteration loops.2,5 To reduce the burden of convergence problems, Montagna11 has suggested solving an integer optimization problem, whereby modified solution paths associated with a minimum number of tears may be obtained. Also, global search strategies based on mathematical programming techniques have been added to sequential modular simulators.9,10 With equation-oriented simulators, such as Fortran solvers, one has to write down the process equations, which will have to be solved by some appropriate method. This is usually performed with Newton-type algorithms, with or without homotopy (continuation) methods for widening the convergence domain toward global convergence, for finding starting guesses close to the desired solution, or for dealing with Jacobian singularities. For a comprehensive discussion on these methods, the reader is referred elsewhere.1-3,6,12,13 The added difficulty of having to write the entire set of equations is offset by the greater flexibility in solving nested recycle loops and in the incorporation of new unit specifications and design constraints.2,5 However, unlike with modular simulators, a good initialization scheme is usually a prerequisite to a successful full-scale solution of the entire system of equations.3 At the optimization level, equation-oriented simulators may be easily interfaced with global search algorithms.14,15 However, the decision variables must be somehow fixed or arbitrated, so that a determinate system is eventually obtained. It is clear that one has the choice of dictating which should be the decision

10.1021/ie9901980 CCC: $18.00 © 1999 American Chemical Society Published on Web 11/16/1999

Ind. Eng. Chem. Res., Vol. 38, No. 12, 1999 4743

variables. With modeling languages such as GAMS16 or AMPL,17 the choice of decision variables is transparent to the user and will depend on the nonlinear solver employed. With Fortran solvers, there are some heuristic rules for choosing these decision variables, such as discreteness, linearization properties, boundness, and obtaining a sequential solution.18 It makes sense to force discrete variables to be decision variables, which will be fixed by the optimization algorithm, rather than allowing them to be output (state) variables. Also, because linearization will make the system of equations easier to solve, then those variables that linearize equations should also be chosen as decision variables. Because the optimizers may require an initial guess for the decision variables, it also makes sense to choose bounded variables as decision variables. Although these criteria are easy to apply, they do not often produce a unique set of decision variables, even for small-scale problems; viz., many sets that meet these criteria are possible. However, some sets produce a sequential solution for the state variables, which is generally, but not always, a better path than forcing a simultaneous solution.12,19 A sequential solution may, however, not be possible, even with small-scale problems, because of a recycled flow of information. It is clear that the structure of the system equations, viz., their ordering or grouping, may have a strong impact on the solution process.3,19-23 Several techniques are available to reduce the effort in solving the nonlinear system, viz., eliminate variables in equations with only one variable and move as many nonlinearities as possible from the equality constraints to the objective function.20,24,25 In any case, the efficiency of these processes is highly dependent on the structure of the system equations. Because this structure may change drastically by specifying different independent variables for the degrees of freedom, this feature may be exploited in order to efficiently solve systems of nonlinear equations, and this is the subject of the present paper. Although this topic has been the subject of significant past studies,19-23 it remains an important one that is, to the authors’ knowledge, not yet resolved satisfactorily in the literature. Tearing/Partitioning Algorithms With sequential modular simulation, there may be recycle streams that have to be identified and computed, usually by fixed-point iteration and acceleration schemes.1-3 To aid in the identification of those recycle streams, tearing/partitioning algorithms were developed. These were aimed at finding the process units that have to be solved as a group (partitioning), their solution ordering (precedence ordering), and the identification of a minimum number of recycle streams (tearing), because the performance of fixed-point algorithms greatly depends on the tear streams. Biegler et al.3 describe a partitioning and precedence ordering algorithm due to Sargent and Westerberg26 and a tearing algorithm, posed as an integer programming problem by Pho and Lapidus,27 that minimizes the number of tear variables, subject to the constraint that all recycles must be broken at least once. Several variants of this algorithm exist, depending on the weights given to the costs of tearing the different streams, but usually several different candidates for the tear streams may give the same cost; viz., different solutions may be equally performant.

With equation-oriented simulation, because a simultaneous strategy is usually employed for its solution, there appears to be less of a need to analyze the structure of the entire system of equations. The effort has been toward the use of efficient decomposition schemes for the sparse Jacobian matrix.3 However, as stated above, SQP-based algorithms are local search methods that employ linearization techniques that may cut off the global optimum.7,8 Thus, they are not adequate for the global optimization of nonconvex functions. Whenever the simulation must be solved at each optimizer step, the problem of the difficulty in solving the nonlinear system of equations arises. Several researchers have developed efficient tearing/partitioning algorithms, which may also profit from the sparse condition of nonlinear systems that usually occur in process design. Christensen20 proposed a tearing method to reduce the number of equality constraints, whereby the number of torn variables is minimized and where the initial partitioning is avoided. This author has shown that, with complete partitioning before tearing, the minimum number of tears may not be obtained. The same has been recognized by Stadtherr et al.,22 who exploit the linearity or solvability of small subsets of equations, obtained by tearing and which they call “allowable subsets”. The basic difference between the work of Stadtherr et al.22 and that of Christensen20 resides in the concept of allowable subsets of rank greater than 1. Ledet and Himmelblau28 have devised a decomposition scheme similar in objectives to the decomposition of Sargent and Westerberg.26 Basically, the idea is to partition the system into the smallest irreducible subsystems and to identify their precedence ordering and a minimum number of tearing variables. Unlike the approach of Christensen,20 the functional form in which the variables appear in the equations can be taken into account as constrained choices of tearing and output (dependent) variables. Thus, it may be fine-tuned to a specific problem in order to minimize convergence problems and unrealistic outputs, i.e., data which exist in tabular form or multiple solutions. To find an output set, i.e., a set of output variables, one from each equation, the algorithm employs Steward’s method for partitioning,29 which is basically a logical sequence of steps to assign variables to equations. This method is also useful to detect errors in formulating a model by identifying sets of equations that contain improperly specified design variables or parameters, or dependent (redundant) equations. It is likely that if an output set cannot be found with Steward’s method, then the system of equations is not determinate.28 Tearing is accomplished by finding the minimum number of tearing variables that reduce the irreducible subsystems obtained by partitioning to subsystems of single equations. The decomposition constraints are added at this step, where certain variables may be prohibited from being tears because of ill-conditioning and/or from being outputs from certain equations, because of discreteness, feasibility to solve the equation for that variable, or presence of multiple solutions. Shacham30 has recognized that tearing may be expensive in terms of computer time, because it requires formula manipulation. Also, the smaller nonlinear systems may be more difficult to solve. This author proposes an alternative tearing method in which new

4744

Ind. Eng. Chem. Res., Vol. 38, No. 12, 1999

variables are introduced in certain equations to remove nonlinear terms and tearing be performed for the smaller linear subsets. Previous Work on the Choice of Decision Variables Several authors have recognized that by specifying different sets for the decision variables, nonlinear sets of different numerical difficulty will occur. The pioneering work of Lee et al.31 for the selection of design variables, using bipartite graph structuring, works very well if an acyclic information structure can be obtained but fails otherwise.19,22 Rudd and Watson18 present this algorithm in a more familiar form of an occurrence matrix manipulation. Book and Ramirez23 developed an algorithm, which produces the entire solution set for the decision variables, but only for systems which may be described by an acyclic information structure. Christensen20 extended the algorithm by Lee et al.31 to handle persistent iteration, eliminating variables in equations with only one variable and moving as many nonlinearities as possible from the equality constraints to the objective function. The same principles are applied by Drud24 in a large-scale nonlinear solver (CONOPT) available in the GAMS environment. However, in Christensen’s algorithm20 a single solution is obtained for the decision variable set, and no information on the numerical properties of the solution procedure is taken into account. Edie and Westerberg21 have recognized the importance of including numerical information on the process of selecting decision variables. Their algorithm ranks decision variables based on avoiding leaving singular or near-singular square systems to solve. Their objective function is quantifiable (minimum eigenvalue principle), but it may produce solution paths that require solution of larger systems than obtained by a minimum tearing criterion.22 It also has the drawback of requiring some meaningful approximation of nonlinear systems, viz., by local linearization around some assumed solution point. In practice, the algorithm is applied in a series of steps, where the designer guesses a structurally valid decision variable set, obtains the solution path, and guesses another set of decision variables if the output subsystems are singular or nearly singular. Thus, this algorithm has the merit of avoiding solution paths, which may be either indeterminate or very ill-conditioned, but on a single pass produces a unique combination of decision and output/tearing variables.19 Stadtherr et al.22 propose an algorithm that extends Christensen’s work20 to allowable subsets of rank greater than 1. However, this algorithm tends to produce implicit nested loops, which are to be avoided by direct substitution.19 Numerical information may be included by prohibiting the solution of an equation for a particular variable. Book and Ramirez19 have proposed that structural as well as numerical (functional) information should be simultaneously considered for the choice of decision variables, where the occurrence matrix is transformed into a functionality matrix, by the identification of key nonlinear terms in the system’s equations. Their objective function is heuristic, based on four criteria which are to be met in sequence: assigning single-valued acyclic output elements (state variables) to single equations; assigning output elements to linear subsets; assigning multiple-valued acyclic output elements to

single equations; and assigning the minimum number of iterative variables to equations which contain those iterates. If the designer is happy with the solution procedure, the algorithm stops; otherwise, a new set of decision variables is chosen by the designer and the process repeated. A New Approach for the Choice of Decision Variables There are several characteristics that should be met by an algorithm for the choice of decision variables. It should (1) include structural information such that some minimization of tears and nested recycles may be achieved, (2) include functional information so that impossible or difficult outputs are avoided, (3) provide for some means of avoiding tears which are not adequate for direct substitution, if that is the solution method sought for any subset of the system’s equations, (4) provide flexibility in fixing a priori a subset for the decision variables, viz., discrete variables or variables that, for some reason, the designer may want to specify, and (5) show a different number of possible solution paths that obey some quantifiable objective function, so that the designer may choose among good alternatives. It is clear that none of the methods proposed so far and described above obeys all of these criteria. The Ledet and Himmelblau28 tearing/partition algorithm described before has some interesting features which make it a good candidate as a starting point for the solution of the proposed problem, because it allows prohibition of unwanted outputs and tears and because the computer code is available. However, this algorithm was devised based on square matrices, although rectangular matrices such as those found in optimization problems may also be used as inputs. The degrees of freedom responsible for the rectangularization are then simply forgotten by the algorithm, because they are assumed to have been fixed in some way (parametrized), and the method simply processes the equations. The general statement for an optimization problem is

minimize F(x,y)

(1)

subject to the constraints hk(x,y) ) 0, k ) 1, 2, ..., m1

(2)

gj(x,y) g 0, j ) 1, 2, ..., m2

(3)

Ri e xi e βi, i ) 1, 2, ..., p

(4)

γi e yi e ξi, i ) p + 1, ..., n

(5)

where x and y are respectively a vector of continuous variables and a vector of discrete variables. The number of degrees of freedom is given by the number of independent variables (n) minus the number of equality constraints (m1) and represents the number of variables that have to be simultaneously arbitrated (or fixed) to perform the simulation step. This general statement is referred to as a MINLP problem, which reduces to a NLP problem if the decision vector is only of continuous variables and to an integer (combinatorial) problem if the decision vector is only of discrete variables. In any case, the application of the tearing/partitioning algorithm of Ledet and Himmelblau28 to the above problem requires the identification of the occurrence matrix, where each variable is numbered from 1 to n and each


equation is numbered from 1 to m1, where n > m1. This matrix is simply an annotation of occurrences of variables within equations. The following steps are then performed: (1) The occurrence matrix is written in an appropriate form (pairs of equation/variable numbers), adding any information regarding variables that cannot be either output variables from the corresponding equations (also in pairs of equation/variable numbers) or tearing variables (in variable numbers). (2) The decomposition computer program (Appendix A in Ledet and Himmelblau28) is applied to the occurrence matrix to obtain an output set. This output set is an information flow matrix that represents the solution sequence of the entire system of equations and identifies which variables should be output from which equation and which irreducible subsets have to be simultaneously solved for. In this case, the output set also identifies the tearing variables, in case fixed-point iteration is employed for convergence of the recycle loops. For any fixed choice of the decision variables, a m1 × m1 matrix is obtained, which may be decomposed by the Ledet/Himmelblau algorithm. There are, however, several cases where decomposition or solution is not possible, i.e., whenever equations end up without variables or if subsets of equations that have fewer variables than equations are obtained. In this case, the solution is to try a new set for the decision variables. A systematization of the possible choice for these decision variables is, however, possible and desirable, because for all but the smallest problems the combinatorial search space quickly becomes intractable. The method devised in this work to implement this approach is the following: (i) Set up the occurrence matrix for the optimization problem, including any functional information, which may be input as prohibited output variable/equation pairs and tearing variables. (ii) Define an objective function, which is somehow related to the difficulty in solving a system of m1 equations. (iii) Interface the Ledet/Himmelblau algorithm with some combinatorial optimizer, which will choose the decision variables according to the minimization criteria defined in ii, constrained by whatever functional constraints and a priori choice of decision variables are appropriate. Point ii is obviously the most difficult to quantify, because there is no a priori criteria that may be universally applicable. To circumvent this shortcoming, it is hypothesized that it is generally easier to solve (a) a system of equations sequentially, (b) small systems rather than larger ones, and (c) disjoint systems rather than nested ones. While many examples may be given to disprove the above claims, they should be of general applicability, in the sense that they correspond probably to most situations encountered in practice. Because the objective function may be constrained by functional constraints, this further adds to the general applicability of the above criteria. An objective function that obeys the above criteria is as follows: SR

Case Studies

NR

ni + ∑nj} ∑ i)1 j)1

Fobj ) Minimize {

tions in simple recycle i, and nj the number of equations in nested recycle j. Although other objective functions are possible, eq 6 is simple and shows that a sequential solution has zero cost and that the cost of solving a system of equations increases with the number of simultaneous subsystems (recycles) and their degree of nesting. Obviously, this only makes sense if functional information is included as constraints to the optimization problem. The units of Fobj are in a number of equivalent equations, corresponding to a single unnested (disjoint) system. To determine these values, code lines were added to the Ledet/Himmelblau algorithm, as described below: (a) For each recycle variable identified by the Ledet/ Himmelblau algorithm, find the corresponding column and line numbers corresponding to the occurrence matrix and the number of equations belonging to the convergence procedure. These are counted from the line number down to the diagonal, because the diagonal element corresponding to the same column number is the output element that will end the recycle loop. (b) For each recycle variable, identify the number of nested recycles and their size, taking into account the column and line number information gathered above. Point iii may be fulfilled with many different algorithms, and in this work the MSGA adaptive random search algorithm15,32 was employed, because it is relatively fast and was found to be robust for various combinatorial problems. For more demanding problems, combinatorial optimizers such as those based on simulated annealing33,34 or faster variants,35 as well as algorithms based on mathematical programming,14,36 may be employed. The MSGA algorithm deals easily with constraints, such as forcing or avoiding certain variables from being decision variables, and penalizes those choices of decision variables that do not produce an output set from the Ledet/Himmelblau algorithm, by simply assigning a large value to the objective function. However, to avoid conflicts between functionality constraints at the simulation and optimization levels, the Ledet and Himmelblau code had to be significantly modified. To have a better knowledge of possible paths to the solution procedure of the simulation step, the optimizer should also rank the different outputs as a sequence of sets of decision variables that correspond to a decreased cost function. Thus, by interfacing the Ledet/Himmelblau algorithm with a combinatorial optimizer, the tearing/partitioning algorithm becomes the simulation step of an optimization procedure with the aim of identifying sets of decision variables that minimize the cost of solving the system of equations. Although Ledet and Himmelblau28 state that their algorithm is sought to be appropriate for systems on the order of 1000 nonlinear equations, the upper limit to the proposed algorithm should be lower than this, because of the significant work that would be necessary to set up the occurrence matrix and to rewrite the equations explicitly for the different output variables.

(6)

where SR (g1) is the number of simple recycles, NR (g2) the number of nested recycles, ni the number of equa-

The cases studied correspond to five problems of increased complexity, from the point of view of the size of the search space for the choice of the decision variables, in the ease in obtaining a sequential solution, and in the functional constraints involved. The first

4746

Ind. Eng. Chem. Res., Vol. 38, No. 12, 1999 Table 1. Occurrence Matrix and Variable Assignment for Case No. 1 a. Occurrence Matrix for Case No. 1 equation

variable

equation

variable

1 2 3 4

1, 2, 3, 4, 9, 10 1, 2, 3, 4, 5, 9 2, 3, 4, 9 1, 3, 4, 6, 9

5 6 7

1, 3, 6, 7, 9 1, 7, 11 1, 2, 8

prohibited output

prohibited output

equation

variable

equation

variable

1 2 3

1, 2 1, 2 2

4 5

1 1, 3, 6

b. Variable Assignment for Case No. 1 1 2 3 4

5 6 7 8

θ1 θ2 a h*

F* b* S* θ0

9 10 11

B0 D d*

Table 2. Occurrence Matrix and Variable Assignment for Case No. 2 a. Occurrence Matrix for Case No. 2 equation

variable

equation

variable

1 2 3 4 5 6

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

7 8 9 10 11 12

9, 10, 11 12, 13 2, 5 9, 10, 11, 14, 15, 16 14, 15 3, 7, 15

prohibited output

prohibited output

equation

variable

equation

variable

3 4

7 7

12

7

b. Variable Assignment for Case No. 2 1 2 3 4 5 6

Figure 1. (a) Lowest cost solution for example 1, constrained by the occurrence matrix in Table 1a. (b) Lowest cost solution for example 1, constrained by occurrence matrix in Table 1a, forcing {D, B0} as decision variables.

example is demonstrative while the remainders correspond to process design problems, of similar size to those that have been analyzed by various authors that have studied the problem of choosing decision variables.19-23 Example 1 This example was taken from Stadtherr et al.22 and represents a seven-equation system with 4 degrees of freedom, obtained in the calculation of forces between two floating cylinders. Thus, there are only 330 potential candidates for a solution strategy. However, not all outputs are allowed, so this is a problem constrained by functionality. Table 1a shows the occurrence matrix with the functional constraints, and Table 1b shows the corresponding variable assignment. Figure 1a shows the

k1 x1 y1 k2 x2 y2

7 8 9 10 11 12

T P L V F z1

13 14 15 16

z2 hL HV hf

output sets obtained by the proposed algorithm. As described by these authors, in conjunction with the functionality constraints, the designer might also want to fix {9, 10} ) {D, B0} as decision variables, in which case no sequential solution could be found (Figure 1b). The solutions were obtained in about 0.3 s CPU time and are similar to those reported. This example shows that the inclusion of functionality constraints for output variables or a priori choice of decision variables is easily done within the proposed framework. Example 2 This example is an adiabatic flash distillation of a benzene/toluene system, taken from Edie and Westerberg.21 Figure 2a shows a schematic diagram for the process. Table 2a shows the occurrence matrix and Table 2b the corresponding variable assignment. The pertinent equations are given in the Appendix. Because there are 12 equations and 16 independent variables, there are 4 degrees of freedom and potentially 1820 candidate solution sets for the decision variables. The solution scheme proposed by Edie and Westerberg21 has {3, 8, 10, 11} ) {y1, P, V, F} as decision variables and

Ind. Eng. Chem. Res., Vol. 38, No. 12, 1999 4747 Table 3. Occurrence Matrix and Variable Assignment for Case No. 3 a. Occurrence Matrix for Case No. 3 equation

variable

equation

variable

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23

1 2 3, 4 4 6, 7 8, 9 10, 49, 50, 51 11, 12, 13 12, 14 13, 15 15, 16, 52 17, 49, 52 4, 13, 18 5, 14, 18 7, 15, 18 18, 50, 52 8, 16, 19 19, 51, 52 9, 16, 20 10, 20, 52 11, 21, 53 12, 22, 23 13, 24, 25

24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45

14, 26, 54 15, 27, 28 16, 29, 55 30, 31, 52 21, 22, 24 22, 26, 32 24, 27, 32 27, 29, 30 13, 33, 34 32, 35, 36 23, 25, 37 28, 38, 54 37, 38, 39 39, 40 36, 37, 38, 40 41, 42, 52 43, 44 33, 44 35, 45 46, 47 41, 47 43, 45, 46, 48

b. Variable Assignment for Case No. 3 1 2 3 4 5 6 7 8 9 10

x11 x22 x13 x23 x24 x15 x25 x36 x46 x47

11 12 13 14 15 16 17 18 19 20

F1 F2 F3 F4 F5 F6 y1 y2 y3 y4

21 22 23 24 25 26 27 28 29 30

z1 z2 T2 z3 T3 z4 z5 T5 z6 z7

31 32 33 34 35 36 37 38 39 40

T7 Q2 V1 t1 A2 ∆Tlm2 W1 W2 W3 W4

41 42 43 44 45 46 47 48 49 50

V3 t3 S1 v1 S2 S3 v3 S x17 x27

51 52 53 54 55

x37 F7 T1 T4 T6

requiring only the solution of two simultaneous equations. However, it is much better to introduce functional constraints to avoid the temperature calculation from either eq 3 or eq 4 (involving logarithms) or eq 12 (involving multiple roots). These constraints are shown in Table 2a, which dictate that T must be a decision variable, because it only occurs in these three equations. With the added constraints, the solution cost rises to four simultaneous (unnested) equations, and several solutions were obtained:

{3, 4, 7, 10} ) {y1, k2, T, V} {1, 2, 7, 10} ) {k1, x1, T, V} {2, 7, 11, 14} ) {x1, T, F, hL} Figure 2. (a) Flash distillation system (Edie and Westerberg21). (b) Original solution for example 2 (Edie and Westerberg21). (c) Lowest cost solution for example 2, constrained by the occurrence matrix in Table 2a.

results in a solvable system of eight simultaneous equations (Figure 2b). The application of the algorithm proposed in this paper (unconstrained) results in several low cost solutions, obtained in less than 1 min CPU time, such as

{1, 2, 10, 11} ) {k1, x1, V, F} {1, 2, 9, 10} ) {k1, x1, L, V} {2, 6, 9, 10} ) {x1, y2, L, V} {6, 10, 11, 12} ) {y2, V, F, z1}

{5, 7, 11, 14} ) {x2, T, F, hL} {1, 7, 11, 14} ) {k1, T, F, hL} {1, 7, 9, 15} ) {k1, T, L, HV} {5, 7, 10, 15} ) {x2, T, HV, V} {7, 8, 10, 14} ) {T, P, V, hL} {5, 6, 7, 11} ) {x2, y2, T, F} {1, 6, 7, 11} ) {k1, y2, T, F} {1, 7, 11, 15} ) {k1, T, F, HV} All of these solutions have in common that the same 4 × 4 linear system has to be solved, namely, eqs {5, 6,

4748


Figure 3. (a) Mixer-exchanger-mixer system (Book and Ramirez19). (b) Lowest cost solution for example 3.


{22, 46, 48} ) {z1, S3, S}

Table 4. Occurrence Matrix and Variable Assignment for Case No. 4

{25, 34, 47} ) {T3, t1, v3}

a. Occurrence Matrix for Case No. 4 equation

variable

equation

variable

{22, 43, 47} ) {z2, S1, v3}

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

2, 3 1, 2, 4, 5, 6 3, 6, 7, 44, 45 3, 8, 9, 10, 45 1, 2, 9, 12 1, 2, 10, 12 7, 8, 11 12, 13 3, 13, 14 3, 14, 45 1, 2, 15, 42 1, 11, 43 1, 16, 17 16, 18 16, 19, 20 17, 20 1, 2, 21, 22 1, 17, 21, 23 21, 24, 25 17, 22, 23, 25

21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40

1, 2, 4, 5 1, 2, 22, 23, 28, 29, 30 29, 31 4, 34 5, 35 15, 35, 36 33, 34, 36 2, 28, 33 28, 37 28, 38, 39 30, 39 29, 32 29, 40, 41 23, 41 5, 23 30 22, 26, 27 4, 27 1, 2, 5, 44 1, 2, 45

{25, 33, 46} ) {T3, V1, S3}

b. Variable Assignment for Case No. 4 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

L D G x1 x2 y2 Nog Hog Hg Hl Z AA DA Gf N

16 17 18 19 20 21 22 23 24 25 26 27 28 29 30

Q4 T5 W4 A4 ∆Tlm4 Q3 Tf T4 A3 ∆Tlm3 q T2 Q1 Q2 T3

31 32 33 34 35 36 37 38 39 40 41 42 43 44 45

AS W2 R Rm Nm X W1 A1 ∆Tlm1 A2 ∆Tlm2 HP1 HP2 x*2 L*

7, 8}. Figure 2c shows one such solution path, and this system, which is not singular or ill conditioned at the solution point of Edie and Westerberg,21 namely, {x2, y2, T, F} ) {0.7, 0.5, 373.2 K, 1 mol/h}, is easily solved by Gauss or LU decomposition. However, if this was not the case, then a higher cost solution would have to be sought. Because the proposed algorithm ranks possible solution sets by a decreased objective function, this task should be relatively easy to perform. Example 3 The following example, taken from Book and Ramirez,19 represents an optimization of a mixer-exchanger-mixer system, with 45 equations (given in the Appendix) in 55 variables. Figure 3a is a schematic diagram of the process. Table 3a shows the occurrence matrix and Table 3b the corresponding variable assignment. These authors have first assumed that 7 out of the 10 degrees of freedom were known, viz., {49, 50, 51, 52, 53, 54, 55} ) {x17, x27, x37, F7, T1, T4, T6}, such that only 3 degrees of freedom remained, with 17 296 possible combinations for the decision vector. They have also assumed an unconstrained problem for the choice of decision variables. The proposed algorithm was applied to this problem, and a minimum cost of three simultaneous equations was obtained in about 4 s CPU time, albeit with many different combinations of decision variables, namely

{25, 46, 48} ) {T3, S3, S} {22, 44, 46} ) {z2, v1, S3} {25, 47, 48} ) {T3, v3, S} Book and Ramirez19 using their algorithm have reached a similar conclusion. These authors have relaxed knowledge of the seven variables described above. A system with 10 degrees of freedom and about 2.9 × 1010 potential solution sets occurs, and a sequential solution was obtained with {8, 10, 23, 34, 42, 49, 52, 53, 54, 55} ) {x36, x47, T2, t1, t3, x17, F7, T1, T4, T6} as decision variables. The application of the proposed algorithm to this problem produced, in about 4 s CPU time, besides the solution reported by these authors, the following alternative sequential solutions:

{4, 18, 20, 24, 43, 48, 52, 53, 54, 55} ) {x23, y2, y4, z3, S1, S, F7, T1, T4, T6} {9, 10, 18, 21, 23, 28, 31, 43, 47, 50} ) {x46, x47, y2, z1, T2, T5, T7, S1, v3, x27} {10, 18, 22, 29, 44, 46, 52, 53, 54, 55} ) {x47, y2, z2, z6, v1, S3, F7, T1, T4, T6} {10, 21, 22, 25, 30, 33, 42, 46, 53, 54} ) {x47, z1, z2, T3, z7, V1, t3, S3, T1, T4} {10, 21, 22, 25, 30, 33, 42, 46, 54, 55} ) {x47, z1, z2, T3, z7, V1, t3, S3, T4, T6} {3, 8, 16, 21, 22, 25, 30, 43, 48, 55} ) {x13, x36, F6, z1, z2, T3, z7, S1, S, T6} {9, 19, 23, 24, 26, 38, 43, 47, 54, 55} ) {x46, y3, T2, z3, z4, W2, S1, v3, T4, T6} {4, 10, 32, 43, 46, 51, 52, 53, 54, 55} ) {x23, x47, Q2, S1, S3, x37, F7, T1, T4, T6} {20, 21, 26, 29, 32, 37, 42, 44, 47, 55} ) {y4, z1, z4, z6, Q2, W1, t3, v1, v3, T6} One such solution can be seen in Figure 3b. Example 4 This example was adapted from Umeda37 and Umeda and Ichikawa38 and represents an absorber-stripping system with heat integration, which is used for solvent recovery. Figure 4a shows a schematic diagram for the process. The pertinent equations are provided in the Appendix, but for a thorough explanation of the equations, the reader is referred to the original papers. The variables specified in the Appendix, together with the system of 40 equations, give 45 unknown variables and 5 degrees of freedom. These correspond to a potential of 1.2 × 106 possible combinations for sets of the decision variables. Table 4a gives the occurrence

4750


Figure 4. (a) Absorber-stripper for solvent recovery (Umeda37 and Umeda and Ichikawa38). (b) Lowest cost solution for the absorberstripper.


matrix and Table 4b the corresponding variable assignment. The algorithm proposed in this work was applied to this example, and several sequential solutions could be obtained in about 3 s CPU time. Umeda37 and Umeda and Ichikawa38 provide a set of decision variables that produces a sequential solution, namely, {1, 2, 4, 26, 33} ) {L, D, x1, q, R}. Without imposing any constraints on the choice of decision variables, the following sets which also produce a sequential solution were obtained by the proposed algorithm:

{2, 13, 20, 23, 36} ) {D, DA, ∆Tlm4, T4, X} {3, 14, 26, 34, 40} ) {G, Gf, q, Rm, A2} {3, 12, 16, 21, 36} ) {G, AA, Q4, Q3, X} {18, 19, 27, 31, 45} ) {W4, A4, T2, As, L*} {3, 26, 34, 40, 45} ) {G, q, Rm, A2, L*} {2, 20, 38, 44, 45} ) {D, ∆Tlm4, A1, x2*, L*} Obviously, some of these sets are better than others, because some variables may be bounded and easier to arbitrate. The inclusion of constraints may, however, produce still better sequential solution paths. In this case, because T2 and T4 have to be computed by some iterative algorithm, it would be better to force them to be decision variables, to linearize eqs 35 and 38. The same does not apply for T3, because it must be an output variable, as it is the only unknown variable in eq 36. Forcing {23, 27} ) {T4, T2} to be decision variables produced the following sets which give a sequential solution:

{3, 26, 36} ) {G, q, X} {15, 18, 38} ) {N, W4, A1} {2, 26, 36} ) {D, q, X} {22, 28, 36} ) {Tf, Q1, X} {19, 20, 36} ) {A4, ∆Tlm4, X} Further, it may be of interest to pose this problem as a MINLP problem, if the number of theoretical stages N is specified as an integer. In this case, forcing {15, 23, 27} ) {N, T4, T2} to be decision variables produced the following sets which give a sequential solution:

{18, 38} ) {W4, A1} {20, 28} ) {∆Tlm4, Q1} {1, 16} ) {L, Q4} where the first solution had already been found. Figure 4b shows one of the sequential solutions obtained, which corresponds to a straightforward simulation. Example 5 The last example was taken from Zaldivar39 and refers to the nitration of toluene in a two-liquid-phase CSTR reactor, according to the following reaction:

H2SO4C7H8 + HNO3 f C7H7NO2 + H2O

The two liquid phases are organic and aqueous, with the reaction proceeding only in the aqueous phase. The organic phase is homogeneous, with a constant bulk concentration. The aqueous phase has a spatial distribution in a liquid film and a constant bulk concentration. The system is isothermal, and all resistance to diffusion is in the liquid film, where the hypotheses of a film model with a stagnant layer apply.40 Furthermore, the product concentration in the aqueous phase is negligible; viz., the reaction in the aqueous phase is followed by an instantaneous diffusion to the organic phase. The mixer geometry and rotation speed are fixed. Figure 5a shows a schematic diagram of the process, and the pertinent equations are given in the Appendix. Table 5a shows the occurrence matrix and Table 5b the corresponding variable assignment. This problem has 45 unknown variables and 36 independent equations, viz., 9 degrees of freedom. There are thus potentially about 9 × 108 candidate sets for assignment to the degrees of freedom. This example was translated into the GAMS16 environment, but convergence could only be achieved if a starting solution close to the optimized solution was provided.41 The Ledet/Himmelblau algorithm was applied to this example, showing that the equation ordering (as entered in the GAMS environment) corresponded to three nested recycles (Figure 5b), with a cost (given by eq 6) corresponding to the solution of 8568 (unnested) equations. The algorithm proposed in this paper was applied (unconstrained) to this example, and several sequential solutions were obtained in about 50 s CPU time. These are

{1, 4, 6, 8, 9, 22, 23, 40, 45} ) {norg, xporg, xpaq, xsaq, xWaq, k2, TR, xAin,org, VR} {4, 6, 7, 8, 10, 15, 29, 36, 37} ) {xporg, xpaq, xNaq, xsaq, Jout,org, v*,org, klA, δ, AP} {3, 6, 8, 9, 10, 24, 26, 34, 37} ) {xAorg, xpaq, xsaq, xWaq, Jout,org, DA, We, Vaq,film, AP} {1, 4, 6, 7, 9, 22, 29, 40, 45} ) {norg, xporg, xpaq, xNaq, xWaq, k2, klA, DA, xAin,org, VR} {4, 5, 7, 8, 9, 10, 17, 36, 45} ) {xPorg, xAaq, xNaq, xSaq, xWaq, Jout,org, CAorg, δ, VR} {5, 9, 10, 15, 17, 22, 24, 36, 45} ) {xAaq, xWaq, Jout,org, V*,org, CAorg, k2, DA, δ, VR} Although it took somewhat long to find sequential solutions, very low-cost solutions, corresponding to solving 2 × 2 systems could be obtained in about 4 s CPU time. The sequential sets above are not good, because they involve the simultaneous independent choice of several compositions belonging to the same stream, and this requires that the stoichiometry be checked and normalized at each simulation. Thus, constraints were added at each stream (aqueous and organic) that there should only be one single composition chosen as the decision variable. With these constraints, low-cost solutions corresponding to solving simultaneously two equations could be obtained, such as

4752



Figure 5. (a) Nitration of toluene (Zaldivar39). (b) Original high-cost solution matrix for the nitration of toluene. (c) Lowest cost solution for the nitration of toluene, constrained by specifying one single composition for each stream as the decision variable.

{3, 7, 19, 21, 23, 24, 30, 40, 45} ) {xAorg, xNaq, CNaq, Vaq, TR, DA, φ, xAin,org, VR} {4, 7, 10, 18, 19, 22, 23, 29, 45} ) {xPorg, xNaq, Jout,org, CAaq, CNaq, k2, TR, klA, VR} The last solution is given in Figure 5c, and the corresponding simulation could be performed without convergence problems, using either the GAMS environment or direct Fortran coding, starting from different points.41 Conclusions For general-purpose process system optimization, SQP-based algorithms are very efficient and do not require that a simulation be performed at each optimization step. However, these are local search algorithms that employ constraint linearizations that may cut off the global optimum.7,8 Thus, they are not adequate for nonconvex global NLP or MINLP optimization. Stochastic algorithms15,32 or mathematical programming algorithms7,8,14,36 are better suited for this purpose, but the prior step of simulation is now required. This may be performed using either the modular or the equation-

oriented approach. While the modular approach has been mostly adopted in industrial environments, the equation-oriented approach has the benefit of a detailed knowledge of the model equations. The added difficulty of having to write the entire set of equations is offset by the greater flexibility in solving nested recycle loops and in the incorporation of new unit specifications and design constraints.2,5 Some modeling languages have been developed, such as GAMS16 or AMPL,17 that take out much of the burden associated with this process. However, unlike with modular simulators, a good initialization scheme is usually a prerequisite to a successful full-scale solution of the entire system of equations.3 For the case of optimization problems requiring the prior step of simulation, the decision variables must be somehow fixed or arbitrated, so that a determinate system is eventually obtained. The structure of the system equations, viz., their ordering or grouping, may however have a strong impact on the solution process. Because this structure may change drastically by specifying different independent variables for the degrees of freedom, this feature may be exploited in order to efficiently solve sparse systems of nonlinear equations. This has been recognized before, and efficient algo-

4754


Table 5. Occurrence Matrix and Variable Assignment for Case No. 5 a. Occurrence Matrix for Case No. 5 equation

variable

equation

variable

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

10, 12, 13, 38 11, 12, 13, 39 3, 12, 13, 38, 40 3, 4 5, 12, 13, 14, 39, 41 5, 6, 7, 8, 9 7, 12, 13, 14, 39, 42 8, 12, 13, 39, 43 9, 12, 13, 14, 39, 44 3, 4, 15, 23 5, 6, 7, 8, 9, 16, 23 3, 15, 17 5, 16, 18 7, 16, 19 1, 15, 20 2, 16, 21 7, 9, 22, 23 1, 2, 10, 38, 39

19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36

7, 8, 9, 23, 24 21, 25, 26, 45 3, 4, 27 5, 6, 7, 8, 9, 28 15, 26, 27 24, 25, 29 19, 22, 24, 29, 30 17, 31 12, 18, 29, 30, 31, 37 18, 29, 30, 31, 32, 37 12, 32, 33 14, 18, 19, 22, 33, 35 34, 36, 37 21, 34, 35 21, 25, 37 24, 29, 36 13, 14 20, 21, 45

b. Variable Assignment for Case No. 5 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

norg naq xAorg xPorg xAaq xPaq xNaq xSaq xWaq Jout,org Jout,aq JAdif,φ JPdif R v*,org

16 17 18 19 20 21 22 23 24 25 26 27 28 29 30

v*,aq CAorg CAaq CNaq Vorg Vaq k2 TR DA d32 We Morg Maq klA φ

31 32 33 34 35 36 37 38 39 40 41 42 43 44 45

CAaq,φ JAdif,δ ∆JAdif Vaq,film Vaq,bulk δ AP Jin,org Jin,aq xAin,org xAin,aq xNin,aq xSin,aq xWin,aq VR

rithms have been devised for systems that may be expressed as acyclic.23,31 For the general case where an acyclic structure might not be obtainable, proposed algorithms19-22 rely on the designer to choose alternative sets of decision variables. In this paper, the simulation step of an optimization problem is posed as an optimization problem itself, whereby a cost function for solving the entire system of equations is associated with the choice of the decision variables. By coupling the tearing/partitioning algorithm of Ledet and Himmelblau28 with a MINLP optimizer,32 a combinatorial problem is formulated whereby minimum cost output sets are obtained. These optimal solutions correspond to sets of decision variables that produce either a sequential solution for the entire system of equations or simulations with a minimum number of nested recycles and simultaneous equations, while obeying all imposed functional constraints and a priori choice of decision variables. The proposed approach was applied to several optimization problems, and the results obtained reveal that it may be a powerful preprocessing step for nonconvex global optimization. The results given in this paper for five selected examples show that functional constraints may be easily added at the simulation level and decision variable constraints at the optimizer level to fine-tune the solution procedure and to take into account problemspecific characteristics. Acknowledgment This work was partially supported by FCT (Portuguese Foundation for Science and Technology), under

Contract PRAXIS 3/3.1/CEG/2641/95, employing the computational facilities of Instituto de Sistemas e Robo´tica (ISR), Porto. Nomenclature Fobj ) objective function gj ) inequality constraint, j ) 1 to m2 hk ) equality constraint, k ) 1 to m1 m1 ) number of equality constraints m2 ) number of inequality constraints n ) number of independent parameters (continuous and discrete) p ) number of independent continuous parameters x ) continuous decision vector y ) discrete decision vector Greek Symbols Ri ) lower bound on continuous parameter xi (i ) 1, p) βi ) upper bound on continuous parameter xi (i ) 1, p) γi ) lower bound on discrete parameter yi (i ) 1, n - p) ξi ) upper bound on discrete parameter yi (i ) 1, n - p) Abbreviations MINLP ) mixed-integer nonlinear programming MSGA ) Minlp SGA (Salcedo-Gonc¸ alves-Azevedo) algorithm NLP ) nonlinear programming SQP ) sequential quadratic programming

Appendix Example 2.

f1 ) k1x1 - y1 ) 0 f2 ) k2x2 - y2 ) 0 A + 2.30B )]/P - k ) 0 [ (-0.120 T -0.120 f ) [exp( A + 2.30B )]/P - k ) 0 T

f3 ) exp 4

1

1

1

2

2

2

f5 ) x1L + y1V - z1F ) 0 f6 ) x2L + y2V - z2F ) 0 f7 ) L + V - F ) 0 f8 ) z1 + z2 - 1 ) 0 f9 ) x1 + x2 - 1 ) 0 f10 ) hLL + HVV - hfF ) 0 f11 ) HV - hL - K ) 0 f12 ) HV - H0V1(y1) - H0V2(1 - y1) [R1y1 + R2(1 - y1)]T + [β1y1 + β2(1 - y1)]T2 ) 0 where xi, yi, and zi are molar fractions (i ) 1, 2), F, V, and L are molar flow rates, HV, hL, and hf are enthalpies, T is temperature, P is pressure, Ki is an equilibrium constant (i ) 1, 2), and H0V1, H0V2, R1, R2, β1, β2, and K are known constants.


Example 3.

f34 ) W1 + T3 - T2 ) 0

Stoichiometry:

f35 ) W2 + T5 - T4 ) 0 f1 ) x11 - 1 ) 0

f36 ) W2W3 - W1 ) 0

f2 ) x22 - 1 ) 0

f37 ) W4 - ln(W3) ) 0

f3 ) x13 + x23 - 1 ) 0

f38 ) W4∆Tlm2 + W2 - W1 ) 0

f4 ) x24 - 1 ) 0

f39 ) V3 - F7t3/F7 ) 0

f5 ) x15 + x25 - 1 ) 0

Capital cost estimation:

f6 ) x36 + x46 - 1 ) 0

f40 ) ln(S1) - 0.515 ln(v1) - 3.354 ) 0

f7 ) x17 + x27 + x37 + x47 - 1 ) 0

f41 ) v1 - V1 ) 0

Flow balances: f8 ) F1 + F2 - F3 ) 0 f9 ) F2 - F4 ) 0

f42 ) ln(S2) - 0.699 ln(A2) - 2.414 ) 0 f43 ) ln(S3) - 0.515 + ln(v3) - 3.354 ) 0 f44 ) v3 - V3 ) 0

f10 ) F3 - F5 ) 0

f45 ) S - S1 - S2 - S3 ) 0

f11 ) F5 + F6 - F7 ) 0 Component balances: f12 ) x17F7 - y1 ) 0 f13 ) x23F3 - y2 ) 0 f14 ) x24F4 - y2 ) 0 f15 ) x25F5 - y2 ) 0 f16 ) x27F7 - y2 ) 0 f17 ) x36F6 - y3 ) 0 f18 ) x37F7 - y3 ) 0 f19 ) x46F6 - y4 ) 0

where v1, V1, v3, and V3 are capacities (m3) and A2 is heat exchanger area (m2). Example 4.

Absorber: f1 ) F - G - Dx3 ) 0 f2 ) Fy1 + Lx2 - (L + D)x1 - Gy2 ) 0 f3 ) Nog -

[ ( ) ]

H y - x* 1 HG 1 P 2 HG + ln 1 )0 HG H PL* PL* 1y2 - x*2 PL* P HGHl )0 f4 ) Hog - Hg + PL*

(

)

f20 ) x47F7 - y4 ) 0 Energy balances:

f5 ) Hg -

f21 ) F1T1 - z1 ) 0 f22 ) F2T2 - z2 ) 0 f23 ) F3T3 - z3 ) 0

f6 ) Hl -

f24 ) F4T4 - z4 ) 0

[

f32 ) V1 - F3t1/F3 ) 0 f33 ) Q2 - U2A2∆Tlm2 ) 0

)0

f8 ) AA - π

f28 ) C1z1 + C2z2 - C3z3 ) 0

Equipment specification:

]

γ

(L + D)Ml δ φ(Scl)0.5 ) 0 AAµl

f27 ) F7T7 - z7 ) 0

f31 ) C5z5 + C6z6 - C7z7 ) 0

]

(L + D)Ml AA

DA2 )0 4

f26 ) F6T6 - z6 ) 0

f30 ) C3z3 + Q2 - C5z5 ) 0

[

FMg β 0.5 Scg AA

f7 ) Z - NogHog ) 0

f25 ) F5T5 - z5 ) 0

f29 ) C3z3 + Q2 - C2z2 ) 0

( )

R

f9 ) DA -

[

x

4GMg

0.75πGf(3600)

]

1 (µ /µ )0.2 f10 ) log Gf2(ap/3) gFlFg l g L*Ml 1.74 GMg

( )

0.25

)0

(Fl/Fg)-0.125 ) 0

Pump no.1 f11 ) HP1 - KP1(L + D)N ) 0 Pump no.2 f12 ) HP2 - KP2LZ ) 0

4756


Heat exchanger no.4 f13 ) Q4 - Lcp(T5 - T1) ) 0 f14 ) Q4 - W4cpw∆Tcw ) 0

exp A′1 -

f15 ) Q4 - U4A4∆Tlm4 ) 0 f16 ) ∆Tlm4 -

(T5 - Tcw0) - (T1 - Tcwi) )0 T5 - Tcw0 ln T1 - Tcwi

(

Px3

(

f36 )

)

(

)

B′1 T2 + Tref

P

f17 ) Q3 - (L + D)cp(Tf - T1) ) 0

(T1 - T5) - (Tf - T4) )0 f20 ) ∆Tlm3 T1 - T5 ln Tf - T 4

)

Stripper and heat exchangers:

f22 ) (L + D)cpTf + Q2 - DcpT3 + LcpT4 + Q1 ) 0

f24 ) Rm -

λMgKxFg(Fl - Fg)

[ [

)0

]

1 x3 a(1 - x3) )0 a - 1 x1 1 - x1

f25 ) Nm -

f26 ) N -

]

x3(1 - x2)

ln

x2(1 - x3) ln(a)

)0

Nm + X )0 1-X

[ (

f27 ) X - 0.75 1 -

) ]

R - Rm R+1

0.5668

)

+

(

)

B′2 T2 + Tref

Makeup f39 ) x*2 -

Lx2

)0

L + D(1 - x3)

f40 ) L* - D(1 - x3) - L ) 0 Known operating conditions:

Known physical properties: phase densities molar weights specific heat (streams w1 and w4) latent heats (heat exchangers 3 and 4) latent heat (from stripper) latent heat (stream w2) phase viscosities solute diffusivity in both phases antoine parameters parameters specific area of fillings relative volatility pump parameters Known global heat-transfer coefficients: Known Henry’s constant: Known reference temperature:

F l, F g Ml, Mg cpw cp λ λw µl, µg Dl, Dg A1′, A2′, B1′, B2′ R, β, γ, δ, φ, K ap/3 a KP1, KP2 Ui (i ) 1, 4) H Tref

Example 5.

)0

f1 ) Jin,org - Jout,org - JAdif,φ + JPdif ) 0 f2 ) Jin,aq - Jout,aq + JAdif,φ - JPdif ) 0

f28 ) Q1 - Dλ(R + 1) ) 0 f29 ) Q1 - W1cpw∆Tcw ) 0

f3 ) Jin,org(xAin,org - xAorg) - xAorg(JPdif - JAdif,φ) JAdif,φ ) 0

f30 ) Q1 - U1A1∆Tlm1 ) 0 (T3 - Tcw0) - (T3 - Tcwi) )0 f31 ) ∆Tlm1 T3 - Tcw0 ln T3 - Tcwi f32 ) Q2 - W2λw ) 0

(

f34 ) ∆Tlm2 - (Ts - T4) ) 0

f35 )

P

B′1 T4 + Tref

)

f4 ) xAorg + xPorg - 1 ) 0 f5 ) Jin,aq(xAin,aq - xAaq) - xAaq(JAdif,φ - JPdif) + JAdif,φ + R ) 0 f6 )

f33 ) Q2 - U2A2∆Tlm2 ) 0

x2 exp A′1 -

∑i xiaq - 1 ) 0;

i ) (A, P, N, S, W)

f7 ) Jin,aq(xNin,aq - xNaq) - xNaq(JAdif,φ - JPdif) in,aq

f8 ) J

+

(

-1)0

F, P, y1, x3, T1, Ti ) Tcwi, T0 ) Tcwo, ∆Tcw, Ts

f21 ) (L + D)x1 - Lx2 - Dx3 ) 0

Q2VM

)

P

f19 ) Q3 - U3A3∆Tlm3 ) 0

f23 ) AS -

(

(1 - x1) exp A′2 -

f18 ) Q3 - Lcp(T4 - T5) ) 0

(

P(1 - x3) B′2 exp A′2 T3 + Tref 1)0

f37 ) Tf - T2 + (q - 1)λ/cp

x1 exp A′1 f38 )

Heat exchanger no.3

B′1 T3 + Tref

+

)

B′2 (1 - x2) exp A′2 T4 + Tref -1)0 P

(xS

in,aq

aq

aq

- xS ) - xS (JA

dif,φ

R)0

dif

- JP ) ) 0

f9 ) Jin,aq(xWin,aq - xWaq) - xWaq(JAdif,φ - JPdif) + R ) 0 f10 ) v*,org - xA*xAorg + vP*xPorg ) 0


f11 ) v*,aq -

∑i vi*xiaq ) 0;

i ) (A, P, N, S, W)

f12 ) xAorg - CAorgv*,org ) 0 f13 ) xAaq - CAaqv*,aq ) 0 f14 ) xNaq - CNaqv*,aq ) 0 f15 ) Vorg - norgv*,org ) 0 f16 ) Vaq - naqv*,aq ) 0 f17 ) ln(k2) + 5.01 - 18.64xWaq + 30.54xNaq + 218.92 47.14xNaqxWaq + (-40.76 + 73.00xWaq TR 32.48xNaq + 84.03xNaqxWaq) ) 0 f18 ) Jout,org(norg + naq) - (Jin,org + Jin,aq)norg ) 0 f19 ) DA - 7.4107 × 10 TRx1.05xN MN + 2.6xW MW + 2.0xS MS ) 0 -9

aq

aq

(

f20 ) d32 - 4.05 × 10-1 1 + 4.47

aq

)

VW |We|-0.6 ) 0 VR

f21 ) Morg - xPorg MP - xAorgMA ) 0

∑i xiaqMi ) 0;

f22 ) Maq -

i ) (A, P, N, S, W)

Morg f23 ) We - 1.7408 × 10-1 ,org ) 0 v* DA f24 ) klA - 179 )0 d32 f25

k2CNaqDA x ) φ - 10 )0 klA

f26 ) CAaq,φ - 3 × 10-4CAorg ) 0 dif,φ

f27 ) JA

( (

) )

CAaq φ aq,φ - 3.6APklA C )0 tanh(φ) A cosh(φ)

CAaq,φ φ f28 ) JAdif,δ - 3.6APklA - CAaq ) 0 tanh(φ) cosh(φ) f29 ) ∆JAdif - (JAdif,φ - JAdif,δ) ) 0 f30 ) R - ∆JAdif - 3.6k2CNaqCAaqVaq,bulk ) 0 f31 ) Vaq,film - APδ ) 0 f32 ) Vaq,bulk - (Vaq - Vaq,film) ) 0 f33 ) AP -

0.6Vaq )0 d32

DA )0 f34 ) δ - 100 klA f35 ) JPdif - R ) 0 f36 ) VR - (Vaq + Vorg) ) 0

where vi* ) c1(1.0-TR/c2)c3/C4 and c1, c2, c3, and c4 are known constants, AP is the phase contact area (m2), Cij(φ,δ) is the concentration (mol/L), d32 is the Sauter mean diameter (cm), DA is molecular diffusivity (cm2/ s), i is component (A, P, N, S, W), j is phases org or aq, Jin,j, Jidif(film,bulk), Jout,j, and ∆Jidif are molar flow rates (kmol/h), k2 is a reaction rate constant (L/mol‚s), klA is the mass-transfer coefficient (mm/s), Mj is the mean molecular weight (g/mol), ni is the number of moles, R is production by reaction (kmol/h), TR is the reactor temperature (K), v*,j is the molar volume of the mixture, phase j (L/mol), vi* is the molar volume of pure component i (L/mol), Vj(film,bulk) is the volume, phase j (L), VR is the reactor volume (L), We is the Weber number (related to the mixing pattern), xiin,j, and xIj are molar fractions, δ is the film thickness (mm), φ is a nondimensional number (related to the reaction/diffusion competition). Literature Cited (1) Himmelblau, D. M. Basic Principles and Calculations in Chemical Engineering, 5th ed.; Prentice-Hall: Englewood Cliffs, NJ, 1989. (2) Seider, W. D.; Brengel, D. D.; Widagdo, S. Nonlinear Analysis in Process Design. AIChE J. 1991, 37 (1), 1. (3) Biegler, L. T.; Grossmann, I. E.; Westerberg, A. W. Systematic Methods of Chemical Process Design; Prentice-Hall: Englewood Cliffs, NJ, 1997. (4) Edgar, T. M.; Himmelblau, D. M. Optimization of Chemical Processes McGraw-Hill: New York, 1989. (5) Biegler, L. T. Chemical Process Simulation. Chem. Eng. Prog. 1989, Oct, 50. (6) Biegler, L. T. On the Simultaneous Solution and Optimization of Large Scale Engineering Systems. Comput. Chem. Eng. 1988, 12 (5), 357. (7) Smith, E. M. B.; Pantelides, C. C. Global Optimization of Nonconvex MINLPs. Comput. Chem. Eng. 1997, 21 (Suppl), S791. (8) Kocis, G. R.; Grossmann, I. E. Global optimization of Nonconvex Mixed-Integer Nonlinear Programming (MINLP) Problems in Process Synthesis. Ind. Eng. Chem. Res. 1988, 27 (8), 1407. (9) Diwekar, U. M.; Grossmann, I. E.; Rubin, E. S. An MINLP Process Synthesizer for a Sequential Modular Simulator. Ind. Eng. Chem. Res. 1992, 31 (1), 313. (10) Diwekar, U. M.; Rubin, E. S. Efficient Handling of the Implicit Constraints Problem for the ASPEN MINLP Synthesizer. Ind. Eng. Chem. Res. 1993, 32 (9), 2006. (11) Montagna, J. M. Problem-oriented modules for Process SimulationsResolution Strategies for Simulation and Optimization. Can. J. Chem. Eng. 1993, 71, 634. (12) Press, W. H.; Flannery, B. P.; Teukolsky, S. A.; Vetterling, W. T. Numerical Recipes in Fortran 77: The art of scientific computing; Cambridge University Press: New York, 1992; Vol. 1. (13) Motard, R. L.; Shacham, M.; Rosen, E. M. Steady-State Chemical Process Simulation. AIChE J. 1975, 21 (3), 417. (14) Grossmann, I. E.; Daichendt, M. M. New trends in optimization-based approaches to process synthesis. In Proceedings of PSE′94, Fifth International Symposium on Process Systems Engineering; En Sup Yoon, Ed.; Korean Institute of Chemical Engineers: Kyongju, Korea, 1994; p 95. (15) Cardoso, M. F.; Salcedo, R. L.; Feyo de Azevedo, S.; Barbosa, D. A Simulated Annealing Approach to the Solution of MINLP Problems. Comput. Chem. Eng. 1997, 21 (12), 1349. (16) Brooke, A.; Kendrick, D.; Meeraus, A. GAMS: A User’s Guide, Release 2.25; The Scientific Press Series, Boyd & Fraser Pub. Co.: Redwood City, CA 1992. (17) Fourer, R.; Gay, D. M.; Kernighan, B. W. AMPL: A Modeling Language for Mathematical Programming; The Scientific Press Series, Boyd & Fraser Pub. Co.: Redwood City, CA, 1993. (18) Rudd, D. F.; Watson, C. C. Strategy of Process Engineering; Wiley International: New York, 1968. (19) Book, N. L.; Ramirez, W. F. Structural Analysis and Solutions of Systems of Algebraic Design Equations. AIChE J. 1984, 30 (4), 609.

4758


(20) Christensen, J. H. The Structuring of Process Optimization. AIChE J. 1970, 16 (2), 177. (21) Edie, F. C.; Westerberg, A. W. Computer-aided Design: Part 3sDecision Variable Selection to Avoid Hidden Singularities in Resulting Recycle Calculations. Chem. Eng. J. 1971, 2, 114. (22) Stadtherr, M. A.; Gifford, W. A.; Scriven, L. E. Efficient Solution of Sparce Sets of Design Equations. Chem. Eng. Sci. 1974, 29 (4-I), 1025. (23) Book, N. L.; Ramirez, W. F. The Selection of Design Variables in Systems of Algebraic Equations. AIChE J. 1976, 22 (1), 55. (24) Drud, A. S. CONOPTsA Large-Scale GRG Code. ORSA J. Comput. 1994, 6 (2), 207. (25) Chen, X.; Rao, K. S.; Yu, J.; Pike, R. Comparison of GAMS, AMPL and MINOS for optimization. Chem. Eng. Educ. 1996, summer, 220. (26) Sargent, R. W. H.; Westerberg, A. W. SPEED-UP in Chemical Engineering Design. Trans. Inst. Chem. Eng. 1964, 42, T190. (27) Pho, T. K.; Lapidus, L. Topics in Computer-Aided Design: Part I. An Optimum Tearing Algorithm for Recycle Systems. AIChE J. 1973, 19 (6), 1170. (28) Ledet, W. P.; Himmelblau, D. M. Decomposition Procedures for the Solving of Large-Scale Systems. Advances in Chemical Engineering. Academic Press: New York, 1970; Vol. 8, p 185. (29) Steward, D. V. Partitioning and Tearing Systems of Equations. J. SIAM Numer. Anal. 1965, 2, Ser. B, 345. (30) Shacham, M. Decomposition of Systems of Nonlinear Algebraic Equations. AIChE J. 1984, 30 (1), 92. (31) Lee, W.; Christensen, J. H.; Rudd, D. F. Design Variable Selection to Simplify Process Calculations. AIChE J. 1966, 12 (6), 1104.

(32) Salcedo, R. L. Solving Nonconvex Nonlinear Programming and Mixed-Integer Nonlinear Programming Problems with Adaptive Random Search. Ind. Eng. Chem. Res. 1992, 31 (1), 262. (33) Kirkpatrick, S.; Gelatt, C. D.; Vechi, M. P. Optimization by Simulated Annealing. Science 1983, 220, 671. (34) Kirkpatrick, S. Optimization by Simulated Annealing: Quantitative Studies. J. Stat. Phys. 1984, 34 (5/6), 975. (35) Cardoso, M. F.; Salcedo, R. L.; Feyo de Azevedo, S. NonEquilibrium Simulated Annealing: A Faster Approach to Combinatorial Minimization. Ind. Eng. Chem. Res. 1994, 33 (8), 1908. (36) Floudas, C. A. Nonlinear and Mixed-Integer Optimizations Fundamentals and Applications; Oxford University Press: Oxford, U.K., 1995. (37) Umeda, T. Optimal Design of an Absorber-Stripper System. Ind. Eng. Chem. Process Des. Dev. 1969, 8 (3), 308. (38) Umeda, T.; Ichikawa, A. A Modified Complex Method for Optimization. Ind. Eng. Chem. Process Des. Dev. 1971, 10 (2), 229. (39) Zaldivar, J. M. Mathematical modelling and numerical simulation of aromatic nitrations by mixed and discontinuous reactors. Ph.D. Dissertation, University Twente, Twente, The Netherlands, 1995. (40) Westerwerp, K. R.; Van Swaaij, W. W. M.; Beenackers, A. A. C. M. Chemical Reactor Design and Operation; John Wiley: New York, 1984. (41) Esparta, R. Institut fur Systemdynamik und Regelungstechnik, University of Stuttgart, personal communication, 1998.

Received for review March 22, 1999 Revised manuscript received September 8, 1999 Accepted September 29, 1999 IE9901980

On the Optimum Choice of Decision Variables for Equation-Oriented

Recommend Documents