M Aucejo

L Maxit

N Totaro

J.-L Guyader

Convergence acceleration using the residual shape technique when solving structure-acoustic coupling with the Patch Transfer Functions method

Keywords: Strong coupling, PTF, FEM, residual modes 1. Introdution

The forced response of the structure-water-filled cavity system is investigated from the Patch Transfer Functions method. In such a case, a poor convergence of the PTF method is observed when using standard mode expansion to build the cavity-PTF. To improve its convergence and maintain the advantages of substructuring, residual shapes are introduced in the cavity-PTF computation, which is the new material of this article. This technique is successfully applied on numerical examples, highlighting the interest of such an approach, especially in heavy fluid.

The forced response of a structure backed by an acoustic cavity has been studied extensively during the last few decades. Classical models are based on the coupling of in vacuo structural modes and rigid-walled cavity modes, providing comprehensive formulations to solve the coupled system [START_REF] Dowell | Acoustoelasticity: General Theory, Acoustic Modes and Forced Response to Sinusoidal Excitation, including Comparisons with Experiment[END_REF][START_REF] Guy | The Transmission of Sound Through a Cavity-Backed Finite Plate[END_REF][START_REF] Gagliardini | The Use of a Functional Basis to Calculate Acoustic Transmission Between Rooms[END_REF]. These models introduce a modal cross-coupling term to describe the interaction between the structure and the cavity. Sum and Pan [START_REF] Sum | On Acoustic and Structural Modal Cross-Couplings in Plate-Cavity Systems[END_REF] have analyzed this interaction term to underline the physical mechanisms of cross-coupling. When considering a geometrically complex system, analytical methods are often replaced by the Finite Element Method (FEM). However, the direct resolution of the coupled problem is time consuming. To overcome this difficulty, model reduction techniques were developed such as the Component Mode Synthesis (CMS) methods [START_REF] Wc Hurty | Dynamic Analysis of Structural Systems using Component Modes[END_REF][START_REF] Tournour | Validation, Performance, Convergence and Application of Free Interface Component Mode Synthesis[END_REF]. These methods, based on the modal approach, allow describing the dynamics of the coupled system with a reduced but sufficient number of degrees of freedom. In the case of strong coupling, the use of uncoupled modal bases can lead to poor convergence if the coupling between the high order modes of a subsystem with low order modes of the other system is not considered in the computation. To improve the convergence of the modal approach, several techniques were developed, such as pseudostatic correction [START_REF] Hjp Morand And R Ohayon | Fluid-Structure Interaction[END_REF][START_REF] Ohayon | Reduced Models for Fluid-Structure Interaction Problems[END_REF][START_REF] Tournour | Pseudostatic Corrections for the Forced Vibroacoustic Response of a Structure-Cavity System[END_REF] and residual modes [START_REF] Qh Tran | Comparison of Robust CMS Methods for Damped Vibroacoustic Problems[END_REF][START_REF] Roy | Efficient Computation of the Radiated Sound Power of Vibrating Structures using a Modal Approach[END_REF][START_REF] Cuppens | Evaluation of the FRF Based Substructuring and Modal Synthesis Technique Applied to Vehicle FE Data[END_REF]. The aim of these techniques is to obtain simple and efficient reduction bases to maintain the advantages of modal superposition while improving its accuracy. In the present paper, the Patch Transfer Functions (PTF) method is used to solve the structure-cavity coupling, when the acoustic cavity is filled with water. This method developed by Maxit et al. [START_REF] Maxit | Airborne Noise Prediction using Patch Acoustic Impedance[END_REF] is a tool for predicting the vibro-acoustic response of complex systems coupled through surfaces divided into elementary areas called patches. This approach, based on substructuring [START_REF] Kergoulay | Model Reduction for Efficient FEM/BEM Coupling[END_REF][START_REF] Tournour | Validation, Performance, Convergence and Application of Free Interface Component Mode Synthesis[END_REF] and the impedance and mobility approach [START_REF] Kim | A Compact Matrix Formulation using the Impedance and Mobility Approach for the Analysis of Structural-Acoustic Systems[END_REF], consists in studying each subsystem independently in order to build a set of transfer functions defined by using mean values on the patches, called Patch Transfer Functions. Then, assembling PTF by using the continuity relations leads to the fast resolution of the coupled problem. However, as presented in this paper, the use of a standard mode expansion to build the cavity-PTF leads to poor convergence of the PTF method in heavy fluid. To alleviate this problem, we introduce the residual mode technique in the present approach. This novelty allows improving the convergence of the coupled system response and maintain the advantages of substructuring, by computing a residual shape basis for the rigid-walled cavity independently of the structure. Classically, these residual shapes, which are attachment modes, correspond to the response of the cavity subjected to structural mode excitation. Nevertheless, in the PTF method, the residual shapes correspond to the response of the cavity excited by a constant normal displacement imposed on a patch at a specific frequency. This ensures the convergence of the proposed method. In the following sections, the theoretical background of the proposed method is first recalled and the PTF computation is then derived from the FE formulation of each subsystem is developed. In sections 5 and 6, the structural-acoustic problem is considered by comparing the forced response of a plate coupled to a water-filled or air-filled cavity, either obtained by a direct resolution using FEM, or by the PTF method including the residual mode technique. The introduction of the residual shapes concept in cavity-PTF computation allows obtaining good results with a coarser interface mesh than that used for standard FEM. Finally, to illustrate the advantage of the present approach on an industrial case, the radiation of a bulkhead two ballasts is considered. The influence of residual shapes is studied by comparing PTF results obtained either from a direct computation of the PTF of both ballasts or a mode expansion with or without the introduction of residual shapes in the calculations.

Description of the vibro-acoustic problem

Let us consider a structure backed by an acoustic cavity with rigid walls and excited by a point force F on its surface S v as presented in Fig. [START_REF] Dowell | Acoustoelasticity: General Theory, Acoustic Modes and Forced Response to Sinusoidal Excitation, including Comparisons with Experiment[END_REF]. The surface S v can also be identified as the coupling surface of the structure-cavity system.

K s -ω 2 M s {U s } = {F } -[A] {P } (1)
where M s and K s are the mass and stiffness matrices of the structure, A the coupling matrix of the structure-cavity system and ω the radian frequency.

The coupled acoustic cavity problem is governed by the homogeneous Helmholtz equation. Rigid-walled boundary conditions are imposed on the rigid surface S r and a source term

Q s = -ω 2 [A]
T {U s } resulting from the vibration of the structure is imposed on S v . The acoustic pressure in the cavity is harmonic {P (M, t)} = {P (M )} e jωt and satisfies the following equation:

K f -ω 2 M f {P } = -ω 2 [A] T {U s } (2)
where M f and K f are the mass and stiffness matrices of the acoustic cavity.

Eqs. ( 1) and (2) leads to the linear matrix system given by the following equation:

K s A 0 K f -ω 2 M s 0 -A T M f U s P = F 0 (3)
Eq. [START_REF] Gagliardini | The Use of a Functional Basis to Calculate Acoustic Transmission Between Rooms[END_REF] shows that for the general structure-acoustic problem, the system to be solved has a non-symmetric mass and stiffness matrices. This lack of symmetry leads to difficulties in computing the dynamic response of the system, since Eq. ( 3) is costly to solve. Likewise, the extraction of eigenvalues for large structure-acoustic systems can be difficult with this non-symmetric formulation. To alleviate these difficulties, several symmetrization procedures have been developed. These procedures are based on different approaches, namely the displacement and pressure (u, p) formulation or displacement, pressure and vorticity moment (u, p, Λ) formulation [START_REF] Kj Bathe | A Mixed Displacement-Based Finite Element Formulation for Acoustic Fluid-Structure Interaction[END_REF][START_REF] Wang | Displacement/Pressure Based Mixed Finite Element Formulations for Acoustic Fluid-Structure Interaction Problems[END_REF], the velocity potential formulation [START_REF] Gc Everstine | A Symmetric Potential Formulation for Fluid-Structure Interactions[END_REF][START_REF] Olson | Analysis of Fluid-Structure Interactions -A Direct Symmetric Coupled Formulation Based on a Fluid Velocity Potential[END_REF] and the displacement formulation [START_REF] Ma Hamdi | A Displacement Method for the Analysis of Vibrations of Coupled Fluid-Structure Systems[END_REF]. Finally, it should be noted that the matrix system given by Eq. ( 3) is not used in the proposed method, since each subsystem is considered independently and coupled by continuity equations.

Patch Transfer Functions method

Description of the method

The PTF method [START_REF] Maxit | Airborne Noise Prediction using Patch Acoustic Impedance[END_REF][START_REF] Ouisse | Patch Transfer Functions as a Tool to Couple Linear Acoustic Problems[END_REF][START_REF] Totaro | Extension of the Patch Transfer Functions Method (PTF Method) to High Frequency Domain (Sub-Cavities Decomposition)[END_REF] is a substructuring approach that allows computing the vibro-acoustic response of complex systems coupled through their common surfaces divided into elementary areas called patches. Then, the Patch Transfer Functions (PTF), corresponding to an impedance or a mobility are computed for each patch of the coupling surface and assembled by using the superposition principle and writing the continuity relations. This technique is efficient and well adapted to complex systems, since each subsystem is defined by its own FE model and solved independently of the others. Consequently, modification of one of the subsystems only leads to a simple update of this subsystem, provided the geometry of the structure-cavity interface remains unchanged. Furthermore, the PTF method leads to small, symmetric systems fully populated with complex coefficients. Moreover, the impedance and mobility matrices, which are the basic quantities of the method, can be obtained by measurements, analytical modelling and numerical modelling (FEM, BEM, etc.). Consequently, parallel computation of the basic quantities for each subsystem is possible, since the basic quantities are established on uncoupled subsystems. This is one of the method's advantages, since different approaches can be combined to obtain the response of a multi-coupled complex system. Moreover, if the PTFs are derived from the FE formulation, the use of incompatible meshes at the interface of the subsystems is possible, since the problem of compatibility is solved by the patch mesh as illustrated by Fig. [START_REF] Guy | The Transmission of Sound Through a Cavity-Backed Finite Plate[END_REF]. The management is indeed permitted because the PTF method does not relate nodes values but values averaged over the surface of a patch. However, as presented in Fig. [START_REF] Guy | The Transmission of Sound Through a Cavity-Backed Finite Plate[END_REF], the mesh of the subsystems can not be highly incompatible at the interface insofar as nodes have to be included in each patch.

Finally, the PTF method differs from standard substructuring methods, e.g. the Craig-Bampton fixed interface approach [START_REF] Craig | Coupling of Substructures for Dynamic Analysis[END_REF], since the subsystems are coupled through patches instead of nodes. It should also be noted that Ouisse et al. [START_REF] Ouisse | Patch Transfer Functions as a Tool to Couple Linear Acoustic Problems[END_REF] and Chazot and Guyader [START_REF] Chazot | Prediction of Transmission Loss of Double Panels with a Patch-Mobility Method[END_REF] have shown that the number of patches is lower than the number of interface nodes. As a consequence, the size of the PTF model is smaller than that obtained with the Craig-Bampton method. Indeed, the size of the reduced model is N patch × N patch (N patch : Number of patches) for the PTF method, while for a Craig-Bampton method the size of the reduced model is (N im + N idof ) × (N im + N idof ) (N im : Number of fixed interface modes chosen for each subsystem, N idof : Number of interface dofs).

PTF definition

In the case of the structure-cavity problem, the coupled system is divided into two subsystems, namely a structure and an acoustic cavity, and PTF are defined for each subsystem on the coupling surface S v divided into patches.

For the structure, the PTF Y s jk between an excited patch k and a receiving patch j is the mobility defined as the ratio of the mean structural velocity vs j on patch j and the mean pressure ps k on patch k (see Eq. ( 4)),

Y s jk = vs j ps k (4)
where •p = 1 Sp Sp •dS p is the space-averaged quantity on the patch p.

For the acoustic cavity, the PTF Z c jk between an excited patch k and a receiving patch j is the impedance defined as the ratio of the mean acoustic pressure pc j on patch j and the mean normal velocity vc k on patch k (see Eq. ( 5)),

Z c jk = pc j vc k (5)

Subsystems coupling

Coupling the subsystems is performed by writing continuity conditions on each patch of the coupling surface S v as expressed in Eq. ( 6) for a patch j.

ps j = pc j = pj vs j = vc j = vj (6)
These coupling conditions are then used to determine the coupling normal velocity vj on each patch j. For this purpose, the linearity of the problem is used. The normal structural velocity vs j on the patch j is thus defined as the superposition of the in vacuo structural velocity vs j (uncoupled structural response) and the normal velocity Y s jk pk due to the pressure caused by the acoustic cavity on the structure (see Eq. ( 7)),

vs j = vs j + Y s jk pk (7)
Similarly, the pressure pc j on patch j of the cavity can be calculated by using linearity, as the superposition of the pressure created by the acoustic source inside the cavity pc j and the pressure due to structural vibrations Z c jk vk (see Eq. ( 8)).

pc j = pc j + Z c jk vk (8)
By introducing the coupling conditions (Eq. ( 6)), simple mathematical operations make it possible to obtain the coupling normal velocity vj as presented in Eq. ( 9) in the absence of acoustic sources in the cavity.

vj = I -Y s ji Z ik -1 vs k ( 9
)
where I is the identity matrix.

PTF computation from a FE formulation and standard mode expansion

Structure-PTF computation

Defining the structure-PTF requires determining the structural velocity when a constant pressure pk is imposed on patch surface S k . Consequently, Eq. ( 10) has to be solved for each patch k of the vibrating structure.

K s -ω 2 M s {U s } = {F k } (10)
In Eq. ( 10), the excitation term F k corresponding to the constant pressure pk imposed on the patch k is given by the following equation:

F k = A k N pk on S k 0 on S v \ S k (11)
where A k is the area of patch k and N is the number of FE nodes included in patch k.

Eq. ( 11) is solved by standard mode expansion. For this purpose, displacement {U s } is expanded on its structural modes as expressed in the following equation:

{U s } = [ψ n ] {u n } (12)
where [ψ n ] is the mode shape matrix and {u n } is the modal amplitude vector.

The generalized mass [M sn ], stiffness [K sn ] and force [F n ] are defined as presented in the following equation:

     [K sn ] = [ψ n ] T [K s ] [ψ n ] [M sn ] = [ψ n ] T [M s ] [ψ n ] {F n } = [ψ n ] T {F f } (13)
The reduced system presented in Eq. ( 14) is therefore solved at each frequency to obtain structural displacement {U s }.

K sn -ω 2 M sn {u n } = {F n } ( 14
)
This procedure is repeated for each excited patch k leading to the structural velocity matrix [V ] of dimensions N dof × N patch , where N dof is the number of degrees of freedom of the structure and N patch the number of patches.

The structure-PTF Y s jk is thus obtained by averaging the normal velocities for the nodes included in the receiving patch j.

Cavity-PTF computation

The pressure response to a unit normal velocity vc p imposed on the patch surface S p has to be computed for each patch p of the coupling surface of the acoustic cavity, leading to the solution of the following equation:

K f -ω 2 M f {P } = {Q p } (15)
In Eq. ( 15), the excitation term Q p corresponding to the constant normal velocity vc p imposed on the patch p is given by the following equation:

Q p = -jω Ap N vc p on S p 0 S v \ S p (16)
where A p is the area of patch p and N is the number of FE nodes included in patch p.

In this case too, Eq. ( 15) is solved by standard mode expansion, meaning that the acoustic pressure {P } is expanded on its rigid-walled cavity modes, as shown in the following equation:

{P } = [φ n ] {p n } (17)
where [φ n ] is the mode shapes matrix and {p n } is the modal amplitudes vector.

We now use the same resolution procedure as that used to compute the structure-PTF Y s jk to obtain the cavity-PTF Z c jk .

Concluding remarks

The equations developed in this section are general, since they are based on FE models. Consequently, the approach can be applied to complex systems. Furthermore, the standard mode expansion procedure has been used up now to compute the PTF without introducing the residual shapes concept. As demonstrated in the following sections, convergence of the method is difficult to achieve for acoustic cavities filled with heavy fluids such as water. In this case, residual shapes are computed (see section 6) to enrich the original cavity modal basis and improve the convergence of the method.

Structure-cavity coupling

In this section, the PTF method is used to compute the forced response of a structure-cavity system and compared with standard FEM calculation.

To achieve this, we consider the academic case of a rectangular simply-supported plate excited by a point force F and coupled to a parallelepiped rigid-walled cavity filled either with air, or water, as shown in Fig. (3). This example is voluntarily simplified in order to facilitate studying the convergence parameters of both methods, namely the interface FE mesh size and the number of modes taken into account in the cavity-PTF computation. The simulation parameters are listed in table 1.

Convergence of FEM

Before using the FE result as a reference, it is important to study the convergence of the method in order to define the optimal interface mesh criterion. To do this, a coincident mesh is used at the structure-cavity interface and a direct resolution of the coupled problem using NASTRAN is performed to avoid modal truncation effects on the results presented. The interface mesh criterion, which provides the size of the finite elements at the structure-cavity interface, is based on the smallest wavelength λ of the uncoupled plate and cavity and computed at the maximal frequency of the range of interest f max (see Eq. ( 18)).

       λ s = √ 2πf -1/2 max 4 Eh 2 12ρ(1-ν 2 ) λ f = c fmax λ = min(λ s , λ f ) (18)
Four interface mesh criteria are used to study the convergence of FEM ( λ 6 , λ 12 , λ 28 and λ 60 ) when the plate is coupled to a water-filled cavity and the mean square velocity of the plate is computed for each of them. This measurement is usually used in vibro-acoustics (see Ref. [START_REF] Tournour | Validation, Performance, Convergence and Application of Free Interface Component Mode Synthesis[END_REF]) and provides an overview of the response of a system. Fig. [START_REF] Sum | On Acoustic and Structural Modal Cross-Couplings in Plate-Cavity Systems[END_REF] shows that convergence in heavy fluid requires a very fine interface mesh to capture the interface dynamics, since convergence occurs for a λ 28 mesh criterion. This difficulty is linked to coupling with heavy fluids such as water. For a cavity filled with air, convergence is much better and the usual λ 6 mesh criterion is adequate to properly describe the response of the coupled system.

The presence of the cavity leads to an added mass effect on the plate that decreases the flexural wave speed and thus the flexural wavelength. As a consequence, the standard λ 6 interface mesh criterion, based on the wavelength of the uncoupled subsystems, leads to an underestimation of the size of the interface finite elements necessary for describing the behavior of the coupled system.

Convergence of PTF method

Light fluid case

Ouisse et al. demonstrated in [START_REF] Ouisse | Patch Transfer Functions as a Tool to Couple Linear Acoustic Problems[END_REF] that a λ 2 patch mesh criterion for the coupling surface was sufficient to obtain good results when two acoustic air-filled cavities are coupled. However, this patch mesh criterion is also valid for the structure-cavity coupling presented in Fig. [START_REF] Wc Hurty | Dynamic Analysis of Structural Systems using Component Modes[END_REF]. In this example, the simulation parameters are given in table 1. Furthermore, the structural modal basis contains 12 modes, corresponding to resonance modes below 508 Hz, while the acoustic basis contains 9 modes, corresponding to resonance frequency modes below 515 Hz.

In Fig. [START_REF] Wc Hurty | Dynamic Analysis of Structural Systems using Component Modes[END_REF], the plate response is presented for a classical FE computation with a λ 6 mesh criterion and a PTF calculation with a patch size of λ 2 that ensures the convergence of the proposed method. Consequently, the PTF method gives good results with a patch mesh coarser than a standard FE interface mesh.

Heavy fluid case

Unlike the light fluid case, in heavy fluid, account has to be taken of the decrease of wavelength described in section 5.1, since the structure-cavity coupling is strengthened. Consequently, we used a λ 6 patch mesh criterion, i.e. 56 patches. To avoid a mode truncation effect, structure-PTF and cavity-PTF are obtained by direct resolution of Eqs. ( 10) and [START_REF] Kim | A Compact Matrix Formulation using the Impedance and Mobility Approach for the Analysis of Structural-Acoustic Systems[END_REF]. Fig. [START_REF] Tournour | Validation, Performance, Convergence and Application of Free Interface Component Mode Synthesis[END_REF] shows that in this case too the PTF method gives satisfying results for a patch mesh coarser than the interface FE mesh. 6. Introduction of residual shapes in the cavity-PTF computation

Influence of the mode truncation effect on the cavity-PTF computation in heavy fluid

As in section 5.2.2, the PTF can be obtained by direct resolution of FE models. However, when the size of the model is large, it is preferable to reduce it by using mode expansion. Unfortunately, in heavy fluid the results obtained are sensitive to the truncation effect of acoustic cavity modes, as shown in Fig. [START_REF] Hjp Morand And R Ohayon | Fluid-Structure Interaction[END_REF]. 200, 300 and 400 acoustic cavity modes, corresponding to resonant frequency modes below 8.8 kHz, 10.5 kHz and 11.7 kHz, were used to compute the cavity-PTF Z c jk , while 5 structural modes, corresponding to resonant frequency modes below 205 Hz, were used to compute the structure-PTF Y s jk . This was sufficient to ensure its convergence since unchanged results were observed with more structural modes. This figure shows that more than 400 acoustic cavity modes are required to ensure the convergence of the PTF method. ), (-.-) PTF method for 300 acoustic cavity modes (patch size: λ 6 ) and (. . . ) PTF method for 400 acoustic cavity modes (patch size: λ 6 )

Introduction of residual shapes -methodology

To achieve convergence of the PTF method in heavy fluid, high-order cavity modes are very important to provide a good description of the coupling mechanism. To overcome this difficulty, the concept of residual modes (see Ref. [START_REF]A Hybrid Method of Component Mode Synthesis[END_REF][START_REF] Rubin | Improved Component-Mode Representation for Structural Dynamic Analysis[END_REF]) is used for the cavity-PTF computation. The residual shapes do not correspond to physical normal modes but allow enriching the modal basis with a set of functions describing the contribution of high-order modes. In the PTF method, the idea consists in computing the residual shapes or attachment modes for the acoustic cavity independently of the structure in order to keep the substructuring aspect of the method. In the standard residual mode technique, the residual functions for an acoustic cavity correspond to the response of the cavity excited by structural normal modes [START_REF] Qh Tran | Comparison of Robust CMS Methods for Damped Vibroacoustic Problems[END_REF][START_REF] Roy | Efficient Computation of the Radiated Sound Power of Vibrating Structures using a Modal Approach[END_REF]. Consequently, a modification of the structure leads to a modification of structural normal modes that requires computing the residual shapes again. In our case, a residual shape corresponds to the response of the cavity excited by a constant normal displacement imposed on a patch p at a specific radian frequency ω c . Consequently, a structural modification does not change the residual cavity shapes as long as the patch mesh remains unchanged. Furthermore, this definition leads to a number of residual shapes equal to the number of patches. This is a limitation of the method in the high frequency range, since the number of patches increases with frequency. Using the notations used in section 4.2, the residual shapes {P r } of the cavity are computed by the direct resolution of the following equation:

K f -ω 2 c M f {P r } = {Q p } (19)
In Eq. ( 19), the excitation term Q p corresponding to the unit normal displacement ūc p imposed on the patch p is given by the following equation:

Q p = ω 2 c Ap N ūc p on S p 0 on S v \ S p (20)
Thus we obtain a new reduction basis containing the N physical modes kept in the original basis, enriched by the N r residual shapes {P r } (see Eq. ( 21)),

[T ] = [φ 1 , . . . , φ N |P 1 , . . . , P Nr ] (21)
The reduction basis [T ] has to be reorthogonalized to maintain the validity of the mass and stiffness orthogonality relations. To do this, the reduced eigenvalue problem given by Eq. ( 22) has to be solved [START_REF] Balmès | Theoretical Course Notes: Methods for Vibration Design and Validation[END_REF],

[T ]

T [K f ] [T ] -ω 2 [T ] T [M f ] [T ] {ϕ R } = {0} (22)
Finally, Eq. ( 22) leads to a new reduction basis {Φ} containing N physical modes and N r residual shapes (see Eq. ( 23)).

{Φ} = [T ] {ϕ R } (23)
An example of the residual shape computed from Eq. ( 23) is plotted in Fig. [START_REF] Ohayon | Reduced Models for Fluid-Structure Interaction Problems[END_REF]. As expected, the residual shape varies slowly along the three directions of the cavity. It can also be seen that the residual pressure shape reaches its maximum in the vicinity of the excited patch whose center coordinates are (x,y,z)=(0.22 m, 0.14 m, 0 m).

Numerical results

In the example presented in this article, the new reduction basis contains 2 physical modes, corresponding to resonant frequency modes below 750 Hz, and 56 residual shapes computed for ω c = 942.5rad.s -1 (150 Hz). When using the residual shapes, the mean square velocity computed by the PTF method with a λ 6 patch mesh is in very good agreement with the reference result obtained by FEM, contrary to that computed by the PTF method without residual shapes, as shown in Fig. [START_REF] Tournour | Pseudostatic Corrections for the Forced Vibroacoustic Response of a Structure-Cavity System[END_REF]. This figure shows the advantage of using residual shapes, since it provides a visibly spectacular improvement of convergence.

Fig. [START_REF] Tournour | Pseudostatic Corrections for the Forced Vibroacoustic Response of a Structure-Cavity System[END_REF] shows that the PTF method gives good results with a coarser interface mesh than that required in standard FEM to ensure the convergence of the method. The better convergence of the PTF method can be explained in part by the independent convergence of each subsystem obtained with a standard λ 6 mesh criterion and the coupling technique that couples subsystems through patches (surfaces) instead of nodes. The PTF method associated with the residual shapes technique therefore provides effective time-saving computation compared to standard FEM.

In the example shown in Fig. [START_REF] Tournour | Pseudostatic Corrections for the Forced Vibroacoustic Response of a Structure-Cavity System[END_REF], the value of ω c corresponds to 20% of the first non-zero cavity mode.

In fact the value that must be used to compute the residual shapes is not clearly established, since for Tran et al. [START_REF] Qh Tran | Comparison of Robust CMS Methods for Damped Vibroacoustic Problems[END_REF] and Cuppens et al. [START_REF] Cuppens | Evaluation of the FRF Based Substructuring and Modal Synthesis Technique Applied to Vehicle FE Data[END_REF], it re-spectively corresponds to 50% and 75% of the maximal frequency of the range of interest. However, to thoroughly study the influence of this parameter, the problem is solved for three values of ω c : 314.1 rad.s -1 (50 Hz), 942.5 rad.s -1 (150 Hz) and 6,283.2 rad.s -1 (1000 Hz). It is important to bear in mind that, in the present case, physical acoustic cavity modes are computed up to 750 Hz. In other words, the first two values of ω c are below the frequency of the last physical mode kept in the modal basis, while in the last case ω c is above this frequency. Fig. [START_REF] Qh Tran | Comparison of Robust CMS Methods for Damped Vibroacoustic Problems[END_REF] shows that the convergence of the PTF method is ensured for ω c = 314.1 rad.s -1 and ω c = 942.5 rad.s -1 , which is not at all the case for ω c = 6283.2 rad.s -1 , since the contribution of modes between 750 and 1000 Hz is absent. Thus, the choice of ω c is not arbitrary. However, convergence is achieved as long as the value of this parameter remains below the frequency of the last physical mode retained to compute the cavity-PTF.

Radiation of a bulkhead into two ballasts

The advantage of the residual shapes concept is illustrated on an industrial example. The proposed case concerns two water-filled ballasts, considered as rigid-walled acoustic volumes, which are separated by a clamped bulkhead mechanically excited, as illustrated in Fig. [START_REF] Roy | Efficient Computation of the Radiated Sound Power of Vibrating Structures using a Modal Approach[END_REF]. The first ballast has dimensions 4 m (width) by 2 m (height) by 2.75 m (depth) approximately, while the second ballast has dimensions 4 m (width) by 2 m (height) by 3.55 m (depth) approximately. Thus, the bulkhead, which is a steel shell of thickness 0.02 m, has dimensions 4 m×2 m approximately. The excitation consists of a prescribed point force acting on the bulkhead along the x-axis.

A FE model of the whole system is used to compute the basic quantities of the PTF method. Table 2 sums up the main features of the FE mesh of the whole system. The aim is to build a PTF model in which both ballasts and the bulkhead are modeled by their PTF. Then, the influence of residual shapes for both ballasts is studied by comparing the pressure level at one point of each ballast, when the cavity-PTF are obtained either from a direct computation using FEM or a mode expansion with or without residual shapes. The structure-PTF are computed from a mode expansion insofar as they converge rapidly.

In the present case, the PTF equations, namely Eqs. ( 6)-( 9), has to be rewritten to take into account the presence of the second ballast, corresponding to a third subsystem. To this end, continuity equations are rewritten as expressed by the following equation:

ps j = pb1 j -pb2 j vs j = vb1 j = -v b2 j = vj (24)
These continuity conditions are used to calculate the coupling normal velocity vj on each patch j. Again, the introduction of continuity conditions in Eq. ( 7) allows obtaining the coupling normal velocity vj given by the following equation:

vj = I -Y s ji Z b1 ik + Z b2 ik -1 vs k ( 25
)
where Z b1 ik is the cavity-PTF of the first ballast and Z b2 ik is the cavity-PTF of the second ballast.

The pressure at point L (resp. M ) belonging to the first ballast (resp. the second one) is thus obtained from the coupling velocity by using the cavity-PTF Z Lj calculated between the point L and the patch j.

p b1 (L) = Z b1 Lj vj p b2 (M ) = Z b2 Lj vj (26)
PTF calculations have been performed for 329 patches. The results obtained from a direct computation of the cavity-PTF are compared in Fig. [START_REF] Cuppens | Evaluation of the FRF Based Substructuring and Modal Synthesis Technique Applied to Vehicle FE Data[END_REF] with results yielded by a mode expansion of the cavity-PTF with or without introducing the residual shapes. It can be also noticed that the mode expansion have been performed using normal modes computed up to 1600 Hz. This corresponds to 168 normal modes for the first ballast and 171 normal modes for the second one. Fig. [START_REF] Cuppens | Evaluation of the FRF Based Substructuring and Modal Synthesis Technique Applied to Vehicle FE Data[END_REF] shows a very good agreement between a direct computation of the cavity-PTF and a modal expansion introducing the residual shapes, while the calculation using resonant modes only does not converge.

The introduction of residual shapes allows decreasing in computational efforts as highlighted in Table 3, where a comparison of total CPU times is performed for a direct computation of the cavity-PTF and a modal calculation with or without the residual shapes. This table shows that the total CPU time for a modal approach introducing the residual shapes is much lower than that of a direct calculation, while the accuracy is equivalent for both approach. It can be noticed that the CPU time for a modal calculation with or without the residual shapes are similar. Consequently, the introduction of residual shapes improves the accuracy of the modal approach, without significantly increasing the computational effort.

Conclusion

The PTF substructuring approach allows computing the response of structurecavity systems. It has the advantage of solving vibro-acoustic problems with reduced computational cost compared to standard FEM. The coupling technique consists in discretizing the coupling surface in elementary areas called patches on which continuity equations are verified. The novelty of the present paper concerns the introduction of the residual shapes concept in order to accelerate the convergence of the PTF method, when solving a structure-acoustic coupling in heavy fluid. The studies presented in this article show that the concept of residual modes associated with interface coupling patches permits substantially improving the convergence of the PTF approach while keeping the uncoupled resolution of equations for each subsystem.

Figure 1 :

 1 Figure 1: Basic vibro-acoustic problem The coupled structural problem is governed by its equation of motion with boundary conditions imposed on the structure boundary ∂S v . The structure is excited by a harmonic point force {F } and the cavity pressure-induced force {F c } = -[A] {P } acting on the structure. The displacement of the structure is harmonic {U s (M, t)} = {U s (M )} e jωt and satisfies the following equation:

Figure 2 :

 2 Figure 2: Use of incompatible mesh -(a) Definition of the PTF model and (b) Detection of FE nodes included in patches: (•) nodes of the first subsystem and (⋄) nodes of the second subsystem

Figure 3 : 1 ×

 31 Figure 3: Definition of the rectangular simply-supported plate excited by a point force F and coupled to a parallelepiped water-filled cavity

1 )Figure 4 :Figure 5 :

 145 Figure 4: Convergence of FEM in heavy fluid -Comparison of the mean square velocity obtained by direct resolution for 4 interface mesh criteria: (-) λ 6 , (--) λ 12 , (-.-) λ 28 and (. . . ) λ 60

1 )Figure 6 :

 16 Figure 6: Convergence of the PTF method in heavy fluid -Comparison of mean square velocity for: (-) FEM solution (interface element size: λ 28 ), (--) PTF method (patch size: λ 6 )

1 )Figure 7 :

 17 Figure 7: Truncation effect of acoustic cavity modes -Comparison of mean square velocity for: (-) FEM solution (interface element size: λ 28 ), (--) PTF method for 200 acoustic cavity modes (patch size: λ 6 ), (-.-) PTF method for 300 acoustic cavity modes (patch size: λ 6 ) and (. . . ) PTF method for 400 acoustic cavity modes (patch size: λ 6 )

Figure 8 : 1 )Figure 9 :

 819 Figure 8: Example of residual pressure shape -Spatial distribution of acoustic pressure (a) on boundaries of the cavity and (b) in the cavity

1 )Figure 10 :

 110 Figure 10: Influence of ωc: (-) FEM solution, (--) ωc = 314.1 rad.s -1 (50 Hz), (-.-) ωc = 942.5 rad.s -1 (150 Hz) and (. . . ) ωc = 6, 283.2 rad.s -1 (1000 Hz)

Figure 11 :

 11 Figure 11: Definition of the FE model

Figure 12 :

 12 Figure 12: Influence of residual shapes -Comparison of the pressure level (a) at the point (-2.55 m, 1.77 m, -0.66 m) of the first ballast and (b) at the point (0.46 m, 0.59 m, -0.82 m) of the second ballast when the cavity-PTF are obtained from (-) a direct calculation, (--) a mode expansion without the residual shapes and ( * ) a mode expansion with residual shapes

Table 2 :

 2 Main features of the FE mesh of the whole system

	System	Number of degrees of freedom
	Ballast 1	27 864 dofs
	Ballast 2	29 147 dofs
	Shell	9 630 dofs

Table 3 :

 3 Comparison of the total CPU times for the calculation of the cavity-PTF

	direct	normal modes normal modes + residual shapes
	Ballast 1 35223 s	1423 s	1456 s
	Ballast 2 38103 s	1456 s	1465 s

Acknowledgment

The authors are grateful to the French Ministry of Defense for funding this work.


Agostino Martinelli
email: agostino.martinelli@inria.fr.

The Unicycle in

published or not. The documents may come   L'archive ouverte pluridisciplinaire

INTRODUCTION

In many mobile robotics applications one fundamental requirement for a mobile robot is its ability to estimate its own configuration during the motion (localization problem).

It has become praxis in robotics to provide an observability analysis prior to solving an estimation problem. Usually this is a proof that the system is locally weakly observable [START_REF] Hermann | Nonlinear Controllability and Observability[END_REF] according to nonlinear observability theory. In the last decade, a great effort has been devoted to analyze the observability properties in several estimation robotics problems (e.g., in SLAM [START_REF] Bryson | Observability Analysis and Active Control for Airbone SLAM[END_REF][START_REF] Huang | Analysis and improvement of the consistency of extended Kalman filter based SLAM[END_REF][START_REF] Lee | On the Observability and Observability Analysis of SLAM[END_REF][START_REF] Perera | On the linear and nonlinear observability analysis of the SLAM problem[END_REF][START_REF] Vidal-Calleja | On the Observability of Bearing-only SLAM[END_REF], in visual-inertial sensor fusion [START_REF] Hesch | Towards consistent vision-aided inertial navigation[END_REF][START_REF] Jones | Visual-inertial navigation, mapping and localization: A scalable real-time causal approach[END_REF][START_REF] Kelly | Visual-inertial simultaneous localization, mapping and sensor-to-sensor selfcalibration[END_REF][START_REF] Kottas | On the consistency of vision-aided inertial navigation[END_REF][START_REF] Li | Improving the accuracy of EKF-based visualinertial odometry[END_REF][START_REF] Martinelli | Vision and IMU data fusion: closedform solutions for attitude, speed, absolute scale and bias determination[END_REF][START_REF] Martinelli | Observability Properties and Deterministic Algorithms in Visual-Inertial Structure from Motion[END_REF][START_REF] Mirzaei | A Kalman filterbased algorithm for IMU-camera calibration: Observability analysis and performance evaluation[END_REF]). Note that in all these works, the dynamics of the considered systems are assumed either unaffected by disturbances or, in presence of disturbances, these are assumed to be constant in time.

Very recently, also the case when some of the inputs are unknown and time dependent has been considered in the framework of mobile robotics [START_REF] Belo | 3 Known Landmarks are Enough for Solving Planar Bearing SLAM and Fully Reconstruct Unknown Inputs[END_REF][START_REF] Martinelli | Observability Properties and Deterministic Algorithms in Visual-Inertial Structure from Motion[END_REF]. In [START_REF] Belo | 3 Known Landmarks are Enough for Solving Planar Bearing SLAM and Fully Reconstruct Unknown Inputs[END_REF] the authors investigate the observability properties of a fundamental problem in the framework of mobile robotics (the bearing SLAM). In [START_REF] Martinelli | Observability Properties and Deterministic Algorithms in Visual-Inertial Structure from Motion[END_REF] the author derives the observability properties in the case of robots equipped with visual and inertial sensors, when some of the inertial sensors are missing. In both these cases, the authors introduced new methods to solve a fundamental problem in control theory that is the Unknown Input Observability problem. The problem of state observability for systems driven by unknown inputs (UI) is a fundamental problem in control theory. This problem was introduced and firstly investigated in the seventies [START_REF] Basile | On the observability of linear, time invariant systems with unknown inputs[END_REF][START_REF] Bhattacharyya | Observer design for linear systems with unknown inputs[END_REF][START_REF] Guidorzi | On Wonham stabilizability condition in the synthesis of observers for unknowninput systems[END_REF][START_REF] Wang | Observing the states of systems with unmeasurable disturbance[END_REF]. A huge effort has then been devoted to design observers for both linear and nonlinear systems in presence of UI. On the other hand, the problem of finding a simple analytic condition that allows us to check the weak local observability of the state is still open in the nonlinear case and for any number of UI. The methods introduced in [START_REF] Belo | 3 Known Landmarks are Enough for Solving Planar Bearing SLAM and Fully Reconstruct Unknown Inputs[END_REF][START_REF] Martinelli | Observability Properties and Deterministic Algorithms in Visual-Inertial Structure from Motion[END_REF] only provide sufficient conditions for the state observability. They are based on a suitable state extension.

Very recently, this fundamental problem has been solved in the general case [START_REF] Martinelli | Nonlinear Unknown Input Observability: The General Analytic Solution[END_REF]. In other words, the analytic criterion that allows us to check the weak local observability of the state in the general nonlinear case and characterized by any number of unknown inputs has been introduced. This analytic tool is based on a simple recursive algorithm that computes the observable codistribution for a nonlinear system whose dynamics are driven by multiple known inputs, a nonlinear drift and multiple unknown inputs (note that the case of a single unknown input was previously dealt in [START_REF] Martinelli | Nonlinear Unknown Input Observability: Analytical expression of the observable codistribution in the case of a single unknown input[END_REF]). In this paper, in section 3, we provide a concise description of the analytic criterion. Then, the paper adopts the analytic criterion in order to investigate the observability properties of a mobile robot that moves on a planar surface by satisfying the unicycle dynamics and that is equipped with exteroceptive sensors. The robot dynamics and observation are defined in section 2. In accordance with the unicycle dynamics, the motion is powered by two independent controls, which are the linear and the angular speed, respectively. We assume that both these speeds are known. In other words, we consider these two speeds as two known inputs of our system. We consider the case when the robot motion is affected by a disturbance that produces an additional (unknown and time-dependent) robot speed along a fixed direction. We consider both the case when this fixed direction is known and unknown. The goal of the paper is to obtain the observability properties of the state that characterizes the robot configuration. Note that the disturbance is an unknown and time dependent input of our system. Hence, our system is characterized by two known inputs and one unknown input. In section 4 we derive the observability properties of the systems defined in section 2, by using the analytic tool provided in section 3. In section 5, we simulate the aforementioned system and we show that, a simple estimator based on an Extended Kalman Filter, provides results that agree with what we could expect from the observability analysis. Finally, conclusions are provided in section 6

The considered system

We consider a vehicle that moves on a 2D-environment. The configuration of the vehicle in a global reference frame can be characterized through the vector [x v , y v , θ v ] T , where x v and y v are the Cartesian vehicle coordinates, and θ v is the vehicle orientation. We assume that the vehicle motion satisfies the unicycle dynamics and we also consider the case when the robot motion is affected by a disturbance that produces an additional (and unknown) robot speed (denoted by w) along a fixed direction (denoted by γ). Hence, the dynamics are characterized by the following differential equations:

(2.1)

   ẋv = v cos θ v + w cos γ ẏv = v sin θ v + w sin γ θv = ω v
where v and ω v are the linear and the rotational speed, respectively, in absence of disturbance (w = 0). We assume that these two speeds are known (we refer to them as to the known inputs), w is unknown (we refer to it as to the unknown input or disturbance) and γ is constant in time and we consider both the case when it is known and unknown. We consider the following three cases of output (see also figure 1 for an illustration):

Figure 1: The vehicle state together with the three considered outputs.

1. the distance from the origin, denoted by r (e.g., a landmark is at the origin and its distance is measured by a range sensor);

2. the bearing of the origin in the local frame, denoted by β (e.g., a landmark is at the origin and its bearing angle is measured by an on-board camera);

3. the bearing of the vehicle in the global frame, denoted by φ (e.g., a camera is placed at the origin).

Extended Observability Rank Condition

This section provides the analytic criterion to investigate the observability properties of a nonlinear system whose dynamics are driven by multiple known inputs a nonlinear drift and multiple unknown inputs (or disturbances). In absence of unknown inputs, the analytic criterion is the observability rank condition introduced by Hermann and Krener [START_REF] Hermann | Nonlinear Controllability and Observability[END_REF]. We call the new criterion, able to account for the presence of unknown inputs, the extended observability rank condition. The derivation of this criterion is available on the book in [START_REF] Martinelli | Nonlinear Unknown Input Observability: The General Analytic Solution[END_REF]. In the driftless case and with a single unknown input, the criterion was introduced in [START_REF] Martinelli | Nonlinear Unknown Input Observability: Analytical expression of the observable codistribution in the case of a single unknown input[END_REF].

We will refer to a nonlinear control system with m u known inputs (u

[u 1 , • • • , u mu ] T ) and m w unknown inputs or disturbances (w [w 1 , • • • , w mw ] T
). The state is the vector x ∈ M , with M an open set of R n . We assume that the dynamics are nonlinear with respect to the state and affine with respect to the inputs (both known and unknown). Finally, for the sake of simplicity, we will refer to the case of a single output y (the extension to the case of multiple outputs is trivial and is given in [START_REF] Martinelli | Nonlinear Unknown Input Observability: The General Analytic Solution[END_REF]). Our system is characterized by the following equations:

(3.2)        ẋ = g 0 (x) + mu i=1 f i (x)u i + mw j=1 g j (x)w j y = h(x)
where g 0 (x), f i (x), i = 1, • • • , m u , and g j (x), j = 1, • • • , m w , are vector fields in M and the function h(x) is a scalar function defined on the open set M . Finally, we assume that the unknown inputs w 1 , • • • , w mw are analytic functions of time.

The system in (2.1) is characterized by m w = 1 and dynamics linear (and not affine) in the inputs (i.e., without the term g 0 ). In other words, the system belongs to the following class of systems:

(3.3)      ẋ = mu i=1 f i (x)u i + g(x)w y = h(x)
These systems are precisely the ones dealt in [START_REF] Martinelli | Nonlinear Unknown Input Observability: Analytical expression of the observable codistribution in the case of a single unknown input[END_REF]. We now outline all the steps to investigate the weak local observability at a given point x 0 of a nonlinear system characterized by (3.3). The reader is addressed to [START_REF] Martinelli | Nonlinear Unknown Input Observability: The General Analytic Solution[END_REF] (sections 3.3 and 4.2) for the general case described by (3.2). Basically, these steps are the steps necessary to compute the observable codistribution (i.e., the steps of algorithms 2 below) and to prove that the differential of a given state component belongs to this codistribution.

In the sequel, we will denote with the symbol D the differential with respect to the state x. For instance, if

x = [x 1 , x 2 ] T and h = x 1 + x 2 2 , we have 1 : Dh = Dx 1 + 2x 2 Dx 2 .
For a given codistribution Ω and a given vector field f = f (x) (both defined on the open set M ), we denote by L f Ω the codistribution whose covectors are the Lie derivatives along f of the covectors in Ω. We remind the reader that the Lie derivative of a scalar function h(x) along the vector field f (x) is defined as follows:

L f h ∂h ∂x f
which is the product of the row vector ∂h ∂x with the column vector f . Hence, it is a scalar function. Additionally, by definition of Lie derivative of covectors, we have:

L f Dh = DL f h. Finally, given two vector spaces V 1 and V 2 , we denote by V 1 + V 2 their sum, i.e.
, the span of all the generators of both V 1 and V 2 . To start we need to introduce the following scalar:

1 The span of the differentials of a set of scalar functions is a codistribution. The reader non familiar with the theory of distributions can simply consider the differential as the gradient operator. The gradient of a scalar function is a line vector. For instance, if x = [x 1 , x 2 ] T and h = x 1 + x 2 2 , we obtain for its gradient the line vector function [1, 2x 2 ]. Later, we adopt this representation. According to this, a codistribuion will be the span of a set of line vectors and a covector (i.e., an element of a codistribution) will be a line vector.

(3.4) L 1 g L g h
The analytic computation is based on the assumption that L 1 g = 0 on a given neighbourhood of x 0 (the reader is addressed to [START_REF] Martinelli | Nonlinear Unknown Input Observability: The General Analytic Solution[END_REF] to see how to proceed in the case when this assumption is not satisfied). Finally, we need to introduce a new set of vector fields i φ m ∈ R n (i = 1, • • • , m u and for any integer m). They are obtained recursively by the following algorithm:

Algorithm 1 1. i φ 0 = f i ; 2. i φ m = [ i φm-1, g] L 1 g
where the parenthesis [•, •] denote the Lie bracket of vector fields, defined as follows:

[a, b] ∂b ∂x a(x) - ∂a ∂x b(x)
In other words, for each i = 1, • • • , m u , we have one new vector field at every step of the algorithm. In the case when m u = 1, we denote by φ m the vector field 1 φ m . The algorithm that generates the entire observable codistribution for the system in (3.3) is the following: Algorithm 2 Observable codistribution in the case m w = 1 and g 0 = 0

1. Ω 0 = span{Dh}; 2. Ω m = Ω m-1 + mu i=1 L f i Ω m-1 + L g L 1 g Ω m-1 + mu i=1 span Li φm-1 Dh
It is possible to prove that this algorithm converges in less than n + 2 steps (n is the dimension of the state). The convergence criterion is trivial in the case dealt by the following lemma (the proof can be found in [START_REF] Martinelli | Nonlinear Unknown Input Observability: The General Analytic Solution[END_REF]).

Lemma 1 Let us denote by Λ i j the distribution generated by i φ 0 , i φ 1 , • • • , i φ j (i.e., the vectors obtained by running algorithm 1 for a given i = 1, • • • , m u ) and by m i (≤ n-1) the smallest integer for which

Λ i mi+1 = Λ i mi (n is the dimension of the state x). If Li φj L 1 g = 0, ∀i = 1, • • • , m u and ∀j = 0, • • • , m i ,
the convergence of algorithm 2 occurs in at most n -1 steps and it occurs at the smallest integer j such that Ω j+1 = Ω j .

The steps to investigate the weak local observability at a given point x 0 of a nonlinear system characterized by (3.3) are reported below. Note that, in the trivial case analyzed by the previous lemma, the criterion provided below simplifies, since we do not need to compute the 2 , and we do not need to check that its differential belongs to the codistribution computed at every step of algorithm 2. In practice, we skip the steps 4 and 5 in the procedure below.

quantity τ = L 2 g h (L 1 g h)
1. For the chosen x 0 , compute L 1 g (= L 1 g h). In the case when L 1 g = 0, choose another function in the space of functions F (defined as the space that contains h and its Lie derivative up to any order along the vector fields f 1 , • • • , f mu ) such that its Lie derivative along g does not vanish 2 .

2. Compute the codistribution Ω 0 and Ω 1 (at x 0 ) by using algorithm 2.

3. Compute the vector fields i φ m (i = 1, • • • , m u ) by using algorithm 1, starting from m = 0, to check if the considered system is in the special case dealt by lemma 1. In this trivial case, set m = 0, use the recursive step of algorithm 2 to build the codistribution Ω m for m ≥ 2, and skip to step 6.

4. Compute τ = L 2 g h (L 1 g ) 2
and Dτ . 7. If the differential of a given state component (

x j , j = 1, • • • , n) belongs to Ω * (namely if Dx j ∈ Ω * )
on a given neighbourhood of x 0 , then x j is weakly locally observable at x 0 . If this holds for all the state components, the state x is weakly locally observable at x 0 . Finally, if the dimension of Ω * is smaller than n on a given neighbourhood of x 0 , then the state is not weakly locally observable at x 0 .

4 Observability properties of the system defined in section 2

In this section we derive the observability properties of the system defined in section 2 by using the analytic criterion presented in the previous section. We deal with both the cases when γ is known and unknown (in section 4.1 and 4.2, respectively). For both cases, we consider the three outputs defined in section 2, separately. For simplicity sake, we actually consider the following three outputs:

y = r 2 = x 2 v + y 2 v instead of y = r, y = tan β = yv-xv tan θv
xv+yv tan θv instead of y = β and y = tan φ = yv xv instead of y = φ. Obviously, the result of the observability analysis does not change. As it will be seen, in all the cases, the systems always meet the assumptions of lemma 1. Hence, we will skip the fourth and the fifth step in the procedure given in section 3. Note that the most of systems investigated in [START_REF] Martinelli | Nonlinear Unknown Input Observability: The General Analytic Solution[END_REF] do not satisfy these assumptions and for them we had to compute the scalar τ and verify that its differential belongs to the codistribution returned by algorithm 2.

Observability properties when γ is known

The state is [x v , y v , θ v ] T and its dynamics are provided by the three equations in (2.1), where γ is a known parameter. From (2.1) and (3.3) we easily obtain:

m u = 2, u 1 = v, u 2 = ω v , f 1 =   cos θ v sin θ v 0   , f 2 =   0 0 1   , g =   cos γ sin γ 0  
We follow the seven steps mentioned in section 3.

y = r 2

First Step We compute L 1 g in (3.4). We obtain:

L 1 g = 2(x v cos γ + y v sin γ) Second Step We have Ω 0 = span{[x v , y v , 0]} and Ω 1 = span{[x v , y v , 0], [cos θ v , sin θ v , y v cos θ v -x v sin θ v ] }.

Third

Step We compute 1 φ 1 and 2 φ 1 through algorithm 1. We obtain:

1 φ 1 = 2 φ 1 =   0 0 0   . As a result,
all the subsequent steps of algorithm 1 provide null vectors. Therefore, the assumptions of lemma 1 are trivially met. We set m = 0 and we skip to the sixth step.

Sixth

Step We use algorithm 2 to compute Ω 2 . We obtain Ω 2 = Ω 1 + dL g L 1 g

L f 1 h. Hence, we obtain that Ω 2 has dimension equal to 3, which is the dimension of the state. Hence, algorithm 2 has converged.

Seventh

Step We conclude that the state is weakly locally observable.

y = tan β

First Step We compute L 1 g in (3.4). We obtain:

L 1 g = - y v cos γ -x v sin γ x 2 v cos 2 θ v + 2 sin θ v cos θ v x v y v -y 2 v cos 2 θ v + y 2 v ,

Second

Step By an explicit computation (by using algorithm 2) we obtain that the dimension of Ω 0 is 1 and the dimension of Ω 1 is 2.

Third

Step We compute 1 φ 1 and 2 φ 1 through algorithm 1. They vanish as in the previous case. As a result, all the subsequent steps of algorithm 1 provide null vectors. Therefore, the assumptions of lemma 1 are trivially met. We set m = 0 and we skip to the sixth step.

Sixth

Step By using algorithm 2 we obtain that Ω 2 has dimension equal to 3, which is the dimension of the state. Hence, algorithm 2 has converged.

Seventh

Step We conclude that the state is weakly locally observable.

y = tan φ

First Step We compute L 1 g in (3.4). We obtain:

L 1 g = - y v cos γ -x v sin γ x 2 v Second Step We have Ω 0 = span{[-y v , x v , 0]}.
Additionally, we obtain that the dimension of Ω 1 is 2.

Third

Step We compute 1 φ 1 and 2 φ 1 through algorithm 1. They vanish as in the previous case. As a result, all the subsequent steps of algorithm 1 provide null vectors. Therefore, the assumptions of lemma 1 are trivially met. We set m = 0 and we skip to the sixth step.

Sixth

Step By using algorithm 2 we obtain that Ω 2 has dimension equal to 3, which is the dimension of the state. Hence, algorithm 2 has converged.

Seventh

Step We conclude that the state is weakly locally observable.

Observability properties when γ is unknown

The state is [x v , y v , θ v , γ] T and its dynamics are provided by the following four equations:

(4.5)      ẋv = v cos θ v + w cos γ ẏv = v sin θ v + w sin γ θv = ω v γ = 0
From (4.5) and (3.3) we easily obtain:

m u = 2, u 1 = v, u 2 = ω v , f 1 =     cos θ v sin θ v 0 0     , f 2 =     0 0 1 0     , g =     cos γ sin γ 0 0    
As in section 4.1, we follow the seven steps mentioned in section 3 and we consider the same expressions for the outputs. y = r 2

First

Step We obviously obtain the same expression for L 1 g as in the case y = r 2 of section 4.1.

Second

Step

We have Ω 0 = span{[x v , y v , 0, 0]}. Additionally, Ω 1 = span{[x v , y v , 0, 0], [cos θ v , sin θ v , y v cos θ v -x v sin θ v , 0] }.

Third

Step We compute 1 φ 1 and 2 φ 1 through algorithm 1. We obtain:

1 φ 1 = 2 φ 1 =     0 0 0 0     . As a result,
all the subsequent steps of algorithm 1 provide null vectors. Therefore, the assumptions of lemma 1 are trivially met. We set m = 0 and we skip to the sixth step.

Sixth

Step By using algorithm 2 we obtain that Ω 2 has dimension equal to 3. Additionally, by a direct computation, it is possible to check that Ω 3 = Ω 2 meaning that m * = 2 and Ω * = Ω 2 .

Seventh

Step The dimension of the observable codistribution is 3 < 4. We conclude that the state is not weakly locally observable. In particular, since the differential of every state component does not belong to Ω 2 , we conclude that no state component is observable.

y = tan β

First Step We obviously obtain the same expression for L 1 g as in the case y = tan β of section 4.1. Second Step We compute Ω 0 and Ω 1 and we obtain that their dimensions are 1 and 2, respectively.

Third

Step We compute 1 φ 1 and 2 φ 1 through algorithm 1. They vanish as in the previous case. As a result, all the subsequent steps of algorithm 1 provide null vectors. Therefore, the assumptions of lemma 1 are trivially met. We set m = 0 and we skip to the sixth step.

Sixth

Step By using algorithm 2 we obtain that Ω 2 has dimension equal to 3. Additionally, by a direct computation, it is possible to check that Ω 3 = Ω 2 meaning that m * = 2 and Ω * = Ω

Seventh

Step The dimension of the observable codistribution is 3 < 4. We conclude that the state is not weakly locally observable. In particular, since the differential of every state component does not belong to Ω 2 , we conclude that no state component is observable.

y = tan φ

First Step We obviously obtain the same expression for L 1 g as in the case y = tan φ of section 4.1. Second Step We compute Ω 0 and Ω 1 and we obtain that their dimensions are 1 and 2, respectively.

Third Step

We compute 1 φ 1 and 2 φ 1 through algorithm 1. They vanish as in the previous case. As a result, all the subsequent steps of algorithm 1 provide null vectors. Therefore, the assumptions of lemma 1 are trivially met. We set m = 0 and we skip to the sixth step.

Sixth

Step By using algorithm 2 we obtain that Ω 2 has dimension equal to 4, which is the dimension of the state. Hence, algorithm 2 has converged.

Seventh

Step We conclude that the state is weakly locally observable.

Table 1 summarizes the results of the observability analysis carried out in this section. We conclude this section by remarking that these results agree with our expectation. By using the observability rank condition in [START_REF] Hermann | Nonlinear Controllability and Observability[END_REF], we easily obtain that, in absence of the disturbance, the dimension of the observable codistribution is 2 for the first two observations (y = r and y = β) and 3 for the last one (y = φ). In particular, for the first two observations, all the initial states rotated around the vertical axis are indistinguishable. In other words, in these two cases, the system exhibits a continuous symmetry (the reader is addressed to [START_REF] Martinelli | State Estimation Based on the Concept of Continuous Symmetry and Observability Analysis: the Case of Calibration[END_REF][START_REF] Martinelli | Vision and IMU data fusion: closedform solutions for attitude, speed, absolute scale and bias determination[END_REF] for the meaning of continuous symmetry). In presence of the disturbance, when γ is known, the aforementioned system invariance is broken and the entire state becomes observable. When γ is unknown, the symmetry still remains (and obviously also concerns the new state component γ). Note also a very important aspect. The presence of a disturbance improves the observability properties of a system (this regards the case when γ is known). In particular, if w = 0 (absence of disturbance), the state becomes unobservable even if it is known that it is zero, while, when w = 0, the state is observable even if it is unknown. Note that the analytic derivation of the analytic criterion holds when the disturbance does not vanish [START_REF] Martinelli | Nonlinear Unknown Input Observability: The General Analytic Solution[END_REF].

Simulations

The scope of this section is to show that, by using a very simple estimator based on an Extended Kalman Filter, we obtain results which agree with our observability analysis. Note that the goal of this paper is not to derive a new observer for a system characterized by an unknown input. On the other hand, the results obtained by using the aforementioned estimator, are useful to further validate the analytic criterion in section 3.

Simulated trajectories and robot sensors

The trajectories are simulated as follows. The equations in (2.1) are discretized with a time step of 5 * 10 -4 s and they last for 10 s. The linear speed, i.e. v, is constant and equal to 0.1 m s -1 . The angular speed, i.e. ω v , is generated randomly. Specifically, its value is settled at each step and follows a Gaussian distribution with mean value 2π 10 rad s -1 and variance 6π 10 2 rad 2 s -2 . A typical example of trajectory, obtained with this setting, is displayed in fig. 2, left side. The disturbance is generated as follows. The parameter γ is set equal to 0.7 π. The unknown input w is generated randomly. Specifically, its value is settled at each step and follows a Gaussian distribution with mean value 0.05 m s -1 and variance (0.005) 2 m 2 s -2 . A typical example of trajectory is displayed in fig. 2, right side. The robot is equipped with proprioceptive sensors able to measure at each time step the linear and the angular speed. These measurements are affected by errors. Specifically, each measurement is generated at every time step of 5 * 10 -4 s by adding to the true value a random error that follows a Gaussian distribution.

The mean value of this error is zero and the standard deviation is 0.01 times the true value for both the linear and the angular speed.

Regarding the exteroceptive measurements, they are generated at a lower frequency. Specifically, the measurements are generated each 2.5 * 10 -2 s. Also these measurements are affected by errors. Specifically, the measurement is generated by adding to the true value a random error that follows a Gaussian distribution. The mean value of this error is zero and the standard deviation is 0.01 m for the range measurements and 1 deg for the two angular measurements (bearing of the origin in the robot frame and bearing of the robot in the global frame).

Estimation results

We adopt an Extended Kalman Filter that estimates an extended state that includes the unknown input together with its first order time derivative. In other words, in the case when γ is known, the estimated state is: [x v , y v , θ v , w, w (1) ]. Its dynamics is obtained from (2.1) and are:

(5.6)

         ẋv = v cos θ v + w cos γ ẏv = v sin θ v + w sin γ θv = ω v ẇ = w (1)
ẇ(1) = w (2) When γ is unknown, the estimated state also includes γ and the dynamics are given by (5.6) and the additional equation γ = 0. In order to implement the prediction phase of our Extended Kalman Filter we have to provide the value of w (2) , which is unknown. We set this quantity to zero. Note that the simulated trajectory does not satisfy this hypothesis since the disturbance is randomly generated. However, the estimator is able to provide good performance as it is shown in figures 3-5.

Figure 3 displays the estimated trajectory in the case when the output is provided by the range sensor (h = r). Dots blue are the ground truth, red circles the estimated trajectory by only using the knowledge of the proprioceptive measurements (i.e., the measurements of v and ω v ), black stars the trajectory estimated by our estimator. Left side is the case when γ is known and right side when it is unknown. In accordance with our observability analysis, the estimator follows the true trajectory only in the case when γ is known. Figure 4 displays the estimated trajectory in the case when the output is provided by the on-board bearing sensor (h = β). Also in this case, the estimator follows the true trajectory only in the case when γ is known, in accordance with our observability analysis. Finally, figure 5 displays the estimated trajectory in the case when the output is provided by the bearing sensor at the origin (h = φ). In this case, the estimator follows the true trajectory both in the case when γ is known and when it is unknown. Also this result agrees with our observability analysis. We conclude this section by remarking that the results of our simulations fully agree with the observability analysis provided in section 4. Note that we showed the result of a single trial. However, we found always similar results by running many trials.

Conclusion

In this paper we investigated the observability properties of the unicycle in presence of a disturbance. This is a nonlinear system whose dynamics are driven by two known inputs (the linear and the angular speed that would result in absence of the disturbance) and an external unknown input. The unknown input consists of an additional unknown and time dependent robot speed, along a fixed direction. We considered both the case when this fixed direction is known and unknown. The goal of the paper was to obtain the observability properties of the state that characterizes the robot configuration. The novelty of this observability analysis was that it accounted for the presence of an unknown and time dependent disturbance. In order to account for the time dependent disturbance, the paper used the analytic criterion recently introduced in [START_REF] Martinelli | Nonlinear Unknown Input Observability: The General Analytic Solution[END_REF]. Note that this analytic criterion is the solution of a fundamental open problem in control theory. The analytic criterion was previously introduced in [START_REF] Martinelli | Nonlinear Unknown Input Observability: Analytical expression of the observable codistribution in the case of a single unknown input[END_REF] for the easier driftless case and characterized by a single unknown input. We showed that the application of this criterion is very simple and can be done automatically. Additionally, we simulated the aforementioned systems. In particular, in the simulations, the disturbance was generated time dependent. We showed that, a simple estimator based on an Extended Kalman Filter, provides results that fully agree with what we could expect from the observability analysis.

5 . 6 .

 56 Use the recursive step of algorithm 2 to build the codistribution Ω m for m ≥ 2, and, for each m, check if Dτ ∈ Ω m . Denote by m the smallest m such that Dτ ∈ Ω m . For each m ≥ m , check if Ω m+1 = Ω m and denote by Ω * = Ω m * , where m * is the smallest integer such that m * ≥ m and Ω m * +1 = Ω m * (note that m * ≤ n + 2).

Figure 2 :

 2 Figure 2: Typical simulated trajectories. Left side without disturbance, right side with disturbance.

Figure 3 :

 3 Figure 3: Estimated trajectory in the case when the output is h = r. Dots blue are the ground truth, red circles the estimated trajectory by only using the knowledge of the proprioceptive measurements, black stars the trajectory estimated by our estimator. On the left side is the case when γ is known and on the right side when it is unknown.

Figure 4 :

 4 Figure 4: The same as in figure 3 in case when the output is h = β.

Figure 5 :

 5 Figure 5: The same as in figure 3 in the case when the output is h = φ.

Table 1 :

 1 [START_REF] Belo | 3 Known Landmarks are Enough for Solving Planar Bearing SLAM and Fully Reconstruct Unknown Inputs[END_REF] . Weak local observability of the state in all the considered scenarios

	γ	Output State observability
	known	y = r	yes
	known	y = β	yes
	known	y = φ	yes
	unknown	y = r	no
	unknown	y = β	no
	unknown	y = φ	yes

If the Lie derivative of any function in F vanishes, it means that the unknown input can be ignored to obtain the observability properties (the system is not canonic with respect to the unknown input, as shown in[START_REF] Martinelli | Nonlinear Unknown Input Observability: The General Analytic Solution[END_REF]).

This work was supported by the French National Research Agency ANR through the project VIMAD.