State Space Kriging model for emulating complex nonlinear dynamical systems under stochastic excitation
Kai Chenga Iason Papaioannoua MengZe Lyub Daniel Straub








organization=Engineering Risk Analysis Group, Technical University of Munich,
 addressline=Theresienstr. 90, 
city= Munich,
 postcode=80333, 
 country=Germany
organization=College of Civil Engineering, Tongji University,
 addressline=Siping Rd. 1239, 
city=Shanghai,
 postcode=200092, 
 country= China
 
§ ABSTRACT 

We present a new surrogate model for emulating the behavior of complex nonlinear dynamical systems with external stochastic excitation. The model represents the system dynamics in state space form through a sparse Kriging model. The resulting surrogate model is termed  state space Kriging (S2K) model. Sparsity in the Kriging model is achieved by  selecting an informative training subset from the observed time histories of the state vector and its derivative with respect to time. We propose a tailored technique for designing the training time histories of state vector and its derivative, aimed at enhancing the robustness of the S2K prediction. We validate the performance of the S2K model with various benchmarks. The results show that S2K yields accurate prediction of complex nonlinear dynamical systems under stochastic excitation with only a few training time histories of state vector.



 Stochastic dynamical system; Surrogate model; Active learning; Gaussian process; Sparse learning.








§ INTRODUCTION


Dynamical systems are widely used in modern engineering and applied science for modeling complex underlying physical phenomena <cit.>. With the increase of computational power, numerical simulation offers a feasible way to study and predict the behavior of complex dynamical systems. However, the response of complex dynamical systems is governed by uncertainties, due to the stochastic external excitations, uncertain boundary conditions, and natural variability of system properties <cit.>. To obtain effective prediction,  these uncertainties must be accounted for. To this end, uncertainty quantification of stochastic dynamical systems has gained particular interest in the last few decades <cit.>.  


In the context of uncertainty quantification, sampling-based simulation methods, i.e., Monte Carlo simulation (MCS) and various variance reduction methods <cit.>, are generally used to propagate uncertainties from system inputs to system response quantities of interest. These methods are robust, but they are infeasible when only a small number of expensive simulations are affordable or available. To address this issue, surrogate models have been widely used to construct computationally efficient approximations of the expensive computational model. Various surrogate modelling techniques have developed, including: Gaussian process regression (aka, Kriging) <cit.>, support vector regression (SVR) <cit.>, polynomial chaos expansion (PCE) <cit.>, neural networks <cit.>. These surrogate modelling techniques are powerful for approximating the behavior of traditional static “black-box”models, but exhibit difficulties when applied to dynamical systems under stochastic excitation <cit.>. In these problems, the number of input parameters due to discretization of the stochastic external excitations, often represented in terms of a white noise process, can be extremely high. Most common surrogate models suffer from the “curse of dimensionality”, since effective learning typically requires at least twice as many samples as the number of input parameters.
To address this issue, a standard technique is to insert existing surrogate modelling techniques into the framework of nonlinear auto-regressive with exogenous input (NARX) modelling <cit.>. The NARX model is a powerful system identification technique, which is established based on the principle of causality, i.e., it assumes that the system response quantity of interest at the current time instant is only affected by its previous several response values and the  current and past several values of external excitations. Based on this cause-consequence effect, the current response quantity of interest is assumed to be a function of the response values at past multiple time instants and the input excitation at the current and previous multiple instants. This function can be emulated through application of  existing surrogate modelling techniques.  Although the NARX model has proved its effectiveness in several structural dynamical problems <cit.>, it is only an empirical model, and there is no rigorous way to select the model hyper-parameters, e.g., the time lags of both input and output. In addition, the NARX model struggles to emulate the response of complex highly nonlinear dynamical problems. Recently, a manifold NARX (mNARX) model <cit.>  has been proposed to address the above problems, in which a transformation function is adopted to map the input into a problem-aware manifold, which is expected to be more suitable for constructing the NARX model than in original input space. However, physical information of the dynamical system under investigation is required to construct such input manifold in mNARX. If no physical information is available, the mNARX degenerates to the traditional NARX. 
 
In the present work, we propose to emulate complex nonlinear dynamical systems under stochastic excitation in their state space form with the Kriging model, termed as the state space Kriging (S2K) model. The state space representation of a stochastic dynamical system can be considered as a multi-input multi-output (MIMO) function, in which the input is the state vector and the external excitation, and the output is the derivative of the state vector. The Kriging model is utilized to learn every component of this MIMO function separately. To overcome the inefficiency of the Kriging surrogate for large training data set, we introduce an active learning algorithm to select a reduced training set from the training time histories, resulting in a sparse Kriging model. By learning the state space form of a dynamical system, the S2K model avoids the curse of dimensionality 
resulting from the discretization of the stochastic external excitation. Numerical examples demonstrate that the S2K model is effective for emulating various complex nonlinear dynamical systems with stochastic excitation with only a few training time histories of the state quantities, and outperforms the NARX model.

The layout of this paper is as follows. In section 2, we review the fundamental definitions of nonlinear dynamical systems under stochastic excitation and the basic idea of the NARX model. The S2K model is presented in Section 3, together with the active learning algorithm and the technique for designing the training time history of state vector. In Section 4, several benchmarks are used to assess the performance of our method. The paper concludes with final remarks in Section 5.



§ BACKGROUND



A dynamical system is generally represented by a high-order differential equation. It can be equivalently transformed into multiple coupled first-order equation systems by introducing new state variables, a form known as the state space representation. 

    dx/dt   = f_1(x) 
    d^2 x/dt   = f_2 (x,dx/dt)
       ⋮
    d^n x/dt   = f_n(x,dx/dt,...,d^n-1x/dt^n-1)

 In the present work, we consider a general nonlinear dynamical system under stochastic excitation in its state space form <cit.>, which can be generally expressed as
 
    (t)  = ((t),(t)),  with(0) = _0,

 where (t)=[X_1(t),...,X_n(t)]^ T∈R^n is the state vector at time t; _0 is the initial condition; (t) is the derivative of (t) with respect to t; (t)=[U_1(t),...,U_m(t)]^ T∈R^m is the external stochastic excitation vector acting on the structure; (·) is the n-dimensional nonlinear vector function . 
 
 For every realization of the time history (t) of (t), the time history of the state vector (t) of Eq. (<ref>) exists and uniquely depends on the initial conditions _0, and one can apply various numerical discretization methods, e.g., the Runge-Kutta method, to find the solutions (t_i)(i=1,...,N_t) at N_t time steps within the time period of interest [0,T]. In general, (·) is a computational expensive  "black-box" model, which makes UQ of complex  dynamical systems under stochastic excitation a computationally demanding task.  

  The NARX model is a popular method for UQ of dynamical system under stochastic excitation <cit.>. It is a system identification technique developed based on the discrete-time step representation of the dynamical system. Given a discrete time history of the state vector (t_1),...,(t_N_t) corresponding to a discrete realization of the stochastic excitation (t), namely, (t_1),...,(t_N_t), the NARX model represents the response quantity of interest  x(t_i)∈R  at the current time instant as a function of its past values and the input excitation values at the current and previous instants as

    x(t_i)= g((t_i),(t_i-1),...,(t_i-n_u),x(t_i-1),...,x(t_i-n_x)) + ϵ_t,

where n_u and n_x represent the maximum excitation and response time lags; ϵ_t is the residual of the NARX model; g(·) is the underlying model to be learned. For dynamical systems, the function g(·) is usually learned by polynomial model or Kriging model <cit.>.  

The NARX model can be interpreted as an empirical model that tries to capture the behaviour of the stochastic dynamical system with a low-dimensional auto-regressive surrogate model. Its accuracy highly depends on the choice of the time lags n_u and n_x, and there is no rigorous way to determine them. For problems with long memory, both n_u and n_x become large, and the NARX model is high-dimensional. Moreover, NARX gives poor prediction for strongly nonlinear dynamical systems, which we also confirm in Section <ref>. 




§ METHODOLOGY


 In this section, we introduce the proposed algorithm for learning  dynamical systems under stochastic excitation in state space form. We also present the active learning algorithm for selecting the informative training set and the technique for designing the training time history of state vector and its derivative. 



 §.§ Kriging model for learning dynamical systems 


By denoting (t)=(t), the original dynamical system in Eq. (<ref>) can be expressed as
 
    (t)  = ((t),(t))
    .


    Note that (·): R^m+n→R^n is a deterministic function, which maps the state vector (t)  ∈R^n and the external excitation (t)  ∈R^m to (t) ∈R^n at time instant t. The input dimension of (·) depends on the dimension of the state vector n of the associated dynamical system and the cardinality of the external excitation vector m. For every realization of the time history (t) of the excitation (t), the time history of the corresponding state vector (t) and its derivative  (t) over the time period of interest [0,T] can be estimated numerically. 

In this work, we employ the Kriging model <cit.> to approximate the state space representation (·). Note that other existing surrogate modelling techniques, such as PCE, SVR and neural network could also be implemented in this context. 
As in the ordinary Kriging model, it is assumed that every component y_i(t) of (t) is a  Gaussian process with constant mean β_i:

    y_i(t) = β_i + σ_i^2 Z((t)),   for i =1,...,n,

where (t)=[(t),(t)]^
T∈R^n+m, σ_i^2 is the process variance, and Z((t)) is a stationary Gaussian process with zero mean and unit variance. The correlation function of Z((t)) is given by
R((t),(t');_i), which describes the correlation between (t) and (t') at time t and t', with the correlation lengths in various coordinate directions being controlled by the hyper-parameter vector _i∈R^m+n. Moreover, we assume that y_i(t) is independent of y_j(t) for all i ≠ j, which allows learning (t) component-wisely.
   
Given the experimental design _t=[(t_1),...,(t_N)]∈R^(n+m)× N, the joint distribution y_i(t) and  _i=[y_i(t_1),...,y_i(t_N)]^
T∈R^N corresponding to _t is multivariate Gaussian given by

    [[ y_i(t);     _i ]] ∼𝒩([[   β_i; β_i F ]],σ_i^2[[         1 r^ T((t));    r((t))         R ]] ),

where =[1,...,1]^ T∈R^N and

    r((t))   := [R((t_i_1),(t);_i)]_i_1∈R^N, 
    R   :=[R((t_i_1),(t_i_2);_i)]_i_1,i_2∈R^N× N.



The predictive distribution of y_i(t) given the observation _i=_i is still Gaussian, i.e.,

    ŷ_i(t)  ∼𝒩(μ_i(t), s_i^2(t)),

where the predictive mean and predictive variance are given by


    μ_i(t) = β_i +r^T((t))R^-1(_i-β_i)  ,

     and
     
    s_i^2(t) = σ_i^2 (1-r^T((t))R^-1r((t))+(1-^TR^-1r((t)))^2/(^TR^-1)).


The predictive mean is the kriging surrogate model prediction, and the predictive variance is used for measuring the predictive uncertainty. 
The prediction (t) of the vector (t) is obtained by stacking ŷ_i(t)(i=1,...,n), namely, (t)=[ŷ_1(t),...,ŷ_n(t)]^ T.  (t) is a Gaussian vector, with mean m(t) = [m_1(t),..., m_n(t)]^ T and diagonal covariance matrix with diagonal vector s^2(t) = [s_1^2(t),...,s_n^2(t)]^ T, namely

    (t) = {[ ŷ_1(t);      ⋮; ŷ_n(t) ]}∼𝒩({[ μ_1(t);      ⋮; μ_n(t) ]}, [[ s_1^2(t)        ⋯        0;        ⋮        ⋱        ⋮;        0        ⋯ s_n^2(t) ]] ).


In the current work, the Matérn-5/2 correlation function is used, which is given by

    R((t),(t');_i) = 
     (1+√(5)r +5/3r^2)  exp(-√(5)r),

where r=√(((t) -(t'))^ TΛ((t) -(t'))), and Λ is a diagonal matrix with diagonal element _i.

The hyper-parameter β_i, σ_i^2 and _i for i-th surrogate model are determined by the maximum likelihood estimation method <cit.>. The likelihood function of the hyper-parameters given data  reads

    f(β_i,σ_i^2,_i)=[ detR (_i)]^-0.5/√((2πσ_i^2)^N) exp(-1/2σ_i^2(_i-β_i)^ TR^-1(_i)(_i-β_i)).

     
    After taking the logarithm of Eq. (<ref>), the optimal values of β_i and σ_i^2, depending on hyper-parameter , can be derived analytically. This yields
     
    β_i(_i)    =    (F^TR^-1(_i)F)^-1F^TR^-1(_i)_i,
         
    σ_i^2(_i)    =   1/N(_i-β_iF)^ TR^-1(_i)(_i-β_iF).

    
 There is no closed-form solution for the optimal hyper-parameter _i, and one has to use numerical optimization algorithms to determine its value. Substituting σ_i^2(_i) in Eq. (<ref>) into the likelihood function in Eq.  (<ref>), we are left with minimizing the following _i-dependent objective function
    
    ℓ(_i)=Nlnσ_i^2(_i) + ln[ detR(_i)].


 The objective function in Eq. (<ref>) is generally highly nonlinear and multimodal <cit.>, and thus a global optimization algorithm is required to find the hyper-parameter. In the present work, we use the multi-starts gradient-free “Hooke & Jeeves” pattern
search method <cit.> for determining the optimal hyper-parameter vector. Note that the single-start “Hooke & Jeeves” pattern
search method has been implemented in the popular DACE toolbox <cit.>, and here we adapt it to multi-starts case for improving the robustness of Kriging model. 
 


 §.§ Prediction with S2K model


Once the S2K model is trained, it can be used to predict the time history of system response (t) for arbitrary realization of the time history  (t) of (t) by solving the following approximate dynamical system in state space form

    (t) = ((t),(t))
     ,

where (·) is the Kriging approximation of (·).

In general, the Runge-Kutta method can be used to estimate the time history of the state vector (t_i)(i=1,...,N_t) at N_t discretized time steps within the time period of interest [0,T]. 
At each time step, the Kriging model provides a  Gaussian prediction (t) given a realization of the current predicted state vector (t) and external excitation (t). However, the whole time-history of (t) and (t) over the time period of interest is non-Gaussian. This is due to the recursive nature of the dynamical system, i.e., the output quantities at current time instant will be used as input to predict the output quantities at the next time step, and carry their predictive uncertainty with them. Since the surrogate model (·) is nonlinear, the predictive distribution of the entire time history is generally 
 intractable. In practice, one can adopt MCS to approximate the predictive distribution of the state vector over the time period of interest given a realization of the excitation. This can be achieved by emulating the time history of state vector N_MC times based on the trained S2K model, and use these time history trajectories to estimate the empirical predictive distribution.

In Fig. <ref>, we show an example of the probabilistic prediction results of the time history of the displacement of the nonlinear Bouc-Wen hysteretic oscillator (see Section <ref>), where MCS with N_MC=50 randomly predicted trajectories of the displacement time histories are used. It is found that the 95% confidence interval 





 §.§ Sparse Kriging with active learning


We now analyze the computational cost of the S2K model.  When training a Kriging model with N training samples, one needs to determine the hyper-parameter _i by minimizing the likelihood function in Eq. (<ref>). To  numerically evaluate this function, one needs to invert the correlation matrix  R∈^N× N. When a Cholesky factorization is used to decompose R, the corresponding computational costs are 𝒪(N^3). The computational costs for training the S2K model are therefore 𝒪(nN^3), which tends to be very time-consuming for large number of samples N and large n.

 
 In the Kriging model, the computational cost for making a single prediction is 𝒪(N) for the mean in Eq. (<ref>), and 𝒪(N^2) for the variance in Eq. (<ref>). 
To estimate the predictive distribution of the entire time history of the state vector, one needs to predict every component of the n-dimensional state vector (t) at N_t time steps,  N_MC times for every realization of the time history (t) of (t). The associated operations are 𝒪(nNN_tN_MC) for the mean and 𝒪(nN^2N_tN_MC) for the variance. In general, the sample size N, the time history length N_t and the MCS sample size N_MC are all quite large, which makes it time-consuming to train the S2K model and to  estimate the predictive distribution of the entire time history of the state vector corresponding to a single realization of the time history (t). 

To improve the training and prediction efficiency, we suggest constructing a sparse Kriging <cit.> model by selecting a training subset of size N_s(N_s ≪ N) from the available training set for every component of (t). This selection is performed adaptively, through first training a Kriging model with a few initial training samples uniformly selected from the whole sample set, and then updating the Kriging model by enriching informative samples selected sequentially through maximization of the following mean square error criterion

    L_i(t) = (m_i(t) - y_i(t))^2 + s_i^2(t),

where the first term represents the square of the predictive bias, and the second term denotes the predictive variance. The bias term prefers points with large local predictive error, and the predictive variance term tends to select points with large global predictive uncertainty. The acquisition function in Eq. (<ref>) therefore accounts for both the bias for local exploitation and the variance for global exploration. 

The proposed method selects one sample in each active learning step until the maximum of L_i(t)/σ^2_y is less than a threshold δ_i (e.g., 10^-5), where σ^2_y is the variance of observed model response in the training set.
 The convergence threshold δ_i is a user-specified parameter, depending on the complexity of the dynamical system to be learned. With a small δ_i, many samples will be selected to construct the S2K model, which improves the prediction accuracy, but decreases the prediction efficiency; with a large δ_i, only a small portion of the training set will be selected to construct the S2K model, which leads to a very sparse model, but the accuracy may be low. 
Based on our numerical experiments, we suggest setting δ_i∈ [10^-3, 10^-7]. A higher threshold values should be chosen when the nonlinearity of f_i(·) is expected to be high, and vice versa. 



 §.§ Design of training time history


The Kriging model is powerful for interpolation, but its performance degenerates when extrapolating away from the training data. To ensure the accuracy of the S2K model, the training data should cover the whole parameter space of both the state vector and the external excitation. However, for a black-box model, the parameter space of state vector is unknown beforehand, since one cannot know the minimum and maximum values of the time histories of the state vector associated with future realizations of external excitation.

To address this issue, one can collect many observed time histories of state vector corresponding to different realizations of external excitation, and select a subset of the time histories (e.g., requiring that the maximum value exceeds a pre-defined threshold) to train the surrogate model, as has been proposed for training the NARX model <cit.>. However, this is infeasible when only a few time histories of the state vector are available. In the current work, we suggest collecting the time histories of a pseudo state vector resulting from pseudo external excitation with magnified variability. By magnifying the variability of the excitation, the state vector will exhibit stronger variability, hence, it is expected that a few pseudo time histories of state vector will sufficiently populate its parameter space under the true excitation. To this end, one can magnify the standard deviation of the stochastic excitation at every time instant by a magnification factor σ . 

In Fig. <ref>, a time history of state variables (displacement and velocity) of the Bouc-Wen hysteretic oscillator (see Section <ref>) with white noise ground acceleration is presented, where the state variables are obtained with both the true excitation and a pseudo excitation with magnified variability (σ=2). It is found that the peaks and the valleys of the pseudo responses (both displacement and velocity) largely exceed the ones of the original response. Consequently, the pseudo time history is more likely to cover the whole parameter space of the state vector under the original excitation. By training the Kriging model with the pseudo training time history, extrapolation in the prediction stage can be potentially avoided, thereby improving the accuracy and robustness of the S2K model. 
Note that the magnification factor σ is a problem-dependent parameter. A small σ cannot guarantee full coverage of the parameter space of the state vector, while a large σ will lead to a less accurate Kriging model, since many samples will be enriched during the active learning procedure to explore the  unnecessary region.  Numerical investigation reported in Section <ref>  show that σ∈ [1.5, 2] is a good choice when only one training time history of the state vector is available. However, when one can afford multiple training time histories, it is beneficial to set different σ values for each training time history. In this work, we suggest setting

    σ_k=1+k-1/n_t-1, 
    k=1,...,n_t

 where n_t(n_t≥ 2) is the number training time histories. 
 In doing so, the mixture of different training time histories of state vector should better cover the entire parameter space, leading to more accurate prediction, as demonstrated in Section <ref>. 
For problem with white noise excitation with intensity S discretized from 0 to T with step length t, namely, f(t_i)=√(2π S/ t)ζ(t_i),   ζ(t_i)∼𝒩(0,1), one can magnify the variability of the excitation by sample from ζ(t_i)∼𝒩(0,σ^2), where σ is the magnification factor.




 §.§ Summary of the S2K algorithm


The S2K algorithm for emulating the response of dynamical systems under stochastic excitation is summarized in Algorithm <ref>.


 






 §.§ Multi-degrees of freedom dynamical system

For problems with many degrees-of-freedom, the state vector is high-dimensional, which hinders the application of Kriging model for learning the state space representation of the dynamical system. To address this issue, we first project the original dynamical system into a latent low-dimensional space identified with proper orthogonal decomposition (POD), and then learn the dynamics in the latent space. 

Denoting the POD basis as , the original state vector can then be expressed as _r(t)=^T(t), and the original dynamical system can be reformulated as

    _r(t)  = ^T(_r(t),(t)),   _r(0) = ^T_0,








§ NUMERICAL INVESTIGATIONS


We investigate the effectiveness of S2K model on four example applications. To assess the accuracy of the S2K model, the relative error  of a whole time history of the state vector corresponding to a realization of the stochastic excitation is used to define the accuracy metric as
 
    ϵ_i =∑_j=1^N_t(x_i(t_j) - x̂_i(t_j) )^2/∑_j=1^N_t(x_i(t_j) - x̅_i )^2,

where x_i(t_j) and x̂_i(t_j) are true and predicted i-th state quantity of interest at time instant t_j, and x̅_i is the mean over the whole time history. In addition, we use the mean value ϵ̅_i of ϵ_i 
corresponding to N_MC different realizations of the excitation to measure the global accuracy of i-th state, namely
 
    ϵ̅_i = 1/N_MC∑_j=1^N_MCϵ_i^(j),

where ϵ_i^(j) is the relative error of i-th state quantity of interest corresponding to j-th realization of the stochastic excitation. We set N_MC=1000.


For the first three example applications, we compare the performance of the S2K model with the NARX model in Eq. (<ref>). In the NARX model, the polynomials are used as bases, and the least angle regression technique <cit.> is used to select the most relevant basis terms. In addition, the maximal time lags of both excitation and response n_u and n_x are chosen equal to twice
the number of degrees of freedom of the considered dynamical system <cit.>. In these examples, both S2K and NARX model are trained 10 times independently, and the relative errors are depicted by boxplots, whereby the central mark indicates the median, and the bottom and
top edges of the box indicate the 25th and 75th percentiles, respectively. 

In both NARX and S2K model, the training time history of the state vector (t) and its derivative are evaluated at equidistant time instants t_i=i t(i=1,...,T/ t) over the time period of interest [0,T] by the Matlab solver ode89, where t is the step length. Note that while equidistant discretization of the training time history is required in the NARX model, it is not essential for the S2K model. However, we apply it here for consistency.
 




 §.§ Quarter car model


 We first consider a quarter car model represented by a
nonlinear two degree-of-freedom system <cit.> depicted in Fig. <ref>. The displacements of
the masses are governed by the following system of ordinary differential equations (ODEs), namely,

    m_s x_1(t)    = -k_s(x_1(t)-x_2(t))^3 - c(ẋ_1(t)-ẋ_2(t)), 
      m_u x_2(t)    = k_s(x_1(t)-x_2(t))^3 - c(ẋ_1(t)-ẋ_2(t)) + k_u(u(t)-x_2(t)),

where the sprung mass m_s=22.7   kg and the unsprung mass m_u=42   kg are connected by a nonlinear spring of stiffness k_s = 1897.02   N/m^3 and a linear damper with damping coefficient c=601.8   N· s/m. An external excitation is applied to m_u through a linear spring of stiffness k_u = 1771.4   N/m. x_1(t) and x_2(t) are the displacements of m_s and m_u, respectively. In this work, the excitation is modeled by a stochastic process as u(t) = A  sin(b t), where the parameters A, b follow uniform distributions, A ∼ U(0.09,0.11) (m) and b ∼ U(1.8π,2.2π) (rad/s).


Denoting (t)=[x_1(t),ẋ_1(t),x_2(t),ẋ_2(t)]^ T and (t)=(t), one can express the ODEs in Eq. (<ref>) in state space form as

    y_1(t)      =ẋ_1(t), 
    y_2(t)     = -k_s/m_s(x_1(t)-x_2(t))^3 - c/m_s(ẋ_1(t)-ẋ_2(t)),  
    y_3(t)      =ẋ_2(t), 
    y_4(t)     = k_s/m_u(x_1(t)-x_2(t))^3 - c/m_u(ẋ_1(t)-ẋ_2(t)) + k_u/m_u(u(t)-x_2(t)).





We use the S2K model to construct a surrogate model (·):  R^5→R^4 for the state space representation in Eq. (<ref>), in which the input is (t) = [(t), u(t)]^ T∈R^5, and the output is (t) ∈R^4. 
In this example, we draw a time history of length 10s of the excitation randomly (A=0.0964,b=6.5248), and the time history of the state vector and its derivative are estimated with time step t = 0.002s.  We use this time history of the state vector and the corresponding derivative to construct the four Kriging models, in which 5 equidistant samples selected from the whole time history are used to train the initial Kriging models, and the active learning algorithm introduced in Section <ref> is utilized to select the informative samples sequentially with the convergence threshold δ=10^-6. The corresponding results are depicted in Fig. <ref>, where a test time history is also presented. The results demonstrate that the S2K model is highly accurate. Indeed, the relative errors of the four Kriging models corresponding to the four state quantities of the test time histories depicted in Fig. <ref> are ϵ_1=6.09× 10^ -4, ϵ_2=4.54× 10^ -4, ϵ_3=6.02× 10^ -4, and ϵ_4=4.50× 10^ -4, respectively. In addition, it is observed that only a small portion of the samples are selected from the whole time history (a time history consists 5001 discrete time instants in total) after the active learning procedure. Specifically, the final sample sizes of the four Kriging models are 11, 24, 7, and 23, respectively. That is, the corresponding degrees of sparsity are 0.22%, 0.48%, 0.14%, and 0.46%.

Fig. <ref> presents the average relative error ϵ̅_1 of x_1(t) predicted by the S2K model with varying the number of time histories and various convergence thresholds. It shows that the S2K model with only one training time history of the state vector already yields very accurate results. As expected, the accuracy can be further improved by increasing the number of the training time histories as well as by decreasing the convergence threshold. In addition, the training sample sizes of the four Kriging models corresponding to the four components of the state space representation in Eq. (<ref>) are shown in Fig. <ref>. It is observed that only a limited number of samples are retained in the four Kriging models for various settings, which confirms the sparsity of the S2K model. 

The comparisons of relative error between the S2K (δ=10^-5) and the NARX model (maximum polynomial order set to 3 <cit.>) is shown 
in Fig. <ref>. It shows that both methods provide very accurate prediction, but the S2K outperforms the NARX model, especially when more training time histories are available. 
















 §.§  Duffing oscillator


The second example is a Duffing oscillator <cit.> subjected to random loading for a duration of T=10s. The oscillator is governed by the following ODE

    x(t) + 2ζω_nẋ(t) + ω^2_nx(t) +β x^3(t) = u(t),

where ζ=0.05,  ω_n=10 rad/s, and β=2000 m^-2s^-2. 
The random loading u(t) is modeled by a white noise ground acceleration discretized
in the frequency domain as <cit.>

    u(t) = √(2Sω)∑_i=1^d/2 [ϑ_i cos(ω_it) + ϑ_d/2+isin(ω_it) ],

where ω_i=i ω with ω=30π/d, and ϑ_i∼𝒩(0,1), for i=1,...,d. In the current work, the spectral intensity of the Gaussian white noise is set to S=0.1 m^2/s^3, and we choose d=150 <cit.> . 

Denoting  (t)=[x(t),ẋ(t)]^ T and (t)=(t), the original differential equation in Eq. (<ref>) can be expressed by its state space representation as

    y_1(t)      =ẋ(t), 
    y_2(t)     =  u(t) -  2ζω_nẋ(t) - ω^2_nx(t) -β x^3(t).


The S2K model is used to learn the dynamics of the Duffing oscillator in Eq. (<ref>). In this example, we use only one time history of the state vector estimated for the total duration T = 10 s with time step t = 0.002 s to train the two Kriging models separately.
The results are depicted in Fig. <ref>, where a test time history is also presented. One can see that the S2K model is highly accurate in emulating the Duffing model, with relative errors ϵ_1=9.68× 10^-6 for x(t) and ϵ_2=9.09× 10^-6 for ẋ(t). Again, only a few informative samples are enriched during the active learning procedure, which leads to highly sparse Kriging models. 

The average relative error ϵ̅_1 of x(t) predicted by the S2K model by varying the training time histories with different variability magnification factors of excitation are presented in Fig. <ref>, in which σ∈ [1,2] signifies that different magnification factors are assigned to different training time histories according to Eq. (<ref>). The variability of the excitation is magnified by amplifying the standard deviation of ϑ_i(i=1,...,150) in Eq. (<ref>).
The S2K model with only one training time history of the state vector already provides prediction  with less than 1% relative error, and its accuracy can be improved by about two orders of magnitude by magnifying the variability of the excitation by σ=1.5 or σ=2. The mixture of 3 training time histories of state vector with different σ values yields the optimal results for n=3.

The training sample sizes of the two Kriging models corresponding to the two components of the state space representation in Eq. (<ref>) after the active learning procedure are depicted in Fig. <ref>. As expected, more samples are selected to train the Kriging model when using a larger magnification factor σ. However, both Kriging models remain  sparse for all settings. 

In Fig. <ref>, we compare the relative errors of x(t) obtained with the S2K model (δ=10^-6,σ=1.5) and the NARX model (maximum polynomial order set as 3 <cit.>) by varying the number of training histories. Due to the strong variability of the state quantity x(t) over the time domain, the NARX model gives poor prediction. 
















 §.§ Nonlinear Bouc-Wen hysteretic oscillator


In the third example, we consider a Bouc-Wen hysteretic oscillator under random external excitation <cit.>, described by the following differential equation

    m x(t) + cẋ(t) + k [α x(t) +(1-α)x_y z(t)] = m u(t), 
    
    ż(t) = 1/x_y [Aẋ(t) -β |ẋ(t)|| z(t)|^d-1z(t) -γẋ(t) |z(t)|^d ],

where the mass, stiffness and damping of the oscillator are 
m=6×10^4kg, k=5×10^6N and c = 2mζ√(k/m), respectively, with ζ=0.05. The degree of
hysteresis is defined by α, which is chosen as 0.5. In addition, we set x_y=0.04m, β =γ = 0.5 and A=1,d=3.


Denoting  (t)=[x(t),ẋ(t),z(t)]^ T and (t)=(t), the above differential equation of Bouc-Wen oscillator can be expressed by its state space form as 

    y_1(t)      =ẋ(t), 
    y_2(t)     = -c/mẋ(t) -k/m [α x(t) +(1-α)x_y z(t)] + u(t),  
    y_3(t)     = 1/x_y [Aẋ(t) -β |ẋ(t)|| z(t)|^n-1z(t) -γẋ(t) |z(t)|^n ].


The excitation is modeled by a white noise ground acceleration discretized in frequency domain as in Eq. (<ref>) with spectral intensity being S=0.05 m^2/s^3. The time period is set to T=8 s <cit.>. 

We first draw one time history of excitation, and the state vector is estimated with  time step t=0.002 s. The corresponding excitation is depicted in Fig. <ref>, together with the nonlinear response of the auxiliary variable z of the Bouc-Wen model. We use this time history to train the S2K model with  convergence threshold δ=10^-4. The results are depicted in Fig. <ref>, in which a test time history of state vector is also presented. Again, it is shown that the S2K model yields accurate prediction, with relative errors corresponding to the test time history of state vector depicted in Fig. <ref> of  ϵ_1=1.44× 10^-4 for x(t), ϵ_2=1.65× 10^-5 for ẋ(t) and ϵ_3=2.47× 10^-5 for z(t). Since the third state equation in Eq. (<ref>) is highly nonlinear, 173 samples (5 initial samples plus 168 enriched samples) are selected from the whole time history of z(t) to construct the third Kriging model. 







The average relative error ϵ̅_1 of x(t) predicted by the S2K model by varying the training time histories with different 
magnification factors is presented in Fig. <ref>, in which σ∈ [1, 2]  signifies that different magnification factors are assigned to different training time histories
according to Eq. (<ref>). One can see that the S2K model with only one training time history of the state vector already yields accurate predictions. Its accuracy can be improved by about one order of magnitude by magnifying the variability of the excitation by σ=1.5 times, but there is no additional improvement if a larger magnification factor σ=2 is used. In addition, the accuracy can be further improved by collecting more time histories of the state quantities. Again, one can see that the mixture of the 3 training time histories of state vector with different σ values yields the most accurate S2K model when n=3.

The training sample size of the three Kriging models corresponding to the three components of the state space representation in Eq. (<ref>) after the active learning process is depicted in Fig. <ref>. It shows the first two Kriging models remain sparse for various parameter settings while more samples are selected to train the third Kriging model. 



 


The comparison of the relative errors of x(t) obtained with both the S2K model (δ=10^-4,σ=1.5) and the NARX model (maximum polynomial order 5) is depicted in Fig. <ref>.
Note that different from <cit.>, we do not use any physical information to construct the NARX model here. For this challenging nonlinear problem, the vanilla NARX model with polynomial basis fails to provide satisfactory prediction. By contrast, the S2K model yields accurate prediction with relative errors that are 4 orders of magnitude lower than that of the NARX model. 





 §.§ Two-story nonlinear hysteretic structure 



 We apply the S2K model to emulate the behavior of a two-story two-span nonlinear hysteretic frame structure subjected to ground motion excitation <cit.>. The corresponding motion equation reads

    M(t) + C(t) +G[(t),(t)] = -MI u(t),

where (t)∈R^2,(t)∈R^2, and (t)∈R^2 represent the displacement, velocity and acceleration vector; M∈R^2×2,  C∈R^2×2 and I∈R^2×2 are the lumped mass, Rayleigh damping and the identity matrices; u(t) is modeled by a  white noise ground acceleration process as that in Eq. (<ref>); G(·) is the restoring force vector characterized by the Bouc-Wen model as 

    G̃_j((t),(t)) = α K_j x̃_j(t) + (1-α)K_j z_j(t),

in which α = 0.04; G̃_j(·), x̃_j(t), and K_j represent the j-th inter-story restoring force, drift, and story initial stiffness, respectively; z_j(t)(j = 1,2) is the j-th auxiliary hysteretic displacement, which is described by

    ż_j(t) = ϖ/1+d_ηϵ_j(t)(1-ξ exp(-{z_j(t) sgn[ẋ̃̇_j(t)]-q/(β+γ)κ/[ψ + d_ψϵ_j(t) ](λ +ζ_sξ )}^2))
    ,

where κ =1+d_νϵ_j(t), ϖ = ẋ̃̇_j(t) - κ[β |ẋ̃̇_j(t)|z_j(t) + γẋ̃̇_j(t) |z_j(t)| ], ξ= ζ_s[1-e^-pϵ_j(t)]; the parameters β and γ control the basic hysteresis shape; d_ν and d_η are the strength and stiffness degradation; ζ_s measures the total slip; q,p,ψ,d_ψ, and λ are the initiation, slope, magnitude, rate, and severity interaction of pinching <cit.>; 
ϵ_j(t) is the j-th story hysteretic-dissipated energy, which is described by

    ϵ̇_j(t) = ẋ̃̇_j(t)z_j(t), for   j=1,2.


In the current work, the system parameters are set as follows: lumped mass M_j=2.6×10^5   kg, initial stiffness K_j = 10^8   N/m, β=15  m^-1, γ=150  m^-1,  d_ν=d_η=p=1000  m^-2, q=0.25, d_ψ= 5   m^-2, λ=0.5, ζ_s=0.99, ψ=0.05  m. In addition, the damping ratios of the first two modes are taken as ζ_1=ζ_2=0.05.
 
The state vector (t)=[x_1(t),x_2(t),ẋ_1(t),ẋ_2(t),z_1(t), z_2(t), ϵ_1(t), ϵ_2(t)]^ T is 8-dimensional, and we construct 8 Kriging models in the S2K model to learn the state space representation of this dynamical system, where the input is [(t)^ T, u(t)]^ T, and the output is (t)=(t). The time history of the state vector over the whole time period of interest T∈ [0,10] s is evaluated with time step t = 0.01 s. The average relative errors ϵ̅_1 of x_1(t) predicted by the S2K model with varying number of training time histories are presented in Fig. <ref>, in which  different magnification factors are assigned to different training time histories
according to Eq. (<ref>). Due to the strong non-linearity and non-smoothness property of this model, at least 5 time histories of the state vector are required to train an accurate S2K model. To illustrate its performance, a specific test time history of the state vector is depicted in Fig. <ref>. It shows that the eight state quantities are well predicted by the S2K model over the whole time period of interest (with 5 training time histories).  In this example, only the last four components of the state space representation are strongly nonlinear, especially  y_6(t)=ż_1(t) and y_7(t)=ż_2(t) in Eq. (<ref>).  We therefore only present the training sample sizes of the last four Kriging models in Fig. <ref>. As expected, as the number of training time histories increases, an increasing number of samples is selected during the active learning procedure to train the Kriging models, and more than 20 % of the total samples from the 5 training time histories are selected to train the two challenging functions: y_6(t)=ż_1(t) and y_7(t)=ż_2(t) in Eq. (<ref>). In this example, application of the NARX model led to meaningless predictions after some time period, it is not suitable for emulating the response of this system due to the high nonlinearity in Eq. (<ref>). Hence we do not report the relative error of the NARX model for comparisons. 




 


 



§ CONCLUDING REMARKS



 In this work, we have presented a novel surrogate modeling framework for emulating dynamical systems under stochastic excitation by learning its state space representation with Kriging (S2K) model. Several conclusions can be drawn from the work:
 
 (1) Learning the state space representation of a dynamical system under stochastic excitation can avoid the curse of dimensionality of surrogate model training due to discretization of the stochastic excitation. 

 (2) The proposed active learning algorithm can select an informative sample subset from the whole training sample set efficiently, resulting in a sparse Kriging model. 
 
 (3) The proposed technique for designing the training time history of the state vector by magnifying the variability of the excitation can improve the accuracy of S2K model. 
 
 (4) Numerical examples demonstrate that the S2K model is powerful for emulating various complex nonlinear dynamical systems under stochastic excitation. It provides highly accurate prediction (Relative error is less than 10^-3) with only a few time histories of the state vector, and it outperforms the NARX model in terms of accuracy and efficiency. 
 
In the future, we plan to adapt this framework to emulating complex nonlinear dynamical systems with both random system parameters and stochastic excitation. Moreover, we will also consider combining the S2K model with model reduction techniques, e.g., proper orthogonal decomposition or auto-encoder, to emulate stochastic dynamical systems with large number of degrees of freedom. 

 


§ ACKNOWLEDGEMENTS

 This work was supported by the Alexander von Humboldt Foundation.


plain



