## 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 [6–8]. 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.

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.

*g*(·) stands for a pychrono simulation model, $u=[u1,u2,\u2026,uN]$ is a vector of

*N*discretized elevation values over a path, $\epsilon $ is a vector of simulation noise values due to numerical discretization and model assumptions, and $a=[a1,a2,\u2026,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.

*ξ*

_{l},

*l*= 1, 2, …,

*N*

_{KL}are independent standard normal random variables,

*N*

_{KL}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

*N*

_{d}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.

*y*

_{i}is the output (state variable) of interest,

*i*is a time index up to

*q*steps or “lags,” and

*c*

_{j},

*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.

*exogenous*terms and can be similarly included by appending a similarly structured sequence of exogenous lags as additional predictands as follows:

*u*

_{i−k},

*k*= 1, …,

*p*is the exogenous input, and

*b*

_{k},

*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.

*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):

*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 *C*_{t} by combining the information inherited from the previous cell state *C*_{t−1} with the new candidate information $C~t$. Finally, the output of an LSTM cell is obtained.

*W*

_{f},

*W*

_{i},

*W*

_{C}, and

*W*

_{o}denote the weight matrices of each unit, and

*b*

_{f},

*b*

_{i},

*b*

_{C}, and

*b*

_{o}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.

*N*

_{p}paths are used to generate synthetic mobility data for the off-road AGV, we have the terrain elevation along the

*N*

_{p}training paths as

*u*

_{i,j},

*i*= 1, 2, …,

*N*

_{p};

*j*= 1, 2, …,

*N*

_{i}, where

*N*

_{i}is the length of the

*i*th training path. The mobility simulation model of Sec. 2.1 is run for each of the

*N*

_{p}training paths to obtain the corresponding vertical acceleration of the AGV as

*a*

_{i,j},

*i*= 1, 2, …,

*N*

_{p};

*j*= 1, 2, …,

*N*

_{i}. According to the NARX scheme of Sec. 2.2, the synthetic mobility data of the

*N*

_{p}paths can be converted into a training data matrix as follows:

*N*

_{T}is the total number of training data collected from the

*N*

_{p}paths. Since each path is used to generate multiple training data point due to the NARX scheme,

*N*

_{T}is much greater than

*N*

_{p}. The vectors $u1=[u11,u12,\u2026,u1p]\u2208R1\xd7p$ and $a1=[a11,a12,\u2026,a1q]\u2208R1\xd7q$ are, respectively, the first

*p*and

*q*time-steps of the first training path, and the corresponding output is

*a*

_{T1}≜

*a*

_{1(q+1)}. The remaining rows of the training data in Eq. (9) can be explained similarly.

*p*+

*q*,

*Q*is the number of Gaussian components,

*w*

_{i}is the weight of the

*i*th Gaussian component, $\u2211i=1Qwi=1$, and $\varphi (xT|\mu i,\Sigma i),i=1,2,\u2026,Q$ is the PDF of the

*i*th multivariate Gaussian distribution:

*N*(·) and

*W*(·) are, respectively, the PDF of a Gaussian and a Wishart distribution, $mi$ is the posterior mean vector of $\mu i$,

*β*

_{i}is the confidence in the estimation of $mi$, and

*v*

_{i}represents the extent of effect on the estimated $Oi$, which is a scaled parameter of the covariance matrix $\Sigma i$.

The *i*th posterior mean and covariance matrix for the *i*th cluster, $\mu i$ and $\Sigma 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 $\mu i$. For a given number of Gaussian components/clusters, the GMM parameters can be estimated using the VI method.

*i*th cluster can be estimated by a famous Chinese restaurant process as follows [47]:

*N*

_{ci}is the number of samples in the

*i*th 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

*α*.

*α*, 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:

*η*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].

*j*th data sample belongs to the

*i*th cluster, which is predicted using the learned GMM, we have $\u2211i=1QP(xT,j,i|\theta )=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*th 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.

*i*th cluster is expressed as follows:

*i*th cluster, we have

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:

*i*th cluster.

*Q*NARX models as follows:

*p*and

*q*lags of the terrain elevation and vehicle acceleration, and $g^a,i(\u22c5)$ is the

*i*th 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).

### 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,\u2026,uNnew]$, where *N*_{new} 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.

*t*

_{k}as an example, the inputs $xT,k=[uk,ak]$ to the NARX models are expressed as follows:

*p*lags of terrain elevation of

*t*

_{k}and $ak$ represents the previous

*q*lags of vehicle vertical acceleration.

*t*

_{k}as follows:

*t*

_{k}are expressed as follows:

*i*th GP-NARX model at

*t*

_{k}.

*t*

_{k}using the ensemble of the NARX models, the prediction is used as input to predict the vehicle mobility at the next step

*t*

_{k+1}as follows:

The updated inputs are then used to determine the ensemble weights for *t*_{k+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.

### 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 *i*th cluster from Sec. 3.2 as follows:

*V*different training datasets by adding more lags to the data matrix of Eq. (27) as follows:

*r*th training dataset of the

*i*th 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

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:

*r*th new dataset generated from the

*i*th cluster.

*V*NARX models for the

*i*th cluster using $xT,i(r),r=1,2,\u2026,V$ and $aT,i$ as follows:

*i*th cluster trained using $xT,i(r)$ and $aT,i$.

*i*th cluster, we compute the following prediction error of the

*r*th NARX model for the validation dataset:

*l*

_{2}norm and $a^V,i(r)$ is the prediction obtained from the

*r*th NARX model as follows:

*i*th cluster after adding

*r*more lags to the original cluster.

*i*-cluster is then determined by

*Q*clusters, the NARX models, each with a different number of lags, are ensembled for prediction at

*t*

_{k}as follows:

## 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.

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.

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 × 10^{7} 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).

*S*

_{i},

*T*

_{i},

*i*= 1, …, 5) of the five paths are as follows:

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.

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.

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.

Method | Number of lags | ||
---|---|---|---|

Cluster 1 | Cluster 2 | Cluster 3 | |

DENA-DIFF-Net | 10 | 15 | 13 |

DENA-DIFF-GP | 10 | 9 | 15 |

Method | Number of lags | ||
---|---|---|---|

Cluster 1 | Cluster 2 | Cluster 3 | |

DENA-DIFF-Net | 10 | 15 | 13 |

DENA-DIFF-GP | 10 | 9 | 15 |

Figures 14–16 compare the prediction of different methods for a test path against the true vertical acceleration in scenario 1. Figures 17–19 show the prediction error of different methods for this particular path.

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.

Metrics | RMSE | Average absolute error | ||
---|---|---|---|---|

Number of training paths | 40 | 80 | 40 | 80 |

LSTM | 4.194 | 3.517 | 3.027 | 2.550 |

GP-NARX | 3.818 | 3.822 | 2.606 | 2.613 |

NARX-Net | 2.938 | 2.036 | 2.116 | 1.402 |

DENA-GP | 2.301 | 1.921 | 1.512 | 1.314 |

DENA-Net | 3.723 | 2.749 | 2.717 | 2.031 |

DENA-DIFF-GP | 2.003 | 1.910 | 1.306 | 1.236 |

DENA-DIFF-Net | 2.942 | 2.563 | 2.717 | 1.926 |

Metrics | RMSE | Average absolute error | ||
---|---|---|---|---|

Number of training paths | 40 | 80 | 40 | 80 |

LSTM | 4.194 | 3.517 | 3.027 | 2.550 |

GP-NARX | 3.818 | 3.822 | 2.606 | 2.613 |

NARX-Net | 2.938 | 2.036 | 2.116 | 1.402 |

DENA-GP | 2.301 | 1.921 | 1.512 | 1.314 |

DENA-Net | 3.723 | 2.749 | 2.717 | 2.031 |

DENA-DIFF-GP | 2.003 | 1.910 | 1.306 | 1.236 |

DENA-DIFF-Net | 2.942 | 2.563 | 2.717 | 1.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.