C V Mai

M. D Spiridonakos

E N Chatzi

B Sudret

Surrogate modelling for stochastic dynamical systems by combining NARX models and polynomial chaos expansions

Keywords:

The application of polynomial chaos expansions (PCEs) to the propagation of uncertainties in stochastic dynamical models is well-known to face challenging issues. The accuracy of PCEs degenerates quickly in time. Thus maintaining a sufficient level of long term accuracy requires the use of high-order polynomials. In numerous cases, it is even infeasible to obtain accurate metamodels with regular PCEs due to the fact that PCEs cannot represent the dynamics. To overcome the problem, an original numerical approach was recently proposed that combines PCEs and non-linear autoregressive with exogenous input (NARX) models, which are a universal tool in the field of system identification. The approach relies on using NARX models to mimic the dynamical behaviour of the system and dealing with the uncertainties using PCEs.

The PC-NARX model was built by means of heuristic genetic algorithms. This paper aims at introducing the least angle regression (LAR) technique for computing PC-NARX models, which consists in solving two linear regression problems. The proposed approach is validated with structural mechanics case studies, in which uncertainties arising from both structures and excitations are taken into account. Comparison with Monte Carlo simulation and regular PCEs is also carried out to demonstrate the effectiveness of the proposed approach.

Introduction

Modern engineering and applied sciences have greatly benefited from the rapid increase of available computational power. Indeed, computational models allow for an accurate representation of complex physical phenomena, e.g. fluid and structural dynamics. Physical systems are in nature prone to aleatory uncertainties, such as the natural variability of system properties, excitations and boundary conditions. In order to obtain reliable predictions from numerical simulations, it is of utmost importance to take into account these uncertainties. In this context, uncertainty quantification has gained particular interest in the last decade. The propagation of uncertainties from the input to the output of the model is commonly associated with sampling-based methods, e.g. Monte Carlo simulation, which are robust, however not suitable when only a small number of simulations is affordable or available. Metamodelling techniques can be employed to circumvent this issue. Polynomial chaos expansions (PCEs) [START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF][START_REF] Soize | Physical systems with random uncertainties: chaos representations with arbitrary probability measure[END_REF]) are a powerful metamodelling technique that is widely used in numerous disciplines. However, it is well-known that PCEs exhibit difficulties when being applied to stochastic dynamical models (Wan and Karniadakis, 2006a;[START_REF] Gerritsma | Time-dependent generalized polynomial chaos[END_REF].

Many efforts have been focused on solving this problem. Wan and Karniadakis (2006a,b) proposed multi-element PCEs, which rely on the decomposition of the random input space into sub-domains and the local computation of PCEs on each sub-element. Note that this approach is carried out in an intrusive manner, i.e. the dynamical mechanism is inherently taken into account by the handling of the set of equations describing the system under investigation. [START_REF] Gerritsma | Time-dependent generalized polynomial chaos[END_REF] introduced the concept of time-dependent PCEs, which consists in adding on-the-fly random variables to the set of uncertain parameters, which are the response quantities at specific time instants. From a similar perspective, [START_REF] Luchtenburg | Long-time uncertainty propagation using generalized polynomial chaos and flow map composition[END_REF] proposed to combine PCEs with flow map composition, which reconstructs response time histories from short-term flow maps, each modelled with PCE. A common feature of the approaches mentioned above is that the dynamical mechanism is mimicked by a specialized tool different from PCE. More precisely, flow map composition and iterative enrichment in time of the random space are used to capture the dynamics of the system.

Recently, [START_REF] Spiridonakos | Metamodeling of structural systems through polynomial chaos ARX models[END_REF][START_REF] Billings | Nonlinear system identification: NARMAX methods in the time, frequency, and spatio-temporal domains[END_REF], 2015b,a) introduced a numerical approach that is based on the combination of two techniques, namely PCEs and nonlinear autoregressive with exogenous input (NARX) modelling, which is a universal tool in the field of system identification. In the proposed approach, the NARX model is used to represent the dynamical behaviour of the system, whereas PCEs tackle the uncertainty part. A two phase scheme is employed. First, a stochastic NARX model is identified to represent the dynamical system. It is characterized by a set of specified NARX model terms and associated random coefficients.

Second, the latter are represented as PCEs of the random input parameters which govern the uncertainties in the considered system. In the two phases, both the NARX terms and the polynomial functions are selected by means of the heuristic genetic algorithm, which evolves randomly generated candidates toward better solutions using techniques inspired by natural evolution, e.g. mutation, crossover. The PC-NARX model is distinguished from conventional deterministic system identification tools in that it allows one to account for uncertainties arising both from the system properties, e.g. stiffness, hysteretic behaviour and energy dissipation, and from the stochastic excitations, e.g. ground motions in structural analysis. The approach proved its effectiveness in several case studies in structural dynamics (Spiridonakos and Chatzi, 2015b,a).

It is worth mentioning that early combinations of system identification tools with polynomial chaos expansions can be found in the literature. [START_REF] Ghanem | Identification and prediction of stochastic dynamical systems in a polynomial chaos basis[END_REF] regressed the restoring force of an oscillator on the Chebychev polynomials of state variables of the system, then used PCEs to represent the polynomial coefficients. [START_REF] Wagner | A polynomial chaos approach to ARIMA modeling and terrain characterization[END_REF] used PC-ARIMA models with a-priori known deterministic coefficients for characterizing terrain topology. Linear ARX-PCE models were also used by [START_REF] Kopsaftopoulos | A functional model based statistical time series method for vibration based damage detection, localization, and magnitude estimation[END_REF], [START_REF] Samara | Aircraft virtual sensor design via a time-dependent functional pooling NARX methodology[END_REF] and [START_REF] Sakellariou | Functionally pooled models for the global identification of stochastic systems under different pseudo-static operating conditions[END_REF]. However, in those studies the input parameters are not characterized by known probability density functions, thus the bases are constructed from arbitrarily selected families of orthogonal polynomials. [START_REF] Jiang | A time varying coefficient vector AR modeling of nonstationary covariance time series[END_REF] and [START_REF] Poulimenos | Parametric time-domain methods for non-stationary random vibration modelling and analysis -A critical survey and comparison[END_REF] used time-dependent autoregressive moving average model with stochastic parameter evolution, in which the model parameters are random variables that change with time under stochastic smoothness constraints. Latterly, [START_REF] Kumar | Robust nonlinear MPC based on Volterra series and polynomial chaos expansions[END_REF] represented the coefficients of a Volterra series with PCEs. The above mentioned approaches were carried out in the time domain. In the frequency domain, [START_REF] Pichler | A mode-based meta-model for the frequency response functions of uncertain structural systems[END_REF] used linear regression-based metamodels to represent the frequency response functions of linear structures. In particular, the same strategy was recently used with PCEs by [START_REF] Yang | Application of polynomial chaos expansion and model order reduction for dynamic analysis of structures with uncertainties[END_REF] and [START_REF] Jacquelin | Polynomial chaos expansion and steady-state response of a class of random dynamical systems[END_REF].

In summary, the PC-NARX model consists of two components, namely a NARX and a PCE model. Spiridonakos and Chatzi (2015a,b) computed the two components with the heuristic genetic algorithm. This paper aims at introducing the so-called least angle regression (LARS) technique [START_REF] Efron | Least angle regression[END_REF] for the computation of both NARX and PCE which are merely linear regression models. Indeed LARS has proven to be efficient in computing adaptive sparse PCEs at a relatively low computational cost [START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on Least Angle Regression[END_REF]. LARS has been recently used for selecting the NARX terms in the context of system identification [START_REF] Zhang | Forward and backward least angle regression for nonlinear system identification[END_REF]. Yet the original contribution of this paper is to use LARS as a selection scheme for both NARX terms and PCE basis terms. This way we provide a new fully non-intrusive surrogate modelling technique that allows to tackle nonlinear dynamical systems with uncertain parameters.

The paper is organized as follows: in the first two sections, the theory of PCEs and NARX model are briefly recalled. Section 4 presents the PC-NARX model and the proposed LARSbased approach for computing it. The methodology is illustrated with three benchmark engineering case studies in Section 5. Detailed discussions on the approach are then given before the final conclusions.

Polynomial chaos expansions

Polynomial chaos expansions

Let us consider the computational model Y = M(Ξ) where Ξ = (Ξ 1 , . . . , Ξ M ) is a Mdimensional input vector of random variables with given joint probability density function f Ξ defined over an underlying probability space (Ω, F, P) and M : ξ ∈ D Ξ ⊂ R M → R is the computational model of interest, where D Ξ is the support of the distribution of Ξ. Without loss of generality, we assume that the input random variables are independent. When dependent input variables are to be considered, an iso-probabilistic transform is used to map the latter to independent auxiliary variables [START_REF] Sudret | Polynomial chaos expansions and stochastic finite element methods[END_REF].

The scalar output Y , which is assumed to be a second order random variable, i.e. E Y 2 < +∞, can be represented using generalized polynomial chaos expansion as follows [START_REF] Xiu | The Wiener-Askey polynomial chaos for stochastic differential equations[END_REF][START_REF] Soize | Physical systems with random uncertainties: chaos representations with arbitrary probability measure[END_REF]:

Y = α∈N M y α ψ α (Ξ), (1)
in which ψ α (Ξ) = M i=1 ψ i αi (Ξ i ) are multivariate orthonormal polynomials obtained by the tensor product of univariate polynomials ψ i αi (Ξ i ), y α are associated deterministic coefficients, α = (α 1 , . . . , α M ) is the multi-index vector with α i , i = 1, . . . , M being the the degree of the univariate polynomial ψ i αi (Ξ i ). The latter constitutes a basis of orthonormal polynomials with respect to the marginal probability measure. For instance, when Ξ i is a uniform (resp. standard normal) random variable, the corresponding basis comprises orthonormal Legendre (resp. Hermite) polynomials [START_REF] Abramowitz | Handbook of mathematical functions[END_REF].

In practice, it is not tractable to use an infinite series expansion. An approximate representation is obtained by means of a truncation:

Y = α∈A y α ψ α (Ξ) + ≡ M P C (Ξ) + , (2)
in which A is the truncation set and is the truncation-induced error. [START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on Least Angle Regression[END_REF] introduced a hyperbolic truncation scheme, which consists in selecting all polynomials satisfying the following criterion:

A M,p q =    α ∈ N M : α q def = M i=1 α q i 1/q ≤ p    , (3)
with p being the highest total polynomial degree, 0 < q ≤ 1 being the parameter determining the hyperbolic truncation surface. To further reduce the number of candidate polynomials, one can additionally apply a low-rank truncation scheme which reads [START_REF] Blatman | Adaptive sparse polynomial chaos expansions for uncertainty propagation and sensitivity analysis[END_REF]:

A M,p q,r = α ∈ N M : α 0 def = M i=1 1 αi>0 ≤ r, α q ≤ p , (4)
where α 0 is the rank of the multivariate polynomial ψ α , defined as the total number of nonzero components α i , i = 1, . . . , M . The prescribed rank r is usually chosen as a small integer value, e.g. r = 2, 3.

Computing the coefficients by least square minimization

The computation of the coefficients in Eq. ( 2) can be conducted by means of intrusive approach (i.e. Galerkin scheme) or non-intrusive approaches (e.g. stochastic collocation, projection, regression methods) [START_REF] Sudret | Uncertainty propagation and sensitivity analysis in mechanical models -Contributions to structural reliability and stochastic spectral methods[END_REF][START_REF] Xiu | Numerical methods for stochastic computations -A spectral method approach[END_REF]. The regression method consists in estimating the set of coefficients ŷα = {y α , α ∈ A} that minimizes the mean square error:

E 2 def = E   Y - α∈A y α ψ α (Ξ) 2   , (5)
which means: ŷα = arg min

y α ∈R cardA E   M(Ξ) - α∈A y α ψ α (Ξ) 2   . (6)
In practice, the coefficients are obtained by minimizing an empirical mean over a sample set: ŷα = arg min

y α ∈R cardA 1 N N i=1 M(ξ (i) ) - α∈A y α ψ α (ξ (i) ) 2 , ( 7
)
where X = ξ (i) , i = 1, . . . , N is an experimental design (ED) obtained with random sampling of the input random vector. The computational model M is run for each point of the ED, yielding the vector of output values Y = y (i) = M(ξ (i) ), i = 1, . . . , N . By evaluating the polynomial basis onto each sample point in the ED, one obtains the information matrix, which is defined as follows:

A def = A ij = ψ j (ξ (i) ), i = 1, . . . , N, j = 1, . . . , cardA , (8)
i.e. the i th row of A is the evaluation of the polynomial basis functions at the point ξ (i) in the ED. Eq. ( 7) basically represents the problem of estimating the parameters of a linear regression model, for which the least squares solution reads:

ŷα = A T A -1 A T Y. (9)

Error estimator

Herein, the accuracy of the truncated expansion (Eq. ( 2)) is estimated by means of the leaveone-out (LOO) cross validation technique, which allows a fair error estimation at an affordable computational cost [START_REF] Blatman | An adaptive algorithm to build up sparse polynomial chaos expansions for stochastic finite element analysis[END_REF]. Assume that one is given an experimental design (ED) X = ξ (i) , i = 1, . . . , N and the associated vector of model responses Y = y (i) , i = 1, . . . , N obtained by running the computational model M. Cross-validation consists in partitioning the ED into two complementary subsets, training a model using one subset, then validating its prediction on the other subset. In this context, the term LOO means that the validation set comprises only one sample and the model is computed using all the N -1 remaining samples. A single point of the experimental design is left out at a time, a PCE M PC\i is built using the remaining points and the corresponding prediction error on the validation point is computed:

∆ i def = M(ξ (i) ) -M PC\i (ξ (i) ). ( 10
)
The LOO error is defined by:

Err LOO = 1 N N i=1 ∆ 2 i = 1 N N i=1 M(ξ (i) ) -M PC\i (ξ (i) ) 2 . ( 11
)
At first glance, computing the LOO error may appear computationally demanding since it requires to train and validate N different PCE models. However, by means of algebraic derivations, one can compute the LOO error Err LOO from a single PCE model M PC (•) built with the full ED X as follows [START_REF] Blatman | Adaptive sparse polynomial chaos expansions for uncertainty propagation and sensitivity analysis[END_REF]:

Err LOO = 1 N N i=1 M(ξ (i) ) -M PC (ξ (i) ) 1 -h i 2 , ( 12
)
in which h i is the i th diagonal term of the matrix A A T A -1

A T .

Adaptive sparse PCEs based on least angle regression

To obtain satisfactory results with the least-square minimization method, the number of model evaluations N is commonly required to be two to three times the cardinality P = cardA of the polynomial chaos basis [START_REF] Blatman | Adaptive sparse polynomial chaos expansions for uncertainty propagation and sensitivity analysis[END_REF]. In case when the problem involves high dimensionality and when high order polynomials must be used, i.e. P is large, the required size of the ED soon becomes excessive. An effective scheme that allows one to obtain accurate polynomial expansions with limited number of model evaluations is of utmost importance.

In the following, we will shortly describe the adaptive sparse PCE scheme proposed by [START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on Least Angle Regression[END_REF] which is a non-intrusive least-square minimization technique based on the least angle regression (LARS) algorithm [START_REF] Efron | Least angle regression[END_REF]. LARS is an efficient numerical technique for variable selection in high dimensional problems. It aims at detecting the most relevant predictors among a large set of candidates given a limited number of observations. In the context of PCEs, LARS allows one to achieve a sparse representation, i.e. the number of retained polynomial functions is small compared to the size of the candidate set. This is done in an adaptive manner, i.e. the candidate functions become active one after another in the descending order of their importance. The relevance of a candidate polynomial is measured by means of its correlation with the current residual of the expansion obtained in the previous iteration. The optimal PCE is chosen so as to minimize the leave-one-out error estimator in Eq. ( 12). In particular, LARS generally requires a relatively small ED, thus being effective in realistic applications. The reader is referred to [START_REF] Efron | Least angle regression[END_REF] for more details on the LARS technique and to [START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on Least Angle Regression[END_REF] for its implementation in the adaptive sparse PCE scheme.

Time-frozen PCEs and associated statistics

In this paper, we focus on stochastic dynamical systems where the random output response is a time dependent quantity Y (t) = M(Ξ, t). In this context, the polynomial chaos representation of the response may be cast as: In such an approach, the metamodel of the response is computed independently at each time instant, hence the name time-frozen PCEs. This naive approach is presented here for the sake of comparison with the PC-NARX method introduced later.

Y (t) = α∈A y α (t)ψ α (Ξ) + (t), (13)
Eq. ( 13) can be used to compute the evolution of the response statistics. The multivariate polynomial chaos functions are orthonormal, i.e. :

E ψ α (Ξ) ψ β (Ξ) def = DΞ ψ α (ξ) ψ β (ξ) f Ξ (ξ) dξ = δ αβ ∀α, β ∈ N M , (14)
in which δ αβ is the Kronecker symbol that is equal to 1 if α = β and equal to 0 otherwise. In particular, each multivariate polynomial is orthonormal to

ψ 0 (Ξ) = 1, which means E [ψ α (Ξ)] = 0 ∀α = 0 and Var [ψ α (Ξ)] = E ψ 2 α (Ξ) = 1 ∀α = 0.
Thus, the time dependent mean and standard deviation of the response can be estimated with no additional cost by means of the post-processing of the truncated PC coefficients in Eq. ( 13) as follows:

μY (t) def = E α∈A y α (t)ψ α (Ξ) = y 0 (t) , (15)
σ2

Y (t) def = Var α∈A y α (t)ψ α (Ξ) = α∈A α =0 y 2 α (t). ( 16
)
3 Nonlinear autoregressive with exogenous input model

Let us consider a computational model y(t) = M(x(t)) where x(t) is the time-dependent input excitation and y(t) is the response time history of interest. System identification aims at building a mathematical model describing M using the observed data of the input and output signals.

In this paper, we focus on system identification in the time domain. One discretizes the time duration under investigation in T discrete instants t = 1, . . . , T . A nonlinear autoregressive with exogenous input (NARX) model allows one to represent the output quantity at a considered time instant as a function of its past values and values of the input excitation at the current or previous instants [START_REF] Chen | Orthogonal least squares methods and their application to non-linear system identification[END_REF][START_REF] Billings | Nonlinear system identification: NARMAX methods in the time, frequency, and spatio-temporal domains[END_REF]:

y(t) = F (z(t)) + (t) = F (x(t), . . . , x(t -n x ), y(t -1), . . . , y(t -n y )) + t , (17)
where F(•) is the underlying mathematical model to be identified, z(t) = (x(t), . . . , x(tn x ), y(t -1), . . . , y(tn y )) T is the vector of current and past values, n x and n y denote the maximum input and output time lags, t ∼ N (0, σ 2 (t)) is the residual error of the NARX model.

In standard NARX models, the residuals are assumed to be independent normal variables with zero mean and variance σ 2 (t). There are multiple options for the mapping function F(•). In the literature, the following linear-in-the-parameters form is commonly used:

y(t) = ng i=1 ϑ i g i (z(t)) + t , (18)
in which n g is the number of model terms g i (z(t)) that are functions of the regression vector z(t) and ϑ i are the coefficients of the NARX model.

Indeed, a NARX model allows ones to capture the dynamical behaviour of the system which follows the principle of causality, i.e. the current output quantity (or state of the system) y(t) is affected by its previous states {y(t -1), . . . , y(tn y )} and the external excitation {x(t), . . . , x(tn x )}. Note that the cause-consequence effect tends to fade away as time evolves, therefore it suffices to consider only a limited number of time lags before the current time instant. It is worth emphasizing that the model terms may be constructed from a variety of global or local basis functions. For instance the use of polynomial NARX model with g i (z(t)) being polynomial functions is relatively popular in the literature.

The identification of a NARX model for a system consists of two major steps. The first one is structure selection, i.e. determining which NARX terms g i (z(t)) are in the model. The second step is parameter estimation, i.e. determining the associated model coefficients. Note that structure selection, particularly for systems involving nonlinearities, is critically important and difficult. Including spurious terms in the model leads to numerical and computational problems [START_REF] Billings | Nonlinear system identification: NARMAX methods in the time, frequency, and spatio-temporal domains[END_REF]. Billings suggests to identify the simplest model to represent the underlying dynamics of the system, which can be achieved by using the orthogonal least squares algorithm and its derivatives to select the relevant model terms one at a time [START_REF] Billings | Nonlinear system identification: NARMAX methods in the time, frequency, and spatio-temporal domains[END_REF]. Different approaches for structure selection include trial and error methods, see e.g. [START_REF] Chen | On-board identification and control performance verification of an MR damper incorporated with structure[END_REF]; [START_REF] Piroddi | Simulation error minimisation methods for NARX model identification[END_REF], and correlation-based methods, see e.g. [START_REF] Billings | An adaptive orthogonal search algorithm for model subset selection and non-linear system identification[END_REF]; [START_REF] Wei | Model structure selection using an integrated forward orthogonal search algorithm assisted by squared correlation and mutual information[END_REF]. There is a rich literature dedicated to this topic, however discussions on those works are not in the scope of the current paper.

The identified model can be used for several purposes. First, it helps the analysts reveal the mechanism and behaviour of the underlying system. Understanding how a system operates offers one the possibility to control it better. Second, the identified mathematical model can be utilized for predicting future responses of the system. From this point of view, it can be considered a metamodel (or approximate model) of the original M.

Polynomial chaos -nonlinear autoregressive with exogenous input model

Consider a computational model y(t, ξ) = M(x(t, ξ x ), ξ s ) where ξ = (ξ x , ξ s ) T is the vector of uncertain parameters, ξ x and ξ s respectively represent the uncertainties in the input excitation

x(t, ξ x ) and in the system itself. For instance, ξ x can contain parameters governing the amplitude and frequency content of the excitation time series, while ξ s can comprise parameters determining the system properties such as geometries, stiffness, damping and hysteretic behaviour. Spiridonakos and Chatzi (2015b,a) proposed a numerical approach based on PCEs and NARX model to identify the metamodel of such a dynamical system with uncertainties arising from both the excitation and the system properties. The time-dependent output quantity is first represented by means of a NARX model:

y(t, ξ) = ng i=1 ϑ i (ξ) g i (z(t)) + g (t, ξ), (19)
in which the model terms g i (z(t)) are functions of the regression vector z(t) = (x(t), . . . , x(t - represented by means of truncated PCEs as follows [START_REF] Soize | Physical systems with random uncertainties: chaos representations with arbitrary probability measure[END_REF]:

n x ), y(t -1), . . . , y(t -n y )) T ,
ϑ i (ξ) = n ψ j=1 ϑ i,j ψ j (ξ) + i , (20)
in which {ψ j (ξ), j = 1, . . . , n ψ } are multivariate orthonormal polynomials of ξ, {ϑ i,j , i = 1, . . . , n g , j = 1, . . . , n ψ } are associated PC coefficients and i is the truncation error. Finally, the PC-NARX model reads:

y(t, ξ) = ng i=1 n ψ j=1 ϑ i,j ψ j (ξ) g i (z(t)) + (t, ξ), (21)
where (t, ξ) is the total error time series due to the truncations of NARX and PCE models.

In the proposed approach, the NARX model is used to capture the dynamics of the considered system, whereas PCEs are used to propagate uncertainties.

Let us discuss the difference between the PC-NARX model and the conventional timedependent PCE formulation in Eq. ( 13). For the sake of clarity, Eq. ( 21) can be rewritten as follows:

y(t, ξ) = n ψ j=1 ng i=1 ϑ i,j g i (z(t)) ψ j (ξ) + (t, ξ). ( 22
)
At a considered instant t, the polynomial coefficients y j (t) def = ng i=1 ϑ i,j g i (z(t)) are represented as functions of the past values of the excitation and the output quantity of interest. Consequently the polynomial coefficients follow certain dynamical behaviours. In the conventional model in Eq. ( 13), the polynomial coefficients at time t are deterministic. Therefore high and increasing polynomial order is required to maintain an accuracy level and properly capture the dynamics as time evolves (Wan and Karniadakis, 2006a). In contrast, when a functional form is used to relate the coefficients y j (t) with the excitation and output time series, constant and low polynomial order suffices. Spiridonakos and Chatzi (2015a) used PC-NARX models with fourth order PCEs to obtain remarkable results in the considered structural dynamics case studies. In the literature, [START_REF] Gerritsma | Time-dependent generalized polynomial chaos[END_REF] showed that when applying time-dependent PCEs, i.e. adding previous responses to the set of random variables to represent current response, low-order polynomials could also be used effectively. From a similar perspective, the PC-flow map composition scheme proposed by [START_REF] Luchtenburg | Long-time uncertainty propagation using generalized polynomial chaos and flow map composition[END_REF] was also proven efficient in solving the problems with low polynomial order, which was impossible with PCEs alone.

Indeed, not all the NARX and PC terms originally specified are relevant, as commonly observed in practice. The use of redundant NARX or PC terms might lead to large inaccuracy.

Therefore, it is of utmost importance to identify the correct structure of NARX and PC models,

i.e. to select appropriate NARX terms and PC bases. To this end, Spiridonakos and Chatzi (2015a) proposed a two-phase approach, in which the NARX terms and PC functions are subsequently selected by means of the genetic algorithm. However, due to the linear-in-parameters formulations of the NARX model (Eq. ( 19)) and the PC expansions (Eq. ( 20)), the question of selecting NARX and PC terms boils down to solving two linear regression problems. To this end, it appears that one can use techniques that are specially designed for linear regression analysis, for instance least angle regression (LARS) [START_REF] Efron | Least angle regression[END_REF]. LARS has been recently used for selecting NARX terms [START_REF] Zhang | Forward and backward least angle regression for nonlinear system identification[END_REF]. The use of LARS in the field of system identification can be classified as a correlation-based method, which selects the NARX terms that make significant contribution to the output using correlation analysis, see e.g. [START_REF] Billings | An adaptive orthogonal search algorithm for model subset selection and non-linear system identification[END_REF]; [START_REF] Wei | Model structure selection using an integrated forward orthogonal search algorithm assisted by squared correlation and mutual information[END_REF]. LARS has also been used in the adaptive sparse PCE scheme and proven great advantages compared to the other predictor selection methods, i.e.

fast convergence and high accuracy with an ED of limited size.

Least angle regression-based approach

In this section, we introduce least angle regression (LARS) for the selection of appropriate NARX and PCE models. A two phase approach is used, which sequentially selects NARX and PCE models as follows:

• Phase 1: Selection of the appropriate NARX model among a set of candidates.

-Step 1.1: One specifies general options for the NARX model (model class and related properties), e.g. type of basis functions (polynomials, neural network, sigmoid functions, etc. ), maximum time lags of input and output, properties of the basis functions (e.g. maximum polynomial order). Note that it is always preferable to start with simple models having a reasonable number of terms. In addition, any available knowledge on the system, e.g. number of degrees of freedom, type of non-linear behaviour, should be used in order to obtain useful options for the general NARX structure. This leads to a full NARX model which usually contains more terms than actually needed for a proper representation of the considered dynamical system. At this stage, one assumes that the specified full NARX model contains the terms that can sufficiently describe the system. This assumption will be verified in the final step of this phase.

- For experiment #k, the one-step-ahead prediction of the response reads:

ŷp (t, ξ k ) = ng i=1 ϑ i (ξ k ) g i (ẑ p (t)), (23)
in which

ẑp (t) = (x(t), . . . , x(t -n x ), y(t -1), . . . , y(t -n y )) T . (24)
It is worth emphasizing that y(t) is the recorded data used for training the NARX model. Denoting φ(t) = {g i (ẑ p (t), i = 1, . . . , n g )} T and ϑ = {ϑ i , i = 1, . . . , n g } T , the residual time series reads:

p (t, ξ k ) = y(t, ξ k ) -ŷp (t, ξ k ) = y(t, ξ k ) -φ T (t) ϑ(ξ k ). (25)
Thus, the sum of squared errors is given by:

T t=1 [ p (t, ξ k )] 2 = T t=1 y(t, ξ k ) -φ T (t) ϑ(ξ k ) 2 . ( 26
)
Using Eq. ( 25) and assembling all time instants in the k-th experiment, one obtains:

     y(1, ξ k ) . . . y(T, ξ k )      =      φ T (1) . . . φ T (T )      ϑ(ξ k ) +      p (1, ξ k ) . . . p (T, ξ k )      (27)
The above equation can be rewritten in matrix notations as follows:

y k = Φ k ϑ(ξ k ) + k , (28)
where y k is the T × 1 vector of output time-series, Φ k is the T × n g information matrix with the i th row containing the evaluations of NARX terms φ(t) at instant t = i and k is the residual vector. This is typically the equation of a linear regression problem, for which the relevant NARX regressors among the NARX candidate terms φ(t) can be selected by LARS.

Note that the same candidate model might be obtained from different selected experiments. In theory all the available experiments can be considered, i.e. the number of candidate models is at most the size of the ED. Herein, we search for the appropriate NARX model among a limited number of experiments which are exhibiting strong non-linearity.

-Step 1.3: For each candidate NARX model, the corresponding NARX coefficients are computed for each of the experiments by ordinary least-squares. The parameters ϑ(ξ k ) minimizing the total errors in Eq. ( 26) is the least-squares solution of Eq. ( 28), i.e. :

ϑ(ξ k ) = arg min ϑ ( T k k ) = Φ T k Φ k -1 Φ T k y k , (29)
with the information matrix Φ k containing only the NARX regressors specified in the NARX candidate. Having at hand the NARX coefficients, the free-run reconstruction for each experiment output is conducted as follows:

ŷs (t, ξ k ) = ng i=1 ϑ i (ξ k ) g i (ẑ s (t)), (30)
in which ẑs (t) = (x(t), . . . , x(tn x ), ŷs (t -1), . . . , ŷs (t

-n y )) T . (31)
It is worth underlining that the free-run reconstruction of the response is obtained using only the excitation time series x(t) and the response initial condition y 0 . The response is reconstructed recursively, i.e. its estimate at one instant is used to predict the response at later instants. This differs from Eq. ( 23) where the recorded response was used in the recursive formulation. The relative error for simulation #k reads:

k = T t=1 (y(t, ξ k ) -ŷs (t, ξ k )) 2 T t=1 (y(t, ξ k ) -ȳ(t, ξ k )) 2 , ( 32
)
in which ŷs (t, ξ k ) is the output trajectory reconstructed by the NARX model and ȳ(t, ξ k ) is the mean value of the response time series y(t, ξ k ).

-Step 1.4: One selects the most appropriate NARX model among the candidates.

Herein, the error criterion of interest is the mean value of the relative errors:

¯ = 1 K K k=1 k . ( 33
)
We propose to choose the model that achieves a sufficiently small overall error on the conducted experiments, e.g. ¯ < 1 × 10 -3 , with the smallest number of NARX terms.

In other words, the appropriate model is the simplest one that allows to capture the system dynamical behaviour, thus following the principle suggested by [START_REF] Billings | Nonlinear system identification: NARMAX methods in the time, frequency, and spatio-temporal domains[END_REF].

To refine the estimated coefficients, a nonlinear optimization for minimizing the simulation error (Eq. ( 32)) may be conducted afterwards Spiridonakos and Chatzi (2015a).

However, this is not used in the current paper due to the fact that LARS allows one to detect the appropriate NARX terms, therefore the models estimated by ordinary least-squares appear sufficiently accurate.

If an appropriate NARX model is not obtained, i.e. the initial assumption that the full NARX model includes an appropriate candidate is not satisfied, the process is re-started from Step 1.1 (choice of model class), when different options for the full NARX model should be considered. For instance, one may use more complex models with larger time lags, different basis functions, etc.

• Phase 2: Representation of the NARX coefficients by means of PCEs using the sparse adaptive PCE scheme which is based on LARS (see section 2.4). The NARX coefficients obtained from Phase 1 are used for training the PC expansion.

For the sake of clarity, the above procedure for computing a PC-NARX model is summarized by the flowchart in Figure 1.

Use of the surrogate model for prediction

The PC-NARX surrogate model can be used for the prediction1 of the response to a set of input parameters ξ . Given the excitation x(t, ξ x ) and the initial conditions of the response y(t = 1, ξ ) = y 0 , the output time history of the system can be recursively obtained as follows: Step 1.1

ŷ(t, ξ ) = ng i=1 n ψ j=1 ϑ i,j ψ j (ξ ) g i (ẑ(t, ξ )), t = 2, . . . , T, (34)
Step 1.2

Step 1.3

Step

ẑ(t, ξ ) = x(t, ξ x ), . . . , x(t -n x , ξ x ), ŷ(t -1, ξ ), . . . , ŷ(t -n y , ξ ) T . (35)
Currently, no close-form formulation for computing the time-dependent statistics of the output quantity is available as opposed to time-frozen PCEs. However, the evolutionary response statistics can be straightforwardly obtained by means of Monte Carlo simulation using the PC-NARX model.

Validation of the surrogate model

The PC-NARX model is computed using an ED of limited size. The validation process is conducted with a validation set of large size which is independent of the ED. A large number, e.g. n val = 10 4 , of input parameters and excitations is generated. One uses the numerical solver to obtain the response time histories sampled at the discrete time instants t = 1, . . . , T . Then PC-NARX model (Eq. ( 34)) is used to predict the time dependent responses to the excitations and uncertain parameters of the validation set. The accuracy of the computed PC-NARX model is validated by means of comparing its predictions with the actual responses in terms of the relative errors and the evolutionary statistics of the response. For prediction #i, the relative error reads:

val,i = T t=1 (y(t, ξ i ) -ŷ(t, ξ i )) 2 T t=1 (y(t, ξ i ) -ȳ(t, ξ i )) 2 , ( 36
)
where ŷ(t, ξ i ) is the output trajectory predicted by PC-NARX and ȳ(t, ξ i ) is the mean value of the actual response time series y(t, ξ i ). The above formula is also used to calculate the accuracy of the time dependent statistics (i.e. mean, standard deviation) predicted by PC-NARX. The mean value of the relative errors over n val predictions reads:

¯ val = 1 n val n val i=1 val,i . (37)
The relative error for a quantity y, e.g. the maximal value of the response (resp. the response at a specified instant) is given by:

val,y = n val i=1 (y i -ŷi ) 2 n val i=1 (y i -ȳ) 2 , ( 38
)
where y i is the actual response, ŷi is the prediction by PC-NARX and ȳ is the mean value defined by ȳ = 1 n val n val i=1 y i .

Numerical applications

The use of LARS-based PC-NARX model is now illustrated with three nonlinear dynamical systems with increasing complexity, namely a quarter car model subject to a stochastic sinusoidal road profile, a single degree-of-freedom (SDOF) Duffing and a SDOF Bouc-Wen oscillator subject to stochastic non-stationary excitation. In all considered numerical examples, uncertainties arising from the system properties and from the excitation are taken into account. PC-NARX models are computed using a small number of numerical simulations as experimental design.

The validation is conducted by comparing their response predictions with the reference values obtained by using Monte Carlo simulation (MCS) on the numerical solvers.

Quarter car model

In the first numerical example, a quarter car model of a vehicle suspension represented by a nonlinear two DOF system [START_REF] Kewlani | A polynomial chaos approach to the analysis of vehicle dynamics under uncertainty[END_REF] (Figure 2) is considered. The displacements of the masses are governed by the following system of ordinary differential equations (ODEs):

   m s ẍ1 (t) = -k s (x 1 (t) -x 2 (t)) 3 -c ( ẋ1 (t) -ẋ2 (t)) m u ẍ2 (t) = k s (x 1 (t) -x 2 (t)) 3 + c ( ẋ1 (t) -ẋ2 (t)) + k u (z(t) -x 2 ) (39)
in which the sprung mass m s and the unsprung mass m u are connected by a nonlinear spring of stiffness k s and a linear damper of damping coefficient c. The forcing function z(t) is applied to m u through a linear spring of stiffness k u . x 1 (t) and x 2 (t) are the displacements of m s and m u respectively. A sinusoidal function road profile with amplitude A and frequency ω is considered:

z(t) = A sin(ω t). (40)
The parameters of the quarter car model and of the excitation ξ = {k s , k u , m s , m u , c, A, ω} are modelled by independent random variables with marginal distributions given in Table 1. Note that the mean value of the parameters are the deterministic values specified in [START_REF] Kewlani | A polynomial chaos approach to the analysis of vehicle dynamics under uncertainty[END_REF]. In addition, Gaussian distributions are used as in [START_REF] Kewlani | A polynomial chaos approach to the analysis of vehicle dynamics under uncertainty[END_REF], although it would be more appropriate to use e.g. lognormal variables to ensure the positivity of mass and stiffness parameters. [START_REF] Kewlani | A polynomial chaos approach to the analysis of vehicle dynamics under uncertainty[END_REF] addressed this numerical example with the multi-element PCE approach.

Figure 2: Quarter car model

We now aim at building the metamodel for representing the displacement time histories x 1 (t)

of the sprung mass m s . For this purpose, N = 100 analyses of the system are conducted with 100 samples of the uncertain input parameters generated by Latin hypercube sampling. The system of ODEs are solved by means of the Matlab solver ode45 (explicit Runge-Kutta method with relative error tolerance 1 × 10 -3 ) for the total duration T = 30 s and the time step dt = 0.01 s.

In the first place, a NARX model structure is chosen, in which the model terms are polynomial functions of past values of the output and excitation

g i (t) = x l 1 (t -j) z m (t -k) with l + m ≤ 3, 0 ≤ l ≤ 3, 0 ≤ m ≤ 1, j = 1, . . . , 4, k = 0, . . . , 4.
The specified full NARX model contains 86 terms. It is worth noticing that the initial choice of the NARX structure was facilitated by the knowledge of the dynamical nonlinear behaviour of the system of interest. For instance, polynomial functions of order up to 3 are used because of the cubic nonlinear behaviour as in Eq. ( 39). As a rule of thumb, the maximum time lags n x = n y = 4 are chosen equal to twice the number of degrees of freedom of the considered system.

Next, the candidate NARX models were computed. To this end, we selected the simulations with maximum displacement exceeding a large threshold, i.e. max(|x 1 (t)|) > 1.2 m and retained 15 experiments. For each selected simulation, LARS was applied to the initial full NARX model to detect the most relevant NARX terms constituting a candidate NARX model. This procedure resulted in 10 different candidate NARX models in total.

For each candidate NARX model, we computed the NARX coefficients for all simulations by ordinary least-squares to minimize the sum of squared errors (Eq. ( 29)). We then reconstructed all the output time histories with the computed coefficients and calculated the relative errors k of the reconstruction (Eq. ( 32)).

Among the 10 candidates, we selected the most appropriate NARX model which achieves a sufficiently small overall error with the smallest number of terms. The selected model results in a mean relative error ¯ = 3.56 × 10 -4 for 100 simulations in the ED and contains 6 terms, namely the constant term, z(t -4), x 1 (t -4), x 1 (t -1), x 3 1 (t -1), x 2 1 (t -4) z(t -4). LARS proves effective in selecting the appropriate NARX model by retaining only 6 among the 86 candidate terms available to describe the system.

In the next step, we expanded the 6 NARX coefficients by adaptive sparse PCEs of order p ≤ 20 with maximum interaction rank r = 2 and truncation parameter q = 1. The NARX coefficients computed in the previous step were used for training the metamodel. PCE models that minimize the LOO errors were selected. This led to LOO errors smaller than 10 -7 . The optimal PCE order selected by the adaptive scheme is up to 6.

For the sake of comparison, we represent the response x 1 (t) by means of time-frozen PCEs.

For this purpose, adaptive sparse PCEs are used with an ED of size N = 500. The best PCE with maximum degree 1 p 20, maximum interaction rank r = 2 and truncation parameter q = 1 is selected. The PC-NARX and time-frozen PCEs are used to predict the output time histories for an independent validation set of size n val = 10 4 which is pre-computed by the numerical solver. The accuracy of the two PCE approaches are compared in the following.

We compare the predictions of the mass displacement at an early time instant t = 5 s and a late instant t = 30 s (Figure 3). The instantaneous responses x 1 (t) predicted by two PCE approaches are plotted in Figure 3(A) and Figure 3(C) versus the actual responses obtained by numerically solving the considered system of ODEs. One observes that PC-NARX outperforms time-frozen PCEs. At t = 30 s, PC-NARX is still capable of providing highly accurate predictions ( val = 4.21 × 10 -3 ) whereas time-frozen PCEs show large inaccuracies already from t = 5 s. In terms of statistics, the probability density functions (PDFs) of the instantaneous responses are depicted in Figure 3(B) and Figure 3(D). The reference PDF is obtained from MCS using the validation set of size n val = 10 4 . One notices that the instantaneous PDFs computed by time-frozen PCEs do not differ significantly from the reference function, whereas the predictions are actually of poor quality. Indeed, it is always possible that two sets of samples provide identical PDFs while their pair-wise comparison shows large discrepancies. This example shows a problem that is overlooked in the literature on surrogate models for dynamical systems, when conclusions are commonly drawn only based on the statistics (for instance, PDFs, mean and standard deviation) of the predicted quantities. This might be misleading for judging the accuracy of the metamodel under investigation in predicting specific output values. The PDFs computed by PC-NARX show perfect agreement with the reference functions, which is obvious because every single PC-NARX-based prediction is highly accurate.

We investigate now the effect of the sample size on the MCS estimates. PC-NARX is used to predict the responses for a different validation set of size n val = 10 6 . Figure 3(B) and Figure 3(D) show that the PDFs obtained with 10 6 runs, which are closer to the true PDFs, differ slightly from the PDFs obtained with 10 4 runs. For instance, the peaks of the 10 6 runs-based PDFs become marginally higher. It is believed that the PDF obtained with 10 6 runs of the PC-NARX surrogate is the most accurate, although this cannot be validated by running 10 6 times the computationally expensive model. In general, it is computational expensive to obtain the true PDFs with MCS on the original numerical model, however it is feasible using PC-NARX model. Around t = 5 s, time-frozen PCEs start showing signs of instability and the predictions become inaccurate. In contrast, PC-NARX provides predictions that are indistinguishable from the actual responses. Over 10 4 validation trajectories, PC-NARX leads to a mean relative error ¯ val = 0.17 × 10 -2 , and only 5 simulations among them exhibit a relative error val,i exceeding 0.1.

The two PCE approaches are now compared in terms of the evolutionary response statistics.

Figure 6 shows that time-frozen PCEs can represent relatively well the mean trajectory, except for the late instants (t > 20 s) when the discrepancy is noticeable. However, time-frozen PCEs fail to capture the oscillatory behaviour of the standard deviation when the prediction starts deviating significantly from the actual trajectory at t = 5 s. The improvement in accuracy of PC-NARX is outstanding, in particular because it can represent in detail the oscillatory response statistics. The relative errors are val,M ean = 0.73 × 10 -2 and val,Std = 0.91 × 10 -2 for the mean and standard deviation time-histories respectively. This numerical case study illustrates the effectiveness of the LARS-based PC-NARX model with respect to the usual PCE approach. It is worth noting that the process was conducted using the uncertainty quantification software UQLab [START_REF] Marelli | UQLab: A framework for uncertainty quantification in Matlab[END_REF], more specifically the polynomial chaos expansion toolbox [START_REF] Marelli | UQLab user manual -Polynomial chaos expansions[END_REF]. The next case studies will be (2082.94, 2187.10, 15.46, 40.08, 564.26, 0.10, 6.41)

Duffing oscillator subject to synthetic ground motions

Let us consider a SDOF Duffing oscillator [START_REF] Kafali | Seismic fragility analysis: Application to simple linear and nonlinear systems[END_REF] subject to stochastic excitation. Note that the vibration of a Duffing oscillator is a commonly used benchmark in the literature, see e.g. Spiridonakos and Chatzi (2015a); [START_REF] Orabi | Nonstationary response analysis of a Duffing oscillator by the Wiener-Hermite expansion method[END_REF]; [START_REF] Lucor | Adaptive generalized polynomial chaos for nonlinear random oscillators[END_REF]; [START_REF] Kougioumtzoglou | An approximate approach for nonlinear system response determination under evolutionary stochastic excitation[END_REF]; [START_REF] Chouvion | Development of error criteria for adaptive multi-element polynomial chaos approaches[END_REF] in which an oscillator subject to either periodic or white-noise external force was considered. The dynamics of the oscillator can be described by the following equation of motion:

ÿ(t) + 2 ζ ω ẏ(t) + ω 2 (y(t) + y(t) 3 ) = -x(t), (41)
in which y(t) is the oscillator displacement, ζ is the damping ratio, ω is the fundamental frequency, is the parameter governing the nonlinear behaviour and x(t) is the excitation.

Herein, the excitation is generated by the probabilistic ground motion model proposed by [START_REF] Rezaeian | Simulation of synthetic ground motions for specified earthquake and site characteristics[END_REF]. The ground motion acceleration is represented as a nonstationary process by means of a modulated filtered white noise process as follows:

x(t) = q(t, α)    1 σ h (t) t -∞ h[t -τ, λ(τ )]ω(τ ) dτ    . (42)
The white-noise process denoted by ω(τ ) passes a filter h[tτ, λ(τ )] which is selected as an impulse-response function:

h[t -τ, λ(τ )] = ω f (τ ) 1 -ζ 2 f (t) exp[-ζ f (τ )ω f (τ )(t -τ )] × sin[ω f (τ ) 1 -ζ 2 f (τ )(t -τ )] for τ ≤ t, h[t -τ, λ(τ )] = 0 for τ > t, (43)
where λ(τ ) = (ω f (τ ), ζ f (τ )) is the vector of time-varying parameters of the filter h. ω f (τ ) and ζ f (τ ) are respectively the filter's frequency and bandwidth at time instant τ . They represent the evolving predominant frequency and bandwidth of the ground motion. A linear model is assumed for ω f (τ ) and ζ f (τ ) is constant during the entire signal duration:

ω f (τ ) = ω mid + ω (τ -t mid ) and ζ f (τ ) = ζ f , (44)
in which t mid is the instant at which 45% of the expected Arias intensity I a is reached, ω mid is the filter frequency at instant t mid and ω is the slope of linear evolution of ω f (τ ). After being normalized by the standard deviation σ h (t), the integral in Eq. ( 42) becomes a unit variance process with time-varying frequency and constant bandwidth, which represents the spectral nonstationarity of the ground motion.

The non-stationarity in intensity is then captured by the modulation function q(t, α). This time-modulating function determines the shape, intensity and duration of the motion as follows:

q(t, α) = α 1 t α2-1 exp(-α 3 t). ( 45
)
The vector of parameters α = (α 1 , α 2 , α 3 ) is directly related to the physical characteristics of the ground motion, namely the expected Arias intensity I a , the time interval D 5-95 between the instants at which the 5% and 95% of I a are reached and the instant t mid .

In the discrete time domain, the synthetic ground motion in Eq. ( 42) becomes:

x(t) = q(t, α) n i=1 s i (t, λ(t i ))U i , (46)
where the standard normal random variable U i represents an impulse at instant t i and s i (t, λ(t i ))

is given by:

s i (t, λ(t i )) = h[t -t i , λ(t i )] k j=1 h 2 [t -t j , λ(t j )] for t i < t k , t k ≤ t < t k+1 , (47)
= 0 for t ≤ t i .

Herein, the parameters ζ and ω of the SDOF oscillator are considered deterministic with values ζ = 0.02 and ω = 5.97 rad/s. The uncertain input vector contains parameters of the oscillator and parameters representing the main intensity and frequency features of the ground motion model ξ = ( , I a , D 5-95 , t mid , ω mid , ω , ζ f ). Table 2 represents the probabilistic distributions associated with the uncertain parameters. The six parameters describing the ground motion are considered dependent with a Nataf distribution (a.k.a Gaussian copula) [START_REF] Liu | Multivariate distribution models with prescribed marginals and covariances[END_REF][START_REF] Lebrun | An innovating analysis of the Nataf transformation from the copula viewpoint[END_REF]. The correlation matrix is given in Table 3. We now aim at building the metamodel for representing the displacement time histories y(t).

For this purpose, N = 200 analyses of the system were conducted with the input parameters generated by Latin hypercube sampling. The equation of motion was solved by means of the Matlab solver ode45 (explicit Runge-Kutta method with relative error tolerance 1 × 10 -3 ). Note that for the sake of consistency, all the synthetic motions were generated with the total duration T = 30 s and time step dt = 0.005 s. In the first place, a NARX model structure was chosen 29)). The responses are reconstructed using the computed coefficients, leading to the errors k (Eq. ( 32)). The best NARX model achieves a mean relative error ¯ = 7.4 × 10 -4

g i (t) = y l (t -j) x m (t -k) with l + m ≤ 3, 0 ≤ l ≤ 3, 0 ≤ m ≤ 1, j = 1, 2, k = 0,
over 200 experiments and contains 7 terms, namely the constant term, x(t -2), x(t), y(t -2), y(t -1), y 2 (t -2), y 3 (t -1).

In the next step, we represented the NARX coefficients by adaptive sparse PCEs of order up to 20 with maximum interaction rank r = 2 and truncation parameter q = 1. The PCEs of the NARX coefficients have LOO errors smaller than 7.34 × 10 -4 . The optimal PCE selected by the adaptive scheme is of total degree 3.

A validation set of size n val = 10 4 was pre-computed by the numerical Matlab solver. The displacements were then predicted by the PC-NARX model. Figure 7 depicts two specific excitations and the corresponding response time histories. PC-NARX provides predictions that are in remarkable agreement with the actual responses. Over 10 4 validation trajectories, the mean relative error is ¯ val = 3.53 × 10 -2 . Less than 5% of those simulations exhibit a relative error val,i exceeding 0.1. Note that predicting the response of a mechanical system subject to nonstationary excitation is never an easy task. From our experience, it is of no interest to apply time-frozen PCEs to this type of problems. informative therefore not presented herein.

Finally, we compare the maximal displacements predicted by PC-NARX with the actual values computed by the numerical solver. Figure 9 shows that overall the predictions are quite accurate with a relative error val,max = 0.7 × 10 -2 . The simulations with large discrepancies between the prediction and the actual values correspond to large displacements, i.e. they belong to the domain of rare events (the upper tail of the distribution) that was not well represented in the small training set of size 200.

Bouc-Wen oscillator subject to synthetic ground motions

In this application, we consider the SDOF Bouc-Wen oscillator [START_REF] Kafali | Seismic fragility analysis: Application to simple linear and nonlinear systems[END_REF] subject to stochastic excitation. The equation of motion of the oscillator reads: sets of parameters. One observes that the velocity trajectories are relatively well predicted by PC-NARX. Indeed, the mean relative error over 10 4 validations is ¯ val = 1.82×10 -2 and only 3% of those time-histories have a relative error val,i exceeding 0.1. Despite the high accuracy of the PC-NARX model for the velocity, the predicted displacements exhibit some slight discrepancies with respect to the actual trajectories. The inaccuracies are slightly large in the later phase of the considered time domain due to the occurrence of residual displacements, while remaining acceptable.

   ÿ(t) + 2 ζ ω ẏ(t) + ω 2 (ρ y(t) + (1 -ρ) z(t)) = -x(t) ż(t) = γ ẏ(t) -α | ẏ(t)| |z(t)| n-1 z(t) -β ẏ(t) |z(t)| n ( 48
Figure 12 compares the time-dependent standard deviation of the two response quantities predicted by PC-NARX with those obtained from Monte Carlo simulation. The high frequency content is observed in the evolutionary response statistics, which was also reported in the literature when applying Monte Carlo simulation to a nonlinear structure [START_REF] Kougioumtzoglou | Nonlinear MDOF system stochastic response determination via a dimension reduction approach[END_REF].The standard deviation of the velocity is remarkably well captured by PC-NARX (relative error val,std = 0.39 × 10 -2 ). For the displacement, the standard deviation tends to increase in time, which is different from the Duffing oscillator that does not exhibit residual displacement.

The discrepancy between the prediction and the actual time histories is also increasing in time.

However, the resulting relative error remains rather small ( val,mean = 1.57 × 10 -2 ). (6.32, 48.24, 0.05, 9.08, 4.53, 19.76, -0.11, 0.28)

Discussion

In this section, we discuss certain features of the PC-NARX approach. It is clear that PC-NARX outperforms traditional time-frozen PCEs. The latter are reported to fail in representing the random response of dynamical systems. The reason is that PCEs are designed for propagating uncertainties, but not for capturing the evolutionary behaviour in time domain. In order to use PCEs in this context, one should utilize a specific tool to tackle the dynamical aspect of the problem. For instance, one can use the available system of equations describing the model in an intrusive way. This approach is, however, also known to deteriorate in long term integration due to the accumulation of error in time (Wan and Karniadakis, 2006a;[START_REF] Ghosh | Applicability of the spectral stochastic finite element method in time-dependent uncertain problems[END_REF]. One also apprehends the remarkable performance of PC-NARX approach, which handles separately the two aspects of the problem, namely dynamics and uncertainties with the specific tools of NARX modelling and PCEs, respectively. PC-NARX can also be considered to belong to a class of approaches where a pre-processing of the response time histories is conducted before applying PCEs. The response trajectories are first projected onto a different space, for instance phase-space [START_REF] Desai | Uncertainty quantification of a nonlinear aeroelastic system using polynomial chaos expansion with constant phase interpolation[END_REF], transformed time space [START_REF] Le Maître | Asynchronous time integration for polynomial chaos expansion of uncertain periodic dynamics[END_REF] or space of NARX functions (Spiridonakos and Chatzi, 2015a).

In the case of PC-NARX, the NARX functions form a "coordinate system" that moves along the trajectory. Due to the projection of the trajectories on the suitably selected set of basis functions, the projection coefficients become smooth functions of the uncertain parameters, therefore they can be represented effectively by low order PCEs.

The PC-NARX approach also has its limitations. It uses data observed in discrete time.

When the recorded data is sparse in the time domain, i.e. the time step is too large, the mechanism that relates the current state with the previous state and excitation cannot be re- vealed. In addition, the recursive form of the PC-NARX formulation renders it more difficult for post-processing. For instance, closed-form solutions for time-dependent statistics and sensitivity analysis are currently not available.

In general, we are convinced that approaches which make use of interdisciplinary tools such as PC-NARX are promising. In other words, PCEs should be combined with tools designed for predicting non-linear time series in uncertainty quantification of stochastic dynamical systems.

Conclusions

Polynomial chaos expansions are a versatile tool for uncertainty quantification, which have proven effective in various fields. However, the propagation of uncertainties in stochastic dynamical systems remains a challenging issue for PCEs as well as all other metamodelling techniques. It is widely known that PCEs fail to capture the long-term dynamics of the underlying system. Therefore, a specialized tool must be used to handle this aspect of the problem. Nonlinear autoregressive with exogenous input (NARX) models are universally used in the field of system identification for revealing the dynamical behaviour from observed time series. The combination of NARX and PCEs have been proposed recently and shown great effectiveness in the context of stochastic dynamics.

In this paper we introduced the least angle regression (LARS) technique for building PC-NARX models of nonlinear systems with uncertainties subject to stochastic excitations. In particular, the approach consists in solving two linear regression problems. LARS proves suitable for selecting both the appropriate NARX and PCE models, which plays a crucial role in the approach. The LARS-based PC-NARX approach is applied to predict the response time histories of three benchmark engineering models subject to non-stationary stochastic excitation with remarkable accuracy. It is important to observe that, to our knowledge, no generic nonintrusive tools based on PCEs were available so far to build surrogate models of dynamical system responses based on a set of computer experiment.

Investigation on the application of the proposed approach to more complex problems is ongoing, with particular interests being focused on the seismic analysis of structures.

  n x and n y denote the maximum input and output time lags, ϑ i (ξ) are the coefficients of the NARX model, g (t, ξ) is the residual error, with zero mean Gaussian distribution and variance σ 2 (t). The proposed NARX model differs from the classical NARX model in the fact that the coefficients ϑ i (ξ) are functions of the uncertain input parameters ξ instead of being deterministic. The stochastic coefficients ϑ i (ξ) of the NARX model are then

Figure 1 :

 1 Figure 1: Computation of the LARS-based PC-NARX model

Figure 3 :

 3 Figure 3: Quarter car model -Instantaneous displacements: comparison of the two approaches.

Figure 5 :

 5 Figure 5: Quarter car model -Two particular trajectories and their predictions by means of timefrozen PCEs and PC-NARX.

Figure 6 :

 6 Figure 6: Quarter car model -Mean and standard deviation of the trajectories: comparison of the two approaches.

Figure 8 Figure 7 :

 87 Figure8represents the evolutionary standard deviation of the displacement. The curve computed by PC-NARX is in excellent agreement with the reference one with a relative error val,Std = 0.5 × 10 -2 . The mean trajectory, which is slightly fluctuating around zero, is not

Figure 8 :

 8 Figure 8: Duffing oscillator -Standard deviation trajectory of the displacement.

Figure 9 :

 9 Figure 9: Duffing oscillator -Maximal displacements.

  ) in which ζ is the damping ratio, ω is the fundamental frequency, ρ is the post-to pre-yield stiffness ratio, γ, α, β, n are parameters governing the hysteretic loops and the excitation x(t) is a ground motion generated by means of the probabilistic model presented in Section 5.2. Deterministic values are used for the following parameters of the Bouc Wen model: ζ = 0.02, ρ = 0, γ = 1, n = 1, β = 0. The remaining parameters are considered indepen-

Figure 10 :

 10 Figure 10: Bouc Wen oscillator -Experiment leading to the selected NARX structure.

Figure 13

 13 Figure 13 compares the maximum values of velocity and displacement predicted by PC-NARX with those values computed by the numerical solver. Despite the complexity of the problem, the predictions are remarkably consistent with the true values, with sufficiently small validation errors of val,max(|v(t)|) = 9 × 10 -3 and val,max(|y(t)|) = 1.9 × 10 -2 . The accurate predictions allow one to obtain the probability density functions of the maximum responses that are in good agreement with the reference functions, as shown in Figure 13(B) and Figure 13(D).

Figure 11 :

 11 Figure 11: Bouc Wen oscillator -Two particular trajectories of velocity v(t) and displacement y(t) and their predictions by means of PC-NARX.

  Probability density function of maximal values of displacement

Figure 13 :

 13 Figure 13: Bouc Wen oscillator -Maximal responses computed by numerical solver and predicted by PC-NARX model.

  Step 1.2: One selects some candidate NARX models being subsets of the specified full model. To this end, one considers the experiments exhibiting a high level of nonlinearity. For instance those experiments can be chosen with measures of nonlinearity or by inspection of the simulations with maximum response values exceeding a specified threshold. For each of the selected experiments, one determines a candidate NARX model containing a subset of the NARX terms specified by the full model. This is done using LARS and the input-output time histories of the considered experiment.

Table 1 :

 1 Parameters of the quarter car model and of the road excitation

	Parameter Distribution Mean & Standard deviation
	k s (N/m 3 )	Gaussian	(2000, 200)
	k u (N/m)	Gaussian	(2000, 200)
	m s (kg)	Gaussian	(20, 2)
	m u (kg)	Gaussian	(40, 4)
	c (N s /m)	Gaussian	(600, 60)
	Parameter Distribution	Support
	A (m)	Uniform	[0.09, 0.11]
	ω (rad/s)	Uniform	[1.8 π, 2.2 π]

Table 2 :

 2 Marginal distributions of the stochastic ground motion model (after Rezaeian and Der Kiureghian (2010)) and of the Duffing oscillator's non linearity .

	Parameter	Distribution	Support Mean Standard deviation
	I a (s.g)	Lognormal	(0, +∞) 0.0468	0.164
	D 5-95 (s)	Beta	[5, 45]	17.3	9.31
	t mid (s)	Beta	[0.5, 40]	12.4	7.44
	ω mid /2π (Hz)	Gamma	(0, +∞)	5.87	3.11
	ω /2π (Hz)	Two-sided exponential [-2, 0.5] -0.089	0.185
	ζ f (.)	Beta	[0.02, 1] 0.213	0.143
	(1/m 2 )	Uniform	[90, 110]	100	5.773

Table 3 :

 3 Correlation matrix of the Nataf distribution of the stochastic ground motion model (after[START_REF] Rezaeian | Simulation of synthetic ground motions for specified earthquake and site characteristics[END_REF]).

		I a	D 5-95 t mid ω mid	ω	ζ
	I a	1	-0.36	0.01 -0.15 0.13 -0.01
	D 5-95 -0.36	1	0.67 -0.13 -0.16 -0.2
	t mid	0.01	0.67	1	-0.28 -0.2 -0.22
	ω mid -0.15 -0.13 -0.28	1	-0.2	0.28
	ω	0.13	-0.16	-0.2	-0.2	1	-0.01
	ζ	-0.01	-0.2	-0.22 0.28 -0.01	1
	in which the model terms are polynomial functions of past values of the output and excitation

  1, 2. The chosen full NARX model contains 10 terms. Next, candidate NARX models were computed. To this end, we selected the simulations with maximum displacement exceeding a large threshold, i.e. max(|y(t)|) > 0.07 m, leading to 19 selected experiments. For each simulation previously selected, LARS was applied to the initial full NARX model to detect the most relevant NARX terms constituting a candidate NARX model. This procedure resulted in 12 candidates in total. The NARX coefficients corresponding to each candidate model are then computed for each simulation in the ED by means of ordinary least squares (Eq. (

In what follows, the term "prediction" is employed to refer to the NARX model's so-called "simulation mode" as addressed in signal processing literature, which stands for the estimation of the response relying only on its initial condition and feedback of the excitation. The term "prediction" is used however because it is the standard wording in the surrogate modelling community.

dent random variables with associated distributions given in Table 4. The uncertainties in the considered system is therefore characterized by means of vector of uncertain parameters ξ = (ω, α, I a , D 5-95 , t mid , ω mid , ω , ζ). (explicit Runge-Kutta method with relative error tolerance 1 × 10 -3 ) with the total duration T = 30 s and time step dt = 0.005 s as in the previous example. In the first place, a NARX model structure is chosen, in which the model terms are

The use of absolute terms has proven effective in capturing the hysteretic behaviour of nonlinear systems in (Spiridonakos and Chatzi, 2015a). The initial NARX model contains 19 terms in total.

Next, the candidate NARX models were computed. For this purpose, we selected the simulations with maximum velocity exceeding a large threshold, i.e. max(|v(t)|) > 0.25 m/s and obtained 15 experiments. LARS was applied to the initial full NARX model to detect the most relevant NARX terms constituting a candidate NARX model from each simulation previously selected. This procedure resulted in 11 candidates in total. OLS (Eq. ( 29)) is used to determine the NARX coefficients corresponding to each NARX candidate model for all the simulations.

To evaluate the accuracy of the NARX candidate, Eq. ( 32) is used to compute the error indicators. The most appropriate NARX model achieves a mean relative error ¯ = 6.27 × 10 -4

over 200 experiments and contains 12 terms, namely constant term, x(t -4), x(t -4) |v(t -1)|, 10 depicts the experiment from which the most appropriate NARX model is selected. Note that the nonlinear behaviour is noticeable and the oscillator exhibits a residual displacement after entering the domain of nonlinearity.

Then we represented the NARX coefficients by sparse PCEs. The optimal polynomial of order p = 2 was found adaptively with maximum interaction order r = 2 and truncation parameter q = 1 so that the resulting PC-NARX model led to the smallest error when reconstructing the responses in the ED. The PCEs of the NARX coefficients have LOO errors smaller than 1.68 × 10 -4 . The PC-NARX model of the velocity was obtained and used for predicting the velocity on the validation set. The displacement time history is then obtained by integration.

Figure 11 depicts two specific velocity and displacement trajectories due to distinct validation