Rachit Agarwal
email: rachit.agarwal@inria.fr

Nikolaos Georgantas
email: nikolaos.georgantas@inria.fr

Valerie Issarny
email: valerie.issarny@inria.fr

Impact of Spatial Scale on Mobility Parameters

• Huge body of work focusing on the understanding of human mobility characteristics [1] utilizing various types of data [2] • How spatial scale affect the human mobility parameters is still the question?

• Positive lognormal distribution fits best (See Figure 2 and Table 1) both jump length and R g • As expected spatial scale affects the distribution parameters • Other results:

• Pause time also follow positive lognormal distribution but bias present due to nature of the dataset (crowdsourced dataset) used. • Distribution of max and min jump length from users also follow positive lognormal distribution with parameters µ=9.95 σ=2.75 and µ= 3.68 σ= 1.70 respectively at fine grain scale.

The findings reported here will be leveraged for revisiting the identified parameters, based on further analysis using complementary datasets.

[1] L.

•

  FIGURE 2: Log-Log plot of best fit distribution for jump-length and R g at different spatial scales with different values of x min . Top row (Jump length): (a) fine grain samples, (b) samples mapped to arrondissements, (c) samples mapped to Ile-de-France cities. Bottom row (R g ): (d) fine grain samples, (e) samples mapped to arrondissements, (f) samples mapped to Ile-de-France cities.

TABLE 1 :

 1 Distribution fit with identified parameters.

TABLE 2 :

 2 KL divergence between identified distributions at different spatial scale.

	x min in meters	power law α D	positive log-normal µ σ D	exponential λ D	truncated power law α β in km D

FIGURE 1: Ambiciti application (formerly called as SoundCity) monitors the environmental pollution the user is exposed to [3].


Yassine Abidi
email: yassine.abidi@enit.rnu.tn

Mourad Bellassoued
email: mourad.bellassoued@enit.utm.tn&moncef.mahjoub@enit.rnu.tn.

Moncef Mahjoub

Nejib Zemzemi
email: nejib.zemzemi@inria.fr.

IONIC PARAMETERS IDENTIFICATION OF AN INVERSE PROBLEM OF STRONGLY COUPLED PDE'S SYSTEM IN CARDIAC ELECTROPHYSIOLOGY USING CARLEMAN ESTIMATES

Keywords: Mathematics Subject Classification. 35Q92, 34A55. November 15, 2018 Lipschitz stability estimate, carleman estimate, cardiac electrophysiology, bidomain system, physiological ionic model, ionic parameters

In this paper, we consider an inverse problem of determining multiple ionic parameters of a 2 × 2 strongly coupled parabolic-elliptic reaction-diffusion system arising in cardiac electrophysiology modelling. We use the bidomain model coupled to an ODE system and we consider a general formalism of physiologicaly-detailed cellular membrane models to describe the ionic exchanges at the microscopic level. Our main result is the uniqueness and a Lipschitz stability estimate of the ion channels conductance parameters of the model using subboundary observations over an interval of time. The key ingredients are a global Carleman-type estimate with a suitable observations acting on a part of the boundary.

Introduction

The electric wave in the heart is governed by a system of reaction-diffusion partial differential equations called the bidomain model. This system is coupled nonlinearly to an ordinary differential equations (ODEs) modeling the cellular membrane dynamics. The cellular membrane electrical activity model is based on the Hodgkin and Huxley (HH) formalism [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF] which has been adapted to cardiac Purkinje cells by Denis Nobel [START_REF] Noble | Cardiac action and pacemaker potentials based on the hodgkin-huxley equations[END_REF]. After that, many other model have been introduced to describe the electrical activity of the cell membrane in the myocardium. In 1977, Beeler and Reuter [START_REF] Beeler | Reconstruction of the action potential of ventricular myocardial fibres[END_REF] introduced a ventricular cells model. Di Francesco and Noble proposed in 1985 a model that takes into account ion pumps, which allows different chemical species such as the sodium potassium and calcium to regain their stable states [START_REF] Francesco | A model of cardiac electrical activity incorporating ionic pumps and concentration changes[END_REF]. Rudy and his collaborators proposed series of models based on the HH formalism which are chronologically more and more complex [START_REF] Li | A model of canine purkinje cell electrophysiology and ca2+ cycling[END_REF][START_REF] Luo | A model of the ventricular cardiac action potential: Depolarization, repolarization, and their interaction[END_REF][START_REF] Luo | A dynamic model of the cardiac ventricular action potential. i. simulations of ionic currents and concentration changes[END_REF][START_REF] O'hara | Simulation of the undiseased human cardiac ventricular action potential: Model formulation and experimental validation[END_REF][START_REF] Shaw | Electrophysiologic effects of accute myocardial ischemia: Theoretical study of altered cell exitability and action potential duration[END_REF] and take into account more and more the physiological behaviors of the ion channels. Other models like [START_REF] Buen-Orovio | Minimal model for human ventricular action potentials in tissue[END_REF][START_REF] Courtemanche | Ionic mechanisms underlying human atrial action potential properties: insights from a mathematical model[END_REF][START_REF] Tusscher | A model for human ventricular tissue[END_REF] have been also extensively used in the computational electrophysiology community. In all these models, the ion channels maximal conductance parameter plays an important role in generating the action potential but also to consider some pathological conditions or when taking into account the effect of drugs. In [START_REF] Mcintyre | Sensitivity analysis of a model of mammalian neural membrane[END_REF], authors studied the sensitivity of the strength and duration of the action potential in the parameters describing the sodium channel. In [START_REF] Zemzemi | Computational assessment of drug-induced effects on the electrocardiogram: from ion channel to body surface potentials[END_REF], authors demonstrate the ability of computational models to simulate the effect of drug action on the electrical activity of the heart, at the level of the ion-channel, cell, heart and ECG body surface potential. In [START_REF] Paci | Computational modelling of lqt1 in human induced pluripotent stem cell derived cardiomyocytes[END_REF][START_REF] Paci | Computational modeling supports induced pluripotent stem cell-derived cardiomyocytes reliability as a model for human lqt3[END_REF], authors showed how, with few changes in the maximal conductances of the ion currents, they reproduced different control action potentials. In [START_REF] Abbate | In silico assessment of the effects of various compounds in mea/hipsc-cm assays: Modeling and numerical simulations[END_REF], the authors propose a mathematical approach for the analysis of drugs effects on the electrical activity of hiPSC-CMs based on multi-electrode array experiments where the drug acts directly on the maximal conductance of the targeted ion channel. All these examples demonstrate the importance of identifying or estimating the maximal conductance parameters when dealing with numerical simulations in multi-scale models. Despite their importance, ion channel maximal conductances have not been subject of abondant theoretical studies. The only theoretical work that we know is [START_REF] Abidi | On the identification of multiple space dependent ionic parameters in cardiac electrophysiology modelling[END_REF], which proves under certain conditions the stability of the ionic conductance parameters identification problem in case of the monodomain approximation. Two other works has been dedicated to the parameters of phenomenological ionic models still in the case of the monodomain approximation. In [START_REF] Boulakia | Stability estimates for some parameters of a reaction-diffusion equation coupled with an ode[END_REF], authors showed Lipschitz stability inequalities for the identification of some parameters of the FitzHugh-Nagumo model from measurements on the cardiac potential and the ionic variable. The paper [START_REF] Lassoued | Stability results for the parameter identification inverse problem in cardiac electrophysiology[END_REF] shows the stability of the identification a reaction parameter in the Mitchell-Schaeffer model for the monodomain system. To the best of our knowledge no work proves the stability of the ion channels conductances in the case of the bidomain model. This would be the subject of our paper.

In the framework of the bidomain model, there some works that have been dedicated to study the stability of the cardiac tissue conductivities identifiability problem. The paper by Ainseba et al. [START_REF] Ainseba | Stability of conductivities in an inverse problem in the reaction-diffusion system in electrocardiology[END_REF] obtains the stability results for the conductivities diffusion coefficients to a strongly reaction-diffusion system modeling electrical activity in the heart using Carleman estimate. But their method is based on an approximation of the monodomain model and establish a Carleman estimate for a family of parabolic equations. Recently, Wu and Yu [START_REF] Wu | Hölder stability of an inverse problem for a strongly coupled reaction-diffusion system[END_REF] convert the strongly coupled terms to the derivatives of the state variables components. The authors establish a Hölder stability result for the inverse conductivities problem for a linearized reaction term in the bidomain model in electrophysiology where the conductivities are scalar functions. However, considering a linear reaction term significantly reduces the complexity of the problem. Authors also consider a linear and simple approximation of the dynamic system governing the cell membrane activity at the microscopic level. Both simplification are non-realistic but also non phenomenologically based models in the sense that one could not generate the shape of an action potential with these models. However, the trick that they proposed to weaken the coupling in the bidomain system is very interesting and gives a good way to deal with a strongly coupled problem.

Strongly coupled systems are notoriously difficult to treat. Many of the standard results are not valid for a strongly coupled system. For example, there is no maximum principle nor DeGiorgi-Nash-Moser type estimates nor Carleman type estimates unless the system has special structures. In literature, it has been mostly assumed that the systems are regular elliptic and bounndedness of solutions are known in the study of the direct problem. Many inverse problems for the parameter identification use Carleman's estimates which are quite difficult to deal with, because the equations are coupled in the highest derivatives terms. Until now, only partial results were available in the literature concerning the Carleman estimates for strongly coupled systems. Different from the existing methods dealing with weakly or strongly coupled system, such as Fan and Chen [START_REF] Fan | Stability estimate for a strongly coupled parabolic system[END_REF] and Bellassoued & Yamamoto [START_REF] Bellassoued | Carleman estimate and inverse source problem for biot's equations describing wave propagation in porous media[END_REF]. Carleman estimates for a weakly coupled parabolic system could be derived by adding up all Carleman estimates corresponding to each equation. The weakly coupled terms are absorbed by the terms on the left-hand side of the obtained inequality in this process, see Benabdallah et al. [START_REF] Benabdallah | Inverse problem for a parabolic system with two components by measurements of one component[END_REF], Cristofol et al. [START_REF] Cristofol | Inverse problems for a 2 × 2 reaction-diffusion system using a carleman estimate with one observation[END_REF]. As for a strongly coupled parabolic or hyperbolic system, the method above does not work, because the strongly coupled terms could not be simultaneously absorbed by adding up each Carleman estimate. However, to our best knowledge, if the matrix (a i,j ) 1≤i,j≤2 formed by the parameters a 1,1 , a 2,2 describing the self diffusion of the state variables and a 1,2 , a 2,1 describing the influence of each component on another, is diagonalizable then the transformed system satisfies a weakly coupled system. In this case, Carleman estimate is obtained similarly to that of the weakly coupled system, see [START_REF] Bellassoued | Carleman estimate and inverse source problem for biot's equations describing wave propagation in porous media[END_REF][START_REF] Fernández-Cara | Controlling linear and semilinear systems formed by one elliptic and two parabolic pdes with one scalar control[END_REF][START_REF] Wu | Determination of an unknown source for a thermoelastic system with a memory effect[END_REF]. But due to the fact that the anisotropy in the bidomain equations depends on the fiber directions, and the fact that the fiber direction are space dependent, one cannot diagonalize the problem as in [START_REF] Bellassoued | Carleman estimate and inverse source problem for biot's equations describing wave propagation in porous media[END_REF][START_REF] Wu | Determination of an unknown source for a thermoelastic system with a memory effect[END_REF].

In this paper, we study the stability of the inverse problem of identification of these conductances parameters for a bidomain model: the reconstruction of these conductances from the measurement of electrical potentials over a part of the space boundary during the time and the state variables at a suitable time t 0 . The stability results is based on a global Carleman estimate for a non-linear coupled system with one observation. Our inverse stability results are new because the bidomain system contains a strong coupling term. The technics we shall discuss are similar to the framework using Carleman estimates for inverse problems but the obtained estimates differs from those of [START_REF] Cristofol | Inverse problems for a 2 × 2 reaction-diffusion system using a carleman estimate with one observation[END_REF][START_REF] Sakthivel | Stability of diffusion coefficients in an inverse problem for the lotka-volterra competition system[END_REF][START_REF] Yuan | Lipschitz stability in the determination of the principal part of a parabolic equation[END_REF] because of the strongly coupled terms. We use the same strategy introduced by Wu and Yu in [START_REF] Wu | Carleman estimate for a linearized bidomain model in electrocardiology and its applications[END_REF] to convert the strongly coupled terms to the gradient of the extra-cellular potential u e . Then, we consider the new version of the bidomain model as a whole to establish Carleman estimate under a weaker coupled equations. The paper is organized as follows: In Section 2, we briefly recall the general structure of cardiac cellular membrane models describing the transmembrane potential and the ionic exchange at the cell membrane. Then, we present the bidomain model describing the electrical wave propagation and recall some existence and uniqueness results that have been shown in [START_REF] Veneroni | Reaction-diffusion systems for the macroscopic bidomain model of the cardiac electric field[END_REF] and we deal with a regularity results of the bidomain solution which will be useful in stability analysis. In section 3, we announce the main stability result including the conditions we need for the identification of multiple parameters. The proof of the main result is divided into two sections. In Section 4, we prove the global Carleman inequality for the reaction-diffusion system. Most of the non-classical parts of the proof of the main result are presented in Section 5 where we prove the stability estimate of conductances parameters.

Mathematical Model

Let Ω ⊂ R n (n ≥ 1) be a bounded connected open set whose boundary Γ = ∂Ω is regular enough, (Ω ⊂ R 3 being the natural domain of the hearth). Let T > 0 and Γ 0 ⊂ Γ be a relatively open subset of Γ. We will use the notation Q = Ω × (0, T ), Σ = Γ × (0, T ) and Σ 0 = Γ 0 × (0, T ).

We introduce a two reaction-diffusion parabolic PDEs system called bidomain model, coupled to a system of ODEs. This model was proposed in the late 1970s by Tung [START_REF] Tung | A bidomain model for describing ischemic myocardial d-c potentials[END_REF] and is now the generally accepted model of electrical behaviour of cardiac tissue (see Henriquez [START_REF] Henriquez | Simulating the electrical behavior of cardiac tissue using the bidomain model[END_REF], Keener and Sneyd [START_REF] Keener | Mathematical Physiology, I: Cellular Physiology[END_REF]), can be written as:

                     c m ∂ t v + I ion (¯ , v, w, z) = div(σ i ∇u i ) + I s i in Q, c m ∂ t v + I ion (¯ , v, w, z) = -div(σ e ∇u e ) -I s e in Q, ∂ t w = F (v, w) in Q, ∂ t z = G(¯ , v, w, z) in Q, σ i ∇u i .ν = σ e ∇u e .ν = 0 on Σ, (2.1)
where u i = u i (x, t) and u e = u e (x, t) represent, respectively, the intracellular and extracellular electric potentials, whose difference v := u i -u e is the transmembrane potential. w := w(x, t) = (w 1 , . . . , w k )(x, t) and z := z(x, t) = (z 1 , . . . , z m )(x, t) are the gating variables and the ionic intracellular concentration variables, respectively. The surface capacitance of the membrane is represented by the constant c m > 0. For simplicity, we shall suppose c m = 1. Here, I s i := I s i (x, t) and I s e := I s e (x, t) are respectively the internal and the external applied current sources.

In an isolated heart, no current flows out of the heart, as expressed by the Neumann conditions. We complete this model with initial data:

v(•, 0) = v 0 , w(•, 0) = w 0 , z(•, 0) = z 0 , x ∈ Ω.
(2.

2)

The ionic current I ion := I ion (x, t) and the functions F and G depends of the considering ionic model.

In the membrane model, the ionic current I ion has the following general structure [START_REF] Franzone | Mathematical Cardiac Electrophysiology[END_REF]:

I ion (¯ , v, w, z) = N i=1 ¯ i y i (v) k j=1 w pj,i j (v -E i (z)), (2.3)
where N is the number of ionic currents, ¯ i := ¯ i (x) is the maximal conductance associated with the i th current, y i is a gating function depending only on the membrane potentiel v, p j,i are positive integers exponents and E i is the reversal potential for the i th current I i , which is the related equilibrium (Nernst) potential and is given by

E i (z) = γ i log z e z i , z = (z 1 , . . . , z m ), (2.4)
where γ i is a constant and z i , i = 1, . . . , m, are the intracellular concentrations. The constant z e denotes an extracellular concentration. The anisotropic properties of the two media are modeled by an intracellular and extracellular conductivity tensors σ i (x) and σ e (x). Generally, the conductivities σ i and σ e are two matrices given by σ j (x) = σ t j (x)I + (σ l j (x) -σ t j (x))a l (x)a T l (x), where σ l j and σ t j , j ∈ {i, e} are the intra-and extracellular conductivities along and transversal to the direction of the fiber (parallel to a l (x)), respectively. In the case of equal anisotropy, i.e the so-called anisotropy ratios σ t i /σ l i = σ t e /σ l e = 1, the σ i and σ e are simplified as σ i (x) = σ i (x)I and σ e (x) = σ e (x)I with σ i (x) = σ t i (x), σ e (x) = σ t e (x), which is the case we discussed. For the ODEs, the dynamics of the gating variable w is described in the Hodgkin-Huxley formalism by a system of ordinary differential equations which when w j is a gating variable (0 ≤ w j ≤ 1) are governed by the following equation,

∂ t w j = F j (v, w j ) := α j (v)(1 -w j ) -β j (v)w j , j = 1, . . . , k, (2.5)
where α j and β j are positive rational functions of exponentials in v. A general expression for both α j and β j is given by

µ 1 e µ2(v-vn) + µ 3 (v -v n ) 1 + µ 4 e µ5(v-vn) , (2.6)
where µ 1 , µ 3 , µ 4 , v n are non-negative constants and µ 2 , µ 5 are positive constants. The dynamics of the ionic concentration variables z is described by the additional system of ODEs:

∂ t z i = G i (¯ , v, w, z) := -J i (¯ , v, w, log z i ) + H i (¯ , v, w, z), i = 1, . . . , m. (2.7)
Assumptions 2.1. We assume that: 1) The conductivities of the intracellular and extracellular σ i , σ e ∈ C 2 ( Ω), such that 0 < σ j < σ j (x) < σ j , j ∈ {i, e}, in Ω. 2) We use the regularized form of the variable y i (v) in hyperbolic functions introduced in [START_REF] Endresen | A theory for the membrane potential of living cells[END_REF]. In this case, y i (v) is a C ∞ function with respect to the variable v for i = 1, . . . , N , and then is locally Lipschitz since v is bounded, similarly for the function y i (v)v. 3) J i and H i are locally Lipschitz continuous functions where:

J i ∈ C 2 (R * + × R × R k × R), 0 < g * (w) ≤ ∂J i ∂τ (¯ , v, w, τ ) ≤ g * (w), ∂J i ∂v (¯ , v, w, 0) ≤ L v (w), (2.8) g * , g * , L v belong to C 1 (R k , R + ), and
H i ∈ C 2 (R * + × R × R k × (0, +∞) m ) ∩ Lip(R * + × R × [0, 1] k × (0, +∞) m ).
(2.9)

Remark 2.2. One could find in the literature many refined models based on Hodgkin-Huxley formalism taking into account different quantities. For example, we recall here the following models: Beeler-Reuter ( [START_REF] Beeler | Reconstruction of the action potential of ventricular myocardial fibres[END_REF], N = 4, k = 6, m = 1), phase-I Luo-Rudy ( [START_REF] Luo | A model of the ventricular cardiac action potential: Depolarization, repolarization, and their interaction[END_REF], N = 6, k = 6, m = 1), phase-II Luo-Rudy ( [START_REF] Luo | A dynamic model of the cardiac ventricular action potential. i. simulations of ionic currents and concentration changes[END_REF], N = 10, k = 6, m = 5). Now, we propose a new approach to the bidomain equations (2.1). Our idea is to use a reformulation of (2.1) as a parabolic PDE coupled to an elliptic one, by replacing u i = v + u e in the first PDE. The boundary condition is also reformulated in terms of v and u e . We have:

                           ∂ t v -div(σ i ∇v) + I ion (¯ , v, w, z) = div(σ i ∇u e ) + I s i in Q, -div(σ e ∇u e ) -I ion (¯ , v, w, z) = ∂ t v + I s e in Q, ∂ t w = F (v, w) in Q, ∂ t z = G(¯ , v, w, z) in Q, σ i ∇v.ν + σ i ∇u e .ν = 0 on Σ, σ e ∇u e .ν = 0 on Σ, (2.10)
By using the second equation in (2.10), we obtain

-div(σ i ∇u e ) = -∇ • (σ i σ -1 e σ e ∇u e ) = -∇(σ i σ -1 e ) • σ e ∇u e -σ i σ -1 e div(σ e ∇u e ) = σ -1 e (σ i ∂ t v + σ i (I ion + I s e )) + A(x) • ∇u e , (2.11)
with A(x) := -∇(σ -1 e σ i )σ e = σ -1 e (σ i ∇σ e -σ e ∇σ i ).

(2.12) Substituting (2.11) into (2.10) yields

                           (1 + σ i σ -1 e )∂ t v -div(σ i ∇v) = -(1 + σ i σ -1 e )I ion (¯ , v, w, z) -A(x) • ∇u e + I s i -σ i σ -1 e I s e in Q, -div(σ e ∇u e ) = ∂ t v + I ion (¯ , v, w, z) + I s e in Q, ∂ t w = F (v, w) in Q, ∂ t z = G(¯ , v, w, z) in Q, σ i ∇v.ν + σ i ∇u e .ν = 0 on Σ,
σ e ∇u e .ν = 0 on Σ, (2.13) Hence we consider the problem of finding unknown functions v, u e , w and z verifying (2.13). Since v = u i -u e , it is natural decompose the initial condition v 0 as v 0 = u i,0 -u e,0 . Now, we state some regularity and estimate results proved in [START_REF] Veneroni | Reaction-diffusion systems for the microscopic cellular model of the cardiac electric field[END_REF][START_REF] Veneroni | Reaction-diffusion systems for the macroscopic bidomain model of the cardiac electric field[END_REF] of a strong solution to problem (2.1).

Lemma 2.3. Assume that (v 0 , w 0 , z 0 ) ∈ H 2 (Ω) × L 2 (Ω) k × L 2 (Ω) m , with log z 0 := (log z 0,1 , . . . , log z 0,m ) ∈ L 2 (Ω) m , and
I s i,e ∈ L p (0, T ; L 2 (Ω)) ∩ H 2 (0, T ; H 2 (Ω)), p > 4, (2.14) then we have 1) u i,e ∈ L p (0, T ; H 2 (Ω)), v ∈ W 1,p (0, T ; L 2 (Ω)) ∩ L 2 (0, T ; H 2 (Ω)) ∩ C 0 ([0, T ], C 0 (Ω)), for p > 4, w : Q → [0, 1] k measurable, z : Q → (0, +∞) m measurable, w j (x, .) ∈ C 1 (0, T ) ∩ C 0 ([0, T ]), for a.e., x ∈ Ω, j = 1, . . . , k, z i (x, .) ∈ C 1 (0, T ) ∩ C 0 ([0, T ]), for a.e., x ∈ Ω, i = 1, . . . , m, z ∈ H 1 (0, T ; L 2 (Ω)) m ∩ L ∞ (Q) m , log z := (log z 1 , . . . , log z m ) ∈ L ∞ (Q) m . 2) There exists a constant C > 0, independant of v, w, z, such that |z(x, t)| ≤ C(1 + |z 0 (x)| + v(x) L 2 (0,t) ), a.e., x ∈ Ω, (2.15
)

and |log z(x, t)| + |∂ t z(x, t)| ≤ C(1 + |z 0 (x)| + v(x) C 0 (0,t) ), a.e., x ∈ Ω, (2.16) ∀t ∈ [0, T ], for a.e. x ∈ Ω.
3) There exists a constant C > 0, depending on m, T such that, ∀i = 1, . . . , m,

z i (x, t) ≥ exp[-C(1 + z 0 L ∞ + v C 0 ([0,T ],C 0 (Ω)) )] > 0, ∀t ∈ [0, T ] , a.e., x ∈ Ω.
(2.17)

4) There exists M ∞ > 0, depending on the data of the problem, such that:

sup{|v(x, t)| : (x, t) ∈ Q} ≤ M ∞ . (2.18)
Now, we will establish regularity results for the solution of the bidomain system (2.13). The aim is to improve the regularity results given in Lemma 2.3 in order to satisfy some assumptions useful for the stability result.

Proposition 2.4. Let (v, u e , w, z) be the solution of bidomain system (2.1), with initial conditions (v 0 , u e,0 , w 0 , z 0 ). andI s i,e verify the regularity (2.14), then

• If v 0 , u e,0 ∈ H 2 (Ω), w 0 ∈ L 2 (Ω) k , z 0 ∈ L 2 (Ω) m ,
v, u e ∈ W 1,∞ (0, T ; H 1 (Ω)) ∩ L 2 (0, T ; H 2 (Ω)) ∩ H 1 (0, T ; H 1 (Ω)), w ∈ W 1,∞ (0, T ; L 2 (Ω)) k , and z ∈ W 1,∞ (0, T ; L 2 (Ω)) m .
(2. [START_REF] Endresen | A theory for the membrane potential of living cells[END_REF])

Moreover if w 0 ∈ H 1 (Ω) k , and z 0 ∈ H 1 (Ω) m , (2.20) then w ∈ W 1,∞ (0, T ; H 1 (Ω)) k , and z ∈ W 1,∞ (0, T ; H 1 (Ω)) m .
(2.21)

• If v 0 , u e,0 ∈ H 4 (Ω), w 0 ∈ H 2 (Ω) k , z 0 ∈ H 2 (Ω) m ,
and I s i,e verify the regularity (2.14). Then the solution of (2.1) satisfies

v, u e ∈ H 1 (0, T ; H 3 (Ω)) ∩ H 2 (0, T ; H 1 (Ω)), w ∈ W 1,∞ (0, T ; H 2 (Ω)) k ∩ H 2 (0, T ; L 2 (Ω)) k , z ∈ W 1,∞ (0, T ; H 2 (Ω)) m ∩ H 2 (0, T ; L 2 (Ω)) m . (2.22) Moreover if w 0 ∈ H 3 (Ω) k , and z 0 ∈ H 3 (Ω) m , (2.23) then w ∈ W 1,∞ (0, T ; H 3 (Ω)) k → H 1 (0, T ; H 3 (Ω)) k , z ∈ W 1,∞ (0, T ; H 3 (Ω)) m → H 1 (0, T ; H 3 (Ω)) m .
(2.24)

The proof of proposition 2.4 is provided in the Appendix 7.

Inverse problem: Main result

Our inverse problem is related to determination of multiple ionic parameter ¯ i , 1 ≤ i ≤ N , of a non linear parabolic reaction diffusion system coupled with an ordinary differential equations from a set of boundary observations (u e, , w , z ) |Γ0×(0,T ) and (u e, (x, t 0 ), w (x, t 0 ), z (x, t 0 )), x ∈ Ω, 1 ≤ ≤ N .

In order to formulate our results, we need to introduce the following notations: For a sequence functions

( v , w , z ) ∈ H 3 (Ω) × C 1 (Ω) k × C 1 (Ω) m , we define the N × N matrix Λ as follows Λ( v (x), w (x), z (x)) =         S 1,1 (x) S 2,1 (x) . . . S N,1 (x) S 1,2 (x) S 2,2 (x) . . . S N,2 (x) . . . . . . . . . . . . S 1,N (x) S 2,N (x) . . . S N,N (x)         , (3.1)
where

S i, (x) = y i ( v (x)) ( v (x) -E i ( z (x))) k j=1 ( w ) pj,i j (x) , 1 ≤ , i ≤ N .
Let us fix constant M 0 > 0. We introduce an admissible set of unknown coefficients vector ¯ by

A = ¯ ∈ H 3 (Ω) N , ¯ l 2 (H 3 (Ω) N ) ≤ M 0 . (3.2)
We obtain the following stability result.

Theorem 3.1. Let t 0 ∈ (0, T ), Γ 0 be a subdomain of ∂Ω and let ¯ (2) ∈ A be arbitrary fixed. We assume that

I s, i,e ∈ L p (0, T ; L 2 (Ω)) ∩ H 1 (0, T ; H 2 (Ω)), p > 4, 1 ≤ ≤ N , satisfy det Λ v (2) (x, t 0 ), w (2) (x, t 0 ), z (2) (x, t 0 ) = 0, ∀x ∈ Ω. (3.3)
Here v (2) , w (2) and z (2) are the state variables of (2.13) with ¯ = ¯ (2) and I s i,e = I s, i,e . Furthermore, we assume that v

(2)

C 0 ([0,T ];C 1 (Ω)) + w (2) C 0 ([0,T ];C 1 (Ω)) k + z (2) C 0 ([0,T ];C 1 (Ω)) m ≤ M, (3.4
) for some positive M . Then there exists a constant C > 0 such that we have:

¯ (1) -¯ (2) l 2 (L 2 (Ω)) N ≤ C N =1 (u (1)
e, -u

e, ) H 2 (0,T ;H 1 (Γ0))

+ (u (1)
e, -u

(2) e, )(•, t 0 ) H 2 (Ω) + (w (1) -w (2) )(•, t 0 ) L 2 (Ω) k + (z (1) -z (2) )(•, t 0 ) L 2 (Ω) m , (3.5)
for all ¯ (1) ∈ A.

Remark 3.2. The condition (3.4) is a straightforward consequence of the Proposition 2.4 and a Sobolev embedding theorem (e.g., Thm. 5.4 in [START_REF] Adams | Sobolev spaces[END_REF]).

H 1 (0, T ; H 3 (Ω)) → C 0 ([0, T ]; C 1 (Ω)), (3.6)
if the initial conditions verify

v (2) (t = 0) ∈ H 4 (Ω), u (2)
e, (t = 0) ∈ H 4 (Ω), w (2)

(t = 0) ∈ H 3 (Ω) k , and z (2) (t = 0) ∈ H 3 (Ω) m . (3.7)
We introduced the condition (3.4), because it may be obtained with a less regularities on the initial conditions.

By Theorem 3.1, we can readily derive the uniqueness in the inverse problem.

Corollary 3.3. Under the same assumptions as in Theorem 3.1 and if

(u (1)
e, (x, t 0 ), w (1)

(x, t 0 ), z (1) (x, t 0 )) = (u (2)
e, (x, t 0 ), w

2) (x, t 0 ), z (2) (x, t 0 )), in x ∈ Ω, (
e, (x, t) = u

(2)

e, (x, t), in Γ 0 × (0, T ), (3.9) for = 1, . . . , N , then ¯ (1) = ¯ (2) in Ω.

Since the number of the unknown coefficients is N , it is natural to expect that N -times observations can yield the Lipschitz stability. Our tool are a weighted L 2 -norm estimates of the solution of a PDE called Carleman estimates, where the weight takes an exponential form. They are an important tool in subjects in analysis of PDEs such as unique continuation [START_REF] Bukhgeim | Global uniqueness of class of multidimensional inverse problems[END_REF][START_REF] Yu | An inverse problem for the dynamical lamé system with two sets of boundary data[END_REF], control theory [START_REF] Chaves-Silva | Uniform null controllability for a degenerating reaction-diffusion system approximating a simplified cardiac model[END_REF][START_REF] Yu | Controllability of parabolic equations[END_REF] and coefficient inverse problems [4, 6, 7, 27-29, 32, 52].

Global Carleman inequality for reaction-diffusion system

In this section, we give Carleman estimate for the reaction-diffusion model with a finite observations acting on a subboundary Σ 0 of ∂Ω on the right-hand side of the estimate. This Carleman estimate would be used later for the stability and uniqueness of the solution of the parameter identification problem. We are interested in identifying the parameters ¯ i , i = 1, . . . , N , where ¯ i is the maximal conductance associated with the i th current.

In order to frame a Carleman type estimate, we shall first introduce a particular type of weight functions.

Weight functions and Carleman estimate

For our Carleman estimate, we need a weight function β with the following properties [START_REF] Yu | Controllability of parabolic equations[END_REF]. Proof 4.1. The proof of the Lemma 4.1 is given in [START_REF] Yu | Lipschitz stability in inverse parabolic problems by the Carleman estimate[END_REF].

We set

ϕ(x, t) = e λβ(x) t(T -t) (4.2) and η(x, t) = e 2λ β ∞ -e λβ(x) t(T -t) (4.3)
where λ > 0 and t ∈ (0, T ). Let us consider a boundary value problem for the parabolic operator:

       Lu(x, t) ≡ ∂ t u - n i,j=1 σ ij (x)∂ i ∂ j u = F, in Q, ∇u.ν = 0, on Σ. (4.4)
We recall that

H 2,1 (Q) = {u ∈ L 2 (Q)| ∂ t u, ∂ i u, ∂ i ∂ j u ∈ L 2 (Q), 1 ≤ i, j ≤ n}, (4.5)
Moreover, we set

γ = n i,j=1 σ ij C 2 ( Ω) , (4.6)
By β constructed in Lemma 4.1, we have the following Lemma:

Lemma 4.2. Let σ(x) ∈ C 2 (Ω) and the function ϕ, η be defined by (4.1)-(4.3). Then there exists a number λ 0 > 0 such that for an arbitrary λ ≥ λ 0 , we can choose a constant s 0 (λ) ≥ 0 satisfying: there exists a constant C = C(s 0 , λ 0 ) > 0 such that:

Q    (sϕ) p-1 (|∂ t u| 2 + n i,j=1 |∂ i ∂ j u| 2 ) + (sϕ) p+1 λ 2 |∇u| 2 + (sϕ) p+3 λ 4 |u| 2    e -2sη dxdt ≤ C Q (sϕ) p |Lu| 2 e -2sη dx dt + Ce C(λ)s Σ0 |∂ t u| 2 + |∇u| 2 + |u| 2 dSdt, (4.7)
for all s > s 0 , p = 0, 1, 2, and all u ∈ H 2,1 (Q). Here the constant C > 0 depends continuously on λ 0 but is independent of s, and λ 0 depends continuously on γ.

Proof 4.2. The proof of this lemma is similar to the proof given in [START_REF] Chae | Exact controllability for semilinear parabolic equations with neumann boundary conditions[END_REF]. We can further refer to [START_REF] Fursikov | Controllability of evolution equations[END_REF].

Global Carleman inequality for bidomain equations

We consider the solutions (v (n) , u

e , w (n) , z (n) ), n = 1, 2, to the following systems

                                 (1 + σ i σ -1 e )∂ t v (n) -div(σ i ∇v (n) ) + (1 + σ i σ -1 e )I ion (v (n) , w (n) , z (n) ) = A(x) • ∇u (n) e + I s i -σ i σ -1 e I s e in Q, -div(σ e ∇u (n) e ) = I ion (v (n) , w (n) , z (n) ) + ∂ t v (n) + I s e in Q, ∂ t w (n) = F (v (n) , w (n) ) in Q, ∂ t z (n) = G(v (n) , w (n) , z (n) ) in Q, σ i ∇v (n) .ν + σ i ∇u (n)
e .ν = 0 on Σ,

σ e ∇u (n) e .ν = 0 on Σ, (4.8)
and we consider the difference 2) , u e = u (1) e -u (2) e , w = w (1) -w (2) , z = z (1) -z (2) , ¯ = ¯ (1) -¯ (2) . (4.9)

v = v (1) -v (
Then, (v, u e , w, z) is solution to the following problem

                           (1 + σ i σ -1 e )∂ t v -div(σ i ∇v) = -(1 + σ i σ -1 e )h -A(x) • ∇u e in Q, -div(σ e ∇u e ) = h + ∂ t v in Q, ∂ t w = Φ in Q, ∂ t z = Ψ in Q, σ i ∇v.ν + σ i ∇u e .ν = 0 on Σ, σ e ∇u e .ν = 0 on Σ, (4.10)
Here h := = I ion (¯ (1) , v (1) , w (1) , z (1) ) -I ion (¯ (2) , v (2) , w (2) , z (2) ) = S (v (2) , w (2) , z (2) ) • + R(v (1) , w (1) , z (1) , v (2) , w (2) , z (2) ), (4.11) where the vector S = (S i ) 1≤i≤N is defined as follow:

S i (v (2) , w (2) , z (2)

) = y i (v (2) )(v (2) -E i (z (2) )) k j=1 (w (2) j ) pj,i , i = 1, . . . , N, (4.12)
and R(v (1) , w (1) , z (1) , v (2) , w (2) , z (2) ) = I ion (¯ (2) , v (1) , w (1) , z (1) ) -I ion (¯ (2) , v (2) , w (2) , z (2) ). (4.13)

In (4.11), represent the transpose of any matrix and the expression S (v (2) , w (2) , z (2) ) • is the Euclidian scalar products of the row vector S and the colon vector := ( i ) 1≤i≤N formed by the ionic parameters i . Finally, the functions Φ and Ψ are respectively given by: Φ = F (v (1) , w (1) ) -F (v (2) , w (2) ), (

Ψ = G(¯ (1) , v (1) , w (1) , z (1) ) -G(¯ (2) , v (2) , w (2) , z (2) ). (4.15)

In the next, we denote R(x, t) = R(v (1) , w (1) , z (1) , v (2) , w (2) , z (2) )(x, t), S(x, t) = S (v (2) , w (2) , z (2) )(x, t), (

and h(x, t) := S (x, t) • (x) + R(x, t). (4.17) Our Carleman estimate for the bidomain equations (4.10) is as follow.

Lemma 4.3. Let the Assumptions 2.1 be fulfilled and the function ϕ, η be defined by (4.1)-(4.3). Then there exists a number λ 0 > 0 such that for an arbitrary λ ≥ λ 0 , we can choose s 0 (λ) ≥ 0 such that for each s ≥ s 0 (λ), a solution (v, u e ) ∈ H 2,1 (Q) of the problem (4.10) satisfies the following inequality:

Q    sϕ(|∂ t v| 2 + n i,j=1 |∂ i ∂ j v| 2 ) + (sϕ)

Stability estimate of conductances parameters

This section is devoted to proof Theorem 3.1. Since the proof technically looks very awkward firstly we demonstrate some preliminary lemma. Lemma 5.1. Under the same hypothesis of Lemma 4.3, we have

Q    ∂ 2 t v 2 + n i,j=1 |∂ i ∂ j ∂ t v| 2 + (sϕ) 2 λ 2 |∇∂ t v| 2 + (sϕ) 4 λ 4 |∂ t v| 2 + (sϕ)(|∂ t v| 2 + n i,j=1 |∂ i ∂ j v| 2 ) +(sϕ) 3 λ 2 |∇v| 2 + (sϕ) 5 λ 4 |v| 2 + (sϕ) -1 n i,j=1 |∂ i ∂ j ∂ t u e | 2 +
(1 + σ i σ -1 e )∂ 2 t v -div(σ i ∇∂ t v) = -A(x) • ∇∂ t u e -(1 + σ i σ -1 e )∂ t h. (5.2)
Applying Lemma 4.2 with p = 1 to the solution ∂ t v of (5.2), we have the estimate

Q    ∂ 2 t v 2 + n i,j=1 |∂ i ∂ j ∂ t v| 2 + (sϕ) 2 λ 2 |∇∂ t v| 2 + (sϕ) 4 λ 4 |∂ t v| 2    e -2sη dx dt ≤ C Q (sϕ)(|∇∂ t u e | 2 + |∂ t h| 2 )e -2sη dx dt + Ce C(λ)s Σ0 ∂ 2 t v 2 + |∇∂ t v| 2 + |∂ t v| 2 dSdt. (5.3)
Taking the time derivative to the second equation of (4.10), we get

-div(σ e ∇∂ t u e ) = ∂ 2 t v + ∂ t h (5.4)
Applying Lemma 4.2 with p = 0 to the solution ∂ t u e of (5.4), we have

Q    (sϕ) -1 n i,j=1 |∂ i ∂ j ∂ t u e | 2 + (sϕ)λ 2 |∇∂ t u e | 2 + (sϕ) 3 λ 4 |∂ t u e | 2    e -2sη dx dt ≤ C Q ( ∂ 2 t v 2 + |∂ t h| 2 )e -2sη dx dt + Ce C(λ)s Σ0 |∇∂ t u e | 2 + |∂ t u e | 2 dSdt. (5.5)
Using ϕ ≥ 4 T 2 , (2C) × (5.3) + (5.5), and λ sufficiently large, we obtain

Q    ∂ 2 t v 2 + n i,j=1 |∂ i ∂ j ∂ t v| 2 + (sϕ) 2 λ 2 |∇∂ t v| 2 + (sϕ) 4 λ 4 |∂ t v| 2 + (sϕ) -1 n i,j=1 |∂ i ∂ j ∂ t u e | 2 +(sϕ)λ 2 |∇∂ t u e | 2 + (sϕ) 3 λ 4 |∂ t u e | 2 e -2sη dx dt ≤ C Q (sϕ) |∂ t h| 2 e -2sη dx dt + C Σ0 ∂ 2 t v 2 + |∇∂ t v| 2 + |∂ t v| 2 + |∇∂ t u e | 2 + |∂ t u e | 2 dSdt. (5.6)
Summing (4.18), (5.6) and for λ sufficiently large, we deduce the Lemma 5.1.

Lemma 5.2. We assume that t 0 = T /2 without loss of generality by changing the scale of t. We get Ω e -2sη(x,t0) S (x, t 0 ) • (x)

2 dx ≤ C Q s -1 | (x)| 2 e -2sη dx dt + v(•, t 0 ) 2 H 2 (Ω) + u e (•, t 0 ) 2 H 2 (Ω) + w(•, t 0 ) 2 L 2 (Ω) k + z(•, t 0 ) 2 L 2 (Ω) m + v 2 H 2 (0,T ;H 1 (Γ0)) + u e 2
H 1 (0,T ;H 1 (Γ0)) . (5.7)

Proof 5.2. We evaluate the first equation of (4.10) at a fixed time t 0 , we have

(1 + σ i σ -1 e )∂ t v(x, t 0 ) -div(σ i ∇v(x, t 0 )) + (1 + σ i σ -1 e )h(x, t 0 ) = -A(x)
• ∇u e (x, t 0 ).

(5.8)

Then, we integrate over Ω the square of (5.8) with the weight function e -2sη(x,t0) , we obtain

Ω e -2sη(x,t0) S (x, t 0 ) • ¯ (x) 2 dx ≤ C Ω e -2sη(x,t0) |∂ t v(x, t 0 )| 2 dx + Ω e -2sη(x,t0) |div(σ i ∇v(x, t 0 ))| 2 dx + Ω e -2sη(x,t0) |∇u e (x, t 0 )| 2 dx + Ω e -2sη(x,t0) |R(x, t 0 )| 2 dx ≤ C Ω e -2sη(x,t0) (|∂ t v(x, t 0 )| 2 + |∇u e (x, t 0 )| 2 + |R(x, t 0 )| 2 )dx + v(•, t 0 ) 2 H 2 (Ω) . (5.9)
Similarly, for the second equation of (4.10) at a fixed time t 0 , we get -div(σ e ∇u e (x, t 0 )) = h(x, t 0 ) + ∂ t v(x, t 0 ).

(5.10)

Then Ω e -2sη(x,t0) S (x, t 0 ) • ¯ (x) 2 dx ≤ C Ω e -2sη(x,t0) |∂ t v(x, t 0 )| 2 dx + Ω e -2sη(x,t0) |div(σ e ∇u e (x, t 0 ))| 2 dx + Ω e -2sη(x,t0) |∇u e (x, t 0 )| 2 dx + Ω e -2sη(x,t0) |R(x, t 0 )| 2 dx . (5.11)
Combining (5.9) and (5.11), we have

Ω e -2sη(x,t0) S (x, t 0 ) • ¯ (x) 2 dx ≤ C Ω e -2sη(x,t0) |∂ t v(x, t 0 )| 2 dx + v(•, t 0 ) 2 H 2 (Ω) + u e (•, t 0 ) 2 H 2 (Ω) + w(•, t 0 ) 2 L 2 (Ω) k + z(•, t 0 ) 2 L 2 (Ω) m . (5.12)
since

Ω e -2sη(x,t0) |R(x, t 0 )| 2 dx ≤ C v(•, t 0 ) 2 L 2 (Ω) + w(•, t 0 ) 2 L 2 (Ω) k + z(•, t 0 ) 2 L 2 (Ω) m . (5.13)
For the first term in the RHS of (5.12), we have

Ω e -2sη(x,t0) |∂ t v(x, t 0 )| 2 dx = t0 0 Ω d dt (e -2sη(x,t) |∂ t v| 2 )dx dt = t0 0 Ω -2s∂ t η(x, t) |∂ t v| 2 + 2p∂ 2 t v e -2sη(x,t) dx dt ≤ C t0 0 Ω s 2 ϕ 2 |∂ t v| 2 + s -2 ϕ -1 ∂ 2 t v 2 e -2sη(x,t) dx dt. (5.14)
Then, we apply the Carleman inequality (5.1) satisfied by ∂ t v, using 4

T 2 ≤ ϕ ≤ ( T 2
) 2 ϕ 2 , we obtain for s and λ sufficiently large

Ω e -2sη(x,t0) |∂ t v(x, t 0 )| 2 dx ≤ C s 2 Q (sϕ |∂ t h| 2 + (sϕ) 2 |h|
2 )e -2sη dx dt

+ C Σ0 ∂ 2 t v 2 + |∂ t v| 2 + |∇∂ t v| 2 + |∇v| 2 + |v| 2 + |∇∂ t u e | 2 + |∂ t u e | 2 + |∇u e | 2 + |u e | 2 dSdt. (5.15)
Using the Lemma 5.1 given in [START_REF] Abidi | On the identification of multiple space dependent ionic parameters in cardiac electrophysiology modelling[END_REF], we have

|h| 2 ≤ C( S • (x) 2 + |v| 2 + |w| 2 + |z| 2 ), (5.16)
and

|∂ t h| 2 ≤ C( ∂ t S • (x) 2 + |∂ t v| 2 + |v| 2 + |∂ t w| 2 + |w| 2 + |∂ t z| 2 + |z| 2 ) ≤ C(| (x)| 2 + |∂ t v| 2 + |v| 2 + |w| 2 + |z| 2 ),
(5.17)

we complete the proof.

Proof of Theorem 3.1.

In order to prove (3.5), we use the following integro-differential equation

v = {div(σ i ∇)} -1 (div(σ i + σ e )∇u e )) := Bu e ,
We can see that the operator B is linear and uniformly continuous [START_REF] Bourgault | Existence and uniqueness of the solution for the bidomain model used in cardiac electrophysiology[END_REF]. Replacing the new expression of v in the right hand side of (5.7) and using the uniform continuity of B, the estimate (5.7) becomes

Ω e -2sη(x,t0) S (x, t 0 ) • (x) 2 dx ≤ C Q s -1 | (x)| 2 e -2sη dx dt + u e (•, t 0 ) 2 H 2 (Ω) + w(•, t 0 ) 2 L 2 (Ω) k + z(•, t 0 ) 2 L 2 (Ω) m + u e 2 H 2 (0,T ;H 1 (Γ0)) . (5.18)
Summing up the above estimate over = 1, . . . , N , we get that

Ω e -2sη(x,t0) |Λ(x) • ¯ (x)| 2 dx ≤ C N s -1 Ω e -2sη(x,t0) |¯ (x)| 2 dx + N =1 u e 2 H 2 (0,T ;H 1 (Γ0)) + u e (•, t 0 ) 2 H 2 (Ω) + w(•, t 0 ) 2 L 2 (Ω) k + z(•, t 0 ) 2 L 2 (Ω) m , (5.19)
where the

N × N real matrix Λ(x) = Λ(v (2) (x, t 0 ), w (2) (x, t 0 ), z (2)
(x, t 0 )), for x ∈ Ω. Notice that we have

Λ(x)ξ R N ≥ α 1 (x) |ξ| R N , x ∈ Ω, ξ ∈ R N , (5.20)
where (α j (x)) 1≤j≤N ⊂ R N + denotes the increasing sequence of the singular values of Λ(x), and |ξ| R N stands for the Euclidian norm of ξ. Moreover, by Sobolev embedding theorem (e.g., Thm. 5.4 in [START_REF] Adams | Sobolev spaces[END_REF], Cor. 9.1, p. 46, inVol. 1 of [START_REF] Lions | Non-homogeneous Boundary Value Problems and Applications[END_REF]), we see that

H 3 (Ω) ⊂ C 1 (Ω), hence v (2) (., t 0 ), w (2)
(., t 0 ) and log z (2) (., t 0 ), = 1, . . . , N , being and Theorem 3.1 follows directly from this and (5.19) by choosing s so large that CN s -1 < α 2 0 .

Discussion and conclusions

Personalizing a numerical model in cardiac electrophysiology simulation allows to assimilate the electrophysiology conditions of a patient, including anatomical information but also functional conditions. These functional conditions are generally observed in the electrical signals measured either invasively in the heart domain or non-invasively on the body surface. The underlying substrate behind the electrical changes could be related to macroscopic parameters like the conductivity distribution of the tissue, but also and most importantly to the microscopic parameters at the cell and sub-cell levels. In particular, the conductances of the ion channels play an important role in the modulation of the transmembrane potential and consequently the ECGs. In this paper, we addressed the problem of identifying ion channels conductances from a set of extracellular potentials measurements on the boundary of the heart. Our result says that it is possibly to uniquely identify the conductances parameters if the mathematical model governing the electrical activity satisfies some hypothesis. Our approach is based on a Carleman inequality for the bidomain reaction diffusion model coupled to a general form of ordinary differential equation system. This result is a step further in the parameters identification problem compared to the existing results [START_REF] Abidi | On the identification of multiple space dependent ionic parameters in cardiac electrophysiology modelling[END_REF][START_REF] Ainseba | Stability of conductivities in an inverse problem in the reaction-diffusion system in electrocardiology[END_REF][START_REF] Boulakia | Stability estimates for some parameters of a reaction-diffusion equation coupled with an ode[END_REF][START_REF] Lassoued | Stability results for the parameter identification inverse problem in cardiac electrophysiology[END_REF][START_REF] Wu | Carleman estimate for a linearized bidomain model in electrocardiology and its applications[END_REF]. There are two novelties in our result compared to the results in those paper. First the cited papers deal with simpler models either by simplifying the ionic model or by simplifying the anisotropy of the bidomain model to reach the condition of the monodomain model approximation. In our case, we only remove the anisotropy but we keep ratio between intra and extracellular conductivities space dependent and we consider physiologically detailed ionic model. Second, we consider that the observations are given on a part of the accessible boundary of the domain which is more realistic compared to observations on a sub-domain of the heart.

Acknowledgements. This work has been supported by EPICARD cooperative research program, funded by INRIA international laboratory LIRIMA. The LAMSIN researcher's work is supported on a regular basis by the Tunisian Ministry of Higher Education, Scientific Research and Technology. This study received also financial support from the French Government as part of the 'Investments of the future' program managed by the National Research Agency (ANR), Grant reference ANR-10-IAHU-04.

Appendix

Proof of Proposition 2.4. First, we propose a parabolic-elliptic formulation of (2.1) given as follow:

                                 ∂ t v -div(σ i ∇v) + I ion (¯ , v, w, z) = div(σ i ∇u e ) + I s i in Q, -div(σ i ∇v + (σ i + σ e )∇u e ) = I s i + I s e in Q, ∂ t w = F (v, w) in Q, ∂ t z = G(¯ , v, w, z) in Q, σ i ∇v.ν + σ i ∇u e .ν = 0 on Σ, σ i ∇v.ν + (σ i + σ e )∇u e .ν = 0 on Σ, v(x, 0) = v 0 (x), u e (x, 0) = u e,0 (x), w(x, 0) = w 0 (x), z(x, 0) = z 0 (x) in Ω. (7.1)
Taking the time derivative of equations system (7.1), we obtain

                           ∂ tt v -div(σ i ∇∂ t v) + ∂ t I ion = div(σ i ∇∂ t u e ) + ∂ t I s i in Q, -div(σ i ∇∂ t v + (σ i + σ e )∇∂ t u e ) = ∂ t I s i + ∂ t I s e in Q, ∂ tt w = ∂ t F (v, w) in Q, ∂ tt z = ∂ t G(¯ , v, w, z) in Q, σ i ∇∂ t v.ν + σ i ∇∂ t u e .ν = 0 on Σ, σ i ∇∂ t v.ν + (σ i + σ e )∇∂ t u e .ν = 0 on Σ, (7.2)
with initial conditions

∂ t v(t = 0) = div(σ i ∇v 0 ) + div(σ i ∇u e (t = 0)) + I s i (t = 0) -I ion (t = 0) in Ω, ∂ t u e (t = 0) = ∂ t u e,0
in Ω,

∂ t w(t = 0) = F (v 0 , w 0 ) in Ω, ∂ t z(t = 0) = G(¯ , v 0 , w 0 , z 0 ) in Ω. (7.3) Using the hypothesis v 0 ∈ H 2 (Ω), u e (t = 0) ∈ H 2 (Ω), w 0 ∈ L 2 (Ω) k and z 0 ∈ L 2 (Ω) m , we deduce that div(σ i ∇v 0 ) ∈ L 2 (Ω), div(σ i ∇u e (t = 0)) ∈ L 2 (Ω), I ion (t = 0) ∈ L 2 (Ω), (7.4)
and then ∂ t v(t = 0) ∈ L 2 (Ω) since I s i,e satisfy the hypothesis (2.14). Also we deduce that

∂ t w(t = 0) ∈ L 2 (Ω) k , and ∂ t z(t = 0) ∈ L 2 (Ω) m . (7.5)
We integrate over Ω the sum of the first equation of (7.2) multiplied by ∂ t v, the 2 nd , 3 rd and the 4 th one multiplied by ∂ t u e , ∂ t w and ∂ t z, respectively, we obtain

1 2 d dt ( ∂ t v 2 L 2 (Ω) + ∂ t w 2 L 2 (Ω) k + ∂ t z 2 L 2 (Ω) m )+ < σ i ∇∂ t (v + u e ),
Ω ∂ t v∂ t I s i dx ≤ 1 2 ( ∂ t v 2 L 2 (Ω) + ∂ t I s i 2 L 2 (Ω) ) (7.9
)

Ω ∂ t u e ∂ t (I s i + I s e )dx ≤ ∂ t u e 2 L 2 (Ω) + C ( ∂ t I s i 2 L 2 (Ω) + ∂ t I s e 2 L 2 (Ω) ), (7.10)
Computing the time derivative of the ionic curren I ion , we get

∂ t I ion = N i=1 ¯ i ∂ t v y i (v) k j=1 w pj,i j (v -E i (z)) + N i=1 ¯ i y i (v)( k l=1 p l,i w p l,i -1 l ∂ t w l j =l w pj,i j )(v -E i (z)) + N i=1 ¯ i y i (v)∂ t v k j=1 w pj,i j + N i=1 ¯ i ∂ t z i z i y i (v) k j=1 w pj,i j . (7.11) Thanks to Lemma 2.3, we have v ∈ L ∞ (Q), log z ∈ L ∞ (Q) m , w ∈ [0, 1]
k and y i is C ∞ . Applying Cauchy-Schwarz inequality, we deduce from (7.9) and (7.11) that

Ω ∂ t v∂ t I s i dx + Ω ∂ t v∂ t I ion dx ≤ C( ∂ t v 2 L 2 (Ω) + ∂ t w 2 L 2 (Ω) k + ∂ t z 2 L 2 (Ω) m + ∂ t I s i 2 L 2 (Ω) ). (7.12)
We recall that the function F j given by (2.5) is C 2 (R 2 ) for j = 1, . . . , k, and then we can write

∂ t F j (v, w j ) = ∂ t v∂ 1 F j + ∂ t w j ∂ 2 F j , (7.13)
where ∂ l is the partial derivative with respect the l th variable, l = 1, 2. Since α j , β j are C ∞ (R), v ∈ L ∞ (Q) and w j ∈ [0, 1], j = 1, . . . , k, we deduce that

∂ 1 F j (v, w j ) = α (v)(1 -w j ) -β (v)w j , and ∂ 2 F j (v, w j ) = -α(v) -β(v) , j = 1, . . . , k, (7.14)
are bounded. Then, by Cauchy-Schwarz inequality and Young inequality, we have

Ω ∂ t w.∂ t F dx ≤ C( ∂ t v 2 L 2 (Ω) + ∂ t w 2 L 2 (Ω) k ). (7.15)
Similarly for the variable z, for i = 1, . . . , m, we can write

∂ t G i (¯ , v, w, z) = -∂ t v∂ 2 J i + k j=1 (∂ t w j ∂ j+2 J i ) + ∂ t log z i ∂ k+3 J i +∂ t v∂ 2 H i + k j=1 ∂ t w j ∂ j+2 H i + m j=1 ∂ t log z j ∂ k+2+j H i . (7.16)
By hypothesis (2.8) and (2.9), we deduce that ∂ l J i , l = 1, . . . , k + 3, and

∂ l H i , l = 1, . . . , m + k + 2, are bounded since v ∈ L ∞ (Q), w ∈ [0, 1] k , z ∈ L ∞ (Q) m . Then, we have Ω ∂ t z.∂ t Gdx ≤ C( ∂ t v 2 L 2 (Ω) + ∂ t w 2 L 2 (Ω) k + ∂ t z 2 L 2 (Ω) m ). ( 7
d dt ( ∂ t v 2 L 2 (Ω) + ∂ t w 2 L 2 (Ω) k + ∂ t z 2 L 2 (Ω) m ) + σ i ∇∂ t v + ∇∂ t u e 2 L 2 (Ω) + 1 2 σ e ∇∂ t u e 2 L 2 (Ω) ≤ C( ∂ t v 2 L 2 (Ω) + ∂ t w 2 L 2 (Ω) k + ∂ t z 2 L 2 (Ω) m + ∂ t I s i 2 L 2 (Ω) + ∂ t I s e 2 L 2 (Ω) ). (7.18)
Applying Gronwall Lemma, we deduce

∂ t v 2 L 2 (Ω) + ∂ t w 2 L 2 (Ω) k + ∂ t z 2 L 2 (Ω) m +σ i ∇∂ t v+∇∂ t u e 2 L 2 (Ω) + 1 2 σ e ∇∂ t u e 2 L 2 (Ω) ≤ C 0 , ∀t ∈ [0, T ], (7.19)
where

C 0 = C( v 0 H 2 (Ω) , w 0 L 2 (Ω) k , z 0 L 2 (Ω) m , I s i,e H 1 (0,T ;L 2 (Ω)) ) > 0. (7.20) Thus, v, u e ∈ W 1,∞ (0, T ; H 1 (Ω)), w ∈ W 1,∞ (0, T ; L 2 (Ω)) k , and z ∈ W 1,∞ (0, T ; L 2 (Ω)) m . (7.21)
Now, for t ∈ (0, T ), we intergrate (7.18) over (0, t), we obtain

∂ t v 2 L 2 (Ω) + ∂ t w 2 L 2 (Ω) k + ∂ t z 2 L 2 (Ω) m + σ i t 0 ∇∂ t v + ∇∂ t u e 2 L 2 (Ω) dτ + 1 2 σ e t 0 ∇∂ t u e 2 L 2 (Ω) dτ ≤ 2C 0 , (7.22)
and then, we conclude v, u e ∈ H 1 (0, T ; H 1 (Ω)). (7.23) On the other hand, from the system (2.10), we have div(σ e ∇u e ) = -∂ t v -I ion -I s e ∈ L 2 (0, T ; L 2 (Ω)), ⇒ u e ∈ L 2 (0, T ; H 2 (Ω)). (7.24) and div(σ i ∇v) = ∂ t v + I ion -div(σ i ∇u e ) -I s i ∈ L 2 (0, T ; L 2 (Ω)), ⇒ v ∈ L 2 (0, T ; H 2 (Ω)). (7.25) Let us now prove the regularities (2.21). Deriving equations (2.5) and (2.7) over the space variable x, we obtain

∂ t (∂ x w j ) = ∂ x v∂ 1 F j + ∂ x w j ∂ 2 F j , j = 1, . . . , k, (7.26
) and

∂ t (∂ x z i ) = N j=1 ∂ x j ∂ j G i + ∂ x v∂ N +1 G i + k j=1 ∂ x w j ∂ j+N +1 G i + m j=1 ∂ x z j ∂ k+N +1+j G i , i = 1, . . . , m. (7.27)
Without loss of generality, ∂ x is the space derivative over one direction (here could be the first, the second or the third dimension of the space R 3 ).

According to Lemma 2.3 and the hypothesis (2.7)-(2.9), there exist a constant C depending on T , such that

∂ 1 F j 2 L 2 (Q) + ∂ 2 F j 2 L 2 (Q) ≤ C, j = 1, . . . , k, (7.28)
and

N j=1 ∂ x j ∂ j G i 2 L 2 (Q) + ∂ N +1 G i 2 L 2 (Q) + k j=1 ∂ j+N +1 G i 2 L 2 (Q) + m j=1 ∂ k+N +1+j G i 2 L 2 (Q) ≤ C, i = 1, . . . , m.
(7.29) Multiplying the equation (7.26) by (∂ x w j ), j = 1, . . . , k, the equation (7.27) by (∂ x z i ), i = 1, . . . , m, and integrating the sum over Ω, and applying the Young inequality of the right hand side, we obtain

d dt ( ∂ x w 2 L 2 (Ω) k + ∂ x z 2 L 2 (Ω) m ) ≤ C( ∂ x v 2 L 2 (Ω) + ∂ x w 2 L 2 (Ω) k + ∂ x z 2 L 2 (Ω) m ). ( 7

.30)

Thanks to the regularity (7.23), we have

∂ x v L 2 (Ω) ∈ C 0 [0, T ].
Applying Gronwall lemma and using hypothesis (2.20), we obtain

∂ x w 2 L 2 (Ω) k + ∂ x z 2 L 2 (Ω) m ≤ C( ∂ x w 0 L 2 (Ω) k , ∂ x z 0 L 2 (Ω) m ), ∀t ∈ [0, T ]. (7.31)
Thus, w ∈ L ∞ (0, T ; H 1 (Ω)) k , and z ∈ L ∞ (0, T ; H 1 (Ω)) m . (7.32) On the other hand, from (7.26) and (7.27), we deduce

w ∈ W 1,∞ (0, T ; H 1 (Ω)) k , and z ∈ W 1,∞ (0, T ; H 1 (Ω)) m . ( 7.33)
• Let us now consider the initial conditions:

• ∂ tt v(0) = div(σ i ∇∂ t v(0)) + div(σ i ∇∂ t u e (0)) + ∂ t I s i (0) -∂ t I ion (0), • ∂ tt w(0) = ∂ t F (v(0), w(0)), • ∂ tt z(0) = ∂ t G( ¯ , v(0), w(0), z(0)). (7.34)
We start by showing that ∂ t v(t = 0) ∈ H 2 (Ω). Using hypothesis (2.14), v 0 ∈ H 4 (Ω), and u e (t = 0) ∈ H 4 (Ω), we have div(σ i ∇v 0 ) ∈ H 2 (Ω), and div(σ i ∇u e (t = 0)) ∈ H 2 (Ω). We also have I s i (t = 0) ∈ H 2 (Ω). It's obvious that

I ion (t = 0) = N i=1 ¯ i y i (v 0 ) k j=1 w pj,i 0,j (v 0 -E i (z 0 )) ∈ L 2 (Ω). (7.35)
Since ¯ is bounded in H 3 (Ω) N , there exists a constant C > 0 such that

|∂ x I ion (t = 0)| 2 ≤ C(|∂ x v 0 | 2 + |∂ x w 0 | 2 + |∂ x z 0 | 2 ), |∂ xy I ion (t = 0)| 2 ≤ C(|∂ x v 0 | 2 + |∂ x w 0 | 2 + |∂ x z 0 | 2 + |∂ y v 0 | 2 + |∂ y w 0 | 2 + |∂ y z 0 | 2 + |∂ xy v 0 | 2 + |∂ xy w 0 | 2 + |∂ xy z 0 | 2 ). (7.36) So, from hypothesis v 0 ∈ H 4 (Ω), w 0 ∈ H 2 (Ω) k , z 0 ∈ H 2 (Ω) m ,
v(t = 0) ∈ L 2 (Ω), ∂ tt w(t = 0) ∈ L 2 (Ω) k , and ∂ tt z(t = 0) ∈ L 2 (Ω) m .
Now, we take the second derivative with respect to time of the system (7.1), we multiply the first (respectively the 2 nd , 3 rd and the 4 th ) equation by ∂ tt v (respectively by ∂ tt u e , ∂ tt w, ∂ tt z), we obtain

• ∂ ttt v∂ tt v -∂ tt vdiv(σ i ∇∂ tt v) = ∂ tt vdiv(σ i ∇∂ tt u e ) + ∂ tt v∂ tt I s i -∂ tt v∂ tt I ion , • -∂ tt u e div(σ i ∇∂ tt v) -∂ tt u e div((σ i + σ e )∇∂ tt u e ) = ∂ tt u e (∂ tt I s i + ∂ tt I s e ), • ∂ ttt w j ∂ tt w j = ∂ tt F j (v, w j )∂ tt w j , j = 1, • • • , k, • ∂ ttt z i ∂ tt z i = ∂ tt G i ( ¯ , v, w, z)∂ tt z i , i = 1, • • • , m.
(7.37) Deriving (7.11) with respect to time, integrating over Ω the sum of the first and second equation of (7.37), using Cauchy Schwarz inequality, we get 1 2

d dt ∂ tt v 2 L 2 (Ω) + 1 2 < σ i ∇∂ tt (v + u e ), ∇∂ tt (v + u e ) > + 1 2 < σ e ∇∂ tt u e , ∇∂ tt u e > ≤ C 1 ( ∂ tt v 2 L 2 (Ω) + ∂ tt u e 2 L 2 (Ω) ) + C 2 ∂ tt w 2 L 2 (Ω) k + C 3 ∂ tt z 2 L 2 (Ω) m +C 4 ( ∂ t v 2 L 2 (Ω) + ∂ t w 2 L 2 (Ω) k + ∂ t z 2 L 2 (Ω) m + ∂ tt I s i 2 L 2 (Ω) + ∂ tt I s e 2 L 2 (Ω) ). (7.38)
Deriving in time (7.13) and (7.16), integrating over Ω the 3 rd and the 4 th equation of (7.37), using Cauchy Schwarz inequality, for j = 1, . . . , k and i = 1, . . . , m, we get We obtain v, u e ∈ H 2 (0, T ; H 1 (Ω)), w ∈ H 2 (0, T ; L 2 (Ω)) k , and z ∈ H 2 (0, T ; L 2 (Ω)) m . (7.43) On the other hand, from the system (2.10), we have div(σ e ∇∂ t u e ) = -∂ tt v -∂ t I ion -∂ t I s e ∈ L 2 (0, T ; L 2 (Ω)), ⇒ u e ∈ H 1 (0, T ; H 2 (Ω)). (7.44) and div(σ i ∇∂ t v) = ∂ tt v + ∂ t I ion -div(σ i ∇∂ t u e ) -∂ t I s i ∈ L 2 (0, T ; L 2 (Ω)), ⇒ v ∈ H 1 (0, T ; H 2 (Ω)). Let us now prove the regularities (2.22). Getting second derivatives of equations (2.5) and (2.7) over the space variable x, multiplying both equations by ∂ xx w j and ∂ xx z i respectively and using that fact that ¯ is bounded in H 3 (Ω) N , integrating the sum over Ω and applying the Young inequality of the right hand side, we obtain

• 1 2 d dt ∂ tt w j 2 L 2 (Ω) + C α,β ∂ tt w j 2 L 2 (Ω) ≤ C 5 ∂ tt v 2 L 2 (Ω) + C 6 ∂ tt w j 2 L 2 (Ω) +C 7 ( ∂ t v 2 L 2 (Ω) + ∂ t w j 2 L 2 (Ω) ), • 1 2 d dt ∂ tt z i 2 L 2 (Ω) + C ∂ tt z i 2 L 2 (Ω) ≤ C 8 ∂ tt v 2 L 2 (Ω) + C 9 ∂ tt w 2 L 2 (Ω) k + C 10 ∂ tt z i 2 L 2 (
+ ∂ tt w 2 L 2 (Ω) k + ∂ tt z 2 L 2 (Ω) m dτ ≤ C ∂ tt I s i,e 2
d dt ( ∂ xx w 2 L 2 (Ω) k + ∂ xx z 2 L 2 (Ω) m ) ≤ C 1 ( ∂ xx w 2 L 2 (Ω) k + ∂ xx z 2 L 2 (Ω) m ) + C 2 , (7.48)
where C 1 and C 2 are two non negative constants. Applying Gronwall Lemma, we obtain

∂ xx w 2 L 2 (Ω) k + ∂ xx z 2 L 2 (Ω) m ≤ C( ∂ xx w 0 L 2 (Ω) k , ∂ xx z 0 L 2 (Ω) m ). (7.49)
Thus, w ∈ L ∞ (0, T ; H 2 (Ω)) k , and z ∈ L ∞ (0, T ; H 2 (Ω)) m . (7.50) Similarly, computing the third space derivatives of (2.5) and (2.7) , multiplying both equations by ∂ xxx w j and ∂ xxx z i respectively and using that fact that ¯ is bounded in H 3 (Ω) N , integrating the sum over Ω and applying the Young inequality of the right hand side, we obtain

d dt ( ∂ xxx w 2 L 2 (Ω) k + ∂ xxx z 2 L 2 (Ω) m ) ≤ C 3 ( ∂ xxx w 2 L 2 (Ω) k + ∂ xxx z 2 L 2 (Ω) m ) + C 4 . (7.51)
since v ∈ L ∞ (0, T ; H 3 (Ω)), where C 3 and C 4 are two non negative constants. Using w 0 ∈ H 3 (Ω) k , z 0 ∈ H 3 (Ω) m and applying Gronwall Lemma, we have

∂ xxx w 2 L 2 (Ω) k + ∂ xxx z 2 L 2 (Ω) m ≤ C( ∂ xxx w 0 L 2 (Ω) k , ∂ xxx z 0 L 2 (Ω) m ). ( 7

Lemma 4 . 1 .

 41 Let Γ 0 = ∅ ⊂ ∂Ω be an arbitrary relatively open subset. Then there exists a function β ∈ C 2 (Ω) such that β(x) > 0 ∀x ∈ Ω, |∇β(x)| > 0 ∀x ∈ Ω, and σ i,e (x)∇β(x) • ν(x) ≤ 0, x ∈ ∂Ω \ Γ 0 . (4.1)

taken in C 1 (

 1 Ω) thanks to Proposition 2.4. Then, α 1 ∈ C 1 (Ω; R + ) from [30] [Thm 6.8 p. 122]. This combined with equation (3.4), yields α 0 := inf x∈Ω α 1 (x) > 0. As a Consequence and by (5.20), we have Ω e -2sη(x,t0) |Λ(x)¯ (x)| 2 dx ≥ α 2 0 Ω e -2sη(x,t0) |¯ (x)| 2 dx, (5.21)

L 2 ( 2 L 2

 222 0,T ;L 2 (Ω)) + C 15 C 0 T + ∂ tt v(t = 0) (Ω) + ∂ tt w(t = 0) 2 L 2 (Ω) k + ∂ tt z(t = 0) 2 L 2 (Ω) m .(7.42)

(7. 45 )

 45 Deriving(7.44) and (7.45) with respect to space, we obtain that∂ x div(σ e ∇∂ t u e ) = -∂ x ∂ tt v -∂ x ∂ t I ion -∂ x ∂ t I s e ∈ L 2 (0, T ; L 2 (Ω)), ⇒ u e ∈ H 1 (0, T ; H 3 (Ω)).

∂

  x div(σ i ∇∂ t v) = ∂ x ∂ tt v + ∂ x ∂ t I ion -∂ x div(σ i ∇∂ t u e ) -∂ x ∂ t I s i ∈ L 2 (0, T ; L 2 (Ω)), ⇒ v ∈ H 1 (0, T ; H 3 (Ω)).

. 52 )

 52 Thenw ∈ L ∞ (0, T ; H 3 (Ω)) k , and z ∈ L ∞ (0, T ; H 3 (Ω)) m . (7.53) Moreover, from the expressions of ∂ t (∂ xxx w) and ∂ t (∂ xxx z), we deducew ∈ W 1,∞ (0, T ; H 3 (Ω)) k ⊂ H 1 (0, T ; H 3 (Ω)) k → C 0 ([0, T ]; C 1 (Ω)) k , z ∈ W 1,∞ (0, T ; H 3 (Ω)) m ⊂ H 1 (0, T ; H 3 (Ω)) m → C 0 ([0, T ]; C 1 (Ω)) m . (7.54)The proof of Proposition 2.4 is finished.

  Similarly, applying estimate (4.7) with p = 1 to the solution u e of (4.10), we get

			 	n	|∂ i ∂ j u e |	2 + (sϕ) 2 λ 2 |∇u e | 2 + (sϕ) 4 λ 4 |u e | 2	 	e -2sη dx dt
		Q	 i,j=1				
								≤ C	(sϕ)(|∂ t v| 2 + |h| 2 )e -2sη dx dt + Ce C(λ)s	|∇u e | 2 + |u e | 2 dSdt. (4.20)
								Q	Σ0
	By	4 T 2 ≤ ϕ ≤ (	T 2	) 2 ϕ 2 , (2C) × (4.19) + (4.20), and λ sufficiently large, we obtain the Lemma 4.3.
			 	sϕ(|∂ t v| 2 +	n	|∂ i ∂ j v|	2 ) + (sϕ) 3 λ 2 |∇v| 2 + (sϕ) 5 λ 4 |v| 2	 	e -2sη dx dt
		Q					i,j=1	
								≤ C	(sϕ) 2 (|∇u e |	2 + |h|	2 )e -2sη dx dt + Ce C(λ)s	|∂ t v|	2 + |∇v| 2 + |v|	2 dSdt. (4.19)
								Q	Σ0

3 λ 2 |∇v| 2 + (sϕ) 5 λ 4 |v| 2 + n i,j=1

|∂ i ∂ j u e | 2 +(sϕ) 2 λ 2 |∇u e | 2 + (sϕ) 4 λ 4 |u e | 2 e -2sη dx dt ≤ C Q (sϕ) 2 |h| 2 e -2sη dx dt + C Σ0 |∂ t v| 2 + |∇v| 2 + |v| 2 + |∇u e | 2 + |u e | 2 dSdt, (4.18)

where the constant C depends continuously on λ 0 .

Proof 4.3. From the boundary conditions of the Problem (4.10), we can write σ i ∇v.ν = 0 on Σ. Then applying estimate (4.7) with p = 2 to the solution v of (4.10), we have the estimate

  ∇∂ t (v + u e ) > + < σ e ∇∂ t u e , ∇∂ t u e >= T ∇ϑdx. From Assumptions 2.1, we consider < σ i ∇∂ t (v + u e ), ∇∂ t (v + u e ) >≥ σ i ∇∂ t v + ∇∂ t u e ∇∂ t u e , ∇∂ t u e >≥ σ e ∇∂ t u e

	∂ t v∂ t I s i dx +	∂ t u e (∂ t I s i + ∂ t I s e )dx +	∂ t v∂ t I ion dx +	∂ t w.∂ t F dx +	∂ t z.∂ t Gdx. (7.6)
	Ω	Ω	Ω	Ω	Ω
	where < σ∇ϑ, ∇ϑ >=
				2 L 2 (Ω) ,	(7.7)
	and
		< σ e 2 L 2 (Ω) .	(7.8)
	Using Cauchy-Schwarz inequality, we have

Ω

(σ∇ϑ)

  we deduce that ∂ x I ion (t = 0) and ∂ xy I ion (t = 0) belong to L 2 (Ω). Thus I ion (t = 0) ∈ H 2 (Ω), and then ∂ t v(t = 0) ∈ H 2 (Ω). From (7.34), we get ∂ tt

  Ω) +C 11 ( ∂ t v 2 L 2 (Ω) + ∂ t w 2 L 2 (Ω) k + ∂ t z i , C 2 , . . . , C 11 , C α,β , C positive constants such that α j (v) + β j (v) ≥ C α,β and 1 z i ∂ k+3 J i ≥ C. Summing L 2 (Ω) + ∂ tt w 2 L 2 (Ω) k + ∂ tt z 2 L 2 (Ω) m ) + 1 2 σ i ∇∂ tt v + ∇∂ tt u e +C α,β ∂ tt w 2 L 2 (Ω) k + C ∂ tt z 2 L 2 (Ω) m ≤ C 12 ( ∂ tt v 2 L 2 (Ω) + ∂ tt u e 2 L 2 (Ω) ) + C 13 ∂ tt w 2 L 2 (Ω) k +C 14 ∂ tt z 2 L 2 (Ω) m + C 15 ( ∂ t v 2 L 2 (Ω) + ∂ t w 2 L 2 (Ω) k + ∂ t z 2 L 2 (Ω) m + ∂ tt I s , with C p is the Poincaré constant depending on Ω, ∂ tt v 2 L 2 (Ω) + ∂ tt w 2 L 2 (Ω) k + ∂ tt z 2 L 2 (Ω) m ) + C( ∇∂ tt v + ∇∂ tt u e 2 L 2 (Ω) + ∇∂ tt u e 2 L 2 (Ω) + ∂ tt w 2 L 2 (Ω) k + ∂ tt z 2 L 2 (Ω) m ) ≤ C 15 ( ∂ t v 2 L 2 (Ω) + ∂ t w 2 L 2 (Ω) k + ∂ t z 2 L 2 (Ω) m +∂ tt I s Using hypothesis (7.34) and inequality (7.19), for t ∈ [0, T ], we integrate over the time interval (0, t), we get∂ tt v 2 L 2 (Ω) + ∂ tt w 2 L 2 (Ω) k + ∂ tt z 2 L 2 (Ω) m + C

	where C 1 all equations in (7.38) and (7.39), we get
		1 2	d dt	( ∂ tt v 2					2 L 2 (Ω) +	1 2	σ e ∇∂ tt u e	2 L 2 (Ω)
										(7.40)
										i	2 L 2 (Ω) + ∂ tt I s e	2 L 2 (Ω) ).
	Applying Poincaré inequality by choosing C 12 = σ i C 13 = 1 4C p 1 2 C α,β , and C 14 = 1 2 C, we obtain
	1 2	d dt	( i	2 L 2 (Ω) + ∂ tt I s e	2 L 2 (Ω) ), (7.41)
	where C = min{	1 2	C α,β ,	1 2	C,	1 4	σ i ,	1 2	σ e }.
										(7.39)
										2 L 2 (Ω) ),

t 0 ∇∂ tt (v + u e ) 2 L 2 (Ω) + ∇∂ tt u e 2 L 2 (Ω)


Amine Lazrak

François Boudehenn

Sylvain Bonnot

Gilles Fraisse

Antoine Leconte

Philippe Papillon

Bernard Souyri

Development of a dynamic artificial neural network model of an absorption chiller and its experimental validation

Keywords: Thermal systems, Absorption chiller, Performance estimation, Dynamic modelling, Artificial neural networks, System testing

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

From more than ten years now, the development of air-conditioning offers more and more comfort and safety to people. Major part (approximately 99%) of commercial air-conditioners use electricity-powered compression cycles.

This system must face a paradox: the more air-conditioners installed in a city, the more heat is released to urban atmosphere, and the more ambient air temperature increases, decreasing the performance efficiency of airconditioners and increasing the cooling load of building [START_REF] Pons | Performance comparison of six solar-powered airconditioners operated in five places[END_REF]. In fine, peak electricity demand for cooling must be tripled [START_REF] Santamouris | On the impact of urban climate on the energy consumption of buildings[END_REF]. According to Pons et al. [START_REF] Pons | Performance comparison of six solar-powered airconditioners operated in five places[END_REF], part of the solution could be the use of thermally driven chillers powered by waste-heat or solar energy.

The major part of thermally driven chillers available on the market is absorption chiller. The basic physical process consists of at least two chemical components, one of them serving as the refrigerant and the other as the sorbent.

The operation of such systems is well documented and is not described here [START_REF] Henning | Solar assisted air conditioning of buildings -an overview[END_REF]. Main advantages of this technology are the continuous thermodynamic cycle and the high thermal coefficient of performance (COP th ), compared to adsorption chiller.

Absorption chillers are available on the market in a wide range of capacities and designed for different applications.

However, only very few systems are available in a range below 100 kW of cooling capacity. 75% of these systems are single effect absorption chillers [4]. Absorption chillers behavior is highly dependent on the climatic conditions and building quality (i.e. on the boundary conditions) [START_REF] Boudéhenn | Dynamical studies with a semi-virtual testing approach for characterization of small scale absorption chiller[END_REF]. This is why systems could have poor energy performance in some environments. For this reason it becomes necessary that users could have reliable information about the long-term performance of the system in the wanted boundary conditions. This requires reliable and faithful models of absorption chillers under reals conditions. Two types of models can be used to predict behavior and performances of absorption chillers: physical and empirical ones. Physical models describe all the thermodynamic cycle of the absorption chiller, heat exchangers performances and compute then heating and cooling capacity and temperatures of external fluids. On the contrary simplified models only focused on external fluids computing temperatures and power by a set of non-physical equations fitted with experimental data. Examples of those types of models are well studied and compared in steady states conditions in [START_REF] Boudéhenn | Comparison of Different Modeling Methods for a Single Effect Water-Lithium Bromide Absorption Chiller[END_REF] and [START_REF] Swider | A comparison of empirically based steady-state models for vapor-compression liquid chillers[END_REF].

Among simplified models, Artificial Neural Networks (ANNs) seems to be the most powerful mathematic tool to solve this modelling problem. In fact, it was shown that ANNs are universal function approximators [START_REF] Cybenko | Approximation by Superpositions of a Sigmoidal Function[END_REF], so they can be used to approximate the system function. ANNs were applied successfully to solve complex, non-linear, dynamic and multivariable problems. They tolerate errors, imprecisions and missing data too [START_REF] Kalogirou | Artificial neural networks in renewable in renewable energy systems application: a review[END_REF]. ANNs were extensively used during the last decade and have been especially used to solve prediction modelling problems in renewable energy thermal systems [START_REF] Yilmaz | Modeling of mechanical cooling system with variable cooling capacity by using artificial neural network[END_REF], [START_REF] Ahmad | Investigating the performance of a combined solar system with heat pump for houses[END_REF], [START_REF] Almonacid | Estimation of the energy of a PV generator using artificial neural network[END_REF], [START_REF] Vig | Neural network Modeling of thermal stratification in a solar DHW storage[END_REF], [START_REF] Kalogirou | Artificial neural networks for the performance prediction of large solar systems[END_REF], [START_REF] Kalogirou | Artificial neural networks and Genetic Algorithms in Energy Applications in Buildings[END_REF]. In the following are presented selected studies that focus on modelling absorption chiller systems:

• In [START_REF] Frey | Artificial Neural Network modelling of sorption chillers[END_REF] authors developed an ANN based model of an adsorption chiller. The ANN has 6 inputs (3 inlet temperatures and 3 flow rates) and 3 outputs (3 outlet temperatures). The model developed is able to predict with a high accuracy, error less than 2°C, the output temperatures. However, because it was trained using a quasi-steady state database, the ANN cannot predict the short-term behavior of the system;

• In [START_REF] Rosiek | Performance study of solar-assisted air-conditioning system provided with storage tanks using artificial neural networks[END_REF] (another similar study is presented in [START_REF] Rosiek | Modelling a solar-assisted air-conditioning system installed in CIESOL building using an artificial neural network[END_REF] as well) authors developed an ANN based model of a solar-driven absorption chiller. The model developed was able to predict with a low error both the COP th and the system cooling capacity. For this study only 5 inputs were relevant to model the whole system e.g. the inlet and outlet temperatures of the evaporator and generator, and the average temperature of the hot storage tank. The inputs of the modelling configuration used are not suitable to evaluate the performance of an absorption chiller when only inlet temperatures and flow rates are available;

• In [START_REF] Manohar | Modelling of steam fired double effect vapour absorption chiller using neural network[END_REF] authors modelled a double effect absorption chiller in steady state using neural networks. The model was validated using experimental data of about 250 samples. The used ANN predicted the performance of the absorption chiller quite accurately (R² greater than 0.99);

• In [START_REF] Labusa | Inverse neural network based control strategy for absorption chillers[END_REF] authors have developed a control system based on an inverse ANN model. The static ANN model, which is an analytical function, leads to a short computing time and as a consequence make this methodology suitable for the on-line control of absorption cooling systems.

Further studies where static ANNs were used to model or optimize absorption chillers could be consulted in [START_REF] Chow | Global optimization of absorption chiller system by genetic algorithm and neural network[END_REF], [START_REF] Congradac | Recognition of the importance of using artificial neural networks and genetic algorithm to optimize chiller operation[END_REF], [START_REF] Şencan | Performance of ammonia-water refrigeration systems using artificial neural networks[END_REF] and [START_REF] Şencan | Thermodynamic analysis of absorption systems using artificial neural network[END_REF].

The presented studies focus only on steady or quasi-steady state behavior of the system and the modelling configuration used in some of them cannot be helpful to predict the long term absorption chiller energy performance. Dynamic simulation plays a very important role in the description of the real performance of an energy conversion system, especially during the activation stage or part-load operation. Such a problem is extremely relevant for absorption chillers, where the high mass of the internal components and the accumulation of the fluids inside the vessels usually make the transient period longer than for mechanical compression chillers [START_REF] Evola | Proposal and validation of a general model for the transient simulation of single-stage LiBr/water absorption chillers[END_REF]. The carried literature shows that there is a need of a complete method to model absorption chillers in order to predict their energy performance. The current work presents a methodology to model absorption chillers in a dynamic way. Indeed, the proposed method does not assume any prior knowledge about the system to be modelled or its components. This makes the method generic and more relevant for a future tool to evaluate absorption chillers energy performance within the framework of a certification context. This methodology is based on a short dynamic test and the identification of a "black box" dynamic model of the system using ANN.

The absorption chiller

The absorption chiller used in this study is a single effect water-lithium bromide chiller manufactured by the German company EAW with a nominal cooling capacity of 15 kW. The nominal thermal COP specified by the manufacturer [START_REF] Eaw | Technische Beschreibung, Absorptionskälteanlagen WEGRACAL SE 15[END_REF] is 0.71. The hydraulic specifications of the chiller are presented in Table 1. In order to be representative, the characterization of the absorption chiller must be realized under operating conditions equivalent to those of a real system. For this, a semi virtual testing method is used. This method is based on the use of a solar cooling system model under TRNSYS, which defines at each time step the operating conditions (temperatures and flow rates) applied to the chiller. Then a dynamical test bench transforms the numerical information of the system model into real physical values at each time step.

The methodology used for the dynamical characterization of the absorption chiller is based on the dynamical study of solar combi-system (SCS) developed at INES. This protocol estimates the SCS performances on 12 days representative of the 12 months of a year [START_REF] Leconte | Global approach test improvement using a neural network model identification to characterise solar combisystem performances[END_REF]. Despite the fact that this methodology is difficult to transfer directly to solar cooling system, it is possible to test the solar cooling system on three representative days of the different conditions met during the summer period and corresponding to a day with a high cooling demand and a medium irradiation (day 1), a day with a high cooling demand and a high irradiation (day 2) and a day with a low cooling demand and a high irradiation (day 3). These three days are chosen between June the 1st and September the 30th in the Carpentras (south of France) weather data file, used for the test of the system. For each test day (8 hours), absorption chiller is tested and the rest of the system is simulated with TRNSYS. During the other period, the complete system, including the absorption chiller, is simulated. One advantage of this methodology is that it is able to define good initial conditions of the system (building, storage, etc.) at the beginning of the physical test of the chiller. Compiling the 3 days, is it possible to have representative experimental data of operating chiller conditions integrated into a solar cooling system.

The semi-virtual test bench consists mainly on a boiler room able to supply hot (54 kW at 180°C) and cold (150 kW at -10°C) water with two distribution loops, and on test modules of 25 or 50 kW able to reproduce the desired dynamical thermal loads within these range of temperature, Fig. 1. This test bench allows testing of various thermal systems for space heating, air conditioning and DHW preparation for different kind of applications (single or multifamily houses, small industry, small tertiary buildings). The operating conditions of each module are representative of all or part of a system and their behavior is controlled by a TRNSYS numerical model linked. With this method, a component or several components are characterized in a semi virtual environment, in fact by emulating some of the components of the system and the climatic conditions to which it is subjected. In the case of application, the climatic conditions, the solar loop, the heat rejection loop and the building are emulated. The TRNSYS absorption solar cooling reference system used include the following loops:

• the solar loop, used for the alimentation of the thermally driven chiller, comprised flat-plate collectors, a plate heat exchanger and a hot water storage. Using a parametric study, [START_REF] Boudéhenn | Dynamical studies with a semi-virtual testing approach for characterization of small scale absorption chiller[END_REF] have defined the following criteria for an optimal design of the solar loop of system : 4.5 m 2 /.kW & 55 l/kW;

• the heat rejection of the chiller is connected to a closed circuit cooling tower, designed at nominal values provided by the chiller manufacturer;

• the cold loop uses a heat exchanger like fan-coils able to transfer the cold water power provided by the chiller into cold air power introduced in the building. In order to adjust the cooling power provided by the chiller to the cooling demand of the building, the cooling power dissipated by the fan-coils into the building (P cool ) is calculated using the maximal cooling demand of the building (P SFH60 ), the nominal cooling power of the chiller (P N ), the effective cooling power of the fan coils (P fc ) and the effectiveness of the fan coils (η fc ), with the following equation:

( ) ( ) N fc fc SFH cool P P P P ⋅ ⋅ = η 60 eq 1
The climatic conditions and the building used for the simulation have been defined by the IEA-SHC Task 32 [START_REF] Heimrath | Project Report A2 of Subtask A: The Reference Heating System, the Template Solar System[END_REF]: Carpentras (FRANCE) and a 60 kWh.m -2 .year On the right picture of Figure 2, the three identified days show different profiles:

• during the first day, the chiller operates during 274 minutes with a decrease of the inlet generator temperature during the day;

• during the second day, the chiller operate during 471 minutes with first an increase of the inlet generator temperature from 80 to 93.7°C before decreasing;

• During the day 3, due to a low cooling demand, the chiller operates in an on/off mode with a frequency close to 30 minutes for each phase. 2. During the operating mode of the absorption chiller, the thermal losses to the environment of the absorption chiller can be ignored regarding the heat transfers between the system and its boundary conditions. For this reason the ambient temperature was not considered in the model input configuration.

The modelling approach

The ANN theory

ANN theory is clearly presented in [START_REF] Dreyfus | Neural networks methodology and applications[END_REF] and [START_REF] Norgaard | Neural networks for modelling and control of dynamic systems[END_REF]. ANNs are parametric analytical functions whose concept takes inspiration from the human central nervous system. A neuron, the basic element of an ANN, can compute values from a weighted summation of its inputs (eq. 2). The summation coefficients are called synaptic weights. The subscript denotes the neuron number. The neural operation is presented in eq. 1. The function is called the neural activation function (AF).

eq. 2

Inter-connected neurons constitute what is commonly called a neural network. There are several network architectures and each one is more suiTable for a specific problem than others. The most famous architecture for prediction and fitting problems is the class of the multi-layer perceptron (MLP). A MLP is a feed-forward network built of neurons, arranged in layers. It has an input layer, one or more hidden layers and an output layer. In Figure 3 a MLP, with inputs, neurons in the hidden layer and outputs, is represented. The th output of the network can be obtained using eq. 3.

eq. 3

Modelling a system using an ANN means building the appropriate network architecture (hidden layers, neurons, AFs) and identifying the different weights. The ANN learning or training is the process of determination of an ensemble of weights so that the underlying function approximates the real system function. In fact, the objective of the training process is to minimize, for a specific ANN structure, a cost function, with respect to weights, knowing a set of data (training data).

Usually, the resulting models fit with high accuracy the learning data set (especially when a lot of neurons are used e.g. a lots of parameters). However, the model prediction precision may be very poor for unseen data (e.g. overfitting). In some cases, this means that the resulting model has learned only characteristics of the training data set and not the system behavior. In other cases probably the optimizer has been stacked in a local minimum. The learning strategies are ways to deal with those optimization algorithms and the neural network structure in order to get a model that has good generalization ability e. g. accurate estimations for unseen data.

One must be aware that the most optimization algorithms commonly used are local methods. In fact, they are valid for a region near the minimum. This is why a good way of preventing the algorithms to be stacked on a local poor minimum is to train the ANN several times with different starting points (weights initialization).

Training eq. 4

Bayesian regularization modifies the objective function by adding an additional term: the sum of squares of the network weights (eq. 5) where is a vector of parameters to be determined (synaptic weights), and are two constant parameters calculated using a Bayesian regularization methodology [START_REF] Forsee | Gauss-Newton approximation to Bayesian learning[END_REF]. By constraining the size of weights the training process produce an ANN with good generalization ability [START_REF] Mackay | Proceedings of the International Joint Conference on Neural Networks[END_REF]. In fact, by keeping the weights small the ANN response will be smooth and so the over-fitting is supposed to be prevented. In this study the objective function minimization is done using the Levenberg-Marquardt algorithm.

eq. 5

Usually, training algorithms (even if coupled to an adequate training strategy) do not necessarily guarantee the production of efficient networks based on row data. This is why it is essential to carry out some data pre-processing before training. By normalizing or standardizing the input and target data vectors, the ANN training will be easier and faster, and all vectors will be equally taken into account during the learning process. Equation 6 will be used to pre-treat the training data in order to fall between and (normalization bounds); represents a vector (input or output of the model) of data over time.

eq. 6 and represent the minimum and maximum values of each variable.

Generally ANNs do not give good results in extrapolation. One way to reduce the effect of this problem is to select carefully the normalization bounds in eq. 6.

The system modelling process

Determination of the modelling input-output configuration is crucial to develop a reliable and exploitable model. The goal of this study is to develop a model able to predict the energy performance of the system (COP) and its cooling capacity. To do so, it is necessary to predict the outlet temperatures at the level of the evaporator and the generator.

These two unknown variables are then used as model outputs and the other ones which are control variables e. g. inlet temperatures (condenser, evaporator and generator) and the three flow rates (condenser, evaporator and generator) are used as model inputs. In order to make the ANN learn the dynamic of the system, delayed outputs must be used as inputs too (by feeding them back). This two outlet temperatures (evaporator and generator), physically, are affected by the outlet temperature of the condenser. This is why it is relevant to add it as output of the model in order to be used as a delayed input. The resulting input-output modelling configuration is presented in The absorption chiller was tested in a dynamic way as described in section 2.

To determine the optimal model that describes best the system dynamic behaviour, several ANN with different number of neurons (from 1 to 15) were trained. All the models contain one hidden layer and tangent hyperbolic or linear function as the AF. Three normalization bounds: ±1, ±0.3 and ±0.8 were investigated as well. As said in section 3.1 in order to prevent the minimization algorithm getting stuck in local minima each network is retrained 15 times using the Nguyen-Widrow algorithm developed in [START_REF] Nguyen | Improving the learning speed of 2-layer neural networks by choosing initial values of the adaptative weights[END_REF]. Because the true output (actual data) is available during the training of the network, the actual outputs are used instead of feeding back the predicted outputs, as shown in Figure . 4. This has the advantage to use an accurate delayed output as an input so the network will learn more efficiently.

Input delays ANN Input delays

Figure. 4 :Open loop architecture

The Bayesian information criterion (BIC) was chosen to select the most relevant network [START_REF] Dreyfus | Neural networks methodology and applications[END_REF], [START_REF] Qia | An investigation of model selection criteria for neural network time series forecasting[END_REF]. It is defined by the following: eq. 7

where is the number of model weights and the size of the learning data.

The BIC is calculated using only the training data and in a closed loop architecture (see Figure 5), in which only available data, e.g. inputs and the initial values of the outputs, are used to predict the outputs (long-term model simulation). The open loop architecture is used only during training. The BIC selects ANNs with a small size. This is advantageous because parsimonious ANNs have a better generalization power as stated in [START_REF] Khosravi | Quantifying uncertainties of neural network-based electricity price forecasts[END_REF].

Input delays ANN Feedback delays 294 The Figure 9 presents a comparison between predictions of the selected ANN and the experimental results for the outlet temperatures of each hydraulic loop. This first comparison is based on the identification data. The absolute relative instantaneous error is calculated equation 8, where is the estimate of the variable . This relative error is calculated only when the absorption chiller is operating e.g. when there is a flow rate in each hydraulic loop, otherwise the error is set to zero.

eq. 8

The ANN predictions are very close to the experimental data. Outlet instantaneous temperature errors are, in the whole, less than 6% for the generator and the condenser-absorber and less than 10% for the evaporator. The learning of the latter seems to be more difficult. The mean error values are presented on the

316

eq. 9

317

The mean error values on temperatures and thermal power are presented on the Table 6.

318 In general the ANN predicts with a good precision degree each temperature for both data sets "validation 1" and "validation 2". In fact, mean error values of the outlet temperature in each hydraulic loop are less than 5%. As expected, the prediction of the outlet evaporator temperature by the ANN is less accurate than the two other outlet temperature variables. However, the corresponding instantaneous errors are low, still less than 10%. The neural model seems reproducing with a very satisfactory precision degree the dynamic of the system, even in the on-off operating mode from 1100min to 1600min.

The figures show that the ANN performances on the "validation 1" data are better than on the "validation 2" data.

This can be explained by the fact that the "validation 1" data are close to the data used during the learning process.

Some oscillations or big differences between the ANN predictions and the test field data can be noticed in the beginning. During this period the system is off, no flow rate in the hydraulic loops, there is no heat transfer between the system and its boundary conditions except the thermal losses to the environment. In the modelling configuration of the system the ambient temperature was not taken into account, this explains why the ANN was not able to predict well the outlet temperatures in this period. However, even when the absorption chiller is off but during the test process the model gives a good estimation of the temperatures. This shows than the ANN could be able to model systems even when a potential influential variable was not included in the inputs configuration set.

The mean error values for the thermal powers predictions are about 9% which is bigger than the values found for the temperature predictions. This is a mathematical amplification which depends on the value of the inlet temperature when calculating the powers.

The absolute relative errors of the transferred energy are given in the Table 7. They were calculated for all the three hydraulic loops (driving loop, heat rejection loop, chilled water loop) of the absorption chiller. Similar results, errors between 0.14% and 2.53% were found in [START_REF] Frey | Artificial Neural Network modelling of sorption chillers[END_REF] but using an adsorption chiller and in steady and quasi-steady state operating conditions.

Conclusion

In the present paper, the results of the development of a methodology to model absorption chillers for building applications are presented. This can be applied to systems as found on the market and there is no need to create a physical model of the system which could be very difficult and time consuming. The dynamic ANN models developed are able to predict, with a good level of precision, outlet temperatures and the transferred energy of all the three hydraulic loops (driving loop, heat rejection loop, and chilled water loop). In fact, the ANN generalization ability makes possible to predict the system's behaviour in various conditions (by changing the loads), different from the ones used during the ANN learning process. Absolute relative prediction errors are almost within the range of 1.1-5% for outlet temperatures and 1-6.6% for the transferred energy. The proposed approach will be helpful in the context of energy performance guarantees. Moreover, ANN simulation takes only a few seconds. This is ideal for engineers and designers to compare different solutions and select the system most suitable for a given building and location.

Although in this study the neural model was used in the simulation mode which means that predictions are done for long-term, there is a good agreement between the neural model predictions and the measured data. Usually, short term predictions are more accurate. Therefore ANN models could be used in advanced control systems (predictive control) in order to asses effectively the absorption system operation.

Figure 1 :

 1 Figure 1: Picture of the EAW absorption chiller connected to the semi-virtual test bench The test bench is equipped with temperature sensors (uncertainty of 0.1K on T and 0.25K on the ∆T) at the inlet and outlet of the tested component and with mass flow rate sensors (uncertainty of 0.2% of the measured value) on each hydraulic loop.

  -1 single family house (SFH 60) with an air change rate of 0.4 volume per hour. Two controls are used, one on the solar heat production (comparison between the solar collector temperature and the storage tank temperature with two hysteresis of +5°C and of -2°C) and one on the internal temperature of the building (26°C with a hysteresis of -1°C). The Figure 2 left gives an overview of the inlet and outlet temperature profiles on the three hydraulic loops for the 3 days compiled and the Figure 2 right presents the variation of the inlet temperature profiles for the three considered sequences: identification and validation 1

Figure 2 :

 2 Figure 2: Operating temperatures on the three hydraulic loops (left) and inlet temperature profile variations for identification and validation test sequences (right)

  Figure xx right. The "Identification" and "Validation 1" temperature profiles are very close, contrary to the Validation 2 temperature profiles lower in terms of inlet temperature at the generator. The temperature and mass flow rate variations for each profile are described in Table2.

Figure. 3 :

 3 Figure. 3: Formal representation of a neuron (left), example of a neural network MLP with one hidden layer (right)

  using Cross-Validation and Early Stopping are examples of the most famous learning strategies that are commonly used to prevent over-fitting. The both methods are based on a data set to do validation which must be taken from data available for training. Physical system tests are expensive and usually time consuming. The short time of the system test, which is long, in this study, of only one day, restricts the amount of available data for training. Therefore it is necessary to use a learning strategy that can use a restricted data set without compromising the generalization ability of the model. For this reason it is relevant to use Bayesian regularization method for the learning process. Typically, training aims to minimize the sum of squared errors ( ) see (eq. 4) with the target data at the sequence time and the ANN output at the same time.

Figure

  Figure. 4.

Figure. 5 : 6 nFigure. 6 :Figure. 7 :Figure. 8 :

 56678 Figure. 5 : Closed loop architectureThe whole training and selection process was developed in MATLAB R2012b and is represented in Figure6

Figure 9 :

 9 Figure 9: Outlet temperature comparisons between numerical and experimental results for each hydraulic loops

Figure 10 :Figure 11 :

 1011 Figure 10: Outlet temperature comparisons between numerical and experimental results for each hydraulic loops

Table 1 : Hydraulic specifications of the 15 kW cooling capacity EAW absorption chiller [26]

 1

		Capacities Inlet/outlet nominal	Inlet min/max	Flowrates	Pressure
	EAW Wegracal SE 15		temperatures	temperatures		drops
		[kW]	[°C]	[°C]	[m 3 /h]	[mbars]
	Evaporator	15	17 / 11	6 / -	1.9	400
	Generator	21	90 / 80	70 / 95	1.8	400
	Absorber-Condenser	35	30 / 36	25 / 40	5	900

Table 2 : Temperature and mass flow rate variations for the different day for each profile

 2

			Identification			Validation 1			Validation 2
			Average	Min.	Max.	Average	Min.	Max.	Average	Min.	Max.
	T E, in	[°C]	22.2	11.4	45.4	22.0	10.2	44.8	22.6	10.6	45.7
	M E	[kg/h] 1900.8	1891.3	1900.6	2191.2	1892.3	2192.5	1900.3	1240.9	2241.8
	T G, in	[°C]	80.9	62.2	101.4	80.6	62.0	101.1	77.0	61.8	90.8
	M G	[kg/h] 1900.0	1855.7	1899.9	1905.9	1895.8	1903.5	1900.0	1886.1	1904.6
	T AC, in	[°C]	29.4	17.6	37.2	29.4	17.6	37.5	28.3	17.6	37.2
	M AC	[kg/h] 5005.9	4986.9	5100.0	5006.1	4981.7	5100.0	5007.8	4982.3	5100.0

Table 3 : Selected ANNs with tangent hyperbolic AF and their corresponding BIC value 295

 3

	Normalization interval	TD Number of neurons	Initialization number	BIC
	[-0.3 0.3]	1	10	10	0.618
	[-0.3 0.3]	2	7	6	0.694
	[-0.3 0.3]	3	6	1	0.625
	[-0.8 0.8]	1	11	3	0.375
	[-0.8 0.8]	2	15	4	1.942
	[-0.8 0.8]	3	14	7	1.388
	[-1 1]	1	9	2	0.358
	[-1 1]	2	15	10	1.723
	[-1 1]	3	5	10	2.648

Table 4 : Selected ANNs with linear AF and their corresponding BIC value

 4

	Normalization interval	TD Number of neurons	Initialization number	BIC
	[-0.3 0.3]	1	4	1 0.965
	[-0.3 0.3]	2	6	15 1.135
	[-0.3 0.3]	3	4	7 0.699
	[-0.8 0.8]	1	5	6 0.790
	[-0.8 0.8]	2	5	9 0.411
	[-0.8 0.8]	3	7	9 0.735
	[-1 1]	1	6	14 0.752
	[-1 1]	2	6	5 0.867
	[-1 1]	3	5	5 0.668

Table 5 ,

 5 they do not include the zero values. The mean values are less than 3% which means that the performance of the ANN training process is good. However, to effectively test the generalization ability of the artificial neural model it is necessary to use a set of experimental unseen data.

Table 5 : Mean relative error values for the identification step

 5

		T E, out	T G, out	T AC, out
	Identification	2.9%	0.9%	0.9%

Table 6 : Mean relative error values 319

 6

		T E, out	P E	T G, out	P G	T AC, out	P AC
	Validation 1	4.9%	8.7%	1.1%	4.4%	1.2%	6.6%
	Validation 2	5.0%	9.3%	2.1%	9.3%	1.6%	7.9%

Table 7 : The absolute relative error of the transferred energy

 7

		Evaporator	Absorber-Condenser	Generator
	Validation 1	6.6%	2.3%	0.2%
	Validation 2	5.1%	2.6%	0.1%

Acknowledgments

This study has been supported and funded by the French Agency for Environment and Energy Management (ADEME), the 490 National Institute of Nuclear Sciences and Techniques (INSTN) and the French National Research Agency (ANR) in the framework of the "ABCLIMSOL" project sponsored by Prebat program 2007.