Biswarup Bhattacharyya
email: biswarup.bhattacharyya@ifsttar.fr

Eric Jacquelin
email: eric.jacquelin@univ-lyon1.fr

Denis Brizard
email: denis.brizard@ifsttar.fr

A Kriging-NARX model for uncertainty quantification of nonlinear stochastic dynamical systems in time domain

Keywords: NARX model, Kriging, Nonlinear stochastic dynamical system, LARS, Uncertainty quantification

A novel approach, referred to as sparse Kriging-NARX (KNARX), is proposed in the present paper for the uncertainty quantification of nonlinear stochastic dynamical systems.

It combines the

Nonlinear AutoRegressive with eXogenous (NARX) input model with the high fidelity surrogate model Kriging. The sparsity in the proposed approach is introduced in the NARX model by reducing the number of polynomial bases using the least angle regression (LARS) algorithm. Sparse KNARX captures the non-linearity of a problem by the NARX model, whereas the uncertain parameters are propagated using the Kriging surrogate model and, further, LARS makes the model efficient. The accuracy and the efficiency of the sparse KNARX is measured through uncertainty quantification applied to three nonlinear stochastic dynamical systems. The time dependent mean and standard deviation are predicted for all the numerical examples. Instantaneous stochastic response characteristics and maximum absolute response are also predicted. All the results are compared with the full scale MCS results and a mean error is calculated for all the numerical problems to measure the accuracy.

INTRODUCTION

Uncertainty is an inherent property for any type of real world engineering problem which may arise due to material heterogeneity, variability in dimensions, external forces. It is evident from various studies [START_REF] Grigoriu | Response of dynamic systems to poisson white noise[END_REF]Lucor et al. 2004;[START_REF] Kundu | Transient Response of Structural Dynamic Systems with Parametric Uncertainty[END_REF][START_REF] Chakraborty | Polynomial Correlated Function Expansion for Nonlinear Stochastic Dynamic Analysis[END_REF][START_REF] Mai | Surrogate models for oscillatory systems using sparse polynomial chaos expansions and stochastic time warping[END_REF]) that these sources of uncertainties are always present in a dynamical system. Now, if the input quantities are uncertain, it is obvious that the response characteristics would also be uncertain. The propagation of those uncertain input parameters through a dynamical system and the prediction of uncertain response characteristics are the main objectives of this paper.

The most used method for uncertainty quantification of any engineering problem is Monte Carlo simulation (MCS) [START_REF] Sheppard | Computer simulation of stochastic processes through model-sampling (Monte Carlo) techniques[END_REF][START_REF] Muscolino | Monte Carlo simulation in the stochastic analysis of non-linear systems under external stationary Poisson white noise input[END_REF]. Mostly, MCS is used to predict the stochastic response characteristics for any type of problem and often result found by MCS is used as the principal result for any problem. However, the accuracy of MCS greatly depends on the number of simulations because MCS requires a large number of model evaluations to predict an accurate result for a complex engineering problem. Thus, for a nonlinear dynamical system, this method is highly time consuming and very much computationally expensive. To overcome the issues of MCS, several surrogate models have been developed by the researchers (Lucor and Karniadakis 2004;Lucor et al. 2004;[START_REF] Gerritsma | Time-dependent generalized polynomial chaos[END_REF][START_REF] Kundu | Transient Response of Structural Dynamic Systems with Parametric Uncertainty[END_REF][START_REF] Luchtenburg | Long-time uncertainty propagation using generalized polynomial chaos and flow map composition[END_REF][START_REF] Mai | Surrogate models for oscillatory systems using sparse polynomial chaos expansions and stochastic time warping[END_REF] in the last few decades for stochastic dynamical systems. These surrogate models are the combination of MCS and some kind of polynomials which maintain a trade-off between efficiency and accu-racy such that the stochastic response behavior can be predicted with higher accuracy and lower computational cost. In this regards, one class of surrogate model widely known as polynomial chaos expansion (PCE) (Xiu and Karniadakis 2002) has been used extensively in the last few decades for the uncertainty propagation in stochastic dynamical systems (Lucor et al. 2004;[START_REF] Wan | An adaptive multi-element generalized polynomial chaos method for stochastic differential equations[END_REF][START_REF] Kundu | Transient Response of Structural Dynamic Systems with Parametric Uncertainty[END_REF][START_REF] Jacquelin | Polynomial chaos expansion in structural dynamics: Accelerating the convergence of the first two statistical moment sequences[END_REF][START_REF] Jacquelin | Polynomial chaos-based extended Padé expansion in structural dynamics[END_REF][START_REF] Ozen | A dynamical polynomial chaos approach for long-time evolution of SPDEs[END_REF].

PCE was proposed only to propagate Gaussian input random variables in [START_REF] Wiener | The homogeneous chaos[END_REF] which was further improved in (Xiu and Karniadakis 2002) to account the other types of random variables. This method is also known as generalized PCE (gPCE). Most of the methods developed till date have utilized gPCE as the basic tool for propagating uncertainties in the stochastic dynamical systems. Therefore, Wan and Karniadakis [START_REF] Wan | An adaptive multi-element generalized polynomial chaos method for stochastic differential equations[END_REF] developed a multi-element gPCE to capture the non-linearity of the problem. This method decomposes the domain into multiple segments (when the variance of the stochastic response becomes too worse) such that the non-linearity and the probability density function (PDF) of the response can be captured efficiently. A similar method has been developed in [START_REF] Gerritsma | Time-dependent generalized polynomial chaos[END_REF]) and [START_REF] Luchtenburg | Long-time uncertainty propagation using generalized polynomial chaos and flow map composition[END_REF]) using some different criteria. However, for a system of ordinary differential equations (ODE) having three random variables, it is shown in [START_REF] Gerritsma | Time-dependent generalized polynomial chaos[END_REF]) that PCE needs to be constructed at every time-step. On the other hand, a different type of development has been made in [START_REF] Maitre | Asynchronous Time Integration for Polynomial Chaos Expansion of Uncertain Periodic Dynamics[END_REF] and [START_REF] Mai | Surrogate models for oscillatory systems using sparse polynomial chaos expansions and stochastic time warping[END_REF]) by rescaling the time domain which is often called time warping PCE. However, it can be seen clearly in [START_REF] Mai | Surrogate models for oscillatory systems using sparse polynomial chaos expansions and stochastic time warping[END_REF] that the time warping PCE cannot predict the variance of the stochastic response properly in the later time for nonlinear dynamical system with a few random variables. In contrast, a surrogate model has been recently proposed [START_REF] Spiridonakos | Metamodeling of dynamic nonlinear structural systems through polynomial chaos NARX models[END_REF][START_REF] Mai | Surrogate modeling for stochastic dynamical systems by combining nonlinear autoregressive with exogenous input models and polynomial chaos expansions[END_REF]) by combining PCE and Nonlinear Auto-Regressive with eXogenous input (NARX) [START_REF] Billings | Identification of MIMO non-linear systems using a forward-regression orthogonal estimator[END_REF][START_REF] Chen | Non-linear system identification using neural networks[END_REF][START_REF] Wei | Improved parameter estimates for non-linear dynamical models using a bootstrap method[END_REF][START_REF] Billings | Nonlinear System Identification: NARMAX Methods in the Time, Frequency, and Spatio-Temporal Domains[END_REF]) model to account for the dynamic behavior and the non-linearity of the problem. In this method, the non-linearity of the problem is taken care by the NARX model while the uncertain parameters of the dynamical system are propagated by PCE. However, the proposed model still uses high degree polynomials even for nonlinear single degree of freedom (SDOF) dynamical systems [START_REF] Mai | Surrogate modeling for stochastic dynamical systems by combining nonlinear autoregressive with exogenous input models and polynomial chaos expansions[END_REF]) which increases the computational cost drastically.

Another class of surrogate model widely known as Gaussian process regression or, Kriging surrogate model [START_REF] Krige | A Statistical Approach to Some Basic Mine Valuation Problems on the Witwatersrand[END_REF][START_REF] Santner | The design and analysis of computer experiments[END_REF] has proved the efficiency in predicting the stochastic response characteristics for high-dimensional problems [START_REF] Mukhopadhyay | A Critical Assessment of Kriging Model Variants for High-Fidelity Uncertainty Quantification in Dynamics of composite Shells[END_REF][START_REF] Bhattacharyya | A Critical Appraisal of Design of Experiments for Uncertainty Quantification[END_REF]. However, if one wants to predict the time-dependent stochastic response behavior of an uncertain dynamical system, a new Kriging model is required to fit for each time-instances and almost seems the limitation like MCS. It is also mentioned in [START_REF] Chakraborty | Polynomial Correlated Function Expansion for Nonlinear Stochastic Dynamic Analysis[END_REF][START_REF] Chatterjee | An efficient sparse Bayesian learning framework for stochastic response analysis[END_REF]) that for a non-linear dynamical systems, Kriging requires a large number of model evaluations. However, the efficiency in predicting the stochastic response behavior for the engineering system having spatial random parameters has been well established through several researches [START_REF] Simpson | Kriging Models for Global Approximation in Simulation-Based Multidisciplinary Design Optimization[END_REF][START_REF] Tong | A hybrid algorithm for reliability analysis combining Kriging and subset simulation importance sampling[END_REF][START_REF] Sugai | Application of Kriging method into practical estimations of earthquake ground motion hazards[END_REF][START_REF] Huang | Assessing small failure probabilities by AKSS: An active learning method combining Kriging and Subset Simulation[END_REF][START_REF] Lu | Improved Kriging with extremum response surface method for structural dynamic reliability and sensitivity analyses[END_REF]), but the use of Kriging surrogate model is still questionable for nonlinear stochastic dynamical systems.

Therefore, several improvements have been made by the researchers [START_REF] Kersaudy | A new surrogate modeling technique combining Kriging and polynomial chaos expansions -Application to 29 uncertainty analysis in computational dosimetry[END_REF][START_REF] Chakraborty | An efficient algorithm for building locally refined hpadaptive H-PCFE: Application to uncertainty quantification[END_REF][START_REF] Bayarri | Computer model validation with functional output[END_REF][START_REF] Chatterjee | Refined sparse Bayesian learning configuration for stochastic response analysis[END_REF] to alleviate the issue of the curse of dimensionality and account for the non-linearity of the response behavior, but still these methods are suffering to describe the time dependent response behavior of random non-linear dynamical systems efficiently. To address this issue, a surrogate model is constructed in the present paper in a similar way to [START_REF] Mai | Surrogate modeling for stochastic dynamical systems by combining nonlinear autoregressive with exogenous input models and polynomial chaos expansions[END_REF] by combining the NARX model with the Kriging surrogate model. Therefore, the issue of capturing the high order non-linearity of a stochastic dynamical system can be resolved by utilizing the NARX model and the uncertainties are propagated by the Kriging. Indeed, NARX model has been investigated by very few researchers in context of stochastic dynamical systems [START_REF] Spiridonakos | Metamodeling of dynamic nonlinear structural systems through polynomial chaos NARX models[END_REF][START_REF] Mai | Surrogate modeling for stochastic dynamical systems by combining nonlinear autoregressive with exogenous input models and polynomial chaos expansions[END_REF]Worden et al. 2018) as compared to system identification tool for deterministic nonlinear dynamical systems [START_REF] Billings | Identification of MIMO non-linear systems using a forward-regression orthogonal estimator[END_REF][START_REF] Wei | Improved parameter estimates for non-linear dynamical models using a bootstrap method[END_REF]Xie et al. 2009;Zhang and Li 2015).

The rest of the paper is organized as follows. A brief description of Kriging and NARX model are given in section 2 and section 3 respectively. Section 4 introduces the proposed model along with an algorithm to implement it for dynamical systems. Further, the applicability of the proposed approach for dynamical systems is illustrated through three typical nonlinear stochastic dynamical systems in section 5 and section 6 describes the important conclusion drawn from the present study.

KRIGING

Kriging is a well established surrogate model technique which combines a regression function with a Gaussian process. Kriging was first introduced in geostatistics [START_REF] Matheron | Principles of geostatistics[END_REF]) and afterwards it was used for the analysis of computer experiments [START_REF] Sacks | Design and Analysis of Computer Experiments[END_REF]).

One of the important characteristics of Kriging is that it predicts the quantity of interest in a finite region by an unbiased estimator. A brief description of Kriging is presented in this section.

Construction of model

Consider Ξ = {ξ 1 , ξ 2 , . . . , ξ d } as the d-dimensional input random variables in the probability space P, where Ξ ∈ D ⊂ R d . For N number of sample points, Y = {y (Ξ 1 ) , y (Ξ 2 ) , . . . , y (Ξ N )} T are the corresponding function evaluations. In Kriging, the performance function for a realization is given by:

M (Ξ) = β T φ (Ξ) + Z (Ξ) (1)
In Equation 1, the first part of the expression β T φ (Ξ) represents the regression part of the model which can also be written as:

β T φ (Ξ) = P i=1 β i φ i (Ξ) (2)
φ i (Ξ) are the polynomial basis functions which are constructed with the normalized input variables and β i are the corresponding coefficients of the basis functions. P represents the total number of terms in the basis function which depends on the degree of the polynomial.

According to the form of the polynomial, Kriging model has several variants in the literature [START_REF] Mukhopadhyay | A Critical Assessment of Kriging Model Variants for High-Fidelity Uncertainty Quantification in Dynamics of composite Shells[END_REF]. A universal Kriging model is utilized in the present work (Equation 2). The second part of Equation 1 defines the Gaussian process with mean zero and the process covariance is

cov [Z (Ξ i ) , Z (Ξ j )] = σ 2 Z R (Ξ i , Ξ j ) ; i, j = 1, 2, . . . , N (3)
where, Ξ i and Ξ j are two different sample points and σ 2 Z is the process variance. R (Ξ i , Ξ j ) is the auto-correlation function between two sample points. A variety of auto-correlation functions have been adopted by the researchers [START_REF] Kaymaz | Application of kriging method to structural reliability problems[END_REF][START_REF] Bhattacharyya | A Critical Appraisal of Design of Experiments for Uncertainty Quantification[END_REF][START_REF] Sacks | Design and Analysis of Computer Experiments[END_REF]) such as linear, exponential or Gaussian. The mostly used Gaussian auto-correlation function is utilized in the present paper which is given by:

R (Ξ i , Ξ j ) = d k=1 exp -θ k (ξ i,k -ξ j,k ) 2 (4)
where θ k is the hyper-parameter of the auto-correlation function. The Kriging model parameters β and σ Z , and the hyper-parameter θ k can be determined by the maximum likelihood estimation (MLE) [START_REF] Sacks | Design and Analysis of Computer Experiments[END_REF]). The procedure of computing the parameters is explained in [START_REF] Sacks | Design and Analysis of Computer Experiments[END_REF][START_REF] Santner | The design and analysis of computer experiments[END_REF][START_REF] Kaymaz | Application of kriging method to structural reliability problems[END_REF].

Prediction

Consider any untried sample point for prediction Ξ 0 ∈ R d . The function value at the new sample Ξ 0 is predicted by the best linear unbiased predictor (BLUP). The BLUP and the predicted variance are given by:

M (Ξ 0 ) = w T (Ξ 0 ) Y (5) = βT φ (Ξ 0 ) + r T (Ξ 0 ) R -1 βT φ (Ξ 0 ) (6) σ 2 M (Ξ 0 ) = σ2 Z    1 -φ T (Ξ 0 ) r T (Ξ 0 )    0 F T F R    -1    φ (Ξ 0 ) r (Ξ 0 )       (7)
where, r

(Ξ 0 ) = {R (Ξ 0 , Ξ 1 ) , R (Ξ 0 , Ξ 2 ) , . . . , R (Ξ 0 , Ξ N )} is the correlation matrix between
the new untried sample point and the initial samples, and F is the matrix of basis function at the initial sample points of dimension N × P . R is the correlation matrix for the initial samples. Further, the unknown parameters β and σ 2 Z are calculated as:

β = F T R -1 F -1 F T R -1 Y (8) σ2 Z = 1 N Y -F β T R -1 Y -F β (9)
Once the unknown parameters using the Kriging model are estimated, the prediction of the response at the untried point can be made easily by utilizing Equation 6.

NARX MODEL

Overview of NARX model

A brief overview of NARX model is presented in this section. The time dependent response characteristics of a dynamical system at current time instance can be predicted by the responses of some previous time instances and the excitation of some previous and current time instances through a suitable NARX model [START_REF] Chen | Modelling and analysis of non-linear time series[END_REF][START_REF] Billings | Nonlinear System Identification: NARMAX Methods in the Time, Frequency, and Spatio-Temporal Domains[END_REF]). If we consider a dynamical system having the time dependent excitation, then according to the NARX model, the dynamical system can be expressed as:

y (t) = F [z (t)] + ε (t) (10)
where, z (t) = x (t) , x (t -∆t) , x (t -2∆t) , . . [START_REF] Cantelmo | Adaptive model selection for polynomial NARX models[END_REF], wavelet [START_REF] Billings | A New Class of Wavelet Networks for Nonlinear System Identification[END_REF], sigmoid function [START_REF] Sjöberg | Nonlinear black-box modeling in system identification: a unified overview[END_REF], RBF [START_REF] Li | A fast nonlinear model identification method[END_REF]) and neural network [START_REF] Lin | Learning long-term dependencies in NARX recurrent neural networks[END_REF]. Out of these, the effectiveness of polynomial function has already been proved in the literature [START_REF] Leontaritis | Input-output parametric models for non-linear systems Part II: stochastic non-linear systems[END_REF][START_REF] Cantelmo | Adaptive model selection for polynomial NARX models[END_REF][START_REF] Cheng | An efficient identification scheme for a nonlinear polynomial NARX model[END_REF]. Therefore, a linear-in-parameter form of polynomial has been used in the present paper which is represented by

F [z (t)] = M i=1 ϕ i ψ i [z (t)] (11)
In the above equation, M is the total number of terms in the polynomial basis function,

Ψ k (t, Ξ k ) = {ψ 1 [z k (t, Ξ k )] , ψ 2 [z k (t, Ξ k )] , . . . , ψ M [z k (t, Ξ k )]} T ; Ψ k ∈ R M ×1 (12)
ϕ k = {ϕ 1 , ϕ 2 , . . . , ϕ M }; ϕ k ∈ R 1×M (13)
Now, if we discretize the total time T in N t number of time instances i.e. t ∈ {t 1 , t 2 , . . . , t Nt }, then for N t number of time-steps, the polynomial basis function matrix can be given by

Ψ k ∈ R M ×Nt .
The procedure of choosing the type of polynomial for F [z (t)] and estimating the corresponding coefficients are described in the next section.

Model formulation and parameter estimation of the NARX model

The polynomial basis function for the NARX model can be formulated by two variables, i.e. the excitation and the response of the dynamical system, either as an independent variable basis function [START_REF] Spiridonakos | Metamodeling of dynamic nonlinear structural systems through polynomial chaos NARX models[END_REF] or as a composition of both variables [START_REF] Mai | Surrogate modeling for stochastic dynamical systems by combining nonlinear autoregressive with exogenous input models and polynomial chaos expansions[END_REF] with certain maximum degree of the polynomial basis. For a nonlinear system, it is always important to use the interaction terms [START_REF] Billings | Nonlinear System Identification: NARMAX Methods in the Time, Frequency, and Spatio-Temporal Domains[END_REF]. Hence, a similar kind of polynomial basis function has been adopted in the present paper. Another important issue with the polynomial basis function is the selection of the maximum time lags n xm and n ym for the excitation and the response respectively. The maximum time lags are selected according to [START_REF] Spiridonakos | Metamodeling of dynamic nonlinear structural systems through polynomial chaos NARX models[END_REF][START_REF] Mai | Surrogate modeling for stochastic dynamical systems by combining nonlinear autoregressive with exogenous input models and polynomial chaos expansions[END_REF]) in the present paper, which are two-times the number of DOF of the system.

On the other hand, the computation of the NARX coefficients is one of the challenging tasks in the construction of the NARX model. The coefficients of the NARX model can be computed easily by ordinary least square method (for the k-th sample point). However, it has been found often in the literature [START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on least angle regression[END_REF][START_REF] Spiridonakos | Metamodeling of dynamic nonlinear structural systems through polynomial chaos NARX models[END_REF] that all the terms in the polynomial basis do not get involved in predicting the response characteristics of the system, and reducing the degree of the polynomial may reduce the accuracy of prediction. Thus, it is important to capture the important terms in the polynomial basis function which are solely responsible for the response behavior of a system. Recently, the important NARX basis terms were identified using genetic algorithm in [START_REF] Spiridonakos | Metamodeling of dynamic nonlinear structural systems through polynomial chaos NARX models[END_REF]. However, due to the form of NARX model, the important terms can also be identified with the popularly used regularized least square method, least angle regression (LARS) [START_REF] Efron | Least angle regression[END_REF]. Consequently, LARS is used in the present paper for selecting the important terms in the polynomial basis.

Sparse NARX model

The sparse NARX model is selected in a similar way to [START_REF] Mai | Surrogate modeling for stochastic dynamical systems by combining nonlinear autoregressive with exogenous input models and polynomial chaos expansions[END_REF] in the present paper. At first, it is required to select the samples exhibiting high non-linearity. The responses having highly non-linear behavior can be selected by the measurement of nonlinearity of the response [START_REF] Spiridonakos | Metamodeling of dynamic nonlinear structural systems through polynomial chaos NARX models[END_REF] or, specifying some threshold value for the response series [START_REF] Mai | Surrogate modeling for stochastic dynamical systems by combining nonlinear autoregressive with exogenous input models and polynomial chaos expansions[END_REF]). In the present paper, a combination of both the mentioned methods [START_REF] Spiridonakos | Metamodeling of dynamic nonlinear structural systems through polynomial chaos NARX models[END_REF][START_REF] Mai | Surrogate modeling for stochastic dynamical systems by combining nonlinear autoregressive with exogenous input models and polynomial chaos expansions[END_REF]) has been utilized.

Firstly, the response versus restoring force is plotted arbitrarily without performing any simulations (keeping other parameters constant at their mean values). The intensity of the response is increased upto a certain limit till the nonlinear behavior is noticed in the force-displacement relationship. From the restoring force curve, the threshold value for a response series can be selected easily by observing the starting point of the non-linearity of the response series. Imposing the threshold value would reduce the number of samples to N 1 < N . Now, for each of the N 1 samples (e.g. for the k-th sample), a full NARX model is formulated as follows:

y (t, Ξ k ) = ϕ k Ψ k (Ξ k ) (14)
where, Ψ k is the matrix having all the M terms of the NARX model polynomial basis for all the time-steps with dimension M × N t which is formulated by the current and previous excitation, and previous responses. y (t, Ξ k ) is the response series for the k-th sample point.

At this step, the coefficient vector ϕ k is the unknown.

Remark 1: One should discretize the total time T with a suitable time-step ∆t (small enough to capture the dynamics correctly). Besides, the time-step plays a vital role in case of NARX model. In the present paper, we have used the same time-step for the time integration and the construction of the NARX model.

The LARS algorithm is then employed to select the important terms for each of the N 1 full NARX models.

Remark 2: It is evident that N 1 different NARX models are found at this step. The most important terms for all the NARX models are found using LARS algorithm. However, it may happen that more than one NARX model contains similar terms in the polynomial basis matrix. Therefore, only the N 2 different NARX models are selected at this step.

The unique NARX models are selected from the N 1 number of NARX models which may further reduce the number of NARX models to N 2 ≤ N 1 . Hence, for each of the N 2 sparse NARX models, the coefficients of all the N initial samples are computed by ordinary least square. The response series of the initial N samples are then re-constructed using the coefficients and the selected polynomial bases, and the error for each sample point is predicted as follows:

k = Nt i=1 [y (t i , Ξ k ) -ŷ (t i , Ξ k )] 2 Nt i=1 [y (t i , Ξ k ) -ȳ (Ξ k )] 2 (15)
where, ȳ (Ξ k ) is the time average mean of the k-th response series which is given by:

ȳ (Ξ k ) = 1 N t Nt i=1 y (t i , Ξ k ) (16)
The mean predicted error for all the sample points by a particular sparse NARX model out of N 2 sparse NARX models is given by:

¯ = 1 N N k=1 k (17)
The finally selected sparse NARX model is the one having the predicted mean error for N sample points lower than some threshold value. In the present paper, the threshold value of the mean error is imposed as 1 × 10 -3 for all the examples.

Remark 3: It is of utter importance to mention that if the predicted mean errors for two different NARX model are found same (or, lower than the threshold value) then the sparse NARX model having less number of terms in the polynomial basis is selected as the final sparse NARX model.

KRIGING-NARX MODEL

Consider a dynamical system having some uncertain input parameters Ξ = {ξ 1 , ξ 2 , . . . , ξ d }, then the time dependent response of the system can be represented by a NARX model as:

y (t, Ξ k ) = M i=1 ϕ i (Ξ k )ψ i [z k (t)]; k = 1, 2, . . . , N (18)
It is observed from Equation 18that the coefficients of the NARX model are dependent on the sample points which does not make the model stochastic in nature. To get the independent coefficients of the NARX model, the NARX coefficients are represented by the Kriging surrogate model as given in Equation 1:

ϕ i (Ξ) = β T i φ (Ξ) + Z i (Ξ) ; i = 1, 2, . . . , M (19)
Further, the model is constructed by the Kriging surrogate model as discussed in section 2.

Therefore, the full KNARX model is expressed as: The response quantity of a dynamical system at some untried sample points Ξ 0 can be predicted by BLUP as mentioned in Equation 6in accordance with the full NARX model in an auto-regressive manner.

y(t, Ξ) = M i=1 β T i φ (Ξ) + Z i (Ξ) ψ i [z (t)] (20
ŷ (t, Ξ 0 ) = M i=1 βT i φ (Ξ 0 ) + r T (Ξ 0 ) R -1 βT i φ (Ξ 0 ) ψ i [z (t)] (21)
Equation 21 helps in predicting the time dependent response characteristics of a dynamical system at large number of sample points having d dimensional random input variables.

It is evident from Remark 4 that the computational cost increases with the increase of number of terms in the NARX polynomial basis. For that reason, the sparse NARX model as discussed in section 3 is used in accordance with the Kriging surrogate model in the present paper. Consequently, the number of terms for the NARX model is decreased to M 1 < M which ultimately reduces the number of Kriging model calibration (M 1 ) and computational cost. The algorithm for constructing the sparse KNARX model is provided in Table 1.

NUMERICAL APPLICATIONS

The sparse KNARX model as described in the previous section has been utilized for uncertainty quantification of three typical nonlinear stochastic dynamical systems. For each of the examples, the accuracy of the method is measured using the predicted mean error as given in Equation 17 It is well known that a suitable sample point (often known as experimental design point) generation strategy is required to use for the generation of initial samples in step 2 of Table 1.

Consequently, one of the most widely used sampling strategies, Latin hypercube sampling (LHS) [START_REF] Mai | Surrogate modeling for stochastic dynamical systems by combining nonlinear autoregressive with exogenous input models and polynomial chaos expansions[END_REF]; Chatterjee and Chowdhury 2017) has been utilized for the initial N number of sample points generation.

Duffing oscillator

A non-linear duffing oscillator is considered as the first example for the illustration of the proposed model. The governing differential equations of the duffing oscillator are given by:

ÿ (t) + 2ζω ẏ (t) + ω 2 y (t) + εy 3 (t) = x (t) (22) x (t) = A sin (ω x t) (23)
In Equation 22, ω represents the natural frequency considering the undamped linear structure (ε = 0, ζ = 0), ζ is the damping ratio and ε controls the non-linearity of the system.

x (t) denotes the excitation part of the dynamical system which is considered sinusoidal in this case (see Equation 23). The initial conditions are y (0) = 0 and ẏ (0) = 0. The numerical integration has been performed for T = 30 s with a time-step of ∆t = 0.01 s through MATLAB solver ode45. All the parameters of the duffing oscillator are considered uncertain (i.e. Ξ = {ω, ζ, ε, A, ω x }) [START_REF] Mai | Polynomial chaos expansions for uncertain dynamical systems Applications in earthquake engineering[END_REF]. The distribution types along with the parameters of the distribution are listed in Table 2.

Time dependent stochastic displacement of the duffing oscillator is predicted by MCS, Kriging, sparse PCE-NARX [START_REF] Mai | Surrogate modeling for stochastic dynamical systems by combining nonlinear autoregressive with exogenous input models and polynomial chaos expansions[END_REF]) and sparse KNARX model. MCS has been performed with 3 × 10 4 sample points and N = 50 samples have been generated by LHS for the prediction of the stochastic response behavior using Kriging. The Kriging model has been calibrated for each time-step (i.e. 3001 times). The step by step procedure of constructing the sparse KNARX model (according to Table 1) is described below:

1. For the duffing oscillator, d = 5 (see Table 2).

2. The sparse KNARX model is constructed using N = 25 LHS points.

3. For the construction of the sparse KNARX and the sparse PCE-NARX model, it is important to specify the threshold value for capturing the nonlinear response series according to step 3 of Table 1. Full NARX models are to be constructed on those samples which satisfy the criterion of threshold value. For that reason, the displacement y (t) versus restoring force f s = ω 2 [y (t) + εy 3 (t)] has been plotted in Figure 1

ψ i [z (t)] = x l i (t -n x i ∆t) y m i (t -n y i ∆t) (24)
where, x and y are the excitation and the response of the duffing oscillator respectively.

6. The maximum time lags for the excitation and the response are taken as n xm = 2 and n ym = 2 respectively due to 1-DOF system. n x ∈ {0, 1, 2} and n y ∈ {1, 2} represent all the time lags for the excitation and the response respectively. l ∈ {0, 1} and m ∈ {0, 1, 2, 3} are the degrees for the excitation and the response respectively, with a maximum degree of the polynomial 3 (i.e. l i + m i ≤ 3), due to the cubic non-linearity of the system. A total 22 number of terms have been obtained in the polynomial basis matrix utilizing all the above-mentioned criterion which depicts M = 22 and i = 1, 2, . . . , M in Equation 24.

7. Therefore, 2 full NARX models are constructed using the basis function as given in Equation 24, and they have 22 terms in the NARX polynomial basis matrix.

8. The sparse NARX models are then constructed using the LARS algorithm by selecting the most important terms in the bases for both the full NARX models. 9. The sparse NARX models are then retained which have unique set of polynomial bases and for the duffing oscillator, both the sparse NARX models are found unique. The polynomials selected by the LARS algorithm for both the NARX models are listed in Table 3.

10. The coefficients for all the 25 samples are obtained with the ordinary least square for both the sparse NARX model. For this step, explicitly N × N 2 number of ordinary least squares is performed.

11. Then, the response series are reconstructed again in a recursive manner using the coefficients and the polynomial bases.

12. ¯ is then computed by both the sparse NARX models for the N samples.

13. ¯ for both the sparse NARX models is found below the threshold value. Hence, the Model 2 of Table 3 is chosen as the best sparse NARX model (according to Remark 3) and the predicted mean error for the N = 25 samples is found as ¯ = 1.68 × 10 -7 .

14. The coefficient vector ϕ i (Ξ) ∈ R N ×1 ; i = 1, . . . , 8 corresponding to each of the polynomials are uncertain. These 8 coefficient vectors are considered as the uncertain response quantities for the Kriging model. Hence, the unknown parameters of the Kriging model β, σ2 Z and θ k are computed by the MLE for all the 8 NARX coefficients separately.

Consequently, 8 Kriging models are calibrated for converting the model stochastic.

15. The same 3×10 4 MCS samples are used here for the prediction by the surrogate model.

16. The 8 NARX coefficient vectors for the 3 × 10 4 samples are then predicted using the 8 calibrated Kriging models by BLUP. 17. Further, the predicted coefficients and the selected polynomials are used to predict the response series auto-regressively for the 3 × 10 4 MCS samples.

The same N = 25 LHS samples are also utilized for uncertainty quantification using the sparse PCE-NARX [START_REF] Mai | Surrogate modeling for stochastic dynamical systems by combining nonlinear autoregressive with exogenous input models and polynomial chaos expansions[END_REF] model. Here, the sparse adaptive PCE model is constructed using the UQLab module [START_REF] Marelli | UQLab user manual -Polynomial chaos expansions[END_REF]. All the parameters of the sparse PCE model are taken similar to [START_REF] Mai | Surrogate modeling for stochastic dynamical systems by combining nonlinear autoregressive with exogenous input models and polynomial chaos expansions[END_REF] for the present computation. Therefore, the maximum interaction is chosen as 2 and the q-norm is taken as 1 with the degree of the polynomials varying between 1 and 20 [START_REF] Marelli | UQLab user manual -Polynomial chaos expansions[END_REF][START_REF] Mai | Surrogate modeling for stochastic dynamical systems by combining nonlinear autoregressive with exogenous input models and polynomial chaos expansions[END_REF]). However, the predicted response is not converged using the sparse PCE-NARX model with the same number of samples. For that reason, the number of samples for the sparse PCE-NARX model is increased to N = 35 and the time dependent responses are predicted with the sparse PCE-NARX model. Here, the final sparse NARX model is found having M 1 = 9 terms in the NARX polynomial basis which are

y (t -∆t) , y (t -2∆t) , y 2 (t -∆t) , y 3 (t -∆t) , y 3 (t -2∆t) , x (t) , x (t -2∆t) , x (t) y (t -∆t) , x (t -2∆t) y 2 (t -2∆t) .
To reduce further the computational cost, another study has been performed using less number of initial samples with N = 21 (for sparse KNARX). In this case, 3 samples are retained initially based on the criterion of threshold value of the response series and the 3 full NARX models are formulated using the basis function as mentioned in Equation 24.

Further, the sparse NARX models are found by applying the LARS algorithm and all the three models are found as the unique sparse NARX models. Therefore, the coefficients for the N = 21 samples are found by the ordinary least square method and then, all the 21 response series are reconstructed using the coefficients of the sparse NARX models. The finally selected sparse NARX model is the Model 1 of Table 3 and

Bouc-Wen oscillator

A non-linear Bouc-Wen oscillator [START_REF] Bouc | Forced vibration of mechanical systems with hysteresis[END_REF][START_REF] Wen | Method for Random Vibration of Hysteretic Systems[END_REF]) is investigated in this example.

The governing differential equation for the Bouc-Wen oscillator is given by:

ÿ (t) + 2ζω ẏ (t) + ω 2 [ρy (t) + (1 -ρ) w (t)] = -A sin (ω x t) ẇ (t) = γ ẏ (t) -α | ẏ (t)| |w (t)| n-1 w (t) -β ẏ (t) |w (t)| n (25)
where, ζ and ω are the damping ratio and the natural frequency of the oscillator respectively. ρ = 0, γ = 1, β = 0 and n = 1 are considered for this problem. The excitation of the oscillator is considered as x (t) = A sin (ω x t) and w (t) is the hysteric displacement as given in [START_REF] Wen | Method for Random Vibration of Hysteretic Systems[END_REF]). The initial conditions at rest are y (0) = 0, ẏ (0) = 0 and w (0) = 0. Similar to [START_REF] Mai | Surrogate models for oscillatory systems using sparse polynomial chaos expansions and stochastic time warping[END_REF], 5 uncertain parameters are considered for the oscillator, which are Ξ = {ζ, ω, α, A, ω x }. The type of distributions for all the uncertain parameters are given in Table 6.

The stochastic response is calculated for the oscillator in the time domain t ∈ [0 s, 30 s] with a time-step of ∆t = 0.005 s. As it has already been illustrated through the previous example, Kriging is unable to predict the stochastic response behavior for the non-linear stochastic duffing oscillator even with higher samples than the sparse KNARX. Therefore, Kriging has not been utilized from now onward. Similar to the previous example, the sparse KNARX model has been constructed with two different sizes of sample points N = 40 and N = 10. A different type of basis function has been considered for the Bouc-Wen oscillator according to [START_REF] Mai | Surrogate modeling for stochastic dynamical systems by combining nonlinear autoregressive with exogenous input models and polynomial chaos expansions[END_REF], which is given by:

ψ i [z (t)] = x l i (t -n x i ∆t) | ẏ (t -∆t)| m i , ẏl i (t -n y i ∆t) | ẏ (t -∆t)| m i (26)
In Equation 26, the basis function relies on the excitation and the velocity of the oscillator.

Thus, the velocity is computed by the sparse KNARX model and further, the displacement of the system is obtained through numerical integration by utilizing the MATLAB solver ode45.

Due to the hysteric displacement component, the bases are computed by two expressions as given in Equation 26which has already been considered in [START_REF] Mai | Surrogate modeling for stochastic dynamical systems by combining nonlinear autoregressive with exogenous input models and polynomial chaos expansions[END_REF]. l ∈ {0, 1} and m ∈ {0, 1} are considered for the oscillator as single degree is noticed for the excitation and the velocity in Equation 25. n xm = n ym = 4 is considered as the problem can be appraised as 2-DOF system due to an extra hysteric displacement part i.e. n x ∈ {0, . . . , 4}

and n y ∈ {1, . . . , 4}. A total 21 number of terms are found in the full NARX polynomial basis matrix by utilizing Equation 26and imposing all the mentioned conditions.

A threshold similar to the previous example was employed to detect the highly nonlinear samples. The threshold value for the response of the Bouc-Wen oscillator is chosen as max (| ẏ (t)|) > 0.3 m s -1 . The threshold value reduces the number of samples from N = 40 to N 1 = 3 which means only 3 full NARX models are required. Thus, 3 full NARX models are formulated utilizing the basis function of Equation 26and further the LARS algorithm has been employed to make the full NARX models sparse. All the three sparse NARX models are found unique in this step. Therefore, the coefficients for all the 40 initial samples are computed by the ordinary least square method for the 3 sparse NARX models and the mean errors (Equation 17) are predicted by computing the time series in a recursive manner. The finally selected sparse NARX model produces a mean error of ¯ = 1.10 × 10 -5 for the 40 samples and contains 8 terms in the polynomial basis matrix which are

| ẏ (t -∆t)| , x (t) , x (t -4∆t) , x (t -4∆t) | ẏ (t -∆t)| , ẏ (t -∆t) , ẏ (t -4∆t) , ẏ (t -∆t) | ẏ (t -∆t)| , ẏ (t -4∆t) | ẏ (t -∆t)| .
A similar procedure is adopted by substantially reducing the initial number of sample points to N = 10. 4 samples are found mostly non-linear based on the previous threshold value. Out of the 4 full NARX models, 2 are found as unique sparse NARX models after applying the LARS for detecting the most important terms. Therefore, one sparse NARX model is selected based on the predicted mean error from both sparse NARX models. The final sparse NARX model is found to have the same terms in the polynomial basis as the previous one with N = 40.

The stochastic responses are also predicted here by the sparse PCE-NARX model [START_REF] Mai | Surrogate modeling for stochastic dynamical systems by combining nonlinear autoregressive with exogenous input models and polynomial chaos expansions[END_REF]) (the sparse PCE model is constructed here with the same conditions as the previous example). However, the sparse PCE-NARX is unable to predict the stochastic response behavior with N = 40 number of model evaluations for the Bouc-Wen oscillator.

For that reason, the number of model evaluations is increased to N = 50 and the stochastic responses are predicted with the same M 1 = 8 terms in the sparse NARX model as in the sparse KNARX model with N = 40.

This study aims at quantifying the uncertainty associated with the response quantity.

The uncertain response characteristics are predicted for the displacement and the velocity of the Bouc-Wen oscillator. The time dependent mean and the standard deviation of y (t) and ẏ (t) are shown in Figure 6 and Figure 7 respectively. The figures show that both the statistical moments predicted by the sparse KNARX are following the MCS results with utmost accuracy. However, the accuracies of the statistical moments are deteriorating from the initial time-steps using the sparse PCE-NARX model even with more samples. The scatter diagrams and the PDFs at three different time instances (10 s, 20 s and 30 s) are also plotted in Figure 8 and Figure 9 for the displacement and the velocity respectively. It is seen from all the figures that the sparse KNARX performs well in all instances with very few samples. The error metrics for all the instantaneous displacements and velocities are given in Table 7. An excellent accuracy in terms of error ( ) and R 2 value is noticed for all the time instances by the sparse KNARX model.

For the prediction of the safety corridor under uncertainty, the scatter plots and the PDFs of the maximum absolute displacement and velocity are plotted in Figure 10 and Figure 11 respectively. The worst prediction is noticed for the PDFs by the sparse PCE-NARX model, whereas an excellent accuracy is observed in both the cases for the predicted maximum responses by the sparse KNARX model. These plots, which represent the uncertain maximum response behavior, can be utilized to measure the safety margin of the system.

To observe the accuracy of the surrogate models, the mean of the displacement and the velocity are plotted in state space by the MCS, the sparse PCE-NARX and the sparse 22 KNARX in Figure 12. It is evident from the figure that the accuracy in predicting the state space behavior of the response is greatly achieved even using much less number of samples (N = 10) by the sparse KNARX model whereas the sparse PCE-NARX predicted result is the worst even with more samples. Similar to the previous example, the accuracy and the efficiency measurement metrics are given in Table 8. The table is showing clearly that the sparse KNARX has predicted results with higher accuracy by minimum computational cost. The efficiency of the proposed sparse KNARX method is also observed over the recently proposed method time warping PCE [START_REF] Mai | Surrogate models for oscillatory systems using sparse polynomial chaos expansions and stochastic time warping[END_REF] which required N = 100 (>> N = 10) samples to predict the stochastic response behavior for the same example.

A 2-DOF dynamical system

Finally, a 2-DOF dynamical system [START_REF] Mai | Surrogate modeling for stochastic dynamical systems by combining nonlinear autoregressive with exogenous input models and polynomial chaos expansions[END_REF]) has been considered for the applicability of the proposed sparse KNARX model. The dynamical system is shown in Figure 13 and the governing differential equation of the problem is:

m s ÿ1 (t) = -k s [y 1 (t) -y 2 (t)] 3 -c [ ẏ1 (t) -ẏ2 (t)] m u ÿ2 (t) = k s [y 1 (t) -y 2 (t)] 3 + c [ ẏ1 (t) -ẏ2 (t)] + k u [x (t) -y 2 (t)] ( 27
)
where m s is the sprung mass which is connected with the mass m u by a non-linear spring with stiffness k s and a damper with damping coefficient c. k u is a linear spring attached to the ground having a sinusoidal displacement function x (t) = A sin (ω x t).

Similar to [START_REF] Mai | Surrogate modeling for stochastic dynamical systems by combining nonlinear autoregressive with exogenous input models and polynomial chaos expansions[END_REF]), all the parameters of the system are considered uncertain i.e. Ξ = {k s , k u , m s , m u , c, A, ω x }. The mean and the standard deviation of all the uncertain parameters are given in Table 9.

The main aim of this example is to predict the uncertain response y 1 (t) of mass m s which is attached with the non-linear spring. This system has been solved by the sparse KNARX model with two different number of sample points (N = 50 and N = 20). As it has already been seen from the previous two examples that the sparse PCE-NARX model is unable to predict the stochastic response behavior of the dynamical systems properly, the sparse PCE-NARX model is not utilized for this problem. It should be noted that the problem has been solved in the time domain t ∈ [0 s, 30 s] with a time-step of ∆t = 0.01 s by ode45. The initial conditions for the system at rest (t = 0 s) are given by:

y 1 (0) = 0 ẏ1 (0) = 0 y 2 (0) = 0 ẏ2 (0) = 0 (28)
For the construction of the NARX model, the NARX basis function is chosen similar to the first example. The NARX basis function for this problem is given by:

ψ i [z (t)] = x l i (t -n x i ∆t) y m i 1 (t -n y i ∆t) (29)
In Equation 29, l ∈ {0, 1} and m ∈ {0, 1, 2, 3} with l i + m i ≤ 3 are chosen due to cubic nonlinearity of the system.

(t -∆t) , y 1 (t -4∆t) , y 3 1 (t -∆t) , x (t -4∆t) , x (t -4∆t) y 2 1 (t -4∆t)}.
Therefore, only 5 number of Kriging models are required to calibrate for converting the sparse NARX model stochastic which is highly efficient.

The same example has also been solved by reducing the initial number of sample points with N = 20. For this case, initially, 3 samples are selected as the measure of non-linearity and 3 full NARX models are constructed using the basis function as mentioned in Equation 29. Further, all the 3 NARX models are found unique sparse NARX model utilizing the LARS. The predicted mean error for the 20 samples using the finally selected sparse NARX model was ¯ = 1.85 × 10 -6 with 5 terms in the polynomial. Here, the 5 terms

are {y 1 (t -∆t) , y 1 (t -4∆t) , x (t) , x (t -4∆t) , x (t -4∆t) y 2 1 (t -4∆t)}.
The same number of Kriging models are calibrated as the previous case to make the sparse NARX model stochastic.

The time dependent mean and standard deviation of the displacement (y 1 ) and velocity 10.

The error for a particular time instance y(t) is computed using Equation 15. The results

show a high accuracy of the sparse KNARX in predicting the instantaneous response characteristics. The accuracy and the efficiency are also noticed for the sparse KNARX at 30 s (N = 50, y 1 (t) = 2.3 × 10 -3 ) as compared to the recently proposed sparse PCE-NARX model (N = 100, y 1 (t) = 4.21 × 10 -3 ) [START_REF] Mai | Surrogate modeling for stochastic dynamical systems by combining nonlinear autoregressive with exogenous input models and polynomial chaos expansions[END_REF]. The accuracy of the stochastic response behavior has also been checked by plotting the absolute maximum displacement and the velocity in Figure 17 and Figure 18 respectively. A very accurate result is noticed using sample points N = 50.

The mean state space behavior of the predicted responses are also plotted in Figure 19.

It is noticeable that the mean state space behavior is predicted well even with low number of sample points (N = 20) by the sparse KNARX. The global mean error of the predicted time series is measured by using Equation 17. Along with this, the accuracy in predicting the maximum absolute displacement using different sample points are calculated using Equation 15 and R 2 value. All these results along with the efficiency measurement metrics of the proposed sparse KNARX are listed in Table 11. Here, the efficiency of the sparse KNARX

DATA AVAILABILITY

Some or all data, models, or code generated or used during the study are available from the corresponding author by request.

• NARX model 1. Declare the d number of uncertain variables with the type of distribution. 2. Generate N number of sample points for the uncertain variables.

3. Draw the restoring force versus response curve to decide the threshold value for capturing nonlinear response series. 4. Get the response of the dynamical system upto time T for each of the N samples. 5. Select the samples and the response series having high order non-linearity using some threshold (according to step 3) on the response series. This step selects N 1 samples (N 1 < N ). 6. Decide the maximum time lags n xm and n ym for the excitation and response quantity, respectively. Decide the type of NARX polynomial basis function along with the maximum degrees also. 7. For each of the N 1 samples, build the NARX model (see Equation 10 and 11) which has M number of terms. 8. Select the most important terms for each of the N 1 NARX models using LARS [START_REF] Efron | Least angle regression[END_REF]). 9. Select the N 2 different sparse NARX models (see Remark 2) having similar terms in the NARX polynomial basis (N 2 ≤ N 1 ). 10. Perform ordinary least square on N samples for the N 2 sparse NARX models to get the NARX coefficients of all the response series. 11. Reconstruct the N response series using the coefficients computed in step 10 by all the N 2 sparse NARX models. 12. Predict the mean error ¯ for all the N 2 number of sparse NARX models using Equation 17.

13. Select the most appropriate sparse NARX model having ¯ lower than some threshold value (1 × 10 -3 in the present paper) and less number of terms (M 1 ) in the polynomial. 14. Calibrate M 1 < M number of Kriging models using the NARX coefficients as the response parameter and the sample points generated in step 2 as the uncertain input quantities. 15. Generate a large number of new untried samples for the prediction. 16. Predict the NARX coefficients for the untried samples using Kriging models generated in step 14. 17. Predict the response series, in a auto-regressive way (see Equation 21), at the untried samples generated in step 15 using the coefficients of step 16 and the M 1 number of NARX polynomial bases selected in step 13.

ψ

  i [z (t)] are the polynomial basis functions and ϕ i are the corresponding coefficients of the NARX model. The polynomial basis function for the NARX model is dependent on the time varying excitation and response. Consequently, the polynomial basis matrix and the coefficient matrix for a particular sample point k are represented by:

)

  The coefficient of the regression part, β T and the Gaussian process part Z (Ξ) are dependent on the number of terms in the NARX model, whereas the polynomial basis function of the Kriging model for all the NARX coefficients remain the same as the basis is purely dependent on the uncertain input parameters.Remark 4: It should be noted that the Kriging model should be calibrated M times (M is the total number of terms for a full NARX model) because the coefficients of the NARX model for N samples corresponding to a polynomial basis act as the single response quantity for the Kriging model.

  [START_REF] Mai | Surrogate modeling for stochastic dynamical systems by combining nonlinear autoregressive with exogenous input models and polynomial chaos expansions[END_REF]) (taking MCS as the reference) and the value of the coefficient of correlation R 2 . Uncertainty quantification for all the problems is made by predicting the time dependent mean and standard deviation of the responses. Further, the PDFs of the responses are also predicted at some time instances. The computational efficiency is measured through the computational cost (CPU time), the number of surrogate model calibration and the number of initial sample points N . The first problem is also solved by the Kriging model to measure the efficiency and the accuracy of the sparse KNARX model over the Kriging. For the first two examples, the sparse KNARX predicted results are also compared with the recently proposed sparse PCE-NARX surrogate model[START_REF] Mai | Surrogate modeling for stochastic dynamical systems by combining nonlinear autoregressive with exogenous input models and polynomial chaos expansions[END_REF]).

  (keeping other parameters constant at their mean values). From Figure 1, it is clear that the displacement behaves almost linearly for y (t) ∈ [-0.045 m, 0.045 m] and beyond this region, the displacement is highly nonlinear. Consequently, to capture the highly nonlinear response series, the threshold value for the nonlinear displacement is set as max (|y (t)|) > 0.045 m. Thus, only those samples will be picked for the construction of NARX model which falls beyond the region y (t) ∈ [-0.045 m, 0.045 m]. 4. The displacement series are then obtained at the initial N = 25 samples. Out of all the response series, the displacement versus restoring force for two different samples (13-th and 25-th sample point) are plotted in Figure 2. It is clearly seen that the 13-th sample point exhibits higher order non-linearity than the 25-th sample point. It should be noted that the region of displacement value is almost restricted in [-0.02 m, 0.02 m] for the 25-th sample point which is far less than the specified threshold value. On the other hand, the displacement intensities are far beyond the threshold value and behaves nonlinearly for the 13-th sample point. 5. After imposing the above-specified threshold criterion, only N 1 = 2 samples are selected as the highly nonlinear response series. Thus, only 2 full NARX models are required for the duffing oscillator using the sparse KNARX model (N = 25). The basis function for the full NARX model is chosen as:

  the predicted mean error of the selected sparse NARX model is found as ¯ = 1.80 × 10 -7 for the 21 samples. Further, the coefficients corresponding to the sparse NARX polynomial bases are modelled using the Kriging surrogate model. On the contrary, the sparse PCE-NARX model is not able to predict the response with less samples. The time varying stochastic response characteristics are predicted by the time dependent mean and standard deviation which are plotted in Figure 3. The figure clearly depicts the efficiency (initial number of sample points N ) and the accuracy of the sparse KNARX model over the Kriging model. The sprase KNARX has predicted better results with fewer number of model evaluations (N = 21) and the sparse PCE-NARX model is incapable to predict a result even using N = 25 samples. To illustrate the accuracy of the sparse KNARX, the instantaneous response characteristics are also plotted. The scatter diagrams and the PDFs of the predicted response are plotted in Figure 4 at three different time instance (10 s, 20 s and 30 s). The accuracy metrics of the instantaneous response characteristics are listed in Table4, they clearly suggest that the sparse KNARX is highly efficient as compared to the Kriging model. The accuracy of the sparse KNARX model is also comparatively higher using less model evaluations than the sparse PCE-NARX model.The stochastic absolute maximum displacement max (|y (t)|), plotted in Figure5, ultimately measures the safety corridor for the dynamical system. It is seen clearly that the sparse KNARX outperforms the Kriging and the sparse PCE-NARX with very less number of sample points in predicting the PDF of max (|y (t)|), and the accuracy of the sparse KNARX is given by the R 2 value in Table5which is very close to 1.0 by N = 25.Further, the accuracy of the overall model is computed by the mean error (Equation17)for the predicted responses and the error for the predicted max (|y (t)|). The efficiency of the sparse KNARX has already been shown by the initial number of sample points. An accurate result is predicted by the sparse KNARX even with much less number of model evaluations (N = 21) as compared to the other methods. The efficiency of the sparse KNARX is also measured by the number of surrogate model calibrations (n K ) for a method and by the CPU time. All the accuracy and the efficiency measurement metrics are reported in

  n x ∈ {0, 1, . . . , 4}, n y ∈ {1, . . . , 4} are chosen because the system is having 2-DOF. The full NARX model is constructed with this NARX basis function which has M = 58 terms. Initially, after investigating the response versus restoring force, a threshold criterion of responses has been decided as max (|y 1 (t)|) > 1.2 m which reduces the number of samples to N 1 = 10. Consequently, only 10 number of full NARX models are constructed by using the basis function as mentioned in Equation 29. The sparsity is introduced in this step by using the LARS algorithm on the 10 full NARX models to get the unique sparse NARX models which further reduces the number of unique sparse NARX models to N 2 = 5. Thus, 5 unique sparse models are used to get the coefficients of the sparse NARX models by ordinary least square for N = 50 samples. Finally, the mean error for each of the sparse NARX models is computed using Equation 17 and the finally selected sparse NARX model has mean error of ¯ = 2.95 × 10 -4 with 5 terms in the NARX polynomial out of 58 which are {y 1

(

  ẏ1 ) are predicted for this problem. The velocities are predicted by numerical differentiation using the Newton's central difference scheme. The time dependent statistical characteristics of displacement and velocity are plotted in Figure 14 and Figure 15 respectively. It is evident from both the figures that the sparse KNARX performs very well even with very few sample points. Similar to the previous examples, the instantaneous predicted displacements by the sparse KNARX are compared with the MCS results at three different time instances. The instantaneous scatter plots and the PDFs are plotted in Figure 16. The error of the predicted responses at instantaneous time instances along with the R 2 values are given in Table

  model has not been emphasized by the number of surrogate model calibration (n K ), instead it is greatly affected by the initial number of model evaluations (N ). All the results show sparse KNARX perform very well in all instances. CONCLUSION A novel method has been proposed in the present paper for uncertainty propagation and quantification of nonlinear stochastic dynamical systems in the time domain. The well known NARX model (a nonlinear system identification technique) has been coupled with Kriging (a high fidelity surrogate model) to propose the method referred as KNARX. Further, the sparsity in the NARX model has been introduced by utilizing the LARS algorithm to formulate the sparse KNARX model. The proposed sparse KNARX model has been used for uncertainty quantification on three typical nonlinear stochastic dynamical systems: a 1-DOF duffing oscillator, a Bouc-Wen oscillator and a 2-DOF dynamical system. Time dependent mean and standard deviation are predicted for all the examples. Along with this, PDFs are plotted for instantaneous response characteristics at some time instances and for maximum absolute response. All the results have been compared with full scale MCS. The first example is also solved with Kriging, whereas the first two examples are computed by the sparse PCE-NARX model (Mai et al. 2016). The accuracy has been measured by a global mean error of the predicted model whereas the efficiency has been measured by number of initial sample points, number of surrogate model evaluations and CPU time. The results for the first example have shown that sparse KNARX outperforms Kriging in all cases even with a very limited number of initial samples. Consequently, Kriging has not been utilized for the last two examples. On the other hand, the sparse PCE-NARX model predicted acceptable results by using a higher number of sample points than the sparse KNARX for the first two examples. Along with this, the efficiency of the sparse KNARX model is higher as compared to the sparse PCE-NARX model considering the CPU time. All the results predicted by the sparse KNARX are in excellent agreement with the full scale MCS results in all instances. The introduction of sparsity has reduced substantially the number of Kriging model evaluations which ultimately enhanced the efficiency of the sparse KNARX model. All the examples have also been solved by reducing number of model evaluations to check the accuracy of the fitted model with less sample points. All the results have been predicted very close to the MCS results even with very few model evaluations for all the examples. This study suggests that reducing the number of model evaluation does not affect the accuracy too much due to the fact that the major non-linearity of the system has been captured by the sparse NARX model. In conclusion, the proposed sparse KNARX outperforms the state-of-the-art methods (Mai and Sudret 2017; Mai et al. 2016) in predicting the time dependent stochastic response characteristics for nonlinear dynamical systems. The accuracy and the efficiency of the proposed method has also been illustrated through several examples. The proposed sparse KNARX model is limited to use for the nonlinear dynamical systems having parametric uncertainties only. Therefore, this model is not useful for the dynamical systems with white noise excitation for the prediction of time dependent statistical responses.

-FIG. 2 :FIG. 3 :FIG. 4 :FIG. 6 :FIG. 7 :FIG. 8 :FIG. 9 :FIG. 10 :FIG. 11 :FIG. 12 :FIG. 14 :FIG. 15 :FIG. 17 :FIG. 18 :FIG. 19 :

 23467891011121415171819 FIG. 1: Displacement versus restoring force plot for the duffing oscillator

Table 5 .

 5 Table 5 is suggesting that the sparse KNARX model outperforms the Kriging and the sparse PCE-NARX in accuracy and efficiency. It is noticeable that the CPU time

is lower with higher value of N for the sparse KNARX model due to less number of Kriging model calibration.

TABLE 2 :

 2 Uncertain parameters for the duffing oscillator

	Variable Distribution type Mean Standard deviation Unit
	ω ζ ε	Uniform Uniform Uniform	2π 0.03 100	π √ 3 0.015 √ 3 10 √ 3	rad s -1 --
	A	Normal	0.6	0.06	N
	ω x	Normal	1	0.1	rad s -1

TABLE 4 :

 4 Accuracy of the surrogate models in predicting instantaneous response characteristics for the duffing oscillator

	Method	N Time instance	y(t)	R 2
	Kriging	50		337.8 × 10 -3 0.6622
	Sparse PCE-NARX 35 Sparse KNARX 21	t = 10 s	9.6 × 10 -3 8.8 × 10 -3	0.9904 0.9912
	Sparse KNARX	25		4.6 × 10 -3	0.9954
	Kriging	50		926.9 × 10 -3 0.0731
	Sparse PCE-NARX 35 Sparse KNARX 21	t = 20 s	4.4 × 10 -3 5.4 × 10 -3	0.9956 0.9946
	Sparse KNARX	25		3.1 × 10 -3	0.9969
	Kriging	50		1701.1 × 10 -3	-
	Sparse PCE-NARX 35 Sparse KNARX 21	t = 30 s	2.3 × 10 -3 3.0 × 10 -3	0.9977 0.9970
	Sparse KNARX	25		2.0 × 10 -3	0.9980

TABLE 5 :

 5 Comparison of accuracy and efficiency of the surrogate models for the duffing oscillator

	Accuracy	Efficiency

TABLE 7 :

 7 Accuracy of the surrogate models in predicting the instantaneous response characteristics for theBouc-Wen oscillator

	Displacement	Velocity

TABLE 10 :

 10 Accuracy of sparse KNARX in predicting instantaneous response characteristics for the 2-DOF dynamical system

	N Time instance	y 1 (t)	R 2
	20 50	t = 10 s	11.6 × 10 -3 0.9884 1.1 × 10 -3 0.9989
	20 50	t = 20 s	16.2 × 10 -3 0.9838 1.9 × 10 -3 0.9981
	20 50	t = 30 s	16.2 × 10 -3 0.9838 2.3 × 10 -3 0.9977

TABLE 11 :

 11 Accuracy and efficiency of sparse KNARX in computing y 1 (t) for the 2-DOF dynamical system

			Accuracy	Efficiency
	Method	¯	max(|y 1 (t)|)	R 2 max(|y 1 (t)|)

Model 2 (M 1 = 8) x (t) x (t) x (t -2∆t) x (t -2∆t) y (t -∆t) y (t -∆t) y (t -2∆t) y (t -2∆t) y 2 (t -∆t) y 2 (t -∆t) y 3 (t -∆t) y 3 (t -∆t) y 3 (t -2∆t) y 3 (t -2∆t) x (t -2∆t) y 2 (t -2∆t) x (t -2∆t) y 2 (t -2∆t) x (t -2∆t) y (t -2∆t)

of Gaussian process NARX models and their higher-order frequency response functions."

Mechanical Systems and Signal Processing, 104, 188-223. Xie, H., Tang, H., and Liao, Y.-H. (2009)