S Samagassi

E Jacquelin

A Khamlichi

M Sylla

Bayesian sparse regularization for multiple force identification and location in time domain

Keywords: Multi-force identification, Bayesian regularization, regularization parameter, time domain, compressed sensing, sparse, inverse problem

In this work, reconstruction and location in time domain of multiple forces acting on a linear elastic structure are achieved through a Bayesian approach to solve an inverse problem. The Bayesian solution of the inverse problem is provided in the form of a posterior probability density function. The unknown forces are determined through Markov chain Monte Carlo (MCMC) method, the Gibbs algorithm.

This posterior density integrating both the likelihood and prior information was considered for the particular case of a linear elastic structure. The measurements are affected by an additive random noise. Two particular cases were analyzed: unperturbed and uncertain model representing the structure. The unperturbed model was used to identify a single force. When the model is uncertain, compressed sensing technique was used to provide an adequate sparse representation of the inverse problem through a wavelet basis.

With this sparse representation, the possibility of achieving automatic location of the forces was investigated.

This requires to identify all the degrees-of-freedom along with the identified actions are not vanishing. Also, the possibility of force identification with less sensors than forces was studied. The proposed approach is illustrated and validated on numerical examples. This proposed approach is compared with classical approach of force identification based on Tihonov regularization associated with the GCV criterion.

Introduction

Force identification has been extensively addressed for the last two decades as a consequence of developing structural health monitoring [START_REF] Klinkov | An updated comparison of the force reconstruction methods[END_REF][START_REF] Zhang | Identification of coexistent load and damage[END_REF]. So many works dealt with the possibility to identify the force and the damage of a structure [START_REF] Zhang | Identification of coexistent load and damage[END_REF][START_REF] Fe Gunawan | Time and frequency domains iterative regularization for inverse analysis of an instrumented one-point bend specimen[END_REF][START_REF] Zhang | Simultaneous identification of excitation time histories and parametrized structural damages[END_REF]. This is a well-known ill-posed problem and the early work of Tikhonov [START_REF] Tikhonov | Solutions of Ill-Posed Problems[END_REF] on regularization techniques happens to be the main solution used to solve it. In general an optimization problem must be solved. A direct solution may be written [START_REF] Hansen | A matlab package for analysis and solution of discrete ill posed problems[END_REF][START_REF] Hansen | Rank-Deficient and Discrete III-Posed Problems[END_REF][START_REF] Hansen | Deconvolution and regularization with toeplitz matrices[END_REF] but some authors preferred specific algorithms like the Levenberg-Marquardt [START_REF] Gunawan | Levenberg-Marquardt iterative regularization for the pulse-type impact-force reconstruction[END_REF] or the particle swarm optimization [START_REF] El-Bakari | Assessing impact force localization by using a particle swarm optimization algorithm[END_REF] to solve it. The close relation between singular value decomposition and regularization was highlighted by Hansen in many papers [START_REF] Hansen | A matlab package for analysis and solution of discrete ill posed problems[END_REF][START_REF] Hansen | Rank-Deficient and Discrete III-Posed Problems[END_REF][START_REF] Hansen | Deconvolution and regularization with toeplitz matrices[END_REF] : in particular the notion of filter factor is a comprehensive approach of the regularization principle including Tikhonov regularization, the truncation method [START_REF] Jacquelin | Force reconstruction analysis and regularization of a deconvolution problem[END_REF]. [START_REF] Lin | A new investigation into regularization techniques for the method of fundamental solutions[END_REF] studies different regularization approaches to investigate the solution stability of the method of fundamental solutions (MFS). Three regularization methods in conjunction with two different regularization parameters to find the optimal stable MFS scheme are illustrated.

The forces may be identified by indirect measurements through a relation between the responses and the solicitations. In time domain, this relationship is established thanks to the transfer function (or the Green's function) [START_REF] Jacquelin | Force reconstruction analysis and regularization of a deconvolution problem[END_REF][START_REF] Sturm | Force reconstruction using an adaptive algorithm in time domain[END_REF] or to the system Markov parameters [START_REF] Mao | A state space force identification method based on Markov parameters precise computation and regularization technique[END_REF][START_REF] Wang | Sensor placement methods for an improved force identification in state space[END_REF]. A modal-like method that involved specific basis functions was developed by Liu [START_REF] Liu | An improved method for the reconstruction of a distributed force acting on a vibrating structure[END_REF]. Chen et al. [START_REF] Chen | Impact source identification in finite isotropic plates using a timereversal method: theoretical study[END_REF][START_REF] Chen | Development of time-reversal method for impact source identification on plate structures[END_REF] used the time-reversal approach to identify the forces. Considering that the model and the signals are uncertain, Zhang et al. [START_REF] Zhang | Bayesian force reconstruction with an uncertain model[END_REF] developed a Bayesian approach: by using a Markov chain Monte Carlo method, the procedure may be viewed as an iterative regularization method. Note that usually the structural behaviour is supposed to be linear even if some researchers [START_REF] Jang | Indirect measurement of the impulsive load to a nonlinear system from dynamic responses: Inverse problem formulation[END_REF] tried to study nonlinear systems.

However, all the methods necessitate knowing the force location [START_REF] Doyle | An experimental method for determining the location and time of initiation of an unknown dispersing pulse[END_REF]- [START_REF] Ahmari | Impact location and load identification through inverse analysis with bounded uncertain measurements[END_REF]. So a lot of works must be done to estimate this location first. Wang et al. [START_REF] Wang | Sensor placement methods for an improved force identification in state space[END_REF] studied the influence of the sensor location on the force identification. Jacquelin et al. [START_REF] Jacquelin | Analysis of a force reconstruction problem[END_REF] suggested that the measurement location has a direct influence on the condition number of the transfer matrix. Samagassi et al. [START_REF] Samagassi | Reconstruction of multiple impact forces by wavelet relevance vector machine approach[END_REF] proposed a method of location through the identification of impact multiple forces. In [START_REF] Samagassi | Reconstruction of multiple impact forces by wavelet relevance vector machine approach[END_REF], reconstruction of forces generated by multiple impacts occurring on linear elastic structures has been achieved through wavelet relevance vector machine approach to inverse problem solution. Bayesian hierarchical modeling according to the relevance vector machine approximation was then applied in order to estimate the forces generated by impact events. The obtained results were remarkably good as the reconstructed forces were found to be very close to the original forces at the system input. The method was found also to provide a way for localization of impact forces. Almost all the works are related to identify one force, except [START_REF] Samagassi | Reconstruction of multiple impact forces by wavelet relevance vector machine approach[END_REF]. Sometimes a distributed force is studied [START_REF] Liu | An improved method for the reconstruction of a distributed force acting on a vibrating structure[END_REF] but the problem came down to identify one force. Multiple-force identification was addressed in [START_REF] Samagassi | Reconstruction of multiple impact forces by wavelet relevance vector machine approach[END_REF], [START_REF] Robert | Multiple force identification for complex structures[END_REF], [START_REF] Sekine | Identification of locations and force histories of multiple point impacts on composite isogrid-stiffened panels[END_REF] and [START_REF] Martin | Impact force location in frame structures[END_REF]: in the latter two publications they obtained very good results without performing any regularization.

Many of these previous studies demonstrate that more sensors are required to determine the impact-load location. In fact, more sensors can also improve the identification accuracy of impact-force. Thus, some researchers have worked on the possibility of identifying impact-forces with fewer sensors than there are forces [START_REF] Tan | Efficient Sparse Bayesian Learning via Gibbs Sampling[END_REF][START_REF] Ji | Bayesian Compressive Sensing[END_REF][START_REF] Cao | Bayesian signal detection with compressed measurements[END_REF][START_REF] Khoo | Impact force identification with pseudo-inverse method on a lightweight plate for under-determined, even-determined and over-determined cases[END_REF][START_REF] Qiao | Impact-force sparse reconstruction from highly incomplete and inaccurate measurements[END_REF][START_REF] Qiao | Sparse regularization for force identification using dictionaries[END_REF]. Known as compressed sensing (CS) theory, Candès and Tao wrote the fundaments of this theory [START_REF] Durand | Matrices aléatoires et norme L1 pour le compressed sensing[END_REF]. The CS considers the sparse character which exists in most of physical data. In particular reconstructing data with much lower sample rate than classical Shannon limit is possible in the CS framework.

In structural dynamics, data acquisition may require several sensors. In situations when the number of mea-surements is limited [START_REF] Earls | Bayesian inference of hidden corrosion in steel bridge connections: Non-contact and sparse contact approaches[END_REF], it is still possible to identify a force in a CS context by considering the sparsity of the force. In many situations, the signal-force is not directly sparse : however, in many cases, it is possible to find a basis in which the force is approximately sparse [START_REF] Tan | Efficient Sparse Bayesian Learning via Gibbs Sampling[END_REF][START_REF] Ji | Bayesian Compressive Sensing[END_REF][START_REF] Hu | Compressive sampling in computed tomography: Method and application[END_REF].

Identifying forces is useful for determining the amount of damage undergone by a structure after a shock event and enables to evaluate through reanalysis of the structure its residual capacity. When the structure is represented by a discrete model, having the general form of a Toeplitz like matrix, the information acquired by sensors which are implemented on the structure does not allow by simple inversion of that matrix to straightforwardly recovering the input forces. On one hand, the problem is habitually ill-conditioned and may be even ill-posed, and on the other hand the information acquired by sensors as well as the Toeplitz matrix can be noisy and/or uncertain. To tackle these issues, which will be considered in this paper, Bayesian approach appeared to be an interesting way of achieving regularization of the problem [START_REF] Zhang | Bayesian force reconstruction with an uncertain model[END_REF][START_REF] Samagassi | Reconstruction of multiple impact forces by wavelet relevance vector machine approach[END_REF].

In this paper, the force identification problem in time domain was investigated. On account of the presence of measurement noise and modeling error, the impact-force identification is tackled within a Bayesian framework. The Bayesian approach considers the unknown quantities of interest as random variables [START_REF] Zhang | Bayesian force reconstruction with an uncertain model[END_REF][START_REF] Samagassi | Reconstruction of multiple impact forces by wavelet relevance vector machine approach[END_REF][START_REF] Lee | Bayesian Statistics : An Introduction[END_REF][START_REF] Tarantola | Inverse Problem Theory and Methods for Model Parameter Estimation[END_REF]. This approach has several advantages. First, it endows the unknown force with prior information in the form of a probability density function (p.d.f.), which naturally imposes an intrinsic regularization. Second, the Bayesian approach provides a rigorous probabilistic framework to account for all possible sources of errors, which participate to the uncertainty of the identified forces including the measurement noise and the model uncertainty.

The main objective of this paper is to study the possibility for the Bayesian approach together with the so-called compressed sensing technique [START_REF] Candés | An introduction to compressive sampling[END_REF][START_REF] Candès | Sparsity and incoherence in compressive sampling[END_REF] to identify multiple forces that act on a linear elastic structure, through indirect measurements. This Bayesian identification of multiple forces is analyzed through the Gibbs sampler.

The rest of the paper is organized as follows. The general description of force reconstruction problem is tackled in section 2. In Section 3, two strategies of regularization including the l 2 -norm regularization and l 1 -norm regularization are illustrated and a Bayesian sparse regularization approach is proposed as sparse deconvolution model of impact-force identification. Section 4 highlights the technique of Bayesian sparse regularization approach for identifying impact-force. Before concluding in Section 6, numerical simulations will be make on a planar beam modeled both by an analytical approach and by the finite element method in Section 5.

Obviously, the load identification issue is not new but it seems that the limitations of such an identification must be highlighted. that's why the identification of force will also be carried out by the classical approach through the regularization of Tikhonov in order to make a comparative study of the two approaches. The quality of the results obtained by the regularization of Tikhonov seems to be related to the choice of the regularization parameter. Several methods exist to reach a "good choice" of the regularization parameter [START_REF] Hansen | A matlab package for analysis and solution of discrete ill posed problems[END_REF] : the L-curve method, the generalized cross validation criterion (GCV), the quasi-optimality criterion and the Morozov's discrepancy principle are among the most used criteria. In the following, the GCV method is used to provide the identified actions because the criterion GCV seems to be the one that would lead to better results [START_REF] Jacquelin | Force reconstruction analysis and regularization of a deconvolution problem[END_REF][START_REF] Lin | A new investigation into regularization techniques for the method of fundamental solutions[END_REF]. However, the L-curve and the quasi-optimality were used to compare the regularization parameters.

General description of the force reconstruction problem

In the following, the structure is supposed to have a linear elastic behaviour and to remain elastic during and after applying forces.

Let G ij (t) be the impulse response between degrees-of-freedom (dof) i and j. This function may be determined thanks to a modal expansion of the response, either analytically (for very simple structures) or from a finite element model. Suppose that n F forces act on the beam along several dofs, and that n m measurements are carried out.

Then response s i l (t) of dof i l is given by the convolution equation:

∀l ∈ 1, • • • , n m , s i l (t) = n F ∑ k=1 ∫ t 0 G i l j k (t -τ ) F j k (τ ) dt) = n F ∑ k=1 G i l ,j k (t) ⋆ F j k (t) (1)
where {j k } k=1•••n F is the set of the dofs along which a force is applied; similarly, {i l } l=1•••nm is the set of the dofs along which a measurement is performed; so index l (resp. k) of i l (resp. j k ) restricts the dof set to the dofs where a signal (resp. a force) is measured (is applied); "⋆" stands for the convolution product.

Convolution equation ( 1) is discretized. This leads to a system of algebraic equations ( see appendix C for more details about the identification of multiple forces):

∀l ∈ 1, • • • , n m , s i l = n F ∑ k=1 G i l ,j k F j k (2)
where:

-∆t is the time step; n t is the number of time sampling, that is to say, the number of time steps used for sampling the continuous response measured by the sensor.

-G i l ,j k is a n t -by-n t transfer matrix:

G i l ,j k = ∆t             G i l ,j k (∆t) 0 0 G i l ,j k (2∆t) G i l ,j k (∆t) . . . G i l ,j k (3∆t) G i l ,j k (2∆t) . . . . . . . . . . . . . . . . . . 0 G i l ,j k (n t ∆t) G i l ,j k ((n t -1)∆t) . . . . . . G i l ,j k (∆t)             -s i l = [s i l (∆t), . . . , s i l (n t ∆t)] t , -F j k = [F j k (0), . . . , F j k ((n t -1)∆t)] t ,
Equations (2) lead to the following problem:

s = G F (3)
where

s = [s i 1 t , . . . , s in m t ] t ; F = [F j 1 t , . . . , F jn F t ] t is the force vector; G is a Toeplitz block G kl matrix such as G kl = G i l ,j k for k = 1 • • • n F and l = 1 • • • n m ; so G is a m-by-n matrix, s is a m-vector, F is a n-vector, where m = n t × n m and n = n t × n F .
However, given the observed response s of structure after the impact events and the transfer matrix G which describes completely the structure 's dynamic characteristic, solving Eq.( 3) for the unknown impct-force vector F is called impact-force deconvolution. An important issue in practice is that actual observations always contain some amount of noise. Noise may arise in particular because of numerical round-off or nonlinearities that have not been taken into account for instrument readings. Taking into account the measurement noise, Eq.(3) becomes :

s = G F + η (4)
where η, m-vector, is a random vector representing measurement noise.

If we consider impact-force signal F in Eqs.( 3) or ( 4) and an associated signal w which is r-sparse in an appropriate basis B, that which means that w contains no more than r non-zero elements or significant components with typically r ≪ n. This property can be expressed under vector-matrix form by a projection taking the following notation :

F = Bw (5)
where w ∈ R n×1 , and B ∈ R n×n . The columns of B define n vectors representing an orthonormal basis.

For most natural force signals F that are piecewise smooth, the wavelet transform such as Daubechies and Symlets wavelets transform has been shown to yield sparse representation [START_REF] Daubechies | Ten Lectures on Wavelets[END_REF][START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF].

B is an orthogonal basis matrix, so B satisfies B t B = BB t = I. Thus, Eq.( 5) can be easily inverted :

w = B t F (6)
where () t represent the transpose operation. By using Eq.(5), Eq.( 4) rewrites

s = GBw + η (7)
or equivalently

s = Aw + η ( 8
)
where A = GB ∈ R m×n .

We can divided the inverse problem of impact-force identification into three categories depending on the total number of unknown impact-sources F and known responses s -Case 1:

Eq.( 8) is a under-determined case : the number of measurements m is less than the number of sources n.

-Case 2:

Eq.( 8) is a even-determined case: the number of measurements m is equal to the number of sources n.

-Case 3:

Eq.( 8) is a over-determined case: the number of measurements m is greater than the number of sources n.

Solving Eq.( 8) in case 1 for the unknown vector w given the observed response s as well as matrix A is a compressed sensing problem [START_REF] Candés | An introduction to compressive sampling[END_REF][START_REF] Candès | Sparsity and incoherence in compressive sampling[END_REF][START_REF] Candès | Decoding by linear programming[END_REF][START_REF] Candès | Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information[END_REF]. Compressed sensing (CS) appears,so, as a simple method for finding the sparsest solution to an under-determined system of linear equations having the form of Eqs. [START_REF] Hansen | Deconvolution and regularization with toeplitz matrices[END_REF] or [START_REF] Hansen | Rank-Deficient and Discrete III-Posed Problems[END_REF]. Also, CS can be see as an effective solution approach particularly in situations where the number of sensors is limited because of implementation constraints or cost and when the sensing process provides a small number of measurements due to a poor sampling rate [START_REF] Daubechies | Ten Lectures on Wavelets[END_REF].

The first two cases have been considered in this paper.

Note that for real systems, the matrix G can also be polluted. There are two ways in order to take into account system matrix uncertainties. In the first approach, one uses model uncertainty propagation based approaches that enable to estimate variability of G following those of the basic design parameters. In the second approach, one apply posterior noise on the matrix G in order to take into account globally and non-parametrically the various physical perturbations, modeling errors as well as material and geometrical variability.

The method based on propagation gives raise to major difficulties and is out of the scope of the present work. Thus, in the following, only the second method is considered and the matrix of system is assumed to be perturbed by a general white Gaussian noise.

Note also that the problem of force identification is naturally ill-posed. The condition number of the transfert matrix G can be very large, which means that the desired solution F through w is very sensitive to the small noise in response s. The direct Least squares (LS) approach that minimize the l 2 -norm of the response residual ∥s -Aw∥ 2 in order to solve Eq.( 8) never provides a satisfactory solution. Thus in order to stabilize the problem, regularization strategies must be developed to reconstruct the impact-force .

3. Strategies of regularization for impact-force reconstruction

Method of regularization based on l 2 -norm

If there is a method of regularization that is most commonly used and well-known for solving various inverse problems, it is probably Tikhonov regularization [START_REF] Hansen | Rank-Deficient and Discrete III-Posed Problems[END_REF][START_REF] Lin | A new investigation into regularization techniques for the method of fundamental solutions[END_REF]. The primary purpose of this regularization is to find a stable solution by minimizing the weighted combination of the residual norm and the solution norm :

minimize ∥s -Aw∥ 2 2 + β∥w∥ 2 2 (9)
where β is a regularization parameter and ∥.∥ 2 is the l 2 -norm. With this term β∥w∥ 2 2 , the ill-posed inverse problem becomes well-posed. A fair balance between minimizing ∥s -Aw| 2 and ∥w∥ 2 is established by β.

Regularization of Tikhonov always has an explicit and unique solution to Eq. ( 9), for any β fixed :

w tik = ( A t A + βI ) -1 A t s ( 10
)
where I is an identity matrix. The quality of the results obtained by this approach seems to be related to the good choice of the regularization parameter β. Several methods exist to reach a "good choice" of the regularization parameter [START_REF] Hansen | Rank-Deficient and Discrete III-Posed Problems[END_REF][START_REF] Lin | A new investigation into regularization techniques for the method of fundamental solutions[END_REF]: the L-curve method, the generalized cross validation criterion, the quasi-optimality criterion and the Morozov's discrepancy principle are among the most used criteria. In the following, the GCV method is used to provide the identified actions. However, the L-curve and the quasi-optimality were used to compare the regularization parameters.

Method of regularization based on l 1 -norm

It is well-known the classical regularization methods based on l2-norm such as Tikhonov regularization fail to solve the case 1, that is to say under-determined case. Inversely, the l1-norm regularization, that is to say sparse regularization, refers to finding the minimum l1-norm solution to an under-determined system.

Recently it has received much attention, particularly motivated by compressive sensing theory under sparsity condition, where the minimum l1-norm solution is also the sparsest solution [START_REF] Candés | An introduction to compressive sampling[END_REF][START_REF] Candès | Sparsity and incoherence in compressive sampling[END_REF][START_REF] Candès | Decoding by linear programming[END_REF][START_REF] Candès | Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information[END_REF]. In fact, when the desired impact-force is known to be sparse in the sense that most components of w are close to zeros, more measurements will not be necessary. Therefore, the sparsity nature of impact-force w fundamentally changes the under-determined problem, making, thus, unique solution possible. The sparse regularized problem associated to the problem Eq.( 8) [START_REF] Samagassi | Reconstruction of multiple impact forces by wavelet relevance vector machine approach[END_REF][START_REF] Qiao | Impact-force sparse reconstruction from highly incomplete and inaccurate measurements[END_REF][START_REF] Qiao | Sparse regularization for force identification using dictionaries[END_REF] can be written as :

minimum { ∥s -Aw| 2 2 + β∥w∥ 1 } ( 11
)
where β is also called parameter of regularization. The sparse solution of Eq.( 8) is constructed by l 1norm term. The problem of optimization in Eq.( 11) is known as the basis pursuit denoising (BPDN)

problem [START_REF] Qiao | Impact-force sparse reconstruction from highly incomplete and inaccurate measurements[END_REF]. Note that the inherent sparsity of impact-event is the prerequisite of sparse regularization.

While the solution of Tikhonov can be explicitly expressed, sparse regularization requires solving the convex optimization problem Eq.( 11) by an iterative approach, that which has no explicit analytic solution. There are much theoretical work in compressing sensing domain that have shown that the so-called restricted isometry property (RIP) of the matrix A guarantees that sparse regularization will yield an accurate and robust solution from the incomplete and noisy measurement s [START_REF] Candés | An introduction to compressive sampling[END_REF][START_REF] Candès | Sparsity and incoherence in compressive sampling[END_REF][START_REF] Candès | Decoding by linear programming[END_REF][START_REF] Candès | Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information[END_REF]. If the matrix A of system verifies the RIP, any small perturbation of measurement resulting from noise will yield a small perturbation on the identified sparse signal. It has been shown in [START_REF] Samagassi | Reconstruction of multiple impact forces by wavelet relevance vector machine approach[END_REF][START_REF] Qiao | Impact-force sparse reconstruction from highly incomplete and inaccurate measurements[END_REF] that the RIP is a sufficient condition, that which means that even though a matrix does not satisfies the restricted isometry property (RIP), it can yet be used as a matrix of measurement. Results obtained in [START_REF] Samagassi | Reconstruction of multiple impact forces by wavelet relevance vector machine approach[END_REF][START_REF] Qiao | Impact-force sparse reconstruction from highly incomplete and inaccurate measurements[END_REF] are very satisfactory.

For our particular problem of force reconstruction, it has not still been proven mathematically that the matrix A satisfy the RIP. Nevertheless, some studies have shown that the verification of the RIP strongly resides in the random nature of the matrix of the system [START_REF] Candès | Sparsity and incoherence in compressive sampling[END_REF][START_REF] Candès | Decoding by linear programming[END_REF]. This is why, we perturbed by a Gaussian term the matrix of the system in order to take into account globally and non-parametrically the various physical perturbations, modeling errors as well as material and geometrical variability. This uncertainty which is propagated to the response is assumed small enough : A + δA. So, robustness can only be verified numerically.

In this paper, iterative approach proposed in order to solve Eq.( 11) is Bayesian hierarchical model. Two algorithms based on Gibbs sampling are proposed.

4. Bayesian sparse regularization for impact-force identification

Bayes'formula

Denoting by x, the vector including the unknown parameters (F (or w), noise η, model uncertainty), the Bayesian solution of equation Eq. [START_REF] Jacquelin | Force reconstruction analysis and regularization of a deconvolution problem[END_REF] or Eq.( 8) is obtained by the a posteriori probability density function (pdf) of x given by the Bayes' formula:

π (x|s, J) = π (s|x, J) π (x, J) π (s, J) ∝ π (s|x, J) π (x, J) (12)
Information J, which will be more detailed later, is introduced to make information a priori explanatory.

The elements of equation ( 12) are listed as follows:

the likelihood function, π (s|x, J), reflects the probability of observing the data s given a set of parameters x and J;

the prior pdf, π (x, J), reflects our knowledge on the parameters before experiments are undertaken.

the marginal function, π (s, J) = ∫ π (s|x, J) π (x, J) dx, is often seen as a normalization constant.

Eq.( 12) provides the updated information on the parameters involved in the inference. The force vector, w, is finally estimated with an iterative process, the Markov chain Monte Carlo, through a probability density function (pdf).

The construction of the likelihood requires two steps. The first one is the definition of the measurement noise;

the second one is the representation of the noise by a probability density function. Reference [START_REF] Kaipio | Statistical and computational inverse problems[END_REF] shows that the pdf of the noise, π noise (η), and the pdf of the measurement s are identical because of the mutual independence between F (or w) and η . In the following, the noise that spoils the data, s, is supposed to be a Gaussian random noise with zero mean: η ∼ N (0, Γ noise ), where Γ noise = σ 2 η I is the covariance matrix and I is the identity matrix.

First algorithm for Bayesian force identification

One of the first impact-force reconstruction algorithms based on Bayesian approach was developed by Zhang and al. [START_REF] Zhang | Bayesian force reconstruction with an uncertain model[END_REF]. This impact-force reconstruction has been discussed in the frequency domain.

In this subsection, impact-force identification is discussed in time domain and we explain a Bayesian iterative method to solve convex optimization problem expressed by Eq.( 11). We consider situation where matrix system is unperturbed, ie A = GB. Thus, the likelihood function reads:

π ( s|w, σ 2 η ) = π noise ( s -Aw|w, σ 2 η ) (13)
This means that: s|w, σ 2 η ∼ N (Aw, Γ noise ).

A conjugate prior1 pdf is attributed to force w :

w ∼ N (w 0 , Γ pr ) (14)
with w 0 the vector of mean values and the covariance matrix Γ pr is assumed to be diagonal:

Γ pr = diag(σ 2 w )
where

σ 2 w = [ σ 2 w (1), σ 2 w (2), . . . , σ 2 w (n t )
] .

In Eqs. ( 13) and ( 14), parameters w 0 , σ 2 w (i), and σ 2 η , called hyper-parameters are unknown. These hyperparameters may be viewed as regularization parameters.

In addition, these hyper-parameters are the parameters of the prior distribution and noise. So they are a source of a priori information [START_REF] Zhang | Bayesian force reconstruction with an uncertain model[END_REF] : w 0 , σ 2 w (i), and σ 2 η are hyper-parameters and are random variables following probability distributions :

w 0 ∼ N (U 0 , C u ) (15) σ -2 η ∼ Γ (k η , β η ) (16) σ -2 w (i) ∼ Γ (k w , β w ) (17)
where Γ (k, β) stands for the gamma distribution and C u = σ 2 u I.

Denote J = { U 0 , σ 2 u , k w , β w , k η , β η }
the set of parameters supposed to be known a priori; then the posterior probability density is:

π ( w, w 0 , σ -2 w (i), σ -2 η |s ) ∝ π noise ( s -Aw|w, σ -2 η ) π pr ( w|w 0 , σ -2 w (i) ) π pr ( w 0 , σ -2 w (i), σ -2 η |J ) (18)
Eq.( 18) is the Bayesian solution of the force identification problem. The determination of inference parameters w, w 0 , σ -2 w (i), σ -2 η can be achieved through the MCMC Gibbs algorithm (Algorithm 1 below) : The maximization of this posterior probability density with respect to w yields the posterior maximum estimator [START_REF] Kaipio | Statistical and computational inverse problems[END_REF] :

w = ( Γ -1 pr + A T Γ -1 noise A ) -1 ( A T Γ -1 noise s + Γ -1 pr w 0 ) (19)
A comparison between the Bayesian approach and classical Tikhonov regularization is given in Appendix A.

The principle of MCMC consists in generating random walks whose stationary pdf coincides with the posterior pdf of interest. Among the MCMC methods [START_REF] Hastings | Monte Carlo sampling methods using Markov chains and their applications[END_REF][START_REF] Geman | Stochastic relaxation, Gibbs distributions and the Bayesian restoration of images[END_REF][START_REF] Gilks | Markov Chain Monte Carlo in Practice[END_REF], the Gibbs sampling algorithm [START_REF] Geman | Stochastic relaxation, Gibbs distributions and the Bayesian restoration of images[END_REF] is used to explore the posterior pdf. The implementation of the Gibbs sampling on force reconstruction is resumed as follows (See Appendix B for details).

-Algorithm 1

Initialize parameters of J = { U 0 , σ 2 u , k w , β w , k η , β η } 1. Draw w 0 ∼ N ( U 0 , σ 2 u I ) σ -2 w (i) ∼ Γ (k w , β w ) σ -2 η (i) ∼ Γ (k η , β η ) 2. Draw w|s, w 0 , σ -2 w (i), σ -2 η ∼ N ( w, Γ post ) with w = ( Γ -1 pr + A T Γ -1 noise A ) -1 ( A T Γ -1 noise s + Γ -1 pr w 0 ) Γ post = ( Γ -1 pr + A T Γ -1 noise A ) -1 3. Draw σ -2 η |w, s ∼ Γ ( kη , βη ) with kη = k η + nt 2 βη = ∥s-Aw∥ 2 2 2 + β η 4. Draw w 0 |w, σ -2 w (i) ∼ N ( Û0 , Ĉu ) with Ĉu = ( Γ -1 pr + σ -2 u I ) -1 Û0 = Ĉu ( Γ -1 pr w + σ -2 u U 0 ) 5. Draw σ -2 w (i) |w, w 0 ∼ Γ ( kw , βw ) with kw = k w + nt 2 βw = ∥w-w 0 ∥ 2 2 2
+ β w 6. go to step 2 and repeat until a sufficiently large sample is collected after the burn-in phase.

In practice the initial values given to J are assigned arbitrarily. However, as the algorithm is iterative, these parameters are updated in the loop and their final values does not depend on their arbitrary initial value.

The Markov chain needs a "heating" stage to reach the stationary target p.d.f; this stage is the so-called burn-in phase, in which the samples are not yet drawn from the target probability distribution. Thus, the first few samples are not taken into account when estimating unknown parameters (ie, w, w 0 , σ 2 η , σ 2 w (i))

of the inference. The estimate of these unknown is made by calculating the average for each unknown parameter without the samples drawn in the burn-in stage. The convergence of a Markov chain can be roughly verified by inspecting its trajectory.

Sparse Bayesian learning via Gibbs sampling for impact-force identification

From the Bayesian point of view, Eq.( 11) is equivalent to the maximum a posteriori method for identifying w, considering that η is a normal random vector and w has a Laplacian a priori distribution. For solving Eq.( 11) many approaches have recently been proposed [START_REF] Tan | Efficient Sparse Bayesian Learning via Gibbs Sampling[END_REF][START_REF] Qiao | Impact-force sparse reconstruction from highly incomplete and inaccurate measurements[END_REF][START_REF] Qiao | Sparse regularization for force identification using dictionaries[END_REF][START_REF] Tan | Compressed sensing via sparse Baysian and Gibbs sampling[END_REF]. It is proposed, here, a novel Gibbs sampling based sparse Bayesian learning (SBL) method in order to solve Eq.( 11).

As mentionned above, for real systems the matrix A can be polluted by applying posterior noise on the matrix G in order to take into account globally and non-parametrically the various physical perturbations, modeling errors as well as material and geometrical variability. Thus, matrix A can be expressed by A =

(G + δG) B.
To solve the problem defined in Eq.( 8) or Eq.( 11), consider the following Bayesian model :

1. the likelihood function is

s|w, σ 2 η ∼ N ( Aw, σ 2 η I ) ,
2. the prior informations are modeled by:

(a) w|σ 2 w ∼ N (0, Γ pr ) (b) σ 2 w (i) ∼ IG (k w , β w ) (c) π ( σ 2 η ) ∝ 1
where

Γ pr = diag ( σ 2 w (i) ) , σ 2 w = [ σ 2 w (1) , σ 2 w (2) , • • • , σ 2 w (n) ]
. The sparsity of w relies heavily in the sparsity of σ 2 w . We consider the inverse gamma prior for σ 2 w (i) because it is known that σ 2 w is sparse and the inverse gamma prior promotes sparsity in the estimate of σ 2 w , when β w is small [START_REF] Tan | Efficient Sparse Bayesian Learning via Gibbs Sampling[END_REF][START_REF] Tan | Compressed sensing via sparse Baysian and Gibbs sampling[END_REF]. No a priori information on noise σ 2 η is available, that is why its prior distribution, π

( σ 2 η )
, is assumed to be uniform. So, conditional distributions can readily be calculated :

the joint posterior pdf:

π ( w, σ 2 w , σ 2 η |s, k w , β w ) = π ( s|w, σ 2 η ) π ( w|σ 2 w ) π ( σ 2 w |k w , β w ) -the posterior pdf: π ( w|s, σ 2 w , σ 2 η , k w , β w ) = π ( s|w, σ 2 η ) π ( w|σ 2 w ) ∝ e -1 2 ((w-µ) t Σ -1 (w-µ))
where Σ = ( (

σ 2 η ) -1 A T A + Γ -1 pr ) -1 and µ = ( σ 2 η ) -1 ΣA T s
updating the noise:

σ 2 η |s, w, k w , β w ∼ IG ( m 2 -1, ∥s -Aw∥ 2 2 2 )
updating the variance:

σ 2 w (i) |s, w, k w , β w ∼ IG ( k w + 1 2 , β w + w 2 i 2 )
Note (.) T represent transpose operation. Using these relationships, we can derive the Gibbs sampling algorithm as follow :

-Algorithm 2

Initialize w, k w and β w Perform the following steps at the t-th iteration :

1. Draw ( σ 2 w (i) ) t |s, w (t-1) , k w , β w ∼ IG ( k w + 1 2 , β w + ( w (t-1) i ) 2
2

)

with i = 1, . . . , n. 2. Draw ( σ 2 η ) t |s, w (t-1) , k w , β w ∼ IG ( m 2 -1, ∥s-Aw (t-1) ∥ 2 2 2 ) 3. Draw (w) t |s, ( σ 2 w ) t , ( σ 2 η ) t , k w , β w ∼ N ( µ t , Σ t )
where

Σ t = [ ( ( σ 2 η ) -1 ) t A T A + ( Γ -1 pr ) t ] -1 Σ t = (Γ pr ) t -(Γ pr ) t A T ( ( σ 2 η ) t I + A (Γ pr ) t A T ) -1 A (Γ pr ) t and µ t = ( ( σ 2 η ) -1 ) t Σ t A T s µ t = (Γ pr ) t A T ( ( σ 2 η ) t I + A (Γ pr ) t A T ) -1 s
4. go to step 1 and repeat until a sufficiently large sample is collected after the burn-in phase.

We recommend CoSaMP developped by Needell and al. [START_REF] Needell | CoSaMP: Iterative signal recovery from incomplete and inaccurate samples[END_REF] in order to obtain the initial estimate of w.

CoSaMP requires the knowledge of the level of sparsity of w, that is to say, the number of non-zeros elements in w (i.e r), which is unknown a priori in most applications. Since r is unknown, we simply suggest to set it to m/2.

An affine transformation of a normal random vector z can be used to generate the normal random vector (w) t :

(w) t = µ t + (ϕ) t z where z ∼ N (0, I) and (ϕ) t ( (ϕ) t ) T = Σ t .
Conventionally, ( (ϕ) t ) can be obtain by the Cholesky decomposition. However, the decomposition of Cholesky is of a complexity of O ( n 3 ) [START_REF] Tan | Efficient Sparse Bayesian Learning via Gibbs Sampling[END_REF][START_REF] Tan | Compressed sensing via sparse Baysian and Gibbs sampling[END_REF], and so it is not suitable for large-scale applications. Instead, we illustrate an efficient way to generate (w) t as follows.

Consider z 1 and z 2 two independent normal random vectors of dimensions m and n, respectively. Then :

z 3 = ( ( σ 2 η ) t ) -1 2 A T z 1 + ( (Γ pr ) t ) -1 2 z 2 (20)
is a normal random vector with 0 as mean and

( Σ t
) -1 as matrix of covariance. The matrix vector product A T z 1 determine the complexity of Eq.( 20). Next, it can be verify that

(w) t = µ t + Σ t z 3 (21)
is a normal random vector where µ t is mean and Σ t matrix of covariance . Using these precedent formulas

Σ t = (Γ pr ) t -(Γ pr ) t A T ( ( σ 2 η ) t I + A (Γ pr ) t A T ) -1 A (Γ pr ) t and µ t = (Γ pr ) t A T ( ( σ 2 η ) t I + A (Γ pr ) t A T ) -1 s
Eq.( 21) can be simplify as follow :

(w) t = (Γ pr ) t z 3 -(Γ pr ) t A T ( ( σ 2 η ) t I + A (Γ pr ) t A T ) -1 E ( 22
) with E = s -A (Γ pr ) t z 3 .
Conjugate gradient (CG) method can be used to compute

( ( σ 2 η ) t I + A (Γ pr ) t A T ) -1 E.
Thus, we can identify the desired impact-force F through Eq.( 5). Algorithm 2 being a Markov chain, it needs a burn-in phase to reach the stationary target p.d.f: the first few samples are not taken into account to estimate the pdf. The estimate of the unknowns is done by calculating the average for each unknown parameter. The convergence of a Markov chain can be roughly verified by inspecting its trajectory, or more reliably by several repetition of simulated Markov chains that should converge to the same probability distribution.

Numerical simulations

The objectives of this section are twofold: firstly the Bayesian probabilistic method is used and its capacity to identify a force that may be identically equal to zero is assessed, and secondly the difficulties are identified and highlighted. These objectives will be illustrated by testing Algorithms 1 and 2.

A pinned-pinned beam is studied to illustrate the identification process proposed in this paper. Two cases are analyzed: in the first application, the beam is analytically modeled whereas in the second application a finite element model of the beam is used. The beam is made of a homogeneous linear elastic material having a uniform rectangular cross section. The geometrical characteristics of the beam are shown in Fig. 2 below.

Forces and moments can be applied to the beam: in the following, an action will denote either a force or a moment. All the actions applied in the following are proportional to:

f (t) = f 0 t α e -βt ( 23
)
where in example 1 f 0 = 2.10 15 N/s 6 , α = 6 and β = 300 s -1 , and f 0 = 2.10 19 N/s 6 , α = 6 and β = 1000 s -1 in example 2. In order to vary the actions applied, the force reconstructed by the regularization of Tihonov has for parameters : f 0 = 3.10 11 N/s 6 , α = 3 and β = 1300 s -1 . This function is depicted in Fig. 1.

An action that is not identically equal to zero will be referred to a non-zero action (NZA); an action that identically equals zero is an identically zero action (IZA). A delay may be introduced in the action to test its influence on the force identification.

To validate the quality of the identification of the NZA, a relative error between the actual action ( A true )

and the identified action ( A id ) is used:

E N ZA (%) = 100 × ∥A true -A id ∥ 2 ∥A true ∥ 2 (24)
For IZA, the error is defined as follows:

E IZA (%) = 100 × ∥A id ∥ 2 ∥A N ZA ∥ 2 (25)
The noise, η, is Gaussian with zero mean and standard deviation proportional to the maximum amplitude of the response. The standard deviation is defined in each of the applications below. The total number of iterations is 10,500. The number of iterations that constitute the burn-in phase is 500: the 500 samples drawn in this phase are not taken into account.

It is important to remember the importance of identifying a null action. The objective of this manuscript is to locate the impact zones through the reconstruction of the impact forces that the structure has undergone.

For a finite element model such as ours, assuming priori that the load to be identified is distributed over all the nodes of the studied structure, one eliminates the problem of localization. Thus, where the identified load will be (almost) zero, we can then conclude that there was no loading applied: it is then possible to locate the loading zone posteriori. The quality of the result is therefore directly related to the ability of the methods used to identify a constantly effort equal to zero and to simultaneously identify multiple actions. Although the method can be applied in general to any elastic structure for which the matrix G in Eq.( 4) is available, the two Gibbs's algorithms will be tested on a planar beam having a symmetric section and loaded orthogonally to its mean fiber, in its plane of symmetry. Fig. (3) depicts the considered beam which is assumed to be made from a homogeneous elastic linear material and having a uniform cross section. The beam section is assumed to be rectangular. This enables to derive simple analytical time transfer functions between any excitation point where a transverse force f (x i , t) is applied and any response point where the longitudinal normal strain s (x j , t) is assumed to be measured. The measurement point is assumed to be located on the upper beam fiber where a strain gauge sensor is placed.

It is straightforward to derive for the simply supported beam considered here that the Toeplitz like matrix G (x i , x j ) giving the discrete time response in terms of longitudinal axial strain of the upper fiber s as function of the discrete force vector F according to Eq.( 3), writes :

G (x i , x j ) = π 2 h ρAL 3 N modes ∑ n=1 n 2 sin ( nπx i L ) sin ( nπx j L ) g (w n , ξ n , (k -l)∆t) ( 26
)
with g (w n , ξ n , (k -l)∆T ) = 0 if k ≺ l, and

g (w n , ξ n , (k -l)∆t) = ∆t exp (-δ n w n (k -l)∆t) w n √ 1 -ξ 2 n sin ( w n √ 1 -ξ 2 n (k -l)∆t ) if k ≥ l ( 27
)
and the rotational frequencies are

w n = n 2 π 2 h 2L 2 √ E 3ρ (28)
where k, l ∈ 1, 2, • • • , N t , E is the Young's modulus, ρ the density, h the beam thickness, L the beam length, ξ n the modal damping associated to mode number n, ∆t the time step and Nmodes the total number of modes selected owing to modal truncation. N modes is determined from a convergence study or can be estimated by considering the highest frequency contained in the impact signal. The time step is chosen in order to satisfy Shannon sampling condition with regards to the maximum frequency contained in the excitation, even if this condition is not necessary.

In this first application, the identification of a single force is considered. The force is applied at s, where ∆t max = 1/(2f 5 ). The applied force spectrum is taken with a maximum frequency which is smaller than f 5 .

x i = L/
The results provided by Algorithm 1 are plotted in Figs. 4 to 6. The values to initialize algorithm 1 are

J = { U 0 = [0], σ 2 u = 100, k w = 10, β w = 15, k η = 10, β η = 0.1 } .
Three noise standard deviations were tested: 0.01 %, 0.1%, and 1% of the response maximum amplitude.

Among these results, only those of Fig. 4 are excellent even if those of Fig. 5, affected by parasitic oscillations, are still acceptable. Note that the delay does not really improve the reconstruction: for a given figure, all errors are almost the same whatever the delay. This means that Algorithm 1 is suitable for noise low level.

Algorithm 2, initialized with k w = 1 and β w = 10 -5 , was also tested with noise standard deviation equal to 1 % (Fig. 7 Note also that the delay in signal start contributed to improve slightly the results.

Thus, Algorithm 2 appears to be more effective compared to Algorithm 1. Indeed, the results obtained are excellent until a noise level about 2 %, which highlights its robustness with respect to Algorithm 1. So, in the following, Algorithm 2 is used to identify all the actions.

Multiple-force identification

The pinned-pinned beam is discretized in n e Bernoulli beam finite elements and then has 2 × n e degrees-offreedom (dof): the dofs are related to translations and rotations. An action is applied along dof i F direction, and a direct problem is solved: then some dofs are numerically measured. These measurements are spoilt by a random noise whose standard deviation equals to 2 % of the response maximum. The forces along dof i F and along another dof are identified.

The objective is then two-fold: testing the ability to identify two actions, and to assess the identification of an identically zero action (IZA). Several cases will be addressed depending on the nature of the action (force or moment) and the measurements (number, translation or rotation).

Identically zero action identification

The beam is discretized in four elements (see Fig. 8); the non zero action (NZA) is applied along one dof (translation or rotation) and must be identified; an IZA along another dof (translation or rotation) is 9: the dof responses related to case 1 were used. The NZA is very well identified. Algorithm 2 is also efficient to identify IZA: the maximum of F 2 is about 25,000 times smaller than the maximum of F 4 . So, compared to F 4 , F 2 may be considered as an IZA.

To have an opinion on the value given by the discrepancies defined in Eqs. [START_REF] Samagassi | Reconstruction of multiple impact forces by wavelet relevance vector machine approach[END_REF][START_REF] Robert | Multiple force identification for complex structures[END_REF], the discrepancy of both forces plotted in figure 9 is equal to E F 2 = 9, 09×10 -4 % and E F 4 = 3, 36% respectively. All the discrepancies are plotted in figure 10 in order to check the influence of the measurements (position and nature) on the results. Figure 10 shows the different discrepancies between the actual actions and the identified actions for different cases listed in Tables 1 and2. All these results are excellent.

For the NZA (Fig. 10 (b)), action pairs 1, 2, and 3 associated with measurement 2, have the highest discrepancies (about 4 %), which is very satisfactory. Figure 10 (a) shows that the results are also very satisfactory for IZA: the highest discrepancy is about 0.02 % and concerns action pair 4 associated with measurement 4. These results (Figure 10) are due to a good estimate of regularization parameters (σ 2 w , σ 2 η ) that are supposed to give a good reconstruction of the actions if they are well estimated. Also, note that a rapid con-vergence of the identification algorithm is related to the initialization of this algorithm and a good estimate of regularization parameters depends on the number of iterations. The parameter of interest (the actions) and the variance σ 2 w are the only parameters to be initialized. The variance σ 2 w is initialized via β w and k w from its prior distribution. By choosing values such as β w be as low as possible (k w = 2, β w = 10 -10 for example) and a large number of iterations (about 20,000), the regularization parameters are well estimated and therefore the actions are well identified.

The quality of the results may be due to three main factors. First how to model the data and the problem, second, the method used to solve it, and third, obtaining a good estimate of regularization parameters.

The difficulty in the Bayesian approach lies in its prior distribution that naturally imposes an intrinsic regularization: a prior well controlled favors a satisfactory regularized solution. A good modelling of this prior is therefore of capital importance. As it is mentioned in [START_REF] Zhang | Bayesian force reconstruction with an uncertain model[END_REF], this prior information can not be modeled exactly: as already mentioned, identifying actions in the CS context enables a better model of the prior information. The quantity of interest is reconstructed through its sparse representation which follows a normal Gaussian distribution. Algorithm 2 uses the Gibbs sampler method, which is very well adapted to identify forces because the laws of probability involved in this algorithm are usual, so samples can be easily drawn.

In the deterministic approach, the reconstruction of actions is much more difficult when at least two actions must be identified [START_REF] Tran | Reconstruction de sollicitations dynamiques par méthodes inverses[END_REF]. This difficulty originates in the choice of the regularization parameter which is determined by classical methods such as the L-curve [START_REF] Tran | Reconstruction de sollicitations dynamiques par méthodes inverses[END_REF]. Indeed, when there are several forces to identify, the number of the L-curve corners increases [START_REF] Tran | Reconstruction de sollicitations dynamiques par méthodes inverses[END_REF]: estimating the regularization parameter is then difficult.

Thus, it is not easy to do a good choice for the regularization parameter. In the Bayesian formulation, regularization parameters (σ 2 w and σ 2 η ) are variable, so they are adaptive. They are determined by an iterative process through Gibbs sampling. As a consequence, the regularization parameters are improved gradually through the MCMC process. Thus, the higher the number of iterations, the better the estimate of the regularization parameters and therefore the quality of identification.

Figure 10 also shows that it is difficult to predict the benefit of a measurement case to another. Overall, the identification is always satisfactory.

Identification of two forces

Two forces were applied along dof 2 (F 2 ) and dof 4 (F 4 ). They are identical and proportional to figure 1.

The identification is achieved with the measurement cases listed in Table 1. The discrepancies are shown in The same conclusion is reached for σ 2 η and σ 2 w (1): both histograms converged to the gamma distribution.

Indeed Fig. 12(d) shows a histogram which looks like to an exponential distribution: it is known that the exponential distribution is a gamma distribution with a shape parameter equal to 1.

To test the robustness of the convergence of the Markov Chain, we doubled the number of points (samples).

The results obtained with 512 points per force are illustrated on the figure.Thus, all of these histograms give us a good sign that the sampling algorithm has converged.

Identification of pressure

A uniform pressure, p(t), is applied on the third element: the shape of the pressure is still given by Fig.

1

. The pressure is projected on the finite element model dofs, so it is transformed in four actions: forces

F 4 = F 6 = p L e /2 and moments M 5 = -M 7 = p L 2 e /12
, where L e is the length of a finite element. The measurements of ν 2 , ν 4 , ν 5 , and ν 6 were used to identify the 4 actions. Indeed, experimentally it is much more easy to measure translations than rotations.

Fig. 14 shows that the identified forces are in a very good agreement with the initial force (2.36% for F 4 and 3.26% for F 6 ). However the moments are not satisfactorily identified since the discrepancies are about 78.29% for M 5 and 94.56% for M 7 .

This moment identifcation problem is certainly due to the difference in nature between force and moment that makes the simultaneous identification difficult.

However, the work provided by the pressure is mainly reflected in the work provided by the forces: the work of the moments is negligible as shown in figure 15. Consequently, this wrong identification of the moments would not have a significant incidence on a design of the beam: it is particularly true when the number of finite element increases. So, it is not necessary to identify the moments. In the following all the moments will be supposed to be IZA and won't be identified: this procedure will be referred as "under-identification".

The results are satisfactory as a discrepancy of 3.5 % (resp. 3.42 %) is obtained for F 4 (resp. F 6 ). It must be emphasized that if the responses are evaluated with the identified actions, they are very close to the actual responses.

Identification of the force location

In this subsection, the purpose is to achieve the location through the identification of multiple forces. So, impact identification involves both force localization and time history reconstruction. The structure is meshed and at each point of the mesh the force is identified: if, at a node of the mesh the force is almost equal to zero (IZA) then it is considered that no force is applied at this point. Thus, the location involves to identify both NZA and IZA: the forces are located at the points of the mesh where the identified force is not equal to zero. Obviously a coarse mesh leads to a coarse estimation of the force location and then there is a trade-off between accuracy and numerical cost.

All the dofs can be considered but it turns out that the results were much better when only translations were used: this simplification is more and more efficient when the number of nodes increases.

The beam is now divided in 7 elements (see Fig. 16). A pressure is applied on the second and the fifth ones:

that is forces are applied along translations 2, 4, 8 and 10, and moments are applied along rotations 3, 5, 9 and 11. As explained previously, the forces are identified but the moments are assumed to be equal to zero.

Also, translations ν 2 , ν 4 , ν 6 , ν 8 , ν 10 , and ν 12 are measured to identify the forces.

The shape of pressure p 1 (t) on the second element is represented in Fig. 1, whereas pressure p 2 (t) applied on the fifth element is equal to zero for a period of time and then equal to a sine signal (see 17(d)).

Note that the accuracy may depend on the size of the mesh, mainly when a distributed force has to be identified as the distributed force is replaced with point forces.

The results are illustrated in Fig. 17. The identification of the NZA is satisfactory. Indeed, the discrepancy between the real forces and the identified forces is 3.7 % for F 2 , 3.7 % for F 4 , 27.4 % for F 8 , and 28.2%

for F 10 . The large discrepancies observed with forces F 8 and F 10 come from the very end of the signal. A boundary effect could explain this feature. Indeed, the procedure identifies one large vector in which all the vectors are gathered. So the procedure must identify the large discontinuity between F 8 (n t ), and M 9 (1).

Similarly there is a large discontinuity between F 10 (n t ) and M 11 [START_REF] Klinkov | An updated comparison of the force reconstruction methods[END_REF]. It turns out that such discontinuities are difficult to identify. However, it may be considered that F 8 and F 10 are well identified except at the very end of the signal.

Identification of forces which are supposed to be identically zero is very satisfactory: F 6 and F 12 can be considered identically zero actions and then, dof 6 and 12 are not in the action location.

Thus, the results show that the Bayesian approach seems to be adapted to locate the action.

The main drawbacks of the procedure is the number of sensors required to measure displacements: 6 sensors were used for identifying 6 forces. However, the CS is claimed to be efficient even when the number of sensors is less than the number of actions. So the forces were identified with only 4 measurements: ν 2 , ν 4 , ν 8 , and ν 10 .

The results plotted in figure 18 are satisfactory and the same conclusions as before can be reached. Table 4: List of the dofs related to the action pairs to be identified

conclusions

A Bayesian method for multiple-force reconstruction has been proposed, discussed and validated on numerical examples. It was shown that identifying the signals in a basis where they are sparse is much robust with respect of the noise level.

In the numerical example, the identification of forces and moments applied on a beam was discussed. To improve the force identification problem, an under-identification was performed: the moments are supposed to be equal to zero. Finally, it was demonstrated it is possible to find the force location by identifying accurately an identically zero action from the identification of all the forces. The results obtained were conclusive. The possibility to identify actions with less sensors than forces was also discussed.

It is noted that the Bayesian approach has several advantages. First, the priori information in the form of a probability density function imposes an intrinsic regularization. Second, the Bayesian approach provides a rigorous probabilistic framework that takes into account all possible sources of errors (noise measurement and modeling error). Finally, it proposes the solution of the inverse problem in the posteriori probability density form from which drawings can be made. One of the key points that makes this approach feasible is its implementation using MCMC methods.

Multi-forces identification by Tikhonov regularization

In order to better control the degrees of freedom influencing the quality of reconstruction, Table 1 and2 have been slightly modified. Thus, we consider the tables 3 and 4 in the rest of the manuscript. Unless specified, the beam in 8 is considered in all that follows. Several measurements were considered. They are listed in Table 3. Similarly several action pairs were identified: they are listed in table 4. Identified forces F 2 and F 4 (force pair 1) are plotted in figure 19: the dof responses related to case 1 were used. The criterion to assess the discrepancy between the identified action (A id ) and the true action (A tr )

for ANN ( F 4 for example ) is rewritten :

D 24 (F 4 ) = ∥F id 4 -F tr 4 ∥ ∥F tr 4 ∥ × 100 (29)
In order to better appreciate discrepancy of IZA, a new criterion is proposed for the IZA:

D 0 24 (F 2 ) = 1 n t ∥F id 2 ∥ max(∥F tr 4 ∥) (30)
In this latter expression max(∥F tr 4 ∥) is a kind of scaling factor. However, if one action is a force and the other one is a moment, then D 0 ij has a unity. That is why it is not very obvious to compare the discrepancies between different IZAs. To make an opinion on these discrepancies, the discrepancy of the forces plotted in -Regularization parameter :

The worst configuration turned out to be for identifying an IZA along the translation at mid-span and a NZA along the rotation at mid-span for measurement case 8. These actions are depicted in figure 21. It is clear that a large spurious oscillation spoils the identified actions: the identification process does not provide a good result. This is due to a bad estimation of the regularisation parameter by the GCV method. Indeed, by using the regularization parameter given by the L-curve criterion, the results are excellent: the discrepancy of the forces is equal to D 0 25 (F 4 ) = 0.4 and D 25 (M 5 ) = 2.7 % respectively. This was observed for each bad result: usually, the GCV criterion gave the best estimation of the regularization parameter except in some cases; in that case, the one given by the L-curve provides a quite good result, that is with a discrepancy less than 15 % for the NZA.

The choice of the regularization parameter is much difficult when two actions have to be identified rather than one. Indeed, if the L-curve criterion, the GCV criterion or the quasi-optimality criterion seem to be suitable for a one load identification it is not always the case for the multiple action identification. So, a maximum of curvature is seek for the L-curve criterion. However, as shown in figure 22(a), several corners arise and the best choice is not always the maximum curvature corner.

It is similar for the GCV and the quasi-optimality criteria where a minimum is sought: several local minima are found (see figures 22(b) and 22(c)) and the minimum-minimorum is not always a good choice. That is an additional difficulty compared to the one action identification: in general, the number of corners or minima raises with the number of actions to be identified. So the advice is to combine the three criteria to make an opinion on the "best" choice, as already mentioned in [START_REF] Jacquelin | Force reconstruction analysis and regularization of a deconvolution problem[END_REF]. It is also of interest to mention the poor results given by the quasi-optimality criterion: in almost all the cases, it provides the worst choice, and it leads quite often to very poor identified actions. The results plotted in figure 20 did not show any influence on the location and the nature (translation or rotation) of the measured dofs, as far as two dofs are used.

Case 9 requires to measure the three translations (dofs 2, 4 and 6) but the results are not much better than case 5 where dofs 2 and 6 are measured. It is probably a consequence of the symmetry of the loading. However, a case involving a force acting along dof 2 only (not presented here) showed that a third measured dof does not improve significantly the identification.

All the simulations showed that the best case is when all the dof are measured: this is in agreement with the study on the sensor location and number [START_REF] Wang | Sensor placement methods for an improved force identification in state space[END_REF][START_REF] Chen | Impact source identification in finite isotropic plates using a timereversal method: theoretical study[END_REF], and this is expected as more information is used. This may be explained by checking the singular values of the transfer function, as explained in [START_REF] Jacquelin | Analysis of a force reconstruction problem[END_REF]. Indeed, when 2 or even with 3 dofs are used, the inverse problem may be not only ill-posed but also rank deficient. This is the case when case 2 of measurement points is used. Figure 23(a) shows a large gap between the first generalized singular value and the next one: so the first singular value is less than the computer precision. On the contrary, when all the dofs are used (case 10) the problem is not rank-deficient: the lower singular value is approximately 10 -3 and the singular value evolution does not show any large gap (see figure 23(b)). Case 10 is not the only case for which the rank is full: this is the case as soon as a measurement is done along each action (e.g. measurement case 1 for action pair 1). However, it was remarked that case 10 corresponds to the best conditioned case: in that specific case the regularization is even not needed as the least-square solution corresponds to a very good solution. Two forces were applied along dof 2 (F 2 ) and dof 4 (F 4 ). They are identical but a delay is introduced in force F 2 (see figure 24).

The errors were evaluated for the different measurement cases: the results are plotted in figure 25 for measurement case 1. No specific dof seems to be better than others. It was not expected that the error may be different for F 2 and F 4 , for the same case: cases 4 and 8 are even better than case 10 (all the dofs are used) to identify F 2 , but the error on F 4 is 4 time larger. The only case with low errors for both identified forces is case 10: once more, the redundancy helps. The pressure is projected on the finite element model dofs, so it is transformed in four actions: forces

F 4 = F 6 = p L e /2 and moments M 5 = -M 7 = p L 2 e /12
, where L e is the length of a finite element. All the dof are used. Note that the regularization requires at least as measurements as identified actions, that is at least four measurements. Case 10 (all the dofs are measured) was used to identify the 4 actions.

The results shown in figure 26 show that the identified forces are quite acceptable (discrepancy approximately 30 %) whereas the moments are not satisfactory (discrepancy approximately 90 %). However it is of interest to be aware that much of the work provided by the pressure is related to the work of the forces: figure 27 shows that the work of the moments is quite low compared to the work of the forces. Also, in this configuration the GCV was not suitable to provide the regularization parameter (see figure 28). Indeed the GCV is much too flat around the minimum and it turns out that the best solution was not for a minimum of the GCV function but for a corner of the L-curve: the latter value was in the flat part of the GCV curve. In this section, "under-identification" means that a NZA occurred but this action is supposed to be an IZA and then is not identified. The idea is to deal with the pressure identification issue and then to identify only the forces as the moments are poorly identified and the related work is low. Further, experimentally it is much more easy to measure translations than rotations. So the forces that come from the projected pressure on the finite element dofs are identified from the translations measurements.

The results are satisfactory as a discrepancy of 14 % (resp. 24 %) is obtained for F 4 (resp. F 6 ) that is even better than the previous ones whereas less measurements are used. It must be emphasized that if the responses are evaluated with the identified actions, they are very close to the actual responses (ie the ones given by the actual actions): this is observed for the response of all the dofs, even the rotations whereas the moments were neglected. Then from a design point of view the identified actions are really satisfactory:

figure 29 shows that the work provided by the identified forces is quite close to the work of the actual actions. As it is possible to identify several NZA and IZA, it may be possible to locate the forces actually applied on a structure, without any additional effort. To address this issue, the beam of Fig. 16 is considered. A pressure is applied on the second and the fifth ones: that is forces are applied along translations 2, 4, 8 and 10 and, moments are applied along rotations 3, 5, 9 and 11. As explained in the previous subsection, the forces are identified but the moments are assumed to be equal to zero.

The pressure on the second element is similar to the one proposed in the last two subsections, whereas the pressure applied on the fifth element is equal to zero for a period of time and then equal to a sine signal (see figures 30(a) and 30(e) for the shapes of the pressure).

The results are given in figure 30. As it is expected from the previous subsections, the NZA identification is satisfactory: the discrepancies are approximately 15 %. Similarly, the responses evaluated from the identified forces are very close to the true responses whereas the moments are neglected. However, the spurious oscillations that appear in the IZA have an amplitude much higher than observed previously: then from these results it is not possible de firmly decide whether F 6 and F 12 are an IZA. As a consequence, it is not obvious from these results that the third and the fifth elements are the only loaded elements.

A calculation was performed to identify all the actions (forces and moments): the results were almost unchanged regarding the forces. The moments were all almost equal to zero.Then it does not seem possible to locate the actions with only a force identification: further work is needed to identify better an IZA, even when several actions have to be identified.

Conclusions

Multi action identification was addressed: the inverse problem was solved with a Tikhonov regularization associated with a GCV criterion to find out the regularization parameter. First, this identification process was assessed by testing the ability to identify an identically zero action. The results were satisfactory. It appears that the quality of the identification does not seem to be correlated to a specific dof (location, translation or rotation), but it is improved when all the possible measurements are used.

Then several NZA were applied and the results were still satisfactory except for identifying moments.

However, it was shown that it is possible to identify a pressure on an element as soon as it is considered as two identical forces located on the nodes of the element and it is assumed that the moments related to the pressure projected on the finite element dofs are zero. However, when several actions that are very different are applied on a structure, it is almost impossible to identify know whether an action is identically zero or non zero.

The key point of a multiple load identification is the determination of the regularization parameter. It appeared that several minima occurred in the GCV function or in the quasi-optimality criterion, and the L-curve has several corners. All these criteria must be used, compared, and analysed in order to find out the optimal parameter. However, when actions of different nature must be identified (force and moment for example), several regularization parameters should be introduced: at least one for the force and another for the moments.

Further efforts may be done to improve the identification of an identically zero force in order to be able to locate the forces without any additional work.

Final conclusions

The main objective of this manuscript has been to study the possibility to locate load location through an identification of multiple forces that act on a linear elastic structure. The Bayesian approach together with the so-called compressed sensing technique as well as Tikhonov regularization associated with the GCV criterion were the two methods used to achieve this purpose.

The possibility to locate load location through an identification of multiple forces was based on an important note. It is important to note that if the structure is studied thanks to a finite element modelling, this includes the distributed force identification, such as a loading is transformed in some point loads (forces, moments) applied at the mesh nodes. Then, if it is possible to identify all the loads acting along all the degrees-offreedom (dof) of the structure, the location of the loads is also identified: a force (almost) equals to zero must be satisfactorily identified where no force is actually applied. This remark was the central theme of this paper and was used to validate (or not) the results. The quality of the result is therefore directly related to the ability of the methods used to identify a constantly effort equal to zero and to simultaneously identify multiple actions.

The results obtained by the Bayesian approach were much more eloquent than those obtained by the regularization of Tikhonv.

where:

π(s | w, w 0 , σ -2 w (i) , σ -2 η ) = π(s|w, w 0 , σ

Figure 1 : 3 Figure 2 :

 132 Figure 1: Shape of the actions

3 and the measurement is carried out at x j = L/ 2 .

 2 The reconstruction will be done on the interval [0 0.1]s with 256 points and the chosen time step value was ∆t = 3.90 10 -4 s.The five first frequencies of the pinned-pinned beam are: f 1 = 11, 68 Hz, f 2 = 46, 72 Hz,f 3 = 105, 1 Hz, f 4 = 186, 9 Hz and f 5 = 292 Hz. The time step was fixed so as to satisfy Shannon's sampling condition which is here equivalent to ∆t ≺ ∆t max = 1.7123 × 10 -3

  (a) and Fig.7(b)) and 2 % (Fig.7(c) and Fig.7(d)) of the response maximum amplitude. The identified forces are in a very good agreement with the initial force.

Figure 3 Figure 4 :Figure 5 :Figure 6 :Figure 7 :

 34567 Figure 3: A pinned-pinned beam

Figure 8 :Figure 9 :Figure 10 :

 8910 Figure 8: nodal action and dof of a 4-finite element discretized beam

Figure 11 . 2 w and σ 2 η

 1122 Figure 11. Measurement 2 associated with force F 4 has a discrepancy of the order of 4 % and discrepancy is about 3 % for all the other cases. Also, note that measurement 10 leads to both identified forces closest to the initial forces.The results given in Fig.12(a) are obtained with measurement case 10. Fig.12illustrates the process of reconstruction of both NZA. The histograms of the identified forces and the regularization parameters (σ 2 w

Figure 11 :Figure 13 :

 1113 Figure 11: NZA discrepancy for both forces: F2; F4

Figure 14 :Figure 15 :

 1415 Figure 14: Identified actions: (a) F4 (b) F6 (c) M5 (d) M7

Figure 16 :Figure 17 :Figure 18 :

 161718 Figure 16: Dofs of a 7-finite element discretized beam

5. 4 . 1 .Figure 19 :

 4119 Figure 19: Identified force along (a) dof 2 (b) dof 4

figure 19 is equal to D 0 24 (

 24 figure 19 is equal to D 0 24 (F 2 ) = 0.2 and D 24 (F 4 ) = 8.95 % respectively. All the discrepancies are plotted in figure 20 in order to check the influence of the measurements (position and nature) on the results.

Figure 20 :

 20 Figure 20: (a) IZA discrepancy (b) NZA discrepancy vs. the measurement cases for the different action pairs; pair 1: 'o'; pair 2: '+'; pair 3: '*'; pair 4: 'x'

Figure 21 :

 21 Figure 21: Identified action along (a) dof 5 (M5) (b) dof 4 (F4) for measurement case 8.

Figure 22 :

 22 Figure 22: Regularization parameter criterion: (a) L-curve, (b) GCV, (c) quasi-optimality

Figure 23 :

 23 Figure 23: Generalized singular values: (a) measurement case 1, (b) measurement case 10

Figure 24 : 2 Figure 25 :

 24225 Figure 24: Actual forces (dotted line) and identified forces (solid line) applied along (a) dof 4 and (b) dof 2

Figure 26 :Figure 27 :Figure 28 :

 262728 Figure 26: Identified actions: (a) F4 (b) F6 (c) M5 (d) M7

Figure 29 :

 29 Figure 29: True work (dotted line) and identified work (solid line) (a) of F4 and (b) F6

Figure 30 :

 30 Figure 30: True force (dotted line) and identified force (solid line) (a) F2, (b) F4, (c) F6, (d) F8, (e) F10, (f) F12

Table 1 :

 1 Measurement case list

	action pair number	1	2	3	4	5
	IZA dof	2	1	5	4	4
	NZA dof	4	4	4	5	2

Table 2 :

 2 List of the dofs related to the action pairs to be identified identified as well.Several measurements were considered. They are listed in Table1. Similarly several action pairs were identified: they are listed in Table2. Identified forces F 2 and F 4 (force pair 1) are plotted in figure

Table 3 :

 3 Measurement case list

	action pair number	1	2	3	4
	IZA dof	2	1	5	4
	NZA dof	4	4	4	5

A conjugate prior pdf is in the same pdf family as the posterior

AppendixA. Comparison between the Bayesian approach and the classical Tikhonov regularization [START_REF] Kaipio | Statistical and computational inverse problems[END_REF], [START_REF] Ayasso | Une approche bayésienne de l'inversion. Application à l'imagerie de diffraction dans les domaines micro-onde et optique[END_REF] Consider the following linear model:

with A ∈ R m×n a known matrix; X ∈ R n , Y ∈ R m and E ∈ R m are random vectors/matrices. E is an additive noise. Suppose further that X and E are mutually independent Gaussian vectors with probability densities:

and

With this information, Bayes' formula provides the posterior distribution of X knowing Y:

If Γ noise = σ 2 N I and e 0 = 0 , then

Therefore, if one chooses a MAP estimator, the optimization problem is equivalent to finding the minimum of the following criteria:

This corresponds to Tikhonov regularization where φ(x)/2 plays the role of the stabilizer functional. The minimum of this criterion is given by the relation:

AppendixB. Detailed description of the conditional pdfs in Algorithm 1

Conditional p.d.f. of the force

Hence:

Where:

with:

AppendixC. Identification of multiple forces When only one force, F, has to be identified, the problem to be solved is:

whereas in the case of multiple forces, it is: