Route Planning for Agricultural Machines with Multiple Depots: Manure Application Case Study

Capacitated field operations involve input/output material flows where there are capacity constraints in the form of a specific load that a vehicle can carry. As such, a specific normal-sized field cannot be covered in one single operation using only one load, and the vehicle needs to get serviced (i.e., refilling) from out-of-field facilities (depot). Although several algorithms have been developed to solve the routing problem of capacitated operations, these algorithms only considered one depot. The general goal of this paper is to develop a route planning tool for agricultural machines with multiple depots. The tool presented consists of two modules: the first one regards the field geometrical representation in which the field is partitioned into tracks and headland passes; the second one regards route optimization that is implemented by the metaheuristic simulated annealing (SA) algorithm. In order to validate the developed tool, a comparison between a well-known route planning approach, namely B-pattern, and the algorithm presented in this study was carried out. The results show that the proposed algorithm outperforms the B-pattern by up to 20.0% in terms of traveled nonworking distance. The applicability of the tool developed was tested in a case study with seven scenarios differing in terms of locations and number of depots. The results of this study illustrated that the location and number of depots significantly affect the total nonworking traversal distance during a field operation.


Introduction
The decreasing marginal income and increasing costs in agriculture require modern agriculture to become increasingly productive [1]. Over the last decades, the production of large and powerful agricultural machinery has been the main engineering focus of manufacturing sectors to increase productivity and efficiency. Although usage of these large machines can increase the unit capacity, it also leads to undesirable and adverse environmental and biological effects such as soil compaction, which potentially reduces the yield [2,3]. Therefore, more efforts are being contributed to the development of advanced information and communication technology (ICT) systems and operation management tools to achieve higher operational efficiency and machinery productivity [4]. These tools/systems have been developed ranging from aiding and supporting navigation efforts to full autosteering systems [5][6][7][8]. Most of these systems can navigate and supervise the operator to complete the field tasks by setting the route for a vehicle either manually or using a predefined fieldwork pattern.
In recent years, a larger number of vehicle route planning algorithms have been developed to be integrated with the aforementioned systems. Commonly, the objectives of these algorithms are twofold: one is creating a geometrical representation of the operational environment in which the field is divided into a set of parallel, straight, or curved tracks and headland passes, and several methods addressing field representation have been developed [9,10]. The other regards the optimization of the routing of vehicles within the geometrically represented field. Concerning this problem, advanced methods based on combinatorial optimization have recently been introduced based on one or multiple optimization criteria that include the minimization of the total nonworking traveled distance [11], total operational time [12], and soil compaction risk [13,14]. All of these algorithms used the metaheuristic algorithms (i.e., Clarke-Wright savings algorithm [15,16], genetic algorithm [17,18], simulated annealing algorithm [6,19], ant colony algorithm [20,21]) to find optimal or near-optimal solutions in a reasonable time. The reason is that the route planning problem is considered as a non-deterministic polynomial (NP-hard) combinatorial optimization problem, and exact algorithms tend to perform poorly in the case of large-size instances [22].
All of the aforementioned algorithms only take into account one depot as the replenishment place during field operations; however, in some operations, multiple depots might be required. For instance, in a capacitated operation such as a manure application, several refillings are needed to complete the entire operation; in a case that the slurry storage facilities are far from the field in order to reduce the total nonworking traveled distance and subsequently increase the field efficiency, farmers often set up multiple slurry buffer tanks (as depots) at different locations for refilling so that each trip can begin and end at different depots. Regarding the multiple depot vehicle route planning problem (MDVRP), it has not received much attention in the agricultural sector. While in the industry sector, this problem has been extensively studied by multiple authors [23][24][25][26]. Due to the specific nature of agricultural operations, existing methods from industry cannot be directly applied. To this end, in this paper, a route planning tool for the MDVRP problem is developed and the optimization problem is solved by using the simulated annealing algorithm. The remainder of this paper is organized as follows. In Section 2, the methodology and mathematical formulation of the problem are explained. The algorithm of the developed tool is validated by comparing it with a well-known route planning method B-pattern [27], and then its applicability is demonstrated in the manure application process considered in Section 3. Finally, in Section 4, the conclusions are drawn.

Problem Definition
Agricultural field operations such as seeding, spraying, and fertilization are executed by a fleet of machines with limited capacity. To achieve the optimal performance of such machines, a route with several trips back and forth between depots needs to be carefully planned. A trip is an operational cycle that consists of suboperations carried out by a machine: 1. filling the manure tank at a depot, 2. leaving that depot to process the field tracks, and 3. driving back to a depot for a new refilling. In addition, as a frequent practice of farmers and machine operators, the distributor can only turn in the headland area in order to reduce the soil compaction and avoid crop damage, and each track is treated only once.
In this study, the manure application process is selected as the study case where multiple refilling depots are considered. The routing problem in our case is formulated as a capacitated vehicle routing problem with multiple depots (MD-CVRP) that has been extensively studied in an industrial section like logistic planning. This problem can be classified into two different categories: fixed and nonfixed destination MD-CVRP. The fixed destination MD-CVRP cannot apply to our research problem, in which each trip begins and ends at the same depot. Therefore, a nonfixed destination MD-CVRP must be developed for this study.

Mathematical Formulation
In this study, all the fieldwork tracks are defined as a set of directed arcs, and each arc has the form of a = v i , v j where v i , v j is called a tail and head vertex, respectively. In this way, each fieldwork track is represented by two arcs with the same length but in opposite directions. These two arcs are called sibling arcs. Therefore, to cover the entire set of tracks, one of these sibling arcs should opt in the optimized solution and a depot is defined as a special arc, which starts and ends at the same point and its sibling is itself. Bochtis and Vougioukas [11] formulated this problem as a graph G = (V, A), , v m+n is the vertex sets corresponding to the ending points of fieldwork tracks. The V t further is classified into two sets: even V te and odd V to vertices. Therefore can be considered as the set of arcs representing tracks, and Therefore, the number of arcs A = A t ∪ A d is equal to m, which is the number of vertices. Each arc a is assigned a value y t , where y t is equal to the total demand of arc a when a ∈ A t , otherwise y t is ∞, meaning that the depot arc has an infinite capacity for the refilling machine unit. The total demand for a field can be calculated by the total summation of values of all arcs' demand divided by two. In addition, a fleet of m homogeneous vehicles of capacity Q is considered in the problem. Each arc a ij has a related traveling cost C ij .
The nonfixed destination MD-CVRP is to consider m vehicle routes with these criteria: • Every path should start and end at the depot; • Every track is covered exactly once; • The total demand for any vehicle path should not exceed the capacity of the vehicle.
For any path, the total nonworking distance can be defined as f (∂), which is the summation of distances from the starting depot to the first fieldwork track in the planned track sequence and all headland turning distances, and finally the movement from the last track to the final depot. The optimal solution which minimizes the cost function is the permutation ∂ * and can be calculated as follows: The decision variable X ij of the formula is defined as follows: X ij = 1 : (i f the vehicle immediately goes f rom vertex i to vertex j) ; 0 : (i f the vehicles need to visit depots f or re f illing) The objective function can be defined as follows: Min : The constraints are as follows: i∈H j∈H Agronomy 2020, 10, 1608 4 of 14 Equation (4) states that only one of the sibling arcs needs to part of the solution. Equations (5) and (6) ensure that each route should be started and finished at a depot. Equation (7) specifies that the total summation of demand related to each trip should not exceed the capacity of the distributor's tank. Equation (8) represents that if the vehicle starts covering a track K, it will leave the track at the end of the process. Equation (9) will remove the subtours from feasible solutions. The last equation, (10), states the type of the variable.

Methodology
The problem explored in this study is to find the optimal traversal sequence of tracks in a field with multiple depots by considering the capacity of the distributor to minimize the total nonworking distance. Therefore, this problem can be classified as a well-known multidepot capacitated vehicle routing problem (MD-CVRP), which is an NP-hard problem and the tracks' arcs are going to visit instead of customers, and to complete the task, the distributor should visit the depots several times for refilling. A metaheuristic algorithm, simulated annealing (SA), was selected for finding the solution for this problem.

Field Representation
Field representation is a process of generating a geometrical representation of a field area to provide a concise representation of the operational environment, which includes three tasks: generation of headland passes (H), determination of the fieldwork tracks (T) that completely cover the field area, and generation of the connecting paths that connect the depots, tracks, and headland passes. Specifically, these connecting paths that demonstrated in Figure 1 are:

•
A path connects fieldwork tracks' ends to adjacent fieldwork tracks' ends (T2T); • A path connecting each depot and to the first headland pass (D2H); • A path connecting each headland pass to its adjacent headland pass (H2H); • A path connecting fieldwork tracks ends at each headland pass (T2H). It is worth noting that the curvature of connecting paths is constrained by the minimum turning radius of the machine for making the paths traversable by machines.

Cost Matrix Generation
The cost matrix gives the transition cost between every pair of vertices of the field graph. The field graph consists of all the edges in the field with their corresponding weight (which is the actual length of the edge). The cost matrix is square (n * n ) with an element C ij which represents the transition cost from vertex i to vertex j, where i j; otherwise, the cost is equal to 0. The shortest path search algorithm was applied for calculating the transition cost between every pair of nodes [28].

Simulated Annealing Algorithm
Simulated annealing (SA) as a stochastic algorithm is used in this study to investigate the near-optimal solution. This method has been implemented in several problems in various fields such as computer design, route planning, and image processing [6,28,29]. The mechanism of this algorithm is based on the simulation of a cooling process called annealing where a solid is gradually cooled. The SA algorithm has some primary features such as the initial temperature, cooling rate, and termination policy, which should be defined carefully in advance since the performance of the algorithm depends on them. For instance, at the beginning when we have a higher temperature, the probability of accepting newer worse solutions is high and with time it will decrease. By this ability, the SA algorithm enables a search of all solution space and an escape from trapping in local minima [29].
All of the SA algorithm steps can be listed as follows: 1.
Define the SA main parameters such as initial temperature = 200, temperature-damping rate = 0.9.

2.
SA main loop a. Generate initial solution I.
Generate the initial solution based on the maximum capacity of the vehicle to determine when it is needed to go to the depot for refilling II.
Calculate the cost of an initial solution by using the cost function that represents the total nonworking distances traveled by the vehicle b. Set the temperature (T) equal to the initial temperature (T 0 ) c.
Update the best solution ever found If the new solution is convincing, then accept it as the best solution d.
If the new solution is not convincing, then there is an opportunity to accept it by the condition based on T e.
Update the best solution f.
Reduce the T based on the damping rate (α).
A flowchart (Figure 2) summarizes all the steps of the SA algorithm:

Initial Solution Generation
Before the initial solution construction, all arcs are clustered based on the distance of its head to its nearest depot. For example, as presented in Figure 3, the nearest depot to arc 3 is depot 1, while arc 4 is depot 2. The construction criteria are: • Each solution should be initiated and ended with a depot that can be the same or different.

•
For simplicity, only one of each sibling arcs is presented in the solution.

•
Based on the solution already constructed, if the remaining capacity of the agricultural machine is not sufficient to process the next track, then the nearest depot has to be selected for refilling service.

Solution Evaluation
In order to evaluate the quality of a solution, a cost function is defined, specifically for the traveling cost. Let C(a i ) be the traveling cost for covering the arc a i . Then for a solution, the total traversal cost of the agricultural vehicle is equal to the summation of the costs C(a i ) corresponding to the arc a i , which is traversed by the vehicle. The total cost function F(a) can be defined as follows: F(a) = i∈(visited arcs list) C(a i ).

Neighborhood Solution Generation
Random swaps, random swaps of subsequences, random insertions, random insertions of subsequences, reversing a subsequence, random swaps of reversed subsequences, and finally random insertions of reversed subsequences are among the neighborhood operators that can be used for making a new solution in the SA algorithm. Table 1 shows how these operators make changes in the initial solution. Table 1. Applied neighborhood operators [6].

Results and Discussion
The SA algorithm was implemented in MATLAB ® technical programming language (The MathWorks Inc., Natwick, MA, USA) on a computer with the following specifications: Intel ® Core ™ i5-3210 M CPU @ 2.50GHZ with 4GB RAM. To investigate the quality of the solutions achieved by the algorithm developed, the results presented by the well-known approach (B-pattern) [11] are used for comparison. Furthermore, [29] also compared their proposed algorithm with the B-pattern approach. The field shapes and machine-related specifications were considered in this work as shown in Table 2, which are used as input data for the optimization algorithm. The capacity of the machine is set to infinite because disk-harrowing was considered as the field operation in [11]. In the B-pattern approach [11], the authors apply their proposed method in four fields. Field numbers 1 to 4 are divided into 8, 20, 12, and 12 tracks, based on the defined operating width of the machine and the size of those fields. The results achieved from the comparison of the proposed optimization algorithm in this study with the solutions provided in [11] and [29] are presented in Table 3. The comparisons illustrate that all the results from the algorithm presented in this study outperformed the solutions provided by [11]. The total nonworking traveled distance in the headland part for the first field, based on the presented method in this study, is 94.439 m, which shows a 1.4% improvement in comparison with the quality of the solution found in [11]. Moreover, the total nonworking traveled distance in the headland area for field numbers 2 to 4, based on the proposed approach in this work, is 232.566, 141.027, and 141.027 m, respectively. The comparison of the results achieved in this work with the solutions presented in [11] shows that our solutions are 1.4%, 20%, and 15.3% more efficient than those presented in [11]. Moreover, the comparison shows that the solutions presented in this work are 1.2%, 3.4%, and 3.1% more efficient than those provided in [29] for field numbers 2, 3, and 4, respectively. Table 3. Comparison of the generated solutions based on the proposed method in this study with the results provided in studies [11] and [29]. (A) Method proposed in [11], (B) method proposed in [29], and (C) method proposed in this study.  Figure 5 shows the coverage of all four mentioned fields with the type of turnings used for traveling from one track to another.  [11], (B) method proposed in [29], and (C) method proposed in this study. The track's index starts from the left and moves to the right in each field.

Proposed
To demonstrate the applicability of this method for the multiple-depot route planning problem, seven scenarios for the manure application process were carried out based on the number of the depots and their locations in the field, as well as the distributor specifications. Specifically, the tested field is selected from the GIS database of the Danish Ministry of Food, Agriculture, and Fisheries, and the information regarding the properties of the field is presented in Table 4.  Figure 6 illustrates the location of possible depots and field partitioning. The technical characteristics of the distributor are presented in Table 5. In this regard, the field can be divided into 25 fieldwork rows. The total demand of the field is equal to 4235 (liter), therefore at least five refilling trips between the field and depot are needed to cover the entire field. A total of seven scenarios are considered where they only differ in the locations and number of depots used.
Agronomy 2020, 10, x; doi: FOR PEER REVIEW www.mdpi.com/journal/agronomy Depot 3 56°36′25.046″ N 10°13′45.1398″ E Figure 6 illustrates the location of possible depots and field partitioning. The technical characteristics of the distributor are presented in Table 5. In this regard, the field can be divided into 25 fieldwork rows. The total demand of the field is equal to 4235 (liter), therefore at least five refilling trips between the field and depot are needed to cover the entire field. A total of seven scenarios are considered where they only differ in the locations and number of depots used.   Table 6 shows the list of tracks generated after partitioning the sample field. Each track has a pair of sibling arcs and in the provided solution for simplicity, only one of these sibling arcs is demonstrated. For instance, arc 1 starts from node 1 and ends at node 2 and arc 2 is from node 2 to node 1.  Table 5. Technical characteristics of the applied manure distributor.

Technical Characteristics of the Distributor Values
Capacity (m 3 ) 25 Minimum turning radius (meter) 9 Operating width (meter) 12 Application rate (m 3 /ha) 17 Table 6 shows the list of tracks generated after partitioning the sample field. Each track has a pair of sibling arcs and in the provided solution for simplicity, only one of these sibling arcs is demonstrated. For instance, arc 1 starts from node 1 and ends at node 2 and arc 2 is from node 2 to node 1.  Table 7 illustrates all the scenarios with their corresponding solution, including the optimal order of arcs to cover. To better compare the defined scenarios, the total traversal nonworking distance in the field is calculated.
The results in Table 7 show that the locations of the depots affect the nonworking traversal distance in the field. By comparing the three scenarios (1,2,3) in which only one depot is considered, it can be seen that scenario 3 produces the best solution. Scenario 3 with 3220.9 (m) nonworking traversal distance has 35.6% and 40% cost savings in comparison with scenarios 1 and 2, respectively. Furthermore, a different combination of depots in scenarios 4 to 6 illustrates how the position and the order of visiting depots can affect the traversal distance. As an example, scenario 6 shows approximately 8% cost saving in comparison with scenario 4. The scenario evaluation indicates that the number of depots can also influence the traversal distance. For instance, scenario 6 with a nonworking distance of 2745.8 (m) indicates 48.9% cost saving in comparison with scenario 2. Moreover, scenario 7 with three depots shows 11.4% and 50.9% cost saving in comparison with scenarios 4 and 2, respectively. Furthermore, the comparison of running time for all the scenarios shows that by increasing the number of depots, the complexity of the problem will increase and as a result, the algorithm takes more time for calculation.
The number of required depots for a field can be determined as follows. First, we need to know the rental/buying price of one slurry buffer tank (depot). Second, the amount of fuel consumption of the vehicle is required. The cost of fuel consumption for the nonworking distance traveled in each scenario can be calculated by multiplying the cost of fuel by the fuel consumption for nonworking traveled distance. The total cost for each scenario is the summation of fuel consumption cost and the price of slurry buffer tanks (depots). The location and the number of depots can affect the nonworking traveled distance and the total cost related to each scenario. The comparison of the total cost for each scenario can show the best scenario.
Since the proposed tool can reduce the in-field traveled distance and traffic intensity, it potentially can reduce the soil compaction risk in the field. Moreover, other factors also influence soil compaction, such as increasing soil organic matter [30,31], adequate crop rotation [2,32], types of soils [33], and contact pressure [14], which should be considered in assessing the soil compaction risk in the field.
Generally speaking, the application of the route planning tools in farm operations is beneficial currently due to their potential ability to reduce the cost units and increase the operational efficiency in the field. As one advantage of the tool presented in this study, restricting all the nonworking traversal distance to the headlands area makes this tool able to control the traffic in the main cropping area of the field.
The novel ICT systems in agricultural fleet management [4,34] can facilitate the implementation of this tool for navigating the vehicles in the field. Especially, the applicability of the proposed tool depends on the integration of in-field navigation technology and navigation aid systems. The optimized sequence of tracks generated by this tool can be used to navigate an agricultural vehicle in the field. In this regard, the operator of the machine can follow the turning paths generated by the field representation module in the headland area.
In the vehicle routing problem, the optimization is accomplished based on the cost matrix. In this paper, the nonworking traversal distance was used to generate the cost matrix. Potentially, the tool presented could be implemented in the fields with obstacles. Several field representation algorithms were developed by other researchers [9,10], which considered obstacles in the field. In such tools, the cost matrix is generated based on the field representation module, and therefore the occurrence of obstacles in the field could simply be reflected in the cost matrix. This approach could also apply to fields with slopes by applying a 3D field representation algorithm [35]. To implement this tool in MVRP (multiple vehicle route planning) problems when heterogeneous vehicles with different capacity or maneuver capabilities are present, an individual cost matrix should be considered for each vehicle.
For future studies, the tool provided can be modified to handle other types of VRP problems. For instance, in the SVRP (stochastic vehicle route planning) problem, the proposed tool can reflect the uncertainty by stochastic costs in the cost matrix [36,37]. Furthermore, the real-time calculation of deviations between the expected and real progress of route planning can be added to this tool, which is useful for illustrating the precision of the applied tool.

Conclusions
In this paper, a tool for optimal route planning for agricultural machines with multiple depots was developed. The optimization model developed was validated by comparing it to the well-known route planning approach, B-pattern. The comparison illustrates that the proposed algorithm outperformed the B-pattern approach up to 20% in terms of reducing the nonworking traveled distance in the field. Moreover, the applicability of the presented tool for multiple depots vehicle routing problem was tested in a case study. Specifically, seven scenarios were defined to consider the different combinations of multiple depots to study the effects of the quantity and the location of the depots on the nonworking traversal distances in the field. The scenario comparisons demonstrate that changing the location of the depots could make up to 40% and increasing the number of depots could bring up to 50.9% reduction in the traveled nonworking distance in the field. The tool developed can be used as part of a decision support tool to provide optimized route planning for field operations by considering field features and machinery specifications. Since this tool can generate curved paths in the field representation module, it can be used for autonomous agricultural vehicles as well as robotics for complete field coverage.