Duc Thach

Son Vu

Ibtihel Ben Gharbia

Mounir Haddou

Quang-Huy Tran

A new approach for solving nonlinear algebraic systems with complementarity conditions. Application to compositional multiphase equilibrium problems

Keywords: complementarity condition, unified formulation, phase equilibrium, Newton's method, interior-point methods Mathematics subject classification 80M25, 90C33, 90C51

We present a new method to solve general systems of equations containing complementarity conditions, with a special focus on those arising in the thermodynamics of multicomponent multiphase mixtures at equilibrium. Indeed, the unified formulation introduced by Lauser et al. [Adv. Water Res. 34 (2011), 957-966] has recently emerged as a promising way to automatically handle the appearance and disappearance of phases in porous media compositional multiphase flows. From a mathematical viewpoint and after discretization in space and time, this leads to a system consisting of algebraic equations and nonlinear complementarity equations. Such a system exhibit serious convergence difficulties for the existing semismooth and smoothing methods. This observation led us to design a new strategy called NPIPM (Non-Parametric Interior-Point Method ). Inspired from interior-point methods in optimization, the technique we propose avoids any parameter management while ensuring good theoretical convergence results. These are validated by extensive numerical tests, in which we compare NPIPM to the Newton-min method.

Introduction 1.Motivation and objectives

In many applied scientific fields such as mechanics, electronics or chemical kinetics [START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems[END_REF][START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems, I[END_REF] we often encounter models of the form ΛpXq 0, R ¡m , (1.1a) minpGpXq, HpXqq 0, R m , (1.1b) * IFP Energies nouvelles, 1 et 4 avenue de Bois Préau, 92852 Rueil-Malmaison Cedex, France. thachsonnt94@gmail.com, ibtihel.ben-gharbia@ifpen.fr, quang-huy.tran@ifpen.fr where the unknown X D is to be found in some open domain D R and where the given functions Λ : D R Ñ R ¡m and G, H : D R Ñ R m , with 0 m ¤ , are assumed to be continuously differentiable on D. The first ¡ m equations (1.1a) are "ordinary" algebraic equations. By contrast, the last m equations (1.1b) are rather "special" in that they involve the componentwise minimum function and are therefore nondifferentiable. They represent the so-called complementarity conditions (or unilateral conditions), the exact significance of which is 0 ¤ GpXq u HpXq ¥ 0, or equivalently, GpXq ¥ 0, HpXq ¥ 0, xGpXq, HpXqy 0.

(1.2)

This name is justified by the observation that for each index α t1, . . . , mu, at least one of the two quantities G α pXq and H α pXq vanishes while the other remains nonnegative. Each complementarity condition minpG α pXq, H α pXqq 0 expresses two possible functioning regimes by a same single equation. The m complementarity conditions thus enable us to conveniently envision the 2 m configurations of the physical system in a unified manner. By "unified" we mean that a fixed set of equations and unknowns is assembled throughout the simulation, namely, F pXq 0, with F pXq ΛpXq minpGpXq, HpXqq & .

(1.3)

Complementarity conditions were introduced in thermodynamics by Lauser et al. [START_REF] Lauser | A new approach for phase transitions in miscible multi-phase flow in porous media[END_REF] for the phase equilibrium problem of compositional multiphase mixtures, with an aim to speed up the simulation of flows in porous media. The difficulty of this problem lies in the appearance and the disappearance of various phases (liquid, gas, oil), the handling of which is quite delicate. Reservoir engineers traditionally use the variable-switching formulation [START_REF] Coats | An equation of state compositional model[END_REF], in which only those unknowns and equations that correspond to a currently present phase are considered. This has the advantage of keeping the system small-sized, but is cumbersome and costly to implement, insofar as "switching" may occur all the time. In this respect, Lauser et al.'s unified formulation is a tremendous progress, not only for the practical comfort it offers but also for the theoretical properties that it encapsulates, as emphasized in [9].

Subsequent works by Ben Gharbia [START_REF] Gharbia | Résolution de problèmes de complémentarité. : Application à un écoulement diphasique dans un milieu poreux[END_REF], Ben Gharbia and Jaffré [10], Masson et al. [START_REF] Masson | Formulations of two phase liquid gas compositional Darcy flows with phase transitions[END_REF][START_REF] Masson | Coupling compositional liquid gas Darcy and free gas flows at porous and free-flow domains interface[END_REF] and Beaude et al. [START_REF] Beaude | Non-isothermal compositional liquid gas Darcy flow: formulation, soil-atmosphere boundary condition and application to high-energy geothermal simulations[END_REF] testified to the superiority of the unified formulation over the variable-switching one regarding computational time in simple cases, i.e., with Henry's law for fugacity coeffients. Another series of works at IFPEN [START_REF] Gharbia | Study of compositional multiphase flow formulation using complementarity conditions[END_REF][START_REF] Ben Gharbia | Study of compositional multi-phase flow formulations with cubic EOS[END_REF][START_REF] Lusetti | Numerical methods for compositional multiphase flow models with cubic EOS[END_REF][START_REF] Peton | Comparaison de plusieurs formulations pour les écoulements multiphasiques et compositionnels en milieu poreux[END_REF] seem to reach the opposite conclusion for realistic fugacity coefficients given by cubic equations of state, such as Peng-Robinson's law. All these works use the Newton-min method [START_REF] Aganagić | Newton's method for linear complementarity problems[END_REF][START_REF] Kräutle | The semismooth Newton method for multicomponent reactive transport with minerals[END_REF] to solve the nondifferentiable algebraic system of equations (1.1). It was observed that Newton-min may suffer from periodic oscillations or converge to a wrong solution. There are two possible explanations to this lack of robustness:

1. System (1.1) is ill-posed for some data and thermodyanimic laws. It may not have a solution or some components of Λ are not well-defined over the whole domain of interest D. This issue pertains to physical modelling and was addressed in [9].

2. Despite its popularity among numericist, thanks to its simplicity and the semismooth local convergence theory of [START_REF] Mifflin | Semismooth and semiconvex functions in constrained optimization[END_REF][START_REF] Qi | A nonsmooth version of Newton's method[END_REF], Newton-min may not be well-suited to the specific nature of our problem. We need another algorithm. This issue is the subject of this article.

This work is a first step in the direction of an alternative numerical method to solve (1.1), with a better guarantee of convergence and a greater robustness with respect to the parameters of the problem. Our quest is deeply rooted in the belief that the unified formulation has a strong potential to improve the performance of compositional multiphase flow simulators. After all, it is already a major advance to have succeeded in recasting the continuous problem under a unified language. It would be a pity not to "convert the try" for want of an adequate numerical method.

Main results and outline of the paper

First of all, it is essential to review the existing classical methods that could be considered for our problem. This is the purpose of §2. In addition to Newton-min ( §2.1), which is part of the semismooth family, we also revisit the smoothing family with the θ-smoothing methods ( §2.2) and the interior-point methods ( §2.3). Although we prefer this second category for its conceptual elegance, we insist on its main shortcoming, which is the absence of a reliable strategy to drive the regularization parameter towards zero.

The identification of this difficulty prompts us to propose a new interior-point method ( §3), in which the smoothing parameter becomes a full-fledged unknown. The latter is governed by a new equation coupling the new variable with the remaining ones. The augmented system is then solved with the ordinary smooth Newton method, with an optional Armijo-type line search to ensure global convergence. Since the parameter is now updated "automatically" and no longer "manually," the new algorithm is given the name NonParametric Interior-Point Method (NPIPM). This is our main contribution.

In preparation for the numerical tests, we consider two models giving rise to systems of the form (1.3). Both are reduced versions of more realistic models for two-phase compositional mixtures. The first one, expounded in §4.1, is the stationary phase equilibrium problem at a given global composition. A thorough analysis of it was conducted in [9], where we highlighted its good mathematical properties and clarified the assumptions ensuring the existence of a solution, detailed in §4.2. In §4.3, we prove the regularity of all solutions barring two degenerate cases, which is a favorable feature for quadratic convergence of NPIPM. This model is deployed in combination with Henry's and Peng-Robinson's fugacity laws. Thanks to the simplicity of this model, we were able to diagnose [9] an intrinsic deficiency of Peng-Robinson's law that makes it impossible for F to be well-defined over the whole domain D, thus jeopardizing the unified formulation. The domain extension procedure proposed in [9] and recalled in §4.4 are therefore systemically applied in order to prevent the numerical methods from stopping prematurely.

The second model, which includes an evolution in time, is built on top of the first one by prescribing an ordinary differential equation on the global mass fraction ( §5.1). Far from being arbitrary, this differential equation can be heuristically derived from a one-dimensional two-phase binary flow model discretized in space by a one-cell mesh. The interest of this second model lies in the new difficulty associated with the time-step. This is why it is deployed in combination with only Henry's law for fugacities, which also allows the exact solution to be determined ( §5.2) and proven to be regular ( §5.3).

Finally, extensive numerical results are provided in §6, where NPIPM is systematically compared to Newton-min. For each model, we sweep over the initial points and the parameters. This test campaign clearly demonstrates NPIPM's excellence in comparison with Newton-min: most of the time, NPIMP achieves a 100% convergence ratio.

Standard approaches for solving systems with complementarity conditions

For a general function F : D R Ñ R that is not differentiable everywhere, the numerical methods available in the literature to solve F pXq 0 can be grossly divided into two categories:

semismooth methods and smoothing methods.

Newton-min, a popular semismooth method

The semismooth approach generalizes the classical smooth Newton method by resorting to a broader notion of Jacobian matrix. Actually, the semismooth Newton theory [START_REF] Mifflin | Semismooth and semiconvex functions in constrained optimization[END_REF][START_REF] Qi | A nonsmooth version of Newton's method[END_REF]

D Ñ R ¢ defined as f B F pXq 2 M R ¢ | hpX n q kN C F , X n Ñ X, ∇F pX n q Ñ M @ . (2.1a)
In other words, the Bouligand subdifferential f B F pXq is the set of all matrices M are the limits of the Frechet differentials ∇F pX n q for a sequence X n converging to X. 2. The C-subdifferential or the generalized Jacobian of F at X is the set-valued mapping fF :

D Ñ R ¢ given by fFpXq convpf B F pXqq. (2.1b)
In other words, the Clarke subdifferential fFpXq is the convex hull of the Bouligand subdifferential f B F pXq. The idea is then to prescribe the local linear approximation d Þ Ñ F pX k q M k d around X k with M k fFpX k q and to look for a zero of this approximation scheme. This defines the sequence

X k 1 X k d k after having solved M k d k ¡FpX k q for d k .
Unfortunately, in general this linear local model does not satisfy the technical conditions required by the nonsmooth framework to ensure convergence [START_REF] Vu | Numerical resolution of algebraic systems with complementarity conditions[END_REF]Definition 4.5]. This is why we have to restrict ourselves to a subclass of those locally Lipschitz functions that comply exactly with these conditions. Definition 2.2 (Semismooth function). Let F : D R Ñ R be a locally Lipschitz-continuous function. We say that F is semismooth at s

X D if lim sup XÑ X M fFpXq F pXq M p s X ¡ Xq ¡ F p s Xq X ¡ s X 0.
(2.

2)

The subclass of semismooth functions is rich enough to cover all functions of interest in real applications. It happens, however, that the generic element of fFpX k q is difficult to identify. This is precisely the case for the functions F having the form (1.3). It is then recommended [START_REF] Izmailov | Newton-Type Methods for Optimization and Variational Problems[END_REF] to pick

M k f B F pX k q instead,
which might be easier to determine. The following statement provides the generic element of f B F pX k q when F is of the type

(1.3). Proposition 2.1. If Λ, G, H : D R Ñ R are continuously differentiable, then F is semis- mooth. Its B-subdifferential consists of all matrices M R ¢ of the form f B F pXq 4 M ∇ΛpXq { ∇ & , { ∇ R m¢ B , (2.3a)
where the α-th row of { ∇ for α t1, . . . , mu is

{ ∇ α 6 9 8 9 7 ∇G α pXq if G α pXq H α pXq, ∇G α pXq or ∇H α pXq if G α pXq H α pXq, ∇H α pXq if G α pXq ¡ H α pXq.
(2.3b)

Algorithm 1 Newton-min algorithm 1. Choose X 0 D R . Set k 0. 2. If F pX k q 0, stop. 3. Select an element M k f B F pX k q as in (2.3). Find a direction d k R such that F pX k q M k d k 0. (2.4) 4. Set X k 1 X k d k and k Ð k 1. Go to step 2.
Proof. The proof makes use of general results on the semismoothness of the componentwise minimum of two semismooth functions and on the B-subdifferential of a vector-valuded function [ Algorithm 1 summarizes the successive steps of the Newton-min method for problem (1.3). Historically, the first report on Newton-min seems to date back to Aganagić [START_REF] Aganagić | Newton's method for linear complementarity problems[END_REF]. Instead of the minimum function, we could have expressed componentwise complementarity by means of another

C-function, that is, a function ψ : R 2 Ñ R such that ψpa, bq 0 ô a ¥ 0, b ¥ 0, ab 0. (2.5)
Besides ψ min pa, bq minpa, bq, two other C-functions are worth knowing of: the Fischer-Burmeister function ψ FB pa, bq c a 2 b 2 ¡ pa bq. This C-function is differentiable everywhere except at p0, 0q. In addition, its square ψ 2 F B pa, bq is continuously differentiable on the entire plane. Introduced in [START_REF] Fischer | A special Newton-type optimization method[END_REF], the Fischer-Burmeister function played a central role in the early development of algorithms. The corresponding semi-smooth method is called Newton-FB.

the Mangasarian function ψ M pa, bq ζp|a¡b|q¡ζpaq¡ζpbq, where ζ : R Ñ R is an increasing function and ζp0q 0. It can be made differentiable everywhere by an appropriate choice of ζ, e.g., ζptq t 3 as advocated in [START_REF] Mangasarian | Equivalence of the complementarity problem to a system of nonlinear equations[END_REF]. The coresponding Newton-M method is nevertheless not very successful, due to a flaw shared by all smooth C-functions, that is, ∇ψp0, 0q p0, 0q.

This makes the Jacobian matrix singular at the solution.

Remark 2.1. It was observed [START_REF] Gharbia | Résolution de problèmes de complémentarité. : Application à un écoulement diphasique dans un milieu poreux[END_REF] that the direction d k computed by (2.4) is not always a descent direction for the least-squares merit function ΘpXq 1 2 F pXq 2 . Consequently, globalization of Newton-min by means of a line search (cf. §3.3) along d k remains a delicate issue.

θ-smoothing, a typical smoothing method

In contrast to semismooth methods, smoothing methods first try to regularize the F function, which introduces a regularization parameter that must be progressively pushed to zero. A regularization of F is the choice of a family of functions 2 r

F p¤; νq : D R Ñ R , ν ¡ 0 @ (2.6)
such that: (i) r

F p¤; νq is a smooth (continuously differentiable) function of X, for all ν ¡ 0; (ii) r F p¤; νq is continuous with respect to ν, in some functional sense; (iii) lim νÓ0 r F p¤; νq F p¤q, in some functional sense. Starting from a current pair of values pX k , ν k q, the overall strategy of a smoothing method is to 1. solve r F pX k 1 ; ν k q 0 in the unknown X k 1 by means of the smooth Newton method, using X k as the initial point; 2. decrease the regularization parameter from ν k to ν k 1 by some "rule of thumb." Start over the process until F pX k 1 q 0. The regularized system in Step 1 is supposedly easier to solve and may benefit from the enhanced properties of the smooth Newton method. Moreover, it is customary to do just one iteration in Step 1, with a view to lower the computational cost. This paradigm, known as the diagonal Newton, induces more approximation error but is of great practical interest.

To smooth out the scalar complementarity equation 0 ¤ a u b ¥ 0, let us first rewrite it under the plainly equivalent form

a ¥ 0, b ¥ 0, Spaq Spbq ¤ 1, (2.7)
where Sptq 5

0 if t 0, 1 if t ¡ 0 (2.8)
is the step function. The latter serves as an indicator of positive arguments t ¡ 0 over R . Since the step function is discontinuous, let us regularize it using the notion of a smoothing function introduced by Haddou and his coauthors [START_REF] Auslender | Asymptotic analysis for penalty and barrier methods in convex and linear programming[END_REF][START_REF] Haddou | A new class of smoothing methods for mathematical programs with equilibrium constraints[END_REF]. (2.9a) Furthermore, if θ can be defined for negative arguments t p¡T, 0q, with T ¡ 0, while remaining continuous, nondecreasing and concave, it is required that θptq 0 for t p¡T, 0q.

(2.9b)

The smoothing function is a "father" function, by compression of which regularized step functions will be generated. The two most common examples of smoothing functions are: the rational function θ 1 : p¡1, Vq Ñ p¡V, 1q defined by

θ 1 ptq t t 1 .
(2.10a) the exponential function θ 2 : R Ñ p¡V, 1q defined by θ 2 ptq 1 ¡ expp¡tq.

(2.10b) Definition 2.4 (θ-smoothing family). Let θ be a θ-smoothing function. The family of functions

2 θ ν ptq : θpt{νq, ν ¡ 0 @ (2.11)
is said to be the θ-smoothing family associated with θ.

Obviously, θ ν is a smooth function of t ¥ 0 for all ν ¡ 0. It is also continuous with respect to ν at each fixed t ¥ 0. From the defining properties (2.9), it can be readily shown that lim νÓ0 θ ν ptq Sptq, dt ¥ 0.

(2.12)

In other words, S is the limit of θ ν in the sense of pointwise convergence. The equivalence between 0 ¤ a u b ¥ 0 and (2.7) suggests us to impose a ¥ 0, b ¥ 0, θ ν paq θ ν pbq 1 (2.13) for ν ¡ 0, as a smooth approximation of (2.7). Replacing S by θ ν in (2.7) is logical. Replacing "¤" by "=" in (2.7) is motivated by the fact that we want an equality to be mounted into the system of equations. Let us examine the impact of (2.13) on the examples (2.10).

For the rational function (2.10a), we can easily prove the remarkable equivalence

θ 1 ν paq θ 1 ν pbq 1 ô ab ν 2 .
(2.14a)

The equality ab ν 2 appears to be a natural relaxation of ab 0. This smoothing paradigm will prevail in interior-point methods, with ν in the right-hand side instead of ν 2 .

For the exponential function (2.10b), we readily check the equivalence

θ 2 ν paq θ 2 ν pbq 1 ô ¡ν lnrexpp¡a{νq expp¡b{νqs 0, (2.14b)
In the left-hand side, the function min ν pa, bq : ¡ν lnrexpp¡a{νq expp¡b{νqs can be regarded as a smooth approximation of minpa, bq. This relaxation of minpa, bq 0 could have been worked out from the more well-known fact [START_REF] Beck | Smoothing and first order methods: A unified framework[END_REF] that max ν pa, bq : ν lnrexppa{νq exppb{νqs is a smooth approximation of maxpa, bq.

For F of the form (1.3), we consider the regularization family t r F p¤, νq, ν ¡ 0u, where

r F pX , νq ΛpXq ν pθ ν pGpXqq θ ν pHpXqq ¡ 1q & . (2.15)
Here, it is understood that θ ν operates componentwise on GpXq and HpXq, while 1 R m is the vector whose entries are all equal to 1. The premultiplication by ν of the second line in (2.15) is aimed at controlling the magnitude of their partial derivatives. Indeed, for all t ¥ 0, θ I ν ptq ν ¡1 θ I pt{νq can be seen to blow up when ν Ó 0, while νθ I ν ptq tends to the finite limit θ I p0q.

To complete the two-step procedure described at the beginning of §2.2, we also need some "rule of thumb" to steer ν to 0. This issue will be discussed at the end of §2.3.

Interior-point, a reference in optimization

Renowned for their efficiency in linear programming thanks to their polynomial complexity, interiorpoint methods [START_REF] Wright | The interior-point revolution in optimization: history, recent developments, and lasting consequences[END_REF] can be interpreted as regularization methods. We are mostly interested in primal-dual methods [START_REF] Wright | Primal-Dual Interior-Point Methods[END_REF], in which primal variables (initial unknowns) and duals (Lagrange multipliers) enjoy the same status. When one dissects a primal-dual inner-point method, one realizes that it is basically a method for solving the algebraic system of Karush-Kuhn-Tucker (KKT) optimality conditions. The fact that the system comes from a constrained minimization problem is ultimately of little importance in the method. This opens up the prospect of transposing these methods to the case of a general system containing complementarity conditions.

Let us consider the family of regularized problems r

F pX; νq 0, with r F pX; νq ΛpXq GpXq d HpXq ¡ ν1 & R , (2.16)
where ν ¥ 0 is the smoothing parameter, 1 R m is the vector whose components are all equal to 1, and d denotes Hadamard's componentwise product. It is usually more convenient to cast the previous sytem under the form FpX ; νq 0, with FpX ; νq

" " ! ΛpXq GpXq ¡ V HpXq ¡ W V d W ¡ ν1 ( 0 0 ) R 2m , (2.17)
where

X rX T ; V T ; W T s T D ¢ R m ¢ R m R 2m
is the augmented unknown and where the slack variables pV, W q R m ¢ R m are subject to the componentwise positivity conditions V ¥ 0, W ¥ 0.

(2.18)

Enlarging the size of the system and the number of unknowns does not change the determinant of the Jacobian matrix at the corresponding solution. Let us state this result for later use. Due to defintions (2.16) and (2.17), the Jacobian matrices ∇ X r F pX; νq and ∇ X FpX; νq do not depend on ν. For short, they will be denoted by ∇ r F pXq and ∇FpXq. Lemma 2.1. Let X D and X rX T ; V T ; W T s T D ¢ R m ¢ R m such that V GpXq and W HpXq. Then, det ∇FpXq det ∇ r F pXq.

(2.19)

Proof. The determinant of the Jacobian matrix of FpX ; νq is equal to

det ∇FpXq ∇ΛpXq 0 0 ∇GpXq ¡I m 0 ∇HpXq 0 ¡I m 0 I m d W I m d V ,
where the Hadamard product between a matrix and a vector is defined as the matrix whose each column is the Hadamard product between the corresponding column of the matrix and the vector. By linear combination of the last (block)-row with the second and third (block)-rows, we obtain

det ∇FpXq ∇ΛpXq 0 0 ∇GpXq ¡I m 0 ∇HpXq 0 ¡I m ∇GpXq d W ∇HpXq d V 0 0 .
By means of 2m row permutations, we end up with det ∇FpXq ∇ΛpXq 0 0

∇GpXq d W ∇HpXq d V 0 0 ∇GpXq ¡I m 0 ∇HpXq 0 ¡I m ∇ΛpXq ∇GpXq d W ∇HpXq d V .
For V GpXq, W HpXq, the last determinant coincides det r F p¤; νq. Note that the Lemma does not require X to be a solution of (2.16).

A primal-dual interior-point method strives to generate a sequence

tX k u kN I : 2 X pX, V, W q R 2m | V ¡ 0, W ¡ 0 @ , (2.20)
as well as an auxiliary sequence tν k u kN R ¦ such that pX k , V k , W k q Ñ p s X, Gp s Xq, Hp s Xqq and

ν k Ñ 0, where s X D is a zero of F r F p¤; 0q. Algorithm 2 describes a single-stage interior-point
method which consists of one Newton iteration (Step 3), followed by a truncation (Step 4) and an update for the regularization parameter (Step 6). Below are a few common empirical ways to progressively drive ν k to 0:

ν k 1 pν k q 2 ;
(2.23a)

ν k 1 pν k q 2 ;
(2.23b)

ν k 1 min 0.5 ν k , pν k q 2 ¨;
(2.23c)

ν k 1 min 0.5 ν k , pν k q 2 , xV k 1 , W k 1 y { m ¨.
(2.23d)

The geometric sequence (2.23a) has the advantage of going slowly to zero, which is recommended when ν k is still large. The power sequence (2.23b) goes quickly to zero, which is relevant when ν k is already small. The hybrid geometric-power sequence (2.23c) combines the advantages of the first two strategies. The hybrid geometric-power sequence involving a duality measure (2.23d)

Algorithm 2 Single-stage interior-point algorithm 1. Choose X 0 I. Set k 0, ν 0 xV 0 , W 0 y { m, γ 0.99. 2. If FpX k ; 0q 0, stop. 3. Find a direction d k pdX k , dV k , dW k q R 2m such that FpX k ; ν k q ∇FpX k q d k 0. (2.21) 4. Compute ς k p0, 1q such that X k ς k d k I by ς k γ ¤ arg max 2 ς r0, 1s | V k ς dV k ¥ 0, W k ς dW k ¥ 0 @ . (2.22) 5. Set X k 1 X k ς k d k . 6. Set ν k 1 one of the heuristic strategies (2.23). 7. Set k Ð k 1. Go to step 2.
allows the sequence to be reconnected to some current "reality" [20, §5]. Unfortunately, there is no universal magic formula to monitor the sequence of regularization parameter tν k u. A heuristic strategy that works fine with one problem may fail with another. We have to try several sequences tν k u before knowing which one is best suited to the problem at hand.

Design of a new nonparametric interior-point method

To compensate for the weakness of smoothing methods regarding the lack of a good strategy to decrease ν to 0, we undertake to look for another way of dealing with the regularization parameter. Our guiding principle is to treat it as an unknown in its own right.

NPIPM in a nutshell

In system (2.17), the status of the parameter ν is very distinct from that of the variable X. While X is computed by a Newton iteration, ν has to be updated in an ad hoc manner. Usually, progress occurs when two objects of different natures are put on an equal footing. This is why we feel that it would be judicious to incorporate the parameter ν into the variables X.

Let us consider the enlarged vector X rX T ; νs

T D ¢ R m ¢ R m ¢ R R 2m 1 . We seek a system of 2m 1 equations FpXq 0, with FpXq FpX ; νq f pX, νq & (3.1)
to be prescribed on X. The last equation f pX, νq 0 must be devised in such a way that every solution of (2.17) with ν 0 is also a solution of the enlarged system (3.1). To enforce a few "desirable" properties that will be enumerated later, we advocate

f pX, νq 1 2 V ¡ 2 1 2 W ¡ 2 u 2m 2 pxV, W y q 2 ην ν 2 , (3.2a)
where u and η are two positive parameters and

V ¡ 2 m α1 pV ¡ α q 2 , W ¡ 2 m α1 pW ¡ α q 2 . (3.2b)
To explain the rationale behind the choice (3.2), let us remind ourselves that our ultimate goal is to make ν equal 0, while ensuring the inequalities (2.18). Thus, it seems really natural to first consider f pX, νq ν. This construction turns out to be too naive. Indeed, if we start from some ν 0 ¡ 0 and solve the smooth system (3.1) by the smooth Newton method, since the last equation is linear, we end up with ν 1 0 at the first iteration. Once the boundary of the interior region is reached, we are "stuck" there.

To prevent ν from rushing to zero in just one iteration, we could set f pX, νq ν 2 , which is equivalent at the continuous level. At the level of Newton iterates, there is still a flaw: since ν 0 is now a double root of the last equation, quadratic convergence will be lost when ν k approaches 0! A remedy to this is to add a small linear term, that is, f pX, νq ην ν 2 , where η ¡ 0 is a small parameter. The price to be paid for recovering quadratic convergence is that there is now a spurious negative solution ν ¡η 0. This should not be a problem, however, if we start from a positive value for ν.

At this stage, system (3.1) is still not yet adequate. Indeed, the last equation is totally decoupled from the others. Everything happens as if ν is governed by a pre-determined sequence, generated a priori by the Newton iterates of the scalar equation ην ν 2 0, regardless of X. It is desirable to couple ν and X in a tighter way. The choice (3.2) does achieve this purpose:

As long as ν ¥ 0, the cancellation of f pX, νq implies V ¡ W ¡ 0, which amounts to saying that V ¥ 0 and W ¥ 0. Should a component of V or W become negative during the iteration, the penalty terms 1 2 V ¡ 2 and 1 2 W ¡ 2 are strong incentives for it to return into the interior domain I.

Thanks to upxV, W y q 2 {2m 2 , where u ¡ 0 is a small parameter, the update of ν depends on V and W . Without this term, we run again into the previous problem of ν "living its own separate life" inside I. Note that at the continuous level, we have xV, W y{m ν as a consequence of V d W ν1, but at the level of iterates, things can be very different.

The idea is now to apply the standard Newton method to the smooth system (3.1)-(3.2), which updates X and ν simultaneously. To enforce a globally convergent behavior, we also opt for Armijo's line search.

Determinant of the Jacobian matrix and parameter increment

Before writing down the new algorithm, let us mention some insightful properties regarding the determinant of the Jacobian matrix and the parameter increment when the current iterate lies in the interior region. Lemma 3.1. Let X I, where I is the interior region defined in (2.20). Let X rX T ; νs T for some ν R. Then, det ∇FpXq pη 2ν uxV, W y{mq det ∇FpXq.

(3.3)
If ν ¡ ¡η{2, the two Jacobian matrices are singular or nonsigular at the same time.

Proof. In the same fashion as in the proof of Lemma 2.1, the Jacobian matrix of the enlarged system (3.1)-(3.2) can be decomposed blockwise as

∇FpX q " " " " ! ∇ΛpXq 0 0 0 ∇GpXq ¡I m 0 0 ∇HpXq 0 ¡I m 0 0 I m d W I m d V ¡1 0 pV ¡ q T uxV, W y W T {m 2 pW ¡ q T uxV, W y V T {m 2 η 2ν ( 0 0 0 0 ) (3.4)
where V ¡ is the vector of components V ¡ α minpV α , 0q and similarly for W ¡ . If X I, then V ¡ W ¡ 0. Subtracting the product of uxV, W y 1 T {m 2 with the fourth (block)-row from the last row, we can eliminate the scalar products xV, W y ¥ 0 in the second and third columns, so that det ∇FpXq The interest of upxV, W y q 2 {2m in (3.2) can also be revealed by inspecting the increment dν k ν k 1 ¡ ν k of the Newton method (without truncation nor line search). The following statement shows that whenever ν k is lower than the central value xV k , W k y{m, either ν k increases or ν k decreases with a smaller magnitude than it would have done without upxV, W y q 2 {2m. Put another way, the parameter ν k cannot go too fast to zero without "waiting" for xV k , W k y{m. Proposition 3.1. For X k I, the Newton increment for the parameter is

∇ΛpXq 0 0 0 ∇GpXq ¡I m 0 0 ∇HpXq 0 ¡I m 0 0 I m d W I m d V ¡1 0 0 0 η 2ν uxV, W y{m . ( 3
dν k ¡ ην k pν k q 2 ¡ upxV k , W k y{2m ¡ ν k q xV k , W k y{m η 2ν k uxV k , W k y{m . (3.6)
In comparison with the value

dν k 0 ¡ ην k pν k q 2 η 2ν k (3.7)
that corresponds to u 0, we have either

dν k 0 ¤ 0 dν k |dν k 0 | or dν k 0 dν k ¤ 0 as soon as ν k xV k , W k y{m.
Proof. Applying the row transformations described the proof of Lemma 3.1 to the linear system ∇FpX k qd k ¡FpX k q and taking into account the right-hand side, we end up with (3.6). When u 0, the increment (3.6) degenerates to (3.7). Let us introduce

d k dX k dν k & ¡ ∇FpX k q ¡f ν F 0 η 2ν k uxV k , W k y{m & ¡1 FpX k ; ν k q ην k pν k q 2 ¡ a k & , where a k upxV k , W k y{2m ¡ ν k q xV k , W k y{m, hence
A k ην k pν k q 2 , a k upxV k , W k y{2m ¡ ν k q xV k , W k y{m, (3.8a) B k η 2ν k , b k uxV k , W k y{m, (3.8b) so that dν k ¡pA k ¡ a k q{pB k b k q and dν k 0 ¡A k {B k . Then, since b k ¡ 0 and B k ¡ 0, the inequality |dν k | |dν k 0 | occurs if and only if pA k ¡ a k qB k A k pB k b k q, which boils down to a k b k ¡ ¡ A k B k . (3.9)
This holds true, e.g., for a k ¥ 0, that is, ν k ¤ xV k , W k y{2m. But this also holds true when a k is negative but "not too much." The exact condition (3.9) reads

xV k , W k y{2m ¡ ν k ¡ ¡ ν k pη ν k q η 2ν k pν k q 2 η 2ν k ¡ ν k and reduces to xV k , W k y{m ¡ 2pν k q 2 {pη 2ν k q. Since 2pν k q 2 {pη 2ν k q ν k , the condition ν k xV k , W k y{m becomes sufficient for |dν k | |dν k 0 | to take place.

Globalized algorithm

Introduce the least-squares potential ΘpX q 1 2 FpXq 2 . A detailed description of NPIPM is given in Algorithm 3. The initial point X 0 pX 0 , ν 0 q must be an interior point, i.e., X 0 I. Furthermore, it is often taken at equilibrium, that is, V 0 GpX 0 q and W 0 HpX 0 q, so that the initial parameter ν 0 xV 0 , W 0 y{m has the correct order of magnitude.

Algorithm 3 Nonparametric interior point algorithm with Armijo line search 1. Choose X 0 pX 0 , ν 0 q, X 0 I, ν 0 xV 0 , W 0 y{m, κ p0, 1{2q, p0, 1q. Set k 0. 2. If FpX k q 0, stop. 3. Find a direction d k R 2m 1 such that FpX k q ∇FpX k qd k 0.

(3.10) 4. Choose ς k j k p0, 1q, where j k N is the smallest integer such that

ΘpX k j k d k q ¤ p1 ¡ 2κ j k q ΘpX k q. (3.11) 5. Set X k 1 X k ς k d k and k Ð k 1. Go to step 2.
Remark 3.1. There is no need to truncate the Newton direction d k to preserve positivity for V k 1 and W k 1 , since nonnegativity is guaranteed at convergence. However, the possibility for the iterates to get out of the interior region makes this method not strictly "interior-point." A truly interior-point variant can be cooked up by carrying out the damping (2.22) before Step 4.

Remark 3.2. The positive parameters η and u are chosen once and for all. They not need to be dynamically adapted in some "smart" way during the iterations. It is in this sense that the adjective nonparametric is to be understood.

The convergence of Algorithm 3 is governed by the general theory for the smooth Newton method [11, §6]. This global convergence theory heavily relies on the regularity of zeros.

Definition 3.1 (Regular zero). Let s X D R be a zero of F , that is, F p s Xq 0. If the Jacobian matrix ∇F p s Xq is nonsingular, s
X is said to be a regular zero of F .

The three items of the upcoming Theorem illustrate the conditions and the qualities of convergence of the algorithm. Item (i) corresponds to the behavior of the algorithm near a regular zero. Item (ii) states the rate of convergence in some particular situations. Item (iii) summarizes all of the possible scenarios when running the algorithm. Theorem 3.1. Let F : R 2m 1 Ñ R 2m 1 be a continuously-differentiable function.

(i) [Local analysis] Let s X be a regular zero of F. If X 0 is close enough to s X, then ς k 1 for all k, and X k Ñ s X superlinearly (and we recover the standard Newton method).

(ii) [Limit point] Let r X be a limit point of sequence tX k u. If ∇Fp r Xq is invertible, then r X is a regular zero of F. If r X is a regular zero of F , then ς k 1 for k big enough and X k Ñ r X superlinearly.

(iii) [General behavior] At least one of three possibilities below holds:

(a) FpX k q Ñ 0. (b) dpX k q is unbounded. (c)
The sequence tX k u converges to r X where ∇Fp r Xq is not invertible.

Proof. See [11, Theorem 6.9] or the condensed exposition of [START_REF] Vu | Numerical resolution of algebraic systems with complementarity conditions[END_REF]Theorem 5.2].

Stationary model for compositional two-phase mixtures

We are going to compare Newton-min and NPIPM on two models involving complementarity conditions for multicomponent two-phase mixtures. Designed as reduced versions of those commonly encountered in reservoir simulations, they enable us to gain valuable insights into what happens to the algorithms and to conduct a thorough numerical study with respect to a full range of parameters ( §6). The first model is the unified formulation for the stationary phase equilibrium problem investigated in [33, §2] and [9].

Phase equilibria in the unified formulation

We consider a mixture consisting of K ¥ 2 distinct components or species, labeled by the elements of the set K 2 I, II, . . . , K

@

. Each component i K may be present under at least one phase but at most two phases. To fix ideas, the phases are labeled by the elements of P tG, Lu, which respectively stand for Gas and Liquid. Let tc i u iK be nonnegative numbers such that °iK c i 1. Each c i represents the global fraction of species i. The vector c pc I , . . . , c K¡1 q Ω, called global composition, lies in

Ω 2 x px I , . . . , x K¡1 q R K¡1 | x I ¡ 0, . . . , x K¡1 ¡ 0, 1 ¡ x I ¡ . . . ¡ x K¡1 ¡ 0 @ . ( 4.1)
For any vector x px I , . . . , x K¡1 q s Ω, it is understood that x K 1 ¡ x I ¡ . . . ¡ x K¡1 r0, 1s is the fraction corresponding to the last component.

Each phase α P is characterized by a fundamental function g α : Ω Ñ R known as the (intensive) Gibbs free energy of the phase. We require g α to be as smooth as necessary in Ω and continuously extendable to fΩ. However, ∇g α may blow up on fΩ. From g α , we define K functions

µ j α : Ω Ñ R, j K, called chemical potentials by µ j α pxq g α pxq x∇g α pxq, δ j ¡ xy (4.2)
for x Ω, where the vector δ j pδ j,1 , δ j,2 , . . . , δ j,K¡1 q R K¡1 is made up of Kronecker's symbols.

The following statement summarizes some helpful identities between g α and µ i α .

Lemma 4.1 (Connection between Gibbs energy and chemical potentials). For all x Ω:

g α pxq K jI x j µ j α pxq; (4.3a) fg α fx j pxq µ j α pxq ¡ µ K α pxq, dj KztKu. (4.3b) 0 K iI x i α ∇ xα µ i α px α q. (4.3c) Proof. See [33, Lemma 2.1].
The first equation (4.3a) relates to Gibbs energy to the potentials. The second formula (4.3b) provides the gradient of the Gibbs energy from the potentials. The last identity (4.3c) is known as the Gibbs-Duhem condition. The potentials usually take the form

µ i α pxq lnpx i Φ i α pxqq, (4.4)
in which Φ i α is called the fugacity coefficient (or activity coefficients) of component i in phase α. Substituting the form (4.4) into (4.3a), we obtain

g α pxq K iI x i ln x i Ψ α pxq, (4.5)
where the first sum °K jI x j ln x j is the ideal part and

Ψ α pxq K iI x i ln Φ i α pxq, (4.6)
is the excess part. The relations between Ψ α and ln Φ i α are similar to those between g α and µ i α .

Lemma 4.2 (Connection between excess energy and fugacity coefficients). For all x Ω: ln Φ j α pxq Ψ α pxq x∇Ψ α pxq, δ j ¡ xy, dj K;

(4.7a) fΨ α fx j α pxq ln Φ j α pxq ¡ ln Φ K α pxq, dj KztKu; (4.7b) 0 K iI x i α ∇ xα tln Φ i α upx α q. (4.7c) Proof. See [33, Lemma 2.2].
A family of positive real-valued functions tΦ i α u pi,αqK¢P is said to be admissible if, for each α P, there exists a Gibbs energy function g α such that they are the fugacity coefficients. This implies, in particular, that the functions Φ i α satisfy the Gibbs-Duhem condition (4.7c). Given c Ω and an admissible family of fuagacity coefficients, the phase equilibrium problem in the unified formulation amounts to solving the nonlinear system

Y ξ i G p1 ¡ Y qξ i L ¡ c i 0, di KztKu, (4.8a)
ξ i G Φ i G px G q ¡ ξ i L Φ i L px L q 0, di K, (4.8b) min Y, 1 ¡ °jK ξ j G ¨ 0, (4.8c) min 1 ¡ Y, 1 ¡ °jK ξ j L ¨ 0, (4.8d) in the unknowns pY , tξ i G u iK , tξ i L u iK q R ¢R K ¢R K . for α tG, Lu.
Here, the fraction Y r0, 1s measures the overall importance of phase G in the mixture. As a result, 1 ¡ Y reflects the global presence of phase L. The quantities ξ α pξ I α , . . . , ξ K α q R K are called extended fractions of the components in phase α. If a phase α tG, Lu is present, that is, Y ¡ 0 or 1 ¡ Y ¡ 0, then the corresponding extended fractions sum to 1, as can be seen from the complementarity conditions (4.8c)-(4.8d). In such a case, they coincide with the familiar partial fractions.

The K ¡ 1 equations (4.8a) are material balance of component i KztKu. It seems that the material balance for the last species K is missing, but in fact it can be readily recovered by combining (4.8a) and (4.8c)-(4.8d). Therefore, it is redundant and must be left out of the problem statement. The K equalities (4.8b) express the extended thermodynamic equilibrium of each component across the phases. In (4.8b), it is very important to be aware of the fact that x α is the renormalized fraction vector defined as

x i α ξ i α °jK ξ α j . (4.8e) System (4.8) is of the form (1.3) with 2K 1, m 2, X pY, ξ I G , . . . , ξ K G , ξ I L , . . . , ξ K L q, ΛpXq " " " " " " " " ! Y ξ I G p1 ¡ Y qξ I L ¡ c I . . . Y ξ K¡1 G p1 ¡ Y qξ K¡1 G ¡ c K¡1 ξ I G Φ I G px G q ¡ ξ I L Φ I L px L q . . . ξ K G Φ K G px G q ¡ ξ K L Φ K L px L q ( 0 0 0 0 0 0 0 0 ) R 2K¡1 (4.9a)
and

GpXq Y 1 ¡ Y & R 2 , HpXq 1 ¡ ξ 1 G ¡ . . . ¡ ξ K G 1 ¡ ξ 1 L ¡ . . . ¡ ξ K L & R 2 , (4.9b)
Remark 4.1. The thermodynamic functions g α , Ψ α and Φ i α also depend on the pressure P and the temperature T. We voluntarily omitted to indicate this dependence, since the phase equilibrium problem (4.8) is set at fixed pP, Tq.

Existence and construction of a solution

We recall that the functions tg α u αP are smooth (say, twice differentiable), take finite values on the boundary fΩ but their gradients blow up there, i.e., lim xÑfΩ ∇g α pxq V. The latter is due to the presence of logarithms in the ideal parts of the Gibbs functions. In [9], we showed that under the following additional hypotheses, an explicit solution of (4.8) can be worked out.

Hypotheses 4.1. The gradient map ∇g α : Ω Ñ R K¡1 is surjective. Moreover, the Gibbs energy g α : Ω Ñ R is strictly convex, that is, it satisfies one of the two conditions below, which are equivalent [START_REF] Boyd | Convex Optimization, Berichte über verteilte messysteme[END_REF] for a twice differentiable function:

(a) For all px, yq Ω ¢ Ω with x $ y, d ∇g α pxq ¡ ∇g α pyq, x ¡ y h ¡ 0.

(4.10) (b) For all x Ω, the Hessian matrix ∇ 2 g α pxq is definite positive.

This solution, which may not be unique, is inspired from Gibbs' geometric construction for the binary case [START_REF] Deiters | High-pressure Fluid Phase Equilibria: Phenomenology and Computation[END_REF] pK 2q, which is depicted in Figure 1. Let g mintg G , g L u and let q g be the lower convex envelope of g on s Ω. By design, q g is a convex and continuous function. It can be shown [9, Lemma 3.2] that under Hypotheses 4.1, q g is differentiable at all interior point.

Thus, for c Ω, it makes sense to speak about the gradient ∇q gpcq and the tangent hyperplane, defined as the graph of the linearized expansion T c q gpxq q gpcq x∇q gpcq, x ¡ cy. We introduce s Γpcq 2 α P | h s x α Ω, g α ps x α q T c q gps x α q, ∇g α ps x α q ∇q gpcq @ (4.11)

as the set of thoses phases whose Gibbs function g α is tangent to the hyperplane T c q g. It can then be proven [9, Lemma 3.3] that: (i) s Γpcq $ r; (ii) for each α s Γpcq, the contact point s x α is unique;

(iii) if P ¤ K, then c convts x α u α s

Γpcq . The last property means that c is a convex combination of the contact points, that is, there exist t s Y α u α s Γpcq ¥ 0 such that °α s Γpcq be defined as above and set

s ξ i α s x i α ,
for pα, iq s Γpcq ¢ K.

(4.12a)

For pβ, iq Pz s Γpcq ¢ K, set s Y β 0, s x β r∇g β s ¡1 p∇q gpcqq, s ξ i β exprT c q gps x β q ¡ g β ps x β qss x i β . (4.12b)
This procedure supplies us with a solution of (4.8).

Proof. See [9, Theorem 3.6]. Note that r∇g β s ¡1 is well-defined thanks to Hypotheses 4.1.

Regularity of zeros

According to Theorem 3.1, the promise of global convergence for the NPIPM algorithm hinges on the regularity of the zeros of the system at hand. Put another way, if we could prove that the Jacobian matrix ∇Fp s Xq at a solution s X is nonsingular, this would be an auspicious sign of the adequacy of the NPIPM algorithm to the problem. To study the regularity of a zero of (4.8), we need to pay attention to two kinds of "singular" solution.

Definition 4.1 (Transition point). A solution p s

Y , s ξ G , s ξ L q R ¢ R K ¢ R K of (4.8) is said
to be a transition point when both arguments of one of the complementarity conditions vanish simultaneously, that is,

4 s Y 0, 1 ¡ iK s ξ i G 0 B or 4 s Y 1, 1 ¡ iK s ξ i L 0 B . (4.13)
In the two-phase framework, such a point marks the change in the nature of the solution, from a two-phase regime to a single-phase regime or vice-versa. Figure 2 illustrates two azeotropic situations in the binary case. Note that c alone is not responsible for azeotropy. It takes the two Gibbs functions to behave in a peculiar way to satisfy (4.14). If azeotropy occurs at c Ω, the solution of (4.8) is not unique: since the two contact points s

x α , α tG, Lu, coincide with each other, s Y could be replaced by any other value Y r0, 1s.

The following Theorem tells us that, under Hypotheses 4.1, a solution of (4.8) is regular except for the two previous pathological cases.

Theorem 4.2. Let s X p s X, s V , W , s
νq R 2K 6 be a solution of (3.1), (3.2) using the functions (4.9). Assume that s ν 0 and that the Gibbs energy functions g G and g L meet Hypotheses 4.1.

Then, s

X is a regular zero if and only if s

X is neither a transition point nor an azeotropic point.

Proof. The proof is based on a brute-force calculation. We first invoke Lemma 3.

Y ξ I G p1 ¡ Y qξ I L ¡ c I . . . Y ξ K¡1 G p1 ¡ Y qξ K¡1 G ¡ c K¡1 lnpξ I G Φ I G px G qq ¡ lnpξ I L Φ I L px L qq . . . lnpξ K G Φ K G px G qq ¡ lnpξ K L Φ K L px L qq ( 0 0 0 0 0 0 0 0 )
. This is because the logarithm is an increasing function and because we are at a solution, which allows us to take out the same positive factor on each line. To compute d

, we go through a long sequence of row and column linear combinations. A lot of simplifications occur thanks to the Gibbs-Duhem condition (4.3c). At the end of the process, we obtain the following cases:

If s Y 1, i.e., the solution is single-phase in G, then Y 0, i.e., the solution is single-phase in L, and the discussion is similar.

d 1 ¡ s σ L ps σ L q K det ∇ 2 g L ps x L q,
If s Y p0, 1q, i.e., the solution is two-phase, then

d s Y p1 ¡ s Y q det p1 ¡ s Y q∇ 2 g G ps x G q s Y ∇ 2 g L ps x L q $ ∆s x T pM G ps x G q M L ps x L qq∆s x, where ∆s x s x G ¡ s
x L and where both symmetric matrices

M G p∇ 2 g G q 1{2 3 I K¡1 ¡ I K¡1 1 ¡ s Y s Y p∇ 2 g G q ¡1{2 ∇ 2 g L p∇ 2 g G q ¡1{2 % ¡1 A p∇ 2 g G q 1{2 M L p∇ 2 g L q 1{2 3 I K¡1 ¡ I K¡1 1 ¡ s Y s Y p∇ 2 g L q ¡1{2 ∇ 2 g G p∇ 2 g L q ¡1{2 % ¡1 A p∇ 2 g L q 1{2
can be shown to be positive definite owing to strict convexity of g α . To cancel d requires ∆s x 0, i.e., s x G s

x L . This is precisely the characterization of an azeotropic solution.

Full details of the calculation can be found in [START_REF] Vu | Numerical resolution of algebraic systems with complementarity conditions[END_REF]Theorem 5.3].

Laws for fugacity coefficients

We now give two examples of fugacity coefficients and their associated Gibbs functions. The first one, due to Henry, is very simple and complies with Hypotheses 4.1. The second one, due to Peng-Robinson, is one of the most physically advanced laws by current standards and does not necessarily satisfy these assumptions.

Henry's law

If phase α tG, Lu is an ideal fluid, for which Ψ α 0, then it can be checked to fulfill Hypotheses 4.1. Henry's law is a straightforward generalization, with

Ψ α pxq K i1 x i ln k i (4.16)
where tk i u iK are positive constants, each of them embodying a property of the corresponding species. The fugacity coefficients calculated by (4.7a) are then ln Φ j α pxq ln k i , for all j K.

(4.17) This is why this law is also called the constant coefficients law. Again, it is not difficult to show [33, Proposition 3.1] that the Gibbs energy function g α associated with Henry's law fulfills Hypotheses 4.1 for all pk I , . . . , k K q pR ¦ q K .

Peng-Robinson's law

In Peng-Robinson's law, the fugacity coefficients of the two phases G, L are coupled with each other and, in a sense, are computed "simultaneously." Each component i K in a pure state is characterized by a pair of positive parameters a i (cohesion term) and b i (covolume). These are highly sophisticated functions of the pressure and the temperature, but at fixed pP, Tq can be considered as constants. This gives rise to a pair of postive dimensionless parameters

A i Pa i pRTq 2 , B i Pb i RT . ( 4

.18)

A multicomponent mixture is supposed to behave approximately as a fictitious pure component endowed with an averaged value for the pair pA, Bq. The latter is computed from the pA i , B i q's and the current partial fractions by means of a mixing rule x Þ Ñ pApxq, Bpxqq. There can be found [START_REF] Orbey | Modeling Vapor-Liquid Equilibria: Cubic Equations of State and Their Mixing Rules[END_REF] a wide variety of mixing rules . We choose the most common one, namely,

Apxq

¢ jK

x j c A j 2 , Bpxq jK x j B j . (4.19)
The next step is to consider the cubic equation

Z 3 pxq pBpxq ¡ 1qZ 2 pxq rApxq ¡ 2Bpxq ¡ 3B 2 pxqsZpxq rB 2 pxq B 3 pxq ¡ ApxqBpxqs 0 (4.20)
in the variable Zpxq. In the most favorable situation, there are three real roots, all greater than Bpxq. These are then named Bpxq Z L pxq Z I pxq Z G pxq. In other words, the smallest root is associated with the liquid phase L, while the largest one is associated with the gas phase G. As for the intermediate root Z I pxq, it does not have any physical meaning. Let α tG, Lu and assume that the real root Z α pxq ¡ Bpxq is well-defined. Then, the excess Gibbs energy of Peng-Robinson's law is defined as

Ψ α pxq Z α pxq ¡ 1 ¡ ln rZ α pxq ¡ Bpxqs ¡ Apxq 2 c 2Bpxq ln Z α pxq p1 c 2qBpxq Z α pxq ¡ p c 2 ¡ 1qBpxq & . (4.21)
From this, the fugacity coefficients can be deduced with the help of (4.7a using the "matrix-vector" product

A i pxq c A i ¡ jK x j c A j © . (4.23)
In the unfavorable situation when equation (4.20) has only on real root greater than Bpxq, two subcases have to be envisaged. If we manage to assign a "natural" phase label α G or L to the real root, then the corresponding excess Gibbs energy Ψ α is defined by (4.21), leaving its counterpart in the other phase undefined. Otherwise, Ψ α is undefined in both phases. This process raises two serious questions:

1. When does the cubic equation has three real roots greater than Bpxq and when does it have only one real root greater than Bpxq?

2. When and how can a "natural" phase label be assigned to the unique real root greater than Bpxq and when is it impossible?

These questions were addressed at length in [33, §3.2.3] (see also [9] for a shorter presentation), where we stressed the critical issue of the Gibbs functions g G , g L not being defined on the whole domain Ω. To sketch out the difficulty, let us consider the binary case where x x p0, 1q Ω. Then, for some data pA I , B I q and pA II , B II q, there exist 0 x S x U 1 such that the quantities Z L pxq, Ψ L pxq, g L pxq are well-defined only for x r0, x U s, while the quantities Z G pxq, Ψ G pxq, g G pxq are well-defined only for x rx S , 1s. Furthermore, since g I G px S q and g I L px ¡ U q are finite, the image sets g I G prx S , 1qq and g I L pp0, x U sq are not equal to R. This prevents the unified formulation from being always able to assigning a well-defined extended fraction to a vanishing phase by inverting g I G and g I L . To circumvent this obstacle, we proposed [33, §3.3] a procedure to extend the domains of defintion of the Gibbs functions to s Ω. When the cubic equation does not have three real roots, the idea is to use the common real part of the two complex conjugate roots, as a "surrogate" of the missing real root. Assume that Z α is the only real root greater than B of Peng-Robinson's cubic equation

Z 3 pB ¡ 1qZ 2 pA ¡ 2B ¡ 3B 2 qZ pB 2 B 3 ¡ ABq 0.
To alleviate notations, we do not explicitly indicate the dependency of A, B and Z on x. Let β be the other phase, that is, β L if α G and β G if α L. Since the sum of the three (complex) roots is equal to 1 ¡ B, the two remaining conjugate roots share the common real part

W β 1 ¡ B ¡ Z α 2 . (4.24)
It can then be proven that for physically reasonable values of B, we have

W β ¡ B and Z α W β if α L, W β Z α if α G.
These properties of W β are the constraints to which Z β would have been subject, had it existed. They speak in favor of the enrollment of W β as a substitute for Z β . Doing so yields the Gibbs function

Ψ β W β ¡ 1 ¡ ln W β ¡ B $ ¡ A 2 c 2B ln W β p c 2 1qB W β ¡ p c 2 ¡ 1qB & . (4.25)
for the missing phase β. By (4.7a), we can derive the corresponding fugacity coefficients. In view of the mixing rule (4.19), these are given by

ln Φ i β pxq B i B rW β ¡ 1s ¡ ln rW β ¡ Bs B i B ¡ 2A i pxq A & A 2 c 2B ln W β p c 2 1qB W β ¡ p c 2 ¡ 1qB & x∇W β , δ i ¡ xy W β ¡ x∇B, δ i ¡ xy B & W β Υ A,B pW β q pW β ¡ BqpW 2 β 2BW β ¡ B 2 q (4.26) for all i K, with Υ A,B pWq W 3 pB ¡ 1qW 2 pA ¡ 2B ¡ 3B 2 qW pB 2 B 3 ¡ ABq.
The last step of this procedure is a "sacrifice" of a tiny portion of the three-root region in order for ∇Z β not to blow up when approaching the one-root region. To this end, we introduce

ϑ Z I ¡ Z L Z G ¡ Z L r0, 1s (4.27)
as an indicator of the closeness to the transition boundary. Indeed, the cubic equation has double roots when ϑ 0 or ϑ 1. Let ε p0, 1{4q be a small threshold. If ϑ r2ε, 1 ¡ 2εs, we apply the usual formulas for the case of three real-roots. If ϑ p1 ¡ 2ε, 1s, the two roots Z I and Z G are close to each other. We keep

Z L but progressively replace Z G by W G 1 2 p1 ¡ B ¡ Z L q 1 2 pZ I Z G q whose gradient is bounded. Instead of Z G , we plug r Z G r1 ¡ ν G pϑqsZ G ν G pϑqW G into (4.21),
where

ν G pϑq 6 9 8 9 7 0 if ϑ ¤ 1 ¡ 2ε, qppϑ ¡ p1 ¡ 2εqq{εq if ϑ p1 ¡ 2ε, 1 ¡ εq, 1 if ϑ ¥ 1 ¡ ε, (4.28)
and qpyq y 2 p3 ¡ 2yq. The rescaled function y Þ Ñ qpy{εq serves as a C 1 step function over the interval r0, εs. We note that qp0q 0, qp1q 1 and q I p0q q I p1q 0. From the modified excess Gibbs energy

Ψ G r Z G ¡ 1 ¡ ln r Z G ¡ B $ ¡ A 2 c 2B ln r Z G p c 2 1qB r Z G ¡ p c 2 ¡ 1qB & (4.29a)
and from the rule (4.7a), the fugacity coefficients are inferred as

ln Φ i G B i B r r Z G ¡ 1s ¡ ln r r Z G ¡ Bs B i B ¡ 2A i pxq A & A 2 c 2B ln r Z G p c 2 1qB r Z G ¡ p c 2 ¡ 1qB & x∇ r Z G , δ i ¡ xy r Z G ¡ x∇B, δ i ¡ xy B & r Z G Υ A,B p r Z G q p r Z G ¡ Bqp r Z 2 G 2B r Z G ¡ B 2 q . ( 4.29b)
If ϑ r0, 2εq, we proceed in a similar and symmetric fashion to replace Z L by r Z L r1¡ν L pϑqsZ L ν L pϑqW L in the expression of Ψ L , while preserving Z G .

To conclude this section, we lay emphasis on the fact the domain extension procedure is aimed at improving the unified formulation's chance of "survival," which will be corroborated by the numerical experiements of §6. We do not strive to fulfill Hypotheses 4.1, since these assumptions may already be violated for the original unextended Gibbs functions.

Evolutionary model for binary two-phase mixtures

The second model on which comparisons between Newton-min and NPIPM are carried out is a simplified system with a time evolution. To build this model, we start from the stationary twophase binary equilibrium problem and impose a differentiel equation for the global fraction on the top of it.

A simplified ODE system

In the binary case of the stationary model (4.8), the global fraction of the first component c is a given data. Let us make it depend on time by considering the algebro-differential system

dc dt ¡ Y p1 ¡ Y qξ I G p1 ¡ 1{k I q 0, (5.1a)
Y ξ I G p1 ¡ Y qξ I G {k I ¡ c 0, (5.1b) minpY, 1 ¡ ξ I G ¡ ξ II G q 0, (5.1c) minp1 ¡ Y, 1 ¡ ξ I G {k I ¡ ξ II G {k II q 0 (5.1d) in the four unknowns pc, Y, ξ I G , ξ II G q, equipped with the intial condition pc, Y, ξ I G , ξ II G qpt 0q pc 0 , Y 0 , pξ I G q 0 , pξ II G q 0 q (5.2)
subject to the equilibrium relations

Y 0 pξ I G q 0 p1 ¡ Y 0 qpξ II G q 0 {k I ¡ c 0 0, (5.3a) minpY 0 , 1 ¡ pξ I G q 0 ¡ pξ II G q 0 q 0, (5.3b) minp1 ¡ Y 0 , 1 ¡ pξ I G q 0 {k I ¡ pξ II G q 0 {k II q 0.
(5.3c)

At fixed c, the last three equations of (5.1) are in fact the stationary equilibrium problem (4.8) with K 2, in which the variables ξ I L and ξ II L have been eliminated with the help of the equilibrum relations ξ I

G k I ξ I L , ξ II G k II ξ II
L . The latter mean that phase G is an ideal gas and phase L is governed by Henry's law.

For simplicity, we restrict ourselves to

k I ¡ 1 ¡ k II ¡ 0. Let K G , K L be the two constants k I p1 ¡ k II q k I ¡ k II : K G ¡ K L : 1 ¡ k II k I ¡ k II ¡ 0.
(5.4)

It can then be proven that the exact solution of (5.1)-( 5.3) is given by cptq 6 9 9 9 8 9 9 9 7

c 0 if c 0 r0, K L s, K G γ 0 expptq K L γ 0 expptq 1 if c 0 pK L , K G q, c 0 if c 0 rK G , 1s, (5.5)
where

γ 0 pc 0 ¡ K L q{pK G ¡ c 0 q.
c, k II p1 ¡ cqq if c r0, K L s, ¢ c ¡ K L K G ¡ K L , K G , k II p1 ¡ K L q if c pK L , K G q, p1, c, 1 ¡ cq if c rK G , 1s, (5.6)
which are none other than the solution of the stationary two-phase binary equilibrium problem (4.8) at fixed c. But our primary interest is the algebraic system that arises when we apply the Euler backward scheme to (5.1) with a time-step τ ¡ 0. This system reads

c ¡ c S ¡ τ ¢ 1 ¡ 1 k I ξ I G Y p1 ¡ Y q 0, (5.7a)
Y ξ I G p1 ¡ Y qξ I G {k I ¡ c 0, (5.7b) minpY, 1 ¡ ξ I G ¡ ξ II G q 0, (5.7c) minp1 ¡ Y, 1 ¡ ξ I G {k I ¡ ξ II G {k II q 0, (5.7d)
where c S denotes the value of c at the previous time-step. In (5.7), c S r0, 1s, τ ¡ 0 and k I ¡ 1 ¡ k II ¡ 0 play the role of parameters. System (5.7) is of the form (1.3) with 4, m 2,

ΛpXq c ¡ c S ¡ τ p1 ¡ 1{k I q ξ I G Y p1 ¡ Y q Y ξ I G p1 ¡ Y q ξ I G {k I ¡ c & , (5.8a)
and

GpXq Y 1 ¡ Y & , HpXq 1 ¡ ξ I G ¡ ξ II G 1 ¡ ξ I G {k I ¡ ξ II G {k II & .
(5.8b)

Reference solution

The upcoming Proposition deals with the exitence and uniqueness of a solution to (5.7).

Proposition 5.1. Let K L , K G be the constants defined by (5.4). Except for the case 3(b) in the enumeration below, system (5.7) has a unique solution ps c, s

Y , s ξ I G , s ξ II G q r0, 1s ¢ r0, 1s ¢ R ¢ R called reference solution.
1. If c S rK G , 1s, then the reference solution is in the G single-phase regime and given by

s c c S , s Y 1, s ξ I G c S , s ξ II G 1 ¡ c S .
(5.9)

2. If c S pK L , K G q, then the reference solution is in the two-phase regime and given by

s c K G K L 2 ¡ K G ¡ K L 2τ 4 1 ¡ 1 ¡ 2 K G K L ¡ 2c S K G ¡ K L τ τ 2 & 1{2 B .
(5.10)

The values of s Y , s ξ I G and s ξ II G are deduced from s c by formulas (5.6).

3. If c S r0, K L s, then the number

τ max K G K L ¡ 2c S K G ¡ K L d ¢ K G K L ¡ 2c S K G ¡ K L 2 ¡ 1 (5.11)
is well-defined and greater than or equal to 1.

(a) For τ τ max , the reference solution is in the L single-phase regime and given by

s c c S , s Y 0, s ξ I G k I c S , s ξ II G k II p1 ¡ c S q;
(5.12) (b) For τ ¥ τ max , in addition to (5.12) that we declare to be the reference solution, there are two spurious solutions (counted with multiplicity).

Proof. The last three equations of model (5.7) are exactly the stationary binary model (4.8).

Therefore, pY, ξ I G , ξ II G q can be expressed as functions of c by means of (5.6). In particular,

Y c ¡ K L K G ¡ K L 1 pK L ,K G q pcq
for all phase regimes, using the characteristic function 1. Inserting this into the first equation (5.7a) and invoking k I K G {K L , we obtain a scalar equation on c, namely,

c ¡ c S τ K G ¡ K L pc ¡ K L qpc ¡ K G q 1 pK L ,K G q pcq 0.
(5.13)

The rest of the proof relies on studying the function representing the left-hand side of the above equation. This part is not difficult and is left to the readers.

Regularity of zeros

The most significant result for this model is that the reference solution corresponds most of the time to a regular zero.

Theorem 5.1. For all τ ¥ 0, the reference solution of (5.7) defined in Proposition 5.1 gives rise to a regular zero for the NPIPM system, except at transitional and azeotropic points. Xq. In this case, we have

d 1 ¡τ∆ s ξ I p1 ¡ 2 s Y q ¡τ∆ s ξ I s Y p1 ¡ s Y q{ s ξ I G 0 ¡1 ∆ s ξ I s Y p1 ¡ s Y q{k I 0 0 1 ¡ s σ G ¡ s Y ¡ s Y 0 ¡1 s σ L p s Y ¡ 1q{k I p s Y ¡ 1q{k II , with ∆ s ξ I s ξ I G ¡ s ξ I L s ξ I G p1 ¡ 1{k I q, s σ G s ξ I G s ξ II G , s σ L s ξ I L s ξ II L .
Expanding the determinant with respect to the first column, we find

d d 0 ¡ τ d 1 , (5.14)
where

d 0 ∆ s ξ I s Y p1 ¡ s Y q{k I 0 1 ¡ s σ G ¡ s Y ¡ s Y ¡1 s σ L p s Y ¡ 1q{k I p s Y ¡ 1q{k II
is the determinant of the stationary binary model and

d 1 ∆ s ξ I p1 ¡ 2 s Y q ∆ s ξ I s Y p1 ¡ s Y q{ s ξ I G 0 1 ¡ s σ G ¡ s Y ¡ s Y ¡1 s σ L p s Y ¡ 1q{k I p s Y ¡ 1q{k II . Assume s Y 1, i.e.,
G s σ L 1, ξ G x G and ξ L x L . We get d 0 ∆s x I ¡ s Y ¡ s Y p s Y ¡ 1q{k I p s
Y ¡ 1q{k II can be expressed as a quadratic form and hence d 0 ¥ 0, with equality if and only if x G x L , namely, at an azeotropic point. Let us compute d 1 . By expanding with respect to its first column and by noticing that the, we obtain

d 1 p1 ¡ 2 s Y q d 0 .
Coming back to (5.14), we get

d d 0 r1 ¡ τ p1 ¡ 2 s Y qs.
From (5.10), we infer by (5.6) that

τ p1 ¡ 2 s Y q 1 ¡ 1 ¡ 2 K G K L ¡ 2c S K G ¡ K L τ τ 2 & 1{2 1.
Consequently, d has the same sign behavior as d 0 .

Numerical results

On the grounds of the previous models, we now compare NPIPM and the Newton-min method. For the stationary model of §4, we consider the binary and ternary cases. For each case, we consider two tests: one with the ideal and Henry's laws, the other with Peng-Robinson's law. The domain extension procedure (4.24)-(4.26) must be activated, otherwise both Newton-min and NPIPM may crash.

For each test, we display two figures: the first one represents the solution for various c Ω using the same initial point in both methods; the second one indicates not only the number of iterations to reach convergence, starting from this initial point, but also the percentage of elements within a generated set of initial points for which convergence occurs.

Unless otherwise specified, the stopping criterion is FpXq 10 ¡7 and the maximum number of iterations to be 50 in all tests. With NPIPM, the parameters are η 0.5, u 1, κ 0.4 and 0.99.

Stationary two-phase binary model

In the two-component case, we write c, x α instead of c I , x I α . Model (4.8a) then becomes α q for α tG, Lu.

Y ξ I G p1 ¡ Y qξ I L ¡ c 0, ξ I G Φ I G px G q ¡ ξ I L Φ I L px L q 0, ξ II G Φ II G px G q ¡ ξ II L Φ II L px L q 0, minpY ; 1 ¡ ξ I G ¡ ξ II G q 0, minp1 ¡ Y ; 1 ¡ ξ I L ¡ ξ II L q 0, (6.1a)

With Henry's law

Phase G is an ideal gas, phase L obeys Henry's law with k I 2, k II 0.5. Thanks to the extreme simplicity of (6.1c), we can even eliminate ξ I L , ξ II L from the equations, as was done in §5.1 for the evolutionary model.

In the first test, starting from the same initial point pY, ξ I G , ξ II G q p0.99, 0.67, 0.327q and sweeping over the grid of parameters c t0.01; 0.02; . . . ; 0.99u, we run the two algorithms and observe the computed solution in Figure 3. Barring a single divergence point (red dot at c 0.4q for Newton-min, the two methods find the same solution, whose Y -component is drawn. The number of iterations required for each method is shown in panel (b) of Figure 4b.

In the second test, we sweep over the grid of parameters c t0.0001, 0.0002, ..., 0.9999u and the set of initial points

D 0 2 pY, ξ I G , ξ II G q 0 M 3 | 1 ¡ pξ I G q 0 ¡ pξ II G q 0 ¡ 0 and 1 ¡ pξ I G q 0 {k I ¡ pξ II G q 0 {k II ¡ 0 @
, where M t0.1; 0.2; . . . ; 0.9u. The number of initial points used for the tests is |D 0 | 216.

For each c, we count the number of initial points for which the method converges and then plot the percentage of success for each algorithm in Figure 4a. The success rate of NPIPM is 100%, while that of Newton-min is around 90%. In this Figure, we also plot the percentage of success corresponding to the θ-smoothing method using θ 1 and the Mehrotra predictor-corrector algorithm. We see that they are significantly worse than Newton-min and NPIPM. This is the reason why we keep only Newton-min as a reference for comparison.

With Peng-Robinson's law

Let us choose pA I , B I q p0.2153, 0.03q and pA II , B II q p0.1861, 0.02q to ensure that the Gibbs functions g G , g L are "visually" strictly convex once extended.

In the first test, we use the same initial point pY, ξ I G , ξ II G , ξ I L , ξ II L q 0 p0.2, 0.4, 0.4, 0.6, 0.2q and sweep over the grid of parameters c t0.01; 0.02; . . . ; 0.99u. The computed solutions are displayed in Figure 5. Note that Newton-min now has many divergence points (red dots) for c ¤ 0.5. The number of iterations required to reach convergence is supplied Figure 6a.

In the second test, we measure the percentage of convergence over the grid of parameters c t0.0001, 0.0002, ..., 0.9999u and the set of initial points

D 0 2 pY, ξ I G , ξ II G , ξ I L , ξ II L q 0 M 5 | 1 ¡ pξ I G q 0 ¡ pξ II G q 0 ¡ 0 and 1 ¡ pξ I L q 0 ¡ pξ II L q 0 ¡ 0 @ ,
where M t0.2; 0.4; 0.6; 0.8u. The number of initial points used for the tests is |D 0 | 144. The success rate for each method shown in Figure 6b as a function of c. Again, NPIPM achieves 100%, while Newton-min culminates at about 85%.

Stationary two-phase ternary model

In the three-component case, model (4.8) reads α q for α tG, Lu, i tI, IIu.

Y ξ I G p1 ¡ Y qξ I L ¡ c I 0, (6.2a) Y ξ II G p1 ¡ Y qξ II L ¡ c II 0, (6.2b) ξ I G Φ I G px I G , x II G q ¡ ξ I L Φ I L px I L , x II L q 0, (6.2c) ξ II G Φ II G px I G , x II G q ¡ ξ II L Φ II L px I L , x II L q 0, (6.2d) ξ III G Φ III G px I G , x II G q ¡ ξ III L Φ III L px I L , x II L q 0, (6.2e) minpY, 1 ¡ ξ I G ¡ ξ II G ¡ ξ III G q 0, (6.2f) minp1 ¡ Y, 1 ¡ ξ I L ¡ ξ II L ¡ ξ III L q 0, ( 6

With Henry's law

Phase G is an ideal gas, phase L obeys Henry's law with k I 0.2, k II 6, k III 2. The stopping criterion is FpXq 10 ¡12 . Starting from the same initial point pY, ξ I G , ξ II G , ξ III G q p0.9, 0.1, 0.7, 0.1q and sweeping over the grid of parameters C tpc I , c II q P 2 | c I c II 1u where P t0.01; 0.02; . . . ; 0.99u, we display the regime type (single-phase G, single-phase L or two-phase) of the computed solution in Figure 7, highlighting divergence in the red color. In Figure 8a, we display the number of iterations after linearly indexing the elements of C.

In the second test, we sweep over the same grid of parameters and the set of initial points

D 0 2 pY, ξ I G , ξ II G , ξ III G q 0 M 4 | 1 ¡ pξ I G q 0 ¡
pξ II G q 0 ¡ pξ III G q 0 ¡ 0 and 1 ¡ pξ I G q 0 {k I ¡ pξ II G q 0 {k II ¡ pξ III G q 0 {k III ¡ 0 @ , where M t0.1; 0.2; . . . ; 0.9u. The number of initial points used for the tests is |D 0 | 252. For each c C, we count the number of initial points for which the method converges and then plot the percentage of success for each algorithm in Figures 8b-8c. Figure 8c demonstrates the efficiency of NPIPM relatively to Newton-min, with 100% of convergence.

With Peng-Robinson's law

The stopping criterion is FpXq 10 ¡10 and η 10 ¡4 in the last equation of the NPIPM system. We select pA I , B I q p0.0883, 0.01q, pA II , B II q p0.1861, 0.02q, and pA III , B III q p0.2153, 0.03q so that g G and g L are "visually" strictly convex once extended.

In the first test, we use the same initial point pY, ξ I G , ξ II G , ξ III G , ξ I L , ξ II L , ξ III L q 0 p0.4, 0.3, 0.5, 0.1, 0.325, 0.2, 0.17q and sweep over the grid of parameters C tpc I , c II q P 2 | c I c II 1u, where P t0.01; 0.02; . . . ; 0.99u.

The regime type (single-phase G, single-phase L or two-phase) of the computed solution is represented in Figure 9, where divergence is reported in red. In Figure 10a, we display the number of iterations.

In the second test, we sweep over the grid of parameters C t pc I , c II q P 2 | c I c II 1 u, where P t0.05; 0.10; . . . ; 0.95u, and the set of initial points D 0 2 pY, ξ I G , ξ II G , ξ III G , ξ I L , ξ II L , ξ III L q 0 M 7 | 1 ¡ pξ I G q 0 ¡ pξ II G q 0 ¡ pξ III G q 0 ¡ 0 and 1 ¡ pξ I L q 0 ¡ pξ II L q 0 ¡ pξ III L q 0 ¡ 0 @ , where M t0.1; 0.2; . . . ; 0.9u. The number of initial points used for the tests is |D 0 | 64. For each c C, we plot the percentage of successful initial points for in Figures 10b-10c. Once again, we notice that NPIPM is much more effective than Newton-min.

Evolutionary two-phase binary model

Setting pk I , k II q p2, 0.5q, we sweep over the two-dimensional grid of parameters pc S , τ q t0.01; 0.02; . . . ; 0.99u¢ 2 0.1; 0.2; . . . ; 10 @ and the set of initial points D 0 2 pY, ξ I G , ξ II G , cq 0 M 4 | 1 ¡ pξ I G q 0 ¡ pξ II G q 0 ¡ 0 and 1 ¡ pξ I G q 0 {k I ¡ pξ II G q 0 {k II ¡ 0 @ where M t0.1; 0.2; . . . ; 0.9u. The number of initial points used for the test is |D 0 | 1944. For each pair pc S , τ q, we plot the percentage of successful initial points. The results are shown in Figure 11. Although NPIPM no longer claims a 100% rate, especially for large τ , it remains obviously better than Newton-min for τ ¤ 6 or c S ¤ 0.35.

Conclusion

In this work, we laid the first stone in the construction of a new method that would be better suited to the resolution of systems containing complementary conditions arising in the unified formulation of thermodynamical equilibria. Numerical results on simple compositional two-phase models revealed an overwhelming superiority of NPIPM over Newton-min. In this sense, they are very promising. The deliberately small size of the systems enabled us to understand in depth the behavior of the algorithms, and most notably the physical obstruction due to Peng-Robinson's law, for which the remedy (4.24)-(4.26) saved us from a compromising situation.

There is still a long way to go. It is of course essential to continue the comparison between NPIPM and Newton-min on the models presented here with a larger number of components, for example a dozen or even a hundred. It is also crucial to try the proposed method on realistic reservoir simulations, in which the system to be solved comes from the finite-volume discretization of the partial differential equations of a porous media flow model. The size of the problem would then be much larger, which could create additional difficulties for the algorithms.

Definition 2 . 3

 23 (θ-smoothing function). A function θ : R Ñ r0, 1q is said to be a θ-smoothing function if it is continuous, nondecreasing, concave, and θp0q 0, lim tÑ V θptq 1.

Remark 2 . 2 .

 22 Most of today's interior-point general-purpose softwares for linear programming are based on Mehrotra's predictor-corrector algorithm[START_REF] Mehrotra | On the implementation of a primal-dual interior point method[END_REF]. This celebrated two-stage algorithm rests upon a subtle tuning of the regularization parameter. It can be similarly extended to an equation solver[33, §4.3.3.3], but the numerical results are not always as good as the other methods [33, §6.1.1].

. 5 )

 5 Expanding(3.5) with respect to the last row yields the desired result. Note that this Lemma does not require X to solve (3.1)-(3.2).

Γpcq s Y α 1

 1 and °α s Γpcq s Y α s x α c. The remaining unknowns are determined as follows.

Figure 1 :

 1 Figure 1: Gibbs' geometric construction for the phase equilibrium of a two-phase binary mixture. Left: two-phase solution. Right: single-phase solution.

Theorem 4 . 1 (

 41 Geometric solution). Assume K ¥ 2, c Ω. Let ts x α u α s Γpcq , t s Y α u α s

Definition 4 . 2 (

 42 Azeotropic composition). A global composition c Ω is said to be azeotropic if the graphs of g G and g L are tangent to each other at c. In other words, g G pcq g L pcq, ∇g G pcq ∇g L pcq.

Figure 2 :

 2 Figure 2: Azeotropic compositions for a two-phase two-component mixture.

Proof.G

  Let X pc, Y, ξ I G , ξ II G q et recall the notations of (5.8). By Lemma 3.1 and Lemma 2.q is the reference solution, instead of the sign of det ∇Fp s Xq or det ∇Fp s

  renormalization x α ξ I α {pξ I α ξ II

Figure 3 :

 3 Figure 3: Henry's law, tested with the same initial point.

  Average number of iterations: 3.88 (Newton-min), 3.80 (NPIPM) Newton-min NPIPM (b) Number of iterations with the same initial points.

Figure 4 :

 4 Figure 4: Henry's law

Figure 5 :

 5 Figure 5: Peng-Robinson's law: one initial point.

  Percentage of convergence over all initial points.

Figure 6 :

 6 Figure 6: Peng-Robinson's law.

Figure 7 :

 7 Figure 7: Henry's law: one initial point.

  Number of iterations with the same initial point. NPIPM: percentage of convergence over all initial points.

Figure 8 :

 8 Figure 8: Henry's law.

  Number of iterations with the same initial point. NPIPM: percentage of convergence over all initial points.

Figure 10 :

 10 Figure 10: Peng-Robinson's law.

  is a concrete embodiment of the nonsmooth Newton theory [17, §7.2], an abstract framework that was developed much earlier for locally Lipschitz functions. By Rademacher's theorem [14, §3.4.1], every locally Lipschitz-continuous function is continuously differentiable almost everywhere. Put another way, the set C F of points X D where ∇F pXq exists in the classical sense is non-empty and its complement DzC F has measure zero. This enables two notions of subdifferential to be introduced. Definition 2.1 (Bouligand and Clarke subdifferentials). Let F : D R Ñ R be a locally Lipschitz-continuous function and C F D be the set of points at which F is differentiable. 1. The B-subdifferential or the limiting Jacobian of F at X is the set-valued mapping f B F :

  21, §1.75, §1.54]. See details in [33, Proposition 4.2].

  The values of Y ptq, ξ I G ptq and ξ II G ptq are deduced from cptq by

	the formulas
	p s Y , s ξ I G , s ξ II G q	9 9 9 9 7 6 9 9 9 9 8	p0, k I

  the solution is in the gas phase. Then, s σ G 1. We obtain that d 0 1¡ s Y p0, 1q, i.e., the solution is in the two-phase regime. Then, s σ

	σ L . σ L ¥ 0. Equality holds at a By a direct computation, we have d 1 0. Therefore, d d 0 1 ¡ s transition point. The other single-phase case s Y 0 is similar.
	Assume now s

[10] I. Ben Gharbia and J. Jaffré, Gas phase appearance and disappearance as a problem with complementarity constraints, Math. Comput. Simul., 99 (2014), pp.