Ke Wu

Gang Zheng

Junfeng Zhang

FEM-based Trajectory Tracking Control of a Soft Trunk Robot ?

de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

As a novel class of robots, soft robots have demonstrated many desirable mechanical properties than traditional rigid robots [START_REF] Trivedi | Soft robotics: Biological inspiration, state of the art, and future research[END_REF] due to their nature of being compliant, flexible and hyper-redundant. For instance, these robots provide great adaptability to complex environments [START_REF] Trivedi | Soft robotics: Biological inspiration, state of the art, and future research[END_REF], safe human robot interaction (HRI) [START_REF] Rus | Design, fabrication and control of soft robots[END_REF], energy-saving actuation [START_REF] Kim | Soft robotics: a bioinspired evolution in robotics[END_REF], the maneuverability to display diverse mechanical properties [START_REF] Trivedi | Soft robotics: Biological inspiration, state of the art, and future research[END_REF][2], a large power-to-weight ratio [START_REF] Andrew D Marchese | Dynamics and trajectory optimization for a soft spatial fluidic elastomer manipulator[END_REF] [START_REF] Morales Bieze | Finite element method-based kinematics and closed-loop control of soft, continuum manipulators[END_REF] and increased fault tolerance [START_REF] William | Continuum robot dynamics utilizing the principle of virtual power[END_REF] [START_REF] Isuru S Godage | Dynamics for variable length multisection continuum arms[END_REF]. Traditional rigid-body robots which consist of rigid links connected at kinematic joints, serve as mechanical devices to transfer motion, force and energy by the synthesized motion of rigid members whereas soft robots are designed to function similarly just through the elastic deformation of their flexible bodies. Due to the redundant degrees of freedom (DoF), these soft robots are more maneuverable than conventional rigid robots in terms of achieving some complex tasks that need unique and more compliant configurations. For example, this type of robots can be well-used in rescue missions [START_REF] Fang | Origami-based earthworm-like locomotion robots[END_REF], pipeline flaw detection [START_REF] Grazioso | A geometrically exact model for soft continuum robots: The finite element deformation space formulation[END_REF][10] [START_REF] Boyer | Macrocontinuous dynamics for hyperredundant robots: application to kinematic locomotion bioinspired by elongated body animals[END_REF], object grasping [START_REF] Li | Design and analysis of a square-shaped continuum robot with better grasping ability[END_REF], medical surgery [START_REF] Burgner-Kahrs | Continuum robots for medical applications: A survey[END_REF], space tasks [START_REF] Manas M Tonapi | A novel continuum robotic cable aimed at applications in space[END_REF] and so on.

On the other hand, the advantageous characteristics of these robots also lead to difficulties in deriving the kinematic and dynamic models of these manipulators [START_REF] Michael | Kinematics and the implementation of an elephant's trunk manipulator and other continuum style robots[END_REF]. Obviously, the control theory developed for rigid robots would hardly be feasible for soft robots, and the research on controlling those robots therefore remains huge possibility for exploration.

In the state-of-art literature, many research efforts for modeling and controlling continuum and soft robots (manipulators) have been reported. The contributors of [START_REF] Michael | Kinematics and the implementation of an elephant's trunk manipulator and other continuum style robots[END_REF] took the lead to propose an effective geometric modeling method called Constant-Curvature-assumed model (CC) for continuum manipulators. Essentially speaking, this method masterly utilizes the nature of Denavit-Hartenberg method and modifies it into an effective method that's able to map the relationship between these robots' configuration space and operational space. Due to this merit, CC model once became a popular choice of designing static controllers for continuum robots [START_REF] Michael | Kinematics and the implementation of an elephant's trunk manipulator and other continuum style robots[END_REF] [START_REF] Bryan | Kinematics for multisection continuum robots[END_REF]. However, in some cases CC assumption becomes less accurate for the whole slender body of a continuum robot. Then, Piecewise Constant Curvature model (see details in [START_REF] Webster | Design and kinematic modeling of constant curvature continuum robots: A review[END_REF]) was proposed to take over the modeling problem. PCC basically assumes the continuum robot can be divided into several independent segments where the bending of each segment is characterized by a circular arc of constant length (CC model applied in each segment). This simplification has made PCC one of the most popular and effective methods for modeling continuum manipulators [START_REF] Webster | Design and kinematic modeling of constant curvature continuum robots: A review[END_REF]. For example, it is not only used in controlling soft robots that is main scope of this paper: kinematic control [START_REF] Bryan | Kinematics for multisection continuum robots[END_REF], feedforward dynamic control [START_REF] Valentin Falkenhahn | Model-based feedforward position control of constant curvature continuum robots using feedback linearization[END_REF], feedback dynamic control [START_REF] Della | Dynamic control of soft robots interacting with the environment[END_REF][20], configuration control [START_REF] Isuru S Godage | Modal kinematics for multisection continuum arms[END_REF], trajectory tracking or path following [START_REF] Mahl | A variable curvature continuum kinematics for kinematic control of the bionic handling assistant[END_REF], but also applied to the cases of mechanical design [START_REF] Runge | A framework for the automated design and modelling of soft robotic systems[END_REF], sensing [START_REF] Kim | Optimizing curvature sensor placement for fast, shape of continuum robots[END_REF] and workspace study [START_REF] Taha Chikhaoui | Kinematics and performance analysis of a novel concentric tube robotic structure with embedded soft micro-actuation[END_REF]. Cosserat theory was also employed as a feasible option to deduce the kinematics of the trunk-like soft robots [START_REF] Jones | Three dimensional statics for continuum robotics[END_REF][27] [START_REF] Zhao | A continuum manipulator for continuously variable stiffness and its stiffness control formulation[END_REF], and it's also extended to capture the dynamics of continuum robots [START_REF] Renda | Dynamic model of a multibending soft robot arm driven by cables[END_REF] [START_REF] Abu Alqumsan | Robust control of continuum robots using cosserat rod theory[END_REF]. Some representative cases of Cosserat-model-based controllers have been reported: robust sliding mode dynamic control [START_REF] Abu Alqumsan | Robust control of continuum robots using cosserat rod theory[END_REF], dynamic configuration tracking [START_REF] Doroudchi | Configuration tracking for soft continuum robotic arms using inverse dynamic control of a cosserat rod model[END_REF], synthesized control for parallel continuum manipulators [START_REF] Till | Efficient computation of multiple coupled cosserat rod models for real-time simulation and control of parallel continuum manipulators[END_REF] and so on. For some slender continuum compliant robots, Euler-Bernoulli beam theory can be used to develop the kinematics [START_REF] Xu | Analytic formulation for kinematics, statics, and shape restoration of multibackbone continuum robots via elliptic integrals[END_REF] and dynamics [START_REF] Yang | A modular approach for dynamic modeling of multisegment continuum robots[END_REF] of the robots, and their static and dynamic controllers can be therefore developed accordingly. Recently, Finite Element Method (FEM) has gradually become a promising option of modeling soft or continuum robots [START_REF] Duriez | Control of elastic soft robots based on real-time finite element method[END_REF] but it's rather computationally expensive due to its high-dimensional nature. Because of this fact, some order reduction methods have been developed to deduce its computational cost [START_REF] Goury | Fast, generic, and reliable control and simulation of soft robots using model order reduction[END_REF]. However, compared to the mentioned methods, FEM does not need the studied deformable object to have a certain regular geometry (such as rod-like or trunk-like ones) since irregular geometries can be discretized into a group of fine mesh elements for further dynamic analysis using Euler-Lagrange method [START_REF] Singiresu | The finite element method in engineering[END_REF]. This nature has made FEM quite a universal strategy to model deformable robots of any arbitrary shape [START_REF] Duriez | Control of elastic soft robots based on real-time finite element method[END_REF].

Likewise, FE model has been proved to be option to develop controllers for soft robots: [START_REF] Zheng | Robust control of a silicone soft robot using neural networks[END_REF] introduced a robust control strategy using neural networks; a static observer-based control method was proposed in [START_REF] Zhang | Kinematic modeling and observer based control of soft robot using real-time finite element method[END_REF]; [START_REF] Zheng | Controllability pre-verification of silicon soft robots based on finite-element method[END_REF] proposed a way of modeling the deformation of soft robots and also checking its controllability; a feedback quasi-static control strategy for soft robots has been realized by applying FEM and the visual tracking technique in [START_REF] Zhang | Vision-Based Sensing of External Forces Acting on Soft Robots Using Finite Element Method[END_REF]. In this paper, we also choose FEM as the modeling tool, and then promote our previous work [START_REF] Wu | Fem-based gain-scheduling control of a soft trunk robot[END_REF] where a FEM-based gainscheduling control strategy for a soft trunk robot. This former contribution focuses on global path following by smoothly switching the gains in each sub-workspace without considering time constraint. However, in this paper, we propose an effective dynamic controller that can realize time-varying and global trajectory-tracking for the same robot. Precisely, we first plan a feasible trajectory for the studied robot in SOFA (a FEM-based simulator) by solving a MPC-based optimization problem. Then, we conduct linearization of this robotic system around the pre-designed trajectory, and then an associated controller can be developed. We validate its feasibility and effectiveness by doing a comparison test between a common PID controller and this proposed one, and afterwards its robustness against temporary and permanent disturbances is proved by experimental testing as well.

Problem statement 2.1. Description of the studied soft trunk-like robot

In this paper, we still focus on the same robot studied in [START_REF] Wu | Fem-based gain-scheduling control of a soft trunk robot[END_REF] (see Fig. 1) to demonstrate the promotion on the control strategy. We fabricated the prototype with a long uniform hole space from its based to its tip where the position sensor can be properly settled at its tip. Four inextensible plastic cables, which are functioned by 4 steppers and mounted through the robot's body, are designed to drive the robot moving but the soft rubber can easily wear if cables repeatedly travels through the rubber. To reduce this undesired friction, we 3D-printed 15 rigid rings to be mounted along the robot's body (see Fig. 1) where the sliding friction happens between the rigid rings and cables instead. The mass of trunk robot m = 0.045 kg, Young's modulus E = 6.5 ù 10 6 Pa and Poisson ratio ⌫ = 0.45.

Control objective

The objective of this paper is to control the end-effector of the trunk-like robot to track a planned trajectory with or without the interference of external disturbances. As stated above, this work is regarded as an extension of our previous work [START_REF] Wu | Fem-based gain-scheduling control of a soft trunk robot[END_REF] to realize global trajectory tracking tasks in its feasible workspace.

Modeling and trajectory generation

Modeling via Finite Element Method

Soft robots are strenuous to model owing to its high-DoF nature and complex deformation of materials [START_REF] Duriez | Control of elastic soft robots based on real-time finite element method[END_REF]. In this paper, FEM is employed to obtain the dynamic model of the soft trunk manipulator on the basis of SOFA Framework, a FEM-based simulator. Therefore, the associated controller for trajectory tracking can be logically developed via synthesizing the information of the FE model. Concerning the derivation of FE dynamical model, the details can be found in [START_REF] Singiresu | The finite element method in engineering[END_REF] and [START_REF] Amehri | Fem based workspace estimation for soft robots: a forward-backward interval analysis approach[END_REF] where the basic idea is similar to the classical one: we first choose the type and the number of fine mesh elements (displayed in Fig. 2) to discretize the geometry of the modeled deformable object (which in this The studied trunk robot [START_REF] Wu | Fem-based gain-scheduling control of a soft trunk robot[END_REF] case is the studied robot); secondly, we calculate the displacement of any point inside the element via the displacement of its nodes through shape functions, and the corresponding strain and stress (might be linear or nonlinear according to the applied constitutive law); thirdly, we use either Euler-Lagrange equation or virtual work principal to derive the FE dynamical model for the studied manipulator. In this paper, we focus on the displacement control of the robot's end-effector (q e À R 3ù1 ), and q e is defined with respect to the inertial frame F I (q e ) as shown in Fig. 3. Therefore, this control problem is essentially an output regulation problem of the following nonlinear system based on the derived FE model:

M(q) á q + D( Ü q, q) Ü q + K(q)q = H(q) T u q e = Cq ( 1
)
where M(q) À R 3nù3n denotes mass matrix; D( Ü q, q) À R 3nù3n denotes damping matrix; K(q) À R 3nù3n denotes stiffness matrix; H(q) À R 4ù3n denotes the input control matrix; u À R 4ù1 denotes the tension forces exerted at the four cables; the friction between the cables and the rubber is ignored in our FE model so H(q) only distributes u to all the elements containing the cable nodes where the cables travel through within the rubber (red nodes displayed in Fig. 2); q À R 3nù1 denotes the displacement of all the n nodes (supposing q i is the ith element of q, q i refers to the displacement of the ith node with respect to its inertial frame F I (q i ) as shown in Fig. 3); logically, Ü q À R 3nù1 and á q À R 3nù1 denote the corresponding velocity and acceleration respectively; the function of matrix C À R 3ù3n is to pick out the end-effector displacement from all n (n = 1484) nodes defined in our derived FE model; the units of q and u are centimeter (cm) and Newton (N) respectively, and the rest units can be deduced accordingly.

In terms of the mesh type, tetrahedron elements (see Fig. 2) are chosen as it is compatible with complex geometry, and the robot's body is spatially discretized into n = 1484 nodes. The size of the mesh elements has been carefully chosen in order to deal with the trade-off between accuracy and computation expense.

Trajectory generation via SOFA

With the deduced dynamical model (1) of the investigated soft trunk robot, this subsection discusses how to generate a feasible end-effector trajectory via the FEM-based simulator (SOFA). Note that FE simulations are conducted in an off-line manner in this paper to complete two tasks (see details in Step 1 and Step 2 respectively in this chapter): numerically defining the workspace of the robot's end effector; numerically generating planned trajectory based on a MPC-based optimization structure. Because these two steps are both done in an off-line manner, the computational expense would not be a critical issue for this methodology. To treat this problem, let us firstly recall the definition of the end-effector workspace [START_REF] Amehri | Fem based workspace estimation for soft robots: a forward-backward interval analysis approach[END_REF].

Definition 1. Given a configuration of a soft manipulator dynamically modeled by (1) with a bounded actuator u À U , the workspace W E of the end-effector q e , is a subspace of R 3

, defined below:

W E = { q e À R 3  q e = Cq, K(q)q * H(q) T u = 0, ≈u À U } Remark 1.
The above definition of the workspace, which is similar to the one defined for a rigid manipulator, implies that W E contains all possible equilibrium positions for which the end-effector q e can reach and stay there. According to the above definition, given two positions in the end-effector's workspace, named as the initial position q init e À W E and the final position q f in e À W E , the trajectory generation problem for such a soft trunk robot is to seek a feasible path starting from q init e and ending with q f in e over a time interval [0, T ], i.e., q e (0) = q init e and q e (T ) = q f in e , and at the same time satisfy the dynamical constraint (1).

In the literature, differential flatness [START_REF] De | Trajectory planning and control for planar robots with passive last joint[END_REF] has been widely used to plan a feasible trajectory for different types of robots. Roughly speaking, the output of a nonlinear dynamical system is called to be flat if its state and input can be written as functions of this output and its derivatives. However, concerning the studied soft robot with dynamics (1), it is clear that the state q and the input u cannot be function of the output q e and its derivatives. Otherwise, it means that the displacement of each node of any mesh can be uniquely determined by the position of the end-effector and its derivatives, which is obviously not true in practice. In other words, the chosen output (the position of end-effector), is not flat, and therefore we cannot apply the existing differential-flatness-based method to generate feasible trajectories for the investigated soft robot.

This paper adopts numerical simulation method to solve this problem for the soft manipulator (as shown in Fig. 1) described in [START_REF] Trivedi | Soft robotics: Biological inspiration, state of the art, and future research[END_REF]. Precisely, such a method contains the following two steps:

Step 1 : Numerically evaluate W e , in which we determine the initial and final positions of the desired trajectory , i.e.,

q init e À W E and q f in e À W E ;

Step 2 : Freely choose a differentiable path in W e linking q init e and q f in e , and then generate the desired trajectory q r (t) for t À [0, T ] where T is the prescribed time interval.

It is worth noting that the above method is to generate the desired trajectory directly in the workspace of the endeffector (i.e., directly planning the end-effector's trajectory q r e in the output space by solving the inverse dynamics to get q r ), and it is possible as well to generate the desired trajectory in the configuration space of q (i.e., planning a trajectory of q r in the configuration space and then map it to get q r e ), like the trajectory generation for rigid manipulator. This paper adopts the direct manner since it is more visible to evaluate the controller performance which will be detailed in Section 5.

Step 1: In this step, in order to evaluate the end-effector's workspace, for all u À U , we solve the static model K(q)q * H(q) T u = 0 to compute the associated configuration q which then yields the associated end-effector position q e = Cq.

Precisely, denote the control input as

u = [u 1 , 5 , u 4 ] T for u i À U i with U i = [u i , u i ]
where u i and u i imply respectively the lower and upper bound of the ith input. Then, we discretize the set U i with a chosen fixed sampling

positive constant h = u i *u i k
where k represents the number of sampled intervals. Following such a discretization, U i has been sampled as

U s i = u i , u i + h, 5 , u i oe U i , i = 1, 2, 3, 4 It is clear that, for any given u À U s = U s 1 ùU s 2 ùU s 3 ùU s 4
, we can solve the algebraic equation K(q)q*H(q) T u = 0 by applying the classical gradient method, which yields

q (j) = q (j*1) * ( q ⌅ K(q)q * H(q) T u ⇧ ⌅ K(q)q * H(q) T u ⇧  q=q (j*1)
where q (j) means the jth iteration of the value q and ( q represents the gradient operator with respect to q. After the convergence (in practice the convergence criteria is set as q (j) * q (j*1)  f ✏ where ✏ is a pre-defined threshold), then we can obtain the corresponding end-effector position via q e = Cq (j) for the given input u.

It is worth noting that, for a given u, the solutions of K(q)q * H(q) T u = 0 might not be globally unique, and they heavily depend on the initial value q (0) of the applied gradient method. To overcome this problem, we start the evaluation procedure from the rest configuration where u i = 0, for which the corresponding value q is known. Then, for each evaluation we modify only one input. By using the final value obtained in the former evaluation as the initial guess of the new evaluation, we can guarantee that the applied gradient method can quickly converge to its local unique solution.

By exhaustively evaluate all inputs in U s , we can obtain the corresponding sampled end-effector workspace W s E oe W E (see Fig. 4b), whose precision can be tuned by changing the sampling step h. Note that the inertial frame of the end effector F I (q e ) is chosen as the global coordinate system of FEM simulations where O is the origin of this coordinate system (see Fig. 4a).

In order to generate a feasible trajectory for the end-effector, it is required that the chosen initial and final positions should satisfy q init e À W s E and q f in e À W s E

Step 2: After having obtained the sampled workspace W s E , we can obtain the boundary of the real workspace, and this enables us to choose a differentiable path, noted as q r e , belonging to W E , i.e.,

q r e oe W E

In the experimental part, we use SOFA to evaluate the workspace W E , and choose a path from the sampled workspace W s E obtained from SOFA between the q init e and q f in e :

q r e = (q r e 1

, q r e 2

, q r e 3

, 5 , q r e n

) Then, we assign the time interval [0, T ] for this path, which logically generates a user-designed trajectory q r e (t) with t = 0, ⌧, 2⌧...T . However, to achieve this planned q r e (t), we need to find the possible input and output of the studied robotic system that is able to approximate q r e (t) as close as possible, which yields the following optimization problem [START_REF] Gritschneder | Fast trajectory planning for automated vehicles using gradient-based nonlinear model predictive control[END_REF]:

min u J = T 0 ⌅ q r e (s) * Cq(s) 2 + u T (s)u(s) ⇧ ds (2) s.t. M á q(t) + D( Ü q(t), q(t)) Ü q(t) + K(q(t))q(t) = H(q(t)) T u(t) Cq(t 0 ) = q r e (t 0 ); C Ü q(t 0 ) = Ü q r e (t 0 ); u(t) À U
Essentially speaking, we use the FEM-based simulator SOFA as the numerical black box to simulate the dynamical constraint M(q) á q +D( Ü q, q) Ü q +K(q)q = H(q) T u, considering J as the objective function and u(t) as the input variables.

This process enables us to calculate the numerical gradient

( u J = ⌧ )J )u T
. Precisely, the following procedure is applied:

Step 1: Divide the time interval [0, T ] into i segments: 0 = t 0 < t 1 < 5 < t i*1 < t i = T with t = T i ; Step 2: Define the input u(t) as following:

u(t) = u k , ≈t À ⌅ t k , t k+1 ⇧ , k = 0, 1, 2 5 i * 1.
where

u k = u(t k ) = ⌅ u k,1 , 5 , u k,4 ⇧ T = ⌅ u 1 (t k ), 5 , u 4 (t k ) ⇧ T is treated as a constant for t À [t k , t k+1 ] with u k,s
being the value of the sth input (0 f s f 4) at the time t k , and construct the input sequence as a vector:

U c = ⌅ u T 0 , 5 , u T i*1 ⇧ T À R 4iù1
and concatenate the corresponding configuration as

Q c = ⌅ q T (t 0 ), 5 , q T (t i ) ⇧ T À R 3n(i+1)ù1
Step 3: For the jth iteration, with the knowledge of U j c , Q j c and Ü Q j c (for the initialization where j = 0 the robot is on its initial static configuration q init with the corresponding input u 0 , therefore

U j c = ⌅ u T 0 , 5 , u T 0 ⇧ T , Q j c = ⌅ [q init ] T , 5 , [q init ] T ⇧ T and Ü Q j c = 0), evaluate its cost function J j ;
Step 4: Based on U j c , generate a new sequence of control input:

É U j c = U j c + [0, 5 , 0, j l , 0, 5 , 0] T
where j l represents the variation of the lth component of U j . Then using SOFA to simulate the dynamics (Eq. ( 1)) to get the new corresponding configuration É Q j c and Ü É Q j c which enables us to compute the new value of cost function É J j ;

Step 5: Calculate the numerical partial derivative of J with respect to the lth component of U c as

É J j *J j j l
, and concatenate all numerical partial derivatives for 1 f l f 4i to obtain the numerical gradient as

( U j c J j = 4 )J j )U j c 5 T .
Step 6: Applying projected gradient descent or projected Newton method to calculate U j+1 c , and use SOFA to calculate

Q j+1 c and Ü Q j+1 c
which enables us to compute the cost function J j+1 ;

Step 7: If J j+1 * J j  > ✏ with ✏ being a pre-defined threshold, then go to Step 3 with j = j + 1 for another round of iteration. Otherwise, the procedure is terminated and the desired feasible trajectory is

Q r = Q j+1 c with Ü Q r = Ü Q j+1 c
and the corresponding open-loop controller is U r = U j+1 c . Note that the obtained Q r , Ü Q r and U r are discrete, based on which a continuous open-loop controller u r (t) À C 0 and a differentiable trajectory q r (t) À C 1 with Ü q r (t) À C 0 need to be generated. For this, interpolation technique can be applied.

Precisely, to construct a continuous u r (t

) from U r = ⌅ u T 0 , 5 , u T i*1
⇧ T , the following linear interpolation is used:

u r (t) = Ç u(t k+1 ) * Ç u(t k ) t (t * t k ) + Ç u(t k ), k = 0, 1, 2 5 i * 1; (3) with Ç u(t 0 ) = u 0 , Ç u(t i ) = u i*1 and Ç u(t k ) = [u k + u k*1 ]_2 for 1 f k f i * 1. Similarly, to compute a differentiable trajectory q r (t) from Q r = ⌅ q T (t 0 ), 5 , q T (t i ) ⇧ T and Ü Q r = ⌅ Ü q T (t 0 ), 5 , Ü q T (t i ) ⇧ T ,
a cubic polynomial is adopted to calculate q r (t) for each interval t À ⌅ t k , t k+1 ⇧ with 0 f k f i*1 since we have position and velocity constraints at t k and t k+1 . Such an interpolation yields:

q r (t) = < b k0 + b k1 t + b k2 t 2 + b k3 t 3 , ≈t À ⌅ t k , t k+1 ⇧ 0, otherwise (4)
with b km for 0 f m f 3 being the 4 unknown polynomial coefficients on the kth interval. In total, we have 4i unknown parameters to be identified. Due to the fact that q r (t) should satisfy the boundary conditions of each interval t À ⌅ t k , t k+1 ⇧ , thus this yields the following 4i linear algebraic equations to be fulfilled:

For t = t k : b k0 + b k1 t k + b k2 t k 2 + b k3 t k 3 = q r (t k ) and b k1 + 2b k2 t k + 3b k3 t 2 k = Ü q r (t k ) For t = t k+1 : b k0 + b k1 t k+1 + b k2 t 2 k+1 + b k3 t 3 k+1 = q r (t k+1 ) and b k1 + 2b k2 t k+1 + 3b k3 t 2 k+1 = Ü q r (t k+1 )
Solving the above linear equation enables us to determine the values of b km , then the differentiable trajectory q r (t) and its associated derivative Ü q r (t) can be obtained via (4).

Controller design via trajectory linearization

Given the planned end-effector trajectory Cq r via (4) and its corresponding control u r defined in (3), the first intuitive manner to design trajectory-tracking controller is the well-known trajectory linearization method. The main idea of this approach is to conduct linearization of the FE model along the pre-obtained trajectory, and then design the controller based on the corresponding u r .

Precisely, for the soft manipulator described in (1), we can rewrite it into the following projected space:

á q e + ↵(q, Ü q) = CM(q) *1 H(q) T u [START_REF] Morales Bieze | Finite element method-based kinematics and closed-loop control of soft, continuum manipulators[END_REF] with ↵(q, Ü q) = CM(q) *1 D( Ü q, q) Ü q + CM(q) *1 K(q)q (6)

Then, we introduce new variables:

X = 4 x 1 x 2 5 = 4 q e Ü q e 5 (7)
so system ( 5) can be transformed into the state-space form:

Ü X = f (X, u) Y = P X (8)
where

f (X, u) = 4 x 2 CM(q) *1 H(q) T u * ↵ 5 (9)
and

P = [I, 0] (10)
Based on the desired trajectory Cq r (t) and its corresponding controller u r (t) obtained via (2) for a pre-defined time interval [0, T ], we can then define

X r (t) = 4 Cq r (t) C Ü q r (t) 5 and Y r (t) = Cq r (t), t À [0, T ] (11)
Then, by noting

X(t) = X r (t) + X(t) and u(t) = u r (t) + u(t) (12)
and

Y (t) = Y r (t) + Y (t) (13)
system ( 8) can be then linearized around the designed trajectory as follows [START_REF] Hespanha | Linear systems theory[END_REF]:

Ü X(t) = A X(t) + B u(t) Y (t) = P X(t) (14) with A = )f (X, u) )X  X=X r ,u=u r (15)
and

B = )f (X, u) )u  X=X r ,u=u r = 4 0 CM(q r ) *1 H(q r ) T 5 (16)
where q r is obtained via (2), representing the corresponding displacement of all nodes along the pre-designed endeffector trajectory.

Since P = [I, 0], then system ( 14) can be written as:

á Y (t) = P A 2 X(t) + P AB u(t) = Ñ ↵(t) + CM(q r ) *1 H(q r ) T u(t) (17)
with

Ñ ↵(t) = P A 2 X(t) (18)
In the above equation ( 17), we can regard Ñ ↵(t) as a vanishing external bounded disturbance, since Ñ ↵(t) = 0 when Y = 0. As a result, we can propose the following robust controller:

u(t) = u r (t) + u(t) (19)
with

u(t) = (CM(q r ) *1 H(q r ) T ) *1 u e (e q e , Ü e q e ) (20)
where the observation error e q e is defined as follows:

e q e = Y r (t) * Y (t) = * Y (t) (21)
and the term u e (e q e , Ü e q e ) has a classic PID form:

u e (e q e , Ü e q e ) = k p e q e + k i t 0 e q e ds + k d Ü e q e [START_REF] Mahl | A variable curvature continuum kinematics for kinematic control of the bionic handling assistant[END_REF] where the diagonal gain matrices k p , k i , k d must have the appropriate dimension. Considering [START_REF] Webster | Design and kinematic modeling of constant curvature continuum robots: A review[END_REF] to ( 22) together, we can finally reach: á e q e + k d Ü e q e + k p e q e + k i t 0 e q e ds = Ñ ↵(t) [START_REF] Runge | A framework for the automated design and modelling of soft robotic systems[END_REF] which is a classical 3-order linear system with a bounded disturbance Ñ ↵(t) which will be vanished when e q e = 0. Therefore, we can always find matrices k p , k i , k d to guarantee that e q e converges to 0.

When the robot is perpetuated by disturbance, then the derived closed-loop system can be written as á e q e + k d Ü e q e + k p e q e + k i t 0 e q e ds = Ñ ↵(t) + d(t) [START_REF] Kim | Optimizing curvature sensor placement for fast, shape of continuum robots[END_REF] where d(t) represents external disturbance. According to [START_REF] Kim | Optimizing curvature sensor placement for fast, shape of continuum robots[END_REF], the proposed controller provides robustness with respect to vanishing or piece-wise constant (or slowly time-varying) disturbance d(t) due to the integral term in the proposed controller and Ñ ↵(t) is vanishing. For other types of disturbance (such as quickly time-varying), if they are bounded, e q e will be bounded as well considering the input-state-stable (ISS) property provided by [START_REF] Kim | Optimizing curvature sensor placement for fast, shape of continuum robots[END_REF].

Experimental testing

The experimental testing regarding control feasibility and robustness was conducted to validate the proposed control strategy. The overall setup of the testing platform (see Fig. 5) mainly has three parts: the mechanical base, the Polhemus magnetic position sensor, the micro communication electronic board and the actuation system. In our work, 4 different tests were implemented: 2 feasibility tests for comparison between a common PID controller and the proposed trajectory tracking controller; 2 robustness analysis tests regarding rejecting temporary and permanent disturbances respectively. Similar to the FEM simulations shown in Fig. 4a, the experimental coordinate system of the position sensor F s is set as F s = F I (q e ). As shown in Fig. 5, O represents the origin of the coordinate system.

Trajectory generation

Following the two steps described in Section 3.2, we chose a sinusoidal-shape path located on the top of the external envelop surface of sampled workspace W s E and added a pre-designed time constraint into the chosen path, which yields the planned trajectory shown in Fig. 6a. Then, the practical tracked trajectory can be found accordingly (see Fig. 6b). This trajectory was used as the targeted trajectory both for the feasibility tests and robustness analysis. The reasons mainly lie in two parts: firstly, it is a typically 3D trajectory; secondly, it contains acceleration and deceleration processes, which provides good indices to evaluate the performance of the proposed controller.

Feasibility tests 5.2.1. Trajectory tracking using a common PID controller

In this test, the trunk robot was expected to track the planned sinusoidal shape trajectory as shown in Fig. 6b using a common PID controller. To develop the PID controller, we linearize the studied system around the initial position to determine the gains of the PID controller. Precisely, for (5), we linearize CM(q) *1 H(q) T around its initial configuration q r (t 0 ), noted as À R 3ù4 , then the designed PID takes the following form:

u = *1 (K d Ü
e q e + K p e q e + K i t 0 e q e (s)ds) + u r (t 0 )

where e q e = q r e * q e , K d , K p , K i À R 3ù3 . Therefore the closed-loop system is ì e q e + K d á e q e + K p Ü e q e + K i e q e = Ü Ñ ↵(t)

where Ñ ↵(t) = á q r e + ↵(t) * ⌅ CM *1 H(q) T * ⇧ u. Finally the gains of PID (i.e., K d , K p and K i ) are determined by guaranteeing that ì e q e + K d á e q e K p Ü e q e + K i e q e = 0 is stable. According to Fig. 7, it can be easily noticed that the robot end-effector was not able to catch up with the trajectory and ended up with a huge amount of displacement delay. The results imply that a common PID controller is not valid in terms of tracking a relatively quick and complex trajectory in this case. Theoretically speaking, the PID controller only works locally around q r (t 0 ) if the linearization matrix is close to CM *1 H(q) T , and will fail if their values are far enough. Similarly, the PID controller will also fail if the value of [ á q r e + ↵(t)] varies fast. In our case, the mentioned two reasons may both partially contribute to the failure of this test.

Trajectory tracking using the proposed controller

In this test, we used the proposed controller instead to track the same sinusoidal shape trajectory (Fig. 6b). The results are displayed in Fig. 8 which demonstrates that the robot end effector effectively tracked the trajectory with very little hysteresis. In other words, the proposed controller has been proved to be feasible to track a relatively quick and complex trajectory. Likewise, the results have also proved that the unknown term Ñ ↵(t) defined in Eq. 23 can be regarded as a bounded disturbance, which agrees with our theoretical assumption.

Robustness Analysis

To validate the robustness of the proposed controller, the process of rejecting temporary disturbance and permanent disturbance was observed in 2 different tests (see Fig. 9).

Temporary disturbance

In this test, the trunk robot was expected to track the same sinusoidal shape trajectory (Fig. 6b) under four temporary disturbances exerted by a pen as shown in Fig. 9a. The corresponding results are displayed in Fig. 10 where the four disturbances are marked as 1 -4: the robot end-effector was able to quickly return to and stably catch up with the trajectory even after being affected by four relatively large temporary disturbances (see Fig. 10a, Fig. 10b and Fig. 10c), which implies that the proposed controller is valid for rejecting vanishing bounded disturbances. This experimental phenomenon explains why we can consider d(t) defined in Eq. 24 as a vanishing disturbance, and this term can be logically eliminated due to the existence of the integral term in Eq. 24.

Permanent disturbance

In this test, the robot was planned to track the sinusoidal shape trajectory (Fig. 6b) while carrying a heavy piece of red rubber from the beginning of the trajectory-tracking process. In particular, the extra piece of material on the robot caused the displacement delay at the beginning but the controller rapidly dealt with this permanent disturbance quickly made up for the hysteresis in a very short time (see Fig. 11a, Fig. 11b and Fig. 11c). Therefore, the proposed controller has demonstrated reliable robustness with respect to bounded permanent disturbances. In this case, the piece of rubber can be considered as a permanent constant disturbance, which essentially agrees with the definition of d(t) defined in Eq. 24. Logically, this term can be eliminated by the integral term defined in Eq. 24 if d(t) is constant or piecewise constant.

CONCLUSIONS

Soft robots have been attracting more and more research attention due to their unique properties that are different from conventional rigid robots. Their compliance and flexibility are considered as promising features for the current robotic systems. However, many scientific challenges regarding dynamically controlling these novel robots still exist mainly because of their inherently complex nature. In this paper, we proposed a dynamic controller for a soft trunk robot to track trajectories. The feasibility has been validated through experimental testing, and the robustness of the proposed method has been illustrated by rejecting both the temporary and permanent external disturbances. Note that the axial elongation of the robot body can not be achieved by the current set up since the needed actuators (such as pneumatic actuators) to stretch the robot have not been integrated in this robotic system. However, this will be further discussed and investigated in our future work. Besides, we will also promote this control strategy on other soft robots as well as optimizing it for better control performance.

Figure 1 :

 1 Figure 1:The studied trunk robot[START_REF] Wu | Fem-based gain-scheduling control of a soft trunk robot[END_REF]

Figure 2 :

 2 Figure 2: Mesh elements of the modeled trunk robot.

Figure 3 :

 3 Figure 3: The inertial frames defined for the studied robot.

Figure 4 :

 4 Figure 4: Workspace W E and path planning.

Figure 5 :

 5 Figure 5: Experimental settings

Figure 6 :

 6 Figure 6: The designed trajectory

Figure 7 :

 7 Figure 7: Trajectory tracking using a common PID controller

Figure 8 :

 8 Figure 8: Trajectory tracking using the proposed controller

Figure 9 :

 9 Figure 9: Experimental settings

Figure 10 :

 10 Figure 10: Temporary disturbance test using the proposed controller

Figure 11 :

 11 Figure 11: Permanent disturbance test using the proposed controller


Mohamed Anis Ghaoui
email: mohamed-anis.ghaoui@universite-paris-saclay.fr

Bastien Vincke

Roger Reynaud

Weinberg Madgwick

Human Motion Likelihood Representation Map-Aided PDR Particle Filter

Keywords: Indoor Localization, Pedestrian Dead Reckoning, Particle Filter, Floor Map, Human Motion, Likelihood Grid. Fusion Primary sources IMU Map Prediction

come

L'archive ouverte

Introduction

Localization is one of the most desired services that has been developing over the last decades due to the increasing demand of Location-Based Services [START_REF] Zafari | A survey of indoor localization systems and technologies[END_REF][START_REF] Gu | Indoor localization improved by spatial context-a survey[END_REF]. Outdoor localization has reached a satisfying level of maturity and consumer-ready solutions exist, thanks to satellite-based solutions like Global Navigation Satellite System. Indoor localization techniques see efforts to enhance their capabilities in terms of accuracy, uncertainty, and price. Since satellite signals cannot be received in indoor environments, other ways to locate a device were sought. There are numerous techniques that deploy an infrastructure using various technologies such as WiFi [START_REF] Caso | Performance comparison of wifi and uwb fingerprinting indoor positioning systems[END_REF], Bluetooth [START_REF] Faragher | Location fingerprinting with bluetooth low energy beacons[END_REF], UWB [START_REF] Caso | Performance comparison of wifi and uwb fingerprinting indoor positioning systems[END_REF], RFID [START_REF] Motroni | Sar-based indoor localization of uhf-rfid tags via mobile robot[END_REF] or vision-based [START_REF] Morar | A comprehensive survey of indoor localization methods based on computer vision[END_REF]. These techniques can have prohibitive deployment and energy costs and are often considered intrusive [START_REF] Yassin | Recent advances in indoor localization: A survey on theoretical approaches and applications[END_REF]. Others worked on self-contained techniques that use embedded sensors when using wearable embedded sensors [START_REF] Hou | Pedestrian dead reckoning with wearable sensors: A systematic review[END_REF].

A phone is the typical consumer device that is ubiquitous and always carries an Inertial Measurement Unit (IMU) which enables a self-contained navigation technique: Inertial Navigation System (INS). INSs are attractive because IMUs are reliable, cheap, omnipresent, have a high data rate and low power.

Pedestrian Dead-Reckoning (PDR) seems a promising technique. It leverages accelerometer, gyroscope, and magnetometer data to reconstruct the displacement of each step of a pedestrian. However, in an indoor environment, the magnetometer data is often considered partially or totally unusable [START_REF] Hu | Improving the heading accuracy in indoor pedestrian navigation based on a decision tree and kalman filter[END_REF], making the PDR a relative positioning system to its starting position.

The Particle Filter (PF) models the dynamics of the localization problem by representing the hypotheses space with individual particles and by adding randomness during particles' propagation. This randomness is related to the IMU sensors noises. The PF is particularly interesting in the case of indoor localization because of the ability to fuse a new type of information [START_REF] Tian | A low-cost ins and uwb fusion pedestrian tracking system[END_REF]: the indoor map. Indoor maps are mathematical objects that represent some geometrical and/or semantic information to describe the indoor environment. They can be represented in various forms and levels of details such as floor maps, grids, and spatial graphs [START_REF] Chen | Indoor cartography[END_REF]. In this paper, we pick up on the works of [START_REF] Lu | Indoor positioning system based on chest-mounted imu[END_REF], analyze it and propose an enhancement. We are particularly interested in leveraging the human motion likelihood extracted from a map in the form of a grid. We gradually explain the steps that led to the final result and the impact of this behavior described by the grid. This grid is different from an occupancy grid used in probabilistic robotic localization. Our grid is a mapping of the likelihood of a step performed by the human where each particle of the PF will draw its weight from. While occupancy grids aim to represent the environment by an evenly spaced field random variable and to indicate the probability of a cell being occupied. This paper is organized as follows: Section II presents the main contributions of this paper. Section III summarizes the state of the art related to our problem and discusses the works of [START_REF] Lu | Indoor positioning system based on chest-mounted imu[END_REF]. Section IV sets the mathematical foundations of the problem. Section V is a description of the system. Section VI presents the experiments. Finally, Section VII analyzes the obtained results and provides a discussion surrounding them.

Contributions

This paper presents the following contributions:

• An extension of [START_REF] Lu | Indoor positioning system based on chest-mounted imu[END_REF] via a robust particle filtering based on floor map usage, and on a particle weighting via a human motion likelihood representation of the indoor environment.

• An evaluation method of the proposed system in terms of accuracy, robustness, and integrity.

• A dataset with IMU (accelerometers and gyroscopes), landmarks data, and the floor map of the building.

• A Free Open Source implementation of the framework. 1 3 Related works

Embedded Sensors Indoor localization

Murata [START_REF] Murata | Smartphone-based indoor localization for blind navigation across building complexes[END_REF] identified six main challenges related to indoor localization which are often overlooked. We noted that researchers worked on each point with different approaches. Localization needs to be continuous and accurate in order to offer reliable guidance. Their preliminary tests consider that two meters accuracy is the limit for their indoor application. Probabilistic approaches are designed to consider the sensors' noise but often incorporate correction techniques that may produce a jump in position estimation which is highly undesirable [START_REF] Murata | Smartphone-based localization for blind navigation in building-scale indoor environments[END_REF].

The techniques relying on understanding human motion [START_REF] Baker | The Conventional Gait Model : The Success and Limitations[END_REF] need to be flexible in order to adapt to the individual's distinctive walking pattern. Such adaptivity is far from reach, and thus motion models are approximated. It is often mentioned that humans have multiple motion patterns among: forward, backward, sidestep, turn, ... and each of them presents specific sensor signatures. That is why human motion recognition is problematic in itself [START_REF] Zhang | Human motion recognition based on svm in vr art media interaction environment[END_REF][START_REF] Lin | Deep heading estimation for pedestrian dead reckoning[END_REF]. In this paper, the user is considered to be a forward walker.

Pedestrian Dead-Reckoning and Particle Filter

Pedestrian Dead-Reckoning (PDR) [START_REF] Hsu | Indoor localization and navigation using smartphone sensory data[END_REF] is a technique that combines data from gyroscopes, accelerometers, and magnetometers to reconstruct a human walking motion. There are two main variations [START_REF] Deng | Heading estimation fusing inertial sensors and landmarks for indoor navigation using a smartphone in the pocket[END_REF]: the strap-down approach and the step-and-heading approach. The latter interests us because it sequentially predicts the current pedestrian position by computing the displacement relative to the previous position. This is done by predicting PDR's three main components: step detection, step length, body heading. [START_REF] Khedr | A smartphone step counter using imu and magnetometer for navigation and health monitoring applications[END_REF] deployed an adaptative step detection and counting using accelerometers and gyroscopes signals. By placing an adaptive low-pass filter and a peaks/valleys detection in the signals. All thresholds used in their work are computed adaptively.

Both [START_REF] Lu | Indoor positioning system based on chest-mounted imu[END_REF] and [START_REF] Zhao | Smartphone-based 3d indoor pedestrian positioning through multi-modal data fusion[END_REF] built a PDR-PF system leveraging IMUs data: accelerators and gyroscopes, and additionally the magnetometer for [START_REF] Zhao | Smartphone-based 3d indoor pedestrian positioning through multi-modal data fusion[END_REF]. Even though their prediction functions differ, they produce the same observation data: step length and step heading. To build on the work of [START_REF] Lu | Indoor positioning system based on chest-mounted imu[END_REF], we will need to use the same or equivalent prediction functions presented in section 5.1. For instance, both Lu [START_REF] Lu | Indoor positioning system based on chest-mounted imu[END_REF] and Zhao [START_REF] Zhao | Smartphone-based 3d indoor pedestrian positioning through multi-modal data fusion[END_REF] used the Madgwick filter [START_REF] Madgwick | Estimation of imu and marg orientation using a gradient descent algorithm[END_REF] for angle prediction, which is known to be fast and robust to sensor error. It also deduces the absolute roll and pitch angles from the vertical acceleration, leaving only the yaw relative to its starting angle if the magnetometer is not used. [START_REF] Lu | Indoor positioning system based on chest-mounted imu[END_REF] used double accelerometer integration for step length prediction. While [START_REF] Zhao | Smartphone-based 3d indoor pedestrian positioning through multi-modal data fusion[END_REF] said: "The human gait is modelled by an inverted pendulum of a knee-less biped on the sagittal plane". [START_REF] Do | Personal dead reckoning using imu mounted on upper torso and inverted pendulum model[END_REF] explains the inverted pendulum model even further which is very similar to the Weinberg approximation [START_REF] Weinberg | Using the adxl202 in pedometer and personal navigation applications[END_REF][START_REF] Jimenez | A comparison of pedestrian dead-reckoning algorithms using a low-cost mems imu[END_REF]. We will use the Madgwick filter for heading prediction without magnetometer data and Weinberg Method for step prediction.

Due to drift phenomena associated with noisy gyroscope data, localization with PDR becomes erroneous. Thus, the addition of a PF on top is often done as seen in the IPIN 2018 [START_REF] Renaudin | Evaluating indoor positioning systems in a shopping mall: The lessons learned from the ipin 2018 competition[END_REF]. This is understandable to leverage the discrete prediction of the PDR associated with the discrete multihypothesis capabilities that the PF performs. Previously cited papers do this to achieve better localization [START_REF] Liang | Indoor mapping and localization for pedestrians using opportunistic sensing with smartphones[END_REF][START_REF] Yu | Map-based indoor pedestrian navigation using an auxiliary particle filter[END_REF][START_REF] Ning | Combining a modified particle filter method and indoor magnetic fingerprint map to assist pedestrian dead reckoning for indoor positioning and navigation[END_REF] or to use of a kind of map [START_REF] Lu | Indoor positioning system based on chest-mounted imu[END_REF][START_REF] Zhao | Smartphone-based 3d indoor pedestrian positioning through multi-modal data fusion[END_REF].

The PF is especially interesting since [START_REF] Lu | Indoor positioning system based on chest-mounted imu[END_REF][START_REF] Zhao | Smartphone-based 3d indoor pedestrian positioning through multi-modal data fusion[END_REF][START_REF] Liang | Indoor mapping and localization for pedestrians using opportunistic sensing with smartphones[END_REF][START_REF] Yu | Map-based indoor pedestrian navigation using an auxiliary particle filter[END_REF][START_REF] Ning | Combining a modified particle filter method and indoor magnetic fingerprint map to assist pedestrian dead reckoning for indoor positioning and navigation[END_REF] all differed in the number of particles used for their PF, respectively 256, 500, 1000, 1000, 1000 particles. This indicates that, for what appears to be the same problem, either not every approach is very efficient with its particles use, or the problem dealt with has a hidden layer of complexity.

Researchers have tried machine learning by feeding IMU data to their models and by training them to perform a robust double integral. In their previous work, [START_REF] Yan | Ridi: Robust imu double integration[END_REF] used a linear model to perform a robust integration. But in [START_REF] Yan | Ronin: Robust neural inertial navigation in the wild: Benchmark, evaluations, and new methods[END_REF], they trained a long short-term memory neural network to perform the robust double integration. The aforementioned neural network named RoNIN was trained on the RoNIN dataset, the most significant IMU dataset (40 hours of records). They then compared their results to the well-known Oxford model IONet [START_REF] Chen | Ionet: Learning to cure the curse of drift in inertial odometry[END_REF] trained on the OxIOD IMU dataset [START_REF] Chen | Oxiod: The dataset for deep inertial odometry[END_REF]. We will also compare our system to RoNIN, because for the same inputs, our results are clearly better as shown in Fig. 15.

Map Usage

According to [START_REF] Chen | Indoor cartography[END_REF], maps come in different types, embedded pieces of information, and levels of detail: floor-plan, Building Information Modelling, point cloud, textured mesh, and even hand-drawn. [START_REF] Luperto | Robot exploration of indoor environments using incomplete and inaccurate prior knowledge[END_REF] and [START_REF] Walker | A framework for multi-agent uav exploration and targetfinding in gps-denied and partially observable environments[END_REF] used incomplete or inaccurate prior knowledge, i.e. maps, to obtain better results in their applications. [START_REF] Yu | Map-based indoor pedestrian navigation using an auxiliary particle filter[END_REF] used a floor plan map to eliminate candidate particles that land behind an impassable wall.

[36] considered the indoor environment as a graph where every corridor is an edge and intersection is a vertex. Then, they tracked the gyroscope for 90 • rotation signature and updated their position in the graph by observing the walked distance.

For estimation, we will stick to the usage of a map-assisted PF. Both [START_REF] Lu | Indoor positioning system based on chest-mounted imu[END_REF][START_REF] Zhao | Smartphone-based 3d indoor pedestrian positioning through multi-modal data fusion[END_REF] used the floor map of the building where the experiment was conducted. They applied particles wallfiltering via the floor map where an updated particle that goes into an inaccessible area would be eliminated; else, it survives. The major flaw in their systems is that their particles weights are binary. Thus, they only use the PF as a multi-hypotheses propagator rather than a complete probabilistic estimation filter. This is sensible since the wallfiltering is already efficient at distinguishing particles that perform an impossible physical movement. However, this also means that anything that is not an impossible physical location is plausible. Which, in our view, is flawed. We believe that humans in an indoor space tend to move according to a motion likelihood related to the indoor environment configuration. We chose to represent this behavior as a grid in section 5.2.2.

4 Problem statement and Methodology

Mathematical Model

Both [START_REF] Lu | Indoor positioning system based on chest-mounted imu[END_REF] and [START_REF] Zhao | Smartphone-based 3d indoor pedestrian positioning through multi-modal data fusion[END_REF] processed the floor transition via barometer tracking. In this paper, we will consider working only on a single floor. This means locating on a 2D plan the phone that a pedestrian is carrying in terms of position (x, y) and heading θ. Thus, the device is represented by the state vector P t at each step of timestamp t:

P t = [x t y t θ t ] T | P t ∈ R 3 (1)
A data collection system on the phone is needed to be able to acquire the accelerometer and gyroscope data. From this data, we produce an observation vector U t described in section 5.1:

U t = [dx t dy t dθ t ] T | U t ∈ R 3 (2)
Where dx t , dy t are the displacements along the x and y axes, respectively, and dθ t is the change in heading between two consecutive steps. Then, the displacement of the pedestrian is deduced by observing the accelerometer and gyroscope data at each step via the following equations:

dx t = SL t × cos(θ t ) dy t = SL t × sin(θ t ) (3)
Where SL t is the step length at step index t. θ and SL are obtained via the predictions explained in Section 5.1. The result is then accumulated over the previous position to compute the new position:

P t = P t-1 + U t (4)   x t y t θ t   =   x t-1 y t-1 θ t-1   +   dx t dy t dθ t   (5)
with

x 0 = X 0 , y 0 = Y 0 , θ 0 = Θ 0 as initial conditions.
The predictions of the PDR are then sequentially fed to the PF to estimate the next position. The PF uses multiple particles representing possible trajectories that the pedestrian could have stepped into [START_REF] Jouin | Particle filterbased prognostics: Review, discussion and perspectives[END_REF]. Each particle is represented by a state vector and weight: P i t , w i t (ith particle at t time index). The PF adds a noise vector N t i to both the displacement (dx i t , dy i t ) and differential heading dθ i t . Equation (5) becomes:

P i t = P i t-1 + U i t + N i t ( 6
)
with N i t = δ i t cos(θ i t ) δ i t sin(θ i t ) α i t T (7)
Where δ i t and α i t are the independent random variables respectively drawn from Gaussian distributions N δ and N α with parameters µ δ , σ δ and µ α , σ α defining them respectively (their empirical values are given in section 7.1). Replacing equations ( 3) and ( 5) in equation ( 6) develops into equation [START_REF] Hou | Pedestrian dead reckoning with wearable sensors: A systematic review[END_REF].

         x i t y i t θ i t   =   x i t-1 y i t-1 θ i t-1   +   SL i t × cos(θ i t ) SL i t × sin(θ i t ) dθ i t   +   δ i t × cos(θ i t ) δ i t × sin(θ i t ) α i t   w i t = w i t-1 × L(x i t , y i t ) (8)
A weight is applied to each particle's state vector. The weight allows to filter some particles according to their likelihood L which is a function of the particle's location is:

(x i t , y i t )
. We compute the particles' weighed average Pt as specified in equation ( 9) to choose the pedestrian's new position. This presents an edge case where the Pt may land in a forbidden area.

Pt = xt = N i w i t × x i t , ȳt = N i w i t × y i t , θt = N i w i t × θ i t T (9)
In the Fig. 1, we noticed that Pt (green particle) might land inside/behind the wall, which is physically impossible. To counter this edge case, after computing Pt , we choose the closest valid particle among all particles. This ensures that Pt is valid without having to perform a domain validity analysis. To our experience, this method does not degrade the estimation's performance and is deemed necessary. This will produce the successive positions that a pedestrian walked for each step. It effectively means that our system is Event Triggered by the steps of the user. Once the positions are obtained, we can compute metrics to evaluate the impact of different parameters.

Coordinate Systems

We define three coordinate systems: Global (GCS), User (UCS), and Device (DCS) according to Fig. 2. We make assumptions about the pedestrian verticality to be the same as the building they are in (z g and z u are the same). So, only a translation T u→g and rotation R u→g around the heading θ are needed to express the change of system. The DCS is presented relative to the user carrying the device. We suppose the user to hold the device horizontally as explained in section 6. This implies that we suppose both ψ and φ meaning that only a translation T d→u (corresponding the user's arm) and rotation R d→u around α (gyroscope zaxis) are needed to express the change of system.

System description

The system is split into two main parts: the prediction and estimation functions. The former are a set of tools that transform the raw synchronized accelerometer and gyroscope data into step events t and observation vector U t . The estimation functions are the set of tools that leverage the prediction result by deploying probabilistic mechanisms. The proposed system is offline processing. The processing time of our system is elaborated in section 7.5.

Prediction

Step detection

Step detection techniques have been studied and tested [START_REF] Khedr | A smartphone step counter using imu and magnetometer for navigation and health monitoring applications[END_REF]. Zero-crossing peak valley detection is a technique where the centered acceleration norm is computed from the 3-axis accelerometer A norm . Then, it is filtered via a low-pass secondorder Butterworth filter with a cut-off frequency of 2 Hz. Then, it finds the moments where A norm crosses 0, respectively, from negative/positive to positive/negative, noted t + i , t - i where i ∈ N. Then, find the peaks/valleys in each interval such as:

t peak = argmax t A norm (t) | t ∈ [t + i , t - i ] (10)
t valley = argmin t A norm (t) | t ∈ [t - i , t + i+1 ] (11)
Since each step is defined as a peak then a valley, every pair (t peak , t valley ) defines a step interval. Finally, only peaks/valleys above/below a manually fixed threshold are kept.

Step length

A well-known method is Weinberg's [START_REF] Weinberg | Using the adxl202 in pedometer and personal navigation applications[END_REF] which considers the step length SL to be proportional to a function of the difference between the peak and valley accelerations:

SL = K × 4 a peak -a valley | K ∈ R + ( 12
)
Where K is the Weinberg gain proportional to the length of the pedestrian's leg.

Heading

The heading is predicted by the Madgwick filter [START_REF] Madgwick | Estimation of imu and marg orientation using a gradient descent algorithm[END_REF] that uses gyroscope and accelerometer data to estimate the attitude's quaternion at each sample. We convert the quaternions to Euler angles. Moreover, using the coordinate system changes explained in section 4.2, we only keep the heading value, θ.

The Madgwick filter returns a continuous array of angles. So, we subsample it by only picking the steps events from the continuous heading prediction. Then, it is differentiated to obtain the change in heading dθ t .

dθ t = θ t -θ t-1 (13)
It is important to note that we consider choosing the heading of a step closer to when the foot landed (closer to t valley ).

Estimation

Particle Filter

A PF mainly consists of 4 processes [START_REF] Jouin | Particle filterbased prognostics: Review, discussion and perspectives[END_REF] : Initialization, Prediction, Update and Resampling. We will see and discuss these below.

Initialization We select a starting region by randomly normally generating particles around the starting position and a heading that is the starting heading with some noise added to it. So, there is already some incertitude around the pedestrian's starting position for robustness's sake. Each particle weighs 1 /N.

Prediction The PF predicts the state of particles every time a step is detected via the state model from equation [START_REF] Hou | Pedestrian dead reckoning with wearable sensors: A systematic review[END_REF]. It receives two informations from the prediction stage (Section 5.1), the step length SL and the differential heading dθ t . Every particle computes its new heading and position with angle and radius additive Gaussian noises according to equation [START_REF] Hou | Pedestrian dead reckoning with wearable sensors: A systematic review[END_REF].

Update

In the case of indoor localization, providing a floor map allows the introduction of filtering as [START_REF] Lu | Indoor positioning system based on chest-mounted imu[END_REF] and [START_REF] Yu | Map-based indoor pedestrian navigation using an auxiliary particle filter[END_REF] did. The updated particles are declared valid if they do not cross a wall. This is a straightforward constraint rule. Nevertheless, it is the most accurate mapping of the physical behavior of the pedestrian by the PF. It is often referred to as Map Aiding. The number of particles after this step is lower than the optimal number that the PF is designed for. A floor map of the experiment building is shown in Fig. 3. Resampling It removes lighter-weight particles and duplicates the heavier-weight ones to avoid what is known as PF degeneracy. This process is deemed of the utmost importance [START_REF] Neil J Gordon | Novel approach to nonlinear/non-gaussian bayesian state estimation[END_REF][START_REF] Bolic | Resampling algorithms and architectures for distributed particle filters[END_REF] and has been well studied [START_REF] Bolic | Resampling algorithms and architectures for distributed particle filters[END_REF][START_REF] Elvira | Adapting the number of particles in sequential monte carlo methods through an online scheme for convergence assessment[END_REF][START_REF] Hol | On resampling algorithms for particle filters[END_REF][START_REF] Murray | Parallel resampling in the particle filter[END_REF]. On the other hand, resampling too often presents the risk of obstructing the hypotheses space because the heaviest particles will cluster together and will shadow every other particle (a blackhole-like effect). This is why Kong [START_REF] Kong | Sequential imputations and bayesian missing data problems[END_REF] and Pham [START_REF] Pham | A kalman-particle kernel filter and its application to terrain navigation[END_REF] each introduced a criterion to know when resampling is needed. After implementing and testing different resampling algorithms and criteria, we noticed that the behavior of PF was not enhanced nor degraded due to already using Filtering and Regeneration. We thus abandon resampling in favor of the following one.

Regeneration Although it was not formally named in [START_REF] Lu | Indoor positioning system based on chest-mounted imu[END_REF] and [START_REF] Yu | Map-based indoor pedestrian navigation using an auxiliary particle filter[END_REF], it consists of regenerating the missing particles to have the entirety of particles best represent the updated hypotheses space. The optimality zone of the PF is sensitive to the particles count that is deemed necessary for the PF to be able to correctly represent the hypotheses space (see section 7.1). Since Filtering deleted some particles, the PF is left with Ñ particles, such as Ñ ≤ N . We need to regenerate the N -Ñ missing particles. The way of regenerating these missing particles affects how the PF paves the hypotheses space and what it "believes" to be the more likely position. If all particles had the same weight, the particles would regenerate according to a majority vote. But, we propose in section 5.2.2 a method to weigh our particles and to introduce a priori likelihood.

The association of Map Filtering and Regeneration emulates a resampling where the particles landing in inaccessible zones (e.g., behind walls) are eliminated, and only the valid particles will have offspring.

Subset Regeneration: We studied more formally the approach [START_REF] Lu | Indoor positioning system based on chest-mounted imu[END_REF] implemented in their open-source code. In order to reach N particles, it considers that the new M particles need to be generated from M randomly selected subsets of N particles defined as N = ss × Ñ where ss ∈]0, 1] is the subset proportional size. Fig. 4 shows a visualization of algorithm 1.

Our contribution is to randomly draw new particles from // check if the candidates do not cross walls 8 valid = filter(candidates, mean particles) 9 return valid the statistical properties (mean and variance) of the subset and to verify that the new proposed particles do not cross a wall of the map. This acts as a random re-sampling based on the paved hypotheses space where each randomly selected subset is a pavement from where the new particle is picked. It emulates a re-sampling algorithm but bases itself on both statistical properties of the particles and physical properties of the indoor environment. After our tests, we empirically selected ss = 5% and preferably Ñ × ss ≤ 25.

Human Motion Likelihood grid

Our final contribution is the creation of a grid that represents the human motion likelihood. To represent the human motion inside a building on a floor map, we set an a priori about the likelihood during motion associated with the different indoor dispositions and configurations. We suppose that a human is more likely to move in the center of a corridor rather than hugging the wall, walk through the center of a door rather than the sides.

We convert the floor map of length L and width W to a grid: G ∈ R l×w of length l = L/scale and width w = W/scale, with each cell g ∈ G | g ∈ [0, 1]. In our case, scale = 5 cm. So the grid is composed of 5 × 5 cm 2 cells. Then, we compute the distance from each cell to the closest wall.

The resulting distances map is converted into a likelihood grid via a membership function f defined as f : N l×w → [0, 1] l×w . The choice of the mapping functions f , shown in Fig. 5 as an example, affects how the grid constrains the particles during their propagation. In every case, the function f creates a centering constraint as it pushes the particles into the center of corridor. It has 2 parameters: min range, that determines the end of the forbidden area; and max range, that determines the end of the transition area. The resulting grid is shown in Fig. 6 and represents our a priori on human motion likelihood. A zoomed n example of a corridor grid is given in Fig. 7a. Black cells are walls with an associated likelihood of 0. Fairer cells have a higher likelihood of a step occurring on them and thus associate a more important weight to the particles landing on them.

Setup and experiment

In order to test different configurations of the system's parameters, we need to set up an experiment that reflects a realistic human motion in indoor space. We use a OnePlus 6T phone to acquire LSM6DS3 IMU data (accelerometer + gyroscope. The summary of this IMU's characteristics 2 in Table 1 shows that the IMU uses 16 signed bits to code its data. The Linear acceleration measurement range is selected because we noticed that the accelerometer data never surpassed 1.5g. Also, human rotation didn't exceed 180 The user holds the phone horizontally, the USB port towards their chest, and the back camera looks at the ground. This ensures that the coordinate system does not change during the experiment and considers Android's Z rotation axis as the pedestrian heading θ.

To produce relevant metrics to measure the performance, we need to set up a pedestrian ground truth. This ground truth must have a satisfying spatial and temporal accuracy and uncertainty. So, we put landmarks in the form of crosses on the ground and ask the pedestrian to step on them and simultaneously press a timestamp button on the phone. This procedure is illustrated in Fig. 7b. We use a high-precision laser telemeter to locate the crosses relative to the building's walls. If synchronization is done correctly, the moment the user presses the button should correspond to when they step on the cross. This reduces the ground truth effort into a synchronization problem. Experiments showed that this evaluation method was satisfactory as we can count how many timestamps the user produced and invalidate the run if any are missing. We also have a worst case error of a step around the cross compared to when the user pressed the timestamp button (the user pressed the button a step before, while on, or after the cross).

While looking for an acquisition tool, we settled for the INRIA-Tyrex team's open-source app: SensLogs [START_REF] Michel | Attitude Estimation for Indoor Navigation and Augmented Reality with Smartphones[END_REF]. The app allows for a plethora of sensor acquisition and different sampling frequencies. We select the maximum frequency our phone allows: 400 Hz. The temporal synchronization is guaranteed because all the sensors are timestamped with the same clock source. Thanks to this, the tedious work of synchronization is abstracted, and we have an all-ready solution.

We define a trajectory to have various indoor cases as shown in Fig. 3. The trajectories ought not to loop around a single place many times as this biases when performing map aiding. This happens if the first loop is correctly performed, then the others are quasi-automatically correct, which makes a system looks more permanent than it is.

We previously mentioned that our system had many parameters. However, two of them are particular because they are specifically related to the user holding the phone: the Weinberg gain from section 5.1.2 and the acceleration threshold from section 5.1.1. We tested our system's invariance to its user by having four different users performing the same trajectory. Then, we fit the aforementioned two parameters to each user.

To measure our system's performance at each user-made timestamp, we compute the Euclidean distance from the average weighed position Pt to the landmark l t for landmarks count, K:

E = K t=1 E t K | E t = l t -Pt 2 ( 14
)
And the Mahalanobis distance of the landmark l t to the N particles' positions P t ∈ R N ×2 with weights W t ∈ R N and weighed average Pt presented in equation ( 9).

M = K t=1 M t K | M t = (l t -Pt ) × S -1 t × (l t -Pt ) T (15)
With S ∈ R 2×2 being the covariance matrix of the particles where each particle has two components: P t = (x t , y t ):

S t = E W T t P t -Pt × W T t P t -Pt T (16)
It is important to notice that the weighted average position Pt differs from the one presented in section 4.1. It is not the closest valid particle but the statistically weighted average of the particles. This makes our measure take into account the worst-case scenario for the PF.

The parametric study must be done over several runs over the same data due to the random nature of the Monte Carlo sequence. Then, these runs are averaged to obtain the curves shown below. We performed an average over ten runs. The success rate of a trajectory is defined as completed runs /total runs with a run considered completed if it has as many user timestamps as landmarks there are, and a return to initial position error defined as E f inal ≤ 2.5m inspired by [START_REF] Zhao | Smartphone-based 3d indoor pedestrian positioning through multi-modal data fusion[END_REF]. An example of a run is given in Fig. 8.

Results, analysis, and discussions

We are interested in showing how the integrity of the system evolves and the impact of the human motion likelihood grid that we conceptualized in section 5.2.2. The available data and code can recreate the results and cover how to obtain them in more details. The first step is to determine Weinberg's gain and acceleration threshold experimentally. This is simply done by making users walk a measured distance, count their steps, and then empirically reverse engineer each parameter for each user.

Optimal particles count

In order to discuss results, we first need to determine how many particles our system requires to model the hypotheses space correctly. For that, we proceed to compute both average E and M for the same trajectory, same data, and same user for a given number of particles: N . The Fig. 9 shows that the system converges for at least 200 particles. Counts beyond this number do not significantly change either metric but increase the necessary computing resources. We will be using 200 particles for the remaining experiments. We also give the particles' propagation noise used in equation (7) in table 2.

Weinberg's functional interval

Starting from the empirical Weinberg's gain value from section 7, we figure out the optimal value for a user. The Fig. 10 shows that both E and M (red and yellow curves) reach a minimum value close to the empirical Weinberg gain's value. Since this value is never supposed to change for a user, we make them walk the exact trajectory but in reverse order to obtain the blue and cyan curves. We also compute the success rate of each trajectory to understand the impact of this parameter. We first notice that the two trajectories do not have optimal Weinberg gain value. Nevertheless, an interval exists where the success rates overlap and keep the errors to a minimal value. Table 3 shows the comparison between the empirically found and statistically computed Weinberg gain's value. This implies that the step noise δ introduced in equation ( 6) corrects the errors in step length predictions.

Grid effect and integrity

By measuring the signed error in the position of all particles relative to a landmark: e l and its standard deviation σ e l , we can observe the effect of the grid on the distribution of particles throughout the whole experiment. The errors in X and Y, e x l and e yx are respectively surrounded by ±3σ ex and ±3σ ey . We compare the performance of the system for the same data and parameters with a scenario where the likelihood grid is enabled and disabled on an average of 10 runs each time. This is equivalent to a comparison between [START_REF] Lu | Indoor positioning system based on chest-mounted imu[END_REF] and our system. In order to do this comparison, we need to find the optimal values of the two parameters min range and max range from section 5.2.2. To find the optimal values of both parameters, we perform a parametric study on the grid. The goal is to maximise the score of the grid version vs no-grid version, defined as:

score = 1 K K t=1
σ e l t no-grid -σ e l t grid [START_REF] Lin | Deep heading estimation for pedestrian dead reckoning[END_REF] Equation ( 17) is analogous to computing the area difference of incertitude ellipsoids at each key-point between the no grid and grid version in Fig. 14 and13. In addition, equation ( 17) allows a balance of the centering constraint from section 5.2.2. We make sure that the system keeps its integrity by ensuring that the ground truth is still bounded by the error i.e. that the 0 y-axis is included in the purple corridor for Fig. 11. This is expressed as:

score t = 0 if (e t + σ t ) × (e t -σ t ) ≥ 0 score t else ; (18)
We compute the score for each combination of the two parameters according to these two equations. We obtain the heat map presented in Fig. 12. We first notice that, if min > max, the system cannot work because the forbidden area includes the allowed area, and all the particles receive a weight of 0. Secondly, the blue area shows that the score was at its maximum for those corresponding values of parameters. Thirdly, We notice that beyond a min value of 14 cells (70 cm), the forbidden area from a wall is too big and cuts the path for the particles. This is explainable because, in our path, some corridors are 150 cm wide. This means that we would have a forbidden area of 70 cm on each side with an error due to the discretization of the map into a grid. We will pick: min range = 11 and max range = 13 (one of the darkest blue cells in the Fig. 12).

In Fig. 11, the 0 on y-axis is always included between ±3σ for both version. This means that in 99.7% of cases, the system includes the ground truth position inside the uncertainty ellipsoid. This expresses its integrity, i.e., how faithfully it can always contain the correct positive in its hypotheses space.

Another aspect to note is the impact of the usage of the grid. We do notice that enabling the grid (with both its parameters, see section 5.2.2) makes the area imprisoned between ±3σ narrower than without the likelihood grid. This means that the likelihood grid, in that configuration, compresses the particles to the center of the corridor, which means that, on average more particles are closer to the correct position (landmark).

If we compare the ellipsoid of the particles' distribution at each landmark for an enabled and disabled grid as shown in Fig. 14. We notice that the ellipsoids are constrained by the walls and rarely go through them when the grid is enabled. This shows that the grid has an ameliorative effect on the localization for these indoor configurations. The colors of the particles indicate their weights. In the case of the disabled grid, all particles have the same weight (color) 1 /N. The system's behavior is interesting at some particular points that we encounter in the building during the walk, namely: 8, 12, and 24. By observing Fig. 8, 11 and 14, we notice that landmark 8 comes after a walk of ≈ 10m in a semi-open space. Thus, non-constrained particles tend to spread and uniformly pave the hypotheses spaces. Landmark 12 shows an interesting compression of the particles' ellipsoid. It comes after a 90 • turn with some semi-open space, yet the grid version achieves an important compression and has no particles in the adjacent rooms like landmark 9 in Fig. 13. Landmark 24 is an edge case as it combines both difficulties of 8 and 12. It comes to a semi-open hyafter a 90 • turn. Particles easily spread to fill the space. Then, they get funneled into the tight corridor. This makes the ellipsoid spread to its maximum and then fold back into the corridor. All of these physical phenomenons are correctly displayed by spikes in error e t and std σ t on both axis X and Y in the Fig. 11.

Discussions

Map Filtering and Aiding

[12] deployed a backtracking procedure to ensure that the particle's life was coherent and did not cross walls over the last 20 steps. This seems redundant to us as our map filtering, presented in section 5.2.1, guarantees by construction that a particle dies if it hits a wall. The comparison to [START_REF] Lu | Indoor positioning system based on chest-mounted imu[END_REF]'s system is made by comparing the presence and absence of our grid in section 7.3. In Fig. 15, we show the comparison of 3 methods: a fully deterministic PDR (no wall filtering), the prediction of the already trained RoNIN [START_REF] Yan | Ronin: Robust neural inertial navigation in the wild: Benchmark, evaluations, and new methods[END_REF] and finally, our system. The resulting trajectory shows that map filtering is mandatory for indoor environments because it maps the physical displacement more faithfully than sheer statistical methods. It is also important to note that RoNIN relies on Android's orientation (Game Rotation Vector) which itself leverages the faulty indoor magnetometer data. To further show the difference between the 3 approaches, Table 4 gives the computed RMSE for each system. We also resolved some edge cases like the weighted average of particles in section 4.1.

System

Human Motion Likelihood Grid

At first, there seemed to be no enhancement in terms of accuracy. However, we showed that the integrity of the system and its ability to map the hypotheses space in the indoor environment more accurately are enhanced. Nevertheless, although the users were told to "walk as they like" and "follow the path of the crosses", they will tend to walk in segments from one landmark to another. Which makes for corridor walks to seem more accurate. However, this isn't true for landmarks 8, 12, 16 and 24 in Fig. 14. We notice that those configurations differ from simple corridors. However, the grid reinforces the system's integrity by preventing the propagation of deceptive particles in the previous steps, making the localization at those difficult spots more reliable.

Nevertheless, this comes at some cost of flexibility. The Fig. 12 showed that there is a balance in how constraining the grid is and how relaxed it is. Thus, the choice of the likelihood function presented in section 5.2.2 has a key role. Our hypothesis remains valid as long as we suppose a single pedestrian is walking alone in a building. It begs how this likelihood evolves if more pedestrians were walking either in the same, opposite direction, or even at a stop.

Processing Time

Our system is offline processing on CPU: Intel i5-9500(6)@4.400GHz. The floor map from Fig. 3 is composed of nearly 300 segments that represent walls. The grid is made of 5x5 cm 2 cells for a total of 1400 × 720 ≈ 10 6 cells (see Fig. 7a), with each cell containing a floating value. This grid is built only once per floor. It is then a read-only memory object.

The particles count was set to be 200 particles in Sec-tion7.1. The implementation is done on Python with an average step processing speed of 4.5 steps/s. With the knowledge that a human can walk between 1 to 3 steps a second, the current offline processing is satisfactory even before optimization.

Optimization can be done in several ways in future studies. The most straightforward one is to change the implementation language from an interpreted one to a compiled one (C++, C#, Java/Kotlin for Android). This change will reduce memory usage and execution time.

Conclusion

This paper presents an indoor localization system based on a phone that the pedestrian carries. We have based our work on [START_REF] Lu | Indoor positioning system based on chest-mounted imu[END_REF]. Then, we have detailed the prediction and estimation functions that we set up. Then, we have introduced the human motion likelihood representation as a grid in the PF. This grid allowed to weigh the particles. We have also defined an experiment and evaluation method with reliable ground truth for the pedestrian's indoor movement. We have proceeded to the analysis of the results of the experiments.

Furthermore, we have given a qualitative appreciation for the metrics we have defined. These metrics showed that the integrity of our system is more important with the enabled grid. We compared our system to a variant of [START_REF] Lu | Indoor positioning system based on chest-mounted imu[END_REF] and to [START_REF] Yan | Ronin: Robust neural inertial navigation in the wild: Benchmark, evaluations, and new methods[END_REF].

In the future, we wish to experiment on multi-agent representation of different behaviors. More semantic elements will be added into the map to enhance indoor localization's integrity and robustness.

Data Transparency and Experiment Reproducibility

We declare that the used data is from real-world subjects and was collected with anonymity and integrity. The building's map is delivered as a numerical object with faithful measured dimensions of a frequently busy physical building. Any sensitive pieces of information that are not required for the experiments have been censored without compromising the integrity of the data or its ground truth.

proach centered on intrinsic autonomy augmented by external assistance. His teaching focuses on instrumentation for physical measurements, signal processing, robotics, and data fusion.

Figure 1 :

 1 Figure 1: Particles around the corner.

Figure 2 :

 2 Figure 2: Coordinate systems changes.

Figure 3 :

 3 Figure 3: A floor-plan of the experiment building: grayed out areas are forbidden. User walks along the landmarks from 0 to 28.

Figure 4 :Algorithm 1 : 6 N

 416 Figure 4: N particles subset regeneration are green.

Figure 5 :

 5 Figure 5: Example of distance to likelihood functions f (d) = L. .

Figure 6 :

 6 Figure 6: Likelihood grid built upon the floor map in Fig. 3: darker areas are less likely.

  (a) Example of corridor grid. (b) Walking over a cross and pressing the timestamp button.

Figure 7 :

 7 Figure 7: Left: grid concept. Right: experiment setup.

Figure 8 :

 8 Figure 8: An estimated trajectory of a run. Red dots are the steps when the user pressed on the timestamp in Fig. 7b.

Figure 9 :

 9 Figure 9: Measured distances as a function of particles count.

Figure 10 :

 10 Figure 10: Measured distances as a function of Weinberg gain.

Figure 11 :

 11 Figure 11: Measured signed errors on each key-points.

Figure 12 :

 12 Figure 12: Score heat map for both grid parameters.

Figure 13 :

 13 Figure 13: Incertitude ellipsoids with likelihood grid disabled.

Figure 14 :

 14 Figure 14: Incertitude ellipsoids with likelihood grid enabled.

Figure 15 :

 15 Figure 15: Comparison between 3 different approaches: Our system, RoNIN [31] and a deterministic PDR on the same trajectory.

Table 1 :

 1 Summary of LSM6DS3 IMU.

	Parameter	values	unit
	Linear acceleration measurement range	± 2	g
	Linear acceleration sensitivity at± 2g	0.061	mg/LSB
	Angular rate measurement range	± 250

• /s. So, an Angular rate measurement range of ±250 • /s seams reasonable. • /s Angular rate sensitivity at ± 250 • /s 8.75 m • /s/LSB

Table 2 :

 2 PF's propagation noise distributions.

Table 3 :

 3 Weinberg gain: empirical vs statistical values.

Table 4 :

 4

		RMSE[m]
	Ours	0.811
	RoNIN	7.534
	deterministic PDR	4.671

RMSE of comparison of Fig.15.

Roger Reynaud is a Pr. in Paris-Saclay University. His research topics deal with: data fusion using probabilistic and belief theories; obstacle detection, tracking, localization; Algorithms Prototyping with "Algorithm Architecture Adequacy" methodology; and Design of quasi-autonomous systems with an ap-