Guillaume Gautier
email: guillaume.gautier@ifsttar.fr

Laurent Mevel

Jean-Mathieu Mencik

Roger Serra

Michael Döhler

Variance analysis for model updating with a finite element based subspace fitting approach

Keywords: Stochastic system identification, Subspace fitting, Uncertainty bounds, Finite element model

Recently, a subspace fitting approach has been proposed for vibration-based finite element model updating. The approach makes use of subspace-based system identification, where the extended observability matrix is estimated from vibration measurements. Finite element model updating is performed by correlating the model-based observability matrix with the estimated one, by using a single set of experimental data. Hence, the updated finite element model only reflects this single test case. However, estimates from vibration measurements are inherently exposed to uncertainty due to unknown excitation, measurement noise and finite data length. In this paper, a covariance estimation procedure for the updated model parameters is proposed, which propagates the data-related covariance to the updated model parameters by considering a first-order sensitivity analysis. In particular, this propagation is performed through each iteration step of the updating minimization problem, by taking into account the covariance between the updated parameters and the data-related quantities. Simulated vibration signals are used to demonstrate the accuracy and practicability of the derived expressions. Furthermore, an application is shown on experimental data of a beam.

Introduction

Linear system identification methods are of interest in mechanical engineering for modal analysis. Using outputonly vibration measurements from structures, Operational Modal Analysis (OMA) has been successfully used as a complementary technique to the traditional Experimental Modal Analysis (EMA) methods [START_REF] Heylen | Modal Analysis Theory and Testing[END_REF][START_REF] Maia | Modal analysis identification techniques[END_REF][START_REF] Peeters | Stochastic system identification for operational modal analysis: a review[END_REF]. With methods originating from stochastic system realization theory for linear systems, estimates of the modal parameters of interest (natural frequencies, damping ratios and observed mode shapes) can be obtained from vibration data. Among these methods, the stochastic subspace identification (SSI) techniques [START_REF] Van Overschee | Subspace Identification for Linear Systems: Theory, Implementation, Applications[END_REF][START_REF] Peeters | Reference-based stochastic subspace identification for output-only modal analysis[END_REF] identify the system matrices of a state-space model, from which the modal parameters are retrieved. Subspace methods are well-suited for the vibration analysis of structures in operation, which is due to the fact that they have excellent theoretical and computational properties, for instance numerical efficiency and robustness, see e.g. [START_REF] Benveniste | Nonstationary consistency of subspace methods[END_REF][START_REF] Döhler | Fast multi-order computation of system matrices in subspace-based system identification[END_REF].

For any system identification method, the estimated parameters are afflicted with variance errors due to finite data length, unknown excitation and measurement noise. The variance of the modal parameter estimates is a most relevant information for assessing their accuracy. It depends on the chosen system identification algorithm. A practical approach for the variance estimation of modal parameters was developed in [START_REF] Pintelon | Uncertainty calculation in (operational) modal analysis[END_REF], where an estimated covariance on the measurements is propagated to the desired parameters by considering a sensitivity analysis. The required sensitivities are derived analytically through first-order perturbation theory, from the data to the identified parameters, and are then computed using the system identification estimates. In [START_REF] Reynders | Uncertainty bounds on modal parameters obtained from stochastic subspace identification[END_REF], the covariance computation scheme for the covariancedriven subspace method (SSI-cov) has been developed, and in [START_REF] Döhler | Efficient multi-order uncertainty computation for stochastic subspace identification[END_REF] a fast and memory efficient implementation of the covariance computation for SSI-cov has been proposed.

The identification of modal parameters or, more generally, the system identification results of a structure in operation find an important application in the calibration of a numerical model of the investigated structure. Finite element (FE) models are used e.g. to verify design specifications, to assess stress fields in structures, to predict vibration levels under prescribed harmonic excitations, to detect abnormal structural behavior in the context of structural health monitoring [START_REF] Friswell | Damage identification using inverse methods[END_REF][START_REF] Simoen | Dealing with uncertainty in model updating for damage assessment: A review[END_REF], and so on. With model updating techniques, the parameters of the FE model are calibrated such that some model properties are close to the truly observed structural properties. Vibration-based FE model updating techniques [START_REF] Friswell | Finite element model updating in structural dynamics[END_REF][START_REF] Mottershead | The sensitivity method in finite element model updating: a tutorial[END_REF] identify model parameters by minimizing a cost function involving the identified and model-based modal parameters (or derived variables thereof). The involved experimental data are subject to uncertainties. In a broad sense, these uncertainties can be classified into two categories of aleatory (irreducible) and epistemic (reducible) uncertainties [START_REF] Der Kiureghian | Aleatory or epistemic? does it matter?[END_REF]. Aleatory uncertainty may result from geometric dimension variability due to manufacturing tolerances or inherent variability of materials such as concrete, while epistemic uncertainty is caused by lack of knowledge (e.g. due to finite number of data samples, undefined measurement noises, unknown excitations, and so on). These uncertainties can be considered in two ways in model updating. First, the uncertainty of the updated parameters can be evaluated based on the uncertainty of the experimental data. Second, the uncertainty of the experimental data can be taken into account in stochastic updating techniques as suggested in [START_REF] Simoen | Dealing with uncertainty in model updating for damage assessment: A review[END_REF]. In this paper, we consider the first way.

While many FE updating methods analyze the impact of aleatory structural uncertainties on the updated parameters (see for instance [START_REF] Khodaparast | Perturbation methods for the estimation of parameter variability in stochastic model updating[END_REF]), we consider the problem of evaluating the uncertainty of the updated FE parameters that result from the uncertainty when estimating parameters from vibration measurements. This statistical uncertainty falls into the category of epistemic uncertainties since it reduces as the number of measurements increases. We consider a recently proposed subspace fitting (SF) approach for FE model updating [START_REF] Gautier | A finite element-based subspace fitting approach for structure identification and damage localization[END_REF], which is a deterministic approach closely linked to subspace identification [START_REF] Swindlehust | A subspace fitting method for identification of linear state-space models, Automatic Control[END_REF]. In this framework, the model parameters of a coarse FE mesh are updated in a minimization problem that consists in correlating a FE-based extended observability matrix with an experimental one that is identified from SSI-cov. Since the experimental observability matrix is estimated from measured vibration data, it is afflicted with uncertainty. In this paper, we propose a covariance analysis of the structural parameters obtained from the SF approach taking into account this uncertainty. The identified observability matrix and its covariance are estimated with SSI-cov [START_REF] Reynders | Uncertainty bounds on modal parameters obtained from stochastic subspace identification[END_REF][START_REF] Döhler | Efficient multi-order uncertainty computation for stochastic subspace identification[END_REF] and propagated to the updated FE parameters by means of a sensitivity analysis. The expressions of the covariances of the FE parameters are an original contribution of the present paper. In particular, we show how the covariance originating from the vibration data is propagated through the iterations of the minimization problem for the solution of the FE parameters, taking into account the covariance between the parameter estimate in each iteration and the data-related quantities. The derived expressions are validated through Monte Carlo simulations. Notice that the covariance estimation from multiple datasets as in Monte Carlo simulations is computationally much more demanding [START_REF] Tondreau | Numerical and experimental analysis of uncertainty on modal parameters estimated with the stochastic subspace method[END_REF], compared to the covariance estimation from a single dataset as in the proposed approach. Also, an application of the approach on a lab experiment is reported.

The paper is organized as follows. In Section 2, the theoretical framework of the FE-based SF method is presented.

The covariance expressions of the updated FE model parameters are derived in Section 3. In Section 4, the covariance estimation procedure is illustrated regarding the vibration data issued from a numerical model of a beam. Results are validated with standard deviations obtained from Monte-Carlo simulations. Finally, in Section 5, the proposed method is applied to model updating using experimental data of a beam.

SF method for model updating

Model updating has the purpose to calibrate the parameters of a FE model such that some model properties are close to the truly observed structural properties. In the proposed SF method, the considered quantity for calibration is the extended observability matrix O of the underlying linear system. The main idea is to correlate -in a least squares sense -the matrix Ô obtained from experimental data with the matrix O h (θ θ θ h ) issued from a FE model of the structure, where θ θ θ h ∈ R n h is the vector of structural parameters of the FE model to be updated.

In this section, the subspace identification method is introduced to obtain Ô from experimental data, and the SF method, based on Ô and O h (θ θ θ h ) for updating the vector of parameters θ θ θ h is presented.

Stochastic subspace identification

The vibration behavior of a linear elastic structure, which is observed at some sensor positions, can be described through the following equations [START_REF] Heylen | Modal Analysis Theory and Testing[END_REF]

           M q(t) + γ q(t) + Kq(t) = v(t), y(t) = H d q(t) + H v q(t) + H d q(t) + w(t), (1)
where M ∈ R n×n , K ∈ R n×n and γ ∈ R n×n are, respectively, the mass, stiffness and viscous damping matrices, t denotes continuous time, q ∈ R n is the vector of displacements, and n is the number of degrees of freedom (DOFs). Boolean matrices H d , H v and H a ∈ R r×n localize the DOFs at which displacements, velocities and accelerations are measured in the vector y ∈ R r . Finally, v ∈ R n and w ∈ R r are white noise vectors of unmeasured excitation forces and measurement noise, respectively.

Sampling system (1) at discrete time instants t = kτ, where τ is the sampling rate, leads to the discrete-time linear time-invariant (LTI) state-space model

           x k+1 = Ax k + v k , y k = Cx k + w k , (2)
where x k = [q(kτ) T q(kτ) T ] T ∈ R 2n is the state vector, y k ∈ R r is the vector of measurement outputs. Also, A ∈ R 2n×2n

and C ∈ R r×2n are state transition and output matrices, respectively, expressed by

A = exp                     0 I n -M -1 K -M -1 γ           τ           , C = H d -H a M -1 K | H v -H a M -1 γ . (3)
The vectors v k ∈ R 2n and w k ∈ R r are zero-mean white noise, i.e., E[v k ] = 0 and E[w k ] = 0, with covariance matrices

[4] E                     w p v p           w q v q           =           Q S S T R           δ pq , (4)
where δ pq is the Kronecker delta. The stochastic process ( 2) is assumed to be stationary, i.e., E[x k+1 x T k ] = Σ where Σ ∈ R 2n×2n does not depend on time k, with zero mean, i.e., E[x k ] = 0. Also, v k and w k are assumed to be independent on the actual state, i.e., E[

x k v T k ] = 0 and E[x k w T k ] = 0. Define the output covariance matrix R i ∈ R r×r as R i = E[y k+i y T k ], (5)
and the state-output covariance matrix G ∈ R 2n×r as

G = E[x k+1 y T k ]. (6)
The output covariance matrices R i can be stacked into a block Hankel matrix H ∈ R (p+1)r×qr , where p and q are chosen such that min(pr, qr) ≥ 2n, as follows

H =                           R 1 R 2 . . . R q R 2 R 3 . . . R q+1 . . . . . . . . . . . . R p+1 R p+2 . . . R p+q                           . ( 7
)
The Hankel matrix H has the factorization property [START_REF] Stoica | Introduction to spectral analysis[END_REF]

H = OC, (8)
where O ∈ R (p+1)r×2n and C ∈ R 2n×qr are the extended observability matrix and controllability matrix, respectively, which are given by

O =                           C CA . . . CA p                           , C = G AG . . . A q-1 G . (9)
Notice that the matrix O is full column rank [START_REF] Van Overschee | Subspace Identification for Linear Systems: Theory, Implementation, Applications[END_REF][START_REF] Gautier | A finite element-based subspace fitting approach for structure identification and damage localization[END_REF]. Within the SSI framework [START_REF] Peeters | Reference-based stochastic subspace identification for output-only modal analysis[END_REF], an estimate Ĥ of the output

Hankel matrix H can be built from N + p + q measurements as

Ĥ = 1 N Y + (Y -) T , (10)
where

Y -=                           y q y q+1 . . . y N+q-1 y q-1 y q . . . y N+q-2 . . . . . . . . . . . . y 1 y 2 . . . y N                           , Y + =                           y q+1 y q+2 . . . y N+q y q+2 y q+3 . . . y N+q+1 . . . . . . . . . . . . y q+p+1 y q+p+2 . . . y N+p+q                           . (11)
From ( 8), it can be readily proved that H and O share the same column space. Hence, an estimate Ô of the extended observability matrix can be obtained from a Singular Value Decomposition (SVD) of Ĥ and its truncation at the order 2n

Ĥ = U∆V T = U 1 U 0           ∆ 1 0 0 ∆ 0                     V T 1 V T 0           ≈ U 1 ∆ 1 V T 1 , (12)
and thus

Ô = U 1 , (13)
where ∆ 1 is the diagonal matrix of the 2n dominant singular values of Ĥ, and U 1 ∈ R (p+1)r×2n is the related matrix of left singular vectors. While estimates ( Ĉ, Â) of the system matrices and subsequently the modal parameters are obtained from Ô in subspace identification, the SF approach requires only the matrix Ô in the following. Thus, any of these state-space systems may be considered. Notice however that the true transformation matrix Ti.e., the one which is considered to express the estimate of O in [START_REF] Friswell | Finite element model updating in structural dynamics[END_REF] in the basis corresponding to (3) -is unknown.

In other words, there exists an unknown matrix T such that ÔT is an estimate of O in (9).

FE-based SF method

In many applications, a numerical FE model is available that approximates the vibration behavior of a structure whose dynamic equilibrium equation is provided by ( 1)-( 2). This model yields an analogous equation of motion

           M h qh (t) + γ h qh (t) + K h q h (t) = v h (t), y(t) = H h d q h (t) + H h v qh (t) + H h d qh (t) + w h (t), (14)
where M h ∈ R n×n , K h ∈ R n×n and γ h ∈ R n×n are, respectively, the mass, stiffness and viscous damping matrices which are issued from the FE modeling of the structure, and the other variables being analogously defined as in [START_REF] Heylen | Modal Analysis Theory and Testing[END_REF]. The matrices M h , K h and γ h depend on a vector

θ θ θ h = [θ h 1 θ h 2 • • • θ h n h ] T ∈ R n h
, which contains the structural parameters that will be updated.

Similarly as in the previous section, sampling [START_REF] Mottershead | The sensitivity method in finite element model updating: a tutorial[END_REF] at rate τ leads to a discrete-time state-space model as in [START_REF] Maia | Modal analysis identification techniques[END_REF] with the system matrices

A h = exp                     0 I n -(M h ) -1 K h -(M h ) -1 γ h           τ           , C h = H h d -H h a (M h ) -1 K h | H h v -H h a (M h ) -1 γ h , (15)
which depend on the parameters {θ h j } j . Also, the related extended observability matrix is expressed by

O h (θ θ θ h ) =                           C h C h A h . . . C h (A h ) p                           . ( 16
)
Outputs y k When the assumed FE model fits the experimental data perfectly, the extended observability matrices O h (θ θ θ h ) issued from the FE model and Ô estimated from experimental data are equal up to a change of basis (see Remark 1), i.e. O h (θ θ θ h ) = ÔT for some matrix T. Then, both matrices share the same column space. However, this is hardly the case in practice and the vector of parameters θ θ θ h needs to be adjusted to satisfy this condition. The key idea behind SF is to adjust the respective subspace defined by O h (θ θ θ h ) such that

Inputs v k

Finite element mesh

θ θ θ h = argmin ||O h (θ θ θ h ) -ÔT|| 2 F , (17)
where || • || F denotes the Frobenius norm.

Assume that the matrix O h depends on a few parameters θ θ θ h . This leads to a minimization problem whose resolution does not require excessive computational times.

The first step in solving this minimization problem is to express the unknown matrix T. It is easy to see that

T = Ô † O h (θ θ θ h ) ( 18
)
is the minimum-norm solution of the linear least squares problem [START_REF] Gautier | A finite element-based subspace fitting approach for structure identification and damage localization[END_REF] for fixed θ θ θ h , where Ô † is the Moore-Penrose pseudoinverse of Ô. By replacing the expression of T in the original problem [START_REF] Gautier | A finite element-based subspace fitting approach for structure identification and damage localization[END_REF], the minimization problem takes the following form [START_REF] Golub | Separable nonlinear least squares: the variable projection method and its applications[END_REF]

θ θ θ h = argmin ||(I (p+1)r -Ô Ô † )O h (θ θ θ h )|| 2 F . (19)
Using vector notations with the column stacking vectorization operator vec{•}, the minimization problem turns out to be

θ θ θ h = argmin ||r|| 2 2 , (20)
where r ∈ R 2(p+1)rn has the meaning of a residual vector, defined as

r = I 2n ⊗ (I (p+1)r -Ô Ô † ) vec O h (θ θ θ h ) , (21)
and || • || 2 is the Euclidean vector norm. One of the most popular and effective algorithms for solving least squares problems like [START_REF] Stoica | Introduction to spectral analysis[END_REF] is the Gauss-Newton method [START_REF] Wang | Gauss-newton method[END_REF]. The method is based on a second-order expansion of the objective function ||r|| 2 2 about some approximated values of the parameters θ θ θ h , and iteration steps so as to find the local extrema of ||r|| 2 2 . Assuming an initial deterministic parameter value θ θ θ h 0 , the k-th iteration of the Gauss-Newton algorithm can be written as

θ θ θ h k = θ θ θ h k-1 -J † r k r k (22) for k ≥ 1, where θ θ θ h k = [θ h 1,k • • • θ h n h ,k ] T and θ θ θ h k-1 = [θ h 1,k-1 • • • θ h n h ,k-1 ]
T are the vectors of structural parameters identified at iterations k and k -1of the Gauss-Newton algorithm, respectively. Also, r k is the residual vector at iteration k

r k = I 2n ⊗ (I (p+1)r -Ô Ô † ) vec O h (θ θ θ h k-1 ) ,
J † r k is the Moore-Penrose pseudoinverse of the Jacobian matrix J r k ∈ R 2(p+1)rn×n h , namely the derivative of r k with respect to the approximated values of θ θ θ h

J r k = I 2n ⊗ (I (p+1)r -Ô Ô † ) J O h ,θ θ θ h k-1 . (23)
In [START_REF] Wang | Gauss-newton method[END_REF],

J O h ,θ θ θ h k-1 is given by J O h ,θ θ θ h k-1 = ∂vec{O h (θ θ θ h )} ∂θ h 1 . . . ∂vec{O h (θ θ θ h )} ∂θ h n h θ θ θ h =θ θ θ h k-1 . (24)
A simple way to calculate ∂vec{O h (θ θ θ h )} ∂θ h j θ θ θ h =θ θ θ h k-1 is by numerical differentiation using the central difference theorem [START_REF] Andersen | Estimation of modal parameters and their uncertainties[END_REF] ∂vec{O

h (θ θ θ h )} ∂θ h j θ θ θ h =θ θ θ h k-1 = 1 2δθ h j,k-1 vec O h (θ θ θ h k-1 + δ j θ θ θ h k-1 ) -vec O h (θ θ θ h k-1 -δ j θ θ θ h k-1 ) , (25)
where δ j θ θ θ h k-1 ∈ R n h is a vector whose components are zero, except for the j-th component which is equal to δθ h j,k-1 .

Mode-based approach

Notice that the computation of the minimization problem can be cumbersome even though a reduced set of parameters θ θ θ h is involved, which is explained by the fact that the size of the matrix O h (θ θ θ h ) is usually large. To solve this issue, a model reduction technique is employed which makes use of the concept of reduced mode expansion. In this framework, the vector of nodal displacements q h of the FE mesh of the structure is expressed in terms of mode shape vectors {φ h j } j . These are solutions of the following generalized eigenproblem

K h φ h j = (ω h j ) 2 M h φ h j , (26)
where {ω h j } j are to be understood as the eigenpulsations of the vibrating structure. The reduced mode expansion hence consists in expressing the vector of nodal displacements q h using a reduced set of mode shape vectors {

φh j } j=1•••n r which is extracted from the full set {φ h j } j=1,••• ,n , i.e., { φh j } j=1,••• ,n r ⊂ {φ h j } j=1,••• ,n
where n r ≪ n. The reduced mode expansion is expressed as follows

q h (t) ≈ n r j=1 φh j qh j (t) = Φh qh (t), (27)
where { qh j } j have the meaning of generalized coordinates; also, Φh and qh (t) are expressed as

Φh = [ φh 1 • • • φh n r ] and qh (t) = [ qh 1 (t) • • • qh n r (t)] T .
As a result, a reduced extended observability matrix Õh ( θ θ θ h ), having a size of (p + 1)r × 2n r , can be expressed as

Õh ( θ θ θ h ) =                           Ch Ch Ãh . . . Ch ( Ãh ) p                           , ( 28
)
where Ãh and Ch are matrices whose determination involves projecting the second-order differential equation [START_REF] Mottershead | The sensitivity method in finite element model updating: a tutorial[END_REF] onto the reduced basis { φh j } j through Galerkin procedure. This is done by left multiplying this equation by ( Φh ) T and invoking the orthogonality properties of the mode shape vectors -i.e., ( Φh ) T M h Φh = I and ( Φh ) T K h Φh = diag{( ωh j ) 2 } j -as well as the fact that the matrix ( Φh ) T γ h Φh is assumed to be diagonal. As a result, one has

qh (t) + diag{2 ξh j ωh j } j qh (t) + diag{( ωh j ) 2 } j qh (t) = ṽh (t), (29)
where the parameters { ξh j } j denote modal damping ratios. Expressions of the matrices Ãh and Ch in the corresponding discrete-time state-space model follow as

Ãh = exp                     0 I -diag{( ω j h ) 2 } j -diag{2 ξh j ωh j } j           τ           , Ch = H h d Φh -H h a Φh diag{( ωh j ) 2 } j | H h v Φh -H h a Φh diag{2 ξh j ωh j } j . (30)
By considering Eqs. ( 16) and ( 28), the minimization problem involving the reduced extended observability matrix Õh (θ θ θ h ) can be expressed as

θ θ θ h = argmin ||r|| 2 2 , (31)
where

r = I 2n r ⊗ (I (p+1)r -Ô Ô † ) vec Õh ( θ θ θ h ) . (32)
Similarly as in ( 22), the Gauss-Newton iteration is defined as

θ θ θ h k = θ θ θ h k-1 -J † rk rk , (33)
where the residual rk writes as

rk = I 2n r ⊗ (I (p+1)r -Ô Ô † ) vec Õh ( θ θ θ h k-1 ) . (34)
It is understood that the parameters θ θ θ h are almost equal to those involved in the original minimization problem (20), i.e., θ θ θ h ≈ θ θ θ h . For simplicity of notation, only θ θ θ h instead of θ θ θ h is used in the following, neglecting the tilde sign.

Uncertainty quantification

Analyzing the statistical uncertainty of identified parameters is mandatory to assess the quality of the estimates from the data. When estimated from a finite number of data samples, not the "true" parameters of the system are obtained, but estimates that are naturally subject to variance errors depending on the data and the estimation method, due to unknown noise inputs, measurement noise and finite data length.

A variance analysis of the system matrices obtained from stochastic subspace identification has been made in [START_REF] Reynders | Uncertainty bounds on modal parameters obtained from stochastic subspace identification[END_REF][START_REF] Döhler | Efficient multi-order uncertainty computation for stochastic subspace identification[END_REF] by computing a sample covariance that is closely linked to the measurements and propagating it to the system matrices through a sensitivity analysis. In this section, the variance analysis of the updated structural parameters from the FE-based SF method of the previous section is made.

The following notation is used. Let Y be a matrix-valued function of X. Expression of its covariance follows from a Taylor approximation

vec{Y( X)} ≈ vec{Y(X)} + J Y,X vec{ X -X} ⇒ cov(Y( X)) ≈ J Y,X ΣX J T Y,X , (35)
where J Y,X is the sensitivity matrix, defined by J Y,X = ∂vec{Y(X)}/∂vec{X}. A consistent estimate is obtained by replacing, in the sensitivity matrix, the theoretical variables X with consistent estimates X issued from data. Using delta notation, theoretical first-order perturbations are defined, yielding vec{∆Y} = J Y,X vec{∆X}.

Using this relationship, a perturbation ∆ Ĥ of the Hankel matrix Ĥ is propagated towards the observability matrix Ô and ultimately the structural parameters through the Gauss-Newton algorithm to obtain the respective sensitivity matrices. Together with a sample covariance ΣH of the Hankel matrix, which can be easily computed on the measurements, the covariance of the structural parameters is obtained through (35). Since the output-covariances of the data and the Hankel matrix estimate are asymptotically Gaussian (for N → ∞) [START_REF] Hannan | Multiple time series (wiley series in probability and mathematical statistics[END_REF], the statistical delta method ensures that each variable in the computation chain -down to the desired structural parameters -is also asymptotically Gaussian [START_REF] Casella | Statistical inference, thomson learning[END_REF], and that the respective covariance estimate in Eq. ( 35) is exact for N → ∞.

The expressions of the covariances of the Hankel matrix and observability matrix have been derived in past works [START_REF] Pintelon | Uncertainty calculation in (operational) modal analysis[END_REF][START_REF] Reynders | Uncertainty bounds on modal parameters obtained from stochastic subspace identification[END_REF] and are recalled here for the sake of clarity.

Covariance of the Hankel matrix

The derivation of an estimate of the covariance of the Hankel matrix, namely ΣH , is achieved by splitting the data matrices Y + and Y -in [START_REF] Friswell | Damage identification using inverse methods[END_REF] into n b blocks, as follows

Y + = Y + 1 Y + 2 . . . Y + n b , Y -= Y - 1 Y - 2 . . . Y - n b . (37)
For each pair of sub-matrices (Y + j , Y - j ), the corresponding Hankel matrix Ĥ j can be estimated as in Eq. ( 10) as

Ĥ j = 1 N b Y + j (Y - j ) T , (38)
where N b = N/n b . As a result, Ĥ = (1/n b ) n b j=1 Ĥ j and the estimate ΣH is given by

ΣH = 1 n b (n b -1) n b j=1 vec{ Ĥ j } -vec{ Ĥ} vec{ Ĥ j } -vec{ Ĥ} T .
(39)

Covariance of the observability matrix

A perturbation ∆H of the Hankel matrix is propagated towards the observability matrix. In order to obtain ∆ Ô from an estimate of O (see Eq. ( 13)), the sensitivities of the left singular vectors U 1 (see Eq. ( 13)) are to be derived, as follows [START_REF] Döhler | Efficient multi-order uncertainty computation for stochastic subspace identification[END_REF] vec

{∆U 1 } = J U 1 ,H vec{∆ Ĥ} =                   B 1 C 1 . . . B 2n C 2n                   vec{∆ Ĥ}, ( 40
)
where

B i =           I (p+1)r + Ĥ σ i D i           ĤT σ i -           0 qr-1,(p+1)r u T i                     Ĥ σ i D i           , D i =           I qr +           0 qr-1,(p+1)r 2v T i           - ĤT Ĥ σ 2 i           -1 (41)
and

C i = 1 σ i           (I (p+1)r -u i u T i )(v T i ⊗ I (p+1)r ) (I qr -v i v T i )(I qr ⊗ u T i )           , ( 42
)
where u i and v i stand for the i-th left and right singular vectors, respectively. Then, an estimate of the covariance matrix of Ô can be derived from (35) as

ΣO = J U 1 ,H ΣH J T U 1 ,H . ( 43
)

Covariance matrix of the structural parameters

The purpose of the present section is to link the parameter perturbation ∆θ θ θ h k with the perturbation of the observability matrix. The main result is that

∆θ θ θ h k = M k vec{∆ Ô}, ( 44
)
where M k ∈ R n h ×2n(p+1)r . It is defined iteratively by

M k = M k-1 -(r T k ⊗ I n h )L k N k -J † rk Q k ( 45
)
where M 0 = 0 n h ,2n(p+1)r , while matrices L k , N k and Q k will be detailed in the following. Hence, an estimate of the covariance matrix of θ θ θ h k is given by Σθ

θ θ h k = M k ΣO M T k , (46)
where Σθ θ θ h k is the parameter covariance matrix for the k-th iteration. In a constructive proof, relation (44) is derived by induction. Obviously, (44) holds for k = 0, since ∆θ θ θ h 0 = 0 n h due to the deterministic nature of the initial value. Then, for a fixed k with k ≥ 1, assume that (44) holds for k -1, i.e.

∆θ θ θ h k-1 = M k-1 vec{∆ Ô}, ( 47
)
and the relation (44) for k is proved as follows.

The starting point is to consider a perturbation on the iteration θ θ θ h k = θ θ θ h k-1 -J † rk rk in (33). Separating zeroth-and first-order terms and using the fact that vec{AB} = (B T ⊗ I a )vec{A} for any matrices A ∈ R a×b and B ∈ R b×c [START_REF] Brewer | Kronecker products and matrix calculus in system theory[END_REF], this yields

∆θ θ θ h k = ∆θ θ θ h k-1 -(r T k ⊗ I n h )vec{∆J † rk } -J † rk ∆r k . (48)
From Appendix B, both vectors vec{∆J † rk } and ∆r k can be written in terms of vec{∆ Ô} as

vec{∆J † rk } = L k N k vec{∆ Ô}, ( 49
)
∆r k = Q k vec{∆ Ô}, (50) where
• Following the calculations in Appendix B with X = Õh (θ θ θ h k-1 ) and Y = rk , Q k is expressed by

Q k = -Ô † Õh (θ θ θ h k-1 ) T ⊗ I (p+1)r -Õh (θ θ θ h k-1 ) T ⊗ Ô K + I 2n r ⊗ (I (p+1)r -Ô Ô † ) J Õh ,θ θ θ h k-1 M k-1 , (51)
where J Õh ,θ θ θ h k-1 is defined as in [START_REF] Hannan | Multiple time series (wiley series in probability and mathematical statistics[END_REF] and matrix K satisfies vec{∆ Ô † } = Kvec{∆ Ô} with

K = {-[( Ô Ô † ) T ⊗ ( ÔT Ô) -1 ] + [I (p+1)r ⊗ ( ÔT Ô) -1 ]}P (p+1)r,2n -[( Ô † ) T ⊗ Ô † ] (52)
following from Appendix A with permutation matrix P (p+1)r,2n defined in (A.3).

• Matrix L k satisfies vec{∆J † rk } = L k vec{∆J rk } with L k = {-[(J rk J † rk ) T ⊗ (J T rk J rk ) -1 ] + [I 2(p+1)rn r ⊗ (J T rk J rk ) -1 ]}P 2(p+1)rn r ,n h -[(J † rk ) T ⊗ J † rk ] (53)
following from Appendix A.

• Matrix N k satisfies vec{∆J rk } = N k vec{∆ Ô} with N k =                           N 1,k N 2,k . . . N n h ,k                           , (54)
where

N j,k is derived from Appendix B with X = J * Õh ,θ h j,k-1
, where J *

Õh ,θ h j,k-1 ∈ R (p+1)r×2n r is such that vec{J * Õh ,θ h j,k-1 } = J Õh ,θ h j,k-1
, and

Y = I 2n r ⊗ (I (p+1)r -Ô Ô † ) J Õh ,θ h j,k-1
being the j-th column of J rk in [START_REF] Andersen | Estimation of modal parameters and their uncertainties[END_REF]. It follows

N j,k = -Ô † J * Õh ,θ h j,k-1 T ⊗ I (p+1)r -J * Õh ,θ h j,k-1 T ⊗ Ô K + I 2n r ⊗ (I (p+1)r -Ô Ô † ) J J * Õh ,θ h j,k-1 ,θ θ θ h k-1 M k-1 . (55) Matrix J J * Õh ,θ h j,k-1 ,θ θ θ h k-1 satisfies vec{∆J * Õh ,θ h j,k-1 } = J J * Õh ,θ h j,k-1 ,θ θ θ h k-1 vec{∆θ θ θ h k-1
} and is obtained by numerical differentiation similar to [START_REF] Casella | Statistical inference, thomson learning[END_REF].

Finally, introducing Eqs. (47), ( 49) and (50) into Eq. (48) proves the assertion for the perturbation on the parameter at iteration k in (44) with M k defined as in Eq. ( 45). The estimate of the covariance matrix of the structural parameters Σθ θ θ h k (see Eq. ( 46)) can thus be obtained for the k-th iteration of the Gauss-Newton algorithm.

Numerical validation

Introduction

The proposed approach is used to estimate the variances of some mechanical parameters of a vibrating cantilever beam whose measurement data result from a fine FE mesh with known characteristics, see Figure 2. The characteristics of the beam are: length of 1m, cross-section of 0.0249m × 0.0053m, Young's modulus of 200GPa, density of 7,850kg/m 3 and Poisson's ratio of 0.3. The damping matrix is defined so that γ γ γ = M + 10 -6 K, where M and K are the mass and stiffness matrices of the beam. The reference FE mesh is composed of 3D tetrahedral elements with three DOFs per node and contains 5,505 DOFs, see Figure 2. The eigenfrequencies of the 3D FE beam are reported in Table 1. The output data (vector y k ) refer to the transverse displacement of the structure which is recorded by means of one displacement sensor located at 1m from the clamped end. The sampling frequency used for measuring these data is 1,280Hz and the number of samples is N = 12,800. A Gaussian random noise is added to the output with a signal to noise ratio of 20dB.

3D beam

1D beam

E = 160GPa E = 200.21GPa Mode f [Hz] f [Hz] Err [%] f [Hz] Err [%] 1

Young's modulus updating and standard deviation estimation

The FE-based SF approach and uncertainties quantification are investigated. In this framework, a coarse FE mesh based on the Euler-Bernoulli beam theory is considered as shown in Figure 3. This mesh is made up of 10 Euler- where ΣE (i) is the estimate of the variance of the Young's modulus estimated from dataset i, see Eq. (46). It is evaluated by using the procedure explained in Section 3, where n b = 100, and its mean is found to be σ E ∆,S F = 0.53GPa at the last iteration of the Gauss-Newton procedure. This value should be close to the sample standard deviation of the Young's modulus obtained from the Monte-Carlo simulations, which is given by

σ E MC,S F = 1 n e -1 n e i=1 (E (i) -E) 2 , (58)
and is found to be σ E MC,S F = 0.61GPa. Indeed, this value is close to the estimated standard deviation with the proposed method. The difference may be explained by the approximation of the covariance estimation in (35), which is getting more precise when the number of samples N in the datasets increases. In fact, N = 12,800 for a sampling frequency of 1,280Hz may be relatively low for convergence.

The distribution of the estimated standard deviations of the Young's modulus in (57) is displayed in Figure 5. The sample standard deviation of the estimated standard deviations can be computed from the Monte-Carlo simulations and is found to be σ σ E ∆,S F = 0.13GPa. This represents about 25% of the estimated standard deviation. In other words, the estimation error on the standard deviation from a single dataset is about 25%, which is a reasonable value for a second-order statistic. Furthermore, considering that E is the only unknown parameter of the FE model, the standard deviations of the eigenfrequencies -namely, f j -can be obtained from the estimate of the variance of the Young's modulus as (see Eq. ( 36))

Σ f j = J f j ,E ΣE J T f j ,E , (59)
where

J f j ,E = ∂ f j ∂E . ( 60
)
Their means are displayed in Table 2. In comparison, the sample standard deviations of the eigenfrequencies can be obtained from the Monte-Carlo simulations as

σ f j MC,S F = 1 n e -1 n e i=1 ( f (i) j -f j ) 2 , f j = 1 n e n e j=1 f (i) j . (61)
They are reported in Table 2. It can be seen that the uncertainty bounds obtained from the proposed approach are close to those provided by Monte-Carlo simulations, as expected.

Notice that the mean of the CPU time involved to solve the minimization problem (56) and to calculate the estimate of the standard deviation of the Young's modulus is 16.25s, which appears to be quite reasonable. Within the Monte Carlo framework, the CPU time involved to solve the 1,000 minimization problems is 16,942s, which is significantly higher than that of the proposed approach.

mode frequency [Hz] σ f

(i) ∆,S F σ f (i)

Experimental application

An application on a real experimental beam is carried out to investigate further the accuracy of the uncertainty estimation. The experimental structure is shown in Figure 6, with dimensions which are similar to those depicted in Section 4. The density, Poisson's ratio and damping of the structure are supposed to be known. The excitation is a random noise which is generated by a shaker located at 3.5cm from the clamped end. The output data are recorded through a mono-axial accelerometer located on the free end of the beam.

Output data are sampled at a rate of 32,768Hz. The number of samples is set to be N = 496,000, which results in a measurement time of approximately 15s. Afterwards, the data are filtered with a low-pass filter with a cut-off frequency of 300Hz and resampled at 512Hz, which reduced the number of samples to N = 7750. The pre-processed data are used to update the Young's modulus of the 1D FE beam (see Figure 3), which is unknown. The parameters used in the updating procedure are p = 20, n = 8 and n b = 100.

As a whole, 70 runs are conducted in order to estimate 70 values of the Young's modulus by means of the proposed approach. The mean of the updated values of the Young's modulus, which results from an initial guess E = 160GPa, is found to be E = 188.50GPa after 30 iterations as shown in Figure 7. The mean of the standard deviation is found to be σ E ∆,S F = 0.13GPa. Again, the Monte-Carlo procedure is carried out in order to validate the proposed approach. The sample standard deviation of the Young's modulus is computed from the 70 estimations and it is found to be σ E MC,S F = 0.15GPa, see Figure 8. As expected, this value appears to be closely linked with the standard deviation issued from the proposed approach.

Concluding remarks

In this paper, a theoretical framework has been proposed to quantify the variance in a SF model updating procedure. In this framework, a strategy has been proposed to calculate covariance estimates when updating the parameters of a FE model from measurement data. The proposed method has been successfully validated through Monte Carlo experiments conducted on numerical and experimental beams. It has been shown that the variance estimates resulting from the proposed method are in good agreement with those resulting from sample statistics of the Monte Carlo procedure.

Variance estimation procedures have already been developed for subspace methods. Such works only addressed the simpler case of global system identification algorithms. In the current paper, first variance results were obtained for the case of iterative optimization algorithms. It is a promising approach that could be applied to other algorithms based on Gauss-Newton optimizations such as maximum likelihood approaches.

Every uncertainty quantification scheme aims at obtaining uncertainty information based on prior information and collected measurements. In a Bayesian model updating scheme as e.g. outlined in [START_REF] Simoen | Dealing with uncertainty in model updating for damage assessment: A review[END_REF], the prior information, i.e. the non-calibrated FE model parameters, are considered to be random variables with a given probability distribution, and the final uncertainty is derived in a likelihood framework based on the Bayes formula in a direct computation. In the proposed approach, only the uncertainty related to the measurement data is considered, which does not require prior information about the FE parameters as in the Bayesian framework. The estimated parameters issued from SSI methods follow an asymptotically Gaussian distribution [START_REF] Bauer | Consistency and asymptotic normality of some subspace algorithms for systems without observed inputs[END_REF], leading to an asymptotically Gaussian distribution of the updated FE model parameters. Then, a prior Gaussian distribution of the initial FE model parameters could easily be incorporated in the proposed approach.

Remark 1 .

 1 The LTI state-space system (2) is defined up to a change of basis. Hence, for any invertible transformation matrix T, the matrix systems defined by (C, A) or (CT, T -1 AT) are equivalent and share the same modal parameters (eigenfrequencies, observed mode shapes and modal dampings), with respective observability matrices O or OT.

Figure 1 :

 1 Figure 1: Vibrating structure (left) and corresponding FE mesh discretization (right).

Figure 2 :

 2 Figure 2: 3D FE model of the beam.

Figure 4 :

 4 Figure 4: Young's modulus distribution obtained from Monte-Carlo simulations.

Figure 5 :

 5 Figure 5: Distribution of the variance of the Young's modulus.

Figure 6 :

 6 Figure 6: Experimental set-up of the beam.

Figure 7 :

 7 Figure 7: Overview of the iterative updating procedure, regarding the Young's modulus of the experimental beam.

Figure 8 :

 8 Figure 8: Young's modulus distribution obtained from Monte-Carlo simulations.

  The impact of uncertainties still constitutes an open and tough challenge in FE model updating. Future works may involve incorporating the parameter covariance matrix, at each iteration of the FE-based SF identification algorithm, in order to regularize the procedure.• The first term in Eq. (B.3) can be developed as[I 2n r ⊗ (∆[ Ô Ô † ])]vec{X} = [I 2n r ⊗ ∆ Ô Ô † ]vec{X} + [I 2n r ⊗ Ô∆ Ô † ]vec{X} (B.4)where -by using the relation (I 2n r ⊗ AB)vec{C} = ((BC) T ⊗ I (p+1)r )vec{A}, where A ∈ R (p+1)r×2n , B ∈ R 2n×(p+1)r and C ∈ R (p+1)r×2n r [26] [I 2n r ⊗ ∆ Ô Ô † ]vec{X} = [( Ô † X) T ⊗ I (p+1)r ]vec{∆ Ô} (B.5)-by using the relation (I 2n r ⊗ AB)vec{C} = (C T ⊗ A)vec{B} [I 2n r ⊗ Ô∆ Ô † ]vec{X} = (X T ⊗ Ô)vec{∆ Ô † } (B.6) From Eq. (A.4) it comes that vec{∆ Ô † } = K vec{∆ Ô}, (B.7)whereK = {-[( Ô Ô † ) T ⊗ ( ÔT Ô) -1 ] + [I (p+1)r ⊗ ( ÔT Ô) -1 ]}P (p+1)r,2n -[( Ô † ) T ⊗ Ô † ]. (B.8) Eq. (B.6) is then formulated as [I 2n r ⊗ Ô∆ Ô † ]vec{X} = (X T ⊗ Ô)Kvec{∆ Ô}. (B.9)• Using Eq. (B.2), the second term in Eq. (B.3) can be developed as[I 2n r ⊗ (I (p+1)r -Ô Ô † )]vec{∆X} = [I 2n r ⊗ (I (p+1)r -Ô Ô † )]J X,θ θ θ h k-1 ∆θ θ θ h k-1 . (B.10)By using Eq. (47), i.e.,∆θ θ θ h k-1 = M k-1 vec ∆ Ô , (B.11)Eq. (B.10) is rewritten as[I 2n r ⊗ (I (p+1)r -Ô Ô † )]vec{∆X} = [I 2n r ⊗ (I (p+1)r -Ô Ô † )]J X,θ θ θ h k-1 M k-1 vec ∆ Ô (B.12) By introducing Eqs. (B.4), (B.5), (B.9) and (B.12) in Eq. (B.3), the relation between vec{∆Y} and vec{∆ Ô} writes as ∆Y = J Y, Ôvec{∆ Ô}, (B.13) where J Y, Ô = -[( Ô † X) T ⊗ I (p+1)r ] -[X T ⊗ Ô]K + [I 2n r ⊗ (I (p+1)r -Ô Ô † )]J X,θ θ θ h k-1 M k-1 . (B.14)

Table 1 :

 1 Comparison between the eigenfrequencies of the 3D FE beam and those of the updated 1D model.

		4.3293	3.8653	10.72	4.3238	0.13
	2	27.1280	24.2241	10.70	27.0975	0.11
	3	75.9513	67.8433	10.66	75.8907	0.08
	4	148.8102	133.0385	10.60	148.8193	0.01
	5	245.9611	220.2667	10.45	246.3943	0.18
	6	367.3624	329.9845	10.17	369.1267	0.48
	7	512.9839	462.9843	9.75	517.9027	0.96

Table 2 :

 2 Standard deviation of the eigenfrequencies.

Bernoulli beam elements of same length with six DOFs per node. As a result, this coarse FE mesh contains n = 60 DOFs, which appears to be small compared to the number of DOFs used to discretize the 3D structure. The number of DOFs involved can be further reduced by considering the model reduction strategy depicted in Section 2.3. In this framework, a reduced number of modes are used to describe the dynamic behavior of the 1D beam. The manner by which these modes are selected follows from the procedure explained in [START_REF] Gawronski | Discrete-time norms of flexible structures[END_REF], which yields n r = 7 modes. The issue here consists in updating the Young's modulus -namely, E -of the coarse 1D FE model, which is assumed to be unknown. This yields a single-parameter SF minimization problem, as follows

where O h is the observability matrix issued from the 1D FE model of the beam (Figure 3). To solve Eq. ( 56), the algorithm procedure depicted in Section 2 can be easily implemented using MATLAB R where p = 50 and 2n r = 14.

In order to validate the proposed method, a comparison is made with Monte-Carlo simulations [START_REF] Mooney | Monte carlo simulation[END_REF], where a large number n e of output-only datasets are generated from random excitation. For each dataset, the parameter i=1 E (i) -, which results from an initial value E = 160GPa, is found to be E = 200.21GPa after 20 iterations as shown in Figure 4. The eigenfrequencies of the 1D beam, which result from the initial and updated Young's modulus, are reported in Table 1 and compared with those of the reference 3D model.

From the uncertainty quantification method described in Section 3, the standard deviation of Young's modulus for each dataset can be estimated by

Appendix A. Vectorization of Moore-Penrose pseudoinverse matrix perturbation Consider a matrix X ∈ R a×b which is assumed to be full column rank. Then the Moore-Penrose pseudoinverse X † of X is expressed as

The first-order perturbation of X † writes Hence: where X ∈ R (p+1)r×2n r , Ô and X are afflicted with uncertainties, and assume

2

The goal is to derive the relationship between vec{∆Y} and vec{∆ Ô}.


Pierre Boulet

Philippe Devienne

Pierre Falez

Guillermo Polito

Mahyar Shahsavari

Pierre Tirilly

N2S3, an Open-Source Scalable Spiking Neuromorphic Hardware Simulator

One of the most promising approaches to overcome the end of Moore's law is neuromorphic computing. Indeed, neural networks already have a great impact on machine learning applications and offer very nice properties to cope with the problems of nanoelectronics manufacturing, such as a good tolerance to device variability and circuit defects, and a low activity, leading to low energy consumption. We present here N2S3 (for Neural Network Scalable Spiking Simulator), an open-source simulator that is built to help design spiking neuromorphic circuits based on nanoelectronics. N2S3 is an event-based simulator and its main properties are flexibility, extensibility, and scalability. One of our goals with the release of N2S3 as open-source software is to promote the reproducibility of research on neuromorphic hardware. We designed N2S3 to be used as a library, to be easily extended with new models and to provide a user-friendly special purpose language to describe the simulations.

I. Introduction

Neuromorphic computing has the potential to bring very low power computation to future computer architectures and embedded systems [START_REF] Merolla | A million spiking-neuron integrated circuit with a scalable communication network and interface[END_REF]. Indeed, parallel neuromorphic computing, by performing computation and storage on the same devices, can overcome the von Neumann bottleneck. Several large projects are based on neuromorphic systems, such as the EU Human Brain Project [START_REF] Markram | The human brain project[END_REF], the DARPA/IBM SYNAPSE project [START_REF] Ananthanarayanan | The cat is out of the bag: cortical simulations with 109 neurons, 1013 synapses[END_REF], and deep learning research led by Google and Facebook, among others.

Recently, emerging nano-scale devices have demonstrated novel properties for producing new memories and unconventional processing units. One of those is the memristor, that was hypothetically presented by Leon Chua in 1971 [START_REF] Chua | Memristor-The missing circuit element[END_REF]; after a few decades, HP was the first to announce a successful memristor fabrication [START_REF] Strukov | The missing memristor found[END_REF]. The unique properties of memristors, such as extreme scalability, flexibility thanks to their analog behavior, and their ability to remember their last state, make memristors very promising candidates to be used as synapses in Spiking Neural Networks (SNN) [START_REF] Jo | Nanoscale Memristor Device as Synapse in Neuromorphic Systems[END_REF].

Given their potential of very low power execution and their capability to handle natural signals, we focus on SNN. A comprehensive introduction and literature review about SNN was published by Paugam-Moisy and Bohte in 2012 [START_REF] Paugam-Moisy | Computing with Spiking Neuron Networks[END_REF]. The authors explore the computational capabilities of SNN, their learning capabilities, and their simulation.

Brette et al. [START_REF] Brette | Simulation of networks of spiking neurons: A review of tools and strategies[END_REF] surveyed and discussed the existing work on SNN simulation in 2007. All the simulators discussed in this report as well as the more recent Brian [START_REF] Goodman | Brian: a simulator for spiking neural networks in Python[END_REF] target the simulation of biological SNN. More recently, Bichler et al. [START_REF] Bichler | Design exploration methodology for memristor-based spiking neuromorphic architectures with the Xnet event-driven simulator[END_REF] proposed Xnet, a C++ event-driven simulator dedicated to the simulation of hardware SNN. In our work, we share the goals of Xnet: "intermediate modeling level, between lowlevel hardware description languages and high-level neural networks simulators used primarily in neurosciences", and "the integration of synaptic memristive device modeling, hardware constraints and any custom features required for the targeted application". In addition to these goals, we put an emphasis on flexibility and usability to allow the study of various kinds of hardware designs (possibly by other researchers than us), scalability to simulate large hardware SNNs, and software engineering best practices (robust and extensible software architecture for maintainability, extensive test suite, continuous integration, open-source distribution). In this report we present N2S3 (Neural Network Scalable Spiking Simulator, pronounced "Nessie", hence its logo), an event-driven simulator dedicated to the exploration of hardware SNN architectures. The internals of N2S3 are based on the exchanges of messages between concurrent actors [START_REF] Wyatt | Akka Concurrency[END_REF], mimicking the exchange of spikes between neurons. N2S3 has been developed from the ground up for extensibility, allowing to model various kinds of neuron and synapse models, various network topologies (especially, it is not restricted to feed-forward networks), various learning procedures, various reporting facilities, and to be user-friendly, with a domain specific language to easily express and run new experiments. It is available as open-source software at https://sourcesup.renater.fr/wiki/n2s3 so that its users can share their models and experimental settings to enable others to reproduce their results. In this spirit, N2S3 is distributed with the implementation of two "classical" experiments: handwritten digit recognition on the MNIST dataset [START_REF] Lecun | Gradient-based learning applied to document recognition[END_REF], [START_REF] Lecun | Convolutional networks for images, speech, and time series[END_REF] and the highway vehicle counting experiment [START_REF] Bichler | Extraction of temporally correlated features from dynamic vision sensors with spike-timing-dependent plasticity[END_REF].

In the remainder of this report, we will first detail the fundamental architectural choices of N2S3, then the models that are already implemented in N2S3, and finally the features that N2S3 provides to implement and run experiments.

II. Event-Driven Simulation Architecture

A. Event-Driven vs Clock-Driven Simulation SNN are essentially defined by standard differential equations, but because of the discontinuities caused by the spikes, designing an efficient simulation of spiking neural networks is a nontrivial problem. There are two families of simulation algorithms: event-based simulators and clock-based ones. Synchronous or "clock-driven" simulation simultaneously updates all the neurons at every tick of a clock, and is easier to code, especially on GPUs, for getting an efficient execution of data-parallel learning algorithms. Event-driven simulation behaves more like hardware, in which conceptually concurrent components are activated by incoming signals (or events).

Event-driven execution is particularly suitable for untethered devices such as neurons and synapses, since the nodes can be put into a sleep mode to preserve energy when no interesting event is happening. Energy-aware simulation needs information about active hardware units and event counters to establish the energy usage of each spike and each component of the neural network. Furthermore, as the learning mechanisms of spiking neural networks are based on the timings of spikes, the choice of the clock period for a clock-based simulation may lead either to imprecision or to a higher computational cost.

There is also a fundamental gap between this event-driven execution model and the clock-based one: the first one is independent on the hardware architecture of computers on which it is running. So, event-driven simulators can naturally run on a grid of computers, with the caveat of synchronization issues in the management of event timings.

B. Technological Choices: Scala and Akka

To address our concurrency and distributability requirements (i.e., ability to scale out a simulation on several computers to handle large networks) we have chosen to use the Scala programming language [START_REF] Odersky | An overview of the Scala programming language (2nd Edition)[END_REF] along with the Akka actor library [START_REF] Wyatt | Akka Concurrency[END_REF].

C. Software Architecture

The architecture of N2S3 is organized in Akka actors. Our choice for an actor model aims mainly for the distribution and scalability of large neural networks. Each entity of the simulation is deployed in a concrete actor that is assigned to a concrete machine at runtime. The core of a typical neural network simulation in N2S3 is composed mainly by the following network entities: a) Containers: Containers are network entities in charge of organizing the network structure. Their main responsibility is to contain network entities and dispatch messages to their children.

b) Neurons: The basic building blocks of a neural network. In N2S3, a neuron is composed of both its nucleus (the soma) and its incoming connections (dentrites and associated synapses).

Each N2S3 actor has a network container capable of containing one or multiple entities. This particularity allows us to control how distribution and parallelism is managed in each simulation. Figure 2 illustrates the architecture with an example network. In this figure, the simulation is made of one input and a neural network organized in two layers. The input of the simulation resides in one actor, the hidden layer is split in one actor per neuron, and all neurons of the output layer reside in a single actor. The reason why one would put several neurons in the same or a different actor is based on the scalability and the synchronization choices for the experiment. The topology of the network (discussed in more details in Section III-C) is one of the major influence on these questions.

The communication between simulation entities requires to identify each object by a URI made of a pair (actor, local_identifier). The component actor is the Akka actor that contains the entity and local_identifier is a path that identifies each object uniquely within its actor.

Since actors are inherently concurrent, one concern is how the temporal order of messages is guaranteed during the simulation. To do so, N2S3 allows to configure several levels of synchronization to be used by the simulation designer. On one end of the spectrum, N2S3 may make use of a unique synchronizer for the simulation, which will ensure that no causality issues happen but can create a bottleneck that will affect the performance of the simulation. On the other end of the spectrum, we can also configure N2S3 to use a synchronization mechanism which is local to each neuron. The latter policy enables better parallelism, but may cause some temporal consistency problems.

A. Neuron Modeling

The neuron is a dynamic element and processing unit that emits output pulses whenever the excitation exceeds some threshold. The resulting sequence of pulses or spikes contains all the information that is transmitted from one neuron to another one. As we have chosen event-driven simulation, the state of a neuron is updated only when an event representing a spike is received.

We provide by default a Leaky-Integrate-and-Fire (LIF) neuron model [START_REF] Maass | Pulsed Neural Networks[END_REF]. There are several reasons for using a LIF model:

• CMOS technology fabrication of this model is available [START_REF] Chicca | A VLSI recurrent network of integrate-and-fire neurons connected by plastic synapses with long-term memory[END_REF], [START_REF] Liu | Temporal coding in a silicon network of integrate-and-fire neurons[END_REF]; • it is fast enough to simulate, and in particular it is effective for large-scale network simulations [START_REF] Brette | Simulation of networks of spiking neurons: A review of tools and strategies[END_REF]. In LIF, neurons integrate the spikes coming from other neurons. These spikes can change the internal potential of the neuron, which is known as the membrane potential or state variable of the neuron. When the membrane potential reaches a given threshold voltage after integrating enough input spikes, an action potential occurs; in other words, the neuron fires, generating an output spike. This is done by updating the membrane potential using the value computed when receiving the last spike and the analytical solution to the differential equations defining the behavior of the LIF neuron.

In order to improve the learning capabilities of neural networks, we provide two refinements of the neuron model: homeostasis and a refractory period. Homeostasis [START_REF] Marder | Variability, compensation and homeostasis in neuron and network function[END_REF] dynamically adapts the threshold of each neuron to the activity of this neuron to prevent a neuron from being over-active or inactive. To allow more specialization of the neurons, we include a refractory period after the firing of a neuron during which it ignores incoming spikes.

B. Synapse modeling and learning

A synapse operates as a plastic controller between two neurons. This plasticity is believed to be the origin of the learning and memorization capabilities of the brain [START_REF] Martin | Synaptic plasticity and memory: an evaluation of the hypothesis[END_REF]. Hence, hardware spiking neural networks usually use some kind of synaptic plasticity to enable learning. In N2S3, we have modeled and simulated hardware synapses and one standard learning behavior, namely Spike Timing-Dependent Plasticity (STDP).

STDP is a local weight modification based on the timing difference between presynaptic spikes and postsynaptic spikes. It increases the connectivity between the neurons when there is a temporal causality between the spikes they emitted. This is implemented by locally remembering the past spikes and updating the synaptic weight according to the chosen STDP rule.

Using a memristor as a nonvolatile synaptic memory has been proposed in previous work [START_REF] Jo | Nanoscale Memristor Device as Synapse in Neuromorphic Systems[END_REF], [START_REF] Alibart | An Organic Nanoparticle Transistor Behaving as a Biological Spiking Synapse[END_REF]- [START_REF] Indiveri | Integration of nanoscale memristor synapses in neuromorphic computing architectures[END_REF]. The basic artificial synapse model that we provide is taken from [START_REF] Querlioz | Immunity to Device Variations in a Spiking Neural Network With Memristive Nanodevices[END_REF], [START_REF] Querlioz | Learning with memristive devices: How should we model their behavior?[END_REF], which describes a nonvolatile memristor suitable for eventdriven and STDP computation. In addition, we designed a new model of synapse which is able to forget unimportant events and remember significant events by combining a nonvolatile memristor and a volatile one [START_REF] Shahsavari | Combining a Volatile and Nonvolatile Memristor in Artificial Synapse to Improve Learning in Spiking Neural Networks[END_REF].

C. Network Topologies

Another purpose of N2S3 is to allow an easy exploration of different neural network configurations. In order to facilitate the network construction, neurons are gathered into groups. Those groups, which usually represent neuron layers, are organized into a specified shape. This information allows the automation of the builder and of the visualization creation.

Inside a group one can specify the use of lateral inhibition to implement the winner-take-all rule [START_REF] Maass | On the Computational Power of Winner-Take-All[END_REF] that states that a spiking neuron inhibits its neighbors during a short time period to allow the specialization of the neurons and hence to enhance unsupervised learning.

Neuron groups can be connected to each other. Connection patterns are managed by different connection policies. By default, N2S3 uses unidirectional full connections (each neuron of the input group is connected to every neuron of the output group). Other policies are also bundled with N2S3, such as the one-to-one policy, which connects each neuron of the input group to only the corresponding neuron of the output group, and a random policy, which creates connections according to a given probability distribution, e.g. to enable Reservoir Computing [START_REF] Verstraeten | An experimental unification of reservoir computing methods[END_REF] modeling.

Several builders have the responsibility to create different topologies. They can use the information of an existing layer (e.g., the shape of the neuron group), if there is one, to construct a new one. Builders bundled with N2S3 include, for instance, builders for the typical layers of Convolutional Neural Networks [START_REF] Lecun | Convolutional networks for images, speech, and time series[END_REF].

IV. N2S3 Features

A. Input Processing

The simulator uses a piped stream system to provide stimuli to the network entities. The input process typically starts by an input reader, which reads data from files or any external source, followed by a number of streams that filter the input data before feeding it to the network.

Input readers provided with N2S3 allow to read data in a variety of format, including standard formats such as Address Event Representation (AER), a data format used by spike-based cameras, or MNIST, used in a standard dataset for handwritten digit recognition (see Section IV-E for details about the data).

Subsequent filter streams available in N2S3 include coding streams, which convert raw numerical data into sequences of spike timings, spike presentation streams (e.g., repeating input spike over a given period, or shuffling spikes), and modifier streams, that alter the input spikes (e.g., by adding noise). Users are free to use one or multiple input readers and to combine any number of filter streams in any order; they may also easily create their own readers and filters.

B. Visualization tools

Users may observe simulation outputs (spikes, weight values…) through network observers. Network observers follow the observer pattern by subscribing to events in the simulation (e.g., when a spike happens), perform some calculations on such events, and make them visible to the user. Examples of such observers range from textual loggers to dynamic visualizations of the spikes of each neuron. Concretely, N2S3 provides a spike activity map of the network, a synaptic weight evolution visualizer, and the calculation of evaluation metrics (e.g. recognition rates, confusion matrices…).

C. Experiment Specification Language

N2S3 includes a dedicated internal Domain Specific Language (DSL) that aims to simplify the creation of simulations. At a higher level of abstraction, users can design experiments (network topology, neuron and synapse properties, observation units…) without having to deal with core features such as synchronization or actor policies. The snippet of code below illustrates (1) the creation of a layer with 18 standard LIF neurons with a 20 mV threshold, (2) its full connection to an existing layer (named outputLayer), and the addition of two observers: (3) a standard weight observer and (4) a custom (i.e., defined by the user) spike observer. The DSL also allows the definition of different stages for the simulation (e.g., splitting the simulation into a training phase and a test phase). At a lower level, the DSL can specify how neurons are organized within actors (e.g., deploying all neurons within the same actor or each neuron in a specific actor). Users may use pre-existing policies or define their own.

D. Software Engineering Practices

To ensure the quality of the simulator's core and its surrounding features, we follow several good practices of sofware engineering. We use automated tests (unit tests and integration tests) to validate the basic behavior of the simulator. These unit tests are run both locally and within a continuous integration environment to detect regressions in the project as soon as possible. N2S3's codebase is open-source, and its developers perform periodic code reviews on it.

E. Standard experiments

N2S3 comes with a set of pre-implemented experiments from the literature. These implementations both demonstrate the features and simulation accuracy of N2S3 and provide code snippets to allow users to understand the ESL and help them designing their own simulations. Specifically, experiments for two standard tasks are provided with N2S3: handwritten digit recognition on the MNIST dataset and car counting on the Freeway dataset. In both cases N2S3 provide simulation results comparable to those reported in the literature.

MNIST [START_REF] Lecun | Gradient-based learning applied to document recognition[END_REF] is a standard dataset for automatic handwritten digit recognition. Since its creation, it has been extensively used to evaluate algorithms for machine learning, computer vision and document recognition. The task is to learn to identify which digit is represented on each image. The dataset includes 60,000 greyscale images of size 128 × 128 pixels, divided into a training set (50,000 images) and a test set (10,000). The implementation of SNN learning on MNIST provided with N2S3 follows the experimental settings from [START_REF] Querlioz | Bioinspired Networks with Nanoscale Memristive Devices that Combine the Unsupervised and Supervised Learning Approaches[END_REF]: 128 × 128 inputs (1/pixel) and one output layer (10 to 300 neurons). Figure 4 shows the result of learning the MNIST database with 30 neurons on the output layer and winner-take-all activated to enhance learning. The initial synaptic weights are random.

Freeway [START_REF]DVS128 Dynamic Vision Sensor Silicon Retina Data[END_REF] is an AER (Address Event Representation) video of a Freeway in Pasadena. It is a standard video demonstrating the capabilities of spike-based cameras. The task here is to count the number of vehicles passing over each of the six lanes of the freeway. The input data is the spikes recorded by a spiking camera of resolution 128 × 128 pixels over a sequence of 78.5 seconds (5.2 millions spikes). The implementation of the Freeway experiment provided with N2S3 reproduces the experimental architecture and settings described in [START_REF] Bichler | Extraction of temporally correlated features from dynamic vision sensors with spike-timing-dependent plasticity[END_REF]: 2 × 128 × 128 inputs, one hidden layer (60 neurones), one output layer (10 neurons), with lateral inhibition on every layer. Figure 5 represents the input of the neural network and figure 6 shows the result of the unsupervised learning process.

V. Conclusion

We have presented in this article a spiking neural network simulator, N2S3. This simulator is dedicated to nanoelectronicsbased hardware neural networks. Its simulation strategy is event-based, for precision and flexibility sakes. As one of the main goals of N2S3 is to promote open research data, i.e. the open distribution of models and experiments for research result reproducibility, it is distributed as open-source software (at https://sourcesup.renater.fr/wiki/n2s3), and a lot of effort has been put in the software architecture, maintainability and enabling features (I/O, graphical outputs, logging, domain specific language). As of version 1.0 (January 2017), N2S3 is distributed with a few neuron and synapse models. It can simulate not only feed-forward networks but also convolutional or even recurrent networks. Finally, it is extensively extensible.

We are planning to periodically release new versions in a timeboxed fashion and to perform more quality measurements in the near future, for example: continuous benchmarking and performance testing to detect non-functional regressions, adding automated code quality and architectural rule enforcement to detect potential bugs, ensuring the respect for architectural and programming idioms, and agilizing the code reviews.

On the functional side, future work includes new models of neurons and synapses, new learning algorithms and procedures, noise modeling, energy modeling, automatic design space exploration, and new applications.

Fig. 1 .

 1 Fig. 1. N2S3 Logo

Fig. 2 .

 2 Fig. 2. N2S3 Architecture. A network is organized in actors that may contain one or more network entities. Such entities could be for example, neurons, inputs or any other.

Fig. 3 .

 3 Fig. 3. N2S3 Packages

Fig. 4 .

 4 Fig. 4. Heat map of the synaptic weights after learning the MNIST data base with 30 neurons on the output layer.

Fig. 5 .

 5 Fig. 5. Input data for the freeway experiment comming from a spike-based camera. The spikes represent a variation of intensity for a given pixel and are generated asynchronously.

Fig. 6 .

 6 Fig. 6. Heatmap showing the reconstruction of the contribution of each input pixel to the activity of the 10 output neurons for the freeway experiment. One can see that some neurons have clearly specialized to detect vehicles on a particular lane.

Acknowledgment

This work has been partly funded by IRCICA (Univ. Lille, CNRS, USR 3380 -IRCICA, F-59000 Lille, France).

The authors would like to thank several colleagues from IEMN for many fruitful discussions about memristor and under the threshold CMOS device and circuit modeling: Fabien Alibart, Alain Cappy, François Danneville, and their teams.