Computers them. Printed
in Great
Engng, Britain.
Vol. 14, No. 8, pp. 889-900, 1990 All rights reserved
DIFFERENTIAL
Copyright
GEOMETRY BASED CONTINUATION
0098-I 354/90 53.00 + 0.00 0 1990 Pergamon Press plc
HOMOTOPY
W. L. RION and V. VAN BRUNT Department (Received
of Chemical
23 August
Engineering,
1989; final revision
University
received
of South Carolina,
2 September
Columbia,
i 990; receivedfor
SC 29208,
publication
U.S.A.
27 March
1990)
Abstract-An improved homotopy continuation algorithm for solving systems of non-linear equations has been developed. The improvement in this algorithm is that in addition to the unit tangent, the curvature and principal unit normal of the solution path Frenet frame are calculated at each continuation step. The curvature and principal unit normal are then incorporated into the prediction and step length control phases to yield a faster more efficient continuation algorithm. Increased computational load over standard procedures is minimized by using the inverse Jacobian matrix required for tangent vector calculation in standard procedures, and a single row/column bordering technique. to obtain the augmented Jacobian matrix inverse required for calculation of the curvature and principal unit normal. The new algorithm is tested on several small, but highly non-linear problems and the results are documented. The new algorithm requires less continuation steps, less Newton corrections, and encounters fewer poorly-chosen step sizes. CPU time requirements are reduced for most problems.
1NTRODUCTION
Many chemical engineering problems can be represented mathematically by a set of non-linear algebraic equations. A generalized set of such equations is given by:
.L(x,,x,,.
. .,-x,,)=O.
(1)
The solution to this problem, written in simplified form as F(X) = 0, requires determination of the vector X such that F(X) = 0. (Here 0 represents the zero vector.) The most common approach to this type of problem is to use Newton’s method which can be written in general form as: xk+‘=xk_
E [
ax
3%. -w
dxl
I 1 dF ax=
8x1
1 F(X), -’
(2)
. . . -aA_ ax,
(3)
With this method X” is the current (or initial, when k = 0) estimate of the solution vector which yields F(X) = 0. Equation (2) is applied with an updated Jacobian matrix (3) at each point until a criterion 889
such as )IX*+‘-XXkII -z L, where L represents the desired solution tolerance, is met. An inherent disadvantage of Newton’s method is that topological phenomena which generate singularity, or near-singularity, of the Jacobian matrix can prevent convergence. Newton’s method is therefore only locally convergent and relies on initial estimates close to the solution for convergence. However, a large or even global convergence domain is often desirable. This is especially true when the algorithm is part of a large simulation or expert system where a good initial estimate may be. difficult to provide and failure to converge would interrupt progress and waster computer time. A method which allows for a greatly expanded convergence domain relative to Newton’s method is differential arclength homotopy continuation. Wayburn and Seader (1987) discuss the conditions under which homotopy continuation is guaranteed to converge, or not converge. Additional advantages of homotopy continuation methods include the ability to systematically search for multiple solutions and to easily determine the dependence of the solution on problem parameters. Because of these features, application of homotopy continuation to engineering and physics is widespread. Seydel and Hlavacek (1987) have recently published a review article on the role of continuation in engineering analysis. Another excellent review of continuation methods and their application to separation problems is given by Wayburn (1987). Several published computer codes exist for the implementation of homotopy continuation to engineering systems. (Rheinboldt and Burkardt, 1983a, b; Kubicek and Marek; 1983; Frantz, 1986; Watson et al., 1987; Morgan, 1987). The algorithm of Frantz,
W.
890
L. RION and V.
which is especially adapted to engineering systems, uses the local continuation variable procedure of Rheinboldt and Burkardt (1983) during the predictive phase of the algorithm. Newton corrections are made in the plane orthogonal to the tangent vector. Quasi-Newton or Broyden corrections may be utilized in place of, or in conjunction with, rigorous Newton corrections. The step-length control is that of Georg (1980). Additional benefits of this algorithm include: automatic incorporation of system equations into common homotopy formulas, use of analytic (through user-supplied subroutines) or numerical (no user input required) derivatives, and ability to continue integration after location of target points or roots. Frantz (1986) and Frantz and Van Brunt (1987) have compared their algorithm to those of Kubicek and Marek (1983) and Rheinboldt and Burkardt (1983a, b) on several test problems. The algorithm of Frantz and Van Brunt performed as well or better on all of the problems considered in terms of number of continuation steps required, number of function evaluations required, and CPU time. In this work the algorithm of Frantz and Van Brunt has been improved by including procedures for calculating the curvature and principal unit normal of the homotopy path Frenet frame from the system equations. This information is then incorporated into the prediction and step-length procedures to yield a faster, more efficient continuation algorithm. The new algorithm requires fewer continuation steps, fewer Newton corrections, and encounters fewer poorly-chosen step sizes. CPU time is also reduced. This is demonstrated by comparison of the new algorithm to the algorithm of Frantz and Van Brunt on several test problems.
VhN
BRUNT
The goal of homotopy continuation is to trace the solution path from X at t = 1 to X at t = 0. Under reasonable assumptions (see Waybum and Seader, 1987) this path will exist. The path is followed by differentiating the equation set with respect to the arclength of the path to obtain an IVP, using Eulertype integration to approximate a forward point on the path, using Newton’s method to correct back to the exact path, and repeating until the end of the path is reached. Euler prediction The set of homotopy equations, formed by transformation of F(X) to H(X, t) can be written in genera1 form as: h,(x,rx~r....x,,,t)=O hr(xI,xrr...rx,,t)=O
h,(x,,x,
,...,
(4)
x,,t)=O.
The parameterization of the equations is most naturally by the Euclidean metric arclength S. To form the IVP, equation (4) is differentiated with respect to the solution arclength to obtain:
1 Bh, ax,
-
dh, ... -
ah,’ -
dx,
at
ds
ax,
ah, . . . ax,
ah, at
dx, ds
dh, ax,
.‘.
dh, ax,,
ah, dt
dx, cl.9
ax,
dt L
ds
= =I 7
0
0
0
(5) REVIEW
The homotopy
OF CONTINUATION
PROCEDURES
map
A homotpy h(X, t) is a blending of two functions, f(X) and g(X), by a parameter t. The homotopies discussed here belong to one of two classes. The first type are the convex linear homotopies by: h(X, 2) = (1 - t)f(x)
+ tg(X).
Here f(X) is the problem whose solution is desired and g(X) is a “simplified” version of the problem with a known solution. At t = 1, h(X, t) =g(X), the easy problem, while at t = 0, h(X, r) =f(X), the real problem. The second type of homotopy considered are those formed by imbedding a parameter into the original equations f(X) or by considering one of the original problem parameters as the imbedded parameter. In the convex linear homotopies, solutions at intermediate values of the parameter r are usually not of interest. However, for problem parameter imbeddings, the entire solution path is of interest.
This system is made determinant by addition of the definitive equation for arclength: ds*=dx;+dx$+
. ..+dx.+dt*,
which can also be written as: * : 2+!?+ . . .+~+!& This yields the following determinant system: ‘-.I ah ax, dh,
ah, ax,
ax,
ax,
2ah ax, dx, ds
dx,
...
ah, dx,
..’
:
.
ah, ax,
...
ah, dx,
dh, ar
dx,, ds
0
dx, -dY
...
dx, ds
dt Z
dr Z
1
ah,
dh; at
ds
dh,
ah,
dx,
ax,
dr
ds
0
0 I
(6)
Differential
This N + 1 x N + 1 system is solved by reducing it to an N x N system using either substitution and elimination or matrix elimination procedures described in Frantz (1986). The last row, and the column corresnonding to- the imbedded 138.rameter t are eliminated- to yield:
-ah -.L
dh --!
ax,
dx,
ah,
ax”
ds
at
.,.
ah, ax,,
dx, ds
_,.
2
ah
dxn
a-r,,
ds
..,
as>
dh2 ax,
dh, dx,
I
-ah,
= -
ah 2 ax, a+ ah
dt
ah,
at dh
_LL
1
=
T-t II ’
where
=
.
(
order
to
avoid
matrix in equation (7),
techniques to reduce the N + 1 x N + 1 system in (6) to an N x N system by eliminating the last row and the column corresponding to the local continuation variable. The local continuation variable is chosen as the one whose component of the tangent vector from the previous iteration has the largest absolute value. Consider the kth variable of equations (4) to be the continuation variable; then (6) can be reduced with respect to the continuation variable to the following N x N system: ah,
Z
ah,
ax, ...
ah2
dh2
ds,
a.Y,
of Newton’s
in one of two ways.
where e’=[O,O,O
. . . . Ii,0 ,...,
ah,,
ah2
dx,
ah, -
ds
ax, =-
ah2
dx,
d++r ds
ah,,
al, _,
)I
ah
axk+l
_..
ah
2
34,
d-x,,+ ,
14 X--F
ah
2
ax,
ds dxI, ds
CACC
0]
and i represents the variable with largest gradient to the homotopy path.
ds
_,
method
: ah,,
(9)
[e’. (X - X’)] = 0,
d-T
ah,
.
ax, ax, .
.
The first method is to hold one of the variables constant. This yields an N x N determinant system. Newton iteration is performed on the remaining variables. Rheinboldt (1980) recommends that the variable with the largest gradient be held constant to prevent numerical instability. This can be accomplished by supplementing the equation set by the equation:
dx,-, ds
.
>
procedures
ax,_, a_~,
dt2 112
dx; “‘ds+Z
Once an estimate to the homotopy path is generated as above, Newton’s method is employed to correct back to the exact homotopy path. The homotopy equations form a set of N equations in N + 1 unknowns. The system is made determinant during application
singularities in the Jacobian as turning points are approached, the local continuation variable technique of Rheinboldt and Burkardt (1983a, b) (and used by Wayburn and Seader; 1987; Frantz and Van Brunt, 1987) can be used. This procedure utilizes algebraic In
. ..!!!I!$+. d-d+,
Correction (7)
IIT-(II
,,,_,,,=(~+~+
The system shown above, along with the initial solution point at I = 1, form the initial value problem desired. The elements of the Jacobian are either approximated numerically or supplied through analytic subroutines. The unit tangent is calculated by inverting the Jacobian matrix in (8) and subsequent use of equations (8) and (9).
at
ds
with I(T_* ([ being the norm of the reduced tangent vector defined by:
+-+ ds
.
_2
891
geometry based homotopy continuation
=
1
IIT-k II’
W. L.
892
In the second method the equation plemented by the equation:
RION
and V.
set is sup-
[T.ul=o, where U is the vector of Newton correctors. This is a geometric constraint that forces corrections to be made in the plane orthogonal to the tangent vector (see Allgower and Georg, 1980). Step -length conrrol
Homotopy continuation procedures are most efficient if step-length control is applied to the predictor. The goal is to take as large a step as possible while remaining in the domain of convergence of the Newton corrector. It is also important not to jump to other paths or to skip portions of the homotopy path that are of interest. Den Heijer and Rheinboldt (1981) showed that the convergence radii of a forward point on the homotopy path cannot be deduced from previous convergence radii. Deuflhard (I 979), Den Heijer and Rheinboldt (1981) and Rheinboldt and Burkardt (1983a, b), proposed step-length controls based on interpolation polynomial approximations to the exact curve. Watson et al. (1987) used an adaptation of the procedure of Rheinboldt and Burkardt. The algorithm of Georg (1980) has been adopted by many researchers including Wayburn and Seader (1987) and Frantz and Van Brunt (1987). In this algorithm the primary step-size determinant is: As,=As,_,p where: As, _ , = step-length used for previous continuation step, Asj = step-length for next continuation step, P = @i&all@? Oldea,= ideal turning angle, user-specified, 0 = local estimate of the turning angle, given by: O=cos-‘[T,_,.T,] and T, = unit tangent T,_ t = unit tangent
VAN BRUNT
the same path, is suspected and the step-length is reduced. Contraction ratio. If the ratio I(F(X’) I(/ 1)F(Xk _ ‘) 11is not less than a user-specified number, then poor convergence behavior is suspected, the step is rejected, and the steplength reduced. Turning angle. If the angle between the tangent vectors at the current point and the previous point is greater than a specified value, then convergence to another path or another portion of the homotopy path is suspected, and the step to the current point is rejected, the step-length reduced, and a smaller step attempted. Maximum and minimum step-length. A maximum limit is placed on the step-length to avoid jumping to other homotopy paths. If the steplength is reduced below a specified value by the other controls, the entire program is re-started with the turning angle and maximum steplength parameters reduced. DIFFERENTIAL
HOMOTOPY
The modified continuation algorithm to be described here is based on the local differential geometry of the solution trajectory. Thus we review the differential geometry concepts used. First the Frenet frame in the classical 3-D case is reviewed (see Stoker, 1969). The extension of the Frenet frame to n dimensions is then discussed. Afterwards it is shown how these concepts are incorporated into the homotopy continuation algorithm. Curvature
Consider
and forsion in R’
a curve C(s) R’(C(s)
in
= [x,(s)> XZ(S),
xj(s)f)
parameterized by arclength. The basic result of the Frenet equations is that there exists an orthonormal frame along C(s), consisting of the mutually orthogonal unit vectors T, N and B which are related by the following equations:
at current step, at previous step.
T’(s) = N’(s)=
The step-length, As, can be limited by the following secondary step-length control parameters: 1. The number of Newton corrections required to converge to the exact solution from the prediction. If the number of iterates exceeds a specified value the step is rejected and the step-length is reduced. This helps prevent prediction beyond the domain of attraction of Newton’s method. 2. Initial size of the corrector step. If the length of the Newton correction vector during the initial step of the correction process is above a specified value, then convergence to another strongly attractive path, or a different part of
GEOMETRY BASED CONTINUATION
KN -KT+
B’(s) =
rB,
--rN,
where; T = unit tangent vector, N = principal unit normal B = unit binormal vector, K = curvature, 5 = torsion. K and
vector,
T completely define the curve to within a proper Euclidean motion. N = C”(s)/ Here T = C’(s), the unit tangent. )(C”(s) (1 is defined as the principal unit normal and
Differential geometry based homotopy continuation is everywhere perpendicular to T since [T . T] = 1, implies lT . T’] = 0. B = T x N is the unit binormal and is everywhere orthogonal to T and N. This shows that T, N and B form an orthonormal frame along C(s). The curvature and torsion can be calculated by the following formulas: K=
11 C”(s) Ij = jlC’(t)
r = -[B’(s)
= [[c’(t)
cos-‘[T(s)
A’+B2
Bt
T
,
ds
and
Fig. 1. Helix curve. Generalization
of the Frenet equations
T’(s)
N(s) = -_
K,(s)
(10)
T’(s) Now
[T.N]=O,
N’= B; = B; =
= K, (s)N(s).
When combined
with the above result: [N’ . T] + K, = 0.
This implies that N’ = - K,T + a vector perpendicular to T and N. The second curvature function is defined by: K2G)
=
II N’(s)
+
K,
(s)‘Vs)
//
and if K2 # 0 for all s we let: B (S) = N’(s)
I
+ K,(s)T(s) K,(s)
Thus the curvature and torsion are constant for the helix. In light of the previous discussions this T’ =
equation is formed:
so [N’.-l-+[T’-w=O.
B
r = A2+B2.
to RN
The Frenet equations are generalized to RN as follows (see Spivak, 1979): Let C(S) : s +R”, be a curve parametrized by arclength. As before we define C’(S) = T(s). Since C is parametrized by S, 11 C’(S) )( = I]T(s) 1)= 1. Thus as before [T.T]= 1, and [T.T’]=O, so T and T’ are perpendicular. The first curvature function is defined by I/C”(s)l( = ((T’(s)(( =K,(s). If K, # 0 for all s we let:
Thus the first “Frenet”
. T(s -t ds)]
i.e. K is the rate at which the angle between tangent vectors is changing with respect to the arclength. The vector B, formed by T x N, is obviously perpendicular to T and N, and is a unit vector. Geometrically the torsion is the rate at which the curve is turning out of the osculating plane. The torsion 5 may he positive or negative. To illustrate these concepts consider a helix curve defined and illustrated in Fig. 1. The curvature and torsion for the helix curve can be computed directly as: A
l
. C”‘(t)]/I/ C’(l) x C(Z) 112.
d -0
K=-
X3(t)
N’(s)]
The expressions above containing s are for unitspeed, or arclength parametrized curves while those containing t are for non-unit-speed, or arbitrarily parametrized curves. Geometrically, T, N and B form a triad of orthogonal unit vectors along C(s). N is a unit vector which points in the same direction as the vector whose components are d2xi/ds2. The length of the vector d2C/ds2 is equal to the curvature. The plane defined by T and N form the osculating plane of the curve. The curvature can be shown to be the radius of curvature of the circle which best fits the projection of C(s) in the osculating plane. By definition the curvature is always positive. The curvature has another important geometric interpretation: K = lim
Xl(t) = A coa(t1
X2(t) - A sin(t)
x C”(t) l]/ljC’(t) lj3,
. N(s)] = [B(s)
x C(l)]
893
.
B, is a unit vector field everywhere perpendicular to T and N. Thus N’(s) = -K, T + K2B,. Continuing in this fashion, the generalized Frenet equations take the form:
K, N, -K,T
+
-K,N
K,B,.
+ -KK,B,
B; = might be intuitive by inspection of the helix illustration. The radius of curvature of the circle approximating the helix in the osculating plane is constant as is the rate at which the curve turns out of the osculating plane.
K,B,, + -K,+iB,
K,B,,
I +
K/+2Bj+I
until N is reached. At N we have: B&-, = -K,y_
,B,_,
+ a vector
perpendicular
to T, N, B, , . . , B, _ ?.
W. L. RION and V. VANBRUNT
894
At this point T, N and B, +B,_ 2 form an orthonormal basis to RN, thus only the zero vector is perpendicular to T, N and B, +BN_ t. So B&_*= -KN_,BN_-) and K,,_, = -[Bh_,-B,_,J= [BN--2.Bp_3] is the form of the final generalized Frenet formula. Thus only the last curvature function is allowed to take on negative values. The above discussion shows that for an N-dimensional problem, N - 1 curvature functions are needed to completely define the curve. It follows that in order to calculate the N - 1 curvature functions the (dC/ds - dNC/dsN) vector fields must be known along C. Incorporation
of curvature
and unit normal into the continuation
algorithm
Calculation of the curvatures and principal unit normal for the system H(X, t) = 0 requires that equations (4) must be differentiated twice with respect to arclength. The second differentiation of hi(X) yields equations of the form:
where x,+ 1= t. Equations (1 I) can be represented in matrix form ai5: dh,
dx,
ah,
ax,
ah,
ah,
ax,
ax,
ah,
... ax,
ah, ax,+,
d’x,
ah,
d2x, ds2
ah, ..’
dx,
Zl
d&r2
ax,,,
=2
= ah”
ah,
dx,
Bx,
ah -2! ax,
-.a
d2x n ds2
ah” ax,+*
=!I
d’x, + , ds2 where zi=c
“+‘“+’ a2hi dx,dx,_ c j=, k=, axjaxk ds ds
To make the above system determinant the system is supplemented with the orthogonality constraint of T and N:
This yields the following, determinant N + 1 x N + 1 system: -ah, ax, ah2
ah, -
ax,
ah,
. . .
dh, -
ax,
ah* ... ax,
-
ah,
ax,+,
ah, ax,+,
ax,
ax,
:
.
ah _! ax,
ax,
..’
G
ax,+,
dx, ds
.
-dx,
-dx,+ I
dx,
_ ds
ZI
d2x, ds2
22
=
ah,
-
d’x --! ds2
ah,
ds
ah”
ds
d*x 2 ds2 d2x, + ,
ds2
Thus:
KN(s) = [Z(s)l[J+(s)l-‘3
.
Z#l
0
.
Differential
geometry
based homotopy
ah,
ah,
ax, ax* . ’ ah,
3x1
ah2
ax, . .
895
partial derivatives appearing in the Jacobian may change at each step as the continuation variable changes. This Jacobian inverse is shown schematically below:
where
ah,
continuation
ah,
ax, ax,+, ah, ah, -ax,, +
ax,,
I
=
[J+(S)]-‘,
hk+,
. .
.
jzn
ai, ah ah,, ai, ax, ax, ... ax, dx,,, dx --! ds Z(s) KN(s)
dx2 ds = rz,
“.
dx,
dt
ds
ds
(Note that the column corresponding to the continuation variable is missing.) The above matrix can be bordered with the column of partial derivatives with respect to the continuation variable and the row to include the orthogonality constraint:
122, . . . , z,,,01=,
= [d2x,/ds2,.
. , d2x,/ds2, d2t/ds*f’.
r jl&-, hk-,
j,k+
1
jzk+,
_i#k
+
,
The principal unit normal vector is found by normalizing the above vector KN. The magnitude of KN is the curvature. Obviously these calculations require prior calculation of the unit tangent vector, the quantities z,, and the inverse of the augmented Jacobian [J’]. The unit tangent is determined as shown previously. The zi are calculated using the components of the unit tangent and either numerically or analytically supplied partial derivatives of the form:
Jacobian
Jl”+l
_izn
hn*,
..
J””
ah#.
dx,
ax,
9
variable
and
forj=l,n+l.
Jt,+ r.1 = ds’
The remaining elements are those of the inverted Jacobian. The bordering method from Faddeev and Fadeeva (1963) can now be employed to find a “scrambled” version of the inverse of [J’]. The elements of [J+]- ’ are given by the following relations:
+‘-L’=ji,k+*i~k/u J‘i.rc
i,k
=l,n,
+(--I) J‘i.n + I = _,+f#;/Q, J..+(--l’ + ,,k=
inverse
--&/a
(i,k = I,n),
j,::::;.!r + I = l/a. The elements vector/matrix ilk
+
,
j2k
+
I
tl/;and ri are calculated multiplications:
.
.
ilk
j21
. .
.
j=-,
j2k
Jut
...
Jnn-I
Jn/i + I
-j,,
.
=
i = 1, n, k = continuation
At each step we have at our disposal the inverse of the reduced Jacobian matrix [J] in equation (8). The
[51,52,---,5,,1=[j,,+1.1
Jnn+ I
/i.k
The inverse of the augmented Jacobian [J’]-’ is calculated with minimal computations by utilizing the previously calculated inverse Jacobian required for unit tangent computation and matrix bordering techniques discussed by Faddeev and Fadeeva (1963).
of augmented
Jl”
. . .
where
as ax,ax;
Cafcuiarion
...
-
I
ilk+
1 +
I
..j.,+,.k-,,j,,+,.k+,...j,,+,.,,+11
by the following
896
W. L. RION and V. VAN BRUNT
The value a is calculated from either of the two following relations (or both as a check on previously calculated 5 and $ values):
I-
I
Correction
I=,
Multiplication of Z by [J+]-’ normal vector shown below:
to other paths (or other parts of the same path), the five secondary step-length control parameters used in the algorithm of Frantz and Van Brunt are used here also.
gives the scrambled
procedure
In the new algorithm Newton corrections are made in the plane orthogonal to the tangent vector as described earlier. APPLICATIONS
= [J+]-’ x Z.
Thus the above method preserves the local continuation variable technique to avoid Jacobian matrix singularities as well as minimizing additional computations for the augmented Jacobian inverse. Incorporation
of curvature
into
predictor
Standard continuation procedures use the unit tangent vector to predict forward points on the homotopy path. Here the unit tangent and principal normal vector are used together to yield a predictor of the form: xi(s + As) = x,(s) + (dx;/ds ), As + (d2x;/ds2)., (As)‘/2, which can be written as: X(s + As) = X(s) + T(s)(As)
+ K(s)N(s)
(As )* -.
This confines the predicted point to the plane formed by T and N and is essentially a parabolic approximation of the actual curve in the osculating plane. The error of the new predictor is of order (As)~ instead of (As)’ as in the case of the tangent-only predictor. Incorporation cedure
of‘ currlature
into
the
step-length
pro-
Since the curvature and turning angle are related by equation (IO) it is natural to define the optimum step-length based on the turning angle desired during the continuation steps. At each step the curvature is calculated and the step-length for the predictor is calculated by:
Oidca,is a user-specified parameter which represents the “ideal” turning angle between tangent vectors, large enough to insure adequate progress along the homotopy path and small enough to prohibit prediction into the attraction of other paths. The above equation defines the primary means of determining the step-length. To avoid predictions beyond the domain of attraction of Newton’s method, or jumps
AND
DEMONSTRATION
The new algorithm incorporating the solution curvature into the predictor and step-length was compared to the algorithm of Frantz and Van Brunt on several small, but highly non-linear test problems. The measures of performance were chosen as: I. The number of continuation steps required to trace the desired continuation path. 2. The number of Newton correctors required to correct the predicted values back to the exact homotopy path. The number of predictor step-length rejections due to large initial contraction size, large contraction ratio or too many iterates during the correction process, or due to large turning angle. Number of Jacobian matrix inversions. Number of function evaluations. CPU time. The first two measures reflect effectiveness of the predictor and step-length. The third is mostly a measure of step-length performance. The Jacobian matrix inversions is the sum of the inversions from tangent vector computations and Newton corrections. The number of function evaluations indicates how many times the system equations were called to calculate numerical partial derivatives and/or equation residuals. When comparing algorithms, tuning parameters such as contraction ratio, ideal turning angle, etc. were optimized for both algorithms so that the measures of performance were minimized for a given problem. The differential geometry based algorithm requires calculation of second-order partial derivatives of the homotopy equations. Numerical computation of these partial derivatives can be costly in terms of overall computational load. Therefore the computer code which incorporates differential geometry utilizes user-specified pairs of problem variables j and k where: t??‘h
---2-
axi dzc,
= 0,
for all i,
to avoid function calls to calculate second-order numerical partial derivatives known to be 0 for all homotopy equations (for example d’hJ&’ = 0 for all convex linear homotopies).
Differential geometry based homotopy continuation
897
tank reactors in series at steady state with recycle and an exothermic, first-order, irreversible reaction. The equations describing this system (from Kubicek, 1976) are: R)(l
f, = 0 = (1 -
-yl)
x exp{lOA/(l fz = 0 = (1 -
x _A4 =
0
-8
RWU
-
10(1 + BMI
fx=O=y~ -8
+ 1OAh)l
-~2+(1
-
lmbedded
Parameter,
The independent
t
y, = reactant
-~2)
~W,/Y)I,
lO(1 +/Y&$2
x exp(lO4J(l
+ lOAh)
+ B,#%lr
+
104,
y,)exp(lO4,l(l
-RN
w41W2/U =
--Ye,
+ (1 -
+ lO&lv)3
variables
R)D(l
-y*)
+ MEZ-
in these equations
conversion
are:
in reactor
i = 1 or 2, Fig. 2. Himmelblau’s
4, = dimensionless
function
temperature
in reactor
i = 1 or 2, Himmefbiau
‘s function
R = recycle
As a first example we consider Himmelblau’s function (from Himmelblau, 1972; Reklaitis et al., 1983). The function is:
The physical parameters defined as follows: y = dimensionless
The
stationary
aflax,
points
= f, = 2x:
are roots
+ 2x, x2 -
21x, + x: -
13x* -
D=
4, = dimensionless
1. Solution
pi= dimensionless
No. No. No. No. No. CPU
steps correctors failures functions Jacobian time (ms)
Numeric 40 127 50 961 222 690
Van
temperature
rise,
temperature,
I- 5, 21,
heat transfer coefficient
reactor
i,
P; E IO, 31. In this example, y = 1000, D = 22, I&, = I#+ = 0, R =O+l. j4=/!&=2, The desired solution diagram is that of reactor 2 temperature $J~vs recycle rate R. Solution of this set of equations is very difficult due to strong attraction to a divergent path with negative recycle coordinates. Solution diagrams and statistics are compared in Fig. 3 and Tabte 2.
comparisons
and
coolant 4, E
CSTRs
Derivatives
adiabatic
001,
D E IO, 601,
The next example (detailed in Kubicek et al., 1980) is a model of two continuous, non-adiabatic, stirred
Frantz
dimensionless
11 = 0.
are
energy,
Y E PO,
7 = 0,
For this problem the fixed-point homotopy was used. The method of Kuno and Seader (1988) was followed to determine the initial points as X, = x2 = -50 so that all roots would be found from one starting point. Solution diagrams generated by the two continuation codes are compared in Fig. 2 and Table 1. For this problem the curvature-based algorithm outperformed Frantz and Van Brunt’s algorithm in every category except function evaluations for numerical partial derivatives.
Table
in the above equations
activation
to the equations:
af/ax, = f* = x: + 2x, x* + 2x: -
fnterfinked
R E [0, I].
ratio,
of
Brunt
Himmelblau’s
algorithm
Analytic 40 127 50 517 222 673
function
Curvature-based NUlTXIiC 34 103 19 1147 160 427
algorithm Analytic 34 103 19 349 160 380
898
W.
(l-R), 0
Frank
Algorithm
and Van Brunt (.tw
Rate
R - Recycle # fn bold)
Fig. 3. Interlinked
*
Curvature
Based
Algorithm. Thie Work ,,t.c. It I” Y”d.,ll”.d I1.,80.,
CSTRs.
As in the previous example the curvature-based algorithm required less continuation steps, less Newton corrections, and encountered less correction failures due to large contraction ratio, etc. CPU time was also reduced. Holodniok
‘s function
The next example is from Kubicek et al. (1981). This problem consists of two equations in two unknowns and is known to have 12 solutions: .f, = 0 = 0.5 sin(x, x2) - x&r - x,/2, .fi=O=(l
- 1/4n)(e2X’-e)+exz/~
-2ex,
All 12 solutions can be found using a Newton homotopy and either of the two sets of initial estimates given by Kubicek et al. (a) x, = -0.25 x2 = -e
or
(b) x, = -0.5 x2 = -2e.
To trace the portion of the path containing the 12 roots, the first pair of initial estimates (a) required the normal initial direction of integration, the direction of decreasing t, while the second initial esti-
Table
2. Solution Frantr
Derivatives No. No. No. No. No. CPU
s,elx
c&rectors
failures functions Jacobians time (ms)
commwison
VAN
of interlinked
and Vzin Brunt algorithm
Numeric
Analytic
29
29
60
60
19 718 II4 783
BRUNT
mates (b) required initial integration in the positive t direction. When comparing algorithms, the step-length control parameters were relaxed to achieve minimal computation time until further relaxation resulted in algorithm failure. In the two previous examples, the step-length control parameters yielding the fastest algorithm performance were identical. In the present example, the tuning parameters yielding optimal performance were different. The difference was in the number of Newton corrections allowed during each continuation step. When more than three Newton corrections per continuation step were allowed, Frantz and Van Brunt’s algorithm encountered divergence of Newton’s method when attempting to follow the homotopy path corresponding to initial estimates (a). However, up to eight corrections per step were allowable for the curvature-based algorithm using initial estimates (a). But when using initial estimates (b), the curvature-based algorithm encountered problems with more than three Newton iterates allowed per step. The curvature-based algorithm would turn back on itself at sharp turning points in the homotopy path and retrace a previous portion of the path in the opposite direction when using initial estimates (b) with more than three corrections per step. Frantz and Van Brunt’s algorithm did not encounter problems with up to eight corrections allowed per step using initial estimates (b). Due to the above difficulties the algorithms were compared in two different ways for this example. The first comparison was made using equal tuning parameters (three Newton corrections per step). Then the best case performances of the two algorithms in finding all 12 roots were compared. Solution diagrams and statistics using equal tuning parameters {with three Newton corrections allowed per step) and initial estimates (a) are shown in Fig. 4 and Table 3. Figure 5 and Table 4 compare solution diagrams and statistics using initial estimates (b) and equal tuning parameters. Table 5 compares best-case performance of the two algorithms using different initial estimates and increased Newton corrections in step. The above results show again the increased efficiency of the curvature-based algorithm. The results also indicate that either of the two algorithms may fail when step-length control parameters are relaxed to achieve absolute minimal computation time.
0.9x
0.976
0.956
0.935
V.
L. RION and
19 148 114 667
CSTR
oroblem
Curvature-based Numeric 24 49 10 1173 a7 740
algorithm Analytic 24 49 10 112 87 590
Differential geometry based homotopy continuation Initial
Eatlmaterr
899
Xl
(a)
tnitiat
Estimates
(b)
1.5
1.5
1
0.5
0
-0.5
1
-0.5
-3
-2
-1
lmbedded
Fig. 4. Holodniok’s
1
0
Parameter,
function,
-2
2
initial estimates
Fig. 5. Holodniok’s
(a).
comparisons
Frantz Derivatives
and
steps No. correctors No. failures No. Jacobians time
(ms)
Table
4.
Solution
No. No.
2127
1885
1887
of Van
Holodniok’s Brunt
function
algorithm
Analvtic
Numeric
steps
Solution
and
95 250
156 693 I396 2053
123 582 3357 2013
Holodniok’s
89 234 108 517 874 1658
initial
estimates
function.
best
(b)
algorithm
95 250 123 582 1039 I797
case
step-length
tuning
parameters Frantz
Derivatives No. No. No. No. No. CPU
steps correctors failures Jacobians functions time (ms)
and Van
Brunt algorithm
Numeric 86 254 82 496 2556 I697 Initial
Analytic
estimate
86 254 82 496 1068 1567 (a)
-_
Analvtic
Numeric
107 304
(a)
algorithm Analvtic
using
156 693 3475 2250
of
estimates
Curvature-based
107 304
comparisons
t
function, initial estimates (b).
89 234 108
275 137
comparisons
initial
Numeric
99
517 2961
correctors
5.
Analytic
603 1185
No. failures No. Jacobiam No. functions CPU time (ms)
Table
using
Curvature-based
603 2994
Fran@ Derivatives
function
algorithm
99 27s 137
functions
CPU
Holodniok’s
Brunt
NUmCriC
No.
No.
of
Van
Parameter,
2
(1987). The new algorithm required less continuation steps, less Newton corrections and less CPU time on several small, but highly non-linear, test problems. The number of step-sizes rejected due to poor Newton convergence characteristics or large turning angle during the correction process was also
An improved homotopy continuation algorithm has been developed by incorporation of differential geometry into the prediction and step-length phases of the algorithm of Frantz and Van Brunt
3. Solution
1
0
lmbedded
t
CONCIAJSIONS
Table
-1
Curvnture-based
algorithm
Numeric 75 214 74 406 2418 I573 Initial
Analytic 75 214 74 406 748 I327
csrimate
(b)
W. L. RION and V. VAN BRUNT
900
reduced. The development of this algorithm is significant since, by incorporation of the intrinsic geometry of the homotopy path, it allows for finding multiple solutions to engineering problems by homotopy continuation with less computation than previous algorithms. NOMENCLATURE B C D f F h H j
= Unitbinormal vector of Frenet frame, B = T x N = Arbitrary curve = Dimensionless adiabatic temperature rise = Function = Function vector = Homotopy function = Homotopy function vector = jik = Element of Jacobian matrix inverse (except when k
corresponds to continuation variable or when i = n + 1) i,., = Jff,/&,, i.e. change in homotopy function i with respect to .%, only when k = continuation variable
j. + ,.( = dx,/d.s, elements of tangent vector corresponding to variable j j+‘-” = Elements of augmented Jacobian inverse [J] = Jacobian marrix [J ‘I= Augmented Jacobian matrix n = Number of problem variables, dimension of F N = Unit normal vector of Frenet frame, N = C”(s)/~/C”(s)~~ R = Recycle ratio s = Arclength of solution trajectory I = Imbedded parameter T = Unit tangent vector to homotopy curve T-, = Reduced tangent vector defined by [d_r,/ds, dx,/d.s, . ,
dx~-,lhd++,
U = Correction
ds, .
, dx.ldr, dx,, ,/dsl
vector generated by Newton’s method x = Independent variable X = Independent variable vector y, = Reactant conversion in reactor i = I or 2 z,= Summation of the form: I$: ;g: 2 = Vector
x52 I
formed by components
2 !, ii
Software
Kubicek
2, 98-107
M.
Bifurcation
and
M.
Theory
(1976).
Marek, Computational Methods in and Dissiparive Structures. Springer-
Verlag, New York (1983). Kubicek M., H. Hofmann, V. Hlavacek and J. Sink& Multiplicity and stability in a sequence of two nonadiabatic non-isothermal CSTR. Chem. Engng Sci. 35, 987-996 (1980). Kubicek M., M. Holodniok and I. Marek, Numerical solutions of non-linear equations by one-parameter imbedding methods. Numer. Funct. Analysis Optimizarion 3(2),
223-264
(1981).
Kuno M. and J. D. Seader, Computing all real solutions to systems of non-linear equations with a global fixed point homotopy. Ind. Engng Chem. Res. 27, 1320-1329 (1988).
Greek B = p, = As = 4 = &= 0 = r #, y I$,
= = = =
Step-length multiplier (no subscript) Dimensionless heat transfer coefficient for reactor i Step-length used during continuation step Column vector of intermediate values used in matrix inversion by bordering technique Row vector of intermediate values used in matrix inversion by bordering technique Turning angle between tangent vectors to homotopy path Torsion Dimensionless temperature in reactor i = I or 2 Dimensionless activation energy Dimensionless coolant temperature
Subscripts ideal = Ideal value of a parameter i,j,k = Arbitrary parameter signification n = n th parameter Superscripts -
REFERENCES
Allgower E. and K. Georg, Simplicial and continuation methods for approximating fixed points and solutions to systems of equations. SIAM Rev. 22( 1), 28-85 (1980). Den Heijer C. and W. C. Rheinboldt, On step-length algorithms for a class of continuation methods. SIAM J. Numer. Analysis W(5), 925-948 (1981). Deuflhard P., A stepsize control for continuation methods and its special application to muitiple shooting techniques. Numer. Mafh. 33, 115-146 (1979). Faddeev D. K. and V. N. Fadeeva, Computational Methods of Linear Algebra, pp. 163-167. Freeman, San Francisco (1963). Frantz R. W.. Development and application of Continuation methods for the solution of engineering systems. M.S. Thesis, University of South Carolina (1986). Frantz R. W. and V. Van Brunt, A differential homotopy continuation method for interlinked solvent extraction cascades. Separal. Sci. Technof. 22(2/3), 243-267 (1987). Georg K., A note on stepsize control for numerical curve following. Proc. NATO Advanced Research Institute on HomoJopy MeJhods and Global Conrwrgence, pp. 145-154. Plenum Press, New York (1980). Himmelblau D. M., Applied Nonlinear Programming, pp. 427428. McGraw-Hill, New York (1972). Kubicek M., Algorithm 502, dependence of solution of Trans. Math. nonlinear systems on a parameter. ACM
I = matrix inversion ’= Differentiation with
*= I’= + = k = T =
respect to s Two differentiations with respect to s Three differentiations with respect to s Member of Augmented Jacobian or its inverse Iteration number Transpose
[A- B] = Inner product or “dot product” of vectors A and B [A x B] = Cross product of vectors A and B /IAl] = Euclidean norm of vector A
Morgan for
A., Solving Polynomial Systems Using Continuatron Engineering and ScienriJic Problems. Prentice-Hall,
New Jersey (1987). Reklaitis G. V., A. Ravindran and K. M. Ragsdell, Engineering Optimization. Wiley, New York (1983). Rheinboldt W. C., Solution fields of nonlinear equations and continuation methods. SIAM J. Numer. Analysis (1980). 17(2), 221-237 Rheinboldt W. C. and J. V. Burkardt, A locally parameterized continuation process. ACM Trans. Math. Software 9(2),
2 I S-253
(1983a).
Rheinboldt program
W. C. and J. V. Burkardt, Algorithm 596. A for locally parameterized continuation process. ACM Trans. Math. Software 9(2), 236241 (1983b). Seydel R. and V. Hlavacek. Role of continuation in engineering analysis. Chem. Engng Sci. 42(6), 1281-1295 (1987). Spivak M., A Comprehensiw Vol. 4. Publish Geomerrv,
Itttroduction
or Perish
to DJrerential
Inc.,
Wilmington Delaware (1979). Stoker J. J., Differential Geometry. Wiley Interscience, New York f 19691. Watson .L. T.1 C. D. Billups and A. P. Morgan, HOMPACK: a suite ofcodes for glob&y convergent homotopy algorithms. ACM Trans. Math. Software 13(3), 2X1-310 (1987). Wayburn T. L., A review of continuation methods and their applications to separations problems. CACHE Newlett. Spring, S-22 (1988). Wayburn T. L. and J. D. Seader, Homotopy continuation methods for computer-aided process design. Comput. them. Engng 11, 7-25 (1987).