M Darche

F Lopez-Caballero

D Aubry

B Tie

B Caicedo

F J Monroy Mendoza

Efficiency of a reduced model of seismic metamaterials

In the present work, we deal with both analytical/numerical and experimental modelings of seismic metamaterials, designed to protect buildings and civil infrastructures from incoming seismic waves. By using Floquet-Bloch transform, it is possible to strongly reduce the size of numerical model compared to a fully simulated model. Indeed, the numerical modeling comes down to the study of the eigenproblems defined on a primitive cell representative of the whole periodic medium. Existence of frequency bandgaps, one of the main issues of metamaterial, is analysed. Moreover, using reduced laboratory samples of soil with periodically distributed holes or concrete inclusions, actual efficiency of the studied seismic metamaterials is tested and compared with the analytical/numerical prediction. The comparison highlights the similarities observed between experimentally recorded transfer functions and numerically predicted bandgaps.

INTRODUCTION

The response of a building facing a seismic signal depends on its own geometric and material properties, but also on both the incident waves and the characteristics of soil in the surrounding area. As the standardization of old structures to the current seismic risk cartography could be tricky, a modification of the soil around the buildings could be a promising alternative to improve their response by rerouting harmful waves.

On the one hand, several studies tended to analyze the effect of barriers or other architectured structures places in the ground to limit the propagation of the waves induced by ground vibrations [START_REF] Avilès | Piles as Barriers for Elastic Waves[END_REF] or railway [START_REF] Van Hoorickx | Double wall barriers for the reduction of ground vibration transmission[END_REF]. Both are based on the previous work of [START_REF] Lysmer | Shear waves in plane infinite structures[END_REF]. [START_REF] Clouteau | Modifications of the ground motion in dense urban areas[END_REF] analyzed how a dense urban area can influence the propagation of seismic waves, especially in the case of periodically placed buildings.

On the other hand, the academic work of [START_REF] Veselago | The electrodynamics of substances with simultaneously negative values of ε and µ[END_REF] about electromagnetism, considering materials with both negative permeability and negative permittivity, have enlightened a new branch of the fundamental physic research about architectured materials. Materials architectured with specific periodic (or non-periodic) sub-wavelength microstructures, called metamaterials, were designed for their unconventional properties. Research advances have been achieved first by the electromagnetic community with for example the use of negative refraction index metamaterials [START_REF] Pendry | Magnetism from conductors and enhanced nonlinear phenomena[END_REF][START_REF] Pendry | Negative Refraction Makes a Perfect Lens[END_REF].

The study and development of metamaterials in the seismic domain to protect structures from incoming waves, especially the surface waves, are more recent. For example, [START_REF] Kim | Seismic waveguide of metamaterials[END_REF], [START_REF] Ungureanu | Auxetic-like metamaterials as novel earthquake protections[END_REF] used metamaterials to control seismic waves. Material and geometrical characteristics of the microstructure of a metamaterial play an essential role on its behaviour with respect to the wave propagation, parametric analyses are therefore necessary and several studies have been performed (e.g. [START_REF] Du | Elastic metamaterial-based seismic shield for both Lamb and surface waves[END_REF]).

Beside the theoretical and numerical investigations, several in-situ experiments have been designed to show the feasibility of seismic protections with, for example, periodically distributed inclusions [START_REF] Brûlé | Experiments on Seismic Metamaterials: Molding Surface Waves[END_REF] or even natural metamaterials like forests [START_REF] Colombi | Forests as a natural seismic metamaterial: Rayleigh wave bandgaps induced by local resonances[END_REF], and to test their efficiency. Contrary to in-situ test, the experimentation on reduced models allows a better control of the different material parameters. For example, Colombi et al.'s study (2016) on the propagation of surface waves through a forest was preceded by laboratory experiments on reduced scale [START_REF] Rupin | Experimental Demonstration of Ordered and Disordered Multiresonant Metamaterials for Lamb Waves[END_REF]. In laboratory models (reduced scale), in order to have an adequate replication of the in-situ stress field (full scale), several works have also been performed using centrifuge modeling. For example, isolation effects of geofoam barriers on wave propagation within soil were investigated using a centrifuge device by [START_REF] Murillo | Ground vibration isolation with geofoam barriers: Centrifuge modeling[END_REF].

In the present work, two types of seismic metamaterials with periodically distributed holes or inclusions are numerically and experimentally investigated. In particular, the existence of frequency bandgaps is analyzed. It is indeed one of the most important issues of using metamaterials to control seismic waves, as it means the non-propagation of waves within some segmented frequency ranges due to the interactions between waves and the periodic microstructure. Based on the previously developed theoretical/numerical studies of periodic media made of beams or plates [START_REF] Tie | Theoretical and numerical investigation of HF elastic wave propagation in two-dimensional periodic beam lattices[END_REF][START_REF] Tie | Theoretical and numerical modeling of membrane and bending elastic wave propagation in honeycomb thin layers and sandwiches[END_REF], new numerical models are built and theoretical frequency bandgaps are computed using Floquet-Bloch theory. Besides, reduced laboratory models of soil are designed and made. Then, their efficiency is tested and their coherence with numerical dispersion curves showing bandgaps is considered.

BLOCH WAVE BASED MODELING OF A PERIODIC MEDIUM

In order to study analytically the existence of frequency bandgaps due to the presence of inclusions or heterogeneities in a periodic media (Fig. 1), the Floquet-Bloch theory is used in this work. Thus, the following elastic wave equation is considered in the frequency domain, with u(x) the displacement field, C(x) the elasticity tensor, ρ(x) the mass density and ω the angular frequency:

div x (C(x) : ε(u(x))) = -ρ(x)ω 2 u(x) (1)
In Equation 1, ε(u) is the infinitesimal strain tensor defined as:

ε(u) = 1 2 ∇ x u + (∇ x u) T (2)
In the case of a periodic medium, we have C(x + γ n ) = C(x) and ρ(x + γ n ) = ρ(x), ∀n ∈ Z dim , with γ n = dim i=1 n i L i , dim the dimension of the physical space, and L i the vectors of periodicity (Fig. 1). Thanks to the periodicity of the periodic medium, the Floquet-Bloch transform allows giving a response on the entire medium with the modal analysis of the only primitive cell [START_REF] Floquet | Sur les équations différentielles linéaires à coefficients périodiques[END_REF][START_REF] Bloch | Quantum mechanics of electrons in crystal lattices[END_REF].

In fact, by the Floquet-Bloch transform a non-periodic field like u(x) defined on a periodic medium can be decomposed into periodic fields u B (x, k), which are called Bloch wave modes and are calculated as follows:

u B (x, k) = n∈Z dim u(x + γ n )e ik • (x+γ n ) (3)
where k is the Bloch wave vector belonging to the Brillouin zone [START_REF] Brillouin | Wave Propagation in Periodic Structures[END_REF]).

According to Equation 3, it is obvious that the principle of the Floquet-Bloch transform is to sum up, over the whole periodic medium, values of the signal taken at all the periodically corresponding points (Fig. 1). The obtained Bloch wave modes are therefore periodic and have the same periodicity than the periodic medium.

The Floquet-Bloch transform applied to the initial problem defined by Equation 1 gives rise to the following eigenvalue problems:

(div x + ik) • C : (ε(u B ) + ik ⊗ s u B ) = -ρω 2 u B (4)
Hence, to solve the wave propagation problem on the whole periodic medium, it is sufficient to study eigenvalues and eigenmodes on the primitive cell. Because the eigenproblem (Equation 4) is too complex to be solved analytically, there is a need to use numerical methods.

Otherwise, according to [START_REF] Brillouin | Wave Propagation in Periodic Structures[END_REF], the eigenanalysis of Bloch wave modes can be limited to the first Brillouin zone, which is in fact the elementary cell in the reciprocal space dual to the primitive cell in the physical space. However, the eigenanalysis is finally done within a smaller zone, called the irreducible Brillouin zone, which is defined by taking into account the Figure 1. Example of periodic medium and its primitive cell symmetries of the first Brillouin zone. In practice, a grid of discret points of k is used in the irreducible Brillouin zone, and for each Bloch wave vector k, the eigenproblem (Equation 4) is solved numerically using the finite element method.

To apply the finite element method to the eigenproblem (Equation 4), its weak form is developed as follows:

Ω C : (ε(u B ) + ik ⊗ s u B ) : (ε (v B ) + ik ⊗ s v B )dΩ = - Ω ρω 2 u • vdΩ (5)
with v virtual displacement field. To complete the definition of the above weak eigenproblem, appropriate periodic boundary conditions have to be applied on the primitive cell.

The primitive cell is discretized by a finite element mesh. In our study, a same finite element mesh is used whatever the value of k.

REDUCED METAMATERIAL MODELS

In order to verify the efficiency of a non-infinite metamaterial compared to the prediction obtained by the Floquet-Bloch theory in an infinite periodic medium, an experimental study is made. The attenuation on the amplitude of some frequencies induced by the presence of inclusions is measured in the frequency domain. For that, a comparison of the obtained response at several points of an homogeneous block of compacted soil and the one with periodically distributed heterogeneities is done. In the experimental model, the heterogeneities are simulated using empty holes and concrete inclusions. It allows to assess the effect of the properties contrast between the inclusions and the soil matrix on both the occurence and the location of frequency bandgaps. In addition, the efficiency of each solution will be evaluated.

Sample preparation and parameters

As previously indicated, by using a reduced laboratory model, the properties of soil are easier to control than in an in-situ test. In the present work, the tested block is made of compacted kaolin soil. The consistency of the soil is chosen to have the water content adjusted around its optimum Proctor (w opt = 26%) to reduce the damping. The Proctor optimum water content corresponds to the highest mass density that can be achieved. Otherwise, to assure the assumption of homogeneity of the metamaterial's matrix adopted in the present work, the hydrated kaolin is passed through a sieve of 4mm before the compaction.

To obtain a sample of size 300 × 300 × 250mm 3 (Fig. 2), the hydrated kaolin is compacted with three successive layers in its depth until a load of 1.1MPa using a hydraulic press. The obtained compact kaolin soil has its density ρ s equal to 1750kg/m 3 , and its elastic modulus E s is measured to be 233MPa, by calculating correlation between time signals recorded at different points (see Section 4) and with a Poisson ratio ν s taken equal to 0.25. From homogeneous samples, metamaterial ones are made in such a way that they are periodic in the horizontal plane on a depth of 8cm by drilling empty holes or by introducing inclusions. For the periodic layer, the primitive cell is designed as a 2D square pattern of dimension c × c = 13 × 13mm 2 with cylindrical empty holes or inclusions with a radius R equal to 4mm (Fig. 3).

Two types of metamaterial samples are made, with respectively empty holes or concrete inclusions as periodic heterogeneities. In the first case, holes are drilled using a broach with a diameter of 8mm. Because the holes are hand-made, there is a large uncertainty in their actual diameters and positions. For the second case, the holes drilled in the first case are filled by quick-setting cement, then tests are performed six days after pouring the concrete. Contrary to originally rigid inclusions such as metallic ones, cast-in-place concrete inclusions have the advantage of fitting the shape of the holes without the use of force, hence, a "prefect" adhesion on matrix/inclusion interfaces can be obtained without local modification of the properties of the compact kaolin soil. After compression tests on samples of concrete poured in the same way, the properties of the concrete inclusions are identified as: E c = 10GPa, ρ c = 2100kg/m 3 and ν c = 0.3.

Boundary condition setting and experimental devices

In order to reduce the influence of reflected waves, the sample is laid on a rubber layer. All its other boundaries are free surfaces.

The incident wave signal is generated with a piezoelectric device at the top horizontal boundary of tested samples (Fig. 2). Especially in our case surface waves are generated, as under an electrical signal, a piezoelectric crystal convert the energy into a stress field and it is therefore possible to control the type of introduced mechanical waves. For the sake of simplicity of experimental For each frequency, the propagation of the signal is measured on the top horizontal surface by means of seven accelerometers (Fig. 2, Fig. 4). The accelerometers have an operational range until 20kHz. The sampling frequency is 100kHz on samples of 2s. As a first analysis of experimental results, only the signals recorded by the accelerometers 1, 3 and 7 (Fig. 4) are considered in this work.

Experimental configurations

Regarding the position of the piezoelectric device, two experimental configurations are tested: one with the source out of the metamaterial barrier and an other one with the source inside (Fig. 4). With the first configuration, it is expected to test the influence of metamaterial barrier on the propagation of the incoming wave, while with the second one, it is expected to investigate the behaviour of the waves trapped by the barrier within it. The same block of compacted soil is used for homogeneous sample and for the improved one. The homogeneous model will be used as a reference case. As the attenuation within metamaterial samples is measured with respect to the reference one, uncertainties on material properties and the lack of knowledge on boundary conditions could be considered as suppressed.

4 NUMERICAL/EXPERIMENTAL COMPARISON 2D numerical models based on the Floquet-Bloch theory presented here above (Section 2) are used to perform numerical/experimental comparison. They represent in fact plane wave propagation in an infinite 2D periodic domain with a primitive cell defined in Figure 3. It is worth noticing that, for the tested samples, neither the vertical direction, in which materials are homogeneous, nor its boundaries are not taken into account by the numerical model. Caution is therefore required in interpreting the comparison presented hereafter.

To set appropriate material properties in the numerical models, especially those of the compact soil, the classical method of cross-correlation [START_REF] Elgamal | Seismic downhole arrays and applications in practice[END_REF]) is used to estimate the Rayleigh wave velocity propagating in it. Cross-correlation function c aiaj between two temporal signals recorded by accelerometers a i (t) and a j (t) is expressed as:

c aiaj (τ ) = 1 N -m N -m n=1 a i (n∆t)a j ((n + m)∆t) 1 N N n=1 a 2 i (n∆t) 1 N N n=1 a 2 j (n∆t) m = 0, 1, 2, . . . , N (6)
with τ the time delay, N the length of the signal and ∆t the time step of the signal. ). The frequency of the signal of 2s is 2MHz; A zoom of the temporal between 0.08s and 0.09s is presented.

For example, the cross-correlation between signals recorded at point 1 and 3 for the frequency equal to 4.25kHz results in a velocity the Rayleigh waves of V R = 206m/s. Then the Young's modulus of the compacted soil is determined using the following equation:

E s = 2ρ s (1 + ν s ) • V 2 S (7)
with V S the shear wave velocity linked to V R by:

V R = 0.862 + 1.14ν s 1 + ν s V S (8)
With the numerical models, eigenvalues for each discret k chosen in the irreducible Brillouin zone are computed by performing a finite element analysis of the primitive cell. Dispersion curves are then plotted and existence of frequency bandgaps can be highlighted.

For the data post-processing of experimental results, the transfer functions α(f ) with respect to the frequency f is defined as the ratio of the signal recorded after the metamaterial barriers to the incident one:

α M (f ) = a 7 (f ) a 1 (f ) ; α (f ) = a 1 (f ) a 7 (f ) (9)
with a i (f ) the Fourier transform of the vertical accelerations measured with the sensor i. In Equation 9, α M gives an estimate of attenuation for the first experimental configuration (see Figure 4(left)) and α C for the second one (see Figure 4(right)). Finally, by calculating the ratio of α(f ) obtained within a metamaterial sample to the one obtained within the homogeneous simple, the attenuation due to the metamaterial barriers can be actually measured without the one caused by the geometrical spreading of wavefront (visible on Fig. 5).

The Figure 6 presents the comparison between numerical predictions and experimental measures for both experimental configurations and for both cases of periodically distributed holes (top line in Fig. 6) or concrete inclusions (bottom line in Fig. 6).

• Figure (a) presents dispersion curves predicted theoretically by solving the eigenproblem (Equation 4) and computed using finite elements. The first Floquet-Bloch modes are shown within the frequency range [0, 20]kHz. In the case of the primitive cell with a hole, the first bandgap is obtained between 9.4kHz and 10.2kHz. In the case of the primitive cell with a concrete inclusion, the first frequency bandgap is obtained between 16.5kHz and 18.6kHz and is much larger. • Figures (b)-(c) present the attenuation levels obtained experimentally in terms of the ratios between transfer functions. We note that, as the valid operating range of the accelerometers are inferior to 20kHz, experimental results at this frequency are questionable. For the metamaterial sample with holes, the experimental results appear coherent with the numerical predictions, as an important attenuation level is observed around 8kHz, close to the frequency range of the first numerically predicted bandgaps. For the metamaterial sample with concrete inclusions, the numerical/experimental comparison is less obvious. However, for the first experimental configuration, a large frequency range of attenuation is observed between 10kHz and 15kHz, which is in agreement with the numerical predictions.

The comparison between numerical and experimental results is encouraging but not conclusive. Several possible causes can explain differences. First of all, the 2D numerical model is not representative of the 3D experiment, because the third vertical direction is ignored and the boundary conditions are not taken into account in the Floquet-Bloch theory. Moreover, there are uncertainties in the identified material properties, such as the value of Rayleigh wave velocity.

CONCLUSIONS

In this work, 2D numerical models based on Floquet-Bloch theory have been developed to study a periodic metamaterial by eigenanalysis of a primitive cell. It was used to locate frequency bandgaps by plotting dispersion curves.

In order to test actual efficiency of metamaterial, an experimental test has been made on a block of compacted soil with periodically distributed empty holes or concrete inclusions. Experimental results showed the existence of frequency ranges of attenuation. However, they did not fit well with the numerically predicted frequency bandgaps, which brought out the limits of the developed numerical model.

Figure 2 .

 2 Figure 2. Presentation of the experimental setting: A homogeneous sample with the electronic devices

Figure 4 .

 4 Figure 4. Top view of two different configurations for the positions of accelerometers (blue circles) and of the piezoelectric device (red)

Figure 5 .

 5 Figure 5. Temporal signals recorded at accelerometers 1 and 3 in the homogeneous sample tested within the first configuration (see Fig. 4 (left)). The frequency of the signal of 2s is 2MHz; A zoom of the temporal between 0.08s and 0.09s is presented.

Figure 6 .

 6 Figure 6. Comparison between the theoretically predicted frequency bandgaps (grey) and experimental transfer functions in the two cases of empty holes (first line) and concrete inclusions (second line). (a) Dispersion curves showing the first Floquet-Bloch modes and the first bandgap in the frequency range [0, 20]kHz. (b)-(c) Transfer functions for the two different configurations (Fig. 4).


Hugo Lhachemi
email: hugo.lhachemi@ucd.ie

Christophe Prieur
email: christophe.prieur@gipsa-lab.fr

Emmanuel Trélat
email: emmanuel.trelat@sorbonne-universite.fr

Neumann trace tracking of a constant reference input for 1-D boundary controlled heat-like equations with delay

Keywords: 1-D reaction-diffusion equation, PI regulation control, Neumann trace, Delay boundary control, Partial Differential Equations (PDEs)

come

Neumann trace tracking of a constant reference input for 1-D boundary controlled heat-like equations with delay

INTRODUCTION

State of the art

Proportional-Integral (PI) regulation control of infinitedimensional systems is an active topic of research. While early works in this area were reported in the 80's for bounded control operators [START_REF] Pohjolainen | Robust multivariable PI-controller for infinite dimensional systems[END_REF][START_REF] Pohjolainen | Robust controller for systems with exponentially stable strongly continuous semigroups[END_REF][START_REF] Xu | A robust PI-controller for infinite-dimensional systems[END_REF], the PI boundary regulation of infinitedimensional systems is more recent. This includes linear hyperbolic systems [START_REF] Santos | Boundary control with integral action for hyperbolic systems of conservation laws: Stability and experiments[END_REF][START_REF] Xu | Multivariable boundary PI control and regulation of a fluid flow system[END_REF][START_REF] Bastin | Stability of linear density-flow hyperbolic systems under PI boundary control[END_REF][START_REF] Lamare | Control of 2×2 linear hyperbolic systems: Backstepping-based trajectory generation and PI-based tracking[END_REF], 1-D nonlinear transport equation [START_REF] Coron | PI controllers for 1-D nonlinear transport equation[END_REF][START_REF] Trinh | Design of integral controllers for nonlinear systems governed by scalar hyperbolic partial differential equations[END_REF], regulation of the downside angular velocity of a drilling string [START_REF] Terrand-Jeanne | Regulation of the downside angular velocity of a drilling string with a PI controller[END_REF], and the regulation of a drilling pipe under friction [START_REF] Barreau | Practical stabilization of a drilling pipe under friction with a PI-controller[END_REF]. The PI regulation of open-loop exponentially stable semigroups with unbounded control operators were reported in (Terrand-Jeanne et al., 2018a, 2019) via a Lyapunov functionalbased design procedure. This paper is focused on the PI regulation control of the left Neumann trace of a one-dimensional reaction-diffusion equation, which might be either open-loop stable or unstable, with a delayed right Dirichlet boundary control. Specifically, we aim at achieving the Neumann trace track-This publication has emanated from research supported in part by a research grant from Science Foundation Ireland (SFI) under grant number 16/RC/3872 and is co-funded under the European Regional Development Fund and by I-Form industry partners.

ing of a constant reference input in spite of the presence of an arbitrarily large constant input delay and a stationary distributed disturbance. It was shown in [START_REF] Guzmán | Stabilization of the linear Kuramoto-Sivashinsky equation with a delayed boundary control[END_REF] that backstepping-based control design can be used to achieve the feedback stabilization of a reaction-diffusion equation in the presence of an arbitrarily large input delay. In this paper, we adopt the approach reported in [START_REF] Prieur | Feedback stabilization of a 1D linear reaction-diffusion equation with delay boundary control[END_REF] which takes advantage of the following design procedure initially reported in [START_REF] Russell | Controllability and stabilizability theory for linear partial differential equations: recent progress and open questions[END_REF]: 1) design of the controller on a finite-dimensional model capturing the unstable modes of the original infinitedimensional system; 2) use of a suitable Lyapunov functional to guarantee the stability of the resulting closedloop infinite-dimensional system. This control design procedure, which was used in (Coron andTrélat, 2004, 2006;[START_REF] Schmidt | Controllability of Couette flows[END_REF] to stabilize semilinear heat, wave or fluid equations via (undelayed) boundary feedback control, was extended in [START_REF] Prieur | Feedback stabilization of a 1D linear reaction-diffusion equation with delay boundary control[END_REF] to the case of delay boundary control of a one-dimensional reaction-diffusion equation in which the contribution of the input-delay was managed by the synthesis of a predictor feedback via the classical Artstein transformation [START_REF] Artstein | Linear systems with delayed controls: a reduction[END_REF][START_REF] Richard | Time-delay systems: an overview of some recent advances and open problems[END_REF][START_REF] Bresch-Pietri | New formulation of predictors for finite-dimensional linear control systems with input delay[END_REF]. This control strategy was first reused in [START_REF] Guzmán | Stabilization of the linear Kuramoto-Sivashinsky equation with a delayed boundary control[END_REF] for the delay boundary feedback stabilization of a linear Kuramoto-Sivashinsky equation and then generalized to the delay boundary feedback stabilization of a class of diagonal infinite-dimensional systems for either a constant (Lhachemi and Prieur, 2020;[START_REF] Lhachemi | Control law realification for the feedback stabilization of a class of diagonal infinite-dimensional systems with delay boundary control[END_REF] or a time-varying (Lhachemi et al., 2019a(Lhachemi et al., , 2020) ) input delay.

Investigated control problem

Let L > 0, let c ∈ L ∞ (0, L), and let D > 0 be arbitrary. We consider the one-dimensional reaction-diffusion equation over (0, L) with delayed Dirichlet boundary control:

y t = y xx + c(x)y + d(x), (t, x) ∈ R * + × (0, L) (1a) y(t, 0) = 0, t ≥ 0 (1b) y(t, L) = u D (t) u(t -D), t ≥ 0 (1c) y(0, x) = y 0 (x),
x ∈ (0, L) (1d) where y(t, •) ∈ L 2 (0, L) is the state at time t, u(t) ∈ R is the control input, D > 0 is the (constant) control input delay, d ∈ L 2 (0, L) is a stationary distributed disturbance, and y 0 ∈ H 2 (0, L) with y 0 (0) = 0 and y 0 (L) = u(-D) is the initial condition.

Our objective is to achieve the PI regulation control of the left Neumann trace y x (t, 0) to some prescribed constant reference input r ∈ R in spite of the stationary distributed disturbance d, i.e., y x (t, 0) → r as t → +∞. Note that an exponentially stabilizing controller for (1a-1d) was reported in [START_REF] Prieur | Feedback stabilization of a 1D linear reaction-diffusion equation with delay boundary control[END_REF] in the disturbancefree case (d = 0) for a system trajectory evaluated in H1 0norm. The control strategy that we develop in the present paper elaborates on the one of [START_REF] Prieur | Feedback stabilization of a 1D linear reaction-diffusion equation with delay boundary control[END_REF], adequately combined with a PI procedure.

CONTROL DESIGN STRATEGY

The sets of nonnegative integers, positive integers, real, nonnegative real, and positive real are denoted by N, N * , R, R + , and R * + , respectively. All the finite-dimensional spaces R p are endowed with the usual Euclidean inner product x, y = x y and the associated 2-norm x =

x, x = √ x x. For any matrix M ∈ R p×q , M stands for the induced norm of M associated with the above 2norms. For a given symmetric matrix P ∈ R p×p , λ m (P ) and λ M (P ) denote its smallest and largest eigenvalues, respectively. In the sequel, the time derivative ∂f /∂t is either denoted by f t or ḟ while the spatial derivative ∂f /∂x is either denoted by f x or f .

Augmented system for PI feedback control

The control design objective is: 1) to stabilize the reactiondiffusion system (1a-1d); 2) to ensure the tracking of the constant reference input r ∈ R by the left Neumann trace y x (t, 0). We address this problem by designing a PI controller. Following the general PI scheme, we introduce a new state z(t) ∈ R taking the form of the integral of the tracking error y x (t, 0) -r:

y t = y xx + c(x)y + d(x), (t, x) ∈ R * + × (0, L) (2a) ż(t) = y x (t, 0) -r, t ≥ 0 (2b) y(t, 0) = 0, t ≥ 0 (2c) y(t, L) = u D (t) u(t -D), t ≥ 0 (2d) y(0, x) = y 0 (x), x ∈ (0, L) (2e) z(0) = z 0 (2f) with z 0 ∈ R the initial condition of the integral component.
As we are only concerned in prescribing the future of the system, we assume that the system is uncontrolled for t < 0, i.e. u(t) = 0 for t < 0. Thus, we assume in the rest of the paper that y 0 ∈ H 2 (0, L) ∩ H 1 0 (0, L).

Spectral reduction

We rewrite (2) as an equivalent homogeneous Dirichlet problem. Assuming 1 that u is continuously differentiable and setting w(t, x) = y(t, x) -x L u D (t), we have

w t = w xx + c(x)w + x L c(x)u D - x L uD (t) + d(x) (3a) ż(t) = w x (t, 0) + 1 L u D (t) -r (3b) w(t, 0) = w(t, L) = 0 (3c) w(0, x) = y 0 (x) - x L u D (0) = y 0 (x) (3d)
z(0) = z 0 (3e) for t > 0 and x ∈ (0, 1). Introducing the real state-space L 2 (0, 1) endowed with its usual inner product f, g

= L 0 f (x)g(x) dx and the operator A = ∂ xx + c id : D(A) ⊂ L 2 (0, L) → L 2 (0, L) defined on the domain D(A) = H 2 (0, L) ∩ H 1 0 (0, L), (3a-3c) can be rewritten as w t (t, •) = Aw(t, •) + a(•)u D (t) + b(•) uD (t) + d(•) (4a) ż(t) = w x (t, 0) + 1 L u D (t) -r (4b)
with a(x) = x L c(x) and b(x) = -x L for every x ∈ (0, L), with initial conditions (3d-3e). Since A is self-adjoint and has compact resolvent, we consider a Hilbert basis (e j ) j≥1 of L 2 (0, L) consisting of eigenfunctions of A associated with the sequence of simple real eigenvalues

-∞ < • • • < λ j < • • • < λ 1 with λ j -→ j→+∞ -∞.
Note that e j (•) ∈ H 1 0 (0, L) ∩ C 1 ([0, L]) for every j ≥ 1 and

e j (0) ∼ 2 L |λ j |, λ j ∼ - π 2 j 2 L 2 , (5)
when j → +∞. Since w(0, •) = y 0 ∈ H 2 (0, L) ∩ H 1 0 (0, L), the classical solution w(t, •) ∈ H 2 (0, L) ∩ H 1 0 (0, L) of (4a) can be expanded as a series in the eigenfunctions e j (•), convergent in H 1 0 (0, L),

w(t, •) = +∞ j=1 w j (t)e j (•). (6)
Thus ( 4) is equivalent to the infinite-dimensional control system:

ẇj (t) = λ j w j (t) + a j u D (t) + b j uD (t) + d j (7a) ż(t) = j≥1 w j (t)e j (0) + 1 L u D (t) -r (7b)
for j ∈ N * , with w j (t) = w(t, •), e j , a j = a, e j , b j = b, e j , and d j = d, e j . Introducing the auxiliary control input v u, and denoting v D (t) v(t -D), (7) can be rewritten as uD

(t) = v D (t) (8a) ẇj (t) = λ j w j (t) + a j u D (t) + b j v D (t) + d j (8b) ż(t) = j≥1 w j (t)e j (0) + 1 L u D (t) -r (8c)
for j ∈ N * . As u(t) = 0 for t < 0, (8a) yields v(t) = 0 for t < 0 and the initial condition u D (0) = 0.

Finite-dimensional truncated model

In what follows, we fix the integer n ∈ N * such that λ n+1 < 0 ≤ λ n . In particular, we have λ j ≥ 0 when 1 ≤ j ≤ n and λ j ≤ λ n+1 < 0 when j ≥ n + 1. Then, introducing as in [START_REF] Prieur | Feedback stabilization of a 1D linear reaction-diffusion equation with delay boundary control[END_REF]: 8a) and the n first equations of (8b) yield the truncated model:

X 1 (t) =     u D (t) w 1 (t) . . . w n (t)     , A 1 =     0 0 • • • 0 a 1 λ 1 • • • 0 . . . . . . . . . . . . a n 0 • • • λ n     , B 1 = (1 b 1 . . . b n ) , D 1 = (0 d 1 . . . d n ) , with X 1 (t) ∈ R n+1 , A 1 ∈ R (n+1)×(n+1) , B 1 ∈ R n+1 , D 1 ∈ R n+1 , (
Ẋ1 (t) = A 1 X 1 (t) + B 1 v D (t) + D 1 . (9)
Now, as

e j (0) λj 2 ∼ 2L π 2 j 2 when j → +∞, we can introduce the change of variable ζ(t) z(t) - j≥n+1 e j (0) λ j w j (t), (10)
whose time derivative is given by

ζ(t) = ż(t) - j≥n+1 e j (0) λ j ẇj (t) = αu D (t) + βv D (t) -γ + n j=1
w j (t)e j (0), where we have used (8b-8c), with

α = 1 L - j≥n+1 e j (0) λ j a j , β = - j≥n+1 e j (0) λ j b j , (11a)
γ = r + j≥n+1 e j (0) λ j d j . (11b)
Then we have

ζ(t) = L 1 X 1 (t) + βv D (t) -γ (12)
with L 1 = α e 1 (0) . . . e n (0) ∈ R 1×(n+1) . Now, defining the augmented state-vector X

(t) = X 1 (t) ζ(t) ∈ R n+2 , the exogenous input Γ = D 1 -γ ∈ R n+2 and the matrices A = A 1 0 L 1 0 ∈ R (n+2)×(n+2) , B = B 1 β ∈ R n+2 ,
(13) we obtain from ( 9) and (12) the control system Ẋ(t) = AX(t) + Bv D (t) + Γ (14) which is the finite-dimensional truncated model capturing the unstable part of the infinite-dimensional system along with the introduced integral component for PI regulation. Putting together the truncated model ( 14) along with (8b) for j ≥ n + 1, we get the final representation used for both control design and stability analyses:

Ẋ(t) = AX(t) + Bv D (t) + Γ (15a) ẇj (t) = λ j w j (t) + a j u D (t) + b j v D (t) + d j (15b) with j ≥ n + 1.

Control design strategy

The adopted control design strategy relies on the use of the classical predictor feedback to stabilize the finitedimensional truncated model (15a). Specifically, introducing the Artstein transformation [START_REF] Artstein | Linear systems with delayed controls: a reduction[END_REF])

Z(t) = X(t) + t t-D e A(t-D-τ ) Bv(τ ) dτ, (16)
we have Ż(t) = AZ(t) + e -DA Bv(t) + Γ.

Since (A, B) satisfies the Kalman condition (the proof of this claim is omitted due to space limitation; details can be found in [START_REF] Lhachemi | PI regulation of a reaction-diffusion equation with delayed boundary control[END_REF]), the pair (A, e -DA B) also satisfies the Kalman condition and we infer the existence of a feedback gain KZ(t) where χ [0,+∞) denotes the characteristic function of the interval [0, +∞), we obtain the stable closed-loop dynamics Ż(t) = A K Z(t) + Γ.

K ∈ R 1×(n+2) such that A K A + e -DA BK is Hurwitz. Setting v(t) = χ [0,+∞) (t)
Remark 1. In original coordinates, the control input v is solution of the fixed point implicit equation

v(t) = χ [0,+∞) (t)KX(t)+K t max(t-D,0) e A(t-D-τ ) Bv(τ ) dτ.
Existence and uniqueness of the solution of the above equation are reported in [START_REF] Bresch-Pietri | New formulation of predictors for finite-dimensional linear control systems with input delay[END_REF]. •

Characterization of the equilibrium point

We characterize the equilibrium point of the closed-loop system associated with a constant reference input r ∈ R and a stationary distributed disturbance d ∈ L 2 (0, L). In the sequel, we denote by a subscript "e" the equilibrium value of the different quantities. Noting that u D,e = u e and v D,e = v e , we have 0 = A K Z e + Γ 0 = λ j w j,e + a j u e + b j v e + d j , j ≥ n + 1

In particular, from v e = KZ e , we have AZ e + e -DA Bv e + Γ = 0. As the first row of A and Γ are null and [1 0 . . . 0] e -DA B = 1, we obtain that v e = 0. Moreover we have Z e = -A -1 K Γ. Then we can set X e = Z e because AX e + Bv D,e + Γ = A K Z e + Γ = 0, which is compatible with the Artstein transformation since v e = 0 implies Z e = X e + t t-D e (t-s-D)A Bv e ds. Then we have

u e = u D,e = [1 0 . . . 0] Z e = -[1 0 . . . 0] A -1
K Γ, and w j,e = -a j λ j u e -d j λ j for j ≥ n + 1. Now, noting that both (w j,e ) j and (λ j w j,e ) j are square-summable sequences, we can define w e j≥1 w j,e e j ∈ D(A) = H 2 (0, L) ∩ H 1 0 (0, L).

Expending the last row of AX e + Γ = 0, we obtain that w e (0) + 1 L u e = r. Consequently, the introduction of y e w e + x L u e yields y e (0) = r. Finally, from Aw e = j≥1 λ j w j,e e j we deduce that Aw e + au D,e + bv D,e + d = 0.

Dynamics of deviations

We now define the deviations of the various quantities with respect to their equilibrium value: ∆X = X -X e , ∆Z = Z -Z e , ∆w = w -w e , ∆w j = w j -w j,e , ∆ζ = ζζ e , ∆u = u -u e , ∆u D = u D -u e , ∆v = v -v e , and ∆v D = v D -v D,e . Then, in original coordinates:

∆w t = A∆w + a∆u D + b∆v D (17) and ∆ Ẋ(t) = A∆X(t) + B∆v D (t)
∆ ẇj (t) = λ j ∆w j (t) + a j ∆u D (t) + b j ∆v D (t) for j ≥ n + 1 with the auxiliary control input ∆v(t) = χ [0,+∞) (t)K∆Z(t) (because v e = KZ e = 0) where

∆Z(t) = ∆X(t) + t t-D e (t-s-D)A B∆v(s) ds. ( 18
)
In Z coordinates, the closed-loop dynamics is given by ∆

Ż(t) = A K ∆Z(t) (19a) ∆ ẇj (t) = λ j ∆w j (t) + a j ∆u D (t) + b j ∆v D (t) (19b) for j ≥ n + 1.

STABILITY ANALYSIS

The stability of the closed-loop infinite-dimensional system is assessed by the following theorem. Theorem 1. There exist κ, C 1 > 0 such that ∆u

D (t) 2 + ∆ζ(t) 2 + ∆w(t) 2 H 1 0 (0,L) (20)
≤ C 1 e -2κt ∆u D (0) 2 + ∆ζ(0) 2 + ∆w(0) 2 H 1 0 (0,L)
. The proof of Theorem 1 relies on the following Lyapunov function:

V (t) = M 2 ∆Z(t) P ∆Z(t) (21) + M 2 t max(t-D,0) ∆Z(s) P ∆Z(s) ds - 1 2 j≥1 λ j ∆w j (t) 2 ,
where P ∈ R (n+2)×(n+2) is the solution of the Lyapunov equation

A K P + P A K = -I and M > 0 is chosen such that M > max γ 1 λ 1 λ m (P ) , 2 γ 1 a 2 + b 2 e -DA K 2 K 2
with γ 1 2 max 1, De 2D A BK 2 . Useful properties of V are stated in the three following lemmas. Due to space limitation, only a sketch of proof is provided. Lemma 1. There exists a constant C 1 > 0 such that

V (t) ≥ C 1 j≥1 (1 + |λ j |)∆w j (t) 2 , (22a)
V (t) ≥ C 1 ∆u D (t) 2 + ∆ζ(t) 2 + ∆w(t) 2 H 1 0 (0,L) , (22b) V (t) ≥ C 1 ∆Z(t) 2 , (22c) for every t ≥ 0.
Sketch of proof Using M > γ1λ1 λm(P ) , the claimed estimates are obtained similarly to the ones reported in [START_REF] Prieur | Feedback stabilization of a 1D linear reaction-diffusion equation with delay boundary control[END_REF]. 2 Lemma 2. There exist κ, C 2 > 0 such that V (t) ≤ e -2κ(t-D) V (D) for every t ≥ D.

Sketch of proof

As A is self-adjoint, we have for t > D,

V (t) = - M 2 ∆Z(t) 2 - M 2 t t-D ∆Z(s) 2 ds -A∆w(t) 2 -A∆w(t), a ∆u D (t) -A∆w(t), b ∆v D (t). The use of Cauchy-Schwarz and Young inequalities show that V (t) ≤ - 1 2 A∆w(t) 2 - γ 2 λ M (P ) ∆Z(t) P ∆Z(t) + t t-D ∆Z(s) P ∆Z(s) ds for all t > D and γ 2 = M/2-γ 1 a 2 + b 2 e -DA K 2 K 2 > 0.
Similarly to [START_REF] Prieur | Feedback stabilization of a 1D linear reaction-diffusion equation with delay boundary control[END_REF], we infer the existence of a constant γ 3 > 0 such that, for all t ≥ 0, -A∆w(t), ∆w(t) ≤ γ 3 A∆w(t) 2 . Consequently, we obtain that V (t) ≤ -2κV (t) for all

t > D with κ = 1 2 min 2γ 2 M λ M (P ) , 1 γ 3 > 0. 2 Lemma 3. There exists C 2 > 0 such that V (t) ≤ C 2 ∆u D (0) 2 + ∆ζ(0) 2 + ∆w(0) 2 H 1 0 (0,L) for all 0 ≤ t ≤ D with ∆u D (0) = -u e .
Sketch of proof Estimations similar to the ones reported in the proof of Lemma 2 show the existence of γ 4 > 0 such that V (t) ≤ γ 4 ∆X(0) 2 for all 0 ≤ t < D. Therefore, V (t) ≤ V (0) + Dγ 4 ∆X(0) 2 for all 0 ≤ t ≤ D. The estimation of V (0) from ( 21) and a direct integration with t ≤ D show the claimed result. 2

The proof of Theorem 1 is now a straightforward combination of the results reported in Lemmas 1, 2 and 3. Recalling that ∆v(t) = K∆Z(t) for t ≥ 0 and ∆v(t) = 0 for t < 0, we also obtain that ∆v D (t) 2 ≤ Ĉ1 e -2κt ∆u D (0) 2 + ∆ζ(0) 2 + ∆w(0) 2

H 1 0 (0,L) (23) for t ≥ 0 with Ĉ1 = K 2 C 1 e 2κD .

SETPOINT REFERENCE TRACKING ANALYSIS

We assess that the tracking of the constant reference input r is achieved in spite of the stationary distributed disturbance d. Theorem 2. Let κ > 0 be provided by Theorem 1. There exists C 2 > 0 such that |y x (t, 0) -r| (24)

≤ C 2 e -κt |∆u D (0)| + |∆ζ(0)| + ∆w(0) H 1 0 (0,L) + A∆w(0) L 2 (0,L) .
Sketch of proof Recalling that w e,x (0) + 1 L u e = r, we have

|y x (t, 0) -r| = w x (t, 0) + 1 L u D (t) -r ≤ |w x (t, 0) -w e,x (0)| + 1 L |∆u D (t)|. (25)
From the exponential convergence of ∆u D (t) to zero provided by (20), it is sufficient to study the term w x (t, 0)w e,x (0) = j≥1 ∆w j (t)e j (0). As e j (0) ∼ 2 L |λ j |, there exists a constant γ 7 > 0 such that |e j (0)| ≤ γ 7 |λ j | for all j ≥ n + 1. Let m ≥ n + 1 be such that η -λ m > κ > 0. Thus λ j ≤ -η < -κ < 0 for all j ≥ m. We have:

|w x (t, 0) -w e,x (0)| ≤ m-1 j=1 |∆w j (t)||e j (0)| + γ 7 j≥m |λ j ||∆w j (t)| ≤ m-1 j=1 e j (0) 2 m-1 j=1 ∆w j (t) 2 (26) + γ 7 j≥m 1 |λ j | j≥m λ 2 j ∆w j (t) 2
where j≥m

1 |λj | < +∞ because |λ j | ∼ π 2 j 2 /L 2 .
Based on (20) and Poincaré inequality, it is sufficient to study the term j≥m λ 2 j ∆w j (t) 2 . To do so, we integrate (19b) for j ≥ m and we use estimates ( 20) and ( 23

λ 2 j ∆w j (t) 2 ≤ 2e -2κt A∆w(0) 2 + 2C 2 3 η 2 (η -κ) 2 e -2κt ∆CI 2 . (27) with C 3 > 0 defined by C 2 3 = j≥m C 2 3,j ≤ 2C 1 a 2 + 2 Ĉ1 b 2 .
Using now (25) along with (26) and estimates ( 20) and ( 27), we obtain the existence of a constant C 2 > 0 such that the claimed estimate (24) holds for all t ≥ 0. 2

NUMERICAL ILLUSTRATION

We take c = 1.25, L = 2π, and D = 1 s. The three first eigenvalues of the open-loop system are λ 1 = 1, λ 2 = 0.25, and λ 2 = -1. Only the two first modes need to be stabilized. Thus we have n = 2 and we compute the feedback gain K ∈ R 1×4 such that the poles of the closed-loop truncated model (capturing the two unstable modes of the infinite-dimensional system plus two integral components, one for the control input and one for reference tracking) are given by -0.5, -0.6, -0.7, and -0.8. The adopted numerical scheme is the modal approximation of the infinite-dimensional system using its first 10 modes. initial condition is set as y 0 (x) = -x L 1 -x L . The obtained simulation results with r = 50 and d(x) = x are depicted in Fig. 1. As expected from the theoretical analysis, the PI controller achieves the stabilization of the reaction-diffusion equation and ensures that the Neumann trace y x (t, 0) tracks the constant reference input r.

CONCLUSION

This paper discussed the PI regulation control of the left Neumann trace of a one-dimensional linear reactiondiffusion equation with delayed right Dirichlet boundary control. The proposed strategy extends to PI control a recently proposed approach for the delay boundary feedback stabilization of infinite-dimensional systems combining spectral reduction and the use of the classical Artstein transformation for handling the delay in the control input. The validity of this control strategy for the tracking of a constant reference input and in the presence of a stationary perturbation was assessed via a Lyapunov-based argument. The extension of these results to the set-point regulation control of a time-varying reference input r(t) and in the presence of a time-varying distributed perturbation d(t, x) can be found in [START_REF] Lhachemi | PI regulation of a reaction-diffusion equation with delayed boundary control[END_REF].

  ), yielding |λ j ∆w j (t)| ≤ e λj t |λ j ∆w j (0j )e λj (t-τ ) {|a j ||∆u D (τ )| + |b j ||∆v D (τ )|} dτ ≤ e -ηt |λ j ∆w j (0)| + C 3,j e λj t t 0 (-λ j )e -λj τ e -κτ dτ ∆CI with ∆CI = ∆u D (0) 2 + ∆ζ(0) 2 + ∆w(0) 2 H 1 0 (0,L) and constant C 3,j = |a j | C 1 + |b j | Ĉ1 . As λ j ≤ -η < -κ for all j ≥ m, we obtain that e λj t t 0 (-λ j )e -λj τ e -κτ dτ = λj λj +κ (e -κt -e λj τ ) ≤ λj λj +κ e -κt ≤ η η -κ e -κt , hence |λ j ∆w j (t)| ≤ e -κt |λ j ∆w j (0)| + C 3,j η η -κ e -κt ∆CI and thus j≥m

  Fig. 1. Time evolution of the closed-loop system

This property will be ensured by the construction carried out in the sequel.


Makia Zmitri

Hassen Fourati

Christophe Prieur

Inertial Velocity Estimation for Indoor Navigation Through Magnetic Gradient-based EKF and LSTM Learning Model

This paper presents a novel method to improve the inertial velocity estimation of a mobile body, for indoor navigation, using solely raw data from a triad of inertial sensors (accelerometer and gyroscope), as well as a determined arrangement of magnetometers array. The key idea of the method is the use of deep neural networks to dynamically tune the measurement covariance matrix of an Extended Kalman Filter (EKF). To do so, a Long Short-Term Memory (LSTM) model is derived to determine a pseudo-measurement of inertial velocity of the target under investigation. This measurement is used afterwords to dynamically adapt the measurement noise parameters of a magnetic field gradient-based EKF. As it was shown in the literature, there is a strong relation between inertial velocity and magnetic field gradient, which is highlighted with the proposed approach in this paper. Its performance is tested on the Openshoe dataset, and the obtained results compete with the INS/ZUPT approach, that unlike the proposed solution, can only be applied on foot-mounted applications and is not adequate to all walking paces.

I. INTRODUCTION

The recent growth of interest in smart systems technology, autonomous robotics and location-aware services generates a strong motivation to develop robust and efficient techniques for tracking the position of a moving body, in both outdoor and indoor environments. Indoor navigation is difficult to deal with because of the limited indoor coverage of Global Positioning System (GPS). Nowadays, different solutions have been proposed to deal with this issue, most of them use a variety of sensors that usually require costly and pre-installed infrastructures, which can be inaccessible in various applications that introduce hostile environments or emergency situations. For this reason, navigation solutions that are independent from pre-equipped buildings are highly preferred. The most used technique in these situations is the one fusing data from low-cost Inertial Measurement Units (IMUs), composed of accelerometers and gyroscopes. Different approaches have been proposed to effectively handle these sensors, as they introduce biases and errors, and a simple integration of their outputs can not provide any accurate position or orientation (attitude) estimations [START_REF] Michel | Attitude estimation for indoor navigation and augmented reality with smartphones[END_REF]. One possible solution is foot-mounted dead-reckoning using the Zero-Velocity Update Technique (ZUPT) [START_REF] Fourati | Heterogeneous data fusion algorithm for pedestrian navigation via foot-mounted inertial measurement unit and complementary filter[END_REF]. However, this approach needs a very specific parameters tuning that depends on the nature of the walking activity, which makes the ZUPT unsuitable to general cases. Another way of exploiting IMUs is by combining them with magnetic sensors. In fact, All authors are with Univ. Grenoble Alpes, CNRS, Grenoble INP, GIPSAlab, F-38000 Grenoble, France. This work has been partially supported by MIAI@Grenoble Alpes, (ANR-19-P3IA-0003).

in indoor environments, important magnetic perturbations are observed [START_REF] Bachmann | An investigation of the effects of magnetic variations on inertial/magnetic orientation sensors[END_REF], that hold valuable information for indoor localization. One of the innovative techniques that takes advantage of these perturbations, has been firstly introduced in [START_REF] Vissière | Using magnetic disturbances to improve IMU-based position estimation[END_REF] and [START_REF] Vissière | Using spatially distributed magnetometers to increase IMU-based velocity estimation in perturbed areas[END_REF]. It states that as long as the magnetic field gradient is non-singular, the velocity can be reconstructed, from only a 3-axis magnetometers array, a 3-axis accelerometer and a 3-axis gyroscope measurements. This innovative technique can be applied on any moving object in indoor environments, and doesn't require any particular sensors location. In [START_REF] Vissière | Using magnetic disturbances to improve IMU-based position estimation[END_REF], authors considered that the magnetic field gradient is unknown and assumed that it satisfies some first order dynamics driven by white noises. In [START_REF] Vissière | Using spatially distributed magnetometers to increase IMU-based velocity estimation in perturbed areas[END_REF], the same authors extended their approach by measuring the magnetic field gradient (using three 3-axis magnetometers) along side with the previous assumption (white noise-based dynamic equation).

To the best of the authors knowledge, this choice is not fit to capture the magnetic field gradient dynamics and can consequently degrade velocity estimation. For this reason, in [START_REF] Zmitri | Improving inertial velocity estimation through magnetic field gradient-based extended kalman filter[END_REF], a magnetic field gradient-based Extended Kalman Filter (EKF) is proposed, that better models the magnetic field gradient dynamics by using a set of spatially distributed magnetometers array, and as a result, velocity estimation is more accurate. One of the well-known difficulties of Kalman filtering is its covariance matrices tuning (both process Q and measurement N covariances). Not providing the adequate values to the different noise parameters in these matrices causes unbounded errors and thus, unexpected estimation results. In fact, there is no general solution to properly tune these covariance matrices and usually, it is based on empirical methods. More particularly, giving fixed noise parameters to the measurement covariance N is unsuitable, as the errors in sensor measurements vary in time, due to different external factors (e.g temperature). In [START_REF] Mahmoud | Tuning of the error covariance parameters in EKF-based INS/GPS systems: A practical approach[END_REF], different tuning methodologies were demonstrated, to accurately handle an EKF-based INS/GPS fusion. In [START_REF] Brossard | AI-IMU dead-reckoning[END_REF], authors used a CNN to tune the noise parameters of an Invariant Extended Kalman Filter (IEKF), for the dead-reckoning of wheeled vehicles. This was done by exploiting the assumptions on the velocity that were presented in the form of pseudo-measurements, and controlling their corresponding covariance matrices. On the basis of this idea, an original architecture is proposed in this paper to improve the performance of the magnetic field gradient-based EKF. The approach dynamically adapts the measurement covariance matrix N of the EKF, using deep learning networks, in the purpose of having a better inertial velocity estimation. To do so, a LSTM model is trained to compute a pseudo-measurement of inertial velocity. The latter is then fed to the EKF and used to properly set the degree of confidence given to it by the filter. This is encoded in its corresponding covariance matrix that is set automatically by the proposed adapter. To evaluate the performance of the proposed approach, a comparison is achieved with the traditional INS/ZUPT method on the open source Openshoe dataset [START_REF] Skog | Foot-mounted INS for everybody -an open-source embedded implementation[END_REF]. First, a notable improvement on the velocity estimation is displayed, compared to when only the magnetic field gradient-based EKF with fixed covariance matrices is used. Then, the estimation results of the proposed approach are shown to be very close to the INS/ZUPT results. This is considered very promising, having in mind that unlike ZUPT, the proposed method is independent from walking paces (i.e. walking slowly or running) and is applicable on any other kind of dataset, which doesn't have to correspond to a foot-mounted framework, and can be extended to other applications.

II. PROBLEM STATEMENT

The proposed work in this paper is related to the navigation of a rigid body in indoor environments, by relying only on measurements from Micro Electro Mechanical Systems (MEMS) inertial sensors, as well as a spatially distributed 3-axis MEMS magnetometer's array. The goal is to design an approach that estimates the attitude matrix R ∈ SO(3) (the group of rotations in 3D Euclidean space), the inertial velocity

v n = [v nx v ny v nz ] ∈ R 3×1 and eventually if possible the position M n = [M nx M ny M nz ] ∈ R 3×1
from a given starting configuration (R 0 , v 0 , M 0 ).

A. Equations of motion

Two frames of coordinates are used to deal with the problem cited above:

• a local inertial frame n fixed to the Earth and its associated orthonormal basis B n = ( i n , j n , k n ); • a body frame b attached to the moving body and its associated orthonormal basis

B b = ( i b , j b , k b ).
Variables expressed in n (resp. b ) are marked by the subscript n (resp. b).

The attitude matrix R represents the rotation between these two frames, from n to b . This matrix can be expressed in terms of quaternion, denoted by q, such as R = R(q) [START_REF] Zmitri | Improving inertial velocity estimation through magnetic field gradient-based extended kalman filter[END_REF]. The quaternion is a hypercomplex number of rank 4 such that q = [q 0 q vect ] , where q 0 is the scalar part and q vect = [q 1 q 2 q 3 ] is the vector part. The reader is invited to refer to [START_REF] Kuipers | Quaternions and rotation sequences[END_REF] for more details about quaternion algebra.

The main equation that describes the variation of quaternion, can be defined from angular velocity measurements ω = [ω x ω y ω z ] ∈ R 3×1 of b with respect to n , given by a 3-axis gyroscope such as,

dq dt = 1 2 [ω q ×]q (1)
where ω q = [0 ω ] and [ω q ×] is its skew-symmetric matrix.

To determine the inertial velocity v n and position M n , the following kinematic model is introduced,

dv n dt = R a b -g ; dM n dt = v n (2)
where a b = [a bx a by a bz ] ∈ R 3×1 is the acceleration in b , measured by a 3-axis accelerometer. In this work the effects of Earth rotation and Coriolis acceleration are ignored, and the gravity vector g = [0 0 g z ] ∈ R 3×1 is considered known.

B. The magnetic field gradient to observe velocity

The following variables can be introduced at first:

• The magnetic field measured in b by a 3-axis magnetometer, which depends on time and space,

B b = [B bx B by B bz ] ∈ R 3×1 ; • The Jacobian matrix ∇B b ∈ R 3×3 , which represents
the gradient of the magnetic field, measured on a fixed point P b = [x bp y bp z bp ] ∈ R 3×1 and defined by

∇B b (P b (t)) = ∂B b (P b (t)) ∂P b (t) (3)
In this work, the magnetic field disturbances, that are often observed indoors [START_REF] Bachmann | An investigation of the effects of magnetic variations on inertial/magnetic orientation sensors[END_REF], are used to generate accurate estimates of velocity. This technique was firstly introduced in [START_REF] Vissière | Using magnetic disturbances to improve IMU-based position estimation[END_REF], and it resides on deriving Maxwell's equations [START_REF] Jackson | Classical Electrodynamics. Third Edition[END_REF] to obtain the following dynamic equation1 of the sensed magnetic field in b ,

dB b dt = -ω × B b + ∇B b v b (4)
where

v b = Rv n ∈ R 3×1 is the velocity in b .
As long as ∇B b is non-singular, v b is observable and can be estimated (see [START_REF] Vissière | Using spatially distributed magnetometers to increase IMU-based velocity estimation in perturbed areas[END_REF] for the observability proof). However, a remaining issue in this case is to reliably measure ∇B b .

In [START_REF] Vissière | Using magnetic disturbances to improve IMU-based position estimation[END_REF] and [START_REF] Vissière | Using spatially distributed magnetometers to increase IMU-based velocity estimation in perturbed areas[END_REF], model ( 4) is considered and ∇B b measurements are computed from a magnetometers array and assumed to satisfy a white noise-based model. In this case, the following points can be mentioned:

• the obtained ∇B b from the magnetometers array may contain a large noise coming first, from each noisy magnetometer measurement, and second, from the approximations taken into account to compute the different Jacobians representing the elements of ∇B b ; • modeling ∇B b components with a white noise is not sufficient to capture full information on their dynamics. These two arguments suggest that an accurate representation of ∇B b needs to be derived, or else, velocity estimation will suffer from unbounded errors especially when ∇B b has low values, as shown in [START_REF] Chesneau | Improving magneto-inertial attitude and position estimation by means of a magnetic heading observer[END_REF]. In [START_REF] Zmitri | Improving inertial velocity estimation through magnetic field gradient-based extended kalman filter[END_REF], a novel method that deals with this issue is proposed, through introducing a new equation that models efficiently the dynamics of ∇B b , using higher order spatial derivatives of the magnetic field,

d∇B b dt = T b v b + ∇B b [ω×] -[ω×]∇B b ( 5
)
where T b ∈ R 3×3×3 is a tensor representing the first spatial derivative of ∇B b in b . The reader can refer to [START_REF] Zmitri | Improving inertial velocity estimation through magnetic field gradient-based extended kalman filter[END_REF] for more information on how T b is measured.

C. Continuous-time dynamic model

To estimate attitude, inertial velocity and eventually position, a magnetic field gradient-based EKF alongside with a LSTM model are proposed. The EKF is established using the continuous-time dynamic model represented as follows:

                         dq dt = 1 2 [ω q ×]q dv n dt = R a b -g dB b dt = -ω × B b + ∇B b Rv n d∇B b dt = T b Rv n + ∇B b [ω×] -[ω×]∇B b (6)
The state vector for this dynamic model is

x = [q v n B b ∇B b ] ∈ R 15×1 , the input vector is u = [ω a b T b ] ∈ R 13×1
, and the output (measurement)

vector is y = [B b ∇B b ] ∈ R 8×1 .
Recall that only 7 elements of T b are sufficient to calculate all the tensor's components, see [START_REF] Zmitri | Improving inertial velocity estimation through magnetic field gradient-based extended kalman filter[END_REF]. This model is used to design the magnetic field gradient-based EKF, however, minor changes are applied when combining the latter with a LSTM model, in particular, a new element is added to the measurement vector. This proposed EKF-LSTM estimation approach is explained further in the next section.

III. EKF-LSTM ESTIMATION APPROACH

In this section, the proposed estimation schema is detailed. A LSTM-based learning model is used to dynamically tune the measurement covariance matrix N ∈ R 11×11 of the magnetic field gradient-based EKF. The approach which consists of three main blocks, is represented in Fig. 1, and is summarized as follows:

• a LSTM-model is trained to give a pseudomeasurement of inertial velocity, denoted

v lstm = [v nx,lstm v ny,lstm v nz,lstm ] ∈ R 3×1 ;
• v lstm is added to the output vector of the magnetic field gradient-based EKF, represented by the continuous-time dynamic model ( 6), such as

y = [B b ∇B b v lstm ] ∈ R 11×1 ;
• the measurement covariance matrix adapter determines the most suitable noise parameters of the measurement covariance matrix N corresponding to v lstm ; • EKF uses the continuous-time dynamic model [START_REF] Hochreiter | GLong short-term memory[END_REF], inertial and magnetic measurements, v lstm and N to better estimate the states, especially the magnetic field gradient and the velocity. In the following subsections a detailed explanation of the two blocks: the LSTM model and the measurement covariance matrix adapter, used in the magnetic field gradient based-EKF, to estimate The use of a deep learning algorithm inside the proposed approach is done to obtain a pseudo-measurement of inertial velocity. This pseudo-measurement is then used in the EKF, and the level of confidence given to this measurement is controlled by N , as it will be explained further in the next subsection. After reviewing the literature, it appears that the LSTM is one of the best learning solutions. This choice can actually be straightforward as LSTM-based models are considered very suitable for time series data, and they do not require any complex implementation. LSTM is an extension of artificial Recurrent Neural Network (RNN) architecture [START_REF] Hochreiter | GLong short-term memory[END_REF]. Its particularity is that it has feedback connections, which distinguish it from standard feedforward neural networks [START_REF] Sazli | A brief review of feed-forward neural networks[END_REF]. These feedback connections qualify the network to not only process single data points, but also an entire sequence of data, by retaining data information over long time periods.

x = [ q v n B b ∇B b ] .

B. Measurement covariance matrix adapter

The pseudo-measurement of inertial velocity v lstm is fed from the trained LSTM-based model to the magnetic field gradient-based EKF. At this point, the role of the measurement covariance matrix adapter, represented with the orange block in Fig. 1, is to dynamically assess the noise parameters in N , that reflect confidence in the new measurement v lstm . The measurement covariance matrix can be represented with N = diag(N 1 , N 2 , N 3 ), where N 1<i≤3 are the measurement covariance matrices corresponding to B b , ∇B b and v lstm respectively. In this work, the main interest is in N 3 ∈ R 3×3 , the covariance matrix of the measurement v lstm . Regarding the other noise parameters (i.e. N 1<i≤2 ), their values are fixed taking into account the noise standard deviation, indicated in the magnetometer datasheet, and their dynamic adaptation is left for future work. The same goes for the process covariance matrix Q ∈ R 15×15 . Algorithm 1 represents the measurement covariance matrix adapter, with α ∈ R >0 a multiplication factor, and ∈ R >0 a chosen threshold that indicates whether the amplitude of the LSTM predicted velocity is sufficiently low. Then, the EKF estimates v n using the magnetic field gradient-based technique, aided by the LSTM estimate v lstm . This is clearly understood when large values are assigned (equal to αI 3 ) to N 3,i . On the other hand, when v lstm has an amplitude that is close to zero, the EKF gives it more EKF initialization; for i = 1 : length (testdata) do N 3,i = αI 3 ; EKF prediction; if v lstm,i ≤ then N 3,i = I 3 ; end EKF update; end Algorithm 1: Measurement covariance matrix adapter confidence by attributing lower values (equal to I 3 ) to N 3,i . The reason for this strategy is the fact that v lstm prediction is more accurate when its 3-axis values are close to zero, than when they are not. This behavior can be explained with the insufficiency of data that corresponds to high amplitudes of inertial velocity in the used training dataset. To back up this assumption, the impact of the high amplitudes training data percentage on the predicted velocity v lstm accuracy is analyzed using the boxplot [START_REF] Williamson | The box plot: A simple visual method to interpret data[END_REF] in Fig. 2.

IV. EXPERIMENTAL SCENARIO AND RESULTS

The performance of the magnetic field gradient-based EKF is presented in this section, and how the inertial velocity estimation results should be improved. Then, the proposed EKF-LSTM approach is applied and it exhibits better performance on the velocity estimation and it impacts the position determination.

A. Magnetic field gradient-based EKF

In order to assess the performance of the magnetic field gradient-based EKF, 50% of the data provided by Openshoe is used2 . Under Matlab, a comparison between the proposed algorithm, and the INS/ZUPT-based velocity estimation [START_REF] Skog | Evaluation of zero-velocity detectors for foot-mounted inertial navigation systems[END_REF] is achieved. In Fig. 3, the estimation results for v n are plotted. The proposed EKF succeeds to estimate accurate inertial velocity at some phases of the walking scenario, without any additional velocity correction as the case for the INS/ZUPTbased method. Despite the intermittent errors that are seen in the comparison above, this first result is very promising, as it doesn't rely on any foot-mounted related technique, and is independent from the walking pace, as the case for ZUPT. To the best of authors knowledge, the proposed method can be extended to any other type of applications (drones, wheeled-robots, etc.), that are in search for reliable methods for inertial velocity estimation.

There are multiple sources of errors that are seen in Fig. 3. For instance, the different approximations taken into account to extract the spatial derivatives of ∇B b (T b for example) can lead to some errors, others can be the result of linearization process of the EKF, etc. One of the additional encountered difficulties when dealing with the Kalman filter, is the tuning of its covariance matrices (i.e. process and measurement covariances). Then, using a fixed N , is not suitable for situations where the error in sensor readings varies in time due to external factors (e.g. temperature). Based on this idea, an EKF-LSTM system is proposed that generates a pseudo-measurement of inertial velocity, feeds it to the magnetic field gradient-based EKF, and uses its dynamic to tune its corresponding covariance matrix.

B. Magnetic field gradient-based EKF aided by LSTM network

The implementation details and evaluation results of the proposed EKF-LSTM approach are detailed in the following, starting from the training process until the estimation one. The result of the reconstructed trajectory corresponding to the estimated velocity is also analyzed and discussed.

For the training process, the Adam optimizer [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF] is used with an initial learning rate of 0.01, and the model is trained for 100 epochs. Before training, a standardization [START_REF] Anysz | The influence of input data standardization method on prediction accuracy of artificial neural networks[END_REF] is applied on the Openshoe dataset, which consists basically in subtracting the arithmetic mean from the given data and dividing it by the standard deviation. Then, 50% of this data are used to form the training set and the other 50% are for the testing. The 50% training portion is equivalent to 7038 samples, which is considered small for a training process. This is the advantage of the proposed approach, that it doesn't require any complex artificial intelligence step that might need a large time computation or memory usage. The result of inertial velocity determination using the LSTM, applied on the testing set is displayed in Fig. 5. As expected, the accuracy of the determined velocity from the LSTM during phases where it has large amplitudes is lower than the phases where the amplitude is close to zero. As demonstrated in Fig. 2, this is due to the usage of a training set that lacks high amplitudes velocity data. 2) EKF-LSTM estimation results: The velocity estimate v lstm is fed to the measurement vector of EKF (as pseudomeasurement). Within its loop, algorithm 1 is ran, to dynamically adapt the corresponding covariance matrix. The threshold is 0.1ms -1 and the multiplication factor α is 10 4 . These values correspond to the tested data, and can change according to the studied dataset. The performance of the proposed EKF-LSTM-based estimation is evaluated on the testing set, by comparing the obtained results to the INS/ZUPT algorithm ones. In Fig. 6, the 3-axis components of inertial velocity v n are plotted. The velocity estimate v n (blue solid line) given by the approach is now very close to the INS/ZUPT estimate (red solid line). The RMSE between the two results is around 0.23ms -1 , which is very promising considering that the proposed method does not use any footmounted related estimation techniques. To analyze the effect of velocity estimation on the trajectory reconstruction, v n is integrated and the 2D estimated position M n is plotted and represented in Fig. 7. The estimated trajectory (blue solide line) is close to the tested one (black solid line), and the same squared shape is maintained. Table I presents the traveled distance, as well as the coordinates of the arrival points for both trajectories. The arrival point of the tested trajectory is different from the starting point, that is the same for both trajectories. The error of the final arrival point between the two trajectories is around 1.8m which competes with top ranked methods in the literature, that require either a costly and heavy preinstalled infrastructure or a very specific sensor location (as the case for ZUPT). It is also important to highlight that in this work, the Openshoe INS/ZUPT estimation results are considered as the ground truth. However, no conclusions can be made on their accuracy, as no actual ground truth is provided. This means that results given by the proposed approach, can be closer to reality than the ones computed from the INS/ZUPT approach.

V. CONCLUSIONS

In this paper, the inertial velocity and position estimation from a magnetic field gradient-based EKF aided by a LSTM algorithm was proposed. The LSTM model was introduced to provide a pseudo-measurement of velocity for the EKF, and a covariance matrix adapter was applied to dynamically tune its corresponding noise parameters. The final inertial velocity estimate from the proposed approach competes with the INS/ZUPT results, as it requires only IMU and magnetic measurements, and doesn't rely on any particular sensor location or walking pace parameterization. It was shown that the reconstructed trajectory from this novel method is very close to the one from INS/ZUPT technique, and small drifts were observed. Other datasets corresponding to different applications (e.g. drones, land robots, etc.) need to be tested to highlight more the efficiency of the proposed solution.

Fig. 1 :

 1 Fig. 1: The proposed EKF-LSTM estimation system

  Fig. 2: Effect of high amplitudes training data percentage on RMSE of predicted velocity. Larger the percentage of training data, the smaller velocity RMSE of the testing data becomes.

1 ]Fig. 3 :

 13 Fig. 3: Estimation of v n with the magnetic field gradientbased EKF and without the LSTM algorithm

1 )Fig. 4 :

 14 Fig. 4: Diagram of the used LSTM network

1 ]Fig. 5 :

 15 Fig. 5: Determination of v n with only the LSTM algorithm

1 ]Fig. 6 :Fig. 7 :

 167 Fig. 6: Estimation of v n with the EKF-LSTM system

  Note that the 30% corresponding to the testing set are always fixed. As the training data percentage increases, the RMSE of the predicted velocity of the testing set gradually decreases (for instance by looking at the median values represented by the red horizontal lines). This result indicates that as more inputs are provided to the LSTM during training for the dynamic phases, the prediction of v lstm is more accurate. Nevertheless, high accuracy of the deep learning network is not mandatory in this work. A small training database is sufficient to detect when v lstm is close to zero. The low accuracy during the dynamic phases is handled through the tuning of the magnetic field gradient-based EKF with the measurement covariance matrix adapter.

	2: Effect of high amplitudes training data percentage
	on RMSE of predicted velocity. Larger the percentage of
	training data, the smaller velocity RMSE of the testing data
	becomes.
	For this experiment, a dataset with 5913 samples is used.
	It is created from the Openshoe dataset and contains sensor
	measurements corresponding to dynamic phases with high
	velocity amplitudes. This dataset is divided into 70% (equiv-
	alent to 4139 samples) for the training, and 30% (equivalent
	to 1774 samples) for the testing. Then the Root Mean Square
	Error (RMSE) between the INS/ZUPT-based velocity and
	the LSTM predicted one (of the testing set) is measured for
	different used percentages of the training dataset. The process
	is repeated 20 times for each training data percentage.

× is the cross product of two vectors in R 3 .

The dataset provides measurements for only one

3-axis magnetometer. The rest of the magnetic array data is simulated.