4296
Ind. Eng. Chem. Res. 2000, 39, 4296-4301
Global Optimization of Modular Process Flowsheets R. P. Byrne† and I. D. L. Bogle* The Centre for Process Systems Engineering, Department of Chemical Engineering, University College London, Torrington Place, London WC1E 7JE, U.K.
This paper presents an approach for applying rigorous global optimization techniques to modular process design approaches. Modular flowsheeting systems are very popular in process engineering and some use local optimization methods. However, even very simple process engineering problems can give rise to nonlinear, nonconvex optimizations with multiple local optima. In these circumstances, local optimization approaches cannot guarantee that the global optima will be found. The paper introduces an approach to process simulation that allows construction of flowsheets in modular approaches that can then be optimized using interval global optimization methods. The modular flowsheets are constructed with generic unit modules that can provide interval bounds, linear bounds, derivatives, and derivative bounds using extended arithmetic types. Extended types are introduced to be used by the generic unit modules. Optimization is achieved by application of global optimization algorithms to the modular flowsheets built from these general models. Using interval analysis and automatic differentiation as the arithmetic types, lower bounding information is used in a branch and bound framework. 1. Introduction Flowsheet optimization can provide significantly better designs at modest cost during the early stages of design. Modular flowsheeting systems dominate the market and are used by much of the design community. Optimization is possible with such systems but, to date, only through local optimization methods. Many process models are nonconvex, so that local optimization approaches cannot guarantee global optimality. To do this, a rigorous global optimization approach must be adopted. Reviews of global optimization methods have been provided by To¨rn and Zˇ ilinskas1 and Horst and Tuy.2 Interest in global optimization has been mirrored in the chemical engineering literature. Floudas3 gave a recent review of the use of the development of global optimization techniques for process engineering problems, and many contributions have been gathered in Grossmann.4 Floudas and Visweswaran5 presented a “primal-dual” decomposition approach based on the generalized benders decomposition. Vaidyanathan and El-Halwagi6 and Byrne and Bogle7 have proposed modifications to the interval method.8,9 Ryoo and Sahinidis10 used convex underestimators to solve a large number of engineering problems. Adjiman et al.11 present a global optimization method for the solution of general chemical process design problems. Quesada and Grossmann developed specialized convex underestimators for the solution of a class of heat exchanger network problems and bilinear mass transport problems.12 All of this previous work has applied global optimization techniques to problems in which the equations are explicit. In this paper, we will address the problem of applying global optimization to modular approaches to flowsheeting. * Author to whom correspondence should be addressed. Tel.: +44 20 7679 3803. Fax: +44 20 7383 2348. E-mail:
[email protected]. † Current Address: Salomon Smith Barney, London SW1, U.K.
There are many ways of formulating flowsheeting problems. Two important distinctions are identified in the literature. In the equation-oriented (EO) formulation, the flowsheet is treated as a set of mass/energy balance equations that are solved simultaneously, whereas the sequential modular (SM) approach views the flowsheet as a set of interconnected black-box models of the unit operations. Comprehensive reviews of flowsheeting techniques are given by Biegler13 and Perkins.14 A recent thorough comparison is provided by Biegler et al.15 Both approaches have their advantages, and consequently, a large body of research has been devoted to the construction of flowsheeting methodologies with the best mixture of properties. This is particularly the case for flowsheet optimization where the sequential component of SM flowsheeting can result in poor performance as the flowsheet is converged at every optimization step. This has led to the simultaneous modular approach that has been used here. Global optimization methods based on interval analysis are discussed in section 2 with an explanation of how they could be implemented with modular approaches. Two simple examples, one without recycle and one with, that demonstrate the approach are presented in section 3 with some numerical results. Section 4 presents some conclusions about the effectiveness of the approach and further developments that may be required to implement this approach on existing systems. 2. A Module Type for Global Optimization To globally optimize a modular simulation, the flowsheet model must provide sufficient information to the global optimization algorithm for it to guarantee that the global minimum is found. For global optimization using interval techniques it is necessary to obtain upper and lower bounds on the problem, not just use the value of the current point. Figure 1 illustrates the problem of obtaining the necessary bounds from black-box models. The full range for y is not given by evaluating y at x1
10.1021/ie990619d CCC: $19.00 © 2000 American Chemical Society Published on Web 10/14/2000
Ind. Eng. Chem. Res., Vol. 39, No. 11, 2000 4297
Figure 2. Generic model.
Figure 1. Black-box models.
and x2. If it is only possible to evaluate the output at fixed points, no information about the behavior of the model between those points is provided. It is not possible to minimize the model rigorously with only end-point information if it is nonconvex. What is needed is a computational representation of a modular flowsheet that can provide the optimization algorithm with the information necessary to guarantee a global solution. To achieve this end, we need to describe the extended arithmetic types used and how they are used with a general module that can deal with interval information. 2.1. Interval Analysis and Extended Arithmetic Types. Interval analysis has been used, to varying degrees, by a number of global optimization approaches, in particular by Vaidyanathan and El-Halwagi,6 Adjiman et al.,11 and Byrne and Bogle.7 An interval is a compound type made from two real numbers, the upper and lower bounds, and is used to represent a range. An extended type is made from a compound type plus a set of rules that define operations on it. Thus, the interval compound, used according to the rules of interval arithmetic,8 is an extended type called an interval. The extended arithmetic type makes it possible to perform calculations using a range of arithmetic operators to obtain information about the calculation. Each type, T, corresponds to performance of the calculation on a different space, T. For example, floating point arithmetic approximates T ) R; interval types perform calculations on T ) I, where I is the set of all possible intervals and so on. Examples of extended types that are useful to optimization are the following: (i) Intervals are an extended type for calculating the results of application of operators to ranges. (ii) Vectors and matrices are extended types built from arrays of scalar types and rules for manipulating the data. (iii) The convex underestimators proposed by McCormick16 provide the rules for addition, subtraction, and multiplication of convex underestimators. A suitable definition of the compound object that maintains the data required at each step makes this an extended type. (iv) The automatic differentiation rules proposed by Rall17 are, when combined with a compound type for the data, an extended type. Of these, intervals and automatic differentiation types are going to be used in the approach presented here. An interval is an ordered pair of real numbers X ) [x,xj], where [x] is the lower bound on x and [xj] the upper bound. It is the set of real numbers {x ∈ R|x e x e xj}. The set of all possible intervals is denoted by I. A series of operations has been defined over sets of intervals8 that are used in calculations involving intervals in this work. A great many interval operations can be defined for general functions. These are called inclusion functions.9,18 Automatic differentiation (AD) is a method for simultaneously computing partial derivatives with a func-
tion’s value by replacing the numbers in a function with a multicomponent object, called a gradient, whose algebraic properties incorporate the chain rule of differentiation. The rules, together with the gradient type, form an extended type that we will call AD. The common arithmetic operations are determined using the chain rule. An advantage of AD is that it can be applied to factorable functions. That is, the function F can be written as a set of procedures, as would be the case in modular flowsheeting systems, and not as a single equation. This approach relies on the fact that every function, no matter how complicated, is executed on a computer as a sequence of elementary operations and elementary functions. By repeated application of the chain rule to the composition of these elementary operations, it is possible to compute derivatives that are accurate to machine precision. The procedures involved may contain arbitrary branches, conditional loops, and subprocedures.17 Griewank et al.19 developed procedures for automatic differentiation of procedures. Unit modules that require iteration can be solved using the interval Newton method proposed by Schnepper and Stadtherr.20 2.2. A Generic Model. For the current generation of modular flowsheeting systems, a module requires values for all of the input streams and unit parameters and calculates the condition of the output streams. Systems use the modules as black-box modules. However, it is possible to create a general module in which more information is passed and retrieved from the module. This more generic module is needed to be able to use interval methods for global optimization where not only single data points but also intervals are required. Rather than only taking one type of scalar input for stream or parameter variables, a generic model is a model that specifies the transformations that need to be applied to some more general variable type, T, to obtain the output, as shown in Figure 2. The model contains the necessary operations, and then the appropriate rules are applied for the type T. For example, a module that adds its inputs should use the interval arithmetic rules if the underlying type is an interval (T ) I) and the rules of automatic differentiation if the underlying type is an AD type as defined above. This is an extension of a common idea. A vector arithmetic is a generic way of applying operations to vectors of some type, T. When two real vectors are added, the elements are added componentwise; when two interval vectors are added, the intervals are added componentwise. The addition of the components differs, but the rule for adding vectors is the same. A form of the generic module that can be applied to any unit operation is shown in Figure 3. Each unit j has n input streams that are the outputs of the preceding units and m output streams that are the feeds for successive units. There are p input parameters, φ, that determine how the unit operates, for example, the split
4298
Ind. Eng. Chem. Res., Vol. 39, No. 11, 2000
Figure 3. Generic unit module.
fraction in a splitter and a cost, γ, that represents some cost associated with the unit. There are r residuals, F, which represent constraints. Once the types are constructed, they can be used with the model to provide the appropriate information. Because the generic model specifies which operations (or transformations) are to be applied and the type defines how those operations are applied to the compound, a generic model can be used with types that are unknown when the model is constructed. Two simple interval unit operations are shown in Figure 4. The equations are defined in terms of extended types. Unit modules consist of a sequence of internal computational operations. For this approach the operations must be cast in terms of the extended type operations. This requires implementation of new operations in the chosen computational language. In this work, the implementation has been done in Matlab. For the many existing modular systems, this will require implementation of such new operations (for example, in Fortran). The extent of further coding modifications required to adapt existing systems is still under investigation. For simple models without iteration existing in some language such as Fortran, these modifications are expected to be modest. Modules with iteration will be solved using the interval Newton method mentioned above.20 For more complex modules and for physical property calculations, further investigation is required. 2.3. Optimization of Extended-Type Flowsheets. Flowsheets can be built from generic units and then evaluated using operations appropriate for some extended type, T, that provides the information necessary for the global optimization algorithm. Consider an algorithm that only uses interval bounding information. The variable type, T, being used is interval. The algorithm is provided with initial bounds on the unit parameters which become the initial bound vector X0. At each iteration k, interval parameters, Xk, are chosen by the algorithm. Each unit calculates the output streams and the cost associated with the unit. The summation of these costs provides the objective function. This overall scheme is shown in Figure 5. The flow of the process streams is shown in solid lines, and the flow of information through the flowsheet is indicated by dashed lines. The design constraints on the products are added to the optimization problem, and the flowsheet is optimized to minimize the sum of the costs subject to these constraints. Interval gradient types have been used, and the global optimization algorithm uses the MV underestimation scheme21 to construct a linear relaxation of the objective
function and constraints in terms of the optimization variables. Solution of this linear relaxation provides the necessary bounds on the optimization problem. The AD type has been used in order that the flowsheet can be solved using an algorithm that uses bounds on the derivatives to construct lower bounding LPs.21 Evaluation of the flowsheet with interval gradients provides all of the information, bounds and derivative bounds, that is required to construct the relaxation. The benefits of using this relaxation are described in ref 21. In the examples used here, the cost of using the interval algorithm without derivative bounding information was excessive. This approach works for other variable types and algorithms in the same way. For example, a flowsheet based on convex underestimators can be optimized by one of the many algorithms based around convex underestimation.10,16,22 Kearfott23 has compared a number of these algorithms on a standard problem set. The optimization algorithm provides the generic unit operation models with information of the appropriate type, which will frequently be a vector of interval parameters, and the information required for the next iteration is calculated through the flowsheet and returned to the algorithm, as shown in Figure 5. 2.4. Solving the Recycle Problem. A technique for dealing with recycles is vital for a modular flowsheeting methodology. To maintain a view consistent with the modular flowsheeting systems, a stream can be torn by introducing a “tear unit” that has an input, an output, and a residual based on equality of the two streams. The tear unit parameters are the output stream values, as shown in Figure 6. Because the units, including the tear unit, may have residuals, this approach allows the problem to be formulated in a number of ways ranging from sequential to simultaneous. Some alternatives are considered next. If the recycle is torn with a tear block, then the range of inputs to the flowsheet, X1, is selected by the interval optimization, and the flowsheet is calculated through the interval units to produce a set of outputs, X2, which bound all of the possible outputs from the flowsheet with parameters determined by X1. When an interval type is used, it is possible to iterate around the recycle, as in a traditional sequential modular system. Given a “guess” for the values in the torn stream, XT, the resulting calculated torn stream, XC, will also be a range that encloses all of the possible values of the recycle stream with XT as an input. This process might be accelerated by application of the interval Newton algorithm. Typically, this involves reducing the range of XT such that it still contains all solutions rather than reducing it to a point. Adopting the sequential modular flowsheet strategy directly, which solves the recycle at each iteration, is difficult with a flowsheet that uses extended-type parameters such as intervals because there may be more than one solution to the recycle problem for a given set of inputs and, to globally optimize the flowsheet, we must not lose any solutions. The simultaneous modular approach corresponds to the use of a recycle module in which the parameters are independent variables from the optimization and the residuals are the residuals of the constraint XT - XC ) 0. In this case, the recycle is converged as the optimization converges, reducing the amount of computation that must be performed per iteration.
Ind. Eng. Chem. Res., Vol. 39, No. 11, 2000 4299
Figure 4. Unit operations (mixer and splitter) based on interval arithmetic.
Figure 5. Optimization of an interval flowsheet.
Figure 6. Tear unit module.
The residuals are of some type, T, and the optimization algorithm is chosen to handle this arithmetic type. This approach can be used with any optimization algorithm defined for type T. 3. Examples Two example problems will be presented. The first is an acyclic problem, the Haverly pooling problem. The second uses a reactor with implicit equations and a recycle. Both problems can be formulated as explicit equations which helps to show how the implementation can be done and to make comparisons easier. The problems are well studied and are known to have multiple optima that have all been established. They are also easily cast as modular problems.
3.1. An Acyclic Problem. Consider the Haverly pooling problem given in Figure 7 as formulated by Quesada and Grossmann.22 The problem is to blend four feeds into two products to minimize the total cost. Each feed has a fixed composition, xA, and cost, c [£/(kmol h)]. Each product has a cost, required composition, and flowrate, F. The feeds are mixed to produce products that satisfy the quality requirements using mixer and splitter units to represent the different blending tanks. The problem is solved here as a modular problem using simple mixer and splitter modules to demonstrate the methodology. An equation-based formulation is also given in Quesada and Grossmann.22 This problem is nonconvex because of the bilinear terms involved in the splitter unit. The problem has an infinite number of local minimizers when the feed flowrates are zero, a strong local minimizer with an objective function value of -100, and a global minimizer at
F ) [0, 100, 0, 0, 0, 100, 100, 200]T and xA ) 0.01 with f(x) ) -400. The problem can be formulated in modular form. All of the mixers are the same, as are all of the splitters. The feeds are unit operations without inputs where the flowrate is a parameter of the unit. The product units represent quality requirements. This is one of many pooling/blending problems that may need to be solved, and it would be advantageous
4300
Ind. Eng. Chem. Res., Vol. 39, No. 11, 2000
Figure 7. Haverly pooling problem.
The reactor model has a single input stream, F2,i, and three parameters, A1, B1, and the extent of reaction, V.
F3,A ) F2,A + V(-A1) F3,B ) F2,B + V(A1 - B1) F3,C ) F2,C + V(A1 - B1) The following constraints are used as residuals:
A1
B1
Figure 8. Problem with recycle.
to be able to formulate the problem as a modular problem: placing units onto the flowsheet and “drawing” the streams to connect them. It is natural to view the problem in terms of interconnected units with each unit performing some transformation of the input stream to provide the output stream. The interval algorithm is applied with the mean value inclusion that utilizes interval gradient types. More details on the optimization algorithms used can be found in Byrne and Bogle.21 The problem was solved in 60 s using an implementation on an Intel 486 with 8M of RAM. Using the reduced interval global relaxation algorithm21 implemented using Matlab, the problem was solved in just under 20 s. 3.2. Recycle Flowsheet. This case study, shown in Figure 8, is taken from a synthesis problem.24 The operating parameters and sizes for one of the synthesis alternatives are optimized using the detailed models and costing information provided.24 Each unit has a capital cost, Cc, and an operating cost, Co, which are incorporated into an objective function through a payout time of 2.5 years. The principle units are a CSTR and two separation columns. The models have been formulated in terms of component flowrates, Fs,i. 3.2.1. Reactor Model. The reactor is a continuous stirred tank reactor which models the reaction between chlorine and benzene (A) to produce monochlorobenzene (B) and dichlorobenzene (C) at constant temperature.
∑i F3,i ) 0.413F3,A ∑i F3,i ) 0.055F3,B
and the capital cost is given by
Cc ) 25 795 + 8178VCo ) 0 3.2.2. Columns. The columns perform sharp splits. The costs are
C1
Cc ) 132 718 + 369F3,A - 1114F3,B Co ) (21.67 + 4.65)39.1F3,A + 10.7F3,B + 3F3,C
C2
Cc ) 25 000 + 6985F4,B - 1114F3,B Co ) (21.67 + 4.65)(55.662F4,B + 26.212F4,C)
The costs of the feeds and products are the purchase price of benzene, $27.98/kmol, and chlorine, $19.88/ kmol, and the sale price of monochlorobenzene, $92.67/ kmol, for which there is a demand of 50 kmol/h. These equations are implemented in Matlab as interval equations with interval operations defined to produce interval outputs. Using the reduced interval global relaxation algorithm, the problem was solved in 299.2 s using a Matlab implementation on the machine described above. 4. Conclusions This paper has shown that, by exploiting extended arithmetic types, in particular those using interval
Ind. Eng. Chem. Res., Vol. 39, No. 11, 2000 4301
types, it is possible to rigorously globally optimize a flowsheet using the simultaneous modular approach. To make this possible, it was necessary to define a generic module type that can handle extended arithmetic types. Such extended types can be implemented using most modern programming languages. The generic modules could be used in the flowsheet optimization. This has the advantage that a single code for each model can be used to generate floating point results, interval results, gradient information and gradient bounding information. The simultaneous modular approach has been used to optimize the two examples. Sequential solution, including the solution of implicitly defined units, is also being investigated. Achievement of the global solution will entail considerably longer times than using local optimization methods, so comparisons between the two are not meaningful. The purpose here is to show that such an approach is possible and that there is a way of formulating the modular approach, with its many advantages, to obtain the globally optimal solution. As the size of the problem increases the complexity of the computational problem will increase at least exponentially. This will make the solution of large problems particularly difficult. The results demonstrate that the approach can be used to globally optimize flowsheets using the modular approach. The extent of modification that might be required for legacy code in existing systems, in particular for physical property prediction subroutines, is currently under investigation. Interval implementations can deal with all expected equations and conditional expressions, and the issue of internal iterations can be solved using the interval Newton algorithm. Acknowledgment The authors acknowledge support from the Engineering and Physical Sciences Research Council of the U.K. and the Centre for Process Systems Engineering. Literature Cited (1) To¨rn, A.; Zˇ ilinskas, A. Global Optimization; Lecture Notes in Computer Science. Springer-Verlag: Berlin, 1989. (2) Horst, R.; Tuy, H. Global Optimization. Deterministic Approaches; Springer-Verlag: Berlin, 1992. (3) Floudas, C. A. Recent Advances in Global Optimization for Process Synthesis, Design and Control: Enclosure of All Solutions. Comput. Chem. Eng. 1999, 23, Supplement S963-S973. (4) Grossmann, I. E. Global Optimization for Engineering Design; Kluwer Academic Publishers: Dordrecht, The Netherlands, 1996. (5) Floudas, C. A.; Visweswaran, V. A Global Optimization Algorithm (GOP) for Certain Classes of Nonconvex NLPs. 1. Theory. Comput. Chem. Eng. 1990, 14 (12), 1397.
(6) Vaidyanathan, R.; El-Halwagi, M. Global Optimization of Nonconvex Programs via Interval Analysis. Comput. Chem. Eng. 1994, 18, 889. (7) Byrne, R. P.; Bogle, I. D. L. Solving Nonconvex Process Optimization Problems Using Interval Subdivision Algorithms. In Global Optimization for Engineering Design; Grossmann I. E., Ed.; Kluwer Academic Publishers: Dordrecht, the Netherlands, 1996. (8) Hansen, E. Global Optimization Using Interval Analysis; Marcel Dekker: New York, 1992. (9) Moore, R. E. Interval Analysis; Prentice-Hall: Englewood Cliffs, NJ, 1966. (10) Ryoo, H. S.; Sahinidis, N. V. Global Optimization of Nonconvex NLPs and MINLPs with Applications in Process Design. Comput. Chem. Eng. 1995, 19 (5), 551. (11) Adjiman, C. S.; Androulakis, I. P.; Maranas, C. D.; Floudas, C. A. A Global Optimization Method, RBB, for Process Design. Comput. Chem. Eng. 1996, 20S, S419. (12) Quesada, I.; Grossmann, I. E. Global Optimization Algorithm for Heat-Exchanger Networks. Ind. Eng. Chem. Res. 1993, 32 (3), 487. (13) Biegler, L. T. Simultaneous Modular Simulation and Optimization. In Proceedings of the Second International Conference on Foundations of Computer-Aided Process Design; CACHE, 1983. (14) Perkins, J. D. Equation Oriented Flowsheeting. In Proceedings of the Second International Conference on Foundations of Computer-Aided Process Design; CACHE Publications: University of Michigan, Ann Arbor, MI, 1983. (15) Biegler, L. T.; Grossmann, I. E.; Westerberg, A. W. Systematic Methods of Chemical Process Design; Prentice Hall: Englewood Cliffs, NJ, 1997. (16) McCormick, G. P. Computability of Global Solutions to Factorable Nonconvex Programs. Part 1. Convex Underestimating Problems. Math. Prog. 1976, 10, 147. (17) Rall, L. B. Automatic Differentiation: Techniques and Applications; Lecture Notes in Computer Science; SpringerVerlag: Berlin, 1981. (18) Ratschek, H.; Rockne, J. New Computer Methods for Global Optimization; Ellis Horwood Ltd.: Chichester, U.K., 1988. (19) Griewank, A.; Juedes, D.; Utke, J. Algorithm 755; ADOLC, a package for the automatic differentiation of algorithms written in C/C++. ACM Trans. Math. Software 1996, 22 (2), 131. (20) Schnepper, C. A.; Stadtherr, M. A. Robust Process Simulation Using Interval Methods. Comput. Chem. Eng. 1996, 20, 187. (21) Byrne, R. P.; Bogle, I. D. L. Global optimization of constrained nonconvex programs using reformulation and interval analysis. Comput. Chem. Eng. 1999, 23, 1341. (22) Quesada, I.; Grossmann, I. E. Global Optimization of Bilinear Process Networks with Multicomponent Flows. Comput. Chem. Eng. 1995, 19 (12), 1219. (23) Kearfott, R. B. Test Results for an Interval Branch & Bound Algorithm for Equality-Constrained Optimization. In State of the Art in Global Optimization; Floudas, C. A., Pardalos, P. M., Eds.; Kluwer Academic Publishers: Dordrecht, The Netherlands, 1996. (24) Floudas, C. A. Nonlinear and Mixed-Integer Optimization. Fundamentals and Applications; Oxford University Press: New York, 1995.
Received for review August 16, 1999 Revised manuscript received June 6, 2000 Accepted August 24, 2000 IE990619D
