Reverse Iteration in Chemical Process Simulation - ACS Publications

4332

Ind. Eng. Chem. Res. 1998, 37, 4332-4340

Reverse Iteration in Chemical Process Simulation Angelo Lucia* and Delong Liu Department of Chemical Engineering, University of Rhode Island, Kingston, Rhode Island 02881-0805

A new technique for generating approximations of Julia sets, called reverse iteration, is proposed. The concept of reverse iteration is simple and based on the idea of solving an iterative map of the general form Zk ) G(Zk-1, p) for Zk-1 instead of Zk. It is shown that if the process of reverse iteration is initiated at any singular point, the collection of inverse images must be members of the Julia set since singular points are in the Julia set and the Julia set is closed. This sequence of reverse iterates, say {Zk-1}, is necessarily distributed throughout the basin boundaries. It is also shown that reverse iteration can have multiple inverse images and a tree structure for the Julia set but that the associated potential combinatorial computational demand is easily resolved by exploiting the fractal nature of any Julia set. From this, practical ways generating initial values that converge to solutions to the given model equations are proposed. Several examples and geometric illustrations are used to elucidate key concepts. 1. Introduction With the recent industrial shift toward equationoriented process simulation, a number of important equation-solving concerns have surfaced. For the most part, Newton or Newton-based methods serve as the engine for many of the current equation-oriented process simulators (e.g., SPEEDUP, HYSIS, PROCESS, and others) and, as a consequence, singularity has become an important practical issue. Surprisingly, even though Newton-like methods have been used for many years in sequential modular architectures (FLOWTRAN, ASPEN Plus, HYSIM, etc.), there is very little information on the impact of singularity in process simulation literature. If a singular region or singular point is encountered during iteration, numerically wild iterates can result and the simulator will usually terminate for reasons of overflow. Furthermore, even if the steps are damped by line searching, trust regions, or step bounding or in some other manner, other related numerical difficulties (e.g., slow convergence, excessive line searching, trust regions that tend to zero, etc.) usually occur. Unlike the situation in the more traditional sequential modular approach, singularity in an equation-oriented architecture is even more difficult to “trap”. This is because the model equations are treated as a single set and, as a result, it is more difficult to isolate the cause of the problem and/or use physical insight to intervene when a singular region or singular point is encountered. When norm-reduction is used to stabilize the equation-solving aspects of any process simulation, the issue of singularity becomes a problem in global optimization, in which solutions to the model equations represent zero-valued global minima and singular points that are attractors are either nonzero-valued local minima or saddle points of the two-norm. Consequently, mixed integer nonlinear programming (MINLP) algorithms as well as other methods of global optimization represent a rigorous way of resolving singularity. However, because many model equations in chemical process simulation are nonconvex, it is unclear whether any rigorous guarantees can be established using MINLP * To whom correspondence should be addressed. Tel.: (401) 874-2814. Fax: (401) 874-4689. E-mail: [email protected].

methods. Other deterministic methods such as the tunneling algorithms of Levy and Montalvo1 also address singularity. However, in our opinion, all tunneling methods are somewhat undesirable because they require tunable parameters such as the pole strength for deflating the function at a point of singularity. Thus they may work in some cases and not in others. Finally, (complex domain) trust region methods of Lucia and coworkers (see, Sridhar and Lucia2) allow convergence to a singular point and then use an eigenvalue-eigenvector decomposition of the Hessian matrix of the two-norm of the model equations to move from a singular point to a solution. However, eigenvalue-eigenvector decomposition can be computationally expensive for large problems if all eigenvalues and eigenvectors are needed and many workers are reluctant to use complex domain process models and numerical methods. In this paper, we take an entirely different approach to resolving the computational difficulties associated with singularity in process simulation. In particular, reverse iteration is used to generate approximations of the boundaries of the basins of attraction, thereby making it possible to choose initial values that are guaranteed to converge to a solution of the model equations. Thus, in section 2, the connection between singularity, Julia sets, and basin boundaries is drawn. The concept of reverse iteration is introduced in section 3 and the associated mathematical formalism is presented in section 4. Numerical examples and geometric illustrations for a continuous stirred tank reactor, a cubic equation of state, and a flash problem are presented in section 5. Concluding remarks and open issues are discussed in section 6. Finally, it is important to stress at the outset that the techniques presented in this paper are equally applicable to real and complex domain process models. 2. Julia Sets and Basin Boundaries Although the theory of Julia sets is quite old, dating back to the early 1900s (see Julia3 and Fatou4), it is only through the publication of Mandelbrot5 that they have resurfaced as a topic of interest (see Barnsley6 and Nusse and Yorke7). In simple terms, the Julia set of a given iterative map can be thought of as the collection

10.1021/ie9802764 CCC: $15.00 © 1998 American Chemical Society Published on Web 10/10/1998

Ind. Eng. Chem. Res., Vol. 37, No. 11, 1998 4333

Figure 1. Basins of attraction for Newton’s method for a CSTR example. Table 1. Numerical Values for Quantities in Equation 4 quantity

numerical value

F (Gmol/cm3) Cp (cal/Gmol K) A (cm3/s) E (cal/Gmol) ∆H (cal/Gmol) C0 (Gmol/cm3) T0 (K) R (cal/Gmol K)

1.0 1.0 4.48 × 109 1.5 × 104 -50 000 0.003 298. 1.987

of initial values (or basin of attraction) that give periodic or aperiodic behavior. Although not dense (or fractal) in a topological sense, the Julia set does lie in the boundary between different basins of attraction of the solutions to a given map. These are key observations that have practical use. Figure 1 is the Julia set for Newton’s method for solving for the temperature roots of the combined mass/energy balance equation for a continuous stirred tank reactor (CSTR) studied in Lucia and Liu,8 which can be found in the section on Numerical Illustrations (i.e., eq 4) along with the associated numerical data (Table 1). Note that Figure 1 is plotted in the complex domain and that the Julia set (i.e., the set of black points) is dispersed throughout the boundaries between the basins of attraction for the real low temperature and complex-valued high temperature roots. It is important to understand that there are members (or equivalently points) of the Julia set in each of the boundaries between regions of different color, even though it is not entirely evident from Figure 1. Figure 2, which is an enlargement of the boundary in the neighborhood of the singular point on the left, clearly shows this. Note again that there are black points dispersed throughout the basin boundaries. Enlargement also shows there are points in the Julia set between the blue and green regions in the imaginary direction of Figure 1, which is why the distinction between complex conjugate roots is made. The basins of attraction for Figure 1, as well as many other related figures in this manuscript, were generated for a combined mass/energy balance equation expressed

in terms of a single complex-valued temperature. What about its real domain counterpart? Figure 1 also clearly shows those results as well. In moving from the left end of Figure 1 to the right along the real line, note there is a single physically meaningful real-valued solution (i.e., a low temperature solution at 298.429 K) on the far left, a basin of attraction for the low temperature solution (i.e., the red region along the real line on the left), spurious real-valued solutions (i.e., the yellow region along the real line), two real-valued singular points, and the Julia set, which is everything on the real line to the right of the yellow region minus those (red) points in that region that jump to the low temperature solution. There is no blue or green region along the real line since those points are associated with the complex-valued solutions. Note also that the Julia set is fractal (i.e., there are points between the black dots that jump to the low temperature solution) and includes points to the right of the singular point on the right. Moreover, because the Julia set contains many points in this region, it is clear that arbitrary perturbation of the singular point would not necessarily result in convergence to the low temperature solution and that some degree of rigor is truly needed to guarantee convergence. Figure 3 shows the basins of attraction for the extended trust region (or dogleg) strategy for the same reactor temperature calculation. See Lucia and Xu9 and Lucia et al.10 for a formal description of complex domain trust region methods. Note the similarity in overall geometric structure and, in particular, that the basin of attraction for the singular points looks very much like the Julia set for Newton’s method shown in Figure 1, except that it is no longer fractal since norm reduction is forced. Although all singular points are saddle points in the complex domain, the singular point at 419.170 K is a nonzero minimum in the two-norm of the function in the real domain and attracts many initial values while the one at 363.031 K is a maximum and attracts only very few initial values (see, p S572 in Lucia and

4334 Ind. Eng. Chem. Res., Vol. 37, No. 11, 1998

Figure 2. Enlarged basins of attraction for Newton’s method for a CSTR example.

Figure 3. Basins of attraction for the extended dogleg method for a CSTR example.

Guo11). All numerical details associated with the generation of Figures 1-3 can be found in Lucia and Liu.8 Three properties of Julia sets are useful for our purposes. First, we have observed, through a number of mathematical and chemical process examples, that all singular points of an iterative map are members of the Julia set. Thus any time an iterative map encounters a point of singularity, we know this singular iterate is in the Julia set; see Figure 1. Second, it is well-known that the Julia set of any iterative map is a closed set (see, Peitgen et al.12). Consequently, once an iterate lands in the Julia set, subsequent iterates of a given map cannot escape. Third, Julia sets are not topologically dense; they are fractal. That is, for any x ∈ J and δ > 0, there are y ∈ B(x, δ), the open ball of radius d centered at x such that y ∈ J. Therefore the one-

dimensional linear subspace containing (or line connecting) any two members of a Julia set must contain points (initial values) from one or more (other) basins of attraction of solutions to the given model equations. 3. The Concept of Reverse Iteration The idea of reverse iteration is simple. Once an attracting singular point is encountered, we exploit the closure property stated briefly in section 2. That is, when a singular iterate occurs, we know that this singular point and all iterates from the initial guess to the singular point are in the Julia set for the given map. Moreover, any further iteration will, in most cases, carry the iterates to infinity. However, because the Julia set is a closed set, it is possible to run the iterative map in


reverse from the singular point and in that way build an approximation of the Julia set. This, in turn, provides insight into the structure of the basins of attraction and their boundaries. With this, a new initial value can be chosen in such a way as to guarantee the following: (i) The new initial value can be chosen from the basin of attraction of at least one solution since the Julia set approximation is not dense and the linear subspace containing any two members of the Julia set must be in some basin of attraction because of the fractal nature of the basin boundaries. Consequently, the calculations must converge to some solution. (ii) Subsequent iterates cannot return to any singular point because the Julia set is a closed set and this new initial value has been deliberately chosen to lie outside the Julia set. Figure 4. Reverse Newton iterates for a CSTR example.

4. Mathematical Formalities Here we state briefly relevant mathematical concepts and sketch one possible way of building an approximation of the Julia set for a nonlinear iterative map. It is important to keep in mind that what follows is equally applicable to either the real or complex domain process models. Consider the system of nonlinear algebraic model equations given by

F(Z,p) ) 0

(1)

where F takes some set D in either R2n or Cn into itself, Z is a vector of unknown variables of appropriate length, and p is some (real or complex) vector of parameters of length m. Furthermore, let Z′ denote a singular point (any iterate for which the Jacobian matrix of F, say J, is rank deficient) and assume that some numerical method for solving eq 1 (e.g., Newton’s method with line searching or a trust region strategy) has converged to this singular point. The Concept of Reverse Iteration Is Straightforward. For example, Newton’s method can be reversed by solving

GRN(Zk-1,p) ) Zk-1 - [J(Zk-1,p)]-1F(Zk-1,p) - Zk ) 0 (2) for Zk-1 instead of Zk, where GRN denotes the reverse Newton map. To do this, we let Zk ) Z′ (i.e., any singular point) and determine a solution to eq 2, say Zk-1. Note that the Jacobian matrix in eq 2 is evaluated at Zk-1 and not at the singular point, thus it is invertible. Then we replace Zk with Zk-1 in eq 2 and solve for a solution, Zk-2. Note we are always determining a solution to the reverse Newton iteration function given by eq 2. We continue to iterate backward, generating in turn solutions to eq 2 Zk-3, Zk-4, Zk-5, ..., until a sufficient sequence of inverse iterates {Zk-1} has been generated to give a rough picture of the Julia set for the given iterative map. Note also that the structure of the Julia set depends on the values of the parameters, and, while any iterative map can, in principle, be reversed (when written in the form Zk-1 ) G(Zk,p)), some maps are more complicated to reverse than others. There are several fundamental issues related to solving eq 2 that need to be mentioned. These include (i) actually solving eq 2, (ii) the existence of multiple inverse images, (iii) the number of inverse images

needed to approximate the Julia set, and (iv) the number of inverse images needed to generate “good” initial values. Solving the Reverse Iteration Problem. We use the complex domain dogleg strategy of Lucia and coworkers (Lucia et al.10) with second-order acceleration (Lucia and Liu8) to solve eq 2 since it is stable with respect to all solutions. The mechanics of using any Newton-like method to solve eq 2 are not entirely transparent, especially in the multivariable case. This is because the first partial derivatives of eq 2 involve the function F(Zk-1,p) and the derivatives of its inverse Jacobian matrix. That is, the Jacobian matrix of reverse Newton’s method, denoted J(GRN), is

J(GRN) ) I - [J(Zk-1,p)]-1J(Zk-1,p) + F(Zk-1,p)T(Zk-1,p) ) F(Zk-1,p)T(Zk-1,p) (3) where T(Zk-1,p) denotes a third-order tensor associated with the derivatives of the inverse Jacobian matrix in all but the single-variable case. Note that the product of F(Zk-1,p) and T(Zk-1,p) is a matrix. While the matrix J(GRN) can be obtained analytically, it is also possible to approximate it (or its inverse) using a quasi-Newton method like the Broyden (or inverse Broyden) update, where the secant information is obtained from differences in eq 2. Other known constraints on J(GRN), in addition to secant conditions, could also be approximated by using the techniques in Venkataraman and Lucia.13 Multiple Inverse Images. Figures 4 and 5 show the functionality of F(Z, p) and the reverse Newton iteration function respectively, for the CSTR example studied in Lucia and Liu.8 In this example, the unknown complex-valued variable Z is the temperature of the reactor, T, and the residence time, θ, is a realvalued parameter. The important point to note in Figure 4 is that the iterate Tk (i.e., the singular point at 419.170 K) has at least three inverse images Tk-1 (i.e., there are at least three solutions to eq 2 in this case). This is an important point because it gives rise to a tree structure for the sequence of inverse iterates, say {Tk-1}. Furthermore, since there are three inverse images for the singular point (i.e., the base node of the tree), there can be nine inverse images at the next level of the tree, 27 at the next, 81 at the next, and so on. This is further illustrated with a numerical example in section 5.


Figure 6. Bifurcation diagram parametrized in residence time.

5. Numerical Illustrations

Figure 5. Reverse Newton function for a CSTR example.

Figure 5, on the other hand, shows the functionality of reverse Newton iteration for the exact same situation. Note that all of the solutions lie to the right of the lefthand asymptote and the two solutions on the right are very close to the right-hand asymptote, making them somewhat sensitive to initial values. Tree Structure of the Julia Set. Because the existence of multiple inverse images can give rise to a tree structure of inverse iterates, construction of the Julia set by reverse iteration can be subject to combinatorial growth. In one aspect of this work, we were interested in generating the complete tree structure to some arbitrary level in order to test the concept of generating Julia sets by reverse iteration. However, for the purpose of generating initial values that converge to a solution to the model equations, this tree structure and the associated combinatorial growth can be avoided. Generating Initial Values. In a problem-solving environment, we are interested in ways of choosing initial values that converge to solutions in a reliable, efficient, and computationally inexpensive way. Thus, the focus for generating one or more initial values to find solutions to the model equations is very different than that for building a “good” approximation of the Julia set. In particular, the number of inverse images required to find initial values that will converge to solutions of the model equations is usually considerably smaller than the number of inverse images needed for an approximation of the Julia set. In fact, sometimes only one inverse image is needed. To find initial values that converge to one or more solutions to the model equations, we exploit the nondense topological (or fractal) structure of Julia sets. In particular, given any two arbitrary members of a Julia set, we know that the one-dimensional linear subspace containing (or line connecting) these points must also contain points in other basins of attraction for one or more solutions to the model equations because of the fractal nature of the basin boundaries. Consequently, initial values can be selected along this line by any number of simple strategies (i.e., bisection, golden section, Fibonacci search, Newton predictor-corrector strategies, etc.) so that convergence to one or more solutions will occur. We note that, unlike homotopycontinuation methods, the proposed strategy has the capability of finding path-disconnected solutions since it is based on the dispersion of points in the Julia set throughout the basin boundaries.

In this section, we provide numerical support for the basic ideas outlined in sections 3 and 4 using models for an adiabatic CSTR, a cubic EOS, and a VT flash problem. 5.1. A CSTR Model. One mathematical model for an adiabatic CSTR with first-order, irreversible, and exothermic kinetics (see, Smith14) is given by the single complex-valued mass/energy balance equation

F(T, θ) ) [FCp/(C0∆H)](T - T0) + θA exp(-E/RT)/ [1 + θA exp(E/RT)] ) 0 (4) where T denotes the reaction temperature in degrees K, θ denotes the residence time in seconds, F is the density, Cp is the specific heat, C0 and T0 are the inlet concentration and temperature, respectively, ∆H is the heat of reaction, E is the activation energy, A is the preexponential factor, and R is the gas constant. As noted earlier, the unknown variable in eq 4 is T, the parameter is the residence time and the particular numerical values for the other quantities that appear in eq 4 are given in Table 1. Figure 6 shows a bifurcation diagram for the solutions to eq 4 parametrized in residence time. Note that at low values of the residence time, in particular for values less than 62.19 s, there is no real-valued high temperature/high conversion solution and it is in this region that Newton’s method behaves periodically/chaotically and trust region methods converge to an attracting singular point (or nonzero-valued minimum in the twonorm of the combined mass/energy balance equation); see Figures 1 and 3. Thus at these residence time values there exists a Julia set for Newton’s method that contains a significant number of points in addition to the singular points. For residence times in which there are three real-valued solutions, the only members of the Julia set for Newton’s method are the singular points and points at infinity. To illustrate the ideas associated with building approximations to the Julia set by reverse iteration we chose a residence time of 62 s, just to the left of the lefthand bifurcation point shown in Figure 6. Again, Figure 1 shows the basins of attraction for Newton’s method for this value of the residence time. Moreover, for this particular value of θ the physically meaningful solutions are T1 ) 298.429 K, T2 ) 419.239 + 2.1406i, and T3 ) 419.239 - 2.1406i, while the singular points are given by T1′ ) 363.031 K and T2′ ) 419.170 K. The important points to note here are that both singular points are real-valued and members of the Julia set and that the singular point at 419.170 K attracts many more initial values than the one at 363.031 K.


Figure 9. Tree structure of the Julia set from reverse Newton iteration for a CSTR example. Figure 7. Julia set for a CSTR example.

Figure 8. Initial values from approximate Julia set for a CSTR example.

The reverse Newton iteration for this example is given by

GRN(Tk-1,θ) ) Tk-1 - F(Tk-1,θ)/F′(Tk-1,θ) - Tk ) 0 (5) where F′(Tk-1,θ) is the first derivative of the combined mass/energy balance with respect to temperature. Note 1/F′(Tk-1,θ) is the inverse Jacobian matrix in eq 2 in this case and because the model for the CSTR is a function of a single complex variable, the Jacobian matrix for the reverse Newton iteration function can be obtained analytically and is given by

J(GRN) ) F(T,θ)F′′(T,θ)/[F′(T,θ)]2

(6)

where F′′(T,θ) is the second derivative of the combined mass/energy balance with respect to temperature. Note the tensor T(T,θ) ) F′′(T,θ)/[F′(T,θ)]2 is a scalar or 1 × 1 matrix for this problem. Figure 7 shows the true Julia set given in Figure 1, which consists of initial values or aperiodic sequences of Newton iterates. Figure 8, on the other hand, is the approximation to the Julia set obtained by repeatedly solving eq 5 starting from the singular point at 419.170 K, using a total of 90 reverse Newton iterations. At each reverse iteration, eq 5 was solved by using a single variable complex domain dogleg strategy (Lucia and Xu;9 Lucia and Liu8) because Newton’s method behaves periodically at some nodes and direct substitution is unstable with respect to some of the fixed points of eq

5. In both Figures 7 and 8, the Julia set and its approximation are denoted by black dots and the singular points by gray dots. Note that the approximate Julia set gives a very reasonable gross approximation of the true Julia set. This is because the true Julia set represents an exhaustive search of the space of initial values, while the approximate Julia set is generated by a finite tree structure. To illustrate the correctness of generating approximate Julia sets by reverse iteration, note that reverse iteration actually finds points in the complex plane along the boundary between the pair of complex conjugate temperature roots (see also Figure 1) and that reverse iteration does not find inverse images to the left of the singular point at 363.031 K (that would be incorrect), both of which are characteristics of the true Julia set. The Tree Structure of the Julia Set. For the base node or the singular point at 419.170 K, there are three inverse images; their actual numerical values are 388.257, 416.875, and 421.194 K. These are the inverse iterates shown in Figure 4. Each of these nodes, in turn, also has multiple inverse images, some real, some complex-valued. For example, there are nine solutions to eq 5 at the second level of the tree. The first three inverse images, 386.519 + 29.849i, 386.519-29.849i, and 419.244 K, are associated with the node at 388.257 K, the second group of three solutions, 390.839, 412.978, and 419.189 K, come from the node at 416.875 K, and the last three inverse images, 386.580, 418.211, and 423.759 K, are associated with the node at 421.194 K. This combinatorial growth of the tree continues as it branches. Several levels of the tree are shown in Figure 9. Generating Initial Values. To avoid many of the computational disadvantages associated with the tree structure of the Julia set (i.e., dependence on starting point, combinatorial demand, etc.), we simply use a few inverse images and the fractal nature of the Julia set in any linear subspace to find initial values that will converge to a solution of the combined mass/energy balance equation. Several linear subspaces containing a few members of the Julia set are shown in Figure 8. The inverse images are numbered in the order in which they were generated starting with the singular point at 419.170 K. If we select the midpoint of the line connecting the singular point, 419.170 K, and the inverse image, 386.519 + 29.849i, as a starting point, it is easily seen from Figure 1 that this initial value is well within the green region or basin of attraction of


Figure 11. Initial values from Julia set for Newton’s method for SRK EOS. Figure 10. Basins of attraction for Newton’s method on SRK EOS.

Table 2. Critical Properties for EOS Example property

N2

C14H30

Tc (K) Pc (bar) ω

126.1 33.496 0.045

694.0 14.192 0.6797

one of the complex-valued high temperature/high conversion solution. From this initial value, T0 ) 402.844 + 14.925i, Newton’s method easily converges the solution, T2 ) 419.239 + 2.1406i, in 10 iterations to a convergence tolerance of 10-8. If instead, we select the midpoint of the line connecting inverse images 1 and 2 (i.e., 419.170 and 388.257 K) as a starting point, Newton’s method converges to the real-valued solution, T1 ) 298.429 K, in 25 iterations to the same convergence tolerance. Homotopy-continuation methods would use far more iterations (or function evaluations) to trace the S-shaped curve shown in Figure 6 and would still only find the real-valued low temperature/low conversion solution. 5.2. Cubic Equation of State. Similar numerical experiments were conducted for the compressibility form of the Soave-Redlich-Kwong (SRK) equation with polynomial mixing rules and a binary interaction parameter of kij ) 0.007 for the binary mixture of 60 mol % nitrogen and 40 mol % C14H30 at 280 K and 3.65 bar. The critical properties used in these calculations are shown in Table 2. For the above conditions, there are three compressibility roots, z1 ) 0.030 12, z2 ) 0.484 94 + 0.098 7i, and z3 ) 0.484 94 - 0.098 7i, two real-valued singular points at z1′ ) 0.193 13 and z2′ ) 0.473 81, and a critical point (i.e., a point at which F′′(z,P) ) 0) at zc ) 0.333 33. Figure 10 shows the basins of attraction for Newton’s method. Note that the Julia set, which is represented by the black dots in this figure, is dispersed throughout the basin boundaries along the real line and in the complex plane. Julia Set Approximation. To generate an approximation of the Julia set for Newton’s method, we used reverse Newton iteration in exactly the same manner as was done for the one-dimensional CSTR example. That is, the reverse Newton iteration was

defined by

GRN(zk-1,P) ) zk-1 - F(zk-1,P)/F′(zk-1,P) - zk ) 0 (7) where the unknown variable is the compressibility, z, and the pressure, P, is the parameter. The corresponding analytical Jacobian matrix for the reverse iteration function given in eq 7 is

J(GRN) ) F(z,P)F′′(z,P)/[F′(z,P)]2

(8)

Figure 11 shows the approximate Julia set generated by solving 30 reverse iteration problems starting from the singular points z2′ ) 0.473 81. Note that the approximate Julia set is a remarkably “good” gross approximation of the true Julia set (by overlaying Figure 11 on Figure 10); it clearly aligns with the basin boundaries throughout the complex plane. From Figure 11, it is certainly possible to get a sense of the structure of the basin boundaries. In particular, the absence of members of the Julia set to the left of the singular point at 0.193 13 as well as the black dots between the two singular points clearly suggests that the liquid root is real-valued and the vapor-like roots are complex-valued. Moreover, the points dispersed in the complex plane that form a shape like the tail of an airplane indicate that this is a region in which there is extreme sensitivity to initial conditions for the Newton’s method. Generating Initial Values. From an initial guess of z0 ) 0.4, our trust region method converges to the singular point z2′. At this point, we use this singular point and reverse Newton iteration to generate, in order, the inverse images 0.473 806, 0.559 756, 0.656 300, 0.415 505, 0.525 287, and 0.298 985 + 0.109575i. These are shown in Figure 11 and labeled points 1, 2, ..., respectively, where point 1 is the attracting singular point z2′. Note that most of the inverse images are on the real line while the last one is complex-valued. The important things to note here are that these few inverse images are dispersed throughout the basin boundaries and that any line connecting any two of these members of the Julia set contains points that are in the basins of attraction for at least one solution to the SRK equation. For example, the line connecting inverse images 1 and 2 contains points that are in the basin of attraction for the real-valued liquid compressibility root z1, while the line connecting inverse images 1 and 6 contains many


points in the basin of attraction for the complex-valued vapor compressibility root. Similar observations can be made for the lines connecting any other two inverse images (i.e., members of the Julia set). Now all that remains is to use a point on the line connecting two members of the approximate Julia set as an initial value for Newton’s method to find a solution. For example, the midpoint between the first and second inverse image gives the initial value z0 ) 0.516 781. Using this initial value, Newton’s method finds the real-valued liquid compressibility root to a convergence tolerance of 10-8 in 11 iterations. If instead we use the midpoint between the first and sixth inverse image, the resulting initial value is z0 ) 0.386 395 5 + 0.054787i and Newton’s method converges to the complex-valued vapor-like compressibility in seven iterations to a convergence tolerance of 10-8. Similar results have been obtained for the midpoint between any two real-valued inverse images. 5.3. VT Flash Problem. Consider the fixed vapor, isothermal (VT) flash problem studied by Lucia and Liu15 involving a mixture of CO2 and n-butane at 344.26 K which contains a retrograde loop in the pressurecomposition (Pxy) diagram, where x and y denote equilibrium liquid and vapor composition, respectively. A complete description of the model equations and unknown variables as well as the Pxy diagram can be found in Lucia and Liu.15 This mixture has a critical point (or turning point), which corresponds to the transition between one- and two-phase behavior, at an overall composition of z1 ) 0.7855 and a pressure P ) 69.40 bar, where the subscript 1 denotes carbon dioxide. From the continuity of first partial derivatives there must be a continuous collection of nondegenerate singular points to the left and right of the turning point that passes through the turning point. These singular points must be nondegenerate by simple tie line arguments. For feed compositions to the right of the critical point, say z1 ) 0.79, and vapor-to-feed and temperature specifications of V/F ) 1 and T ) 344.26 K, respectively, the true solution is single phase behavior. Consequently, the two-phase flash model is infeasible and exhibits a real-valued local minimum (or nondegenerate singular point) at y1 ) 0.79, x1 ) 0.5921, and P ) 71.67 bar and a pair of complex-valued solutions at y1 ) 0.79, y2 ) 0.21, x1 ) 0.5634 ( 0.0775i, x2 ) 0.4366 ( 0.0775i, L/F ) 0, and P ) 70.101 ( 7.303i bar. It is important to understand that there are no real solutions (not even a family of trivial solutions) when complex domain methods and the compressibility root assignments suggested in Gow et al.16 are used. This is because these root assignment rules force the vapor and liquid to have distinctly different compressibilities (and hence equilibrium ratios unequal to one) at each iteration and annihilate trivial solutions. Thus for these specifications, the only stationary point in the real domain is a local minimum in the two-norm and the only solutions are complex-valued solutions! For many real and complex starting points, pure complex domain Newton’s method fails to solve this problem and settles into a stable complex-valued periodic cycle. For example, from a starting point of y1 ) 0.79, y2 ) 0.21, x1 ) 0.55, x2 ) 0.45, L/F ) 0, and P ) 61.25 bar, pure Newton’s method is attracted to a stable complex-valued four-cycle in 52 iterations. On the other hand, our trust region method requires just five iterations from the same starting point to find the

nondegenerate singular point. Real domain homotopycontinuation methods must trace out the entire phase envelop only to conclude that the solution is in the single phase region. Julia Set Approximation. To generate an approximation of the Julia set for Newton’s method, we used multivariable reverse Newton iteration defined by

GRN(Zk-1,p) ) Zk-1 - [J(Zk-1,p)]-1F(Zk-1,p) - Zk ) 0 (9) where the unknown variable vector Z ) (x1, x2, y1, y2, L/F, P) and the feed composition, z, is the parameter. Because this is a multivariable problem, the Jacobian matrix, J(GRN), defined in eq 3 was approximated by using a Broyden (or Schubert) update. Starting from the singular point at y1 ) 0.79, x1 ) 0.5921, and P ) 71.67 bar, reverse Newton iteration finds both real- and complex-valued inverse images. Because the component mass balances and the vapor flow specification equation are linear, y1 ) 0.79 and L/F ) 0 for all reverse Newton iterations and we can denote the inverse images by the shorthand notation Zk-1 ) (x1, P). The first portion of the sequence of real-valued inverse images calculated by reverse iteration is {Zk-1} ) {Z-1, Z-2, Z-3, ...} ) {(0.491 959, 63.5863), (0.396 183, 54.0112), (0.580 01, 71.6012), ...}. The corresponding sequence of complexvalued inverse images is {Zk-1} ) {Z-1, Z-2, Z-3, ...} ) {(0.664 720 + 0.0717037i, 78.8369 + 5.07706i), (0.639 501 + 0.0664465i, 83.4921 + 8.06133i), (0.750 010 + 0.10013i, 93.0021 + 0.50020i), ...}. It is easy to see that these inverse images are dispersed throughout the basin boundaries. It is also important to understand that any reverse iteration problem can be initialized to deliberately look for both real- or complex-valued inverse images. Continuing backward for as many reverse iterations as desired, we generate a reasonably good gross approximation of the Julia set. Generating Initial Values. Starting from y1 ) 0.79, y2 ) 0.21, x1 ) 0.55, x2 ) 0.45, L/F ) 0, and P ) 61.25 bar, our trust region strategy finds the real-valued attracting singular point at y1 ) 0.79, x1 ) 0.5921, and P ) 71.67 bar in just five iterations. From here, reverse iteration generates the pair of inverse images Z-1 ) (0.491 959, 63.5863) and Z-1 ) (0.664 720 + 0.0717037i, 78.8369 + 5.07706i). Using the midpoint between the singular point and each of these inverse images as starting points, Newton’s method converges to the complex-valued solution in eight and seven iterations, respectively, to a tolerance of 10-6. 6. Concluding Remarks It was illustrated that singular points and points in the Julia set occupy the boundaries between different basins of attraction of a given iterative map. The closure of Julia sets was used to develop a method for building approximations to the Julia set and basin boundaries for a given map in process simulation by reverse iteration. It was shown that reverse iteration is capable of generating very good gross approximations to the Julia set and basin boundaries, from which initial values that converge to a desired solution can be chosen. It was also shown that multiple inverse images to the reverse iteration problem can exist and that this can result in a tree of inverse images and computational demand that is combinatorial. A simple strategy for choosing initial values that converge to a solution from


a small set of inverse images was proposed and shown to be effective in finding solutions to process model equations. Several chemical process examples and associated geometric illustrations were used to illustrate the underlying concepts of reverse iteration for generating both Julia set approximations and initial values that converge to solutions to the model equations. Other strategies for determining good initial values can also be used. In particular, we have recently switched to a fractal interval method on the line connecting an attracting singular point and any inverse image in conjunction with a predictor-corrector strategy to improve the quality of initial values generated from reverse iteration. Finally, we note that, unlike homotopy-continuation methods, the proposed strategy for generating initial values is capable of finding path-disconnected solutions such as those described in Russo and Bequette17 in which turning points can exist outside the (feasible) region for which the model is defined or on isola that exist in the neighborhood of an S-shaped curve. Literature Cited (1) Levy, A. V.; Montalvo, A. The Tunneling Algorithm for the Global Minimization of Functions. SIAM J. Sci. Stat. Comput. 1985, 6, 15. (2) Sridhar, L. N.; Lucia, A. Process Analysis in the Complex Domain. AIChE J. 1995, 39, 461. (3) Julia, G. Memoir sur l′Iteration des Function Rationelles. J. Math Pures Appl. 1918, 4, 47. (4) Fatou, P. Sur les Equations Fonctionelle. Bull. Soc. Math France 1919, 47, 161.

(5) Mandelbrot, B. B. Fractal Aspects of z ) λz(1 - z) for complex λ and z. Ann. N.Y. Acad. Sci. 1980, 357, 249. (6) Barnsley, M. F. Fractals Everywhere; Academic Press: San Diego, CA, 1988. (7) Nusse, H. E.; Yorke, J. A. Basins of Attraction. Science 1996, 271, 1376. (8) Lucia, A.; Liu, D. The Dynamics of Equation-Solving. In Fractals and Chaos in Chemical Engineering; Giona, M., Biardi, G., Eds.; World Scientific: London, 1997. (9) 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. (10) Lucia, A.; Guo, X.; Wang, X. Process Simulation in the Complex Domain. AIChE J. 1993, 39, 461. (11) Lucia, A.; Guo, X. Feasibility in Chemical Process Design and Operations. Comput. Chem. Eng. 1996, 20, S569. (12) Peitgen, H. O.; Saupe D.; Haeseler, F. v. Cayley’s Problem and Julia Sets. Math. Intelligencer 1984, 6, 11. (13) Venkataraman, S.; Lucia, A. Exploiting the Gibbs-Duhem Equation in Separation Calculations. AIChE J. 1986, 32, 1057. (14) Smith, J. M. Chemical Engineering Kinetics; McGrawHill: New York, 1981. (15) Lucia, A.; Liu, D. An Acceleration Method for Dogleg Methods in Simple Singular Regions. Ind. Eng. Chem. Res. 1998, 37, 1358. (16) Gow, A. S.; Guo, X.; Liu, D.; Lucia, A. Simulation of Refrigerant Phase Equilibria. Ind. Eng. Chem. Res. 1997, 36, 2841. (17) Russo, L. P.; Bequette, B. W. Impact of Process Design on the Multiplicity Behavior of a Jacketed Exothermic CSTR. AIChE J. 1995, 41, 135.

Received for review April 30, 1998 Revised manuscript received August 6, 1998 Accepted August 7, 1998 IE9802764

Reverse Iteration in Chemical Process Simulation - ACS Publications

Recommend Documents