Tracking Control of a Four-Wheeled Skid-Steered Robot with Slip Compensation and Application of the Drive Unit Model

Abstract

This article focuses on trajectory tracking control of a four-wheeled mobile robot, with non-steered wheels. The issues in terms of robot kinematics are discussed and a dynamics model is derived, which additionally took into account the drive unit model. This paper analyses four versions of the control system, which take into account the possibility of compensating for wheel slip and non-linearities resulting from the drive unit model. It is assumed that the wheel-slip compensation is based on the measurement of the actual robot’s motion parameters. The linear and angular motion parameters of the robot’s mobile platform are taken into account, which allows for the estimation of the wheel slip velocities. The results of the simulation studies are presented, consisting of the evaluation of individual control system solutions in terms of achieving the highest possible accuracy in executing a prescribed trajectory. Additionally, the impact of the investigated control strategies on electric power demand and electric energy consumption by the robot’s drives is analyzed. In order to quantitatively assess the control system solutions, quality indexes were adopted, focusing on tracking accuracy and energy efficiency. The research results indicate that incorporating wheel-slip compensation into the control system enables high accuracy to be achieved in terms of trajectory tracking. In turn, the use of the drive unit model within the control system leads to an increase in the accuracy of the robot’s wheel movements, which does not ultimately result in an increase in the accuracy of the motion of the robot’s mobile platform due to the slipping of the wheels. It was also observed that improving the trajectory tracking accuracy leads to an increase in the maximum electric power demand and electric energy consumption by the robot’s drives.

1. Introduction

Wheeled mobile robots are widely used in various fields, both in open terrain and enclosed facilities. These applications include defense and security, industrial transport and logistics, inspection and services, agriculture, extraterrestrial exploration, education and entertainment, as well as medical assistance.

The use of wheeled robots in such a wide range of areas requires a broad spectrum of solutions, particularly adapted to the type of ground surface, working environment, and the specific nature of the tasks performed. Consequently, various kinematic structures for wheeled robots have been developed. In regard to indoor applications, these kinematic structures primarily include differentially driven and omnidirectional robots, while in regard to outdoor applications, they often involve car-like and rover-type robots.

There are also examples of wheeled robot kinematic structures that perform well in diverse environments, a good example of which are the rovers that are the subject of the work by Iagnemma and Dubowsky, 2004 [1]. These structures also include skid-steered robots, equipped exclusively with non-steered wheels, as well as various designs of robots involving hybrid locomotion systems. Selected solutions in terms of hybrid robots, along with an example of the modeling process, are discussed in a paper by Trojnacki, 2010 [2], while an article by Ni et al., 2020 [3], focuses on the problem of motion control in terms of a four-wheel-legged robot.

Skid-steered robots are characterized by their simple design, good mobility, good maneuverability, and good motion stability. However, the use of non-steered wheels results in wheel slip during in-place turning and rotating, leading to poor dead-reckoning accuracy, and also causes faster tire wear, as noted in a paper by Khan et al., 2020 [4].

Due to the wheel slip feature in regard to skid-steered robots, which is an integral feature of their motion, they present a relatively challenging subject for dynamics modeling and trajectory tracking control.

When it comes to tracking control of skid-steered robots, as follows from, for e.g., the work by Zhang et al., 2020 [5] and Khan et al., 2020 [4], the key challenge is to ensure high motion accuracy under varying operating conditions, while maintaining high energy efficiency. This raises the question of whether it is possible to compensate for wheel slip in a simple and effective manner to achieve high motion accuracy, and how such compensation affects energy efficiency.

Another important aspect of robot motion control is the use of a robot model within the control system. For this purpose, a robot dynamics model is typically used, which allows the control system to compensate for the non-linearity of the controlled object. Much less common is the use of the drive unit model of the robot itself, which can also potentially improve tracking control accuracy. In addition, such a solution may be characterized by greater simplicity in terms of the control system in relation to a case involving the use of a robot dynamics model. The simplicity and effectiveness of the solution are associated with the ease of its implementation, which potentially enables its wide application.

Considering the aspects discussed, it was assumed that the aim of this study was therefore to investigate the influence of wheel-slip compensation and the application of the drive unit model on the trajectory tracking accuracy of a skid-steered robot, as well as the impact of such a solution on the energy efficiency of the robot. This study focuses on a lightweight four-wheeled skid-steered robot design, due to the high popularity of this kind of configuration.

To achieve the aim of the article, it is necessary to answer the following key research questions:

• Does wheel-slip compensation significantly improve tracking control accuracy?
• Can incorporating the drive unit model into the control system significantly enhance the accuracy of tracking control?
• What is the correlation between the accuracy of tracking control and energy efficiency?

To answer these research questions, it is necessary to consider the robot model, taking into account its kinematics, dynamics, and the characteristics of the drive units, as well as to develop various versions of the control system.

The research methodology adopted in this paper includes the following:

• Proposes the desired trajectory, which for the analyzed robot does not lead to excessive wheel slip and is executable, considering the limitations of the drive units on the electric power, maximum current, and maximum rotational speed;
• Develops and implements the robot model, as well as individual versions of the control system in a simulation environment, taking into account wheel-slip compensation, based on the measurement of the actual motion parameters of the robot’s mobile platform, and the inclusion of drive unit models;
• Introduces quality indexes to evaluate the individual control system solutions;
• Conducts simulation studies for selected motion trajectories and the analyzed control system solutions;
• Analyzes and discusses the simulation results, considering the obtained quality indexes, motion paths, and time histories of the selected physical quantities;
• Formulates conclusions based on the conducted studies, provides responses to the research questions, and indicates potential directions for future research work.

This article focuses on simulation studies involving various versions of the motion control system for the PIAP GRANITE robot, with non-steered wheels. The chosen robot design is well-known to the author from previous research, enabling an accurate representation of the robot’s dynamic properties and the inclusion of parameters that are close to those of the real robot.

The main contributions of this article include the following:

• The derivation of a comprehensive robot model is carried out, considering kinematics, dynamics, and Dugoff’s tire model, taking into account wheel slip phenomena, as well as the properties of drive units;
• The formulation of various versions of wheel controllers is achieved, enabling wheel-slip compensation or incorporating drive unit models, and also a discussion on the stability analysis of the proposed solutions for the analyzed robot is presented. A definition of quality indexes is provided, enabling a comprehensive assessment of individual control system solutions, taking into account the issue of motion accuracy and energy efficiency;
• A discussion is presented on how the wheel slip velocity can be estimated, based on the measurement of the linear velocity of a selected point of the robot and the angular velocity of its turning;
• Simulation studies are carried out using the adopted model of the robot and various versions of the control system, providing answers to the key research questions regarding the significance of wheel-slip compensation and the inclusion of the drive unit model within the control system structure, as well as the correlation between motion accuracy and energy efficiency in regard to the analyzed cases.

The contents of the subsequent sections in the article are the following:

• Section 2 focuses on a state-of-the-art review, on the modeling and trajectory tracking control of wheeled mobile robots, with particular attention paid to the distinctive features of skid-steered robots and the wheel slip phenomenon;
• Section 3 introduces the robot model, which encompasses the kinematics and dynamics of the robot, including the so-called tire model, which takes into account the wheel slip phenomenon, as well as the drive unit model, enabling an accurate representation of the real robot’s operation;
• Section 4 defines the various versions of the control system, enabling wheel-slip compensation, based on the measurement of the actual motion parameters of the robot’s mobile platform, or incorporating the drive unit model, as well as generally discussing the issue in terms of control system stability analysis;
• Section 5 focuses on simulation studies involving the robot control system, including presenting the assumed robot’s and environment-related parameters, the desired motion parameters and trajectory, the settings for the individual control system solutions, and the quality indexes used for the quantitative evaluation of these solutions, as well as the simulation results presented in terms of the quality indexes, motion paths, and time histories of the selected physical quantities;
• Section 6 provides the conclusion of the article, formulates answers to the research questions, and indicates potential directions for future research.

2. State-of-the-Art Review

The modeling of wheeled robots requires a description of their kinematics and dynamics. These issues are discussed, for instance, in the article by Alexa et al., 2023 [6], which focuses on wheeled and tracked robots, as well as in the publication by Moshayedi et al., 2022 [7], which examines various configurations of service robots.

Modeling the dynamics of wheeled robots typically employs Lagrange equations, as reflected in works such as the articles by Chakraborty and Ghosal, 2005 [8], Hendzel and Rykała, 2017 [9], and Locardi, 2020 [10]. Another common approach is the use of Newton–Euler equations, an example of which is the article by Khadr, 2024 [11]. Less frequently, alternative methods such as Kane’s approach are applied, as in the study by Thanjavur and Rajagopalan, 1997 [12].

The approach used for modeling wheeled mobile robots depends on the kinematic structure of the robot. For some structures, their dynamics can be modeled with sufficient accuracy, neglecting the wheel slip phenomenon.

Examples of modeling and simulating the motion of wheeled mobile robots, while neglecting wheel slip, are presented in the study by Chakraborty and Ghosal, 2005 [8], which examines a three-wheeled robot with torus-shaped wheels, and in the publication by Hendzel and Rykała, 2017 [9], which focuses on an omnidirectional robot, equipped with Mecanum wheels. In such cases, it is possible to achieve relatively accurate control of the robot based on solving the inverse kinematics problem, which involves determining the angular velocities of the wheels for the desired generalized velocities of the mobile platform.

It should be noted that all wheeled robots experience some degree of wheel slip, for instance, due to the torsional compliance of the tires. The wheel slip aspect is considered, for e.g., in the work by Khadr, 2024 [11], for a three-wheeled robot with a differentially driven kinematic structure. However, there are some mobile robots according to which wheel slip is an inherent feature, which has a significant impact on their motion. This phenomenon must, therefore, be incorporated into the dynamics model. An example of such robots is those with only non-steered wheels, which move during slipping conditions, when turning or rotating in place; therefore, they are called skid-steered robots. The challenges associated with these robots are discussed in the study by Khan et al., 2020 [4].

The wheel slip phenomena involving wheeled vehicles are extensively studied in the literature, especially concerning automotive vehicles. A key reference in this area is the book by Pacejka, 2005 [13], which covers, among other topics, the modeling of tire characteristics, based on empirical tests. For off-road vehicles, the foundational work is the book by Bekker, 1960 [14], whose author initiated a new field of knowledge, known as terramechanics. On the other hand, the book by Wong, 2001 [15], discusses modeling issues for both types of wheeled vehicles, as well as for tracked vehicles.

However, knowledge regarding the modeling of motor vehicles cannot be directly transferred to the modeling of lightweight wheeled mobile robots. This is due to, among others, the following reasons: this type of mobile robot has a specific design in terms of the wheels, which are often non-pneumatic; the forces resulting from the wheel–ground interaction are of a different order; the drive units are of a different type; such vehicles can move on more diversified ground surfaces (and each wheel can even move on different ground surfaces and in different conditions); and for skid-steered robots, the range of lateral slip is much wider than that in the typical operation of motor vehicles. In the case of motor vehicles, a 2 DoF single-track vehicle model (also called a bicycle model) is often used, as in the work by Tufano et al., 2024 [16]. In this type of model, the two left and right wheels are represented by one single wheel. Unfortunately, this approach cannot be used for lightweight wheeled mobile robots with non-steered wheels, which require an individual approach to modeling. The work by Dąbek and Trojnacki, 2016 [17], describes, in more detail, the differences in the structure of the wheels in typical mobile robots in relation to motor vehicles and the requirements for tire modeling of lightweight wheeled robots. This paper tries to fill the gap in regard to the small amount of research that comprehensively considers the characteristics of skid-steered robots, including adequate modeling of the wheel slip phenomenon.

In the case of wheeled mobile robots, their motion depends not only on how the individual wheels are controlled, which determines the kinematics of the locomotion system, but also on the forces occurring in the area of contact between the wheels and the ground. When it comes to modeling the motion of skid-steered robots, it is especially important to know the tire–ground interaction forces at play. In particular, the lateral forces that arise during the turning of such robots can be comparable in magnitude to longitudinal or normal forces, significantly affecting the robot’s motion. These forces are also crucial during the design phase of the robot, as this type of robot requires particularly robust mounting of the wheels to the axles. Modeling the dynamics of such vehicles requires the implementation of a so-called tire model, which accurately describes the wheel–ground interaction, including the wheel slip phenomenon. Examples of dynamics modeling of skid-steered robots that consider the wheel slip phenomenon can be found, for e.g., in the work by Zhang et al., 2021 [18], for a six-wheeled robot, and Ordonez et al., 2012 [19], for a four-wheeled robot, moving on sloped terrain.

When it comes to tracking control methods for wheeled mobile robots, the chosen approach should be based on the robot’s kinematic structure. For example, in the case of differentially driven robots, where appropriately selected motion parameters result in negligible wheel slip, controlling their motion may consist only in controlling the angular velocities of the driven wheels. It is then necessary to measure the angular velocity of the spin of the wheels and it is crucial to ensure high accuracy in terms of wheel control. Such an approach to the motion control of a four-wheeled omnidirectional robot was implemented, for example, in the study by Szeremeta and Szuster, 2022 [20].

In contrast, for robots that typically operate in wheel-slip conditions, such as skid-steered robots, which are the subject of this study, it is necessary to either compensate for wheel slip through the control system or employ control strategies based on the robot’s actual velocities or pose. Therefore, it is necessary to measure the actual motion parameters of the robot. One of the methods of measuring such parameters involves the use of inertial measurement units (IMUs), which is discussed in the work by Trojnacki and Dąbek, 2015 [21].

Another attractive method is the estimation of the robot’s motion parameters, in particular using various variants of the Kalman filter. For example, in the article by Tsunashima et al., 2006 [22], an extended Kalman filter and a multiple model approach were used, which considers various modes of the tire system, depending on the road friction. While the multiple model approach works well for motor vehicles on selected types of ground surface (and in typical temperature and humidity conditions), it cannot be applied to skid-steered wheeled robots that may navigate within a much more diverse environment, including those typical of road and off-road vehicles, as well as on ground surfaces found in buildings.

Therefore, this work uses inertial measurement units to estimate the robot’s motion parameters, which enables independence from the type of ground surface to be achieved when estimating, for e.g., the wheel slip velocity.

The fact is, however, that the estimation of the velocities of the robot is subject to error, which can lead to the accumulation of errors over time. Therefore, in order to avoid this problem, it is worth taking into account the robot’s actual pose when controlling its movement. The application of a pose controller within the control system enables high accuracy in of terms of the robot’s movement to be achieved, which was demonstrated, for e.g., in the work by Cao, 2010 [23], as well as Hassani and Rekik, 2023 [24], where the authors propose the use of such a method of motion control for a planetary exploration rover.

Based on the review of the available publications, it can be observed that there is a lack of studies where the wheel slip velocity is determined and used as the basis for improving the accuracy of tracking control. This is, therefore, one of the issues that is the subject of this article.

Various control methods can be used for trajectory tracking control of wheeled mobile robots, depending on the robot’s kinematic structure and its operating conditions. In some cases, it may be sufficient to use simple control systems based on linear PID controllers, which was, among others, the subject of the MSc thesis by Locardi, 2020 [10]. It is also possible to use the sliding mode control technique with a cascade of PID controllers, as in the article by Xiong et al., 2022 [25].

If the robot’s working conditions are variable, for example, due to changes in the type of ground or the weight of the transported load, it is necessary to use more advanced control methods that enable, for example, compensation for the non-linearity of the robot’s dynamics. For this purpose, adaptive control strategies can be applied, as demonstrated in the work by Cui et al., 2023 [26], Ye and Wang, 2020 [27], as well as Liang et al., 2023 [28], where an adaptive robust control approach was used. Another attractive technique may be model predictive control, which is the subject of the MSc thesis by Dorbetkhany, 2023 [29], papers by Wang et al., 2024 [30] and Li et al., 2023 [31], as well as being discussed in the review article by Harasim and Trojnacki, 2016 [32].

Artificial intelligence techniques can also be used for the tracking control of wheeled robots. They have an additional advantage over classical methods in that they can provide resistance to non-parametric disturbances resulting from inaccuracies in the structure of the robot model. One of the artificial intelligence techniques applied to motion control of wheeled mobile robots is reinforcement learning, which is the subject of the work by Rybczak et al., 2024 [33], and Zhang et al., 2022 [34]. Other possible techniques include fuzzy logic systems, which were used in the work by Pandiaraj and Muralidharan, 2023 [35], Hendzel and Trojnacki, 2023 [36], and also Mehmood, 2021 [37], where a robust fuzzy sliding mode controller was used. In turn, tracking control using artificial neural networks was proposed in the paper by Hendzel and Trojnacki, 2015 [38], for a four-wheeled robot, and in the paper by Hassan and Saleem, 2022 [39], for controlling a differentially driven robot. It is also possible to use various combinations of artificial intelligence techniques or to combine them with classical methods. For instance, in the article by Yue et al., 2022 [40], fuzzy logic systems and a model predictive control strategy were used, while in the publication by Zhao et al., 2023 [41], a self-organizing fuzzy neural network was used for trajectory tracking control of an omnidirectional mobile robot.

The above examples indicate that in order to obtain high accuracy in term of the trajectory tracking control of skid-steered wheeled mobile robots, moving in wheel-slip conditions, it is necessary to use advanced control methods, including those based on the robot dynamics model. It is also noticed that there is a lack of research works in which only the drive unit model of a robot is used within the control system.

Therefore, the question that this article aims to answer is whether it is possible to achieve a significant improvement in the motion accuracy of a skid-steered robot using simpler control system solutions. In particular, it is essential to address whether this can be achieved by implementing wheel velocity control, through the application of wheel-slip compensation, or by incorporating a drive unit model into the robot’s control system.

Most publications on the trajectory tracking control of wheeled mobile robots focus on the issue of motion accuracy. However, there are fewer publications on the issue of energy efficiency. The issue of energy efficiency in regard to mobile robots is the subject of review articles by Wu et al., 2023 [42], and Zhang et al., 2020 [5]. In the publication by Effati et al., 2023 [43], the issue of energy-optimal trajectories for skid-steered rovers is discussed. In the paper by Guo et al., 2024 [44], the issue of path tracking and an energy efficiency coordination control strategy for a skid-steered mobile robot are analyzed.

The conducted review of the state of the art indicates that there are a small number of publications that discuss both the issue of accuracy of the trajectory tracking control of wheeled mobile robots and their energy efficiency. Therefore, the correlation between the accuracy of trajectory tracking control and energy efficiency should be investigated.

3. Model of the Robot

The subject of this study is a mobile robot with non-steered wheels, called PIAP GRANITE, as shown in Figure 1a. A characteristic feature of this type of robot is that wheel slip always occurs when it turns or rotates in place. Modeling of this robot requires the analysis of both its kinematics and dynamics. In order to faithfully reflect the robot’s dynamic properties, it is also advisable to include a model of its drive units.

3.1. Kinematics of the Robot

The kinematic structure of the robot is illustrated in Figure 1b. The robot includes the following basic components: 0—mobile platform, 1–4—wheels. Additionally, the following designations are used for the ith wheel: $A_i$—center, $r_i$ (m)—radius, $L$ (m)—the distance between the front and back axles, $W$ (m)—track width, and ${\theta }_i$ (rad)—spin angle. A moving coordinate system $\left\{R\right\}$ is introduced for the robot, with its origin located at the midpoint of the distance between the center of the wheels and which is marked as point $R$.

The robot’s motion takes place on a horizontal plane, associated with a stationary coordinate system $\left\{O\right\}$. The robot is assumed to be moving in the ${}^{O}x{}^{O}y$ plane, meaning that the roll and pitch of the mobile platform during its motion are neglected.

The robot’s motion can be described using the vectors of generalized velocities, respectively, in regard to $\left\{O\right\}$ and $\left\{R\right\}$ coordinate systems:

$${}_{ }{}^O\dot{\mathit{q}}={\left[{}_{ }{}^O{v}_{Rx},{}_{ }{}^O{v}_{Ry}, {}_{ }{}^O{\omega }_{0z}\right]}^T, {}_{ }{}^R\dot{\mathit{q}}={\left[{}_{ }{}^R{v}_{Rx},{}_{ }{}^R{v}_{Ry}, {}_{ }{}^R{\omega }_{0z}\right]}^T,$$

(1)

where ${}_{ }{}^O{v}_{Rx}={}_{ }{}^O{\dot{x}}_R (\mathrm{m}/\mathrm{s}),{}_{ }{}^O{v}_{Ry}={}_{ }{}^O{\dot{y}}_R$ (m/s) are coordinates of the linear velocity vector of point $R$ of the robot and ${}_{ }{}^O{\omega }_{0z}={}_{ }{}^O{\dot{\phi }}_{0z}$ (rad/s) is the coordinate of the angular velocity vector of the robot, both within the stationary coordinate system $\left\{O\right\}$; ${}_{ }{}^R{v}_{Rx} (\mathrm{m}/\mathrm{s}),{}_{ }{}^R{v}_{Ry}$ (m/s), and ${}_{ }{}^R{\omega }_{0z}$ (rad/s) are similar coordinates of the velocity vectors, but are described in regard to the moving coordinate system $\left\{R\right\}$ associated with the robot.

These two vectors of generalized velocities meet the following relationship:

$${}_{ }{}^O\dot{\mathit{q}}={}_R{}^O\mathit{J} {}_{ }{}^R\dot{\mathit{q}},$$

(2)

where matrix ${}_R{}^O\mathit{J}$ has the following form:

$${}_R{}^O\mathit{J}=\left[\begin{array}{ccc}\cos \left({}^O{\phi }_{0z}\right)& -\sin \left({}^O{\phi }_{0z}\right)& 0\\ \sin \left({}^O{\phi }_{0z}\right)& \cos \left({}^O{\phi }_{0z}\right)& 0\\ 0& 0& 1\end{array}\right].$$

(3)

The velocity vectors of the points $A_i$, which are the centers of the robot’s wheels, can be calculated from the distribution of the velocity vectors. An example of the distribution of the velocity vectors of these points during the turning of the robot is illustrated in Figure 1b. These vectors can be determined based on the known linear velocity vector ${}^R{\mathit{v}}_R$ and the angular velocity vector ${}^R{\mathit{\omega }}_0$ from the equation:

$${}^R{\mathit{v}}_{Ai}={}^R{\mathit{v}}_R+{}^R{\mathit{\omega }}_0\times {}^R{\mathit{r}}_{Ai},$$

(4)

where ${}^R{\mathit{r}}_{Ai}={\left[{}^R{x}_{Ai},{}^R{y}_{Ai},{}^R{z}_{Ai}\right]}^T$ describes the position of point $A_i$ in the robot’s coordinate system $\left\{R\right\}$, $i=1,\dots ,n$ and $n=4$ represents the number of wheels.

Assuming that the robot’s mobile platform is in plane motion, the relationships between the projections of the velocities of point $A_i$ on the ${}_{ }{}^{R}x$ and ${}_{ }{}^{R}y$ axes can be described as a function of the generalized velocities within the robot’s coordinate system $\left\{R\right\}$ in the form:

$$\left[\begin{array}{c}\begin{array}{c}{}^R{v}_{Alx}\\ {}^R{v}_{Arx}\\ {}^R{v}_{Afy}\end{array}\\ {}^R{v}_{Aby}\end{array}\right]=\left[\begin{array}{c}\begin{array}{cc}1& -W/2\\ 1& W/2\\ 0& L/2-{}^R{x}_C\end{array}\\ \begin{array}{cc}0& -L/2-{}^R{x}_C\end{array}\end{array}\right]\left[\begin{array}{c}{}^R{v}_{Rx}\\ {}^R{\omega }_{0z}\end{array}\right]=\left[\begin{array}{cc}\begin{array}{c}0\\ 0\\ 1\\ 1\end{array}& \begin{array}{c}-W/2+{}^R{y}_C\\ W/2+{}^R{y}_C\\ L/2\\ -L/2\end{array}\end{array}\right]\left[\begin{array}{c}{}^R{v}_{Ry}\\ {}^R{\omega }_{0z}\end{array}\right],$$

(5)

where the designations for the individual wheels mean: $l$—left, $r$—right, $f$—front, and $b$—back.

Moreover, it is assumed that $l\in {\mathbb{L}}_w$, $r\in {\mathbb{R}}_w$, $f\in {\mathbb{F}}_w$, and $b\in {\mathbb{B}}_w$, where the sets ${\mathbb{L}}_w=\{1, 3\}$, ${\mathbb{R}}_w=\{2, 4\}$, ${\mathbb{F}}_w=\{1, 2\}$, and ${\mathbb{B}}_w=\{3, 4\}$ introduce indexes for the individual wheels, which relates to the left, right, front, and back, respectively.

From the perspective of robot motion control, the longitudinal velocities of the centers of the wheels are crucial, because they are the basis for controlling the motion of the individual wheels. Therefore, the key relationship is:

$$\left[\begin{array}{c}{}_{ }{}^R{v}_{Alx}\\ {}_{ }{}^R{v}_{Arx}\end{array}\right]=\left[\begin{array}{cc}1& -W/2\\ 1& W/2\end{array}\right]\left[\begin{array}{c}{}_{ }{}^R{v}_{Rx}\\ {}_{ }{}^R{\omega }_{0z}\end{array}\right],$$

(6)

which can alternatively be written as:

$${}_{ }{}^R{v}_{Aix}={}_{ }{}^R{v}_{Rx}-{}_{ }{}^R{\omega }_{0z}\mathrm{s}\mathrm{g}\mathrm{n}({}_{ }{}^R{y}_{Ai})W/2.$$

(7)

In general, the robot moves in wheel-slip conditions, primarily due to the design of the kinematic structure. Namely, employing only non-steered wheels in the robot’s design leads to slippage when turning or rotating in place.

Thus, the velocities of the centers of the wheels arise from the angular velocities of their spin and the slip velocities, according to the relationship:

$${}^R{\mathit{v}}_{Ai}={\left[{}^R{v}_{Aix},{}^R{v}_{Aiy},{}^R{v}_{Aiz}\right]}^T={\left[{\omega }_i{r}_i+{}^R{v}_{Six},{}^R{v}_{Siy},0\right]}^T,$$

(8)

where ${}_{ }{}^R{v}_{Six}$ (m/s) and ${}_{ }{}^R{v}_{Siy}$ (m/s) are the longitudinal and lateral coordinates of the wheel slip velocity vector within the moving coordinate system $\left\{R\right\}$, respectively.

In turn, the slip velocity vector for the ith wheel is equal to:

$${{}^R\mathit{v}}_{Si}={\left[{}^R{v}_{Six},{}^R{v}_{Siy},{}^R{v}_{Siz}\right]}^T={\left[{}^R{v}_{Aix}-{\omega }_i{r}_i,{}^R{v}_{Aiy},0\right]}^T.$$

(9)

The measure of the longitudinal slip for the ith wheel is the longitudinal slip ratio ${\lambda }_i$ (%), which is defined in various ways in the literature. In particular, the actual value of the longitudinal slip ratio ${\lambda }_i$ (%) is determined as the ratio of the longitudinal slip velocity component ${}^R{v}_{Six}$ (m/s) to the circumferential velocity ${\omega }_i{r}_i$ (m/s) during acceleration, or to the longitudinal velocity of the wheel’s geometric center ${}^R{v}_{Aix}$ (m/s) during braking. This leads to the following relationship:

$${\lambda }_i=\left\{\begin{array}{cc}0& \mathrm{i}\mathrm{f} \max (\left|{}^R{v}_{Aix}\right|,\left|{\omega }_i{r}_i\right|)=0,\hfill \\ -{}^R{v}_{Six}/\max (\left|{}^R{v}_{Aix}\right|,\left|{\omega }_i{r}_i\right|)& \mathrm{otherwise}.\hfill \end{array}\right.$$

(10)

It is assumed that the longitudinal slip ratio ${\lambda }_i$ (%) can take values within the range $〈-100\%,100\%〉$ or equivalently $〈-\mathrm{1,1}〉$.

In turn, the angle ${\alpha }_i$ (rad) between the velocity vector ${}^R{\mathit{v}}_{Ai}$ and the plane of the ith wheel is a measure of the lateral slip. Assuming the coordinate system according to the ISO convention, the actual value of the lateral slip angle ${\alpha }_i$ (rad) is calculated from the formula:

$${\alpha }_i=\left\{\begin{array}{cc}0& \mathrm{i}\mathrm{f}\left|{}^R{v}_{Aiy}\right|=0,\hfill \\ \mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}2\left({}^R{v}_{Aiy},\left|{}^R{v}_{Aix}\right|\right)& \mathrm{otherwise}.\hfill \end{array}\right.$$

(11)

In order to determine the angular velocities of the wheels ${\omega }_i$ (rad/s), the inverse kinematics problem is solved, which is applicable within the robot’s motion control system. On the other hand, having measured the angular velocities of the wheels ${\omega }_i$ (rad/s), one can calculate the velocities ${}_{ }{}^R{v}_{Rx}$ (m/s) and ${}_{ }{}^R{\omega }_{0z}$ (rad/s) for the mobile platform, that is solve the forward kinematics problem.

Both problems can be uniquely solved in the case of the robot’s motion without wheel slip. In the case of motion under wheel-slip conditions, the conversion between the velocities ${}_{ }{}^R{v}_{Rx}$ (m/s) and ${}_{ }{}^R{\omega }_{0z}$ (rad/s) of the mobile platform and the angular velocity of the wheels ${\omega }_i$ (rad/s), or vice versa, is often performed while neglecting the impact of wheel slip. This might result from the measurement of only the angular velocity of the wheels ${\omega }_i$ (rad/s), i.e., from the lack of measurement of the actual robot’s motion velocities, which is necessary to determine wheel slip.

3.2. Dynamics of the Robot

It is assumed that the robot is subject to the force of gravity ${}_{ }{}^R\mathit{G}={\left[{{}_{ }{}^{R}G}_x, {{}_{ }{}^{R}G}_y, {{}_{ }{}^{R}G}_z\right]}^T$, as well as the forces and moments of force resulting from the wheel–ground interaction. These reaction forces when reduced to the centers of the wheels are equal to the vectors ${{}_{ }{}^R\mathit{F}}_{Ai}={\left[{{}_{ }{}^{R}F}_{Aix}, {{}_{ }{}^{R}F}_{Aiy}, {{}_{ }{}^{R}F}_{Aiz}\right]}^T$ and ${{}_{ }{}^R\mathit{T}}_{Ai}={\left[{{}_{ }{}^{R}T}_{Aix}, {{}_{ }{}^{R}T}_{Aiy}, {{}_{ }{}^{R}T}_{Aiz}\right]}^T$, respectively.

The force of gravity is equal to ${}_{ }{}^R\mathit{G}=m_R{}_{ }{}^R\mathit{g}$, where $m_R$ (kg) is the total mass of the robot and ${}_{ }{}^R\mathit{g}={\left[{{}_{ }{}^{R}g}_x, {{}_{ }{}^{R}g}_y, {{}_{ }{}^{R}g}_z\right]}^T$ is a gravitational acceleration vector. The magnitude of this vector is $‖{}_{ }{}^R\mathit{g}‖=g=9.81$ (m/s²), where $g$ denotes the gravitational acceleration.

In turn, the location of the robot’s center of mass, to which the gravitational force vector ${}_{ }{}^R\mathit{G}$ is applied, is described by the position vector: ${}_{ }{}^R{\mathit{r}}_{CM}={\left[{{}_{ }{}^{R}x}_{CM}, {{}_{ }{}^{R}y}_{CM}, {{}_{ }{}^{R}z}_{CM}\right]}^T$.

The components of the reaction forces ${{}_{ }{}^R\mathit{F}}_{Aix}$ and ${{}_{ }{}^R\mathit{F}}_{Aiy}$, for $i=1, \dots , n$, acting in the plane of wheel–ground contact, are shown in Figure 2.

The Newton–Euler formalism can be used to describe the robot’s dynamics. The dynamic equations of motion for the discussed force system are written in terms of the coordinate system associated with the robot $\left\{R\right\}$, in the following form:

$$m_R{}^R{\mathit{a}}_{CM}={\sum }_{i=1}^n{}^R{\mathit{F}}_{Ai}+m_R{}_{ }{}^R\mathit{g},$$

(12)

$${\mathit{I}}_R{}^R{\mathit{\epsilon }}_0+{}^R{\mathit{\omega }}_0\times {\mathit{I}}_R{}^R{\mathit{\omega }}_0={\sum }_{i=1}^n\left({{}_{ }{}^R\mathit{T}}_{Ai}+\left({{}_{ }{}^R\mathit{r}}_{Ai}-{{}_{ }{}^R\mathit{r}}_{CM}\right)\times {{}_{ }{}^R\mathit{F}}_{Ai}\right),$$

(13)

where ${\mathit{I}}_R$ is an inertia tensor in terms of the robot, ${}^R{\mathit{a}}_{CM}$ is an acceleration vector in terms of the robot’s center of mass, and ${}^R{\mathit{\epsilon }}_0$ is an angular acceleration vector in terms of the robot’s mobile platform.

It is assumed that the mobile platform of the robot is a rigid body, the ground on which the robot moves is horizontal (i.e., ${}_{ }{}^R\mathit{g}={\left[0, 0,-g\right]}^T$), and the contact points of the wheels with the ground lie in one plane. Thus, the normal components of the ground reaction forces can be determined for the analyzed robot based on the following relationships:

$$\begin{array}{c}{{}^{R}F}_{Aiz}=m_R\left(g{}_{ }{}^R{y}_{CM}-{}^R{a}_{CMy}\left(r+{}^R{z}_{CM}\right)\right)\mathrm{s}\mathrm{g}\mathrm{n}\left({}^R{y}_{Ai}\right)/(2\mathrm{W})+\\ +m_R\left(g{}^R{x}_{CM}-{}^R{a}_{CMx}\left(r+{}^R{z}_{CM}\right)\right)\mathrm{s}\mathrm{g}\mathrm{n}\left({}^R{x}_{Ai}\right)/(2\mathrm{L})+m_{R}g/4.\end{array}$$

(14)

In turn, the normal components of the ground reaction forces ${{}^{R}F}_{Aiz}$ (N) are related to the radial deformation of the tire $∆r_i$ (m), according to the following dependence:

$${{}^{R}F}_{Aiz}=k_r∆r_i=> ∆r_i={{}^{R}F}_{Aiz}/k_r,$$

(15)

where $k_r$ (N/m) represents the radial stiffness of the tire.

In order to determine the tangential components of the ground reaction forces resulting from the interaction of the rolling wheels with the ground, i.e., ${{}_{ }{}^{R}F}_{Aix}$ (N) and ${{}_{ }{}^{R}F}_{Aiy}$ (N), the tire model should be used.

These components depend on the current values of the longitudinal slip ratio, ${\lambda }_i$ (%), and the lateral slip angle, ${\alpha }_i$ (rad).

During the turning of the robot, significant wheel slip may occur, including so-called combined slip, which involves both longitudinal and lateral slip. Particularly high values in terms of lateral slip can occur when the robot rotates in place, as the lateral slip angle ${\alpha }_i$ (rad) may approach values close to ±π/4 rad (±45 deg). Therefore, to properly model the conditions of the wheel–ground interaction, it is necessary to use a tire model for the combined slip case.

One of the popular tire models used for simulating the motion of wheeled vehicles is the Dugoff model, which is discussed, for example, in the paper by Bhoraskar and Sakthivel, 2017 [45], where it is compared to the empirical model developed by Pacejka [13]. This model is both relatively simple and capable of effectively representing the longitudinal and lateral force characteristics, including scenarios with combined slip.

According to the Dugoff tire model and for the purpose of the ISO coordinate system, the tangential coordinates of the ground reaction forces are determined according to the relationships:

$${{}_{ }{}^{R}F}_{Aix}={\displaystyle \frac{C_a{\lambda }_i}{1-{\lambda }_i}}f(z_i), {{}_{ }{}^{R}F}_{Aiy}=-{\displaystyle \frac{C_b\tan ({\alpha }_i)}{1-{\lambda }_i}}f(z_i),$$

(16)

$$f(z_i)=\left\{\begin{array}{cc}{z}_i(2-z_i)& \mathrm{i}\mathrm{f} z_i<1,\hfill \\ 1& \mathrm{otherwise}.\hfill \end{array}\right., z_i={\displaystyle \frac{{\mu }_i{{}^{R}F}_{Aiz}\left(1-{\lambda }_i\right)}{2\sqrt{{\left(C_a{\lambda }_i\right)}^2+{\left(C_b\tan ({\alpha }_i)\right)}^2}}},$$

(17)

$${\mu }_i={\mu }_0\left(1-A_S{V}_{Si}\right),V_{Si}={}^R{v}_{Aix}\sqrt{{\lambda }_i^2+{\tan }^2({\alpha }_i)},$$

(18)

where $C_a$ (N/rad) and $C_b$ (N/rad) are the longitudinal and lateral stiffness of the tires, respectively, $A_S$ (s/m) and $V_{Si}$ (m/s) are the friction reduction factor and magnitude, respectively, and ${\mu }_i$ (–) is the maximum friction coefficient for the ith wheel.

Furthermore, it is assumed that the maximum friction scaling coefficient, ${\mu }_0$ (–), is equal to the coefficient of the static friction, ${\mu }_s$ (–), between the tire and the ground.

The adopted parameters for the tire model should ensure that for the longitudinal slip coefficient ${\lambda }_i={\lambda }_p$ (%) the maximum value of the longitudinal component of the ground reaction force ${{}_{ }{}^{R}F}_{Aix}$ (N) is close to ${\mu }_p{{}^{R}F}_{Aiz}$, where ${\lambda }_p$ (%) denotes the value of the longitudinal slip ratio corresponding to the value of the maximum tire adhesion coefficient, ${\mu }_p$ (–).

Moreover, to achieve the maximum value of longitudinal slip, i.e., $\left|{\lambda }_i\right|={\lambda }_{max}=100\%$, the force component ${{}_{ }{}^{R}F}_{Aix}$ (N) should reach a value close to do ${\mu }_k{{}^{R}F}_{Aiz}$, where ${\mu }_k$ (–) is a coefficient of the kinetic friction between the tire and the ground.

When the actual values of the tangential coordinates of the ground reaction forces are known, that is ${{}_{ }{}^{R}F}_{Aix}$ (N) and ${{}_{ }{}^{R}F}_{Aiy}$ (N), it is possible to calculate the tire adhesion coefficients in the longitudinal and lateral directions, that is ${\mu }_{ix}$ (–) and ${\mu }_{iy}$ (–), respectively. To do this, the following formulas can be used:

$${\mu }_{ix}={{}_{ }{}^{R}F}_{Aix}/{{}_{ }{}^{R}F}_{Aiz}, {\mu }_{iy}={{}_{ }{}^{R}F}_{Aiy}/{{}_{ }{}^{R}F}_{Aiz}.$$

(19)

When the robot is turning or rotating in place, an additional torque component due to friction forces arises, that is ${{}^{R}T}_{Aiz}$ (Nm). This component is proportional to the friction coefficient, the normal coordinates of the ground reaction forces ${{}^{R}F}_{Aiz}$ (N), and the size of the tire–ground contact area.

Let us assume that the tire–ground contact area has a rectangular shape, with dimensions $a_i\times b_i$.

The length $a_i$ (m), measured along the tire, results from its radial deformation $∆r_i$ (m) and can be calculated using the following equations:

$$a_i=2 r\sin ({\vartheta }_i), {\vartheta }_i=\arccos ((r-∆r_i)/r).$$

(20)

In turn, $b_i=b$ (m) corresponds to the width of the tire tread.

Based on this, the average distance from the center of the rectangle can be calculated using the formula:

$$R_{ci}=\sqrt{a_i^2+b_i^2}/4.$$

(21)

Then, the vertical component of the friction torque can be computed based on the approximate relationship, assuming the kinetic friction coefficient ${\mu }_k$ (–) for the calculation:

$${{}^{R}T}_{Aiz}=-{\mu }_k{{}_{ }{}^{R}F}_{Aiz}{R}_{ci}\mathrm{s}\mathrm{g}\mathrm{n}({}^R{\omega }_{0z}),$$

(22)

where the function $\mathrm{s}\mathrm{g}\mathrm{n}({}^R{\omega }_{0z})$ can be approximated by the expression $\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{h}(k_{\omega }{}^R{\omega }_{0z})$ for $k_{\omega }>0$.

When the components of the forces and moments of the force acting on the wheels of the robot are known, as well as the gravitational force, it is possible to calculate the linear and angular accelerations related to the robot’s motion, based on the dynamic equations of motion.

When describing the robot’s motion in the ${}_{ }{}^{O}x{}_{ }{}^{O}y$ plane in terms of the fixed coordinate system $\left\{O\right\}$, the following dependencies are crucial:

$${}^R{a}_{CMx}={\sum }_{i=1}^n{{}^{R}F}_{Aix}/m_R, {}^R{a}_{CMy}={\sum }_{i=1}^n{{}^{R}F}_{Aiy}/m_R,$$

(23)

$${}^R{\epsilon }_{0z}=({\sum }_{r\in {\mathbb{R}}_w}{}^R{F}_{Arx}{w}_r-{\sum }_{l\in {\mathbb{L}}_w}{{}^{R}F}_{Alx}{w}_l+{\sum }_{f\in {\mathbb{F}}_w}{{}^{R}F}_{Afy}{l}_f-{\sum }_{b\in {\mathbb{B}}_w}{{}^{R}F}_{Aby}{l}_b+{\sum }_{i=1}^n{{}^{R}T}_{Aiz})/I_{Rz},$$

(24)

where ${\mathbb{L}}_w=\{1,3\}$, ${\mathbb{R}}_w=\{2,4\}$, ${\mathbb{F}}_w=\{1,2\}$, and ${\mathbb{B}}_w=\{3,4\}$.

The distance in terms of the center of these wheels from the robot’s center of mass, measured along the ${}_{ }{}^{R}x$ and ${}_{ }{}^{R}y$ axes, are equal to, respectively: $l_f=(L/2-{{}_{ }{}^{R}x}_{CM})$, $l_b=(L/2+{{}_{ }{}^{R}x}_{CM})$, $w_r=\left(W/2+{{}_{ }{}^{R}y}_{CM}\right)$, and $w_l=\left(W/2-{{}_{ }{}^{R}y}_{CM}\right)$.

In turn, the spin of the wheels can be described using the following dynamic equation:

$$I_{Wy}\hspace{0.33em}{\epsilon }_i={\tau }_i-{{}_{ }{}^{R}F}_{Aix}r-{{}_{ }{}^{R}F}_{Aiz}r{f}_r\mathrm{sgn}({\omega }_i),$$

(25)

where $I_{Wy}$ (kg·m²) is a mass moment of inertia of the wheel relative to its axis of rotation, ${\tau }_i$ (Nm) is the driving torque, $f_r$ (–) is a coefficient of rolling resistance, while ${\omega }_i$ (rad/s) and ${\epsilon }_i$ (rad/s²) are the angular velocity and acceleration of this wheel, respectively.

Based on Equation (25), the angular accelerations of the wheel’s spin, ${\epsilon }_i$, can be determined (forward dynamics problem) as follows:

$${\epsilon }_i=\left({\tau }_i-{{}_{ }{}^{R}F}_{Aix}{r}_i-{{}_{ }{}^{R}F}_{Aiz}{r}_i{f}_r\mathrm{sgn}({\omega }_i)\right)/I_{Wy},$$

(26)

or by using the driving torques, ${\tau }_i$ (inverse dynamics problem):

$${\tau }_i=I_{Wy}\hspace{0.33em}{\epsilon }_i+{{}_{ }{}^{R}F}_{Aix}r+{{}_{ }{}^{R}F}_{Aiz}r{f}_r\mathrm{sgn}({\omega }_i).$$

(27)

Instead of the function $\mathrm{sgn}({\omega }_i)$, one can use, for example, the function $\tanh (\beta {\omega }_i)$ for $\beta >0$ (–), which enables a gradual change in the rolling resistance moment.

3.3. The Robot’s Drive Units

The previously described dynamics model of the robot can be extended by incorporating a drive unit model, enabling a more precise reflection of all the aspects influencing the robot’s operation.

The analyzed robot is equipped with drive units that enable the wheels to be driven independently. Each drive unit consists of an identical DC motor, encoder, and gear system. The maximum motor input voltage for the drive units is $U_{max}$ (V), which corresponds to the maximum motor rotational speed, $n_{max}$ (rpm), and the resulting maximum angular velocity of the wheel, ${\omega }_{max}$ (rad/s).

In the case of drive units, both a forward and an inverse model can be distinguished, which have analogous meanings, as in the case of a model of robot dynamics.

The forward model of the drive units can be combined with the forward model of the robot dynamics, enabling the determination of the driving torques ${\tau }_i$ (Nm) acting on the robot, based on the motor voltage inputs $U_i$ (V) generated by the control system.

A forward model of the drive units can be described by the equations:

$$d{I}_i/dt=\left(U_i-k_e{n}_d{\omega }_i-R_d{I}_i\right)/L_d, {\tau }_i={\eta }_d{n}_d{k}_m{I}_i,$$

(28)

where $I_i$ (A) is the rotor current, $L_d$ (mH) and $R_d$ (Ω) are the inductance and the resistance of the rotor, respectively, $k_e$ (Vs/rad) is the electromotive force constant, $k_m$ (Nm/A) is a motor torque coefficient, while $n_d$ (–) and ${\eta }_d$ (–) are the gear ratio and an efficiency factor of the transmission system, respectively.

The inverse model of the drive units can be used in regard to the control system, enabling the determination of the values of the rotor current $I_i$ (A) and the motor voltage inputs $U_i$ (V) necessary for the assumed robot motion, based on the desired driving torques, ${\tau }_{di}$ (Nm).

The inverse model of the drive units can be written in the form of the equations:

$$I_i={\tau }_{di}/({\eta }_d{n}_d{k}_m), U_i=L_d\mathrm{d}{I}_i/\mathrm{d}t+k_e{n}_d{\omega }_{id}+R_d{I}_i.$$

(29)

In both cases, for evaluating the electric energy demand and potential optimization of energy efficiency, it is possible to calculate the electric power, $p_i$ (W), necessary to rotate the ith wheel with the angular velocity, ${\omega }_i$ (rad/s), and to generate the driving torque ${\tau }_i$ (Nm) for that wheel by the drive unit. When the input voltage $U_i$ (V) and the rotor current $I_i$ (A) are known, the absolute electric power $p_i$ (W) can be determined from the formula:

$$p_i=\left|U_i{I}_i\right|.$$

(30)

4. Motion Control System

4.1. Desired Motion of the Robot and the Tracking Control Problem

The analyzed robot is a nonholonomic system, which in practice means that it cannot move in the lateral direction, except for movement in regard to wheel-slip conditions. For this reason, ${{}_{ }{}^{R}v}_{Ryd}=0$ is assumed.

Therefore, the desired robot’s motion can be represented as a vector of the desired generalized velocities:

$${\mathit{v}}_d=[{{}^{R}v}_{Rxd},{{}^R\omega }_{0zd}{]}^T=[v_{Rd},{\omega }_{0zd}{]}^T,$$

(31)

where ${{}_{ }{}^{R}v}_{Rxd}=v_{Rd}$ (m/s) and ${{}_{ }{}^R\omega }_{0zd}={\omega }_{0zd}$ (rad/s) are the values of the desired linear velocity of point $R$ in the robot and the desired angular velocity of its mobile platform is $0$, both related to the moving coordinate system $\left\{R\right\}$.

In the general case of the plane motion of the mobile platform of the robot, the desired velocity $v_{Rd}$ (m/s) depends on the desired angular velocity ${\omega }_{0zd}$ (rad/s) and the radius of the curvature $R_z$ (m), according to the relationship:

$$v_{Rd}={\omega }_{0zd}{R}_z.$$

(32)

It is assumed that a wheel controller is used for trajectory tracing control of the robot. In such a case, the trajectory tracking control problem is related to the determination of the control vector, $\mathit{u}$, such that $\mathit{\omega }\to {\mathit{\omega }}_d$ for time $t\to \infty$, where ${\mathit{\omega }}_d={\left[{\omega }_{1d},\dots ,{\omega }_{nd}\right]}^T$ are the current, desired angular velocities of the wheels; and $\mathit{\omega }={\left[{\omega }_1,\dots ,{\omega }_n\right]}^T$ are the current, measured values of those velocities.

To control the robot’s wheels, it is necessary to transform the vector of the desired generalized velocities of the robot into the angular velocities of the wheels.

If the issue of wheel-slip compensation is omitted and assuming that all the robot’s wheels have the same radii, i.e., $r_i=r$ (m), the desired angular velocities of the wheels, ${\omega }_{id}$ (rad/s), can be determined using the relationship:

$$\left[\begin{array}{c}{\omega }_{ld}\\ {\omega }_{rd}\end{array}\right]={\displaystyle \frac{1}{r}}\left[\begin{array}{cc}1& -W/2\\ 1& W/2\end{array}\right]\left[\begin{array}{c}{}_{ }{}^R{v}_{Rxd}\\ {}_{ }{}^R{\omega }_{0zd}\end{array}\right], {\omega }_{id}=\left\{\begin{array}{cc}{\omega }_{ld}& \mathrm{f}\mathrm{o}\mathrm{r} i\in {\mathbb{L}}_w,\\ {\omega }_{rd}& \mathrm{f}\mathrm{o}\mathrm{r} i\in {\mathbb{R}}_w,\end{array}\right.$$

(33)

where ${\mathbb{L}}_w=\{1, 3\}$, ${\mathbb{R}}_w=\{2, 4\}$ are the indices for the left and right wheels, respectively.

The relationship (33) results from the solution of the inverse kinematics problem for the robot and leads to inaccuracy due to the occurrence of wheel slip. This inaccuracy can be compensated by taking into account the wheel slip velocities within the robot’s control system or by additionally introducing a pose controller (see for e.g., the paper by Cui et al., 2023 [26]), if the actual pose is known (e.g., measured using satellite navigation systems). Otherwise, the use of the relationship (33) can lead to tracking errors that accumulate over time.

4.2. Wheel Controller

The simplest form of a robot’s motion control system can be a linear PID controller, used to control the angular velocities of the wheels. The basis for trajectory tracking control is then the relationship between the desired generalized velocities of the robot, i.e., ${}_{ }{}^R{v}_{Rxd}$ (m/s) and ${}_{ }{}^R{\omega }_{0zd}$ (rad/s), and the desired angular velocities of the wheels, ${\omega }_{id}$ (rad/s).

4.2.1. Typical Linear Controller (Solution A)

The measure of tracking control accuracy is the error vector of the angular velocities of the wheels ${\mathit{e}}_{\omega }={\left[e_{\omega 1},\dots ,e_{\omega n}\right]}^T$, which results from the relationship:

$${\mathit{e}}_{\omega }={\mathit{\omega }}_d-\mathit{\omega }.$$

(34)

In addition, the following error tracking function can be defined:

$$\mathit{s}={\mathit{e}}_{\omega }+\Lambda {\mathit{e}}_{\theta },$$

(35)

where the error vector of the self-rotation angles of the wheels ${\mathit{e}}_{\theta }={\left[e_{\theta 1},\dots ,e_{\theta n}\right]}^T$ is obtained by integrating the error vector of the angular velocities of the wheels, ${\mathit{e}}_{\omega }$.

Then, the control law of the linear controller takes the form:

$$\mathit{u}=\mathrm{s}\mathrm{a}\mathrm{t}\left(K_u\mathit{s}\left|\mathrm{s}\mathrm{g}\mathrm{n}\left({\mathit{\omega }}_d\right)\right|,U_{max}\right),$$

(36)

where $\mathit{u}={\left[U_1,\dots ,U_n\right]}^T$ is the control vector, $K_u$ (V·s/rad) and $\Lambda$ (1/s) are the controller gains, and $U_i$ (V) is the input voltage of the motor for the ith wheel, for $i=1,\dots ,n$.

The control vector $\mathit{u}$ is limited using the saturation function, so that the individual control signals do not exceed $\pm U_{max}$ (V), i.e., the maximum input voltage of the motor.

The term $\left|\mathrm{s}\mathrm{g}\mathrm{n}\left({\mathit{\omega }}_d\right)\right|$ is introduced to force the robot to stop in case the desired angular velocities ${\mathit{\omega }}_d$ are zero and the tracking errors ${\mathit{e}}_{\omega }$ and ${\mathit{e}}_{\theta }$ are different to zero. The function $\left|\mathrm{s}\mathrm{g}\mathrm{n}\left({\mathit{\omega }}_d\right)\right|$ may lead to an unwanted chattering phenomenon. Therefore, it can be replaced by the continuous function $\left|\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{h}\left(\gamma {\mathit{\omega }}_d\right)\right|$, with an appropriate value of parameter $\gamma$.

4.2.2. Controller with Wheel-Slip Compensation (Solution B)

If wheel-slip compensation is possible, then when calculating the error vector of the angular velocities of the spin of the wheels ${\mathit{e}}_{\omega }$, the longitudinal wheel slip velocity vector ${\mathit{v}}_S={\left[{}^R{v}_{S1x},\dots ,{}^R{v}_{Snx}\right]}^T$ can additionally be taken into account, in accordance with the formula:

$${\mathit{e}}_{\omega }={\mathit{\omega }}_d-\mathit{\omega }-{\displaystyle \frac{1}{r}}{\mathit{v}}_S.$$

(37)

The relationship (37) takes into account the longitudinal component of the wheel slip velocities, ${}^R{v}_{Six}$ (m/s). It can be determined by measuring the actual longitudinal velocities of the wheels ${}^R{v}_{Aix}$ (m/s) and the angular velocities of the spin of those wheels ${\omega }_i$ (rad/s). For this purpose, inertial measurement units can be used, for e.g., according to the methodology described by Trojnacki and Dąbek, 2015 [21], in which the linear and angular velocities of the robot’s mobile platform are estimated. This approach has the advantage that it is independent of the type of ground on which the robot is moving, so it does not require the estimation of tire–ground friction coefficients, as is the case in the work by Tufano et al., 2024 [16].

The control law of the controller takes a similar form to the previous one. However, the control system shown in Figure 3 incorporates additional feedback from the longitudinal wheel slip velocity vector ${\mathit{v}}_S$, which is taken into account when calculating the tracking error, ${\mathit{e}}_{\omega }$.

It should be noted that the determination of the longitudinal slip velocities of the wheels based on Formula (9) requires the estimation of the longitudinal component of the velocity of the wheel center, ${}^R{v}_{Aix}$ (m/s).

When using inertial measurement units, the estimation of the velocities of the wheel centers can be accomplished, for e.g., by placing MEMS accelerometers above the wheels, and in terms of the determination of the velocities, by integrating the tangential components of the accelerations. It is also possible to use a single set of inertial measurement units in the central part of the robot (around point $R$, according to Figure 1b), including three-axis MEMS accelerometers and gyroscopes. Then, the velocities of the wheel centers can be calculated from Formula (4), based on the determined velocity vectors ${}^R{\mathit{v}}_R$ of point $R$ in the robot and the angular velocity of the mobile platform, ${}^R{\mathit{\omega }}_0$. If only the longitudinal component of the velocity is needed, then relationship (7) can also be used for this purpose.

4.2.3. Controller with Drive Unit Model (Solution C)

In a case where the control system includes the inverse model of the drive units but does not compensate for wheel slip, the error vector of the angular velocities of the spin of the wheels ${\mathit{e}}_{\omega }$ is defined according to Equation (34), while the error tracking function is derived from Formula (35). In this case, first the desired driving torques ${\mathit{\tau }}_d={\left[{\tau }_{1d},\dots ,{\tau }_{nd}\right]}^T$ are determined based on the following equation:

$${\mathit{\tau }}_d=K_{\tau }\mathit{s}\left|\mathrm{s}\mathrm{g}\mathrm{n}\left({\mathit{\omega }}_d\right)\right|,$$

(38)

and then, using the inverse model of the drive units, the control vector $\mathit{u}={\left[u_1,\dots ,u_n\right]}^T$ for the drive units, given in volts (V), is derived based on the following equations:

$${\mathit{I}}_d=k_d{\mathit{\tau }}_d, k_d=1/({\eta }_d{n}_d{k}_m),$$

(39)

$$\mathit{u}=\mathrm{s}\mathrm{a}\mathrm{t}\left(k_e{n}_d{\mathit{\omega }}_d+L_d{\dot{\mathit{I}}}_d+R_d{\mathit{I}}_d,U_{max}\right),$$

(40)

where the gains used in this case have the respective units $K_{\tau }$ (Nm·s/rad) and $\Lambda$ (1/s).

Equations (38) and (39) can optionally include constraints ${\tau }_{max}$ (Nm) and $I_{max}$ (A) for the maximum absolute values of the driving torque and current, respectively. For this purpose, they can be supplemented with appropriate saturation functions.

It should be noted that due to the introduction of the inverse drive unit model within the controller structure, the applied controller has a non-linear structure.

It is also worth emphasizing that this solution requires knowledge of the drive unit model and its approximate parameters. Such parameters can be obtained from the specifications of the drive units, or they can be gained as the result of an identification process.

4.2.4. Controller with Wheel-Slip Compensation and Drive Unit Model (Solution D)

The last variant of the controller incorporates both wheel-slip compensation and the inverse model of the drive units. Thus, the error vector of the angular velocities of the wheels, ${\mathit{e}}_{\omega }$, is defined according to Equation (37), an error tracking function by Dependence (35), and the control law is described by Formula (40).

In this case, as before, the desired driving torques ${\mathit{\tau }}_d$ are calculated using Equation (38) and the inverse model of the drive units, incorporated in Equations (39) and (40), is applied. As in the previous case, the proposed controller has a non-linear structure, because it incorporates the inverse model of the drive units.

4.3. Stability Analysis of the Control System

In the previous section, four versions of the control system for a skid-steered robot were presented. An important element of the control system design is the analysis of its stability. Therefore, in the case of successful verification of the selected solution during simulation tests, an analysis of the stability of the target solution should be performed before its implementation in a real robot. The final solution, depending on the available sensors, may be more advanced than the discussed options. For example, it may include a pose controller or an additional wheel control system component that allows compensation for the object’s dynamic non-linearity and provides resistance to parametric disturbances, as is the case with adaptive control systems.

The basis for the analysis of the system’s stability can be the robot model described by means of the dynamic equations of motion in the matrix-vector form, proposed in the article by Hendzel and Trojnacki, 2023 [36]:

$$\mathit{M}{\dot{\mathit{u}}}_R+{\mathit{F}}_R({\mathit{u}}_R)+{\mathit{\tau }}_z=\mathit{\tau },$$

(41)

where ${\mathit{u}}_R={\left[{}^R{v}_{Rx},{}^R{\omega }_{0\mathrm{z}}\right]}^T$ is a vector of the control signals, $\mathit{M}$ is a constant positive-definite inertia matrix, ${\mathit{F}}_R({\mathit{u}}_R)$ denotes a vector describing the robot’s dynamic non-linearities, $\mathit{\tau }$ is a vector of the driving torque, and ${\mathit{\tau }}_z$ is a vector of bounded unknown disturbances.

The adopted form of the vector ${\mathit{u}}_R$ results from the fact that the analyzed robot is a nonholonomic system, which means that it cannot perform the desired motion in the lateral direction. The matrix $\mathit{M}$ and the vector ${\mathit{F}}_R({\mathit{u}}_R)$ result from the previously described model of the robot’s dynamics and, importantly, the vector ${\mathit{F}}_R({\mathit{u}}_R)$ is a function of the wheel slip.

Moreover, by analogy to dependency (33), one can represent the angular velocities of the driven wheels as functions of the control signals ${\mathit{u}}_R$, which perform the desired trajectory, using the formula:

$$\left[\begin{array}{c}{\omega }_l\\ {\omega }_r\end{array}\right]={\displaystyle \frac{1}{r}}\left[\begin{array}{cc}1& -W/2\\ 1& W/2\end{array}\right]{\mathit{u}}_R, {\omega }_i=\left\{\begin{array}{cc}{\omega }_l& \mathrm{f}\mathrm{o}\mathrm{r} i\in {\mathbb{L}}_w,\\ {\omega }_r& \mathrm{f}\mathrm{o}\mathrm{r} i\in {\mathbb{R}}_w,\end{array}\right.$$

(42)

where ${\mathbb{L}}_w=\{1, 3\}$ and ${\mathbb{R}}_w=\{2, 4\}$ are the indices for the left and right wheels, respectively.

The above dependence is a result of the fact that in the analyzed cases the robot motion control consists in controlling the angular velocities of the wheels and it can be substituted into the dynamics model described by Formula (41).

The stability analysis for the selected wheel controller should be performed using Lyapunov stability theory, which is an ideal tool for non-linear systems. During the stability analysis, the dynamic equations of motion should be written as a function of the velocity error vector, $\dot{\mathit{s}}$.

In order to analyze the stability of the control system using the drive unit model, it is necessary to switch to voltage control. For this purpose, based on Equations (39) and (40), a simplified relationship can be assumed, in the form:

$$\mathit{\tau }=({\eta }_d{n}_d{k}_m)\left(\mathit{u}-k_e{n}_d\mathit{\omega }\right)/R_d,$$

(43)

in which, similarly to the work by An and Eun, 2022 [46], the term $L_d{\dot{\mathit{I}}}_d$ is omitted.

5. Simulation Studies

5.1. Robot and Environment-Related Parameters

For the purpose of the simulation studies using the MATLAB R2024b/Simulink package, the previously described model of the robot and the various versions of the trajectory tracking control system were implemented.

Table 1. Basic parameters of the PIAP GRANITE robot adopted for the simulation studies.

Type of the Parameters	Values of the Parameters
Dimensions	$L=0.425 \mathrm{m}$, $W=0.52 \mathrm{m}$, $r=0.0965 \mathrm{m}$, $b=0.064 \mathrm{m}$
Mass of the bodies	$m_i=1 \mathrm{k}\mathrm{g}$, $m_R=43.4 \mathrm{k}\mathrm{g}$
Center of mass coordinates	${}_{ }{}^R{x}_{CM}=-0.063 \mathrm{m}$, ${}_{ }{}^R{y}_{CM}=0 \mathrm{m}$, ${}_{ }{}^R{z}_{CM}=0.081 \mathrm{m}$
Mass moments of inertia	$I_{Wy}=0.01 \mathrm{k}\mathrm{g}·{\mathrm{m}}^2$, $I_{Rz}=2.8 \mathrm{k}\mathrm{g}·{\mathrm{m}}^2$
Tire parameters	$C_a=C_b=1000 \mathrm{N}/\mathrm{r}\mathrm{a}\mathrm{d}$, $A_S=0.4 \mathrm{s}/\mathrm{m}$, $k_r=\mathrm{18,000} \mathrm{N}/\mathrm{m}$
Drive unit parameters	$L_d=0.0823 \mathrm{m}\mathrm{H}$, $R_d=0.317 \mathsf{\Omega }$, $k_e=0.0301 \mathrm{V}\mathrm{s}/\mathrm{r}\mathrm{a}\mathrm{d}$, $k_m=0.0302 \mathrm{N}\mathrm{m}/\mathrm{A}$, $n_d=53$, ${\eta }_d=0.8$, $U_{max}=32 \mathrm{V}$, ${\omega }_{max}=15.8 \mathrm{r}\mathrm{a}\mathrm{d}/\mathrm{s}$

presents the basic parameters of the PIAP GRANITE robot used in the simulation studies.

The PIAP GRANITE research robot was designed as part of a project, the results of which are described, among others, in the work by Hendzel and Trojnacki, 2023 [36].

This robot is, therefore, well-known to the author of this paper. The adopted robot parameters are the result of the CAD model, measurements, estimations, and experimental tests, as well as from the catalogs of the drive units. The robot is driven by high-quality integrated Maxon drive units, including DC motors, planetary gears, and optical encoders. The robot is also equipped with IMU and GNSS systems, as well as a system for measuring the current consumed by the drives and the forces acting on the wheel axles.

For the simulation studies, it is assumed that the robot moves on concrete ground, which represents the typical conditions for the analyzed robot. The parameters adopted for the tire interaction with this type of ground are presented in

Table 2. Tire–ground contact parameters adopted for the simulation studies.

Ground Type	Contact Parameters
	${\mathit{\mu }}_{\mathit{s}} (–)$	${\mathit{\mu }}_{\mathit{p}} (–)$	${\mathit{\mu }}_{\mathit{k}} (–)$	${\mathit{f}}_{\mathit{r}} (–)$
Concrete	1.00	0.85	0.75	0.015

.

Furthermore, the coefficients $\beta =2$ and $k_{\omega }=2$, approximating the signum functions, are assumed.

5.2. Desired Motion and Initial Conditions

The simulation studies analyze the longitudinal motion and turning of the robot. Therefore, the desired motion path includes:

• A straight-line segment, $L_r$ (m), including the initial part related to acceleration, $l_r$ (m);
• A circular arc of radius $R_z$ (m);
• A straight-line segment, $L_h$ (m), including the final deceleration part, $l_h$ (m).

It is also assumed that due to the desired linear acceleration, the robot’s motion is divided into the following:

• Acceleration with the maximum absolute acceleration, $a_{Rmax}$ (m/s²);
• Steady motion with constant velocity, $v_{Ru}$ (m/s);
• Braking with the maximum absolute acceleration, $a_{Rmax}$ (m/s²).

By analogy, it is assumed that due to the desired angular acceleration, the robot’s turning maneuver consists of the following:

• Acceleration with the maximum absolute angular acceleration, ${\epsilon }_{0zmax}$ (rad/s²);
• Steady turning with the angular velocity ${\omega }_{0zu}=v_{Ru}/R_z$ (rad/s) (where a positive value of the turning radius $R_z$ (m) means a left turn and a negative one means a right turn);
• Deceleration with the maximum absolute angular acceleration, ${\epsilon }_{0zmax}$ (rad/s²).

As a result of the maneuver, the robot should turn according to the angle, ${\phi }_{0zmax}$ (rad).

Table 3. Desired motion parameters of the robot adopted for the simulation studies.

Desired Motion Parameters of the Robot
${\mathit{v}}_{\mathit{R}\mathit{u}}$ (m/s)	${\mathit{a}}_{\mathit{R}\mathit{m}\mathit{a}\mathit{x}}$ (m/s2)	${\mathit{L}}_{\mathit{r}}$ (m)	${\mathit{L}}_{\mathit{h}}$ (m)	${\mathit{l}}_{\mathit{r}},{\mathit{l}}_{\mathit{h}}$ (m)	${\mathit{R}}_{\mathit{z}}$ (m)	${\mathit{\omega }}_{0\mathit{z}\mathit{u}}$ (rad/s)	${\mathit{\epsilon }}_{0\mathit{z}\mathit{m}\mathit{a}\mathit{x}}$ (rad/s2)	${\mathit{\phi }}_{0\mathit{z}\mathit{m}\mathit{a}\mathit{x}}$ (rad)
0.9	6.075	1.0	2.0	0.1	−0.6	−1.5	23.37	−2/3 π

presents the key desired motion parameters of the robot that were adopted for the simulation studies.

The assumed desired motion parameters were selected based on good knowledge of the robot’s capabilities and previous experience, in such a way that the desired motion trajectory was achievable for the robot’s drive units, i.e., that they did not introduce additional errors due to limited electric power, current, or rotational speed.

Moreover, the maximum accelerations were chosen to ensure that, for the analyzed ground, they do not cause a large longitudinal slip of the wheels. This requirement is met when the tangential force is less than the available friction force, i.e., $\left|{{}^{R}F}_{Aixy}\right|<{\mu }_s{{}^{R}F}_{Aiz}$.

Figure 4 shows the time histories of the desired velocities and accelerations of the robot that were assumed for the simulation studies.

5.3. Controller Settings

The wheel controller gains were selected to achieve small tracking errors and maintain control system stability. The parameters of the individual control system solutions, adopted in order to achieve these aims, are presented in

Table 4. Parameters of the individual solutions in terms of the control system assumed for the simulation studies.

Id	Controller Description	Controller Gains
A	Typical linear controller	$K_u=10$$\mathrm{V}·\mathrm{s}/\mathrm{r}\mathrm{a}\mathrm{d}$$, \Lambda =3$${\mathrm{s}}^{-1}$
B	Controller with wheel-slip compensation	$K_u=10$$\mathrm{V}·\mathrm{s}/\mathrm{r}\mathrm{a}\mathrm{d}$$, \Lambda =3$${\mathrm{s}}^{-1}$
C	Controller with drive unit model	$K_{\tau }=10$$\mathrm{N}\mathrm{m}·\mathrm{s}/\mathrm{r}\mathrm{a}\mathrm{d}$$, \Lambda =3$${\mathrm{s}}^{-1}$
D	Controller with wheel-slip compensation and drive unit model	$K_{\tau }=10$$\mathrm{N}\mathrm{m}·\mathrm{s}/\mathrm{r}\mathrm{a}\mathrm{d}$$, \Lambda =3$${\mathrm{s}}^{-1}$

.

It should be noted that analogous gain values were assumed for the individual control systems. They differ only in the units of the gains, due to the fact that solutions C and D in terms of the control system incorporate the inverse model of the drive units.

5.4. Quality Indexes

In order to quantitatively evaluate the results of the simulation studies for the individual control system solutions, it is advisable to introduce quality indexes.

The accuracy of the desired motion execution by the wheel controller can be assessed using quality indexes related to the maximum absolute angular velocity error of the individual wheels and their average over all $n$ wheels, using the following formulas:

$$e_{\omega imax}=\underset{t\in ⟨0,\hspace{0.33em}T⟩}{\mathrm{m}\mathrm{a}\mathrm{x}}\left(\left|e_{\omega i}\right|\right), e_{\omega max}={\sum }_{i=1}^n{e}_{\omega imax}/n.$$

(44)

Cumulative quality indexes can also be calculated, referring to the entire simulation period $T$ (s), that is:

$$E_{\omega i}=\sqrt{{\displaystyle \frac{1}{T}}{\int }_0^T{\left(e_{\omega i}\right)}^2\mathrm{d}t}, E_{\omega }={\sum }_{i=1}^n{E}_{\omega i}/n.$$

(45)

In order to assess the accuracy of the desired generalized velocities of the mobile platform, namely the linear velocity ${}^R{v}_{Rx}$ (m/s) and angular velocity ${}^R{\omega }_{0z}$ (rad/s), the quality indexes defined by the following formulas are employed:

$$e_{kmax}=\underset{t\in ⟨0,\hspace{0.33em}T⟩}{\mathrm{m}\mathrm{a}\mathrm{x}}\left(\left|e_k\right|\right), E_k=\sqrt{{\displaystyle \frac{1}{T}}{\int }_0^T{\left(e_k\right)}^2\mathrm{d}t},$$

(46)

where $k\in \left\{vR,\omega 0\right\}$.

The quality indexes related to the generalized velocities can enable an assessment of the degree of wheel-slip compensation.

In turn, in order to assess the accuracy of executing the desired motion path of the characteristic point $R$ in the robot, the following quality indexes are used:

$$e_{pmax}=\underset{t\in ⟨0,\hspace{0.33em}T⟩}{\mathrm{m}\mathrm{a}\mathrm{x}}\left(\left|e_p\right|\right), E_p=\sqrt{{\displaystyle \frac{1}{T}}{\int }_0^T{\left(e_p\right)}^2\mathrm{d}t},$$

(47)

assuming that the position error is calculated as the Euclidean distance from the actual to the desired position of point $R$, according to the formula:

$$e_p=\sqrt{{\left({{}_{ }{}^{O}e}_x\right)}^2+{({{}_{ }{}^{O}e}_y)}^2}, {{}_{ }{}^{O}e}_c={{}_{ }{}^{O}c}_{Rd}-{{}_{ }{}^{O}c}_R,$$

(48)

where $c=\{x,y\}$.

In addition, in order to assess the energy efficiency of the robot, the following quality indexes are introduced relating to the electric power and electric energy consumed by the robot’s drive units:

$$p_{imax}=\underset{t\in ⟨0,\hspace{0.33em}T⟩}{\mathrm{m}\mathrm{a}\mathrm{x}}\left(p_i\right), p_{max}={\sum }_{i=1}^n{p}_{imax}/n,$$

(49)

$$E_{di}={\int }_0^T{p}_i\mathrm{d}t, E_d={\sum }_{i=1}^n{E}_{di}/n,$$

(50)

where $p_i$ (W) is the electric power needed for the ith drive unit, determined according to Relationship (30), and $E_{di}$ (J) is the electric energy consumed by the same drive.

Similar to the quality indexes related to the tracking error of the wheels, the quality indexes related to energy efficiency also include those expressing the maximum and cumulative values for a single wheel, that is $p_{imax}$ (W) and $E_{di}$ (J), as well as the average for $n$ wheels, that is $p_{max}$ (W) and $E_d$ (J).

In order to make the evaluation of the control system performance independent of the time of motion, a quality index can also be introduced, expressing the average electric power required for the robot’s movement in the assumed period $T$ (s), according to the following equation:

$$p_{avg}=E_d/T.$$

(51)

5.5. The Results of the Simulation Studies

Table 5. Quality indexes obtained for the analyzed solutions in terms of the control system.

Controller Solution	Wheel-Slip Compens	Drive Unit Model	Quality Indexes Obtained for Particular Solutions of the Controller
			${\mathbf{e}}_{\mathit{\omega }\mathit{m}\mathit{a}\mathit{x}}$ (rad/s)	${\mathit{E}}_{\mathit{\omega }}$ (rad/s)	${\mathit{e}}_{\mathit{v}\mathit{R}\mathit{m}\mathit{a}\mathit{x}}$ (m/s)	${\mathit{E}}_{\mathit{v}\mathit{R}}$ (m/s)	${\mathit{e}}_{\mathit{\omega }0\mathit{m}\mathit{a}\mathit{x}}$ (rad/s)	${\mathit{E}}_{\mathit{\omega }0}$ (rad/s)	${\mathit{e}}_{\mathit{p}\mathit{m}\mathit{a}\mathit{x}}$ (m)	${\mathit{E}}_{\mathit{p}}$ (m)	${\mathit{p}}_{\mathit{m}\mathit{a}\mathit{x}}$ (W)	${\mathit{E}}_{\mathit{d}}$ (J)
A	No	No	1.04	0.306	0.140	0.0304	0.853	0.314	2.090	1.073	89.4	116.8
B	Yes	No	1.16	0.363	0.102	0.0255	0.454	0.096	0.172	0.108	121.7	135.2
C	No	Yes	0.48	0.122	0.100	0.0169	0.826	0.306	2.121	1.094	98.8	114.2
D	Yes	Yes	1.60	0.381	0.076	0.0129	0.598	0.135	0.380	0.237	119.7	149.2

presents the quality indexes obtained from the simulation studies for all the analyzed control system solutions and the desired motion parameters illustrated in Figure 4. In terms of compensation for longitudinal wheel slip, the determination of the longitudinal wheel slip velocities ${\mathit{v}}_S={\left[{}^R{v}_{S1x},\dots ,{}^R{v}_{Snx}\right]}^T$, based on the estimation of the linear and angular velocity vectors of the mobile platform, are assumed.

By analyzing the obtained values in regard to the quality indexes for individual solutions in terms of the control system for a four-wheeled skid-steered robot, the following conclusions can be reached:

• The highest accuracy in regard to controlling the angular velocities of the robot’s wheels was achieved with the solution that incorporated the drive unit model in the control system structure; in this case, quality indexes $e_{\omega max}$ and $E_{\omega }$ recorded the lowest values for solution C. In turn, the highest values in terms of these quality indexes occurred when wheel-slip compensation was added to the control system, that is in the case of solutions B and D;
• It should be noted that quality indexes associated with the control of the wheels, that is $e_{\omega max}$ and $E_{\omega }$, are less significant compared to those related to the robot’s actual motion velocities, in particular $e_{vRmax}$, $E_{vR}$, $e_{\omega 0max}$, and $E_{\omega 0}$. This is due to the fact that, in wheel-slip conditions, perfect control of wheel rotation does not guarantee perfect control of the robot’s motion;
• The highest accuracy in regard to controlling the robot’s longitudinal velocity was achieved with solution D, where the control system incorporated both wheel-slip compensation and the drive unit model, as evidenced by the quality indexes $e_{vRmax}$ and $E_{vR}$. A slightly lower performance in terms of longitudinal velocity control accuracy was observed for solution C, which utilized the drive unit model within the control system;
• The highest accuracy in regard to controlling the angular velocity during the robot’s turning was achieved with the solutions that included wheel-slip compensation, as indicated by the quality indexes $e_{\omega 0max}$ and $E_{\omega 0}$, which resulted in the lowest values for solutions B and D;
• The largest deviation from the desired motion path of point $R$ in the robot occurred for solutions A and C, as evidenced by the highest values in terms of the quality indexes, $e_{pmax}$ and $E_p$. This was primarily associated with the wheel slip that occurred during turning. In turn, the smallest deviation from the desired motion path occurred for solutions B and D, thanks to the use of effective wheel-slip compensation when the robot was turning;
• The highest energy efficiency in regard to the robot’s movement was achieved for solution A, with the typical control system, and for solution C, which included the drive unit model, as evidenced by the values of the quality indexes, $p_{max}$ and $E_d$. On the contrary, the lowest energy efficiency was observed for solutions B and D, which incorporated wheel-slip compensation;
• Taking into account the quality indexes related to the accuracy of the tracking control, the best performance was achieved by solutions B and D, i.e., the control systems with wheel-slip compensation. These solutions achieved first or second place in the ranking of the quality indexes related to the accuracy of executing the desired motion. At the same time, these solutions are the worst in terms of energy efficiency, so they should be used only in a situation where energy efficiency is less important than the accuracy of the tracking control.

Figure 5 illustrates the aggregated results for the motion path in terms of the robot’s characteristic point $R$ in regard to the stationary coordinate system $\left\{O\right\}$.

The graph depicts both the desired motion path, marked as ${}_{ }{}^O{\mathit{r}}_{Rd}$, and the motion paths achieved for each of the evaluated solutions in terms of the control system, denoted as ${}_{ }{}^O{\mathit{r}}_R \left(j\right)$, where $j=\{\mathrm{A}, \mathrm{B}, \mathrm{C}, \mathrm{D}\}$ refers to the solution being analyzed.

Analyzing the motion paths of point $R$ in the robot illustrated in Figure 5, it can be observed that the best result was obtained by solution B in terms of the control system, with wheel-slip compensation, as the resulting motion path is the closest to the desired one. Slightly worse results were achieved for solution D, which also included wheel-slip compensation, but with the addition of the drive unit model. This leads to the conclusion that the most decisive factor in achieving the desired motion path is primarily wheel-slip compensation. The worst results in terms of the motion path were obtained for solutions A and C, where no wheel-slip compensation was applied. Additionally, it can be seen that the inclusion of the drive unit model does not contribute to the improvement of the motion path. Therefore, the benefit of applying it within the control system is more evident from the perspective of improving the robot’s motion parameters over time, as indicated by the previously discussed quality indexes.

The next part of this article presents the time histories in terms of the selected physical quantities for the four analyzed control system solutions.

5.5.1. Typical Linear Controller (Solution A)

Figure 6 presents the simulation results obtained for solution A in terms of the control system, that is, without wheel-slip compensation and without the use of the drive unit model. The following aspects are shown:

(a)

The tracking errors in regard to the angular velocities of the robot’s wheels ${\mathrm{e}}_{\omega i}={\omega }_{id}-{\omega }_i$, for two selected wheels ($i=\left\{\mathrm{3,4}\right\}$), i.e., the left and right ones on the robot’s rear axle (the time histories of the errors for the other wheels are similar, as the corresponding left and right ones have the same desired angular velocities, namely ${\omega }_{ld}$ and ${\omega }_{rd}$, respectively);

(b)

The errors in regard to the generalized velocities of the robot, i.e., the longitudinal velocity $e_{vR}={}_{ }{}^R{v}_{Rxd}-{}_{ }{}^R{v}_{Rx}$ and the angular velocity of turning $e_{\omega 0}={}_{ }{}^R{\omega }_{0zd}-{}_{ }{}^R{\omega }_{0z}$;

(c)

The control signals for the robot’s drives $u_i$, for two selected wheels in the robot;

(d)

The driving torques ${\tau }_i$ for all the wheels;

(e)

The normal components of the ground reaction forces ${{}_{ }{}^{R}F}_{Aiz}$ acting on the wheels;

(f)

The longitudinal slip velocities ${}^R{v}_{Six}$ for the wheels.

Analyzing the obtained time histories shown in Figure 6a,b, it can be observed that the highest tracking errors occur at moments of change in the robot’s linear velocity (when the highest acceleration is present during acceleration or deceleration) and in the angular velocity (when the direction of the robot’s motion changes, that is during the transition from longitudinal motion to turning, or vice versa). Moreover, the angular velocity error during the turning of the robot remains consistently high throughout the entire maneuver. The time histories of the control signals for the drives, $u_i$, presented in Figure 6c, are similar to the time histories of the angular velocities of the wheels, ${\omega }_i$ (in steady motion $u_i=k_e{n}_d{\omega }_{id}$). Similarly, the time histories of the driving torque ${\tau }_i$ (illustrated in Figure 6d) are similar to the time histories of the current consumption by the drives $I_i$, which results from the dependence $I_i={\tau }_i/({\eta }_d{n}_d{k}_m)$. Regarding the time histories of the driving torque ${\tau }_i$ (see Figure 6d), the highest values occurred at moments of maximum acceleration and also remained high during the turning of the robot.

One can notice differences in the values of the driving torque ${\tau }_i$ for the front and rear wheels during the starting of the movement and turning of the robot. This is due to the differences in the values of the normal components of the ground reaction forces ${{}_{ }{}^{R}F}_{Aiz}$ for the individual wheels (see Figure 6e). These differences are related to the position of the robot’s center of mass, which is shifted backwards, and result from the accelerations occurring when movement starts and during the turning of the robot. Regarding the longitudinal slip velocities of the wheels ${}^R{v}_{Six}$, illustrated in Figure 6f, they become apparent during acceleration and braking and, most notably, during turning, when they reach significantly higher values.

5.5.2. Controller with Wheel-Slip Compensation (Solution B)

Figure 7 presents the analogous results to those shown previously, but for the control system with wheel-slip compensation and also without the drive unit model (solution B).

Regarding wheel velocity control, there are no major differences when compared to the previously discussed control system solution (see Figure 6a and Figure 7a).

In this case, slightly higher values in terms of these errors occur when the robot changes its direction of motion. However, a noticeable difference can be observed in the robot’s generalized velocity errors (compare Figure 6b with Figure 7b), which now achieve significantly lower maximum values, particularly during steady turning. It is important to emphasize that these errors are much more critical in terms of the robot’s motion control in wheel-slip conditions and play a key role in determining how accurately the robot’s desired trajectory is executed. This is clearly evident when considering the previously discussed quality indexes for generalized velocities (see Table 5) and the illustration of the desired and obtained motion paths in terms of point $R$ in the robot (see Figure 5). When it comes to the driving torque, small differences are noticeable, especially regarding the slightly higher maximum values (compare Figure 6d with Figure 7d). The difference in energy efficiency is more noticeable when comparing the quality indexes for electric power and electric energy for these two cases (see Table 5). As previously noted, solution B in terms of the control system has lower energy efficiency compared to solution A.

For this and subsequent control system solutions, graphs for the normal components of the ground reaction forces ${{}_{ }{}^{R}F}_{Aiz}$ and the longitudinal slip velocities of the wheels ${}^R{v}_{Six}$ are omitted. This is due to the fact that their presentation in regard to the example involving solution A was only intended to explain the differences in the values of the driving torque ${\tau }_i$ for the front and rear wheels of the robot.

To summarize, the application of wheel-slip compensation results in a significant improvement in the accuracy of the tracking control at the cost of a slight decrease in energy efficiency.

5.5.3. Controller with Drive Unit Model (Solution C)

Figure 8 shows the research results for the robot’s control system without wheel-slip compensation and with the drive unit model (solution C).

When it comes to the control of the wheel’s angular velocity, a comparison with the previous results shows that incorporating the drive unit model into the control system significantly reduces tracking errors (compare Figure 8a with Figure 6a and Figure 7a). In this case, the maximum values of these errors are 2–4 times smaller compared to the previously discussed solutions. Due to wheel slip occurring during the robot’s motion, this solution ultimately does not lead to a significant reduction in the generalized velocity errors in terms of the robot compared to the control system without the drive unit model (compare Figure 8b with Figure 6b).

However, the results are comparable to the case where wheel-slip compensation was applied in terms of the linear velocity, but are worse regarding the angular velocity during turning (compare Figure 8b with Figure 7b), especially during steady turning. The time histories of driving torques do not differ significantly compared to the previously discussed versions of the control system (compare Figure 8d with Figure 6d and Figure 7d).

To conclude, the greatest benefit of incorporating the drive unit model into the control system lies in improving the control of the angular velocities of the wheels. Therefore, this solution may be advantageous when measurements of the robot’s actual motion parameters are unavailable and when the control system cannot compensate for the existing wheel slip.

5.5.4. Controller with Wheel-Slip Compensation and Drive Unit Model (Solution D)

Figure 9 presents the results of the simulation studies for the last analyzed case, i.e., the robot control system with wheel-slip compensation and the drive unit model (solution D).

Analyzing the tracking errors in regard to the angular velocities of the robot’s wheels, it can be observed that they are apparently higher than in the previous case (compare Figure 8a and Figure 9a). However, they do not play a decisive role in the accuracy of the robot’s motion execution when wheel-slip compensation is applied. The discussed solution allows for a reduction in generalized velocity errors of the robot in critical motion phases (compare Figure 6b, Figure 7b, Figure 8b and Figure 9b), thus proving to be effective in regard to this aspect, especially in terms of longitudinal velocity. The maximum values of the driving torque are higher compared to the previously analyzed cases (compare Figure 6d, Figure 7d, Figure 8d and Figure 9d). This results in increased electric energy consumption by the robot’s drives during the analyzed maneuver, as reflected by the highest values in terms of the quality indexes connected with energy efficiency (see Table 5).

In summary, incorporating both wheel-slip compensation and the drive unit model into the control system enables high accuracy in executing the desired motion to be achieved in terms of generalized velocities, although at the expense of a decrease in energy efficiency.

6. Conclusions and Directions for Future Works

The subject of this paper was the analysis of the influence of wheel-slip compensation and the use of a drive unit model within the control system structure on the accuracy of the trajectory tracking control and energy efficiency of a four-wheeled skid-steered robot.

The conducted research enabled the following answers to be obtained in response to the key research questions:

• Wheel-slip compensation by the control system increases the accuracy of the tracking control in terms of the robot’s generalized velocities, especially during turning. A more accurate execution of the desired motion path for the selected point of the robot is also noticeable;
• Incorporating the drive unit model into the control system improves the accuracy of the tracking control in regard to the angular velocity of the robot’s wheels. This approach, however, does not increase the accuracy of the execution of the desired motion path. High accuracy in terms of the realization of the robot’s generalized velocities is achieved when the control system incorporates both wheel-slip compensation and the drive unit model;
• Increasing the accuracy of motion results in decreased energy efficiency. This correlation occurs mainly in the case of wheel-slip compensation. It should be emphasized, however, that the improvement in tracking control accuracy is greater than the decrease in energy efficiency. The control system that simultaneously incorporates wheel-slip compensation and the drive unit model demonstrates the lowest energy efficiency.

The conducted research indicates that each of the analyzed solutions in terms of the control system can be advantageous in specific situations, depending on whether motion accuracy or energy efficiency is more critical, as well as what measurement system the robot is equipped with.

It should be emphasized, however, that the object of the research was a four-wheeled skid-steered robot and that the presented conclusions refer only to such a kinematic structure. In the case of robots with a different kinematic structure, analogous studies should be carried out and different conclusions may be reached.

In the presented simulation studies, a desired motion trajectory was adopted that did not lead to large values in terms of wheel slip during the robot’s acceleration and braking, a trajectory which was also feasible for the robot’s drive units. If acceleration values or motion parameters were adopted that were too high, which would not be feasible due to the limitations of the drive units (electric power, current, rotational speed), the analyzed solutions would achieve different values in terms of the quality indexes. In such a case, taking into account the drive unit model, for example, may play a more important role.

Directions for further research work may cover the following aspects:

• Conducting empirical research using a real robot for particular solutions in terms of the control system analyzed in this article and performing a comparative analysis of the results;
• Carrying out simulations or empirical research on scenarios where higher accelerations lead to large wheel slip, or when the desired trajectory cannot be fully achieved due to the limitations of the robot’s drive units;
• Applying the discussed solutions within more advanced control systems, including model-based systems, such as robust, adaptive, or predictive control systems;
• Consideration of the occurrence of parametric disturbances, related to changes in the type of ground and other factors influencing the robot’s movement, such as varying load mass.