## Abstract

Mobility prediction of off-road autonomous ground vehicles (AGV) in uncertain environments is essential for their model-based mission planning, especially in the early design stage. While surrogate modeling methods have been developed to overcome the computational challenge in simulation-based mobility prediction, it is very challenging for a single surrogate model to accurately capture the complicated vehicle dynamics. With a focus on vertical acceleration of an AGV under off-road conditions, this article proposes a surrogate modeling approach for AGV mobility prediction using a dynamic ensemble of nonlinear autoregressive models with exogenous inputs (NARX) over time. Synthetic vehicle mobility data of an AGV are first collected using a limited number of high-fidelity simulations. The data are then partitioned into different segments using a variational Gaussian mixture model to represent different vehicle dynamic behaviors. Based on the partitioned data, multiple surrogate models are constructed under the NARX framework with different numbers of lags. The NARX models are then assembled together dynamically over time to predict the mobility of the AGV under new conditions. A case study demonstrates the advantages of the proposed method over the classical NARX models for AGV mobility prediction.

## 1 Introduction

As an emerging technique, autonomous ground vehicles (AGVs) have gained much attention in different applications in recent years, since they can increase traffic efficiency, reduce the probability of accidents, and replace human in harsh/boring working environments [1,2]. As one of the most important applications of AGVs, the off-road AGV replaces among others, the human in rough terrain conditions and battlefields to transport materials autonomously or explore unknown territory for knowledge discovery, such as the exploration of Mars [3]. In addition, it plays an important role in relieving human effort in the agricultural field or in reducing injuries on the battlefield [4].

To support the design, development, and operation of off-road AGVs, modeling and simulation (M&S) of vehicle mobility is very important [5]. An effective M&S tool can not only predict the performance of a vehicle operating in a certain environment but also can provide insights for decision makers to design the vehicle and plan the mission [6]. Motivated by predicting vehicle mobility under off-road conditions, the U.S. Army Tank Automotive Research, Development, and Engineering Center (TARDEC), now called Ground Vehicle Systems Center, has developed the NATO Reference Mobility Model (NRMM) in the 1960s [68]. To overcome the limitation of NRMM, the next-generation NRMM (NG-NRMM) is being developed to make the prediction of mobility models more accurate and efficient using multibody dynamic modeling and new simulation techniques [6,9,10].

In addition to the NRMM and NG-NRMM, many other mobility analysis models have also been developed in the past decades to predict off-road vehicle performance under different off-road conditions. For instance, Recuero et al. [11] introduced a simulation approach of off-road vehicles using nonlinear finite element tires on granular material. Laughery et al. [12] suggested Bekker's derived terramechanics model using Bekker’s equations to evaluate the off-road mobility. Hetherington [13] proposed a model using mean maximum pressure to specify the off-road performance of vehicles. On the basis of Bekker–Wong theory, Senatore and Sandu [14] designed a 14-degrees-of-freedom vehicle model to study the torque distribution of wheeled vehicles. Serban et al. [15] introduced a co-simulation framework for vehicle–terrain interaction simulation.

Since the high-fidelity vehicle simulation models usually require high computational costs and the semi-empirical models are not accurate enough, surrogate models based on data of high-fidelity simulations or tests are a promising solution to balance the model accuracy and computational efficiency for off-road AGV mobility prediction. They have attracted much attention in the off-road AGV community in recent years. For example, González et al. [16] used the Kriging surrogate modeling method to model terrain elevation in predicting off-road mobility over large spatial regions. A Kriging model is used in Ref. [17] to generate the mobility map for a given area. Choi et al. [18] developed a dynamic Kriging-based method to improve the effectiveness of surrogate modeling in generating stochastic off-road mobility map. Liu et al. [19] developed a single-loop Kriging surrogate modeling method to analyze the mission mobility reliability of the off-road AGVs.

The aforementioned surrogate modeling approaches use speed-made-good, a quasi-static quantity, as the mobility criterion of an AGV. However, a vehicle is a very complicated dynamic system with many other attributes, which can be used to represent different aspects of its mobility, such as energy consumption and vertical acceleration. Mobility quantification using vertical acceleration exhibits complex dynamic behaviors under different off-road conditions. Thus, more sophisticated surrogate modeling methods are needed to construct surrogate models.

Aiming to construct surrogate models for dynamic responses of vehicles, many surrogate modeling approaches have been proposed using machine learning methods in recent years. For example, Hsu et al. [20] proposed an end-to-end method using a deep convolution neural network to construct a mobility model that predicts vehicle steering and acceleration. Ye et al. [21] introduced a deep learning network model called MBSNet to perform multibody simulation analysis for a vehicle-track system. Pan et al. [22] used a deep neural network to predict the vehicle lateral-longitudinal dynamics. Similarly, Paparusso et al. [23] proposed a framework using a deep recurrent neural network based on long short-term memory (LSTM) autoencoders to predict the driver-vehicle system dynamics.

Although great progress has been made in the surrogate modeling of vehicle dynamics, it is observed in our study that it is very difficult to accurately predict the complicated vehicle behavior using a single surrogate model even with the LSTM method. This is attributed to the fact that the vehicle dynamics is highly nonstationary under rough off-road conditions. For instance, the vibration of an AGV could be of high frequency in some regions and low frequency in some other regions. In addition, the complicated vehicle–terrain interaction adds another layer of complexity to dynamic behavior, which makes the surrogate modeling challenging. This article aims to overcome the challenges in surrogate modeling of vehicle mobility with nonstationary dynamic behaviors by developing a dynamic ensemble framework of nonlinear autoregressive models with exogenous inputs (NARX) using a variable number of lags. Multiple NARX models are trained for different dynamic behaviors based on clustering of the training data. The NARX models are then assembled together dynamically over time to predict the dynamic prediction of AGV mobility over a long time period. The accuracy of the proposed framework is investigated using multiple machine learning models including Gaussian process regression and neural networks. A case study demonstrates the advantages of the proposed surrogate modeling method using a dynamic ensemble of NARX models (DENA) over several classical surrogate modeling methods.

The remainder of this article is organized as follows: Sec. 2 introduces the background of vehicle dynamics simulation using pychrono and NARX models for surrogate modeling of dynamic systems. Section 3 describes the proposed dynamic ensemble of NARX models with a variable number of lags. Section 4 uses case studies to demonstrate the proposed approach. Finally, Sec. 5 draws conclusions.

## 2 Background

### 2.1 Off-Road Autonomous Ground Vehicles Dynamic Model

#### 2.1.1 The Need of a Data-Driven Dynamic Model.

The importance of understanding and being able to determine the response of complex dynamic systems has long been recognized [24]. The analysis traditionally begins with mathematical relations and physical laws that work as descriptions and representations of the various behaviors of interacting components within systems. However, in modeling the dynamic behavior of systems, the governing physical laws are not always available, and explicit analytical solutions are not always tractable. This has motivated an increasing amount of research toward data-driven methodologies using of data-driven analysis and modeling [25].

The off-road AGV of Fig. 1 is an example of a complicated dynamic problem, involving complex geometries, nonlinear interactions within components, and often stochastic interactions with the environment. While advancements in software and computation have promoted the use of techniques such as finite element modeling and other numerical simulations, problems in this field are commonly broken into parts and addressed in a limited sense using simplified math models. This is often done by allowing key system characteristics to be represented with simple continuously deformable bodies, concentrated masses, and idealized springs and dampers, such that the often nonlinear dynamics are modeled and analyzed [24]. Many vehicle dynamics simulation codes are constructed and operated in this manner.

Fig. 1
Fig. 1
Close modal

This article utilizes one such software, a python version of chrono: an open-source multiphysics dynamics engine, to simulate vehicular dynamics in an off-road environment [26]. Using a collection of information and data pertaining to the terrain morphology, the vehicle motion is provided as output of the software. We use the latter to build a data-driven approach to meta-model the complex dynamics of a high mobility multipurpose wheeled vehicle (HMMWV) traversing over an off-road terrain. The data-driven model drastically reduces the required computational effort for AGV mobility prediction and thus supports reliability analysis and AGV mission planning under uncertainty [27].

#### 2.1.2 Off-Road Autonomous Ground Vehicles Mobility Simulation Using pychrono.

pychrono (a python version of chrono) is an open-source multiphysics dynamics engine. The software specializes in handling complex and large dynamic systems with many rigid bodies that interact through various contact methods [26,28]. chrono is capable of solving problems involving multiple physics regimes simultaneously. For example, it can compute the dynamics of interconnected rigid bodies (such as moving parts of a mechanism) while computing the dynamics of deformable bodies (such as terrain-mechanics). For this reason, it has been widely employed for vehicle simulation research. Multiple “demo vehicles” and a “demo terrain” are currently included within its installation files. This allows users to easily begin simulating vehicles in a stochastic environment and obtain synthetic mobility data.

A user begins a vehicle simulation by specifying several inputs: a vehicle model type, vehicle parameters (suspension, driveline type, etc.), a terrain height map, a terrain type (rigid or deformable), a contact method for the vehicle–terrain interaction, initial location and orientation of the vehicle on the terrain, and several additional inputs on how the vehicle is powered and controlled throughout the simulation. The latter input can include, for example, a fixed or variable throttle, manual controls, or a speed and/or steering controller to follow a predefined path on the terrain. In addition, parameters to control the numerical computations, such as mesh/grid size of the deformable terrain, and simulation time-step size, are inputs by the user.

A large selection of output data can be extracted using the pychrono vehicle simulation. Essential kinematic information can be recorded for multiple parts of the vehicle such as the global position of the vehicle’s center of mass, position of the driver, vehicle chassis velocity, center of mass (CoM) velocity, vehicle driver’s acceleration and CoM acceleration, wheel velocity, and vehicle rotation. Additional information can be captured involving the vehicle mechanisms such as steering, driveshaft speed, axle torque, wheel force, and powertrain output. For our research, the vehicle’s vertical acceleration is obtained and used as an output.

As indicated in Fig. 2, pychrono uses a terrain map (Fig. 2(a)) and a vehicle model (Fig. 2(b)) to predict various vehicle mobility quantities, such as vertical acceleration (Fig. 2(c)) for a given path.

Fig. 2
Fig. 2
Close modal
Mathematically, a pychrono simulation model can be represented as follows:
$a=g(u)+ε$
(1)
where g(·) stands for a pychrono simulation model, $u=[u1,u2,…,uN]$ is a vector of N discretized elevation values over a path, $ε$ is a vector of simulation noise values due to numerical discretization and model assumptions, and $a=[a1,a2,…,aN]$ is a vector of the predicted vehicle vertical accelerations. Note that the length of $u$ is different for different paths, and therefore, N varies with the path.
In this article, the terrain elevation is modeled as a random field using the Karhunen–Loève (KL) expansion method [29], since it is a widely used method to represent natural variability over space [30] and has been used for stochastic modeling of terrain profiles [31]. For a given coordinate $d$ on the map (Fig. 2(a)), the terrain elevation is represented by
$Yu(d)=μu(d)+σu(d)∑l=1NKLλlξlηl(d)$
(2)
where $μu(d)$ and $σu(d)$ denote the mean and standard deviation of the terrain elevation at $d$, respectively, ξl, l = 1, 2, …, NKL are independent standard normal random variables, NKL is the number of expansion terms, and λl and ηl are the eigenvalues and eigenvectors, respectively, of a correlation matrix whose elements are given by
$ρ(di,dj)=exp[−β1(d1i−d1j)2−β2(d2i−x2j)2],∀i,j=1,2,…,Nd$
(3)
where $di=[d1i,d2i]$ are the spatial coordinates at location $di$, and Nd is the total number of discrete spatial points in the map used for KL expansion.

Even though the pychrono simulation can effectively predict various vehicle mobility quantities, the computation is expensive if a very fine mesh is used. The high computational cost hinders its application in reliability analysis and mission planning under uncertainty. A data-driven model must be constructed to replace the high-fidelity simulation model.

### 2.2 Nonlinear Auto Regressive eXogenous Model.

Many dynamic systems evolve over space and time. To describe and analyze such systems, appropriate methods must be employed. Most dynamic systems are modeled using a state variable representation providing the evolution of state variables through time. This implies utilizing a collection of various time series representing the system input and output.

Time series prediction and forecasting are often done using autoregressive modeling, which accounts for potential for relationship (correlation) within the sequence of values. An autoregressive model uses a basis where the predictands are previous occurrences of the predictor. If the basis is a linear combination of previous outputs, it is called an autoregressive (AR) model and is given by
$yi=c1yi−1+c2yi−2+…+cqyi−q$
(4)
or
$yi=∑j=1qcjyi−j$
(5)
where yi is the output (state variable) of interest, i is a time index up to q steps or “lags,” and cj, j = 1, …, q are the coefficient terms to fit the model. A simple example is the rule to generate the Fibonacci sequence where any term is computed as the sum of the two previous terms.

If the model parameters are adequately fit (typically done by minimizing an error metric using a data training set), the model can be used for prediction. The value at a future time-step can be then computed using the values of previous steps. A computer code can iterate this process allowing the model to march forward in time for forecasting. It should be mentioned that the prediction phase may need to assume prior values (up to q lags) to initiate the simulation and continue with forward propagation.

Such models can also incorporate external terms as input if they affect the state variable of interest. These are called exogenous terms and can be similarly included by appending a similarly structured sequence of exogenous lags as additional predictands as follows:
$yi=∑j=1qcjyi−j+∑k=1pbkui−k$
(6)
in which uik, k = 1, …, p is the exogenous input, and bk, k = 1, …, p are input coefficients to fit the model. Note that if the present time exogenous input (i.e., k = 0) is known, it can also be included as a predictor.
Finally, if a nonlinear basis of terms is used for prediction, or if u and y are used as arguments of a nonlinear function, the method (and naturally the name) becomes nonlinear auto regressive exogenous (NARX) modeling as shown in Eq. (7):
$yi=F(ui−1,…,ui−p;yi−1,…,yi−q)$
(7)
where the function F(·) is a nonlinear mapping from lags of previous outputs and inputs to the current output. This nonlinear function can be a polynomial, a Gaussian process regressor (GP-NARX) [32], or a feed forward artificial neural network (NARX-Net) [33].

Modeling of this kind is akin to finite difference approximations of differentiation, and numerical methods in solving differential equations. For this reason, NARX models are particularly preferred if there are underlying or otherwise undiscovered differential relations governing the behavior of the system, which is often the case for physics-based dynamic problems. Other related techniques have also been developed and have shown success in time series modeling of dynamic systems. Examples are the dynamic mode decomposition, the popular SINDy algorithm, system identification methods, and a host of many different recurrent neural network architectures (LSTM, GRU, etc.) [25,34]. All of these techniques share the common presumption that there exists a tractable relationship between a variable’s future value in time and its previous values in time. Some other methods are mathematically similar [35,36]. Great success in modeling even complex and highly nonlinear dynamical systems has been achieved using techniques such as given in Refs. [25,34,37]. However, extrapolating to forecast significantly far into the future of a system’s evolution, and still retaining good prediction accuracy, is still an active research topic.

In this article, both the GP-NARX and NARX-Net frameworks will be investigated in modeling and forecasting the vertical acceleration of an AGV as it traverses over an off-road terrain. In addition to GP-NARX and NARX-net, the proposed method is also compared with LSTM, which has gained much attention in recent years for modeling various dynamic systems.

### 2.3 Long Short-Term Memory.

LSTM is a variant of a recurrent neural network. It has attracted much interest in dynamic system modeling since it allows a model to sequentially learn and determine whether to forget previous states and/or update the states. A LSTM module first determines what previous information to forget using a “forget gate” [38,39] and then determines what new information $C~t$ to store using an “input gate” and a tanh layer. Based on this, an LSTM cell updates the new cell state Ct by combining the information inherited from the previous cell state Ct−1 with the new candidate information $C~t$. Finally, the output of an LSTM cell is obtained.

The four-step operation of an LSTM cell is given as follows [38,39]:
$Step1:ft=σ(Wf[ht−1,xt]+bf);Step2:it=σ(Wi[ht−1,xt]+bi)C~t=tanh(WC[ht−1,xt]+bC);Step3:Ct=ft×Ct−1+it×C~t;Step4:ot=σ(Wo[ht−1,xt]+bo)ht=ot×tanh(Ct)$
(8)
where Wf, Wi, WC, and Wo denote the weight matrices of each unit, and bf, bi, bC, and bo are the bias terms associated with each unit.

## 3 Dynamic Ensemble of NARX Models With Different Numbers of Lags for Vehicle Mobility Prediction

As shown in Fig. 2(c), the vertical acceleration of an AGV in off-road environment exhibits a nonstationary behavior. This poses challenges for surrogate modeling of vehicle mobility, resulting in a poor performance for GP-NARX, NARX-net, and LSTM as indicated in the case study in Sec. 4.

Different approaches have been developed in the past decade and recently to deal with nonstationary response in surrogate modeling, such as a nonstationary Gaussian process regression model [40], hierarchical multiscale surrogate models [41,42], and neural network-based methods [43]. The current methods, however, are not generalizable and are difficult to implement in practice. For instance, the approach of nonstationary Gaussian process cannot be applied to neural network-based models. Moreover, most of current methods focus on static, quasi-static, or univariate time series modeling. In this section, we propose an alternative approach to the existing methods. The proposed dynamic ensemble of NARX models is easy to implement and generalizable to different types of models including both GP-based and neural network-based methods.

The basic idea of the proposed method is to first partition the prediction domain into different segments. This is achieved through clustering techniques using a Gaussian mixture model (GMM). Using the partition of the prediction domain, different NARX models are constructed and are then ensembled together dynamically over time to predict vehicle mobility. Two strategies are proposed under this dynamic ensemble framework: (1) a dynamic ensemble of NARX models with a fixed number of lags, and (2) a dynamic ensemble of NARX models with a DIFFerent number of lags (DENA-DIFF). In what follows, we explain the proposed framework in detail.

### 3.1 Gaussian Mixture Model for Partition of Dynamic Domains.

If Np paths are used to generate synthetic mobility data for the off-road AGV, we have the terrain elevation along the Np training paths as ui,j, i = 1, 2, …, Np; j = 1, 2, …, Ni, where Ni is the length of the ith training path. The mobility simulation model of Sec. 2.1 is run for each of the Np training paths to obtain the corresponding vertical acceleration of the AGV as ai,j, i = 1, 2, …, Np; j = 1, 2, …, Ni. According to the NARX scheme of Sec. 2.2, the synthetic mobility data of the Np paths can be converted into a training data matrix as follows:
$xT=[u1u2⋮uNTa1a2⋮aNT],aT=[aT1aT2⋮aTNT]$
(9)
where $xT$ is the input training data matrix including previous lags of both terrain elevation and vehicle accelerations, $aT$ is the output training data matrix (i.e., vertical acceleration of AGV), and NT is the total number of training data collected from the Np paths. Since each path is used to generate multiple training data point due to the NARX scheme, NT is much greater than Np. The vectors $u1=[u11,u12,…,u1p]∈R1×p$ and $a1=[a11,a12,…,a1q]∈R1×q$ are, respectively, the first p and q time-steps of the first training path, and the corresponding output is aT1a1(q+1). The remaining rows of the training data in Eq. (9) can be explained similarly.
Since the vehicle has a different dynamic behavior in different segments of a path, we partition the training data of Eq. (9) into different segments using a GMM. A GMM is a unsupervised machine learning method, which approximates the joint probability density function (PDF) of multivariate data as a weighted sum of multivariate Gaussian components as follows [44]:
$f(xT)=∑i=1Qwiϕ(xT|μi,Σi)$
(10)
where $xT$ is the vector of input variables (terrain elevation and vehicle acceleration) with a dimension of p + q, Q is the number of Gaussian components, wi is the weight of the ith Gaussian component, $∑i=1Qwi=1$, and $ϕ(xT|μi,Σi),i=1,2,…,Q$ is the PDF of the ith multivariate Gaussian distribution:
$ϕ(xT|μi,Σi)=1(2π)N/2|Σi|1/2exp{−12(xT−μi)′Σi−1(xT−μi)}(11)$
(11)
where $μi$ is the mean vector and $Σi$ is the covariance matrix.
The application of GMM entails the estimation of (1) the GMM parameters and (2) the optimal number of components/clusters. A commonly used and matured technique for the GMM parameter estimation is the estimation-maximization (EM) algorithm [45]. However, the EM algorithm is sensitive to the starting point and may converge to a local optimum. In this article, we use the variational inference (VI) method to overcome this limitation. In VI, the variational posterior distribution of GMM parameters is expressed as a multiplication of a Gaussian and a Wishart distribution as follows [46]:
$f(μi,Σi)=N(μi|mi,(βiΣi)−1)W(Σi|Oi,vi),i=1,…,Q(12)$
(12)
where N(·) and W(·) are, respectively, the PDF of a Gaussian and a Wishart distribution, $mi$ is the posterior mean vector of $μi$, βi is the confidence in the estimation of $mi$, and vi represents the extent of effect on the estimated $Oi$, which is a scaled parameter of the covariance matrix $Σi$.

The ith posterior mean and covariance matrix for the ith cluster, $μi$ and $Σi$, can be estimated based on the conjugacy properties of the normal inverse Wishart prior distribution [46]. As a result, $mi$ can be interchanged with $μi$. For a given number of Gaussian components/clusters, the GMM parameters can be estimated using the VI method.

To address the second issue related to the optimal number of clusters, the Dirichlet process (DP) is employed in this article to identify the cluster number in an iterative way based on a DP prior distribution. The probability of assigning data samples to the ith cluster can be estimated by a famous Chinese restaurant process as follows [47]:
$p(i|α)={Nci+αNT+α−1,ifi
(13)
where Nci is the number of samples in the ith cluster and α is a parameter that governs the cluster behavior and reflects the growth rate of clusters in observed data samples. It implies that the total number of clusters increases with α.
To select an optimal α, a prior PDF is usually used and then α is calibrated with a Gibbs sampler based on data samples. The PDF of α is derived by Escobar and West [48] as follows:
$p(α|η,i)=RηG(α+i,b−log(η))+(1−Rη)Γ(α+i−1,b−log(η))$
(14)
where $Rη$ is a weighting factor defined as $Rη/(1−Rη)=(α+i−1)/(NT(b−log(η)))$, η is a continuous quantity between 0 and 1 following a beta distribution p(η|α, i) = B(α + 1, η), and Γ( ·, · ) stands for a Gamma distribution. By using Eqs. (13) and (14), the optimal number of clusters Q and optimal α can be estimated through the Gibbs sampling method [49].
If $θ≜{wi,μi,Σi,i=1,2,…,Q}$ includes the GMM parameters and $P(xT,j,i|θ)$ is the probability that the jth data sample belongs to the ith cluster, which is predicted using the learned GMM, we have $∑i=1QP(xT,j,i|θ)=1$. To partition the data matrix of Eq. (9) into different segments while maintaining continuity of the surrogate model prediction, we also impose the following condition during partition
$I(xT,j,i)={1,ifP(xT,j,i|θ)≥0.05Q−1,0,otherwise,$
(15)
where $I(xT,j,i)=1$ means $xT,j$ is classified with the ith cluster. A user-specified value of (0.05/Q − 1) ensures continuity of the response. The case study presented in Sec. 4 shows that this value works well.
Based on Eq. (15), the input training data of the ith cluster is expressed as follows:
$xc,i={x|I(x,i)=1,∀x∈xT},i=1,…,Q$
(16)
According to Eqs. (15) and (16), it is quite possible that a sample may belong to multiple clusters, indicating that such samples will be used as training data for the training of multiple NARX models. If $ac,i$ is the corresponding output training data of the ith cluster, we have
$xT={xc,1∪⋯∪xc,Q},aT={ac,1∪⋯∪ac,Q}$
(17)

Moreover, if the number of samples in a certain cluster is too low, the NARX model prediction will be of very low quality. If the low-quality NARX model is then ensembled with the other models, it may ruin the overall prediction. To avoid this situation, another condition during the partition of the data is set as follows:

$L(xc,i)>NTQ$
(18)
where $L(xc,i)$ represents the number of samples in the ith cluster.
If the aforementioned condition is violated during the data partition, samples from the nearby clusters will be added into the corresponding cluster to ensure the condition of Eq. (18) is satisfied. After the partition of the data matrix of Eq. (9), the data are grouped into different clusters to represent different dynamic behaviors of the AGV in the off-road environment. Finally, we train Q NARX models as follows:
$a=g^a,i(up,aq),i=1,2,…,Q$
(19)
where $up=[u1,u2,…,up]$ and $aq=[a1,a2,…,aq]$ represent the p and q lags of the terrain elevation and vehicle acceleration, and $g^a,i(⋅)$ is the ith NARX model trained using data $xc,1$ and $ac,1$.

The NARX models in Eq. (19) can be of GP-NARX or NARX-Net type. Other types of models can also be used, such as mixed GP-NARX, NARX-net, and LSTM models. Figure 3 summarizes the overall procedure for the partition of dynamic domains using GMM and subsequent training of NARX models. When a NARX-Net is employed, a neural network is trained to capture the nonlinear dynamic behavior of the vehicle over time. If a GP-NARX is employed, a Gaussian process regression model is trained in the NARX scheme to capture the nonlinear function relationship (see Eq. (7) in Sec. 2.2).

Fig. 3
Fig. 3
Close modal

### 3.2 Dynamic Ensemble of NARX Models.

After the Q NARX models are trained, they must be dynamically assembled together over time to predict the vehicle vertical acceleration over any given new path. If the terrain elevation of a new path be $unew=[u1,u2,…,uNnew]$, where Nnew is the number of steps on the new path, the task is to use the NARX models to predict the vehicle mobility without running the high-fidelity mobility simulation.

Taking the prediction at step tk as an example, the inputs $xT,k=[uk,ak]$ to the NARX models are expressed as follows:
$uk=[uk−p−1,uk−p+1,…,uk],ak=[ak−q,ak−q+1,…,ak−1]$
(20)
where $uk$ are the previous p lags of terrain elevation of tk and $ak$ represents the previous q lags of vehicle vertical acceleration.
Using the inputs of Eq. (20), the probability that the inputs/dynamic behavior $xT,k$ at tk belongs to the ith cluster is predicted using the trained GMM presented in Sec. 3.1 as follows:
$pki=P(xT,k,i|θ),i=1,2,…,Q$
(21)
where $∑i=1Qpki=1$.
Using the predicted probabilities, the NARX predictions can be then ensembled together to obtain the predicted vehicle acceleration at tk as follows:
$a^k=∑i=1Qw(xT,k,i)g^a,i(uk,ak)$
(22)
where $w(xT,k,i)=p⌢ki/∑i=1Qp⌢ki$ with $p⌢ki$ given by
$p⌢ki={pki,ifpki≥0.05Q−1.0,otherwise$
(23)
Equation (23) indicates that if the probability that the input/dynamic behavior belongs to a certain cluster is very low, the sample may be at the boundary of that cluster. For such cases, the corresponding NARX model of that cluster usually has a relatively poor performance. Removing the corresponding NARX model from the ensemble avoids contamination of the final ensembled prediction by the NARX model. As shown in Eq. (22), the ensemble weights $w(xT,k,i),i=1,…,Q$ depend on the input $xT,k$. This implies that the ensemble weights will adjust dynamically over time according to the behavior of the response, which enables us to dynamically ensemble the NARX models over time. In addition, as mentioned earlier, the NARX models can be of any type (e.g., GP-NARX or NARX-Net). If GP-NARX models are used, the prediction mean and standard deviation at step tk are expressed as follows:
$μa,i=∑i=1Qw(xT,k,i)μa,i(uk,ak)$
(24)
$σa,i=∑i=1Qw2(xT,k,i)σa,i2(uk,ak)$
(25)
where $μa,i(uk,ak)$ and $σa,i2(uk,ak)$ are the mean and variance of the prediction from the ith GP-NARX model at tk.
After we obtain the prediction of vehicle acceleration at tk using the ensemble of the NARX models, the prediction is used as input to predict the vehicle mobility at the next step tk+1 as follows:
$uk+1=[uk+1−p,uk+2−p,…,uk],a~k+1=[ak+1−q,ak+2−q,…,a^k⏟predicted]$
(26)

The updated inputs are then used to determine the ensemble weights for tk+1 using Eqs. (22) and (23). The aforementioned process is repeated recursively over time allowing us to predict the vehicle mobility over a long time period. This is crucial for mobility prediction of off-road AGVs since the length of paths may be different. Figure 4 shows the flowchart of the proposed dynamic ensemble of NARX models based on the partition of dynamic domains using GMM.

Fig. 4
Fig. 4
Close modal

### 3.3 Dynamic Ensemble of NARX Models With Different Numbers of Lags.

In the DENA method presented in Secs. 3.1 and 3.2, the number of lags is the same across different NARX models. Considering that the response over time may vary for different dynamic behaviors, a constant number of lags for each NARX model may not be optimal. In the ideal scenario, the number of lags in the NARX model should be adaptable to different regions of the nonlinear dynamic response that the NARX model tries to predict. For instance, the number of lags can be large for high-frequency response exhibiting many peaks in a short time and small for low-frequency response. Motived by this observation, this subsection develops the DENA-DIFF method based on the dynamic ensemble method of Sec. 3.2. The basic idea is to use different number of lags for different clusters, where each cluster represents a different dynamic behavior and has its own NARX model.

To explain the proposed DENA-DIFF method, we represent the training data of the ith cluster from Sec. 3.2 as follows:

$xc,i=[ui1ui2⋮uiNTai1ai2⋮aiNT],i=1,…,Q$
(27)
where $ui1=[uk−p,uk−p+1,…,uk−1]∈R1×p$ and $ai1=[ak−q,$$ak−q+1,…,ak−1]∈R1×q$. We then generate V different training datasets by adding more lags to the data matrix of Eq. (27) as follows:
$xc,i(r)=[u~i1(r)u~i2(r)⋮u~iNT(r)a~i1(r)a~i2(r)⋮a~iNT(r)],i=1,…,Q;r=1,…,V$
(28)
where $xc,i(r)$ is the rth training dataset of the ith cluster. The $u~i1(r)$ and $a~i1(r)$ vectors of its first row are, respectively, the updated lags for terrain elevations and vertical accelerations by adding r additional lags to the original training dataset of $xc,i$. They are given by
$u~i1(r)=[uk−p−r,…,uk−p−1⏟Addedrelements,ui1]∈R1×(p+r),a~i1(r)=[ak−q−r,…,ak−q−1⏟Addedrelements,ai1]∈R1×(q+r)$
(29)

Similar expressions hold for the other rows of $xc,i(r)$.

For each of the V different training datasets generated from $xc,i$, we divide the input training data and the corresponding output data into a training set and a validation set as follows:

$xc,i(r)={xT,i(r)∪xVa,i(r)},ac,i={aT,i∪aVa,i},∀r=1,…,V$
(30)
where $xT,i(r)$ and $xVa,i(r)$ are, respectively, the training and validation sets of the rth new dataset generated from the ith cluster.
Subsequently, we train V NARX models for the ith cluster using $xT,i(r),r=1,2,…,V$ and $aT,i$ as follows:
$a=g^i,r(up+r,aq+r),∀i=1,2,…,Q;r=1,…,V$
(31)
where $g^i,r(up+r,aq+r),r=1,2,…,V$ represents the NARX model of the ith cluster trained using $xT,i(r)$ and $aT,i$.
To determine the optimal number of lags for the ith cluster, we compute the following prediction error of the rth NARX model for the validation dataset:
$εV,i(r)=‖aV,i−a^Va,i(r)‖2,r=1,…,V$
(32)
where $‖⋅‖2$ is the l2 norm and $a^V,i(r)$ is the prediction obtained from the rth NARX model as follows:
$a^V,i(r)=g^a,i(xVa,i(r),r)$
(33)
where $xVa,i(r)=[u~Va,i(r),a~Va,i(r)]$, with $u~Va,i(r)$ and $a~Va,i(r)$ being the $u~i1(r)$ and $a~i1(r)$ in Eq. (29) for the validation set of the ith cluster after adding r more lags to the original cluster.
The optimal number of lags for the i-cluster is then determined by
$ri*=argminr=1,…,V{εV,i(r)}$
(34)
in which $εV,i(r)$ is the error given in Eq. (32).
The aforementioned procedures are implemented for all clusters. Equations (28)(34) allows us to use different number of lags for different clusters, each representing a different dynamic behavior. After the optimal number of lags is identified for all Q clusters, the NARX models, each with a different number of lags, are ensembled for prediction at tk as follows:
$a^k=∑i=1Qw(xT,k,i)g^i,ri*(up+ri*,aq+ri*)$
(35)
where $w(xT,k,i)$ are the weights determined similarly as in the DENA approach of Sec. 3.2. The remaining recursive implementation procedure of DENA-DIFF is the same with that of the DENA method. Figure 5 summarizes the overall procedure of the DENA-DIFF method.
Fig. 5
Fig. 5
Close modal

## 4 Case Study

In this section, we use an off-road AGV example to demonstrate the proposed DENA and DENA-DIFF methods. Both methods are implemented using GP-NARX and NARX-Net. The DENA method is referred as “DENA-GP” and “DENA-Net,” respectively, to indicate its implementation with GP-NARX and NARX-Net. Similarly, the implementations of DENA-DIFF with GP-NARX and NARX-net are referred as “DENA-DIFF-Net” and “DENA-DIFF-GP,” respectively. The DENA-GP, DENA-Net, DENA-DIFF-GP, and DENA-DIFF-Net methods are compared with the classical GP-NARX and NARX-Net methods without ensemble of NARX models and the LSTM model with extensive tuning. The LSTM model’s performance is optimized by tuning the number of layers, the number of neurons, the dropout rate, and the learning rate. First, we present in Sec. 4.1 how the synthetic AGV mobility data are collected for the training and testing of different methods. Then, Sec. 4.2 compares the performance of different methods.

### 4.1 Generation of Synthetic Mobility Data.

The off-road AGV mobility simulation model discussed in Sec. 2.1.2 is employed to generate synthetic mobility data. The vehicle module is the pychrono’s build-in HMMWV as shown in Fig. 6, representing an off-road military Humvee.

Fig. 6
Fig. 6
Close modal

The key vehicle parameters are given as follows:

• Vehicle mass = 2422.74 kg

• Wheelbase = 3.378 m

• Wheel track = 1.819 m

• Driveline Type: all-wheel drive

During the simulation, the vehicle uses a proportional–integral–derivative controller with preset gains to keep it on the specified path and a fixed throttle value of 100%. The terrain is described as a Gaussian random field with a mean of 0 m, a standard deviation of 0.26 m, and a correlation length of 0.5 m. These parameters generate a random field that closely resembles realistic off-road terrains in terms of variation of hills and bumps in relation to the terrain size. Figure 7 presents ten different random terrain maps generated using the KL expansion method (Eqs. (2) and (3)). A zoomed-in view of one of the terrain maps is shown in Fig. 8, indicating a maximum variability of the vertical height of about 2 m.

Fig. 7
Fig. 7
Close modal
Fig. 8
Fig. 8
Close modal

The terrain is assumed rigid, and the smooth (penalty) contact method is used in the vehicle mobility simulation with terrain friction = 0.9, terrain restitution = 0.01, and terrain Young’s modulus = 2 × 107 Pa. To investigate the impact of the number of training paths on the performance of different methods, two sets of synthetic vehicle mobility data are generated.

• In the first set of data, five paths are generated from each terrain map resulting in 50 paths in total (referred as Scenario 1).

• For the second set of data, nine paths are generated from each map generating therefore 90 paths in total (referred as Scenario 2).

For mobility data collection, the vehicle is guided by its controller during the simulation as it follows a specified path. The key metrics of interest recorded during the simulation at each time-step are the terrain elevation and the vehicle’s vertical acceleration at the center of mass. Figure 9 shows the selected five paths to collect synthetic mobility data for scenario 1. The starting and end points (Si, Ti, i = 1, …, 5) of the five paths are as follows:
${Si=(−140,d2)Ti=(140,d2),d2=±140,0$
(36)
Fig. 9
Fig. 9
Close modal

The nine paths of each terrain map are generated similarly for scenario 2. Figure 10 shows an example (path 1) of the generated synthetic data, indicating that the vertical acceleration fluctuates considerably over the path.

Fig. 10
Fig. 10
Close modal

Section 4.2 compares the accuracy of different methods including GP-NARX, NARX-Net, DENA-GP, DENA-Net, DENA-DIFF-GP, DENA-DIFF-Net, and LSTM. For scenario 1, 40 of the 50 paths are randomly selected for training and the remaining 10 paths are used for testing to compare the performance of different methods. For scenario 2, 80 of 90 are randomly selected for training and 10 are used for testing.

### 4.2 Comparison of Mobility Prediction Using Different Surrogate Models.

For the training of GP-NARX, NARX-Net, DENA-GP, DENA-Net, and LSTM, the number of lags is optimized for each model to obtain the best performance. For DENA-DIFF-GP and DENA-DIFF-Net, the number of lags is optimized for each cluster. As discussed in Sec. 3.1, the Dirichlet process is employed to determine the optimal number of clusters in DENA and DENA-DIFF. Figure 11 shows the clustering iterative process for scenario 1 with 40 training paths and the resulting optimal estimate of the number of clusters, Q = 3. Following that, Fig. 12 presents the clusters of training data for terrain elevation and vehicle vertical acceleration for scenario 1 with 40 training paths and 10 test paths. As shown in this figure, even though there is an overlapping between different clusters, the shapes of the clusters are different, indicating the difference in local dynamic behaviors of the vehicle. Figure 12 is only a 2D plot of 2 of the 40 inputs comprising 20 lags of the terrain elevation and 20 lags of the vehicle acceleration. To better illustrate the difference of clusters, Fig. 13 shows the clustering of different regions for a test path in the dynamic ensemble. It indicates that the prediction regions are classified into different regions representing different local dynamic behaviors.

Fig. 11
Fig. 11
Close modal
Fig. 12
Fig. 12
Close modal
Fig. 13
Fig. 13
Close modal

Table 1 lists the optimal number of lags for DENA-DIFF-GP and DENA-DIFF-Net for scenario 1. The optimal number of lags is different for different clusters and different NARX models.

Table 1

Optimal number of lags for DENA-DIFF (scenario 1)

MethodNumber of lags
Cluster 1Cluster 2Cluster 3
DENA-DIFF-Net101513
DENA-DIFF-GP10915
MethodNumber of lags
Cluster 1Cluster 2Cluster 3
DENA-DIFF-Net101513
DENA-DIFF-GP10915

Figures 1416 compare the prediction of different methods for a test path against the true vertical acceleration in scenario 1. Figures 1719 show the prediction error of different methods for this particular path.

Fig. 14
Fig. 14
Close modal
Fig. 15
Fig. 15
Close modal
Fig. 16
Fig. 16
Close modal
Fig. 17
Fig. 17
Close modal
Fig. 18
Fig. 18
Close modal
Fig. 19
Fig. 19
Close modal

Figure 20 compares the prediction error distribution for different methods for this path. It shows that the DENA-DIFF-GP method performs the best. Subsequently, Table 2 summarizes the root mean square error (RMSE) and the average absolute error of different methods for ten test paths in the studied two scenarios (i.e., one with 40 training paths and the other one with 80 training paths). The results show that the prediction accuracy of NARX model can be improved using the proposed dynamic ensemble of different NARX models. For both scenarios (i.e., 40 and 80 training paths), the accuracy of DENA-GP-DIFF is the best with DENA-GP and NARX-Net. In addition, as the number of training paths increases from 40 to 80, the performance of all methods improves. The accuracy of different methods also gets closer to each other as the number of training paths becomes high. This implies that the ensemble method may not be necessary for vehicle mobility prediction if there is a large number of training paths. A direct NARX-Net model may satisfy the decision maker’s accuracy requirement. The proposed DENA methods are more suitable for applications with a smaller number of available training paths.

Fig. 20
Fig. 20
Close modal
Table 2

Comparison of RMSE and average absolute error of different approaches for ten testing paths using 40 and 80 training paths

MetricsRMSEAverage absolute error
Number of training paths40804080
LSTM4.1943.5173.0272.550
GP-NARX3.8183.8222.6062.613
NARX-Net2.9382.0362.1161.402
DENA-GP2.3011.9211.5121.314
DENA-Net3.7232.7492.7172.031
DENA-DIFF-GP2.0031.9101.3061.236
DENA-DIFF-Net2.9422.5632.7171.926
MetricsRMSEAverage absolute error
Number of training paths40804080
LSTM4.1943.5173.0272.550
GP-NARX3.8183.8222.6062.613
NARX-Net2.9382.0362.1161.402
DENA-GP2.3011.9211.5121.314
DENA-Net3.7232.7492.7172.031
DENA-DIFF-GP2.0031.9101.3061.236
DENA-DIFF-Net2.9422.5632.7171.926

## 5 Conclusion

Building a cheap surrogate model to substitute the original computationally expensive off-road AGV simulation model plays an important role in supporting real-time decision-making and mission planning under uncertainty. It is observed that a single surrogate model cannot accurately capture the complicated nonlinear dynamics of an AGV over time. This article overcomes this challenge by developing a dynamic ensemble framework for different NARX models. The prediction domain is first partitioned into different dynamic regions using a Gaussian mixture model-based clustering method. Based on the partition, different NARX models are constructed and then dynamically ensembled together over time. To account for the dynamic behavior being potentially different at different time-steps, the dynamic ensemble framework is also extended to include NARX models with different number of lags. The results of a case study showed that the proposed DENA framework can improve the prediction accuracy of NARX models for mobility prediction of off-road AGVs in an off-road environment. Even though the developed approach is for off-road AGVs, it is also applicable to other dynamic systems, which have complicated nonlinear dynamics.

This article mainly explores the dynamic ensemble of a single type of the NARX model, namely, GP-NARX or NARX-Net. It can also be extended to dynamic ensemble of different types of models in the NARX format such as the dynamic ensemble of NARX with LSTM. In this article, the terrain profile is modeled as a Gaussian random field. The application of the proposed methods to mobility prediction of off-road AGVs under more complicated terrain represented by a non-Gaussian random field is worth investigating in the future. In addition, this article mainly focuses on simulation-based vehicle mobility prediction. To use the developed framework for reliability-based off-road ground vehicle design optimization, design-related parameters, such as vehicle parameters, must be included into the proposed surrogate modeling framework. Incorporating these design-related parameters in the proposed framework is an interesting research topic that is worth studying.

## Acknowledgment

This work was supported in part by the Automotive Research Center (ARC) in accordance with Cooperative Agreement W56HZV-04-2-0001 U.S. Army CCDC Ground Vehicle Systems Center (GVSC), Warren, MI. The support is gratefully acknowledged. DISTRIBUTION A. Approved for public release; distribution unlimited. OPSEC # 6140.

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

No data, models, or code were generated or used for this paper.

## References

1.
Li
,
M.
,
Song
,
X.
,
Cao
,
H.
, and
Huang
,
Z.
,
2019
, “
Shared Steering Control Combined With Driving Intention for Vehicle Obstacle Avoidance
,”
Proc. Inst. Mech. Eng. Part D J. Automob. Eng.
,
233
(
11
), pp.
2791
2808
.
2.
Wang
,
B.
,
Gong
,
J.
, and
Chen
,
H.
,
2019
, “
Motion Primitives Representation, Extraction and Connection for Automated Vehicle Motion Planning Applications
,”
IEEE Trans. Intell. Transp. Syst.
,
21
(
9
), pp.
3931
3945
.
3.
Liu
,
Q.
,
Zhao
,
L.
,
Tan
,
Z.
, and
Chen
,
W.
,
2017
, “
Global Path Planning for Autonomous Vehicles in Off-Road Environment via an A-Star Algorithm
,”
Int. J. Veh. Auton. Syst.
,
13
(
4
), pp.
330
339
.
4.
Oksanen
,
T.
, and
Visala
,
A.
,
2009
, “
Coverage Path Planning Algorithms for Agricultural Field Machines
,”
J. Field Rob.
,
26
(
8
), pp.
651
668
.
5.
Korlath
,
G.
,
2007
, “
Mobility Analysis of Off-Road Vehicles: Benefits for Development, Procurement and Operation
,”
J. Terramech.
,
44
(
5
), pp.
383
393
.
6.
McCullough
,
M.
,
Jayakumar
,
P.
,
Dasch
,
J.
, and
Gorsich
,
D.
,
2017
, “
The Next Generation NATO Reference Mobility Model Development
,”
J. Terramech.
,
73
(
1
), pp.
49
60
.
7.
Singh
,
A.
,
2012
, “
Mobility of Military Vehicles at TARDEC
,”
Army Tank Automotive Research Development and Engineering Center
, Warren MI.
8.
Ciobotaru
,
T.
,
2009
, “
Semi-Empiric Algorithm for Assessment of the Vehicle Mobility
,”
Leonardo Electron. J. Pract. Technol.
,
8
(
15
), pp.
19
30
.
9.
,
M.
,
Dasch
,
J.
,
Gonzalez
,
R.
,
Hodges
,
H.
,
Jain
,
A.
,
Iagnemma
,
K.
,
Letherwood
,
M.
,
Mccullough
,
M.
,
Priddy
,
J.
, and
And Wojtysiak
,
B.
,
2016
, “
Next-Generation NATO Reference Mobility Model (NG-NRMM)
,”
Tank Automotive Research, Development and Engineering Center (TARDEC)
, Warren, United States.
10.
McCullough
,
M.
,
Jayakumar
,
P.
,
Dasch
,
J.
, and
Gorsich
,
D.
,
2016
, “
Developing the Next Generation NATO Reference Mobility Model
,”
US Army TARDEC, Warren
, MI.
11.
Recuero
,
A.
,
Serban
,
R.
,
Peterson
,
B.
,
Sugiyama
,
H.
,
Jayakumar
,
P.
, and
Negrut
,
D.
,
2017
, “
A High-Fidelity Approach for Vehicle Mobility Simulation: Nonlinear Finite Element Tires Operating on Granular Material
,”
J. Terramech.
,
72
(
1
), pp.
39
54
.
12.
Laughery
,
S.
,
Gerhart
,
G.
, and
Goetz
,
R.
,
1990
, “
Bekker’s Terramechanics Model for Off-Road Vehicle Research
,” TACOM Research Development and Engineering Center, Warren, MI.
13.
Hetherington
,
J.
,
2001
, “
The Applicability of the MMP Concept in Specifying Off-Road Mobility for Wheeled and Tracked Vehicles
,”
J. Terramech.
,
38
(
2
), pp.
63
70
.
14.
Senatore
,
C.
, and
Sandu
,
C.
,
2011
, “
Torque Distribution Influence on Tractive Efficiency and Mobility of Off-Road Wheeled Vehicles
,”
J. Terramech.
,
48
(
5
), pp.
372
383
.
15.
Serban
,
R.
,
Olsen
,
N.
,
Negrut
,
D.
,
Recuero
,
A.
, and
Jayakumar
,
P.
,
2017
, “
A Co-Simulation Framework for High-Performance, High-Fidelity Simulation of Ground Vehicle-Terrain Interaction
,”
Proceedings of the Conference NATO AVT-265 Specialists’ Meeting
,
Vilnius, Lithuania
,
May 12–19
, Vol. 23, p.
24
.
16.
González
,
R.
,
Jayakumar
,
P.
, and
Iagnemma
,
K.
,
2017
, “
Stochastic Mobility Prediction of Ground Vehicles Over Large Spatial Regions: A Geostatistical Approach
,”
Auton. Robots
,
41
(
2
), pp.
311
331
.
17.
Gonzalez
,
R.
,
Jayakumar
,
P.
, and
Iagnemma
,
K.
,
2017
, “
Generation of Stochastic Mobility Maps for Large-Scale Route Planning of Ground Vehicles: A Case Study
,”
J. Terramech.
,
69
(
1
), pp.
1
11
.
18.
Choi
,
K.
,
Jayakumar
,
P.
,
Funk
,
M.
,
Gaul
,
N.
, and
Wasfy
,
T. M.
,
2019
, “
Framework of Reliability-Based Stochastic Mobility Map for Next Generation Nato Reference Mobility Model
,”
ASME J. Comput. Nonlinear Dyn.
,
14
(
2
), p.
021012
.
19.
Liu
,
Y.
,
Jiang
,
C.
,
Mourelatos
,
Z. P.
,
Gorsich
,
D.
,
Jayakumar
,
P.
,
Fu
,
Y.
,
Majcher
,
M.
, and
Hu
,
Z.
,
2021
, “
Simulation-Based Mission Mobility Reliability Analysis of off-Road Ground Vehicles
,”
ASME J. Mech. Des.
,
143
(
3
), p.
031701
.
20.
Hsu
,
T.-M.
,
Wang
,
C.-H.
, and
Chen
,
Y.-R.
,
2018
, “
End-to-End Deep Learning for Autonomous Longitudinal and Lateral Control Based on Vehicle Dynamics
,”
Proceedings of the 2018 International Conference on Artificial Intelligence and Virtual Reality
,
Nov.,
pp.
111
114
.
21.
Ye
,
Y.
,
Huang
,
P.
,
Sun
,
Y.
, and
Shi
,
D.
,
2021
, “
MBSNet: A Deep Learning Model for Multibody Dynamics Simulation and Its Application to a Vehicle-Track System
,”
Mech. Syst. Signal Process
,
157
(
1
), p.
107716
.
22.
Pan
,
Y.
,
Nie
,
X.
,
Dai
,
W.
,
Xu
,
F.
, and
Li
,
Z.
,
2021
, “
A Deep Neural Network-Based Nonlinear Dynamics Model for the Prediction of Lateral-Longitudinal Vehicle Dynamics
,”
Nonlienar Dyn.
,
1
(
1
). DOI: 10.21203/rs.3.rs-629814/v1
23.
Paparusso
,
L.
,
Melzi
,
S.
, and
Braghin
,
F.
,
2021
, “
Real-Time Forecasting of Driver-Vehicle Dynamics on 3D Roads: A Deep-Learning Framework Leveraging Bayesian Optimization
,”
arXiv preprint
. https://arxiv.org/abs/2103.03825v3
24.
Close
,
C. M.
,
Frederick
,
D. K.
, and
Newell
,
J. C.
,
2001
,
Modeling and Analysis of Dynamic Systems
,
John Wiley & Sons
,
New York
.
25.
Brunton
,
S. L.
, and
Kutz
,
J. N.
,
2019
,
Data-Driven Science and Engineering: Machine Learning, Dynamical Systems, and Control
,
Cambridge University Press
,
New York
.
26.
Tasora
,
A.
,
Serban
,
R.
,
Mazhar
,
H.
,
Pazouki
,
A.
,
Melanz
,
D.
,
Fleischmann
,
J.
,
Taylor
,
M.
,
Sugiyama
,
H.
, and
Negrut
,
D.
,
2015
, “
Chrono: An Open Source Multi-Physics Dynamics Engine
,”
Proceedings of the International Conference on High Performance Computing in Science and Engineering
,
Czech Republic
,
May 25–28
, Springer, Cham, pp.
19
49
.
27.
Liu
,
Y.
,
Jiang
,
C.
,
Zhang
,
X.
,
Mourelatos
,
Z. P.
,
Barthlow
,
D.
,
Gorsich
,
D.
,
Singh
,
A.
, and
Hu
,
Z.
,
2022
, “
Reliability-Based Multi-Vehicle Path Planning Under Uncertainty Using a Bio-Inspired Approach
,”
ASME J. Mech. Des.
,
144
(
9
), p.
091701
.
28.
Mazhar
,
H.
,
Heyn
,
T.
,
Pazouki
,
A.
,
Melanz
,
D.
,
Seidl
,
A.
,
Bartholomew
,
A.
,
Tasora
,
A.
, and
Negrut
,
D.
,
2013
, “
Chrono: A Parallel Multi-Physics Library for Rigid-Body, Flexible-Body, and Fluid Dynamics
,”
Mech. Sci.
,
4
(
1
), pp.
49
64
.
29.
Huang
,
S.
,
,
S.
, and
Rebba
,
R.
,
2007
, “
Collocation-Based Stochastic Finite Element Analysis for Random Field Problems
,”
Probabilistic Eng. Mech.
,
22
(
2
), pp.
194
205
.
30.
Zhang
,
J.
, and
Ellingwood
,
B.
,
1994
, “
Orthogonal Series Expansions of Random Fields in Reliability Analysis
,”
J. Eng. Mech.
,
120
(
12
), pp.
2660
2677
.
31.
Sandu
,
C.
,
Sandu
,
A.
, and
Li
,
L.
,
2005
, “
Stochastic Modeling of Terrain Profiles and Soil Parameters
,”
SAE Int. J. Commer. Veh.
, 1, pp.
211
220
.
32.
Worden
,
K.
,
Becker
,
W. E.
,
Rogers
,
T. J.
, and
Cross
,
E. J.
,
2018
, “
On the Confidence Bounds of Gaussian Process NARX Models and Their Higher-Order Frequency Response Functions
,”
Mech. Syst. Signal Process
,
104
(
1
), pp.
188
223
.
33.
Khalid
,
A.
,
Sundararajan
,
A.
, and
Sarwat
,
A. I.
,
2019
, “
An Arima-Narx Model to Predict Li-Ion State of Charge for Unknown Charge/Discharge Rates
,”
Proceedings of the 2019 IEEE Transportation Electrification Conference (ITEC-India)
,
Bengaluru, India
,
Dec. 17–19
, IEEE, pp.
1
4
.
34.
Li
,
M.
, and
Wang
,
Z.
,
2022
, “
LSTM-Augmented Deep Networks for Time-Variant Reliability Assessment of Dynamic Systems
,”
Reliab. Eng. Syst. Saf.
,
217
(
1
), p.
108014
.
35.
Siegelmann
,
H. T.
,
Horne
,
B. G.
, and
Giles
,
C. L.
,
1997
, “
Computational Capabilities of Recurrent NARX Neural Networks
,”
IEEE Trans. Syst. Man Cybern. Part B Cybern.
,
27
(
2
), pp.
208
215
.
36.
Sum
,
J.
,
Kan
,
W.-K.
, and
Young
,
G. H.
,
1999
, “
A Note on the Equivalence of NARX and RNN
,”
Neural Comput. Appl.
,
8
(
1
), pp.
33
39
.
37.
Menezes
,
J. M. P.
, Jr
, and
Barreto
,
G. A.
,
2008
, “
Long-Term Time Series Prediction With the NARX Network: An Empirical Evaluation
,”
Neurocomputing
,
71
(
16–18
), pp.
3335
3343
.
38.
Olah
,
C.
,
2015
, “
Understanding LSTM Networks
,” http://colah.github.io/posts/2015-08-Understanding-LSTMs/, Accessed February 2022.
39.
Zhang
,
X.
, and
,
S.
,
2020
, “
Bayesian Neural Networks for Flight Trajectory Prediction and Safety Assessment
,”
Decis. Support Syst.
,
131
(
1
), p.
113246
.
40.
Xiong
,
Y.
,
Chen
,
W.
,
Apley
,
D.
, and
Ding
,
X.
,
2007
, “
A Non-Stationary Covariance-Based Kriging Method for Metamodelling in Engineering Design
,”
Int. J. Numer. Methods Eng.
,
71
(
6
), pp.
733
756
.
41.
Chung
,
J.
,
Ahn
,
S.
, and
Bengio
,
Y.
,
2016
, “
Hierarchical Multiscale Recurrent Neural Networks
,”
arXiv preprint
. https://arxiv.org/abs/1609.01704v7
42.
Hu
,
Z.
, and
Mourelatos
,
Z.
,
2018
, “
Efficient Global Surrogate Modeling Based on Multi-Layer Sampling
,”
SAE Int. J. Mater. Manuf.
,
11
(
4
), pp.
385
400
.
43.
Ghazali
,
R.
,
Hussain
,
A. J.
,
Nawi
,
N. M.
, and
,
B.
,
2009
, “
Non-Stationary and Stationary Prediction of Financial Time Series Using Dynamic Ridge Polynomial Neural Network
,”
Neurocomputing
,
72
(
10–12
), pp.
2359
2367
.
44.
Reynolds
,
D. A.
,
2009
, “
Gaussian Mixture Models
,”
Encycl. Biometrics
,
741
(
1
), pp.
659
663
.
45.
,
D. I.
,
Mourelatos
,
Z. P.
, and
Hu
,
Z.
,
2019
, “
Reliability Analysis Using Second-Order Saddlepoint Approximation and Mixture Distributions
,”
ASME J. Mech. Des.
,
141
(
2
), p.
021401
.
46.
Murty
,
M. N.
, and
Devi
,
V. S.
,
2015
,
Introduction to Pattern Recognition and Machine Learning
,
World Scientific
,
Singapore
.
47.
Blei
,
D. M.
,
Griffiths
,
T. L.
,
Jordan
,
M. I.
, and
Tenenbaum
,
J. B.
,
2003
, “
Hierarchical Topic Models and the Nested Chinese Restaurant Process
,”
Proceedings of the NIPS
,
,
Dec. 8–13
.
48.
Escobar
,
M. D.
, and
West
,
M.
,
1995
, “
Bayesian Density Estimation and Inference Using Mixtures
,”
J. Am. Stat. Assoc.
,
90
(
430
), pp.
577
588
.
49.
Murphy
,
K. P.
,
2012
,
Machine Learning: A Probabilistic Perspective
,
MIT Press
,
Cambridge, MA
.