Differential geometry based homotopy continuation

Differential geometry based homotopy continuation

Computers them. Printed in Great Engng, Britain. Vol. 14, No. 8, pp. 889-900, 1990 All rights reserved DIFFERENTIAL Copyright GEOMETRY BASED CON...

1MB Sizes 0 Downloads 75 Views

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).