Surrogate Motion Resistance Models for the Design of Omnidirectional Tracked Robots

1. Introduction

Mobile robots have become one of the key drivers behind the development of modern automation systems. Advances in sensing technologies, embedded computing, autonomous navigation, and intelligent control have significantly expanded their capabilities, leading to widespread adoption across numerous industrial and service applications [1,2,3,4,5,6,7,8].

Despite these developments, conventional mobile robots based on differential drive, Ackermann steering, or standard tracked locomotion remain subject to non-holonomic motion constraints. Such objects cannot generate lateral motion without changing their orientation, often requiring complex maneuvers and longer trajectories [9]. These limitations become particularly significant in applications involving narrow aisles, densely occupied workspaces, or frequent changes in travelling direction, where high maneuverability and precise positioning are essential [10]. To overcome these limitations, omnidirectional mobile robots have been developed. By employing specialized wheels or tracks, these objects can perform independent translational and rotational motion, allowing movement in any planar direction without prior reorientation. Such holonomic mobility enables more efficient operation not only in constrained environments [11], but also across challenging and uneven terrains [12]. Most omnidirectional mobile robots feature omni wheels or Mecanum wheels, whose kinematic and dynamic characteristics have been extensively investigated over the past decades [13,14,15]. However, wheeled omnidirectional robots exhibit several limitations when operating under increased loads or on surfaces requiring higher traction, such as greater sensitivity to surface irregularities, reduced load distribution, or motion disturbances [16].

An alternative solution is provided by omnidirectional tracked mobile robots, which combine the maneuverability of omnidirectional motion with the inherent advantages of tracked locomotion. An example of such a vehicle is shown in Figure 1.

Owing to their continuous contact with the supporting surface, tracked locomotion provides improved traction, enhanced load distribution, reduced ground pressure, and superior obstacle negotiation. These features render them highly effective for heavy-duty material handling, outdoor inspection, rescue operations, and other demanding applications that require both omnidirectional mobility and greater terrain adaptability[17]. Despite these advantages, the performance of an omnidirectional tracked robot is fundamentally governed by the complex interaction between its individual rollers and the supporting surface. The motion resistance generated at the roller-ground interface directly affects the driving torque requirements, energy consumption, traction characteristics, and overall dynamic behavior. Consequently, accurate determination of longitudinal and lateral motion resistance constitutes an important prerequisite for drivetrain sizing, energy-efficient motion planning, and optimization of the robot’s mechanical design. Since the global behaviour of the vehicle orginates from local contact occurring at the roller level, understanding these interactions is essential for developing reliable design and modeling.

Over the last decade, numerous studies have addressed design challenges in omnidirectional mobility. Table 1 summarizes representative literature in this field, considering locomotion type, level of analysis, used method, and primary research scope.

A critical analysis of the literature reveals three major research gaps. First, existing studies focus overwhelmingly on wheeled configurations, leaving omnidirectional tracked robots considerably less investigated. Second, the majority of research treats the robot as a complete macro-system, rather than a fundamental analysis of the mechanical behavior at the individual component level (i.e., the roller). Third, current modeling approaches rely primarily on computationally expensive physics-based models (e.g., full Multibody Dynamics), which are impractical for rapid iterative design or real-time control. In this context, surrogate models offer an alternative means of representing the response of computationally demanding simulations using simplified mathematical models. Their use has become increasingly common in robotics [33,34,35,36], particularly in applications that require repeated model evaluations, as they can substantially reduce computational time while maintaining satisfactory prediction accuracy.

Therefore, this study aims to develop computationally efficient surrogate models for determining the longitudinal and lateral motion resistance forces generated by an omnidirectional tracked roller as functions of its operating conditions and design parameters. First, a multibody dynamics model of a single roller is developed to reproduce the contact interaction between the roller and the supporting surface under different loading conditions, orientation angles, ground surface types, and bearing types. The numerical model is subsequently validated experimentally. The validated simulation model is then employed to generate a comprehensive dataset that describes the influence of selected design and operating parameters on motion resistance. Finally, multiple regression models are developed to provide a rapid determination of the resistance forces and to support the design process of omnidirectional tracked mobile robots. The main contributions of this study are the following:

• development and experimental validation of a multibody dynamics model of an omnidirectional tracked roller for motion resistance determination;
• quantitative identification of the influence of orientation angle, applied mass, surface type (through coefficient of friction consideration), and bearing type (through bearing friction torque consideration) on longitudinal and lateral motion resistance;
• formulation of computationally efficient surrogate models that enable rapid determination of motion resistance forces to directly support the design and real-time control of omnidirectional tracked mobile robots.

2. Materials and Methods

The adopted methodology combines physics-based modelling with data-driven surrogate modelling. A multibody dynamics model is first established to reproduce the roller-ground interaction and experimentally validated to ensure its predictive capability. The validated model is subsequently used to generate simulation data covering different combinations of design and operating parameters, which serve as the basis for constructing surrogate models of motion resistance. The overall workflow is summarized in Figure 2.

2.1. Multibody Dynamics Model of Motion Resistance

To investigate the motion resistance of the omni-track roller under various operating conditions, a three-dimensional multibody dynamics (MBD) model was developed in Hexagon Adams 2022. The model reproduces the rolling-sliding interaction between the roller and the ground surface, allowing the motion resistance forces (${\mathrm{F}}_{\mathrm{x}}$ and ${\mathrm{F}}_{\mathrm{y}}$) to be determined for different combinations of applied mass $m$, roller orientation angle $r_{\alpha }$(the angle between the roller’s rotational axis and the velocity vector), ground surface type (through coefficient of friction $\mu$ consideration), and bearing type (through bearing friction torque $T_b$ consideration).

The numerical model consists of three rigid bodies: a fixed ground surface, the investigated omni-track roller, and a roller carrier to which the roller is connected through a frictionless revolute joint (Figure 3). The revolute joint enables free rotation of the roller about its longitudinal axis while constraining the remaining relative motions.

A cylindrical roller geometry was considered in the simulations with a diameter of 40mm and a height of 68 mm. The geometry of the investigated roller is shown in Figure 4.

The model is defined in a global Cartesian coordinate system. The roller carrier is prescribed to move at a constant velocity along the global X-axis, representing the translational motion of the robot. Lateral translation along the Y-axis is constrained, whereas vertical translation along the Z-axis remains unconstrained, allowing the roller to adjust its vertical position in response to the applied mass and contact forces. Rotations about the X- and Y-axes are constrained, while the rotation about the Z-axis defines the roller orientation angle $r_{\alpha }$ with respect to the direction of motion.

The interaction between the roller and the ground surface is modelled using the Contact Force formulation available in the simulation environment. This contact model accounts for the normal contact stiffness and damping as well as friction, thereby reproducing the rolling and sliding phenomena at the roller-ground interface. The contact parameters adopted in the simulations are summarized in Table 2.

The coefficient of friction between the roller and the ground surface is treated as a variable simulation parameter to account for different surface conditions. Values of 0.60, 0.65, and 0.70 were selected to represent relatively high-friction contact conditions that may occur for elastomeric rollers operating on dry, relatively rough surfaces, such as concrete, asphalt, or textured industrial flooring.

The normal load acting on the roller is controlled by the applied mass, while the roller orientation angle varies to represent different design and operating conditions. In addition, the bearing friction torque is introduced through the revolute joint and varied according to the investigated bearing type. Four torque levels were considered to represent different levels of bearing resistance associated with the sealing configuration. A value of 0.2 Nmm was assigned to an open bearing, while 1 Nmm was used for a metal-shielded bearing (ZZ). For sealed bearings with contact rubber seals (2RS), two friction torque levels were considered: 3 Nmm for a light seal contact and 6 Nmm for a higher seal contact pressure. The simulation input parameters and their investigated ranges are summarized in Table 3.

For each combination of applied mass, roller orientation angle, coefficients of friction, and bearing friction torque, the model was simulated over a prescribed translational distance of 5mm at a constant velocity 10mm/s.

2.2. Test stand for the Experimental Validation of the Physical Model

Experimental tests were conducted to validate the developed multibody dynamics model and to assess its ability to reproduce the longitudinal and lateral motion resistance forces. The experiments were performed under controlled laboratory conditions. Three applied mass levels (0.385kg, 0.65kg, and 0.915kg) and five roller orientation angles, ranging from 15 to 75 in 15 increments, were investigated.

The experiments were conducted using a test stand (Figure 5) that enabled simultaneous measurement of the longitudinal force, acting along the direction of motion, and the lateral force, acting perpendicular to the direction of motion, generated during the translation of a passive roller over the ground surface.

The measuring system consists of a linearly moving aluminium frame driven by a lead screw and guided by linear guides. The frame provides a maximum measurement distance of 400 mm and a maximum translational velocity of 20 mm/s. The longitudinal and lateral forces are measured using strain-gauge force sensors. An MAVIN NS6 sensor with a measurement range of 0-6 N is used to measure the longitudinal force, while an MAVIN NA27 sensor with a measurement range of 0-20 N is used to measure the lateral force. The signals from both sensors are acquired using an ADT42-U measuring unit at a sampling frequency of 80 Hz.

A measuring carriage is connected to the longitudinal strain-gauge sensor by means of an inextensible cable. The carriage contains the roller under test and the applied mass and is designed to allow adjustment of the roller orientation angle. A schematic representation of the test stand and its main components is presented in Figure 6.

To ensure that the measured forces represented the steady-state motion resistance rather than transient effects associated with the initiation of motion, the data recorded during the initial acceleration phase were excluded from the analysis.

2.3. Surrogate Modelling

A multibody dynamics model provides an accurate determination of motion resistance but requires repeated numerical simulations for every combination of design and operating parameters. Such an approach becomes computationally demanding when large numbers of design variants must be evaluated. To enable rapid prediction of motion resistance while preserving satisfactory accuracy, surrogate modelling was employed.

In surrogate modelling, the results obtained from a computationally expensive model are used to construct a computationally efficient approximation, referred to as a surrogate model. The surrogate reproduces the relationship between the model inputs and the corresponding outputs while substantially reducing the computational effort required for subsequent predictions. This reduction in computational time and resources is one of the main advantages of surrogate modelling and makes it particularly useful for engineering problems involving complex and nonlinear systems [37,38].

The general procedure for constructing a surrogate model consists of three main stages [39]:

I. 

Data preparation and selection of the modelling approach

The first stage involves identifying the relevant input variables, defining their ranges, and generating a set of input-output observations that adequately represents the investigated design space. Based on the characteristics of the investigated problem, an appropriate surrogate modelling approach is then selected. Surrogate modelling methods can be broadly divided into polynomial-based and machine learning-based approaches [40].

In this study, the input variables correspond to the design and operating parameters summarized in Table 3. Since the motion resistance is characterized by two response variables, namely the longitudinal force and the lateral force, two independent surrogate models were developed.

Among the various modelling techniques, Polynomial Response Surface (PRS) models were selected due to their simple mathematical formulation, low computational cost, and straightforward interpretation of the influence of individual input variables on the response [41]. In addition, PRS models provide explicit analytical equations that can be directly implemented in design optimization and real-time control algorithms.

The PRS employs regression analysis to establish the relationship between input variables and responses. The general form of a PRS model can be expressed as follows [42]:

$$f(x)={\beta }_0+{\displaystyle \sum _{i=1}^m}{\beta }_i{x}_i+{\displaystyle \sum _{i=1}^m}{\displaystyle \sum _{j\ge i}^m}{\beta }_{ij}{x}_i{x}_j+\dots +\epsilon$$

(1)

where $\epsilon$ represents statistical error, $x_i$ denotes the i-th predictor variable, and β are the regression coefficients determined during the model fitting process.

The above formulation represents the surrogate as a linear combination of polynomial functions of the input variables, with linear and quadratic terms being most commonly used in the PRS approach.

II. 

Parameter estimation and model fitting

In the second stage, the parameters of the selected surrogate model are estimated using the generated data. Given the selection of a polynomial response surface as the modelling approach, the model parameters were estimated using the ordinary least squares method [38].

III. 

Model testing

Finally, the model is tested on data not used during fitting to assess its predictive capability and verify its accuracy. The model performance is evaluated by comparing its predictions with the corresponding results of the original computational model using the adjusted coefficient of determination ($R_{adj}^2$) and error metrics, including the mean absolute error (MAE), root mean square error (RMSE), and mean absolute percentage error (MAPE), defined as follows [42]

$$R^2={\left({\displaystyle \frac{{\displaystyle {\sum }_{i=1}^n}\left(x_i-\overline{x}\right)\left(y_i-\overline{y}\right)}{\sqrt{{\displaystyle {\sum }_{i=1}^n}{\left(x_i-\overline{x}\right)}^2{\displaystyle {\sum }_{i=1}^n}{\left(y_i-\overline{y}\right)}^2}}}\right)}^2$$

(2)

$$R_{adj}^2=1-(1-R^2)\left({\displaystyle \frac{n-1}{n-p-1}}\right)$$

(3)

$$RMSE=\sqrt{{\displaystyle \frac{{\displaystyle {\sum }_{i=1}^n}{\left(x_i-\overline{x}\right)}^2}{n}}}$$

(4)

$$MAE={\displaystyle \frac{{\displaystyle {\sum }_{i=1}^n}\left|x_i-y_i\right|}{n}}$$

(5)

where x-actual values, y-predicted values, p-number of predictors, and n-number of observations.

To assess the predictive capability of the surrogate models, the complete MBD dataset was divided into a model-development dataset and an independent validation dataset. Rather than applying a random split, the validation set was defined by a combination of operating conditions that was not used during model estimation. The model-development dataset comprised 3240 simulation cases corresponding to all 90 investigated orientation angles and six mass levels, combined with two friction coefficients (μ=0.6 and 0.65) and three bearing friction torque levels (0.2, 1, and 3 Nmm). The independent validation dataset comprised 540 cases covering the same range of orientation angles and mass levels, but with the friction coefficient and bearing friction torque simultaneously set to values excluded from model estimation, namely $\mu$ = 0.7 and $F_b$= 6 N. This resulted in a validation set representing an unseen combination of operating conditions rather than a random subset of the simulation data. This approach was adopted to evaluate whether the surrogate models could generalize to previously unseen combinations of friction and bearing friction torque while retaining coverage of the full ranges of the investigated roller orientation angles and applied masses.

3. Results

The obtained results are presented in three stages, following the workflow adopted in the methodology. First, the results of the multibody dynamics simulations are presented, focusing on how the considered design and operating parameters affect the motion resistance of the omni-track roller. Next, the numerical model is experimentally validated using laboratory measurements obtained under selected operating conditions. Finally, based on the simulation results, surrogate models are developed to represent the motion resistance as a function of the investigated parameters.

3.1. Multibody Dynamics Simulation Results

The multibody dynamics model was used to investigate the influence of roller orientation angle, applied mass, bearing friction torque, and friction coefficient on the motion resistance of the omni-track roller. The simulations were performed for roller orientation angles ranging from 0° to 89° in 1° increments, six applied mass levels, four bearing friction torque variants, and three coefficients of friction. In total, 6480 simulation cases were analysed, corresponding to all combinations of the investigated parameters. The complete simulation campaign required approximately five hours of computation.

Figure 7 presents the variation in longitudinal and lateral resistance forces as a function of the roller orientation angle for different applied mass levels, using a bearing friction torque of 1N and a friction coefficient of 0.6.

As the orientation angle increases, the longitudinal resistance force decreases continuously, whereas the lateral resistance force increases. The opposite trends of ${\mathrm{F}}_{\mathrm{x}}$ and ${\mathrm{F}}_{\mathrm{y}}$ indicate that the orientation angle primarily determines the distribution of the resistance force between the longitudinal and lateral directions. The dependence on the orientation angle is clearly non-linear. The decrease in ${\mathrm{F}}_{\mathrm{x}}$ and increase in ${\mathrm{F}}_{\mathrm{y}}$ are not proportional to the angle itself, indicating that a linear representation of $r_{\alpha }$ would not adequately describe the observed behaviour. The trends are consistent with the geometric decomposition of the resistance force into longitudinal and lateral components. In particular, the longitudinal and lateral components exhibit behaviour qualitatively corresponding to the $cos(r_{\alpha }$) and ${sin(r}_{\alpha }$) functions, respectively. These trends can be explained by the change in the direction of the contact reaction as the roller orientation angle increases. A larger orientation angle reduces the longitudinal component of the contact force while increasing its lateral component, resulting in a progressive redistribution of the resistance force between the ${\mathrm{F}}_{\mathrm{x}}$ and ${\mathrm{F}}_{\mathrm{y}}$ directions.

In addition to the angular dependence, both resistance force components increased as the applied mass increased. For a given roller orientation angle, higher applied masses resulted in higher values of both ${\mathrm{F}}_{\mathrm{x}}$ and ${\mathrm{F}}_{\mathrm{y}}$ (Figure 8). The only exception occurred at a roller orientation angle of 0°, where the lateral force remained equal to zero regardless of the applied load due to the absence of a lateral component of the contact force. The increase was approximately linear over the investigated range of applied masses. This behaviour was observed throughout the investigated angular range, although the magnitude of each force component depended strongly on the roller orientation angle. Thus, the applied mass primarily affected the magnitude of the resistance force, whereas the orientation angle determined its distribution between the longitudinal and lateral directions.

The observed linear dependence results from the adopted Coulomb friction model, in which the contact force is proportional to the normal reaction. Consequently, increasing the applied mass proportionally increases the longitudinal and lateral force components, while the roller orientation angle determines the distribution of the motion resistance force between its longitudinal ${\mathrm{F}}_{\mathrm{x}}$ and lateral ${\mathrm{F}}_{\mathrm{y}}$ components. At a roller orientation angle of 0°, the contact force acts exclusively in the longitudinal direction, resulting in a zero lateral force component. The combined effect of applied mass and roller orientation angle indicates that their effects cannot be considered completely independently. An increase in applied mass scales the magnitude of the resistance force, while the orientation angle determines the proportion of this force acting in the longitudinal and lateral directions.

The coefficient of friction also exhibited a systematic influence on the resistance forces. Increasing the friction coefficient resulted in an approximately proportional increase in the force components for a given applied mass and roller orientation angle. This effect is consistent with the frictional nature of the contact resistance, for which the magnitude of the friction force depends on both the applied mass and the coefficient of friction.

The bearing friction torque also affected the predicted motion resistance (Figure 9), although its contribution was substantially smaller than that associated with the applied mass and friction coefficient. For a given orientation angle, changes in bearing friction torque resulted in systematic changes in both force components. Importantly, the effect of bearing friction torque also depended on the roller orientation angle, indicating that its contribution to longitudinal and lateral resistance should be considered alongside the angular transformation.

The preliminary analysis of the MBD results, therefore, identified two distinct groups of effects. The applied mass, coefficient of friction, and bearing friction torque primarily influenced the magnitude of the resistance forces, whereas the roller orientation angle governed their distribution between the longitudinal and lateral directions. The angular dependence was non-linear and was qualitatively consistent with cosine and sine functions for ${\mathrm{F}}_{\mathrm{x}}$ and ${\mathrm{F}}_{\mathrm{y}}$, respectively. Moreover, the results indicated that the effects of applied mass and friction coefficient should be considered jointly, as their combined contribution determines the magnitude of the contact resistance. These observations were used to define the functional form of the subsequent surrogate models. Prior to this stage, the multibody dynamics model was experimentally validated to verify its accuracy.

3.2. Experimental Validation Results

To validate the developed multibody dynamics model, the simulation results were compared with the measurements obtained using the test stand described in Section 2.2. The numerical model was configured to reproduce the experimental conditions by considering five roller orientation angles (15°, 30°, 45°, 60°, and 75°) and three applied mass levels (0.385 kg, 0.65 kg, and 0.915 kg). The coefficient of friction between the roller and the ground surface was set to 0.6, based on the materials used in the experimental setup. Furthermore, the roller was modelled with bearings providing a constant friction torque of 1Nmm, in accordance with the manufacturer’s specifications. A quantitative assessment of the numerical model accuracy was conducted based on the MAE, RMSE, and MAPE metrics determined for the ${\mathrm{F}}_{\mathrm{x}}$ and ${\mathrm{F}}_{\mathrm{y}}$ forces (Table 4).

The obtained mean absolute percentage errors (MAPE) were 11% for the longitudinal resistance force ${\mathrm{F}}_{\mathrm{x}}$ and 8% for the lateral resistance force ${\mathrm{F}}_{\mathrm{y}}$, respectively. The mean absolute errors (MAE) did not exceed 0.24 N. The minor difference between the RMSE and MAE values confirms the absence of significant outliers and highlights the numerical stability of the developed model. The obtained metrics clearly demonstrate a high level of agreement between the model and experimental results, fully justifying its application to generate the dataset for surrogate modeling.

The scatter plots comparing the experimental and numerical model values for ${\mathrm{F}}_{\mathrm{x}}$ and ${\mathrm{F}}_{\mathrm{y}}$ demonstrate strong agreement across the analyzed operational range (Figure 10). The reference line represents the identity line (y = x), corresponding to ideal model accuracy.

For both force components, the data points cluster tightly around the identity line, particularly in the low-to-medium force range (0-2.5N), indicating high model accuracy. At higher load levels (> 3N), a slight systematic deviation is observed, with points lying marginally below the identity line. This reflects the model’s minor overestimation of forces, previously identified in the residual analysis. Overall, the minimal dispersion around the line of identity confirms a strong linear correlation between the model and physical measurements, validating the model for surrogate data generation.

An analysis of the residual plots (Figure 11) for the longitudinal (${\mathrm{F}}_{\mathrm{x}}$) and lateral (${\mathrm{F}}_{\mathrm{y}}$) motion resistance forces reveals two key characteristics of the developed model. First, the positive sign of most residuals indicates a slight, systematic overestimation of force values by the model compared to experimental testing. Second, the residuals exhibit a proportional behavior. Their absolute value increases with load, while the relative error remains constant (MAPE at 8-11%). This phenomenon stems from the idealization of contact geometry and the constant friction coefficient assumed in the simulation environment, whereas real-world conditions involve local micro-slip and nonlinear elastic deformations of the roller material.

The results presented above confirm a high degree of agreement between the numerical model and the experimental data. Therefore, the validated model serves as a reliable tool for generating datasets to construct surrogate models.

3.3. Surrogate Motion Resistence Models

The surrogate modelling procedure comprised models development based on a dedicated subset of the MBD simulation data and subsequent model validation using independent MBD observations.

3.3.1. Models Development

The surrogate models were developed using a dedicated subset of 3240 observations from the MBD simulation dataset, separated from the data used for subsequent model validation. Based on the relationships identified in the MBD simulations, two surrogate models were developed to predict the longitudinal and lateral motion resistance forces, (${\mathrm{F}}_{\mathrm{x}}$) and (${\mathrm{F}}_{\mathrm{y}}$), respectively. Multiple linear regression was employed to construct the surrogate models.

The preliminary analysis presented in Section 3.1 indicated that the relationship between the roller orientation angle and the resistance forces was nonlinear. In particular, the variation of ${\mathrm{F}}_{\mathrm{x}}$ and ${\mathrm{F}}_{\mathrm{y}}$ with the roller orientation angle followed trends consistent with cosine- and sine-type dependencies, respectively. The analysis also indicated that the effect of the orientation angle should not be considered independently of the other model inputs, as the influence of the applied mass and friction coefficient changed with the roller orientation. Therefore, both nonlinear angular effects and interactions between the input variables were considered during model development. To account for the nonlinear angular dependence while retaining a linear regression framework, the roller orientation angle was transformed using the cosine function for the longitudinal force model and the sine function for the lateral force model. This transformation was motivated by the observed force decomposition: ${\mathrm{F}}_{\mathrm{x}}$ decreased towards zero as the roller orientation angle approached 90°, whereas ${\mathrm{F}}_{\mathrm{y}}$ increased towards its maximum. Accordingly, $cos(r_{\alpha }$) was used to represent the angular dependence of ${\mathrm{F}}_{\mathrm{x}}$, while $sin(r_{\alpha }$) was used for ${\mathrm{F}}_{\mathrm{y}}$.

In addition to the trigonometric transformation, interaction terms were introduced to account for the combined influence of the applied mass, friction coefficient, and roller orientation angle. For the longitudinal force, the transformed variables considered during model development included:

$$X_1^{F_x}=cos\left(r_{\alpha }\right)$$

(6)

$$X_2^{F_x}=m\mu cos\left(r_{\alpha }\right)$$

(7)

$$X_3^{F_x}=T_{b}cos\left(r_{\alpha }\right)$$

(8)

For the lateral force, the corresponding variables were defined as:

$$X_1^{F_y}=sin\left(r_{\alpha }\right)$$

(9)

$$X_2^{F_y}=m\mu sin\left(r_{\alpha }\right)$$

(10)

$$X_3^{F_y}=T_{b}sin\left(r_{\alpha }\right)$$

(11)

The final model was obtained using the ordinary least-squares method combined with backward variable selection. Equation (12) presents the surrogate model developed for predicting the longitudinal motion resistance ${\mathrm{F}}_{\mathrm{x}}$:

$$F_x=0.135-0.172cos\left(r_{\alpha }\right)+9.82m\mu cos\left(r_{\alpha }\right)+0.024T_{b}cos\left(r_{\alpha }\right)$$

(12)

The model consists of an intercept and three physically interpretable terms. The first angular term, $cos\left(r_{\alpha }\right)$, represents the direct contribution associated with the orientation of the roller. The $m\mu cos\left(r_{\alpha }\right)$term describes the combined effect of the applied mass and friction coefficient on the longitudinal resistance, while the $F_{b}cos\left(r_{\alpha }\right)$ term represents the contribution of bearing friction torque projected onto the longitudinal direction. The positive coefficient of the $m\mu cos\left(r_{\alpha }\right)$ term indicates that increasing either the applied mass or the friction coefficient increases the longitudinal resistance. The positive coefficient associated with bearing friction torque indicates an additional increase in ${\mathrm{F}}_{\mathrm{x}}$ with increasing bearing friction torque.

Equation (13) presents the surrogate model developed for predicting the lateral motion resistance $F_y$:

$$F_y=-0.66+0.106sin\left(r_{\alpha }\right)+9.69m\mu sin\left(r_{\alpha }\right)-0.026T_{b}sin\left(r_{\alpha }\right)$$

(13)

The structure of the lateral force model is analogous to that of the longitudinal force model, but the cosine transformation is replaced by the sine transformation to reflect the opposite angular dependence of the two force components. The $m\mu sin\left(r_{\alpha }\right)$ term represents the combined effect of the applied mass and friction coefficient on the lateral resistance, whereas $F_{b}sin\left(r_{\alpha }\right)$ accounts for the contribution of bearing friction torque to the lateral force.

Regression coefficients in both models were statistically significant at (p<0.001) when evaluated using heteroskedasticity-consistent standard errors. The models also showed a very high goodness of fit to the development data. For the $F_x$ model, $R_{adj}^2$=0.999871, with a residual standard error of 0.031. For the $F_y$ model, the corresponding values were $R_{adj}^2$=0.999785, and a residual standard error of 0.041. The overall regression models were statistically significant (p<0.001).

Multicollinearity was assessed using variance inflation factors (VIF). The VIF values were low for both models, ranging from 1.261 to 2.238 for $F_x$ and from 1.275 to 2.307 for $F_y$, indicating no substantial multicollinearity among the retained predictors. The Breusch–Pagan test indicated statistically significant heteroskedasticity for both models. This indicates that the residual variance was not constant across the investigated operating conditions. Therefore, heteroskedasticity-consistent HC1 standard errors were used for statistical inference.

The models were subsequently evaluated using the independent MBD observations excluded from model development. The results of the validation procedures are presented in Section 3.3.2.

3.3.2. Models Validation

The predictive performance of the developed surrogate models was evaluated using the remaining 540 MBD observations that were not included in the model development process. These observations represented combinations of operating conditions excluded from model fitting, including a different bearing friction torque and coefficient of friction from those used for model development. The previously estimated regression equations were applied directly to these observations without any further parameter estimation or model adjustment.

The validation results are summarized in Table 5.

Both surrogate models showed a close agreement with the independent MBD results, with $R_{adj}^2$ values exceeding 0.99. The model for ${\mathrm{F}}_{\mathrm{x}}$ achieved a higher predictive accuracy, with an MAE of 0.087 and an RMSE of 0.101, compared with 0.157 and 0.186, respectively, for ${\mathrm{F}}_{\mathrm{y}}$. The maximum absolute residuals were 0.250 for ${\mathrm{F}}_{\mathrm{x}}$ and 0.363 for ${\mathrm{F}}_{\mathrm{y}}$, indicating that the lateral force was more difficult to approximate than the longitudinal force. Nevertheless, the obtained $R_{adj}^2$ values of 0.9989 for ${\mathrm{F}}_{\mathrm{x}}$ and 0.9965 for ${\mathrm{F}}_{\mathrm{y}}$ indicate that the surrogate models reproduced the variation of the corresponding MBD force components with high accuracy

The validation results demonstrate a close agreement between the surrogate predictions and the independent MBD results for both force components. The slightly higher prediction errors observed for ${\mathrm{F}}_{\mathrm{y}}$ indicate that the lateral force was more difficult to approximate than the longitudinal force, although the coefficient of determination remained above 0.99. The comparison between predicted and reference values is presented in Figure 12 as scatter plots, where the proximity of observations to the 1:1 line provides a visual assessment of prediction accuracy.

The close clustering of the observations around the 1:1 line confirms the high agreement between the surrogate predictions and the independent MBD results.

4. Discussion

The results indicate that the proposed surrogate models provide a sufficiently accurate representation of the motion resistance predicted by the MBD model while substantially reducing the computational effort of repeated evaluations. The close agreement between the surrogate predictions and both independent MBD data and experimental measurements suggests that the models can serve as a practical approximation of the motion resistance forces within the investigated operating range.

The developed models have direct potential for application in the design of mobile robotic platforms equipped with omnidirectional tracks. In particular, the predicted longitudinal resistance can be incorporated into the estimation of the required drive torque and motor power. The model can therefore support motor and drivetrain sizing by enabling evaluation of the resistance force under different roller orientations, loads, friction conditions, and bearing friction torque levels without repeatedly performing computationally expensive MBD simulations. The model can also support the selection and comparison of bearing solutions, as bearing friction torque constitutes an explicit input to the surrogate formulation. This may be particularly useful at the preliminary design stage, where many drivetrain configurations must be evaluated.

From a system design perspective, the model also enables consideration of the effects of individual design choices jointly. For example, reducing bearing friction torque may decrease the total motion resistance and consequently reduce the required drive torque, whereas changes in contact conditions or vehicle loading may have a substantially larger effect on the required driving force. The surrogate formulation therefore provides a computationally inexpensive means of assessing such design trade-offs before detailed mechanical design or full system simulation.

The low computational complexity of analytical surrogate models also enables their use in online control. Unlike an MBD model, the proposed equations require only a small number of arithmetic operations. This makes them suitable for implementation in embedded control systems, where computational resources and response time are limited. For example, the predicted longitudinal and lateral resistance could be incorporated into motion controllers to compensate for changes in roller orientation, payload, or estimated contact conditions. Such a model could potentially be used for feed forward compensation, torque allocation, or real-time estimation of the force demand associated with the omnidirectional drive. However, implementation in an online control system would require additional validation under dynamically varying operating conditions, including changes in velocity, acceleration, payload, and contact conditions.

Several limitations of the present study should be acknowledged. First, the models were developed for a single roller geometry and within a predefined range of operating conditions. The angular domain was limited to 0°≤$r_{\alpha }$≤89°. Although the observed dependence of ${\mathrm{F}}_{\mathrm{x}}$ on the orientation angle is well represented by the cosine transformation within this range, the behaviour at exactly 90° is discontinuous, with the longitudinal resistance force dropping to zero. Consequently, the trigonometric formulation adopted in the present models should not be extrapolated to 90°, and the models should be regarded as valid only within the investigated angular range. Second, the investigation focused on a defined range of operating parameters. The experimental validation of the MBD model was performed for three mass levels, five roller orientation angles, one bearing friction torque level, and one coefficient of friction. Therefore, although the MBD model was validated against experimental measurements, the experimental validation did not cover the entire parameter space used for surrogate model development. The additional validation of the surrogate model was performed using independent MBD simulations and therefore primarily assessed its ability to reproduce the numerical model outside the training dataset. Finally, the present models do not explicitly account for parameters such as traveling velocity, acceleration, dynamic load transfer, surface irregularities, roller deformation, or changes in contact conditions during operation. These factors may become relevant when the model is transferred from controlled simulation conditions to a complete mobile robotic platform operating in realistic environments.

Future research should extend the proposed approach from an individual omni-track roller to the complete mobile robotic system. In a full-robot model, the interactions between multiple rollers, track geometry, load distribution, and the drive system could be considered simultaneously. Such an extension would allow the surrogate approach to be used not only for predicting the resistance of an individual roller but also for estimating the overall motion resistance and energy demand of the robotic platform. Further research could also investigate the influence of other drivetrain components and design parameters, including roller geometry, bearing characteristics, track material, contact stiffness, and surface properties. Extending the experimental campaign to a broader range of operating conditions would additionally enable more comprehensive validation of the surrogate models and assessment of their applicability under realistic operating scenarios. Ultimately, integrating such surrogate models with robot-level simulation and control frameworks could provide a computationally efficient link between mechanical design, energy aware motion planning, and real-time control of omnidirectional tracked robots.

5. Conclusions

This study presented a surrogate modelling approach for estimating the motion resistance forces of an omnidirectional tracked roller under varying operating conditions and design parameters. A multibody dynamics model was developed in Hexagon Adams to investigate the effects of roller orientation angle, applied mass, coefficient of friction, and bearing friction torque on the longitudinal and lateral force components. Based on the MBD simulation data, Polynomial Response Surface surrogate models were developed to provide a computationally efficient representation of the resistance forces. The surrogate models achieved a close agreement with the reference MBD results, demonstrating that the main relationships governing the force response can be captured without performing a full multibody simulation for each new combination of input parameters. Independent comparison with experimental measurements further supported the applicability of the modelling approach within the investigated operating conditions.

The proposed surrogate models provide a practical tool for rapid estimation of motion resistance during the preliminary design and analysis of omnidirectional tracked systems. Their use can facilitate repeated evaluations of different operating conditions and support the selection of drive and mechanical components without the computational cost associated with repeated MBD simulations. The approach can also serve as a basis for integrating resistance estimation into higher level robot design, optimization, and control.

Future research should extend the experimental validation to a broader range of operating conditions and investigate the applicability of the proposed approach to complete omnidirectional tracked robot systems.

The main contribution of the study is not only the development of the surrogate models, but also the demonstration that a computationally expensive MBD-based motion resistance analysis can be transformed into a fast predictive tool suitable for engineering applications.