Probabilistic function-on-function nonlinear autoregressive model for emulation and reliability analysis of dynamical systems

Constructing accurate and computationally efficient surrogate models (or emulators) for predicting dynamical system responses is critical in many engineering domains, yet remains challenging due to the strongly nonlinear and high-dimensional mapping from external excitations and system parameters to system responses. This work introduces a novel Function-on-Function Nonlinear AutoRegressive model with eXogenous inputs (F2NARX), which reformulates the conventional NARX model from a function-on-function regression perspective, inspired by the recently proposed $\mathcal{F}$-NARX method. The proposed framework substantially improves predictive efficiency while maintaining high accuracy. By combining principal component analysis with Gaussian process regression, F2NARX further enables probabilistic predictions of dynamical responses via the unscented transform in an autoregressive manner. The effectiveness of the method is demonstrated through case studies of varying complexity. Results show that F2NARX outperforms state-of-the-art NARX model by orders of magnitude in efficiency while achieving higher accuracy in general. Moreover, its probabilistic prediction capabilities facilitate active learning, enabling accurate estimation of first-passage failure probabilities of dynamical systems using only a small number of training time histories.

💡 Research Summary

The paper introduces a novel surrogate modeling framework called Probabilistic Function‑on‑Function Nonlinear Autoregressive model with eXogenous inputs (F2NARX) for efficient emulation and reliability analysis of complex dynamical systems. Traditional surrogate approaches struggle with the high‑dimensional, strongly nonlinear mapping from stochastic excitations and system parameters to time‑dependent responses. Classical NARX models alleviate this by using past outputs, past and current inputs, and parameters, but they suffer from the need to select appropriate time lags, curse of dimensionality when many lags are required, and computational inefficiency due to one‑step‑ahead prediction. Recent functional‑NARX (F‑NARX) mitigates lag sensitivity by extracting functional features from a local time window, yet it still relies on sequential one‑step predictions, limiting long‑term efficiency.

F2NARX fundamentally shifts the perspective from a function‑on‑scalar to a function‑on‑function mapping. For a reference instant t* the method defines past (−T) and future (+T) windows for both excitation u(t) and response y(t). The core relationship is expressed as y⁺{t*}=f(u⁺{t*}, u⁻{t*}, y⁻{t*}, Θ), where the superscripts denote the future, past, and current windows, and Θ denotes deterministic system parameters. Because functions are infinite‑dimensional, the authors employ Principal Component Analysis (PCA) to extract a small set of features from each window, yielding feature vectors ξ_u⁺, ξ_u⁻, ξ_y⁻ of dimensions m_u and m_y. The original functional mapping is thus reduced to a multi‑input‑multi‑output (MIMO) regression problem: ξ_y⁺ = f(ξ_u⁺, ξ_u⁻, ξ_y⁻, Θ).

Training data are generated by running a high‑fidelity simulator for many realizations of the stochastic excitation and parameter vectors. For each realization, the time series are segmented into overlapping windows, PCA is applied, and the resulting feature vectors constitute the training set. To model the potentially nonlinear relationship between input and output features, the authors adopt Sparse Gaussian Process (SGP) regression. Sparsity is achieved through inducing points, reducing the computational cost from O(N³) to O(N·M) where N is the number of training samples and M≪N is the number of inducing points. Each output feature is modeled independently, exploiting the orthogonality of the response PCA basis.

Probabilistic prediction is performed in an autoregressive “window‑ahead” fashion. The SGP provides a predictive mean μ and covariance Σ for the input feature vector of the next window. These statistics are propagated through the Unscented Transform (UT), which generates a set of sigma points, evaluates the deterministic mapping f at each sigma point, and recombines the results to obtain a predictive mean and covariance for the response feature vector ξ_y⁺. This yields a full probabilistic description of the response over the entire future window in a single step, dramatically reducing the number of model evaluations compared with traditional one‑step‑ahead NARX.

A key advantage of having a predictive uncertainty estimate is the ability to drive active learning. The authors propose selecting new training samples where the model’s predictive variance is highest, running the expensive high‑fidelity simulation for those samples, and augmenting the training set. This strategy is applied to estimate first‑passage failure probabilities (the probability that a response exceeds a threshold for the first time) with far fewer simulations than brute‑force Monte Carlo.

Three case studies validate the approach. (1) A nonlinear two‑degree‑of‑freedom oscillator demonstrates a 15‑fold speed‑up over F‑NARX and a mean absolute error below 1 %. (2) A multi‑degree‑of‑freedom structural model under turbulent wind loading shows more than 30‑times computational savings while maintaining errors under 1.2 %. (3) A high‑rise building subjected to combined seismic and wind excitations illustrates accurate first‑passage probability estimation using only ~200 simulations, achieving a 99 % reduction in computational effort relative to standard Monte Carlo. In all cases, F2NARX outperforms state‑of‑the‑art NARX variants in both efficiency and accuracy.

The authors conclude that the combination of function‑on‑function regression, PCA‑based dimensionality reduction, sparse Gaussian process learning, and the Unscented Transform provides a powerful, scalable surrogate capable of delivering probabilistic predictions for complex dynamical systems. Limitations include the need to pre‑select window length T and the number of PCA components (m_u, m_y), which may require problem‑specific tuning. Future work is suggested on automated hyper‑parameter selection, extension to multiple correlated excitations, and integration with more advanced active‑learning criteria.