An Acceleration Method for Dogleg Methods in Simple Singular Regions

strategy over Newton's method, which can behave periodically or aperiodically, and Powell's original dogleg strategy, which can terminate at a singula...
0 downloads 0 Views 101KB Size
Download PDF
1358

Ind. Eng. Chem. Res. 1998, 37, 1358-1363

An Acceleration Method for Dogleg Methods in Simple Singular Regions Angelo Lucia* and Delong Liu Department of Chemical Engineering, University of Rhode Island, Kingston, Rhode Island 02881-0805

The behavior of dogleg methods in singular regions that have a one-dimensional null space is studied. A two-tier approach of identifying singular regions and accelerating convergence to a singular point is proposed. It is shown that singular regions are easily identified using a ratio of the two-norm of the Newton step to the two-norm of the Cauchy step since Newton steps tend to infinity and Cauchy steps tend to zero as a singular point is approached. Convergence acceleration is accomplished by bracketing the singular point using a projection of the gradient of the two-norm of the process model functions onto the normalized Newton direction in conjunction with bisection, thus preserving the global convergence properties of the dogleg method. Numerical examples for a continuous-stirred tank reactor and vapor-liquid equilibrium flash are used to illustrate the reliability and effectiveness of the proposed approach. Several geometric illustrations are presented. 1. Introduction In a sequence of recent papers, Lucia and co-workers (Lucia and Xu, 1992; Lucia et al., 1993; Sridhar and Lucia, 1995) have studied the behavior of trust region methods in the real domain and presented an extension of the dogleg method to the complex domain. In particular, Lucia and Xu (1992) studied the behavior of the dogleg strategy and extended it to the complex domain. Lucia et al. (1993) showed the advantages in numerical performance of the complex domain dogleg strategy over Newton’s method, which can behave periodically or aperiodically, and Powell’s original dogleg strategy, which can terminate at a singular point (or local minimum or saddle point in the two-norm of the function in the real domain). They also proposed a singular-point perturbation for moving iterates from a singular point to a solution. Sridhar and Lucia (1995) gave an analysis for the extended dogleg method showing that all nondegenerate singular points are saddle points of the two-norm of the model functions in the complex domain. This analysis also showed that an eigenvalue-eigenvector decomposition of the Hessian matrix of the two-norm of the function can be used to construct a path from a singular point to a solution. A summary of this equation-solving methodology can be found in the recent paper by Gow et al. (p 2843, 1997). In this paper, we concentrate on methods for the rapid and reliable convergence to singular points. There is considerable literature in the area of applied mathematics on the computation of singular points and related convergence acceleration methods. See, for example, Keller (1977), Rheinboldt (1978), Abbott (1978), Moore and Spence (1980), Decker and Kelley (1980), Griewank (1980), Griewank and Osborne (1981), Doedel (1981), Georg (1981), Decker and Kelley (1981), Kearfott (1983), Griewank and Reddien (1984), and others. However, most of these papers are concerned with the calculation of singular solutions (or turning and bifurcation points) and not necessarily singular points that * To whom correspondence should be addressed. Tel.: (401) 874-2814. Fax: (401) 874-4689. E-mail: [email protected].

do not represent solutions to the given set of nonlinear algebraic equations. Moreover, most of these algorithms are based on Newton’s method and do not consider trust region methods, and none address the issue of complexvalued points of singularity. In fact, there has been no systematic study addressing the performance of dogleg strategies in singular regions. Therefore, the primary objective of this research is to address the difficulties associated with slow convergence and termination in singular regions by developing a convergence acceleration method for the dogleg method. Accordingly, this paper is organized in the following way. In section 2 we present new methodologies for improving the numerical performance of the dogleg algorithm. These improvements consist of new computer tools for (1) the identification of singular regions and (2) convergence acceleration in the Newton direction while locating simple singularities. The ideas developed in section 2 are tested on a two-dimensional continuousstirred tank reactor (CSTR) problem and a multidimensional vapor-liquid equilibrium (VLE) problem in section 3. Conclusions based on the numerical results are presented in section 4. 2. Numerical Behavior of the Dogleg Method in Singular Regions The dogleg algorithm of Powell (1970) chooses between Newton, steepest descent (or Cauchy), and dogleg steps. At any iteration, the desired step is selected using a linear combination of the steepest descent step, µ, in the least-squares function and the Newton step, γ, constrained to lie within a trust region, to produce a monotonically decreasing sequence of least-squares function values that converge to a solution. Because the size of the trust region radius is also adjusted iteratively and because the rules of adjustment require that the trust region radius decrease when norm reduction does not occur, steepest descent steps are guaranteed to be selected in regions of difficulty and herein lies the global convergence characteristics of the dogleg strategy, which guarantees convergence to either a solution or a singular point of f(Z,p). We refer the

S0888-5885(97)00681-7 CCC: $15.00 © 1998 American Chemical Society Published on Web 02/25/1998

Ind. Eng. Chem. Res., Vol. 37, No. 4, 1998 1359

reader to the original paper of Powell (1970) and to those of Lucia and co-workers for the mathematical details of the dogleg strategy and its extension to the complex domain. In the material that follows, we present a strategy for improving the numerical performance of the dogleg method in simple singular regions. By a simple singular region, we mean a region which contains a singular point at which the null space of the Jacobian matrix of the model functions has dimension one. The proposed strategy is comprised of two parts: (1) the identification of a singular region by a ratio test; (2) convergence acceleration of iterates in the Newton direction. 2.1. Identification of Singular Regions. The normal convergence condition that ||f|| be less than a prescribed tolerance will not usually be satisfied in a region of singularity, especially if ||f|| is bounded well away from zero in this region. Therefore, the identification of singular regions is the first step in improving the performance of the dogleg method. When a local minimum of ||f|| (or equivalently a singular point) is being approached, any trust region algorithm will usually take Cauchy (or steepest descent) steps and converge linearly to this stationary point in the twonorm of a function. Furthermore, this convergence will usually be slow. The original condition for terminating the calculations proposed by Powell (1970) is

F(Z,p)k > M||g(Z,p)k||

(1)

where M is an arbitrarily chosen positive number governing the conditions for finishing the iterative process, usually set to an overestimate of the distance from Zk to the solution. Here, Z represents a vector of unknown variables, p is a vector of parameters, F(Z,p) ) ||f(Z,p)||2 ) fTf, where f(Z,p) is a vector of process model functions, ||g(Z,p)|| is the two-norm of the gradient of ||f||2 at (Z,p), where g(Z,p) ) JTf(Z,p), and J is the Jacobian matrix of f. The disadvantages of this original convergence condition are that it does not clearly characterize a region of singularity so that the number M is difficult to select and an excessive number of Cauchy steps (or trust region reductions) is usually required to meet condition 1. This also means that it requires many function evaluations. Moreover, sometimes the calculations terminate prior to reaching a stationary point or fail because the size of the trust region becomes extremely small (i.e., ||g(Z,p)k|| and

(2)

||γk|| gM ||µk||

(3)

where again M is some positive number and γk and µk are the Newton step and Cauchy step, respectively, on the kth iteration. We call condition 3 a ratio test since it is based on the ratio of the two-norm of the Newton and Cauchy steps. When conditions 2 and 3 are satisfied, the singular region can contain either a nondegenerate or degenerate singular point. The reasons for preferring conditions 2 and 3 are due to the facts that in a region of singularity the ratio of the twonorm of the Newton step to the two-norm of the Cauchy step tends to be a large number, and at a singular point, this ratio is infinite. In contrast, in nonsingular regions, the ratio is bounded. Compared with other methods, such as using a condition number for detecting singularity, conditions 2 and 3 capture the main characteristics of any singularity and are easily implemented within the dogleg algorithm. 2.2. Convergence Acceleration in Singular Regions. Once a region of singularity has been identified by conditions 2 and 3, a method of acceleration to move iterates quickly toward the singular point is needed. The goal of the acceleration step is to reduce the two-norm of a gradient of ||f(Z,p)||2 to an acceptable tolerance in a given direction by an efficient and reliable means without sacrificing the global convergence characteristics of the dogleg method. To achieve this goal in a singular region, higher order derivatives and/or geometric information associated with the function f or ||f||2 must be utilized. However, for a set of nonlinear equations in many variables, it can be difficult to obtain the second-order information of f and ||f||2 analytically. Here, we present a geometric approach for locating a simple singular point more efficiently on the basis of the reduction of ||g|| in the normalized Newton direction in a region of singularity. Figure 1 shows the underlying geometry for a two-dimensional problem. The basic idea of this method is to reduce the magnitude of the directional derivative of ||f||2 given by

DD(Z,p) ) gTγ ) [g1 g2]T[γ1 γ2]

(4)

in the normalized Newton direction to a desired tolerance, where gi is the ith component of the gradient of the function ||f||2 with respect to Zi and where γi is the ith component of the normalized Newton step. This is because at the projected minimum (i.e., singular point), the Newton step aligns with the null space of the Jacobian matrix of f(Z,p) and the gradient of ||f||2 is zero. This is a fact that we have known for some time (see Venkataraman and Lucia (1986)). One way of locating a projected local minimum is by maintaining a change in the sign of the directional derivative of ||f||2 in a bracketed region. To do this, the singular point is bracketed by placing an iterate on the trust region boundary in the Newton direction until a change in sign of the directional derivative occurs and the bracketed region is reduced in size by discarding the current end point whose directional derivative is the same as that at the bisection midpoint. This procedure is repeated until ||DD|| is reduced to a small number, say 10-7. The advantages of this approach are that it is reliable, efficient, and geometrically straightforward in accelerating convergence to a singular point. In particular,

1360 Ind. Eng. Chem. Res., Vol. 37, No. 4, 1998

Figure 2. Complex bifurcation diagram for a CSTR. Figure 1. Projection of the gradient of ||f(Z,p)||2 onto the Newton direction. Table 1. Data for a CSTR Example F (mol/mL) Cp (cal/(mol‚K)) C0 (mol/mL) T0 (K) θ (s) A (1/s) E (cal/mol) ∆H (cal/mol)

1.0 1.0 0.003 298.0 300 4.48 × 106 1.50 × 104 -5.0 × 104

bisection requires only one function evaluation per iteration once the singular region has been bracketed. Note, however, we have limited our attention to problems with a one-dimensional null space. For problems with higher-order singularities (e.g., critical points of mixtures, azeotropes, cusps, etc.) a similar approach can be used but additional research is needed to identify useful ways of accelerating the calculations since the null space has a dimension greater than one. 3. Numerical Examples In this section, we demonstrate the usefulness of the techniques developed in section 2 using an adiabatic continuous-stirred tank reactor (CSTR) with first-order, irreversible kinetics and the vapor-liquid equilibrium flash of a binary mixture in the retrograde region. 3.1. A CSTR Example. One mathematical model for an adiabatic CSTR derived from the conservation of mass and energy is given by

E RT X) E 1 + θA exp RT

(

θA exp -

(

)

)

(5)

and

X)

FCp(T - T0) C0(-∆H)

(6)

where T denotes the reaction temperature in degrees Kelvin, θ denotes the residence time in seconds, X is the conversion of the chemical reactant, F is the density of the fluid, Cp is the specific heat, C0 and T0 are the inlet concentration and temperature, respectively, ∆H denotes the heat of reaction, E is activation energy, A denotes the pre-exponential factor, and R is the universal gas constant. In the illustrations in this section, all other quantities were held fixed and their particular numerical values and units are given in Table 1. Figure 2 shows the complex bifurcation diagram for the physi-

Table 2. Numerical Results for a CSTR Example without Acceleration ni (nr)a

X

T (K)

||f||

||g||

1 (0) 2 (0) 3 (0) 4 (0) 5 (0) 6 (0) 7 (0) 8 (0) 9 (3) 10 (5) 11 (6) 12 (7) 13 (10) 14 (16) 15 (18) 16 (20) 17 (21)

0.700 00 1.045 98 1.012 40 0.979 48 0.946 38 0.913 17 0.879 89 0.846 58 0.813 26 0.804 92 0.807 00 0.805 96 0.806 48 0.806 42 0.806 42 0.806 42 0.806 42

450.000 454.988 449.988 444.988 439.988 434.988 429.988 424.989 419.989 418.739 419.051 418.895 418.973 418.963 418.964 418.964 418.964

47.1285 1.817 91 1.191 91 0.745 86 0.437 73 0.235 77 0.113 81 0.050 56 0.028 98 0.028 45 0.028 42 0.028 42 0.028 42 0.028 42 0.028 42 0.028 42 0.028 42

710.207 18.6311 1.477 62 0.561 18 0.089 46 0.110 24 0.155 71 0.133 31 0.084 81 0.001 84 0.001 78 0.000 70 0.000 40 0.000 05 7.49152 × 10-7 1.87266 × 10-7 9.36373 × 10-8

a n ) number of iterations, n ) number of trust region i r reductions. M ) 10 for condition 1 and 50 000 for condition 3.

cally meaningful solutions to eqs 5 and 6 parametrized in residence time from 0 to 800 s. Note that there are regions to the left and right of the two turning points in which there are real- and complex-valued solutions and a region between the two turning points in which there are three real-valued solutions. However, only the real-valued high-temperature (and corresponding highconversion) solution represents the desired operation of the reactor. The low-temperature solution corresponds to essentially no reaction and the intermediate solution is an unstable steady state. Moreover, the left-hand turning point, which occurs at 62.19 s, represents a lower limit of feasibility for the high-temperature/highconversion solution and consequently places limitations on the feed flow and/or volume that result in proper operation of the reactor. Now we study the behavior of the multivariable acceleration method on this two-variable CSTR problem, for which eqs 5 and 6 are solved simultaneously for the conversion, X, and temperature, T. However, first we geometrically demonstrate the numerical behavior of the trust region method in a singular region without acceleration. For initial values (X, T) ) (0.7, 450 K) with a specified residence time of 60 s and an initial trust region radius of 5, the values of X and T at each iteration are shown in Table 2. The associated level curves of ||f|| are shown in Figure 3. At a residence time of 60 s, only one real-valued solution at (X, T) ) (0.00277, 298.415 K) exists since we are to the left of the turning point at 62.19 s. See Figure 2. Note also that iterates 1-6 lie outside the region shown in Figure 3; only iterates 7-12 are shown. Observe that all iterates are located in a narrow strip,

Ind. Eng. Chem. Res., Vol. 37, No. 4, 1998 1361

Figure 3. Level curve of ||f|| for CSTR model equations.

and iterates 9-11 are alternately distributed on both sides of the 12th iterate as ||f|| decreases on each step. On iteration 12, ||f|| ) 2.84218 × 10-2, ||g|| ) 7.04209 × 10-4, and ||γ||/||µ|| ) 5.46 × 104. However, neither the two-norm of the gradient of ||f||2 nor the two-norm of f is close to vanishing and the calculations require increasingly excessive trust region reductions (and therefore function evaluation) in the Cauchy direction on this and all subsequent iterations. In fact, from iteration 9-17, 106 additional function calls are required by trust region reduction to reach a value of ||g|| < 10-7, as shown in Table 2. We can still, however, recognize that there is a singular point (a local minimum in the real domain but saddle point in the complex domain) near (X, T) ) (0.80596, 418.895 K). Note also that with M ) 10 in condition 1, a value suggested by Powell, the calculations terminate on iteration 14, very close to the singular point. However, 47 additional function calls are needed to reach this point. Now we apply acceleration to this example to see if the singular point shown in Table 2 (i.e., iterate 17) can be located more efficiently. On the 10th iteration, the two-norm of the gradient is 1.8419 × 10-3, the two-norm of f(Z,p) is 2.8445 × 10-2, and the ratio of the two-norm of the Newton step to the two-norm of the Cauchy step is 6.7722 × 104. Accordingly, conditions 2 and 3 indicate a region of singularity is nearby and acceleration is automatically invoked on iteration 10sfour iterations sooner than Powell’s condition. Following the identification of this singular region, the singular point is bracketed and reduction of the directional derivative of ||f||2 in the normalized Newton direction (i.e., (∆X, ∆T) ) (0.00667, 0.99998 K)) is accomplished by bisection. The bisection results are shown in Tables 3 and 4. The value of the directional derivative on the ninth bisection iteration clearly indicates that a singular point is found at (X, T) ) (0.80642, 418.963 K). Note that the calculated singular point is, indeed, reasonably close to the 10th-12th iterates shown in Table 2, and that the line connecting the 10th through the 12th iterates can be approximately considered the null space of the Jacobian matrix of f(X,T). As a consequence, we can

Table 3. Numerical Results for a CSTR Example with Acceleration bisection step

X

T (K)

directional derivatives

1 2 3 4 5 6 7 8 9

0.821 59 0.813 25 0.809 09 0.807 00 0.805 96 0.806 48 0.806 22 0.806 35 0.806 42

421.239 419.989 419.364 419.051 418.895 418.973 418.934 418.954 418.963

8.04000 × 10-5 3.18320 × 10-5 1.18966 × 10-5 2.56577 × 10-6 -1.98890 × 10-6 2.80545 × 10-7 8.56060 × 10-7 -2.88237 × 10-7 -3.96758 × 10-9

use any of these normalized Newton steps to approximate the null space of the Jacobian matrix of f(X,T). In addition, we note that each bisection iteration requires only one function evaluation, and thus we have reduced the number of function evaluations required to find this singular point from 44 required by Powell’s condition to 9 needed by our ratio test. Finally, for a residence time of 750 s, with initial values of (X, T) ) (0.20, 350 K), a real-valued local minimum (or singular point) is found at (X, T) ) (0.09372, 312.059 K) with five acceleration steps reaching a final directional derivative value of 5.13290 × 10-9. Note again that this singular point is a saddle point in the complex domain. 3.2. A Vapor-Liquid Equilibrium (VLE) Flash Example. Figure 4 shows the experimental (Olds et al., 1949) and calculated VLE for carbon dioxide and n-butane at 344.26 K using the Soave-Redlich-Kwong (SRK) equation of state. The calculated results in this figure were generated by solving the component mass balances,

fi - li - vi ) 0,

i ) 1, 2, ..., nc

(7)

the phase equilibrium equations,

Ki(li/

∑lj) - (vi/∑vj) ) 0,

the vapor specification equations,

i, j ) 1, 2, ..., nc (8)

1362 Ind. Eng. Chem. Res., Vol. 37, No. 4, 1998 Table 4. Numerical Results for Bisection for a CSTR Example bisection step

mid point

end points

1 2 3 4 5 6 7 8 9

(0.821 59, 421.239) (0.813 25, 419.989) (0.809 09, 419.364) (0.807 00, 419.051) (0.805 96, 418.895) (0.806 48, 418.973) (0.806 22, 418.934) (0.806 35, 418.954) (0.806 42, 418.963)

[(0.804 92, 418.739), (0.838 25, 423.748)] [(0.804 92, 418.739), (0.821 59, 421.239)] [(0.804 92, 418.739), (0.813 25, 419.989)] [(0.804 92, 418.739), (0.809 09, 419.364)] [(0.804 92, 418.739), (0.807 00, 419.051)] [(0.805 96, 418.895), (0.807 00, 419.051)] [(0.805 96, 418.895), (0.806 48, 418.973)] [(0.806 22, 418.934), (0.806 48, 418.973)] [(0.806 35, 418.953), (0.806 48, 418.973)]

Figure 4. Isothermal VLE for CO2(1)-n-C4H10(2) at 344.26 K using the SRK EOS.

(

∑vi - V)/F ) 0

(9)

V/F - O ) 0

(10)

and

to point out that this singular point is a saddle point in the complex domain, that 74 of the 100 iterations without acceleration are Cauchy steps, that many of these Cauchy steps require trust region reductions resulting in more than 322 function evaluations to reach a value of ||g|| < 10-5 and that condition 1 still cannot identify the singular point. In our opinion, this is computationally wasteful. If, on the other hand, the parameter M in eq 3 is set to 107 and the acceleration procedure in section 2 is used, the ratio test invokes acceleration on the 27th iteration, for which the normalized Newton direction is (∆l1, ∆l2, ∆v1, ∆v2, ∆V/F, ∆P) ) (1.2187 × 10-2, -1.2187 × 10-2, 6.769 × 10-7, -6.645 × 10-7, -2.2601 × 10-7, 0.99985) and represents a very good approximation to the null space of the Jacobian matrix. After this, 12 bisection iterations are needed to find the approximate singular point at which the value of ||g|| is 1.9244 × 10-4. Thus, a total of 39 function evaluations (instead of more than 322) from start to finish using our acceleration procedure are required and this, in our opinion, represents a considerable computational savings. 4. Conclusions

at fixed temperature and feed conditions, for li, vi, V/F, and P, where fi, li, and vi are the ith component molar flow rates of the feed, liquid and vapor respectively, Ki is the ith component equilibrium ratio, F is the total feed flow, V/F is the vapor-to-feed ratio, which is specified at a value of O. This mixture exhibits retrograde behavior and has a turning point at yCO2 ) 0.7855 and P ) 69.40 bar. Moreover, because of continuity of first partial derivatives, there are singular points which are not solutions for feed specifications either just to the left or just to the right of the turning point. Consider then 1 mol of feed with a composition of z ) (z1, z2) ) (0.79, 0.21), which is just to the right of the turning point and for which the two-phase solution is complex-valued, indicating that the only real-valued solution is a single-phase solution. See Figure 4. For these feed specifications, the extended dogleg method without acceleration finds a point at Z ) (l1, l2, v1, v2, V/F, P) ) (0.5921, 0.4075, 0.790, 0.2102, -1.318 × 10-4, 71.67) in 100 iterations, at which F ) fTf is 1.1 × 10-6 and ||g|| is 8.5 × 10-6, where the subscripts 1 and 2 denote carbon dioxide and n-butane, respectively. However, Powell’s condition 1 with M ) 10 is still not satisfied even after 100 iterations. The eigenvalues of the Hessian matrix of the two-norm of the model functions at the singular point are (0.02122, 0.0119, 0.0008, 0.00009, 0.00002, 0.0) which clearly show that the two-norm is convex in this region and that the singular point is truly a local minimum in the two-norm when the calculations are restricted to or remain in the real domain. This means that Cauchy steps will converge to the singular point. However, it is important

A strategy comprised of the identification of singular regions by a ratio test and an acceleration method using the reduction of the gradient of the square of the twonorm of the model functions in the normalized Newton direction in a singular region was proposed and tested on two process engineering examples, a CSTR example, and a VLE problem. Numerical results showed that the proposed strategy avoids excessive Cauchy steps or trust region radius reduction failures in singular regions. It was also shown that the normalized Newton direction in a singular region can be considered a reasonable approximation of the null space of the Jacobian matrix and that reducing the projection of the gradient of ||f(Z,p)||2 in that direction is a reliable means for determining a singular point. The concepts that underlie this strategy are simple and its implementation within a dogleg method is straightforward. Also, this strategy is suitable for use in other types of trust region algorithms. Furthermore, the enhanced dogleg algorithm can be used for locating simple local minima of ||f(Z,p)||2. With some modifications, it can also be used to determine turning points in a given domain. Acknowledgment This work was supported by the National Science Foundation under Grant Nos. CTS-9409158 and CTS9696121. Nomenclature A ) frequency factor (s-1) C0 ) inlet concentration (mol mL-1)

Ind. Eng. Chem. Res., Vol. 37, No. 4, 1998 1363 Cp ) specific heat (cal/mol‚K) DD ) directional derivative of gradient of ||f(Z,p)||2 in a normalized Newton direction E ) activation energy (cal mol-1) fi ) component i feed molar flow rate to flash (mol s-1) f(Z,p) ) model functions F ) total feed molar flow rate (mol s-1) F(Z,p) ) ||f||2, square of two-norm of model functions F(Z,p)k ) two-norm of model functions at kth iterate gi ) ith component of the gradient of F(Z,p) g(Z,p)k ) gradient of F(Z,p) at kth iterate evaluated at (Z,p) ∆H ) heat of reaction (cal mol-1) J ) Jacobian matrix J Ki ) component i equilibrium ratio li ) component i liquid molar flow rate from flash (mol s-1) M ) a positive number nc ) number of components R ) universal gas constant (cal mol-1‚K-1) p ) a real-valued parameter vector P ) pressure (bar) X ) conversion of reactant T ) temperature (K) T0 ) inlet temperature (K) vi ) component i vapor molar flow rate from flash (mol s-1) V ) total overhead vapor molar flow rate from flash (mol s-1) z ) feed composition Z ) vector of unknown variables Greek Letters γ ) normalized Newton step φ ) specified vapor-to-feed ratio F ) density of reactant (mol ml-1) θ ) residence time (sec) µ ) Cauchy step Subscripts i ) ith component Superscripts k ) iteration counter n ) dimension of vector of unknowns

Literature Cited Abbott, J. P. An Efficient Algorithm for the Determination of Certain Bifurcation Points. J. Comput. Appl. Math. 1978, 4, 19. Decker, D. W.; Kelley, C. T. Newton’s Method at Singular Points. SIAM J. Numer. Anal. 1980, 17, 66.

Decker, D. W.; Kelley, C. T. Convergence Acceleration for Newton’s Method at Singular Points. SIAM J. Numer. Anal. 1981, 19, 219. Doedel, E. J. Auto: A Program for the Automatic Bifurcation Analysis of Autonomous Systems. Congr. Numer. 1981, 30, 265. Georg, K. On Tracing an Implicitly Defined Curve by QuasiNewton Steps and Calculating Bifurcation by Local Perturbations. SIAM J. Sci. Stat. Comput. 1981, 2, 35. Gow, A. S.; Guo, X.; Liu, D.; Lucia, A. Simulation of Refrigerant Phase Equilibria. Ind. Eng. Chem. Res. 1997, 36, 2841. Griewank, A. O. Starlike Domains of Convergence for Newton’s Method at Singular Points. Numer. Math. 1980, 35, 95. Griewank, A. O.; Osborne, M. R. Newton’s Method for Singular Problems When the Dimension of the Null Space is > 1. SIAM J. Numer. Anal. 1981, 18, 145. Griewank, A. O.; Reddien, G. W. Characterization and Computation of Generalized Turning Points. SIAM J. Numer. Anal. 1984, 21, 176. Kearfott, R. B. Some General Derivative-Free Bifurcation Techniques. SIAM J. Sci. Stat. Comput. 1983, 4, 52. Keller, H. B. Numerical Solution of Bifurcation and Nonlinear Eigenvalue Problems. In Applications of Bifurcation Theory; Rabinowitz, P. H., Ed.; Academic Press: New York, 1977. Lucia, A.; Xu, J. Global Minima in Root Finding. In Recent Advances in Global Optimization; Floudas, C. A., Pardalos, P. M., Eds.; Princeton University Press: Princeton, NJ, 1992. Lucia, A.; Guo, X.; Richey, P. J.; Derebail, R. Simple Process Equations, Fixed Point Methods, and Chaos. AIChE J. 1990, 36, 641. Lucia, A.; Guo, X.; Wang, X. Process Simulation in the Complex Domain. AIChE J. 1993, 39, 461. Moore, G.; Spence, A. The Calculation of Turning Points of Nonlinear Equations. SIAM J. Numer. Anal. 1980, 17, 567. Olds, R. H.; Reamer, H. H.; Sage, B. H.; Lacey, W. N. Phase Equilibrium in Hydrocarbon System. The n-Butane-Carbon Dioxide System. Ind. Eng. Chem. Res. 1949, 41, 475. Powell, M. J. D. A Hybrid Method for Nonlinear Equation. In Numerical Methods for Nonlinear Algebraic Equations; Rabinowitz, P., Ed.; Gordon and Breach: London: 1970. Rheinboldt, W. C. Numerical Methods for a Class of Finite Dimensional Bifurcation Problems. SIAM J. Numer. Anal. 1978, 15, 1. Sridhar, L. N.; Lucia, A. Process Analysis in the Complex Domain. AIChE J. 1995, 41, 585. Venkataraman, S.; Lucia, A. Exploiting the Gibbs-Duhem Equations in Separation Calculations. AIChE J. 1986, 32, 1057.

Received for review September 24, 1997 Revised manuscript received January 14, 1998 Accepted January 15, 1998 IE970681F