M Aucejo

O De Smet

On a general Iteratively Reweighted algorithm for solving force reconstruction problems

Keywords: Linear inverse problem, Force reconstruction, Multiplicative regularization, Iteratively Reweighted algorithm

The multiplicative q -regularization has been recently introduced in structural dynamics for solving force reconstruction problems. Practically, the resolution of this regularization strategy requires the implementation of an iterative procedure. To this end, an Iteratively Reweighted Least-Squares algorithm has been originally implemented. The core idea of this algorithm

is to replace the direct resolution of the inverse problem by an equivalent iterative procedure having an explicit and unique solution at each iteration.

However, the exploitation of this very general idea allows defining other Iteratively Reweighted schemes. The present paper aims at comparing the overall performances of three particular Iteratively Reweighted algorithms for solving force reconstruction problems derived from a more general and original iterative procedure. The numerical applications proposed in this contribution highlight the ability of the considered algorithms in providing consistent regularized solutions with respect to various parameters such as the measurement noise level or the tolerance chosen to stop the iterative process.

Introduction

In structural dynamics, force identification problems are generally solved using additive regularization strategies. The most popular approach is certainly the Tikhonov regularization (a.k.a. 2 -regularization) [START_REF] Tikhonov | Regularization of incorrectly posed problems[END_REF][START_REF] Jacquelin | Force reconstruction: analysis and regularization of a deconvolution[END_REF][START_REF] Leclere | Indirect measurement of main bearing loads in an operating diesel engine[END_REF][START_REF] Liu | Reducing the impact of measurement errors when reconstructing dynamic forces[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] Mao | Experimental study of hammer impact identification on a steel cantilever beam[END_REF][START_REF] Turco | Tools for the numerical solution of inverse problems in structural mechanics: review and research perspectives[END_REF]. Although widely used, it is theoretically mainly applicable to the identification of rather smooth excitation fields or excitation signals [START_REF] Boyd | Convex optimization[END_REF], which is not a desirable effect when a localized source or an impulsive excitation signal have to be identified.

To remedy this problem, LASSO regularization (a.k.a. 1 -regularization) has been developed to promote the sparsity of the regularized solutions, while keeping the inverse problem convex. For this particular reason, it has attracted many attention in the recent years [START_REF] Ginsberg | New approach for impact detection by finding sparse solution[END_REF][START_REF] Qiao | Sparse regularization for force identification using dictionaries[END_REF][START_REF] Qiao | Impact-force sparse reconstruction from highly incomplete and inaccurate measurements[END_REF][START_REF] Pan | Moving force identification based on redundant concatenated dictionary and weighted 1 -norm regularization[END_REF]. Aucejo has generalized the previous approaches by introducing the q -regularization in the context of force reconstruction. One of the merits of this regularization strategy is to allow the resolution of non-convex sparse optimization problems, while including Tikhonov and LASSO regularization as special cases [START_REF] Aucejo | Structural source identification using a generalized Tikhonov regularization[END_REF]. However, all the procedures described above consider a global a priori on the sources to identify, meaning that poor reconstructions can be obtained, in situations where a structure excited by several sources having different spatial distributions, because the a priori, encoded in the regularization term, has to reflect a compromise between contradictory distributions. A possible solution to this issue has been proposed by Aucejo and De Smet to take advantage of prior local information available on the sources to identify [START_REF] Aucejo | Bayesian source identification using local priors[END_REF]. Practically, this approach consists in defining several identification regions for which a local q -regularization term is set. Although they have proved their efficiency, all these techniques requires a proper tuning of the regularization parameter, since it conditions the quality of the identified solution. That is why, several automatic selection methods have been developed, such as the Morozov's discrepancy principle [START_REF] Scherzer | The use of Morozov's discrepancy principle for Tikhonov regularization for solving nonlinear ill-posed problems[END_REF], the Generalized Cross Validation [START_REF] Golub | Generalized cross-validation as a method for choosing a good ridge parameter[END_REF], the Unbiased Predictive Risk Estimator [START_REF] Lin | Application of the UPRE Method to Optimal Parameter Selection for Large Scale Regularization Problems[END_REF], the Bayesian estimator [START_REF] Antoni | A Bayesian approach to sound source reconstruction: Optimal basis, regularization, and focusing[END_REF][START_REF] Pereira | Acoustic imaging in enclosed spaces[END_REF] or the L-curve principle [START_REF] Hansen | Rank-Deficient and Discrete Ill-Posed Problems: Numerical Aspects of Linear Inversion[END_REF]. It can be noted that the aforementioned methods generally requires intensive computations, since they are based either on the calculation of the root of some equation, the minimization of some functional or the determination of the maximum curvature of a certain curve. For this particular reason, Aucejo and De Smet have recently developed the multiplicative q -regularization for reconstructing mechanical sources [START_REF] Aucejo | A multiplicative regularization for force reconstruction[END_REF]. This regularization strategy, originally proposed by Van den Berg et al. [START_REF] Van Den Berg | Multiplicative regularization for constrast profile inversion[END_REF], consists in including the regularization term as a multiplicative constraint in the formulation of the inverse problem, which alleviates the need for the preliminary definition of any optimal regularization parameter. From a practical standpoint, the multiplicative regularization is solved in an iterative manner from an Iteratively Reweighted (IR) Least-Squares algorithm, a.k.a. IR 2 algorithm. The core idea of this algorithm is to replace the direct resolution of the inverse problem by an equivalent iterative process having an explicit and unique solution at each iteration by recasting the q -regularization term into a weighted 2 -regularization term. However, the IR 2 algorithm is not the sole Iteratively Reweighted procedure that can be implemented, since IR 1 algorithms have been proposed by several researchers for solving sparse reconstruction problems [START_REF] Candès | Enhancing sparsity by Reweithed 1 minization[END_REF][START_REF] Chartrand | Iteratively reweighted algorithms for compressive sensing[END_REF][START_REF] Lyu | A comparison of typical p minimization algorithms[END_REF]. As an example, Chartrand and Yin have compared in Ref. [START_REF] Chartrand | Iteratively reweighted algorithms for compressive sensing[END_REF] the performances of IR 2 and IR 1 algorithms. Because the difference between the resulting regularized solutions is negligible, they conclude that IR 2 is the better approach due to its computational simplicity.

In the present paper, an original and general IR p algorithm is introduced for solving the force identification problems. More specifically, it is intended to compare the overall performances of the IR p algorithm for three particular values of p. To properly highlight the main features of the proposed approach, this article is divided into three parts. For the sake of completeness, the multiplicative regularization proposed in Ref. [START_REF] Aucejo | A multiplicative regularization for force reconstruction[END_REF] is first recalled in section 2. Then, section 3 details the implementation of the proposed IR p algorithm. Finally, the overall performances of the IR p algorithm is illustrated numerically in section 4, where sparse and locally sparse sources are reconstructed in the frequency domain. In particular, the numerical applications proposed in this contribution aims at comparing the ability of the three particular algorithms, deriving from the IR p scheme, in providing consistent regularized solutions with respect to various parameters such as the measurement noise level or the tolerance chosen to stop the iterative process.

From the conclusions resulting from the considered applications, a set of recommendation is proposed to properly define the main parameters of the IR p algorithm.

Multiplicative regularization

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 in the frequency domain is completely determined by the transfer functions matrix H, relating the vibration field X to the excitation field F at a particular frequency, so that:

X = HF + N, (1)
where N is an additive noise accounting for measurement and modeling errors.

Reconstructing excitation forces is far from an easy task, since the direct inversion of the previous relation leads generally to inconsistent results, because the identified solution is biased by the term H -1 N. A classical approach to circumvent this problem consists in including in the formulation of the inverse problem some prior information on the forces to identify [START_REF] Tikhonov | Regularization of incorrectly posed problems[END_REF][START_REF] Aucejo | Bayesian source identification using local priors[END_REF][START_REF] Tibshirani | Regression shrinkage and selection via the Lasso[END_REF]. To this end, we have recently proposed to apply the multiplicative regularization, which is mathematically expressed as:

F = argmin F\{0} F(X -HF) • R(F), (2)
where:

• F(X -HF) is the data-fidelity term which controls the a priori on the noise corrupting the data [START_REF] Green | Statistics of Images, the TV Algorithm of Rudin-Osher-Fatemi for Image Denoising and an Improved Denoising Algorithm[END_REF][START_REF] Nikolova | A variational approach to remove outliers and impulse noise[END_REF][START_REF] Le | A variational approach to reconstructing images corrupted by Poisson noise[END_REF].

• R(F) is the regularization term that encodes prior information on the excitation field F [START_REF] Tibshirani | Regression shrinkage and selection via the Lasso[END_REF][START_REF] Chambolle | Total Variation Minimization and a Class of Binary MRF Models[END_REF][START_REF] Hintermuller | Nonconvex TV q -models in image restoration: Analysis and a trust-region regularization based superlineraly convergennt solver[END_REF].

At this stage, it is clear that the quality of the reconstructed excitation field is conditioned to the adequacy of the data-fidelity and regularization terms with the actual noise and the actual source characteristics. Regarding the definition of the data-fidelity term, it is classically assumed that the vibration field X is corrupted by an additive Gaussian white noise. As a result, the data-fidelity term can be defined such that [START_REF] Antoni | A Bayesian approach to sound source reconstruction: Optimal basis, regularization, and focusing[END_REF][START_REF] Faure | Empirical and fully Bayesian approaches for the identification of vibration sources from transverse displacement measurements[END_REF][START_REF] Aucejo | On a full Bayesian inference for force reconstruction problems[END_REF]:

F(X -HF) = X -HF 2 2 . (3)
It should be noted here that if the actual noise is additive but not strictly Gaussian (e.g. impulsive noise), then previous data-fidelity term can be more generally written as [START_REF] Aucejo | Structural source identification using a generalized Tikhonov regularization[END_REF][START_REF] Nikolova | A variational approach to remove outliers and impulse noise[END_REF]:

F(X -HF) = X -HF s s , (4)
• s is the s -norm or quasi-norm. The introduction of such a data-fidelity term in the formulation of the regularization problem only leads to minor changes (not described in this paper) in the algorithm described in section 3.

The definition of the regularization term requires more attention, since it has to reflect one's prior knowledge of the sources to identify. Generally, forces, whose nature and locations are practically roughly known, can simultaneously excite a structure. That is why, it is supposed that the structure is excited in R different regions by local excitation fields F r of various types (localized or distributed, for instance). This assumption allows expressing the regularization term R(F) as a sum of local regularization terms, that is:

R(F) = R r=1 R(F r ), (5)
where R(F r ) is the regularization term associated to the zone r.

In Ref. [START_REF] Aucejo | Bayesian source identification using local priors[END_REF][START_REF] Aucejo | A multiplicative regularization for force reconstruction[END_REF], the local regularization terms are written:

R(F r ) = L r F r qr qr , (6)
where q r is the norm parameter belonging to R + * , while L r is a smoothing operator controlling the regularity of the solution in the zone r. Some examples of the structure of the smoothing operator can be found in Refs. [START_REF] Aucejo | Bayesian source identification using local priors[END_REF][START_REF] Aucejo | A space-frequency multiplicative regularization for force reconstruction problems[END_REF].

Finally, by introducing Eqs. ( 6) and (3) into Eq. ( 2), one obtains the general expression of the proposed multiplicative regularization:

F = argmin F\{0} X -HF 2 2 • R r=1 L r F r qr qr . (7)

IR p algorithm

The resolution of the reconstruction problem defined by the multiplicative regularization presented in section 2 requires the implementation of an iterative procedure, since the solution of a reconstruction problem involving qr -norms has generally no closed-form expression.

In Ref. [START_REF] Aucejo | A multiplicative regularization for force reconstruction[END_REF], the present multiplicative regularization is solved by an IR 2 algorithm. In this section, we propose to generalize this approach by deriving an IR p algorithm, which allows a greater resolution flexibility thanks to the choice of the value of norm parameter p.

General principle

The basic idea behind the proposed IR p algorithm is to replace the direct resolution of the minimization problem, given by Eq. ( 7), by an equivalent iterative process obtained by recasting the qr -regularization term into a weighted p , where p ∈ R + * .

The direct application of this idea to the proposed regularization problem leads to define the estimated force vector

F (k+1)
at iteration k+1 of the iterative procedure as the solution of the following minimization problem:

F (k+1) = argmin F\{0} X -HF 2 2 • R r=1 W (k) r 1/p L r F r p p , (8)
where W (k) r is a local weighting definite positive diagonal matrix.

As classically done in IR-type algorithm, the weighting matrices W (k) r are defined such that:

W (k) r = diag w (k) r,1 , . . . , w (k)
r,i , . . . , w

(k) r,Nr (9)
with:

w (k) r,i = max (k) r , f (k) i qr-p , (10)
where N r is the number of identification points in the zone r, f

(k) i is the i th component of the vector L r F (k)
r and r is a small real positive number acting as a damping parameter. It allows avoiding infinite weights when f (k) i → 0 and q r < p. Practically, the damping parameter is automatically from the cumulative histogram of L r F (k) r . More precisely, its value is chosen so that 5% of the values of L r F (k) r are less than or equal to r [START_REF] Aucejo | Structural source identification using a generalized Tikhonov regularization[END_REF][START_REF] Rodriguez | Efficient Minimization Method for a Generalized Total Variation Functional[END_REF].

For implementation convenience, Eq. ( 8) is expressed in a more compact form by introducing two global matrices, W (k) and L, defined as follows:

W (k) = diag W (k) 1 , . . . , W (k) R and L = diag (L 1 , . . . , L R ) . (11)
These global matrices allow actually writing the minimization problem given by Eq. ( 8) under the following generic form:

F (k+1) = argmin F\{0} X -HF 2 2 • W (k) 1/p LF p p . ( 12
)
The operational form of the previous minimization problem is obtained by applying the first-order optimality condition. In doing so, it readily comes that:

F (k+1) = H H H + α (k+1) L (k)H Ω (k) L (k) -1 H H X. (13)
where α (k+1) is the adaptive regularization parameter, defined such that:

α (k+1) = X -H F (k) 2 2 L (k) F (k) p p with L (k) = W (k) 1/p L, (14)
while Ω (k) is a diagonal weighting matrix whose components ω (k) i are given by:

ω (k) i = p 2 max (k) , f (k) i p-2 , (15)
where

f (k) i is the i th component of the vector L (k) F (k)
and is the damping parameter that allows avoiding infinite weights when

f (k) i → 0 and p < 2.
Its value is chosen so that 5% of the values of L

(k) F (k)
are less than or equal to (k) .

Choice of the initial guess and stopping criterion

The proposed resolution algorithm being iterative, the present section aims at introducing the choices of the initial guess and the stopping criterion.

Choice of the initial guess

It is clear from section 2 that the proposed multiplicative formulation is generally non-convex, meaning that the existence of a unique minimizer is not guaranteed. Consequently, the choice of a good initial guess

F (0)
is crucial for a successful reconstruction.

A good initial guess can be defined as a coarse solution of the problem, easy to compute, but sufficiently close to the solution to ensure the convergence of the iterative process. Such a requirement is fulfilled by the solution of the standard Tikhonov-like regularization, that is:

F (0) = H H H + α (0) L H L -1 H H X, (16)
where α (0) is a rough estimate of the converged value of the adaptive regularization parameter or, equivalently, the value of the regularization parameter λ picked by the L-curve method.

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. This task proves to be difficult in practice, because the order of magnitude of the optimal regularization parameter is a priori unknown from the data only. To avoid using automatic selection procedure, some heuristics can be used, as the one proposed in

Ref. [START_REF] Aucejo | A multiplicative regularization for force reconstruction[END_REF]. Another possibility originates from studies on 1 -regularization [START_REF] Kim | An Interior-Point Method for Large-Scale 1 -Regularized Least Squares[END_REF]. Practically, the initial adaptive regularization parameter is defined as a fraction of regularization parameter above which the optimal solution of the 1 -regularization is the zero vector. From a series of numerical experiments, its value has been chosen such that:

α (0) = 10 -3 2 H H X ∞ . ( 17
)
The influence of the heuristic procedure used for computing the initial solution will be investigated in section 4.

Choice of the stopping criterion

The proposed IR p 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 paper, the relative variation δ of the adaptive regularization parameter is defined such that:

δ = α (k+1) -α (k) α (k) . ( 18
)
As classically done in the literature, the iterative process is stopped when the relative variation δ is less than or equal to some tolerance. The influence of the tolerance chosen to stop the iterative process will be examined in section 4.

Summary and comments

To clearly highlight each step of the proposed IR p algorithm, an overview of the resolution procedure is given in Fig. 1. From the analysis of this figure and the theoretical derivation proposed in the previous section, several comments can be made. First of all, the proposed IR p algorithm includes the IRLS algorithm (p = 2) as a special case. Furthermore, when only one identification region is considered, the case q = p leads to the iterative procedure that would be obtained by applying directly the first-order optimality condition to Eq. ( 7), which is another example of IR algorithm. Consequently, the IR p can actually be considered as a generalization of a certain class of IR algorithms.

Initialization: Compute α (0) and F

, Set p and q

Compute W (k) and Ω (k)

Convergence ?

Yes No F Compute α (k+1) and F (k+1)
Eqs. ( 16) and ( 17)

Eqs. ( 10) and ( 15)

Eqs. ( 13) and ( 14)

Eq. ( 18)

Figure 1: Overview of the IR p algorithm Furthermore, it is worth mentioning that the proposed IR p can be directly applied to additive regularization strategies, since one just has to replace the adaptive regularization parameter α (k) by the regularization parameter λ (k) . The only difficulty in this case is to properly estimate the optimal value of the regularization parameter at each iteration. For the IR 2 algorithm (i.e. IRLS algorithm), the regularization parameter can be estimated from existing automatic selection procedures. On the contrary, for the case p = 2, this problem seems to remain an open issue.

Finally, one should notice that IR algorithms have attracted the attention of the scientific community and remains an active research topic regarding, in particular, the convergence and stability of such algorithms. To the best of our knowledge, some theoretical results related to the convergence and stability of IR p algorithm can be found for p = 2 and p = 1 in Refs. [START_REF] Ba | Convergence and Stability of Iteratively Re-weighted Least Squares Algorithms[END_REF][START_REF] Chen | Convergence of the reweigthed 1 minimization algorithm for 2p minimization[END_REF][START_REF] Ochs | On Iteratively Reweighted Algorithms for Nonsmooth Nonconvex Optimization in Computer Vision[END_REF][START_REF] Lu | Iterative reweighted minimization methods for p regularized unconstrained nonlinear programming[END_REF][START_REF] Wipf | Iterative reweighted and methods for finding sparse solutions[END_REF].

For the case p ≥ 1, convergence results can be found in Ref. [START_REF] Lu | Iterative reweighted minimization methods for p regularized unconstrained nonlinear programming[END_REF]. Unfortunately, for p ≤ 1, it seems that no theoretical results demonstrating the convergence of the IR p algorithm currently exist and this question remains an open issue. However, the results presented in section 4 show that the proposed algorithm seems effective for solving force reconstruction problems.

Numerical applications

This section intends to investigate, through two frequency-domain applications, the behavior and the performances of the IR p algorithm with respect to several parameters of the problem such as the measurement noise level, the considered initial solution or the tolerance used to judge the convergence of the algorithm. In the present numerical applications, only the reconstruction of the spatial distributions of sources associated to a harmonic point force exciting the studied structures at a particular frequency is considered for p = {0.5, 1, 2}. This is not a limitation of the method, since the force spectrum or history can be reconstructed for any value of p from the proposed approach provided that the system matrix H is established accordingly. Finally, it is also worth noting that all the reconstructions are performed outside resonance frequencies. Indeed, at resonance frequencies, the resolution of the inverse problem is known to be extremely difficult. Because the system matrix is singular, the uniqueness of the solution is not ensured and the identified solution is generally not satisfying whatever the regularization approach used.

1D structure

In the present application, the studied structure is a simply-supported steel beam with dimensions 1 m×0.03 m×0.01 m excited by a harmonic unit point force at 350 Hz. The coordinate of the point force, measured from the left end of the beam, is x 0 = 0.6 m.

In the following sections, the reconstructions are performed from the displacement field of the simply-supported beam X computed at 350 Hz, assuming a structural damping ratio of 0.01. To synthesize the measured vibration field, two numerical steps have been implemented. First, a finite element model of the beam made up with 20 plane beam elements has been used to compute the exact displacement field. Then, a Gaussian white noise with a prescribed SNR has been added to the exact data to simulate measurement errors, related to the transducers quality. Regarding the transfer functions matrix H of the system, it has been computed from a FE model of the free-free beam by assuming that only bending motions are measurable. In other words, the transfer functions matrix have been dynamically condensed over the measurable dofs, meaning that the reconstruction problem is fully determined. The main interest in using free boundary conditions to model the dynamic behavior of the beam is to allow the identification of external excitations acting on the structure as well as reaction forces at boundaries [START_REF] Renzi | Vibratory source identification by using the finite element model of a subdomain of a flexural beam[END_REF].

To finalize the definition of the multiplicative regularization, it remains to determine the number of identification regions and the value of the related norm parameters. The description of the mechanical problem suggests that only one identification region needs to be defined. To properly choose the value of the related norm parameter q, it should be noted that smooth solutions are promoted for q = 2 [8, 42], while sparse excitation fields are favored with q ≤ 1 [START_REF] Aucejo | Structural source identification using a generalized Tikhonov regularization[END_REF][START_REF] Aucejo | Bayesian source identification using local priors[END_REF][START_REF] Tibshirani | Regression shrinkage and selection via the Lasso[END_REF][START_REF] Chartrand | Nonconvex regularization for image segmentation[END_REF][START_REF] Fu | Efficient minimization methods of mixed 2 -1 and 1 -1 norms for image restoration[END_REF][START_REF] Grasmair | Non-convex sparse regularization[END_REF][START_REF] Chu | Robust Bayesian superresolution approach via sparsity enforcing a priori for near-field aeroacoustic source imaging[END_REF][START_REF] Li | A hierarchical Bayesian method for vibration-based time domain force reconstruction problems[END_REF]. In the light of this observation, it is reasonable to set q = 0.5, because the target excitation field has theoretically only three non-zero components.

Furthermore, to quantify the accuracy of the identified solutions with respect to the considered operating conditions, the relative error (RE) and the peak error (PE) are evaluated. Formally, the relative error is a global indicator of the reconstruction quality, which is formally defined by the relation:

RE = F ref -F 1 F ref 1 . (19)
Contrary to the relative error, the peak error is a local indicator that allows assessing the reconstruction quality of the point force amplitude. Mathematically, it is defined such that:

PE = F p -F ref p F ref p , (20)
where F ref p is the point force amplitude associated to the reference force vector F ref , while F p is the point force amplitude associated to the identified solution F.

Finally, in absence of contradictory information, the SNR is set to 35 dB, the tolerance used for stopping the iterative is set to 10 -4 , while the initial solution

F (0)
is computed from the initial adaptive regularization parameter α (0) estimated using the heuristic procedure described in Ref. [START_REF] Aucejo | A multiplicative regularization for force reconstruction[END_REF].

Sensitivity of the algorithm to the initial solution

In section 3.2.1, it has be noted that only the definition of the initial adaptive parameter α (0) is required to initialize the algorithm, since the initial force vector F (0) is computed accordingly from Eq. ( 16). In this context, it is interesting to assess the behavior of the IR p algorithm with respect to the heuristics used the estimate the value of the initial adaptive regularization parameter. In the following, H 1 refers to the heuristic procedure described in Ref. [START_REF] Aucejo | A multiplicative regularization for force reconstruction[END_REF], while H 2 designates the estimation provided by Eq. [START_REF] Lin | Application of the UPRE Method to Optimal Parameter Selection for Large Scale Regularization Problems[END_REF]. Fig. 2 presents the excitation fields reconstructed from the IR p algorithm for p = {0.5, 1, 2} depending on the heuristics used to estimate the initial adaptive regularization parameter. The analysis of this figure suggests that all the identified excitation fields are qualitatively in good agreement with the reference one. This observation is confirmed by the results listed in Table 1.

The main conclusion of this study is that the both heuristics considered in this paper can be indistinctly used for initializing the IR p algorithm. Finally, one has to mention that the same conclusions have been drawn for the 2D structure presented in section 4.2. However, these results have not been reproduced here for the sake of conciseness. The results presented in Fig. 3 and listed in Table 2 shows that for high and moderate SNRs, namely 35 dB and 25 dB respectively, the three considered IR algorithms perform equally. However, in case of low SNR, typically 15 dB, all the algorithms greatly underestimate the point force amplitude.

In this case, however, it should be noted that IR 0.5 and IR 1 give similar reconstructions, while IR 2 exhibits a larger peak error. For the sake of completeness, the results obtained when α (0) is estimated from H 2 are presented in Appendix A. Our results suggest that the standard IR 2 algorithm is less effective than the other IR algorithms considered when SNR gets lower for reconstructing sparse excitation fields. As a side note, it should however be noted that, below 10 dB, all the versions of all the IR algorithms considered fail in providing consistent reconstruction, since the solution obtained is the zero vector.

Influence of the choice of the convergence tolerance

In previous sections, all the reconstructions have been performed for a tolerance set to 10 -4 . However, as any iterative procedure, the solution accuracy is partly conditioned to the value of the tolerance used to stop the iterations. For this reason, it is interesting to assess the performance of the IR p algorithm with respect to the value of tolerance defined by the user.

Table 3 gathers the results obtained with the different IR algorithms for tolerance values ranging from 10 -2 to 10 -8 and α (0) estimated from H 1 . It worth noting here that the three considered IR algorithms give similar regularized solutions whatever the value of the tolerance. However, it is important to note that IR 1 is subject to convergence issues, as indicated by the number of iterations made by the algorithm, which dramatically increases as the tolerance decreases. It should be noticed here that the convergence difficulty of IR 1 observed in the case of a high SNR value (here 35 dB) disappears when the SNR gets lower [see Appendix B]. All these results suggest to define the tolerance so as to obtain a reasonable compromise between the solution accuracy and the computational efficiency.

Recommendations

In the light of the previous results, the following recommendations can be made to wisely use the IR p algorithm:

1. Whatever the value of p, the proposed IR p algorithm is almost insensitive to the choice of the heuristics considered in this paper to estimate the initial adaptive regularization parameter α (0) . Consequently, it is recommended to apply the procedure H 2 , since it is requires fewer operations than H 1 to obtain a reasonable initial adaptive regularization parameter and is, de facto, computationally more effective.

2. When the target excitation field is sparse, one suggests to use the IR p algorithm for p ≤ 1, because it leads to more consistent reconstructions in case of very noisy data. prone to think that it is not a wise practice, especially when dealing with non-convex inverse problems. Nevertheless, one suggest to set the tolerance of the iterative process to 10 -4 in order to define a reasonable compromise between the solution accuracy and the computational efficiency.

2D structure

This second application aims at confirming the conclusions drawn in the previous section and potentially review the proposed recommendations by implementing the IR p algorithm on a more complex case study. However, to limit the size of the paper, the whole study carried out in the previous section is not reproduced here. This application is actually focused on the most critical points that allows completing the analysis of the IR p algorithm.

Problem description

In the present application, one seeks to identify a harmonic point force of unit amplitude acting at 350 Hz on a thin simply supported steel plate with dimensions 0.6 m×0.4 m×0.005 m. The coordinates of the point force, measured from the lower left corner of the plate, are (x 0 , y 0 ) = (0.42 m, 0.25 m).

Practically, this configuration allows studying the influence of the definition of local regularization terms, since the present excitation field exhibits two different spatial distributions over the structure, namely a smooth distribution of the reaction forces at boundaries and a singular distribution around the location of the point force.

In this example, the reconstructions are performed from the velocity field X computed at 350 Hz, assuming a structural damping ratio equal to 0.01.

To properly simulate experimental measurements, the exact vibration displacement field X exact is first computed from a FE mesh of the plate made up with 187 shell elements, assuming that only bending motions are measurable. Then, a Gaussian white noise, with a prescribed SNR is added to the exact data to simulate measurement errors.

Finally, the transfer functions matrix H is computed, as in section 4.1, from the FE model of the structure with free boundary conditions, assuming that only normal velocities are measurable, meaning that the reconstruction problem is fully determined.

Application

To numerically validate any force reconstruction strategy, it is first necessary to define the reference force vector F ref that could serve as a proper benchmark. This reference force vector is computed from the transfer functions matrix H and the exact displacement field X exact thanks to the following relation:

F ref = H -1 X exact . (21)
As expected, the reference force vector at 350 Hz corresponds to the description of the test case given in the very beginning of this section, since it exhibits smooth reaction forces at boundaries of the plate as well as a unit point force F 0 at (x 0 , y 0 ) = (0.42 m, 0.25 m).

Consequently, the analysis of the reference force vector suggests that two identification regions can be defined: (i) a central region associated to the shape parameter q 1 and containing the point force only and (ii) a region associated to the shape parameter q 2 and corresponding to the boundaries of the plate [see Fig. 5]. Given this partition of the structure and prior information on the sources to identify associated to each region, it is reasonable to set (q 1 , q 2 ) = (0.5, 2) [START_REF] Boyd | Convex optimization[END_REF][START_REF] Aucejo | A multiplicative regularization for force reconstruction[END_REF].

Following the recommendations given in section 4.1.4, the proposed reconstructions are performed using α (0) estimated from the heuristics H with the reference one. In particular, it should be pointed out that IR 0.5 variant performs better for quantifying the point force amplitude. However, it is not the better approach in this case to properly reconstruct the reaction forces, since IR 1 and IR 2 provide better solutions in this case. Perhaps more interesting is the behavior of the IR p algorithm in case of low SNR. Indeed, when the SNR is set to 15 dB, the three IR algorithms behave rather differently as observed in Fig. 7 and Table 5. In the present situation, only the IR 0.5 and IR 1 algorithms are able to identify the point force quite satisfactorily, while all the IR algorithms provide a reasonable reconstruction of the reaction forces. All in all, the global analysis of obtained results suggests that the IR 1 provide the better compromise between the quantification of the point force amplitude and the overall accuracy of the reconstructed excitation field. Consequently, in the light of the results presented in this section, it seems that the IR p algorithm for p ≤ 1, and especially, p = 1, is the better option, when the excitation field to identify is sparse or locally sparse. This observation suggests that the choice of p is related to the overall distribution of the excitation sources over the structure, i.e. mainly sparse or distributed.

Conclusion

In the present paper, a general and original IR p algorithm has been introduced to solve a force reconstruction problem formulated as a multiplicative The results listed in Table B.7 are in line with those presented in Table 3.

The major noticeable difference lies in the number of iterations made by the algorithm to reach the convergence. Indeed, while IR 1 exhibits convergence issues for high SNR values, it is no more the case when the SNR gets lower.

Figure 2 :

 2 Figure 2: Reconstructed force vector F at 350 Hz with respect to the heuristic used to estimate the initial adaptive regularization parameter α (0) -(a) H 1 and (b) H 2 -(-) Reference, (--) IR 2 , (-• -) IR 1 and (• • • ) IR 0.5

Figure 3 :

 3 Figure 3: Reconstructed force vector F at 350 Hz from a vibration field having an SNR equal to (a) 25 dB and (b) 15 dB -(-) Reference, (--) IR 2 , (-• -) IR 1 and (• • • ) IR 0.5

Figure 4 :

 4 Figure 4: Reference force vector F ref at 350 Hz

2 and a tolerance set to 10 - 4 .Figure 5 : 2 (

 10452 Figure 5: Definition of the identification regions -(•) region 1 (Point force), (×) region 2 (Reaction forces) and (♦) location of the point force

Figure 6 :

 6 Figure 6: Reconstructed force vector F at 350 Hz from IR p algorithm for a SNR equal to 25 dB. (a) Surface plot -IR 2 , (b) Surface plot -IR 1 , (c) Surface plot -IR 0.5 and (d) Section view at y 0 = 0.25 m -(-) Reference, (--) IR 2 , (-• -) IR 1 and (• • • ) IR 0.5

Figure 7 :

 7 Figure 7: Reconstructed force vector F at 350 Hz from IR p algorithm for a SNR equal to 15 dB. (a) Surface plot -IR 2 , (b) Surface plot -IR 1 , (c) Surface plot -IR 0.5 and (d) Section view at y 0 = 0.25 m -(-) Reference, (--) IR 2 , (-• -) IR 1 and (• • • ) IR 0.5

Table 1 :

 1 Performances of the IR p algorithm with respect to the initial estimation of α (0) -N it : Number of iterations of the algorithm Procedure Indicator IR 0.5 IR 1 IR 2After applying the IR p algorithm for a relatively high SNR value, it is legitimate to wonder whether the quality of the identified solutions remains the same when the SNR gets lower whatever the value of p. To answer this question, the IR p is applied to vibration data having SNR equal to 15 dB, 25 dB and 35 dB. For this numerical experiment, α (0) is estimated from H 1 .

		PE (%)	-2.49 -2.79 -2.47
	H 1	RE (%)	1.33	1.91 1.20
		N it	16	19	9
		PE (%)	-2.49 -2.69 -2.47
	H 2	RE (%)	1.33	1.79 1.19
		N it	27	16	9

Table 2 :

 2 Performances of the IR p algorithm with respect to the measurement noise level corrupting the data -N it : Number of iterations of the algorithm initialized using H 1

	SNR (dB) Indicator IR 0.5 IR 1 IR 2
		PE (%)	-2.49	-2.79 -2.47
	35	RE (%)	1.33	1.91	2.47
		N it	16	19	9
		PE (%)	-7.62	-7.69 -7.87
	25	RE (%)	3.70	3.75	3.83
		N it	13	20	9
		PE (%)	-25.11 -25.55 -34.10
	15	RE (%)	12.07 12.31 15.32
		N it	14	15	24

Table 3 :

 3 Performances of the IR p algorithm with respect to the tolerance fixed to stop the iterative process -N it : Number of iterations of the algorithm initialized using H 1 Tolerance Indicator IR 0.5 IR 1 IR 2 Our results suggest that a large tolerance can be used, but we are

		PE (%)	-2.49 -2.84 -2.47
	10 -2	RE (%)	1.33	1.98 1.20
		N it	14	11	8
		PE (%)	-2.49 -2.79 -2.47
	10 -4	RE (%)	1.33	1.91 1.20
		N it	16	19	9
		PE (%)	-2.49 -2.73 -2.47
	10 -6	RE (%)	1.33	1.93 1.20
		N it	18	195	10
		PE (%)	-2.49 -2.73 -2.47
	10 -8	RE (%)	1.33	1.93 1.20
		N it	21	200	11
		21

Table 4 :

 4 Performances of the IR p algorithm for a SNR equal to 25 dB -N it : Number of iterations of the algorithm Indicator IR 0.5 IR 1 IR 2

	PE (%)	-4.33 -5.04 -5.10
	RE (%)	24.18 19.05 19.02

Table 5 :

 5 Performances of the IR p algorithm for a SNR equal to 15 dB -N it : Number of iterations of the algorithm Indicator IR 0.5 IR 1 IR 2

	PE (%)	-9.25 -11.70 -99.99
	RE (%)	30.39 26.33 34.93
	N it	21	31	15

  Table A.6: Performances of the IR p algorithm with respect to the measurement noise level corrupting the data when choosing α (0) from H 2 -N it : Number of iterations of the algorithm SNR (dB) Indicator IR 0.5 IR 1 IR 2

	to 15 dB.
		PE (%)	-2.49	-2.69 -2.47
	35	RE (%)	1.33	1.79	1.19
		N it	27	16	9
		PE (%)	-7.62	-7.70 -7.98
	25	RE (%)	3.70	3.79	4.11
		N it	13	32	9
		PE (%)	-25.11 -25.55 -33.72
	15	RE (%)	12.07 12.30 14.83
		N it	14	15	25

Table B .

 B 7: Performances of the IR p algorithm with respect to the tolerance fixed to stop the iterative process for SNR = 15 dB -N it : Number of iterations of the algorithm Tolerance Indicator IR 0.5 IR 1 IR 2

		PE (%)	-2.49	-2.79 -2.47
	10 -2	RE (%)	1.33	1.91	1.20
		N it	10	11	8
		PE (%)	-25.11 -25.55 -34.10
	10 -4	RE (%)	12.07 12.31 15.32
		N it	14	15	24
		PE (%)	-25.11 -25.55 -34.15
	10 -6	RE (%)	12.07 12.31 15.33
		N it	17	19	43
		PE (%)	-25.11 -25.55 -34.15
	10 -8	RE (%)	12.07 12.31 15.33
		N it	31	23	62

Appendix A. Influence of the measurement noise level when choosing α (0) from H 2

In this appendix, we propose to the reader a presentation of the results associated to the numerical experiments described in section 4, where the performances of the IR p algorithm with respect to the noise level is examined.

In particular, this appendix presents the results obtained when the initial adaptive regularization parameter α (0) is estimated from H 2 . The comparison of the results listed in Tables 2 andA