Ind. Eng. Chem. Res. 1992,31, 262-273
262
Serth, R. W.;Heenan, W.A. Gross Error Detection and Data Reconciliation in Steam-Metering Systems. AZChE J. 1986, 32, 733-742. Sidak, Z. Rectangular Confidence Regions for the Means of Multivariate Normal Distributions. J. Am. Stat. Assoc. 1967, 62, 626-633. Stengel, R.F. Stochustic Optimal Control: Theory and Application; Wiley: New York, 1986. Tamhane, A. C.; Mah, R. S.H. Data Reconciliationand Gross Error Detection in Chemical Process Networks. Technometrics 1986,
27, 409-422. Watanabe, K.; Himmelblau, D. M. Instrument Fault Detection in Systems with Uncertainties. Znt. J. Syst. Sci. 1982,13,137-168. Willsky, A. S.;Jones, H. L. A General Likelihood Ratio Approach to State Estimation in Linear System Subject to Abrupt Chauges. R o c . ZEEE Conf. Decision Control 1974,846-853.
Received for review January 11, 1991 Revised manuscript received July 29,1991 Accepted August 15, 1991
Solving Nonconvex Nonlinear Programming and Mixed-Integer Nonlinear Programming Problems with Adaptive Random Search R. L. Salcedo Departamento de Engenharia Quimica, Faculdade de Engenharia da Universidade do Porto, Rua dos Bragas, 4099 Porto, Portugal
An adaptive random-search optimization algorithm is presented, which is found to be efficient in dealing with nonconvex nonlinear programming (NLP) problems and mixed-integer nonlinear programming (MINLP) problems. The algorithm is an extension of a previously reported algorithm which successfully solved for the global optimum of difficult multimodal unconstrained and constrained NLP problems. The proposed algorithm treats integer variables as discrete integer Variables, thus avoiding approximations which could produce either infeasible results or a wrong solution to the optimization problem. The algorithm does not require any problem decomposition, neither identification or elimination of nonconvexities, prior to its application. It was tested with several nonconvex NLP and MINLP problems published in the literature, and the results obtained reveal its adequacy for the optimization of nonconvex NLP and small to medium scale MINLP problems encountered in chemical engineering practice. A simple two-phase relaxation strategy is proposed which overcomes some of the difficulties associated with the solution of ill-conditioned or larger scale MINLP problems.
Introduction A great deal of attention has been recently drawn to the problem of efficiently solving nonconvex nonlinear programming (NLP) problems and mixed-integer nonlinear programming (MINLP) problems (Grossmann and Sargent, 1979; Duran and Grossmann, 1986; Kocis and Grossmann, 1987,1988;D o h et al., 1989;Floudas et al., 1989;Floudas and Ciric, 1989;Kokwis and Floudas, 1990, Ciric and Floudas, 1990;Cerdd et al., 1990;Ostrovsky et al., 1990). Nonconvex NLP problems arise in many situations in chemical engineering optimization, due either to the objective function or to the constraints imposed on the particular problem. Although a large number of variables and constraints may be present in NLP problems, these tend to have a moderate number of degrees of freedom or decision variables (generally below 10). However, they may be difficult to solve for the global optimum, due to the presence of several local optima. The currently available algorithms for nonconvex NLP problems can lead to suboptimal solutions since nonconvexities may cut off the global optimum from the current search regions (Wang and Luus, 1978; Kocis and Grossmann, 1988;Floudas et al., 1989). MINLP problems occurring in process synthesis, viz., network-type problems and design of multiproduct batch plants, are usually large-scale problems which can have tens of degrees of freedom, since the set of discrete variables can be very large (Kokossis and Floudas, 1990;Ciric and Floudas, 1990). "heae variables are commonly discrete because of the requirements of a set of standard sizes and commonly integer because an undetermined number of
identical elements (processing units) are to be used in the design (Siddall, 1982; Edgar and Himmelblau, 1988). Smaller scale MINLP problems, having usually less than about 20 degrees of freedom, may also occur in chemical engineering practice (Groasmann and Sargent, 1979;Cerd6 et al., 1990). For the solution of MINLP problems, there seems to be no single procedure which can claim a uniform advantage over all others for all problems, or even claim to be generally effective (Ray and Szekely, 1973;Grossmann and Sargent, 1979;Edgar and Himmelblau, 1988). Recently, most of the work related to the solution of nonconvex MINLP problems has been based either on the generalized Benders' decomposition method (Geoffrion, 1972)or the outer approximation of Duran and Grossmann (1986).In both methods the original problem is partitioned into two subproblems (the primal and the master), and iteration between the subproblems is carried out to arrive at the final solution. The primal is the original NLP problem solved for a fixed set of the integer variables y, and the master (a mixed-integer linear problem or MILP) is solved for the continuous variables x to provide new values for the integer variables. In the generalized Benders' decomposition method, the master problem is based on the dualization of the original problem, while in the outer approximation it is based on linearization of the objective function and constraints. In any case,partitioning is made such that the integer variables y may be solved for independently of the continuous variables x. Also, since a number of constraints must usually be evaluated prior to the solution of the problem, the master problem is solved as a series of relaxed subproblems. Kocis and Grossmann
0888-5885/ 9212631-0262$03.oO/0 0 1992 American Chemical Society
Ind. Eng. Chem. Res., Vol. 31, No. 1, 1992 263
(1987)have further adapted the OA algorithm to deal explicitly with equality constraints (the OA/ER algorithm). For the global solution of nonconvex NLP and MINLP problems, Kocis and Grossmann (1988)have applied the OA/ER algorithm in a two-phase strategy; in the first phase, nonconvexities are identified by local/global tests while in the second phase a relaxation is systematically imposed on the invalid solutions (outer approximations) in an attempt to locate an improved solution. This approach may fail to identify nonconvexities in the first phase, although in many cases the two-phase strategy can identify the global MINLP optimum. Floudas et al. (1989)have suggested a different approach, in which partitioning of the original problem is made in such a way as to ensure that the global solutions to the primal and master subproblems are attained at each iteration. This is performed by first identifying the sources of nonconvexities by partitioning the variable set into complicating (viz., responsible for nonconvexities) and noncomplicating variables and by partitioning the nonconvex constraint set. Formulation of the master problem through dualization of the primal problem and iterating between these subproblems follows (viz., through the application of the generalized Benders’ decomposition method), until no further improvement can be obtained. Even though the proposed approach cannot guarantee the identification of the global optimum, Floudas et al. (1989) and Floudas and Ciric (1989)have obtained the global solution for a number of nonconvex NLP and MINLP test problems, starting from different search vectors. Dolan et al. (1989)have used simulated annealing to optimize process network design problems that have combinatorially large numbers of feasible designs. Based on the theory of Markov chains, simulated annealing seems capable of surpassing local optima and arriving close to the global optimum. Algorithms based on simulated annealing have been successfully applied to the optimization of pressure relief header networks and heat exchanger networks (Dolan et al., 1989)and coupled with heuristic search to the traveling salesman problem (Press et al., 1986). However, the success of simulated annealing depends strongly on the cooling schedule and on the problem formulation, and only asymptotically is there guarantee of convergence to the global optimum ( A m and Korst (1989)). Grossmann and Sargent (1979)and Ostrovsky et al. (1990)have utilized the branch and bound (BB) method usually applied to the solution of mixed-integer linear programming (MILP) problems to solve medium-scale MINLP problems. This is a partitioning multistep procedure involving the solution of consecutive NLP subproblems with appropriate branching criteria. The BB method requires efficient branching to be effective, which is problem dependent (Edgar and Himmelblau, 1988)and may become computationally intensive (Grossmann and Sargent, 1979). Motivation The general objective of MINLP problems is to maximize or minimize a function of the form
F = f(x,y)
(1)
subject to the constraints gj(x,y) 1 0, j = 1, 2, ..., m ar I xj I
yi I yi 5
pj, i = 1, 2, ..., p &, i = p + 1, ..., n
(2)
(3) (4)
where x is a vector of continuous variables, y is a vector of integer variables, (ai,Bi) are bounds on parameter lci, and (Ti,ti) are bounds on parameter yi, e.g., {O,l)for binary variables. Although it is tempting to apply standard nonlinear optimization algorithms to the solution of MINLP problems by rounding the components of the solution vector y to the nearest integer values, this should be avoided, since it may lead to solutions that are not feasible or otherwise do not correspond to the optimum (Siddall, 1982; Edgar and Himmelblau, 1988;Ostrovsky et al., 1990). Rounding becomes critical for problems where the discrete Variables can only take few values, typically below about 8 (Siddall, 1982). The discrete variables can usually be treated as continuous only for many values, typically above 20 (Siddall, 1982;Ostrovsky et al., 1990). Apart from the algorithms based on simulated annealing, the various methods described above share one common characteristic, which is the need to solve a series of subproblems obtained from an appropriate decomposition of the original problem. The algorithms based on the generalized Benders’ decomposition or the outer approximation also share the necessity to identify the sources of nonconvexities, and this task may not be easy. Nonconvexities for continuous variables may be determined by solving a perturbed NLP problem around a local solution, by performing local/global tests or by exploiting the available theoretical properties for the form of the objective function and constraints (Kocis and Grossmann, 1988; Floudas et al., 1989). It may be possible, under certain circumstances, to eliminate nonconvexities by appropriate transformations (viz., logarithmic transformations; see Kocis and Grossmann, 1988). In any case, the presence of nonconvexities and the necessity of their identification or elimination may considerably aid in the work involved in solving a specific NLP or MINLP problem. The application of random search to the optimization of nonconvex NLP problems has been described in the literature (Luus and Jaakola, 1973a,b;Luus, 1975;Gaines and Gaddy, 1976;Heuckroth et al., 1976;Wang and Luus, 1977,1978;Price, 1978;Torn, 1978;Martin and Gaddy, 1982;Fey0 de h v e d o et al., 1988;Luus and Brenek, 1989; Salcedo et al., 1990),among others. Although random search has been coupled with heuristic algorithms to provide initial estimates for the solution of combinatorial minimization problems (see for example Reiter and Sherman (1965)and Lin and Kernighan (1973)),to the author’s knowledge the application of direct or adaptive random-search methods per se to the solution of MINLP problems has not been presented. This paper describes the application of an adaptive N L P random-search optimization algorithm (the MSGA algorithm) to the solution of nonconvex NLP and MINLP problems, with the emphasis placed on the latter. The proposed algorithm does not require the prior step of identification or elimination of the sources of nonconvexities, neither the decomposition of the original problem into subproblems which have to be iteratively solved for. The algorithm was tested with several nonconvex NLP and small to medium scale MINLP problems given in the literature, showing a good reliability in attaining the global optimum with a high accuracy. A simple two-phase relaxation strategy extends the applicability of the randomsearch algorithm to the solution of ill-conditioned and larger scale MINLP problems. The MSGA Algorithm For many years the efficiency of optimization routines has been closely linked to the number of iterations involved
264 Ind. Eng. Chem. Res., Vol. 31, No. 1, 1992
in reaching the optimum within a desired accuracy. Today, with the growing availability of inexpensive powerful computational means, the attention is more and more on robustness, i.e., the ability of routines to overcome local optima and their sensitivity to starting conditions. Fey0 de Azevedo et al. (1988) and Salcedo et al. (1990) have presented an adaptive random-search algorithm (the SGA algorithm) which is based on the random-search strategy first proposed by Luus and Jaakola (1973a,b) and subsequentlydeveloped by Gaines and Gaddy (1976). The main features of the SGA algorithm are its robustness and reliability in attaining the global optimum. This is due to ita high insensitivity to the initial estimates of the parameters, search intervals, and random number sequence, viz., the optimization path. The main objective of this work is to provide an extension of the algorithm to include discrete decision variables. Basically, the algorithm generates a sequence of random test points, which may be continuous (with NLP problems), integer (discrete), or of both types (with MINLP problems). The continuous test points are generated centered at their current optimum, whereas the only requirement for the generation of the integer test points is that they are in the range of the allowable values (given by eq 4). During the course of the optimization the search regions for the continuous variables decrease after generation of a predetermined number of test points (referred as one global iteration). The search regions for the integer variables remain constant and equal to their initial values. The algorithm ends after a predetermined number of global iterations are performed (depending on the desired accuracy) or whenever built-in convergence criteria for the independent parameters and the objective function are met. The proposed algorithm (named MSGA for MINLP SGA) can be summarized in the following steps. (a) Iteration j is executed where Pl trials are carried out, with the continuous test points x established by
x = x* + ZrG)
(5)
and the integer test points y established by y =y
+ INT {Ur(l)+ v)
where x* is the current optimum for the continuous variables, y is the lower bound vector for the integer variables (see eq 4), Z and U are square matrices of dimensions p and n - p, respectively, consisting of random numbers in the intervals ]-1/2,1/2[ and ]0,1[, rv) is the search region ..., 1/2] is a constant vector of for iteration j, vT = dimension n - p, and INT denotes the integer operator, viz., returning the largest integer less or equal to the operand. It is important to notice that the search regions for the integer variables are invariant during the course of the optimization, since r(l)in eq 6 refers to the initial search region vector. This provides an opportunity for the integer variables to span their range of possible values while the search regions for the continuous variables decrease. To generate discrete (not necessarily integer) values, eq 6 can still be employed, provided a one-to-one correspondence is made between the integers generated by this equation and a table of the allowable discrete values. (b) After each success, (x*,y*)are replaced by the best feasible values, and, at the end of PI trials, the search region is contracted P2times for the continuous variables only, following ,iv+l)
=
- ci)
(7)
where ri represents the search domain for parameter i and ei its component of the compression vector e. In detail, the overall proposed algorithm may be described as follows. (i) It offers a choice of three “levels of accuracy” (1,2, or 3) which correspond respectively to P2 = 100, 200, or 400 (maximum) global iterations (interval reductions), where in each iteration Pl = 100 trials are typically executed (see (vi) below). (ii) For constrained problems, it incorporates an initial search on a coarse grid for the choice of a feasible “bad” starting point, to minimize the probability of starting near a local optimum. If a feasible starting point cannot be obtained, it starts the search from the midpoint of the search intervals. A feasible or infeasible initial estimate may also be provided. (iii) It defines automatically the search intervals for contrained problems by the bounds on the independent parameters (eq 3 and 41, or they may be chosen by the user for unconstrained problems. (iv) It updates each component of the compression vector t (initial values ci = 0.05) at the end of each global iteration, assuming a value of the discrete domain (0.1,0.05,0.025,0.02,0.01]. This is a function of the relative variation of the optimization parameters between two consecutive successes and of the number of iterations without successes. The key idea is to decrease the search regions for parameters which do not vary much between consecutive successes at a faster rate, since this increases the probability (for a fixed number of Pl random trials), to escape the current solution vector, i.e., optima are systematically assumed as local. For details on the selective reducing schedule employed for t, the reader is referred elsewhere (Salcedo et al., 1990). (v) It shifts the current solution vector if after a number of interval reductions (currently set at 20 for the initial estimate and at 50 for an improved solution) there is no improvement in the objective function. This avoids locking the search around highly stable local optima. The search is restarted from a random estimate, where for unconstrained problems both the search intervals and compression factors are doubled from their original values. A “local character” is assigned to the current optimum, and this “local flag” is carried throughout the entire optimization until a better optimum is obtained. This strategy allows for “wrong-way” moves (i.e., “uphill” for minimization or “downhill”for maximization) which take the trajectory away from the local optima. (vi) “Level of accuracy 1” should lead to a coarse determination of the global optimum. “Level of accuracy 2” provides a better estimate of the objective function determined at the end of level l, through some expansion of the search intervals (Heuckrothet al., 1976; Salcedo et al., 1990) followed by further interval reductions. “Level of accuracy 3” starta by assuming that the optimum at the end of level 2 is a local optimum and as such assigns a poor objective function to the local solution vector. The algorithm then proceeds for another 100 iteractions. However, if the local optimum flag of level l is activated, then the compression factors are decreased by 50% of their current values, with a lower bound of 0.01, to provide a more careful search. This strategy is extremely helpful for the optimization of objective functions which exhibit many local optima (viz., nonconvex problems). At the end of 300 global iterations the current optimum is compared to the local optimum obtained at the end of level 2, and the best optimum is retained for further interval reductions, up to a maximum of 400.
Ind. Eng. Chem. Res., Vol. 31, No. 1, 1992 265 (vii) It performs a test at the extremes of the search intervals at the end of each level of accuracy, since the optimum for constrained problems is often found at one or more of the interval extremes. (viii) It incorporates a dynamic relaxation strategy for the inequality constraints, whereby the relaxations are decreased during the c o m e of the optimization. This provides a method for the algorithm to proceed toward feasibility whenever a feasible solution cannot be determined easily (viz., with heavily constrained problems). A similar relaxation method is also employed by Kokossis and Floudas (1990) as a first step toward the initialization of large-scale problems. The proposed relaxation strategy is described in detail below, as applied to the solution of examples 8 and 9. (ix) Finally, it includes stopping criteria based on convergence factors for the objective function and for the parameters. To decrease the probability of convergence on a local optimum, these criteria must be simultaneously met and are only operational for level 3, after 300 interval reductions, if a reasonable number of global iterations (20) is counted between two consecutive successes. Overall Test Conditions The MSGA algorithm was written in Fortran 77 (the code is available from the author) and run in an Apple Macintosh SE microcomputer with a Radius board (Motorola 68020 CPU and 68881 FPU both at 25 MHz, from Radius Inc., San Jose, CA). All floating point calculations were performed in extended precision (80 bits total corresponding to 64 bits mantissa). Since direct random-search methods can be sensitive to the sequence of pseudo random numbers (Martin and Gaddy, 1982),the pseudo random number generators employed within the computer programs that implement these methods should be uniform and independent, so that deficiencies in the random-search algorithms may not somehow be compensated by nonuniformity. In this work, a Lehmer linear congruential generator (Shedler, 1983) was used. This generator was tested for multidimensional uniformity and found to be satisfactory (Salcedo et al., 1990). To evaluate the performance of the MSGA algorithm, several nonconvex NLP and MINLP problems were selected from the literature. Since the MSGA algorithm is a random-search algorithm, its reliability has to be ascertained on a statistical basis. For this purpose, sets of 100 different runs with different seeds, taken from Shedler (1983), were employed for each problem and each starting solution vector. Since the test functions examined are all constrained, the initial search regions are automatically defied by the bounds on the independent parameters. The initial estimates of the parameters are those proposed by independent authors plus the optional default value from the SGA algorithm whenever this originates a supplementary starting point. Hence, no bias in performance should be attached to the choice of the starting vector or the search intervals. As the extremes of the intervals are automatically searched for in the MSGA algorithm, the global optimum may be found for some problems irrespective of the effectiveness of the random search. Thus, the test on the extremes was deactivated for all problems to provide a correct appraisal of the reliability of the random-search step of the MSGA algorithm. The maximum number of interval reductions was set to 400 (this corresponds to level of accuracy 3 for the SGA/MSGA algorithms), and the number of random trials
Table I. Initial Solution Vectors and Results for Example 1 X1
F
0 -10 0.5 -15.6094 1 -12 2 infeasible
BUCCWMS
x2
F
100 100 100 100
1.5 1.875 0 -30
-16.735 -15.6094 -12 infeasible
successes 100 100 100 100
typical CPU time per run = 11 s
per interval reduction was set to 100. The convergence factors for the relative error in the objective function and independent parameters were set respectively to and to (root mean square deviation). These values were the ones employed previously with the SGA algorithm (Salcedo et al., 1990) for comparison with other randomsearch algorithms. Results and Discussion Three different types of problems were employed to test the reliability of the MSGA algorithm. The first set consists of three nonconvex NLP problems, and the second set consists of four small-scale (less than 10 degrees of freedom) MINLP problems. All these problems are taken from the literature, and all have been used by other authors to test different types of algorithms. Whenever starting points were given by the different authors, these were retained for use with the MSGA algorithm. All problems in these two sets were run with 100 different seeds from each starting point, and the number of successes refer to the percentage of runs that successfully attained the global optimum. Whenever appropriate, a measure of accuracy (defined here as the relative deviation of the objective function from the global optimum) is also provided, giving either the lowest accuracy achieved at any trial or the percentage of trials within a certain accuracy, whichever the case. The third set consists of two larger scale MINLP problems (respectivelywith 10 and 22 degrees of freedom). This serves the purposes of identifying sources of possible problems when a random search is applied to more complex situations, of establishing a strategy to overcome these problems, and possibly of understanding when random search is not adequate. Nonconvex NLP Problems. Example 1. Consider the following simple problem, taken from Soland (1971):
min -12x1 - 7x2 + x22 5.t.
-2x14 + 2 - x 2 = o 04x112 05x243
This problem has been recently used by Floudas et al. (1989) to demonstrate an efficient partitioning method in the search for the global optimum of nonconvex NLP and MINLP problems. The only source of nonconvexity is the 4th order term in the equality constraint. Since this problem has 1degree of freedom, choosing x1 as the decision (independent) variable will eliminate the nonconvexity. The global optimum is at (xl,x2;FJ = {0.718,1.47;-16.7389). The MSGA algorithm was run starting from four different solution vectors (the default from the automatic search built into the algorithm plus three starting (xl) points employed by Floudas et al. (1989)),which both x1 and x2 (in turn) as decision variables. Table I shows that the global optimum was obtained in all cases.
266 Ind. Eng. Chem. Res., Vol. 31, No. 1, 1992 2
min CaA;
(8)
i=l
Figure 1. Superstructure for example 3. Table 11. Initial Solution Vecton and Reroltr for Example 2 x1 *2 4 F SUCCBLU)BB 46 accuracy 0 0 1 3 100 0.4 0 4 0 -3.03594 100 0.4 4 4 0 infeasible 100 0.3 1 3 1 infeasible 100 0.6 3 4 0 infeasible 100 0.5 4 0 1 infeasible 100 0.5 typical CPU time per run = 36 s
Example 2. Consider the following example from Stephanopoulos and Westerberg (1975):
min
+ xp6- 621 - 4ul + 3Uz 32
8.t.
- 321 - 3Ul = 0
+ 2241 I 4 2 2 + 2uz I4 21
05x153 OIUZI1 22,
u1 1 0
This problem has 3 degrees of freedom, has a nonconvex objective function, and exhibits a dual gap (Floudas et al., 1989). Thus, solution to the primal cannot be obtained through dualization (Gro~smannand Sargent, 1979). An obvious choice of decision variables is the set (xl,z2,uz) where 0 I x z I4. The global optimum is at (xl,rz,uz;F(1 = (4/s,4.0;-4.5142). Table I1 shows that the MSGA algorithm also solved this problem for the global optimum in every trial from all starting points. Example 3. We turn now our attention to a heat-exchanger network synthesis example taken from Floudas and Ciric (1989) and Floudas et al. (1989). A stream superstructure of one cold stream and two hot streams for two heat exchangers is given in Figure 1. The problem consists of obtaining the minimum investment cost network confiiation, given a fixed set of process stream matches and a fixed number of units. This superstructure contains all the poesible splittar/mixer flow configurations. The objective function, subject to constraints for maas and heat balances, is Table 111. Input Data for Exomple 3 stream
Ti,, K
s.t. T! ITjITju;j = 1, 4 (9) where Ai = Qi/Ui(LMTDi)is the heat-transfer area of are exchanger i, (a,@)are cost coefficients, and Tj'and Tju appropriate lower and upper limits for the temperature of stream j. This objective function is nonconvex since j3 = 0.6. The problem equations are the stream mass balances at the splitter at the inlet of the superstructure, around the mixers before the heat exchangers,and around the splitters after the heat exchangers and the energy balances at the mixers before the exchangers and at the exchangers, with a total of 9 equations in 18 variables. With a known feed (C,Th) and since Q and U are known, this problem has 3 degrees of freedom. However, it does not have a serial solution, and the four energy equations must be simultaneously solved at each step of the optimization. Also, since the energy equations at the exchangers involve bilinear terms (in flow rate and temperature) of opposite sign, two of the constraints are nonconvex as well. To solve this optimization problem, the flow rates (F2$3$4) were chosen as decision variables and the four simultaneous linear equations were solved by standard LU decomposition (Press et al., 1986). Although this choice of decision variables eliminates nonconvexities in the constraints, nonconvexities in the objective function still remain and local optima may occur. Table III gives the input stream and match data for the two cases studied, taken respectively from Floudas and Ciric (1989) and Floudas et al. (1989). Parallel configurations were used as starting points to the optimization process. These correspond respectively to the solution vectors (F2,F3,F4;FJ= (7.059,2.941,7.059;81266) and (8,4,8;67320) for case 1, and (FZ,F3,F&F) = (5.556,4.444,5.556;13874) and (4,4,4;17297) for case 2. The inequality constraints of this problem are the temperature constraints (eq 9) with a minimum temperature difference of 10 K (Floudas and Ciric, 1989) first case 150 I Ti I240 Ti + 10 IT2 I 490 150 IT3 I190 T3 + 10 ITd I340
second case 100 I TI I 290 Ti + 10 I Tz I 310 100 I T3 I 290 T3 + 10 I TI I 330
and the logical bounds on the flow rates OIFiIC; i=l,8 To compare the results with those from Floudas and Ciric (1989) and Floudas et al. (19891,the LMTD was estimated with Paterson's approximation (Paterson, 1984). The global optima correspond to (FZ$3$#j = (10,10,10;56826) and (10,10,10;12293) respectively for cases 1 and 2. Table IV shows the results obtained from the two different s t a r t i n g points (the initial solutions suggested by Floudaa and Chic (1989) and Floudas et al. (1989) plus the two default solutions from the MSGA algorithm) for each
T-t,K
FiC,, kW K-'
match
U,kW m-*K-'
4 4 10
HIC H2C
0.05 0.05
40 25 10
HIC HZC
2.5 0.2
case 1 a = $1300
H1
500
H 2
C
350 150
H1 H2
320 340 100
250 200 310
case 2 a = $1200
C
300 300 280
Ind. Eng. Chem. Res., Vol. 31, No. 1, 1992 267 Table V. Initial Solution Vectors and Results for Example 4 _____ XI X2 Y F successes % accuracy 1.4 0 3.15667 100 0.004 0.356675 0 3.15667 100 0.001 ~~
~
typical CPU time per run = 25 s Table VI. Initial Solution Vectors and Results for Examule 5 _ _ _ _ _ _ _ _ ~ ~ ~ Y1 YZ YS Y4 F successes 1 0 60 100 1 1 0 0 1 1 -6 100 1 0 1 1 -3 100 0 1 -1 100 0 1 0 1 0 100 1 0 1 1 0 1 2 100 0 1 1 1 8 100 20 100 1 1 1 1 ~
~---.-..-..--....-.--.Y 310K,
typical CPU time per run = 32 s
vious work for a more complete appraisal of the MSGA algorithm in solving nonconvex NLP problems. Small-Scale MINLP Problems. Example 4. This example is taken from Kocis and Grossmann (1987) and involves a small MINLP problem with 2 degrees of freedom and an equality constraint: min -y + 2x1 x 2 x1 - 2 exp(-x2) = 0 s.t. -x1 + x2 y I0 0.5 I~1 I1.4
+
',ooKp J +
1
.-....- -.
$00
K 4 kW K - '
4.021 KW K-'
'
1 0 k W K-l 150K
j c ; i
$--.-
..
...._..__._.._.___ Y 310K,
y..-~------.-.-.-......-...
3
Y E (091) It was solved with the MSGA algorithm with both (xl,yJ and (x2,y)as decision variables. Table V shows that the Optimum ( x ~ , x ~ ~ = ; F{1.375,0.375,1;2.124) J was obtained in every trial starting from the default solution vectors of the MSGA algorithm. Example 5. The next example, which represents a quadratic capital budgeting problem (Kocis and Grossmann, 1988), features only integer variables. min (Y1+ 2 ~ + 2 3 ~ -3y4)(2y1 + ~ Y + Z 3 ~ -3 6 ~ 4 ) 5.t. Y 1 + 2Y2 + Y3 + 3Y4 2 4 Y E (0,114
This is a nonconvex problem with 4 degrees of freedom, where the objective function contains bilinear terms of opposite sign. Table VI shows that the MSGA algorithm found the optimum (y1,y2,y3,y4fl= (0,0,1,1;-6] in every case. Example 6. The following problem was taken from Kocis and Grossmann (1988) and has three continuous variables and two binary variables. Due to the equality constraints, which are the sources of nonconvexities (Floudas et al., 1989), this problem has only 3 degrees of freedom. min 2x1 3x2 + 1 . 5 + ~ 2y2 ~ - 0.5y3 x12 + y1 = 1.25 5.t.
+
+ 1.5y2 = 3.00 x1 + y1 I 1.60 1.333X2 + y2 I3.00 X21'5
-Y1
- Y2 + Y3 I 0 XLO
Y E (0,113
268 Ind. Eng. Chem. Res., Vol. 31, No. 1,1992 Table VII. Initial Solution Vectors and &SUlts for Example 6 Y2 Y3 F successes Y1 1 0 0 8.74025 100 0 1 1 7.667 18 100 0 1 0 8.167 18 100 0 0 0 8.476 32 100 0 0 1 infeasible 100 1 1 1 7.931 11 100 1 0 1 8.24025 100 1 1 0 8.431 11 100 x1 x2 Y3 F successes % accuracy 0.001 0 infeasible 100 0.64 1.8008 0.01 1 infeasible 100 0.80 1.1260 0.015 0 infeasible 100 1.60 2.251 0.03 1 infeasible 100 0 0
The MSGA algorithm was applied to this example with two sets of decision variables, respectively ( y 1 ~ 2 ~ and 3 ) ( ~ l ~ z ~ For 3 ) .this last case, the continuous variables were searched in the interval, 0 Ix1 I1.60 and 0 I x 2 I 2.251, as can easily be derived from the constraints. Table VI1 shows that the algorithm was capable of obtaining the global Optimum ( x ~ , z ~ J ~ ~ J ~ J = ~ ;{F1.118,1.310,0,1,1~7.667) J in every situation and with a high accuracy. It should be pointed that the OA/ER algorithm applied with the two-phase strategy of Kocis and Grossmann (1988) failed to obtain the global optimum starting from ( y y 1 , ~ 2 , ~ 3= ) (l,O,O), whereas the decomposition method of Floudas et al. (1989) succeeded. Example 7. As a last example of a small-scaleMINLP problem, consider the problem taken from Yuan et al. (1989) which was successfully solved by Floudas et al. (1989). It features nonlinearities in both continuous and integer variables and has 7 degrees of freedom. min (yl + (yz - 2)2 + (y3 - - In (y4 + 1) + ( x 1 - 1 ) Z + ( x 2 - 2)Z + ( x 3 - 3)2 y1
+ yz + y 3 + x1 + 2.2 + x 3 I 5 y32
+ + I5.5 + I1.2 y2 + ~2 I 1.8 + I2.5 + I 1.2 yz2 + x z 2 I1.64 y32+ x32 I4.25 yZ2+ x32 I4.64
+
XI2
x22
y1
x1
y3
x3
y4
x1
x32
x10 Y E {0,114 The search intervals for the continuous variables were derived from the problem constraints to be 0 Ix l I1.2, Table VIII. Initial Solution Vectors and Results for Example 7 Xl 3c2 x3 Y1 Yz 0
0
0
0
0 0 0 3 1
0 0 0 3 1
0 0 0
1 1 0 0 1
3 1
Point 2 0
Point 3
+
Point 4 Point 5
0
! . ....... . ...... . .....
.oo1
I
1
.o1
.1
1.
7
1
.
.
.......I
10
..m
1 0
95 Accuracy
Figure 3. Comparison of the performances of different starting points for the optimization of example 7.
typical CPU time per run = 32 s
s.t.
d
0 0 1 1 0 1
0 Ix 2 I 1.281, and 0 I x 3 I 2.062. Table VI11 shows the starting solution vectors and the results obtained from the application of the MSGA algorithm. The two runs without success gave the same local optimum which correspondsto the best solution reported by = Yuan et al. (1989), with (x~,xz,~~,Y~,Yz,Y~,Y~;FJ (0.2,0.8,1.908,0,1,0,1;5.5796) and is about 22% larger than ~4.~ the true global optimum of ( x 1 , 3 ~ 2 , x 3 , y l , y 2 , ~ 3 , = (0.2,0.8,1.908,1,1,0,1;4.5796). It is interesting to note that both the local and global optima share the same values for the continuous variables. Figure 3 shows the number of trials which reached a value of the objective function within a percentage of the global optimum for the six starting points of Table VIII. The fist point in Table VI11 corresponds to the default from the MSGA algorithm and in this case shows a somewhat greater accuracy. Although the default gives usually good reaulta (see Salcedoet al., 19901,the better behavior shown in Figure 3 for the default solution vector should not be generalized to other problems. Larger Scale MINLP Problems. The MSGA algorithm is now tested with larger scale and more complex MINLP problems. The proposed problems were taken from Grossmann and Sargent (1979) and Kocis and Grossmann (1988) and refer to the optimization of a multiproduct batch plant. The plant consists of M processing stages in series where fixed amounts Qi of N products have to be manufactured. The objective is to determine for each stage j the number of parallel units Nj and their sizes and for each product i the corresponding batch sizes Bi and cycle times Tu The problem data are the horizon time H, the size factors Sij and processing times tij of product i in stage j , the required productions Qi, and appropriate cost coefficients aj and & The mathematical formulation of the optimum design of multiproduct batch plants is as follows (Grossmann and Sargent, 1979; Kocis and Grossmann, 1988) M
min C a j N j V p j=1
Y*
Y4
0 1 1 1 0 1
0 0 1 0 0 1
typical CPU time per run = 1 min 8 s
F 20 18 14.3069 16 infeasible infeasible
successes 100 100 100 100 99 99
Ind. Eng. Chem. Res., Vol. 31, No. 1, 1992 269 set.
Vj 1 SijBi NjTLi 1 tij 1 INj INj"
vi' I vj Iv j u TL) I TLi ITLr B) I Bi IBi" N j integer The bounds Nju,V), and Vjuare specified by the problem and appropriate bounds for TLi and Bi can be determined as follows
Qi B) = TL)
H
Biu = min
The optimization of a multiproduct batch plant with the data structure as formulated above has 2(M + N) degrees of freedom, 4(M + N) bound constraints, 1(I inequality ) constraint and 2MN (1) inequality constraints. The objective function is nonconvex as long as the cost factors pi are less than unity. The horizon time constraint and several inequalities (see Kocis and Grossmann (1988)) are also nonconvex. Thus, depending on the values of M and N , a fairly large nonconvex MINLP problem may result. Also, the possible range for N j together with the number of stages M will determine the number of NLP subproblems embedded in the original MINLP formulation. Table IX gives the data chosen to represent two examples of increasing complexity, taken from Grossmann and Sargent (1979) and Kocis and Grossmann (1988). In all cases, H = 6000 h, CY, = $250, and 8. = 0.6 0' = 1, M). Example 8. This problem has 10 degrees of freedom, 20 bound constraints, and 13 inequality constraints. A systematic search over the space of the integer variables N j would require solving 27 NLP subproblems embedded in the original MINLP formulation. Table X shows the default solution vector employed as initial estimate and the results obtained by running the MSGA algorithm with 100 different seeds. Also shown is the global optimum as obtained by Grossmann and Sargent (1979). The global optimum was never obtained, and instead three local optima were obtained. These local optima are similar from the point of view of the objective function, and the best local optimum was obtained most of the time. The MSGA algorithm was rerun from another starting point and similar results were obtained, viz.,the global optimum always escaped detection. To solve this problem and understand the observed behavior, a two-phase relaxation strategy is imposed on the MSGA algorithm. Phase I is an iterative procedure with the objective of obtaining successive infeasible solutions which are closer and closer to feasibility, thus identifying candidates for the global optimum. Phase I1 simply solves the corresponding NLP subproblems. This twophase strategy can be summarized as follows. Phase I. (1) This phase completely relaxes the ineand 2)to provide a relaxed solution quality constraints (I to the problem (this relaxed solution is not necessarily the
same as the unconstrained solution since the bound constraints are always considered active). (2) It determines the maximum violations on both typea of inequality constraints and considers these as relaxations to the following run. (3) Since one wants to proceed toward feasibility, it reduces the relaxations during the optimization process by a fixed amount, e.g., 1% per interval reduction, which is automatically performed by the MSGA algorithm. It iterates between steps 2 and 3 until a feasible solution or until no further improvement is obtained. (4) It repeats steps 1-3 with different optimization paths (only necessary for very difficult problems). Phase 11. This phase solves the corresponding NLP subproblems and considers the best feasible solution as the global optimum of the original MINLP problem. In this work two relaxation reduction schedules were employed, viz., 1.145 and 0.346% per interval reduction (global iteration), giving overall maximum reductions with the MSGA algorithm of 90 and 50%, respectively. The identification process of the candidates for the global optimum (phase I) is straightforward as long as feasible solutions have been obtained. However, if it terminates with infeasible solutions, then the possible candidates for the global optimum have to be estimated with a more robust process. This involves providing different optimization paths to the MSGA algorithm, starting from the least infeasible solution obtained. It should be pointed out that phase I should always proceed with the original search regions for the independent parameters, so that feasibility regions where the global optimum might be are not prematurely cut off from the search space. Table XI shows the results from the application of the relaxation strategy described above to example 8, employing a relaxation reduction schedule of 1.145% per interval reduction and three different optimization paths (these correspond to the three first seeds provided by Shedler (1983) for the pseudo random number generator employed in thiswork). The objective of employing three different seeds is to provide some indication of the sensitivity of the relaxation strategy with different optimization paths. Table XI shows that feasibility was always obtained with the set of integer variables (2,2,1),in agreement with the results of Table X. Phase I suggests the following can(2,2,1),and (1,2,1); didates for the global optimum, (l,l,l), i.e., in this case all integer sets obtained are candidates. Phase I1 solved the corresponding NLP subproblems, giving the global optimum at N = (l,l,liT,V = (480,720,960)T,B = (240,120)T,and TL = (20,16jT,with F = 38499.8. These values agree with those found by Grossmann and Sargent (1979) and correspond to an extremelly ill-conditioned optimum, since variations as small as 0.01 % in any of ita independent continuous parameters produces infeasibility. This is the reason why the MSGA algorithm did not directly find this optimum. The relaxation strategy, however, seems adequate to find such illconditioned points,where the global optimum was detected from all three seeds. A similar problem with 14 degrees of freedom corresponding to 81 NLP subproblems (Groeemann and Sargent, 1979)was also sumeadully solved for the global optimum (againan extremely ill-conditioned point) with the MSGA algorithm and the two-phase relaxation strategy described above. Example 9. This last problem has 22 degrees of freedom, 44 bound constraints, and 61 inequality constraints. A systematic search over the space of the integer variables N j would require solving 4096 NLP subproblems embed-
270 Ind. Eng. Chem. Res., Vol. 31, No. 1, 1992 Table IX. Input Data for the the Multiproduct Batch Plant
v,.
example 8
M 3
N
NY
2
3
Vi 250
2500
9
6
5
4
300
3000
example 8 9
4j
Sij
2 3 4 4 6 3 7.9 2.0 0.7 0.8 0.7 2.6 4.7 2.3 1.2 3.6
5.2 0.9 1.6 1.6 2.4
4.9 3.4 3.6 2.7 4.5
8 20 8 16 4 4 6.4 4.7 6.8 6.4 1.0 6.3 3.2 3.0 2.1 2.5
6.1 4.2 2.1 2.5 3.2 2.9 1.2 2.5 1.6 2.1
8.3 6.5 5.4 3.5 4.2
3.9 4.4 11.9 3.3 3.6
2.1 2.3 5.7 2.8 3.7
81
Q2
83
84
Qs
40 000 250 000
20 000 150000
18OOOO
160OOO
120000
Table X. Default Initial Vector, Global Optimum, and Results for Example 8 Nj v, Bi default 3 2500 334.722 3 2500 217.222 3 2500 solution 1 2 250 120 2 360 60 1 480 solution 2 1 407.273 140 2 610.909 101.818 1 560 solution 3 2 373.333 186.667 1 560 93.333 1 746.667 global optimum (Grossmann and Sargent, 1979) 1 480 240 1 720 120 1 960
TL*
F
1.2 3.2 6.2 3.4 2.2
occurrences
6.66667 10.6667 246 007.0 10 8 40977.5
95
43 813.3
4
41 844.4
1
38 499.8
0
10 16 20 8 20 16
typical CPU time per run = 1 min 36 s
Table XI. Results from Phase I of the Relaxation Strategy for Example 8 seed run N violation on I violation on 2 377 003 613 1 (1JJI 0 1204 2 {1,1,1l 0 628.4 3 11,1,1l 0 336.0 4 I1JJl 0 139.3 5 I1JJl 0 34.56 6 {1,1,1l 0 9.543 7 !1,1,1l 0 5.168 8 {1,1,1l 0 0.8829 9 12,2911 0 0.064884 10 (2,2,1l 0 0 711072 531 1 UJJI 0 1204 2 {1,1,1l 0 251.8 3 MA 0 34.14 4 I1JA 0 10.66 5 11,1,1l 0 7.692 6 11,1,1l 0 3.396 7 I1ml 0 1.614 8 L2JI 0 0.3443 9 12,2,1l 0 0 1752629 996 1 I1JA 0 1204 2 11,lJI 0 707.4 3 (1,lJI 0 115.3 4 {1,1,1l 0 10.96 5 l1,1,11 0 1.357 6 12,2,1l 0 0
ded in the original MINLP formulation. This corresponds to a fairly large and heavily constrained MINLP problem. The default from the MSGA algorithm corresponds to the midpoint of the search intervals (an infeasible point). The MSGA algorithm was run from this initial solution vector with 100 different seeds, and no feasible solutions were obtained. Kocis and Grossmann (1988) have solved this problem with the OA/ER algorithm after convexifying both the objective function and constraints, using logarithmic transformations. This was also tried with the MSGA algorithm, but again no feasible solution could be obtained. Floudas et al. (1989) have also solved this
Table XII. Results from Phase I of the Relaxation Strategy for Example 9 schedule violation violation seed % run N on I on 1 377003613 1.145 1 ~1.1.2.1.1.1l 10270 2420 2414 1815 512.7 893.1 598.1 0 293.3 0 0 6.339 0 3.282 0.346 10270 2420 8421 1945 5293 1809 2482 1032 1499 619.2 489.1 357.1 123.2 85.54 0 7.647 0 5.053 0 2.405 0 1.760 711 072 531 1.145 3494 997.0 2451 598.0 189.8 537.7 20.38 7.806 0 5.218 0 3.211 1752629 996 1.145 18037 1471 7784 678.7 3924 480.3 2665 251.7 719.6 7.303 0 3.533 0 3.056 0 1.448 0 1.072 typical CPU time per run = 3 min 37 s
problem with the generalized Benders’ decomposition without needing to eliminate nonconvexities. The relaxation strategy was thus applied to the MSGA algorithm,
Ind. Ehg. Chem. Res., Vol. 31, No. 1, 1992 271
Table XIV. Results from Phase I1 of the Relaxation Strategy for Example 9 solution Vi TLi 3000 379.7467 1st set 1891.640 770.3064 1974.683 727.5089 2619.195 638.2978 2328.100 525.4531 2109.797
Bi 3.2 3.4 6.2 3.4 3.7
2nd set
3000 1891,046 1974.705 2617.244 2327.021 2109.145
379.7468 769.7316 727.0072 638.2978 525.1368
3.2 3.4 6.191913 3.396574 3.700610
3rd set
2360 1500 1553.376 2079.653 1822.230 1622.898
29817261 611.6624 559.6196 415.9394 416.3328
3.2 3.4 4.0 1.7 1.9
global optimum (KO& and Grossman, 1988)
3000.0 1892.0 1975.0 2619.0 2328.0 2110.0
380.0 770.0 730.0 638.0 525.0
3.2 3.4 6.2 3.4 3.7
giving the results of Table XII. For the first seed two different reduction schedules were now employed, in order to obtain some information on the relevance of this parameter. Table XI1 shows that feasibility was not obtained with any of the four optimization paths. As expected, more subproblems had to be solved with a smaller reduction schedule. Considering only as candidates those results which showed feasibility on the horizon time constraint
N
F 285 610 308 633
909236 311 700 313 618 316 004 334 823 335 427 339 809 342 194 310 148 333 266 334 831 336 331 338 244 340642 369 449 360058 361 014 366 826 389 943 391 508 315 513 335 631 336 605 338 123 341 918 361 080 362 036 362 463 406 738 385 506
(the only Iinequality besides the bound constraints), Table 2UI suggesta 12 seta of integer variables, where only the set {l,l,l,l,l,2)appsars to be repeated. This situation is clearly different from the previous problems, where feasible solutions could be obtained and where the same global optimum occurred with different optimization paths. A robust method is thus needed to ensure that the correct candidatefor the global optimum is not overhked. The method employed consista of providing different o p
272 Ind. Eng. Chem. Res., Vol. 31, No. 1, 1992
timization paths to the MSGA algorithm,starting from the least infeasible solution obtained. Table XI11 shows the results obtained by running the MSGA algorithm with 100 different seeds for each of the four end points given in Table XII (with a reduction schedule of 0.346% per global iteration for all cases). The starting points in Table XI11 are organized in order of decreasing maximum violation (viz., increased feasibility), and the results in order of decreasing total frequency of occurrence. Apart from the 20 sets shown in Table XIII, 51 distinct sets occurring only once or twice were also obtained. These 71 NLP subproblems were solved with the MSGA algorithm. Although this is a large number, it constitutes only about 1.7% of the total number of subproblems embedded in the original MINLP formulation. Table XIV shows the results of the 31 feasible solutions encountered, together with the global optimum given by Kocis and Grwmann (1988). It is interesting to note that all feasible solutions correspond to only three sets of continuous variables (a behavior similar to the one observed before in example 7). Table XIV also shows that, of the five most highly occurring integer sets (i.e., greater than 10% from any starting point), k., @,2,3,3,2,2),(2,2,3,2,2,21,(2,2,2,2,1,11, (2,2,3,2,1,1},and (2,2,3,2,1,2),three are local optima and one is the global optimum. The optima belonging to the first set are extremelly ill-conditioned. The optima of the second set are also ill-conditioned, but they are also feasible for the first set of continuous variables. The optima of the third set are the most stable, since they are feasible over a wide range of the continuous variables. These can be regarded as “operating optima” in contrast to the first two sets where it would not be reasonable to operate the system. The global optimum determined with the MSGA algorithm agreea well with the global optimum reported by Kocis and Grwmann (1988) except for the parameter TL3.However, the solution truncated as reported by these authors violates five inequality constraints (1)by as much as 0.34%. Conclusions An optimization algorithm (the MSGA algorithm) is proposed, which efficiently solves nonconvex NLP and small to medium scale MINLP problems. The MSGA algorithm is an extension of an adaptive random-search NLP solver (the SGA algorithm). The main characteristics of the NLP solver are that it employs a variable, parameter-dependent compression vector for the contraction of the search regions and that it incorporates shifting strategies allowing for wrong-way moves to be produced. These features coupled with the random-search concept are responsible for the good robustness in escaping local optima. The MINLP solver (the MSGA algorithm)treats integer variables as discrete integer variables and does not require any problem decomposition, neither identification or elimination of nonconvexities, prior to its application. It was tested with several nonconvex MINLP problems, and the results obtained reveal ita adequacy for the optimization of small to medium scale MINLP problems encountered in chemical engineeringpractice. However, due to the random nature of the search space, the fact that the MSGA algorithm showed good convergence to the proposed problems does not imply global proof of convergence. Direct application of the algorithm does not seem to provide any information on the presence of very ill-conditioned optima. A simple two-phase relaxation strategy was deacribed which allows the detection of ill-conditioned optima (global and local) or the solution of MINLP problems whenever feasible solutions cannot be obtained.
The objective of phase I is to obtain successive less infeasible solutions from which possible candidates for the global optimum are identified. Phase I1 solves the corresponding NLP subproblems, and the best feasible result is retained as the global optimum. This relaxation strategy becomes especially relevant for larger scale MINLP problems where the number of degrees of freedom can be very large. The relaxation strategy described is not restricted to the SGA NLP solver and can in principle be applied with any direct random-search method, including random-search clustering algorithms (Price, 1978; Torn, 1978). The results presented in this work show that although it may be possible to solve large-scale MINLP problems with direct random search coupled with a relaxation strategy, the work involved in the identification of the candidate integer sets for the global optimum and the corresponding solution of the various NLP subproblems may be substantial. Thus, for larger problems more efficient (but also less straightforward) methods may be preferable, such as the generalized Benders’ decomposition or the outer approximation/equality relaxation algorithms. Acknowledgment
I thank Dr.S. Fey0 de Azevedo for providing most useful comments and for reviewing this paper. Nomenclature Ai = heat-transfer area of heat exchanger i Bi = batch size (kg) C = input mass flow rate (kg/h) to heat-exchanger network F = objective function Fi= mass flow rate of stream i (kg/h) in heat-exchanger network H = horizon time (h) M = number of processing stages in series N = number of products N j = number of parallel units for stage j (j = 1,M) Qi= required production (kg) of product i ( i = 1, N) gj = inequality constraint (j = 1, m) m = number of inequality constraints n = number of independent parameters (continuous and discrete) P1= number of random trials per interval reduction (set at 100) Pz = number of interval reductions (1400) p = number of independent continuous parameters r = search interval vector Sij = size factor (l/kg) of product i in stage j (i = 1, N; j = 1,M) T = absolute temperature TLi= cycle time (h) of product i (i = 1, N) tij = processing time (h) of product i in stage j (i = 1,N; j = 1, M) U = square {(n- p ) X (n - p ) ) matrix of random numbers between 0 and 1 or overall heat-transfer coefficient matrix Vi = size (1)of parallel unit in processing stage j (j = 1, M) x = continuous variable vector x* = continuous variable current optimum vector y = discrete (integer) variable vector y* = discrete (integer) variable current optimum vector Z = square (p x p) matrix of random numbers between -l/2 and +‘/z Greek Symbols a = cost parameter of capital cost of heat exchanger ai = lower bound on continuous parameter xi (i = 1, p) aj = cost coefficient ($) of stage j (j = 1, M) B = cost parameter of capital cost of heat exchanger Pi = upper bound on continuous parameter xi (i = 1,p)
Ind. Eng. Chem. Res., Vol. 31, No. 1, 1992 273
flj = cost coefficient (dimensionless) of stage j 0’ = 1, M) y = lower bound vector for integer parameters yi = lower bound on integer parameter yi ( i = p + 1, n) f i = upper bound on integer parameter yi (i = p + 1, n) ei = compression factor for continuous variable x i Abbreviations
BB = branch and bound method LMTD = log mean temperature difference IMILP = mixed-integer linear programming MINLP = mixed-integer nonlinear programming NLP = nonlinear programming SGA = Salcedo-Gongalves-Azevedo algorithm MSGA = MINLP Salcedo-Gonplves-Azevedo algorithm OA = outer approximation algorithm OA/ER = outer approximation with equality relaxation algorithm Literature Cited Aarta, E.; Korst, J. Simulated Annealing and Boltzmann Machines-A Stochastic Approach to Combinatorial Optimization and Neural Computers; Wiley: New York, 1989;Chapters 3-5. CerdH, J.; Vicente, M.; GutiBrrez, J.; Eaplugas, S.; Mata, J. Optimal Production Strategy and Design of Multiproduct Batch Plants. Znd. Eng. Chem. Res. 1990,29,590-600. Ciric, A. R.; Floudas, C. A. A Mixed Integer Nonlinear Programming Model for Retrofitting Heat-Exchanger Networks. Ind. Eng. Chem. Res. 1990,29,239-251. D o h , W. B.; Cummings, P. T.;Levan, M. D. Process Optimization via Simulated Annealing: Application to Network Design. AIChE J. 1989,35 (5),725-736. Duran, M. A.; Grossmann, I. E. A Mixed Integer Nonlinear Programming Approach for Process Systems Synthesis. AZChE J. 1986,32(4),592-606. Edgar, T.F.; Himmelblau, D. M. Optimization of Chemical Processes; McGraw-Hilk New York, 1988; pp 422-423. Fey0 de Azevedo, S.; Salcedo, R.; Gonqalves, M. J. Optimization Strategim in Chemical Engineering. Recent Prog. Genie Procedes 1988,2 (6),199-204 (Lavoisier, Paris). Floudas, C. A.; Ciric, A. R. Strategies for overcoming uncertainties in heat-exchanger network synthesis. Comput. Chem. Eng. 1989, 13 (10). 1133-1152. Floudaa, C. A.; Aggarwal, A.; Ciric, A. R. Global optimum search for nonconvex NLP and MINLP problems. Comput. Chem. Eng. 1989,13(lo), 1117-1132. Gaines, L. D.; Gaddy, J. L. Process optimization by flow sheet optimization. Znd. Eng. Chem. Process Des. Dev. 1976, 15, 206. Geoffrion, A. M. Generalized Benders’ Decomposition. J. Optim. Theory Appl. 1972,10,237. Grosemann, I. E.; Sargent, R. W. H. Optimal design of multipurpose chelqical plants. Ind. Eng. Chem. Process Des. Dev. 1979,18, 343-348. Heuckroth, M. W.; Gaddy, J. L.; Gaines, L. D. An examination of the adaptive random search technique. AIChE J. 1976,22 (41,744. Kocis, G. R;Grosemam, I. E. Relaxation Strategy for the Structural Optimization of Process Flow Sheets. Znd. Eng. Chem. Res. 1987, 26,1869-1880. Kocis, G. R.; Grossmann, I. E. Global Optimization of Nonconvex Mixed-Integer Nonlinear Programming (MINLP) Problems in Process Synthesis. Znd. Eng. Chem. Res. 1988,27, 1407-1421.
Kokossis, A. C.; Floudas, C. A. Optimization of complex reactor networks. I. Isothermal Operation. Chem. Eng. Sci. 1990,5 (3), 595-614. Lin, S.; Kernighan, B. W. An effective heuristic algorithm for the traveling salesman problem. Oper. Res. 1973,21 (2),498-516. Luus, R. R. Optimization of multistage recycle systems by direct search. Can. J. Chem. Eng. 1975,53,217. Luus, R. R.; Jaakola, T.H. I. A direct approach to the optimization of a complex system. AZChE J. 1973a,19 (3),645. Luus, R. R.; Jaakola, T. H. I. Optimization by direct search and systematic reduction of the size of search region. AIChE J. 1973b, 19 (4),760. Luus, R. R;Brenek, P. Incorporation of gradient into random search optimization. Chem. Eng. Technol. 1989,12,309-318. Martin, D. L.;Gaddy, J. L. Process optimization with the adaptive randomly directed search. AIChE Symp. Ser. 1982,78 (214),99. Ostrovsky, G.M.; Ostrovsky, M. G.; Mikhailow, G. W. Discrete optimization of chemical processes. Comput. Chem. Eng. J. 1990, 14 (l), 111-117. Paterson, W. R. A replacement for the logarithmic mean. Chem. Eng. Sci. 1984,39,1635-1636. Press, W. H.; Flannery, B. P.; Teukolsky, S. A.; Vetterling, W. T. Numerical Recipes: The Art of Scientific Computing; Cambridge University Press: New York, 1986;pp 326-334. Price, W. L. A Controlled Random Search Procedure for Global Optimization. In Towards Global Optimization 2; Dixon, L. C. W., Szego, G. P., Eds.; North-Holland Amsterdam, 1978; pp 71-84. Ray, W. H.; Szekely, J. Process Optimization-With Applications in Metallurgy and Chemical Engineering; Wiley: New York, 1973; pp 151-152. Reiter, S.; Sherman, G. Discrete Optimizing. J. SOC.Znd. Appl. Math. 1965,13,214. Salcedo, R.; GonMves, M. J.; Fey0 de Azevedo, S. An improved random-search algorithm for nonlinear optimization. Comput. Chem. Eng. 1990,14(lo), 1111-1126. Shedler, G. S. Generation methods for discrete event simulation. Computer Performance Modeling Handbook; Academic Press: New York, 1983; pp 227-251. Siddall, J. N.Optimal Engineering Design-Principles and Applications; Dekker: New York, 1982;p 223. Soland, R. M. An algorithm for separable nonconvex programming problem II: nonconvex constraints. Manage. Sci. 1971,17,759. Stanley, R.; Sherman, G. Discrete optimizing. J. SOC.Znd. Appl. Math. 1965,13,864-889. Stephanopoulce,G.; Westerberg,A. W. The Use of Hestenes’ method of multipliers to resolve dual gaps in engineering system optimization. J. Optim. Theory Appl. 1975,15,285. TBrn, A. A. A Search-Clustering Approach to Global Optimization. In Towards Global Optimization 2; Dixon, L. C. W., Szego, G. P., Eds.; North-Holland Amsterdam, 1978; pp 49-62. Wang, B. C.; Luus, R. R. Optimization of nonunimodal systems. Znt. J. Numer. Methods Eng. 1977,11,1235. Wang, B. C.; Luus, R. R. Reliability of optimization procedures for obtaining global optimum. AZChE J. 1978,24 (4),619. Yuan, X.;Zhang, S.; Pibouleau, L.; Domenech, S. Une methode #optimization nonlineaire en variables mixtea pour la conception de procedes. RAZRO-Oper. Res. 1989.
Received for review January 4,1991 Revised manuscript received September 11,1991 Accepted September 26,1991
