Development of Online Adaptive Traction Control for Electric Robotic Tractors

Abstract

Estimation and control of wheel slip is a critical consideration in preventing loss of traction, minimizing power consumptions, and reducing soil disturbance. An approach to wheel slip estimation and control, which is robust to sensor noises and modeling imperfection, has been investigated in this study. The proposed method uses a simplified form of wheels longitudinal dynamic and the measurement of wheel and vehicle speeds to estimate and control the optimum slip. The longitudinal wheel forces were estimated using a robust sliding mode observer. A straightforward and simple interpolation method, which involves the use of Burckhardt tire model, instantaneous values of wheel slip, and the estimate of longitudinal force, was used to determine the optimum slip ratio that guarantees maximum friction coefficient between the wheel and the road surface. An integral sliding mode control strategy was also developed to force the wheel slip to track the desired optimum value. The algorithm was tested in Matlab/Simulink environment and later implemented on an autonomous electric vehicle test platform developed by the Nanjing agricultural university. Results from simulation and field tests on surfaces with different friction coefficients (μ) have proved that the algorithm can detect an abrupt change in terrain friction coefficient; it can also estimate and track the optimum slip. More so, the result has shown that the algorithm is robust to bounded variations on the weight on the wheels and rolling resistance. During simulation and field test, the system reduced the slip from non-optimal values of about 0.8 to optimal values of less than 0.2. The algorithm achieved a reduction in slip ratio by reducing the torque delivery to the wheel, which invariably leads to a reduction in wheel velocity.

1. Introduction

In precision agriculture, the primary objectives are to improve production efficiency; reduce environmental degradation; and, more importantly, optimize the use of available resources such as energy, water, seeds, and fertilizer. In line with these objectives, electric tractors (ET) were introduced into agricultural productions. In ET, electric motors are used for traction, while chemical batteries, fuel cells, or ultracapacitors are used as energy sources. Compared to internal combustion engines, electric tractors have the advantages of having fewer components, being emission-free, having high efficiency, being independent from petroleum, and having quiet and smooth operations [1]. ET is also easy to handle because its dynamic can be easily realized by controlling the electric motor. ET’s advantages are realized because electric motors, which are its prime components, have a quick torque response and can be easily controlled by manipulating voltage [2]. To leverage on the advantages of electric vehicle technologies, electric tractors were developed for pesticide applications and other agricultural purposes. The main advantages are that such operation does not require high power; as such, their power demand can be met by the electric batteries [3,4,5,6,7,8]. Though much progress has been achieved on improving batteries’ performance, its energy density is still lower than that of conventional fuel, and the cost is still high [9,10,11,12,13]. Studies have been conducted with the sole aim of increasing the operating range of electric vehicles, and these studies have led to the development of various techniques for reducing power consumption. Online traction control and optimization is among the most promising researched areas.

A critical role is played by traction control in vehicles’ dynamic control. It improves drive efficiency, increases safety, and makes the system more stable [14]. Traction is the force that propels the vehicles in the direction of motion; it is generated by the friction between the wheels and the road/terrain surface. Friction is dependent on many variables, including road type, road conditions, wheel types, wheel slip, vertical load, tire wear state, and temperature [14,15]. Additionally, it is affected by operation parameters such as velocity and acceleration. Researchers in vehicles’ dynamic controls have lumped all these factors affecting the friction into a single parameter call friction coefficient. This parameter is found to be dependent on the wheel slip. There exists an optimum value of slip, for each road type and set of conditions, at which the friction coefficient is high. Since the slip affects friction coefficient and invariably the traction force, it is used as the control variable in most of the traction control. Therefore, traction controls are designed to regulate slip to the optimum value. By so doing, degradation of driving performances resulting during sudden acceleration in poor road conditions is eradicated. One more benefit is that the power consumption is also reduced. Unfortunately, traction control is generally difficult because of the uncertainty and nonlinearity associated with the friction between the wheel and road surface [16].

The earliest studies on traction control were based on maximum transmissible torque estimation (MTTE), model following control (MFC), and velocity feedback. However, the most recent studies were based on slip feedback. This change in focus is connected to the fact that there is a strong relationship between slip rate and traction force, and the ease of slip estimation due to the development and availability of improved sensors [17,18,19,20,21,22].

Previous studies were mainly directed to address the problems of estimating the velocity of the vehicle. Difficulties, such as drift caused by DC offset while using accelerometers, high cost of optical sensors, and their susceptibility to contamination by mud or water necessitated the development of traction control that does not depend on vehicle velocity estimation [23,24,25,26,27]. Model following control (MFC) does not need information on vehicle velocity or acceleration; it is based on only wheel torque and wheel velocity. However, to accommodate worst-case scenario in a real field environment, high control gain is required in tuning the controllers, and this constrains the performance of the method [2,28,29,30]. MTTE is another means of traction control without the need to estimate vehicle speed. It is an improved form of MFC, and it centers on comparing the maximum torque developed based on the vehicle dynamic model and the one measured using sensors. When wheels slip, a difference is obtained between the two torques, and this difference is feedback to the controller. The method’s major drawbacks are high sensitivity to variation in wheel inertia and steering difficulty in situations where the wheels have different friction coefficients [29,31,32].

Traction control that is more robust to changes in terrain types, the unknown wheel-terrain conditions, and does not have the limitations of MFC and MTTE was achieved through slip feedback control. However, it is only applicable where a reliable estimate of vehicle velocity is available. Because it was found experimentally that the traction is highest at the slip range 15–30% for many terrains, earlier studies on slip feedback pegged the optimum slip within this range [33]. However, it was found that the optimum value of slip is dependent on tires, operation parameters, and road conditions. These have led to the development of slip controls based on a real-time estimate of the optimum slip, which correspond to high friction between tire and road surfaces. Lookup tables were used to estimate optimum slip, and gave a satisfactory performance on the various roads. Because lookup-table-based systems are detuned to consider the worst scenario, such as wet or icy roads, a suboptimal controller is usually obtained [16].

The most recent research trend is the online estimation of optimum slip for different terrains in various conditions. In [33], the performances of fuzzy logic and sliding mode observer were tested for optimum slip estimation. The study established that observer techniques are appropriate for online estimation of peak adhesion coefficients for different tire-road surface conditions. In [34], the optimal slip value was obtained by taking the first derivative of the slip–adhesion function, the point at which the derivative is equal to zero is taken as the point of optimum slip. The main problem of the method is that it gives an average value of the optimum slip. Furthermore, studies have shown that it mainly estimates high slip changes such as the change from dry road to gravel [35]. Additionally, machine learning applications such as the support vector machine (SVM) and artificial neural network (ANN) have been applied in estimating optimum slip, especially in a situation of low changes in slip value (same road and condition). However, the problem of a tradeoff between generality and estimation speed remains [36].

Due to its ability to perform good unbiased estimation, the recursive least square algorithm (RLS) has been employed to estimate road friction conditions [37]. The algorithm also has a fast response speed, strong anti-jamming capability, and can deal with nonlinearities associated with complex tire models [38,39].

This work’s main objective is to develop an algorithm for the online traction control of electric vehicles used in agricultural fields. The original contribution lies in combining the robustness of sliding mode controller and observer, with the simplicity of the Burckhardt model, to accurately estimate the optimum maximum friction coefficient and the slip ratio at which it occurs. The use of sliding mode observers to estimate the longitudinal wheel force has eliminated the use of expensive sensors in measuring the force. In essence, the study includes developing an optimum slip estimation algorithm by using the Burckhardt traction model, and the design of integral sliding mode controller for tracking the optimum slip. The systems’ performance was evaluated in the Matlab/Simulink environment. Finally, the system was implemented on an electric vehicle test platform and tested in different terrain types.

2. Materials and Methods

2.1. Vehicle Modelling

Four degrees of freedom vehicle model was established by considering the motion along the x, y, and z-axis, also considered is the yaw moment along the z-axis. As depicted in Figure 1, the vehicle has a rear-wheel steering system. The x and y axes reside in the ground plane, while the z-axis is perpendicular to the plane of x-y.

Because low speed of operations and flat terrain are considered, the effect of aerodynamic force and grade resistance were neglected to simplify the model. The machine is designed to work in unstructured terrain where soil under the wheels may not be the same. The following are the equations necessary for development and analysis of traction control.

Longitudinal force on the vehicle:

$$M\left(a_x-\dot{\gamma }\ast V_y\right)=F_{xfr}+F_{xfl}+F_{xrr}+F_{xrl}-(F_{roll, fr}+F_{roll, fl}+F_{roll, rr}\cos \delta +F_{roll,rl}\cos \delta )$$

(1)

Lateral force on the vehicle:

$$M\left(a_y-\dot{\gamma }\ast V_x\right)=F_{yfr}+F_{yfl}+F_{yrr}+F_{yrl}-(F_{roll,rr}\sin \delta +F_{roll,rl}\sin \delta )$$

(2)

Yaw moment of the vehicle:

$$I_z=\left(F_{xfr}+F_{xrr}\right)w_r-\left(F_{xfl}+F_{xrl}\right)w_f+\left(F_{yrr}+F_{yrl}\right)l_r-\left(F_{yfr}+F_{yfl}\right)l_f$$

(3)

where $M$ denotes the mass of the vehicle at center of gravity (CG); $a_x$ denotes the longitudinal acceleration of the vehicle; $a_y$ denotes the lateral acceleration of the vehicle; $V_x$ denotes the longitudinal velocity of the vehicle; $V_y$ denotes the lateral velocity of the vehicle; $\dot{\gamma }$ denotes the yaw rate of the vehicle; and $\delta$ denotes the steering angle.

$F_{xi}$ stands for the longitudinal tire-road friction force for the wheels; $F_{yi}$ stands for the lateral tire-road friction force for the wheels; $F_{roll,i}$ stands for the rolling resistance force on wheels direction of travel; the subscript, $i$, in $F_{xi}$, $F_{yi}$, and $F_{roll,i}$ stands for $fl,fr,rl,$ and $rr$, respectively, with $fl,$ as front left; $fr,$ front right; $rl,$ rear left, and $rr$, rear right, respectively.

$l_f$ stands for distance from the vehicle center of gravity to axis that passed through the center of the front wheels; $l_r$, distance from the vehicle center of gravity to the axis that passed through the center of the rear wheels; $W_f$ and $W_r,$ front and rear wheelbase (equal length); $H_{CG}$, height of the center of gravity; and $CG$, the center of gravity.

Vertical forces acting on each wheel are defined as follows [40,41]:

$$\left[\begin{array}{c}{F}_{zfr}\\ \\ F_{zfl}\\ \\ F_{zrf}\\ \\ F_{zrr}\end{array}\right]=M·\left[\begin{array}{c}\left(\frac{l_f\times g+H_{cg}\times a_x}{L}\right)\left(\frac{1}{2}-\frac{H_{cg}}{W_f}\times \frac{a_y}{g}\right)\\ \left(\frac{l_f\times g+H_{cg}\times a_x}{L}\right)\left(\frac{1}{2}+\frac{H_{cg}}{W_f}\times \frac{a_y}{g}\right)\\ \left(\frac{l_r\times g+H_{cg}\times a_x}{L}\right)\left(\frac{1}{2}-\frac{H_{cg}}{W_r}\times \frac{a_y}{g}\right)\\ \left(\frac{l_r\times g+H_{cg}\times a_x}{L}\right)\left(\frac{1}{2}+\frac{H_{cg}}{W_r}\times \frac{a_y}{g}\right)\end{array}\right]$$

(4)

where $F_{zfr}$, $F_{zfl}$, $F_{zrl}$, and $F_{zrr}$ are the vertical force on the front right, front left, rear left, and rear right wheel, respectively; $g$, is acceleration due to gravity; and $L=l_f+l_r$.

The dynamic of each wheel model can be described by the torque balance of the system in Figure 2, as shown in Equation (5):

$${\dot{\omega }}_i=\frac{T_{mi}}{J}-\frac{r·F_{xi}}{J}$$

(5)

where $J$ stands for the wheel inertia including motor and reducer; $T_{mi},$ the torque transmitted to the wheels; ${\dot{\omega }}_i$, wheels rotational acceleration; and $r$, wheel radius.

Wheels longitudinal slip is defined in Equation (6). To avoid a situation whereby the denominator of the slip function becomes zero, which will lead to an infinite value for the slip, a small constant value of 0.05 m/s was set as the minimum value for the slip function denominator.

$${\lambda }_i=\frac{{\omega }_i·r-V_x}{{\omega }_i·r}$$

(6)

where ${\lambda }_i$ stands for wheel longitudinal slip and $\omega$ stands for wheel rotational speed.

The road friction coefficient ${\mu }_i\left(\lambda \right)$, which is a nonlinear function of many physical variables including the velocity of the vehicle, the road surface conditions, and slip ratio, is defined as

$${\mu }_i\left(\lambda \right)=\frac{F_{xi}}{F_{zi}}$$

(7)

2.2. Slip Dynamic Based on a Simplified Vehicle Model

The slip dynamic, which is the basis for formulating traction control algorithm, is developed in this section. In order to reduce computational demand on controller, the following simplifications were made in the slip dynamics: it is assumed that viscous torque developed at the wheels’ axles is negligible compared to the torque developed by the longitudinal force, and the yaw rate and lateral acceleration were also set to zero. Due to weight transfer, the tractive force that develops at each of the wheels is assumed to be different. Furthermore, it is assumed that the wheels longitudinal velocity ($V_{wi}={\omega }_i·r$), and the vehicle longitudinal velocity are greater than or equal to zero ($V_{wi}\ge 0$ and $V_x\ge 0$). To control traction or speed of the vehicle, the torque delivered to the wheels by the electric motor is the parameter to be manipulated; therefore, the torque is the control input. Because the structures of the traction models of each of the wheels is similar, method for controlling the tractive force of only one wheel is going to be described. These simplifications may cause some uncertainties; however, real systems are far from being idealistic, uncertainties and unmodelled dynamics are inevitable, and robust controllers are usually used to take charge of these systems’ limitations [42].

The simplified longitudinal acceleration model based on one-wheel dynamics, derived from Equation (1), is given.

$$a_x=\frac{{{\displaystyle \sum }}^{ }{F}_{xi}-{{\displaystyle \sum }}^{ }{F}_{roll}}{M}$$

(8)

From Equation (6), the wheels rotational speed as a function of slip ratio and vehicle speed can be written as

$${\omega }_i=\frac{V_x}{r}\left(\frac{1}{1-{\lambda }_i}\right)$$

(9)

When the slip function in Equation (6) is differentiated, the following function is obtained:

$${\dot{\lambda }}_i=\frac{V_x}{{\omega }_i^2·r}{\dot{\omega }}_i-\frac{1}{{\omega }_i·r}{a}_x$$

(10)

While the state variables in Equation (10) are $\dot{\omega }$, the wheel rotational acceleration and $a_x$, the vehicle longitudinal acceleration, the variable of interest in the design of the slip controller is the slip ratio. Therefore, the slip dynamic model in (10) was manipulated to make the slip ratio and vehicle velocities the state variables.

By substituting Equations (5) and (9) into Equation (10), the following model for slip dynamic function is obtained:

$$\dot{{\lambda }_i}=-\frac{1}{V_x}\left(\frac{1-{\lambda }_i}{M}\left({{\displaystyle \sum }}^{ }{F}_{xi}-{{\displaystyle \sum }}^{ }{F}_{roll}\right)+\frac{r^2{\left(1-{\lambda }_i\right)}^2}{J}·F_{xi}\right)+\frac{r{\left(1-{\lambda }_i\right)}^2}{J·V_x}{T}_{mi}$$

(11)

For the sake of simplicity purpose, we take

$$f{f}_i=-\frac{1}{V_x}\left(\frac{1-{\lambda }_i}{M}\left({{\displaystyle \sum }}^{ }{F}_{xi}-{{\displaystyle \sum }}^{ }{F}_{roll,i}\right)+\frac{r^2{\left(1-{\lambda }_i\right)}^2}{J}·F_{xi}\right)$$

(12)

And

$$g{g}_i=\frac{r·{\left(1-{\lambda }_i\right)}^2}{J·V_x}$$

(13)

Therefore, the slip dynamic in Equation (11) can be written as

$$\dot{{\lambda }_i}=f{f}_i+g{g}_i·T_{mi}$$

(14)

2.3. Wheel Longitudinal Force Observer

In this section, an observer is proposed to identify the longitudinal force of each wheel. Due to the high cost of commercially available sensors for online measurement of tire longitudinal force, $F_{xi},$ efficient technique of estimation for these forces is strongly recommended. The fast dynamics of wheels’ vertical forces make the longitudinal force, which is dependent on the wheel’s vertical force, unsuitable for adaptive compensation. Additionally, the large uncertainty in the friction coefficient between the terrain and the wheel makes it unsuitable for robust compensation because of the need for high loop gain [43]. Therefore, an observer was designed to estimate the longitudinal force. Due to finite convergence; insensitivity to modeling errors; modeling uncertainties; and measurement noise, a sliding mode observer was chosen [44,45,46,47,48,49,50]. Values of vehicle mass and wheel inertia, accurate measurement of wheel torque, and speed were used in the observer models.

Design of the Observer

Because of the lack of cost-effective sensors for measuring longitudinal force, observers are used to estimate the force from the measurement of other variables. The use of the observer also solves the problem of noise in the measuring sensors. An observer was designed to estimate the friction force between the tire and the soil using wheel angular velocity, torque, and wheel dynamic model. In the observer’s design and implementation, numerical differentiation to obtain angular acceleration was avoided because of its sensitivity to measurement noise and quantization, as reported in [43]. Under the assumption that wheel inertia ($J$) is known, first-order sliding mode observer was designed as follows.

First, based on the procedure described in [3,43,51], the wheel torque balance in Equation (5) is estimated as

$$J·{\dot{\widehat{\omega }}}_i=T_{mi}-r·L\left({\omega }_i,{\widehat{\omega }}_i\right)$$

(15)

where ${\dot{\widehat{\omega }}}_i$ is the state estimate of the ‘wheel’s angular acceleration; ${\widehat{\omega }}_i$ is the state estimate of the ‘wheel’s angular velocity; and $L\left({\omega }_i,{\widehat{\omega }}_i\right)$ is a discontinuous input signal for the observer.

The discontinuous input signal for the observer in Equation (15) is defined as

$$L\left({\omega }_i,{\widehat{\omega }}_i\right)={\Omega }_i·sign\left(s_{\omega i}\right)$$

(16)

$$s_{\omega i}={\omega }_i-{\widehat{\omega }}_i$$

(17)

where ${\Omega }_{\mathrm{i}}$ is the observer gain, and $s_{\omega i}$ is the the sliding variable, which is the difference between the measured and estimated wheel rotational speed.

The discontinuous signum function in Equation (16) is defined as

$$sign\left(s_{\omega i}\right)=\{\begin{array}{c}+1, s_{\omega i}\le 0\\ -1, s_{\omega i}>0\end{array}$$

(18)

By differentiating Equation (17), and substituting Equations (5), (15) and (16), the expression in Equation (19) is obtained

$${\dot{s}}_{\omega i}={\dot{\omega }}_i-{\dot{\widehat{\omega }}}_i=\frac{-r}{J}\left(F_{xi}+\Omega ·sign\left(s_{\omega i}\right)\right)$$

(19)

As the dynamic function in Equation (19) approaches zero, the discontinuous input signal $L\left({\omega }_i,{\widehat{\omega }}_i\right)$ will approach $F_{xi}$, thereby getting an estimate of $F_{xi}$. The function $L\left({\omega }_i,{\widehat{\omega }}_i\right)$ can only take the values of $\pm \Omega$ and zero, so it determines the sign of $F_{xi}.$

The condition for the sliding variable $s_{\omega i}$ to converge to zero in finite time was obtained by formulating the Lyapunov function, $V_1$, [52]. The function is defined as the positive definite function of the estimation error, $s_{\omega i}$:

$$V_1=\frac{1}{2}{s}_{\omega i}^2$$

(20)

The sliding condition that guarantees the attractivity of $s_{\omega i}$ is obtained by differentiating Equation (20) as

$${\dot{V}}_1=s_{\omega i}{\dot{s}}_{\omega i}$$

(21)

By substituting the error dynamic in Equation (19) into Equation (21), the following expression is obtained:

$${\dot{V}}_1=s_{\omega i}·\left(-\frac{r}{J}\left(F_{xi}+{\Omega }_{\mathrm{i}}·sign\left(s_{\omega i}\right)\right)\right)$$

(22)

This implies that

$${\dot{V}}_1=-\frac{r}{J}\left(F_{xi}·s_{\omega i}+{\Omega }_i·\left|s_{\omega i}\right|\right)$$

(23)

Because $s_{\omega i}$ is $\pm$ and the $\left|s_{\omega i}\right|$ is always $+1,$ if $\Omega$ is chosen such that the following inequality is satisfied

$${\Omega }_i>\max \left(\left|F_{xi}\right|\right)$$

(24)

It will give

$${\dot{V}}_1<0$$

(25)

Satisfying the conditions in Equation (20) and the inequality in (24) proves that the sliding variable is asymptotically stable.

Theoretically, the input signal of the sliding mode observer switches at an infinitely high frequency. However, in reality, it is impossible to have such an ideal input signal, so the system oscillates at a limited frequency, with both high- and low-frequency components [53]. Therefore, an equivalent value of the frictional force ($F_{xi}$) was obtained from the low-frequency component of ${\dot{S}}_{wi}$ by using lowpass filter. So, only the convergence of the low-frequency component of discontinuous signal was required to estimate the longitudinal force [43,51,54].

2.4. Wheel Terrain Interaction and Optimum Slip Estimation

It has been found that the value of the frictional force between the wheel and the terrain varied with wheel slip, and a particular slip at which the frictional force is high exists for all the terrains. As a result of this, a lot of traction prediction models that relate to wheel slip, wheel longitudinal force, and other wheel-terrain parameters have been developed. Some of the developed models are purely theoretical, while others are empirical or semi-empirical. The semi-empirical models are computationally faster than the theoretical ones. They are more general in terms of applicability in comparison to purely empirical ones [8]. Because the tractor is designed to work on concrete, road, firm soil, and agricultural farms, a model with broad applicability as well as simplicity is required. Due to its relatively simple structure, high accuracy in fitting traction data, and broad applicability in many driving conditions, the Burckhardt tire model was used to characterize this work’s traction force [55,56]. Compared to the magic formula and other models, which are more accurate, the Burckhardt tire model provided a good tradeoff between modeling accuracy and computational cost [50].

The model equation as a function of lambda ($\lambda$) is described as

$$\mu {\left(\lambda \right)}_i=c_1\left(1-e^{-c_2{\lambda }_i}\right)-c_3{\lambda }_i$$

(26)

where the parameters $c_1, c_2, \mathrm{and} c_3$ are constant variables that are different for different terrain. The values for the parameters for some terrains, including high and low friction coefficient surfaces such as dry concrete and snow, are presented in

Table 1. Parameters of road surface adhesion coefficient [57].

Road Surface	${\mathit{c}}_1$	${\mathit{c}}_2$	${\mathit{c}}_3$
Dry asphalt	1.28	23.99	0.52
Wet asphalt	0.86	33.82	0.35
Wet pebbles	0.40	33.71	0.12
Snow	0.19	94.13	0.06
Ice	0.05	306.4	0.001

.

The main problem is that the values of the model parameters are not available for all the diverse groups of terrains in the real world, and the available ones are not fixed; even for a particular terrain, they change with the change in surface condition. Moreover, as shown in Figure 3, the friction coefficient changes with the change in the slip ratio, and it reaches a peak point at a particular value of slip (optimum), beyond which the value starts to decline.

For high traction, the longitudinal force needs to be controlled at that optimum value of slip. A simple and feasible interpolation method was put forward for estimating the instantaneous values of the parameters, maximum friction coefficient, and the optimum slip for different terrains in previous studies [58,59]. The principle of operation of the method is shown in Figure 4, and the process is as follows.

Step 1. 

By estimating the frictional force as in Section 2.3 and substituting the value equation (7), the instantaneous value of friction coefficient is obtained, while the instantaneous value of slip was obtained from Equation (6).

Step 2. 

Using the instantaneous value of slip (${\lambda }_i$) and the model parameters of the five different terrain types shown in Table 1, five sets of transient values for friction coefficient (${\mu }_i$) are calculated by using Equation (26).

Step 3. 

Then, five sets of error values (${\Delta }_i$) are generated, by subtracting the estimate friction coefficient $\left({\widehat{\mu }}_i\right)$ from the calculated ones (${\mu }_i$).

Step 4. 

The two smallest error values (${\mathsf{\Delta }}_1$ and ${\mathsf{\Delta }}_2$) from step 3 are selected.

Step 5. 

Using Equation (26) and the values of soil parameters ($c_1, c_2 \mathrm{and} c_3$) at the two points where ${\mathsf{\Delta }}_1$ and ${\mathsf{\Delta }}_2$ are obtained, two functions are generated:

$${\mu }_{x1}{\left(\lambda \right)}_i=c_1\left(1-e^{-c_2{\lambda }_i}\right)-c_3{\lambda }_i$$

(27)

$${\mu }_{x2}{\left(\lambda \right)}_i=c_1\left(1-e^{-c_2{\lambda }_i}\right)-c_3{\lambda }_i$$

(28)

Step 6. 

Using ${\mathsf{\Delta }}_1$ and ${\mathsf{\Delta }}_2$ from step 4 and the functions ${\mu }_{x1}{\left(\lambda \right)}_i$ and ${\mu }_{x2}{\left(\lambda \right)}_i$ from step 5, a new function is reconstructed to estimate the maximum friction coefficient and the corresponding slip values at which it occurs (Figure 4 and Equation (29)).

$$\mu {\left(\lambda \right)}_i=\frac{{\Delta }_1}{{\Delta }_2+{\Delta }_1}{\mu }_{x1}{\left(\lambda \right)}_i+\frac{{\Delta }_2}{{\Delta }_2+{\Delta }_1}{\mu }_{x2}{\left(\lambda \right)}_i$$

(29)

Step 7. 

By differentiating Equation (29) and setting the value to zero, the optimum slip (${\lambda }_{opt,i}$), at which the friction coefficient $\mu {\left(\lambda \right)}_i$ is maximum is obtained. The maximum friction coefficient ${\mu }_{max,i}$ can be obtained as in Equation (30)

$${\mu }_{\max ,i}=f\left({\lambda }_{opt,i}\right)$$

(30)

where ${\mu }_{\max ,i}$ is the maximum tire-road friction coefficient and ${\lambda }_{opt,i}$ is the corresponding slip ratio of ${\mu }_{max,i}$.

3. Traction Control Design

The procedure of developing a controller for achieving the system’s objective of forcing the system to follow the desired optimum slip is described in this section.

3.1. Controller Design and Analysis

Due to the need of having high robustness to bounded parameter variations and the modeling imperfections that exist in wheel-terrain interaction systems, and also the issue of unmodeled dynamics due to pitch and roll moments, a sliding mode control (SMC) technique was proposed in this study. SMC has shown superior performance to other controllers such as PI in our previous study [3]. The controller was first designed without considering parameter variations, and later it was made robust to bounded parametric variations. As in all SMC design processes, the sliding surface was designed first, then a control law, that forced and retained the system at the sliding surface was designed.

3.1.1. The Sliding Surfaces Design

Designing a sliding surface for a second and higher-order systems results in exponential convergence of the state errors to zero, while in the first-order systems it resulted in a vertical line, which does not yield exponential convergence, in the phase plane. Furthermore, the dynamics outside the sliding surface become difficult to control [43]. In such cases, integral sliding mode control is used to ensure the state error’s exponential convergence, even when the system dynamic is not equal to zero. One more advantage of integral sliding mode control is that it eliminates the reaching phase completely. An auxiliary sliding surface was designed as follows [51,60]:

$$S_{aux,i}={\tilde{\lambda }}_i\left(t\right)+\mathsf{\eta }{{\displaystyle \int }}_0^t{\tilde{\lambda }}_i\left(\tau \right)d\tau -{\tilde{\lambda }}_i\left(0\right)$$

(31)

$${\tilde{\lambda }}_i={\lambda }_i-{\lambda }_i^{*}$$

(32)

where $S_{aux,i}$ represents the auxiliary sliding surface; ${\lambda }_i^{*}$ the optimum slip ratio; $\mathsf{\eta }$ a positive constant that determines the convergence rate, $\tau$ a dummy variable used for integration, and ${\tilde{\lambda }}_i\left(0\right)$ a constant value that shifts the sliding surface in the phase plane to ensure that the auxiliary sliding mode starts from the beginning without any reaching phase, for all initial conditions (that is, $S_{aux}\left(0\right)=0$).

The value ${\tilde{\lambda }}_i\left(0\right)$ was calculated as described in [43,60]; first, the sliding surface was set to zero:

$$S_{aux}\left(0\right)=0$$

(33)

Equation (31) becomes

$${\tilde{\lambda }}_i\left(t\right)+\mathsf{\eta }{{\displaystyle \int }}_0^t{\tilde{\lambda }}_i\left(\tau \right)d\tau ={\tilde{\lambda }}_i\left(0\right)$$

(34)

The dynamic of the system in (34) was solved by change of variables, such that

$$z{\left(t\right)}_i={{\displaystyle \int }}_0^t{\tilde{\lambda }}_i\left(\tau \right)d\tau$$

(35)

So, if $z{\left(t\right)}_i={{\displaystyle \int }}_0^t{\tilde{\lambda }}_i\left(\tau \right)d\tau$, then ${\tilde{\lambda }}_i\left(t\right)=\dot{z}{\left(t\right)}_i$ and ${\tilde{\lambda }}_{ }\left(0\right)=\dot{z}\left(0\right)$; hence, by substituting these values, Equation (34) becomes

$$\dot{z}{\left(t\right)}_i+\mathsf{\eta }·z{\left(t\right)}_i=\dot{z}{\left(0\right)}_i$$

(36)

As given in [60], the solution the inhomogeneous ordinary differential equation (ODE) gives in Equation (36) is then the dynamic of $z\left(t\right)$, and it is given as

$$z{\left(t\right)}_i=z{\left(0\right)}_i·e^{-\mathsf{\eta }·t}+\frac{1}{\mathsf{\eta }}·\dot{z}{\left(0\right)}_i·\left(1-e^{-\mathsf{\eta }·t}\right)$$

(37)

By expressing Equation (37) in the original variable and simplifying, the following value was obtained:

$${{\displaystyle \int }}_0^t{\tilde{\lambda }}_i\left(\tau \right)d\tau =\frac{1}{\mathsf{\eta }}{\tilde{\lambda }}_i\left(0\right)·\left(1-e^{-\mathsf{\eta }·\mathrm{t}}\right)$$

(38)

Equation (38) has shown that for any nonzero initial slip error (${\tilde{\lambda }}_i\left(0\right)$), the integrated error converges to a nonzero value. However, the variable of interest is the slip error, not the integrated error; therefore, Equation (38) was differentiated once to the slip error as follows:

$${\tilde{\lambda }}_i\left(t\right)={\tilde{\lambda }}_i\left(0\right)·e^{-\mathsf{\eta }·\mathrm{t}}$$

(39)

So, it can be seen from (39) that the error ${\tilde{\lambda }}_i$ approaches zero as time ($t$) approaches infinity, with an exponential rate of $\mathsf{\eta }$.

3.1.2. Control Law Design

By assuming that the dynamic model’s initial conditions are known, asymptotic output tracking can be achieved by developing a sliding mode control law that compensates for the bounded disturbance in the slip dynamic and drives the sliding variable to zero as time increases.

By taking the derivative of the auxiliary sliding surface in Equation (33), we get

$${\dot{S}}_{aux,i}={\dot{\tilde{\lambda }}}_i+\mathsf{\eta }·{\tilde{\lambda }}_i$$

(40)

Substituting Equation (32) into (40), and assuming that the derivative of the reference slip is zero ($\dot{{\lambda }_i^{\ast }}=0$), the following expression was obtained:

$${\dot{S}}_{aux,i}={\dot{\lambda }}_i+\mathsf{\eta }·{\tilde{\lambda }}_i$$

(41)

The slip dynamic $\dot{\lambda }$ from Equation (14) is substituted into (41), and it becomes:

$${\dot{S}}_{aux}=f{f}_i+g{g}_i·T_{mi}+\mathsf{\eta }·{\tilde{\lambda }}_i$$

(42)

So, the problem of stabilizing the slip dynamics is now reduced to that of forcing auxiliary sliding variable in Equation (42) to zero, and this was achieved by designing a stabilizing control input $T_{mi}$.

In order to design a stabilizing control input $T_{mi}$, a Lyapunov function ($V_2$) was formulated, and the control input ($T_{mi}$), was chosen such that ($\dot{V_2}$) is negative definite. The following positive definite function ($V_2$) was chosen as the candidate Lyapunov function:

$$V_2=\frac{1}{2}{S}_{aux,i}^2$$

(43)

The time derivative of $V$ is given as

$${\dot{V}}_2=S_{aux,i}·{\dot{S}}_{aux,i}$$

(44)

Substituting (42) into (44) gives

$${\dot{V}}_2=S_{aux}·\left(f{f}_i+g{g}_i·T_{mi}+\mathsf{\eta }·{\tilde{\lambda }}_i\right)$$

(45)

In order to make (45) negative definite, the control input T_mi should be chosen as

$$T_{mi}=-\frac{f{f}_i+\mathsf{\eta }·{\tilde{\lambda }}_i+{\mathsf{\Phi }}_i·sign\left(S_{aux,i}\right)}{gg}$$

(46)

With the parameter ${\mathsf{\Phi }}_i$ as positive integer, and substituting (46) in (45), the following is obtained:

$${\dot{V}}_2=-{\mathsf{\Phi }}_i·sign\left(S_{aux,i}\right)·S_{aux,i}=-{\mathsf{\Phi }}_i\left|S_{aux,i}\right|$$

(47)

Because Equation (47) is negative definite, it is proved that (43) is a Lyapunov function, and $S_{aux}=0$ is an asymptotically stable equilibrium. Therefore, the system trajectory will converge to $S_{aux}=0$ as time approaches infinity, and the sign function in Equation (46) is defined as

$$sign\left(S_{aux,i}\right)=\{\begin{array}{c}\begin{array}{c}-1, S_{aux,i}<0\\ 0, S_{aux,i}=0\end{array}\\ +1, S_{aux,i}>0\end{array}$$

(48)

3.2. Robust Controller Design

Due to the uncertainty associated with the functions $f{f}_i and g{g}_i,$ their exact value is not known, only their estimate ${\widehat{ff}}_i$ and ${\widehat{gg}}_i$ is known. As reported in [52], the estimation errors on $ff$ and $gg$ were assumed to be bounded by the following:

$$\left|\widehat{f{f}_i}-f{f}_i\right|\le F_i$$

(49)

$$0<g{g}_{min,i}\le gg\le g{g}_{max,i}$$

(50)

where $F_i$ is a known function bounding $f{f}_i$; $g{g}_{min,i}$, minimum value of the function $g{g}_i$; and $g{g}_{max,i}$, the maximum values of the function $g{g}_i$.

Because the control input ($T_{mi}$) multiplies the function $\left(g{g}_i\right)$, so the estimates of $gg$ were chosen as the geometric mean of it bound:

$${\widehat{gg}}_i=\sqrt{g{g}_{min,i}·g{g}_{max,i}}$$

(51)

The bound on the estimation error of $g{g}_i$ in Equation (53) was also expressed in another form below:

$${\beta }^{-1}\le \frac{{\widehat{gg}}_i}{g{g}_i}\le \beta$$

(52)

The parameter $\beta$ will be used later in calculating the controller gain (${\mathsf{\Phi }}_i$) and it is given as:

$$\beta =\sqrt{\frac{g{g}_{max,i}}{g{g}_{min,i}}}$$

(53)

Vertical load on the wheels, wheels radius, and rolling resistance were assumed to be of unknown exact value; this was due to the weight transfer that is caused by longitudinal acceleration, roll and pitch acceleration, changes in the wheels rolling radius, and the effect of rolling resistance. Therefore, the functions $f{f}_i$ and $g{g}_i$ were modified to consider the bounded disturbances on these vehicle parameters as follows [43,52]:

$${\widehat{ff}}_i=-\frac{1}{V_x}\left(\left(\frac{1-{\lambda }_i}{\widehat{M}}\left({{\displaystyle \sum }}^{ }{F}_{xi}-{{\displaystyle \sum }}^{ }{\widehat{F}}_{roll}\right)\right)+\left(\frac{{\widehat{r}}^2·F_{xi}·{\left(1-{\lambda }_i\right)}^2}{J}\right)\right)$$

(54)

$${\widehat{gg}}_i=\left(\sqrt{r_{min}·r_{max}}\right)\left(\frac{{\left(1-{\lambda }_i\right)}^2}{J·V_x}\right)$$

(55)

$${\widehat{F}}_{roll}=\widehat{M}·g·\widehat{C}$$

(56)

The bounded interval that is used in calculating controller gain (${\mathsf{\Phi }}_i$) is also described in terms of vehicle system parameters follows [43,52]:

$$F_i=\frac{1}{V_x}·\left(\left|1-{\lambda }_i\right|\right)\left({{\displaystyle \sum }}^{ }\left|F_{xi}\right|·\frac{\overline{M}}{\left(\widehat{M}·M_{min}\right)}+g·{\overline{C}}_{ }\right)$$

(57)

where $M_{min}$ is the minimum weight on a wheel; $\widehat{M}$ the estimated weight; $\overline{M}$ the maximum deviation of weight on a wheel; $r_{min}$ the minimum wheel rolling radius; $r_{max}$ the maximum wheel rolling radius; $\widehat{C}$ the estimated value of rolling resistance coefficient; and $\overline{C}$ the maximum deviation of rolling resistance coefficient.

The modified control law, which is robust to bounded disturbance, is estimated as follows:

$$T_{mi}=\frac{1}{{\widehat{gg}}_i}\left(-{\widehat{ff}}_i-\mathsf{\eta }·{\tilde{\lambda }}_i\right)+\frac{1}{{\widehat{gg}}_i}\left({\mathsf{\Phi }}_i·sgn\left(S_{aux,i}\right)\right)$$

(58)

With this modified control law in (58), the Lyapunov function will be satisfied when the condition in the inequality in (59), which is used to determine the controller gain (${\mathsf{\Phi }}_i$), is met.

$${\mathsf{\Phi }}_i\ge \beta ·F_i+\left(\beta -1\right)\left|{\widehat{ff}}_i+\mathsf{\eta }·{\tilde{\lambda }}_i\right|$$

(59)

3.3. Chattering Elimination

The discontinuous sign function $\left(sign\left(S_{aux,i}\right)\right)$ was replaced with a hyperbolic tangent function $\left(tanh\left(S_{aux}\right)\right)$, so that the chattering of the system around $S_{aux}=0,$ which leads to the demand for high control effort, was eliminated [51,52]. The function is defined as follows:

$$\tanh \left(\frac{S_{aux,i}}{\epsilon }\right)=\frac{e^{\frac{S_{aux,i}}{\epsilon }}-e^{-\frac{S_{aux,i}}{\epsilon }}}{e^{\frac{S_{aux,i}}{\epsilon }}+e^{-\frac{S_{aux,i}}{\epsilon }}}$$

(60)

where $\epsilon >0;$ the steepness of the hyperbolic tangent function is determined by ε value.

The control was then modified to

$${\widehat{T}}_{mi}=-\frac{{\widehat{ff}}_i+\mathsf{\eta }·{\tilde{\lambda }}_i+{\mathsf{\Phi }}_i·tanh\left(\frac{S_{aux,i}}{\epsilon }\right)}{{\widehat{gg}}_i}$$

(61)

4. Matlab Simulations

To verify the validity of the control strategy and to tune the controller gains for satisfactory performance before the actual implementation on a real vehicle, a MATLAB/SIMULINK based model of the electric vehicle was built. The values of fixed and bounded parameters required to simulate the algorithm are shown in

Table 2. Parameters of the vehicle’s major components.

Parameters	Value
Vehicle weight (kg)	1200
Width (m)	2.3
Height (m)	2.5
Wheelbase (m)	2.2
Ground clearance (m)	2.2
Rated power (kW)	4.5
Rated torque (Nm)	21.5
Maximum torque (Nm)	54.9
Rated speed (rpm)	2000
Maximum speed (rpm)	3000
Motor speed reduction ratio	32:1
Total wheel inertia (Kgms2)	1.95
Wheel radius (M)	0.326
Center of Gravity (CG)	1.1

and

Table 3. Bounded model parameters.

Parameter/Value	Minimum Value	Estimate	Maximum Value	Maximum Deviation
Mass on the wheels (kg)	200	375	600	225
Rolling resistance coefficient ($C_{roll}$)	0.1	0.15	0.2	0.05

.

The developed wheel slip controller, which is based on integral sliding mode control, the friction force observer, and the optimum slip estimation algorithms, was also implemented in the environment. The diagram in Figure 5 shows the complete Simulink model.

To evaluate the control functionality and to find an optimal controller setup, several test scenarios that include surfaces with different friction conditions, different speeds, and different longitudinal accelerations were considered during the simulation. The controller gains were tuned to get satisfactory tracking performance.

4.1. Tracking of Friction Coefficient and Robustness Performance

Figure 6 shows the optimum friction coefficient estimation algorithm’s performance on terrains with friction coefficient values of 0.2, 0.5, and 0.18, respectively. The algorithm could track an abrupt change in friction coefficient from 0.2 to 0.5 and 0.5 to 0.18, respectively. It converges to a narrow domain, which contains the optimum value within 10 milliseconds (ms).

Figure 7 shows the robustness of the control algorithms to the bounded disturbances in the vertical load on the wheels and the changes in rolling resistance of the terrains, respectively. In Figure 7A, the rolling resistance coefficients used for the simulation were 0.1 (C250), 0.15 (C350), and 0.20 (C450), respectively. The corresponding effect on the tracking ability of the control algorithms was negligible as the graph for overlap each other. Figure 7B, the response of control algorithm due change in the weight on the wheel is shown. Three weights of 250 kg (M250), 350 kg (M350), and 450 kg (M450) were used during the simulation. As in the case of rolling resistance, only a little change is observed in the controller response speed. However, the graphs do not overlap. This is an indication that the system is robust to bounded disturbance.

The developed system was simulated for acceleration performances, with and without the traction control algorithm, respectively, in a low μ surface (0.2). In each case, a step input was applied at the accelerator pedal at 2 s; as a result, the wheel speed increased rapidly and diverged from that of the vehicle when the controller was not used (Figure 8A). At 5 s, the wheel speed reached over 11 m/s, while the vehicle speed was just 1.3 m/s. In the second instance of the simulation, in which the controller was used to the system, the wheel and vehicle velocity were regulated to about 6.4 and 6.0 m/s in 5 s, respectively (Figure 8B).

Furthermore, accelerating the vehicle without the controller leads to an increase in the slip from the values around 0.1 before 2 s to a value up to 0.80 after 2-s of the simulation. However, when the controller was used, the slip was regulated within the range of the optimum value of 0.08–0.11 (Figure 9A). As shown in Figure 9B, the torque supplied by the motor was regulated to less than 40 Nm, unlike in the first instance when it reached as high as 100 Nm despite not accelerating the vehicle to a value comparable to the increase in wheel velocity. This shows that the increased torque was wasted in slipping the wheels.

4.2. μ Transition Simulation

Another simulation was conducted on a variable μ surface; the value of μ was gradually reduced from that of dry concrete surface to wet concrete, ice, and finally snow surface (0.8–0.18); this was to simulate the worst case of μ transition in real field situations. The vehicle model was steadily accelerated in the test from 0 to 5 s by applying a ramp input. Like the previous case, the wheel velocity diverged to more than twice the value of vehicle velocity, as shown in Figure 10A, which resulted in the increased wheel slip (Figure 11A), thereby not accelerating the vehicle significantly.

However, when the traction controller was applied, the slip ratio was regulated within the optimum range. The wheel velocity did not diverge as in the first case (Figure 10A and Figure 11A).

As shown in Figure 11B, the torque supplied by the motor when the traction control was not implemented reached up to 90 Nm before settling down at about 27 Nm. However, when TC was implemented, it increased steadily to a maximum of about 30 Nm, and the fluctuation was less. The increase in torque when TC was not used did not result in increased vehicle velocity, as depicted in Figure 10B, while in the previous instance, it was used up, to increase slip the wheels.

4.3. Field Experiment

In order to verify the performance of the proposed algorithms, experiments were designed and conducted using an experimental platform developed by Nanjing agricultural university (Figure 12), at the Jiangsu Agricultural Machinery testing station. The vehicle was equipped with four in-wheel electric motors, and its parameters were the same as those used during the simulations. The terrain surfaces were dry concrete and moist soil with estimated friction coefficients of 0.8 and 0.2–0.3, respectively.

The velocity and drive torque of each wheel was measured using Panasonic A52 series servo motor driver. The model number and precision of the servo driver were MFME454G1G and 0.1 m/s, respectively. The vehicle velocity and acceleration were measured by P3 Beidou reference station receiver with a speed precision of 0.03 m/s. The control algorithm simulates in Section 3 were converted and uploaded to S7-200 Microcontroller, which is available on the experimental platform. The platform was tested for acceleration on a low μ surface (wet grassy clayey soil surface). Two tests were conducted, one with the traction control engage and the other without engaging the algorithm. In each of the cases, the vehicle was run on a straight line for 35 s at a speed range of about 0–2 m/s. The vehicle was accelerated and decelerated from the beginning to the end of field during the test. The values of wheel torque, slip ratios, wheel speed, and vehicle speed were recorded during the test.

As shown in Figure 13A,B, the slip ratio reached up to 0.8 value when the controller was not engaged, and the value was regulated at less than or equal to 0.15 when the controller was engaged.

By controlling the slip, the wheel speed was maintained close to the vehicle speed as depicted in Figure 14A. In contrast, the wheel speed diverges from the vehicle speed significantly in Figure 14B, during the test without the controller.

The torque performances of the system are shown in Figure 15; it is observed that the torque delivered in the controlled case is about a maximum value of 19 Nm, while in the uncontrolled case, it reached up to 24 Nm; this indicates that the motor delivered more torque in the uncontrolled test, even though the vehicle accelerated more in the test when the controller was applied. Therefore, the controller assists in reducing torque wastage.

On a general note, comparing the controller’s simulated and field response, a similar trend was observed during acceleration and transition from one surface to another. The controller reduces the torque delivered to the wheels, which decreases the wheel speed, thereby suppressing the slip ratio to a lower value that cannot cause soil disturbance and unnecessary loss of energy. However, it is observed that the controller gave more torque (up to 30 Nm) during the simulation, in comparison to the field test, which was about 19 Nm (Figure 11B and Figure 15A). In both simulation and field tests, the trend of the results for slip ratio, torque, and velocity was similar. It can be observed that at some point in the field test results, the torque and slip ratio were negative; this was due to sudden application of brake while maneuvering the vehicle in the field.

5. Conclusions

The study proposed and developed an efficient and feasible means of estimating an optimum slip that corresponds to maximum traction and low soil disturbances. A robust integral sliding mode control scheme was applied to regulate the slip at the optimum estimated value that guaranteed high traction force. Using an instrumented electric vehicle test platform developed by Nanjing agricultural university, the proposed system was verified to be effective during simulation and field experiments. The controller was able to quickly estimate and track the optimum slip despite the modelling uncertainties and changes in terms of weight on wheels and the rolling resistance. The algorithm assigned an optimal torque to the wheels during acceleration on low μ surfaces, and in the transition from high to low μ surfaces. It suppressed the slip ratio from the non-optimal value of about 0.8 and above, to optimal values of less than 0.2. Therefore, the controller is recommended for use in traction control for maximum traction, reduced soil disturbance, and lower power consumption.