Slip estimation model for traversability-based motion planning of cargo rover on extraterrestrial surface

As part of the robotics technologies required for In-situ resource utilization (ISRU), the development of cargo rovers for transporting resources is needed. However, these cargo rovers have unique technical challenges that differ from conventional exploration rovers, including the need to traverse rough terrains with their varying mass due to transporting payloads. Moreover, research addressing these challenges has been limited, and the relevant technologies have not been fully established. To address these challenges, this paper proposes a parametric model for estimating wheel slippage. The model is formulated as a function of four input parameters: slope angle, rover heading angle, payload mass, and wheel angular velocity, and is applicable to resource-transporting rovers with varying mass. Additionally, the use of a parametric model reduces computational load, which offers advantages for onboard implementation. The proposed estimation model was quantitatively evaluated by comparing datasets obtained from multi-body dynamics analysis. This paper also introduces a new traversability assessment model which incorporates the proposed slip estimation model. We demonstrated the proposed model by integrating it into a sampling based motion planning. The simulation result of the motion planning show that the planner with our model can generate safer motions and enables the rover to reach the target regardless of the cargo payload.

1 Introduction

In recent years, there has been active discussion regarding long-term exploration activities of planets and moons. The Artemis program, for example, aims to establish technologies for sustainable human presence in space, with the Moon as a base for future Mars missions. A key component of Artemis is In-Situ Resource Utilization (ISRU), which seeks to use local resources for producing water, fuel, and construction materials (Smith et al., 2020; Sanders and Kleinhenz, 2024; Werkheiser et al., 2024). Achieving ISRU will require advanced robotics, particularly rovers capable of transporting resources to designated locations. Additionally, future operations may involve multiple rovers, and there is an increasing expectation for rovers to move autonomously while considering traversability. However, implementing resource-transporting rovers presents several technical challenges not encountered with conventional exploration rovers. First, there are differences in the environmental conditions in which rovers are operated. One of the most critical risks in rover operations is the vehicle getting stuck in sand, which is extremely difficult to recover from via remote control, and in some cases, this could directly lead to mission failure. Conventional exploration rovers have been operated conservatively, avoiding areas prone to slipping, due to the unique circumstances of space missions, such as being one-time missions and the fact that human repairs are not feasible. However, when considering applications for resource transportation, the roads in excavation areas are unpaved, and repeated rover movements and excavation activities can alter the terrain, impeding the rover mobility. Autonomous systems capable of safely operating in such harsh environments have not yet been demonstrated in space, making this a highly challenging technical task. Next, transporting resources means that the rover’s weight varies dynamically. The mass properties are closely linked to the rover dynamics, making the rover traversability significantly varies.

Experiments have been conducted to investigate the traversability of rovers, revealing that slip characteristics are influenced by a wide range of factors (Ishigami et al., 2007; Chhaniyara et al., 2012; Sutoh et al., 2010). Based on the result, numerous non-parametric models have been proposed for estimating slip characteristics (Sevastopoulos and Konstantopoulos, 2022). Non-parametric models have the advantage of being able to flexibly represent complex real-world phenomena. Gonzalez et al. (2018) proposed a machine learning algorithm to estimate slip characteristics using proprioceptive sensors such as RTK-GPS and IMUs, without relying on environmental information such as terrain slope. A major advantage of this method, which solely uses proprioceptive sensors, is that it does not require computationally expensive processes like Visual Odometry and is not affected by environmental conditions such as optical conditions. Another end-to-end approach involves a neural network model that estimates slip ratio from three parameters: translational speed, slope angle, and heading angle (Sakayori and Ishigami, 2021). This model is designed to output both the slip ratio and energy consumption simultaneously. One key feature of neural network models is their ability to integrate various types of information into a single architecture for estimation. However, as the phenomena become more complex and the dimensionality of the model increases, a large amount of training data is required. This raises concerns, especially in the context of planetary exploration, about whether sufficient datasets can be prepared. Furthermore, the low interpretability of black-box models can lead to issues such as overfitting. Consequently, ensuring that these models do not exhibit unexpected behavior requires exhaustive validation across a wide range of scenarios, which makes their application to real missions challenging.

The main contributions of this paper are as follows:

• A new parametric slip estimation model that incorporates parameters such as mass variation and wheel velocity.

• A new traversability assessment model that integrates the above slip estimation model.

In this paper, we propose a new slip estimation model using parametric representation. The model is not only interpretable, but is formulated based on terramechanics theory to minimize deviation from the actual environment. The impact of rover mass on slip characteristics is also considered, and incorporates this factor into the model. As a result, the model is expected to be used for applications where mass variations are anticipated, such as construction and resource transport rovers on the lunar surface. Moreover, by formulating the model in a parametric manner, the computational load is reduced compared to non-parametric machine learning-based methods. This advantage facilitates implementation on space-grade CPUs with limited computational resources. Additionally, we propose a new traversability assessment model that uses the developed slip model. This model incorporates uncertainties, such as model errors, into its parameters. We further demonstrate the motion planning that explicitly considers the slip risk using the proposed model.

The outline of this paper is as follows. In Section 2, we propose a parametric slip estimation model. Multi-body dynamics simulations were conducted to identify the key factors influencing slip, and the identified factors were then incorporated into the proposed model. Section 3 introduces the traversability assessment model incorporating the proposed parametric models. Section 4 presents the results of simulation verification that the proposed risk metrics can be used to move to the target point while minimizing the risk of wheel slip.

2 Slip estimation model

2.1 Slip metrics definition

First, the concept for evaluating slippage is explained. When a rover traverses on loose soil, the wheels of a rover can slip, resulting in a difference between the desired and actual velocity vectors (Wong, 1978). Slip in the longitudinal direction can be quantified by the slip ratio s, as shown in Equation 1:

$s=1-\frac{v_x}{v_{\text{ref}}},(1)$

where $v_x$ is the linear velocity in the longitudinal direction. The $v_{\text{ref}}$ is the desired velocity vector and is expressed as $v_{\text{ref}}=r\omega$ using the wheel radius $r$ and the angular velocity of the wheel $\omega$. When the wheel does not slip and the desired translational velocity is generated, $s=0$, and when the wheel is stuck and the velocity is zero, $s=1$. There are various definitions of lateral slip, but in this study, we use the slip angle as defined in Equation 2.

$\beta ={\tan }^{-1}\left(\frac{v_y}{v_x}\right),(2)$

where $v_y$ is the linear velocity in the lateral direction. The slip angle is an index that expresses the degree of lateral velocity relative to longitudinal velocity in terms of an angle.

2.2 Parametric slip model

Slippage is closely related to various environmental and rover system conditions. When estimating slip ratios and slip angles, one approach would be to consider as many of these factors as possible and build nonparametric models using techniques such as machine learning. However, while onboard sensors such as LiDAR can measure the terrain geometry, it can hardly directly identify soil properties such as cohesion and friction. Therefore, in this study, we extract parameters that are relatively easy to measure and have a significant impact on slip metrics, and approximate them using a parametric model for computational simplicity. One advantage of the parametric model is that it reduces computational load.

Investigating the relationship between various parameters and slip characteristics through experiments is not practical due to the enormous number of test cases required and accurately reproducing space conditions on Earth is difficult. Therefore, simulating low-gravity conditions and analyzing rover behavior is also a valuable approach for space applications. In fact, studies have analyzed the slip characteristics of NASA’s VIPER using a multi-body dynamics simulator (Hu et al., 2025). In this study, the open-source software Project Chrono was used to create datasets (Tasora et al., 2016; Mazhar et al., 2013). Considering the trade-off between computation time and model accuracy, the Soil Contact Model (SCM) was adopted (Krenn and Hirzinger, 2009). SCM represents the wheel as a polygonal mesh and performs contact detection with the ground for each grid, calculating the pressure acting on the wheel based on the Bekker-Wong theory (Bekker, 1969). This calculation is very simple, light on computational load, because assuming ideal contact conditions. SCM provides results that are reasonably close to actual phenomena, but the soil on actual planetary surfaces may consist of mixed soils with different properties and may include coarse materials such as gravel and rocks. As another approach, Discrete Element Method (DEM) simulations provide results that are more accurate than SCM, but it is still known to deviate from real-world behavior. In addition, DEM requires long computation times, making it unsuitable for applications such as generating large datasets, as in this study. Figure 1 shows a simulation of a 4-wheeled rover in the Project Chrono environment. The area where the wheels are in contact with the soil and pressure is generated is visualized in different colors.

Figure 1

Figure 1. Multi-body dynamics simulation with project chrono.

The parameter conditions for obtaining the dataset are shown in Table 1. The definitions of slope angle $\gamma$ and heading angle $\psi$ are illustrated in Figure 2. Slope angle, heading angle, total mass, and wheel angular velocity were selected as the parameters with particularly large effects on slippage, and data were collected for a total of 840 cases. The simulation was performed with parameter conditions within the range where the analysis was stable. For example, when the slope angle exceeds 25° or the translational velocity reaches 0.5 m/s, inertial and other dynamic effects become significant, causing the assumptions of the quasi-static motion model to break down and reducing the model’s reproducibility. However, such extreme parameter regions are not included in the applicable range of the slip risk metric introduced later in this study, because in practice, operation under these conditions would be excluded already at the motion planning stage. Therefore, these extreme conditions do not need to be explicitly considered for the purposes of this study. Constructing a dataset to evaluate the robustness of the proposed model is impractical because it would require handling a large number of parameters and an enormous number of analysis cases. Instead, it is more reasonable to assess robustness through experimental validation, which we consider as future work, while focusing here on the model’s nominal fitting ability. In addition, developing online parameter adaptation strategies for slip estimation, given the large parameter space, is a substantial challenge beyond the scope of this study and is also identified as future work. The slip ratio and slip angle are calculated based on measurements from each wheel and subsequently aggregated into representative values for the rover as a whole. The mass range was determined by assuming that the rover itself weighs 20 kg and that the payload can be up to 80 kg. In the case of a lunar rover application, the vehicle’s mass and payload are expected to be larger. However, the proposed slip estimation model described later is parametric, allowing it to be formulated in the same manner. Although the behavior should change depending on soil conditions, these parameters are not directly used to estimate the slip ratio, but rather other parameters are adjusted to simulate differences in characteristics depending on soil conditions.

Table 1

Table 1. SCM simulation conditions.

Figure 2

Figure 2. Definition of the slope angle $\gamma$, heading angle $\psi$ and reference velocity $v_{\text{ref}}$.

Figure 3 shows the relationship between the slip characteristics and the three parameters of heading angle, slope angle and linear velocity command. Here, the mass is fixed at 60 kg. Regardless of the other parameters, the slip ratio tends to increase as the slope angle increases. The slip ratio increases as the heading angle approaches zero, while the slip angle tends to decrease. The slip ratio and slip angle also tended to change more sensitively with the slope angle as the translational speed increased. In addition, while linear velocity had a large effect on the slippage when the velocity is small, the slippage do not change much when the velocity increases to a certain degree.

Figure 3

Figure 3. Multi-body dynamics simulation results: (a) slip ratio vs. linear velocity, (b) slip angle vs. linear velocity.

Figure 4 shows the results of the verification of how the slip characteristics change depending on the three parameters of heading angle, slope angle, and mass, with linear velocity fixed at 0.2 m/s. The same relationship between the slope angle, heading angle, and slip characteristics can be seen in the figures. As the mass decreases, the slip ratio becomes more sensitive to the slope angle. A similar trend was observed in the relationship between slip angle and mass.

Figure 4

Figure 4. Multi-body dynamics simulation results: (a) slip ratio vs. total mass, (b) slip angle vs. total mass.

Based on the trends of various parameters and slip characteristics obtained from the dataset, the results of formulating the slip ratio $s$ and slip angle $\beta$ as a parametric model are as follows.

$s\left(\gamma ,\psi ,\omega ,k_m\right)=a_1\tan \left(a_2\gamma \right)\cdot \cos \left(\psi \right)\cdot \tanh \left(a_3\omega \right)\cdot e^{-a_4{k}_m}+a_5,(3)$

$\beta \left(\gamma ,\psi ,\omega ,k_m\right)=b_1\tan \left(b_2\gamma \right)\cdot \sin \left(\psi \right)\cdot \tanh \left(b_3\omega \right)\cdot e^{-b_4{k}_m}+b_5,(4)$

where $a_1\dots a_5$ and $b_1\dots b_5$ are hyperparameters, which is tuned according to environmental conditions and rover configurations. $k_m$ is the ratio of the total mass to the reference mass $m_0$. The slip ratio and slip angle have the same equations except where the effect of the heading angle is expressed. The influence of each parameter on slip characteristics is formulated individually, and the overall model is expressed as a product of these terms. For example, when the slope angle $\gamma =0$, the rover is assumed to be traveling on flat terrain, and the first term on the right-hand side becomes zero. The proposed slip estimation model is constructed using a dataset generated by a multibody dynamics simulator that incorporates uncertainties such as rover system noises and environmental disturbances. As a result, the proposed model is formulated with explicit consideration of these uncertainties. In other words, rather than isolating and addressing each source of uncertainty individually, the proposed approach adopts a data-driven strategy to verify that the model functions robustly under the presence of various uncertainties. The following subsection qualitatively discusses the meaning of each term in the proposed slip estimation model.

2.3 Discussion on slip model

The parametric model of slip ratio and slip angle proposed in this study is discussed from a terramechanics perspective. Let ${\mu }_s$ be the coefficient of static friction, the conditions for the rover to slip against the ground are given by Equation 5

${\mu }_s>\tan \left(\gamma \right).(5)$

In other words, whether the rover slips or not is determined by the tangent of the slope angle. While the slip behavior on soft soil is not strictly deterministic, we simplify the model by assuming that the slip ratio and the slip angle, which indicate how prone the rover is to slipping, are proportional to $\tan \gamma$, as shown in the first factors of the first product terms in Equations 3, 4, respectively.

The direction of the drawbar pull that the rover is subjected to due to gravity is determined by the heading angle $\psi$. For example, when $\psi =0$, the rover is facing upslope, and the slip ratio is the highest at this time. Since the drawbar pull force acting on the rover is proportional to $\cos (\psi )$, the slip ratio is also modeled to be proportional to $\cos (\psi )$, as shown in the second factor of the first product term in Equation 3.

The normal and shear stresses beneath a wheel on the loose soil can be modeled as shown in Figure 5. In Equations 3, 4, the effect of angular velocity is formulated using hyperbolic tangents. If the angular velocity is sufficiently low, the rover behavior follows the terramechanics theory under static conditions. Equation 6 represents the shear stress $\tau$ (Janosi and Hanamoto, 1961).

${\tau }_i=\left(c+p\left(\theta \right)\tan \left(\phi \right)\right)\left[1-e^{-j_i/k_i}\right]\left(i=x,y\right),(6)$

where $c$ represents the cohesion stress, $\phi$ is the internal friction angle, $j_x$ and $j_y$ is the total soil deformation. $p(\theta )$ is the normal stress. $k_x$ and $k_y$ is the shear deformation modules. $j_x$ and $j_y$ are the total soil deformations and are expressed as follows (Wong and Reece, 1967; Yoshida and Ishigami, 2004):

$j_x\left(\theta \right)=r\left[{\theta }_f-\theta -\left(1-s\right)\left(\sin \left({\theta }_f\right)-\sin \left(\theta \right)\right)\right],(7)$

$j_y\left(\theta \right)=r\left(1-s\right)\left({\theta }_f-\theta \right)\cdot \tan \left(\beta \right),(8)$

where ${\theta }_f$ is the entry angle. As the angular velocity increases, the wheels move more soil, which can leads to an increase in $j$. Additionally, from Equations 8, 9, it is clear that as the slip ratio $s$ increases, $j_x$ also increases. Therefore, the simulation results, which show that the slip ratio $s$ increases with angular velocity, are consistent with this. Furthermore, $j_y$ increases as the slip angle $\beta$ increases, indicating that as angular velocity increases, the slip angle $\beta$ also increases. These characteristics are incorporated into the third factors in the first product terms of Equations 3, 4, respectively.

Figure 5

Figure 5. Wheel-soil contact model.

The fourth factors in the first product terms of Equations 3, 4 corresponds to the mass-dependent term. As the mass increases, the slip ratio becomes less sensitive to the slope angle, and a similar trend is observed for the slip angle. This is likely because the increased mass causes the wheels to sink into the ground, increasing the contact area between the wheel’s surface and the soil, which in turn increases the shear stress.

2.4 Hyperparameters optimization

In the previous section, the relationship between various parameters and slippage was investigated through multi-body dynamics simulations, and a new parametric models of slip ratio and slip angle were proposed. Here, we discuss the results of hyperparameter tuning to investigate how well the proposed model can simulate the simulation results. By solving the following optimization problem, the hyperparameters of the slip ratio model, $\mathit{a}≔[a_1,\dots ,a_5]$ are adjusted.

${\mathit{a}}^{*}=\underset{\mathit{a}}{\mathrm{arg\; min}}{\displaystyle \sum _{\gamma }}{\displaystyle \sum _{\psi }}{\displaystyle \sum _{\omega }}{\displaystyle \sum _{k_m}}{\left[s_{\text{SCM}}\left(\gamma ,\psi ,\omega ,k_m\right)-s\left(\gamma ,\psi ,\omega ,k_m\right)\right]}^2.(9)$

The hyperparameters of the slip angle model $\mathit{b}≔[b_1,\dots ,b_5]$ are obtained by solving the following equation.

${\mathit{b}}^{*}=\underset{\mathit{b}}{\mathrm{arg\; min}}{\displaystyle \sum _{\gamma }}{\displaystyle \sum _{\psi }}{\displaystyle \sum _{\omega }}{\displaystyle \sum _{k_m}}{\left[{\beta }_{\text{SCM}}\left(\gamma ,\psi ,\omega ,k_m\right)-\beta \left(\gamma ,\psi ,\omega ,k_m\right)\right]}^2.(10)$

The functions $s_{\text{SCM}}$ and ${\beta }_{\text{SCM}}$ correspond to values contained in the simulation dataset. The values of the hyperparameters obtained from Equations 9, 10 are shown in Table 2. For the optimization algorithm, we used Sequential Quadratic Programming (SQP). The error between the parametric model and the simulation data was evaluated using the Root Mean Squared Error (RMSE). As a result, the RMSE for the slip ratio was 0.0293, and the RMSE for the slip angle was 0.0279 radians. Figure 6 shows the results of comparing parametric slip models and simulations with fixed mass parameters. The primary purpose of calculating RMSE here is to provide reference values when setting the modeling errors in the traversability assessment model, rather than to evaluate the model’s goodness of fit. On the other hand, a previous study examining the relationship between slip ratio, slip angle, and wheel-generated drawbar pull reported that, despite estimation errors, the resulting drawbar pull remains small and does not exhibit large variations Ishigami et al. (2007). The dashed lines represent the simulation results, while the solid lines indicate the parametric model. It was confirmed that the proposed model is expressive enough to effectively reproduce the simulation results. Figure 7 is a graph comparing the slip model and simulation results with constant linear velocity. Similar to the previous results, it can be seen that the proposed model shows similar characteristics to the multi-body dynamics simulation.

Table 2

Table 2. Hyperparameter optimization results.

Figure 6

Figure 6. Model fitting results (solid lines: proposed parametric model; dotted lines: simulation data): (a) slip ratio vs. linear velocity, (b) slip angle vs. linear velocity.

Figure 7

Figure 7. Model fitting results (solid lines: proposed parametric model; dotted lines: simulation data): (a) slip ratio vs. total mass, (b) slip angle vs. total mass.

3 Traversability assessment model

In this section, we construct a traversability assessment model that serves as a guideline for rover motion planning by utilizing the proposed parametric model of slip characteristics. By referring to this assessment model during motion planning, the rover is expected to avoid the risk of the vehicle getting stuck and achieve safe traversal toward its target location. Furthermore, since the model adopts a parametric form, it also offers implementation advantages in terms of computational efficiency on onboard processors. Since slip ratio and slip angle represent different concepts, they are nondimensionalized to convert risk factor. When $s<0$, the actual translational velocity exceeds the commanded velocity, indicating that the rover is slipping in the direction of travel. Conversely, $s=1$ represents a state where the wheels are completely stuck, while $s>1$ indicates that the rover is moving backward. Since positive and negative values of $s$ correspond to different slip state, separate threshold values should be defined to determine hazardous slip conditions in each case. Thus, the risk factor $R_s$ for slip ratio is computed as follows:

$R_s=\left\{\begin{array}{cc}|s/s_{\text{th,up}}|\hfill & \mathrm{i}\mathrm{f}s\ge 0,\hfill \\ |s/s_{\text{th,low}}|\hfill & \mathrm{i}\mathrm{f}x<0,\hfill \end{array}\right.(11)$

where $s_{\text{th,up}}>0$ and $s_{\text{th,low}}>0$. The risk factor $R_{\beta }$ for slip angle is written by

$R_{\beta }=|\beta /{\beta }_{\text{th}}|,(12)$

where ${\beta }_{\text{th}}>0$. Figure 8 shows the normalization process of $R_s$ and $R_{\beta }$, as defined in Equations 11, 12. $s_{\text{th,low}}$, $s_{\text{th,up}}$ and ${\beta }_{\text{th}}$ are design parameters and are set to values so that each risk factor exceeds 1 when the slip ratio and slip angle exceed the mission tolerance.

Figure 8

Figure 8. Normalization process. (a) Risk factor $R_s$ and the slip ratio $s$, (b) the risk factor $R_{\beta }$ and the slip angle $\beta$.

To account for uncertainties, such as modeling errors, each risk factor is treated as a stochastic variable. The risk associated with the slip ratio is assumed to follow the distribution ${\mathcal{R}}_s\sim \mathcal{N}(R_s,{\sigma }_s^2)$. Similarly, the risk with respect to the slip angle is assumed to follow ${\mathcal{R}}_{\beta }\sim \mathcal{N}(R_{\beta },{\sigma }_{\beta }^2)$. The conversion from random variables to risk metrics uses Conditional Value-at-Risk (CVaR) (Majumdar and Pavone, 2020; Dixit et al., 2024). CVaR represents the expected value of a random variable when it exceeds a certain threshold. When $\mathcal{R}$ follows a normal distribution, CVaR can be obtained analytically as follows:

$\rho \left(\mathcal{R}\right)=\mu +\sigma \frac{\varphi \left({\mathrm{\Phi }}^{-1}\left(\alpha \right)\right)}{1-\alpha },(13)$

where $\varphi$ and $\mathrm{\Phi }$ is the probability density function and cumulative distribution function of a normal distribution, respectively. $\alpha \in \left.\left(\mathrm{0,1}\right.\right]$ is the risk level and the larger the value, the more conservative the risk metric. Based on Equation 13, $\rho ({\mathcal{R}}_s)$ and $\rho ({\mathcal{R}}_{\beta })$ can be calculated, and the risk is determined by whether each value exceeds 1.

4 Performance verification

4.1 Robot kinematic model

Simulation analysis were conducted to validate the effectiveness of the proposed traversability assessment model. In this study, we assume a differential drive rover as shown in Figure 9. The two-dimensional position and orientation of the rover’s body coordinate system relative to the inertial coordinate system is $\mathit{x}≔{[x,y,\theta ]}^T$, and the rover is assumed to move according to the following kinematics:

$\mathit{x}=\left[\begin{array}{c}\hfill \dot{x}\hfill \\ \hfill \dot{y}\hfill \\ \hfill \dot{\theta }\hfill \end{array}\right]=\left[\begin{array}{cc}\hfill r\cdot \cos \left(\theta \right)/2\hfill & \hfill r\cdot \cos \left(\theta \right)/2\hfill \\ \hfill r\cdot \sin \left(\theta \right)/2\hfill & \hfill r\cdot \sin \left(\theta \right)/2\hfill \\ \hfill -r/d\hfill & \hfill r/d\hfill \end{array}\right]\left[\begin{array}{c}\hfill {\omega }_l\hfill \\ \hfill {\omega }_r\hfill \end{array}\right],(14)$

where ${\omega }_l$ and ${\omega }_r$ are the angular velocity command values for the left and right wheels, respectively.

Figure 9

Figure 9. Coordinate frames and rover configuration parameters.

The differential two-wheel is a nonholonomic system and has some limitations, which means that the rover cannot move directly sideways and is constrained in the heading direction and the direction of travel. Therefore, in order to ensure that the rover can follow the generated path properly, a motion planning method has to guarantee that the kinematics is satisfied.

4.2 Motion planning algorithm

The proposed traversability assessment model should be set as a hard constraint for motion planning because it is a metric for avoiding vehicle getting stuck and any deviation from the specified value may lead to mission failure. As a result, a risk-constrained kinodynamic Rapidly-exploring Random Tree (RRT) algorithm was adopted for performance evaluation in this study. Algorithm 1 represents the algorithm of risk-constrained kinodynamic RRT. The initial state of the rover is ${\mathit{x}}_{\mathit{init}}$ and the target state is ${\mathit{x}}_{\mathit{goal}}$. In this case, these state vector consists of position and orientation $[x,y,\theta ]$ in a 2D plane. The $z$ can be uniquely obtained from terrain elevation data and 2D location. $K$ is the maximum number of iterations. $T$ is the tree. $C$ is a set of parameters used to calculate traversability, such as the mass $m_0$.

Algorithm 1

Algorithm 1. Risk-Constrained Kinodynamic RRT.

The algorithm is an improvement of the basic RRT (LaValle, 1998) to handle both kinematic constraints and traversability constraints (LaValle and Kuffner, 2001). The exploration starts with the tree initialized with the rover’s initial state ${\mathit{x}}_{\mathit{init}}$. In each iteration, a random sample of the 2D position and orientation ${\mathit{x}}_{\mathit{rand}}$ is selected from the exploration space (line 3). The function $\text{NearestVertex}(\mathit{x},T)$ searches for the vertex that is closest to the input $\mathit{x}$ (line 4). In this process, relative attitude is not considered; instead, the nearest neighbor is determined by computing the Euclidean distance based on relative position. ${\mathit{u}}_{\mathit{rand}}$ is a 2D state variable consisting of the rotational velocities of the rover’s wheels. A value is randomly sampled from the possible range of angular velocities (line 5). $\text{Propagation}(\mathit{x},\mathit{u})$ on line 6 is a function to compute the kinematics based on the state $\mathit{x}$ and control input $\mathit{u}$ using Equation 14. Substituting ${\mathit{x}}_{\mathit{near}}$ and ${\mathit{u}}_{\mathit{rand}}$ into the function, the rover’s state at the next time step is calculated and stored as ${\mathit{x}}_{\mathit{new}}$. Instead of directly sampling the state variables and adding them to the tree, the control inputs are sampled and propagated through time, ensuring that the generated random state variables are dynamically and kinematically valid. The risks based on the generated state ${\mathit{x}}_{\mathit{new}}$, control input ${\mathit{u}}_{\mathit{rand}}$, and parameter $C$ are calculated (lines 7–8). If the risks exceed a threshold, the state is rejected as it presents a high risk of slip. Otherwise, the new vertex and edge information is added to the tree $T$. Since the algorithm is sampling-based method, it easily handles constraints by simply checking if the new state satisfies the constraints. If ${\mathit{x}}_{\mathit{new}}$ is sufficiently close to the goal, the exploration is terminated; otherwise, it continues (lines 9–15).

In this study, the effectiveness of the traversability metric is validated using a risk-constrained kinodynamic RRT. The primary objective of the simulation is to evaluate the performance of the modeling method; therefore, a basic RRT was intentionally used as the motion planning algorithm. Improvements to the motion planning algorithm are considered beyond the scope for this work. Note that the proposed risk metric is not limited to this approach and can be applied to other motion planning techniques as well. Numerous constraint-aware motion planning methods exist, and they can be flexibly selected based on the specific application.

4.3 Simulation results

The terrain map used for validation is shown in Figure 10. The starting point is set at $[x,y]=[\mathrm{10,3}]$ and the goal point at $[x,y]=[\mathrm{10,17}]$. The rover is assumed to initially face the slope direction, $\psi =0$. The heading angle of the rover at the goal point is not specified. Both points are placed on a flat surface with no slope. The slope between the starting point and the goal point was set to an inclination of 15°. When torque limitations are present, slopes beyond a certain angle cannot be climbed, which can affect operational strategies. However, torque limitations are hardware constraints of the rover itself and are not directly effected to the slip estimation model or the traversability assessment model proposed in this study. Therefore, the effects of torque limitations are not considered in this simulation. In addition, performance comparisons on actual hardware are out of scope in this paper and the feasibility of implementing the algorithm on a space-grade CPU is not also evaluated.

Figure 10

Figure 10. Terrain for motion planning simulations.

Figure 11a shows a graph comparing how the generated paths change depending on the risk level. In this analysis, wheel radius $r=0.168$ m, the distance between wheels $d=0.5$ m, and mass $m=50$ kg were assumed. The hyperparameters for slip ratio and slip angle were obtained through optimization calculations, using the values from Table 2. The parameter settings for the traversability assessment model are shown in Table 3. The risk levels for slip ratio and slip angle, denoted by $\alpha$, are set to the same values, and in Figure 11a, they are set to (i) 0.5, (ii) 0.7, (iii) 0.9, and (iv) 0.95, respectively. By adopting the Kinodynamic RRT, smooth paths that satisfy the kinematic constraints were successfully generated. In all conditions, paths were generated that climbed diagonally in sloped areas. The risk metric for slip ratio decreases as the heading angle increases when climbing the slope, whereas the risk metric for slip angle decreases as the heading angle decreases, as shown in Equations 3, 4. Therefore, it can be considered that the path climbs the slope at a balanced angle that reduces both risk metrics.

Figure 11

Figure 11. Comparison of motion planning results under different risk levels $\alpha$: (a) generated paths, (b) time history of the linear velocity, and (c) time history of the risk metrics. The risk levels were set to (i) $\alpha =0.5$, (ii) $\alpha =0.7$, (iii) $\alpha =0.9$, and (iv) $\alpha =0.95$.

Table 3

Table 3. Parameter settings for traversability assessment model.

The comparison of the time history of the rover’s linear velocity is shown in Figure 11b. During periods when the rover is not climbing a slope, the linear velocity reaches nearly 0.3 m/s, but it decreases during slope climbing. Specifically, in Figure 11b, i the velocity during climbing remains below 0.2 m/s, while in Figure 11b, iv, it decreases to below 0.1 m/s. Thus, it was confirmed that the linear velocity during slope climbing tends to become lower as the value of $\alpha$ increases. Equations 3, 4 show that both the slip ratio and slip angle risk metrics decrease as the linear velocity decreases. These results demonstrate that the conservativeness of the motion planning can be adjusted by the risk level $\alpha$.

Figure 11c is the comparison of the time history of the risk metrics. From this figure, we can observe that the both risks of slip ratio and slip angle increased during the slope-climbing phase, and by reducing the linear velocity, the rover was able to adjust the risks to avoid exceeding the threshold of 1 while still moving. Even when the rover is located in an area where the slope angle is 0°, $\rho ({\mathcal{R}}_s)$ and $\rho ({\mathcal{R}}_{\beta })$ were not zero, indicating that a positive bias exists. This bias is caused by the second term in Equation 13, and it becomes larger as the risk level $\alpha$ increases. This explains why the conservativeness of the motion planning varies depending on the value of $\alpha$.

Next, we evaluated how the motion planning results varied with mass. Figure 12a shows the generated paths compared by different masses. The hyperparameters of the slip estimation model were set based on the values listed in Table 2. Parameters for the traversability assessment model were taken from Table 3, and the risk levels $\alpha$ for both slip ratio and slip angle were set to 0.7. Regarding the variances, they were set to include a margin, considering the modeling errors of the slip estimation model: the RMSE of slip ratio was 0.0293, and the RMSE of slip angle was 0.0279. In all cases, the generated path tended to climb the slope diagonally. In Figures 12a iii, iv there are sections where the path climbs more directly uphill.

Figure 12

Figure 12. Comparison of motion planning results under different total mass $m$: (a) generated paths, (b) time history of the linear velocity, and (c) time history of the risk metrics. The total masses were set to (i) $m=30$ kg, (ii) $m=50$ kg, (iii) $m=70$ kg, and (iv) $m=90$ kg.

This is because the increased mass results in a lower slip ratio, as described by Equation 3, which makes it less risky to move more directly in the slope direction.

The comparison of the time history of the rover’s linear velocity is shown in Figure 12b. As the mass increases, there is a tendency for the linear velocity while climbing a slope to become faster. Both the slip ratio and slip angle decrease with increasing mass, as shown in Equations 3, 4. Conversely, higher linear velocity tend to increase the slip ratio and slip angle. In other words, increasing the mass reduces the risk of slipping, which allows for higher linear velocity during slope climbing.

Figure 12c shows the time history of the risk metric for each case. The risk metrics remained below 1 in all cases. It was demonstrated that by using the RRT algorithm, the rover can safely reach the target position even under conditions of increased slip risk due to mass parameter settings, by appropriately adjusting its linear velocity and trajectory.

5 Conclusion

In this study, we first proposed parametric models of slip ratio and slip angle, which use four input parameters: slope angle, heading angle, mass, and wheel angular velocity. By comparing the results with multi-body dynamics simulations, it was confirmed that the proposed model has sufficient expressive power to estimate slip. We also proposed a traversability assessment model that incorporates the parametric models for slip ratio and slip angle. This assessment model is designed to allow tuning, taking into account uncertainties such as modeling errors. To evaluate the effectiveness of the proposed method, we conducted a simulation evaluation. As an example of a motion planning method, we introduced the risk-constrained kinodynamic RRT algorithm, which incorporates the proposed assessment model. Through the simulation results, we analyzed how the coefficients of the slip parametric model affect the generated paths. For future work, it would be useful to verify the effectiveness of the proposed slip estimation model through real-world tests. The test results would be an important step toward practical application. Additionally, we intend to model the estimation errors of the four input parameters when calculated onboard and evaluate their impact on the robustness of the motion planning approach.

Data availability statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

Author contributions

TN: Conceptualization, Data curation, Investigation, Methodology, Software, Validation, Visualization, Writing – original draft, Writing – review and editing. GI: Methodology, Supervision, Writing – review and editing.

Funding

The author(s) declare that no financial support was received for the research and/or publication of this article.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Bekker, M. G. (1969). Introduction to terrain-vehicle systems. Ann Arbor, Michigan: The University of Michigan Press.

Google Scholar

Chhaniyara, S., Brunskill, C., Yeomans, B., Matthews, M., Saaj, C., Ransom, S., et al. (2012). Terrain trafficability analysis and soil mechanical property identification for planetary rovers: a survey. J. Terramechanics 49, 115–128. doi:10.1016/j.jterra.2012.01.001

CrossRef Full Text | Google Scholar

Dixit, A., Fan, D. D., Otsu, K., Dey, S., Agha-Mohammadi, A.-A., and Burdick, J. W. (2025). STEP: stochastic traversability evaluation and planning for risk-aware navigation; results from the DARPA subterranean challenge. IEEE Transactions on Field Robotics 2, 81–99. doi:10.1109/TFR.2024.3512433

CrossRef Full Text | Google Scholar

Gonzalez, R., Apostolopoulos, D., and Iagnemma, K. (2018). Slippage and immobilization detection for planetary exploration rovers via machine learning and proprioceptive sensing. J. Field Robotics 35, 231–247. doi:10.1002/rob.21736

CrossRef Full Text | Google Scholar

Hu, W., Li, P., Rogg, A., Schepelmann, A., Chandler, S., Kamrin, K., et al. (2025). A study demonstrating that using gravitational offset to prepare extraterrestrial mobility missions is misleading. J. Field Robotics 42, 3772–3794. doi:10.1002/rob.22597

CrossRef Full Text | Google Scholar

Ishigami, G., Miwa, A., Nagatani, K., and Yoshida, K. (2007). Terramechanics-based model for steering maneuver of planetary exploration rovers on loose soil. J. Field robotics 24, 233–250. doi:10.1002/rob.20187

CrossRef Full Text | Google Scholar

Janosi, Z., and Hanamoto, B. (1961). “The analytical determination of drawbar pull as a function of slip for tracked vehicles in deformable soils,” in Int. Conf. on terrain-vehicle systems.

Google Scholar

Krenn, R., and Hirzinger, G. (2009). “SCM - a soil contact model for multi-body system simulations,” in 11th European regional conference of the international society for terrain-vehicle systems.

Google Scholar

LaValle, S. (1998). Rapidly-exploring random trees: a new tool for path planning. Res. Rep., 9811.

Google Scholar

LaValle, S. M., and Kuffner, J. J. (2001). Randomized kinodynamic planning. Int. J. robotics Res. 20, 378–400. doi:10.1177/02783640122067453

CrossRef Full Text | Google Scholar

Majumdar, A., and Pavone, M. (2020). “How should a robot assess risk? towards an axiomatic theory of risk in robotics,” in Robotics research (Springer), 75–84.

Google Scholar

Mazhar, H., Heyn, T., Pazouki, A., Melanz, D., Seidl, A., Bartholomew, A., et al. (2013). Chrono: a parallel multi-physics library for rigid-body, flexible-body, and fluid dynamics. Mech. Sci. 4, 49–64. doi:10.5194/ms-4-49-2013

CrossRef Full Text | Google Scholar

Sakayori, G., and Ishigami, G. (2021). Energy-aware trajectory planning for planetary rovers. Adv. Robot. 35, 1302–1316. doi:10.1080/01691864.2021.1959396

CrossRef Full Text | Google Scholar

Sanders, G. B., and Kleinhenz, J. E. (2024). “Update on NASA’s ISRU development and mission plans for the artemis program,” in 75th international astronautical congress (IAC).

Google Scholar

Sevastopoulos, C., and Konstantopoulos, S. (2022). A survey of traversability estimation for mobile robots. IEEE Access 10, 96331–96347. doi:10.1109/access.2022.3202545

CrossRef Full Text | Google Scholar

Smith, M., Craig, D., Herrmann, N., Mahoney, E., Krezel, J., McIntyre, N., et al. (2020). “The artemis program: an overview of NASA’s activities to return humans to the moon,” in IEEE aerospace conference.

Google Scholar

Sutoh, M., Yusa, J., Nagatani, K., and Yoshida, K. (2010). Traveling performance evaluation for planetary rovers on weak soil. J. Field Robotics, 546–551.

Google Scholar

Tasora, A., Serban, R., Mazhar, H., Pazouki, A., Melanz, D., Fleischmann, J., et al. (2016). “Chrono: an open source multi-physics dynamics engine,” in High performance computing in science and engineering (Springer), 19–49.

Google Scholar

Werkheiser, N., Ching, M., Galica, C., Sanchez, A. J., and Payne, S. (2024). “NASA’s lunar surface innovation initiative: ensuring a cohesive, executable strategy for technology development,” in AIAA aviation forum and ascend 2024.

Google Scholar

Wong, J. Y. (1978). Theory of ground vehicles. John Wiley and Sons.

Google Scholar

Wong, J. Y., and Reece, A. (1967). Prediction of rigid wheel performance based on the analysis of soil-wheel stresses part i. performance of driven rigid wheels. J. Terramechanics 4, 81–98. doi:10.1016/0022-4898(67)90105-x

CrossRef Full Text | Google Scholar

Yoshida, K., and Ishigami, G. (2004). Steering characteristics of a rigid wheel for exploration on loose soil in IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), Sendai, Japan, 28 September 2004 - 02 October 2004 (IEEE).

CrossRef Full Text | Google Scholar