Nonlinear Homotopic Continuation Methods: A Chemical Engineering

Sep 19, 2013 - ... Cited-by Linking service. For a more comprehensive list of citations to this article, users are encouraged to perform a search inSc...
0 downloads 0 Views 773KB Size
Review pubs.acs.org/IECR

Nonlinear Homotopic Continuation Methods: A Chemical Engineering Perspective Review Hugo Jiménez-Islas, Gloria M. Martínez-González, José L. Navarrete-Bolaños, José E. Botello-Á lvarez, and J. Manuel Oliveros-Muñoz* Departamento de Ingeniería Química-Bioquímica, Instituto Tecnológico de Celaya Avenida Tecnológico y Antonio García Cubas s/n, Celaya, Gto. C.P. 38010, México S Supporting Information *

ABSTRACT: It is difficult to imagine mathematics or engineering to date without nonlinear algebraic systems; these topics have deserved more diverse and vigorous interest over the years. In fact, one of the most important modern branches of mathematics originated in this context, topology, which has been inseparable from homotopy theory, describes a powerful intuitive tool to solve nonlinear algebraic systems. In the present work, we attempt to briefly summarize the theoretical beginnings of homotopy, its eventual transformation into a numerical method, its most reiterative applications in chemical engineering, and some current challenges in this field from an engineering-friendly perspective, especially approachable for chemical engineers.

1. INTRODUCTION Algebraic systems are one of the most powerful resources in mathematics; these tools represent the abstraction of several simultaneous truths. Philosophically, this representation is an analytic perspective of any real phenomenon that considers only practical aspects to solve algebraic systems while neglecting rigorous mathematical definitions. One of these considerations is a root or solution vector, which is defined by the calculation of the set or sets of values for the variables at which the equalities of the systems become true. By this definition, the calculation of a zero corresponds to the system solution when each of the system equations equals zero. Furthermore, nonlinear algebraic systems have additional drawbacks when compared to linear systems that are attributable to multiple solutions. In other words, more than one combination of values for the nonlinear system variables may satisfy the problem in many cases. This phenomenon is known as “multiplicity of solutions” and has been controversial in many fields of applied science because many systems with multiple solutions have only one thermodynamically feasible solution vector, and the rest are only mathematical solutions. However, some engineering problems require that all solutions are known to choose the most appropriate solution. Phase separation problems, the modeling of real gases, chemical equilibrium calculations, heat and mass transfer problems, parameter estimation, the control of microbiological processes, and predicting the three-dimensional structure of macromolecules constitute just a few examples of such problems in chemical and biochemical engineering. Furthermore, beyond being a mathematical curiosity, the multiplicity of solutions in algebraic systems has been experimentally demonstrated in a variety of physical phenomena.1−6 Thus, the multiplicity of solutions will always have practical importance from an engineering analysis perspective. Furthermore, determining the solution vectors for nonlinear algebraic systems harkens back to the solution methods themselves. The more sophisticated and robust methods to © 2013 American Chemical Society

solve this class of problems are iterative methods, which are numerical strategies to modify an initial value for the variables of the system, known as an initial guess, until a certain tolerance value is reached.7,8 Deterministic solution methods for nonlinear algebraic systems can be grouped into two general categories: the local methods, which only converge from numerically close initial guesses to any solution vector, and the global methods, which converge from any initial guess irrespective of its proximity to the solution vectors.7−9 In the present work, the so-called global methods are of particular interest. Specifically, this work focuses on “numerical homotopic methods” because they not only converge from any initial guess but also theoretically allow this feature from all locations of the solution vectors of any system of algebraic equations, regardless of how intricate its topology is.8−10 We present an intuitive straightforward historical analysis of homotopic methods with a chemical engineering perspective from its theoretical origins to the present date. This paper is organized as follows. In Section 2, we provide a brief review of the homotopy theory and homotopic continuation algorithms. Section 2.1 presents one geometrical analogy instead of rigorous mathematical formalisms; this facilitates the comprehension of the connection among some important abstract concepts in the homotopy theory and almost the numerical methods described here. Section 2.2 provides the generic sequence of homotopic continuation methods as a basis for supporting the development of the subsequent sections. A historical review of continuation algorithms development, in the chemical engineering literature mainly, is presented in Sections 2.3−2.7. In Section 2.8, we list some of the most cited homotopy forms and their particular Received: Revised: Accepted: Published: 14729

July 29, 2013 September 12, 2013 September 19, 2013 September 19, 2013 dx.doi.org/10.1021/ie402418e | Ind. Eng. Chem. Res. 2013, 52, 14729−14742

Industrial & Engineering Chemistry Research

Review

would represent to the simpler algebraic system easy to solve (at the t value of 0). All the stages of the homotopic deformation process (the application of the function H on the meta-domain) has its geometrical analogy with the intuitive process of compression. 2.2. The Numerical Method. After the contributions of Poincaré (1854−1912), other authors extended the concept of his primitive algebraic topology to include other vector spaces that did not necessarily have a geometrical representation.12,13 Algebraic equations and systems of algebraic equations were treated with this approach to generate several formulations. Thus, it was deduced that any system of algebraic equations can be deformed into another system of the same dimensions under certain conditions.14 Finally, the mathematical theory was applied to a practical algorithm first proposed by Lahaye,15,16 refined by Ficken,17 and adapted to the most known approach by Davidenko18 and Haselgrove.19 The basic ideas of this algorithm are presented below: (1) To solve any algebraic system (F(x)) defined in the field of real numbers, a simpler system (E(x)) can be proposed with only one condition: E(x) must have the same dimensions as F(x). (2) A continuous topological deformation path relates the original complex system and the easier system. The homotopic path is described by

advantages and drawbacks. Section 2.9 presents a group of strategies that bounds homotopic paths. Section 3 provides a historical summary of the algebraic system solving approach in phase separation problems and the hot topic of semianalytical homotopic methods for transport problems. Finally, in Section 4, we conclude the review and mention some perspectives for homotopic methods, especially in chemical engineering. Additionally, we include Appendix A (Supporting Information) with a detailed solution of a two-equation algebraic system; all the nonarbitrary graphics of the paper are referred to in this algebraic system.

2. HOMOTOPY THEORY The theory of homotopy is a current topic of vigorous research in pure mathematics. From an engineering perspective, the basic concepts of this theory are sufficient to understand the philosophy behind even the current methods. 2.1. First Elements. The first notions of homotopy theory are, in fact, the basis of the mathematical branch of topology. The first homotopy theory proponent was Henri Poincaré (1854−1912), who established the connection between two topological bodies through his intuitive ideas and translated these ideas into algebraic topology as a symbolic and ascertainable formalization of the emerging theory of homotopy. The central concept at this step of topology was “the homotopic deformation”,11 which can simply be interpreted as a “continuous mapping from an N dimensional space to an N + 1 dimensional space, where the extra variable can take values from 0 to 1 (eq 2) relating two functions defined in the N dimensional space (X and Y)”. A rigorous definition of “homotopic deformation” can be found in numerous textbooks and research and review articles. H: X × [0, 1] → Y

H(x , t ) = tF(x) + (1 − t )E(x) = 0

(2)

where F(x) is the original system to solve, E(x) is the simpler system proposed at the beginning of the method, and t is a parameter used to effect the homotopic deformation. The topological path known as the “homotopic deformation path” can be computed as a set of discrete points (sequence of deformation steps) that can be calculated by a classical Newton method. The homotopic path can be viewed as an Ndimensional graph that facilitates its management and comprehension. Any system of two equations would have a 3-D curve like those shown in Figure 3. Note that if the homotopic parameter t = 0 is substituted in the mentioned equation, H(x, t) = E(x). Furthermore, t = 1 yields H(x, t) = F(x). The N-dimensional trajectory, where H(x, t) = 0 is the homotopic path H−1 and its spatial location is not known a priori, in principle, connects the initial approximation to the solution vectors located on the hyperplane t = 1. The choice of E(x) and the method that varies t from 0 to 1 are nontrivial aspects because they determine the number of solutions found and the sensitivity of the method to the initial guesses (Appendix A, Supporting Information, shows a detailed example for the solution of a system of two equations). Although a search of the concurrent literature reveals applications of “Newton homotopy” (eq 3), the first simpler system proposed within homotopic methods may have been the “fixed-point homotopy” (eq 4) because one of the previous mathematical steps in the works of Lahaye15,16 was the “fixed-point theorem”.12,13 This “fixed-point theorem” is philosophically and theoretically linked to “fixed-point homotopy”. Nevertheless, both the Newton and fixed-point approaches present numerical setbacks for certain problems such that they would derive in infinite cycles for the variable values that diverge if the homotopic parameter does not reach the value of 1 or localize repeated roots when t = 1. Modern homotopic methods result from the evolution of the three listed points and have spawned several other methods, such as piecewise linear or simplical methods. However, we

(1)

To clarify the concept of homotopic deformation, consider the mathematical abstraction (coordinates) of some points that constitute a cube in 3-D space. The cube can be compressed to obtain a squared plane (see Figure 1), so the 3-D coordinates

Figure 1. Diagram of the compression and expansion of a cube/plane.

can be mapped to mathematically describe the plane. In other words, an intuitive and geometrical meaningful process is transformed into an abstract mapping between two vector spaces (a 2-D and 3-D spaces). Because any group of points defined in a coordinate system forms a vector space, the geometrical analogy for its homotopic deformations is obvious. In Figure 2, the homotopic mapping is presented as a function applied on a meta-domain (named A), which not only has the original 3-D coordinates of the cube in its range but also the compressed one point with 2-D coordinates of the square. In the solution of a nonlinear algebraic system, the 3-D coordinate points would represent to the complex original algebraic system (at the t value of 1), and the 2-D coordinate points eventually 14730

dx.doi.org/10.1021/ie402418e | Ind. Eng. Chem. Res. 2013, 52, 14729−14742

Industrial & Engineering Chemistry Research

Review

Figure 2. Schematic diagram of the homotopic deformation between two geometric solids and its abstraction through vector mapping.

Figure 3. Homotopic paths for the same nonlinear algebraic system presented in Appendix A (Supporting Information). (a) Using Newton homotopy. (b) Using fixed-point homotopy.

engineering field spanning from last year to the current date is presented. 2.3. First Challenges. Researchers encountered several difficulties after the first algorithms were proposed. For example, the primitive methods tend to fail in some systems when the homotopic path presents a behavior called “return”. A practical definition of a homotopic path “return” would be a numerical return in the homotopic parameter values (t in eq 2) due to the local convergence character for each point of the homotopic path. The fact that two geometrically noncontiguous steps in the homotopic deformation process coincide in at least one homotopic parameter value (the planes

only refer to predictor−corrector methods in this work. Complete reviews of all current algorithms and their historical development in the context of a rigorous mathematical discussion have been put forth by Seydel,20 Rheinboldt,10 and Allgower and Georg.7 More engineering-friendly approaches that focus on “predictor−corrector methods” have been published by Garcia and Zangwill,7 Seydel and Hlavacek,21 Wayburn and Seader,22 and Seader et al.,23 among others. Thus, we present a brief summary designed to contextualize that ranges from 1961 to 2003 and focuses on computational applications. Furthermore, an original review in chemical 14731

dx.doi.org/10.1021/ie402418e | Ind. Eng. Chem. Res. 2013, 52, 14729−14742

Industrial & Engineering Chemistry Research

Review

Figure 4. Homotopic path (x1 ⊥ t ⊥ x2) for the example of Appendix A (Supporting Information) employing an initial guess of x0 = [2, 2]t.

Figure 5. Disadvantages of horizontal predictions in homotopic path tracking methods. Example of an arbitrary homotopic path with (a) divergence of the correction steps and (b) erroneous path constructed due to bad predictions.

Figure 6. Advantages of the differential predictions against horizontal prediction in a typical homotopic path followed by a predictor−corrector. (a) Differential predictors and (b) horizontal predictors and the loss of two hypothetic solution vectors.

For example, simply augmenting the homotopic parameter value (horizontal predictors) would have serious disadvantages for certain problems. Consider the homotopic path shown in Figure 5. A prediction located outside of the convergence zone

shown in Figure 4 intersects the homotopic path at two different points) was a significant challenge for the first algorithms. 14732

dx.doi.org/10.1021/ie402418e | Ind. Eng. Chem. Res. 2013, 52, 14729−14742

Industrial & Engineering Chemistry Research

Review

Figure 7. Homotopic path of example of Appendix A (Supporting Information). (a) Correct path. (b) Correct sequence of construction from t = 0. (c) Incorrect path and its slopes. (d) Wrong sequence of path tracking due to the plane homotopic slope making a false return.

2.4. The Bifurcation Approach. After the pioneering work of Haselgrove,19 a number of authors began to use “differential predictions” that are numerically very similar to the Jacobian predictions in the N-dimensional Newton method. It is not surprising that some numerical errors emerge when the homotopic slope becomes planar (homotopic Jacobian approaches a singularity), for a geometrical interpretation, consider the correct construction of a return (Figure 7a). Because of the tangent at every point in the homotopic path is required to predict the next point, the presence of a return (horizontal tangent in Figure 7c) would make difficult the prediction of the point that truly follows in the correct sequence of the homotopic path construction (Figure 7b). If the appropriate strategies are not implemented, the backward reconstruction of the homotopic path becomes probably (Figure 7d). A second, similar proposal was made within the differential predictor approach. The pure mathematical theory of topology

surrounding the homotopic path results in one of two possibilities: the correction steps could either diverge (because they are iterations of local methods) or converge to another section of the homotopic path, which could result in a mistake commonly named “curve-jumping”24−26 (see Figure 5). The contributions of Haselgrove19 were made in the above context. His proposal constitutes a differential prediction, which is a more robust approach used today. This differential prediction was so beneficial to the homotopic methods that some classical authors have referred to it as “the stabilization of the approach of Davidenko” (see Figure 6). Typical examples of returns and their disadvantages for simple differential and horizontal predictors are shown in Figures 5 and 6. Figure 6 indicates that the close proximity of all predictions to the correct section of the homotopic path (Figure 6a) is very important. If this condition is violated (Figure 6b), some solution vectors may be lost. 14733

dx.doi.org/10.1021/ie402418e | Ind. Eng. Chem. Res. 2013, 52, 14729−14742

Industrial & Engineering Chemistry Research

Review

that the challenging problems pointed by Choi and Book43 can be solved without losing solution vectors if the homotopic parameter is also capable of taking on complex values. According to Seader et al.,23 there is one additional possibility for complex homotopic continuation; the entire path can be tracked using complex arithmetic, which locates both the complex and real roots. This approach has been considered disadvantageous in the literature because the majority of chemical engineering problems are not modeled by algebraic systems with physical meanings for the complex values of its variables. Furthermore, the solution method tends to be slow and challenging for large systems. However, this method is an active topic in the literature.45−50 Many of the important computational features in the complex domain have found applications in control theory, fluid mechanics, and chemical reaction systems, which are mentioned in the contemporary applications section. If the reader is interested in this subject, we suggest the work of Chow et al.51 2.5. Normal Flow, Differential Equation-Based, and Augmented Jacobian Strategies. The classical work of Watson52 mentions four steps for the “prediction−correction methods”: prediction, correction, step size estimation, and computation of each solution. There are three general options for the prediction−correction steps. The first option is the “normal flow approach” (Figure 9a), which yields linear predictions (horizontal or differential increases in the homotopic parameter; see Figures 5 and 6) and the Newton method, which constitutes its correction steps. This method has been questioned because it requires a large number of calculations,52 which would have been especially disadvantageous in the past when computers were not powerful. The traditional literature presents an alternative to this calculationintensive method, which is well-accepted even today: the “quasi-Newton augmented Jacobian” method, which introduces a polynomial function to predict values rather than the lineal ones used by normal flow algorithms (see Figure 9b). Furthermore, the Jacobian matrix that maps the homotopy is augmented with a tangent vector and updated QR-factorizations and quasi-Newton calculations constitute the correction steps, rather than the complete Newton method used by the normal flow strategies. This method has a particular advantage that is very useful for certain modern methods:50,53,54 it permits a tangent vector location when the homotopy parameter value equals 1, which avoids numerical instabilities for the Jacobian matrix singularities at the solution.

was again used to analyze the method, and the concept of bifurcation was introduced to detect critical points in the homotopic path where the returns occur.27−36 For simplicity, one type of a bifurcation could be interpreted as a set of values for the variables of the algebraic system that correspond to a singularity in the homotopic Jacobian matrix.7,37 However, it is important to note that bifurcation is a coarser, deeper, and more complicated mathematical concept. If the reader is interested in a description of this topic, we suggest the work of Allgower and Georg7 to introduce formal definitions in the context of prediction−correction methods and the papers of Seydel38 and Meijer and Kalmár-Nagy39 for discussions of the dynamical analysis context. Cayley40 and Brolin41 studied the behavior of complex variable iterative maps in mathematical physics. Many important contributions to the theory of nonlinear systems, such as bifurcation theory, have stemmed from their studies. Within the context of bifurcations, a second type of homotopic path tracking algorithm, “the complex domain continuation”, was generated. Because singularities in the homotopic Jacobian (turning points) occur in many cases, the homotopic path takes on complete complex values, which implies that the imaginary part of the homotopic curve does not equal zero (Figure 8).

Figure 8. Example of a turning point with an obvious horizontal slope coinciding with a complex bifurcation (modified from Wayburn and Seader22).

Some authors have suggested that complex homotopic branches could connect isolated segments of the homotopic path or opened paths to infinity with bounded segments that connect real solution vectors22,42 while only permitting imaginary numbers for the variables of the nonlinear system. Choi and Book43 demonstrated that some systems have some “unreachable roots” when this approach is used, even if the path is later tracked in complex space. Gustafson and Beris44 showed

Figure 9. Schematic representation of (a) the augmented Jacobian strategies and (b) normal flow methods. 14734

dx.doi.org/10.1021/ie402418e | Ind. Eng. Chem. Res. 2013, 52, 14729−14742

Industrial & Engineering Chemistry Research

Review

space (x1, x2, and the homotopic parameter t for a two equation example) with a common vertex located in the initial guess. The increments can be mapped in an auxiliary space (P ⊥ t) where they have a more manageable visualization when considered as a side hick of a right triangle used to measure every segment of the arc-length in a discrete way. Thence, eq 2, which is defined in terms of the original problem variables (x ∈ Rn in vector form) and the homotopic parameter t, is reexpressed in terms of an auxiliary vector that consists of a grouping of these variables, y ∈ Rn+1; thus, the homotopic path is a function of its own continuation process, which ensures that every path tracking point differs. This feature prevents the possible false infinite cycles from succumbing to turning points and thus adequately predicts values to construct complete and real paths. Recently, the contributions and applications of these ideas have been noted. Specifically, Dhooge et al.54,55 have developed the convenient and widely available toolbox CL_MATCONT in MATLAB, which is a modern suite that primarily consists of the augmented Jacobian technique. However, modern methods are the result of the combination of the three or more approaches. 2.6. Different Predictors. Furthermore, additional, differing predictors were published concurrently with the pseudoparameterization of the arc length. For example, polynomial predictors were inspired in the sequential numerical construction of the homotopic path. Thus, we have every previous point of the homotopic curve, and polynomials, which are functions capable of adjusting known data to extrapolate one or more predictions (see Figure 11a). This characteristic can become problematic when estimating the parameters of cubic, Hermite60 or Langrage61 polynomials, and many other available algebraic forms. Closed predictors also have also been described in the literature. This approach was first described by Lyness and McHugh,62 who used hypercubes to achieve progressive integration. Nevertheless, not until the work of Jiménez-Islas9 was the first algorithm for closed hypersurfaces (see Figure 11b) proposed in the context of homotopic methods. This strategy uses classical differential predictors and some auxiliary equations, like a generalization of the sphere equation to Ndimensions to facilitate the Newton correction steps. Also, thorough comparison with other closed surfaces shapes have been performed.9 Other authors have published several works using the hyperspherical path tracking methodology in diverse science fields.62−65 Furthermore, recently published works

A third strategy that has often been reported in the literature 9,22,52,55−59 is the “differential equation based approach”. In this method, an artificial additional parameter P is introduced and can be geometrically interpreted as the arclength of the homotopic path. According to this scheme, the path is calculated and analyzed as a 2-D curve (see Figure 10)

Figure 10. Homotopic path for the example of Appendix A (Supporting Information) and its corresponding arc-length parametrized path (P ⊥ t).

composed of ordered pairs (P ⊥ t) where P is an advance parameter. For practical purposes, P is ever referred to the initial guess as sum of discrete increments. Every increment would represent the base of a non-right triangle in the N + 1

Figure 11. Schematic diagram of (a) a generic polynomial predictor and (b) closed hypersphere predictor. 14735

dx.doi.org/10.1021/ie402418e | Ind. Eng. Chem. Res. 2013, 52, 14729−14742

Industrial & Engineering Chemistry Research

Review

homotopic path tracking or on the variation of the homotopic slope.9 2.8. Different Homotopies. As discussed in The Numerical Method section, the selection of a simple system of algebraic equations that will be deformed by the method is not a trivial task. Many authors have used distinct algebraic forms, and very few guidelines have been published to select or construct these forms. Historically, the relevant publications that have been most cited discuss the two well-known homotopy forms, Newton and fixed-point. Other alternative algebraic forms that have since been proposed are listed below.

feature graphical analyses and demonstrations of this technique.68 2.7. The Step-Length in the Continuation Process. The literature has not solely focused on the predictor method; the step size of the prediction has also been an active topic of discussion. This step size is entirely linked to the prediction method; the reason for this connection is obvious and presented in Figure 5. A variety of strategies were proposed, which can be classified into two categories. The first such category is a rigorous approach that is based on the famous theorem by Newton-Kantorovich.60,61,69,70 This approach can be explained as follows: consider that a ball exists around a discrete point of the homotopic curve within which all possible predictions converge at the center (see Figure 12). Once the

E(x) = F(x) − F(x 0)

(3)

2.8.1. Newton Homotopy. The Newton homotopy (eq 3) yields bounded homotopic paths in finite domains.23 However, as with any other homotopic form, some paths have closed sections in the real domain (isolas) and opened paths that connect some solutions after virtually infinite tracking distances. See in Figure 13 an example where the homotopic path is constituted by an isola (shown as a segmented line) with four solution vectors and an opened section (shown with a continuous line) that effectively contain five additional solutions. Some publications have recently begun to analyze and characterize the demeanor of paths for every homotopy. A high frequency of homotopic paths with isolas has been identified for the Newton method.69 To illustrate the presence of opened homotopic paths, consider the simple system of two nonlinear algebraic equations shown in Appendix A (Supporting Information). This system was solved by Newton homotopy with an initial guess of x0 = [6, −5.2]t (Figure 3). The two segments shown are not connected to each other, and it is not possible to continue the process after the homotopic parameter reaches a value of 242182.7841, which indicates bifurcation. This bifurcation may be complex, but a continuous real homotopic path for this problem is not available irrespective of this complexity. Figure 13 shows that homotopic paths can be defined by more than one continuous segment. It is not easy to construct such paths and is only presented for descriptive purposes

Figure 12. Geometrical representation of the Newton-Kantorovich theorem describing a ball of convergence around one point of an arbitrary homotopic path.

radius of the ball of convergence is known, simple geometric deductions can be made to derive expressions that are very useful to predict the maximum radius length of the prediction capable of triggering convergence. The second strategy proposed by some authors consists of complete heuristic methods to predict the step length; these were based on the variation of the variable values along the

Figure 13. The two segments of the complete homotopic path for a system of two equations shown in Appendix A (Supporting Information). The isola (segmented line) contains four roots not presented in the opened path. A magnification of the graph shows five roots in this last one. 14736

dx.doi.org/10.1021/ie402418e | Ind. Eng. Chem. Res. 2013, 52, 14729−14742

Industrial & Engineering Chemistry Research

Review

2.8.6. Multiparameter Homotopies. Multiparameter homotopy formulations are a group of strategies that consider more than one homotopy parameter,74 unlike traditional approaches (eq 2). The behavior of these mappings is very similar to the ideas noted in the previous sections because the equilibrium equation (eq 8) has one trivial solution at t1 = 0, t2 = 0, t3 = 0, ..., tn = 0, and the sought solution vectors are reached at t1 = 1, t2 = 1, t3 = 1,..., tn = 1 where n is the total number of homotopic parameters used.

because segments that do not connect with the initial guess have to be calculated by initializing the path tracking process from known roots. E (x ) = x − x 0

(4)

2.8.2. Fixed-Point Homotopy. Even though fixed-point homotopy (eq 4) has drawbacks that are similar to that of the Newton form (discontinuous paths), some heuristic criteria42 are available to choose initial guesses that allow the localization of all the roots for many systems and only one continuous homotopic path. However, several authors have reported the failure of these criteria.43 As such, Kangas et al.71 found that the fixed-point homotopy yields paths that converge to only certain solution vectors, which means that there are discontinuous homotopic paths in real domain. At this point, it is necessary to say that in other engineering fields Affine homotopy (Section 2.8.3) is placed in the same category of fixed-point homotopy.72 Semantically this is not a mistake. However, the chemical engineering literature handles the terminology outlined in this article, so that if the reader is interested in the solution and understanding of some study cases, we suggest the separation of the two mentioned categories. E(x) = A(x − x 0)

H(x , t1 , t 2 , ..., tn) = 0

There are several examples of multiparameter homotopies.65,66 However, the presentation of the easiest case (biparameter homotopy) is sufficient to understand and appreciate the advantages and the philosophy of these techniques from a comprehensive perspective. Two simultaneous deformations are produced from two parameters (t1 and t2), one in the function of the original algebraic system F(x) and another in the homotopic mapping function H(x, t). The initial guess is defined by [x, t1, t2]t = [x0, 0, 0]t (named p1 from now), and when [x, t1, t2]t = [x*, 1, 1]t, we arrive at one nontrivial solution vector of F(x) (named p2 from now). Thus, we aim to graphically derive an interpolation expression that satisfies the two described conditions (see Figure 14).

(5)

2.8.3. Affine Homotopy. Affine homotopy (eq 5) presents a significant improvement in that it combats some of the disadvantages of the fixed-point form. This method defines a weight matrix, A, which reduces numerical difficulties caused by the disparity of scales between the variables of a system.7 This problem is very common in systems with multiple nonlinear algebraic operators of different types. Kangas et al.71 identify the frequency of unfeasible roots for this homotopy and distinguish this situation from the isolas for the Newton form. However, no behavior that is distinct to isolas or infinity-opened paths has been reported as a cause of unfeasible roots. E(x) = x − x 0 + F(x) − F(x 0)

Figure 14. Geometrical interpretation for the interpolation function of the biparametric homotopy.

(6)

The equation that can describe the path shown is named the parametric function M(t1, t2). This equation traverses three points p1, p2, and p3, where p3 is an arbitrarily selected point. M can be derived in several ways; one of them is the eq 9.64,74,75

2.8.4. Combined Fixed-Point Newton Homotopy (FPN). The combined FPN homotopy (eq 6) is a linear superposition of the Newton homotopy and fixed-point form that has the advantages of both methods. Specifically, this approach features criteria to choose the initial approximations and trigger the construction of finite, bounded, and continuous paths that connect all the roots of one nonlinear equation.49 This homotopy was extrapolated to systems of nonlinear algebraic equations and is a form that is free of isolas with multiple advantageous characteristics, including the use of bifurcation points (as opposed to a drawback), the existence of a limit point (which can be employed to implement formal stop criteria), and an extensive number of transcendental systems that have been successfully solved.50 The only drawback for this system is that one root of the original system is required to initialize the method. However, the authors proposed a very efficient method to localize the initial guesses of their methods.

E (x ) =

1 || x − x 0 ||2 2

(8)

M(t1 , t 2) =

(t

2

+

B(A − 1) AB + 1 − 2A

B(A − 1) 2 AB + 1 − 2A



)

2A − B − 1 AB + 1 − 2A

− t1 (9)

Equation 9 necessarily participates in the definition of homotopy. Furthermore, certain “stop criteria” can be applied because we know that the path has finished at t2 = 1.67 This finding is probably the most important contribution of the multiparameter homotopies; the fact that we can impose limits to the homotopic path, eliminating opened homotopic paths that contain solutions after practically infinite path tracking distances (a problem addressed in Section 2.8 and graphically shown in Figure 13). To date, no homotopy or tracking method ensures the complete construction of a homotopic path. In other words, there is no way to ensure that all of the roots of a problem have been located. In homotopy theory, the only premise is that the homotopic path connects all the solutions of a system, but there is no reason to suppose that a homotopic path is continuous and finite, and the selection of the algebraic form

(7)

2.8.5. Probability One Homotopy. In the case of the probability one homotopy (eq 7), a single root can assuredly be located that represents the global optimum for optimization.73 14737

dx.doi.org/10.1021/ie402418e | Ind. Eng. Chem. Res. 2013, 52, 14729−14742

Industrial & Engineering Chemistry Research

Review

from an N-dimensional generalization of the sphere equation) in the homotopic system in every correction step.

for E(x) continues to be a heuristic and crucial aspect of homotopic methods. 2.9. Mapping of Homotopic Paths. A more recent topic of discussion has been the approach of imposing limits on any of the homotopic tracking variables. The strategy has certain advantages. The first advantage is that homotopic paths can be tracked at will in selected domains,76−78 which is very useful for certain problems where the variables should be physically meaningful, and the solutions must fall within this restriction. The second advantage is that the largest possible number of roots can be discovered because the delimitation of the homotopic parameter is a robust alternative that circumvents certain difficulties such as the isolas and opened homotopic paths that must be built after calculating an infinite number of steps in the tracking process. The bounded homotopies for small systems are mathematically defined as the mapping of tracking variables through the linear superposition of two functions. One function is defined by the product of a penalty function and the common homotopy function, which annihilates the value of the path H−1 when the variables exceed the limits. The other function has an offsetting effect that allows the tracking process when the variables are defined outside the domain and are annihilated.74 Although bounded homotopic paths have some advantages, they also have some drawbacks. For example, bounding maps produce compressed topologies, which complicate traditional predictor−corrector iterations. Furthermore, the homotopic path suffers abrupt changes in its slopes near the frontiers of the selected domain. These changes trigger numerical difficulties because the prediction lies outside the mathematical definition of the problem (a geometrical abstraction of this effect is depicted in Figure 15). Three alternatives have been proposed

3. CONTEMPORARY APPLICATIONS With the development of global methods to solve nonlinear systems of algebraic equations that include homotopic methods, different mathematical tools have become necessary, and different analytic approaches have emerged in the literature. In chemical engineering and related fields, phase separation problems have likely been most studied by homotopic methods. A brief and analytical overview with a historical perspective is presented below. 3.1. Phase Separation Problems. Many chemical process design problems, such as the design of distillation columns or alternative separation systems, involve phase equilibrium. For nonreacting multicomponent systems, the problem is the determination of the number of phases that exist at equilibrium and the composition and quantity of each phase. Methods for calculating the equilibrium composition of a mixture of chemical species can be divided into two groups: 1. Equilibrium-constant methods 2. Free-energy-minimization methods The advantages and disadvantages of both approaches have been continually discussed since 1960. The preferential approach depends on the time period and scientific field in question. Prior to 1970, chemical engineers preferred equilibrium constant methods, while the rocket industry preferred free-energy-minimization methods.79 Over time, the second method was extended to chemical engineering and other disciplines.80,81 With respect to the tools that manage and solve mathematical models derived from equilibrium phase methods, several strategies used for group 1 solve the phase equilibrium equations. In one approach, Chavez et al.56 applied homotopic continuation from multiple starting points to find all real positive solutions to three examples of interlinked distillation systems where a physical equilibrium was assumed at each stage. The work was extended by using a global homotopy continuation to find all real positive solutions from just one starting point by continuing the path until all roots were found. To date, the majority of efforts have been directed to construct correct homotopic paths by accurately calculating it in discrete points. The mentioned works constitute a paradigm shift in homotopic methods. Almost all of the reports that have been published since the 1980s have focused on finding a good number of solution vectors.46,82−85 Some attempts were made to ensure the localization of all the zeros of particular nonlinear algebraic systems.50,57,86−88 Optimization schemes for group 2 have received more attention, which has significant implications. The nonlinear objective function is the Gibbs free energy, which is subjected to a set of linear and nonlinear constraints (mass balances and non-negativity of moles). The mathematical properties of the nonlinear objective function depend completely on the mathematical form of the equation of state (EOS) or fugacity correlations chosen to represent each of the phases that may exist at equilibrium. Nonconvexities can appear due to the form of the objective function and can lead to metastable equilibria or points whose multiphase solutions converge to a trivial solution. The global optimal solution of the mathematical formulation corresponds to the true physical equilibrium solution. Therefore, equilibrium compositions are given at the

Figure 15. Schematic representation of an arbitrary bounded homotopic path of one nonlinear algebraic equation (modified from Malinen and Tanskanen57).

to date to overcome these situations. An obvious first solution is to maintain a very short step size of the predictor to ensure that predictions fall sufficiently close to the next homotopic step.76,77 Another possible solution, proposed by Malinen and Tanskanen58 and Malinen,101 consists of a second mapping (the bounding approach is a mapping itself) to follow the unbounded homotopic paths in bounded artificial domains, which allows the construction of correct compressed homotopic paths. The third solution, proposed by OliverosMuñoz and Jiménez-Islas,68 implements a special path tracking methodology named “the hyperspherical path tracking methodology”. This approach includes an auxiliary equation (derived 14738

dx.doi.org/10.1021/ie402418e | Ind. Eng. Chem. Res. 2013, 52, 14729−14742

Industrial & Engineering Chemistry Research

Review

Table 1. Applications of Homotopic Methods in Chemical Engineering application fluid mechanics

contribution application of the Yorke algorithm on FORTRAN a known code named HOMPACK 77

elliptic porous slider compression of a fluid between two plates deceleration of rotating disks for viscous fluids flow in a porous channel of a rotational system optimization (nonlinear complementary problems) absorption columns different configurations in distillation columns interconnected

homotopic first tracking algorithms first applications of the differential tracking method with respect to the arc length of the homotopic path reaction in a catalytic particle, diffusion problem with auto catalysis and implementation of tracking algorithms with a focus on the management of hypersonic flow branches. engineering flowcharts affine homotopy to prevent cyclic paths resolution of equations diffusion type in semiconductors combination of finite element method with tracking arc length parametrization of the homotopic path numerical solution of mass balances for continuous stirred reactors application of transformation methods for selecting the initial guesses used in with radical intermediates fixed-point homotopy calculating azeotropes in more than one phase for multicomponent application of methods of detection of bifurcations mixtures analysis of reactive multiphase systems application tracking algorithm with forks and branches handling complex homotopy path solution of mathematical models (partial differential equations) steady solution of systems resulting from the discretization of differential equations state a solution of reduced model for controlling an activated sludge process. application of homotopic methods in biological reacting systems parameter estimation for elliptic PDE combination of monitoring homotopic to the method of Tikonov regularization simulated distillation columns with calculation of the equilibrium application of penalty functions for parameter homotopy of homotopies defined curves for ternary mixtures.

ref 52

103 56 21 22 104 105 106 107 108 109,110 111 58

achieved by an explicit expression to describe complex systems. Some examples of the applications are the scale-up of mass and heat transfer phenomena in biofilms,94 water transport in porous media,95 diffusion reaction models,96 inverse heat conduction problems,97,98 general Sturm-Liouville problems,99 and many others. If the reader is interested in the analytical solution of differential equations that use this approach, the work of He100 is suggested. Finally, more applications are available in the literature, and a short list of the most cited works is shown in Table 1.

global minimum of the Gibbs free energy surface. However, poor initial guesses for the phase distribution and composition can lead to local minima. Constrained minima occur when too few phases are assumed or when the correct number of phases is assumed but the composition guesses are poor. Furthermore, local minima can occur when too many phases are assumed. The extremely difficult task of finding the multiple minima for nonideal systems was attempted by Jalali et al.46 through some contributions to the application of the homotopic methods for constructing paths that combine real and imaginary spaces, a work proposed first by Kuno and Seader.42 Furthermore, the most novel approaches have been used in phase separation problems. Kangas et al.71 used the bounded homotopy strategies of Malinen and Tanskanen57 to develop a robust method to calculate global minima. The general demeanor of the three classical homotopy forms was characterized based on their studies (Newton, FPN). They observed that for Newton homotopy, isolas topologies are very common, while fixed-point and affine homotopy feature unfeasible roots. 3.2. Transport Phenomena Problems. The aim of this paper is to analyze the numerical methods grouped by the predictor−corrector homotopic continuations with a modest but novel perspective. However, the emergence of a vigorous topic necessitates the dedication of a few brief lines to their effect. In recent years, an increasing number of articles have been published that use homotopy theory in a way that is similar to the predictor−corrector methods used to solve differential equations. The pioneering works of Liao,88−91 Liao and Cheung,92 and He93 exploit the topologic connection between a very difficult problem and one that is easier to develop with two analytical strategies: the homotopic analysis and the perturbation analysis. To date, many papers have been generated by a significant number of authors that address a variety of physical problems with the analytical perspective

4. CONCLUSIONS AND FUTURE WORK The real problems in chemical engineering have been addressed by many approaches. One of the more robust approaches is the application of “the homotopic methods”, which are a very simple and intuitive interpretation of a powerful idea based on the topological connection between two solid bodies by a deformation process. In an abstract way, we can observe the same relationship between two algebraic forms. Thus, to solve equations or systems of equations with the appropriate proposition of one simple algebraic system or algebraic equation, we can deform the system to obtain the very complex system that motivated our interest. In the literature, many problems whose solutions emerge in the development of numerical tools are related to phase equilibriums calculations. Thus, a short review of contributions to this field was presented that focused on the mathematical aspects of the topic. Additional, similarly important efforts from authors in the chemical engineering literature were also listed. No one method to date is totally secure and rigorously demonstrable in finding all the solution vectors for all the nonlinear algebraic systems. The more successful attempts range from the homotopic continuation in the complex domain to the bounding of the homotopic parameter and other 14739

dx.doi.org/10.1021/ie402418e | Ind. Eng. Chem. Res. 2013, 52, 14729−14742

Industrial & Engineering Chemistry Research

Review

(9) Jiménez-Islas, H. (In Spanish) SEHPE: Programa para la solucion de sistemas de ecuaciones no lineales mediante metodo homotopico ́ 1996, 6 (2), 174−179. con seguimiento hiperesferico. Av. Ing. Quim. (10) Rheinboldt, W. C. Numerical Continuation Methods: a perspective. J. Comput. Appl. Math. 2000, 124, 229−244. (11) Stochastikson [Online]; Encyclopedia; Biography of Jules Henri Poincaré. http://encyclopedia.stochastikon.com/ (accesed Sep 10, 2009) (12) Brouwer, L. E. J. (In German), Beweis der Invarianz der Dimensionenzahl. Math. Ann. 1911, 70, 161−165. (13) Brouwer, L. E. J. (In German), Ü ber Abbildung von Mannigfaltigkeiten. Math. Ann. 1912, 71, 97−115. (14) Leray, J.; Schauder, J. (In French), Topologie et équations fonctionnelles. Ann. Sci. de L’É.N.S. 1934, 3 (51), 45−78. (15) Lahaye, E. (In French), Une méthode de resolutio d’une catégorie d’équations transcendentes, C. R. Acad. Sci. Paris 1934, 198, 1840−1842. (16) Lahaye, E. (In French), Sur la résolution des systémes d’équations trascendantes. Acad. R. Belg. Bull. Cl. Sci. 1948, 34 (5), 809−827. (17) Ficken, F. A. The Continuation Method for Functional Equations. Commun. Pure Appl. Math. 1951, 4, 435−456. (18) Davidenko, D. On a new method of numerical solution of systems of nonlinear equations (in Russian). Dokl. Akad. Nauk USSR 1953, 88, 601−602. (19) Haselgrove, C. B. The Solution of Nonlinear Equations and of Differential Equations with Two-Point Boundary Conditions. Computer. J. 1961, 4 (3), 255−259. (20) Seydel, R. From Equilibrium to Chaos: Practical Bifurcation and Stability Analysis, 1st ed.; Elsevier Science Ltd.: New York, 1988. (21) Seydel, R.; Hlavacek, V. Role of Continuation in Engineering Analysis. Chem. Eng. Sci. 1986, 42, 1281−1295. (22) Wayburn, T. L.; Seader, J. D. Homotopy Continuation Methods for Computer-Aided Process Design. Comput. Chem. Eng. 1987, 11, 7−25. (23) Seader, J. D.; Kuno, M.; Lin, W.; Johnson, S. A.; Unsworth, K.; Wiskin, J. W. Mapped continuation methods for computing all solutions to general systems of nonlinear equations. Comput. Chem. Eng. 1990, 14, 71−85. (24) Bates, D. J.; Hauenstein, J. D.; Sommese, A. J.; Wampler, C. W. Adaptive multiprecision path tracking. SIAM J. Numer. Anal. 2008, 46 (2), 722−746. (25) Bates, D. J.; Hauenstein, J. D.; Sommese, A. J.; Wampler, C. W. Stepsize control for adaptive multiprecision path tracking. Contemp. Math. 2009, 496, 21−31. (26) Vazquez-Leal, H.; Marin-Hernandez, H.; Khan, Y.; Yildirim, A.; Filobello-Nino, U.; Castaneda-Sheissa, R.; Jimenez-Fernandez, V. M. Exploring collision-free path planning by using homotopy continuation methods. Appl. Math. Comput. 2013, 219, 7514−7532. (27) Klopfenstein, R. W. Zeros of nonlinear functions. J.A.C.M. 1961, 8, 336−373. (28) Brezzi, F.; Descloux, J.; Rappaz, J.; Zwahlen, B. On the rotating beam: Some theoretical and numerical results. Calcolo 1984, 21, 345− 367. (29) Brezzi, F.; Rappaz, J.; Raviart, P. A. Finite dimensional approximation of nonlinear problems. Part 1: Branches of nonsingular solutions. Num. Math. 1980, 36 (1), 1−25. (30) Brezzi, F.; Rappaz, J.; Raviart, P. A. Finite dimensional approximation of nonlinear problems. Part 2: Limit points. Num. Math. 1980, 37 (1), 1−28. (31) Brezzi, F.; Rappaz, J.; Raviart, P. A. Finite dimensional approximation of nonlinear problems. Part 3: Simple bifurcation points. Num. Math. 1981, 38 (1), 1−30. (32) Golubitsky, M.; Schaeffer, D. G. Singularities and Groups in Bifurcation Theory, Vol. I; Springer Verlag: New York, 1985. (33) Golubitsky, M.; Stewart, I.; Schaeffer, D. G. Singularities and Groups in Bifurcation Theory, Vol. II; Springer Verlag: New York, 1988. (34) Spence, A.; Jepson, A. D. The numerical calculation of cusps, bifurcation points and isola formation points in two parameter

mappings. Some of these attempts are very recent and should be explored further. For example, the combined bounding approaches of Malinen101 and Paloschi76−78 to obtain bounded paths both for the homotopic parameter and for the original system variables as well as the extensive validation and application of the promising homotopic form of Khaleghi Rahimian50 both should be explored further. Furthermore, some novel applications are emerging; the global optimization techniques and semianalytic solutions of partial differential equations, which are more developed techniques in other science and engineering fields,102,103 where relevant modifications to the original perturbation methods have been proposed.



ASSOCIATED CONTENT

S Supporting Information *

Appendix A. This material is available free of charge via the Internet at http://pubs.acs.org.



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Tel: +52 (461) 611-7575. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS We gratefully acknowledge the financial support of SEPCONACYT, Mexico, via grant Basic Science CB-2011-01167095 and DGEST Grant 5037.13-P.



NOMENCLATURE H(x,t) = homotopic mapping F(x) = system of algebraic equations for vector of variable x E(x) = simpler system derived from the original problem t = homotopic parameter x = vector of the original variables x* = solution vector A = weight matrix for the affine homotopy x0 = vector of the initial guess



REFERENCES

(1) Kovach, J. W., III. Heterogeneous azeotropic distillation: homotopy-continuation methods. Comput. Chem. Eng. 1987, 11 (6), 593−605. (2) Müller, D.; Marquardt, W. Experimental Verification of Multiple Steady States in Heterogeneous Azeotropic Distillation. Ind. Eng. Chem. Res. 1997, 36, 5410−5418. (3) Fussmann, G. F.; Elmer, S. P.; Shertzer, K. W.; Hairston, N. G. Crossing the Hopf bifurcation in a live predator−prey system. Science 2000, 290, 1358−1360. (4) James, A.; Brindley, J. Stability of Multiple Steady States of Catalytic Combustion. Combust. Flame 2002, 130, 137−146. (5) Hao, W.; Hauenstein, J. D.; Hu, B.; McCoy, T.; Sommese, A. J. Computing steady-state solutions for a free boundary problem modeling tumor growth by Stokes equation. J. Comput. Appl. Math. 2013, 237 (1), 326−334. (6) Tolsma, J. E.; Barton, P. I. In Software Architectures and Tools for Computer Aided Process Engineering; Braunschweig, B., and Gani, R., Eds.; Amsterdam, The Netherlands: Elsevier Science B.V., 2002; Chapter 3.2 Numerical Solvers, pp 127−162. (7) Allgower, E. L.; Georg, K. Introduction to numerical continuation methods. Classics in Applied Mathematics. Soc. Ind. Appl. Math. 2003. (8) García, C. B.; Zangwill, W. I. Pathways to Solution Fixed Point in Equilibria, 1st ed.; Prentice Hall: New York, 1981. 14740

dx.doi.org/10.1021/ie402418e | Ind. Eng. Chem. Res. 2013, 52, 14729−14742

Industrial & Engineering Chemistry Research

Review

problems. In Numerical Methods for Bifurcation Problems; Küpper, T.; Mittelmann, H. D.; Weber, H., Eds.; Birkhäuser: Verlag, Basel, 1984; pp 502−514. (35) Cliffe, K. A.; Spence, A. In Numerical Methods for Fluid Dynamics, Numerical Calculation of Bifurcations in the Finite Taylor Problem, Morton, K. W., II.; Baines, M. J., Eds., Oxford University Press: Oxford, UK, 1985; pp 177−197. (36) Cliffe, K. A.; Winters, K. H. The use of symmetry in bifurcation calculations and its application to the Bénard problem. J. Comput. Phys. 1986, 67, 310−326. (37) Seydel, R. Nonlinear computation. J. Franklin Inst. 1997, 334B (5/6), 1015−1047. (38) Seydel, R. In Interdisciplinary Applied Mathematics, Practical Bifurcation and Stability Analysis, From Equilibrium to Chaos, 2nd ed.; Springer: New York, 1994; Vol. 5. (39) Meijer, H. G. E.; Kálmar-Nagy, T. The Hopf-van der Pol System: Failure of a Homotopy Method. Differential Equations Dynam. Syst. 2012, 20 (3), 323−328. (40) Cayley, A. The Newton-Fourier Imaginary Problem. Am. J. Math. II 1879, 97. (41) Brolin, H. Invariant Sets under Iteration of Rational Functions. Arkiv Mat. 1966, 6, 103−144. (42) Kuno, M.; Seader, J. D. Computing all real solutions systems of nonlinear equations with a global fixed-point homotopy. Ind. Eng. Chem. Res. 1988, 27, 1320−1329. (43) Choi, S. H.; Book, N. L. Unreachable roots for global homotopy continuation methods. AIChE J. 1991, 37, 1093−1098. (44) Gustafson, J. B.; Beris, A. N. Evaluating all real roots of nonlinear equations using a global fixed-point homotopy method. AIChE J. 1991, 37, 1749−1752. (45) Jalali, F. Process simulation using continuation method in complex domain. Comput. Chem. Eng. 1998, 22 (1), 943−946. (46) Jalali, F.; Seader, J. D.; Khaleghi Rahimian, S. Global solution approaches in equilibrium and stability analysis using homotopy continuation in the complex domain. Comput. Chem. Eng. 2008, 32, 2333−2345. (47) Nikkhah-Bahrami, M.; Oftadeh, R. An effective iterative method for computing real and complex roots of systems of nonlinear equations. Appl. Math. Comput. 2009, 215, 1813−1820. (48) Oftadeh, R.; Nikkhah-Bahrami, M.; Najafi, A. A novel cubically convergent iterative method for computing complex roots of nonlinear equations. Appl. Math. Comput. 2010, 217, 2608−2618. (49) Khaleghi Rahimian, S.; Jalali, F.; Seader, J. D.; White, R. E. A new homotopy for seeking all real roots of a nonlinear equation. Comput. Chem. Eng. 2011, 35, 403−411. (50) Khaleghi Rahimian, S.; Jalali, F.; Seader, J. D.; White, R. E. A Robust Homotopy Continuation Method for Seeking All Real Roots of Unconstrained Systems of Nonlinear Algebraic and Transcendental Equations. Ind. Eng. Chem. Res. 2011, 50, 8892−8900. (51) Chow, S.; Mallet-Paret, J.; Yorke, J. A. Finding zeroes of maps: homotopy methods that are constructive with probability one. Math. Comput. 1978, 32, 887−899. (52) Watson, L. T.; Billups, S. C.; Morgan, A. P. ALGORITHM 652 HOMPACK 90: A Suite of Codes for Globally Convergent Homotopy Algorithms; GM Research Laboratories: Detroit, MI, 1987. (53) Rheinboldt, W. C.; Burkardt, J. V. A Locally Parameterized Continuation Process. ACM Transact. Math. Software 1983, 9, 2. (54) Dhooge, A.; Govaerts, W.; Kuznetsov, Yu. A. MatCont: A Matlab package for numerical bifurcation analysis of ODEs. ACM Transact. Math. Software 2003, 29 (2), 141−164. (55) Dhooge, A.; Govaerts, W.; Kuznetsov, Y. A.; Mestrom, W.; Riet, A. M.; Sautois, B. MATCONT and CL MATCONT: Continuation Toolboxes in MATLAB, http://www.matcont.ugent.be/manual.pdf, 2006. (Accessed September 9, 2013). (56) Chavez, C. R.; Seader, J. D.; Wayburn, T. L. Multiple steadystate solutions for interlinked separation systems. Ind. Eng. Chem. Fundamentals 1986, 25, 566−576.

(57) Gritton, K. S.; Seader, J. D.; Lin, J. W. Global homotopy continuation procedures for seeking all roots of a nonlinear equation. Comput. Chem. Eng. 2001, 25, 1003−1019. (58) Malinen, I.; Tanskanen, J. Modified bounded homotopies to enable a narrow bounding zone. Chem. Eng. Sci. 2008, 63, 3419−3430. (59) Malinen, I.; Tanskanen, J. Homotopy parameter bounding in increasing the robustness of homotopy continuation methods in multiplicity studies. Comput. Chem. Eng. 2010, 34, 1761−1774. (60) Rheinboldt, W. C. An adaptive continuation process for solving systems of nonlinear equations. Banchach Center Publications 1975, 3, 129−142. (61) Deuflhard, P. A Step Size Control for Continuation Methods and its Special Application to Multiple Shooting Techniques. Num. Math. Springer-Verlag: New York, 1979. (62) Lyness, J. N.; McHugh, B. J. J. Integration over multidimensional hypercubes I. A progressive procedure. Reprint of Classics; http://comjnl.oxfordjournals.org (July 9, 2013). (63) Yamamura, K. Simple Algorithms for Tracing Solution Curves. IEEE Trans. Circuits Systems-I: Fundam. Theory Appl. 1993, 40 (8), 537 −541. (64) Inoue, Y. A practical algorithm for DC operating-point analysis of large-scale circuits. Electron. Commun. Jpn. (Part III: Fundam. Electron. Sci.) 1994, 77 (10), 49−62. (65) Soler-López, M. A.; Barajas-Férnandez, J. (In spanish) Seguimiento de curvas homotópicas. Rev. Cienc. Básicas UJAT 2008, 7 (1), 26−34. (66) Vazquez-Leal, H.; Castaneda-Sheissa, R.; Rábago-Bernal, F.; Hernández-Martínez, L.; Sarmiento-Reyes, A.; Filobello-Niño, U. Powering Multiparameter Homotopy-Based Simulation with a Fast Path-Following Technique. International Scholarly Research Network, ISRN Appl. Math. 2011, No. 610637. (67) Vazquez-Leal, H.; Hernandez-Martinez, L.; Sarmiento-Reyes, A.; Castaneda-Sheissa, R.; Gallardo-Del-Angel, A. Homotopy method with a formal stop criterion applied to circuit simulation. IEICE Electron. Express 2011, 8 (21), 1808−1815. (68) Oliveros-Muñoz, J. M.; Jiménez-Islas, H. Hyperspherical path tracking methodology as correction step in homotopic continuation methods. Chem. Eng. Sci. 2013, 97, 413−429. (69) Wassestrom, E. Numerical Solutions by the Continuation Methods, SIAM Review. Col. 15. No. 1, 1973. (70) Georg, K. S. On Tracing an Implicitly Defined Curve by cuasiNewton Steps and Calculating Bifurcation by Local Perturbations. SIAM J. 1981, 2, 1. (71) Kangas, J.; Malinen, I.; Tanskanen, J. Modified Bounded Homotopies in the Solving of Phase Stability Problems for Liquid− Liquid Phase-Splitting Calculations. Ind. Eng. Chem. Res. 2011, 50 (11), 7003−7018. (72) Dunlavy, D.; Ó leary, D. P.; Thirumalai, D. HOPE: A Homotopy Optimization Method for Protein Structure Prediction. J. Comput. Biol. 2005, 12 (10), 1275−1288. (73) Yamamura, K. A Fixed-Point Homotopy Method for Solving Modified Nodal Equations. IEEE Transact. Circuits and Systems-I: Fundam. Theory Appl. 1999, 46, 6. (74) Wolf, D. M.; Sander, S. R. Multiparameter, Homotopy Methods for finding DC operating points of Nonlinear Circuits. IEEE Transact. Circuits Systems-I: Fundam. Theory Appl. 1996, 43 (10), 824−837. (75) Roychowdhury, J.; Melville, R. Delivering global DC convergence for large mixed-signal circuits via homotopy/continuation methods. IEEE Transact. Comput.-Aided Des. Integr. Circuits Syst. 2006, 25 (1), 66−78. (76) Paloschi, J. R. Bounded homotopies to solve systems of algebraic nonlinear equations. Comput. Chem. Eng. 1995, 19, 1243− 1254. (77) Paloschi, J. R. Bounded homotopies to solve systems of sparse algebra in nonlinear equations. Comput. Chem. Eng. 1997, 21, 531− 541. (78) Paloschi, J. R. Using sparse bounded homotopies in the SPEEDUP simulation package. Comput. Chem. Eng. 1998, 22, 1181− 1187. 14741

dx.doi.org/10.1021/ie402418e | Ind. Eng. Chem. Res. 2013, 52, 14729−14742

Industrial & Engineering Chemistry Research

Review

Inspired by Homotopy Continuation Methods. Hindawi Publishing Corporation. Math. Problems Eng. 2012, DOI: 10.1155/2012/309123. (103) Vázquez-Leal, H. Rational Homotopy Perturbation Method. Hindawi Publishing Corporation. Math. Problems Eng. 2012, DOI: 10.1155/2012/490342. (104) Salgovic, A.; Hlavacek, H.; Ilavskýa, J. Global simulation of countercurrent separation processes via one-parameter imbedding techniques. Chem. Eng. Sci. 1981, 36, 1599−1604. (105) El Boukili, A.; Marroco, A. Arclength continuation methods and applications to 2D drift-diffusion semiconductor equations. Institut National de Recherche en Informatique et en Automatique. Report de Recherhche No. 2546, 1995. (106) Dovi, V. G. Numerical solution of mass balance equations in continuous stirred reactors with radical intermediates. Chem. Eng. Sci. 1997, 53, 175−181. (107) Eckert, E.; Kubicek, M. Computing heterogeneous azeotropes in multicomponent mixtures. Comput. Chem. Eng. 1997, 21 (3), 347− 350. (108) Jalali, F.; Seader, J. D. Use of homotopy-continuation method in stability analysis of multiphase, reacting systems. Comput. Chem. Eng. 2000, 24, 1997−2008. (109) Deuflhard, P. Adaptive Pseudo-transient Continuation for Nonlinear Steady State Problems. Konrad-Zuse-Zentrum für Informationstechnik Berlin, ZIB-Report 02, 2002. (110) Mulas, M.; Tronci, S.; Baratti, R. Development of a 4measurable states activated sludge process model deduced from the asm1. 8th International IFAC Symposium on Dynamic Control of Process System, Preprints 1, Cancún, Mexico, 2007. (111) Liu, Y.; Han, B.; Dou, Y. X. A homotopy-projection method for the parameter estimation problems. Inverse Probl. Sci. Eng. 2008, 16 (2), 141−157.

(79) Zeleznik, F. J.; Gordon, S. Calculation of complex chemical equilibria. Ind. Eng. Chem. 1968, 60, 27. (80) Gordon, S.; McBride, B. J. Computer Program for Calculation of Complex Chemical Equilibrium Compositions, Rocket Performance, Incident and Reflected Shocks, and Chapman-Jouguet Detonations, NASA SP-273, 1971. (81) Gautam, R.; Seider, W. D. Calculation of phase and chemical equilibria, Part I, local and constrained minima in Gibbs free energy. AIChE J. 1979, 25. (82) Savageau, M. A. Finding Multiple Roots of Nonlinear Algebraic Equations Using S-System Methodology. Appl. Math. Comput. 1993, 55, 187−199. (83) Choi, S. H.; Harney, D. A.; Book, N. L. A robust path tracking algorithm for homotopy continuation. Comput. Chem. Eng. 1996, 20 (6−7), 647−655. (84) McDonald, C. M.; Floudas, C. A. GLOPEQ: A new computational tool for the phase and chemical equilibrium problem. Comput. Chem. Eng. 1997, 21 (1), 1−23. (85) Laiadi, D.; Hasseine, A.; Merzougui, A. Homotopy method to predict liquid-liquid equilibria for ternary mixtures of (water + carboxylic + organic solvent). Fluid Phase Equilib. 2012, 313, 114−120. (86) Dai, Y.; Kim, S.; Kojima, M. Computing all nonsingular solutions of cyclic-n polynomial using polyhedral homotopy continuation methods. J. Comput. Appl. Math. 2003, 152, 83−97. (87) Zwolak, J. W.; Tyson, J. J.; Watson, L. T. Finding all steady state solutions of chemical kinetic models. Nonlinear Analysis: Real World Appl. 2004, 5 (5), 801−814. (88) Liao, S. J. An approximate solution technique not depending on small parameters: a special example. Int. J. Nonlinear Mech. 1995, 303, 371−80. (89) Liao, S. J. Boundary element method for general nonlinear differential operators. Eng. Anal. Boundary Elements 1997, 202, 91−99. (90) Liao, S. J. Beyond Perturbation: Introduction to the Homotopy Analysis Method; Chapman & Hall/CRC Press: Boca Raton, FL, 2003. (91) Liao, S. J. On the Homotopy Analysis Method for nonlinear problems. Appl. Math. Comput. 2004, 47 (2), 499−513. (92) Liao, S. J.; Cheung, K. F. Homotopy analysis of nonlinear progressive waves in deep water. J. Eng. Math. 2003, 45 (2), 103−116. (93) He, J. H. Approximate solution for nonlinear differential equations with convolution product nonlinearities. Comput. Methods Appl. Mech. Eng. 1998, 167, 69−73. (94) Srivastava, V.; Rai, N.; Das, S. Analytical approach to micro scale bio-film heat transport using homotopy perturbation method. Int. J. Appl. Math. Comput. 2009, 1 (3), 148−158. (95) Hossein, J.; Firoozjaee, M. A. Application of Homotopy Analysis Method for Water Transport in Unsaturated Porous Media Studs. Nonlinear Sci. 2010, 1 (1), 08−13. (96) Sami-Bataineh, A.; Noorani, M. S. M.; Hashim, I. The homotopy analysis method for Cauchy reaction-diffusion problems. Phys. Lett. 2008, A 372, 613−618. (97) Shidfar, A.; Molabahrami, A. A weighted algorithm based on the homotopy analysis method: Application to inverse heat conduction problems. Commun. Nonlinear Sci. Num. Simul. 2010, 15 (10), 2908− 2915. (98) Domairry, G.; Bararnia, H. An approximation of the analytic solution of some nonlinear heat transfer equations: A Survey by using Homotopy analysis method. Adv. Stud. Theor. Phys. 2008, 2, 507−518. (99) Abbasbandy, S.; Shirzadi, A. MLPG method for two-dimensional diffusion equation with Neumann’s and non-classical boundary conditions. Appl. Num. Math. 2011, 61 (2), 170−180. (100) He, J. H. Comparison of homotopy perturbation method and homotopy analysis method. Appl. Math. Comput. 2004, 156 (2), 527− 539. (101) Malinen, I. Improving the robustness with modified bounded homotopies and problem-tailored solving procedures. Ph.D. Thesis; University of Oulu: Finland, 2010. (102) Vázquez-Leal, H.; Filobello-Niño, U.; Castañeda-Sheissa, R.; Hernández-Martínez, L.; Sarmiento-Reyes, A. Modified HPMs 14742

dx.doi.org/10.1021/ie402418e | Ind. Eng. Chem. Res. 2013, 52, 14729−14742