# Application of Soft Computing Techniques for the Analysis of Tractive Properties of a Low-Power Agricultural Tractor under Various Soil Conditions

Authors: Katarzyna Pentoś; Krzysztof Pieczarka; Krzysztof Lejman
Journal: Complexity
 (2020)

---

## Abstract

Considering the fuel consumption and soil compaction, optimization of the performance of tractors is crucial for modern agricultural practices. The tractive performance is influenced by many factors, making it difficult to be modeled. In this work, the traction force and tractive efficiency of a low-power tractor, as affected by soil coefficient, vertical load, horizontal deformation, soil compaction, and soil moisture, were studied. The optimal work of a tractor is a compromise between the maximum traction force and the maximum tractive efficiency. Optimizing these factors is complex and requires accurate models. To this end, the performances of soft computing approaches, including neural networks, genetic algorithms, and adaptive network fuzzy inference system, were evaluated. The optimal performance was realized by neural networks trained by backpropagation as well as backpropagation combined with a genetic algorithm, with a coefficient of determination of 0.955 for the traction force and 0.954 for the tractive efficiency. Based on models with the best accuracy, a sensitivity analysis was performed. The results showed that the traction performance is mainly influenced by the soil type; nevertheless, the vertical load and soil moisture also exhibited a relatively strong influence.


---

## Body

## 1. Introduction

The mechanization of agricultural operations is essential for modern agricultural practices. Most agricultural operations are conducted using tractors that exert a high traction force (e.g., during tilling, cultivation, and seeding) or a light traction force (e.g., during harvesting and haymaking). Studies have found that 20–55% of the available tractor power is lost because of the interaction between the tires and topsoil [1]. The traction performance of a tractor significantly affects the fuel consumption and field performance resulting from soil compaction. Therefore, optimizing the performance is crucial for tillage management. The parameters significantly influencing the performance of drive wheels include the tire inflation pressures, wheel slip, and vertical wheel loads [2–5].A mathematical modeling of the soil-tire interaction process can help improve the tractor design and minimize fuel consumption. Such models can be based on empirical, semiempirical, or analytical methods [6]. Analytical models have been developed for predicting traction parameters [7, 8]. However, Tiwari et al. emphasized some of the difficulties limiting the widespread use of analytical models, including the complex tire-soil interaction [6]. Semiempirical models are based on the vertical deformation of the soil and the shear deformation of the soil under a traction device. Rosca et al. proposed using experimentally derived parameters in semianalytical models to predict the traction performance of a driving tractor tire [9]. Empirical models are simpler than analytical and semiempirical models; however, their applicability is limited to cases in which the service and experimental conditions used to develop the model are similar.The complex soil-tire interaction and the limitations in implementing the mathematical models described above have encouraged researchers to develop models based on soft computing techniques. Artificial neural networks (ANNs) have been employed for predicting draught requirement of tillage implements under sandy clay loam soil conditions [10], prognosticating the energy efficiency indices of driven wheels [11], and modeling the relationship between travel reduction-to-net traction ratio and tractive efficiency [12]. Although an ANN is a powerful tool for solving stochastic and complex problems and can produce highly accurate prediction models, more sophisticated (hybrid) soft computing techniques have also been employed [13]. An ANN-genetic algorithm (GA) has been used to model the power of agricultural tractors as a function of the wheel load, slip, and speed [14] and to model the dynamic characteristics of a tractor on sloping terrain [15]. Another hybrid method used for modeling complex relationships in the agricultural field is the adaptive network fuzzy inference system (ANFIS), which combines the advantages of neural networks and fuzzy logic [16, 17].Most studies on topsoil-tire interactions were conducted for high-power tractors (engine power > 73 kW). Studies concerning low-power tractors (engine power < 15 kW) are lacking. A high-power tractor is not always necessary. For example, in horticulture, vineyards, maintenance of green spaces, or in foothill areas on sloping terrain, the use of big, high-power tractors may be difficult. Furthermore, using low-power tractors is likely to reduce machinery costs and decrease soil compaction as a result of low tractor weight [18]. When high traction force is not required, the use of low-power tractors will result in lower fuel consumption, exhaust emissions reduction, and reduction of environmental pollution.The tractive properties of agricultural tractors need to be optimized to minimize fuel consumption. However, the optimization process requires an accurate objective function. Nevertheless, based on measurement data, it is relatively easy to develop an empirical model of soil-tire interaction to be used as the objective function for an optimization algorithm. Considering the nonlinearity and complexity of a soil-tire interaction, soft computing techniques can be used to develop highly accurate models.The objective of this study was to develop ANN, ANN trained by GA, ANN trained by backpropagation (BP) and GA, and ANFIS models to estimate the traction force or tractive efficiency of a low-power tractor as a function of the soil coefficient (which is indicative of soil texture as defined in equation (1)), vertical load, horizontal deformation, soil compaction, and soil moisture. Additionally, various modeling methodologies were compared to determine their accuracy in predicting the experimentally measured traction force or tractive efficiency. An analysis on the importance of predictor variables was conducted with highly accurate models. The parameters that significantly affect the tractive properties and can be varied by the tractor operator to optimize the traction force and tractive efficiency were found.


## 2. Materials and Methods

### 2.1. Experimental Data Acquisition

The tests were conducted on the following types of soil: sand, fine sandy loam, sandy loam, and silty clay loam. To obtain a universal model that can be used for different soil types, the soil texture was determined according to the USDA Comprehensive Soil Survey System [19]. The soil coefficient was calculated as follows:(1)Sc=c1+c2+c3100,where Sc is the soil coefficient [−], c1 is the proportion of medium silt in the test sample [%], c2 is the proportion of fine silt in the test sample [%], and c3 is the proportion of clay in the test sample [%].Based on equation (1), the values of Sc for sand, fine sandy loam, sandy loam, and silty clay loam were calculated to be 0.08, 0.21, 0.33, and 0.68, respectively. Table 1 lists the values of soil moisture corresponding to various soil types, field capacity, beginning of plant growth inhibition, and strong inhibition of plant growth.Table 1
Values of soil moisture (%).
1.25 field water capacity
Field water capacity
Beginning of plant growth inhibition
Strong inhibition of plant growth
Sand
15.31
12.25
10.00
7.00
Fine sandy loam
24.38
19.50
16.00
11.20
Sandy loam
25.63
20.50
18.00
12.60
Silty clay loam
37.81
30.25
27.00
18.90A soil bin testing facility (see Figure1) was used in the laboratory to provide a controlled environment for evaluating the interaction between the soil and the tire. The tractive force was measured using a load cell with a measurement range of up to 1 kN and a precision of 1 N.Figure 1
Schematic of a soil bin testing facility: 1, linear potentiometer; 2, wheel, 3, box with soil; 4, load cell; 5, data agent.Low-power tractors are usually equipped with wheels having a rim diameter of 10 inches. On this type of rim, tires with different overall widths can be mounted, e.g., 4.00–10, 4.50–10, and 5.00–10. In the case of 5.00–10 tires, the contact area is the largest, resulting in a low-unit pressure and less soil compaction. Therefore, this tire is particularly advantageous for soils of low compaction because the larger tread blocks are unlikely to penetrate a highly compact soil. Conversely, in the case of highly compact soils, a 4.00–10 tire would be most beneficial, as the narrow tread blocks will be able to extend deeper into the soil ensuring an optimum traction force. In this work, experiments were conducted on soils of both low and high compactions. Therefore, the drive tire Kenda 4.50–10, type K365 for small agricultural tractors was used for the measurements. Considering that tractor tires must to some extent be universal, the results obtained for Kenda 4.50–10 can be extended to other tires manufactured for low-power tractors.During the measurement, the peripheral speed of the wheel was constant and low (0.3 rad/s). The decision to use a low speed value was caused by the need to eliminate the impact of dynamic phenomena that could affect the results when the tests were performed on soils with different textures. The following five vertical load values were used in this work: 375, 490, 638, 785, and 932 N. The values of the vertical load were chosen on the basis of technical documentation provided by the manufacturers of the low-power tractors. The average vertical load recommended for low-power tractors is in the range of 800–900 N. The aim of this research was to develop mathematical models of the relationships under study, which can subsequently be used for optimizing the selected operating parameters in order to minimize soil compaction. Therefore, the authors selected three values of the vertical load, all of which were lower than that recommended for low-power tractors. The influence of tire inflation pressure on the tractive properties was not a subject of this research. Thus, all the measurements were conducted using a constant pressure of 0.16 MPa (as recommended by the tire manufacturer).For measuring the static-loaded radius, the tire deflection under certain vertical load and tire inflation pressures was determined. The measurement method is detailed in Figure2(a).Figure 2
Method of static loaded radius measurement (a) and graphical representation of the parameters used for horizontal deformation calculation (b).
(a)
(b)The static loaded radius was calculated as follows:(2)rs=r−e,where rs is the static loaded radius, r is the tire radius when the vertical load is equal to 0, and e is the tire deflection.Based on the static loaded radius values and wheel rotation angle, the horizontal deformation can be calculated as follows:(3)j=α⋅rs⋅π180,where j is the horizontal deformation, α is the wheel rotation angle, and rs is the static loaded radius (see Figure 2(b)).The result of the measurements conducted using the soil bin testing facility (see Figure1) was the traction force as a function of the rotation angle of the wheel. For each combination of the independent variables, the traction force increases with an increase in the horizontal deformation. Consequently, the maximum traction force could not be determined. Therefore, the horizontal deformation that produced a wheel slip of 20% (the slip limit accepted in the case of agricultural tractors) was determined analytically, and the traction force corresponding to this horizontal deformation was considered the maximum. The wheel slip depends on both the horizontal and vertical deformations. Hence, the horizontal deformation that produces a wheel slip of 20% can be determined only for a certain value of the vertical deformation (affected by vertical load, soil type, soil moisture, and soil compaction). In this research, the vertical deformation was calculated as an arithmetic mean of the maximum values measured for each soil type. The horizontal deformation corresponding to this vertical deformation and under a slip of 20% is found to be 0.05 m. Therefore, the traction force measured for a horizontal deformation of 0.05 m was considered the maximum. To determine the relationship between the traction force and the horizontal deformation, measurements were also taken for horizontal deformations of 0.02, 0.03, and 0.04 m (corresponding to slips of 5, 10, and 15%, respectively).The tractive efficiency was calculated as follows:(4)η=∫0jPTjdj∫0jPTjdj+G⋅h,where η is the tractive efficiency, PT is the traction force (N), j is the horizontal deformation (m), G is the vertical load of the wheel (N), and h is the rut depth (m).Table2 lists the statistics of the experimental data.Table 2
Statistics of experimental data by the soil type.
The parameter
Minimum
Maximum
Mean
Standard deviation
Horizontal deformation (m)
0.01
0.05
0.03
0.01
Vertical load (N)
375.00
932.00
644.00
199.53
Sand
Soil compaction (kPa)
201.04
571.85
386.45
124.51
Soil moisture (%)
7.00
15.00
11.00
2.91
Traction force (N)
65.81
596.56
284.91
114.27
Tractive efficiency (%)
5.04
85.55
33.61
17.43
Fine sandy loam
Soil compaction (kPa)
102.69
498.11
300.85
125.92
Soil moisture (%)
11.00
24.00
17.76
4.82
Traction force (N)
99.78
787.63
339.94
140.15
Tractive efficiency (%)
4.45
85.71
33.44
17.39
Sandy loam
Soil compaction (kPa)
90.88
486.30
288.59
125.84
Soil moisture (%)
13.00
26.00
19.50
4.72
Traction force (N)
36.09
683.61
326.27
133.30
Tractive efficiency (%)
0.99
64.98
27.97
14.38
Silty clay loam
Soil compaction (kPa)
95.89
520.83
308.36
128.13
Soil moisture (%)
19.00
38.00
28.50
6.81
Traction force (N)
48.83
770.65
362.24
150.15
Tractive efficiency (%)
2.05
73.77
31.80
14.91The 1600 datasets (the vectors of measured parameters) obtained during the measurement process were randomly separated into training (80%) and validation (20%) sets. Table2 lists the minimum and maximum input and output parameter values. Prior to utilizing the dataset for model development, the data were normalized to a range of 0-1 using the following equation:(5)NV=V−VminVmax−Vmin,where NV is the input or output normalized vector, V is the input or output data, Vmax is the maximum of the input or output data, and Vmin is the minimum of the input or output data. We followed the methods of Pentoś and Pieczarka [20].


### 2.2. Artificial Neural Networks

An ANN is a highly simplified model of the biological structure of neurons in the human nervous system. An ANN is considered an effective substitute for the empirical and statistical process modeling techniques and is widely used in agricultural applications. In this work, a feed-forward neural network, namely, multilayer perceptron (MLP), was used. The ability of an ANN strongly depends on its topology, i.e., the number of hidden layers and the number of neurons in these layers. The optimal topology and learning parameters are usually determined by a trial and error method, requiring many simulations. For this study, an MLP with a single hidden layer was chosen as the ANN architecture. The input layer was composed of five nodes (soil coefficient, vertical load, horizontal deformation, soil compaction, and soil moisture). The number of neurons in the hidden layer was set to a range of 10–40, and nonlinear sigmoid neurons were implemented in this layer. There was one neuron in the output layer, producing the predicted value of the traction force or tractive efficiency. For each ANN architecture, 10 simulations were performed, and as a result, 310 ANNs were trained for each output parameter. The following were the three training methods used for the simulations: resilient backpropagation with and without weight backtracking and a modified globally convergent algorithm. The resilient backpropagation algorithm is based on the traditional backpropagation; however, in this algorithm, a separate learning rateηk is used for each weight in the network and can be changed during the training process. Contrary to the traditional backpropagation, in resilient backpropagation, only the sign of the partial derivatives is used to indicate the direction of weight updation. The weights were modified using the following equation [21]:(6)wkt+1=wkt−ηkt⋅sign∂Et∂wkt,where wk is the kth connection weight and E is the error function. Weight backtracking implies weight update reversal when the sign of the partial derivative changes [22]:(7)Δwkt=−Δwkt−1,if∂Et−1∂wkt⋅∂Et∂wkt<0.The modified globally convergent algorithm presented by Anastasiadis et al. is based on resilient backpropagation. A new modification to the learning rate is proposed [23]:(8)ηit=−∑k;k≠inηkt⋅∂Et/∂wkt+δ∂Et/∂wit,where ∂Et/∂wit≠0 and 0 < δ << ∞. This modification improves the convergence speed and stability of the learning algorithm.The “neuralnet” package, version 1.44.2, for the R environment (R Foundation for Statistical Computing (https://www.r-project.org/)), was used for the simulations [24].


### 2.3. Artificial Neural Network Combined with Genetic Algorithm

A typical problem with the ANN trained by an algorithm based on backpropagation is the possibility of falling into a local minimum of the error function, resulting in a slow convergence. Therefore, this technique needs to be improved by hybridizing the ANN with an optimization tool such as the GA. The GA, suggested by Holland, can be a pragmatic alternative to conventional local search methods [25]. In this study, two hybridizations of the ANN combined with GA were used. The first one (ANN + GA) utilizes a GA to realize the optimal allocation for the given network weights and biases starting from the random values. The second one (ANN_BP + GA) uses one of the algorithms based on backpropagation (resilient backpropagation with and without weight backtracking and modified globally convergent algorithm) for initial ANN training. The GA is then employed for the final optimization starting from the initial chromosome population produced by the ANN training. The chromosome in this work is the vector of genes, real numbers representing the ANN weights and biases. The following operations are performed on chromosome population during GA: selection and genetic operations (crossover and mutation). During the selection, the population of the chromosomes is chosen based on the fitness function for genetic operations. In the selection procedure, the probability that an individual can become a parent should be higher for higher fitness functions. A crossover operation is when two individuals (parents) exchange genes with each other. A crossover is performed with probability Pc, which is usually high. Mutation is a small random tweak in the chromosome that leads to a new individual. It is performed with probability Pm, which is usually low. The function of a mutation operation is to ensure a higher population diversity, consequently preventing the GA from falling into local extremes. After reaching the maximum generation, the GA converges to produce the best chromosome, which represents an optimal or near-optimal solution. In this work, the population size was set to 100 chromosomes, the crossover probability was set to 0.8, and the mutation probability was set to 0.01. As a fitness function, the root mean square error (RMSE) of the ANN was used for the validation dataset. The following were the three selection methods used: roulette wheel, tournament selection, and fitness proportional selection with fitness linear scaling. The genetic operations were performed with local arithmetic crossover and uniform random mutation. An R package, “GA” version 3.0.2, for optimization with the GA was used for the simulations [26].


### 2.4. Adaptive Network Fuzzy Inference System

The ANFIS is a global search soft computing technique, which combines the advantages of fuzzy logic and ANN. It is based on the first-order Takagi–Sugeno fuzzy inference system introduced by Jang [27]. This machine-learning technique generates fuzzy rules from a given input/output dataset and can adjust the membership function parameters directly from the data during the training process. The membership function parameters are adjusted using a combination of gradient descent and the least squares method. The typical fuzzy IF-THEN rules for the first-order Takagi–Sugeno fuzzy model are as follows:(9)IFx1=A1andx2=A3THENf=p11x1+p21x2+r1,IFx1=A2andx2=A4THENf=p12x1+p22x2+r2.The ANFIS architecture contains a five-layer feed-forward neural network. Figure3 shows the architecture used in the present work.Figure 3
Adaptive neurofuzzy inference system structure.Layer 1 is the fuzzification layer. Each node in this layer represents a membership function (MF) and defines the membership grades for each set of input. There are various types of MFs. Because a normalized Gaussian function is used in this study, the output of this layer is given by(10)O1,n=μAnx=exp−x−cn22σn2,where cn and σn are the parameters that make up a premise set.Layer 2 is a multiplicative layer. Each node uses a multiplication operator and calculates the firing strength of the rule as a product of the previous membership grades. For instance, for the first node, this is given by(11)O2,1=w1=μA1x1μA4x2μA7x3μA10x4μA13x5.Layer 3 normalizes the firing strength of the rules:(12)O3,n=w¯n=wn∑jwj.Layer 4 consists adaptive nodes that compute a linear function in which the parameterspn and rn are adapted using the error function of the feed-forward neural network:(13)O4,n=w¯n⋅fn=w¯n⋅∑kpknxk+rn.Layer 5 has a single node and produces the output signal of the ANFIS, which is the sum of the outputs of the nodes from layer 4:(14)O5,1=∑nO4,n.The package “anfis”—Adaptive Neurofuzzy Inference System in R, version 0.99.1, was used for the simulations. Excessive membership functions in the ANFIS model is not appropriate because many parameters need to be predicted. Therefore, the number of membership functions was set to 3, and the number of iterations was set to 20.


### 2.5. Comparison Criteria

The accuracy of the models were compared in terms of the mean absolute error (MAE), RMSE, and coefficient of determination (R2), which are expressed as follows:(15)MAE=1n∑i=1nYpredicted−Ymeasured,RMSE=1n∑i=1nYpredicted−Ymeasured2,R2=∑Ymeasured−Ymean_measYpredicted−Ymean_predict​∑Ymeasured−Ymean_meas2∑Ypredicted−Ymean_predict22,where Ypredicted and Ymean_predict are the absolute and average predicted values and Ymeasured and Ymean_meas are the absolute and average measured values, respectively.The MAE, using the absolute values of the differences between the measured and predicted values, rates over- and underestimations equally. The RMSE is an indicator of the distribution of positive and negative errors of the estimated values.R2 estimates the strength of the relationship between the output value calculated using a model and the expected value. Therefore, the closer MAE and RMSE are to 0 and the closer R2 is to 1, the better is the accuracy of the model in estimating the dependent variable.


### 2.6. Sensitivity Analysis

Based on the mathematical model, the importance of the independent variables can be estimated. Many methods have been proposed for model sensitivity analysis. The relevance of the method depends on the characteristics of the particular model. In this research, the partial derivatives method, dedicated for a neural network model, was used, whereby the contribution of the input variables is determined based on the connection weights and a bias matrix [28]. It is difficult to select an optimal ANN model architecture; thus, the contribution of the predictor variables should be determined based on a group of ANN models [29]. In the present research, a group of twenty ANN models with the highest R2 values and the lowest MAE and RMSE values was selected. As the final result for each dependent variable (traction force and tractive efficiency), the arithmetical mean of the results produced by the twenty ANNs was calculated.


## 2.1. Experimental Data Acquisition

The tests were conducted on the following types of soil: sand, fine sandy loam, sandy loam, and silty clay loam. To obtain a universal model that can be used for different soil types, the soil texture was determined according to the USDA Comprehensive Soil Survey System [19]. The soil coefficient was calculated as follows:(1)Sc=c1+c2+c3100,where Sc is the soil coefficient [−], c1 is the proportion of medium silt in the test sample [%], c2 is the proportion of fine silt in the test sample [%], and c3 is the proportion of clay in the test sample [%].Based on equation (1), the values of Sc for sand, fine sandy loam, sandy loam, and silty clay loam were calculated to be 0.08, 0.21, 0.33, and 0.68, respectively. Table 1 lists the values of soil moisture corresponding to various soil types, field capacity, beginning of plant growth inhibition, and strong inhibition of plant growth.Table 1
Values of soil moisture (%).
1.25 field water capacity
Field water capacity
Beginning of plant growth inhibition
Strong inhibition of plant growth
Sand
15.31
12.25
10.00
7.00
Fine sandy loam
24.38
19.50
16.00
11.20
Sandy loam
25.63
20.50
18.00
12.60
Silty clay loam
37.81
30.25
27.00
18.90A soil bin testing facility (see Figure1) was used in the laboratory to provide a controlled environment for evaluating the interaction between the soil and the tire. The tractive force was measured using a load cell with a measurement range of up to 1 kN and a precision of 1 N.Figure 1
Schematic of a soil bin testing facility: 1, linear potentiometer; 2, wheel, 3, box with soil; 4, load cell; 5, data agent.Low-power tractors are usually equipped with wheels having a rim diameter of 10 inches. On this type of rim, tires with different overall widths can be mounted, e.g., 4.00–10, 4.50–10, and 5.00–10. In the case of 5.00–10 tires, the contact area is the largest, resulting in a low-unit pressure and less soil compaction. Therefore, this tire is particularly advantageous for soils of low compaction because the larger tread blocks are unlikely to penetrate a highly compact soil. Conversely, in the case of highly compact soils, a 4.00–10 tire would be most beneficial, as the narrow tread blocks will be able to extend deeper into the soil ensuring an optimum traction force. In this work, experiments were conducted on soils of both low and high compactions. Therefore, the drive tire Kenda 4.50–10, type K365 for small agricultural tractors was used for the measurements. Considering that tractor tires must to some extent be universal, the results obtained for Kenda 4.50–10 can be extended to other tires manufactured for low-power tractors.During the measurement, the peripheral speed of the wheel was constant and low (0.3 rad/s). The decision to use a low speed value was caused by the need to eliminate the impact of dynamic phenomena that could affect the results when the tests were performed on soils with different textures. The following five vertical load values were used in this work: 375, 490, 638, 785, and 932 N. The values of the vertical load were chosen on the basis of technical documentation provided by the manufacturers of the low-power tractors. The average vertical load recommended for low-power tractors is in the range of 800–900 N. The aim of this research was to develop mathematical models of the relationships under study, which can subsequently be used for optimizing the selected operating parameters in order to minimize soil compaction. Therefore, the authors selected three values of the vertical load, all of which were lower than that recommended for low-power tractors. The influence of tire inflation pressure on the tractive properties was not a subject of this research. Thus, all the measurements were conducted using a constant pressure of 0.16 MPa (as recommended by the tire manufacturer).For measuring the static-loaded radius, the tire deflection under certain vertical load and tire inflation pressures was determined. The measurement method is detailed in Figure2(a).Figure 2
Method of static loaded radius measurement (a) and graphical representation of the parameters used for horizontal deformation calculation (b).
(a)
(b)The static loaded radius was calculated as follows:(2)rs=r−e,where rs is the static loaded radius, r is the tire radius when the vertical load is equal to 0, and e is the tire deflection.Based on the static loaded radius values and wheel rotation angle, the horizontal deformation can be calculated as follows:(3)j=α⋅rs⋅π180,where j is the horizontal deformation, α is the wheel rotation angle, and rs is the static loaded radius (see Figure 2(b)).The result of the measurements conducted using the soil bin testing facility (see Figure1) was the traction force as a function of the rotation angle of the wheel. For each combination of the independent variables, the traction force increases with an increase in the horizontal deformation. Consequently, the maximum traction force could not be determined. Therefore, the horizontal deformation that produced a wheel slip of 20% (the slip limit accepted in the case of agricultural tractors) was determined analytically, and the traction force corresponding to this horizontal deformation was considered the maximum. The wheel slip depends on both the horizontal and vertical deformations. Hence, the horizontal deformation that produces a wheel slip of 20% can be determined only for a certain value of the vertical deformation (affected by vertical load, soil type, soil moisture, and soil compaction). In this research, the vertical deformation was calculated as an arithmetic mean of the maximum values measured for each soil type. The horizontal deformation corresponding to this vertical deformation and under a slip of 20% is found to be 0.05 m. Therefore, the traction force measured for a horizontal deformation of 0.05 m was considered the maximum. To determine the relationship between the traction force and the horizontal deformation, measurements were also taken for horizontal deformations of 0.02, 0.03, and 0.04 m (corresponding to slips of 5, 10, and 15%, respectively).The tractive efficiency was calculated as follows:(4)η=∫0jPTjdj∫0jPTjdj+G⋅h,where η is the tractive efficiency, PT is the traction force (N), j is the horizontal deformation (m), G is the vertical load of the wheel (N), and h is the rut depth (m).Table2 lists the statistics of the experimental data.Table 2
Statistics of experimental data by the soil type.
The parameter
Minimum
Maximum
Mean
Standard deviation
Horizontal deformation (m)
0.01
0.05
0.03
0.01
Vertical load (N)
375.00
932.00
644.00
199.53
Sand
Soil compaction (kPa)
201.04
571.85
386.45
124.51
Soil moisture (%)
7.00
15.00
11.00
2.91
Traction force (N)
65.81
596.56
284.91
114.27
Tractive efficiency (%)
5.04
85.55
33.61
17.43
Fine sandy loam
Soil compaction (kPa)
102.69
498.11
300.85
125.92
Soil moisture (%)
11.00
24.00
17.76
4.82
Traction force (N)
99.78
787.63
339.94
140.15
Tractive efficiency (%)
4.45
85.71
33.44
17.39
Sandy loam
Soil compaction (kPa)
90.88
486.30
288.59
125.84
Soil moisture (%)
13.00
26.00
19.50
4.72
Traction force (N)
36.09
683.61
326.27
133.30
Tractive efficiency (%)
0.99
64.98
27.97
14.38
Silty clay loam
Soil compaction (kPa)
95.89
520.83
308.36
128.13
Soil moisture (%)
19.00
38.00
28.50
6.81
Traction force (N)
48.83
770.65
362.24
150.15
Tractive efficiency (%)
2.05
73.77
31.80
14.91The 1600 datasets (the vectors of measured parameters) obtained during the measurement process were randomly separated into training (80%) and validation (20%) sets. Table2 lists the minimum and maximum input and output parameter values. Prior to utilizing the dataset for model development, the data were normalized to a range of 0-1 using the following equation:(5)NV=V−VminVmax−Vmin,where NV is the input or output normalized vector, V is the input or output data, Vmax is the maximum of the input or output data, and Vmin is the minimum of the input or output data. We followed the methods of Pentoś and Pieczarka [20].


## 2.2. Artificial Neural Networks

An ANN is a highly simplified model of the biological structure of neurons in the human nervous system. An ANN is considered an effective substitute for the empirical and statistical process modeling techniques and is widely used in agricultural applications. In this work, a feed-forward neural network, namely, multilayer perceptron (MLP), was used. The ability of an ANN strongly depends on its topology, i.e., the number of hidden layers and the number of neurons in these layers. The optimal topology and learning parameters are usually determined by a trial and error method, requiring many simulations. For this study, an MLP with a single hidden layer was chosen as the ANN architecture. The input layer was composed of five nodes (soil coefficient, vertical load, horizontal deformation, soil compaction, and soil moisture). The number of neurons in the hidden layer was set to a range of 10–40, and nonlinear sigmoid neurons were implemented in this layer. There was one neuron in the output layer, producing the predicted value of the traction force or tractive efficiency. For each ANN architecture, 10 simulations were performed, and as a result, 310 ANNs were trained for each output parameter. The following were the three training methods used for the simulations: resilient backpropagation with and without weight backtracking and a modified globally convergent algorithm. The resilient backpropagation algorithm is based on the traditional backpropagation; however, in this algorithm, a separate learning rateηk is used for each weight in the network and can be changed during the training process. Contrary to the traditional backpropagation, in resilient backpropagation, only the sign of the partial derivatives is used to indicate the direction of weight updation. The weights were modified using the following equation [21]:(6)wkt+1=wkt−ηkt⋅sign∂Et∂wkt,where wk is the kth connection weight and E is the error function. Weight backtracking implies weight update reversal when the sign of the partial derivative changes [22]:(7)Δwkt=−Δwkt−1,if∂Et−1∂wkt⋅∂Et∂wkt<0.The modified globally convergent algorithm presented by Anastasiadis et al. is based on resilient backpropagation. A new modification to the learning rate is proposed [23]:(8)ηit=−∑k;k≠inηkt⋅∂Et/∂wkt+δ∂Et/∂wit,where ∂Et/∂wit≠0 and 0 < δ << ∞. This modification improves the convergence speed and stability of the learning algorithm.The “neuralnet” package, version 1.44.2, for the R environment (R Foundation for Statistical Computing (https://www.r-project.org/)), was used for the simulations [24].


## 2.3. Artificial Neural Network Combined with Genetic Algorithm

A typical problem with the ANN trained by an algorithm based on backpropagation is the possibility of falling into a local minimum of the error function, resulting in a slow convergence. Therefore, this technique needs to be improved by hybridizing the ANN with an optimization tool such as the GA. The GA, suggested by Holland, can be a pragmatic alternative to conventional local search methods [25]. In this study, two hybridizations of the ANN combined with GA were used. The first one (ANN + GA) utilizes a GA to realize the optimal allocation for the given network weights and biases starting from the random values. The second one (ANN_BP + GA) uses one of the algorithms based on backpropagation (resilient backpropagation with and without weight backtracking and modified globally convergent algorithm) for initial ANN training. The GA is then employed for the final optimization starting from the initial chromosome population produced by the ANN training. The chromosome in this work is the vector of genes, real numbers representing the ANN weights and biases. The following operations are performed on chromosome population during GA: selection and genetic operations (crossover and mutation). During the selection, the population of the chromosomes is chosen based on the fitness function for genetic operations. In the selection procedure, the probability that an individual can become a parent should be higher for higher fitness functions. A crossover operation is when two individuals (parents) exchange genes with each other. A crossover is performed with probability Pc, which is usually high. Mutation is a small random tweak in the chromosome that leads to a new individual. It is performed with probability Pm, which is usually low. The function of a mutation operation is to ensure a higher population diversity, consequently preventing the GA from falling into local extremes. After reaching the maximum generation, the GA converges to produce the best chromosome, which represents an optimal or near-optimal solution. In this work, the population size was set to 100 chromosomes, the crossover probability was set to 0.8, and the mutation probability was set to 0.01. As a fitness function, the root mean square error (RMSE) of the ANN was used for the validation dataset. The following were the three selection methods used: roulette wheel, tournament selection, and fitness proportional selection with fitness linear scaling. The genetic operations were performed with local arithmetic crossover and uniform random mutation. An R package, “GA” version 3.0.2, for optimization with the GA was used for the simulations [26].


## 2.4. Adaptive Network Fuzzy Inference System

The ANFIS is a global search soft computing technique, which combines the advantages of fuzzy logic and ANN. It is based on the first-order Takagi–Sugeno fuzzy inference system introduced by Jang [27]. This machine-learning technique generates fuzzy rules from a given input/output dataset and can adjust the membership function parameters directly from the data during the training process. The membership function parameters are adjusted using a combination of gradient descent and the least squares method. The typical fuzzy IF-THEN rules for the first-order Takagi–Sugeno fuzzy model are as follows:(9)IFx1=A1andx2=A3THENf=p11x1+p21x2+r1,IFx1=A2andx2=A4THENf=p12x1+p22x2+r2.The ANFIS architecture contains a five-layer feed-forward neural network. Figure3 shows the architecture used in the present work.Figure 3
Adaptive neurofuzzy inference system structure.Layer 1 is the fuzzification layer. Each node in this layer represents a membership function (MF) and defines the membership grades for each set of input. There are various types of MFs. Because a normalized Gaussian function is used in this study, the output of this layer is given by(10)O1,n=μAnx=exp−x−cn22σn2,where cn and σn are the parameters that make up a premise set.Layer 2 is a multiplicative layer. Each node uses a multiplication operator and calculates the firing strength of the rule as a product of the previous membership grades. For instance, for the first node, this is given by(11)O2,1=w1=μA1x1μA4x2μA7x3μA10x4μA13x5.Layer 3 normalizes the firing strength of the rules:(12)O3,n=w¯n=wn∑jwj.Layer 4 consists adaptive nodes that compute a linear function in which the parameterspn and rn are adapted using the error function of the feed-forward neural network:(13)O4,n=w¯n⋅fn=w¯n⋅∑kpknxk+rn.Layer 5 has a single node and produces the output signal of the ANFIS, which is the sum of the outputs of the nodes from layer 4:(14)O5,1=∑nO4,n.The package “anfis”—Adaptive Neurofuzzy Inference System in R, version 0.99.1, was used for the simulations. Excessive membership functions in the ANFIS model is not appropriate because many parameters need to be predicted. Therefore, the number of membership functions was set to 3, and the number of iterations was set to 20.


## 2.5. Comparison Criteria

The accuracy of the models were compared in terms of the mean absolute error (MAE), RMSE, and coefficient of determination (R2), which are expressed as follows:(15)MAE=1n∑i=1nYpredicted−Ymeasured,RMSE=1n∑i=1nYpredicted−Ymeasured2,R2=∑Ymeasured−Ymean_measYpredicted−Ymean_predict​∑Ymeasured−Ymean_meas2∑Ypredicted−Ymean_predict22,where Ypredicted and Ymean_predict are the absolute and average predicted values and Ymeasured and Ymean_meas are the absolute and average measured values, respectively.The MAE, using the absolute values of the differences between the measured and predicted values, rates over- and underestimations equally. The RMSE is an indicator of the distribution of positive and negative errors of the estimated values.R2 estimates the strength of the relationship between the output value calculated using a model and the expected value. Therefore, the closer MAE and RMSE are to 0 and the closer R2 is to 1, the better is the accuracy of the model in estimating the dependent variable.


## 2.6. Sensitivity Analysis

Based on the mathematical model, the importance of the independent variables can be estimated. Many methods have been proposed for model sensitivity analysis. The relevance of the method depends on the characteristics of the particular model. In this research, the partial derivatives method, dedicated for a neural network model, was used, whereby the contribution of the input variables is determined based on the connection weights and a bias matrix [28]. It is difficult to select an optimal ANN model architecture; thus, the contribution of the predictor variables should be determined based on a group of ANN models [29]. In the present research, a group of twenty ANN models with the highest R2 values and the lowest MAE and RMSE values was selected. As the final result for each dependent variable (traction force and tractive efficiency), the arithmetical mean of the results produced by the twenty ANNs was calculated.


## 3. Results and Discussion

### 3.1. Soft Computing Models

The optimal solution, according to energy savings, is to have high values of both traction force and tractive efficiency. Figure4 shows the dependence of the traction force and tractive efficiency on the vertical load and soil moisture measured on sand, for a horizontal deformation of 0.05 m.Figure 4
3D surface curves of traction force (a) and tractive efficiency (b) as affected by the interactions of vertical load and soil moisture.
(a)
(b)As shown in Figure4, the optimal traction force is produced under a high vertical load, whereas an optimal tractive efficiency can be achieved under a rather low vertical load. Therefore, achieving an optimal balance between the traction force and the tractive efficiency is very difficult and requires accurate mathematical models of the tractive properties.Before the development of the models based on the ANN, linearly dependent predictor variables must be removed from the dataset. Therefore, Pearson’s correlation coefficients between the explanatory variables were calculated. A high positive correlation was observed only between the soil coefficient and the soil moisture (r = 0.767). However, after analyzing the experimental procedures, it was evident that both soil coefficient and soil moisture must be considered as input variables and used for model development. The correlation coefficient is high because the soil moisture range (Table 1) was determined depending on the soil texture.Table3 lists the statistical parameters, namely, R2, RMSE, and MAE, for the best configurations of the mentioned models. The calculations were performed using normalized data (equation (5)).Table 3
Error metrics of best model performances.
Model
Train
Validation
MAE
RMSE
R
2
MAE
RMSE
R
2
Traction force
ANN
0.026
0.037
0.957
0.029
0.040
0.954
ANN + GA
0.064
0.084
0.808
0.063
0.087
0.820
ANN_BP + GA
0.026
0.037
0.958
0.028
0.040
0.955
ANFIS
0.042
0.061
0.892
0.045
0.064
0.888
Tractive efficiency
ANN
0.022
0.030
0.975
0.024
0.037
0.954
ANN + GA
0.076
0.102
0.789
0.068
0.096
0.752
ANN_BP + GA
0.022
0.030
0.975
0.024
0.037
0.954
ANFIS
0.038
0.061
0.883
0.040
0.064
0.872In the case of the traction force model, the best architecture of the ANN was a network with 28 neurons in the hidden layer, and this network was chosen for weights optimization by the GA. The best ANN architecture trained by the GA starting from random values (ANN + GA) was a network with 14 neurons in the hidden layer. The best neural network for the tractive efficiency model had 26 neurons in the hidden layer, and the best ANN + GA model contained 17 neurons in the hidden layer.Figures5 and 6 show the measured and predicted traction force and tractive efficiency values for all the predictive models on the validation dataset. These graphs show the number of data points located very close to the diagonal line, thus facilitating the assessment of the model accuracy.Figure 5
Scatterplot of model-predicted values versus actual values of traction force: (a) ANN, (b) ANN + GA, (c) ANN_BP + GA, and (d) ANFIS.
(a)
(b)
(c)
(d)Figure 6
Scatterplot of model-predicted values versus actual values of tractive efficiency: (a) ANN, (b) ANN + GA, (c) ANN_BP + GA, and (d) ANFIS.
(a)
(b)
(c)
(d)As listed in Table3 and shown in Figures 5 and 6, for both output model parameters, among all the computational models, the ANN and ANN_BP + GA models exhibit the best performance, as indicated by high values of R2 (0.954 and 0.955 for traction force and 0.954 for tractive efficiency) and low values of MAE and RMSE for the validation dataset. The use of the GA for optimizing the weights and biases adjusted by the BP algorithm produced slightly better accuracy in the case of the traction force model. The accuracy of the ANFIS model was lower than those of the ANN and ANN_BP + GA models, with R2 values below 0.9. Additionally, the computational time required for the calculations during the ANFIS model development was significantly higher than that required in the case of models based on MLP. The ANN + GA technique seems to be unsuitable, exhibiting a low accuracy in estimating the traction force and tractive efficiency (R2 = 0.820 and 0.752 for the validation dataset, respectively). Generally, it can be stated that, in agriculture, mathematical models (also based on machine learning) with coefficient of determination (R2) exceeding 0.9 are useful for real life applications [30].Neural networks and hybrid methods were also used by other researchers to model the behavior of agricultural tractors. The ANFIS-based modeling was found to be a promising technique for prognosticating the traction coefficient and tractive power efficiency, withR2 values of 0.98 and 0.97, respectively [17] and for prognosticating the drawbar pull energy of tractor driving wheels, with MSE and R2 values of 0.00236 and 0.995, respectively [16]. In the case of ANN combined with a GA, Taghavifar et al. demonstrated that this method drastically decreased the error and increased the performance of the model of power provided by agricultural tractors as affected by wheel load, slip, and speed [14]. They obtained high values for the coefficient of determination for the ANN + GA model: 0.9696 for the training dataset and 0.9672 for validation dataset.Comparing the current results with those presented by other researchers, it is unclear which technique most accurately models the nonlinear and complex relationships such as the ones investigated in this study. Similar results were obtained by other researchers. Johann et al. compared computational models based on ANN and ANFIS in estimating the soil moisture from the stochastic information of the horizontal and vertical forces acting on a no-till chisel opener using autoregressive error function parameters [31]. The ANN model (R2 = 0.79 and RMSE = 1.27) outperformed the ANFIS model (R2 = 0.69 and RMSE = 1.62) in the test phase. Citakoglu applied ANN and ANFIS for estimating the solar radiation in Turkey using the calendar month number and pertinent meteorological data and obtained a higher accuracy when using the ANN (R2 = 0.930 and RMSE = 1.650) in comparison with using the ANFIS (R2 = 0.926 and RMSE = 1.691) [32]. In contrast, the ANFIS was found to be more suitable than the ANN for estimating the soil cation exchange capacity as affected by clay, silt, sand, organic carbon, and pH in arid rangeland ecosystems and for estimating the oxidation parameters of Kilka oil [33, 34]. Based on relevant error (RE) values, Ping and Fei showed that the accuracy of an ANN combined with a GA (RE = 1.48%) is better than that of a traditional ANN model (RE = 3.91%) for Guangdong port throughput forecasting [35]. Similarly, Srinivasulu and Jain found that the predictive capability of an ANN combined with GA rainfall-runoff models is better than that trained using a BP algorithm [36].


### 3.2. Sensitivity Analysis

A highly accurate mathematical model can give additional information about the relationships under study. A sensitivity analysis has been performed to determine the contribution of independent variables in black box data mining models. The neural network trained by backpropagation combined with a GA was found to be the best model for the relationships analyzed in the present research. This soft computing technique was used for developing the models for the sensitivity analysis. For each dependent variable (traction force and tractive efficiency), a group of twenty ANN models was developed. Table4 lists the parameters of the models.Table 4
Statistics of neural model architectures used for sensitivity analysis.
Dependent variable
The range of number of neurons in hidden layer
The range ofR2 for validation dataset
The range of MAE for validation dataset
The range of RMSE for validation dataset
Traction force
16–40
0.944–0.955
0.028–0.032
0.040–0.045
Tractive efficiency
16–38
0.928–0.954
0.024–0.032
0.037–0.047The results of the relative importance of the input parameters for each model were determined as the arithmetical mean of the results produced by the group of twenty ANN models. The results revealed that the traction force and tractive efficiency are most affected by the soil type (58.3 and 74.5%, respectively). This is in agreement with the results reported by other authors [37, 38]. The two additional parameters that significantly influenced the traction force and tractive efficiency are the vertical load (18.3 and 10.1%, respectively) and soil moisture (19.8 and 10.3%, respectively). The significant effects of these parameters on the tractive performance have been highlighted in other studies as well [11, 39]. It is worth emphasizing that the vertical load is one of the more easily manageable parameters during traction performance optimization. Soil moisture can also be varied in a certain range, as the operator can advance or delay agricultural operation depending on the weather condition. The influence of both horizontal deformation and soil compaction on the traction performance is very low (does not exceed 4%).


## 3.1. Soft Computing Models

The optimal solution, according to energy savings, is to have high values of both traction force and tractive efficiency. Figure4 shows the dependence of the traction force and tractive efficiency on the vertical load and soil moisture measured on sand, for a horizontal deformation of 0.05 m.Figure 4
3D surface curves of traction force (a) and tractive efficiency (b) as affected by the interactions of vertical load and soil moisture.
(a)
(b)As shown in Figure4, the optimal traction force is produced under a high vertical load, whereas an optimal tractive efficiency can be achieved under a rather low vertical load. Therefore, achieving an optimal balance between the traction force and the tractive efficiency is very difficult and requires accurate mathematical models of the tractive properties.Before the development of the models based on the ANN, linearly dependent predictor variables must be removed from the dataset. Therefore, Pearson’s correlation coefficients between the explanatory variables were calculated. A high positive correlation was observed only between the soil coefficient and the soil moisture (r = 0.767). However, after analyzing the experimental procedures, it was evident that both soil coefficient and soil moisture must be considered as input variables and used for model development. The correlation coefficient is high because the soil moisture range (Table 1) was determined depending on the soil texture.Table3 lists the statistical parameters, namely, R2, RMSE, and MAE, for the best configurations of the mentioned models. The calculations were performed using normalized data (equation (5)).Table 3
Error metrics of best model performances.
Model
Train
Validation
MAE
RMSE
R
2
MAE
RMSE
R
2
Traction force
ANN
0.026
0.037
0.957
0.029
0.040
0.954
ANN + GA
0.064
0.084
0.808
0.063
0.087
0.820
ANN_BP + GA
0.026
0.037
0.958
0.028
0.040
0.955
ANFIS
0.042
0.061
0.892
0.045
0.064
0.888
Tractive efficiency
ANN
0.022
0.030
0.975
0.024
0.037
0.954
ANN + GA
0.076
0.102
0.789
0.068
0.096
0.752
ANN_BP + GA
0.022
0.030
0.975
0.024
0.037
0.954
ANFIS
0.038
0.061
0.883
0.040
0.064
0.872In the case of the traction force model, the best architecture of the ANN was a network with 28 neurons in the hidden layer, and this network was chosen for weights optimization by the GA. The best ANN architecture trained by the GA starting from random values (ANN + GA) was a network with 14 neurons in the hidden layer. The best neural network for the tractive efficiency model had 26 neurons in the hidden layer, and the best ANN + GA model contained 17 neurons in the hidden layer.Figures5 and 6 show the measured and predicted traction force and tractive efficiency values for all the predictive models on the validation dataset. These graphs show the number of data points located very close to the diagonal line, thus facilitating the assessment of the model accuracy.Figure 5
Scatterplot of model-predicted values versus actual values of traction force: (a) ANN, (b) ANN + GA, (c) ANN_BP + GA, and (d) ANFIS.
(a)
(b)
(c)
(d)Figure 6
Scatterplot of model-predicted values versus actual values of tractive efficiency: (a) ANN, (b) ANN + GA, (c) ANN_BP + GA, and (d) ANFIS.
(a)
(b)
(c)
(d)As listed in Table3 and shown in Figures 5 and 6, for both output model parameters, among all the computational models, the ANN and ANN_BP + GA models exhibit the best performance, as indicated by high values of R2 (0.954 and 0.955 for traction force and 0.954 for tractive efficiency) and low values of MAE and RMSE for the validation dataset. The use of the GA for optimizing the weights and biases adjusted by the BP algorithm produced slightly better accuracy in the case of the traction force model. The accuracy of the ANFIS model was lower than those of the ANN and ANN_BP + GA models, with R2 values below 0.9. Additionally, the computational time required for the calculations during the ANFIS model development was significantly higher than that required in the case of models based on MLP. The ANN + GA technique seems to be unsuitable, exhibiting a low accuracy in estimating the traction force and tractive efficiency (R2 = 0.820 and 0.752 for the validation dataset, respectively). Generally, it can be stated that, in agriculture, mathematical models (also based on machine learning) with coefficient of determination (R2) exceeding 0.9 are useful for real life applications [30].Neural networks and hybrid methods were also used by other researchers to model the behavior of agricultural tractors. The ANFIS-based modeling was found to be a promising technique for prognosticating the traction coefficient and tractive power efficiency, withR2 values of 0.98 and 0.97, respectively [17] and for prognosticating the drawbar pull energy of tractor driving wheels, with MSE and R2 values of 0.00236 and 0.995, respectively [16]. In the case of ANN combined with a GA, Taghavifar et al. demonstrated that this method drastically decreased the error and increased the performance of the model of power provided by agricultural tractors as affected by wheel load, slip, and speed [14]. They obtained high values for the coefficient of determination for the ANN + GA model: 0.9696 for the training dataset and 0.9672 for validation dataset.Comparing the current results with those presented by other researchers, it is unclear which technique most accurately models the nonlinear and complex relationships such as the ones investigated in this study. Similar results were obtained by other researchers. Johann et al. compared computational models based on ANN and ANFIS in estimating the soil moisture from the stochastic information of the horizontal and vertical forces acting on a no-till chisel opener using autoregressive error function parameters [31]. The ANN model (R2 = 0.79 and RMSE = 1.27) outperformed the ANFIS model (R2 = 0.69 and RMSE = 1.62) in the test phase. Citakoglu applied ANN and ANFIS for estimating the solar radiation in Turkey using the calendar month number and pertinent meteorological data and obtained a higher accuracy when using the ANN (R2 = 0.930 and RMSE = 1.650) in comparison with using the ANFIS (R2 = 0.926 and RMSE = 1.691) [32]. In contrast, the ANFIS was found to be more suitable than the ANN for estimating the soil cation exchange capacity as affected by clay, silt, sand, organic carbon, and pH in arid rangeland ecosystems and for estimating the oxidation parameters of Kilka oil [33, 34]. Based on relevant error (RE) values, Ping and Fei showed that the accuracy of an ANN combined with a GA (RE = 1.48%) is better than that of a traditional ANN model (RE = 3.91%) for Guangdong port throughput forecasting [35]. Similarly, Srinivasulu and Jain found that the predictive capability of an ANN combined with GA rainfall-runoff models is better than that trained using a BP algorithm [36].


## 3.2. Sensitivity Analysis

A highly accurate mathematical model can give additional information about the relationships under study. A sensitivity analysis has been performed to determine the contribution of independent variables in black box data mining models. The neural network trained by backpropagation combined with a GA was found to be the best model for the relationships analyzed in the present research. This soft computing technique was used for developing the models for the sensitivity analysis. For each dependent variable (traction force and tractive efficiency), a group of twenty ANN models was developed. Table4 lists the parameters of the models.Table 4
Statistics of neural model architectures used for sensitivity analysis.
Dependent variable
The range of number of neurons in hidden layer
The range ofR2 for validation dataset
The range of MAE for validation dataset
The range of RMSE for validation dataset
Traction force
16–40
0.944–0.955
0.028–0.032
0.040–0.045
Tractive efficiency
16–38
0.928–0.954
0.024–0.032
0.037–0.047The results of the relative importance of the input parameters for each model were determined as the arithmetical mean of the results produced by the group of twenty ANN models. The results revealed that the traction force and tractive efficiency are most affected by the soil type (58.3 and 74.5%, respectively). This is in agreement with the results reported by other authors [37, 38]. The two additional parameters that significantly influenced the traction force and tractive efficiency are the vertical load (18.3 and 10.1%, respectively) and soil moisture (19.8 and 10.3%, respectively). The significant effects of these parameters on the tractive performance have been highlighted in other studies as well [11, 39]. It is worth emphasizing that the vertical load is one of the more easily manageable parameters during traction performance optimization. Soil moisture can also be varied in a certain range, as the operator can advance or delay agricultural operation depending on the weather condition. The influence of both horizontal deformation and soil compaction on the traction performance is very low (does not exceed 4%).


## 4. Conclusions

The optimization of the traction performance of an agricultural tractor is essential considering fuel economy. The development of highly accurate mathematical models that describe the tractive properties is an integral part of the optimization process. In this work, four soft computing techniques were used for predicting the traction force and tractive efficiency of a low-power tractor as affected by the soil type (expressed as soil coefficient), vertical load, horizontal deformation, soil compaction, and soil moisture. Comparisons of the error statistics revealed that the neural network model trained by a traditional BP algorithm or by a combination of BP and GA performs better in estimating both the traction force and tractive efficiency than an ANFIS model or an ANN trained by only a GA. An ANN structure with 28 neurons in the hidden layer produced the best model of the traction force, with anR2 value of 0.954, a mean absolute error of 0.029, and an RMSE of 0.040. Similarly, an ANN with 26 neurons in the hidden layer was found to be the best structure for the tractive efficiency model, with R2 = 0.954, MAE = 0.024, and RMSE = 0.037. Using GA for optimizing the weights and biases in the ANN model trained by BP led to a slight improvement in model accuracy. Considering the results presented by other authors, it can be stated that the potential usability of a certain technique depends strongly on the data characteristics. Moreover, the behavior of each machine-learning algorithm is affected by its parameters. Thus, for improving the optimization process, different techniques should be employed, and the model with the highest accuracy should be chosen. Considering the computational time required for ANFIS model development, the neural network trained by the backpropagation algorithm seems to be the best soft computing technique. The results of a sensitivity analysis conducted on a group of models with the highest accuracy showed that the soil type is the parameter most affecting the traction performance of a low-power tractor. A relatively strong influence was also found for the vertical load and soil moisture, which can be varied by the tractor operator to optimize the traction performance.The results of this research are expected to be useful in saving energy in agricultural production systems. However, it should be noted that the application of the empirical models obtained by the authors is limited to conditions similar to those present during the measurements.



---

*Source: 7607545-2020-01-10.xml*