Global Optimization for the Chemical and Phase Equilibrium Problem

using Interval Analysis

K.I.M. McKinnon, C.G. Millar^{†}, M. Mongeau^{∗}

MS 95-01

January 1995

†: Supported by EPSRC research grant RR31817/18 (ECOSSE project)

∗: Laboratoire Approximation et Optimization

Université Paul Sabatier, 118 route de Narbonne - 31062 Toulouse cedex - France Report also appears as LAO 95-02 and ECOSSE 1995-01

(Submitted to the International Conference on:

State of the Art in Global Optimization: Computational Methods and Applications, Princeton University, April 28-30, 1995)

Department of Mathematics and Statistics

University of Edinburgh, The King’s Buildings, Edinburgh EH9 3JZ Tel. (33) 131 650 5042 E-Mail : ken@maths.ed.ac.uk

Abstract

This paper addresses the problem of minimizing the Gibbs free energy in the n_{c}- component, multi-phase chemical and phase equilibrium problem involving different thermodynamic models. The algorithmic approach used is based on the tangent-plane criterion of Gibbs: the global optimization problem considered, which involves a search space of n_{c}(n_{c} − 1) dimensions, is reduced to a finite sequence of local optimization steps in K(n_{c} − 1)-space, K ≤ n_{c}, and global optimization steps in n_{c} − 1-space.

We describe an algorithm performing the global optimization step involved in this lower-dimensional search space using techniques from interval analysis. We report promising numerical results on instances of the Gibbs free energy minimization problem.

Key words:

Global optimization, interval analysis, tangent-plane criterion, Gibbs free energy, chemical and phase equilibrium, non-convex optimization.

1 Introduction

The problem of phase equilibrium is: Given the total quantities (in number of moles) of n_{c} components, Z_{i} > 0, i = 1, …, n_{c}, determine the number, K, and type (e.g. the state; liquid, vapour) of phases existing at equilibrium and the composition (number of moles of each component), {z_{ik}}_{i=1}^{n_{c}}, of each phase k, 1 ≤ k ≤ K. In the problem of chemical phase equilibrium (i.e. involving reaction), the n_{c} components are constituted of m chemical elements, and the {Z_{i}}_{i=1}^{m} are the total quantities of each of the chemical elements [12].

These equilibrium quantities can be obtained by minimizing globally a thermodynamic function, the Gibbs free energy (GFE). This is a problem which is known to be numerically difficult to solve even in the case where the number of phases K present at equilibrium is known, and even for instances involving as few as n_{c} = 2 components and K = 2 phases [13].

We consider n_{c}-component chemical and phase equilibrium problems which involve different thermodynamic functions to model the different types of phase that may be present at equilibrium. Also, this will require us to make use of an extension, presented in [8], of the tangent-plane criterion [1, 2, 3, 5]. The algorithmic approach we use is also taken from [8]: the problem is reduced to a finite sequence of local optimization steps, and of global optimization steps in a lower-dimensional space. These results are summarized in section 2.

Related work has been done by Michelsen [9, 10] who uses alternately a local method to obtain a stationary point of the GFE corresponding to a given number of phases, and a stability test derived from the tangent-plane-criterion. He however uses a local method to attempt to determine whether a phase is stable, making use of several initial estimates. Also, Michelsen concentrates on phase equilibrium problems in which all phases are described with the same thermodynamic model. In our approach, different thermodynamic models for the different types of phase present at equilibrium are allowed. Moreover, the way we make use of the tangent-plane criterion is different. We use the tangent-plane criterion either to add an extra phase, to improve the Gibbs free energy without changing the number of phases, or to establish global optimality.

The unspecified global optimization subroutine of the generic algorithm, GILO (Global Improvement - Local Optimization), in [8] for the global minimization step in the lowerdimensional space will here be implemented using techniques from interval analysis [4, 11]. We survey, in section 3, the techniques we use for the global step. Section 4 reports numerical results obtained from the implementation of our algorithm on 2-component and 3-component phase equilibrium problems, and we compare our results with those obtained by McDonald & Floudas in [6]. We present some conclusions in Section 5.

We consider the following global optimization problem, which we shall refer to as problem P:

inf F (K, t, Y) ≡ ∑_{k=1}^{K} y_{k}f_{t_{k}}(x_{k})

subject to ∑_{k=1}^{K} y_{k}Ax_{k} = b  (1)

∑_{k=1}^{K} y_{k} = 1
y_{k} ≥ 0, k = 1, 2, …, K
x_{k} ∈ S, k = 1, 2, …, K
t_{k} ∈ τ, k = 1, 2, …, K
K ∈ IN^{+},

where Y = (y_{1}, …, y_{K}, x_{1}^{T} , …, x_{K}^{T}) ∈ IR^{K(n+1)}, A is an m × n real matrix of rank m ≤ n, 0 ≠ b ∈ IR^{m}, S is an open convex subset of (IR^{+})^{n}, τ is a finite index set, and, for k = 1, 2, …, K, t_{k} ∈ τ, y_{k} ∈ IR, x_{k} ∈ IR^{n}, and f_{t_{k}} is a differentiable real-valued function defined over S satisfying

lim_{x→x̄} d^{T} ∇f(x) = −∞,  (2)

for any x̄ in the boundary of X, where d is a direction from x̄ into the interior of S. Note that in the above optimization problem we are interested in a solution for which K is smallest, amongst all the possible global minima (to avoid obtaining a solution with several phases having identical compositions).

In what follows N is equal to n_{c} in the case of the phase equilibrium problem, and N is equal to m + 1 in the case of the chemical phase equilibrium problem. Let us denote the molar fraction of “component” (“chemical element”, in the case of the chemical phase equilibrium problem) i in phase k by

x_{ik} ≡ \frac{z_{ik}}{∑_{i=1}^{n_{c}} z_{ik}} i = 1, 2, …, n_{c}; k = 1, 2, …, K.  (3)

Let

Z = ∑_{i=1}^{N} Z_{i},  (4)

n = N − 1, and

S = {x_{k} ∈ IR^{n} : ∑_{i=1}^{n} x_{ik} < 1, x_{ik} > 0},

(where x_{ik} denotes the ith component of vector x_{k}). The n_{c}-component phase equilibrium problem (with m = n, A = I) and the n_{c}-component, m-element chemical equilibrium problem, are special cases of problem P [with A being a matrix related to the formula matrix whose entries give the number of atoms of each chemical element per molecule of each component, and b having Z_{i}/Z as ith entry, i = 1, …, N (see [5] and [8])].

We assume that the problem has at least one feasible solution. For the chemical and phase equilibrium problem, we have that (2) holds. This implies we can assume [5] that any local infimum of P is a local minimum.

2 GILO: A Generic Global Minimization Algorithm for the Chemical and Phase Equilibrium Problem

In this section, we summarize the results of [8] which produced the GILO (Global Improvement - Local Optimization) algorithm.

The number of phases, K, that may exists at equilibrium is unknown a priori. Let P K̃ ; t̃ denote problem P in which K is fixed to a positive integer K̃ and t is fixed to t̃∈ τ K̃, and let F K̃ ; t̃ : IR^{K̃(n+1)} → IR be the corresponding objective function.

Lemma 1 If Y^{∗} is a Kuhn-Tucker point for P_{t^{∗}}^{K^{∗}}, for some positive integer K^{∗} and some t^{∗} ∈ τ^{K^{∗}}, then there exists an m-dimensional hyperplane Θ^{∗} which is tangent to f_{t_{k}^{∗}} at the points {x_{k}^{∗} : 1 ≤ k ≤ K^{∗} and y_{k}^{∗} > 0}.

Let f_{M} : S → IR be the (nonsmooth) function defined by

f_{M}(x) ≡ min_{t∈τ} f_{t}(x).

Lemma 2 Let (K^{+}, t^{+}, Y^{+}), be a global minimum for P, K^{+} ∈ IN^{+}, t^{+} ∈ τ^{K^{+}}, and Y^{+} ∈ IR^{K^{+}(n+1)}, and Θ^{+} be the corresponding tangent hyperplane described by Lemma 1. Then we have that, at the points {x_{k}^{+} : 1 ≤ k ≤ K^{+} and y_{k}^{+} > 0}, f_{M} is differentiable and Θ^{+} is tangent to f_{M}.

We have the following characterization of global minima for P in terms of the function f_{M}.

Theorem 1 (Tangent-Plane Criterion) Let K^{+} ∈ IN^{+}, t^{+} ∈ τ^{K^{+}}, and Y^{+} ∈ IR^{K^{+}(n+1)}.

If (K^{+}, t^{+}, Y^{+}) is a global minimum for P, then Y^{+} is a Kuhn-Tucker point for P_{t^{+}}^{K^{+}} and

f_{M}(x) ≥ Θ_{M}^{+}(x), for all x ∈ S,  (5)

where Θ_{M}^{+} is the tangent hyperplane corresponding to (K^{+}, t^{+}, Y^{+}), as described by Lemma 2.

Conversely, if Y^{+} is a Kuhn-Tucker point for P_{t^{+}}^{K^{+}} and

f_{M}(x) ≥ Θ^{+}(x), for all x ∈ S,  (6)

where Θ^{+} is the tangent hyperplane corresponding to (K^{+}, t^{+}, Y^{+}) as described by Lemma 1, then (K^{+}, t^{+}, Y^{+}) is a global minimum for P and Θ^{+} = Θ_{M}^{+}.

We know from the following proposition that we do not need to go beyond K = m + 1 phases in order to solve P.

Proposition 1 If Y ∈ IR^{K(n+1)} is a Kuhn-Tucker point for P_{t}^{K}, for some fixed positive integer K and some t ∈ τ^{K}, then there exists a solution Y^{m+1} ∈ IR^{(m+1)(n+1)} for P_{t^{m+1}}^{m+1}, for some t^{m+1} ∈ τ^{m+1}, such that F_{t^{m+1}}^{m+1}(Y^{m+1}) = F_{t}^{K}(Y).

Let LO be a local optimization subroutine which, given the current feasible solution, (K, t, Y), for P, returns a local minimum, Y^{∗} := LO(K, t, Y) ∈ IR^{K(n+1)}, for P_{t}^{K}. Let GI be a global optimization subroutine (for a lower-dimensional search space) which, given a linear function, Θ^{∗} : IR^{n} → IR, returns (t_{0}, x_{0}) := GI(Θ^{∗}), where x_{0} is a global minimum for the subproblem

min_{x∈S} f_{M}(x) − Θ^{∗}(x),  (7)

and where t_{0} ∈ τ is such that f_{t_{0}}(x_{0}) = f_{M}(x_{0}).

GILO Algorithm

Step 0: [Initialization] Let K^{∗} = 1, and let t^{∗} be an arbitrary element of τ. Let Y^{∗} be a feasible solution to P_{t^{∗}}^{K^{∗}}. Step 1: [Local Optimization] Y^{∗} := LO(K, t, Y). t_{k}^{∗} := t_{k}, 1 ≤ k ≤ K. K^{∗} := number of non-zero distinct phases (with the rest renumbered if necessary). Step 2: [Global Improvement Subproblem] Let 1 ≤ j ≤ K^{∗} be arbitrary. Let Θ^{∗} : IR^{n} → IR be defined by Θ^{∗}(x) := ∇f_{t_{j}^{∗}}(x^{∗})^{T} (x − x_{j}^{∗}) + f_{t_{j}^{∗}}(x_{j}^{∗}), (t_{0}, x_{0}) := GI(Θ^{∗}). Step 3: [Global Optimality Condition] If f_{M}(x_{0}) ≥ Θ^{∗}(x_{0}), then stop: (K^{∗}, t^{∗}, Y^{∗}) is a global minimum for P. Step 4: [Force improvement] Solve

∑_{k=1}^{K^{∗}} α_{k} = 1, ∑_{k=1}^{K^{∗}} α_{k}Ax_{k}^{∗} = Ax_{0}  (8)

for {α_{k}}_{k=1}^{K^{∗}}.

• (i) [Interchange a phase]

If (8) has a solution, {α_{k}}_{k=1}^{K^{∗}}:

The new feasible solution to P is (K, t, Y) given by

t_{i} := t_{0},
y_{i} := \frac{y_{i}^{∗}}{α_{i}} , x_{i} := x_{0},
t_{k} := t_{k}^{∗}, y_{k} := y_{k}^{∗} − \frac{y_{i}^{∗}}{α_{i}} α_{k} , x_{k} := x_{k}^{∗}, 1 ≤ k ≤ K^{∗}, k ≠ i,

where

i ≡ argmin\{\frac{y_{k}^{∗}}{α_{k}} : α_{k} > 0, 1 ≤ k ≤ K^{∗}\}.

K := number of phases with non-zero y_{k}^{∗} (with the rest renumbered if necessary).

• (ii) [Split a phase] Otherwise (the system (8) has no solution): Let ϵ_{0} > 0 be such that the new solution (K, t, Y) given below is feasible for P.

K := K^{∗} + 1 , t_{K^{∗}+1} := t_{0},
y_{K^{∗}+1} := \frac{ϵ_{0}}{ϵ_{0} + 1}y_{j}^{∗} , x_{K^{∗}+1} := x_{0},
y_{j} := \frac{1}{ϵ_{0} + 1}y_{j}^{∗} , x_{j} := ϵ_{0}(x_{j}^{∗} − x_{0}) + x_{j}^{∗},
t_{k} := t_{k}^{∗} , 1 ≤ k ≤ K^{∗},
y_{k} := y_{k}^{∗} , x_{k} := x_{k}^{∗}, 1 ≤ k ≤ K^{∗}, k ≠ j.

Go to Step 1.

Note that if m ≥ n, then finding a single-phase starting solution Y^{∗} (Step 0), reduces to solving the equations to find the unique feasible solution (in the case with no reaction [m = n and A = I], Y^{∗} is straightforwardly given by y_{1}^{∗} = 1 and x_{1}^{∗} = b). This starting solution, Y^{∗}, will then be a local minimum of problem P_{t^{∗}}^{K^{∗}}, and Step 1 does not need to be performed in the first iteration.

The phase-stability test (Steps 2 and 3) uses the tangent-plane criterion to determine whether the current solution is optimal, and, if it is not, it finds an improved feasible solution either with the same number or fewer phases [Step 4 (i)], or with one added phase [Step 4 (ii)], and it also determines what type of phase is to be added. The number of phases may indeed be reduced at the end of Step 4 (i) (or at the end of Step 1) by eliminating the phases corresponding to null y_{k}’s (and renumbering the remaining phases). Also, at the end of Step 1, spurious phases (i.e the proportions of the components of two phases of same type are identical), are merged into one phase, and the number of phases is hence decreased. Step 4 (ii) corresponds physically to splitting a phase (phase j) into two phases (the new phase j and phase K^{∗} + 1). Hence, after performing this step, K is always increased by one. In Step 4 (i) we interchange the phase x_{i} for the new phase x_{0}. In both cases, a strict decrease of F is obtained [8] and we leave the current local minimum.

Theorem 2 (Convergence) The GILO Algorithm converges to a global minimum of P, from any point (K^{0}, t^{0}, Y^{0}) which is feasible for P and such that K^{0} = 1, in a finite number of steps, provided that P_{t}^{K} has a finite number of distinct local minimum objective function values, for any 1 ≤ K ≤ m + 1 and any t ∈ τ^{K}.

This convergence theorem remains valid in the case where the subroutine GI (Global Improvement) simply returns any point x_{0} such that

f_{t}(x_{0}) < Θ^{∗}(x_{0}),

for some t ∈ τ (or such that f_{M}(x_{0}) < Θ^{∗}(x_{0})). In other words, in order to obtain from Step 4 a strict decrease in F, one does not necessarily need the point x_{0} that most violates the tangent-plane criterion.

3 Global Optimization Using Interval Analysis

In this section we describe a deterministic method for global optimization, guaranteed to find an ϵ-global solution (i.e. whose objective function value is within ϵ > 0 of a globally optimal objective function value), which is based on interval analysis. With this approach, one recursively partitions the search space until each subdomain either can be shown not to contain a global solution, or contain only ϵ-global solutions. To achieve this aim, one makes use of bounds on the function, gradient and/or Hessian values over the subdomains within a branch-and-bound (B&B) scheme.

The interval-analysis approach is an attractive method to address the global optimization subproblem (7) for many reasons. Firstly, as we mentioned above, Step 2 of the GILO Algorithm does not need to perform a thorough global search in order to escape the current local optimum: we can use interval-analysis techniques only until we obtain any point x_{0} violating the tangent-plane criterion. We do not need a global minimum for the subproblem (7), but only any point such that the objective function value of this problem is negative. Since the objective function takes the value zero at each of the current tangent points, the B&B search can be started with an upper bound of zero on the optimal solution. Any subdomain over which this objective function is shown to be non-negative will then be rejected at once. Secondly, since the search space for the global optimization subproblem is low dimensional (n-space versus (m + 1)n-space for the original problem), the subdividing step for its recursive partition should not be prohibitive. Thirdly, most of the work in the B&B is the same from one global optimization subproblem to another (as the corresponding objective functions differ only, through the tangent planes, by an affine amount): in subsequent searches the gradients differ by a constant amount, and the Hessians are the same. This feature could be exploited by storing these quantities from one B&B search to another.

Let I be a bounded interval. The lower and upper bounds on I are denoted by Ī and Ī, so I = [Ī, Ī]. We call a vector of intervals a box. Let f be a continuous real-valued function defined over IR, and X ⊂ IR^{n} be a box. We define the range of f for x ∈ X, denoted by f^{R}(X), to be the interval

f^{R}(X) ≡ [ inf_{x∈X} f(x), sup_{x∈X} f(x)]

(where this is understood as a semi-infinite interval or as an infinite interval, according to whether one or both bounds are ±∞). When using the interval arithmetic [4, 11] (e.g. if X and Y are bounded intervals, X − Y = [X̄ − Ȳ , X̄ − Ȳ]), we always obtain a valid bound, B, for a function f in a box X, i.e such that f^{R}(X) ⊆ B. However, in general we do not necessarily get sharp bounds (we do not necessarily have f^{R}(X) = B). We call an interval-valued function F having as argument a box of IR^{n}, an (interval-)extension of f if F ([x_{1}, x_{1}], …, [x_{n}, x_{n}]) = f(x_{1}, …, x_{n}). The extension is not unique: for example both

F_{1} = X^{2} − 2X and F_{2}(X) = X(X − 2)

are extensions of f(x) = x^{2} − 2x. Also, depending on which interval extension of a function one chooses, different bounds can be obtained for f over X. In this example, with X = [0, 3], we have

F_{1}([0, 3]) = [0, 9] − [0, 6] = [−6, 9]

(i.e. x ∈ [0, 3] implies f(x) ∈ [−6, 9]), whereas

F_{2}([0, 3]) = [0, 3].[−2, 1] = [−6, 3],

which yields a better bound for f over [0, 3]. Note that f^{R}([0, 3]) = [−1, 3].

Different techniques can be used in order to obtain as tight bounds as possible, the aim being to find where is achieved a global minimum, x^{+}, of f over some specified domain:

1. A reformulation of the function (see the example above) can yield better bounds (another example of a reformulation is to use a Taylor expansion, for instance around the mid point of the box).

2. In the case where f is differentiable, one can use bounds on the gradient of f in order to eliminate subboxes. If we denote by F′(X^{k}) a bound obtained on the gradient of f over some subbox X^{k}, and that we have 0 ∉ F′(X^{k}), then we can eliminate this subbox (recall that from (2), a global minimum of the GFE cannot occur on the boundary of the feasible domain). We shall refer to this technique as a gradient test.

3. In the twice-differentiable case, a box can be eliminated if it can be shown that the Hessian is not positive semi-definite anywhere in the box (Hessian test). We use the computationally cheap condition for this of testing the diagonal Hessian elements. If any of these are negative through the box, the Hessian is nowhere positive semi-definite in the box so the box can be eliminated.

4. The structure of the objective function can often be exploited. Consider for example the univariate case. We have

\overline{f^{R}}(X) = max\{f(X̄), f(X̄)\},

if f is convex over X, and

\underline{f^{R}}(X) = min\{f(X̄), f(X̄)\},

if f is concave over X (for a general box, these bounds are attained at an extreme point of X).

5. The interval Newton step. Given an initial box subdomain containing a global optimum, the interval version of Newton’s method produces a smaller subdomain which also contains a global optimum, reducing all dimensions ultimately at a quadratic rate. This involves the solution of a system of equations with interval coefficients. The (real-valued) Hessian at the mid point of the box is found. This is used as a preconditioner for the interval Newton equations, and one sweep of interval Gauss-Seidel is then applied in an attempt to obtain tighter bounds [4].

6. Interval subdivision involves a recursive partition of the domain into subboxes over each of which bounds are computed. This being too costly in high-dimensional spaces, we integrate this into a B&B scheme. Any evaluation of f at a feasible point gives a lower bound for f over the original domain. In the next section, we shall refer to mid-point evaluations when we systematically evaluate f at the center point of each subbox obtained by subdivision. The smallest lower (respectively upper) bound obtained over all subboxes resulting from the subdivision process yields a lower (respectively upper) bounds for f(x^{+}). The aim of the overall B&B process is to eliminate subboxes which can be shown not to contain a global minimum of f. In addition to the use of bounds, one can also apply techniques 1, 2, 3, 4, and 5 above to eliminate a given subbox.

The reader is referred to [4] and [11] for a more complete description of global optimization using interval analysis.

We shall stop this global optimization process whenever either each remaining box has width less than ϵ_{w}, or a lower bound for the objective function is within ϵ_{f} of the upper bound (over the remaining boxes) for the objective function, where ϵ_{w} and ϵ_{f} are sufficiently small positive constants. Moreover, in the current GFE minimization application, we also stop whenever the mid point of a box is found that lies at a distance of ϵ_{θ} or more below the tangent plane [i.e whenever we find a point yielding a sufficiently negative objective function value for problem (7)], for some small positive constant ϵ_{θ}.

4 Computational Experiments

In this section we report numerical results obtained on two instances of the GFE global minimization problem without reaction (i.e. here A = I). The objective function we are using models the GFE for ideal vapour phases, and for liquid phases which can be described by the NRTL equation [6]. More precisely, the objective function has the following form:

G(K, t, x, y) = Z ∑_{k=1}^{K} y_{k}G_{t_{k}}^{k}(x),  (9)

where 1 ≤ K ≤ N, t_{k} ∈ τ, k = 1, 2, …, K, τ is either {L} (Problem 1) or {L, V } (Problem 2), Z is defined by (4), y_{k} is the molar fraction of phase k, given by

y_{k} ≡ \frac{∑_{i=1}^{n_{c}} z_{ik}}{Z} , k = 1, 2, …, K,

and x_{ik} is defined by (3). In the case of an instance involving only liquid phases (τ = {L}), the function G_{L}^{k} is defined by

G_{L}^{k} = ∑_{i=1}^{n_{c}} x_{ik}(∑_{j=1}^{n_{c}} \frac{x_{jk}m_{ij}}{∑_{l=1}^{n_{c}} x_{lk}G_{lj}} + ln x_{ik}),  (10)

where m_{ij} and G_{ij} are given data. When both liquid and vapour phases are allowed (τ = {L, V }), the functions G_{L}^{k} and G_{V}^{k} are defined as follows:

G_{L}^{k} = ∑_{i=1}^{n_{c}} x_{ik}(∑_{j=1}^{n_{c}} \frac{x_{jk}m_{ij}}{∑_{l=1}^{n_{c}} x_{lk}G_{lj}} + ln x_{ik} + ln P sat_{i}),  (11)

and

G_{V}^{k} = ∑_{i=1}^{n_{c}} x_{ik}(ln x_{ik} + ln P_{i}),  (12)

where m_{ij}, G_{ij}, P sat_{i}, and P_{i} are given data.

We shall describe below the behaviour of the our algorithm on two instances with data taken from [6]. The following equations are required to obtain the data involved in the formulation of both problems:

G_{ij} = exp(−α_{ij}t_{ij}),  (13)
m_{ij} = t_{ij}G_{ij},  (14)
α_{ij} = α_{ji},  (15)
α_{ii} = 0,  (16)
t_{ii} = 0,  (17)

1 ≤ i, j ≤ n_{c}.

The local optimization subroutine LO (referred to in the GILO algorithm) we used is the NAG E04UCF subroutine (sequential quadratic programming).

The global optimization subroutine GI, of the GILO algorithm will here be the implementation of the interval-analysis techniques described in section 3. This yields an algorithm which we call Interval-GILO. The specific interval-analysis techniques implemented, within the interval-subdivision B&B process, are the following. Firstly, we reformulated the objective function of the global optimization subproblem (7) (the tangent-plane distance function) using a first-order Taylor expansion around the mid point of each subbox, of width smaller than ϵ_{T} ≡ 0.05, obtained during the course B&B process. We did the same for the gradient of this objective function. We also obtained valid bounds on the objective function and on its gradient without using the Taylor expansion for subboxes of width larger than ϵT̄ ≡ 0.0005 (when both interval extensions—with and without Taylor expansion—are used, the best bound is used). Secondly, gradient tests were used in order to eliminate some subboxes. Use of the Hessian test did not improve the results. Thirdly, we used the interval version of the Newton step each time the algorithm encountered a subbox of width less than ϵ_{N} ≡ 0.0005. Mid-point evaluations were used in the first N_{m} ≡ 100 iterations. We used ϵ_{w} = ϵ_{f} = ϵ_{θ} = 10^{−5} and ϵ_{0} = 10^{−6} in all our tests. Breadth-first search is used in the B&B implementation.

The global search takes place in the reduced space of the variables x_{1k}, …, x_{n_{c}−1,k}. In order words, we deal with the equality constraints as follows. Bounds on the remaining variables x_{n_{c}k} are obtained at each GI iteration from the intersection of its defined bounds, [0, 1], with bounds obtained from the relations

x_{n_{c}k} ≡ 1 − ∑_{i=1}^{n_{c}−1} x_{ik}, k = 1, …, K,

and the bounds on the reduced-space variables. If the resulting interval is empty, then the current box has no feasible points and can be deleted.

As interval methods are easily programmed in Fortran 90, our techniques were implemented with this new standard for Fortran. All CPU time figures reported in the next sections were obtained from runs on a Sun Sparc5 70MHz workstation. The number of local-optimization iterations we shall report refers to the number of iterations needed for the NAG subroutine to converge. The number of global-optimization iterations will correspond to the total number of times an interval has been divided in two in the B&B interval-analysis process.

4.1 Problem 1: LL Equilibrium

This problem corresponds to the Illustrative Example (Liquid-Liquid Equilibrium n-Butyl-acetate - Water) of [6].

For this problem, n_{c} = 2 and τ = {L}. In this case the simple structure of the objective function can be exploited so as to produce a reformulation which gives tighter bounds for both the objective function and for its gradient in the interval-analysis global optimization step. More specifically, we reformulated the G_{L}^{k} part of the objective function (9) as follows. We first replace x_{2k} by 1 − x_{1k} (recall that by definition of molar fractions, we have x_{1k} + x_{2k} = 1), and we let u ≡ x_{1k} for the sake of notational simplicity. One can then show that, for k = 1, 2, the function G_{L}^{k} can be written in the form:

G_{L}(u) = h_{1}(u) + h_{2}(u) + h_{3}(u),

where

h_{1}(u) = \frac{m_{12}}{\frac{G_{12}}{1−u} + \frac{1}{u}} ,
h_{2}(u) = \frac{m_{21}}{\frac{G_{21}}{u} + \frac{1}{1−u}} ,
h_{3}(u) = u ln u + (1 − u) ln(1 − u).

The denominators of h_{1} and h_{2}, and also the function h_{3} being convex in u, we can obtain sharp bounds for each of the functions h_{1}, h_{2}, and h_{3} by using the corresponding interval extension of each of these functions. The fact that the bound for G_{L}^{k} we compute is the sum of three sharp bounds will yield better bounds for the objective function. Similarly for the derivative of G_{L}(u), one can show that we have:

G′_{L}(u) = l_{1}(u) + l_{2}(u) + l_{3}(u),

i j  t_{ij}  G_{ij}  α_{ij}  Z_{i}

(mol)

1 2 3.00498 0.30794 0.39196 0.5
2 1 4.69071 0.15904 0.39196 0.5

Table 1: Data for Problem 1

t_{1}^{+} = L  t_{2}^{+} = L
i z_{i1}^{+} z_{i2}^{+} G^{+}
1 7.13609 ×10^{−4} 0.499286 -2.01983 ×10^{−2}
2  0.155879  0.344121

Table 2: Solution for Problem 1

where

l_{1}(u) = m_{12} \frac{(1/u − 1)^{2} − G_{12}}{[(1/u − 1) + G_{12}]^{2}} ,
l_{2}(u) = m_{21} \frac{G_{21}(1/u − 1)^{2} − 1}{[1 + G_{21}(1/u − 1)]^{2}} ,
l_{3}(u) = ln(\frac{1}{1/u − 1}).

Again sharp bounds can be obtained for each of the three functions l_{1}, l_{2}, and l_{3}, and this will yield better bounds for the gradient of the objective function. Using these better valid bounds for the objective function and its gradient, we derived the Taylor expansion reformulation of the objective function and of its gradient. We used a second-order Taylor expansion of the objective function, and a first-order Taylor expansion of its gradient.

The necessary data for Problem 1 are obtained from Table 1 together with equations (13) to (17).

The program was started with the only feasible point corresponding to a single phase. The algorithm performed successively a global optimization step, a local optimization step, and a second global optimization step. The first global optimization step yielded the introduction of a second liquid phase, while the second global step simply proved that the iterate found by the local optimization step was (ϵ-)globally optimal. The global minimum, (K^{+}, t^{+}, z^{+}), found by the algorithm corresponds to two liquid phases, and it is given by Table 2, where K^{+} = 2 and the z_{ik}’s are in moles. Table 3 reports numbers of iterations and CPU times required to obtain convergence. McDonald & Floudas [6] report 1.23 seconds of CPU time to converge to the same solution, with ϵ_{f} = 0.0005. The final global step of the Interval-GILO ended with both stopping criteria ϵ_{f} and ϵ_{w} satisfied (i.e. the last iterate, say z̄, is within distance ϵ_{w} of the global minimum, z^{+}, and the objective function evaluated at z̄ is within ϵ_{f} of the global minimum function value), with ϵ_{f} = ϵ_{w} = 10^{−5}. When we do not use the reformulation of G_{L}^{k} described above

# of iterations CPU time (sec.)
1st global step  8  0.03
1st local step  9  0.07
2nd global step  123  0.17
entire program  0.42

Table 3: Results for Problem 1

# of iterations CPU time (sec.)
1st global step  8  0.07
1st local step  9  0.09
2nd global step  190  1.47
entire program  1.83

Table 4: Results for Problem 1: Original Formulation of G_{L}^{k}

and use a first-order Taylor expansion for the objective function, a substantially larger number of global-step iterations is required to obtain convergence, and the total CPU time is much increased (see Table 4).

4.2 Problem 2: LLV Equilibrium

This problem corresponds to Example 6A (Liquid-Liquid-Vapour Equilibrium for Benzene - Acetonitrile - Water) of [6]. For this problem, n_{c} = 3 and τ = {L, V }. The necessary data are obtained from Table 5 together with equations (13) to (17) and:

t_{ij} = \frac{g_{ij} − g_{jj}}{RT} ,

R = 1.9872 cal/K/mol, T = 333 K.

If there is a vapour phase present in the current local minimum solution then the global step will not be looking for a vapour phase. This is because we assume the vapour

i j g_{ij} − g_{jj} g_{ji} − g_{ii}  α_{ij}  i  P sat_{i}  P_{i}  Z_{i}
(atm)  (atm) (mol)
1 2 998.2  65.74 0.88577 1 0.512337709 0.769 0.34483
1 3 3883.2 3849.57 0.24698 2 0.485404650 0.769 0.31034
2 3 363.57 1262.4 0.3565 3 0.195255785 0.769 0.34843

Table 5: Data for Problem 2

t_{1}^{+} = V  t_{2}^{+} = L  t_{3}^{+} = L
i  z_{i1}^{+}  z_{i2}^{+}  z_{i3}^{+}  G^{+}
1  0.104643  7.31183 ×10^{−4}  0.239456
2 6.16355 ×10^{−2} 2.16922 ×10^{−2}  0.227012  -1.40852
3 5.24255 ×10^{−2}  0.262353  3.36516 ×10^{−2}

Table 6: Solution for Problem 2

# of iterations CPU time (sec.)
1st global step  2  0.08
1st local step  16  0.15
2nd global step  1  0.04
2nd local step  33  0.40
3rd global step  2504  52.24
entire program  53.03

Table 7: Results for Problem 2

phase to be ideal. That is to say, G_{V}^{k} is convex and hence one can show that no more than one phase of this type (vapour) needs to occur at any minimum of P.

The program was started with the only feasible point corresponding to a single vapour phase. The algorithm performed successively a global optimization step, a local optimization step, a second global step, a second local step, and a third global step. Again, each global step introduced a new liquid phase, except for the last global step which simply established the global optimality of the local minimum found by the second local step. The global minimum, (K^{+}, t^{+}, z^{+}), found by the algorithm corresponds to one vapour phase and two liquid phases. It is given by Table 6, where K^{+} = 3 and the z_{ik}’s are in moles. Again, the final global step of the Interval-GILO ended with both stopping criteria ϵ_{f} and ϵ_{w} satisfied, with ϵ_{f} = ϵ_{w} = 10^{−5}. Table 7 displays numbers of iterations and CPU times required to obtain convergence. These results, obtained on a Sun Sparc5 workstation, compare again favourably with those obtained by McDonald & Floudas’ GOP algorithm, for which they report, in [6], 766 seconds of CPU time on a Hewlett Packard 9000/730 machine to converge to the same solution, with relative convergence tolerance of 0.001. Note that in a recently submitted paper [7], McDonald & Floudas report a CPU time of only 0.76 seconds to verify the global optimality of the global minimum displayed on Table 6.

5 Conclusions

This paper describes the implementation of the finite generic algorithm (the GILO algorithm) for phase and chemical equilibria with multiple phase models presented by McKinnon & Mongeau in [8]. The contribution of the current paper is the implementation of interval-analysis techniques in order to solve the global optimization subproblem involved in the GILO algorithm: this yielded the Interval-GILO algorithm. We obtained promising numerical results on the two instances of the GFE minimization problem on which the Interval-GILO algorithm was tested.

We intend to study in greater detail the impact of the different possible combinations of the interval-analysis techniques we use here, to experiment on a wider variety of problems and thermodynamic models, and to investigate the effects of different function reformulations.

In this work, we made use of an unmodified NAG subroutine in order to solve the local optimization subproblem involved in the GILO algorithm of [8]. Future work should attempt to solve this numerically difficult local optimization problem more efficiently, for example along the lines proposed by Trangenstein in [13] for phase equilibrium problem involving two phases. Trangenstein adapts a modified Newton’s method so as taking advantage of the particular structure of the phase equilibrium problem.

Although it is known that every component is present in every phase, there are cases where the amounts of some components in some phases are so low that it is numerically simpler to treat them as being absent. This is the case for example in the modeling of petroleum reservoir fluid flow where one often wishes to consider one phase which contains only water and at least one phase which cannot contain water. Also, further work should address the issue of considering thermodynamic models which allow for some components to be absent from some of the phases.

The generic GILO algorithm can be used with any global optimization method and it would be interesting to test the global stability test of McDonald & Floudas’ GLOPEQ package with it. The Interval-GILO algorithm has the great advantage of generality: it does not rely on having a particular functional form for the thermodynamic model of the GFE function. All that is needed is an algorithm to evaluate the function values. Optionally, algorithms for evaluating higher-order derivatives can also be used to improve efficiency. The NRTL model is a severe test for interval methods as it contains an entropy term which is often almost exactly cancelled by a molecular interaction term. This is a situation where it is difficult to obtain good bounds from interval methods and more work needs to be done in the area.

References

[1] L. Baker, A. Pierce, and K. Luks. Gibbs energy analysis of phase equilibria. Society of Petroleum Engineers Journal, 22:731–742, 1982.

[2] J. Gibbs. Graphical methods in thermodynamics of fluids. Trans. Connecticut Acad, 2:311, 1873.

[3] J. Gibbs. A method of geometrical representation of the thermodynamic properties of substances by means of surfaces. Trans. Connecticut Acad, 2:382, 1873.

[4] E. Hansen. Global Optimization Using Interval Analysis, volume 165 of Pure and Applied Mathematics. Marcel Dekker, New York, 1992.

[5] Y. Jiang, W. Smith, and G. Chapman. Global optimality conditions and their geometric interpretation for the chemical and phase equilibrium problem. SIAM Journal on Optimization, 1994. To appear.

[6] C. McDonald and C. Floudas. Global optimization for the phase and chemical equilibrium problem: Application to the NRTL equation. Computers & Chemical Engineering, 1994. To appear.

[7] C. McDonald and C. Floudas. GLOPEQ: A new computational tool for the phase and chemical equilibrium problem. 1994. Submitted.

[8] K. I. McKinnon and M. Mongeau. A generic global optimization algorithm for the chemical and phase equilibrium. Technical Report MS94/1, Dept. of Mathematics and Statistics, University of Edinburgh, U.K., 1994.

[9] M. L. Michelsen. The isothermal flash problem. Part I. Stability. Fluid Phase Equilibria, 9:1–19, 1982.

[10] M. L. Michelsen. The isothermal flash problem. Part II. Phase-split calculation. Fluid Phase Equilibria, 9:21–40, 1982.

[11] H. Ratschek and J. Rokne. New Computer Methods for Global Optimization. Wiley, New York, 1982.

[12] W. Smith and R. Missen. Chemical Reaction Equilibrium Analysis: Theory and Algorithms. Wiley & Sons, 1982.

[13] J. A. Trangenstein. Customized minimization techniques for phase equilibrium computations in reservoir simulation. Chemical Engineering Science, 42(12):2847–2863, 1987.