Alexandre Ern
email: ern@cermics.enpc.fr

Rémi Joubaud
email: joubaudr@cermics.enpc.fr

Tony Lelièvre
email: lelievre@cermics.enpc.fr

Mathematical study of non-ideal electrostatic correlations in equilibrium electrolytes

We undertake the mathematical analysis of a model describing equilibrium binary electrolytes surrounded by charged solid walls. The problem is formulated in terms of the electrostatic potential and the ionic concentrations which have prescribed spatial mean values. The free energy of the system is decomposed as the difference of the internal energy and entropy functionals. The entropy functional is the sum of an ideal entropy and an excess entropy, the latter taking into account nonideality due to electrostatic correlations at low ionic concentrations and steric exclusion effects at high ionic concentrations. We derive sufficient conditions to achieve convexity of the entropy functional, yielding a convex-concave free energy functional. Our main result is the existence and uniqueness of the saddle point of the free energy functional and its characterization as a solution of the original model problem. The proof hinges on positive uniform lower bounds for the ionic concentrations and uniform upper bounds for the ionic concentrations and the electrostatic potential. Some numerical experiments are presented in the case where the excess entropy is evaluated using the Mean Spherical Approximation.

Introduction

The motivation for this work is the underground storage of radioactive waste in clay host rocks, whereby a thorough understanding of the clay behavior is of paramount importance. Clays are complex multiscale media. At the scale of one to several nanometers, they can be viewed as a porous medium in which a negatively charged solid matrix confines an electrolyte consisting of a solvent (water) and dissolved ionic species. Confined electrolytes are encountered in many other fields of application with significant interest, e.g., semiconductor devices and chemical engineering. In this work, we consider equilibrium electrolytes, that is, at steady state and no flow. The properties of such systems constitute the first step to understand more complex chemical and mechanical behaviors. For a recent numerical analysis of the coupling with the flow equations, we refer to Prohl and Schmuck [START_REF] Prohl | Convergent finite element discretizations of the Navier-Stokes-Nernst-Planck-Poisson system[END_REF].

In continuum models, equilibrium electrolytes can be described by the electrostatic potential and the ionic concentrations. For nanometric confinements, the classical Poisson-Boltzmann theory, where the electrostatic potential ψ solves a Poisson equation while the ionic concentrations (c i ) 1≤i≤M of the M species follow the Boltzmann distribution, is not sufficient. A first possibility is to resort to direct molecular simulations (see, e.g, Dufrêche et al. [START_REF] Dufrêche | Molecular hydrodynamics for electro-osmosis in clays: from Kubo to Smoluchowski[END_REF]) with suitable interaction potentials, but, at nanometric scales and beyond, this approach entails substantial computational costs. An alternative approach is to enrich the Poisson-Boltzmann theory by a more elaborate microscopic description. In what follows, we focus more specifically on ion-ion interactions. In the so-called primitive model, ions are described by charged hard-spheres of diameter (σ i ) 1≤i≤M , and the correlations between the ions are evaluated using the Ornstein-Zernike integral equation linking the direct and indirect correlation functions (see, e.g, the reference textbook [START_REF] Hansen | Theory of simple liquids[END_REF]). This allows one to derive the excess free-energy of the equilibrium system. The Ornstein-Zernike integral equation needs a closure relation to be solved. Various approaches can be adopted. Without being exhaustive, we mention the Hyper-Netted Chain approximation (HNC), which is a numerical approach, and the Mean Spherical Approximation (MSA) derived by Lebowitz et al. [START_REF] Lebowitz | Mean spherical model integral equations for charged hard spheres. I. Method of solution[END_REF][START_REF] Lebowitz | Mean spherical model integral equations for charged hard spheres. II. Results[END_REF] and Blum [START_REF] Blum | Mean spherical model for asymmetric electrolytes[END_REF], which is analytical in some physical cases. The interesting feature of the MSA is that it generalizes the Debye-Hückel theory of electrolytes to the case of nonzero ion diameters. Brownian dynamics simulations (where the solvent is a continuum but the ions are still charged explicit particles) within the primitive model have been performed by Jardat et al. [START_REF] Jardat | Salt exclusion in charged porous media: a coarse-graining strategy in the case of montmorillonite clays[END_REF] and compared to direct molecular simulations and to continuum models within the MSA.

The goal of the present work is to analyze mathematically equilibrium electrolytes including non-ideal effects due to electrostatics correlations at low to intermediate ionic concentrations and steric exclusion effects at high ionic concentrations. We consider binary electrolytes consisting of a counter-ion (cation) and a co-ion (anion); we denote the respective concentrations by c := (c + , c -). Non-ideal effects are described by activity coefficients in the electrochemical potentials. Electrostatic correlations are modelled through an activity coefficient log γ 0 which depends on a so-called screening length which itself can be evaluated in terms of the ionic concentrations. This model includes in particular the case where the MSA is used to evaluate the screening length (see §2.2 below for further description). Additionally, we account for a hard-sphere contribution which is modelled here as the lowest order term in the Percus-Yervick approximation, leading to a quadratic term in the ionic concentrations.

The mathematical study of critical points of the free-energy of equilibrium electrolytes has been addressed previously in the literature, either in the Poisson-Boltzmann setting or some of its generalizations. In [START_REF] Li | Minimization of electrostratic free energy and the Poisson-Boltzmann equation for molecular solvation with implicit solvent[END_REF], Li analyzes the generalized Poisson-Boltzmann theory with implicit solvent. This formulation, also considered by Borukhov, Andelman, and Orland [START_REF] Borukhov | Steric effects in electrolytes: A modified Poisson-Boltzmann equation[END_REF], accounts for steric exclusion effects, but not for electrostatic correlations. For an electrolyte with M species, the solvent concentration c 0 is introduced such that

σ 3 c 0 := 1 - M i=1 σ 3 c i ,
where σ > 0 represents the mean ionic diameter, and the ideal contribution of the solvent concentration, c 0 (log(σ 3 c 0 ) -1), is included in the entropy functional. The mathematical analysis has been extended to different ionic diameters by Li [START_REF] Li | Continuum electrostatics for ionic solutions with non-uniform ionic sizes[END_REF]. The critical points of the free-energy are sought in a convex set enforcing c i ≥ 0 for all 0 ≤ i ≤ M , so that the ionic concentrations are a priori bounded from above. One important result in the analysis of [START_REF] Li | Continuum electrostatics for ionic solutions with non-uniform ionic sizes[END_REF][START_REF] Li | Minimization of electrostratic free energy and the Poisson-Boltzmann equation for molecular solvation with implicit solvent[END_REF] is the proof that these constraints are not active, that is, that all the c i 's, 0 ≤ i ≤ M , are bounded uniformly away from zero. The technique of proof, which consists in further optimizing the free-energy by modifying the ionic concentrations at extreme values if the above abounds are not satisfied, will be adapted here to the present context. This extension is not straightforward owing to the coupling of all ionic concentrations in the non-ideal terms. Furthermore, the work of Carlen et al. [START_REF] Carlen | Free energy minimizers for a two-species model with segregation and liquid-vapour transition[END_REF] addresses a similar model using the Green operator for the Poisson equation. The analysis is concerned with nonlinearities of the same kind as those described by the implicit solvent. The mathematical analysis shows that L ∞ -bounds on the concentrations can be enforced in the convex set where the minimization is taken, but also shows that under a doubling condition on both the ideal and steric exclusion terms, L ∞ -bounds can be achieved [START_REF] Carlen | Free energy minimizers for a two-species model with segregation and liquid-vapour transition[END_REF]Theorem B.1]. Moreover, the work of Looker treats the case of the Poisson-Boltzmann free-energy in the ψ-formulation [START_REF] Looker | Semilinear elliptic neumann problems with rapid growth in the nonlinearity[END_REF]. The same model is considered by Allaire, Mikelić, and Piatniski [START_REF] Allaire | Homogenization of the linearized ionic transport equations in rigid periodic porous media[END_REF] in view of homogenization of confined electrolytes coupled with Stokes flow.

Our main result, Theorem 1 below, states that there is a unique saddle point (ψ, c) of the free-energy for an equilibrium binary electrolyte confined by charged walls. The electrostatic potential ψ is sought in the Sobolev space H 1 with zero mean-value, while the ionic concentrations c = (c + , c -) are sought in the closed convex subset of L 2 × L 2 consisting of nonnegative ionic concentrations with prescribed meanvalue. Moreover we prove that ψ and c ± are in L ∞ (Ω), and that c ± are uniformly bounded away from zero. Theorem 1 is established under four main assumptions stated in §3.1. The first two assumptions are the classical bulk electro-neutrality condition on the prescribed mean ionic concentrations and elliptic regularity for the Poisson problem governing the electrostatic potential (with non-homogeneous Neumann boundary conditions). The other two conditions are formulated in an abstract setting for the activity coefficient log γ 0 describing non-ideal electrostatic correlations. These conditions, which in particular encompass the MSA setting, require a sublinear growth condition at large ionic concentrations for log γ 0 , and a lower bound on the derivative of log γ 0 with respect to the ionic strength. This last condition is important to assert the convexity of the entropy functional. Indeed, the ideal entropy and the hard-sphere contribution, which dominate respectively at very low and large ionic concentrations, are both convex, but the electrostatic correlations lead to a nonconvex contribution to the entropy functional. Note that we do not need to enforce a priori L ∞ -bounds on the concentrations. To the best of our knowledge, the present analysis, together with the ongoing work [START_REF] Allaire | Role of non-ideality for the ion transport in porous media: derivation of the macroscopic equations using upscaling[END_REF] investigating the role of non-ideality for homogenized ion transport in porous media, is the first to address mathematically the critical points of the free-energy in the presence of electrostatic correlations.

This work is organized as follows. In §2, we give some details on the electrochemical model, including the modelling of non-ideality, the non-dimensionalization of the governing equations, and the free-energy and entropy functionals. Then, in §3, we state the mathematical assumptions together with our main result, and we present the main steps of its proof. We devote §4 to the proof of various technical lemmas, dealing in particular with the convexity of the entropy functional and the a priori bounds on the electrostatic potential and ionic concentrations. Finally, some numerical illustrations are presented in §5; a more detailed discussion of the electrochemical aspects of these results will be reported elsewhere.

Electrochemical model

Geometry and conservation equations

We consider a periodic setting with elementary cell [0, L * ] d , d ∈ {2, 3}, with length scale L * (expressed in m). The elementary cell contains inclusions Ω S whose boundary ∂Ω S contains negative charges with surface density Σ S (expressed in Cm -2 ) (see Figure 1 left). Our approach also applies to other settings, e.g., confined electrolytes in nanochannels (see Figure 1 right). The problem is posed in the domain Ω := [0, L * ] d \ Ω S and consists in finding the electrostatic potential ψ (expressed in V) and the ionic concentrations c = (c + , c -) (each concentration is expressed in m -3 ) such that where e is the elementary charge (expressed in C), ε = ε 0 ε r the solvent permittivity with ε 0 the vacuum permittivity (expressed in CV -1 m -1 ) and ε r the solvent relative permittivity (dimensionless), Z i the valence (dimensionless) of species i = ±, while the electrochemical potentials µ ± (ψ, c) are given by

-∆ψ = e ε i=± Z i c i in Ω, (1a)
µ + (ψ, c) and µ -(ψ, c) are constant in Ω, (1b)
µ i (ψ, c) = k B T log(σ 3 c i ) + k B T log(γ i (c)) + Z i eψ, i = ±, (2)
where k B is the Boltzmann constant (expressed in JK -1 ), T the temperature (expressed in K), σ the mean ionic diameter (expressed in m), and γ i (c) the activity coefficient (dimensionless) of species i accounting for non-ideal behavior (γ i (c) ≡ 1, i = ±, in the ideal case). Boundary conditions enforce that

ψ is periodic on ∂Ω \ ∂Ω S , (3a)
∇ψ • n = - 1 ε Σ S on ∂Ω S , (3b)
where n denotes the unit outward normal to ∂Ω. Additionally, we prescribe the mean ionic concentrations in the form

c i Ω = c 0 i , i = ±, (4)
where c 0 + and c 0 -are given positive real numbers and where, for any function f ∈ L 1 (Ω),

f Ω := 1 |Ω| Ω f. (5)
We assume that c 0 ± satisfy the bulk electro-neutrality condition

i=± Z i c 0 i = 1 |Ω| ∂Ω 1 e Σ S , (6)
which is a necessary and sufficient condition for the solvability of the Poisson equation (1a) together with the boundary conditions (3) for the electrostatic potential ψ given the ionic concentrations c. Finally, since the electrostatic potential is determined up to an additive constant, we require that ψ Ω = 0.

Modelling of non-ideality

For each species i = ±, its activity coefficient γ i (c), which depends on both ionic concentrations c ± , is split into two parts in such a way that

log(γ i (c)) = log(γ Coul i (c)) + log(γ HS (c)), i = ±. (7)
The first contribution accounts for Coulomb interactions, while the second is a hardsphere term introducing steric effects which dominate at high ionic concentrations.

For the Coulomb term, we consider the mean spherical approximation (MSA) hinging on the screening parameter Γ MSA (expressed in m -1 ) such that

Γ 2 MSA = πL B i=± Z 2 i c i (1 + σ i Γ MSA ) 2 , ( 8
)
where L B is the Bjerrum length defined by L B := e 2 4πεk B T

. Recalling that, in the present work, we consider the mean ion diameter approximation σ ± = σ, equation ( 8) can be solved explicitly. For convenience, we introduce the ionic strength function

I : R 2 ≥0 → R (expressed in m -3 ) such that, for all c ∈ R 2 ≥0 , I(c) := i=± η i c i , η i := 1 2 Z 2 i . (9)
Then, the screening parameter Γ MSA solving ( 8) is given by Γ MSA = Υ MSA (I(c)) with the function Υ MSA : R ≥0 → R ≥0 such that, for all θ ∈ R ≥0 ,

Υ MSA (θ) := 1 2σ 2σ(4πL B ) 1/2 (2θ) 1/2 + 1 -1 . (10)
The activity coefficients γ Coul i (c), i = ±, depend on c only through the ionic strength I(c) and are given by

log(γ Coul i (c)) = -Z 2 i L B Υ MSA (I(c)) 1 + σΥ MSA (I(c)) , i = ±. (11)
It is convenient to introduce the function

γ 0 : R ≥0 → R ≥0 such that, for all θ ∈ R ≥0 , log(γ 0 (θ)) := - L B Υ MSA (θ) 1 + σΥ MSA (θ) , (12)
so that log(γ Coul i (c)) = Z 2 i log(γ 0 (I(c))), i = ±. ( 13
)
The hard-sphere contribution stems from the lowest order term in the Percus-Yervick approximation, that is,

log(γ HS (c)) := 4πσ 3 3 (c + + c -). ( 14
)
Such expressions have been derived, e.g., by Waisman and Lebowitz [START_REF] Lebowitz | Mean spherical model integral equations for charged hard spheres. I. Method of solution[END_REF][START_REF] Lebowitz | Mean spherical model integral equations for charged hard spheres. II. Results[END_REF] and in the form presented above by Blum [START_REF] Blum | Mean spherical model for asymmetric electrolytes[END_REF] who treated the case of asymmetric electrolytes for the unrestricted primitive model.

Remark 1 (Debye-Hückel limit) For extremely low values of the mean ionic diameter, σ → 0, we recover for log γ 0 the expression derived in the Debye-Hückel theory, namely log γ 0

(θ) = -(2πL 3 B θ) 1/2 .

Scaling and non-dimensionalization

There are three length scales in the problem: the characteristic size of the elementary cell L * , the mean ionic diameter σ, and the Bjerrum length L B . It is convenient to introduce the Debye length L D and the non-dimensional ratio λ such that

L D := L 3 * 4πL B , λ := L D L * 2 = L * 4πL B . ( 15
)
The Debye length represents the scale over which the ions screen out the electric field.

In what follows, we use L * as the characteristic length. Moreover, the reference electrostatic potential is the so-called Zeta potential ψ * := k B T /e, the reference electrochemical potential is µ * := k B T , the reference surface charge density is Σ S * := k B T ε/(eL * ), and the reference ionic concentration is c * := L -3 * . Typical values for these quantities are provided in §5. With these reference values, the governing equations are recast into non-dimensional form, and to alleviate the notation, we use the same symbols for non-dimensional quantities.

The Poisson problem for the electrostatic potential ψ takes the form

-λ∆ψ = i=± Z i c i in Ω, (16a)
ψ is periodic on ∂Ω \ ∂Ω S , (16b)
∇ψ • n = -Σ S on ∂Ω S , (16c)
ψ Ω = 0. ( 16d
)
The mean ionic concentrations c 0 ± satisfy the bulk electro-neutrality condition

i=± Z i c 0 i = λ |Ω| ∂Ω Σ S , (17)
and the electrochemical potentials are such that

µ i (ψ, c) = log(σ 3 c i ) + log(γ i (c)) + Z i ψ, i = ±. (18)
The activity coefficients γ ± (c) are decomposed as in [START_REF] Ekeland | Analyse convexe et problèmes variationnels[END_REF]. The Coulomb term is still given by ( 13) with the ionic strength evaluated as in ( 9) and the function

γ 0 : R ≥0 → R ≥0 now such that, for all θ ∈ R ≥0 , log(γ 0 (θ)) := - 1 4πλ Υ MSA (θ) 1 + σΥ MSA (θ) , (19)
with the function Υ MSA : R ≥0 → R ≥0 such that, for all θ ∈ R ≥0 ,

Υ MSA (θ) := 1 2σ 2σλ -1/2 (2θ) 1/2 + 1 -1 . (20)
Finally, the hard-sphere term is still given by ( 14).

Entropy and Free-Energy functionals

We define the ideal entropy s id : R ≥0 → R such that, for all u ∈ R ≥0 ,

s id (u) := u(log(σ 3 u) -1), u > 0, 0, u = 0, (21)
and the excess entropy s ex : R 2 ≥0 → R (resulting from non-ideality) such that, for all c = (c

+ , c -) ∈ R 2 ≥0 , s ex (c) := 2Γ 0 (I(c)) + 2πσ 3 3 (c + + c -) 2 , (22)
where, for all θ ∈ R ≥0 ,

Γ 0 (θ) := θ 0 log(γ 0 (θ ′ ))dθ ′ = - 1 4πσλ θ - 8λσ 3 (Υ MSA (θ)) 3 -2λ (Υ MSA (θ)) 2 . ( 23
)
The ideal entropy is continuous in R ≥0 and continuously differentiable in R >0 , while the excess entropy is continuously differentiable in R 2 ≥0 with

∂s ex ∂c + (c) = log(γ + (c)), ∂s ex ∂c - (c) = log(γ -(c)). (24)
The entropy functional S is such that

S(c) = S id (c) + S ex (c), (25)
with

S id (c) := i=± Ω s id (c i ), S ex (c) := Ω s ex (c). ( 26
)
Finally, the free-energy functional E is such that

E(ψ, c) = U(ψ) -S(c) -B(ψ, c), (27)
where

B(ψ, c) = i=± Ω Z i c i ψ, (28a)
U(ψ) = λ 2 Ω |∇ψ| 2 + λ ∂Ω S Σ S ψ. ( 28b
)
To motivate the free-energy functional E, we observe formally that, if (ψ, c) is a saddle point of E under the constraints ψ Ω = 0 and c i Ω = c 0 i , i = ±, and if E is differentiable at (ψ, c), then (ψ, c) is a critical point of E. Differentiating E with respect to its first argument shows that ψ solves the Poisson problem [START_REF] Li | Minimization of electrostratic free energy and the Poisson-Boltzmann equation for molecular solvation with implicit solvent[END_REF], while differentiating E with respect to its second argument shows that the electrochemical potentials µ ± (ψ, c) are constant in Ω with µ i (ψ, c) = µ 0 i , i = ±, where the µ 0 i 's are the Lagrange multipliers associated with the constraints on the mean ionic concentrations. The purpose of the next section is to give a precise mathematical setting for this result.

Mathematical analysis

In this section, we introduce a mathematical framework and prove that, under the assumptions stated below, the free-energy functional E admits a unique saddle point (ψ, c). Moreover, we establish some a priori bounds on (ψ, c) and show that ψ solves the Poisson problem ( 16) while the electrochemical potentials µ i (ψ, c), i = ±, defined by [START_REF] Prohl | Convergent finite element discretizations of the Navier-Stokes-Nernst-Planck-Poisson system[END_REF] are constant in Ω. In what follows, we consider an abstract setting for the one-real variable function log γ 0 which can be chosen arbitrarily provided assumptions (H γ {1,2}) below are satisfied. The case of the MSA where log γ 0 is defined by ( 19) is a special case of application, for which the verification of assumptions (H γ {1,2}) is discussed in §5.

Assumptions

We consider the sets

H := φ ∈ H 1 per (Ω), φ Ω = 0 , (29)
K := c = (c + , c -) ∈ [L 2 (Ω)] 2 , c ± ≥ 0 a.e. in Ω, c ± Ω = c 0 ± , (30)
where the functional spaces H 1 per (Ω) and L 2 (Ω) are, respectively, the closure of C ∞ per (Ω) for the canonical norms

• H 1 (Ω) and • L 2 (Ω) . It is clear that H is a closed subspace of H 1 per (Ω) and that K is a closed convex subset of [L 2 (Ω)] 2 .
In what follows, we make the following assumptions:

• (H c 0 ) Σ S ∈ L 1 (∂Ω S ), the real numbers c 0 ± are positive and satisfy the bulk electroneutrality condition [START_REF] Looker | Semilinear elliptic neumann problems with rapid growth in the nonlinearity[END_REF].

• (H Ω ) Σ S ∈ H 1/2 (∂Ω S ) and the affine operator L :

L 2 (Ω) → H such that, for all f ∈ L 2 (Ω), L(f ) ∈ H solves -∆L(f ) = f -f Ω + |Ω| -1 ∂Ω S Σ S in Ω with the Neumann boundary condition ∇L(f )•n = -Σ S on ∂Ω S and L(f ) Ω = 0, is bounded from L 2 (Ω) to H 2 (Ω). • (H γ 1) The function θ → log(γ 0 (θ)) is continuous on R ≥0 ; moreover, there is β ∈ [0, 1) and (C 1 , C 2 ) ∈ R 2 ≥0 such that ∀θ ∈ R ≥0 , | log(γ 0 (θ))| ≤ C 1 + C 2 θ β . (31)
• (H γ 2) The function θ → log(γ 0 (θ)) is non-increasing and continuously differentiable on R >0 and there holds, for all θ ∈ R >0 ,

η ♯ θ + 4πσ 3 3 + 2η 2 ♯ + 2πσ 3 3η ♭ θ(η ♯ -η ♭ ) 2 (log γ 0 ) ′ (θ) > 0, ( 32
)
where η ♯ := max(η + , η -) and η ♭ := min(η + , η -).

An important consequence of assumption (H γ 1) is that the excess entropy s ex is a nondecreasing function of both its arguments if at least one ionic concentration is large enough.

Lemma 1 Assume (H γ 1). Then, there is κ γ ∈ R ≥0 such that, for all c ∈ R 2 ≥0 satisfying c + ≥ κ γ or c -≥ κ γ , there holds

∂s ex ∂c + (c) ≥ 0, ∂s ex ∂c - (c) ≥ 0. ( 33
)
Proof. Recall that

∂s ex ∂c ± (c) = log(γ ± (c)) = 2η ± log(γ 0 (I(c)) + 4π 3 σ 3 (c + + c -).
Invoking assumption (H γ 1), we infer that for suitable constants C 3 and C 4 , there holds

∂s ex ∂c ± (c) ≥ 4π 3 σ 3 (c + + c -) -C 3 -C 4 (c β + + c β -).
Since β ∈ [0, 1), the conclusion is straightforward. ♦

Main result

It is readily verified that the functional E maps H × K to R. We say that (ψ, c)

∈ H × K is a saddle point of E if ∀c ∈ K, E(ψ, c) ≤ E(ψ, c) ≤ E(φ, c), ∀φ ∈ H. (34)
We can now state the main result of this work.

Theorem 1 Assume (H c 0 ), (H Ω ), and (H γ {1,2}). Then, the functional E has a unique saddle point (ψ, c) ∈ H×K. Moreover, ψ ∈ L ∞ (Ω) and there are 0 < c ♭ ≤ c ♯ < +∞ such that, for a.e. x ∈ Ω, c ♭ ≤ c i (x) ≤ c ♯ , i = ±. Finally, ψ solves the Poisson problem [START_REF] Li | Minimization of electrostratic free energy and the Poisson-Boltzmann equation for molecular solvation with implicit solvent[END_REF], and the electrochemical potentials µ i (ψ, c), i = ±, defined by [START_REF] Prohl | Convergent finite element discretizations of the Navier-Stokes-Nernst-Planck-Poisson system[END_REF] are constant in Ω.

Proof. The proof is decomposed into several steps.

Existence of a saddle point For any c ∈ K, the functional H ∋ ψ → E(ψ, c) ∈ R is strictly convex, continuous, and satisfies, for all ψ ∈ H with ψ H 1 (Ω) → +∞, E(ψ, c) → +∞). Furthermore, in Lemma 2 below, we use assumption (H γ 2) to prove that the functional S is convex on K, while in Lemma 3, we use assumption (H γ 1) to prove that the functional S is continuous on K. Hence, for any ψ ∈ H, the functional Characterization and bound on electrostatic potential Let (ψ, c) ∈ H × K be a saddle point of the functional E. Since E is differentiable with respect to ψ and since H is a vector space, there holds

K ∋ c → E( ψ, c) ∈ R
∂ ψ E(ψ, c), φ = λ Ω ∇ψ•∇φ + λ ∂Ω S Σ S φ - i=± Ω Z i c i φ = 0, ∀φ ∈ H.
This shows that ψ solves the Poisson problem [START_REF] Li | Minimization of electrostratic free energy and the Poisson-Boltzmann equation for molecular solvation with implicit solvent[END_REF]. Moreover, recalling the affine operator L : L 2 (Ω) → H introduced in assumption (H Ω ) and using the bulk electroneutrality condition [START_REF] Looker | Semilinear elliptic neumann problems with rapid growth in the nonlinearity[END_REF], we infer

ψ = L 1 λ i=± Z i c i .
As a result, ψ ∈ H 2 (Ω), and owing to the Sobolev embedding theorem, ψ ∈ L ∞ (Ω).

Characterization and bound on ionic concentrations Using assumption (H γ 1), we prove in Lemmas 4 and 5 below that there are 0 < c ♭ ≤ c ♯ < +∞ such that, for a.e. x ∈ Ω, c ♭ ≤ c i (x) ≤ c ♯ , i = ±. Finally, owing again to the uniform lower bound on the ionic concentrations and using Lemma 3, we infer that the functional E is Gâteauxdifferentiable at (ψ, c) with respect to c ± along any direction v ∈ C ∞ per (Ω) with v Ω = 0, and there holds

∂ c ± E(ψ, c), v = Ω µ ± (ψ, c)v,
where µ ± (ψ, c) are the electrochemical potentials defined by [START_REF] Prohl | Convergent finite element discretizations of the Navier-Stokes-Nernst-Planck-Poisson system[END_REF]. Let now v ∈ C ∞ per (Ω). Since (ψ, c) is a saddle point, there holds ∂ c ± E(ψ, c), ṽ = 0 with ṽ = vv Ω , whence

Ω {µ ± (ψ, c) -µ ± (ψ, c) Ω }v = 0. By density of C ∞ per (Ω) in L 2 (Ω)
, this shows that the electrochemical potentials are constant in Ω.

Uniqueness

The functional E is strictly convex in its first argument. Moreover, owing to the uniform lower bound on the ionic concentrations, entropy is strictly convex in c (see, again, Lemma 2 below). This yields uniqueness of the saddle point. ♦ Remark 2 (Critical points are saddle points) It is readily verified using convexity arguments that if (ψ, c) (with c uniformly bounded from above and below) is a critical point of the free-energy functional E, then (ψ, c) is a saddle point of E. In §5, we will compute an approximation of the saddle point by solving the Euler-Lagrange equations satisfied by the critical point.

Technical results

In this section, we state and prove the various lemmas invoked in the proof of our main result, Theorem 1. Our first result concerns the convexity of the entropy functional S.

Lemma 2 (Convexity of entropy) Assume (H γ 2). Then, the functional S is convex on K and strictly convex on the subset K >0 := {c ∈ K; c + > 0 and c -> 0 a.e. in Ω}.

Proof.

It suffices to show that the entropy s : R 2 ≥0 → R such that, for all c = (c

+ , c -) ∈ R 2 ≥0 , s(c) := s id (c + ) + s id (c -) + s ex (c)
is convex on R 2 ≥0 and strictly convex on R 2 >0 . Moreover, the convexity of s on R 2 ≥0 follows from the strict convexity of s on R 2 >0 and the continuity of s. Hence, it is sufficient to address the strict convexity of s on R 2 >0 . Set η ♯ := max(η + , η -), η ♭ := min(η + , η -), and η := 1 2 (η + + η -). The Hessian associated with s is given by Hess(s) =

1 c + + A + Bη 2 + A + Bη + η - A + Bη + η - 1 c -+ A + Bη 2 - ,
where A := 4π 3 σ 3 and B := 2(log γ 0 ) ′ (θ) with θ = I(c). We verify that the trace and the determinant of H(s) are positive. We obtain

c + c -det Hess(s) = 1 + A(c + + c -) + η 2 + c + + η 2 -c -+ Ac + c -(η + -η -) 2 B. Since c + + c -≥ θ/η ♯ , η 2 + c + + η 2 -c -≤ η ♯ θ, c + c -≤ θ 2 /(4η -η + )
, and B ≤ 0, we infer that, under condition (32), that is,

η ♯ θ + A + η 2 ♯ + A 4η ♭ θ(η ♯ -η ♭ ) 2 B > 0, ( 35
)
there holds det Hess(s) > 0. Furthermore, 1 2 tr Hess(s) = 1 2

1 c + + 1 c - + A + 1 2 (η 2 + + η 2 -)B.
Since (c -1 + + c -1 -)/2 ≥ η/θ, we infer that, under the condition

η θ + A + 1 2 (η 2 ♯ + η 2 ♭ )B > 0, ( 36
)
there holds tr Hess(s) > 0. Finally, it is readily verified that, for all θ ∈ R >0 , (

η ♯ θ + A) 1 2 (η 2 + + η 2 -) ≤ ( η θ + A)(η 2 ♯ + A 4η ♭ θ(η ♯ -η ♭ ) 2
), so that (35) implies (36). ♦ Remark 3 (Necessary convexity condition for symmetric electrolytes) In the case of symmetric electrolytes, that is, η + = η -= η, condition (32) reduces to

η θ + 4πσ 3 3 + 2η 2 (log γ 0 ) ′ (θ) > 0, ( 37
)
and it is readily verified that this condition is also necessary for the convexity of the entropy.

Our second result deals with the continuity and Gâteaux-differentiability of the entropy functional S.

Lemma 3 (Continuity and differentiability of entropy) Assume (H γ 1). Then, the entropy functional S is continuous on [L 2 (Ω)] 2 . Moreover, for all c ∈ [L 2 (Ω)] 2 such that there is c ♭ > 0 with c i (x) ≥ c ♭ for a.e. x ∈ Ω and all i = ±, S is Gâteauxdifferentiable at c along any direction v ∈ C ∞ per (Ω), and there holds

∂ c ± S(c), v = Ω {log(σ 3 c ± ) + log(γ ± (c))}v. Proof. A
∃a, b ∈ R, ∀u ∈ R, |F (u)| ≤ a + b|u| p/q ,
with 1 ≤ p, q < +∞ then, for all v ∈ L p (Ω), there holds F (v) ∈ L q (Ω), and the superposition operator L p (Ω) ∋ v → F (v) ∈ L q (Ω) is continuous. This result can be applied to the function

F 1 (x) = s id (x) (extended by zero for x ≤ 0) showing that S id (c) is continuous from [L 2 (Ω)] 2 to L 1 (Ω).
The result can also be applied to the function F 2 (x) = Γ 0 (x) (extended by zero for x ≤ 0) since, owing to assumption (H γ 1), F 2 satisfies the above growth condition with exponent 1 + β < 2. Since the function c → I(c) maps continuously [L 2 (Ω)] 2 to L 2 (Ω), we infer that the functional c → 2Γ 0 (I(c)) is continuous from [L 2 (Ω)] 2 to L 1 (Ω). Finally, the continuity of the hard-sphere contribution c

→ 2πσ 3 3 (c + + c -) 2 from [L 2 (Ω)] 2 to L 1 (Ω) is obvious.
Gâteaux-differentiability Let now c ∈ [L 2 (Ω)] 2 such that there is c ♭ > 0 with c i (x) ≥ c ♭ for a.e. x ∈ Ω and all i = ±. We treat the partial derivative with respect to c + ; the other case is treated similarly. Let v ∈ C ∞ per (Ω) with v Ω = 0. There is t 0 > 0 such that, for all t ∈ [-t 0 , t 0 ], c + (x) + tv(x) ≥ 1 2 c ♭ for a.e. x ∈ Ω. Since the function F 3 (x) = log(σ 3 x) for x ≥ 1 2 c ♭ extended by the constant value log( 1 2 σ 3 c ♭ ) for x ≤ 1 2 c ♭ satisfies a linear growth condition. Using the proof of [START_REF] Kavian | Introduction à la théorie des points critiques[END_REF]Lemma 17.1,p. 64], we can prove that the ideal entropy is Gâteaux-differentiable with

∂ c + S id (c), v = Ω log(σ 3 c + )v.
Gâteaux-differentiating the excess entropy is simpler. We use the growth condition in assumption (H γ 1) to infer

∂ c + Ω 2Γ 0 (I(c)) , v = Ω 2η + log(γ 0 (I(c))v = Ω log(γ Coul + (c))v,
while Gâteaux-differentiating the hard-sphere term 2πσ 3 3 (c + + c -) 2 is straightforward. ♦ Our third result delivers an a priori L ∞ -bound on the ionic concentrations c = (c + , c -).

Lemma 4 (Upper bound on c) Assume (H γ 1). Let (ψ, c) ∈ H × K be a saddle point of the functional E. Then, there is c ♯ < +∞ such that, for a.e. x ∈ Ω, c i (x) ≤ c ♯ for all i = ±.

Proof. Let (ψ, c) ∈ H × K be a saddle point of the functional E. Proceeding by contradiction, we assume that there is i = ± such that, for all n ∈ N, the set

A n i := {x ∈ Ω; c i (x) > 2 n
} has positive measure. We then construct modified ionic concentrations c ∈ K such that E(ψ, c) > E(ψ, c), thereby providing the desired contradiction with (34). As a result, for all i = ±, there is n i ∈ N such that the set A n i i has zero measure, yielding the statement of Lemma 4 with c ♯ = max i=± 2 n i . The principle of the construction is that, by diminishing the ionic concentration where it is very large, the entropy can be decreased, and thus the free-energy increased. Both the ideal term and the steric exclusion term are large enough at high concentrations to lead to an entropy decrease. We choose to work with the ideal term since, in the proof of Lemma 5 below, the ideal term is the only one leading to the entropy decrease at small concentrations; thus, the two proofs are similar. Furthermore, we observe that some care is needed when perturbing the ionic concentrations since it is necessary to preserve their mean values. Without loss of generality, we assume that, for all n ∈ N, the set A n + has positive measure. We first observe that there is k ∈ N such that the set

Ω k + := {x ∈ Ω; 2 -k ≤ c + (x) ≤ 2 k } (38)
has positive measure (otherwise, c + is zero or infinity a.e. in Ω which contradicts the fact that c + Ω = c 0 + > 0). In what follows, we fix such k ∈ N, and, without loss of generality, we assume that n ≥ k so that the sets A n + and Ω k + are disjoint. Moreover, since c -Ω = c 0 -and c -(x) ≥ 0 for a.e. x ∈ Ω, we infer that, for all m ∈ N, c 0

-≥ 1 |Ω| 2 m |A m -|, which shows that |A m -| → 0 as m → ∞. As a result, there is m ∈ N such that the set Ω k,m := Ω k + ∩ (Ω \ A m -) has positive measure.
In what follows, we fix such m ∈ N. We observe that in Ω k,m , both ionic concentrations c ± are bounded by C k,m := max(2 k , 2 m ). We now define the function cn + as follows:

cn + (x) =        0, x ∈ A n + , c + (x) + δ n , x ∈ Ω k,m , c + (x), x ∈ Ω \ (A n + ∪ Ω k,m ). with δ n = 1 |Ω k,m | A n + c + .

It is readily verified that cn

+ Ω = c 0 + so that cn := (c n + , c -) ∈ K. We observe that the real number δ n is uniformly bounded since

δ n ≤ δ := 1 |Ω k,m | |Ω|c 0 + for all n ∈ N.
It is important to modify c + only in the set Ω k,m to preserve its mean value (and not in the larger set Ω k + ); indeed, a bound on c -is needed to control the variation of the non-ideal terms between c and cn . To conclude the proof, we show that it is possible to choose n large enough so that ∆E := E(ψ, cn ) -E(ψ, c) > 0.

There holds ∆E = -∆S id -∆S ex -∆B with ∆S id := S id (c n ) -S id (c), ∆S ex := S ex (c n ) -S ex (c), and ∆B := B(ψ, cn ) -B(ψ, c).

We estimate the three terms separately. Since

∆B = -Z + A n + c + ψ + Z + Ω k,m δ n ψ,
we infer, since ψ ∈ L ∞ (Ω) owing to the second step in the proof of Theorem 1, that

|∆B| ≤ 2Z + ψ L ∞ (Ω) A n + c + . Since ∆S id = A n + -s id (c + ) + Ω k,m {s id (c + + δ n ) -s id (c + )} =: T 1 + T 2 ,
we infer that

T 1 ≤ -(log(σ 3 2 n ) -1) A n + c + , |T 2 | ≤ δ n Ω k,m m(c + , c + + δ n ) ≤ m(2 -k , 2 k + δ) A n + c + ,
where we have used the fact that c + ≥ 0, δ n ≤ δ, and that, for b ∈ R >0 and a ∈ R Turning next to ∆S ex , recall that the excess entropy s ex is continuously differentiable in R 2 ≥0 and that ∂sex ∂c + (u) ≥ 0 for all u = (u + , u -) ∈ R 2 ≥0 such that u + ≥ κ γ or u -≥ κ γ as shown in Lemma 1. Let C γ := max u∈Kγ | ∂sex ∂c + (u)| with the compact set

K γ := [0, max(2 k + δ, κ γ )] × [0, max(2 m , κ γ )]. We decompose ∆S ex into ∆S ex = A n + {s ex (c n + , c -) -s ex (c + , c -)} + Ω k,m {s ex (c n + , c -) -s ex (c + , c -)} =: T 3 + T 4 .
Observing that

s ex (c n + , c -) -s ex (c + , c -) = cn + c + ∂sex ∂c + (u + , c -)du + , we obtain |T 4 | ≤ Ω k,m C γ δ n = C γ A n + c + ,
since for all x ∈ Ω k,m and for all u + ∈ c + (x), cn + (x) , (u + , c -(x)) ∈ K γ . Moreover, owing to Lemma 1,

T 3 ≤ A n + {s ex (0, c -) -s ex (min(κ γ , c + ), c -)} ≤ A n + ∩{c -≤κγ } {s ex (0, c -) -s ex (min(κ γ , c + ), c -)} =: T ′ 3 ,
since for c -> κ γ , s ex (0, c -)s ex (min(κ γ , c + ), c -) ≤ 0. Moreover,

|T ′ 3 | ≤ A n + ∩{c -≤κγ } min(κγ ,c + ) 0 ∂s ex ∂c + (u + , c -) du + ≤ A n + ∩{c -≤κγ } C γ min(κ γ , c + ) ≤ C γ A n + c + .
Collecting the above bounds, we infer

∆E ≥ (log(σ 3 2 n ) -1 -C) A n + c + , with C = 2Z + ψ L ∞ (Ω) + m(2 -k , 2 k + δ) + 2C γ .
Taking n large enough so that log(σ 3 2 n ) ≥ 1 + C and since A n + has positive measure, we infer ∆E > 0. ♦ Our last result delivers a uniform positive lower bound on the ionic concentrations.

Lemma 5 (Uniform positive lower bound on c) Assume (H γ 1). Let (ψ, c) ∈ H × K be a saddle point of the functional E. Then, there is c ♭ > 0 such that, for a.e. x ∈ Ω, c i (x) ≥ c ♭ for all i = ±.

Proof. Let (ψ, c) ∈ H × K be a saddle point of the functional E. The structure of the proof is similar to that of Lemma 4, though a bit simpler since we already have upper bounds on c ± . Proceeding by contradiction, we assume that there is i = ± such that, for all n ∈ N, the set B n i := {x ∈ Ω; c i (x) < 2 -n } has positive measure. We then construct modified ionic concentrations c ∈ K such that E(ψ, c) > E(ψ, c), thereby providing the desired contradiction with (34). As a result, for all i = ±, there is n i ∈ N such that the set B n i i has zero measure, yielding the statement of Lemma 5 with c ♭ = min i=± 2 -n i . Without loss of generality, we assume that, for all n ∈ N, the set B n + has positive measure. We fix k ∈ N such that the set Ω k + defined by (38) has positive measure, and, without loss of generality, we assume that n is large enough so that 2

-n ≤ 1 |Ω| |Ω k + |2 -k . Since this implies 2 -n ≤ 2 -k , the sets B n
+ and Ω k + are disjoint. We now define the function cn + as follows:

cn + (x) =        c + (x) + 2 -n , x ∈ B n + , c + (x) -δ n , x ∈ Ω k + , c + (x), x ∈ Ω \ (B n + ∪ Ω k + ). with δ n = 1 |Ω k + | |B n + |2 -n . It is readily verified that cn + Ω = c 0 + . Moreover, since |B n + | < |Ω|, δ n < 1 |Ω k + | |Ω|2 -n ≤ 2 -k
, so that cn + ≥ 0 in Ω. Hence, cn := (c n + , c -) ∈ K. To conclude the proof, we show that it is possible to choose n large enough so that ∆E := E(ψ, cn ) -E(ψ, c) > 0.

As in the proof of Lemma 4, we write ∆E = -∆S id -∆S ex -∆B. Since ψ ∈ L ∞ (Ω) owing to the second step in the proof of Theorem 1, we infer

|∆B| ≤ 2Z + ψ L ∞ (Ω) 2 -n |B n + |.
Moreover,

∆S id = B n + {s id (c + + 2 -n ) -s id (c + )} + Ω k + {s id (c + -δ n ) -s id (c + )} =: T 1 + T 2 ,
with

T 1 ≤ B n + 2 -n log(σ 3 2 -n+1 ) = log(σ 3 2 -n+1 )2 -n |B n + |, |T 2 | ≤ Ω k + δ n m(c + -δ n , c + ) ≤ m(2 -k -δ, 2 k )2 -n |B n + |. Finally, let C ′ γ := max u∈K ′ γ | ∂sex ∂c + (u)| with the compact set K ′ γ := [0, 2 k ] × [0, c ♯ ]. We obtain |∆S ex | ≤ B n + ∪Ω k + cn + c + ∂s ex ∂c + (u + , c -) du + ≤ C ′ γ B n + ∪Ω k + |c n + -c + | ≤ 2C ′ γ 2 -n |B n + |, since for all x ∈ B n + ∪ Ω k + and for all u + ∈ c + (x), cn + (x) (u + , c -(x)) ∈ K ′ γ . Collecting the above bounds, we infer ∆E ≥ (log(σ -3 2 n-1 ) -C)2 -n |B n + |, with C = 2Z + ψ L ∞ (Ω) + m(2 -k -δ, 2 k ) + 2C ′ γ .
Taking n large enough so that log(σ -3 2 n-1 ) ≥ C and since B n + has positive measure, we infer ∆E > 0. ♦

Numerical illustration

This section presents a numerical experiment in the case where the activity coefficient log γ 0 is evaluated using the MSA.

Verifications of assumptions for MSA

We verify the abstract assumptions (H γ {1,2}) in the context of the MSA, that is, when log γ 0 is defined by (19) with the screening parameter Υ MSA defined by (20). Assumption (H γ 1) is straightforward to verify, so that we focus on (H γ 2). For all θ > 0, setting y := 2σλ -1/2 (2θ) 1/2 + 1, we obtain

(log γ 0 ) ′ (θ) = - σ πλ 2 1 y(y + 1) 2 (y 2 -1)
, and a simple calculation shows that (H γ 2) is equivalent to the fact that the sixth degree polynomial

P (y) = y (y -1) 2 (y + 1) 3 - (η ♯ -η ♭ ) 2 16πλση ♭ y 2 -1 2 (y -1) + 6η ♯ πλσ y (y + 1) - 3 2 η ♯ πλσ 2 (y -1)
takes positive values for all y > 1. This condition, in turn, can be checked numerically.

It holds true under a condition of the form σλ > υ 0 , where the threshold υ 0 , which depends on the species valences Z ± , is reported in the second column of Table 1 for several values of Z ± . Recalling the definition (15) of the non-dimensional parameter λ and reverting to dimensional length scales, the above condition can be expressed as σ/L B > 4πυ 0 . Using the values e = 1.60 × 10 -19 C, ε 0 = 8.85 × 10 -12 CV -1 m -1 , ε r = 78.3, and k B = 1.38 × 10 -23 J K -1 , the Bjerrum length can be evaluated as a function of the temperature T (yielding, e.g., L B = 7.1

• A for T = 300 K) and a minimal value σ 0 for the mean ionic diameter can be computed. This value is reported in the the third and fourth columns of Table 1 for T = 300 and T = 350 K respectively. The condition σ > σ 0 is easier to fulfill when the temperature increases. This condition is also more stringent for 1:2 and 2:2 electrolytes than for 1:1 electrolytes. Interestingly, the condition σ > σ 0 shows that the mean ionic diameter cannot take extremely low values within the present physical model (in particular, the Debye-Hückel limit σ → 0 does not yield here a convex entropy).

Z + : Z - υ 0 σ 0 ( • A) - - T = 300 K T =

A periodic medium with charged disks

Referring to §2, we consider a two-dimensional setting where Ω S is a disk of radius R = 0.3L * whose center coincides with that of the elementary cell [0, L * ] 2 . We take L * ∈ {1, 10} nm, T = 300 K, and Σ S = 0.13 Cm -2 . With these values, the Debye length is in the range [3.345, 105.8] • A. We consider a 1:1 electrolyte. The mean ionic concentrations (in dimensional form) are specified as

c 0 + = 1 |Ω| ∂Ω 1 e Σ S + c salt , c 0 -= c salt ,
and we set the concentration of added salt to c salt = 0.15 mol/l. In particular, since the disk is negatively charged, there is always an excess of counter-ions to ensure the global electro-neutrality of the system. It is readily seen that assumptions (H c 0 ) and (H Ω ) hold true. Moreover, the mean ionic diameter is set to σ ∈ {3, 4, 5}

• A, so that assumptions (H γ {1,2}) also hold true. Figure 2 (left) depicts the activity coefficient log γ 0 as a function of ionic strength for the various values of the parameter σ, whereas Figure 2 (right) depicts the hard-sphere activity coefficient log γ HS as a function of total concentration. This figure illustrates that electrostatic correlations have a more sizable effect for small values of σ, whereas the opposite effect is observed for steric exclusion effects. The saddle point of the free-energy functional is computed numerically by solving the associated Euler-Lagrange conditions approximately using conforming finite elements for space discretization (using the FreeFEM++ library [START_REF] Hecht | FreeFem++ documentation[END_REF]) in conjunction with a Newton-Raphson algorithm to solve the discrete set of coupled nonlinear equations. The constraints on the mean value of the electrostatic potential and of the ionic concentrations are conveniently handled using three additional Lagrange multipliers. More details concerning the numerical approach can be found in [11]. Figure 3 depicts iso-values of the counter-ion concentration for the two values of the reference length L * and for a mean ionic diameter equal to 3

• A (the elementary cell is rescaled to [0, 1] 2 in the figure). We observe that for large L * , boundary layers appear near the charged walls: counter-ion concentrations exhibit a steeper gradient close to the charged surface and take almost constant values in the bulk region far from the disk.

To gain further insight, we compare the solutions obtained for the three values of the parameter σ ∈ {3, 4, 5}

Figure 1 .

 1 Figure 1. Geometries for Ω: periodic media with negatively charged inclusions (left); nanochannel with negatively charged walls (right)

  ≥0 , there holds |s id (b)s id (a)| ≤ |b -a|m(a, b) with m(a, b) := max(| log(σ 3 a)|, | log(σ 3 b)|).

Figure 2 .

 2 Figure 2. Left: activity coefficient log γ 0 as a function of ionic strength I(c) (mol/l); Right: hard-sphere activity coefficient log γ HS as a function of total concentration (c + + c -) (mol/l)

•

  A and for the two values of L * ∈ {1, 10} nm. We focus on the values obtained on the horizontal line [0.8L * , L * ] × {0.5L * } joining the rightmost part of the disk to the right vertical side of the elementary cell.

  Figure 4 depicts the values of the activity coefficient log(γ ± (c)) (note that log(γ + (c)) = log(γ -(c)) for a symmetric electrolyte), while Figure5depicts those of the counter-ion concentration;

  classical result of nonlinear analysis [12, Lemma 16.2, p. 61] states that, if F : R → R is a continuous function satisfying the growth condition

Table 1 .

 1 Threshold values for which assumption (H γ 2) holds true.

	350 K

Acknowledgments

Partial support by ANDRA (French National Agency for the Management of Nuclear Waste) is gratefully acknowledged; in particular, the second author is supported through a PhD fellowship by ANDRA. The authors are thankful to O. Bernard,

for completeness, values obtained in the ideal Poisson-Boltzmann case (γ ± (c) = 1) are also reported. We observe three different behaviors in Figure 4: electrostatic correlations dominate for σ = 3

• A (log γ ± (c) > 0), steric exclusion effects dominate for σ = 5

• A (log γ ± (c) < 0), or both effects play a role for σ = 4

• A. In Figure 5, we can observe the influence of the non-ideality on the counter-ion concentration close to the charged disk, especially when comparing the concentrations to those obtained within the Poisson-Boltzmann theory. The main effect of non-ideality is to lower the counter-ion concentration close to the charged surface. Interestingly, the ideal predictions are more accurate for larger cell sizes (L * = 10 nm).


Roland Tóth
email: r.toth@tudelft.nl

Vincent Laurain
email: vincent.laurain@cran.uhp-nancy.fr.

Marion Gilson
email: marion.gilson@cran.uhp-nancy.fr.

Hugues Garnier
email: hugues.garnier@cran.uhp-nancy.fr.

On the closed loop identification of LPV models using instrumental variables

Keywords: LPV models, system identification, Refined Instrumental Variable, Box-Jenkins models, Input/Output, Transfer function

Identification of real life systems is often applied in closed-loop due to stability, performance or safety constraints. However, in case of Linear Parameter-Varying (LPV) systems, closed-loop identification is not well-established despite the recent advances in prediction-error identification. Building on the available results, the paper proposes the closed-loop generalization of a recently introduced instrumental variable scheme for the identification of LPV-IO models with Box-Jenkins type of noise models. Estimation under closed-loop conditions is analyzed from the stochastic point of view with the proposed approach and the performance of the method is demonstrated through a representative simulation example.

INTRODUCTION

Identification of physical or chemical systems is often restricted to data gathered during closed-loop operation due to stability, performance/economical or safety constraints. The fact that the control loop itself introduces correlation between the disturbances and the control signal yields that the statistically optimal estimation of the parameters of a chosen model structure w.r.t. the data-generating system is an essentially different problem than in the open-loop setting. Hence in the identification literature, especially in the linear time-invariant (LTI) case, many approaches have been proposed to provide well-applicable solutions in this problem setting. An overview of the developments in the LTI case can be found in [START_REF] Van Den Hof | Closed-loop issues in system identification[END_REF] and [START_REF] Forssell | Closed-loop identification revisited[END_REF].

Identification of linear parameter-varying (LPV) systems has recently seen significant improvements in terms of open-loop methods, however closed loop-identification has remained sparsely studied (see [START_REF] Casella | LPV/LFT modelling and identification: overview, synergies and a case study[END_REF]; [START_REF] Tóth | Modeling and Identification of Linear Parameter-Varying Systems[END_REF]). The main difficulty in this system class in terms of identification is that even if the signal relations between the inputs u and outputs y of the system are linear, these relations are allowed to change over time as a function of a so called scheduling variable p : Z → P with P ⊆ R np being the so called scheduling space. This allows to describe a large class of nonlinear/time-varying systems in an attractive structure allowing the use of well-established LPV control-synthesis approaches, e.g. [START_REF] Packard | Gain scheduling via linear fractional transformations[END_REF]; [START_REF] Scherer | Mixed H 2 /H ∞ control for time-varying and linear parametrically-varying systems[END_REF]. On the other hand, this varying-relation prevents the use of crucial system theoretical relations, like transfer functions and commutativity of operators, which are often used in LTI closed-loop identification approaches. As a consequence, only preliminary closed-loop approaches has been proposed in the literature without being able to exploit the existing tools and knowledge available in the LTI case. In van [START_REF] Van Wingerden | Subspace identification of bilinear and LPV systems for open-and closed-loop data[END_REF], an approximation based LPV extension of a predictor subspace approach (PSBID) has been proposed which is also applicable in a closed-loop setting, while in [START_REF] Boonto | Closed-loop system identification of LPV input-output models -application to an arm-driven pendulum[END_REF] also an approximation based LPV extension of the CLOE algorithm (see [START_REF] Landau | Recursive algorithms for identification in closed loop: A unified approach and evaluation[END_REF]) has been investigated w.r.t. LPV outputerror (OE) type of models. In [START_REF] Abbas | An instrumental variable technique for open-loop and closed-loop identification of input-output LPV models[END_REF] and [START_REF] Butcher | On the consistency of certain identification methods for linear parameter varying systems[END_REF] a basic instrumental variable method has been introduced. It is a clear indicator of the immature state of this research direction that the stochastic properties of the estimation has been analyzed only in [START_REF] Butcher | On the consistency of certain identification methods for linear parameter varying systems[END_REF]. Unfortunately in [START_REF] Laurain | Refined instrumental variable methods for identification of LPV Box-Jenkins models[END_REF], it has been shown recently that the formulation of the one-step IV approach proposed in [START_REF] Butcher | On the consistency of certain identification methods for linear parameter varying systems[END_REF] does not allow to reach statistically optimal estimates due no-commutativity of certain filtering operations. This highlights that currently no method has been established which allows stochastically optimal estimation under a closed-loop setting especially with general noise models like Box-Jenkins (BJ).

Based upon the available results of [START_REF] Abbas | An instrumental variable technique for open-loop and closed-loop identification of input-output LPV models[END_REF]; [START_REF] Butcher | On the consistency of certain identification methods for linear parameter varying systems[END_REF] and using the recent theoretical advancements of the LPV prediction error-framework introduced in Tóth (2010), we propose in this paper a closed-loop extension of a recently developed IV approach for LPV-BJ models, see [START_REF] Laurain | Refined instrumental variable methods for identification of LPV Box-Jenkins models[END_REF]. In this setting, the noise v o affecting the sampled output measurement y(t k ) is assumed to have a rational spectral density which is not restricted to depend on the actual process dynamics of the data-generating system (general noise structure). As a first step towards the case of a pdependent noise, it is also assumed that this rational spectral density does not depend on p. A stochastic analysis of the proposed closed-loop approach is provided, exploring the limitations and the benefits of this estimation strategy. The performance of the algorithm is also demonstrated on a representative simulation example.

PROBLEM DESCRIPTION

System description

Consider the data-generating LPV system S o given in a closed-loop setting of Fig. 1, where S o is defined in the following form:

A o (p k , q -1 )χ o (t k ) = B o (p k , q -1 )u(t k-d ), (1a)
y(t k ) = χ o (t k ) + v o (t k ), (1b)
Here u(t k ) is the input of the plant, p k is the value of the scheduling variable p at sample time t k , χ o is the noisefree output, v o is the additive noise with bounded spectral density, y is the noisy output of the system, d is the delay, and q is the time-shift operator, i.e. q -i u(t k ) = u(t k-i ).

A o (p k , q -1 ) and B o (p k , q -1 ) are polynomials in q -1 of degree n a and n b respectively:

A o (p k , q -1 ) = 1 + na i=1 a o i (p k )q -i , (2a)
B o (p k , q -1 ) = n b j=0 b o j (p k )q -j . ( 2b
)
Here the coefficients a i and b j are real meromorphic functions 1 with static dependence on p, i.e. dependence only on the instantaneous value of p at time t k . It is assumed that each a i and b j is non-singular on P, thus the solutions of S o are well-defined and the process part G o , defined by (1a), is completely characterized by the coefficient functions {a o i } na i=1 and {b o j } n b j=0 . For the sake of simplicity, the noise v o is assumed to be independent from p. In case of rational spectral density of v o , such a noise process is represented by a discrete-time auto-regressive moving average (ARMA) model:

v o (t k ) = H o (q)e o (t k ) = C o (q -1 ) D o (q -1 ) e o (t k ), (3)
where C o (q -1 ) and D o (q -1 ) = 0 are monic polynomials with constant coefficients and with respective degree n c and n d . The corresponding proper transfer function H o (q) is assumed to be stable and to have a stable inverse. In case C o (q -1 ) = D o (q -1 ) = 1, (3) defines an OE noise model, however with C o (q -1 ) = D o (q -1 ), ( 3) is general enough to represent BJ-type of noise models.

r 2 (t k ) G o C o r 1 (t k ) u(t k ) p k e o (t k ) v o (t k ) χ o (t k ) y(t k ) H o (q -1 ) Fig. 1. Closed-loop LPV system
In terms of the closed-loop setting of Fig. 1,

u(t k ) = r 1 (t k ) + C o (r 2 (t k ) -y(t k ))
, where C o is the operator form of the controller. The controller can be any LTI, nonlinear or LPV controller, under the assumptions that

• C o is a priori known;
• C o ensures BIBO stability of the closed-loop system S o for any p ∈ P Z .

Model considered

Next we introduce a model structure and parameterization for the identification of S o , where, according to datagenerating system, i.e. (1a-b) and ( 3), the noise model and the process model are parameterized separately. The proposed LPV-BJ model, denoted in the sequel as M θ , is defined as:

A(p k , q -1 , ρ)χ(t k ) = B(p k , q -1 , ρ)u(t k-d ), (4a)
D(q -1 , η)v(t k ) = C(q -1 , η)e(t k ), (4b)
y(t k ) = χ(t k ) + v(t k ), (4c)
with parameters θ = [ ρ η ] and with u(

t k ) = r 1 (t k ) + C o (r 2 (t k ) -y(t k )).
The process model part of M θ , denoted by G ρ , is defined in terms of the LPV-IO representation (4a) where A and B are polynomials with order n a and n b respectively and with p-dependent coefficients a i and b j parameterized as

a i (p k ) = a i,0 + nα l=1 a i,l f l (p k ) i = 1, . . . , n a (5a) b j (p k ) = b j,0 + n β l=1 b j,l g l (p k ) j = 0, . . . , n b (5b)
In this parametrization, {f l } nα l=1 and {g l } n β l=1 are a priori chosen meromorphic functions of p, with static dependence, allowing the identifiability of the model (linearly independent functions on P for example). The associated model parameters ρ are stacked column-wise:

ρ = [ a 1 . . . a na b 0 . . . b n b ] ∈ Θ ρ ⊂ R nρ , (6)
where

a i = [ a i,0 a i,1 . . . a i,nα ], b j = [ b j,0 b j,1 . . . b j,n β ], n ρ = n a (n α + 1) + (n b + 1)(n β + 1) andΘ ρ ⊂ R nρ is such that (4a) is BIBO stable for all ρ ∈ Θ ρ and p ∈ P Z .
The noise-model part of M θ , denoted by H η , is defined in terms of (4b) where C and D are monic polynomials with order n c and n d respectively and with constant coefficients c i and d j . These model parameters are stacked columnwise in the parameter vector,

η = [ c 1 . . . c nc d 1 . . . d n d ] ∈ Θ η ⊂ R nη , (7) where n η = n c + n d . As (4b) is LTI, it can be represented by a transfer function H(q,η ) = C(q -1 ,η) D(q -1 ,η) . It is further assumed that η ∈ Θ η ⊂ R nη such that H(z,η
) is stable (analytic in the exterior of the unit circle) and has a stable inverse.

Introduce also

G = {G ρ | ρ ∈ Θ ρ } and H = {H η | η ∈ Θ η },
as the collection of all process and noise models in the form of (4a) and (4b). Then the model set, denoted as M, takes the form

M = {(G ρ , H η ) | col(ρ,η ) = θ ∈ Θ ρ × Θ η } , (8)
corresponding to the set of candidate models in which we seek the model that explains data gathered from S o the best under a given identification criterion.

Identification problem statement

Denote D N = {y(t k ), u(t k ), p(t k )} N
k=1 a data sequence of S o . In [START_REF] Laurain | Refined instrumental variable methods for identification of LPV Box-Jenkins models[END_REF], it has been shown, that a onestep-ahead output predictor ŷθ (t k |t k-1 ) can be formulated w.r.t. the considered model structure (4a-c) under the commonly used assumption that noise-free observation of the sequence {p k , p k-1 , . . .} is available. Recently it has been proved that using estimated moments, such a predictor can be formulated if p is observed up to an additive white noise independent from v o , see [START_REF] Tóth | A prediction-error identification framework for linear parameter-varying systems[END_REF], however for the sake of simplicity we will only consider the previous case. These results allow to formulate (in either case) the estimation of θ, as the minimization of the one-step-ahead prediction error:

ε θ = y(t k ) - ŷθ (t k |t k-1
). This minimization is formulated in terms of an identification criterion W (D N , θ), like the least squares criterion

W (D N , θ) = 1 N N k=1 ε 2 θ (t k ) = 1 N ε θ (t k ) 2 2 , (9)
such that the parameter estimate is θN = arg min θ∈Θρ×Θη W (D N , θ).

(10)

Based on the previous considerations, the identification problem addressed in the sequel is defined as follows: Problem 1. Given a discrete-time LPV data-generating system S o in the closed-loop setting of Fig. 1 with a priori known stabilizing controller C o . Based on the LPV-BJ model structure M θ defined by (4a-c) and a data set D N collected from S o , estimate the parameter vector θ as the minimization of (9) under the following assumptions:

A1 S o ∈ M. A2 In the parametrization (5a-b) of A and B, {f l } nα l=1 and {g l } n β l=1 are chosen such that M θ is identifiable. A3 D N is informative w.r.t. M θ , i.e. if θ 1 = θ 2 then M θ1
and M θ2 lead to different prediction errors given D N . A4 S o is BIBO stable, i.e. for any bounded p ∈ P Z and u ∈ R Z , the output of S o is bounded Tóth (2010). A5 The closed-loop system with the controller C o and model M θ is BIBO stable for any θ ∈ Θ ρ × Θ η , i.e. the controller stabilizes S o and also any model in the considered model set. A6 p is treated as an external variable of the system.

In [START_REF] Laurain | Refined instrumental variable methods for identification of LPV Box-Jenkins models[END_REF], a refined instrumental variable (RIV) based approach has been introduced to solve the estimation problem of M θ in an open-loop setting. However, a major difference between the open-loop setting considered in [START_REF] Laurain | Refined instrumental variable methods for identification of LPV Box-Jenkins models[END_REF] and the closed-loop identification setting of Problem 1 is the correlation of the signal u with the noise e o due to the feedback loop. Thus in the next section, we focus on the problem how to determine a suitable instrument that allows the extension of the LPV-RIV approach to the closed-loop case.

REFORMULATION OF THE IDENTIFICATION PROBLEM

In [START_REF] Laurain | Refined instrumental variable methods for identification of LPV Box-Jenkins models[END_REF] it has been shown that reformulation of (4a-c) as a multiple-input single-output (MISO) LTI model allows an elegant solution of the filtering problem associated with LPV-IV approaches. This reformulation is necessary to avoid the problem of non-commutativity of parameter-varying filters and hence required to solve Problem 1. Consequently, (4a-c) are rewritten as

χ(t k ) + na i=1 a i,0 χ(t k-i ) F (q -1 )χ(t k ) + na i=1 nα l=1 a i,l f l (p k )χ(t k-i ) χ i,l (t k ) = n b j=0 n β l=0 b j,l g l (p k )u(t k-d-j ) u j,l (t k ) (11a) v(t k ) = C(q -1 , η) D(q -1 , η) e(t k ) (11b) y(t k ) = χ(t k ) + v(t k ) (11c)
where F (q -1 ) = 1 + na i=1 a i,0 q -i and g 0 ( ) = 1. Note that in this way, the LPV-BJ model is rewritten as a MISO system with (n b + 1)(n β + 1) + n a n α inputs {χ i,l } na,nα i=1,l=1 and {u j,l } n b ,n β j=0,l=0 . F (q -1 ) does not depend on p k , and (11a) can be rewritten as the following LTI model:

y(t k ) = - na i=1 nα l=1 a i,l F (q -1 ) χ i,l (t k ) + n b j=0 n β l=0 b j,l F (q -1 ) u k,j (t k ) + H(q)e(t k ), (12)
Note that ( 12) is an equivalent form of (4a-c), but it is not a representation of the associated LPV system as it includes lumped output variables in the signals {χ i,l } na,nα i=1,l=1 .

RIV APPROACH FOR CLOSED-LOOP LPV SYSTEMS

Using the MISO-LTI formulation (12), the estimation problem of the parameters can be formulated as a linear regression which allows optimal solution of (9) w.r.t. LPV-BJ models. To achieve this solution, one possible way is to develop an extension of the refined instrumental variable (RIV) approach of the closed-loop LTI identification framework. Next we derive this extension, which, as we will see, provides an easily implementable iterative estimation scheme.

The refined instrumental variable estimate

Using (12), y(t k ) can be written in a linear regression form:

y(t k ) = ϕ (t k )ρ + ṽ(t k ), (13)
where ṽ(t k ) = F (q -1 , ρ)v(t k ) and

r 2 (t k ) C o r 1 (t k ) ů(t k ) p k χ(t k ) G o Fig. 2. Auxiliary model ϕ(t k ) = [ -y(t k-1 ) . . . -y(t k-na ) -χ 1,1 (t k ) . . . -χ na,nα (t k ) u 0,0 (t k ) . . . u n b ,n β (t k ) ] , ρ = [ a 1,0 . . . a na,0 a 1,1 . . . a na,nα b 0,0 . . . b n b ,n β ] .
Two difficulties remain to obtain the minimum of (9) based on (13): the regressor ϕ(t k ) contains unknown terms {χ i,l (t k )} na,nα i=1,l=0 and all of its elements are corrupted with the colored noise v(t k ). To resolve this problem, an appropriate instrumental variable ζ(t k ) can be introduced such that the (extended-IV) estimate of ρ can be given as [START_REF] Söderström | Instrumental Variable Methods for System Identification[END_REF]:

ρXIV (N ) = arg min ρ∈R nρ 1 N N k=1 L(q)ζ(t k )L(q)ϕ (t k ) ρ - 1 N N t=1 L(q)ζ(t k )L(q)y(t k ) 2 W , (14)
where x 2 W = x T W x, with W a positive definite weighting matrix and L(q) is a stable prefilter. If G o ∈ G, the estimate ( 14) is consistent under the following well-known conditions 2 :

C1 Ē{L(q)ζ(t k )L(q)ϕ (t k )} is full column rank. C2 Ē{L(q)ζ(t k )L(q)ṽ(t k )} = 0.
Moreover it has been shown in [START_REF] Söderström | Instrumental Variable Methods for System Identification[END_REF] and [START_REF] Young | Recursive Estimation and Time-Series Analysis[END_REF] that in general for IV methods, based on ( 14), the minimum variance estimator can be achieved if:

C3 W = I. C4 ζ(t k
) is chosen as the noise-free version of (13) and is therefore defined in the present LPV case as:

ζ(t k ) = [ -χ(t k-1 ) . . . -χ(t k-na ) -χ 1,1 (t k ) . . . -χ na,nα (t k ) ů0,0 (t k ) . . . ůn b ,n β (t k )
where ů and χ are the signals from the auxiliary system as presented in Fig. 2 and where the index notations are the same as presented in ( 12). C5 G o ∈ G and n ρ is equal to the minimal number of parameters required to represent G o with the considered model structure. C6 L(q) is chosen as in [START_REF] Laurain | Refined instrumental variable methods for identification of LPV Box-Jenkins models[END_REF]:

L(q) = D(q -1 ) F (q -1 )C(q -1 ) . ( 15
)
Full column rank of Ē{L(q)ϕ(t k )L(q)ϕ (t k )} is ensured by the informativity of the data set (A3). Note that conditions about informativity w.r. E{.} is adopted from the prediction error framework of [START_REF] Ljung | System Identification, theory for the user[END_REF].

(2002); [START_REF] Tóth | Modeling and Identification of Linear Parameter-Varying Systems[END_REF]; [START_REF] Wei | On persistent excitation for parameter estimation of quasi-LPV systems and its application in modeling of diesel engine torque[END_REF] for results in the ARX case). To fulfill C1 under A3, the discussion can be found in [START_REF] Söderström | Instrumental Variable Methods for System Identification[END_REF]. It is also important to note that in a practical situation none of F (q -1 , ρ), C(q -1 , η), D(q -1 , η) or {a i,l (ρ)} na,nα i=1,l=0 , {b j,l (ρ)} n b ,n β j=0,l=0 is known a priori. Therefore, the RIV estimation normally involves an iterative (or relaxation) algorithm in which, at each iteration, an 'auxiliary model' is used to generate the instrumental variables (which guarantees C2), as well as the associated prefilters. Convergence of such an iterative RIV algorithm has not been proved so far and is only empirically assumed [START_REF] Young | The refined instrumental variable method: Unified estimation of discrete and continuous-time transfer function models[END_REF]. Furthermore, the considered LPV model can be reformulated in a LTI-MISO form only under the condition that the noise-free output terms are a priori known. Therefore, even if the presented method considerably lowers the variance in the estimated parameters, the optimality cannot be guaranteed.

Iterative closed-loop LPV-RIV Algorithm

Based on the previous considerations, the iterative scheme of the RIV algorithm can be extended to the closed-loop LPV case as follows. Algorithm 1. (LPV-RIV).

Step 1 Assume that as an initialization, an ARX estimate of M θ is available by the LS approach, i.e. θ(0

) = [ (ρ (0) ) (η (0) ) ] is given. Set τ = 0. Step 2 Compute an estimate of χ(t k ) via A(p k , q -1 , ρ(τ) ) χ(t k ) = B(p k , q -1 , ρ(τ) )u(t k-d ),
where ρ(τ) is estimated in the previous iteration and it is assumed to be stable in analog with A5. Based on M θ(τ) , deduce { χi,l (t k )} na,n α i=1,l=0 as given in (11a). It is important to note that χ(t k ) is correlated with the noise as it is simulated using the measured signal u.

Step 3 Compute the estimated filter:

L(q -1 , θ(τ) ) = D(q -1 , η(τ) ) C(q -1 , η(τ) )F (q -1 , ρ(τ) )
and the associated filtered signals {u

f j,l (t k )} n b ,n β j=0,l=0 , y f (t k ) and { χf i,l (t k )} na,nα i=1,l=0 .
Step 4 Build the noisy filtered estimated regressor φf (t k ) as:

φf (t k ) = [ -y f (t k-1 ) . . . -y f (t k-na ) -χf 1,1 (t k ) . . . -χf na,nα (t k ) u f 0,0 (t k ) . . . u f n b ,n β (t k )
] and compute the filtered instrument ζf (t k ) by simulating the estimated auxiliary model according to Fig. 2:

ζf (t k ) = [ -χf (t k-1 ) . . . -χf (t k-na ) -χf 1,1 (t k ) . . . -χf na,nα (t k ) ûf 0,0 (t k ) . . . ûf n b ,n β (t k ) ]
Note that this instrument fulfills C1 and C2.

Step 5 The IV optimization problem can now be stated in the form ρ(τ+1) (N ) = arg min

ρ∈R nρ 1 N N k=1 ζf (t k ) φ f (t k ) ρ - 1 N N k=1 ζf (t k )y f (t k ) 2 (16)
where the solution is obtained as

ρ(τ+1) (N )= N k=1 ζf (t k ) φ f (t k ) -1 N k=1 ζf (t k )y f (t k ).
The resulting ρ(τ+1) (N ) is the IV estimate of the process model associated parameter vector at iteration τ + 1 based on the prefiltered input/output data.

Step 6 An estimate of the noise signal v is obtained as

v(t k ) = y(t k ) -χ(t k , ρ(τ) ). ( 17
)
Based on v, the estimation of the noise model parameter vector η(τ+1) follows, using e.g. ARMA estimation algorithm of the MATLAB identification toolbox (an IV approach can also be used for this purpose, see [START_REF] Young | The refined instrumental variable method: Unified estimation of discrete and continuous-time transfer function models[END_REF]).

Step 7 If θ (τ +1) has converged or the maximum number of iterations is reached, then stop, else increase τ by 1 and go to Step 2.

Based on a similar concept, the so-called simplified LPV-RIV (LPV-SRIV) method, can also be developed for the estimation of LPV-OE models. This method is based on a model structure (4a-c) with C(q -1 , η) = D(q -1 , η) = 1 and consequently, Step 6 of Algorithm 1 can be skipped.

SIMULATION EXAMPLE

Next, the performance of the proposed RIV method with its SRIV variant are compared via a simulation example.

The considered data-generating system is described as

S o        A o (q, p k ) = 1 + a o 1 (p k )q -1 + a o 2 (p k )q -2 B o (q, p k ) = b o 0 (p k )q -1 + b o 1 (p k )q -2 H o (q) = 1 1 -q -1 + 0.2q -2 (18)
where v o (t k ) = H o (q)e o (t k ) and

a o 1 (p k ) = 1 -0.5p k -0.1p 2 k , (19a)
a o 2 (p k ) = 0.5 -0.7p k -0.1p 2 k , (19b)
b o 0 (p k ) = 0.5 -0.4p k + 0.01p 2 k , (19c) b o 1 (p k ) = 0.2 -0.3p k -0.02p 2 k . (19d)
According to Fig. 1, S o is considered with a stabilizing feedback controller given in this case as

u(t k ) = 1 + 0.5q -1 1 -0.85q -1 (r 2 (t k ) -y(t k )) + r 1 (t k ), (20)
To generate data sets D N of S o for identification and validation, the closed-loop system is simulated using p k = 0.5 sin(0.35πk) + 0.5 and r 1 (t k ) taken as a white noise with a uniform distribution U(-1, 1) and with length N = 4000. In the simulation, the white noise disturbance e o (t k ) ∈ N (0, σ 2 ) is considered with different σ 2 > 0 to investigate low and high noise conditions.

In the sequel, the conventional LPV Least Square (LS) method [START_REF] Bamieh | Identification of linear parameter varying models[END_REF] is compared to the proposed IV approaches. The least square method is not suited for OE models and may not be a good candidate for comparison. Nonetheless, this method is known to result in parameter estimates with low variance and is therefore a good indicator for the efficiency of the RIV estimator.

In order to the demonstrate the efficiency of the proposed RIV estimators without effect of under-modeling and overparametrization, full structural knowledge about system is assumed. Thus the proposed LPV Closed Loop Refined Instrumental Variable method (LPV-CLRIV) is considered with the following LPV-BJ model:

M LPV-RIV θ        A(p k , q -1 , ρ) = 1 + a 1 (p k )q -1 + a 2 (p k )q -2 B(p k , q -1 , ρ) = b 0 (p k )q -1 + b 1 (p k )q -2 H(p k , q,η ) = 1 1 + d 1 q -1 + d 2 q -2
which represents the situation S o ∈ M. In analog, the LPV Simplified Closed-Loop Refined Instrumental Variable approach (LPV-SCLRIV) is used with the following LPV-OE model:

M LPV-SRIV θ    A(p k , q -1 , ρ) = 1 + a 1 (p k )q -1 + a 2 (p k )q -2 B(p k , q -1 , ρ) = b 0 (p k )q -1 + b 1 (p k )q -2 H(p k , q,η ) = 1 which represents the case when G o ∈ G, H o / ∈ H.
In both cases, full parametrization of the associated coefficient functions is considered:

a 1 (p k ) = a 1,0 + a 1,1 p k + a 1,2 p 2 k (21a) a 2 (p k ) = a 2,0 + a 2,1 p k + a 2,2 p 2 k (21b) b 0 (p k ) = b 0,0 + b 0,1 p k + b 0,2 p 2 k (21c) b 1 (p k ) = b 1,0 + b 1,1 p k + b 1,2 p 2 k (21d)
To provide representative results, a Monte-Carlo simulation of N MC = 100 runs is accomplished at different noise levels: 20dB, 15dB, 10dB and 5dB. In Table 1, the norm of the bias (BN) ||ρ o -Ē(ρ)|| 2 and variance (VN) || Ē(ρ -Ē(ρ))|| 2 norms of the resulting estimated parameter vectors are given, where Ē is the mean operator over the Monte-Carlo simulation. The number of iterations (Nit) needed for convergence is also given. Table 1 shows that contrary to the LS method the IV methods are unbiased which is in line with the theoretical results. For SNR down to 5dB, the LPV-CLRIV produces a parameter variance very close to the one obtained with the LS method. Although the statistical optimality of the algorithm cannot be proved, this latter result shows on this example, that the LPV-CLRIV algorithm results in accurate estimates. The suboptimal LPV-CLSRIV method offers satisfying results, considering that the used noise model cannot capture the dynamics of v o . The variance in the estimated parameters is twice as much as in the LPV-CLRIV case but it is still close to the variance of the LS method.

CONCLUSION

In this paper, the extension of an LPV-RIV approach has been introduced to provide an efficient solution for the closed-loop identification of LPV systems with Box-Jenkins type of noise models. The approach has been formulated with the assumption of prior knowledge of the controller but without any restriction on its structure like LTI, LPV, nonlinear, etc. It has been shown that under given conditions the proposed approach provides consistent estimates and its performance has been demonstrated on a representative example including a comparison to the least-squares approach. The presented example has shown that the proposed procedure is robust to noise and that the obtained estimates are accurate.

Table 1 .

 1 2 -norm of the bias and variance of the estimated parameters at different SNR levels.

	Method		20dB	15dB	10dB	5dB
	LS	BN	1.5736	2.3922	2.0812 2.0908
		VN	0.0092	0.0140	0.0253 0.0326
	LPV	BN	0.0067	0.0145	0.0125 0.0052
	-SRIV	VN	0.0011	0.0035	0.0110 0.0319
		Nit	16	16	17	19
	LPV	BN	0.0060	0.0139	0.0185 0.0146
	-RIV	VN	0.000623	0.0018	0.0068 0.0206
		Nit	19	19	20	22

f : R n → R is a real meromorphic function if f = g/h with g, h analytic and h = 0.

Supported by the Netherlands Organization for Scientific Research (grant no. 680-50-0927).


Keywords: Computational modeling, Heart Failure, Cardiac resynchronization therapy

The analysis of follow-up data from patients suffering from heart failure is a difficult task, due to the complex and multifactorial nature of this pathology. In this paper, we present a coupled model, integrating a pulsatile heart into a model of the short to long-term regulations of the cardiovascular system. An interface method is proposed to couple these models, which present significantly different time scales. Results from a sensitivity analysis of the original and integrated models are proposed, with simulations reproducing the main effects of the short and long-term responses of an acute decompensated heart failure episode on a patient undergoing cardiac resynchronization therapy.

Index Terms-

II. INTEGRATION OF PULSATILE VENTRICLES

The original Guyton model is built around a central 'Circulatory Dynamics' module, which interacts with 17 peripheral modules corresponding to various physiological functions (autonomic control, renal function, renine angiotensine system,…). Ventricular function is represented by a static algebraic equation, providing mean ventricular outputs (QLO and QRO for the left and right Manuscript received April 11, 2011. This work was partly supported by the ANR BioSys SAPHIR project. V. Le Rolle, D. Ojeda and A.I. Hernández are with INSERM U642, Rennes, F35000, France and with Université de Rennes 1, LTSI, Rennes, F35000, France. (email: virginie.lerolle, david.ojeda, alfredo.hernandez@univ-rennes1.fr).

Embedding a cardiac pulsatile model into an integrated model of the cardiovascular regulation for heart failure follow-up

Virginie Le Rolle, David Ojeda and Alfredo I. Hernández H ventricles, respectively), which are computed as the product of the 'baseline' ventricular outflow and various other parameters, including the mean arterial pressure (PA), pulmonary pressure (PPA) and the autonomic effect on cardiac contractility (AUH 0 ). In order to integrate a pulsatile heart, the left (LV) and right (RV) ventricles of the Guyton model were substituted with pulsatile ventricular models. To our knowledge, the first attempt to integrate a pulsatile heart into the Guyton model was proposed by Werner et al. [START_REF] Werner | Simulation and prediction of cardiotherapeutical phenomena from a pulsatile model coupled to the Guyton circulatory model[END_REF]. However, their work focused on analysing the short-term response of the system, and no details on coupling with all the Guyton components were given. In addition, simulation results were not compared to Guyton's results. This substitution process is not simple, for it requires i) the definition of coupling transformations in order to preserve the numerical and physiological properties of the original model, ii) parameter identification for the proposed coupling transformations and iii) a sensitivity analysis providing information on the impact of integrating the new pulsatile model.

A. Coupling the Guyton and the Pulsatile Models

The first coupling interface concerns the hemodynamic variables. This coupling can be performed directly (without any transformation), since the physical units used in both models for these variables are the same. Atrial and arterial pressures of the Guyton model are connected as inputs to the pulsatile models and trans-valvular flows obtained from the pulsatile model are connected to the Guyton model. In order to obtain pulsatile variables, a time-varying elastance formalism, including ventricular interaction as proposed by Smith [START_REF] Smith | Simulation of cardiovascular diseases by including the autonomic nervous system into a minimal model[END_REF] was used. Ventricular elastances vary between values obtained from the End Systolic Pressure-Volume Relationship and the End Diastolic Pressure-Volume Relationship (EDPVR). End systolic (P es ) and end diastolic (P ed ) pressures are defined as:

P es (V)=E es (V-V d ) (1) P ed (V)= P 0 (exp(λ(V-V 0 )-1))
(2) where E es is the end systolic elastance; V d is the volume at zero pressure; P 0 , λ, and V 0 are the parameters defining the EDPVR. The pressure-volume relationship of each ventricle is calculated by: P(V) = e(t) P es (V)+(1-e(t)) P ed (V) (3) where e(t) represents the elastance function that will be defined later on. The septum is represented by a flexible common wall between the LV and the RV. The LV free wall volume (V lvf ) and the RV free wall volume (V rvf ) are defined as:

V lvf = V lv -V spt and V rvf = V rv -V spt (4) where V spt , V lv and V rv are respectively the septum, LV and RV volumes. The computation of the septum volume is the solution of the equation linking the septum pressure to the difference between left and right ventricular pressures:

P sept = P lv -P rv (5) P sept =e(t)E es,sept (V sept -V d,sept )+(1-e(t))P 0,sept (e λ(V-Vo) -1) (6)
The heart valves are represented by modulated resistances.

A second coupling interface deals with the modulation of the cardiac activity through continuous variables of the Guyton model representing the autonomic control of the chronotropic and inotropic effects (AUR 0 and AUH 0 respectively). Since AUR 0 and AUH 0 are dimensionless variables, the following linear transformations are applied: AUR = S AUR (AUR 0 -1) + B AUR (7) AUH = S AUH (AUH 0 -1) + B AUH [START_REF] Smith | Simulation of cardiovascular diseases by including the autonomic nervous system into a minimal model[END_REF] where S AUR and S AUH are sensitivity controllers and B AUR and B AUH are baseline controllers. These controller parameters have to be tuned to adjust the level of autonomic regulation.

A transformation based on an Integral Pulse Frequency Modulation (IPFM) model [START_REF] Rompelman | Measurement in heart-rate variability: Part1 -Comparative study of hear-rate variability analysis methods[END_REF] has been further defined to convert AUR into a series of pulses that will activate ventricular elastances. Each emitted pulse of the IPFM generates a variation of the ventricular elastance, which depends on AUR as follows:

e(t)=A•exp(B(t•AUR-C) 2 ) ( 9
)
where t is the time elapsed since the last activation pulse and A = 1; B = 80s -2 and C = 0.27s are the elastance parameters proposed in [START_REF] Smith | Simulation of cardiovascular diseases by including the autonomic nervous system into a minimal model[END_REF]. Finally, the end-systolic elastance E es is modulated by: E es =AUH•E es0 [START_REF] Thomas | SAPHIR: a physiome core model of body fluid homeostasis and blood pressure regulation[END_REF] where E es0 is the basal value for the end-systolic elastance.

B. Identification of the controller parameters

Controller parameters P=[S AUR , B AUR , S AUH , B AUH ] were identified by comparing the simulations obtained from the original Guyton model with those obtained from the proposed integrated, pulsatile model, during the 5 minute-simulation of a sudden severe muscle exercise, which is an original experiment described in the 1972 Guyton et al. paper [START_REF] Guyton | Circulation: overall regulation[END_REF] and presented in detail in [START_REF] Thomas | SAPHIR: a physiome core model of body fluid homeostasis and blood pressure regulation[END_REF]. This simulation is obtained by changing the exercise parameter (EXC) to 60 times its normal value at t = 30 sec. This effect generates approximately a 15-fold increase in the whole-body metabolic rate, triggering a regulatory response. Then, at t=2 min, the value of EXC is reset to normal.

The error function ε, which is minimized during the identification process, is computed as:

€ ε = Y j pulsatile (n) -Y j original (n) n=1 N ∑ j=1 6 ∑ ( 11
)
where n is the sample index, N is the number of simulated samples (equivalent to 5 min at a sampling period of 10 -2 min) and variables and correspond to detrended and scaled versions of the j-th output variable, obtained from the original and pulsatile versions of the model, respectively. The six output variables presented in Figure 1, which were validated against published data in [START_REF] Guyton | Circulation: overall regulation[END_REF], have been selected to calculate ε. In order to identify P, an evolutionary algorithm (EA) has been applied, as in our previous works [START_REF] Le Rolle | Model-Based Analysis of Myocardial Strain Data acquired by Tissue Doppler Imaging[END_REF]. EA are stochastic search methods, inspired by the theories of evolution and natural selection, which can be employed to find an optimal configuration for a given system [START_REF] Michalewicz | Genetic algorithms + data structures = evolution programs[END_REF]. The repeatability of the obtained optimal parameters was assessed applying the identification method four times, with different initial conditions.

C. Sensitivity Analysis

In order to assess the impact of integrating pulsatile ventricles, an input/output sensitivity analysis of the 'Circulatory Dynamics' module was performed. The screening method of Morris [START_REF] Morris | Factorial Sampling Plans for Preliminary Computational Experiments[END_REF] was chosen because it provides information on nonlinearities and interactions between variables, with limited computational costs. In this method, each model output y can be expressed as a function of k factors,

, where x i is the i-th input, which has been scaled in the interval [0,1]. Several realizations m=1,..,M, of x i ( ) can be defined, taking values in {1,1/(p-1),2/(p-1),…,1-Δ}, with equal probability, where

Δ=p/2(p-1) and p is the number of levels dividing the input-variable space. The Elementary Effect ( ) of y, with respect to is defined as:

€ EE i, m = y x 1 ,..., x i, m + Δ,..., x k ( ) -y x 1 ,..., x i, m, ..., x k ( ) Δ (12)
The mean (µ i ) and standard deviation (σ i ) of (for all m) are calculated. A large value of µ i indicates a significant overall effect of x i on y, while a large value of σ i indicates either non-linear or strong interactions with other inputs. All EE i for a given y can be analyzed in the σ i -µ i plane.

III. RESULTS

Table 1 shows the mean values and the standard deviations of the identified parameters obtained with the EA. Identified values for the sensitivity (S) and baseline (B) controllers for four realizations of the identification algorithm.

Using these mean parameter values, Figure 1 shows the comparison of the output of the pulsatile and original models for the simulation of a sudden severe muscle exercise, used during the identification process. In order to facilitate the comparison between both model outputs, pulsatile variables have been low-pass filtered. A close match is observed between both simulations. The mean relative root mean squared error (rRMSE) equals 0.0025. The proposed model provides simulations of pulsatile pressures and volumes for RV and LV. These variables do not exist in the original model. An example of these pulsatile variables is presented in Figure 2, with the simulation of ventricular Pressure-Volume (PV) loops, obtained by changing systemic resistance. The end-systolic PV relation has been found to be linear. These simulations are consistent with clinical observations [START_REF] Aroney | Linearity of the left ventricular end-systolic pressure-volume relation in patients with severe heart failure[END_REF]. Figure 3 shows the Morris input/output sensitivity results on the mean PA with p = 20 and M = 5⋅k (k = 16 and k = 17 respectively for the original and pulsatile models). In both cases, the most influential inputs are the plasma volume (VP), the autonomic regulation of vasoconstriction on arteries (AUM) and the vascular volume caused by relaxation (VVR). A slightly higher sensitivity to inputs that modulate the systemic resistance (ANM, ARM and AMM) is observed on the pulsatile model. These factors are more influential than AUH 0 , which is on the same sensitivity level in models. This is mainly due to the more realistic response of the pulsatile model to changes in afterload. IV. SIMULATION OF AN ADHF Figure 4 presents the main hemodynamic and regulatory variables represented in the proposed pulsatile model for the simulation of a stable HF state. At t = 24 hours of simulated time, parameter values E es0 and V d where reduced (E es0,LV = 0.7 mmHg/ml, E es0,RV = 0.5 mmHg/ml and V d,LV = 20 ml) to correspond to those observed from HF patients [START_REF] Aroney | Linearity of the left ventricular end-systolic pressure-volume relation in patients with severe heart failure[END_REF]. This reduced ventricular function causes a sudden decrease of PA and cardiac output and a significant increase in atrial pressures and ventricular preload, leading to an accumulation of fluid on the systemic and pulmonary spaces. In order to compensate for this hemodynamic response, neurohumoral regulations are initiated with a fast autonomic modulation (AU), combined with the slower response of the renin-angiotensin system (ANM). A stable HF state is reached, with an increased sympathetic tone, fluid retention and reduced PA and cardiac output. Finally, the simulation of an ADHF event is shown in Figure 5. The simulation starts from a stable HF state obtained, for example, after implanting a CRT device. A sudden desynchronization of both ventricles (inter-ventricular delay = 200 ms) is simulated at time t = 5 min. This event reproduces a sudden loss of capture of the LV lead that can be observed on CRT patients due to lead displacement. The model response presents a decreased PA and a regulatory response which are in accordance with clinical observations [START_REF] Whinnett | Haemodynamic effects of changes in atrioventricular and interventricular delay in cardiac resynchronisation therapy show a consistent pattern: analysis of shape, magnitude and relative importance of atrioventricular and interventricular delay[END_REF].

A) Original Guyton Model B) Pulsatile Model

Fig. 5: Simulation of a HF state. PA (mean arterial pressure in mmHg), QLO (cardiac output in l.min -1 ), AU (autonomic activity, ratio to normal).

V. DISCUSSION AND CONCLUSION This paper presents an example of temporal multiscale integration in which the non-pulsatile ventricles of the original Guyton model are replaced by a pulsatile, elastance-based model of the cardiac function. In order to perform this integration, an interfacing method applied to couple these heterogeneous models is presented. Although the proposed model is already useful for analyzing the main interaction effects that may be considered for the development of new ADHF detection methods, it still has to be improved and validated. To that end, an initial validation of some model components, such as the reproduction of patientspecific trans-valvular flows for different CRT pacing configurations using the time-varying elastance model with pulsatile atria and ventricles [START_REF] Ve Koon | Atrioventricular Delay Optimization in Cardiac Resynchronization Therapy Assessed by a Computer Model[END_REF], has been undertaken. Moreover, the proposed interfacing approach has been applied to integrate improved versions of other important components of the model, such as the renin angiotensin system [START_REF] Guillaud | A computational model of the circulating renin-angiotensin system and blood pressure regulation[END_REF]. The validation of the global, interconnected model is a challenging task, mainly because of observability limitations in long-term monitoring. The data captured from new-generation CRT devices will be useful to tackle this issue.

Fig. 1 :

 1 Fig. 1: Comparison of the output of the pulsatile model (black curves) with the original Guyton model (dashed curves) during a 5 minute-simulation of sudden severe muscle exercise. PVO (muscle venous oxygen pressure in mmHg), PMO (muscle cell oxygen pressure in mmHg), PA (mean arterial pressure in mmHg), AUP (sympathetic stimulation, ratio to normal), QLO (cardiac output in l.min -1 ) and BFM (muscle blood flow in l.min -1 ). Black lines are filtered versions of the pulsatile signals, obtained by integrating these signals on each cardiac cycle and dividing by the cardiac period.

Fig. 2 :

 2 Fig. 2: Simulations obtained from the pulsatile model on LV and RV pressures, systemic (PA) and pulmonary (PPA) arterial pressures. PV loops are simulated using different values for the systemic resistance of the Guyton model (R1= 11, R2 = 23 and R3 = 46 mmHg.min/l).

Fig. 3 :

 3 Fig.3: Morris sensitivity results for the arterial pressure obtained by the original (A) and pulsatile (B) models. The inputs are: HPL (hypertrophy effect on LV), HPR (hypertrophy effect on RV), VP (plasma volume), VVR (basic venous volume), AUH 0 (autonomic effect on heart strength), AUM (sympathetic vasoconstrictor effect on arteries), ARM (non-muscle global autoregulation multiplier), AMM (muscle autoregulation multiplier), RBF (renal blood flow), VRC (volume of red blood cells), VV7 (vascular volume due to short-term stress relaxation), ANM (general angiotensin multiplier effect), PC (capillary pressure), AUY (sensitivity of sympathetic control of veins), VIM (blood viscosity effect on resistance), HMD (cardiac depressant effect of hypoxia)

Fig. 4 :

 4 Fig.4: Simulation of a heart failure state. PA (arterial pressure in mmHg), PRA (mean right atrial pressure in mmHg), QLO (cardiac output in l.min -1 ), VEC (extracellular fluid volume in l), AU (autonomic activity, ratio to normal), ANM (angiotensin multiplier effect, ratio to normal).

TABLE 1

 1

	S AUR	0.71 (+/-0.017)	B AUR	0.87 (+/-0.009)
	S AUH	0.33 (+/-0.028)	B AUH	0.29 (+/-0.006)