Mathieu Aucejo
email: mathieu.aucejo@cnam.fr

Olivier De Smet

O De Smet

J.-F Deü

Practical issues on the applicability of Kalman filtering for reconstructing mechanical sources in structural dynamics

Keywords: Inverse problem, Force reconstruction, Augmented Kalman Filter, Structural dynamics

Practical issues on the applicability of Kalman filtering

for reconstructing mechanical sources in structural dynamics

Introduction

Joint input-state estimation based on Kalman-type filtering is a trending topic in the system identification community as evidenced by the numerous papers published in the last two decades. Nowadays, two classes of filters are available for dealing with such a problem. The first one considers the inputstate estimation problem in a separate manner. In this class, the input-state estimation can be achieved either by deriving an optimal recursive filter in an unbiased minimium-variance sense [START_REF] Borisov | Optimal filtering in stochastic discrete-time systems with unknown inputs[END_REF][START_REF] Darouach | Kalman filtering with unknown inputs via optimal state estimation of singular systems[END_REF][START_REF] Chen | Optimal filtering and robust fault diagnosis of stochastic systems with unknown disturbances[END_REF][START_REF] Darouach | Unbiased minimum variance estimation for systems with unknown exogenous inputs[END_REF][START_REF] Hou | Optimal filtering for systems with unknown inputs[END_REF][START_REF] Valcher | State observers for discrete-time linear systems with unknown inputs[END_REF][START_REF] Hsieh | Robust two-stage kalman filters for systems with unknown inputs[END_REF][START_REF] Darouach | Extension of minimum variance estimation for systems with unknown inputs[END_REF][START_REF] Hsieh | Unbiased minimum-variance input and state estimation for systems with unknown inputs: A system reformation approach[END_REF][START_REF] Maes | Joint input-state estimation in structural dynamics[END_REF] or by combining two separate Kalman filters as in the Dual Kalman Filter proposed by Eftekhar Azam et al. [START_REF] Azam | A dual kalman filter approach for state estimation via output-only acceleration measurements[END_REF][START_REF] Azam | Experimental validation of the Kalman-type filters for online and real-time state and input estimation[END_REF]. In the second class of filters, the input-state estimation is performed simultaneously by considering an augmented state vector composed by both the state and input vectors [START_REF] Hsieh | Robust two-stage kalman filters for systems with unknown inputs[END_REF]. In doing so, the resulting input-state estimation problem is solved from a conventional Kalman filter algorithm [START_REF] Lourens | An augmented kalman filter for force identification in structural dynamics[END_REF][START_REF] Naets | Stable force identification in structural dynamics using kalman filtering and dummy-measurements[END_REF][START_REF] Risaliti | A state-input estimation approach for force identification on an automotive suspension component[END_REF].

In structural dynamics, Lourens et al. [START_REF] Lourens | Joint input-response estimation for structural systems based on reduced-order models and vibration data from a limited number of sensors[END_REF], Naets et al. [START_REF] Naets | Stable force identification in structural dynamics using kalman filtering and dummy-measurements[END_REF] and, more particularly, Maes et al. [START_REF] Maes | Design of sensor networks for instantaneous inversion of modally reduced order models in structural dynamics[END_REF] have produced a set of guidelines derived from control theory that guarantee the identifiability, stability and uniqueness of estimated quantities. From the practitioner standpoint, however, these conditions, related to model and sensor configurations, are necessary, but not sufficient to ensure the success of the joint estimation in the presence of noise, because they establish the feasibility of the inversion process for a given experimental set-up from an algorithmic perspective only. Consequently, some practical issues remain partially documented or not clearly illustrated, even if a careful reading and analysis of the extensive literature available provides some rather useful indications.

At this stage, it is important to mention that the Kalman-type strategies evoked above are the most commonly used currently in the structural dynamics community to jointly estimate the state and the external excitation sources applied to a structure. Of course, more sophisticated and state-ofthe-art approaches, such as cubature Kalman Filters [START_REF] Arasaratnam | Cubature kalman smoothers[END_REF] or NIRK-based extended Kalman filters [START_REF] Kulikova | Nirk-based accurate continuous-discrete 60 extended kalman filters for estimating continuous-time stochastic target tracking models[END_REF], could be used for that purpose. This lies, however, outside the scope of this paper, which intends to assess the applicability of the Kalman-type filters to deal with force reconstruction problems in some particular practical situations. This paper is more particularly focused on the applicability of the Augmented Kalman Filter (AKF) proposed by Lourens et al. [START_REF] Lourens | An augmented kalman filter for force identification in structural dynamics[END_REF]. This choice is made because it allows performing joint input state estimation simultaneously from a standard Kalman filter scheme, which has solid theoretical foundations. Because the present paper is focused on purely input estimation problems, this means that state estimation is not of particular interest here. In this regard, several issues are of concern when dealing with the practical implementation of such an identification strategy. More specifically, it is proposed to thoroughly illustrate the behavior of the AKF regarding the type of kinematic data used as output data, the possibility of real-time reconstruction or the influence of sensors and reconstruction con-figurations. To this end, the state-space representation used to implement the AKF is based on a state-space representation built from the generalizedα method, which is unconditionnally stable and second order accurate. To clearly introduce the theoretical background of the filtering strategy and illustrate each of the practical issues evoked previously, the paper is divided into three parts. In section 2, the derivation of the state-space representation of a dynamical system based on the generalized-α method is introduced, as well as the main equation of the AKF. Section 3 is dedicated to the tuning of the Kalman Filter. Indeed, it has been shown in [START_REF] Azam | A dual kalman filter approach for state estimation via output-only acceleration measurements[END_REF][START_REF] Petersen | Estimation of the fullfield dynamic response of a floating bridge using kalman-type filtering algorithms[END_REF] that the innovation norm can be used to properly tune Kalman-type filter. More precisely, best estimations are obtained for tuning parameters minimizing the mean square innovation norm. However, to the best of our knowledge, the relevance of this indicator has not been explained from theoretical arguments. Consequently, this particular point is at core of section 3. Finally, section 4 aims at illustrating the suitability of Kalman-type filters in some particular situations that frequently arise in real-world applications.

Kalman Filter for input estimation

Kalman-type filter relies on the discretized state-space representation of the considered dynamical system. The most widespread approach to discretize the state-space representation of a dynamical system is the zero-orderhold (ZOH) sampling technique, which consists in assuming that the input of the system u(t) is constant within a time step [START_REF] Lourens | An augmented kalman filter for force identification in structural dynamics[END_REF]. Formally, the corresponding discretized state-space representation, composed of the state and output equations, is written:

     x k+1 = A x k + B u k + w x k y k = O x k + D u k + v k , (1)
where x k , u k and y k are the state, input and output vectors at sample k, while A, B, O and D are, respectively, the discretized state, input, output and feedthrough matrices. Here, w x k denotes the Gaussian process noise with zero mean and covariance matrix Q x and v k is the Gaussian measurement noise with zero mean and covariance matrix R.

Although widely used, this model is conditionally stable and its performances strongly depends on the sampling rate. That is why, Liu et al. proposed a state-space model established from the Newmark integration scheme [START_REF] Liu | Explicit form of an implicit method for inverse force identification[END_REF]. However, despite Newmark approach being unconditionally stable for a well-chosen set of parameters (namely γ = 1/2 and β = 1/4) [START_REF] Géradin | Mechanical Vibrations -Theory and Application to Structural Dynamics[END_REF], it does not allow introducing numerical damping on the response of the high frequency modes present in the model. Consequently, if those modes are not properly discretized in space and time, numerical errors can affect the computed mechanical fields [START_REF] Chung | A Time Integration Algorithm for Structural Dynamics With Improved Numerical Dissipation: The Generalizedα Method[END_REF]. A practical way of introducing numerical damping in the Newmark method consists in using the HHT (Hilber-Hughes-Taylor) scheme by choosing the parameters γ and β such that [START_REF] Belytschko | Computational methods for transient analysis[END_REF]:

γ = 1 2 + α and β = 1 4 1 + α 2 with 0 ≤ α < 1. (2)
As appealing it seems, the resulting Newmark scheme is first-order accurate only. To preserve the second-order accuracy, some other numerical integration methods have been introduced. One of the most general approaches is the generalized-α method [START_REF] Chung | A Time Integration Algorithm for Structural Dynamics With Improved Numerical Dissipation: The Generalizedα Method[END_REF], which includes the Newmark method, the HHT scheme [START_REF] Belytschko | Computational methods for transient analysis[END_REF] and the WBZ (Wood-Bossak-Zienkiewicz) approach [START_REF] Wood | An alpha modification of Newmark's method[END_REF].

For a proper choice of parameters, the generalized-α method is unconditionally stable and second-order accurate and minimizes the numerical damping at low frequencies for a given high-frequency damping level. That is why, we propose to extend the approach introduced by Liu et al. in [START_REF] Liu | Explicit form of an implicit method for inverse force identification[END_REF] to derive a state-space model based on the generalized-α method.

State-space representation of the dynamical system

This section intends to provide a careful derivation of the proposed statespace representation of a dynamical system from the generalized-α method.

Practically, it consists in defining a state equation, describing the evolution of the system state over the time, and an output equation, relating the measured data to the system state and input.

State equation

Let us consider a regular partition of the time domain such that

t 0 < • • • < t k < • • • < t f (
t 0 and t f : initial and final times) and let h = t k+1 -t k denote the time step size. To derive the state equation from the generalized-α method for a mechanical system described by its mass, stiffness and damping matrices (M, K, C), one has to start from the following equations [START_REF] Erlicher | The analysis of the Generalizedα method for non-linear dynamic problems[END_REF]:

M .. d k+1-αm + C . d k+1-α f + K d k+1-α f = S u u k+1-α f , (3a)
.

d k+1 = . d k + (1 -γ) h .. d k + γ h .. d k+1 , (3b)
d k+1 = d k + h . d k + 1 2 -β h 2 .. d k + β h 2 .. d k+1 , (3c)
where γ and β are the parameters of the Newmark method, α m and α f are two averaging parameters associated to the inertia and internal/external forces, while d k = d(t k ) and u k = u(t k ) are, respectively, the displacement and excitation vectors at time t k , S u being the selection matrix of the input degrees of freedom. In the previous equation, the quantity d k+1-α must be read as:

d k+1-α = (1 -α) d k+1 + α d k . (4)
For linear systems, it is a common practice to use a modally reduced order model to limit the computational costs. The resulting modal equations are obtained by first expanding the displacement vector on the modal basis of the considered structure, meaning that:

d k = Φ n q k , (5)
where q k is the generalized displacement and Φ n is the matrix of the massnormalized mode shapes. After making use of the orthogonality relations, the previous set of equations becomes:

..

q k+1-αm + Z n . q k+1-α f + Ω 2 n q k+1-α f = Φ T n S u u k+1-α f , (6a)
.

q k+1 = . q k + (1 -γ) h .. q k + γ h .. q k+1 , (6b)
q k+1 = q k + h . q k + 1 2 -β h 2 .. q k + β h 2 .. q k+1 , (6c)
where

Ω n = diag(ω 1 , . . . , ω n ) and Z n = diag(2ξ 1 ω 1 , . . . , 2ξ n ω n ),
where ω n and ξ n are respectively the modal angular frequency and the modal damping ratio of the mode n.

The first step to derive the state equation consists in expressing acceleration and velocity vectors at instant t k+1 in terms of acceleration, velocity and displacement vectors at instant t k and the displacement vector at instant t k+1 [START_REF] Liu | Explicit form of an implicit method for inverse force identification[END_REF]. After some manipulations of Eqs. (6b) et (6c), one obtains:

.

q k+1 = γ β h q k+1 -q k + 1 - γ β . q k + 1 - γ 2β h .. q k , (7a)
..

q k+1 = 1 β h 2 q k+1 -q k - 1 β h . q k - 1 2β -1 .. q k . ( 7b
)
By introducing the last two relations into Eq. (6a), one gets:

q k+1 = L d Φ T n S u [(1 -α f ) u k+1 + α f u k ] + K d q k + C d . q k + M d .. q k , (8)
where

L d = 1 -α m β h 2 I + (1 -α f ) γ β h Z n + (1 -α f ) Ω 2 n -1
,

K d = L d 1 -α m β h 2 I + (1 -α f ) γ β h Z n -α f Ω 2 n , C d = L d 1 -α m β h I + (1 -α f ) γ β -1 Z n -α f Z n , M d = L d (1 -α m ) 1 2β -1 I + (1 -α f ) γ 2β -1 h Z n -α m I .
Replacing now Eq.( 8) in Eq. (7a), one finds:

.

q k+1 = L v Φ T n S u [(1 -α f ) u k+1 + α f u k ] + K v q k + C v . q k + M v .. q k , (10)
where

L v = γ β h L d , K v = γ β h [K d -I] , C v = γ β h C d + 1 -γ β I and M v = γ β h M d + 1 -γ 2β h I.
Similarly, after introducing Eq.(8) in Eq. (7b), one readily obtains:

..

q k+1 = L a Φ T n S u [(1 -α f ) u k+1 + α f u k ] + K a q k + C a . q k + M a .. q k , (11)
where

L a = 1 β h 2 L d , K a = 1 β h 2 [K d -I] , C a = 1 β h 2 [C d -h I] and M a = 1 β h 2 M d - 1 2β -1 I.
All things considered, if the state vector x k is defined such that

x k = [q T k . q T k ..

q

T k ] T , then Eqs. ( 8), [START_REF] Maes | Joint input-state estimation in structural dynamics[END_REF] and [START_REF] Azam | A dual kalman filter approach for state estimation via output-only acceleration measurements[END_REF] can be gathered to form the desired state equation:

x k+1 = A x k + B + u k+1 + B -u k , (12)
where

A =      K d C d M d K v C v M v K a C a M a      , B + = (1-α f )      L d L v L a      Φ T n S u , B -= α f      L d L v L a      Φ T n S u .
Finally, it should be noted that the proposed state equation based on the generalized-α method is unconditionally stable and second-order accurate provided that [START_REF] Géradin | Mechanical Vibrations -Theory and Application to Structural Dynamics[END_REF][START_REF] Chung | A Time Integration Algorithm for Structural Dynamics With Improved Numerical Dissipation: The Generalizedα Method[END_REF]:

α f = ρ ∞ 1 + ρ ∞ , α m = 3 α f -1, γ = 1 2 + α f -α m and β = 1 4 (1 + α f -α m ) 2 , (13)
where ρ ∞ is the spectral radius belonging to the interval [0, 1]. In the present paper, ρ ∞ is set to 1.

Output equation

Considering the previously defined state vector, the output equation can be expressed as:

y k = S d d k + S v . d k + S a .. d k = S d Φ n q k + S v Φ n . q k + S a Φ n .. q k = O x k , (14)
where S d , S v and S a are, respectively, the selection matrices of the displacement, velocity and acceleration data measured on the structure and

O = [S d , S v , S a ]Φ n .

State-space model

From the state and output equations presented above, one can form the state-space representation of the considered dynamical system derived from the generalized-α method, namely:

     x k+1 = A x k + B + u k+1 + B -u k + w x k y k = O x k + v k , (15)
where the process noise w x k and the measurement noise v k has been added to include uncertainties in the modeling and the measurement processes.

In the context of force reconstruction, this formulation of the dynamic problem is non-standard, because the calculation of the state vector at instant k + 1 requires the knowledge of the input vector at times k and k + 1.

It is however possible to reduce to classical state space representation by introducing the reduced state x k+1 defined such that [START_REF] Nordström | A dynamic programming algorithm for input estimation on linear time-variant systems[END_REF]:

x k+1 = A x k + B -u k + w x k . ( 16
)
As a result, the state equation writes at time k:

x k = x k + B + u k . ( 17
)
Introducing now the previous relation in Eqs. ( 16) and ( 14), one finally gets:

     x k+1 = A x k + B u k + w x k y k = O x k + D u k + v k , (18)
where

B = A B + + B -and D = O B + .
In the next parts of this paper, the previous state-space model will be referred to as G-α model.

Augmented Kalman Filter

In the present paper, an Augmented Kalman Filter (AKF) [START_REF] Lourens | An augmented kalman filter for force identification in structural dynamics[END_REF] is implemented. This approach has been chosen because it allows performing joint input-state estimation from a standard Kalman filter scheme [START_REF] Kalman | A New Approach to Linear Filtering and Prediction Problems[END_REF].

The AKF scheme consists in including the input in the state vector, in order to build an augmented state. For that purpose, a fictitious state equation for the input must be introduced. Generally, this equation writes [START_REF] Lourens | An augmented kalman filter for force identification in structural dynamics[END_REF]:

u k+1 = u k + w u k , (19)
where w u k denotes the related Gaussian process noise with zero mean and covariance matrix Q u .

Then, if one defines the augmented state as:

x a k =   x k u k   , (20)
the state-space representation of the system becomes:

     x a k+1 = A a x a k + w a k y k = O a x a k + v k , (21)
with

A a =   A B 0 I   , O a = O D , w a k =   w x k w u k   . (22)
The previous state-space representation allows performing the joint inputstate estimation from a standard Kalman filter. Practically, the filtering algorithm is divided into two main steps:

1. Initialization at t 0

• Estimation of the initial state x a 0 and the corresponding covariance matrix P a 0 ;

• Definition of the covariance matrices Q a et R associated to process and measurement noises w a k and v k respectively;

2. Filtering at t k for k = 1, . . . , n t

• Prediction of the current state x a- k and the corresponding covariance matrix P a-

k x a- k = A a x a k-1 , P a- k = A a P a k-1 A T a + Q a . (23)
• Update of the current state x a k and the related covariance matrix P a k from measurements

K k = P a- k O T a O a P a- k O T a + R -1 , x a k = x a- k + K k y k -O a x a- k , P a k = (I -K k O a ) P a- k . (24)
The implementation of a Kalman filter requires some comments, essentially related to the initialization step. First, the estimation of the initial state x a 0 requires the knowledge of the initial reduced state and input. Practically, the initial reduced state is generally unknown. That is why, one assumes that x a 0 = 0, implying that the structure is initially at rest, i.e. x 0 = 0 et u 0 = 0. Second, the covariance matrix associated to the initial state P a 0 and reflecting our confidence in the estimation of the initial state must be defined. In absence of any prior knowledge, this matrix is chosen such that:

P a 0 = σ 2 0 I. ( 25
)
The success of the filtering process also depends on the definition of the covariance matrices Q a and R related to the process and measurement noises.

These matrices reflect our confidence in the modeling and measurement processes respectively. Generally, process and measurement noises are supposed uncorrelated and white. Consequently, it is supposed that the covariance matrix R writes:

R = σ 2 r I, (26)
where σ 2 r can be seen as the variance of the measurement noise. The components of Q a are more difficult to assess a priori. For this reason, the covariance matrix Q a is expressed as:

Q a =   σ 2 x I 0 0 σ 2 u I   . (27)
Here, σ 2

x reflects the confidence in the evolution model associated to the state x k [see. Eqs. ( 15) and ( 18)], while σ 2 u is a tuning parameter, acting as a regularization parameter that limits the variation in the time history of the input. Its value significantly affects the quality of estimated solutions. In the present paper, it is determined from a selection procedure detailed in the next section. Regarding the choice of σ 2

x et σ 2 r , it will be shown in section 4 that they can be chosen in a large range of values, provided that σ 2 u is chosen accordingly.

Kalman filter tuning -MINOR criterion

A Kalman filter is a Bayesian filter. As a result, it is possible to interpret it in probabilistic terms. The description of the Kalman filtering algorithm shows that the current state update is directly related to the innovation

i k = y k -O a x a- k .
From the Bayesian perspective, it can be proved that the innovation have to theoretically exhibit a random sequence with zero mean and covariance matrix

S k = O a P a- k O T a + R [29, 30].
Formally, the innovation can be seen either as the prediction error of the state given the measurements or as a measure of the information brought by a new measurement. This explains why the innovation can serve as a tool for tuning a Kalman filter.

The question that arises now is related to the existence of an a priori optimal value of σ 2 u . Practically, this issue is difficult to answer, since this value, if exists, depends of the values of σ 2

x and σ 2 r . One can, however, consider the implementation of an automatic selection procedure like the Lcurve principle, for which various values of the regularization parameter are tested and only the value corresponding to the maximum curvature of the L-curve is selected. Here, the idea is to seek the optimal value σ 2 u such that:

σ 2 u = argmin σ 2 u J(σ 2 u |σ 2 x , σ 2 r ), (28)
where J(σ 2 u |σ 2 x , σ 2 r ) is some functional.

At this stage, the goal is to find a metric, based on the innovation, satisfying the previous relation. Practically, the functional J(σ 2 u |σ 2 x , σ 2 r ) is defined such that [START_REF] Azam | A dual kalman filter approach for state estimation via output-only acceleration measurements[END_REF][START_REF] Petersen | Estimation of the fullfield dynamic response of a floating bridge using kalman-type filtering algorithms[END_REF]:

J(σ 2 u |σ 2 x , σ 2 r ) = 1 n t nt k=1 y k -O a x a k -2 2 . ( 29
)
It results from what precedes, that the selection procedure consists in finding the tuning parameter minimizing the information brought by a new measurement over the identification duration. This selection procedure is referred to as MINOR (for Minimum Innovation Norm) criterion in the rest of the paper.

From a numerical point of view, the implementation of this selection procedure requires the calculation of the functional J(σ 2 u |σ 2 x , σ 2 r ) for each possible σ 2 u . To render this approach computationally effective, the calculations of J(σ 2 u |σ 2 x , σ 2 r ) are parallelized.

Application

From existing literature, the feasibility of the inversion is subordinate to the following conditions [START_REF] Lourens | An augmented kalman filter for force identification in structural dynamics[END_REF][START_REF] Naets | Stable force identification in structural dynamics using kalman filtering and dummy-measurements[END_REF][START_REF] Maes | Design of sensor networks for instantaneous inversion of modally reduced order models in structural dynamics[END_REF][START_REF] Massey | Inverses of Linear Sequential Circuits[END_REF][START_REF] Moylan | Stable Inversion of Linear Systems[END_REF][START_REF] Kailath | Linear systems[END_REF]]:

• Observability. This condition is satisfied if all the modes considered in the model are captured by at least one output sensor or contribute to the measured output. Formally, the system is observable provided that the matrix (S d + S v + S a )Φ n does not contain zero columns.

• Controllability. This condition is fulfilled if all the modes are excited by as least one input. Practically, the system is controllable if and only if the matrix S T u Φ n does not contain zero columns.

• Direct invertibility. This condition is required to estimate the system input without time delay. Regarding the state space model given by Eq. ( 18), this condition implies that the number of displacement, velocity or acceleration data, n d , must be greater than or equal to the number of unknown inputs, n u , and that the number of inputs cannot exceed the number of modes n m in the model. Mathematically, the direct invertibility condition is satisfied if and only if rank(D) = n u .

• Stability. This condition is related to behavior of the Kalman filter over the long term and the so-called drift effect. It results from the literature that if only acceleration and/or velocity measurements are used, the system is marginally stable leading to a drift observed on the estimated inputs over the long term. On the contrary, this effect disappears in case of displacement and/or strain measurement only provided that their number is larger than the number of inputs. To verify whether a system is stable, all the finite transmission zeros λ i of the system 1 must be such that |λ i | < 1. If at least one finite transmission zero is such that |λ i | = 1, then the system is marginally stable. In all other cases, the system is unstable [START_REF] Levine | Control System fundamentals[END_REF]. All this is crucial when instantaneous inversion is aimed. However, introducing a time delay as in [START_REF] Hsieh | Optimal time-delayed joint input and state estimation for systems with unknown inputs[END_REF] allows solving this issue at a higher computational cost that possibly prevent real-time applications. It should be said in passing that regularization strategies are based on time-delayed inversion, which explains why this problem is generally not raised in the dedicated literature.

• Uniqueness. This condition is related to the stability condition, since a system with at least one finite transmission zero cannot be uniquely identified. In all the subsequent applications, the uniqueness condition is never met, meaning that in any case the solution is not unique.

1 λ i is a finite transmission zero if rank     A -λ i I B O D     < n s + min(n u , n d ) (n s :
number of states defined in the state vector x k ) [START_REF] Kailath | Linear systems[END_REF].

This section aims at analyzing the applicability of the AKF in some situations arising frequently when dealing with source identification problems only. In particular, it is proposed to answer the following questions:

1. Are all the data equivalent for reconstructing excitation sources in time domain?

2. What is the influence of the considered state-space representation?

3. Is real-time reconstructions possible?

4. What happens over the long term?

5. Should measurements be collocated with excitations?

6. What happens in an under-determined configuration?

As it can be deduced from a careful reading of the feasibility conditions described at the very beginning of this section, some of the answers to these questions have been addressed by the literature or can at least be inferred.

However, this is from an algorithmic perspective only, since the immediate effect of a violation of these necessary conditions in terms of filter tunability or result quality is rarely illustrated. To this end, a numerical experiment is implemented.

Description of the test case

The studied structure is an IPN simply supported steel beam of length 3 m, cross-sectional area 1060 mm 2 and second moment of area 171 mm 4 [START_REF] Ginsberg | New approach for impact detection by finding sparse solution[END_REF].

The beam is excited by a hammer impact at position x 0 = 0.9 m, measured from the left end of the beam.

To synthesize this kind of excitation, it is supposed that the structure is impacted by a hammer equipped with a soft rubber tip, so as to fix the cutoff frequency of the excitation around 500 Hz. The analysis of real impact signals leads us to model such an excitation as:

f (t) = A t p-1 e -t/θ [(p -1)θ] p-1 e 1-p , ( 30
)
where p and θ are, respectively, the shape and scale parameters of the hammer impact excitation and A is the maximum impact amplitude.

In the present example, A = 15 N, while the parameters p and θ are adjusted with respect to the desired cutoff excitation frequency. This requirement leads us to set p to 9.7 and θ to 0.6 ms. Finally, it is supposed that a pre-trigger of 8 ms is applied to the excitation signal to simulate realworld experimental conditions [see Fig. 1].

On the other hand, it is supposed that displacement, velocity or acceleration data are available. To synthesize these experimental data, the equation of motion of the beam is first expanded on the corresponding modal basis, which is computed analytically. Then, the resulting modally reduced system, containing the first 53 analytical bending modes (resonance frequencies up to 1 kHz), is solved using the Newmark method for (γ, β) = (1/2, 1/4) and h = 0.1 ms (sampling frequency: 10 kHz). Here, the choice of the time step h has been made following the criterion defined in Appendix A. To be exhaustive, a modal damping ratio equal to 1% for all the modes has been imposed. In a last step, the obtained output fields are corrupted by an additive Gaussian white noise with a SNR equal to 25 dB. It should also be Finally, in absence of contradictory information, it is supposed that x a 0 = 0, while the variances σ 2 0 , σ 2 x et σ 2 r associated to the covariance matrices P a 0 , Q x and R are set to 10 -20 , 10 -6 and 10 -2 respectively. In the frequency domain, kinematic data like displacement, velocity or acceleration can be used indistinctly to solve source reconstruction problems.

But, is it the case for time domain applications? Actually, the dedicated literature indicates that it is not the case. However, what is exactly the influence of the data type on the reconstructed inputs?

In the next of this section, the reconstructions are performed on the first 0.1 s which is sufficient to avoid stability problems on estimated inputs.

Acceleration data

In the present reconstruction configuration, the system is observable, controllable, directly invertible and marginally stable. Consequently, all the necessary feasibility requirements of the system inversion are not fulfilled, which is consistent with existing literature. However, as shown in Fig. 3, the excitation field reconstructed, using the optimal value of σ 2 u obtained from the MINOR criterion and estimated to 1.08 × 10 27 [see Fig. 4], is in a very good agreement with the target excitation field. Furthermore, it should be noted that when using acceleration data the plot MINOR criterion exhibits a Z-shape. This actually indicates that any value of σ 2 u comprised between 10 11 and 10 40 can lead to a consistent reconstruction.

Velocity data

For velocity-based reconstruction, the system is observable, controllable, directly invertible and marginally stable. Unfortunately, contrary to the previous case, the identified excitation field is far from the target excitation field 5]. Furthermore, it should be noted that in this case, the plot of MINOR criterion has a clear global minimum, meaning that choosing, a priori, a relevant value for σ 2 u is far from an easy task [See Fig. 6].

This result is quite surprising when looking at the conclusions of the feasibility conditions. Indeed, the analysis of the reconstructed excitation clearly leads to think that the direct invertibility condition is not satisfied.

Actually, the reduced state-space model given by Eq. ( 18) respects the direct invertibility condition, which is not the case for the original state-space model given by Eq. [START_REF] Risaliti | A state-input estimation approach for force identification on an automotive suspension component[END_REF]. Furthermore, this is in fact exactly the same for the acceleration-based identification. In the latter case, however, no particular problems had been encountered and the conclusions of the feasibility study are in line with existing literature [START_REF] Naets | Stable force identification in structural dynamics using kalman filtering and dummy-measurements[END_REF][START_REF] Maes | Design of sensor networks for instantaneous inversion of modally reduced order models in structural dynamics[END_REF]. Consequently, as clearly written by Maes in [START_REF] Maes | Design of sensor networks for instantaneous inversion of modally reduced order models in structural dynamics[END_REF], the feasibility condition are not sufficient to guarantee that the quantities of interest can be identified correctly in the presence of noise. Finally, all the previous remarks lead us to wonder about the influence of the state-space model used to perform the reconstruction. This problem is discussed in section 4.3. In particular, it will be shown that the results obtained with the G-α model are actually better than those obtained with a more classical state-space representation based on the ZOH assumption.

Displacement data

For displacement-based reconstruction, the system is observable, controllable, directly invertible and marginally stable. As for velocity-based reconstruction, the results are disappointing [see Figs. 7 and8] and exactly the same observation can be made regarding the conclusions of the feasibility conditions.

Summary

The reconstructions presented previously clearly show that the nature of the data used for estimating the excitation sources has a strong influence on the obtained results. This leads us to conclude that if one has to use only one type of data, then acceleration measurements are certainly the best choice. This result is reassuring, since accelerometers are commonly used in structural dynamics for their ability to accurately measure vibrations at relatively low cost. This question deserved, however, to be illustrated in view of emerging or less standard measurement techniques, such that Digital Image Correlation for measuring displacement [START_REF] Schreier | Image correlation for shape, motion and deformation measurements: Basic concepts, theory and applications[END_REF] or multipoint laser vibrometry [START_REF] Macpherson | Multipoint laser vibrometer for modal analysis[END_REF][START_REF] Yang | A multipoint laser Doppler vibrometer with fiber-based configuration[END_REF] for measuring velocities synchronously.

As a side note, it should be mentioned here that the previous study has also been conducted for the joint input-state estimator [START_REF] Lourens | Joint input-response estimation for structural systems based on reduced-order models and vibration data from a limited number of sensors[END_REF] and the Dual Kalman Filter (DKF) [START_REF] Azam | A dual kalman filter approach for state estimation via output-only acceleration measurements[END_REF]. Unfortunately, some problems have been encoun-tered during their application. For the joint input-state estimator, singularity of the filtering matrices arose when applying it to displacement and velocity data. As for DKF, an extreme sensivity to the value of the variances σ 2 0 and σ 2

x has been observed, making the tuning of the Kalman Filter highly difficult. Consequently, further studies need to be carried out to precisely determine to what extent the conclusions drawn in the present paper are applicable to a larger class Kalman-type filters.

What is the influence of the considered state-space representation?

The previous application has clearly illustrated the influence of the data type on the reconstructed solutions. However, it has also highlighted some potential contradictions between the conclusions of the feasibility conditions and the effective results. As underlined above, this leads us to wonder about the influence of the state-space model used to perform the reconstruction, since the feasibility conditions are based on it.

To properly study this problem, all the reconstructions presented in the previous section are performed using an Augmented Kalman Filter built from a state-space model based on the ZOH assumption [see Eq. ( 1)]. Formally, this model, classically used in structural dynamics and referred to as ZOH model in the rest of this paper, is given by:

     x k+1 = A x k + B u k y k = O x k + D u k , (31)
where

A = e Ach , B = [A -I]A -1 c B c , D = S a Φ n Φ T n S u , O = [S d Φ n -S a Φ n Ω n , S v Φ n -S a Φ n Z n ] , A c =   0 I -Ω 2 n -Z n   , B c =   0 M -1 n Φ T n S u   . (32)
For this state-space model, the system is:

• observable, controllable, directly invertible and marginally stable for acceleration-based reconstruction;

• observable, controllable, marginally stable, but not directly invertible for velocity-based reconstruction;

• observable, controllable, stable, but not directly invertible for displacementbased reconstruction.

As it can be seen, the feasibility conditions applied to the ZOH model give rather different conclusions. However, the quality of the estimated inputs resulting from the ZOH model from acceleration or displacement data is globally similar to that obtained with the G-α model [see Figs. 9, 10, 14 and 15]. It should even be noted that for acceleration-based reconstruction, the G-α model provides a better estimation. On the contrary, for velocitybased estimation, the results are significantly different. Indeed, the input estimated from the ZOH model is extremely noisy, which is not the case with the G-α model that leads to a very smooth solution [see Fig. 11]. This result seems inconsistent given the acceleration and displacement counterparts. The explanation of this surprising result can be found in the shape of the plot of the MINOR criterion [see Fig. 12]. Indeed, the analysis of this figure shows that the functional J(σ 2 u |σ 2 x , σ 2 r ) presents a clear local minimum. When choosing the value of σ 2 u corresponding to this local minimum, here σ 2 u = 1.06 × 10 5 , the reconstruction obtained from ZOH model becomes consistent with that obtains with the G-α model [see Fig. 13]. Consequently, this demonstrates that the G-α model allows an easier and more robust tuning of the AKF than the ZOH model. In this case, however, the system is marginally stable, meaning that on the long term the estimated quantities will diverge [START_REF] Naets | Stable force identification in structural dynamics using kalman filtering and dummy-measurements[END_REF].

This question is addressed in section 4.5.

Finally, another interesting consideration is the influence of the statespace representation with respect to the the time step size. To this end, data generation and input reconstruction have been performed for several time step sizes (i.e. for different sampling frequencies), corresponding to one twentieth (h = 0.1 ms), one tenth (h = 0.2 ms), one fifth (h = 0.4 ms) and one half (h = 10 ms) of the oscillation period of the last mode excited by the hammer impact (around 500 Hz in the present application). To quantify the accuracy of the reconstructed solution from acceleration data with respect to the state-space representation and the time step size, the relative error is evaluated. Formally, the relative error is defined such that:

E = u ref -u id 2 2 u ref 2 2 , ( 33
)
where u ref is the reference excitation field, while u id is the identified excitation field.

The results presented in Table 1 strengthen the motivation behind the development of the G-α model, since, in any case, reconstructions based on the G-α model are far more accurate than those based on the ZOH model. In the light of these results, input reconstruction presented in the next sections will be based on the state-space model introduced in section 2 and acceleration measurements as output data.

Is real-time reconstructions possible?

Possible real-time estimation is one of the most appealing features of the AKF (and Kalman-type filtering in general). However, the previous sections demonstrate that the quality of the reconstructed excitation field is strongly dependent on the value of the variance σ 2 u and cannot consequently be chosen arbitrarily. This is particularly true for displacement and velocity-based reconstructions from the G-α model, for which the plot of MINOR criterion exhibits a clear global minimum. Unfortunately, the value σ 2 u is itself related to the values of σ 2

x and σ 2 r . In this situation, it is extremely difficult or even impossible to infer a priori a value of σ 2 u . As a result, when displacement or velocity data are used, consistent real-time reconstructions seem hard to achieve. For acceleration-based reconstructions, the situation is rather different. Indeed, as shown in Fig. 4, the plot of MINOR criterion has a Z-shape, meaning that for a given couple (σ 2 x , σ 2 r ), σ 2 u can be chosen in a wide range of values to provide similar results. This observation has been stated in passing by Lourens et al. in [START_REF] Lourens | An augmented kalman filter for force identification in structural dynamics[END_REF], but no further indications or illustrations have been given to support this statement. In the present case, the problem is to determine the relation that must satisfy the variances σ 2

x , σ 2 r and σ 2 u to expect almost insensitive reconstructions.

To determine this relation, the optimal triplets ( σ 2

x , σ 2 r , σ 2 u ) have been sought using the MINOR criterion. More precisely, they have been sought such that:

( σ 2 x , σ 2 r , σ 2 u ) = argmin (σ 2 x ,σ 2 r ,σ 2 u ) J(σ 2 x , σ 2 r , σ 2 u ), (34)
where the functional J(σ 2 x , σ 2 r , σ 2 u ) is given by:

J(σ 2 x , σ 2 r , σ 2 u ) = 1 n t nt k=1 y k -O a x a k -2 2 . (35)
For the acceleration data used previously, one gets σ 2 x = 1.83 × 10 5 , σ 2 r = 6.95 × 10 36 and σ 2 u = 6.95 × 10 36 . With this optimal set, the reconstructed excitation field agrees very well with the target one [see Fig. 16].

Combining this result with that obtained in section 4.2.1, one can conclude that the estimation of the excitation field is almost insensitive to the choice of the variance parameters provided that: In other words, if the previous condition is fulfilled, real-time reconstructions can be achieved. To examplify this condition, let us take, for instance,

σ 2 x σ 2 r ≤ σ 2 u . (36)
σ 2
x = 10 -20 , σ 2 r = 10 15 and σ 2 u = 10 20 . As expected, Fig. 17 shows that when the relation [START_REF] Ginsberg | New approach for impact detection by finding sparse solution[END_REF] is fulfilled, then the estimation of the excitation field is very satisfying.

Actually, the condition given by Eq. ( 36) is quite logical, since σ 2 r must be less than or equal to σ 2 u to avoid a too strong filtering of the estimated excitation field. Indeed, having σ 2 u < σ 2 r implies that one has a greater confidence in Eq. ( 19) than in the output equation. In other words, the filter will try to respect the relation u k+1 ≈ u k as best as possible. On the other hand, σ 2

x must necessarily be less than σ 2 r , because it is well-known that inverse problems are very sensitive to noise. Consequently, having a greater that a reliable model of the structure can be obtained, which is the case in the present application. In other words, if the level of modeling error is high, which arises in many situations, Eq. ( 36) may not hold. It should however be noticed that in such a case, the probability of good reconstruction is weak.

Finally, for the sake of completeness, displacement and velocity-based reconstructions obtained from the optimal triplet ( σ 2

x , σ 2 r , σ 2 u ) determined from the MINOR criterion are presented in Appendix B for the G-α and ZOH models. In addition, this appendix provides further insights into the influence of the optimal variance parameters selected from the MINOR criterion and the difficulty of properly tuning a Kalman filter when velocity-based reconstruction are to be carried out.

What happens over the long term?

In the literature, it is well-known that the inputs estimated from acceleration data suffer from the so-called drift effect related to the marginal stability of the system [START_REF] Naets | Stable force identification in structural dynamics using kalman filtering and dummy-measurements[END_REF]. This divergence issue, that occurs in such a situation, is illustrated in Fig. 18 presenting the excitation field reconstructed over 0.5 s. It should be mentioned that σ 2 u has been arbitrarily set to 10 20 , which respects the condition given by Eq. ( 36). This undesirable effect can be avoided by complementing acceleration data with collocated displacement data. To prevent collocated displacement measurements, Naets et al. [START_REF] Naets | Stable force identification in structural dynamics using kalman filtering and dummy-measurements[END_REF] proposed to use fictitious displacement measurements collated with the real measurement positions. This is done by introducing a fictitious output equation to the original state-space model [see Eq. ( 14)], defined such that:

0 = O f x a k + v f k , (37)
where O f = [S a Φ n , 0, 0] and v f k is the Gaussian fictitious measurement noise with zero mean and covariance matrix R f . The addition of the previous relation to the G-α model allows obtaining the following state-space model:

     x a k+1 = A a x a k + w a k y a k = O f a x a k + v a k , (38)
with

y a k =   y k 0   , O f a =   O a O f   , v a k =   v k v f k   , (39)
where v a k is a Gaussian white noise with zero mean and covariance matrix:

R =   R 0 0 R f   avec R f = σ 2 f I, (40)
where σ 2 f is the variance parameter associated to fictitious measurements.

One of the key point of this strategy lies in the choice of σ 2 f . In the literature, it is suggested to base our choice on the order of magnitude of the deformation of the structure [START_REF] Naets | Stable force identification in structural dynamics using kalman filtering and dummy-measurements[END_REF], but this possibly requires numerical calculations. To prevent this, the MINOR criterion is used once again. Consequently, the optimal value of σ 2

f is sought such that:

σ 2 f = argmin σ 2 f J(σ 2 f |σ 2 x , σ 2 r , σ 2 u ), (41)
where the functional J(σ 2 f |σ 2 x , σ 2 r , σ 2 u ) is defined as in Eq. [START_REF] Kalman | A New Approach to Linear Filtering and Prediction Problems[END_REF].

At this stage, it is important to note that to limit the influence of the drift effect using dummy measurements, the condition provided by Eq. ( 36) must be restricted. Indeed, it has been found that it is necessary to impose σ 2 r of the order of magnitude of σ 2 u to obtain consistent estimations.

Consequently, by setting σ 2 x = 10 -6 and σ 2 r = σ 2 u = 10 20 , the application of the MINOR criterion allows finding σ 2 f = 3.78 × 10 15 . In this situation, the system is observable, controllable, directly invertible and marginally stable.

The resulting estimated excitation field is presented in Fig. 19. The analysis of this figure shows, as indicated by Naets et al. [START_REF] Naets | Stable force identification in structural dynamics using kalman filtering and dummy-measurements[END_REF], that complementing acceleration data by a set of collocated dummy measurements allows limiting the drift effect over the long term.

One can, however, wonder if this procedure (addition of dummy measurements + selection of σ 2 f from MINOR criterion) doesn ot definitively prevent any chance of real-time reconstruction. Actually, everything depends on the shape of the plot of the MINOR criterion. If its shape is similar to that presented in Fig. 4, then a real-time reconstruction over the long term is possible. The answer to this question is provided in Fig. 20. This figure shows that the plot of the MINOR criterion has the desired shape, which make possible the application of this strategy for real-time estimations.

However, as for the other variance parameters, the value of σ 2 f cannot 36). More specifically, if has been found from a set of numerical experiments that:

σ 2 x σ 2 f < σ 2 r ≤ σ 2 u , (42)
provided that a reliable model of the structure (i.e. with a low modeling error level) could be obtained.

The fact that σ 2 f must be less than σ 2 r is quite logical, because if its not the case, this means that we trust in the measurements more than in the dummy output equation which aims at stabilizing the inversion. As a result, if this condition is not satisfied, there is no reason for the problem to be stabilized. As an example of Eq. ( 42), let us consider σ 2

x = 10 -6 , σ 2 r = 10 19 , σ 2 u = 10 20 and σ 2 f = 10 14 . For this set of parameters, the corresponding estimated excitation field is presented in Fig. 21.

Finally, for the sake of completeness, reconstructions obtained using real displacement data are presented in Appendix C. This appendix shows surprisingly that a better stabilization of the reconstructed excitation field is obtained using dummy measurements.

Should measurements be collocated with excitations?

This question is very important, since it is not always possible to measure output data at excitation points. Actually, in the vast majority of the articles published in the literature, at least one output data is collocated with or is in the close vicinity of the actual excitations. However, as observed by Lourens et al. in [START_REF] Lourens | An augmented kalman filter for force identification in structural dynamics[END_REF], it seems that the violation of this implicit requirement leads to unsatisfying reconstructions. To illustrate this, let us consider In this configuration, the system is observable, controllable, directly invertible and marginally stable when using acceleration data only. When performing the reconstruction on the first 0.1 s, a time delay is observed on the reconstructed excitation field [see Fig. 23]. Furthermore, this non collocated configurations does not allow real-time reconstructions, since the plot of the MINOR criterion has not a characteristic Z-shape [see Fig. 24].

Even with an optimal triplet (σ 2

x , σ 2 r , σ 2 u ) selected from the MINOR criterion, the situation remains the same, although the estimated excitation field is less noisy as illustrated in Fig. 25. ). However, it is clear that the absolute value of the condition number of the observation matrix indicates that the system is extremely ill-conditioned, which is sufficient to infer that a perturbation in the reconstruction parameters can strongly affect the quality of the estimated quantities, which is the case here, since the sensor configuration is not optimal for an inverse problem. As a consequence, it is strongly recommended to work with a fully collated reconstruction configuration for maximizing the chance of success. This conclusion is applicable to other Kalman-type filters such as the joint-input estimator. If it is not possible, it is recommended to adopt a Kalman smoothing approach as the one proposed recently by Maes et al. [START_REF] Maes | A smoothing algorithm for joint input-state estimation in structural dynamics[END_REF].

What happens in an under-determined configuration?

This problem arises when one suspects that excitation sources can act in unmeasured regions of the structure. Such a situation is known to be difficult to deal with, because the problem cannot have a unique solution or even no solution, if it is inconsistent. To analyze the behavior of the Kalman filter in such a situation, the reconstruction is performed in an a priori favorable configuration, since one of the accelerometers is collocated with the actual excitation [see Fig. 26]. In this situation, the system is observable, controllable, marginally stable, but not directly invertible. Consequently, a poor reconstruction is expected. For the sake of completeness is the MINOR criterion is presented in Fig. 28.

A practical way to deal with under-determined problems is to add constraints to the problem. Unfortunately, Kalman-type filters don't assume any spatial distribution on the sources a priori or any a particular form of the excitation signals. Consequently, in under-determined configuration, it is hard to expect reliable reconstructions. In particular, it has been shown that the G-α model provides reconstructions more robust than those obtained from the ZOH model, classically used for solving estimation problems. Then, specific conditions have been derived

to make real-time reconstructions from acceleration data only possible, even over the long term thanks to the use of dummy displacement measurements.

It has finally been shown that, for improving the chance of good estimation, it is necessary to solve the problem from collocated configurations. Finally, this study has highlighted the fact that the feasibility conditions proposed in the literature have to be used as an indicator, since the violations of some of these conditions doesn't necessarily imply the failure of the estimation process in operating conditions (presence of noisy data for instance). Of course, many issues remain to be studied, such as the influence of modeling errors or the reconstruction of random excitation or moving loads, but these questions deserve specific studies. Finally, it should be noted that, at this stage, the conclusions drawn in this paper are mainly applicable to the AKF. The generalization of conclusions obtained on the basis of the AKF to other Kalman filter based input reconstruction techniques must be studied in details.

Figure 1 :

 1 Figure 1: Synthesized hammer impact excitation signal -(a) Time domain representation and (b) Frequency domain representation

Figure 2 :

 2 Figure 2: Locations of the measurement and reconstruction points all along the beam -( ) Ends of the beam, ( ) Measurement points and (•) Reconstruction points

Figure 3 :

 3 Figure 3: Reconstruction of the excitation field from acceleration data for σ 2 u = 1.08×10 27 -(a) Waterfall representation and (b) Identified time signal at excitation point -(-) Reference signal and (--) Reconstructed signal

Figure 4 :

 4 Figure 4: Plot of the functional J(σ 2 u |σ 2 x , σ 2 r ) for acceleration data -(-) J(σ 2 u |σ 2 x , σ 2 r ) and (•) Location of the minimum of the functional

Figure 5 :Figure 6 :

 56 Figure 5: Reconstruction of the excitation field from velocity data for σ 2 u = 5.29×10 4 -(a) Waterfall representation and (b) Identified time signal at excitation point -(-) Reference signal and (--) Reconstructed signal

Figure 7 :Figure 8 :

 78 Figure 7: Reconstruction of the excitation field from displacement data for σ 2 u = 1.70×10 6 -(a) Waterfall representation and (b) Identified time signal at excitation point -(-) Reference signal and (--) Reconstructed signal

Figure 9 :Figure 10 :Figure 11 :Figure 12 :Figure 13 :Figure 14 :Figure 15 :

 9101112131415 Figure 9: Reconstruction of the excitation field from acceleration data for σ 2 u = 1.70×10 33 with the ZOH model -(a) Waterfall representation and (b) Identified time signal at excitation point -(-) Reference signal, (--) Reconstructed signal from the G-α model and (-• -) Reconstructed signal from the ZOH model

Figure 16 :

 16 Figure 16: Reconstruction of the excitation field from acceleration data for σ 2 x = 1.83×10 5 , σ 2 r = 6.95 × 10 36 and σ 2 u = 6.95 × 10 36 -(a) Waterfall representation and (b) Identified time signal at excitation point -(-) Reference signal, (--) Reconstructed signal

Figure 17 :

 17 Figure 17: Reconstruction of the excitation field from acceleration data for σ 2 x = 10 -20 , σ 2 r = 10 15 and σ 2 u = 10 20 -(a) Waterfall representation and (b) Identified time signal at excitation point -(-) Reference signal, (--) Reconstructed signal

Figure 18 :

 18 Figure 18: Reconstruction of the excitation field from acceleration data over 0.5 s for σ 2 u = 10 20 -(a) Waterfall representation and (b) Identified time signal at excitation point -(-) Reference signal and (--) Reconstructed signal

Figure 19 :Figure 20 :

 1920 Figure 19: Reconstruction of the excitation field from acceleration data and dummy measurements over 0.5 s forσ 2 x = 10 -6 , σ 2 r = σ 2 u = 10 20 and σ 2 f = 3.78 × 10 15 -(a) Waterfall representation and (b) Identified time signal at excitation point -(-) Reference signal and (--) Reconstructed signal

Figure 21 :

 21 Figure 21: Reconstruction of the excitation field from acceleration data and dummy measurements over 0.5 s for σ 2 x = 10 -6 , σ 2 r = 10 19 , σ 2 u = 10 20 and σ 2 f = 10 14 -(a) Waterfall representation and (b) Identified time signal at excitation point -(-) Reference signal and (--) Reconstructed signal

Figure 22 :(Figure 23 :Figure 24 :Figure 25 :

 22232425 Figure 22: Locations of the measurement and reconstruction points all along the beam in not fully collocated configuration -( ) Ends of the beam, ( ) Measurement points and (•) Reconstruction points

Figure 26 :

 26 Figure 26: Locations of the measurement and reconstruction points all along the beam in under-determined collocated configuration -( ) Ends of the beam, ( ) Measurement points and (•) Reconstruction points

Figure 27 :Figure 28 :

 2728 Figure 27: Reconstruction of the excitation field from acceleration data in underdetermined configuration over 0.1 s for σ 2 x = 10 -6 , σ 2 r = 10 -2 and σ 2 u = 1.55 × 10 38 -(a) Waterfall representation and (b) Identified time signal at excitation point -(-) Reference signal and (--) Reconstructed signal

Figure B. 30 :Figure B. 31 :

 3031 Figure B.30: Reconstruction of the excitation field from displacement data from the G-α model and σ 2 x = 10 -20 , σ 2 r = 1.83 × 10 5 and σ 2 u = 5.45 × 10 14 -(a) Waterfall representation and (b) Identified time signal at excitation point -(-) Reference signal, (--) Reconstructed signal

Table 1 :

 1 Relative error (%) obtained from acceleration data for different time step sizes

	with respect to the state-space representation
	h (ms) G-α model ZOH model
	0.1	1.78	3.77
	0.2	4.05	7.90
	0.4	6.56	17.69
	10	9.72	53.67

Acknowledgements

The work presented in this paper has been performed within the research collaboration framework agreement CT-413-1601-CNES "Experimental reconstruction of transient mechanical sources" with the Centre National d'Étude Spatiales (Cnes). The authors are grateful to the Cnes for funding this research.

Appendix A. Estimation of the time step

This appendix aims at providing a simple criterion for properly choosing the time step for G-α and ZOH models and more generally for discretizing mechanical problems in time domain. To this end, let us consider an undamped mechanical system described by its mass and stiffness matrices (M, K). When the structure is linear and time-invariant, the equation of motion at any time t k is classically written:

This equation can be discretized using a finite difference scheme. In doing so, the previous relation becomes:

To assess the discretization error, the Fourier transform of Eqs. (A.1) and (A.2) is computed. In doing so, one obtains:

with:

Consequently, to impose the equality of relations (A.3a) and (A.3b), one necessarily has:

The only possible solution is ωh = 0. Nevertheless, to determine a reasonable value of the step size, one can weaken condition (A.5) by satisfying this condition within some tolerance . Mathematically, the weak condition is expressed as:

By computing the Taylor expansion of the cosine function up to the fifth order, one infers:

Practically, one can set ω = ω max , where ω max is the cut-off frequency of the excitation. For the application presented in this paper, the cut-off frequency is equal to 500 Hz. It results that if the approximation error is set to = 1%, then h ≤ 0.11 ms. Here, the time step has been set to h = 0.1 ms, which approximately corresponds to one twentieth of the oscillation period of the last mode excited by the hammer impact.

Appendix B. Displacement and velocity based reconstructions from

This appendix presents the excitation fields reconstructed from the opti- As in section 4.3, the results obtained from the G-α model and the ZOH model are very similar when displacement data are used. But probably more interesting, the same observation can be made for velocity data, which was not the case in section 4.3. Consequently, when combining this result whith those obtained in sections 4.2 and 4.3, it is clear that the tuning of the AKF is quite challenging when only velocity data are available. In addition, the latter result clearly shown that the G-α model is more robust than the ZOH model for solving purely input estimation problems.

Appendix C. Real vs. dummy displacement measurements

In this appendix the reconstruction is performed from real displacement data instead of dummy measurements. Consequently, the output matrix O is expressed as:

For the G-α model, the system is observable, controllable, directly invert- In both cases, it can be seen that even if a significant stabilization of the results is observed, the drift effect is still present. The comparison with the reconstruction from dummy measurements clearly emphasizes the superiority of the latter approach over a stabilization strategy based on real displacement measurements. It should however be noted that when the actual input has a trend/static component, then only the use of real displacement data allows retrieving it [START_REF] Naets | Stable force identification in structural dynamics using kalman filtering and dummy-measurements[END_REF][START_REF] Maes | Design of sensor networks for instantaneous inversion of modally reduced order models in structural dynamics[END_REF].


M Aucejo

O De Smet

A STATIONARY ITERATED MULTIPLICATIVE REGULARIZATION FOR FORCE RECONSTRUCTION PROBLEMS

Iterated Tikhonov regularization can be seen as an iterative refinement of an initial solution obtained by applying the standard Tikhonov regularization. This kind of approach is particularly suited when noisy data with low signal-to-noise ratio (SNR) are to be used. However, to apply such an approach, an initial regularization parameter has to be estimated. Generally, it is obtained numerically through automatic selection procedures that can require intensive computations. In this respect, the multiplicative regularization is an interesting alternative, since it does not require any selection procedure for choosing the regularization parameter. In return, the amount of regularization is automatically adjusted throughout an iterative resolution process. To combine the best of both worlds for solving force reconstruction problems, we propose a stationary iterated multiplicative regularization that allows using prior information on the type of sources acting on the structure. The validity and the suitability of the proposed approach are illustrated numerically on a free-free beam excited by a point force.

Introduction

Iterated Tikhonov regularization has been extensively studied by the inverse problems community [START_REF] King | Approximation of generalized inverses by iterated regularization[END_REF][START_REF] Engl | On the choice of the regularization parameter for Iterated Tikhonov regularization of ill-posed problems[END_REF][START_REF] Hanke | Nonstationary Iterated Tikhonov Regularization[END_REF][START_REF] Huang | Projected nonstationary iterated tikhonov regularization[END_REF][START_REF] Buccini | Iterated Tikhonov regularization with a general penalty term[END_REF]. It can practically be viewed as an iterative refinement of an initial solution obtained by applying the standard Tikhonov regularization. Surprisingly, this regularization strategy seems to be seldom applied in mechanics [START_REF] Cimitière | Solution of the Cauchy problem using iterated Tikhonov regularization[END_REF][START_REF] Jang | Numerical expremients on determination of spatially concentrated timevarying loads on a beam: an iterative regularization method[END_REF] despite its ability to compute improved identified solutions. In the context of force reconstruction problems, iterated Tikhonov regularization has some drawbacks. First, it requires the estimation of a regularization parameter, that can be numerically computed through 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] or the Generalized Cross-Validation [START_REF] Golub | Generalized cross-validation as a method for choosing a good ridge parameter[END_REF]. However, when the reconstruction problem is large, this calculation can be very expensive. Second, the regularization term promotes distributed solutions, which is not desirable when the structure is actually excited by localized sources. In this respect, the multiplicative regularization [START_REF] Aucejo | A multiplicative regularization for force reconstruction[END_REF] is an interesting alternative, since it eliminates the need for the selection of the regularization parameter while taking into account one's prior knowledge on the sources to identify. However, when this strategy is applied to very noisy data, obtained reconstructions are generally disappointing. In such a situation, an iterated approach is particularly well suited. That is why, to combine the best of both worlds, a stationary iterated multiplicative regularization is proposed. In this contribution, the main features of the proposed regularization strategy are set out, while its validity is assessed through a numerical experiment performed on a free-free beam.

The need for another regularization strategy

Classically, the stationary iterated Tikhonov regularization is defined as follows:

   F (0) = 0 F (k) = argmin F X -HF 2 2 + λ F -F (k-1) 2 2 for k ≥ 1 , (1)
where in the context of force reconstruction:

• X is the measured vibration field;

• H is the transfer functions matrix of the structure describing its dynamic behavior;

• F (k)
is the force vector identified at iteration k;

• λ is the regularization parameter.

When applying this formulation to the identification of a point force exciting a free-free steel beam at 350 Hz from a measured vibration field having an SNR equal to 30 dB [see section 4 for a detailed description of the test case], it is clear that the reconstructions resulting from the application of the standard Tikhonov regularization or its stationary iterated version are too smooth compared to the target excitation fields [see Figure 1]. As explained in the introduction, the previous formulation has two drawbacks, since it does not exploit one's prior knowledge of the sources to identify and it requires the calculation of a regularization parameter, which can be computationally intensive depending on the size of the reconstruction problem. To bypass these potential difficulties, a multiplicative regularization has been introduced [START_REF] Aucejo | A multiplicative regularization for force reconstruction[END_REF][START_REF] Aucejo | Further investigation on "a multiplicative regularization for force reconstruction[END_REF]. Formally, it consists in solving the following minimization problem:

F = argmin F\{0} X -HF 2 2 • F q q , ( 2
)
where q is a tuning parameter included in the interval ]0, +∞[, • q is the q -norm or quasi-norm. Practically, q is set to 2 if the solution is supposed to be distributed, while q can be chosen equal to or less than 1 if the solution is supposed to be localized [START_REF] Aucejo | Bayesian source identification using local priors[END_REF].

The application of the multiplicative regularization from the measured vibration field used previously gives, when q = 0.5, the reconstruction field presented in Figure 2. This figure shows that the reference excitation field and the reconstructed one are very close together when the measured data are little noisy (high SNR). Unfortunately, when the data are very noisy (low SNR), a degradation of the regularized solution accuracy is observed [see Figure 2].

In the light of these observations, it seems interesting to develop an iterated multiplicative regularization, combining the benefits of the multiplicative regularization and the iterated Tikhonov regularization. A first approach to this problem consists in developing a stationary iterated multiplicative regularization, which allows obtaining a refined regularized solution, while preserving the computational efficiency of the regularization strategy.

Stationary iterated multiplicative regularization

This section aims at introducing the proposed stationary iterated multiplicative regularization. To render the presentation more didactic, it is focused on the main concepts and features of the method.

Basic principle

The proposed stationary iterated multiplicative regularization is actually a simplified version of the following non-stationary iterated minimization problem: Practically, the resolution of the minimization problem at iteration k ≥ 1 requires the implementation of an iterative procedure. In the present contribution, an Iteratively Reweighted Least Squares (IRLS) algorithm has been implemented. Basically, it consists in iteratively computing the solution of the problem by recasting the q -norm into a weighted 2-norm. As a result, at iteration k of the proposed iterated regularization, the estimated force vector F (k,j) at iteration j of the IRLS algorithm is explicitly given by [START_REF] Aucejo | A multiplicative regularization for force reconstruction[END_REF]:

     F (0) = 0 F (k) = argmin F\{ F (k-1) } X -HF 2 2 • F -F (k-1) q q for k ≥ 1 . (3)
F (k,j) = H H H + α (k,j) W (k,j-1) -1 H H X + α (k,j) W (k,j-1) F (k-1) , (4)
where:

• W (k,j-1) is a diagonal weighting matrix depending on the tuning parameter q and the solution at iteration j -1, namely F (k,j-1)

[10, 12, 13, 14];

• α (k,j) is the adaptive regularization parameter 1 , defined such that:

α (k,j) := X -H F (k,j-1) 2 2 W (k,j-1) 1/2 F (k,j-1) -F (k-1) 2 2 . (5)
When the IRLS algorithm has converged after Nj iterations, the solution vector F

, the weighting matrix W (k) and the adaptive regularization parameter α (k) are returned by the algorithm and are defined such that:

F (k) := F (k,N j ) , W (k) := W (k,N j ) and α (k) := α (k,N j ) . (6)
Consequently, the previous non-stationary iterated multiplicative algorithm is composed of a main (outer) iteration corresponding to the refinement stage and a nested (inner) iterative procedure to compute F

. From there, the stationary iterated multiplicative algorithm is derived by fixing once for all the values of the weighting matrix W (k) and the adaptive regularization parameter α (k) to their values obtained after the first main iteration (i.e. k = 1) of the non-stationary algorithm, that is:

W = W (k) := W (1,N j ) and α = α (k) := α (1,N j ) . ( 7
)
In other words, the stationary multiplicative regularization can be expressed as:

             F (0) = 0 F (1)
, W, α = argmin

F\{ F (0) } X -HF 2 2 • F -F (0) q q for k = 1 F (k) = argmin F X -HF 2 2 + α W 1/2 F -F (k-1) 2 2 for k ≥ 2 . ( 8
)
Finally, it is important to note that from k = 2, the explicit expression of the force vector

F (k)
is:

F (k) = H H H + α W -1 H H X + α W F (k-1)
.

Consequently, only the matrix A = H H H + α W has to be inverted to compute

F (k)
for any k ≥ 2. This allows deriving a computationally efficient algorithm, since if A is factorized using a LU decomposition for instance, the overall computation of

F (k)
for k ≥ 2 is almost as expensive as than the computation of F (1) .

Stopping criterion

As explained in the previous section, the proposed stationary iterated multiplicative regularization first requires the resolution of a multiplicative regularization at iteration k = 1 from the IRLS procedure described in Ref. [START_REF] Aucejo | A multiplicative regularization for force reconstruction[END_REF], defining thoroughly the initial solution and the stopping criterion. Here, the question is when to stop the main iteration loop. In this contribution, the stopping criterion is related to the relative variation of the functional:

J( F (k) ) = X -H F (k) 2 2 + α W 1/2 F (k) -F (k-1) 2 2 (10)
between two successive iterations. Practically, the algorithm is automatically stopped when a prescribed tolerance defined by the user is reached. Here, the tolerance is set to 10 -2 .

Numerical validation

In the present numerical validation, the studied structure is a free-free steel beam with dimensions 1×0.03×0.01 m 3 excited by a unit point force at 350 Hz. The coordinate of the point force, measured from the left end of the beam, is x0 = 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 also worth mentioning that an additive Gaussian white noise with a prescribed SNR has been added to the simulated data in order to synthesize the measured vibration field. Finally, the FE model of the structure is 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.

To assess the validity of the proposed iterated multiplicative regularization strategy, it is compared to the ordinary multiplicative regularization. For this purpose, the choice of the tuning parameter q is crucial. The analysis of the numerical test case shows that the beam is only excited by a point force. In this context, it is reasonable to set q = 0.5 [START_REF] Aucejo | A multiplicative regularization for force reconstruction[END_REF]. The resulting reconstructions proposed in Figure 3 point out that for measured data with a high SNR, both regularization strategy provide, as expected, similar results. Here, the overall computation time is equal to 8 ms (14 iterations) for the ordinary multiplicative regularization and 12 ms (14 inner iterations for k = 1 and 4 outer iterations) for the stationary iterated multiplicative regularization. But more interestingly, for measured data with a low SNR, the iterated multiplicative regularization allows obtaining a far more better estimation of the point force amplitude than the ordinary multiplicative regularization, while preserving the quality of the reconstructed excitation field. In this situation, the overall computation time is equal to 34 ms (132 iterations) for the ordinary multiplicative regularization and 43 ms (132 inner iterations for k = 1 and 5 outer iterations) for the stationary iterated multiplicative regularization.

Conclusion

In the present contribution, a stationary iterated multiplicative regularization has been presented for identifying mechanical sources acting on a structure. The underlying idea behind this regularization strategy was to combine the advantages of the iterated Tikhonov regularization, which can be thought as an iterative refinement of an initial solution obtained using the standard Tikhonov regularization, and the multiplicative regularization that eliminates the need for the selection of the regularization parameter while taking into account one's prior knowledge on the sources to identify. The resulting regularization strategy is consequently flexible and computationally efficient. The main conclusion of this study is that the proposed iterated multiplicative regularization should be used only when the measured data are very noisy. In all other cases, the application of the ordinary multiplicative regularization is generally sufficient to obtain consistent reconstructions.

Figure 1 :

 1 Figure 1: Reconstruction of the excitation field at 350 Hz from a measured vibration field with high SNR (30 dB) -(-) Reference, (--) Standard Tikhonov regularization and (-• -) Stationary iterated Tikhonov regularization.

Figure 2 :

 2 Figure 2: Reconstruction of the excitation field at 350 Hz from a measured vibration field -(-) Reference, (--) Multiplicative regularization for q = 0.5 and SNR = 30 dB and (-• -) Multiplicative regularization for q = 0.5 and SNR = 13 dB.

Figure 3 :

 3 Figure 3: Spatial reconstruction of the excitation field at 350 Hz from a vibration field corrupted by an additive Gaussian white noise with (a) a SNR of 30 dB and (b) a SNR of 13 dB -(-) Reference, (--) Multiplicative regularization for q = 0.5 and (-• -) Stationary iterated multiplicative regularization for q = 0.5.

because its value is updated at each iteration of the IRLS