Homotopy parameter bounding in increasing the robustness of homotopy continuation methods in multiplicity studies

Computers and Chemical Engineering 34 (2010) 1761–1774

Contents lists available at ScienceDirect

Computers and Chemical Engineering journal homepage: www.elsevier.com/locate/compchemeng

Homotopy parameter bounding in increasing the robustness of homotopy continuation methods in multiplicity studies Ilkka Malinen, Juha Tanskanen ∗ Department of Process and Environmental Engineering, PO Box 4300, FI-90014 University of Oulu, Finland

a r t i c l e

i n f o

Article history: Received 20 May 2009 Received in revised form 5 February 2010 Accepted 24 March 2010 Available online 31 March 2010 Keywords: Homotopy continuation methods Bounded homotopies Path tracking Multiplicity studies Starting point multiplicity Distillation

a b s t r a c t Homotopy continuation methods are globally convergent methods, which can also be utilized in multiplicity studies. However, when the starting point and/or solution multiplicities lie on separate homotopy path branches, one or more of the solutions may be missed. This is due to the absence of real space connections between separate homotopy path branches, thus preventing multiple solutions being reached from a single starting point. In this paper, a concept is presented that enables a tracking starting point and solution multiplicities in cases where the standard problem-independent homotopy method fails. The concept is based on homotopy parameter bounding and enables the connection of separate homotopy path branches. The concept performance is examined using distillation column examples. In the examined cases the concept is found to improve robustness by establishing a path in real space such that solutions are approached that would be unattainable using the standard homotopy method. © 2010 Elsevier Ltd. All rights reserved.

1. Introduction In chemical engineering, process units are commonly described with mathematical expressions. Often, these models involve linear algebraic, nonlinear algebraic and transcendental equations. When these equations are collected together, a set of nonlinear equations f(x) = 0

(1)

is formed. This kind of set of equations is usually utilized to describe a steady state situation prevailing in the system, i.e. a situation that is independent of time. To find a solution for the equation set, a widely utilized strategy is to use Newton–Raphson based solving methods. These methods are based on the successive sequence of iterations according to the equation xk+1 = xk − ˛[f (xk )] f (xk )

−1

f(xk ),

(2) f(xk )

xk ,

is the Jacobian matrix of evaluated at and in which ˛ is a suitable parameter (0 < ˛ < 1) that is utilized to facilitate the convergence. Because the equation set may include thousands of equations, the size of the Jacobian matrix f(xk ) is often substantial. Frequently, Jacobian matrix structure is also sparse including a significant amount of zero elements. In order to reduce the number of function evaluations required in numerical Jacobian matrix

∗ Corresponding author. Tel.: +358 8 553 2340; fax: +358 8 553 2304. E-mail address: juha.tanskanen@oulu.fi (J. Tanskanen). 0098-1354/$ – see front matter © 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.compchemeng.2010.03.013

approximation and thus rationalize the computation, the Jacobian matrix sparsity pattern should be taken into consideration. The Newton–Raphson based solving methods are often feasible methods that may find a solution fast. The implementation of the Newton–Raphson based solving methods is also relatively simple and flexible. These methods require, however, that a favorable starting guess close enough to the root is offered. Therefore, the Newton–Raphson based methods (e.g. Newton’s method, quasiNewton methods) are called local or locally convergent methods. Introduction to the properties and use of the Newton–Raphson based solving methods can be found, e.g. in Bogle and Perkins (1988), Dennis and Schnabel (1983), Ortega and Rheinboldt (1970), Rabinowitz (1970) and Westerberg, Hutchison, Motard, and Winter (1979). Even though the convergence of the Newton–Raphson based solving methods can be enhanced with a good initialization by utilizing a succession of approximate models to provide guess values in the vicinity of the solution (Seider, Brengel, & Widagdo, 1991), the fundamental weakness in robustness of local solving methods cannot be fully overcome. In addition, local solving methods can at best find just one solution for each given starting guess. Since in chemical engineering, the model equations are often very nonlinear and process flowsheets frequently complex in character, multiple solutions are possible. Therefore, it is highly desirable to find all real space solutions for the equation set, especially those lying inside a feasible problem domain, or that no solutions exist with the specifications utilized. Therefore, globally convergent homotopy continuation methods have been adopted to increase robustness in

1762

I. Malinen, J. Tanskanen / Computers and Chemical Engineering 34 (2010) 1761–1774 2

Nomenclature b e f f g h M n x

domain boundary n × 1 vector where every element has the value one equation set Jacobian matrix of f auxiliary function in the homotopy equation homotopy function coefficient dimension variable vector

Greek letters ˛ parameter  homotopy parameter  penalty function ␷ auxiliary function in the bounded homotopy equation

2

into account the arc-length relation (dx1 /ds) + · · · + (dxn /ds) + 2 (d/ds) = 1 and an initial condition h(x0 , 0) = 0, an initial value problem is obtained. Path tracking based on the initial value problem is carried out with a predictor–corrector algorithm to obtain a solution of the original equation set. Different alternatives of combining the basic steps (predictor, parameterization, corrector, step control) may lead to a vast number of possible continuation methods. These issues have been discussed in an engineeringfriendly way, e.g. in Choi, Harney, and Book (1996), Kovach and Seider (1987), Seader et al. (1990), Seydel and Hlavacek (1987) and Wayburn and Seader (1987). A more mathematical discussion can be found, e.g. from Allgower and Georg (2003) and Seydel (1988). Depending on the choice of the function, g(x), the following homotopies, introduced, e.g. in Wayburn and Seader (1987), are obtained: Newton homotopy: h(x, ) = f(x) + (1 − )(f(x) − f(x0 ));

(4)

fixed-point homotopy: Subscripts b bounded i ith  homotopy parameter

h(x, ) = f(x) + (1 − )(x − x0 ); scale-invariant Affine homotopy: h(x, ) = f(x) + (1 − )f (x0 )(x − x0 ).

Superscripts b bounded inf infinity k iteration round max maximum min minimum  homotopy parameter 0 starting point * solution point  boundary

problem solving and to improve the ability to find multiple solutions. Homotopy continuation methods are sophisticated approaches that extend the domain of convergence, and thus, the solution may be achieved from a starting point that is unfavorable for locally convergent solving methods. Homotopy continuation methods have been introduced in an engineering-friendly way in several papers, e.g. in Kuno and Seader (1988), Seader et al. (1990), Seydel and Hlavacek (1987), and Wayburn and Seader (1987). The main idea in homotopy continuation methods is to gradually reach a root, x*, of the equation set, f(x), from a starting point, x0 , which satisfies another (simpler) system of equations, g(x). Both of these systems of equations can be embedded into a so-called homotopy function. Broadly defined, homotopy methods refer to the general form of homotopy equation h(x, ) = f (x) + (1 − )g(x) = 0,

(5)

(3)

which is convex with respect to the homotopy parameter. By gradually increasing the value of the homotopy parameter, , from zero to one, the weight of the auxiliary function, g(x), is diminished and the magnitude of actual nonlinear equation set, f(x), increased. The locus of solutions defines the homotopy path that is tracked with some continuation method. Consequently, homotopy continuation methods consist not only of the homotopy equation itself, but also the homotopy path tracking method, i.e. of some continuation strategy. Homotopy continuation methods are usually based upon differentiating the homotopy equation (Eq. (3)) with respect to the arc length s to obtain (∂h/∂x)(dx/ds) + (∂h/∂)(d/ds) = 0. By taking

(6)

As summarized in Seader et al. (1990), the Newton homotopy has the advantage that closed paths for x on  are possible within a finite domain. Because h(x, ) = 0 may have multiple roots for f(x) = f(x0 ) at  = 0, more than one homotopy path branch may exist. In that case, it cannot be guaranteed that all solutions will be obtained if the homotopy path is tracked from just one starting point. The fixed-point and scale-invariant Affine homotopies can be selected instead of the Newton homotopy to avoid the problem of starting point multiplicity. These methods have the advantage that h(x, ) = 0 is satisfied by only one root at  = 0. The fixed-point homotopy has a poor scaling property that originates from the connection of the dependent variables with the independent variables. Both the Affine and the Newton homotopies are, however, scale invariant. As already mentioned, the homotopy continuation methods are superior because of their convergence character compared with the traditional Newton–Raphson based solving methods. The homotopy methods are, however, time-consuming and require considerable computational effort. It is also notable that in spite of the global convergence character, the homotopy continuation methods cannot guarantee that the solutions can always be approached only based on computations in real space arithmetic. Because the Newton homotopy method may have starting point multiplicities, the homotopy path may return to other solutions of g(x) without passing through a solution of f(x) = 0. The existence of starting point multiplicities also increases the possibility for separate homotopy path branches. In this case, some or all of the solutions may be missed because there is not a homotopy path that could be tracked from a single starting point to the roots of the original problem f(x) = 0. In addition to the Newton homotopy, separate homotopy path branches may also exist with other homotopies. This means that even though fixed-point and Affine homotopies can be utilized to avoid the starting point multiplicity problem, the problem of multiple solutions on separate homotopy path branches still exist. Kuno and Seader (1988) found that all real roots could be detected on the same homotopy path starting from a single starting point provided that a criterion is utilized for a starting point selection. The criterion states that the starting point should minimize the number of real roots of a global fixed-point homotopy at an infinite value of the homotopy parameter. Later, Seader et al. (1990) pre-

I. Malinen, J. Tanskanen / Computers and Chemical Engineering 34 (2010) 1761–1774

sented mapping functions (toroidal and boomerang mappings) on x and/or , which can be utilized to transform the infinite path to a finite domain. By utilizing the mapped variables in path tracking, it was assumed that all the solutions of a general system of nonlinear equations were achieved with a fixed-point homotopy. Wayburn and Seader (1987) and Kuno and Seader (1988) provided examples of fixed-point homotopy paths for which all real roots are reachable through a complex space. In that case the real domain isola containing all real roots can be reached by tracking the homotopy path in a complex domain. Once the first root is reached, the remaining roots can be obtained by calculations in the real domain. The conjecture that all real roots are reachable from an arbitrary starting point if the fixed-point homotopy paths are tracked in a complex variable space was, however, disproved by Choi and Book (1991). They provided examples that demonstrate that some real roots are unreachable by tracking real or complex paths connected to the starting point. Gustafson and Beris (1991) showed, however, that by allowing the homotopy parameter also to take complex values the counterexample of Choi and Book (1991) can be solved. In this case, one of the variables, xi , is utilized as a continuation parameter and the rest of the variables as well as the homotopy parameter are allowed to be complex and/or infinite. According to Seader et al. (1990), there are two options for dealing with problems in a complex space. Firstly, the entire complex computation can be done directly in complex arithmetic, or secondly, the original equations can be put into real form (thus generating twice as many equations) and solved for both the real and imaginary parts of the roots, both being treated as real variables. Even though there are some references (e.g. Jalali, 1998; Jalali, Seader, & Khaleghi, 2008) of the utilization of complex arithmetic in chemical engineering problem solving, the complex arithmetic considerably complicates the solving task. Because flowsheeting problems in chemical engineering are often large in their size and relations for physical and chemical properties seldom tolerate complex variable values, the use of complex variable space can be concluded to be an unfeasible method for general use. The development of solving methods that are able to determine real space solutions without the use of complex variable space is therefore highly desirable. In order to reach this target, this paper proposes a concept based on homotopy parameter bounding to assist standard problem-independent homotopy continuation methods in the determination of multiple solutions that lie on separate homotopy path branches. The proposed concept also enlarges the area of convergence by offering a tool for determining starting point multiplicities for the Newton homotopy method.

2. Multiplicities in distillation In general, distillation characterizes well the modeling aspects of chemical engineering as a whole; the equation sets describing distillation systems are nonlinear, the equation sets are large and sparse, and interactions between distillation units may be tight. Distillation also illustrates the challenges encountered in multiplicity studies. It is an interesting detail that for a very long time separation systems were assumed to have only one unique steady state solution when the degrees of freedom for the system were satisfied. This was because multiple solutions had not been observed in separation systems (Chavez, Seader, & Wayburn, 1986). But as it is known nowadays, the existence of multiple steady states is indeed possible in separation systems, also in distillation. Multiple steady state (MSS) solutions as well as unstable operation points have been observed and reported to take place in distillation in several situations from the binary distillation (Jacobsen & Skogestad, 1991;

1763

Kienle, Groebel, & Gilles, 1995) to more complex types of homogeneous azeotropic distillation (Bekiaris, Meski, Radu, & Morari, 1993; Bekiaris, Meski, Radu, & Morari, 1994; Kannan, Joshi, Reddy, & Shah, 2005; Kienle, Gilles, & Marquardt, 1994), heterogeneous azeotropic distillation (Bekiaris, Meski, & Morari, 1996; Bekiaris, Güttinger, & Morari, 2000; Kovach & Seider, 1987), reactive distillation (Chen, Huss, Doherty, & Malone, 2002; Ciric & Miao, 1994; Güttinger & Morari, 1997, 1999a, 1999b; Kumar & Daoutidis, 1999; Singh, Singh, Kumar, & Kaistha, 2005), and interlinked column systems (Chavez et al., 1986; Lin, Seader, & Wayburn, 1987; Wayburn & Seader, 1984). In addition to the graphical and process simulation based studies, experimental investigation and verification of multiple steady states, as well as dynamic behavior studies, have been performed especially on homogeneous and heterogeneous azeotropic distillation (e.g. Dorn et al., 1998; Güttinger, Dorn, & Morari, 1997; Müller & Marquardt, 1997; Wang et al., 1997). Jacobsen and Skogestad (1991) were the first to present examples of steady state multiplicity in distillation columns with ideal thermodynamics. Two fundamentally different types of multiplicities were presented in their paper: 1a: Multiplicity in input transformations, which may appear especially when liquid streams are given on a mass or volume basis (instead of a molar basis). 1b: Multiplicity for molar inputs, where, in ideal distillation, the multiplicity may appear for specification of molar reflux and reboil rates. The third and fourth types of distillation multiplicities are presented, e.g. in Güttinger and Morari (1999b) and Bekiaris et al. (2000). The multiplicity of type (2a) is caused by the underlying reactive or non-reactive VL(L)E (non-ideal thermodynamics). It is essential that the condition of type (2a) multiplicity prevails in the entire column, i.e. the column is either non-reactive or non-hybrid reactive. In contrast, the fourth type of multiplicity (2b) caused by the underlying reactive and non-reactive VL(L)E occurs only in reactive hybrid column, i.e. in the column having both reactive and non-reactive sections. Based on the classification presented in Gani and Jørgensen (1994), multiplicity types can be defined as follows (Vadapalli & Seader, 2001): Output multiplicity refers to the case where all input variables (operating parameters) are specified and two or more sets of output variables are found. Respectively, input multiplicity refers to the case where one or more output variables are specified and multiple solutions are found for the unknown input variables (operating parameters). In addition to two main types of multiplicities, internal state multiplicity refers to the case, where multiple sets of internal conditions or profiles are found for the same values of the input and output variables. Bekiaris et al. (2000) defined the internal state multiplicity as a type of multiplicity where more than one solution with identical product compositions exists for the same inputs and operating parameters. 2.1. Significance of multiplicities on distillation design and operation The consideration of causes of multiple steady states is a relevant issue when considering distillation processes and operation strategies. Gani and Jørgensen (1994) showed, for example, that the number of solutions may be related to the number of column stages while the number of acceptable solutions depends on the set of specified variables. Therefore, an excess number of

1764

I. Malinen, J. Tanskanen / Computers and Chemical Engineering 34 (2010) 1761–1774

column stages (over-design) increases the possibility of multiplicity, but on the other hand, increases the flexibility of operation. As Vadapalli and Seader (2001) stated, multiplicities usually occur only for certain sets of specifications and only over restricted ranges of certain parameters. Often, one of the stable solutions represents the desired and most cost-effective operation point. However, from the operability point of view all the stable solutions need to be considered. Bekiaris and Morari (1996) discussed the various aspects how multiple steady states affect distillation design, synthesis and simulation, and thus for example on the selection of the entrainer, equipment, and separation scheme. Multiple steady states must also be taken into consideration when the controllability of distillation is considered; one design alternative may be much more difficult to control than another. In some cases multiplicities may even pose the risk of column profile jumping (or shifting) from a desirable steady state to an undesirable steady state. For design calculations (to size the towers properly), it is important to recognize if multiple solutions are possible, to locate them when they exist, to check for stability, and to perform dynamic simulations (Seider et al., 1991). The dynamics, operation, control, and stability issues of distillation systems with multiple solutions have been considered, e.g. in Jacobsen and Skogestad (1994), Kiva and Alukhanova (1997), Kumar and Daoutidis (1999) and Lee, Dorn, Meski, and Morari (1999). 2.2. Methods utilized to determine distillation multiplicities Since multiplicities influence in numerous ways the distillation simulation, design, operation, and control, it is important to find multiple steady states as robustly as possible. The detection of only one solution may result in misleading conclusions regarding the separation. Thus, some competitive alternatives may be missed and wrong conclusions drawn if the information of the existence of the multiple steady states is incomplete. Therefore, methods and procedures have been proposed during the years. The methods and procedures include both graphical and simulation approaches. 2.2.1. Graphical approach As Kiva and Alukhanova (1997) and Dorn et al. (1998) have pointed out, Petlyuk and Avetyan (1971) were the first who suggested the possible occurrence of multiple steady states via ∞/∞ analysis. Later, an intensive research and development work on the use of the ∞/∞ analysis in studying multiple steady states in ternary homogeneous azeotropic distillation (Bekiaris et al., 1993, 1994), heterogeneous azeotropic distillation (Bekiaris et al., 1996), quaternary systems (Bekiaris & Morari, 1996), homogeneous sequences consisting of two or more distillation columns, i.e. columns that form a separation train (Güttinger & Morari, 1996), and on reactive distillation (Güttinger & Morari, 1997, 1999a, 1999b) have been carried out. Bekiaris et al. (1993) derived a necessary and sufficient condition for the existence of multiple steady states based on the geometry of the distillation region boundaries. The ∞/∞ analysis is based on the limiting case of an infinite column length (infinite number of trays) and infinite reflux flow. One important feature is that very little information is needed (pure-component and azeotropic boiling points and compositions, and curvature of residue curve boundaries) to determine the feed region that leads to multiple steady states. Multiple steady states can be determined by constructing a bifurcation diagram, i.e. continuation of solutions, based on the graphical presentation using one of the product flow rates (e.g. distillate rate) as a bifurcation parameter. As Bekiaris and Morari (1996) have listed, the ∞/∞ analysis can be utilized to obtain a wide variety of information: the existence of multiplicities, the feed composition region that leads to multiplici-

ties, the feasible product regions, bifurcation diagram construction, the location of turning (limit) points, the multiplicity range for product flow rate, and the column composition profile. Since the method is graphical, all of this information can be obtained from residue curve maps involving the accurate curvature of residue curve boundaries. Because the orientation and the curvature of distillation boundaries depend on a specific thermodynamic model, switching from one VLE model to another may affect the existence of multiplicities quantitatively, and even qualitatively. It is, however, important to note that, as Bekiaris and Morari (1996) stated, the geometrical condition is only a sufficient (not necessary) condition for the existence of multiplicities for columns operating at finite reflux and having a finite number of trays. There may be mixtures that exhibit multiple steady states in finite columns but not in columns operating at the ∞/∞ conditions. 2.2.2. Simulation approach State-of-the-art commercial simulators cannot determine multiple steady states automatically. One strategy for determining multiple solutions is the trial and error approach, where an initial guess is varied to reveal possible multiple steady states. Convergence can be facilitated by generating proper estimates for composition and temperature profiles. This can be done for example based on the ∞/∞ analysis. The generated temperature and composition estimates (user supplied initial profile) may greatly facilitate the attainment of multiple steady state solutions, both stable and unstable, with (commercial) process simulation programs (e.g. Bekiaris & Morari, 1996; Kannan et al., 2005). In addition to the ∞/∞ analysis, dynamic simulation with stabilizing control has been utilized to overcome the shortcomings of conventional solvers especially in tracking unstable steady states (e.g. Kannan et al., 2005). Systematic methods and procedures are highly valued as they enable multiplicities to be found from only one initial point rather than the tedious strategy of performing several different trials that must be done when utilizing local solving methods. To avoid the utilization of numerous starting points in multiplicity studies, several strategies and solving tools have been proposed. Dalal and Malik (2003) introduced an optimization-based approach to track multiple solutions for a distillation column. However, algorithms based on homotopy continuation are more widely utilized (e.g. Chavez et al., 1986; Chen et al., 2002; Ciric & Miao, 1994; Han, Lee, & Yoon, 1993; Kovach & Seider, 1987; Lin et al., 1987; Singh et al., 2005; Wayburn & Seader, 1984). Problem-specific (problem-dependent or parametric continuation) homotopies have been utilized in the latest of these studies. Among others, liquid holdup on reactive trays, tray catalyst mass, column input (column operating conditions), reboil and reflux ratios, and Damköhler (Da) number have been applied as a continuation parameter. To enable multiplicity determination with commercial process simulators, Vadapalli and Seader (2001) presented a bifurcation technique using an arc-length continuation that can be incorporated as an add-in subroutine to a simulation program. The subroutine enables the automatic tracing of a solution path, including turning points, to obtain multiple solutions with respect to a user-selected parameter. The technique was illustrated with the ASPEN PLUS. Perhaps the most sophisticated option in the determination of multiple steady states is the so-called terrain methodology, introduced by Lucia and Feng (2002, 2003). Terrain methodology has been reported to be a very reliable, efficient, and global solving tool for multivariable process engineering simulation and optimization problems. Lucia and Feng (2004) demonstrated its applicability when solving steady state distillation examples with multiple solutions.

I. Malinen, J. Tanskanen / Computers and Chemical Engineering 34 (2010) 1761–1774

2.3. Problem-independent homotopy continuation methods in multiplicity studies It is impractical to carry out an exhaustive and systematic multiplicity study for every single simulation task separately. A more practical approach would be to improve the robustness, i.e. the reliability, of the solving methods in general. Thus, the capability of determining not only one real root, but also other real roots from an arbitrary starting point, would be increased. The problem-independent homotopy continuation methods are good candidates when developing the solving tools for determination of multiplicities. Unfortunately, none of them can guarantee that multiple steady state solutions can be determined from a single starting point by tracking the homotopy path just in real space. To facilitate the capability of problem-independent homotopy methods to determine multiple solutions lying on separate homotopy path branches, a homotopy path bounding concept, adopted outside the predefined homotopy parameter space [0 1] (see Eq. (3)), is proposed in the following chapters. The reason why the homotopy parameter space [0 1] has been selected is that besides  = 1 including the solutions of the original equation system and  = 0 including the starting point, the values between  = 0 and  = 1 are important in order to be able to track a continuous homotopy path from the starting point to the solution based only on computations in real space arithmetic. 3. Homotopy continuation methods with a bounded homotopy parameter In several cases the homotopy path runs outside the homotopy parameter space [0 1]. This may become problematic especially if the nonlinear equation set has several solutions. Even though the homotopy path runs through the roots, it may make long unnecessary curves outside the homotopy parameter space of interest. In some cases it may also appear that the homotopy path branches have no real space connections between them. In those cases the roots on different branches cannot be achieved from a single starting point with the standard homotopy continuation methods that rely only on computations in real space arithmetic. One alternative proposed here to prevent the unboundedness of the homotopy parameter is to utilize the concept of homotopy parameter bounding. In this way, the homotopy path can be forced to stay inside the predefined homotopy parameter space. The strategy extends the variety of the bounded homotopy methods originally proposed by Paloschi (1995, 1997, 1998) for restricting the homotopy path to stay inside the problem domain with respect to the variables of the problem, x. Later, Malinen and Tanskanen (2007a, 2008) introduced modifications, with the purpose of increasing the usability and robustness of bounded homotopies in chemical engineering problem solving. The modified bounded homotopies, together with variables mapping, force the homotopy path to be tightly bounded even when the bounding zone is narrow, and also enable flexible and accurate homotopy path tracking inside the narrow bounding zone. Even though the homotopy parameter-bounding strategy proposed here originates from work done on bounded homotopies, the basic idea of bounding the homotopy path with respect to the homotopy parameter has not been proposed previously. In general, the homotopies that include homotopy parameter bounding can be formulated as hb (x, ) = ()h(x, ) + ␷ () − ␷ ( b ) = 0.

(7)

In Eq. (7) the magnitude of the selected standard homotopy function h(x, ) is annihilated with the penalty function, (), whenever the homotopy path runs outside the predefined homo-

1765

topy parameter space [0 1]. Otherwise, the bounded homotopy method, hb (x, ), coincides with standard homotopy. Because the homotopy parameter, , exists in every equation of the standard homotopy function it is reasonable to focus both the annihilation effect and the effect of the auxiliary functions, ␷ () and ␷ ( b ), onto every (row) equation of the bounded homotopy method. When the equations in the original nonlinear equation set, f(x), are adequately scaled, i.e. the function values have the same order of magnitude, the auxiliary function as the form ␷ () = Me

(8)

has been found to be sufficient. In Eq. (8) the n × 1 vector, e, in which every element has the value one, is multiplied by the homotopy parameter and the proper coefficient, M. It is notable that the proper numeric value of the coefficient, M ∈ [−∞ +∞], depends to a certain extent on the problem that is being solved. When Eqs. (3), (7) and (8) are put together, the following expression is achieved hb (x, ) = ()(f (x) + (1 − )g(x)) + M( −  b )e = 0.

(9)

It is known that the Jacobian matrix structure (sparsity) of the standard Newton and Affine homotopy methods is the same as the Jacobian matrix structure of the original equation set. It is worth noting that this is also valid when considering the Jacobian matrix structure (sparsity) of the modified bounded homotopies. In the case of Newton and Affine homotopies, the Jacobian matrix of the bounded homotopies expressed as Eq. (9), have exactly the same zero elements, i.e. the same sparsity pattern, as the matrix f (x). In the case of problem-independent homotopy methods, the homotopy parameter is only an artificial parameter without any physical meaning. Thus, the bounding zone width and location may be freely specified. Since there is no need to define the bounding zone of the homotopy parameter to be narrow, the external boundaries are defined here to receive the values −1 and 2. The homotopy parameter bounding is forced to occur always when the homotopy path runs outside the homotopy parameter space [0 1]. The equations utilized to induce the homotopy parameter-bounding effect according to the desired specifications are expressed as follows. The equations shown here are based on the general equations presented in Paloschi (1995). In this study, the parameter,  b , expressed in Eqs. (7) and (9) has been defined as 







 b =−(6 · |− b |5 −15 · |− b |4 +10 · |− b |3 ) · [− b ],

(10)

and the penalty function, (), as () = 1 − (6 · | −  b |5 − 15 · | −  b |4 + 10 · | −  b |3 ).

(11)

b

The parameter,  , in Eq. (10) is defined as





b

=

0, , 1,

<0 0≤≤1 >1

(12)

With these assumptions, 0 ≤ () ≤ 1, () = 1 inside ˝ and  () = 0 outside ˝ . The parameter  b is such that  b =  inside ˝ ,   b   ∈ ı˝ for  outside ˝ , and  b ∈ ˝ for any . 4. Examples As the following test examples show, homotopy parameter bounding is a useful tool that can be utilized when connections between the solutions on the homotopy path branches are either formed, or broken. This enhances the probability of being able to determine multiple solutions, both starting point and solution multiplicities, from a single starting point without the need to utilize

1766

I. Malinen, J. Tanskanen / Computers and Chemical Engineering 34 (2010) 1761–1774

Fig. 1. The Newton homotopy paths for Eq. (13) with different coefficient M values when (a) the starting point is 1.25 (or 1.75) and (b) when the starting point is 0.75 (or 2.25). Starting point multiplicities () and multiple solutions (×).

several different starting points. In addition, the continuous homotopy path can be tracked in real variable space without the need for complex arithmetic. 4.1. Single function example Herein a simple and analytically solvable test function is utilized as an example to illustrate the performance of the homotopy parameter-bounding Eq. (9). The same second order polynomial equation written in the form f (x) = x2 − 3x + 2 = 0

(13)

has previously been studied, e.g. in Christiansen, Morud, and Skogestad (1996), Lin et al. (1987), Kuno and Seader (1988), and Seader et al. (1990). Eq. (13) has two roots: x* = 1 and x* = 2. The solving of the equation will be studied first with Newton homotopy and then with fixed-point homotopy. 4.1.1. Illustration with a Newton homotopy As illustrated in Fig. 1, when the standard Newton homotopy (Eq. (4)) is applied to Eq. (13), the homotopy path forms a parabola. The parabola opens either to the left or right depending on the selected starting point. Clearly, when the coefficient M in Eq. (9) is zero, i.e. in the case of the unbounded Newton homotopy method, the homotopy path runs through two points at  = 0 (through starting point multiplicities) and also through two actual solutions of Eq. (13) at  = 1 (through solution multiplicities). By changing the value of the coefficient M, different homotopy paths are formed. As shown in Fig. 1, the sign of coefficient M either strengthens or weakens the character of the locus of the standard (unbounded) homotopy path. When the coefficient M strengthens the character, the homotopy path to be tracked is shortened on the closed side of the parabola. At the same time the homotopy path branches at the open side of the parabola bend farther from each other, strengthening the unclosing character. Correspondingly, when the coefficient weakens the character of the unbounded homotopy path, it aims to open the closed side of parabola and close the open side of parabola. Thus, the opposite coefficient signs have a reverse effect on the locus of the bounded homotopy path. Fig. 2 illustrates the significance of the absolute numerical value of coefficient M on the form of the Newton homotopy path. When M = 0 (Fig. 2b), the unbounded Newton homotopy path forms a parabola, which opens on the right. If the negative coefficient M value is selected (Fig. 2a), the characteristic of further opening the open side of the parabola strengthens. Correspondingly, when the

positive coefficient M value is utilized, the open side of the parabola is closed. The characteristic of closing the open side of the parabola is strengthened when the value of coefficient M is raised (Fig. 2c–f). The situation is reverse on the closed side of the unbounded Newton homotopy path. A negative coefficient M value tends to shorten the closed arch of the homotopy path further (Fig. 2a). Correspondingly, as illustrated in Fig. 2c–f, raising the positive value of coefficient M strengthens the characteristic of opening the closed side of the parabola. At some ‘critical’ coefficient value a bifurcation point is reached, where two separate homotopy path branches touch each other (Fig. 2e). By further increasing the value of coefficient M the unclosing characteristic is strengthened (Fig. 2f) and the originally closed arch of the unbounded Newton homotopy path is opened. In addition to the numerical value of coefficient M, the starting point x0 also has an effect on the form of the bounded homotopy path. As represented in Fig. 3, starting point x0 has a strong effect on which coefficient M value is required to be able to close the open side and open the closed side of the unbounded Newton homotopy path. Fig. 3 also illustrates how the required coefficient value may vary considerably in both homotopy parameter-bounding zones. On the basis of Figs. 1–3, it can be concluded that:

1. The shape of the unbounded homotopy path depends on the original equation set, f(x), the standard homotopy method utilized, h(x, ), and the starting point, x0 . 2. In general, the effect of the coefficient M in Eq. (9) on the shape of the bounded homotopy path depends on the sign of the coefficient. 3. The absolute numerical value of the coefficient M has an effect on how strong the aim is either to form or break the homotopy path connections between multiple solutions.

4.1.2. Illustration with a fixed-point homotopy When the unbounded fixed-point homotopy method is applied to Eq. (13), several homotopy path branches are formed when x0 ≤ 2. When x0 < 1, three homotopy path branches exist (Fig. 4a). In this case the two roots of Eq. (13) are located on one branch, but the roots cannot be reached from the starting point based only on computations in real space arithmetic. When homotopy parameter-bounding Eq. (9) is applied with the positive coefficient M, two separate homotopy path branches can be connected in the homotopy parameter-bounding zone [−1 0]. However, the branch on which the actual roots exist remains unconnected.

I. Malinen, J. Tanskanen / Computers and Chemical Engineering 34 (2010) 1761–1774

1767

Fig. 2. The Newton homotopy paths for Eq. (13) with different coefficient M values. The starting point multiplicities () are 1.25 and 1.75, and multiple solutions (×) 1 and 2.

When 1 < x0 ≤ 2, three unbounded homotopy path branches exist. As illustrated in Fig. 4c and d, even though the starting point and one root lie on the same branch, the other root on another branch cannot be approached from the selected start-

ing point. This is because no connection exists between the two branches. By adopting homotopy parameter bounding for Eq. (9) with a positive coefficient M, the three separate homotopy path branches can be connected in such a way that both

Fig. 3. Significance of coefficient M on the character of a Newton homotopy path when applied to Eq. (13) with different starting point x0 values (a) in the case  < 0 and (b) in the case  > 1. ‘Open’ refers to the open parabola and ‘closed’ to the closed parabola.

1768

I. Malinen, J. Tanskanen / Computers and Chemical Engineering 34 (2010) 1761–1774

Fig. 4. The unbounded (· · ·) and bounded (–) fixed-point homotopy paths for Eq. (13) with various starting point x0 and coefficient M values. The starting point () and multiple solutions (×).

roots of Eq. (13) are approached from the selected starting point. When x0 = 1, both roots can be approached from the starting point (Fig. 4b). In order to achieve this, the bifurcation point at (x, ) = (1, 0.5) must be determined where the homotopy branches intersect. By tracking the branches that fork at the bifurcation point, the root x* = 2 can be approached. As illustrated in Fig. 4b, by adopting homotopy parameter bounding in Eq. (9) both roots for Eq. (13) can be achieved without the need for intensive bifurcation analysis.

When x0 > 2, a single unbounded homotopy path branch is obtained that runs through the starting point and multiple roots (Fig. 4e). Thus, in this case homotopy parameter-bounding Eq. (9) does not bring any benefit to the solving of the roots of Eq. (13). Herein, the homotopy path branches for the bounded fixed-point homotopy method with the coefficient M = −20 have been represented for illustrative purposes only. Fig. 5 has been applied in order to determine a coefficient value that is able to break the arches in both homotopy parameter-bounding zones, i.e. [−1 0] and [1 2].

I. Malinen, J. Tanskanen / Computers and Chemical Engineering 34 (2010) 1761–1774

1769

Fig. 5. Significance of coefficient M on the character of the fixed-point homotopy path when applied to Eq. (13) with different starting point x0 values (a) in the case  < 0 and (b) in the case  > 1. ‘Open’ refers to the unconnected homotopy path branches and ‘closed’ to the connected homotopy path branches.

4.2. Equation set example The equation set consisting of two nonlinear equations 4x13 + 4x1 x2 + 2x22 − 42x1 − 14 = 0

(14a)

4x23 + 4x1 x2 + 2x12 − 26x2 − 22 = 0

(14b)

has been studied earlier, e.g. in Dalal and Malik (2003). The equation set has altogether nine solutions, which are situated in −4 ≤ x1 ≤ 4 and −4 ≤ x2 ≤ 4. The solutions are shown in Table 1. As the demonstrations with both the Newton and fixed-point homotopies illustrate, the proposed strategy of bounding the homotopy path to run inside the limited homotopy parameter space enables approaching multiplicity solutions, which are not achievable from a single starting point with either the standard homotopies nor the original bounded homotopies proposed by Paloschi. Even though the homotopy path bounding against the artificial homotopy parameter (instead of the problem variables) makes the homotopy path unbounded against the problem variables, this does not necessarily prevent the successful solving of a multiplicity problem. Everything depends on the course of the homotopy path, i.e. whether the path runs through unfeasible problem variable values or not before the multiplicity solutions have been tracked. 4.2.1. Illustration with a Newton homotopy As illustrated in Fig. 6, when the starting point is (−2, −2), a standard Newton homotopy does not allow the approaching of all of the solutions. Four solutions lie on an isola that neither runs through the starting point nor the plane  = 0. When the homotopy parameter-bounding Eq. (9) is adopted with M = 100, all the nine solutions (including the isola solutions) are approached from the selected starting point by tracking the continuous homotopy path in the real space. 4.2.2. Illustration with a fixed-point homotopy As shown in Fig. 7, homotopy parameter bounding is also useful when the fixed-point homotopy method is utilized in solving the equation set example. When the starting point is (−2, −2), there are altogether five separate homotopy path branches, which have no real space connections between them. Only one solution can be

Fig. 6. Newton homotopy paths for the equation set example. The unbounded Newton homotopy paths (—) and the bounded Newton homotopy path (—). The coefficient value utilized in bounding was 100. Starting point (x1 , x2 )0 = (−2, −2) (䊉), starting point multiplicities () and solution multiplicities (×).

found on the same homotopy path branch with the selected starting point. The other eight solutions lie on four separate branches, which each run through two solutions and towards the positive homotopy parameter infinity. When the homotopy parameter-bounding Eq. (9) with M = −100 is adopted, the separate homotopy path branches are connected and all nine solutions can be approached from a single starting point on the same continuous homotopy path. 5. Distillation column examples Two distillation column examples are considered here, illustrating how homotopy parameter bounding may increase the robustness of the standard homotopy continuation methods when solving chemical engineering problems. The examples illustrate how the homotopy parameter bounding can be utilized to determine nontrivial starting points and solution multiplicities. The Newton homotopy method is utilized as a standard homotopy method in both example cases. In the cases, exact mole fraction specifications are set for the product flows. In addition, the full set of MESH equations and non-ideal thermodynamics are utilized. Since the dimension of the enthalpy equations is 1 J s−1 , they are multiplied by the value 1 × 10−6 to keep the values of the MESH equations

Table 1 Solutions of the equation set example.

x1 x2

1

2

3

4

5

6

7

8

9

−3.7793 −3.2832

−3.0730 −0.0814

−2.8051 3.1313

−0.2708 −0.9230

−0.1280 −1.9537

0.0867 2.8843

3.0000 2.0000

3.3852 0.0739

3.5844 −1.8481

1770

I. Malinen, J. Tanskanen / Computers and Chemical Engineering 34 (2010) 1761–1774

Table 2 Thermodynamic parameters for acetone/chloroform/benzene mixture. Parameters are from an HYSYS database. Binary interaction parameters for the Wilson equation (cal mol−1 ) A11 = 0 A21 = −506.852 A31 = −243.965 Component

A12 = 116.117 A22 = 0 A32 = −11.823 A

B

Constants for Antoine equation and molar volumes (dm3 kmol−1 ) Acetone 71.3031 −5952 Chloroform 73.7058 −6055.6 Benzene 169.65 −10314.8

A13 = 682.406 A23 = −71.811 A33 = 0 C

D

E

F

Molar volumes

0 0 0

−8.53128 −8.9189 −23.5895

7.82393e−6 7.74407e−6 2.09442e−5

2 2 2

74.4768 80.7314 89.5551

ln(Psat ) = A + B/(T + C) + D ln(T) + E(TF ); P is pressure (kPa) and T is temperature (K).

roughly on the same scale. Thermodynamic models and parameters are based on the Wilson thermodynamic package and on the database of the commercial simulation program HYSYS. The acetone/chloroform/benzene component set is studied in both column cases. The thermodynamic data used is shown in Table 2. The column models and problem-independent equationoriented solving approach have been implemented in a MATLAB environment. The solving method is based on Newton homotopy and arc-length path tracking strategy (continuation method). The pressure drop throughout the column is ignored and the atmospheric pressure is assumed in both distillation column cases. Liquid and vapor flows of 1 mol s−1 in every column stage, reboil and condenser duties of 20 kW, reflux and reboil ratios of 1, and mole fraction and temperature values of the column feed in every stage of the column have been set as the trivial starting point for the homotopy path tracking. Because the starting point selection is straightforward and no special attention is paid to generating a favorable one, standard homotopy methods may fail in cases where multiplicities exist. The mapped variables are used here instead of unmapped ones to facilitate homotopy path tracking. Variables mapping is use-

Fig. 7. Fixed-point homotopy paths for an equation set example. (a) Unbounded fixed-point homotopy paths and (b) bounded fixed-point homotopy path. The coefficient M value utilized in bounding was −100. Starting point (x1 , x2 )0 = (−2, −2) (䊉) and solution multiplicities (×).

Fig. 8. The liquid sidestream column configuration with the exact mole fraction specifications studied in Section 5.1.

I. Malinen, J. Tanskanen / Computers and Chemical Engineering 34 (2010) 1761–1774

1771

Fig. 9. The Newton homotopy paths for a liquid sidestream with and without homotopy parameter bounding. The starting point (), starting point multiplicities (䊉), and solution (). Table 3 The selected variables characterizing the actual solution of the case studied in Section 5.1. Process variable

Value

Composition of distillatea Composition of sidestreama Composition of bottom producta Reflux ratio Reboil ratio Reboiler duty (kJ mol−1 feed) Condenser duty (kJ mol−1 feed)

0.8000/0.1882/0.0118 0.3883/0.5500/0.0617 0.0078/0.1922/0.8000 13.928 4.440 101.06 100.68

The mole fraction values in bold describe the specifications utilized in solving. a Compositions in mole fractions (acetone/chloroform/benzene).

ful especially when the homotopy path runs close to the problem domain boundaries. In distillation, this regularly arises with mole fractions when mole fraction values close to 0 and 1 are encountered. With proper variables mapping, homotopy path tracking

Fig. 11. The simple distillation column configuration with the exact mole fraction specifications studied in Section 5.2.

becomes more flexible and unnecessary homotopy step predictions outside the meaningful variable space are avoided. The homotopy path is tracked in a mapped infinite variable space, but the variables are still mapped into a finite variable space before substituting them into the equations of the nonlinear equation set, f(x). The variables mapping is carried out as introduced in Malinen and Tanskanen (2007a, 2008) based on the following equations: Variable xi mapping from a finite space into an infinite space:



xiinf = log10

 Fig. 10. Liquid mole fraction profile fulfilling the state distribution specifications illustrated in Fig. 8.

xiinf = log10

2 · (xi − bmin ) i bmax − bmin i i



,

0.5 · (bmax − bmin ) i i bmax − xi i

when xi < 0.5 · (bmax + bmin ) i i

(15)

 ,

when xi ≥ 0.5 · (bmax + bmin ) i i

(16)

1772

I. Malinen, J. Tanskanen / Computers and Chemical Engineering 34 (2010) 1761–1774

Fig. 12. Mole fraction of acetone in the bottom product as a function of the reboil ratio. The acetone mole fraction in the distillate is kept at a constant 0.95.

Correspondingly, variable xiinf mapping from an infinite space into a finite space: xi = bmin + 0.5 · (bmax − bmin ) · 10 i i i xi =

bmax i

− 0.5 ·

(bmax − bmin ) i i xinf

10

i

,

xinf i

,

when

when xiinf < 0 xiinf

≥ 0.

(17) (18)

To enable utilization of the mapping equations, the maximum bmax and minimum bmin values of the variables must be specified. i i In the distillation column examples studied, every variable in the vector, x, has been mapped. 5.1. Starting point multiplicity of the liquid sidestream column The starting point multiplicity discovered in Malinen and Tanskanen (2007b) is studied further in this case. The liquid sidestream column under consideration consists of 28 column stages, a condenser, and a reboiler. The boiling acetone/chloroform/benzene liquid feed is led to stage 20 and

the liquid sidestream is drawn from stage 10 according to Fig. 8. As shown in Fig. 9, when the standard Newton homotopy method is used, the starting point multiplicity is encountered at  = 0. In addition, the selected starting point lies on a separate homotopy path branch to the other two points of the starting point multiplicity. When the homotopy path is tracked to the positive direction from the selected starting point, the homotopy path does not run through the problem solution, but runs outside the problem domain. The same happens when the homotopy path is tracked to a negative direction. Therefore, it is not possible to track a continuous homotopy path from the selected starting point to the solution. By utilizing a homotopy parameter bounding equipped with a proper negative coefficient M value, it is possible to force the separate homotopy path branches to be connected. In this case, the auxiliary function part in Eq. (9) bends the Newton homotopy path so it becomes possible to track a continuous homotopy path from the trivial starting point to the actual solution. It is notable that the absolute value of the coefficient M must be high enough to enable the connection to form between the two branches on where the selected starting point and the actual solution lie. If the absolute value of the coefficient M is too low, the auxiliary function part in Eq. (9) cannot force the two homotopy path branches to be connected, in which case the solution cannot be achieved. The actual solution for the problem is illustrated in Table 3 and the liquid mole fraction profile fulfilling the mole fraction specifications set for product flow purities is shown in Fig. 10.

5.2. Input multiplicity of a simple distillation column The simple distillation column configuration shown in Fig. 11 is studied in this case. The case resembles the multiplicity illustrated earlier in Tanskanen and Malinen (2005), where the effect of the distillation boundary on the mole fractions of the product flows through the reflux and reboil ratio values was illustrated. The column studied here has 98 stages, a condenser, and a reboiler. The boiling acetone/chloroform/benzene liquid feed is routed onto stage 50. As illustrated in Fig. 12, when the mole fraction of acetone in the distillate is maintained as 0.95, and the reboil ratio is increased, the mole fraction of acetone in the bottom product flow first decreases, and after achieving the turning point, starts to increase. As is seen, a multiplicity is formed. For example, an acetone mole fraction value of 0.054 in the bottom product flow is achieved with two different reboil ratio values.

Fig. 13. The Newton homotopy paths for the reboil ratio with and without homotopy parameter bounding. The reboil ratio is 1 at the starting point (䊉) and 5.197 and 9.742 in multiple solutions #1 and #2 ().

I. Malinen, J. Tanskanen / Computers and Chemical Engineering 34 (2010) 1761–1774 Table 4 The selected variables characterizing two solutions of the case studied in Section 5.2. Process variable

Solution #1

Solution #2

Composition of distillatea Composition of bottom producta Reflux ratio Reboil ratio Reboiler duty (kJ mol−1 feed) Condenser duty (kJ mol−1 feed)

0.9500/0.0401/0.0099 0.0540/0.3585/0.5875

0.9500/0.0496/0.0004 0.0540/0.3491/0.5969

4.403 5.197 78.67

9.146 9.742 147.54

78.55

147.41

The mole fraction values in bold describe the specifications utilized in solving. a Compositions in mole fractions (acetone/chloroform/benzene).

Fig. 14. Two liquid mole fraction profiles fulfilling the mole fraction specifications of Fig. 11. Mole fraction profile #1 (䊉) and #2 ().

As illustrated in Fig. 13, the homotopy parameter bounding with the proper coefficient M forces Newton homotopy path bending, so that the path runs from solution #1 to solution #2 without crossing the problem domain boundaries. The coefficient M value has an influence on how strongly the homotopy path is bent. If the absolute value of coefficient M is too small, the bounding property is not achieved and the homotopy path crosses the problem domain. However, if the absolute value of coefficient M is increased, the desired kind of bounding property is achieved and both solutions shown in Table 4 can be tracked. As shown in Fig. 14, two different state distributions fulfill the mole fraction specifications that have been set for the product flows of the studied column configuration. 6. Conclusions Solving methods based on homotopy continuation are globally convergent. In addition, these methods can be utilized to track multiple solutions for a nonlinear equation set. There are, however, also some deficiencies, which decrease the robustness of the problemindependent homotopy continuation methods in tracking multiple solutions. For example, in the case of starting point multiplicity of a Newton homotopy there may not exist a homotopy path in real space from the selected starting point to the actual solution of a nonlinear equation set. From the point of view of multiplicity study, when multiple starting points exist, it is highly recommended that every branch running through the different starting points should be tracked.

1773

In addition, in the case of solution multiplicity, the different homotopy path branches running through the different roots are problematic. When the branches do not have real space connections between them, it is a challenge to obtain the roots that lie on the separate homotopy path branches from the starting point. Therefore, when real space arithmetic is favored (as is often the situation in chemical engineering), a method capable of forcing the homotopy path branches to be connected would be very welcome. In this paper, a concept has been proposed that is capable of forcing separate homotopy path branches to be connected. The proposed concept is based on the bounding of the homotopy parameter. The concept is able to connect separate homotopy path branches, assuming that a proper value for the auxiliary function coefficient M is defined that forces the homotopy path to bend strongly enough outside the homotopy parameter space [0 1]. Correspondingly, the closed homotopy path arcs between multiple solutions can be disconnected with the appropriate auxiliary function coefficient. The performance of the proposed concept was illustrated with distillation column cases, where it was found to improve the robustness of the Newton homotopy method both in starting point and solution multiplicity cases. In the cases examined, the concept enabled approaching solutions that are unattainable with standard problem-independent homotopy continuation method, based only on computations in real space arithmetic. So far, no systematic strategy has been proposed for choosing an appropriate coefficient value M for an auxiliary function to manage robustly in multiplicity cases. It is clearly crucial to formulate the findings introduced in this paper as a reliable algorithm, possessing the capability of adaptively determining appropriate auxiliary function coefficient values at specific points of the homotopy path, i.e. at  = 0 and  = 1. At these particular points, first and second order derivatives may be essential in obtaining a figure of the behavior of the function values relative to the variables, x, and the homotopy parameter, . The concept introduced in this paper is not sufficient when solving isola problems, where no real space homotopy path exists from  = 0 to any of the solutions at  = 1. In fact, solving this kind of isola problem relying only on computations in real space arithmetic is in general an open question that needs to be tackled. The concept introduced in this paper could be one important element on the way to developing a systematic problem-independent method for multiplicity examinations. Acknowledgements The postgraduate program Graduate School in Chemical Engineering (GSCE) and the financial support from Tauno Tönning and Emil Aaltonen Foundations are gratefully acknowledged. The helpful comments of Jani Kangas, University of Oulu, are highly appreciated. References Allgower, E. L., & Georg, K. (2003). Introduction to numerical continuation methods. Philadelphia: SIAM. Bekiaris, N., Güttinger, T. E., & Morari, M. (2000). Multiple steady states in distillation: Effect of VL(L)E inaccuracies. AIChE Journal, 46, 955. Bekiaris, N., Meski, G. A., Radu, C. M., & Morari, M. (1993). Multiple steady states in homogeneous azeotropic distillation. Industrial and Engineering Chemistry Research, 32, 2023. Bekiaris, N., Meski, G. A., Radu, C. M., & Morari, M. (1994). Design and control of homogeneous azeotropic distillation columns. Computers and Chemical Engineering, 18, S15. Bekiaris, N., Meski, G. A., & Morari, M. (1996). Multiple steady states in heterogeneous azeotropic distillation. Industrial and Engineering Chemistry Research, 35, 207. Bekiaris, N., & Morari, M. (1996). Multiple steady states in distillation: ∞/∞ predictions, extensions, and implications for design, synthesis, and simulation. Industrial and Engineering Chemistry Research, 35, 4264.

1774

I. Malinen, J. Tanskanen / Computers and Chemical Engineering 34 (2010) 1761–1774

Bogle, I. D. L., & Perkins, J. D. (1988). Sparse Newton-like methods in equation oriented flowsheeting. Computers and Chemical Engineering, 12, 791. Chavez, C. R., Seader, J. D., & Wayburn, T. L. (1986). Multiple steady-state solutions for interlinked separation systems. Industrial and Engineering Chemistry Fundamentals, 25, 566. Chen, F., Huss, R. S., Doherty, M. F., & Malone, M. F. (2002). Multiple steady states in reactive distillation: Kinetic effects. Computers and Chemical Engineering, 26, 81. Choi, S. H., & Book, N. L. (1991). Unreachable roots for global homotopy continuation methods. AIChE Journal, 37, 1093. Choi, S. H., Harney, D. A., & Book, N. L. (1996). A robust path tracking algorithm for homotopy continuation. Computers and Chemical Engineering, 20, 647. Christiansen, A. C., Morud, J., & Skogestad, S. (1996). A comparative analysis of methods for solving systems of nonlinear algebraic equations. In Proceedings of the 38th SIMS simulation conference Trondheim, (pp. 217–230). Ciric, A. R., & Miao, P. (1994). Steady state multiplicities in an ethylene glycol reactive distillation column. Industrial and Engineering Chemistry Research, 33, 2738. Dalal, N. M., & Malik, R. K. (2003). Solution multiplicity in multicomponent distillation. A computational study. In A. Kraslawski, & I. Turunen (Eds.), Proceedings of ESCAPE13 (pp. 617–622). Elsevier. Dennis, J. E., & Schnabel, R. B. (1983). Numerical methods for unconstrained optimization and nonlinear equations. Englewood Cliffs, NJ: Prentice Hall. Dorn, C., Güttinger, T. E., Wells, G. J., Morari, M., Kienle, A., Klein, E., et al. (1998). Stabilization of an unstable distillation column. Industrial and Engineering Chemistry Research, 37, 506. Gani, R., & Jørgensen, S. B. (1994). Multiplicity in numerical solution of non-linear models: Separation processes. Computers and Chemical Engineering, 18, S55. Gustafson, J. B., & Beris, A. N. (1991). Evaluating all real roots of nonlinear equations using a global fixed-point homotopy method. AIChE Journal, 37, 1749. Güttinger, T. E., Dorn, C., & Morari, M. (1997). Experimental study of multiple steady states in homogeneous azeotropic distillation. Industrial and Engineering Chemistry Research, 36, 794. Güttinger, T. E., & Morari, M. (1996). Multiple steady states in homogeneous separation sequences. Industrial and Engineering Chemistry Research, 35, 4597. Güttinger, T. E., & Morari, M. (1997). Predicting multiple steady states in distillation: Singularity analysis and reactive systems. Computers and Chemical Engineering, 21(SS), S995. Güttinger, T. E., & Morari, M. (1999a). Predicting multiple steady states in equilibrium reactive distillation. 1. Analysis of nonhybrid systems. Industrial and Engineering Chemistry Research, 38, 1633. Güttinger, T. E., & Morari, M. (1999b). Predicting multiple steady states in equilibrium reactive distillation. 2. Analysis of hybrid systems. Industrial and Engineering Chemistry Research, 38, 1649. Han, K. T., Lee, K. J., & Yoon, E. S. (1993). Simulation and analysis of complex distillation columns using continuation methods. Computers and Chemical Engineering, 17, S479. Jacobsen, E. W., & Skogestad, S. (1991). Multiple steady states in ideal two-product distillation. AIChE Journal, 37, 499. Jacobsen, E. W., & Skogestad, S. (1994). Instability of distillation columns. AIChE Journal, 40, 1466. Jalali, F. (1998). Process simulation using continuation method in complex domain. Computers and Chemical Engineering, 22, S943. Jalali, F., Seader, J. D., & Khaleghi, S. (2008). Global solution approaches in equilibrium and stability analysis using homotopy continuation in the complex domain. Computers and Chemical Engineering, 32, 2333. Kannan, A., Joshi, M. R., Reddy, G. R., & Shah, D. M. (2005). Multiple-steady-states identification in homogeneous azeotropic distillation using a process simulator. Industrial and Engineering Chemistry Research, 44, 4386. Kienle, A., Gilles, E. D., & Marquardt, W. (1994). Computing multiple steady states in homogeneous azeotropic distillation processes. Computers and Chemical Engineering, 18, S37. Kienle, A., Groebel, M., & Gilles, E. D. (1995). Multiple steady states in binary distillation—Theoretical and experimental results. Chemical Engineering Science, 50, 2691. Kiva, V. N., & Alukhanova, B. M. (1997). Multiple steady states of distillation and its realisation. Computers and Chemical Engineering, 21(SS), S541. Kovach, J. W., & Seider, W. D. (1987). Heterogeneous azeotropic distillation— Homotopy-continuation methods. Computers and Chemical Engineering, 11, 593.

Kumar, A., & Daoutidis, P. (1999). Modeling, analysis and control of ethylene glycol reactive distillation column. AIChE Journal, 45, 51. Kuno, M., & Seader, J. D. (1988). Computing all real solutions to systems of nonlinear equations with a global fixed-point homotopy. Industrial and Engineering Chemistry Research, 27, 1320. Lee, M., Dorn, C., Meski, G. A., & Morari, M. (1999). Limit cycles in homogeneous azeotropic distillation. Industrial and Engineering Chemistry Research, 38, 2021. Lin, W.-J., Seader, J. D., & Wayburn, T. L. (1987). Computing multiple solutions to systems of interlinked separation columns. AIChE Journal, 33, 886. Lucia, A., & Feng, Y. (2002). Global terrain methods. Computers and Chemical Engineering, 26, 529. Lucia, A., & Feng, Y. (2003). Multivariable terrain methods. AIChE Journal, 49, 2553. Lucia, A., & Feng, Y. (2004). Solving distillation problems by terrain methods. Computers and Chemical Engineering, 28, 2541. Malinen, I., & Tanskanen, J. (2007a). A rigorous minimum energy calculation method for a fully thermally coupled distillation system. Chemical Engineering Research and Design, 85, 502. Malinen, I., & Tanskanen, J. (2007b). Modified bounded Newton homotopy method in solving sidestream column configurations. In V. Ples¸u, & P. S¸. Agachi (Eds.), Proceeding of ESCAPE17. Elsevier [CD-ROM]. Malinen, I., & Tanskanen, J. (2008). Modified bounded homotopies to enable a narrow bounding zone. Chemical Engineering Science, 63, 3419. Müller, D., & Marquardt, W. (1997). Experimental verification of multiple steady states in heterogeneous azeotropic distillation. Industrial and Engineering Chemistry Research, 36, 5410. Ortega, J. M., & Rheinboldt, W. C. (1970). Iterative solution of nonlinear equations in several variables. New York: Academic Press. Paloschi, J. R. (1995). Bounded homotopies to solve systems of algebraic nonlinear equations. Computers and Chemical Engineering, 19, 1243. Paloschi, J. R. (1997). Bounded homotopies to solve systems of sparse algebraic nonlinear equations. Computers and Chemical Engineering, 21, 531. Paloschi, J. R. (1998). Using sparse bounded homotopies in the SPEEDUP simulation package. Computers and Chemical Engineering, 22, 1181. Petlyuk, F. B., & Avetyan, V. S. (1971). Investigation of three component distillation at infinite reflux. Theoretical Foundations of Chemical Engineering, 5, 499 [in Russian]. Rabinowitz, P. (1970). Numerical methods for nonlinear algebraic equations. London: Gordon and Breach Science Publishers. Seader, J. D., Kuno, M., Lin, W.-J., Johnson, S. A., Unsworth, K., & Wiskin, J. W. (1990). Mapped continuation methods for computing all solutions to general systems of nonlinear equations. Computers and Chemical Engineering, 14, 71. Seider, W. D., Brengel, D. D., & Widagdo, S. (1991). Nonlinear analysis in process design. AIChE Journal, 37, 1. Seydel, R. (1988). From equilibrium to chaos. Practical bifurcation and stability analysis. New York: Elsevier. Seydel, R., & Hlavacek, V. (1987). Role of continuation in engineering analysis. Chemical Engineering Science, 42, 1281. Singh, B. P., Singh, R., Kumar, M. V. P., & Kaistha, N. (2005). Steady state analysis of reactive distillation using homotopy continuation. Chemical Engineering Research and Design, 83, 959. Tanskanen, J., & Malinen, I. (2005). A rigorous calculation method for the minimum reflux of nonideal homogeneous multicomponent distillation. In S. Pierucci (Ed.), Chemical engineering transactions (pp. 281–286). Vadapalli, A., & Seader, J. D. (2001). A generalized framework for computing bifurcation diagrams using process simulation programs. Computers and Chemical Engineering, 25, 445. Wang, C. J., Wong, D. S. H., Chien, I.-L., Shih, R. F., Wang, S. J., & Tsai, C. S. (1997). Experimental investigation of multiple steady states and parametric sensitivity in azeotropic distillation. Computers and Chemical Engineering, 21(SS), S535. Wayburn, T. L., & Seader, J. D. (1984). Solution of systems of interlinked distillation columns by differential homotopy-continuations methods. In A. W. Westerberg, & H. H. Chien (Eds.), Proceedings of the Second International Conference on Foundations of Computer-Aided Process Design (pp. 765–862) [CACHE]. Wayburn, T. L., & Seader, J. D. (1987). Homotopy continuation methods for computeraided process design. Computers and Chemical Engineering, 11, 7–25. Westerberg, A. W., Hutchison, H. P., Motard, R. L., & Winter, P. (1979). Process flowsheeting. Cambridge University Press.