Using a Bayesian-Inference Approach to Calibrating Models for Simulation in Robotics

This contribution is motivated by the increasing use of simulation in the design of autonomous robots and vehicles. Simulation reduces design time, cuts costs, and leads to a safer development cycle and a more robust engineered solution [1]. For instance, one might be interested in choosing the model parameters that would allow simple two-degree-of-freedom (DOF) or eight-DOF vehicle models to be used effectively for model predictive control in an autonomy stack [2]. This poses a parameter calibration problem—data obtained from a vehicle is available and the task is to calibrate an expeditious model that is a good proxy for the dynamics of the actual vehicle. Traditionally, the parameter calibration problem for multibody systems has been posed as an optimization problem by coupling the integrator with an optimization procedure [3–7]. With the exception of Ref. [7], these parameter estimation approaches provide single point estimates for the parameters. They require direct access to the model being calibrated since either forward sensitivities or associated adjoint coefficients need to be estimated in processes that require derivative information and might have large memory footprint. Finally, there might also be an expectation of calibration data being smooth and clean, a requirement that might not be met. An approach that takes a different tack is presented in Ref. [7], where inertial parameters of a delta robot are calibrated using a least squares estimator. Its limitation is that the parameters sought must appear linearly in the model, which happens to work for inertial parameters in multibody dynamics. Against this backdrop, the parameter calibration problem is approached herein in a Bayesian framework, a choice motivated by three observations: it works with noisy data; it is not intrusive, in that no access to the equations of the model or their sensitivities is required; and it produces confidence bounds for the results.

it is often difficult to directly evaluate that formula, particularly for the evidence term in the denominator when the posterior distribution is on a high-dimensional space. Therefore, Markov-chain Monte Carlo (MCMC) methods [9–11] have been widely used to sample from the posterior distribution. Then, one does not have a closed form for the posterior, but the samples obtained from MCMC can be used to approximately evaluate any expectation with respect to the posterior distribution, which is after all the task of practical interest.

The MCMC algorithm is setup to produce samples of the distribution. It is suitable for the task of Bayesian inference, as MCMC algorithm can be used when the distribution is only known up to a normalizing constant. This is done by trying to see whether a particular choice of parameters adequately explains the observed data, which calls for running a simulation with the candidate set of parameters. Sequential Monte Carlo (SMC) samplers [12,13] have been explored as an efficient alternative to classical MCMC. In general, SMC samples from a sequence of distributions ${πn(θ)}n=1N$⁠, which could arise from the conditional distributions $πn(θ)=p(θ|y1,…,yn)$ with a sequence of data ${yn}n=1N$⁠, or from a sequence of artificial intermediate distributions [14] that starts from an easy-to-sample distribution, e.g., the prior distribution, and gradually converges to the posterior distribution. Compared to MCMC, SMC is easier to parallelize and less likely to be trapped in local modes.

This contribution is organized as follows. Section 2 covers the fundamentals of the Bayesian inference process used. It discusses the Bayesian framework; the SMC sampler that helps one draw samples from the posterior distribution; and the probabilistic programing framework PyMC [15] used to code in Python the described inference framework. In Sec. 3, the Bayesian machinery is used to answer this question: for simulating vehicle dynamics over short time intervals, can a high-fidelity high-mobility multipurpose wheeled vehicle (HMMWV) Chrono model [16,17] be replaced by a simple model, with eight degrees-of-freedom, that runs faster given that its evolution is captured by a small set of differential equations? A discussion of the Bayesian calibration approach embraced makes up Sec. 4, in which we highlight limitations of the methodology and our use of it. Directions of future work and a series of conclusion wrap up the paper.

where $p(θ)$ is the prior distribution that carries some existing knowledge about the unknown parameters θ. We use SMC to draw samples from the posterior distribution $p(θ|y)$⁠.

where ${ϕn}$ denotes a sequence of hyperparameters that starts from 0 and eventually converges to 1. In practice, SMC employs sequential importance sampling to sample each p_n based on an importance distribution $ηn(θ)$⁠. More specifically, it often starts with $η0(θ)=p(θ)$ and then sequentially updates η_n according to $ηn−1$ and the local random walk kernel $Kn(θ,θ′)$ for $n≥1$ [13]. Compared to the traditional MCMC approach, SMC has a few advantages, including better sampling of multimodal distributions and more straightforward parallel implementation that improves the computational efficiency.

The workflow for Bayesian calibration of the vehicle model is presented in Fig. 1. Within this workflow, a Python package called PyMC [15] is used to sample from the posterior distribution.

Due to the lack of experimental data, we use Chrono and its Chrono::Vehicle module [17] to generate synthetic high-fidelity “experimental” data. Since real data collected in physical experiments is noisy, we add Gaussian white noise on top of the data produced by the simulator. More details about the data generation process are provided in Sec. 3.1. We also define an 8DOF model whose unknown parameters $θ$ we aim to calibrate are based on the synthetic data. The goal is to identify a choice of 8DOF model parameters that makes its dynamics a good proxy for that of the “real” vehicle in Chrono.

We use a Gaussian likelihood function, i.e., the squared exponential term in Eq. (3). The likelihood function demands the evaluation of $G(θ)$⁠, which involves running the 8DOF model and gathering model responses. For the prior distributions of the unknown parameters θ, we assume a uniform distribution within an empirically chosen range for each unknown parameter. The samples obtained from PyMC, along with various sampler statistics, can be used with ArviZ [18] to do exploratory analysis of PyMC results.

In robotics, one needs an expeditious model for tasks such as state estimation, path planning, or model predictive control. While typically task specific, these models are similar in one regard: they are expected to run much faster than real-time. These observations motivate the case study in this section: producing a simple vehicle model that runs fast and whose time evolution (dynamics) is a good proxy for the dynamics of a complex vehicle of interest. The vehicle of interest is simulated in Chrono to produce synthetic training/ground truth data and is subsequently used to calibrate a simplified vehicle model. Upon completing the calibration process, we assess how close the dynamics of the simple model comes to that of the vehicle of interest.

In this contribution, the simple model is an 8DOF model, while the vehicle of interest is a U.S. Army HMMWV. The Chrono high fidelity model used for simulating this HMMWV captures all relevant dynamics for a wide range of on- and off-road mobility applications. Constructed using the template-based approach enabled by Chrono::Vehicle [17], this model includes accurate models of the HMMWV independent double A-arm suspensions, a Pitman-arm or rack-and-pinion steering mechanism, a full dynamic model of the driveline (including models of the central and axle differentials), as well as models of the engine (based on power-torque curves including engine losses) and of the automatic transmission (torque converter and gear box). The main components of the resulting multibody system, which includes both kinematic joints and bushings, are shown in the simulation snapshots of Fig. 2. Chrono::Vehicle offers a variety of different tire models, from the commonly used semi-empirical Pacejka, Fiala, and TMeasy, to full FEA models using ANCF or corotational elements.

where $hrc≡hrcfb+hrcraa+b$ and $Ĵx≡(Jx+mhrc2)$⁠.

In these equations, h_rcf and h_rcr are the vertical distances of the front and rear roll centers below the sprung mass center of mass (CM), and thus h_rc is the vertical distance from the sprung mass CM to the vehicle roll center. It should be noted that Eq. (7) for the roll degree-of-freedom is written by considering moments acting about the vehicle roll center rather than the sprung mass CM, and thus the roll inertia of the sprung mass about the vehicle roll center $(Jx+mhrc2)$ is considered in Eq. (7) [20,22]. The forces F_xgij and F_ygij are the global longitudinal and lateral forces at the four tire contact patches and the subscript “ij” denotes lf (left-front), rf (right-front), lr (left-rear), and rr (right-rear).

where $δ$ is the steering angle at the tire (only the front wheels can be steered), $ugij$ and $vgij$ are the global longitudinal and lateral velocities of the tire contact patch, $ωij$ is the angular velocity of the tire, and $rij$ is the instantaneous tire radius. For instance, the velocities of the tire contact patch at the right front contact patch are given by $ugrf=u+ωzcf/2$ and $vgrf=v+ωza$⁠.

where $ff≡mufhufcf+mb(hcg−hrcf)cf(a+b)(v˙+ωzu)$⁠, and $fr≡(murhurcr+ma(hcg−hrcr)cr(a+b))(v˙+ωzu)$⁠, and $fc≡(mhcg+mufhuf+murhur)(u˙−ωzv)2(a+b)$⁠.

Finally, using Eqs. (8) and (9), with the comprehensive slip $sα$ and the coefficient of friction $U$ as defined in conjunction with the tire vertical forces $Fz$⁠, the critical longitudinal slip ratio $scritical$ is calculated as in Ref. [23]. Specifically, if the absolute value of the longitudinal slip ratio $|s|≤scritical$⁠, then the tire is in elastic deformation state and the local longitudinal tire force is given by $Fxtij=Cxjsij$⁠. Otherwise, the tire is in the complete sliding state and the local longitudinal tire force is $Fxtij=sgn(s)(Fx1ij−Fx2ij)$⁠, where $Fx1ij=Uij|Fzij|, Fx2ij= |(UijFzij)24sijCxj|$⁠, and $Cxj$ is the longitudinal tire stiffness of the front/rear wheels.

When the tire is in the elastic deformation state, the lateral tire force is given as (dropping tire index ij) $Fy=−U|Fz|(1−H3)sgn(α)$⁠, where $H=1−Cy| tan(α)|3U|Fz|$⁠. In the complete sliding state, the lateral tire force is given as $Fy=−UFzsgn(α)$⁠, where $Cy$ is the lateral tire stiffness, and $Fxg$ and $Fyg$ are obtained from $Fxt$ and $Fyt$ using a coordinate transformation.

The parameters of the 8DOF vehicle model are summarized in Table 1. Therein, a distinction is made between a tire and a wheel: the former is regarded as a force elements, the latter combines a tire and a wheel hub.

The calibration of the 8DOF model using HMMWV synthetic data is split into three stages which concentrate, respectively, on the longitudinal dynamics, lateral dynamics, and tire rolling resistance. The longitudinal dynamics stage of the calibration process uses the vehicle longitudinal velocity along with the front and rear angular velocities of the wheels. The lateral dynamics stage uses the lateral velocity of the vehicle, as well as its yaw rate, roll angle, and roll rate. The tire rolling resistance stage of the calibration uses the vehicle's longitudinal velocity, but with a different vehicle maneuver.

A Gaussian white noise is added to the synthetic data. Although this is not a prerequisite to using Bayesian calibration approaches, noise is added to demonstrate that the data used need not be smooth, which is the case when dealing with real world scenarios. The standard deviation of the noise for each data vector is $σu=0.1,σωf=1, σωr=1, σv=0.05, σψ˙=0.02, σϕ=0.005$ and, $σϕ˙=0.002$⁠.

In order to calibrate the 8DOF model, we set out to first identify the parameters of interest that, in practice, are unknown or difficult to quantify using simple experiments. Measurable parameters such as the sprung/unsprung mass of the vehicle (⁠$m,muf,mur$⁠), the position of the CM (⁠$h,a,b$⁠) and the roll center (⁠$hrcf,hrcr$⁠), the track width (⁠$cf,cr$⁠) and the nominal tire radius (⁠$r0$⁠) are not required to be calibrated. Similarly, sprung mass inertias (⁠$Jx,Jz,Jxz$⁠), wheel inertias (⁠$Jw$⁠) and the tire vertical stiffness (⁠$ktf,ktr$⁠) are obtainable through fairly simple experiments. Thus, obtaining distributions over these parameters is not of much interest as they can be quantified fairly accurately without the need of data from vehicle maneuvers.

The parameters in Table 2 are considered known parameters, where the HMMWV model in the Chrono::Vehicle simulation and the 8DOF model use the same set of values. This leads to the following unknown parameters for which we obtain posterior distributions: C_xf, C_xr, C_yf, C_yr, $kϕf, kϕr, bϕf, bϕr$⁠, rr.

We first calibrate the 8DOF model for longitudinal dynamics, i.e., we aim to make the 8DOF model perform like the HMMWV for longitudinal acceleration and deceleration vehicle maneuvers.

To generate calibration data, a simple straight line acceleration test is used. A normalized throttle input increasing from 0 to 0.5, see Fig. 4, is applied to the HMMWV vehicle. The HMMWV vehicle model uses a simple powertrain with a single gear ratio along with a maximum torque and speed map (much like a DC motor). The powertrain is connected to a four-wheel driveline which specifies the intermediate shaft inertias and the gear ratios of the differentials. The 8DOF model uses an identical powertrain and thus takes the same normalized throttle input. However, it does not have a driveline and so the torque from the powertrain is directly applied to the four wheels. It should be noted however that additional gear ratios of the differentials in the HMMWV vehicle are taken into account in the 8DOF model. Also, the 8DOF uses a simple tire model that has poor behavior when the vehicle starts from rest. Thus, the simulation of the 8DOF model starts from the state reached by the HMMWV at time 1.0 s into its simulation; at that point, the state of the 8DOF model (initial conditions) is made to match the state of the Chrono::Vehicle (red vertical line in Fig. 4). Although this does not play a role in the calibration effort, it remains a problem to be addressed in the future.

Out of the unknown parameters, the ones identifiable through a longitudinal acceleration experiment are the front and rear longitudinal tire stiffness $Cxf$ and $Cxr$⁠. According to the empirical knowledge about $Cxf$ and $Cxr$⁠, we provide uniform-distributed prior distributions within a range, i.e., $Cxf∼U(1000,50000)$ and $Cxr∼U(1000, 50000)$⁠. The values 1000 and 50,000 are chosen to cover a reasonable range where we expect $Cxf$ and $Cxr$ to lie.

In addition to $Cxf$ and $Cxr$⁠, we also sample the standard deviations of the noise in the three data vectors of interest: the longitudinal velocity of the vehicle, the left front wheel angular velocity and the left rear wheel angular velocity. Although in this study the standard deviation of the noise added is already known and can be directly used to evaluate the Gaussian likelihood, this may not be the case for real world data. Thus, it makes sense to sample the standard deviations together with other unknown model parameters as it can be particularly useful to build measurement error models when the noise is in fact unknown [25,26]. We choose a half-normal prior distribution as this incorporates our prior knowledge of the noise being small and always positive. The priors are of the form $σu∼H(0.1), σωf∼H(1)$⁠, and $σωr∼H(1)$⁠.

The SMC sampler from PyMC is used for the longitudinal dynamics calibration. Eight SMC chains are run in parallel with each chain drawing 1000 samples. The SMC kernel used is Independent Metropolis Hastings. This means that in each chain, at each stage, the required (algorithm determined) steps of random walk (using Metropolis Hastings kernel) are run from each sample of SMC. The acceptance probability for each Metropolis Hastings chain is set to 0.9. Running multiple chains is essential for obtaining valid inferences from iterative simulations of finite length [27]. Having multiple chains is also critical in gauging chain convergence using diagnostics such as split-$R̂$ [28]. In addition, multiple chains are more likely to reveal multimodality and poor adaptation or mixing [28]—scenarios that are particularly common for engineering applications. The minimum number of chains recommended in practice is four [28].

In practice, it is recommended that the Effective Sample Size (ESS) be greater than 400 [28]. It should be noted that ESS is not the same as the number of draws but rather a quantity of interest that captures how many independent draws contain the same amount of information as the dependent samples obtained by the Markov chains [28]. In this work, we have 1000 draws per SMC chain to ensure that ESS is greater than 400.

The posterior distribution obtained for the parameters $Cxf$ and $Cxr$ is plotted using ArviZ in Fig. 5. The posterior distribution for the standard deviation of the synthetic Gaussian white noise added to the Chrono::Vehicle data is shown in Fig. 6. The posteriors obtained for the standard deviation of the noise are in the range of the values seen in Sec. 3.1. It is noted that the standard deviation obtained through the posterior is slightly larger than that supplied due to the difference in physics captured by the 8DOF and HMMWV models. Although these posteriors are useful as they show a distribution over the parameters, what is more insightful is to analyze the response of the 8DOF model with samples drawn from these posteriors—this will show one how the fact that we are not certain in the value of the calibration parameters reflects in the response (dynamics) of the 8DOF model.

We show the improvement of the posterior compared to the prior, by drawing samples of $Cxf$ and $Cxr$ from them, and assessing the fit between the 8-DOF response and that of the HMMWV. We first draw samples from the prior distribution, which represents our noninformative belief about the 8-DOF model prior to seeing any HMMWV data. Therefore, the 8DOF model can assume very different trajectories, see Figs. 7(a)–9(a). These figures are obtained by taking 100 random samples out of the prior distribution and plotting the 8DOF response. Bringing in the HMMWV data shapes our belief about the parameter values and leads to a posterior distribution. Note that drawing samples from the posterior is as straightforward as picking random samples from the Markov Chain that was used to generate the posterior. Thus, we pick 100 random samples of $Cxf$ and $Cxr$ from their joint posterior distribution obtained from the trace of the chain (to be discussed shortly) and plot 100 responses of the 8DOF model, see blue line in Figs. 7(b)–9(b).

These figures show that drawing samples from the posterior produces a response that is much narrower and closer to the noisy HMMWV data. Plotted along with the 100 responses are also the posterior/prior expectation and the posterior/prior mean. The posterior/prior expectation is obtained by taking the mean of the 100 responses whereas the posterior/prior mean is the response obtained by plugging in the mean values of the posterior/prior distributions of the calibrated parameters. The latter two lines are very close to each other, which suggests that the posterior mean is a good single point estimate of the posterior as it produces similar results to the expectation of the model response [29].

where $G(θ)i$ is the 8DOF model response at time-step i, $yi$ is the noisy data at each time-step and $nt$ is the number of time steps. The mean RMSE for both the prior and the posterior is reported in Table 3. This shows that the response of the 8DOF model with parameters drawn from the posterior produces a much tighter fit than with parameters drawn from the priors, a confirmation that the Bayesian framework uses the HMMWV data meaningfully when it comes to improving the quality of the 8DOF model.

“While MCMC, as well as more general iterative simulation algorithms, can usually be proven to converge to the target distribution as the number of draws approaches infinity, there are rarely strong guarantees about their behavior after finite time” [28]. In order to draw conclusions from the posterior distribution, it is important to first assess how well the sampler is approximating the specified posterior distribution. One way to do this is by visually inspecting the sample paths of the eight chains via trace plots. The trace plots for the longitudinal dynamics calibration can be seen in Fig. 10 (right). Therein, it can be seen that the eight chains, which start from random points in prior parameter space, converge to a similar posterior. This shows that the chains are relatively insensitive to the starting point, thus showing geometrical ergodicity, a critical property for the central limit theorem to hold for approximate posterior expectations [28]. Generally, however, a qualitative assessment of chain convergence using visual inspection is often times considered “notoriously unreliable” [30]. Thus, we also calculate the split-$R̂$ metric using ArviZ as this is the most commonly used convergence diagnostic in practice [28]. “Roughly speaking, the split $R̂$ metric compares between and within chain estimates for model parameters” [28]. If the chain has not mixed well, the variance of all the chains mixed together will be higher than the variance of individual chains and the split-$R̂$ will be larger than 1. It is recommended to only use samples if the split-$R̂<1.01$ [28]. Table 4 shows that the split-$R̂$ for all our parameters is less than 1.01. In addition to this, we also calculate the mean and standard deviation of the Monte Carlo Standard Error (MCSE), which is a statistic that helps quantify the sampler error [28,31]. It can be seen from Table 4 that the mean and standard deviation of the MCSE for all our parameters calculated using ArviZ is much smaller than the scale of the parameters, which is reassuring.

Further, we also plot the ESS evolution plot [28] as it gives a “scale-free” measure of information and is particularly useful when diagnosing sampling efficiency. In particular, we plot the evolution of both the bulk-ESS—the effective sample sizes in the “bulk” (5% to 95% quantiles) of the posterior, as well as the tail-ESS—the effective sample sizes in the “tails” of the posterior (minimum of 5% and 95% quantiles). “For a well explored distribution, we expect both the bulk-ESS and tail-ESS measures to grow linearly with the total number of draws S, or, equivalently, that the relative efficiency (ESS divided S) is approximately constant for different values of S” [28]. Figure 11 indicates that both the bulk-ESS and the tail-ESS for $Cxf$ and $Cxr$ increase linearly with the number of draws and reach a much higher value than the recommended value of 400 [28]. To conclude, the SMC sampler approximates the posterior distribution well.

To calibrate the parameters that dictate the lateral dynamics of the 8DOF model, the following Chrono HMMWV simulation is run to generate synthetic data: first, the vehicle is accelerated for 7 s to bring it to a speed of approximately 40 mph. For the next 3.7 s, the normalized steering input at the steering wheel is increased from 0 to 0.2, linearly (positive value corresponds to left steer). The throttle and steering input can be seen in Fig. 12. The HMMWV uses a rack-and-pinion steering mechanism. The 8DOF model is provided the state of the HMMWV at 7 s and is then given the same normalized steering input. Only data from 7 to 10.7 s is used for the calibration. This is done to separate the longitudinal dynamics from the lateral dynamics during the training stage. Unlike the HMMWV, the 8DOF model does not include a model of the steering mechanism. To ensure a similar steering angle at the front wheels for both the 8DOF and HMMWV for the same normalized steering input at the steering wheel, we take the average of the maximum front wheel steering angles of the HMMWV and set it as the maximum steering angles of the front wheels of the 8DOF model. A left/right average is used to accommodate Ackerman steering which is captured by the Chrono model, but not the 8DOF model. With this, the same normalized steering wheel angle is supplied to both the 8DOF model and HMMWV resulting in similar steering angles at the wheels.

The parameters calibrated during this stage of the process are $Cyf,Cyr,kϕf,kϕr,bϕf,bϕr$⁠, along with the standard deviation of the noise of each of the data vectors of interest—$σv,σψ˙,σϕ,σϕ˙$ (corresponding to lateral velocity, yaw angle, roll angle and roll rate). For the 8DOF model parameters, the same priors from Eq. (14) for $Cyf,Cyr,kϕf$ and $kϕr$ are used. However, for $bϕf$ and $bϕr$ we sample a single vehicle damping coefficient $bϕ$ and then set $bϕf=bϕr=0.5 bϕ$⁠. This choice is further explained in Sec. 4.2. The mean of the posterior obtained in the longitudinal dynamics calibration for $Cxf$ and $Cxr$ is used. A half normal prior was chosen for the noise parameters. The scale of the prior is chosen according to the scale of the data as $σv∼H(0.05), σψ˙∼H(0.05), σϕ∼H(0.005)$⁠, and $σϕ˙∼H(0.005)$⁠.

The SMC sampler described for longitudinal dynamics is used herein; also, the SMC settings stay the same.

The posterior distribution obtained for the parameters of interest is shown in Fig. 13. For the noise parameters, see Fig. 14. The posteriors are very narrow with a concentrated 94% HDI, which is highly desirable. It can also be seen that $bϕf$ and $bϕr$ have identical posteriors, which is expected based the way they are evaluated via $0.5 bϕ$⁠. Also, the standard deviation of the noise obtained is in the same range as those supplied during data generation (see Sec. 3.1). Similar to the longitudinal dynamics calibration, we will first discuss the 8DOF response using the posterior distribution and then discuss chain convergence statistics.

A set of 100 responses using the prior distribution are summarized in Figs. 15(a)–18(a). We pick 100 random samples of $Cyf,Cyr,kϕf,kϕr,bϕf$ and $bϕr$ from the joint posterior distribution to plot 100 responses of the 8DOF model after bringing to bear the HMMWV synthetic data. The 100 responses are seen in Figs. 15(b)–18(b). To quantify the posterior fit, we use again the mean of the RMSE between the 100 response lines and the noisy data. The results in Table 5 confirm that the average RMSE of the responses using samples from the posterior is approximately 10 times smaller than using samples from the prior across all data vectors.

An inspection of the trace plots for the eight chains in Fig. 19 reveals that the eight chains converge to a similar posterior. Table 6 shows that the split-$R̂$ is less than 1.01 for all parameters. The mean and standard deviation of the MCSE is also small compared to the scale of each parameter and is thus acceptable. The bulk-ESS and tail-ESS evolution plots are shown in Fig. 20. The values obtained are well above the required minimum of 400 and also grow linearly with the number of draws.

The SMC sampler is thus able to approximate the posterior distribution well. The posteriors over the parameters $Cyf,Cyr,kϕf,kϕr,bϕf$ and $bϕr$ also produce a well calibrated lateral dynamics response. We now finally calibrate the rolling resistance of the tires using a similar approach. We then test our calibrated 8DOF model against different combinations of vehicle maneuvers that were not seen in the “training” data.

This stage of the calibration process draws on a basic maneuver shown in Fig. 21. We first accelerate the HMMWV in Chrono to a speed of about 90 km/h. We then let the vehicle coast. Note that there is no aerodynamic drag included in the Chrono HMMWV model for the purpose of these simulations and therefore the slow down comes from the rolling resistance of the tires. The 8DOF model is set to the state of the HMMWV at the 5.6 s time-mark as shown by the dotted red line in Fig. 21. The Longitudinal velocity is selected as the data vector of interest that is used to drive the calibration process.

Since we expect the rolling resistance coefficient values of both the 8DOF tire model and the HMMWV tire model to be similar, we can provide a narrow prior around the HMMWV's value (Eq. (14)). We also provide a prior over the standard deviation of the longitudinal velocity noise as $σu∼H(0.1)$⁠. This is the same prior for the noise used in the longitudinal dynamics component of the calibration effort.

The SMC sampler described earlier, with the same settings, is re-used for the rolling resistance calibration.

The posterior distribution obtained for the rolling resistance coefficient is shown in Fig. 22. For the posterior distribution of the standard deviation of the longitudinal velocity noise, see Fig. 23. The posteriors visually look narrow with a thin 95% high density interval (HDI).

We perform a similar exercise to evaluate the response of the 8DOF postcalibration—draw 100 samples and plot the 8DOF system response, see Figs. 24(a) and 24(b). The RMSE results are provided in Table 7. Compared to the lateral calibration, the agreement between sampled simulations is less prominent, nonetheless the posterior mean and expectation show significant improvements over the prior.

We once again evaluate the trace (Fig. 25), the split-$R̂$ (Table 8) and the ESS (Fig. 26). All the convergence diagnostics suggest a successful calibration process, as the SMC sampler is able to approximate the posterior well.

The calibration steps described in Subsections 3.3–3.5, which drew on data collected from the HMMWV over a small number of maneuvers, have lead to a calibrated 8DOF model. The question is whether the 8DOF model has become a good proxy for the HMMWV. To answer this, we subjected both the 8DOF model and HMMWV to a set of maneuvers that were not used during the calibration stage. We use the posteriors from the previous three calibration stages and plot 100 responses of the 8DOF model (Figs. 28–35). The control inputs chosen are illustrated in Fig. 27. We first accelerate the vehicle by applying a steeper ramp throttle from 0 to 3.5 s. We then let go of the throttle and apply a left step steer simultaneously. Then we apply a right ramp steer for approximately 3 s. We then again bring back the steer to 0 and apply brakes for the last 2 s. We set the state of the 8DOF model as the state of the Chrono::Vehicle at the 1 s time-mark of the simulation. This is because, as mentioned before, the 8DOF has difficulties starting from rest.

Similarly we also calculate the mean-RMSE for each of the data vectors (Table 9). Overall, the responses of the 8DOF model is much more concentrated in and around the Chrono data postcalibration.

Figure 32 indicates that the 8DOF model roll angle during the first 3.5 s is relatively constant and close to zero. The evolution of the roll angle in the data is however decreasing. This is because the HMMWV has a small amount of roll even during straight-line acceleration due to the reaction torque of the engine and transmission on the chassis. This cannot be captured by the 8DOF model through any combination of parameters. Similarly, the 8DOF model is also not able to match the transient roll rate during the step steer maneuver (Fig. 33). The 8DOF roll rate does however stabilize to a value much closer to the data.

The 11.6 s scenario described above will be used for a time comparison to gain a rough understanding of the 8DOF model execution speed. A time-step of $5×10−3$ s is used for the 8DOF model; the HMMWV simulation in Chrono uses a time-step of $2×10−3$ s and with an Euler linearized implicit integrator [16,32]. The 8DOF model uses a half-implicit integrator [33]. The simulations were performed on a Intel(R) Core(TM) i7-7500 U CPU @ 2.70 GHz and the results are provided in Table 10. Note that the HMMWV is simulated at a real-time factor (RTF) of about 0.41. Unsurprisingly, the 8DOF model runs more than an order of magnitude faster than the HMMWV, at an RTF of 0.03. The 8DOF simulation can be further accelerated by re-implementing the model in C/C++, a future task that falls outside the scope of this work.

There are multiple samplers available to solve a Bayesian calibration problem. Herein, three of them are employed to understand (a) how close they are in terms of the posterior distributions produced and (b) how efficient they are. Besides SMC, we considered the Metropolis-Hastings sampler [34], which is one of the first algorithms employed for MCMC, and the No-U-Turn Sampler (NUTS) [35], which works using the Hamiltonian Monte Carlo algorithm [10].

The Metropolis-Hastings algorithm is an adaptation of a random walk with an acceptance/rejection rule to converge to the specified target distribution [30]. While simple to implement, it is known to not perform very well in high dimensional parameter spaces due to its random-walk style of exploration [36].

Hamiltonian Monte Carlo (HMC) is a MCMC algorithm that avoids the random walk behavior and sensitivity to correlated parameters that plague many MCMC methods. It accomplishes this taking a series of steps informed by first-order gradient information [35]. This allows the HMC methods to explore high dimensional spaces much more efficiently by producing less correlated samples. NUTS is an extension of HMC. Among other things, it eliminates the need to set the step size and the number of steps, two crucial parameters for the computational efficiency of the HMC algorithm.

To gain insights into their performance, we perform the lateral dynamics calibration once again using all three samplers. The same priors are used for the parameters, and the sample acceptance probability is set to 0.9 for all samplers. NUTS and M-H differ from SMC in their need to have tuning samples. Tuning is the process of discarding the initial few samples drawn by the sampler since they typically do not belong to the converged target distribution [29]. The default number of tuning samples in PyMC is 500 and we use the same in this study. NUTS requires us to take the first order gradient of the likelihood functions; we do this via finite differences.

An effective sample size (ESS) per second metric is used to measure sampler efficiency. The split-$R̂$ obtained for the different parameters is used to compare chain convergence and chain autocorrelation. Finally, it is insightful to check whether all three samplers arrive at a similar posterior albeit through different means.

The posterior obtained using the NUTS sampler can be seen in Figs. 36 and 37. The posterior is very similar to that produced by SMC in Figs. 13 and 14. However, the M-H algorithm does not produce a converged chain after 1000 draws (Fig. 38). We thus increase the number of draws per chain from 1000 to 10,000. Comparing results from 10,000 draws of M-H and 1000 draws of NUTS and SMC is still fair as the metric used for the comparison is scale free. With 10000 draws per chain, the M-H algorithm produces a converged posterior, see Figs. 39 and 40. This posterior is, again, almost identical to the posteriors produced using the NUTS and SMC sampler.

Next, rather than comparing the ESS per second and split-$R̂$ of all the parameters individually between the three methods, we take the average across all the parameters. Table 11 shows that SMC outperforms NUTS and M-H in all the metrics. Although we ran M-H for 10000 draws per chain, we only obtained a total of 5155 effective samples which shows that the M-H algorithm struggles to produce noncorrelated, independent samples. The NUTS algorithm does well to produce independent samples. However, as it requires the gradient of the likelihood for each draw, it is severely slowed down. The performance of NUTS could thus significantly improve if the model was implemented in a way that allows for techniques such as automatic differentiation (AD) to take derivatives rather than having to rely on finite differences.

This subsection highlights one typical scenario in which the Bayesian approach described cannot help. While calibrating the 8DOF model with the HMMWV model in the lateral direction, we had to reduce two of the model parameters, namely, the front (⁠$bϕf$⁠) and rear (⁠$bϕr$⁠) roll damping coefficient into a single vehicle roll damping coefficient $bϕ$ by setting $bϕf=bϕr=0.5 bϕ$⁠. This was done because initial calibration attempts showed that $bϕf$ and $bϕr$ are not uniquely identifiable using the HMMWV synthetic data available. This can be seen in the posterior pairwise plot between the two parameters (plotted after performing Kernel Density Estimation), see Fig. 41, as well as in the marginal posterior distribution which is not very sharp, see Fig. 42. Kernel Density Estimation (KDE) is a nonparametric way to estimate the probability density of a random variable based on kernels as weights. In Fig. 41, we provide a cross section of the joint—probability distribution of our two parameters where the more yellowish the hue, the higher the probability.

This can be explained using Eq. (7) where $bϕf$ and $bϕr$ occur together, as a summation, to provide the coefficient to the roll rate. Although in the vertical load transfer equations (Eqs. (10)), $bϕf$ and $bϕr$ occur independently, their low magnitude when combined with the roll rate is almost negligible in comparison with the rest of the terms. This also explains why the problem is not seen with $kϕf$ and $kϕr$⁠, which have a sufficiently distinguishable effect on the vertical forces. This served as a reminder of a classical scenario in which the Bayesian calibration is not capable of providing accurate estimation of each individual parameter—when several parameters are nonidentifiable based on the data. Such kind of nonidentifiability would be manifested in the strong correlation between certain dimensions of the joint posterior distribution as illustrated in Fig. 41.

This contribution outlines a Bayesian methodology to calibrate models used in robotics and autonomous vehicle applications, wherein there is a need for expeditious models to estimate state, plan motion, or determine control policies. The calibration approach discussed: (i) is robust in that it handles noisy data; and, (ii) it provides levels of confidence in the model parameters obtained via probability distributions. In other words, the approach does not produce one parameter set, but rather an infinite number of possible parameter sets, each labeled by a probability of explaining the experimental data observed. We relied on a Sequential Monte Carlo technique to draw samples from the joint posterior distribution of the model parameters, a process that in combination with classical Monte Carlo estimation opens the door to predicting with quantified uncertainties other quantities of interest tied to the dynamic evolution of the simplified model, e.g., most likely trajectory, most likely speed evolution.

Looking ahead, we plan to apply the same Bayesian inference framework to terramechanics calibration, wherein the goal is to use a high-fidelity approach to generate data subsequently used to calibrate a lower-fidelity but faster terramechanics model. Two cases are of interest—using a continuum representation of the terrain [37] to calibrate a lower-fidelity Chrono terrain [38,39]; or using a fully resolved, discrete representation of the terrain [40], to produce calibration data for a lower-fidelity continuum representation of the terrain [37]. Second, we plan to employ automatic differentiation to better gauge the potential of the Hamiltonian Monte Carlo approach, which is currently hindered by the lack of gradient information. Finally, we are in the process of using a real vehicle to generate calibration data, instead of relying on synthetic data generated in Chrono.

The software used to generate the results reported herein is available on GitHub as open source and distributed under a permissive BSD3 license [19].