Mathieu Aucejo
email: mathieu.aucejo@lecnam.net

Olivier De Smet
email: olivier.desmet@lecnam.net

BAYESIAN FORMULATIONS FOR FORCE RECONSTRUCTION PROBLEMS

Keywords: Force reconstruction, Bayesian inference, Generalized Gaussian priors, MCMC

To identify mechanical sources acting on a structure, Tikhonov-like regularizations are generally used. These approaches, however, only provide point estimates, meaning that the uncertainty about the regularized solution is not quantified. In practice, such information is essential to guarantee the quality of reconstructed sources. In this contribution, three possible Bayesian formulations of the source identification problem are presented and their limitations discussed. To assess the posterior uncertainty on the parameters appearing in each formulation given a simulated vibration field and a mechanical model, a Gibbs sampler is implemented. The proposed numerical validations highlight the practical interest of these formulations in terms of parameters estimations and posterior uncertainty quantification.

INTRODUCTION

In structural dynamics, most of the research papers are generally focused on the modeling of the dynamic response of structures subject to perfectly determined excitation sources. However, the latter are only roughly or partially known in practice. As a result, an error in the definition of the excitation vector is propagated to the dynamic response through the model and can have a significant impact on the subsequent mechanical analysis. However, direct measurement of excitation sources can be practically unfeasible. A possible alternative is to perform indirect measurements using a model of the dynamic behavior of the studied structure and accessible quantities such as displacement or acceleration fields. Unfortunately, the reconstruction of mechanical sources from vibration measurements is an ill-posed inverse problem. A classical approach to bypass this difficulty consists in constraining the space of solutions by using prior information on the noise and the sources to reconstruct. A convenient and efficient way to deal with such prior information is the Bayesian framework, because it allows combining both probabilistic and mechanical data. The most widespread and popular approaches deriving from Bayesian statistics are certainly Tikhonov-like regularizations [START_REF] Tikhonov | Regularization of incorrectly posed problems[END_REF][START_REF] Thite | The quantification of structure-borne transmission paths by inverse methods. Part 2 : Use of regularization techniques[END_REF][START_REF] Aucejo | Bayesian source identification using local priors[END_REF]. Although widely used and deeply studied, these methods generally provide only point estimates. In other words, there is no information on the uncertainty about the regularized solution given the measured data and the mechanical model. However, such information is of primary interest for industrial applications, in which it is essential to guarantee the quality of obtained results. In the present contribution, three possible Bayesian formulations of the source identification problem, based on the use of Generalized Gaussian distributions, Gamma and truncated Gamma distributions, are proposed. The limitations of each formulation are discussed from the identification of two close point forces acting on a free-free beam. To assess the posterior uncertainty on the parameters appearing in each formulation given a simulated vibration field and a mechanical model, a Gibbs sampler [START_REF] Geman | Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images[END_REF], including Hamiltonian Monte Carlo [START_REF] Duane | Hybrid Monte Carlo[END_REF] updates, is implemented to perform the inference. The proposed MCMC procedure is detailed in this contribution. The proposed numerical validations highlight the practical interest of the proposed formulations in terms of parameters estimations and posterior uncertainty quantification.

BAYESIAN FORMULATIONS OF THE RECONSTRUCTION PROBLEM

This section aims at introducing the three formulations of the Bayesian force reconstruction problem as well as the related parameters. To render this section more didactic, each formulation is applied on an academic test case in order to better highlight its advantages and limitations. The MCMC algorithm used to perform the inferences will be detailed in section 3.

Description of the benchmark test case

Before presenting the three Bayesian formulations of the reconstruction problem, we propose to introduce the test case that will serve as a benchmark to analyze and compare each of them. The studied structure is a free-free steel beam with dimensions 1×0.03×0.01 m 3 excited by two point forces of unit amplitude at 350 Hz. The coordinate of the point forces, measured from the left end of the beam, are x 1 = 0.6 m and x 2 = 0.7 m [see Fig. 1].

To perform the reconstruction, a model has to be derived to relate the measured vibration field and the excitation field. For this purpose, let us consider the practical situation where the vibration field X, measured over the surface of a structure, is caused by an unknown excitation field F. If the structure is linear and time invariant, its dynamic behavior is completely deter- mined by the transfer functions matrix H, relating the vibration field X to the excitation field F, so that:

X = HF + N, (1)
where N is the measurement noise vector.

A finite element model of the beam made up with 20 plane beam elements has been used to compute the transfer functions matrix H by assuming that only bending motions are measurable. In other words, the transfer functions matrix is dynamically condensed over the measurable dofs, corresponding here to the bending motions. To simulate the measured vibration field X, Eq. ( 1) is applied using a noise vector N corresponding to a Gaussian white noise, which is calculated assuming a signal-to-noise ratio equal to 34 dB. Finally, it is worth mentioning that the model defined in Eq. ( 1) is the basis of the proposed Bayesian formulations, meaning that the transfer functions matrix H is also used in the inverse problem. In other words, modeling errors are not considered here.

Standard Bayesian formulation

Formally, the Bayesian paradigm considers all the parameters of the problem as random variables. Consequently, the uncertainty on each parameter is modeled by a probability distribution, describing the state of knowledge or the prior on this parameter. From a mathematical standpoint, the Bayesian reconstruction problem relies on the Bayes' rule:

p(F|X) ∝ p(X|F) p(F), (2)
where:

• p(F|X) is the posterior probability distribution, representing the probability of observing F given a vibration field X. It defines what it is known about the excitation field F after making vibration measurements;

• p(X|F) is the likelihood function, representing the probability of measuring X given an excitation field F. It reflects the uncertainty related to the measurement of the vibration field X;

• p(F) is the prior probability distribution, representing our knowledge on the unknown excitation field F before measuring the vibration field X.

Generally, the quality of the force reconstruction strongly depends on the choice of the likelihood function and the prior probability distribution. That is why, the choices made in this contribution have to be carefully explained.

Choice of the likelihood function

The likelihood function reflects the uncertainty related to vibration measurements. By definition, this uncertainty is mainly related to the measurement noise N. Consequently, the likelihood function p(X|F) can be written under the following form:

p(X|F) = p(X -HF|N), (3)
representing the probability of obtaining X -HF = 0 given the measurement noise N.

If the noise is supposed to be due to multiple independent causes, then the likelihood function can be represented by a complex multivariate normal distribution with zero mean and precision τ n :

p(X|F, τ n ) = τ n π N exp -τ n X -HF 2 2 , (4)
where N is the number of measurement points.

Choice of the prior probability distribution

The prior probability distribution reflects the uncertainty related to the unknown excitation field F. Actually, it can be seen as a measure of the a priori knowledge of the experimenter on the sources to identify.

For practical reasons, the excitation field F is supposed to be a real random vector, whose components are independent and identically distributed random variables following a Generalized Gaussian distribution. As a result, the prior probability distribution is written:

p(F|τ s , q) = q 2 Γ(1/q) M τ M q s exp -τ s F q q , (5)
where:

• q is the shape parameter of the distribution. Its value is defined in the interval ]0, +∞[; • • q is the q -norm or quasi-norm, if q ≥ 1 and q < 1 respectively;

• τ s is the scale parameter of the distribution, which can be viewed as a generalized measure of the precision of the distribution;

• M is the number of reconstruction points;

• Γ(x) = +∞ 0 t x-1 e -t dt is the gamma function.
It should be noted that the choice of a multivariate generalized Gaussian distribution allows a great flexibility for describing prior knowledge of the sources to identify. Formally, sparse excitation fields are promoted for q ≤ 1, while distributed excitation fields are favored for q = 2 [START_REF] Boyd | Convex optimization[END_REF]. In practical situations, the possible values of the shape parameter are bounded in the interval ]0, 2].

Summary and application

From the explanations given above, the standard Bayesian formulation of the reconstruction problem finally writes:

p(F|X, τ n , τ s , q) ∝ p(X|F, τ n ) p(F|τ s , q). ( 6
)
This formulation is said standard, because it leads to Tikhonov-like regularizations, which correspond to the MAP estimate of Eq. ( 6). To explore the posterior probability distribution, 10 5 samples are drawn from Eq. ( 6). Fig. 2 and Tables 1 and2 present the excitation fields and the estimated values of the forces F 1 and F 2 obtained using either with a proper choice of q (q = 0.5) [see Fig. 2a] or with a poor choice of q (q = 2) [see Fig. 2b]. It should be noted that optimal values of τ n and τ s are strongly related to the value of q and are computed accordingly as explained in section 3. Here, (τ n , τ s ) = (2.77 × 10 15 , 1.48) for q = 0.5 and (τ n , τ s ) = (2.87 × 10 15 , 12.44) for q = 2. Table 1: Standard Bayesian formulation -Summary of the inference result on the model parameters for q = 0.5

Several interesting conclusions can be drawn at the light of the results presented in Fig. 2 and Tables 1 and2. First of all, the choice of the supposedly known parameters have a significant impact on the quality of reconstruction. Indeed, a proper choice of the shape parameter leads to a very good agreement of the reference field with the median of the samples, associated with sharp posterior uncertainty region characterizes by the related 95% credible interval. On the contrary, a poor choice of the shape parameter leads to a disappointing result, since the shape of the median of the samples does not allow to discriminate the two point forces. Furthermore, the posterior uncertainty region is quite large, which is closely related to the value of the shape parameter q. All things being equal, the larger the value of the shape parameter is, the larger is the area of the posterior uncertainty region.

The standard formulation gives information on the credibility of the identified excitation field given the measured vibration field X, the precisions τ n and τ s and the shape parameter q. Actually, this is the main drawback of the standard formulation, because the quality of the inference is conditioned to the knowledge of the precisions and the shape parameter. If their values are poorly chosen, then the resulting inference won't be representative of the actual distribution. As a consequence, it is compulsory to determine near-optimal values of q, τ n and τ s if one wants to perform a relevant statistical inference.

Extended Bayesian formulation

To alleviate the limitations of the standard formulation, the approach generally adopted in the literature consists in considering the precisions τ n and τ s as random variables, while letting fixed the shape parameters q. If we further consider the precisions as independent variables, the following extended formulation is obtained:

p(F, τ n , τ s |X, q) ∝ p(X|F, τ n ) p(F|τ s , q) p(τ n ) p(τ s ), (7)
where p(τ n ) and p(τ s ) are the prior probability distributions of the precisions τ n and τ s respectively.

Choice of the prior probability distribution of the precisions

The choice of the priori probability distributions p(τ n ) and p(τ s ) is first limited to distribution having a strictly positive support, because the precisions τ n and τ s are real positive numbers. The common choice, made in the literature, is the Gamma distribution. The reason for this is rather clear, since the conjugate prior for the precision of a generalized Gaussian distribution is a Gamma distribution [START_REF] Bishop | Pattern recognition and machine learning[END_REF]. Practically, the Gamma distibution is defined by:

G(τ |α, β) = β α Γ(α) τ α-1 exp(-β τ ) with α > 0, β > 0, (8)
where α and β are respectively the scale parameter and the rate parameter of the distribution.

However, the use of a Gamma distribution is questionable, since it has been chosen for mathematical convenience and does not reflect any real prior information on the precisions, except their positiveness. That is why, the prior distribution on τ n and τ s should be as minimally informative as possible [START_REF] Hoff | A first course in Bayesian statistical methods[END_REF]. To this end, one sets α n = α s = 1 and β n = β s → 0.

Summary and application

Considering the previous choice, the extended Bayesian formulation of the reconstruction problem is given by:

p(F, τ n , τ s |X, q) ∝ p(X|F, τ n ) p(F|τ s , q) p(τ n |α n , β n ) p(τ s |α s , β s ), (9)
where (α n , β n ) are the hyperparameters related to the precision τ n , while (α s , β s ) are the hyperparameters related to the precision τ s .

This extended formulation has given rise, when q = 2, to the augmented Tikhonov regularization [START_REF] Bin | Augmented Tikhonov regularization[END_REF]. This method provides a point estimate corresponding to a critical point of the opposite of the logarithm of the posterior probability distribution. The main advantage of this approach is to determine the regularized solution and the precision simultaneously using an iterative process.

As previously done, 10 5 samples have been drawn from the MCMC algorithm presented in section 3 to explore the posterior probability distribution. Fig. 3 presents the excitation fields obtained using either a proper choice of q (q = 0.5) [see Fig. 3a] or with a poor choice of q (q = 2) [see Fig. 3b], while Tables 3 and4 summarize the inference results on each parameter of the model for q = 0.5 and q = 2 respectively. Figure 3: Extended Bayesian formulation -Real part of the reconstructed force vector (a) with a proper choice of q (q = 0.5) and (b) a poor choice of q (q = 2) -(-) Reference, (--) Median of the samples and ( a ) 95% credible interval

Obtained results clearly show that the quality of the inference strongly depends on a proper choice of the shape parameter q. Consequently, when setting q = 2, as classically done in the literature to perform the inference, one takes the risk to draw erroneous conclusions if the structure is actually excited by localized sources, even if the precisions are estimated from Bayesian inference.

Complete Bayesian formulation

The applications of standard and extended Bayesian formulations have pointed out the need for properly defining the value of the shape parameter q. However, choosing a priori relevant value is far from an easy task for non-experienced user. That is why, it is interesting to infer the shape parameter from a Bayesian analysis. Practically, this is done by considering this parameter as a random variable. In doing so, one obtains the complete Bayesian formulation:

p(F, τ n , τ s , q|X) ∝ p(X|F, τ n ) p(F|τ s , q) p(τ n |α n , β n ) p(τ s |α s , β s ) p(q), (10)
where p(q) is the prior probability distribution of the shape parameter q.

Choice of the prior probability distribution of the shape parameter

The only available information is that the value of the shape parameter is bounded and positive. In absence of more precise knowledge on this parameter, the probability distribution is not only chosen to reflect the available information but also for its mathematical tractability. A probability distribution that meets these requirements is the truncated Gamma distribution defined by:

G T (q|α q , β q , l b , u b ) = Γ(α q ) γ(α q , β q u b ) -γ(α q , β q l b ) G(q|α q , β q ) I [l b ,u b ] (q), (11)
where:

• G(q|α q , β q ) is the Gamma distribution defined in Eq. ( 8);

• I [l b ,u b ] (q) is the truncation function defined between the lower bound l b and the upper bound u b . More precisely, this function simply writes:

I [l b ,u b ] (q) = 1 if q ∈ [l b , u b ] 0 otherwise ; (12)
• γ(s, x) = x 0 t s-1 exp(-t) dt is the lower incomplete Gamma function. Here, the choice of the truncated Gamma distribution has been made for mathematical convenience, because other continuous truncated distribution could have theoretically been used. To avoid biasing the inference, the shape of the prior distribution needs to be weakly informative. Here, this means that the hyperparameters should be defined such that α q = 1 and β q → 0. On the other hand, even if the lower and upper bounds l b and u b can theoretically take any positive value, one knows that the value of q practically lies in the interval ]0, 2]. For this particular reason, we set l b = 0.05 and u b = 2.05.

Summary and application

From the above considerations, the complete Bayesian formulation is given by:

p(F, τ n , τ s , q|X) ∝ p(X|F, τ n ) p(F|τ s , q) p(τ n |α n , β n ) p(τ s |α s , β s ) p(q|α q , β q , l b , u b ). ( 13
)
To explore the posterior probability distribution, 10 5 samples have been drawn from the MCMC algorithm presented in section 3. Fig. 4 presents the excitation field obtained after performing the inference, while Table 5 summarizes the inference results on the model parameters. As expected, obtained results clearly show that the inference performed from the complete Bayesian formulation is able to provide parameters estimates as well as a quantification of the posterior uncertainty on those parameters. In particular, it can be stressed that the median of the reconstructed excitation field agrees very well with the reference one. Furthermore, it could be added that the size of the associated 95% credible interval is small, meaning that the reconstructed excitation field is highly probable given the mechanical model and the measured vibration field.

BAYESIAN INFERENCE -MCMC ALGORITHM

In the previous section, the main features of each Bayesian formulation have been studied without describing the algorithm behind the results presented. Actually, a Gibbs sampler have been implemented to explore the posterior probability distribution of each Bayesian formulations. To render the presentation of the proposed sampler more concise, the proposed Gibbs sampler is derived for the complete Bayesian formulation, since the samplers associated to the standard and the extended Bayesian formulations are, in fact, only particular cases of this more general sampler.

The implementation of a Gibbs sampler requires the knowledge of the full conditional probability distributions. For the complete Bayesian formulation, one has for:

• the shape parameter q: p(q|X, F, τ n , τ s ) ∝ τ M q s Γ(1/q) M q αq+M -1 exp -β q q -τ s F q q I [l b ,u b ] (q); (14)
• the precision τ s :

p(τ s |X, F, τ n , q) ∝ G τ s α s + M q , β s + F q q ; (15)
• the precision τ n :

p(τ n |X, F, τ s , q) ∝ G τ n α n + N, β n + X -HF 2 2 ; (16)
• the force vector F:

p(F|X, τ n , τ s , q) ∝ exp -τ n X -HF 2 2 -τ s F q q , (17)
which corresponds to the standard Bayesian formulation.

From the previous full conditional probability distributions, we can derive the following Gibbs sampler:

1. Set k = 0 and initialize q (0) , τ

(0) s , τ (0)
n and F (0) ; 2. Draw N s samples from full conditional distributions

for k = 1 : N s a. draw q (k) ∼ p q|X, F (k-1) , τ (k-1) n , τ (k-1) s b. draw τ (k) s ∼ p τ s |X, F (k-1) , τ (k-1) n , q (k) c. draw τ (k) n ∼ p τ n |X, F (k-1) , τ (k) s , q (k) d. draw F (k) ∼ p F|X, τ (k) n , τ (k) s , q (k)
end for

Monitor the convergence of the Markov chains

In the next of this section, each step of the Gibbs sampler is detailed.

Initialization of the sampler

In the present contribution, the sampler is initialized from a starting point having a reasonably high probability. Here, the initial force vector F (0) is obtained from the MAP estimate of the standard Bayesian formulation, which is defined by [see Ref. [START_REF] Aucejo | Bayesian source identification using local priors[END_REF] for more details]:

F (0) = argmin F X -HF 2 2 + λ (0) F q (0) q (0) , (18)
where

λ (0) = τ (0) s /τ (0)
n is the regularization parameter.

To obtain a relevant initial force vector, it is necessary to determine reasonable values of the shape parameters q (0) and the precisions τ (0) s and τ (0) n . Practically, the values of the shape parameter q (0) can be chosen without any calculation using the indications given in section 2.2.2. However, because q (0) can take any value in the range ]0,2], the solution of the previous optimization problem has generally no closed-form expression and hence the optimal value of the regularization parameter λ (0) can not be directly computed from automatic selection procedures such as the L-curve principle [START_REF] Hansen | Rank-Deficient and Discrete Ill-Posed Problems: Numerical Aspects of Linear Inversion[END_REF].

Practically, the optimization problem is solved using an Iteratively Reweighted Least Squares (IRLS) algorithm [START_REF] Rodriguez | An Iteratively Weighted Norm Algorithm for Total Variation Regularization[END_REF]. From a Bayesian standpoint, the IRLS procedure can be viewed as a transformation of the Generalized Gaussian prior into a multivariate Gaussian-like prior. More precisely, one has:

p(F|X, τ (0) s , q (0) ) = p(F|τ (0) s , q (0) ) ∝ p(F|W, τ (0) s ) ∝ exp -τ (0) s W 1/2 F 2 2 , ( 19
)
where W is a diagonal global weighting matrix depending explicitly on F and q (0) . This explains why Eq. (18) can only be solved using an iterative procedure.

At each iteration j of the IRLS algorithm, one has to solve:

F (0,j) = argmin F X -HF 2 2 + λ (0,j) W (j) 1/2 F 2 2 , (20)
until the convergence is reached.

After convergence of the iterative process, the optimal force vector F (0) , the global weighting matrix W, as well as the optimal value of the regularization parameter λ (0) , that has been updated at each iteration using automatic selection procedure such as the L-curve principle, are obtained. Consequently, to complete the initialization step it remains to determine the values of τ (0) n and τ (0) s . Here, we follow the approach proposed by Pereira et al. [START_REF] Pereira | Empirical Bayesian regularization of the inverse acoustic problem[END_REF] consisting in finding the optimal values of τ (0) n and τ (0) s given the measured vibration field X only using a Bayesian approach. With the proposed parametrization, one has:

τ (0) s = N X H λ (0) I + HW -1 H H -1 X and τ (0) n = τ (0) s λ (0) . ( 21
)

Drawing samples for the shape parameter

The Gibbs sampler is generally the first choice for conditionally conjugate models, where samples can be drawn directly from each conditional probability distribution. Unfortunately, this is not the case for the shape parameters q, for which the conditional distribution is not standard. To bypass this difficulty, we propose to update the value of the shape parameter using a Hamiltonian Monte Carlo (HMC) step, because it limit the random walk behavior of the sampler.

In HMC, the Hamiltonian function can be written as follows:

H(q, s) = U (q) + K(s), (22)
where U (q) is called the potential energy and is defined as:

U (q) = -log p q|X, F (k-1) , τ (k-1) n , τ (k-1) s , (23)
while K(s) is called the kinetic energy and is chosen for mathematical convenience. Usually, it is defined such that:

K(s) = -log [N (s|0, 1)] = 1 2 s T s + const. (24)
In this framework, q and s define the state space of the dynamical system and can be viewed as the position and the momentum respectively. The evolution of the state (q, s) over the time t is governed by the following system of equations:

     dq dt = ∂H(q, s) ∂s ds dt = - ∂H(q, s) ∂q . (25)
Practically, Eq. ( 25) is solved using the leapfrog method. Starting from an initial state (q 0 , s 0 ), it allows determining a proposal state (q , s ), while taking into account the constraints on the parameter q. Once a proposal state is obtained, one has to decide to accept it as a new state with probability:

P = min (1, exp [H(q 0 , s 0 ) -H(q , s )]) , (26)
or to reject it with probability 1 -P .

Drawing samples for the precisions

The conditional probability distribution for the precisions τ n and τ s are actually gamma distributions. Consequently, samples can be easily drawn from standard statistical packages.

Drawing samples for the force vector

At first sight, drawing samples from p F|X, τ

(k) s , τ (k)
s , q (k) seems a difficult task. However, the IRLS algorithm allows writing when the iterative process has converged:

p F|X, τ (k) n , τ (k) s , q (k) ∝ exp -τ (k) n X -HF 2 2 -τ (k) s W 1/2 F 2 2 ∝ N c (F|µ F , Σ F ), (27)
where N c (z|µ, Σ) is the (circurlarly-symmetric) complex multivariate Gaussian distribution with mean µ and covariance matrix Σ. Here, the mean µ F and the covariance matrix Σ F are expressed as:

µ F = τ (k) n Σ F H H X and Σ F = τ (k) n H H H + τ (k) s W -1 . (28)
Once µ F and the covariance matrix Σ F are obtained, samples of the force vector F (k) can be easily drawn. A widely used method for drawing samples from a complex multivariate Gaussian distribution consists in:

1. generating a complex random vector v ∼ N c (v|0, I) from two real normal random vectors

v 1 ∼ N (v 1 |0, I) and v 2 ∼ N (v 2 |0, I): v = v 1 + j v 2 √ 2 , (29)
where j = √ -1 is the imaginary unit;

2. determining any matrix L satisfying the relation LL H = Σ F . This can be typically done using either a Cholesky decomposition or a spectral decomposition;

3. computing the new sample F (k) = µ F + L v.

Convergence diagnostics

In the present contribution, all the inferences have been performed from one single chain for each monitored parameter. More precisely, we have monitored the convergence of the precisions (τ n , τ s ) and the shape parameter q. The procedure used for convergence monitoring is based on the combination of two complementary diagnostics:

1. the Raftery-Lewis diagnostic [START_REF] Raftery | One long run with diagnostics: Implementation strategies for Markov chain Monte Carlo[END_REF] that estimates the total run length and the burn-in period;

2. the Geweke diagnostic [START_REF] Geweke | Evaluating the Accuracy of Sampling-Based Approaches to Calculating Posterior Moments[END_REF] that tests the null hypothesis that the Markov chain is in the equilibrium distribution and produces z-statistics for each estimated parameter.

Practically, the Raftery-Lewis is first applied to the chains to monitor. Then, the Geweke diagnostic is applied to the resulting chains.

Samplers for the standard and the extended Bayesian formulations

The MCMC algorithm for the complete Bayesian formulation can be used to derive the MCMC algorithms corresponding to the standard and the extended Bayesian formulations. Indeed, the MCMC algorithm related to the extended Bayesian formulation is obtained by not considering the step 2a of the Gibbs sampler (i.e. by fixing the values of q), while the sampler associated to the standard Bayesian formulation is derived from the general Gibbs sampler by computing steps 1 and 2c only (i.e. by considering fixed values for τ n , τ s and q).

CONCLUSIONS

In structural dynamics, regularization approaches can be used to deal with source identification problems. Although widely used and deeply studied, these methods generally provide only point estimates. In other words, there is no information on the uncertainty about the regularized solution given the measured data and the mechanical model. However, such information is of primary interest for industrial applications, in which it is essential to guarantee the quality of obtained results. In this contribution, three Bayesian formulations of the force reconstruction problem have been proposed to tackle this issue and an original hybrid Gibbs sampler has been implemented to perform the inferences. One of the merits of proposed sampler is to be flexible enough to deal with each formulation and provide estimates and uncertainty quantification of all the parameters of the formulation considered. The proposed numerical example clearly highlight the benefits and the limitations of the proposed formulations. In particular, it has been shown that the standard and extended formulations should be used when precise information about the problem is available. On the contrary, the complete formulation should be used in case of too vague information about the problem. In this case, however, the solution is obtained at a higher computational cost.

Figure 1 :

 1 Figure 1: Reference force field

Figure 2 :

 2 Figure2: Standard Bayesian formulation -Real part of the reconstructed force vector (a) with a proper choice of q (q = 0.5) and (b) a poor choice of q (q = 2) -(-) Reference, (--) Median of the samples and ( a ) 95% credible interval

Figure 4 :

 4 Figure 4: Complete Bayesian formulation -Real part of the reconstructed force vector -(-) Reference, (--) Median of the samples and ( a ) 95% credible interval

Table 2 :

 2 Standard Bayesian formulation -Summary of the inference result on the model parameters for q = 2

Table 3 :

 3 Extended Bayesian formulation -Summary of the inference result on the model parameters for q = 0.5

	Parameter	Median	Mode	95% CI
	F 1	0.565	0.565	[0.253, 0.889]
	F 2	0.497	0.497	[0.181, 0.817]
	τ n	4.29 × 10 15 4.04 × 10 15 [2.17, 7.80] × 10 15
	τ s	10.37	9.81	[5.38, 17.88]

Table 4 :

 4 Extended Bayesian formulation -Summary of the inference result on the model parameters for q = 2

Table 5 :

 5 Complete Bayesian formulation -Summary of the inference result on the model parameters


Mathieu Aucejo

Olivier De Smet

On a space-frequency regularization for source reconstruction

To identify mechanical sources acting on a structure, Tikhonov-like regularizations are generally used. These approaches, referred to as additive regularizations, require the calculation of a regularization parameter from adapted selection procedures such as the Lcurve method. However, such selection procedures can be computationally intensive. In this contribution, a space-frequency multiplicative regularization is introduced. The proposed strategy has the merit of avoiding the need for the determination of a regularization parameter beforehand, while taking advantage of one's prior knowledge of the type of the sources as well as the nature of the excitation signal. By construction, the regularized solution is computed in an iterative manner, which allows adapting the importance of the regularization term all along the resolution process. The validity of the proposed approach is illustrated numerically on a simply supported beam.

Introduction

The knowledge of external loads is essential to properly design a mechanical system or control the vibration levels. The direct measurement of excitation sources is generally the first-line approach. However, such strategy is intrusive and potentially time-comsuming to implement and calibrate. Furthermore, direct measurements can be impossible in practical situations. To circumvent these undesirable features, a possible alternative is to perform indirect measurements using accessible quantities such as displacement or acceleration fields and a model of the dynamic behaviour of the studied structure. Unfortunately, the reconstruction of mechanical sources from vibration measurements is an ill-posed inverse problem. A classical approach to bypass this difficulty consists in including in the formulation of the inverse problem a regularization term reflecting one's prior knowledge on the spatial distribution of sources or the type of excitation signal to constrain the space of solutions. This idea is at core of Tikhonov-like regularization methods [1,2,3], that are widely used to tackle the inverse problem in the frequency domain. However, the reconstruction problem is generally solved frequency by frequency, meaning that only the spatial prior information on the sources to identify is exploited. From a spectral point of view, this is equivalent to suppose that the frequency spectrum of the identified sources is discontinuous. This lack of continuity can induce potential inaccuracies in the reconstructed excitation fields if the sources are broadband [4]. Consequently, it appears that the vast majority of the frequency methods proposed in the literature are unable to simultaneously exploit both the spatial and spectral features at the same time. When one seeks to localize excitation sources and determine their frequency content, it is crucial to take advantage of both information to aid the reconstruction process in converging to a consistent regularized solution. These requirements are actually satisfied by regularization terms derived from a mixed p,q -norm. In the context of force reconstruction, Rezayat et al. first derive a Tikhonov-like approach using a regularization term based on a mixed 2,1 -norm in order to identify broadband point forces [4]. From a technical standpoint, the problem is then solved by a generalization of the Fast Iterative Shrinkage-Thresholding Algorithm (FISTA) proposed by Beck and Teboulle [5].

In the present contribution, an original regularization strategy is proposed to solve both the localization and the spectral reconstruction problem within a unique framework. The proposed approach first relies on the definition of a regularization term constructed from a mixed 2,q -norm. Such a space-frequency regularization term has the merit to be flexible enough to reflect one's prior knowledge on the type of the excitation sources and the nature of the excitation signal. Then, this regularization term is included in the formulation of the inverse problem as a multiplicative constraint. This regularization strategy, originally developed by Van den Berg et al. [6], has several advantages compared to the more classical Tikhonov-like regularizations (a.k.a. additive regularizations). In particular, it avoids the use of expensive selection procedures such as the L-curve principle [START_REF] Hansen | Rank-Deficient and Discrette Ill-Posed Problems: Numerical Aspects of Linear Inversion[END_REF], since, by construction, the regularized solution has to be computed in an iterative manner. This specific feature allows adapting the importance of the regularization term all along the resolution process. Practically, the proposed inverse formulation is solved from an adapted Iteratively Reweighted Least-Squares (IRLS) algorithm. The ability of the proposed space-frequency approach in providing consistent reconstructions is illustrated on a simply supported beam excited by a broadband point force using synthetic data.

Reconstruction model

The implementation of the space-frequency regularization requires the definition of a reconstruction model describing the relation between the measured vibration field and the excitation field to identify. If the structure is linear and time-invariant, its dynamic behaviour at a particular frequency f j is completely defined by the transfer functions matrix H(f j ) relating the vibration field X(f j ) to the unknown excitation field F(f j ). Now, let us suppose that the vibration field is corrupted by a measurement noise N(f j ). In such a situation, the reconstruction model simply writes:

X(f j ) = H(f j )F(f j ) + N(f j ), ∀ j = 1, . . . N, (1)
where N is the number of studied frequencies.

From the previous equation, it is possible to construct a global model by collecting all the studied frequencies together. In doing so, the global reconstruction model is:

X = H F + N, (2)
where the global quantities are given by:

H =    H(f 1 ) 0 . . . 0 H(f N )    , X =    X(f 1 ) . . . X(f N )    , F =    F(f 1 ) . . . F(f N )    and N =    N(f 1 ) . . . N(f N )    . (3)

Proposed regularization strategy

In the present contribution, it is proposed to seek the unknown global excitation field F as the solution of the following minimization problem:

F m = argmin F F X -H F • R F , (4)
where:

• F X -H F is the data-fidelity term that encodes prior information on the noise N corrupting the data [START_REF] Nikolova | [END_REF];

• R F is the regularization term that controls the a priori on the unknown excitation field

F [9].
To obtain consistent reconstructions, a proper choice of the data-fidelity and regularization terms is crucial, since the more prior information on the noise and the sources is meaningful, the more the confidence in the reconstruction is high [3].

Definition of the data-fidelity term

As evoked previously, the data-fidelity term reflects prior information on the noise corrupting the data. A common assumption consists in considering that the vibration field is corrupted at each frequency by an additive Gaussian white noise. In this case, the data-fidelity term can be expressed as:

F X -H F = X -H F 2 2 .
(5)

Definition of the regularization term

The main objective of this contribution is to simultaneously exploit one's prior knowledge of the type of sources (localized or distributed) and the nature of the excitation signal (broadband or narrowband). To this end, it is relevant to define a regularization term based on the mixed p,q . Indeed, if one represents the global unknown force vector F as a matrix, where the rows correspond to the force spectrum at a particular location and the columns to the force field at a specific frequency, it comes:

F = [F(f 1 ) • • • F(f j ) • • • F(f N )] =        F 1 (f 1 ) • • • F 1 (f j ) • • • F 1 (f N ) . . . . . . . . . F i (f 1 ) • • • F i (f j ) • • • F i (f N ) . . . . . . . . . F M (f 1 ) • • • F M (f j ) • • • F M (f N )        , ( 6
)
where M is the number of spatial reconstruction points.

Furthermore, if one reminds that the mixed p,q -norm (or quasinorm) is defined by:

F p,q =    M i=1   N j=1 |F ij | p   q p    1 q , ∀ (p, q) ∈ ]0, +∞[ 2 , (7)
it becomes clear that such a regularization term introduce explicitly a coupling between the coefficients of the matrix. Consequently, a mixed norm allows promoting some structures observed in real signals [10]. To illustrate this particular property of mixed norms, let us consider the case (p, q) = (2, 1). In this situation, the global force F is supposed to be sparse along the lines (space) and full along the rows (frequency). In other words, one promotes the spatial sparsity of the excitation field (localized sources) and the continuity of its frequency spectrum (broadband source signal).

From the previous considerations, the proposed regularization term is such that:

R F = F q 2,q = M i=1 F[i] q 2 , (8)
where F[i] is the i th row of the global force vector written under the matrix form [see Eq. ( 6)].

Consequently, we restrict ourselves to localized (q ≤ 1) or distributed (q = 2) broadband (p = 2) excitation sources.

Practical form of the identification problem

From the foregoing, the generic form of the proposed space-frequency multiplicative regularization is given by:

F m = argmin F X -H F 2 2 • F q 2,q . (9)

Resolution of the identification problem

The solution of the minimization problem given in the previous section can only be found in an iterative manner, insofar as no explicit solution exists for Eq. (9). For this particular reason, the regularized solution is computed from an adapted Iteratively Reweighted Least-Squares (IRLS) algorithm [START_REF] Rodriguez | Proc. of the 2006 Asilomar Conference on Signals, Systems and Computers[END_REF].

General principle

The core idea of the IRLS algorithm is to replace the direct resolution of the minimization problem by an equivalent iterative process having an explicit unique solution at each iteration.

To this end, the 2,q -norm is replace by a weighted 2 -norm, that is:

∀ q, F q 2,q = W(F) 1/2 F 2 2 , (10)
where W(F) is the global weighting matrix defined as a function of F.

As part of an iterative scheme, one tries to find the solution F ), so as to recover the equality (10) when the convergence is reached. Here, the excitation field F (k+1) m at iteration k+1 is the solution of the following minimization problem:

F (k+1) m = argmin F X -H F 2 2 • W (k) 1/2 F 2 2 . (11)
In the previous relation, the global weighting matrix W

(k)
is such that:

W (k) = I N ⊗ W (k) , (12)
where I N is the identity matrix of dimension N and W (k) is the reduced weighting matrix.

By construction, the reduced weighting matrix writes:

W (k) = diag W (k) 1 , . . . , W (k) i , . . . , W (k) M , (13)
where:

W (k) i = max , F (k) m [i] 2 q-2 . ( 14
)
In the previous relation, is a damping parameter that allows avoiding infinite weights when

F (k) m [i]
2 → 0 and q < 2. Practically, the damping parameter is automatically selected once for all during the initialization of the iterative process so as to 5% of the values of

F (0)
m are less than or equal to .

At this stage, the explicit form of Eq. ( 11) should be given in order to highlight the regularization properties of the proposed multiplicative strategy. After some simple calculations, it readily comes:

F (k+1) m = H H H + α (k+1) W (k) -1 H H X, (15)
where α (k+1) is an adaptive regularization parameter, defined such that:

α (k+1) = X -H F (k) m 2 2 W (k) 1/2 F (k) m 2 2 . ( 16
)
The previous equation clearly shows that the proposed multiplicative strategy avoids the selection of an optimal regularization parameter at preliminary stages of the iterative process, since the adaptive regularization parameter automatically adjusts the amount of regularization throughout the resolution procedure. This particular feature constitutes the definite advantage of the proposed approach over its additive counterpart (i.e. Tikhonov-like regularization).

Choice of the initial guess

Since the resolution algorithm is iterative, another critical issue is the choice of the initial solution. Choosing a good initial guess is a key point of the convergence of the algorithm, since the functional to minimize is non-convex when q ≤ 1. The question that arises here is: What is a good initial guess? Actually, it is a coarse solution of the problem, easy to calculate, but sufficiently close to the final solution to ensure the convergence of the iterative process. Such a requirement is fulfilled by the solution of a standard Tikhonov-like regularization, that is:

F (0) m = H H H + α (0) I -1 H H X, ( 17
)
where I is the identity matrix of dimension M • N and α (0) is a rough estimate of the converged value of the adaptive regularization parameter.

The parameter α (0) has to be ideally determined without using any selection procedures or large computational efforts in order to preserve the advantage of the multiplicative strategy. To this end, one has to notice that the optimal regularization parameter is generally comprised between the smallest and the largest singular values of H H H. From this observation and a series of numerical experiments, we propose a heuristic rule for determining α (0) , that limits the computational efforts and leads to consistent identified solutions. The proposed estimation procedure is divided into three steps:

(i) Find estimates of the largest and the smallest singular values of H H H, noted σ 1 and σ n respectively.

The estimate of the largest singular value is given by the upper bound of σ 1 , namely [START_REF] Horn | Topics in Matrix Analysis[END_REF]:

σ 1 H H H = H H H ∞ H H H 1 . (18)
The estimation of the smallest singular value is obtained from σ 1 and an estimate κ of the condition number of H H namely:

σ n H H H = σ 1 H H H κ H H H . (19)
(ii) Define a set S α 0 of possible values of α (0) ∈ [ σ n , σ 1 ] using a constant logarithmic spacing to take into account the decrease of the singular values.

(iii) Choose α (0) = median (S α 0 ).

Because this estimation procedure is heuristic, it may sometimes fail to give a good starting point for the iterative process. In such a situation, it is always possible to choose α (0) as the regularization parameter picked by the L-curve. Incidentally, the computational efficiency of the overall procedure is affected in proportion to the size of the global transfer functions matrix.

Choice of the stopping criterion

The proposed iterative algorithm offers a natural definition of the stopping criterion, based on the relative variation of the adaptive regularization parameter between two successive iterations. In the present contribution, the relative variation δ of the adaptive regularization parameter is defined such that:

δ = α (k+1) -α (k) α (k) . ( 20
)
As classically done in the literature, the iterative process is stopped when the relative variation δ is less than or equal to some tolerance. Experimentally, it has been found that setting the tolerance to 10 -8 allows obtaining consistent reconstructions.

Numerical validation

In the present numerical validation, the studied structure is a simply supported steel beam with dimensions 1×0.03×0.01 m 3 excited by point force of unit amplitude from 50 Hz to 500 Hz. The coordinate of the point force, measured from the left end of the beam, is x 0 = 0.6 m. Moreover, to simulate the vibration displacement field X, a finite element model of the beam made up with 20 plane beam elements has been used. It is worth mentioning that an additive Gaussian white noise has been added to the simulated data in order to synthesize the measured vibration field. The measurement noise has been computed to obtain a signal-to-noise ration equal to 30 dB. Finally, a FE model of the structure with free boundary conditions is used to compute the transfer functions matrix H by assuming that only bending motions are measurable. In other words, the transfer functions matrices H(f j ) are dynamically condensed over the measurable dofs. The main interest in using free boundary conditions to model the dynamic behaviour of the beam is to allow the identification of the point force acting on the structure as well as reaction forces at boundaries [3,[START_REF] Renzi | [END_REF].

To better highlight the advantage of the proposed regularization term [see Eq. ( 8)], the spacefrequency (SF) regularization is compared with the corresponding frequency-by-frequency (FbF) regularization defined such that:

F m (f j ) = argmin F(f j ) X(f j ) -H(f j )F(f j ) 2 2 • F(f j ) q q , ∀ j = 1, . . . N. (21)
To compare both approaches properly, one has to assess their ability in identifying the mechanical sources acting on the structure. For this purpose, the definition of the regularization term is crucial. The analysis of the numerical test case shows that the beam is only excited by broadband point forces. In this context, it is reasonable to set q = 0.5 [see section 3.2]. The reconstructions proposed in Fig. 1 point out that both FbF and SF approaches provide a consistent identification of the reaction forces at boundaries. However, it should be noted that the FbF regularization fails in localizing the point force at one of the resonance frequencies of the beam, which is not the case with the SF strategy [see Fig. 1b]. This observation is confirmed by the analysis of the spectrum of the point force reconstructed by the FbF or the SF approach and presented in Fig. 2. Indeed, the force spectrum obtained using the FbF strategy exhibits large reconstruction errors around the resonance frequencies (94 Hz, 211 Hz and 375 Hz in the frequency range of interest). On the contrary, the reconstruction error is limited to 1.6 dB at most over all the frequency range with the SF strategy. Finally, one has to notice that the proposed regularization strategy is well adapted to solve large reconstruction problem, since the calculation of a SVD, that is generally used to compute the optimal regularization parameter, is avoided. This is all the more interesting than the calculation of the SVD can be time-consuming or even impossible for large-scale systems on a personal computer.

Conclusion

In the present study, the initial motivation was to propose a formulation of the force reconstruction problem able to fully exploit information available a priori on the type of the sources and the nature of the excitation signals. To this end, a space-frequency multiplicative regularization has been introduced. This formulation is quite flexible, since it can be used to identify localized or distributed broadband sources. Practically, the regularization problem is solved from an adapted IRLS algorithm. The numerical results are encouraging, since the proposed strategy allows obtaining consistent reconstructions even at structural resonance frequencies.

Figure 1 :

 1 Figure 1: Spatial reconstructions of the excitation field at (a) a non-resonant frequency and (b) a resonant frequency of the beam -(-) Reference, (--) Space-frequency regularization and (-• -) Frequency-by-frequency regularization

Figure 2 :

 2 Figure 2: Spectrum of the reconstructed point force -(-) Reconstructed force spectrum and (--) Reference force spectrum