Mohamed Kara
email: mohamed.kara@univ-setif.dz

Salim Mesbahi
email: salimbra@gmail.com

Philippe Angot
email: philippe.angot@univ-amu.fr

The Fictitious Domain Method with Sharp Interface for Elasticity Systems with General Jump Embedded Boundary Conditions

Keywords: AMS subject classifications: 65M55, 65N08, 65N85 Fictitious domain method, sharp interface, elasticity system, jump embedded boundary conditions, finite volume method

In framework of the fictitious domain methods with immersed interfaces for the elasticity problem, the present contribution is to study and numerically validate the jump-integrated boundary conditions method with sharp interface for the vector elasticity system discretized by a proposed finite volume method. The main idea of the fictitious domain approach consists in embedding the original domain of study into a geometrically larger and simpler one called the fictitious domain. Here, we present a cell-centered finite volume method to discretize the fictitious domain problem. The proposed method is numerically validated for different test cases. This work can be considered as a first step before more challenging problems such as fluid-structure interactions or moving interface problems.

Introduction

The conventional grid method approach to simulating transmission problems in heterogeneous elastic body with complex boundaries is to discretize the governing equations on a curvilinear mesh that conforms to the boundaries. The main advantages of this approach are that imposition of boundary conditions is greatly simplified and the scheme can be easily designed so as to maintain sufficient accuracy and conservation property. However, the geometrical complexity of immersed boundaries can adversely impact the stability, the convergence and the operation count of the solver. An alternative method approach which is gaining popularity in recent years is the fictitious domain methods approach where the governing equations are discretized on a uniform structure Cartesian mesh which does not conform to the immersed boundaries. This greatly simplifies mesh generation and retains the relative simplicity of the governing equations in Cartesian coordinates, (see [1-8, 10, 11]). In addition, since the Cartesian mesh scheme does not depend on the location of the immersed boundary, there is no need for remesching strategies and, thus, has a significant advantage in simulating problems with moving boundaries and complicated shapes. However, the difficulty in using the fictitious domain method in conjunction with a uniform structured Cartesian mesh scheme is in the position of the original boundary conditions at the immersed boundaries. Thus, to avoid this complication, an approximate interface Σ h of the original Σ is constructed by a series of cell-sides that are cut by the boundary on which we apply an algebraic transmission boundary conditions proposed in our fictitious domain model-problem.

The numerical validation of a new fictitious domain method associated with general Jump Embedded Boundary Conditions (J.E.B.C) is proposed in the present paper. To solve the linear elasticity system governed by the given problem ( P) in the original domain Ω, a fictitious domain technique is used. The key issue here is to construct an extended imperfect transmission problem (P) of ( P) defined in the extended domain Ω of the original physical domain Ω in which its geometric shape is simpler than that of Ω such that

Ω = Ω∪Σ∪Ω e ,
where Ω e is the external domain and Σ the common interface between Ω and Ω e , see Fig. 1. With an appropriate choice of the data and the transmission conditions in the auxiliary domain Ω e and on the interface Σ respectively, the transmission problem (P) will be well-posed and the two problems are equivalents in the following sense : If u is a solution of the fictitious problem (P) defined in the fictitious domain Ω, then the restriction of the fictitious solution ũ = u| Ω or, at least, that u η = u| Ω is a solution of the physical problem (original) ( P).

Our objective is to use a simple structured mesh in Ω, e.g., a uniform Cartesian grid, independent of the shape of the immersed interface Σ, instead of using an adaptive mesh to which it becomes difficult to find fast and efficient solvers. The method is applied to several test cases for which an analytical and finite volume solution exists. Comparisons of the numerical and analytical results show a very good performance of the method. Our quest is to solve, with a fictitious domain method in Ω, the following problem originally defined in Ω ⊂ Ω with a general boundary condition on Σ : Find the displacement vector ũ and the symmetric stress tensor σ = (σ i,j ) d i,j=1 such that

     -∇•σ( ũ) = f in Ω, ũ = 0 on Γ, ũ = u D or -σ( ũ)•n = A ũ+q on Σ, (1.1)
where f ∈L 2 ( Ω) d is the body force, u D ∈ H 1/2 (Σ) d is the Dirichlet boundary condition, A ∈ L ∞ ( Σ) d×d a uniformly positive matrix and q ∈H -1/2 (Σ) d are given on Σ. Here λ and µ ∈ L ∞ ( Ω) are the Lamé coefficients, which we assume that inf

Ω φ(x) = φ -> 0, φ + = sup Ω φ(x),
where φ stands for λ, µ or ρ, and the constitutive equation in an isotropic linear elastic solid may be written as

σ( ũ) = λ∇• ũI+2µε( ũ), here ε( ũ) = 1 2 [∇ ũ+(∇ ũ) t ],
is the linear strain tensor, I is the identity tensor, the Lamé coefficients λ and µ are related to the its Young's modulus E and the Poisson's ratio ν by the following relations:

λ = Eν (1-2ν)(1+ν) , µ = E 2(1+ν) .
In the next section the governing equations of the fictitious domain problem together with the treatment of the embedded original boundary conditions on the common interface Σ are recalled. This is followed by a description of the finite volume discretization scheme and the solution algorithm. Finally, the method's capabilities are demonstrated by applying it to a number of test cases. Moreover, a particular scalar test problem with singular solution is also simulated.

Fictitious model with jump embedded transmission conditions on Σ

In this section, we present the governing equations of the fictitious model, together with the weak form and treatment of the embedded original boundary conditions. Ω is defined by ∂Ω = Γ∪Γ e . Let n be, either the outward unit normal vector on Γ e , or the outward unit normal vector on Σ oriented from Ω to Ω e . For a function ψ ∈ H 1 ( Ω∪Ω e ), let ψ -and ψ + be the traces of ψ| Ω and ψ| Ω e on each side of Σ respectively, ψ| Σ = 1 2 (ψ -+ ψ + ) the arithmetic mean of trace of ψ, and ψ Σ = (ψ + -ψ -) the jump of traces ψ on Σ oriented by n.

The strong form

For the data f ∈ L 2 (Ω) d , g and h given in H -1/2 (Σ) d , we consider the elasticity system problem for the real-valued vector u = (u 1 ,••• ,u d ) defined in Ω and includes immersed transmission conditions on Σ links the trace of jumps of both the normal stress vector σ(u)•n = (λ(∇•u)I+2µ ε(u))•n and the traces of solution u through the interface Σ :

           -∇•σ(u) = f in Ω = Ω∪Ω e , u = 0 on Γ∪Γ e , σ(u)•n Σ = Mu| Σ -h on Σ, σ(u)•n | Σ = S u Σ -g on Σ, (2.1)
where the Lamé coefficients µ, λ and the transfer matrices S and M in the (J.E.B.C.) third and fourth equations of (2.1) on Σ are measurable and bounded functions verifying ellipticity assumptions

(A1) µ, λ ∈ L ∞ (Ω), ∃ µ -, λ -> 0, µ > µ -, λ > λ - a.e. in Ω, (A2) M, S ∈ L ∞ (Σ) d×d , ∀ξ ∈ R d , ξ t •M(x)•ξ > 0, ξ t •S(x)•ξ > 0 a.e. on Σ.
Moreover, in order to recover the original problem (1.1) in Ω, we choose µ and λ a L ∞extension of μ and λ respectively in Ω and f a L 2 -extension of f in Ω such that

µ = μ in Ω, µ e in Ω e , λ = λ in Ω, λ e in Ω e , f = f in Ω, f e in Ω e ,
as well as, the boundary condition the second equation in (2.1) on Γ e is chosen to ensure the solvability of the fictitious domain problem (2.1).

The weak form

We now introduce the Hilbert space :

W(Ω) = v ∈ L 2 (Ω) d , v| Ω ∈ H 1 ( Ω)
∈ W such that a(u,v) = l(v), ∀v ∈ W, (2.2)
with

a(u,v) = Ω σ(u):∇vdx+ Σ Mu| Σ v| Σ ds+ Σ S u Σ v Σ ds, l(v) = Ω fvdx+ h,v| Σ -1 2 , 1 2 + g, v Σ -1 2 ,
u W ≤ 1 α 0 ( f L 2 (Ω) +c(Ω,Ω e )( g -1 2 ,Σ + h -1 2 ,Σ )).
Proof. We begin by deriving the weak form of the problem (2.1). With the first and second equations in (2.1) and using the Green-Stokes formula, •,• -1 2 ,Σ being the duality pairing between H -1 2 (Σ) d and H 1 2 (Σ) d , we get respectively over Ω and Ω e :

Ω σ(u) : ∇vdx-σ(u) -•n,v - -1 2 ,Σ = Ω fvdx, ∀v ∈ W, v| Ω e = 0, Ω e σ(u) : ∇vdx+ σ(u) + .n,v + -1 2 ,Σ = Ω e fvdx, ∀v ∈ W, v| Ω = 0.
Summing now the two previous equations yields :

Ω σ(u) : ∇vdx+ σ(u) + •n,v + -1 2 ,Σ -σ(u) -•n,v - -1 2 ,Σ = Ω fvdx, ∀v ∈ W.
Then, noticing that for any bilinear form •,• -1 2 ,Σ defined on Σ, we have the equality below :

u + ,v + Σ -u -,v -Σ = u Σ ,v Σ Σ + u Σ , v Σ Σ , ∀u,v,
we obtain the following weak form in Ω :

Ω σ(u) : ∇vdx+ σ(u)•n Σ ,v Σ -1 2 ,Σ + σ(u)•n Σ , v Σ -1 2 ,Σ = Ω fvdx, ∀v ∈ W.
Then, using the jump transmission conditions i.e., the third and fourth equations in (2.1) on Σ, we get the nice weak formulation below: Find u ∈ W, sush that ∀v ∈ W,

Ω σ(u) : ∇vdx+ Σ Mu Σ v Σ + Σ S u Σ v Σ = Ω fvdx+ h,v Σ -1 2 ,Σ + g, v Σ -1 2 ,Σ .
With the ellipticity assumptions (A1), (A2), it is now easy to verify using the Korn inequality, in Ω, Ω e and standard trace lemmas that the left-hand side of the above nice weak formulation is a bilinear continuous and coercive form in W×W, whereas the righthand is a linear continuous form in W. Hence, by the Lax-Milgram theorem, we have existence and uniqueness of the weak solution u in W.

Remark 2.1. Since ∇•σ(u) ∈ L 2 (Ω) d , then we can define σ(u) ± Σ in H -1/2 (Σ) d .

Treatment of the original boundary condition on Σ

In this section, we will use for each kind of desired original boundary condition imposed on the immersed boundary Σ, only the particular variant with no exterior or surface control for both Fourier or Dirichlet boundary conditions respectively. These conditions will be enforced using a thin approximation of the immersed interface Σ, either the so-called Exterior interface Σ Ext. or Cut interface Σ Cut. by surface penalty of the Dirichlet boundary conditions and by choosing the transfer matrices in the first equation of system (2.3) which satisfies Eq. (2.5) or the third equation in (1.1) in the original problem, see [START_REF] Girault | Error analysis of a fictitious domain method applied to a Dirichlet problem[END_REF][START_REF] Khadra | Fictitious domain approach for numerical modilling of Navier-Stokes equations[END_REF][START_REF] Lions | Problèmes aux Limites dans les équations aux Dérivées Partielles[END_REF][START_REF] Mesbahi | Mathematical analysis of a reaction diffusion model for image restoration[END_REF][START_REF] Ère | A fictitious domain approach with spread interface for elliptic problems with general boundary condition[END_REF][START_REF] Ère | A general fictitious domain method with immersed jumps and multilevel structured meshes[END_REF][START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF][START_REF] Sarthou | The sub-mesh penalty method, Finite Volumes for Complex Applications V[END_REF][START_REF] Ye | An accurate Cartesian grid method for viscous incompressible flow with complex immersed boundaries[END_REF].

Embedded Fourier or Neumann boundary condition on Σ

Let -σ(u) - Σ •n and -σ(u) + Σ •n be the traces of the normal stress on each side of Σ. In the previous fictitious domain model, we use the algebraic transmission conditions the third and fourth equations in (2.1) on Σ. the two quantities of the traces can be written in the following way :

     -σ(u) - Σ •n = S+ 1 4 M u - Σ -S- 1 4 M u + Σ +g- h 2 on Σ, -σ(u) + Σ •n = -S+ 1 4 M u + Σ + S- 1 4 M u - Σ +g+ h 2 on Σ.
(2.3) λ|

Ω e = µ| Ω e = 1, f| Ω e = 0 M = 4S = 2A, h 2 -g = q Dirichlet (D1)
λ|

Ω e = µ| Ω e = 1, f| Ω e = 0 M = 4S = 2 η , h 2 -g = 1 η u D Dirichlet (D2)
λ|

Ω e = µ| Ω e = 1 η , f| Ω e = 1 η u e S = 1 η , M = g = h = 0
Then, with the particular choice S = 1 4 M the system of Eq. (2.3) yields the Fourier boundary condition below, in this case the exterior fictitious domain problem and the interior (physical) one are independent : -

     -σ(u) - Σ •n = 1 2 Mu - Σ +g- h 2 on Σ, -σ(u) + Σ •n = - 1 2 Mu + Σ +g+ h 2 on Σ.
σ(u) - Σ •n = Au - Σ +q on Σ. (2.5)
This formulation enables us to deduce the transmission coefficients S, M, g and h on Σ and the data µ, λ and f in Ω e as well as on Σ, which are presented in the Table 1.

Other variants are proposed there require either u + Σ or -σ(u) + Σ •n to be controlled by L 2 or H 1 volume penalty methods [START_REF] Ph | A unified fictitious domain model for general embedded boundary conditions[END_REF][START_REF] Ph | A model of fracture for elliptic problems with flux and solution jumps[END_REF][START_REF] Ph | On the well-posed coupling between free fluid and porous viscous flows[END_REF][START_REF] Ph | A fictitious domain model for the Stokes/Brinkman problem with jump embedded boundary conditions[END_REF][START_REF] Sarthou | The sub-mesh penalty method, Finite Volumes for Complex Applications V[END_REF] performed with the parameters in Ω e for the convection-diffusion problem.

Embedded Dirichlet boundary condition on Σ

The Dirichlet boundary condition u - Σ =u D of the problem (1.1) is treated as a penalization of the previous Fourier boundary condition in the system of Eqs. (2.3) or (2.4). Indeed, the Fourier condition in the system of Eqs. (2.3) or (2.4) can be penalized by a surface penalty on Σ with

M = 2 η I -→ +∞ and h 2 -g = 1 η u D ,
when η -→ 0, where η is a real penalty parameter, (see (D1)) in the Table 1. The second variant consists in using a L 2 or H 1 volume penalty in the exterior domain Ω e such that lim

η-→0 u + ηΣ = u D for the model with S = 1 η I -→ +∞, M = 0, h = 0,
and thus u Σ -→ 0, φ•n Σ = 0, as proposed in [START_REF] Delcourte | Développement de Méthodes de Volumes Finis Pour la Mécanique des Fluides[END_REF].

Remark 2.2. In the practice, the penalty parameter η is chosen small enough, typically of the order of 10 -7 such that the error due to the fictitious domain modeling in Ω for the original problem Ω remains smaller than the error due to the numerical approximation.

Characterization of the fictitious domain (original and auxiliary)

For the numerical approximation of the fictitious domain problem (2.1), we use the finite volume methods because they are robust, relatively cheap, inherently conservative, and able to handle strong discontinuities in the differential operators. We consider the two-dimensional case (d = 2) in the Cartesian coordinates (x,y), when the domain Ω is rectangular for a fictitious domain model. The discretization procedure is separated in two parts: discretization of the computational domain and equations discretization, see [5-8, 10-12, 16, 18, 19].

In this part, we use the level set method which is a numerical technique for tracking moving interface, it is based upon the idea of representing the interface as a level set curve of higher-dimensional function φ(x,t). In keeping with the applications pursued in this paper, we only use the level set theory for representing the static interface Σ.

We decompose the computational domain (fictitious domain) Ω into two subdomains Ω and Ω e which are respectively associated to the part of the physical domain and auxiliary one. The principle is to define a distance function in the computational domain. If the interface of the physical domain is given by Σ, we can define by ϕ(x,0) =

x-x Σ in Ω, -x-x Σ in Ω e , the characteristic function of the physical domain is easily obtained by looking at the sign of ϕ(x i ,0) in the nodes x i of the mesh :

χ ph (x i ) = 1, if ϕ(x i ,0) > 0, 0, if ϕ(x i ,0) < 0,
finally, the characteristic function from the distance function in the control volume K can be given by

χ ph | K = Σ + K |Σ| K ,
where Σ + K is the sum of the nodes of element K which are ϕ(x i ,0) >0, and |Σ| K is the sum of absolute values of ϕ(x i ,0). the value of function χ s | K is 1 or 0 depending on whether all nodes of the element K are positive or negative, but for the control volumes that are crossed by the embedding interface Σ, they, have a value between 0 and 1.

Numerical results

In order to illustrate the accuracy of the fictitious domain method for the elasticity problems, we have considered some problems with exact solutions. In that follows, we focus on 2-D problems. We present the results obtained from numerical simulations carried out on homogeneous and heterogeneous isotropic media. For the sake of simplicity, we carried out the calculations for the model where the Lamé Coefficients are constant, see [1-4, 9, 16] and [START_REF] Ère | A general fictitious domain method with immersed jumps and multilevel structured meshes[END_REF].

In the computation, the material properties are chosen as : Young's modulus E=1.542 and poisson's ratio ν = 0.285. These correspond to the Lamé's coefficients λ = 8000 and µ = 6000. The errors between the numerical and analytic solutions to the following test problems can be appreciated by the calculation of either the L 2 -norm or the L ∞ -norm in Ω h .

Validation of the Finite Volume Scheme (F.V.S.) for the elasticity problem

We begin by the validation of the finite volume scheme in the areas where the computational domain coincides with the physical one. For this purpose, the physical domain is the unit square Ω = [0,1]×[0,1]. The domain Ω is meshed by uniform square cells K with a grid step varying from h = 1 10 to 1 120 and Σ = {(x,y) ∈ Ω such that x = 1}. The first test used to demonstrate the convergence of the proposed scheme, the nonhomogeneous Dirichlet and Fourier problems are considered and the Fig. 2 illustrates the relative L 2 -norm errors versus a discretisation step h which shows the second order accuracy of this scheme.

Validation of the fictitious domain method for two-dimensional linear elasticity system

The conventional way of validating the efficient accuracy of the present scheme is to compute an elastic body which has both a L-shaped or a curved immersed boundary Σ for which an analytical solution exists. In the numerical computations, the fictitious domain is the unit square, a convergence study is carried out using uniform rectangular meshes, the tests are performed for two different domains : L-shaped domain and a quarter disk domain. In the figures of errors, the plot of error in the energy norm versus the mesh spacing is shown on a log-log plot.

First test problem: A L-shaped domain

Dirichlet problem. First, we consider the following homogenous Dirichelet problem :

P DH    -∇•σ( ũ) = f in Ω, ũ = 0 on Γ, ũ = 0 on Σ(u D = 0), 10 -1 10 -1.3 10 -1.6 10 -2 10 -4 10 -3 10 -2 log(discretization step h) log( u-u h / u )

Fourier (FP) and homogeneous Dirichlet problems (DP)

Log(er F ). Log(h 2 ) Log(er D ). where Σ = Σ x ∪Σ y , such that

Σ x = (x,y) x = 1 2 ; 1 2 ≤ y ≤ 1 and Σ y = (x,y) y = 1 2 ; 1 2 ≤ x ≤ 1 ,
which has the analytic solution ũ(x,y) = ṽ(x,y) = sin(2πx)sin(2πy) in Ω, for the appropriate data for f . Fig. 3 gives an example of the uniform grid and the immersed interface. This geometry is used for the test problems presented in this section.

The fictitious domain problem is solved in Ω with the J.E.B.C. method (D1), using the surface penalty on Σ h without exterior control in Ω e . In this case, the equation solved in Ω e,h has no influence on the solution obtained in the physical domain. So, the equation parameters in Ω e,h are arbitrary extensions of the original coefficients of physical problem. The Table 2 synthesizes the results obtained with the J.E.B.C. method (D1) and the parameter of penalization η = o(10 -7 ). Note that with a Cartesian mesh on Ω, the approximate interface Σ h of Σ is exactly defined. So, the global error is defined as the sum of a modelling error (due to the penalization in the Dirichlet case) and numerical scheme error.

Fig. 4 shows the convergence of the numerical error with respect to the discretization step h by J.E.B.C. method (D1) in the L-shaped domain Ω.

In order to validate the J.E.B.C. method (D1) for any Dirichlet boundary condition,

0 1 1 Γ σ(u)•n = 0 Γ e σ(u)•n = 0 o ù u = 0 Σ u = u D Ω , P , ũ Ωe
( P DNH )      -∇•σ( ũ) = f in Ω, ũ = 0 on Γ, ũ = u D on Σ,
where

u D =      3 16 sin(πy), 3 16 sin(πy), 3/4 ≤ y ≤ 1, x = 3/4, u D = x(x-1)sin(3π/4), x(x-1)sin(3π/4), 3/4 ≤ x ≤ 1, y = 3/4, 10 -1 10 -1.3 10 -1.6 10 -2 10 -2 10 -1 log(discretization step h) log( u-u h / u ) Penalization parameter η = 10 -7
Log(er). Log(h 1 ). which has the analytic solution ũ(x,y) = ṽ(x,y) = x(x-1)sin(πy) in Ω, the table below gives the relative error norm for non homogeneous Dirichlet problem in the L-Shaped Domain (see Table 3).

The fictitious problem over Ω is solved using J.E.B.C. method (D1). As Fig. 5 shows, the first order a accuracy is reached for L 2 -norm for exterior approximate interface.

Robin problem. We now consider the Robin problem with an appropriate data of f such that :

P N      -∇•σ( ũ) = f in Ω, ũ = 0 on Γ, σ( ũ)•n = ũ+g on Σ (A = αI = I, q = g), 10 -1 10 -1.3 10 -1.6 10 -2 10 -3 10 -2 10 -1 log(discretization step h) log( u-u h / u ) Penalization parameter η = 10 -7
Log(er1). Log(h 1.2 )

on the Σ x = {x = 1, 3/4 ≤ y ≤ 1} and g = (-µ(2x-1)+αx(x- 1
))sin(3π/4)-µπx(x-1)cos(3π/4), (-λ(2x-1)+αx(x-1))sin(3π/4)-π(λ+2µ)x(x-1)cos(3π/4),

on the Σ y = {3/4 ≤ x ≤ 1, y = 1} , Σ x ∪Σ y = Σ.
The results obtained are illustrated in the Table 4. Fig. 6 shows a second-order accuracy of the J.E.B.C. method (F), for the Robin condition.

10 -1 10 -1.3 10 -1.6 10 -1.9 10 -2. 15 10 -4 10 -3

10 -2 log(discretization step h ) log( u ex -u h / u ex ) Fourier problem in L-Shaped domain, Σ = Σ h Log(error).
Log(h 2 ) The approximate solution obtained with the Fourier approach and the analytic one at the diagonal square section is represented in Fig. 7. It can be seen that the behavior of the approximate solution at the corner (x = 3 4 , y = 3 4 ) of the physical domain Ω.

Second test problem. A quarter disk domain

We now consider a quarter disk domain Ω immersed in the unit square Ω = [0,1]×[0,1].

That defines an immersed interface Σ for quarter disk, see Fig. 8. As in the L-shaped case, the fictitious square domain is meshed with uniform grid step h. The approximate interface Σ h , lying on sides of the mesh, is chosen such that it crosses the physical immersed interface Σ. Fig. 8 describes of the quarter disk domain, associated fictitious domain mesh and approximate immersed interface.

Homogeneous Dirichlet case. We consider the below homogeneous Dirichlet Problem with the adapted right-hand side f which has the below analytic solution ũ(x,y) = ṽ(x,y) = 4xy 1-exp 1 4 (r 2 -1)

P DH      -∇•σ( ũ) = f in Ω, ũ = 0 on Γ, ũ = u D = 0 on Σ,
in Ω, which r = x 2 +y 2 , 6.504512e-2 2.863497e-2 9.679815e-3 3.856029e-3 3.238966e-3 where (r,θ) are the polar coordinates at the origin. Table 5 gives the relative error norms for homogeneous Dirichlet problem in square domain, with η = h 2(λ+2µ) , using J.E.B.C. method (D1) with the cut approximate interface Σ cut h . Now, we use the following variant g = q = 0 and M = S = 1 η , (see D2 in Table 1). The following computations are performed with η = 1e-10 to get the modeling error negligible compared to the discretization error using J.E.B.C. method (D2) with cut approximate interface Σ cut (see Table 6).

P DNH        -∇•σ( ũ) = f in Ω, ũ = 0 on Γ, ũ = u D = 3 2 sin(2θ) on Σ,
where θ = arctan( y x ), has the analytic solution ũ(x,y) = ṽ(x,y) = xy(4-r 2 ) in Ω, which r = x 2 +y 2 .

The computations are illustrated in Table 7. In Fig. 10 we show a log-log plot of the errors both displacement components u and v versus the discretization step h. Also shown is a straight with slope of 1 which corresponds to first order accuracy for the method (D1).

The approximate solution obtained with the J.E.B.C. method (D1) and the analytic one at the diagonal square section is represented in Fig. 11.

Robin or non homogeneous Neumann case.

When the interface Σ is approximated by the mesh volume edges into Σ h , the non-homogeneous Neumann (or Robin) embedded boundary condition case requires the surface correction in the scalar case, as shows the test below. But in the vectorial one, the surface correction is not required.

P N      -∇•σ( ũ) = f in Ω, ũ = 0 on Γ, σ( ũ)•n = q on Σ.
We have chosen f such that the problem ( P N ) has the analytical solution

ũ(x,y) = ṽ(x,y) = 4xy 1-exp 1 4 (x 2 +x 2 -1) in Ω, where q = -2((λ+2µ)y(1-y 2 )+λy 2 1-y 2 )n x -2µ(x(1-x 2 )+x 2 √ 1-x 2 )n y , -2µ(y(1-y 2 )+y 2 1-y 2 )n x -2((λ+2µ)x(1-x 2 )+λx 2 √ 1-x 2 )n y .
10 -1 10 -1. The fictitious domain is solved in Ω with the J.E.B.C. method (F) without exterior control, as described in the

P F      -∇•σ( ũ) = f in Ω, ũ = 0 on Γ, σ( ũ)•n = u+q on Σ, which has the analytic solution ũ(x,y) = ṽ(x,y) = xy(4-(x 2 +y 2 )) in Ω.
The fictitious domain problem is solved in Ω with the J.E.B.C. method (F) without exterior control, as described in Table 1. We investigate in Fig. 13 the effect of the cut approximate interface Σ cut h and exterior one Σ ext h . We observe the first-order accuracy for the L 2 -norm for both approximated immersed interfaces Σ cut h and Σ ext h , since the method involving the cut interface is more accurate than the exterior one (see Fig. 13). M. Kara, S. Mesbahi and P. Angot / Adv. Appl. Math. Mech., 12 (2020), pp. 1-22 10 -1 10 -1.3 10 -1.6 10 -1.9 10 -2.15

Conclusions and perspectives

An approach of the fictitious domain method to solve linear elasticity problems has been introduced which is based on a sharp interface approach and treats boundary conditions General immersed (Dirichlet, Neumann and Robin). A numerical finite volume scheme with jumps of flux and solution through the immersed interface is used to calculate the solution. The main advantage of this method is its low cost. Indeed, the numerical solution is calculated on a simple Cartesian grid without locally modifying the numerical scheme or introducing an unknown local. This method is a first-order accuracy.

In the next step of this work, we shall introduce another fictitious domain method based on a diffuse interface approach and we will use a numerical finite element scheme to deal with general immersed boundary conditions. These two fictitious domain methods remain of the first order. Then, each method will be combined with a local multi-level algorithm refinement mesh to increase the accuracy of the solution.

The other perspective will be to extend these methods in the fictitious domain elasticity problems with moving boundary conditions. These fictitious domain methods allow to simulate the moving boundaries with a low computational cost and require no remeshing of the domain.

Conclusions and perspectives

An approach of the fictitious domain method to solve linear elasticity problems has been introduced which is based on a sharp interface approach and treats boundary conditions General immersed (Dirichlet, Neumann and Robin). A numerical finite volume scheme with jumps of flux and solution through the immersed interface is used to calculate the solution. The main advantage of this method is its low cost. Indeed, the numerical solution is calculated on a simple Cartesian grid without locally modifying the numerical scheme or introducing an unknown local. This method is a first-order accuracy.

In the next step of this work, we shall introduce another fictitious domain method based on a diffuse interface approach and we will use a numerical finite element scheme to deal with general immersed boundary conditions. These two fictitious domain methods remain of the first order. Then, each method will be combined with a local multi-level algorithm refinement mesh to increase the accuracy of the solution.

The other perspective will be to extend these methods in the fictitious domain elasticity problems with moving boundary conditions. These fictitious domain methods allow to simulate the moving boundaries with a low computational cost and require no remeshing of the domain.

Figure 1 :

 1 Figure 1: Two examples of original domain (physical) Ω embedded in a fictitious domain Ω = Ω∪Ω e .

(2. 4 )

 4 When the Fourier or Neumann boundary condition in the third equation in (1.1) is desired for the original problem (1.1) in Ω, the following immersed boundary condition on Σ must be satisfied by the solution of the fictitious domain problem in Ω :

Figure 2 :

 2 Figure 2: The second order accuracy of the schema for a (FP) and (DP) problems.

, Pe , ue x y Figure 3 :

 y3 Figure 3: L-shaped physical domain immersed inside a fictitious domain Ω = Ω∪Σ∪Ω e = [0,1]×[0,1] in an uniform grid.

Figure 4 :

 4 Figure 4: First order accuracy of the J.E.B.C. method (D1) for P DH problem.

Figure 5 :

 5 Figure 5: Order of the convergence of the J.E.B.C. method (D1) for P DNH problem.

Figure 6 :

 6 Figure 6: Order of the convergence of the J.E.B.C. method (F) for P N problem.

Figure 7 :

 7 Figure 7: Approximate solution and analytic one in the Fourier case at a diagonal square section.

Figure 8 :

 8 Figure 8: Approximate interface Σ h and approximated domain Ω h . (a) Cut approximate interface Σ Cut h , (b) Exterior approximate interface Σ Ext h .

Fig. 9 .

 9 Convergence for the L 2 -norm of the discretization error for the homogeneous 10

Figure 9 :

 9 Figure 9: The first order accuracy for (D1) and (D2) methods with cut interface in a quarter disk domain.

Figure 10 :

 10 Figure 10: The first order accuracy for (D1) J.E.B.C. method for non homogeneous Dirichlet problem.

Figure 11 :

 11 Figure 11: Approximate solution and analytical one in the non homogeneous Dirichlet case at a diagonal square section.

Figure 12 :

 12 Figure 12: The first order accuracy of the J.E.B.C. Method (F).

Figure 13 :

 13 Figure 13: Convergence for the L 2 -norm of the discretization error with h of the J.E.B.C. method (F) for the Fourier problem with Σ cut h or Σ ext h .

  d and v| Ω e ∈ H 1 (Ω e ) d ; v = 0 on Γ∪Γ e equipped with the natural inner product and associated norm inH 1 (Ω) d .Let u ∈ W satisfy the first equation in (2.1) and f ∈ L 2 (Ω) d , the weak form of the problem (2.1) can be written as : Find u

unique weak solution u ∈ W(Ω) satisfying (2.2) for all v ∈ W(Ω), such that ∃α 0

	1 2	.
	Theorem 2.1 (Global solvability of the fictitious domain model with J.E.B.C.). If the ellipticity assumptions (A1), (A2) hold, the problem (2.1) with f ∈ L 2 (Ω) d , and g, h ∈ L 2 (Σ) d has a

(Ω,Ω e ,µ, S , M ) > 0,

Table 1 :

 1 Summary of the data in Ω e and on Σ for the J.E.B.C. Method.

	J.E.B.C. method	Parameters in Ω e	Parameters on Σ
	Fourier (F)

Table 2 :

 2 Variation of relative L 2 error norms according to the number of discretization points for Dirichlet condition in L-shaped domain.

	Mesh step h u ex -u h u ex	1/10 7.267270e-2 3.516438e-2 1.696973e-2 8.292075e-3 4.728605e-3 1/20 1/40 1/80 1/140
	we now study a non-homogeneous Dirichlet problem

Table 3 :

 3 The relative L 2 error norms according to the number of discretization points for non homogeneous Dirichlet problem in L-shaped domain.

	Mesh step h u ex -u h u ex	1/16 1.075866e-2 4.604882e-3 2.164175e-3 1.032435e-3 9.361342e-4 1/32 1/64 1/128 1/140

Table 4 :

 4 Variation of relative L 2 error norms according to the number of discretization points for Fourier condition in L-shaped domain.

	Mesh step h u ex -u h u ex	1/16 4.914116e-3 1.181075e-3 2.825945e-4 6.513587e-5 5.352028e-5 1/32 1/64 1/128 1/140
	where
		g =	    	--	1 2 1 2	(λ+2µ)+ µ+ 3α 16 sin(πy)+ 3α 16 sin(πy)+ 3 16 πµcos(πy), 3 λπcos(πy), 16

Table 5 :

 5 Variation of relative L 2 error norms according to parameter h for homogeneous Dirichlet problem in the quarter disk domain using (D1).

	Mesh step h |u ex -u h | |u ex |	1 10	1 40	1 80	1 120	1 140

Table 6 :

 6 Variation of relative L 2 error norms according to parameter h for homogeneous Dirichlet problem in the quarter disk domain using (D2).

	Mesh step h |u ex -u h | |u ex |	1 10 6.494105e-2 2.470404e-2 1.024406e-2 8.240813e-3 5.464180e-3 1 40 1 60 1 80 1 120

Table 7 :

 7 Variation of relative L 2 error norms according to parameter h for non homogeneous Dirichlet problem in the quarter disk domain using (D1).

	Mesh step h |u ex -u h | |u ex |	1 20 3.026736e-1 1.168293e-1 9.872092e-2 6.252881e-2 3.941031e-2 1 40 1 60 1 80 1 120
	Dirichlet elasticity problem with Σ cut h and using (D1) and D2 J.E.B.C. methods.
	Non homogeneous Dirichlet case. We consider the below non homogeneous Dirichlet Problem defined in Ω

Table

  1, the results for the relative L 2 -norm error are registered in the

	The J.E.B.C. method (F) of Neumann problem with Cut interface Σ cut h
		10 0
	ex )
	u
	10 -1 -u h /
	ex
	u
	log(	10 -2	J.E.B.C. method, Log(h 1 )
		10 -2	10 -1
			log(discretization step h )

Table 8 :

 8 Variation of relative L 2 error norms according to parameter h for a Fourier problem in the quarter disk domain using (F).

	Mesh step h u ex -u h u ex	1 10 9.777809e-1 3.026692e-1 1.168142e-1 6.250331e-2 3.349663e-2 1 20 1 40 1 80 1 140

Table 8 .

 8 Finally, In order to validate the J.E.B.C. method (F) for any Robin boundary conditions, we now study a Fourier problem

  Figure 13: Convergence for the L 2 -norm of the discretization error with h of the J.E.B.C. method (F) for the Fourier problem with Σ cut h or Σ ext h .

			The J.E.B.C. method (F) of Fourier problem
	ex )	10 -1
	u
	-u h /
	u ex log(	10 -2	Cut inter. Exterior inter.
			Log(h 1.2 )
			log(discretization step h )


André Monin
email: monin@laas..fr

Gérard Salut

Multi-emission discrimination using IIR Volterra Filtering

Introduction.

The aim of this paper is to build a test algorithm devoted to detect when several independent emissions are using the same carrier and, in such a situation, to determine the number of them.

Let us consider the case of phasis or frequency modulation. If there is only one modulation in the observed signal, its module remains constant with respect to time and the fluctuations of the module of the observed signal is only due to additive white noise. If now, one consider the sum of two modulations, as the phasis of each modulations are assumed to be independents, the module of the observed signal will be fluctuating. As a consequence of this remark, the algorithm which is proposed is based on the prediction of the module of the observed signal conditionally to the number of modulation supposed. The innovations involved are then used to construct a sufficient statistic for the detection in the sense that its expectation maximize some contrast criterium.

Statistics of a multi-emission.

In our problem, the hypothesis to be tested is the number of emissions in the observed complex signal obtain after two channels demodulation. It can be written:

(2.1) where:

(2.2)
is the sum of independent modulations and is a complex white noise with real and imaginary assumed to be independent and of variance . The moments of the white noise are all zero except the following:

(2.3) with:

n y t y t s t n v t + = s t n x t i i 1 = n ∑ = n v t v t 1 jv t 2 + = R v t E v t 2k [ ] R k C k p α p α k p - p 0 = k ∑ = and
(2.4)

Let us examine the statistics of multi-emission in the particular case of a digital phasis modulation. Each signal can be represented, after demodulation, as the complex signal:

(2.5) with where is the (constant) amplitude of the signal and the initial phasis.

is the current phasis carrying the information to be send. Let be the sample period. Between to samples, the evolution of the signal is defined as:

(2.6) where designates the phasis increment. As the transmitted information is unknown, may be viewed as a white noise with characteristics depending of the kind of modulation used. Let . For an analogical modulation, the density of probability of may be assumed gaussian. In such a situation, if the variance of is , the moments of are defined by:

(2.7)

where is complex conjugate of . For a digital modulation, the phasis jumps time of occurrence may be modelised with a Poisson process of frequency and the amplitude uniformly distributed on points on the unit circle. In that case, the density of probability of is:

(2.8)

where is the Dirac distribution centered on . The resulting moments of are then defined by:

(2.9)

One may of course consider other kind of modulations but it is important to note that, for all types of modulations, one has always the same property:

(2.10)

with . Let us now analyze the evolution of moments of one modulation. According to (2.10), one has:

(2.11)

As a consequence, because

, it is clear that all moments of are asymptotically zero expect the set which are constant.

C k p k! p! k p - ( )! --------------------------- = α p 1 3 … 2 p 1 - ( ) × × × = x t x 0 jϕ t ( ) exp = x 0 A 0 jϕ 0 ( ) exp = A 0 ϕ 0 ϕ t t ∆ x t t ∆ + x t j ϕ t ∆ ( ) exp = ϕ t ∆ ϕ t ∆ w t j ϕ t ∆ ( ) exp = ϕ t ∆ ϕ t ∆ Q w t k l E w t k w t * ( ) l [ ] , ≠ ∀ 1 2 ---k l - ( ) 2 Q -     exp = k ∀ E w t k w t * ( ) k [ ] , 1 = w * w λ m ϕ t ∆ p ϕ t ∆ ( ) 1 λ t ∆ - ( )δ 0 ϕ t ∆ ( ) λ t ∆ m 1 - -------------- δ 2πp m ---------- ϕ t ∆ ( ) p 1 = m 1 - ∑ + = δ a ϕ ∆ ( ) a w t k l E w t k w t * ( ) l [ ] , ≠ ∀ 1 m m 1 - --------------λ t ∆ - = k ∀ E w t k w t * ( ) k [ ] , 1 = k l E w t k w t * ( ) l [ ] , ≠ ∀ a = k ∀ E w t k w t * ( ) k [ ] , 1 = a 1 < k l E x t t ∆ + k x t t ∆ + * ( ) l [ ] , ≠ ∀ aE x t k x t * ( ) l [ ] = k ∀ E x t t ∆ + k x t t ∆ + * ( ) k [ ] , E x t k x t * ( ) k [ ] = a 1 < x t E x t k x t * ( ) k [ ] k 1 ≥ , { }
The first moment of a multi-emission can now be computed using (2.11). One may derive that all moments of degree are asymptotically zero except the following:

(2.12)

(2.13)

Let us compare these moments when the global power is the same but when the number of modulation is different. On has:

(2.14) and the fourth order moment is:

(2.15)
This relation means that the number of signals is detectable as the ratio within the fourth order moment and the square of the second order moment depends on . This property will be exploited in the sequel to construct the detection test based on predictors of the signal square module.

The IIR Volterra quadratic predictor.

Let us applied the quadratic Volterra predictor as defined in [2] to our problem. When one deals with real numbers, the purely quadratic IIR predictor of any function of the state has the following recursive structure:

(3.1)
where is the output and are time variant parameters computed to minimize the mean square prediction error, that is:

(3.2)
The optimal parameters are computed using the following orthogonalities:

(3.3)
which leads to the following linear system: --------------------------------E s

(3.4) 4 ≤ E s t n s t n * [ ] E x t k x t k * [ ] k 1 = n ∑ = E s t n s t n * ( ) 2 [ ] E x t k x t k * ( ) 2 [ ] k 1 = n ∑ 2 E x t k x t k * [ ]E x t l x t l * [ ] k l ≠ k l , 1 = n ∑ + = P k ∀ E x t k x t k * [ ] , P n ⁄ = E s t n s t n * ( ) 2 [ ] 1 2 n 1 - ( ) + n -
t n s t n * [ ] ( ) 2 = n φ t p t t ∆ + A t 1 p t A t 0 y t + = z t t ∆ + A t 2 z t p t t ∆ + y t + = z t A t i { } J t 1 2 ---E φ t z t - 2 [ ] = E φ t t ∆ + z t t ∆ + - ( ) y t 2 p t T y t z t T 0 = A t 2 E z t y t ( ) 2 [ ] A t 1 E p t y t ( ) 3 [ ] A t 0 E y t ( ) 4 [ ] + + E φ t t ∆ + y t ( ) 2 [ ] = A t 2 E z t p t T y t [ ] A t 1 E z t 1 p t T y t ( ) 2 [ ] A t 0 E p t T y t ( ) 3 [ ] + + E φ t t ∆ + p t T y t [ ] = A t 2 E z t z t T [ ] A t 1 E p t z t T y t ( ) 2 [ ] A t 0 E z t T y t ( ) 2 [ ] + + E φ t t ∆ + z t T [ ] =
The optimal parameters are then computable if all the expectations in (3.4) are computable. Note that this condition is achieved when the prediction system of is bilinear. Let us generalize this estimator to the complex number case. Let to be estimated. It is easy to verify that, according to (2.13), the only non-zero quadratic projection of is the projection over the term . As a consequence, the structure of the predictor must be as following:

(3.5)

with

. The term makes sure that the predictor is not biased. The stochastic observability concept introduced in [2] is used here to determine the dimension of the predictor, that is the dimension of . In the general case, the variable is stochastically observable if the following relation:

(3.6)
is one to one. To verify this property, let us compute first the conditional expectation:

(3.7)
According to (2.11), one may write:

(3.8) In the same way, one has: Let:

(3.9)
It is clear that the relation is one to one. As a consequence, the dimension of must be at least two. Note that it is equivalent to consider and the more natural following set:

(3.10)

Note moreover that if , the two components of are identical and the dimension of the quadratic predictor is reduced to .

+ + + = z t R ∈ B t z t Π t t ∆ † + Π t t ∆ + E y τ 1 y τ 2 * Π t t ∆ + [ ] τ 1 ∀ τ 2 t ≤ , , ( )       → E y t y t * Π t t ∆ + [ ] E x t k x t l * Π t t ∆ + [ ] k l , 1 = n ∑ 2 R + = E x t k x t k * Π t t ∆ + [ ] k 1 = n ∑ E x t k x t l * Π t t ∆ + [ ] k l ≠ n ∑ 2 R + + = E y t y t * Π t t ∆ + [ ] E x t t ∆ + k x t t ∆ † + k* Π t t ∆ + [ ] k 1 = n ∑ 1 a --- E x t t ∆ + k x t t ∆ † + l* Π t t ∆ + [ ] k l ≠ n ∑ 2 R + + = E y t t ∆ - y t t ∆ - * Π t t ∆ + [ ] E x t t ∆ + k x t t ∆ † + k* Π t t ∆ + [ ] k 1 = n ∑ 1 a 2 ----- E x t t ∆ + k x t t ∆ † + l* Π t t ∆ + [ ] k l ≠ n ∑ 2 R + + = Π t x t k x t k * k 1 = n ∑ x t k x t l* k l ≠ n ∑ , = Π t t ∆ + E y t y t * Π t t ∆ + [ ] E y t t ∆ - y t t ∆ - * Π t t ∆ + [ ] , { } → z t Π t Ψ t x t k x t k * k 1 = n ∑ x t k k 1 = n ∑ 2 , = n 1 = Ψ t 1
4 Detection test.

The detection test used here takes one's inspiration from the linear gaussian case. Recall that in that classical case, one constructs the optimal predictor of the signal under each hypothesis . One then computes the likelihood of each hypothesis that has the following expression:

(4.1)
The decision process is then defined as: if , the hypothesis is chosen.

Our test is based on the innovation processes generated by the signal square module prediction. Let be the one step predictor of under hypothesis " modulations in the signal". The test variable is imposed to have the following form:

(4.2) where and are parameters devoted to normalize the contrast within the two hypothesis. A natural way to satisfied this condition is to impose:

(4.3)
assuring that the expectation of the distance within the test variable under each hypothesis increases linearly with the number of data treated. The test is then reduced to:

and

(4.4)

The constraints (4.3) leads to the following values of :

(4.5)

where:

(4.6)
Note that all the variances may be computed a priori since the system is bilinear.

Simulation results.

In all the following simulations. the signal to noise ratio is fixed to .

x ˆn k 1

+ ( ) t ∆ k t ∆ H n H n L t n y k 1 + ( ) t ∆ x ˆn k 1 + ( ) t ∆ k t ∆ - ( ) 2 E y k 1 + ( ) t ∆ x ˆn k 1 + ( ) t ∆ k t ∆ - ( ) 2 [ ] --------------------------------------------------------------------------------- E y k 1 + ( ) t ∆ x ˆn k 1 + ( ) t ∆ k t ∆ - ( ) 2 [ ] ( ) log + k 0 = t t ∆ ⁄ ∑ = L t n L t m 2 < n φ ˆnt t t ∆ + φ t n S t S t y k 1 + ( ) t ∆ 2 φ ˆn1 k 1 + ( ) t ∆ k t ∆ - ( ) 2 α k y k 1 + ( ) t ∆ 2 φ ˆn2 k 1 + ( ) t ∆ k t ∆ - ( ) 2 β k + + k 0 = t t ∆ ⁄ ∑ = α k β k E S t H n 1 [ ] t = E S t H n 2 [ ] t - = S t 0 < n 1 → S t 0 > n 2 → α k β k , α k 2 J k 23 - J k 22 + J k 33 J k 32 - ----------------------------------- = β k 1 J k 23 - α k J k 33 - = J k ij E y k 1 + ( ) t ∆ 2 φ ˆni k 1 + ( ) t ∆ k t ∆ - ( ) 2 H n j [ ] = J k ij 10dB
5.1 Detection within 1 and 2 modulations.

The first step of construction of the test is the computation of the predictor of the signal square module. One shows on figure 1 the estimation of the signal square module under each hypothesis when only one signal is present with power equal one. Il appears clearly that the estimator which deals with , the good hypothesis, is much more accurate than the estimator which deals with . The situation is inverted when two signals are present as it appears on figure 2. The test variable as defined by (4.2) is presented on figure 3. The dotted lines are the expectations of under each hypothesis. One may oberve that the detection is very fast.

S t S t

  φ

figure 1 :

 1 figure 1: Module estimation -1 modulation.


A Monin
email: monin@laas.fr

Submarine Floating Antenna Model for LORAN-C Signal Processing

The goal of this paper is to propose an electromagnetic model of the Àoating antenna used by submarines for LORAN-C radionavigation and Very Low Frequency (VLF) communications. We present a description of the Àoating antenna and analyze its electric characteristics. The current in the horizontal wire induced by the lateral wave is then computed and we suggest a model that yields to a non-rational transfer function. A parameter estimation algorithm based on parallel-extended Kalman ¿lters is then proposed and validated on real data.

Introduction.

The antenna used by submarines, for LORAN-C radionavigation and Very Low Frequency (VLF) communications, is made of a long insulated wire ( ' .ff 4) dragged by the submarine when it is submerged. A signi¿cant part of this wire is located near the surface. For VLF signals, the wave draws into seawater and is picked up by the antenna. The behavior of such an antenna is often reduced, in practice, to a simple model. One generally just considers that the current phase is inverted when the wave path comes from the front to the back of the carrier. Moreover, since the phase center location of the antenna (virtual point where the wave is supposed to be picked up) is not accurately known, the antenna is only used for radio-navigation when the incident angle of the path wave with the axis of the antenna is near zero (in practice D ). In order to improve the location accuracy with LORAN-C signal processing, it is necessary to develop a relevant model of such an antenna. To achieve this goal, it is necessary to deal with the electromagnetic properties of these kinds of antenna.

The electric ¿eld induced in such an antenna is related to the lateral electromagnetic waves theory (see [START_REF] Ronold | Lateral Electromagnetic Waves -Theory and Applications to Communications[END_REF] for a complete bibliography). It may be viewed as a diffraction effect due to the differences between the wave velocity in the air and seawater and may be analyzed using optical Fresnel formula [START_REF] Andréani-Corbabé | Etude d'Antennes sous-marines pour la réception de signaux electromagnétique de très basse fréquence[END_REF]. This kind of waves has been fully studied and the results have been applied to many ¿elds including geophysical exploration, communications and remote sensing.

The paper is structured as follows. After some electromagnetic preliminaries dealing with the propagation properties of classic electromagnetic waves, we shall present a physical description of the Àoating antenna and analyze its electric characteristics. We then evaluate the electric ¿eld properties of a vertical dipole transmitting a wave travelling along the boundary between sea and air in reference to the LORAN-C transmitter. Note that these ¿rst paragraphs are essentially derived from [START_REF] Ronold | Lateral Electromagnetic Waves -Theory and Applications to Communications[END_REF]. We then compute the main current in the horizontal wire induced by the lateral waves. It is shown that the computation of this current is obtained by integration of elementary currents along the axis of the wire. The model involved exhibits a non-rational transfer function between the current induced in the transmitter and the current in the receiver. We then develop a parameter estimation algorithm based on the use of parallel-extended Kalman ¿lters. Finally, the model is validated with real data.

Electromagnetic preliminaries.

Let us recall that the Fourier transform of the electric ¿eld HE{> w in all homogeneous isotropic media has the following expression:

e HE$> { ' e DE$h nE$ ' $ u E% m $ ( 1
)
where $ stands for the Fourier variable, nE$ is the wave number, is the permeability, % is the permittivity and is the conductivity of the medium (m ' ). e D is an integration constant resulting of sources. Performing the inverse Fourier transform, one gets the general solution of the Maxwell equations, that is:

HEw> { ' 2 ] e DE$h g$
Let us consider a vertical wire source with current L along the axis of the wire. If the wavelength is much smaller than the wire length o, one may assume that the antenna is equivalent to a dipole with electric moment L o [3] [4]. When observed in a direction that is orthogonal to the antenna axis, the Fourier transform of the electric ¿eld radiated at range , in the far-¿eld region, can be written:

e HE$> ' m$ 2 h e LE$o (2)
where e LE$ stands for the current Fourier transform. Conversely, for reception, the current induced in a wire antenna of length k by the impressed electric ¿eld e HE$> can be represented by:

e LE$ ' 2 m$ e HE$> k (3)
The characteristic impedance of such an antenna, de¿ned by ] ' Ho@L, is:

] ' m$ 2
3. Floating antenna.

Description.

The Àoating antenna behavior is described by the theory of lateral electromagnetic waves [START_REF] Ronold | Lateral Electromagnetic Waves -Theory and Applications to Communications[END_REF]. Suppose that an electromagnetic source is lying on the boundary between air and seawater. Three waves are then generated. In both media, there is a direct wave travelling at velocities depending on their electromagnetic characteristics. In air, this direct wave is the only one that reaches the receiver. In seawater and for frequencies near ff !O3, as it is the case for LORAN-C signal, this direct wave is extremely attenuated. Its power is divided by two after a course of Hf U4 long. A third wave, called lateral wave, travels along the surface of the water in air which continuously generates waves that propagate into the water at the critical angle determined by the ratio of velocities of each medium [START_REF] Andréani-Corbabé | Etude d'Antennes sous-marines pour la réception de signaux electromagnétique de très basse fréquence[END_REF]. This ratio is generally suf¿ciently great so that the angle of refraction in water is practically bf . The receiver in the water observes mainly this lateral wave since it travels a long distance in air and then a short distance in water.

The Àoating antenna used by submarines, for Very Low Frequency (VLF) telecommunications and LORAN-C transmitting, consists of a conductor of length k enclosed in an insulating sheath located at a small depth g in seawater [START_REF] Ronold | Lateral Electromagnetic Waves -Theory and Applications to Communications[END_REF](¿gure 1). It is driven at one end against the impedance ] of the receiver. At the other end, the insulated conductor is terminated in second impedance ] . Let ] be its characteristic impedance. In practice, a maximum current is maintained along the antenna when it is terminated either in a low impedance (] ] ) or in its characteristic impedance (] ' ] ).

The current along the insulated horizontal wire is a generalized travelling wave with wave number n depending on the antenna properties. Let n and n be the waves number of medium and 2 respectively, de¿ned according to [START_REF] Ronold | Lateral Electromagnetic Waves -Theory and Applications to Communications[END_REF]. In region (air), the medium is lossless and its wave number is de¿ned by:

n ' $ s % ' $ f ( 4
)
where f is the light velocity in air. In region 2 (seawater), the conductivity of the medium is not zero and the wave number has the following regular expression:

n ' n u % m % $ (5)
Note that the relative permittivity of seawater % is near Hf when its conductivity is approximately =b 5@ 4. For frequencies near ff !O3, the % term is negligible in the square root of (5). In fact, @E% $ * .=f2 f which is much smaller than Hf. Therefore, the wave number in region can be reduced to:

n ' s h s $ The wave number of the Àoating antenna has the following general expression [START_REF] King | Antennas in Matter: Fundamentals, Theory and Applications[END_REF]:

n ' s } | with: | ' m2n $ *?Ee@d } ' } m$ 2 *? e d n } n }
where n ' $ s % is the wave number of the dielectric sheath (outer radius e), d is the radius of the insulated wire, } is the internal impedance per unit length, } is the series external impedance per unit length and } is the series mutual impedance per unit length. The expression of } can be approximated, depending on the value of n ' s m$ > the wave number of the conductor:

mn dm ? 2 ', } * d mn dm A D ', } * mn 2d In our case, for copper ( ' D=SD f 5@ 4), at frequency i ' ff !O3 and with a conductor radius of 44, one has mn dm ' S=SH A D. Consequently, contrary to the case studied in [START_REF] Ronold | Lateral Electromagnetic Waves -Theory and Applications to Communications[END_REF], the leading approximation of } is:

} * E m 2d
u $ In the same way, if mn em ? and mn gm ? , } and } can be approximated by:

} * $ e m$ 2 E*?E 2 n e } * $ e m$ 2 E*?E n g
with ' f=D..2. Using all these expressions, the wave number of the insulated antenna is:

n ' n v n *?E E d s $ n *?E 2 mn em n *?E mn gm 2 n mE d s $ n
and its characteristic impedance is:

] ' u } | ' $ n 2n *?E e d (6)

Evaluating the incident ¿eld

The transmitting LORAN-C antenna is equivalent to a vertical dipole in air as described in Section 2. The wave travels along the sea surface acts as a waveguide since its conductivity is not negligible. Therefore, the axial component of the electric ¿eld near the boundary of the two media is not quite zero as it is the case in homogeneous media. The transverse and axial components of the far ¿eld have the following expressions (the subscript is composed of one number related to the region and of one letter, for the axial component and } for the transverse one):

H E ' n I E> n > n H E ' n JE> n > n ( 7
)
where stands for the distance between the source and the receiver. Functions I and J will be de¿ned in the sequel. In water, at depth }, these components are:

H E ' H Eh H E ' n n H Eh
Note that since n @n * f in practice, it is clear that the transverse component in water (H E) is quite negligible with respect to its value in air, as describe in section .1. Conversely, the axial component in sea is the main component of the entire ¿eld, only attenuated due to the imaginary part of n . Functions I and J are de¿ned by:

I E> n > n ' $ 2 Ei Eh n i Eh JE> n > n ' $ 2 Ej Eh n i Eh
with i E ' mn @ @ , j E ' mn @ @ nmn @ , i E ' n @n s @En IEn @E2n and n ' n E n @E2n . The function I is the Fresnel integral de¿ned by:

IEx ' 2 E m ] h s 2w gw
Functions i E (solid line), j E (dashed line) and i E (dash-dot line) are shown in ¿gure 2. Clearly, for ranges smaller than Dfff !4, the magnitude of i E is quite negligible with respect to other terms. Moreover, for such ranges, i E * j E. Note that for greater frequencies, say such that n mn m @n , the term @ in the expression of i E is dominant since the term @ is cancelled by the Fresnel integral. Furthermore, for ranges greater than f !4, the term @ is also negligible with respect to the term @. Finally, for frequencies near ff !O3 and for ranges within f !4 and Dfff !4, the actual formula for functions I and J are: Using the expression of n in (5), the axial component of the incident electric ¿eld in water can be written:

I E> n > n * JE> n > n * m$ n 2
H E> } * h 2 s % $ h h (8)
The transverse component in air is:

H E * m$ 2
h Note that H has the same expression as the electric ¿eld radiated by a wire source in a homogeneous medium [START_REF] Andréani-Corbabé | Etude d'Antennes sous-marines pour la réception de signaux electromagnétique de très basse fréquence[END_REF]. Finally, the ratio between the transverse component in air and the axial component in water at depth } is:

H E> } H E ' s $ s % h h (9)
5. Evaluating the radiated ¿eld.

Current along the insulated wire.

Let us de¿ne Y the electromotive force of the generator. It can be shown that the current is distributed like the one in a generalized transmission line. Its Fourier transform has the following expression:

e LE{ ' m e Y E] n ] t?En Ek { m
ULtEn k m where is the terminal function de¿ned by:

' UL| E ] ] (10)
The current in the load ] ({ ' f) is related to the current at { by: Let us note ! the angle between the propagation direction and the {-axis of the antenna as illustrated ¿gure 3. The radial and transverse components of the radiated electric ¿eld due to the horizontal dipole g e LE{ in the sea (medium 1) at distance from the origin are:

g e H E$> > ! ' m$ 2 n n h ULtE!g e LE{ g e H E$> > ! ' m$ 2 E n n h ULtE!g e LE{
where n and n are de¿ned according to ( 5) and ( 4). Note that, with frequency near ff !O3, one has:

n n ' u % $ ' =b f
This leads to consider only the transverse component of the ¿eld in the air. Let be the distance to the receiver referred to the point of reference { ' f. Assume that k. At distance { from the origin, the distance to the receiver is:

' { ULtE!
The contribution of the elementary moment e LE{g{ to the entire ¿eld is then:

g e H E$> {> ! ' e H E$> > !h e LE{ e LEf g{ ( 14
)
where e H E$> > ! is the electric ¿eld due to an horizontal electric dipole with moment e LEf located at { ' f, that is:

e H E$> > ! / m$ 2 n n h ULtE! e LEf ( 15
)
Let us introduce the effective length of the antenna de¿ned by:

k E$> ! / ] h e LE{ e LEf g{
The integration of (14) along the axis of the antenna leads to:

e H E$> > ! ' k E$> ! e H E$> > ! ( 16
)
since it is assumed that the amplitude of the incident ¿eld is approximately constant over the length k of the wire.

Formula (12) and ( 13) allow the computation of the effective lengths, that is: Matched load.

k E$> ! ' ] h g{ ' mEn n ULtE! E h Low-impedance termination. k E$> ! ' ] ULtEn Ek { ULtEn k h g{ ' En n ULt E! ULtEn k Emn ULtE!Eh
ULtEn k n n t?En k 5.3. Current in the wire for reception.

Let us consider the LORAN-C transmitting antenna for reception. Let e H E$> > ! be the incident ¿eld de¿ned by ( 16). According to (3), the current in the load is:

e L E$ ' 2 m$ o e H E$> > ! ( 17
)
where o stands for the receiving antenna length. Using (15), ( 16) and ( 17), one has:

e L E$ ' o n n h ULtE!k E$> ! e LEf
Applying the reciprocal theorem to this couple of antennas, if the current in the LORAN-C transmitting antenna is noted L Ew, the current in the receiving buoyant wire antenna is then:

e L E$> ! ' o n n h ULtE!k E$> ! e L E$
where:

L Ew ' DE w w h t?E2i w ( 18
)
with w ' SDv and i ' ffnK}.

In conclusion, with ( 5) and ( 4), omitting the attenuating factor @ and the pure propagation delay h , the transfer function between the transmitted and received current is: Matched load.

KE$> ! ' s $h ULtE! En n ULtE! E h (19)
Low-impedance termination.

KE$> ! ' s $h ULtE! En n ULt E! ULtEn k En ULtE!EULtEn k h n mn t?En k (20)
Note that these expressions are only de¿ned for $ A f. As the response of the system is supposed to belong to U, one must have, for all $ ? f, KE$ ' K E $, where K stands for the complex conjugate of K. This completes the de¿nition of KE$> ! for all $ 5 U.

Identifying the antenna parameters.

Problem setting.

In order to calculate high quality parameter estimates, it is necessary to combine measurements from several transmitter stations or several measurements from a single transmitter. The variables in this context are: antenna characteristics (see ¿gure q may be viewed as the antenna index near the nominal frequency (ff !O3) and as an exponential attenuation coef¿cient. We have veri¿ed that with such an approximation, the distortion induced in the waveform is lower than I. Moreover, if the deployed length of the antenna is a parameter which is theoretically known under the operational conditions, the part Àoating on the surface of water is a priori dif¿cult to determine. Therefore, the working length k will be considered as an unknown parameter. The angle ! is directly computed from ! ' where stands for the submarine course and is the azimuth of the transmitter n.

In short, the state variable to be estimated is reduced to:

{ ' dD > ===> D > o
and the parameters are: ' dq > > ko

Modelling.

The reference signal can be represented in the following way:

uEw> > ! ' I EKE$> > !IEvEw> w 5 U (22)
where I stands for the Fourier transform, vEw is the ideal LORAN-C wave de¿ned by (18) and KE$ is the de¿ned by ( 19) or (20). As only one station is observed at each time, it is convenient to consider that the sampled observation process | > w 5 Q is a vector of U de¿ned by: | ' where {w is the sampling period, y stands for the output white noise and is the propagation time delay between the station n and the receiver. It is computed from the GPS data and the station location using the so-called Salt formula [START_REF]LORAN-C User Handbook[END_REF]:

' G f n 2=S f G n f=fe G e= f
where G is the geodesic distance from the station n. Note that we do not take into account here the sky wave interference [START_REF] Last | Effect of skywave interference on the coverage of LORAN-c[END_REF] because we use only stations close to the receiver in the identi¿cation procedure (G * fff !4).

As the signal amplitude can Àuctuate with time, their variations can be represented by:

D ' D n D ( 24
)
where D is a white noise (possibly Gaussian). Similarly, the Àuctuation of the ampli¿er delay is introduced in view to compensate possible errors and variations of GPS location data:

' n ( 25
)
where is also a Gaussian white noise.

Parameter estimation.

The estimation algorithm proposed here uses an approximation of the a posteriori probability density function based on a sample of the parameters to be identi¿ed. Let us consider the following discrete time dynamic model:

{ ' i Ew> { > n jEw> { > z | ' kEw> { > n y
where { is the state vector, is a time-invariant vector of parameters, | is the output and z and y are independent Gaussian white noises. If this system is locally linear with respect to { , for each value of , one can build the Extended Kalman Filter (EKF) of { . Recall that the EKF use is relevant if [START_REF] Andrew | Stochastic Processes and Filtering Theory[END_REF]: i) the initial uncertainty of the state is small enough compared to the ¿eld of linearization validity ii) the noises have standard deviations small enough to ensure that the standard deviation of the Kalman ¿lter lies on this ¿eld of validity.

This leads to the following approximation of the probability density function:

sE{ m| > ' JE{ ( { E> S E
where { E is the EKF of { and S E is its covariance matrix. Here, JE{( {> S stands for the standard Gaussian distribution with mean { and covariance S . Suppose that at time w , the conditional probability density function for is approximated by: sEm| '

[ E ( 26
)
where E stands for the centered on . Using the Bayes rule, one has:

sEm| ' sE| m| > sEm| (27)
where is a normalizing constant independent of . Replacing sEm| in ( 27) by its value de¿ned in ( 26), one has:

sEm| '

[ sE| m| > E
Because | is supposed to be "almost" Gaussian, one has:

sE| m| > ' JE| _ E > > _ { E S E KEw> > _ { E n U ( 28
)
where KEw> > { ' and U is output noise variance. As a consequence, at time w, the conditional probability density function for has the following approximation:

sEm] ' [ E
The normalized weights are de¿ned recursively using (28):

' sE| m| > S sE| m| > (29)
The minimum variance parameters estimation is then computed using:

e ' [ (30)
and the covariance matrix by: S '

[ E e E e (31)
Recall that this approximation is relevant only if standard deviations t S E are suf¿ciently small with respect to the ¿eld of validity of the linearization.

In fact, this ¿ltering procedure can be applied to our case. Indeed, the dynamic model (24) (25) is clearly linear Gaussian and the amplitudes D appear linearly in the observation model (23). Moreover, the observation model can be linearized with respect to as long as the standard deviation is low compared to the signal period (f t).

The identi¿cation is carried out by positioning Q samples of on a ¿xed grid. Each one of its nodes represents a set of parameter ' d > q > k o. For each sample , one computes the EKF that estimates the state vector { ' dD > ===> D > o using the following output matrix:

K E ' CkEw> { > C{ ' 5 7 u Ew{w > === f D m === === === === f == u Ew{w > D m 6 8
The identi¿ed parameters are then computed using (28), ( 29) and (30). The evaluation of the standard deviation of the estimate given by ( 31) is just used to survey the consistency of the identi¿ed parameters.

7. Experimental results.

Data set.

We use a set that has been collected by a civilian boat freighted by the General Army Direction of France (DGA). The LORAN-C signal was sampled at eff !O3 and time-stamped to a very high accuracy (a few nanoseconds), thanks to the use of a cesium clock. The position reference of the carrier was acquired using a classical GPS receiver. Course and speed are also recorded. The duration of each record is f 4?. The identi¿cation had been which the carrier (distance from transmitters near fff !4). The distances and incidences of the transmitters used are presented in table 1. In order to reduce the quantity of calculations necessary to the identi¿cation, it is convenient to initially carry out the summation of the samples on a burst of eight pulses, the period of the summation being 4t. The signal obtained can then be added up again over several Group Repetition Interval (GRI) periods (* f= t). This reduces the signal-to-noise ratio (division of the variance of the additive noise by the number of pulses added up). Note that, as the carrier location is known and the signal dating very accurate, it is possible and appropriate to dismiss all pulses that are supposed to be spoiled by some other LORAN-C transmitter. It should not however be forgotten that, in the calculation, the movement cannot be compensated, which yields to deform the signal. This deformation is connected with an integration of wave for one length of time corresponding to the distance covered during the summation. As an example, for a carrier navigating at D 4@ , after D t, it traversed approximately D 4, which corresponds to a time of integration of approximately f=fD t. If we compare the starting signal with the integrated one, it is observed that the induced distortion is lower than f=fI. For theses reasons, this duration of summation (D t) had been appointed. At this stage, each D t, we have a pulse at one's disposal which looks like the one shown ¿gure 5.

Algorithm.

Step 1 Sampling the parameter space. Step 2 Reference signals.

For each sample, one calculates the signal of reference as being the result of the ¿ltering of the LORAN-C signal through the de¿ned transfer function (20). The computation is achieved using numerical Fourier transform.

Step 3 Signal summation.

One carries out here the summation of the impulses of the signal, within a GRI, then between GRI as described in .2. The impulses re-entering in conÀict with others are eliminated on a power criterion. The number of pulses proceeded is collected. The variance of the additive noise divided by this number of summations gives the new variance of noise used by the EKFs. Note that this noise variance can be easily derived by estimating the power of the signal non-contaminated by LORAN-C pulses.

Step 4 Ampli¿er delay and amplitudes estimation.

The EKFs estimates and iD j, conditionally with the value of the parameters . The skywave is not taken into account in this estimate, the studied stations being suf¿ciently close to the receiver (time of arrival of the sky wave A SD t) [START_REF] Last | Effect of skywave interference on the coverage of LORAN-c[END_REF].

Step 5 Weight computation (29)

Step 6 Parameter estimation (30).

Step 7 Convergence test (31).

If the signal-to-noise ratio is such that the estimate has not the accuracy required, one carries out a translation of the point of reference in conformity with dynamics of the carrier (given course and speed) and one turns over in 3. The variances of the Kalman ¿lters are brought up to date in accordance with the prediction errors (variances of and D in the dynamic model ( 25) and (24). Otherwise, the algorithm provides the estimated parameters.

Identi¿cation.

Experimental results have been achieved using data described in Section .1. In view to test the algorithm, the parameter estimation algorithm used records 1 and 2 and transmitters 2 (incidences closest to f and Hf ), 1+2 and 1+2 +3. The results are presented in table 2.The induced radiation pattern (computed from (20) at ff !O3) and the incidences of transmitters used for identi¿cation are shown ¿gure 6. In ¿gure 7, we show the wave obtained with our identi¿ed model for an incidence angle near (transmitter 2 -record 1). It is compared to the wave obtained with the whip aerial model. It appears clearly that our model is able to ¿t the data with much more accuracy than the whip aerial model does. Note that, however, with this antenna, the error involved in the location is near Hff 4. On the other hand, as shown in ¿gure 8, the whip aerial model remains accurate when used with incidence angle near f (transmitter 2 -record 2). Indeed, the wave deformation induced by the antenna is quite negligible in that case, which explains that this Àoating antenna is only used in practice in such a situation. The magnitude of the current in the wire de¿ned by ( 20) is theoretically very dependent on the incidence ! of the wave, as it appears on ¿gure 6. Table 3 presents the comparison between the magnitude of the signal estimated for some records and their theoretical values computed from [START_REF]Groundwave propagation curves[END_REF]. In the third column, the theoretical attenuations at ff !O3 (nKE$> !n with $ ' 2i ) are displayed illustrating variations of the antenna gain for some incidences. In the last column, the difference between the amplitudes estimated and their theoretical values are shown. It turns out that these differences are lower than 2dB in practice. One may conclude that the radiation pattern produced by our model is close to reality.

Model validation.

Positioning accuracy.

Finally, the relevance of our model has been tested by the comparison of the propagation delay estimated with both models (Àoating antenna and whip aerial). To achieve this goal, we had used several records with distances greater than those had used for identi¿cation. The results are presented in table 4. It appears clearly that the whip aerial model is only relevant for incidences lower than f when our model is accurate for all incidences.

Conclusion.

We have presented a new way for modelling the behavior of the wire antenna used by submarines. The model leads to an important variation of the waveform when used with different incidence angles. It has been shown that the modelled waveforms match real data with a very good accuracy. Moreover, the variations of the antenna gain with the incidences are well treated by such a model.

¿gure 1 :

 1 Floating antenna model 3.2. Wave numbers.

¿gure 2 :

 2 Variation of attenuation terms with distance from the source.

e 11 )

 11 LE{ ' e LEf t?En Ek { m t?En k m (The behavior of the antenna depends on the value of the terminal function : Matched load. If ] ' ] , then ] @] ' . According to (10), one has ' 4. The current behaves like a travelling wave and (11) is simpli¿ed as follows: e LE{ ' e LEfh (12) Low-impedance termination. If ] ] > one has: ' u n mE 2 n { with mu m and m{ m . In that case, the wave is partially reÀected and the current in the antenna becomes: e LE{ ' e LEf ULtEn Ek { ULtEn k (13) 5.2. Radiated ¿eld. The antenna may be viewed as a continuous set of dipoles with elementary electric moments g e LE{ ' e LE{g{. ¿gure 3: Top view of the antenna

  4): d: wire radius e: sheath radius g: submersion depth -% : dielectric permittivity k: antenna length -! : angle between the propagation direction from the transmitter n and the axis of the antenna D : amplitude of the signal from the transmitter n ' ===Q : ampli¿er delay ¿gure 4: Insulated antenna Concerning the antenna characteristics, in order to reduce the number of parameters, it is convenient to consider approximation of n E$ the actual frequency of the signal treated: n E$ * Eq m

u

  Ew> ' uEw{w> > !

5 :

 5 Signal after summationOne de¿nes the searching interval for all parameters, that is dq > q o, [ > o and dk > k o. A mesh-grid is then constructed with steps de¿ned by: {{ ' E{ { s Q where Q denotes the number of samples used. The standard deviations for is set to t ( f t) making sure the convergence of EKFs. The initial means of f D are set to their theoretical values computed using[START_REF]Groundwave propagation curves[END_REF]. The standard deviations of the amplitudes D are set to half of their theoretical values.

7 :

 7 Comparison wire antenna / whip aerial -! ' .2 Consequently, it allows extending the use of transmitters with any azimuth, without change of course, as it is necessary at present time. Moreover, one can use at the same time several transmitting stations, which is known to improve the global location accuracy [10] [11].

8 :

 8 Comparison wire antenna / whip aerial -! ' 2

table 1 :

 1

	Record 1	Record 2
	Transmitter 1 1090 km / 156.1	1040 km / 28.6
	Transmitter 2 975 km / 172.2	930 km / 11.6
	Transmitter 3 1225 km / 140.2	1170 km / 45.2

Distances and incidence angles 7.2. Preprocessing.

table 2 :

 2 Parameters estimated

	Stations	q		k
	Transmitter 2	4=68 3=34 3=6; 3=34 9;5 p 9 p 3=4: v 3=37 v
	Transmitter 1 + 2	4=69 3=34 3=73 3=339 9;4 p 6 p 3=4< v 3=36 v
	Transmitter 1 + 2 + 3 4=6: 3=34 3=73 3=34 9;4 p 5 p 3=4: v 3=36 v
	Mean	4=69	3=6<	9;4 p	3=4; v

Acknowledgements.

This work was supported by the DGA (General Army Direction -France) and the DIGINEXT company (France).


André Monin

ARMAX Identi¿cation Via Hereditary Algorithm

Keywords:

In this paper, we limit ourselves to the most commonly models used, say the Output Error (OE) model and Auto Regressive Moving Average with exogenous input (ARMAX). In both cases, we treat the case multi-input, single-output (MISO). Our aim is to give the parameters of the model under lattice form, in reference to its numerical robustness allowing to deal with great state space dimension (A 2f). This calculation is achieved considering time varying predictor of the output based on the predictor derived at the previous step. It leads to an hereditary algorithm introduced in previous works [1].

I. INTRODUCTION.

The paper concerns the problem of linear system identi¿cation. The stochastic linear identi¿cation is an issue that has been fully investigated for many years and many algorithms have been designed to this purpose (see [START_REF] Ljung | cation: Theory For The User[END_REF] for bibliography), the most recent being the subspace identi¿cation method [START_REF] Van Overchee | Subspace Identi cation of LInear Systems: Theory, Implementation and Applications[END_REF]. In some previous work [START_REF] Monin | Exact ARMA lattice predictors from autocorrelations functions[END_REF] [1], we had derived some new realization and identi¿cation algorithms based on a hereditary computation in the single output AutoRegressive Moving Average (ARMA) case. We extend here these algorithms when the system has exogenous non-stationary multi-input. The main difference lies in the fact that the lattice parameters, involved in the computation of the system parameters, are no more stationary, as it was the case in [START_REF] Monin | ARMA lattice identi cation: A new hereditary algorithm[END_REF].

A. Problem setting.

The aim of linear identi¿cation is to ¿t some linear stochastic realization to an experimental set measurement of the input/output of some unknown system. We restrict ourselves here to the set of AutoRegressive Moving Average with exogenous input (ARMAX) models:

| @ [ d | . [ [ e x . [ f y . y (1)
where | stands for the (measured) output, ix > n @ 4===qj stands for the (measured) inputs, y stands for a white noise with unknown variance and id > l @ 4===qj ^ie > l @ 3===q 4> n @ 4===sj ^if > l @ 4===qj stands for the set of unknown parameters to be identi¿ed. Note that, for convenience, we ¿x the same order q for the AR, the MA and the X part.

Referring to the Gaussian case, the accuracy of the model derived is commonly measured through the one-step-aheadprediction error of the output, that is:

M @ 4 w [ | e | (2)
A. Monin is with the LAAS-CNRS, Toulouse, France. E-mail: monin@laas.fr.

where e | stands for the predictor of the output. If one may assume that the zeros of the MA part are all inside the unit circle, the output predictor e | is [START_REF] Ljung | cation: Theory For The User[END_REF]:

e | @ [ d e | . [ [ e x . [ n h | (3) where h | @ | e |
is the prediction error and n @ d f . Note that it is assumed that the whole data record over the interval ` 4> w`is available. Obviously, in practice, only data over the interval ^4> w`are known. The common way to get round this dif¿culty is to use windowed data. Most often, one just assumes that the past data, for w ? 4, is zero. This is equivalent to consider that the initial condition, in a state space representation, is zero, which is not true in most cases. The main drawback of this assumption is that, for short experiments, it is more dif¿cult to ¿t data to the model with non-zero initial conditions. The true model for ¿nite data record would be:

;w @ 4===q> | @ { . [ d | . [ [ e x . [ f y . y
where i{ > w @ 4===qj stands for the q initial conditions of the system. For w A q, 1 is then fully de¿ned. To represent the initial condition effect, it is convenient to introduce an extra virtual impulsive input x de¿ned by x @ 4 and x @ 3> ;w 9 @ 4. The complete model is then:

| @ [ d | . [ [ e x . [ f y . y
with, e @ { > ;l @ 3===q 4. With this modi¿ed model, all data must be assumed zero ;w ? 4.

B. Computing the parameters via hereditary projection.

Recall that the criterion 2 is not quadratic with respect to the parameters, as it is the case for ARX models (f @ 3> l @ 4===q). Indeed, if one develops the dependence of the criterion with respect to the parameters d for instance, the AR part appearing in 3 involves an in¿nite development, which is not explicitly tractable. Many approximating algorithms have been developed to compute the model parameters as the so-called Extended Least Square [START_REF] Benveniste | AR and ARMA algorithms of levinson type: An innovation approach[END_REF] [6] [START_REF] Ljung | cation: Theory For The User[END_REF] or Gauss-Newton algorithms [START_REF] Shynk | Adaptive IRR ltering[END_REF] [6] [START_REF] Michaut | Méthodes Adaptatives Pour Le Signal[END_REF]. In none of these approaches, transient behavior is optimized. In this paper, we deal with a time varying approach preserving a transient convergence similar to that of [START_REF] Monin | ARMA lattice identi cation: A new hereditary algorithm[END_REF]. Let us denote e | > @ 4===w 4 the optimal predictor that minimizes the criterion M . The parameters have been obtained from the knowledge of the data record i| > x > @ 4===w 4> n @ 3===sj. At time w, the new predictor trajectory is computed as follows:

; @ 4===w> e | @ [ d e | .
[

[ e x

That means that the whole trajectory from 4 to w is recomputed at each step. In such a time varying scheme, the criterion is quadratic with respect to id > e j (but not with respect to id > e > @ 4===w 4j previously computed). Therefore, the optimal parameters can be computed using classical projection formula. On the other hand, it leads to a hereditary algorithm. That means that the memory needed for computation is linearly growing with time due to the updating of the overall least-square criterion as data accumulate. In other words, the non-linearity complexity is handled via hereditary computation.

II. GENERAL LEVINSON FORMULA.

Let us note the experimental mean of a process K 5 U de¿ned over the interval ^4> w`as H ^K `@ S K . De¿nition 1: The projection of a process [ 5 U > @ 4===w over the process \ 5 U > @ 4===w, is the process ] / S m\ `@ N \ where

N H ^[ \ `+H ^\ \ `, .
It is clear that this projection yields the process that minimizes the mean square error H ^n[ ] n `.

Let us de¿ne ] 5 U

as:

] @ ^x > ===> x > e | > h | `> @ 4===w (4)
Assume that ; ? 4, ] @ 3. The linear predictor e | is then computed using the projection of de¿nition 1:

e | @ S ^| m] > ===> ] ` (5)
Suppose that ; ? 4, [ @ \ @ 3. One then has: Proposition 2: S ^[ m\ `@ S ^[ m\ Às it appears in 5, the linear predictor is computed as a linear projection over the space spanned by i] > ===> ] j. The Levinson formulae are devoted to recursively compute an orthogonal basis of this space according to the inner product de¿ned by k[ m\ l / H ^[ \ `. Let us de¿ne the family iY > ===> Y j (called backward residuals) to be the orthonormalized basis of the space spanned by i] > ===> ] j. According to this de¿nition, one has:

H ^Y +Y , `@ 3> ;l 9 @ m H^Y +Y , `@ L> ;l @ 3===q 4

where L stands for the identity matrix. One has: Proposition 3: Backward and forward normalized residuals obey the following relations:

; @ 4===w> ; A A ? A A = Y @ + O , ] Y @ k +L V +V , l +Y V Z ,> ;l @ 4===q 4 (6)
; @ 4===w>

; A ? A = Z @ +O , ] Z @ L V +V , +Z +V , Y ,> ;l @ 4===q 4
where:

O @ H ^] +] , ` (7)
V @ H ^Y +Z , `> ;l @ 3===q 5 [START_REF] Michaut | Méthodes Adaptatives Pour Le Signal[END_REF] Moreover, thanks to the orthonormality property of the residuals, the one-step-ahead predictor has the following expression:

e | @ [ h P Y > ; @ 4===w (9)
where h stands for the vector of U de¿ned by: h @ 4, ;l 9 @ n> h @ 3 and:

P @ H ^] +Y , ` (10)
See appendix VIII for proof.

III. HEREDITARY ALGORITHM.

Suppose that, at previous step (w 4), the whole sequences iY > l @ 3===q 4> @ 4===w 4j (backward residual trajectories) are available. One then computes the matrices iP > l @ 3===q 4j and O de¿ned by 10 and 7 performing a summation over the interval ^4> w`. Note that although this calculus seems implicit, it can be performed using the orthogonality property of the innovation process H ^+| e | ,Y `@ 3. One may then compute the matrices iV > l @ 3===q 4j using recursively equations 18 and 25 of appendix VIII. It remains then to update the residual trajectories. Here is the hereditary part of the algorithm:

For @ 4===w, do:

e | @ S h P Y ] @ x === x e | +| e | , Z @ +O , ] For l @ 4===q 4, Z @ L V +V , +Z +V , Y , For l @ q 4===4> Y @ k +L V +V , l +Y V Z , Y @ + O , ]
Obviously, during the ¿rst steps of the algorithm, the covariance matrices, such as V , are singular since all data prior to @ 4 are assumed to be zero. Therefore, the coef¿cients related to these data are undetermined. In fact, only one new non-zero coef¿cient is derived at each step. That means that one needs at least q +s . 6, steps to assure the regularity of all matrices. There are two ways to solve this problem.

One may use pseudo-inverse procedure when one needs to inverse a singular matrix.

One may also compute previously the rank of matrices by analyzing the algorithm. For example, one has to compute the inverse of the square root of O . It is clear that, if the input-output sequence is not singular, the rank of O 5 U is equal to w, since w 9 +s . 6,. Moreover, one has to compute the inverse of the square root of the matrices D @ +L +V , V ,. If u denotes the rank of V > the rank of D is obviously s . 6 u. Because only one parameter is derived at each step, one needs to wait +s . 6,+l . 4, steps to assure that V has a rank equal to s . 6. More precisely, one has:

-If w 9 +s . 6,+l 4, . 4 or w A +s . 6,l, then un+D , @ s . 6 -If w 5 ^+s . 6,+l 4, . 4> +s . 6,l`, then un+D , @ +s . 6,l w . 4 -If w 5 ^+s . 6,l . 4> +s . 6,+l . 4,`, then un+D , @ w +s . 6,l where un+D, stands for the rank of the matrix D. Note that these procedures do not harm the optimality of the algorithm since it just provides a particular set of parameters among all possible solutions.

IV. COMPUTING THE TRANSFER FUNCTION COEFFICIENTS

For some purposes, it could be convenient to deal with the transfer function related to the linear system involved under lattice form. To achieve this goal, it is necessary to take into account that the lattice parameters are not stationary, since the inputs are not, in general, a perfect white noise process. Indeed, some parameters are resulting from the projection of the inputs over their own past. First note that, at time w, the state of the linear one-step-ahead predictor can be reduced to the sets:

iY > ===> Y j @ i] > ===> ] j
thanks to formula 6. According to 3, if one sets all components of the state to zero except one term to transfer as the output of the ¿lter. More precisely, if one sets: ] @ 3 == 3 4 3 === 3 ] @ 3> ;m 9 @ l one obtains d for m @ s.5, g for m @ s.6, e for m @ n (m 9 @ 3) and { , the l-th component of the initial condition, for m @ 3. The algorithm is then structured as follows:

For l @ 4===q, do Y @ + O , ] where ] is initialized as described above.

For l @ 4===q, do Z @ t Y For n @ 4===q 4 do for @ 4===q 4 do:

-Y @ k +L V +V , l +Y V Z , -Z @ L +V , V +Z +V , Y , @ S h P Y
where stands for the parameter to be computed, that is d , e or { , depending on the value of m.

V. EXTENSIONS TO MORE GENERAL LINEAR SYSTEMS.

In the introduction, we stated that, for convenience, we ¿xed the same order q for the numerators and denominator of the transfer function. In fact, there is no trouble to deal with more general transfer functions.

A. Introduction of delays on the inputs.

If one wishes to introduce time delays on the input, one just has to replace ] de¿ned in 4 with ] @ ^x > x ===> x > e | > h | `> @ 4===w where stands for the delay ¿xed for the n-th input. Note that the ¿rst input x being the virtual impulsive input intended to represent the initial condition effect, it remains unchanged.

B. Dealing with different orders for numerators and denominator.

If one wishes to deal with different orders for numerators and denominator, ¿rst replace q with the maximum of the orders desired. The de¿nition of the predictor e | is then modi¿ed as e | @ S ^| mH ] > ===> H ] `where H is a square matrix that sets to zero the appropriate components of ] . All computations of the algorithms being linear operations, there is no dif¿culty to modify the algorithm with this change.

C. Time-varying systems.

In the case of time-varying systems, if W stands for the tracking-time constant, one just has to modify the de¿nition of the experimental mean as follows:

H ^[ `@ 4 w [ [ > ;w ? W H ^[ `@ 4 W [ [ > ;w W VI. EXPERIMENTAL RESULTS.
As the main goal of this identi¿cation algorithm is to optimize the transient (short experiments), the ¿rst test we made was proceeded in the noise free case. We have veri¿ed, in simulation, that for q-dimensional linear systems with s inputs, the parameters are obtained exactly (in fact with a precision lower that 43 ) after q+s . 6, steps, which is obviously the minimal number of steps required.

We used for experiments some datasets coming from DAISY [START_REF] Moor | DaISy: Database for the Identi cation of Systems[END_REF]. We have compared the performances of our hereditary algorithm with the results obtained by the standard MATLAB identi¿cation toolbox functions: n4sid.m (subspace identi¿cation), armax.m and oe.m (Newton like algorithms). For the ARMAX and OE functions, we added the extra impulsive input described in I, in view to deal with an unknown initial condition An example of results got with the dataset [96-011] with q @ 7 is shown ¿gure 1. We have represented the relative output error t S +| e | , @ S | , computed over the whole dataset (W @ 49;3 samples), when the parameters are computed over only a smaller varying data length (from 43 to 833 per step of 43). In view to clarify the graphical representation, we set this error to be zero when the system derived was unstable. It appears clearly that, for this example, the output error of the hereditary algorithm has almost reached its minimum after only 483 steps. Note that this output error is no more Àuctuating for longest datasets, contrary to those obtained with standard algorithms. Moreover, all systems obtained with our algorithm are stable for a dataset of length A 63 and leads to an error 9 58(.

In view to illustrate the behavior of the algorithm, we have simulated a 9-order linear system with zeros located at i} @ 4> } @ 3=< h{s+ l@7,> } @ 3=< h{s+ 6l@;,j and poles located at is @ 3=< h{s+ l@;,> s @ 3=< h{s+ l@6,> s @ 3=8 h{s+ 46l@65,j. The input is a Gaussian white noise and signal to noise ratio is 45gE. The hereditary estimates of the ¿rst two parameters of the denominator are plotted ¿gure 2 and compared to the results obtained with the standard oe.m MATLAB identi¿cation toolbox function. It appears clearly that the algorithm that uses a classical Newton type method falls on several occasions in local minima. We have derived a new identi¿cation algorithm for ARMAX systems based on a hereditary computation. The main bene¿t of this algorithm is that it is optimal in the transient. We have shown that the introduction of an extra virtual impulsive input could be used to extract the estimation of initial conditions. Moreover, it appears that no tuning parameter is needed when the input-output is not singular. When it is not the case, for binary signals for example, the algorithm uses pseudo-inverses since the dataset remain singular, which does not harm its ef¿ciency.

VIII. PROOF OF PROPOSITION 3.

Let us de¿ne the following sets: Backward residuals:

; @ 4===w> [ @ ] [ @ ] S ^] m] > ===> ] `> ;l @ 4===q 4 (11) Forward residuals: ; @ 4===w> ] @ ] ] @ ] S ^] m] > ===> ] `> ;l @ 4===q 4
With these de¿nitions, one clearly has:

i] > ===> ] j @ i[ j === i[ j (12) i] > ===> ] j @ i] j === i] j (13)
Using 13, the backward residuals can be rewritten as follows:

[ @ ] S ^] m] > ===> ] ` S ^] m] (14)
But, at previous step w 4, according to 11, one had:

[ @ ] S ^] m] > ===> ] @ ] S ^] m] > ===> ] àccording to property 2: [ @ ] S ^] m] > ===> ]
Therefore, formula 14 takes the following form:

[ @ [ S ^[ m] ` (15)
In the same way, using 12 in the expression of forward residuals, one has:

] @ ] S ^] m[ ` === S ^] m[ @ ] S ^] m[
The aim is now to construct the set of residuals iY > l @ 3===q 4j, where Y @ Q [ , such that H ^Y +Y , `@ L. Note that, using the assumption that Y @ 3> ; ? 4, this is equivalent to ensure that:

Assumption 16 leads to:

Q +Q , @ w w . 4 +L V +V , , (19)
Let Z @ +O , ] be the normalized form of the forward residuals. Expression 6 is a directly involves from 17 and 19. Note that, according this de¿nition, one has V @ H ^Y +Z , `. In the same way, the recursive equation of the forward residuals can be written:

Reporting 24 in 23 leads to the result.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Heat ow density through a two layer wall dataset

H

  ^Y +Y , `@ w . 4 w L (16) Let us de¿ne O / H ^] +] , `> ;l @ 3===q 4 and P / H ^] +Y , `> ;l @ 3===q 4. According to 15, one has: Y @ Q +Q , The autocorrelation matrix of Y can then be written as follows: H ^Y +Y , `@ Q +Q , L +P , +O , P +Q , +Q , where D stands for the inverse of transposed matrix D. Let us de¿ne: V / +P , +O ,

  to the normalization procedure, one has H ^Z +Z , `@ L. The computation of the experimental covariance of Z leads to:L @ +O , +O , L +V , V +O , +O ,(24) which leads to: +O , +O , @ L +V , V