Development of a responsive optimisation framework for decision-making in precision agriculture

Journal Pre-proof

Development of a responsive optimisation framework for decision-making in precision agriculture Qingyuan Kong , Kamal Kuriyan , Nilay Shah , Miao Guo PII: DOI: Reference:

S0098-1354(19)30595-2 https://doi.org/10.1016/j.compchemeng.2019.106585 CACE 106585

To appear in:

Computers and Chemical Engineering

Received date: Revised date: Accepted date:

4 June 2019 25 August 2019 30 September 2019

Please cite this article as: Qingyuan Kong , Kamal Kuriyan , Nilay Shah , Miao Guo , Development of a responsive optimisation framework for decision-making in precision agriculture, Computers and Chemical Engineering (2019), doi: https://doi.org/10.1016/j.compchemeng.2019.106585

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Ltd.

Highlights     

Responsive optimisation framework developed for precision agriculture Mathematical optimisation to bring real-time data into precision decision-making Meta-heuristic algorithms proposed to reduce computational time by over 95% Hybrid metaheuristic algorithm for trade-off between responsiveness and optimality Case study to demonstrate precision farming for sugarcane harvesting in South Africa

1

Development of a responsive optimisation framework for decision-making in precision agriculture Qingyuan Kong1, Kamal Kuriyan1, Nilay Shah1, Miao Guo1* 1Center

for Process Systems Engineering, Department of Chemical Engineering, Imperial College London, SW7 2AZ Corresponding: [email protected]

Abstract Emerging digital technologies and data advances (e.g. smart machinery, remote sensing) not only enable Agriculture 4.0 to envisage interconnected agro-ecosystems and precision agriculture but also demand responsive decision-making. This study presents a mathematical optimisation model to bring real-time data and information to precision decision-support and to optimise short-term farming operation. To achieve responsive decision-support, we proposed two meta-heuristic algorithms i.e. a tailored genetic algorithm and a hybrid genetic-tabu search algorithm for solving the deterministic optimisation. The developed responsive optimisation framework has been applied to a hypothetical case study to optimise sugarcane harvesting in the KwaZulu Natal region in South Africa. In comparison with the optimal solutions derived from the exact algorithm, the proposed meta-heuristic methods lead to near optimal solutions (less than 5% from optimality) and significantly reduced computational time by over 95%. Our results suggest that the tailored genetic algorithm enables rapid solution searching but the solution quality on sugarcane harvesting cannot compete with the exact method. The hybrid genetic-tabu search algorithm achieved a good trade-off between computational time reduction and solution optimality, demonstrating the potential to enhance responsive decision making in precision sugarcane farming. Our research highlights the development of the responsive optimisation framework combining mixed integer linear programming and hybrid metaheuristic search algorithms and its applications in real-time decision-making under Agriculture 4.0 vision. Keywords: event-based optimisation, MILP, heuristic optimisation, precision agriculture, responsive decision-making, harvesting

1. Introduction Industry 4.0 ("Recommendations for implementing the strategic initiative INDUSTRIE 4.0," 2013) envisages industries and supply chains where resources, products and machines are connected via the internet to achieve smart data analyses and coordinated processes. An Industry 4.0 co-concept in the agriculture sector has emerged, which is termed as Agriculture 4.0. Building on digital technologies and data advances (Internet of Things (IoT), sensor data and remote sensing technologies), Agriculture 4.0 has the potential to enhance precision farming and improve farm system responsive performance and precise decision-making in response to operational uncertainties and real-time data updates. As the part of IoT driven precision solutions, agricultural machinery is expected to catalyse the agricultural sector transformation by shifting isolated farms towards performance-optimised and vehicleinterconnected agro-ecosystems. Across such agro-ecosystems, different decision levels can be 2

grouped as long/mid-term planning and short-term operation and scheduling based on temporal scales. Despite the importance, development and application of real-time decisionmaking tools to optimise precision planning and operation of interconnected smart machinery remain largely unexplored. This study will demonstrate optimisation research under the Agriculture 4.0 vision with a particular focus on harvest planning and operation for a perennial crop, sugarcane. Sugarcane is one of the most important commodities globally (Jena & Poggi, 2013); its production has increased fourfold since 1965, reaching 1842 million metric tons in 2017(Statistica, 2018). Sugarcane, widely cultivated in tropical and subtropical countries in Africa, Asia and Latin America and the Caribbean (Nations., 2018), not only dominates global sugar production (86%) but is also a primary bioethanol feedstock, contributing to approximately 30% of global ethanol supply (30 billion litre) (Nations., 2018). BRICS (refers to five emerging national economies i.e. Brazil, Russian Federation, India, China and South Africa) countries representing the growing economy also play significant roles in the global sugarcane market and along with EU nations dominate the global sugar market. Sugarcane contributes substantially to national development in BRICS countries e.g. GDP in South Africa (Farzad et al., 2017) and Brazil. Sugarcane represents a highly perishable agro-product due to cane quality deterioration which can be affected by variety, harvesting, and pre- or post-harvest operations (Solomon, 2009). Thus sugarcane farming increases the system complexity and requires coordinated planning and scheduling (Kusumastuti et al., 2016). As highlighted by Florentino et al. (2018), decision support to optimise operational planning and to estimate the harvest process quality will be of benefit to the sugarcane and wider agricultural sectors. Sugarcane harvesting has evolved from manual operation to mechanical harvesting over the past decades(Ma et al., 2014). A few studies have been undertaken to improve the efficiency and reduce the cost involved in the sugarcane harvesting process. They can be broadly classified into three levels, namely the planning level, short-term operational level and in-field routing optimisation problems. The planning horizon typically spreads over an entire harvesting period, which can range from one to several months (9 months in South Africa)(Ramburan, 2011). At this level, the primary goal for farmers is to harvest each field as close as possible to its peak maturation period to maximise the sugar outputs, while simultaneously satisfying the demands for sugarcane and minimising the usage of harvesters. For instance, Florentino et al. (2018) developed a multiobjective mixed integer linear programming model that minimises the time deviation between the ideal and the actual harvesting time of each field, while maintaining the number of employed harvesters to a minimum. A large-scale integer programming model was proposed by A. J. Higgins (1999) to optimise the decisions on harvest date, crop cycle length, and whether to fallow for all paddocks within a mill region. Stray et al. (2012) built an optimisation-based decision making tool to inform the scheduling tasks of sugarcane harvesting operations in South Africa. da Silva et al. (2015) introduced a multi-choice goal programming (RMCGP-LHS) model to deal with uncertainty in harvest scheduling for the sugar and ethanol milling industry with an objective to harvest sugarcane plots in the time window of the highest sucrose levels, while also minimising agro-industrial costs. Ramos (2016) developed a nonlinear optimisation model to plan sugarcane planting and harvesting over a period of five years, with the objective to optimise overall sugarcane production.

3

At the second level, which has a relatively shorter time window compared with the planning stage, key decisions such as the harvesting sequences and schedules are made. Karlsson et al. (2003) devised a superstructure of harvesting sequences and road networks, from which a mixed integer programming (MIP) model is applied to determine the optimal sequence. Sethanan and Neungmatcha (2016) focused on the path planning of the harvester, involving the constraints on moving direction and field accessibility to simultaneously minimise harvester moving distance and maximise the sugarcane yields. Junqueira and Morabito (2017) developed a mathematical optimisation approach to support decision-making on the scheduling of harvest operations. The work proposed by Jena and Poggi (2013) presents a generalised MIP model which can be adapted for both planning and operational stages of sugarcane harvesting. Another aspect in agricultural operations is the in-field operation of machinery (Zhou et al., 2014), where field geometry has a significant impact on the efficiency of agricultural machines. For the harvesting of sugarcane, Santoro et al. (2017) proposed an integer programming model which minimizes the manoeuvring time of the harvesting machine and, consequently, reduces operational costs e.g. fuel and labour costs. However, almost all the mathematical models developed for the planning and scheduling of sugarcane harvesting are based on discretised time intervals, which only provide a reference grid of time for resource-competing operations (e.g. harvesting which competes on shared harvesters). Although various constraints in the scheduling problems can be formulated in a relatively straightforward manner, the discrete-time formulations are by definition only approximations of the actual problems (Floudas & Lin, 2004). Furthermore, one of the key issues in such approaches is the arbitrary discretisation of the time horizon with uniform time grid intervals, which presents a trade-off between the solution quality and the computational time. For harvesting operations, the varying geometries of different fields lead to variations in harvesting time; however, the discrete-time approaches can only provide approximate descriptions of the actual process, which can deviate substantially from the true solutions (Floudas & Lin, 2004). The real-time decision making and accurate optimisation solutions over continuous time horizon are of particular importance for coordinated operations and interconnected smart machinery with real-time communication (Agriculture 4.0 vision). Due to the inherent limitations of the discrete-time approaches, significant research efforts have been devoted to develop continuous-time representations, notably in the field of Process Systems Engineering (PSE) for process scheduling (Floudas & Lin, 2004). However, this research challenge remains open in the field of agriculture planning and scheduling e.g. sugarcane harvesting. Limited previous research has been published on the application of continuous time formulations to sugarcane harvesting and processing which has a particular focus on minimising the cut-to-crush times and the number of trucks required in the operation (A. Higgins, 2006; Lamsal et al., 2016a, 2016b). A modelling gap exists for continuous-time algorithm development to provide more accurate optimisation solutions for the operation and scheduling of sugarcane harvesting, which are operated in a coordinated manner through smart harvesters. Inspired by the optimisation advances in the PSE field (Floudas & Lin, 2004; Kondili et al., 1993), this study presents a harvest scheduling optimisation model based on the eventbased continuous-time algorithms. Another challenge lies in responsive decision-making based on deterministic optimisation methods. Different solution searching methods have been proposed, including exact algorithms, heuristic and metaheuristic algorithms. Exact approaches to derive optimal solutions mainly 4

adopt tree-search algorithms such as branch and bound and branch and cut algorithms (Belotti et al., 2013). However, finding global optimal solutions using exact algorithms can be challenging in particular for large-scale optimisation problems, which require responsive decision-making as the computational time increases exponentially with the optimisation problem complexity. Alternatively, heuristics and metaheuristics may offer faster responsive approaches in generating near-optimal solutions. The most widely adopted (meta) heuristic algorithms to search solution spaces include genetic algorithms, differential evolution, ant colony optimization, particle swarm optimization (PSO), simulated annealing and tabu search or hybrid approaches e.g. genetic-Tabu hybrid algorithms (Y. P. Zhang et al., 2019). Genetic algorithms are often employed for their search ability however, their performance may be hindered by “premature convergence” and the poor hill-climbing performance (Y. P. Zhang et al., 2019). Tabu search (F. Glover, 1990) methods implement an efficient algorithm with hillclimbing ability in solution searching; however, their searching capability may be limited by the quality of the initial solution. Hybrid approach enables algorithms to collectively and cooperatively solve a predefined problem. Thus hybrid metaheuristic algorithms have the potential to enable responsive decisions and reduce the computational time of large-scale deterministic optimisation. Studies attempted to classify the taxonomy of hybrid metaheuristic algorithms (Blum et al., 2011; O Ting et al., 2015; Parejo et al., 2012; Talbi, 2002) but such classification efforts were hindered by the fact that these algorithms possess overlapping characteristics (Lopes Silva et al., 2018). With the increasing interests in hybridisation of metaheuristics, a number of algorithms have been developed such as genetic algorithm hybrids (Al Chami et al., 2018; Liu & Kozan, 2016; Nouri et al., 2016; Yanik et al., 2014), PSO hybrids (Duan et al., 2013; Ebadi & Navimipour, 2019; Gao et al., 2014; Sahin & Kellegoz, 2019; Vahdani et al., 2012). They have been applied to solve large-scale complex optimisation problems e.g. assembly line optimisation (Şahin & Kellegöz, 2019),real-time task assignment problem (Marimuthu et al., 2018), closed-loop supply chain design(Devika et al., 2014), electrical vehicle charging scheduling (Schneider et al., 2014). Genetic-tabu search approaches were proposed in operational research such as applications in job shop scheduling optimisation (Q. Zhang et al., 2012), vehicle routing with time windows constraint (Y. P. Zhang et al., 2019). Hybrid approaches based on exact, genetic, tabu search algorithms were also presented to solve bi-level optimisation problem to design efficient transportation systems (Parvasi et al., 2017). However, hybrid metaheuristic application in agricultural operational and scheduling problems remains largely unexplored. This study explores different solution searching methods and compares exact algorithms and hybrid metaheuristic approaches in solving harvest scheduling optimisation problems. This study proposes an event-based continuous-time formulation to optimise short-term sugarcane harvesting operations. To demonstrate the model functionality, we present a hypothetical case study based on sugarcane operation using smart harvesters in the Kwazulu Natal region in South Africa, where the optimisation solution approaches will be compared.

2. Problem Statement The harvesting optimisation problem can be stated as follows: given a set of sugarcane fields with known maturity curves, a set of harvesting crews and harvesting windows for each field, to determine the schedules of each harvesting crew that maximise the total profits of harvesting operations. All fields are ready for harvesting; whereas each field has a fixed time horizon

5

within which harvesting must be completed. The number of harvesting crews is given, which bounds the resource availability for harvesting operations.

3. Model Formulation A mixed integer linear programming (MILP) optimisation model has been formulated, which is presented in this section.

3.1.

Nomenclature

Sets

Parameters

Definitions Sets of harvesting crews Sets of all fields The warehouse ‘field’ Set of all event points Definitions The operation time horizon (hr) A large number representing the harvesting capacity (ha) Harvesting speed expressed as the amount of sugarcane harvested per unit time (ha/hr) Unit time per unit sugarcane harvested (hr/ha) Initial amount of sugarcane ready for harvesting in each field (ha) Lower bound of the hard harvesting window for field (hr) Upper bound of the hard harvesting window for field (hr) A large number representing the maximum allowed idle time (hr) Minimum idle time (hr) Time needed for harvesters to travel from one field to another (hr) A large arbitrary constant (Big M) Lower bound of the soft harvesting window for field (hr) Upper bound of the soft harvesting window for field (hr) Penalty cost coefficient for violating the ideal harvesting window($/hr) Revenue generated from selling a unit of harvested sugarcane ($/ha) Cost coefficient for harvesting per unit sugarcane ($/ha) Cost coefficient for harvester movement ($/hr) Penalty cost coefficient for the idle time of harvester crews($/hr)

Variables

Definitions The objective value

Binary Variables if crew starts travel from field to

at event point

if crew is at field at event point if crew starts harvesting in field at event point if crew becomes idle in field at event point 6

Positive Variables The start time of event of crew The finish time of event of crew Amount of sugarcane harvested by crew in field at event point Idle time for each idle event

of crew

Total revenue generated by the harvested sugarcane Total harvesting cost Total costs for harvester movement between fields Total penalty costs for harvester being idle Total penalty costs for harvesting operation outside the ideal harvesting window

3.2.

Objective Function

As given in Eq. (1), the objective function is to maximise the overall profit, which is determined by : 1) the revenue from harvested sugarcane ( ) defined by Eq.(2); 2) the operational costs induced by harvesting ( -Eq.(3)) ; 3) the cost incurred by harvesters travelling between fields ( -Eq.(4)); 4) the penalty costs for harvester idle time in fields ( - Eq.(5); 5) the penalty costs for operation outside the ideal harvesting window ( -Eq.(6)), and this is discussed in greater detail in section 3.7. (

) ∑

(2)

∑

(3)

∑

(4)

∑

(5) ∑

3.3.

(1)

(6)

Crew Location Constraints

As defined in Eq. (7), each harvester crew can only be present in one field at each event point. ∑

(7)

At event point , all harvesters are located in the warehouse, which is defined as the warehouse field . The harvester crew that is present in each field at the beginning of each event point follows the ‘resource balance’ equality constraint Eq. (8): ∑

(

)

(8)

7

where at the event point , the status of crew at field is determined by its status at the previous event point ( ) and its traveling pattern in the previous event (i.e. the harvester crew travelling to or leaving field j). The movement of each harvester is subject to the following constraints (Eq. (9) and Eq. (10)): ∑

(9) ∑

(10)

Eq. (9) constrains the maximum number of movements between different fields for each harvester crew over the entire harvesting period, i.e. harvester crew can only move from field to field up to once during the entire harvesting time window. Eq. (10) defines the upper bound for the total number of movements for crew at each event point . Thereby, within each event and field, a crew is limited to three operational options: harvesting, departure from the field or being idle. These three options are mutually exclusive and are subject to the existence of the crew in the field, which can be expressed as Eq. (11). ∑

3.4.

(11)

Harvesting Constraints

Harvesting is assumed to be conducted by each crew in a continuous fashion, which means that each crew only passes a given field once throughout the harvesting season and does not leave the field until its scheduled harvesting operation is completed. The mathematical expression is given in Eq. (12). ∑

(12)

The amount of sugarcane harvested in field by each crew at event point crew harvesting capacity (Eq. (13)).

is bounded by the

(13) The total amount of sugarcane harvested from a given field during the entire harvesting season must not exceed the maximum sugarcane (ha) cultivated in the field (Eq.(14)). ∑

3.5.

(14)

Crew Idle Time Constraint

The constraint for the idle time (

) of each idle event is defined by Eq. (15). (15)

3.6.

Time and Sequence Constraints

A start time and a finish time are assigned to each event point of for each crew and the finish time is equal to the start time plus the event duration. The event durations, as given in Eq. (16), are determined by three variables namely the time spent by crew in moving between 8

fields (∑

), the harvesting time of each crew (∑

) and the idle time of

each crew (∑ ). Since events are mutually exclusive within a given field for each crew, only one of the abovementioned variables would be non-zero and contribute to the finish time of each event. ∑

∑

∑

(16)

A time sequencing constraint Eq. (17) is introduced to ensure the start time of an event ( is equal to the finish time of the previous event .

)

(17) For a given field with a strict harvesting time window (i.e. hard time window), the finish time of harvesting operation must be within the upper bound of the time window (Eq. (18)) and the start time is constrained by the lower bound of the strict harvesting time window (Eq. (19)). ∑

∑

∑

∑ Eq. (18) is a redundant constraint if

3.7.

(18) (19)

or

is equal to 1 and it is active only if

.

Ideal Harvesting Period

Apart from the hard time window defined above, the concept of soft harvesting window is introduced for each field, defined as the optimal time period at which the peak maturation of the sugarcane is reached. The soft time window for each field is used to quantify the amount of sugar loss due to the harvesting outside this window, which is penalised in the objective function. To assess the penalty term, the relative positions of the harvesting window and the actual harvesting time are categorised into six different scenarios illustrated in Figure 1. A unified equation is used to calculate the time of each harvesting event that lies outside the ideal harvesting window ( ): (20) where the expressions on the right hand side of the equation are defined in Table 1. Table 1: the expression on the right hand side of equation (20)

The actual harvesting duration of crew at event The absolute time difference between the finish time of the harvesting event and the upper bound of the ideal harvesting window (only active if ) – (Figure 1d) The absolute time difference between the start time of the harvesting event and the lower bound of the ideal harvesting window (only active if ) – (Figure 1c) The duration of the ideal time window of each harvesting event The absolute time difference between the finish time of the harvesting event and the lower bound of the ideal harvesting window (only active if ) (Figure 1e) The absolute time difference between the start time of the harvesting event and the 9

upper bound of the ideal harvesting window (only active if

) (Figure 1f)

Figure 1: six scenarios for the relative positions of the actual and the ideal harvesting windows for each harvesting event

As defined in Eq. (20), the variable enables the objective function (Eq. (1)) and Eq.(6) to account for the loss of sugar content (caused by the harvesting outside ideal time window). The right hand side of Eq. (20) varies with six scenarios given in Figure 1, which are formulated using the generalised disjunctive programming approach and are summarised in Table 2. Table 2: the value/representation of right hand side of Eq. (20).

a)

b)

c)

d)

e)

f)

∑

∑

∑

∑

∑

∑

10

As defined in Table 3, for a harvesting event, (

) binary variables are used to determine

which scenarios are enforced. Given the binary variables, the disjunctive constraints are formulated in Eq. (21). Table 3: Definition of the binary variables used to categorise the six scenarios

( ( ( (

;0

)

if otherwise

;0

)

if otherwise

;0

)

if otherwise

;0

)

if othersise

(21)

3.8.

Inequalities based on Assumptions about Harvesting Operations

In this model, it is assumed that each harvesting crew’s relocation is always followed by an idle or a harvesting event, and this leads to inequality constraint given in Eq. (22). ∑

(22)

Additionally, it is assumed that if the starting time of the ideal harvesting window of field is earlier than that of the field ( )and both of which are harvested by the same crew , the event

of harvesting must not occur after the event {

3.9.

|

of harvesting }

(Eq. (23)). (23)

Other Constraints to initialise operational time horizon

Upper bounds have been introduced for the starting and finish times of each event (Eq. (24) and Eq. (25)). (24) (25)

11

The first event ( ) of each crew is initialised at ( ( crew must end with the end of the operational time horizon Eq. (25).

)), whereas the last event of each and is already constrained by

(26) As formulated in Eq. (27) and Eq. (28), at the beginning and the end ( ) of each harvesting time horizon, the location of all harvester crews are limited to the warehouse ( ). Eq. (29) constrains the harvesting event, which is not allowed at the warehouse. (27) (28) (29) Finally, movement within the same field is not considered as an event (Eq. (30)) (30) In summary, by implementing constraints from Eq. (2) to Eq. (30), the optimisation problem is formulated as a MILP model with an objective function to maximise the operational profit given in Eq. (1).

4. Solution Strategy The mathematical model presented above is formulated in Python 3.7 using Pyomo 5.6.2 and is solved using CPLEX 12.8.0 on a 12-Core Intel(R) Xeon(R) CPU E5-2697 v2 @ 2.70GHz with 96GB of memory. The exact solution of the operation schedule for a given case study is presented in section 5. However, with recent data advances, responsive decision-making in response to real-time data collection and analyses is essential, which requires a trade-off between computational time and solution quality. Metaheuristic approaches such as genetic algorithm and tabu search offer promising solution searching methods in that regard. Genetic algorithm delivers good performance in generating feasible solutions and searching for optimal solutions within significantly reduced computational time for large scale optimisation problems. However, genetic algorithm performance is hindered by its poor hill-climbing ability, which inclines to lead to local optimal solutions in particular for complex optimisation problems. In contrast, tabu search algorithm, which features hill-climbing capability in solution searching by moving to immediate neighbours to find better solutions, enables escaping of the local optimality. Tabu search performance is limited by its initial solutions. A genetic-tabu search hybrid solution would combine both capabilities and overcome performance limitations. Therefore, by adapting the original mathematical model, a heuristics-based model is built and solved first using genetic algorithm, which is detailed in section 4.1 and further refined using tabu search given in section 4.2.

4.1.

Genetic Algorithm

Genetic algorithm (GA), is a type of evolutionary meta-heuristic algorithm representing an intelligent exploitation of a random search applied to solve optimisation problems. This section presents the GA application in the sequencing and scheduling of harvesting crews. 4.1.1. Solution Representation

12

A solution for the harvesting problem is treated as a chromosome represented as an array such as given in Figure 2. The solution contains the detailed information regarding the operation mode (harvesting or idle), field and the time of operation of each crew at each event point. Note here that the moving mode ( ) is not included as a decision variable in GA, but simply treated as a dependent variable that is based on the decision made at each event point. For illustration purposes, each cell contains the information of the crew at that event point

Figure 2: GA solution representation of each chromosome

including the operation mode, the field number and the operational time. 4.1.2. Initial Population The initial population is comprised of a number of chromosomes and each chromosome contains the event information of each crew as demonstrated in Figure 2. For each event, the operation mode is sampled from two choices (harvesting or idling), the operation field is sampled from the given field choice and the operational time is sampled between 0 to a suitable number. All sampling procedures are uniform. Two modifications are introduced to the population generation procedure to produce a schedule that ensures all fields are harvested (as the 8-field case study shown in Figure 2). The first modification is to introduce ‘main’ and ‘side’ event points. The ‘main’ events are constrained so that they can only be harvesting events, the number of which is matched with the number of fields. Therefore, all fields are harvested once by the ‘main’ events. The ‘side’ events represent the operations with random operational mode, field, and operational duration. The second modification is to introduce a sampling constraint which ensures that the quantity of the sugarcane harvested by each crew from a given field is equal to the initial amount of cane ready for harvesting in that field. These modifications ensure that the generated solutions account for both randomness and feasibility. 4.1.3. Fitness The pseudocode given in Table 4 is used calculate the harvesting schedule for each crew ( ) at each event point ( ), if the task details such as the operational field ( ), operational mode (H/I) and operational duration ( ) are given.

13

Table 4: Pseudocode I for calculating the fitness of each solution

1 2 3 4 5

FOR each crew FOR each event IF is the starting event SET to warehouse SET to

6 7

SET

8 9 10

SET ELSE SET SET

11

SET

12

to to (

)

to to to (

13 14 15

SET END IF END FOR SET to

16

SET

)

to (

)

to

17 SET to 18 END FOR where denotes the field from which crew travelled at event point ; is the travel time of crew moving from previous field at event point ; is the travel time of crew returning from its final task to the mill warehouse; represents the time point, when crew finishes all its tasks; and defines the idle time crew spent at the warehouse. Additionally, to examine the time spent by crew at event that is outside the soft harvesting window ( ), the algorithm in Table 5 has been introduced. Table 5: Pseudocode II for calculating the fitness of each solution

1 2 3

FOR each crew For each event IF

4 5

SET ELIF

to

6

SET

to

7 8

AND

ELIF SET

AND to

9

ELIF

AND

10 11

SET ELIF

to

12 13 14

SET ELSE SET

to

15 16

END IF END FOR

AND

to

(

)

14

17 END FOR The above information is then used to construct the fitness function (Eq. (31)) ∑

∑

∑

∑

∑

(31)

where parameters are cost coefficients and and are penalty coefficients. In Eq. (31) the total travelling cost is determined by the travelling time of each crew at each event ∑ point ( ). represents the penalty costs induced by idling crews , which are determined by the idle time of each crew at each event point. The penalty term ( ∑ ) accounts for the idle time of each crew at the warehouse at the end of the ∑ harvesting operation whereas is calculated by the total harvesting time that sits ∑ outside the ideal harvesting window of the harvesting field. is introduced to penalise the harvesting time crews spend outside of the time horizon which ensures no such solution is passed to the next generation of population. The last term of Eq. (31) is used to force crews to attend each field once at most (which is generally the case in real world practise), where represents sum of the number of extra times each crew harvesting the same field more than once. The quantities and are positive and negative idle times defined in equations (32) and (33), respectively. (32) (33) 4.1.4. Crossover Crossover is an operator that combines two parents’ genetic information and passes it to new offspring. To preserve good characteristics, solutions with better fitness are more likely to undergo crossover recombination. This is implemented by ranking all solutions in the parental generation based on their fitness values, and the solutions with higher fitness will get higher probabilities of being chosen to breed. In this study, a tailor-made crossover operator is designed to preserve as much of the parents’ good characteristics as possible. The steps are detailed below. 1. As illustrated in Figure 3, two solutions are chosen as parent A and B, where the main and side events are highlighted by orange and white cells, respectively. Parent A Schedule

Parent B Schedule

Figure 3: example solution of parent A (left) and parent B (right)

15

2. The solution array of each parent is flattened to a one-dimensional array, and all the ‘Side’ events are removed from crossover. The positional index of each element of the flattened array in the original array is also stored (Figure 4) 3. The crossover operator used in this paper is the order 1 crossover, which is a fast operator as it requires virtually no overhead operations (Interactive, 2019). During the crossover, a random swath of consecutive alleles from parent A is first selected. Secondly, the swath is dropped to child 1 at the same position and the corresponding elements (based on the harvesting fields) in parent B are removed. Finally, starting from left, the empty cells in child A are filled with elements from parent B that are not removed. A detailed illustration is given in Figure 5. 4. Finally, the elements in the one-dimensional array of child A are used to replace the ‘main’ events in parent A one by one using their positional indices (Figure 6). Child B can be produced by flipping the starting position of parent A and B.

Parent A

Child A

Parent B

Child A

Parent A Schedule flattened

Figure 5: Order 1 crossover in detail

Parent B Schedule flattened 16

Figure 4: Flattened solutions of parent A and B

Child A

Figure 6: Child A from crossover

This particular type of crossover is especially applicable in this instance as it enables the new solution to inherit characteristics from both parents in a systematic way while ensuring all fields are still being harvested. However, in the child solutions, infeasible solutions may occur in the case when the maximum harvestable amount of a particular field is exceeded by the sum of the ‘main’ and ‘side’ events harvesting that field. To check and remove this infeasibility, a feasibility-check procedure is devised: the amount of sugarcane harvested from a given field is checked against the initial sugarcane inventory. If the total harvested amount exceeds the field’s initial inventory, the ‘main’ event harvesting that field is kept the same while all the ‘side’ events are resampled based on the procedure in 4.1.2 so that the solution becomes feasible. 4.1.5. Mutation After crossover, the new born child solutions are subject to a mutation operator, which aims to avoid premature convergence of the GA algorithm. While the crossover operator operates only on the ‘main’ events of parent solutions, the mutation operator is designed to alter only the ‘side’ events of the selected solutions. Each new solution has a pre-defined probability to mutate ( in this study). When mutation occurs, each of the ‘side’ events has a probability of 0.5 to alter their event type, operational field and event duration through resampling. In the case of infeasible solution, the feasibility-check procedure is implemented to repair the solution. 4.1.6. Next Generation The child solutions, mutated or not, become the initial population of the next generation. Secondly, the top solutions of the parental population are migrated to the next generation without any modification to preserve the convergence of the algorithm. Thirdly, to further introduce randomness and expand the search space, new solutions are generated by the operator discussed in section 4.1.2. 4.1.7. Solution Update and Stopping Criteria The fitness value of the best solution generated from current generation is compared with the best solution stored so far. If the new solution outperforms the stored solution, the solution will be updated. In this study, the stopping criteria for the GA is when algorithm reaches a predefined number of generations. The pseudocode for the GA is summarised in Supplementary Information (SI Table S2).

17

4.2.

Tabu Search

After the GA, the obtained solution is subject to tabu search to identify optimal harvesting sequences for each crew. Tabu search is a metaheuristic algorithm for solving optimisation problems, which uses a neighbourhood search procedure to iteratively move from one potential solution to a solution in the neighbourhood, until a stopping criterion satisfied. It explicitly makes use of memory structures to guide a hill-descending heuristic to continue exploration without being trapped by the absence of improvement movements. The capability of the tabu search to escape out of local optimums can lead to an efficient near-optimal solution. The following sections present the basic steps of tabu search applied to the harvesting scheduling problem in our study including diversification, which is used to search unexplored solution space respectively (Fred Glover & Laguna, 1997). 4.2.1. Initial Solution and Representation The solution from GA is basically a timetable represented as a matrix where each row represents the complete schedule for a given crew and each cell gives the detail operation of a given crew at a particular event. This solution serves as the initial solution to the tabu search. An example initial solution is shown in Figure 7, with each incident assigned a numbering tag (1-6).

Solution from GA

1

2

3

4

5

6

Figure 7: representation of the example initial solution

4.2.2. Neighbourhood Definition A neighbourhood to the given solution is defined as any other solutions that are obtained by a pair wise exchange of any two cells in the solution. Figure 8 exemplifies a neighbourhood solution obtained by swapping incidents 2 and 5 of the initial solution. Starting from an initial solution, the neighbourhood with the best fitness value is selected and a new current solution is

1

2

3

4

5

6

Figure 8: Neighbourhood solution obtained by swapping the order of incident 2 and 5 of the example solution

18

produced. 4.2.3. Tabu List One of the characteristics that allows tabu search to enhance its performance is that neighbourhood movements deteriorating the fitness function are also permitted if no improving move is available. Thus, to avoid cycling solutions, a mechanism called short-term memory is employed. These memory structures form the tabu list, a set of forbidden movements used to Recency: Short-term Memory

1

2

3

4

5

Recency: Short-term Memory

6

1 2

1

1

2

2

3

3

4

4

5

5

6

6 Frequency: Long-term Memory

3

4

5

6

n

1 Frequency: Long-term Memory

Figure 9: tabu structure of the example solution: initial tabu (left); after the 2-5 swap (right)

filter solutions that are recently visited from the search. In another word, any movements that have been made recently (less than iterations ago, also known as tabu tenure) will not be revisited. In this work, the tabu list is stored in the upper triangle of the tabu structure matrix given in Figure 9. When the swap is made between cell No.2 and No.5, the tabu tenure of that particular move is updated, changing to , which represents the number of iterations this move will be forbidden for. For each iteration that follows, the previous tabu tenure will drop by 1 ( ) and the tabu tenure of the swap of that iteration will be added. Once the tabu tenure of a movement drops to 0, the ban of this movement is lifted. 4.2.4. Aspiration criterion The short-term memory strategy can be restrictive, which may prohibit attractive moves, even when there is no danger of cycling, or lead to an overall stagnation of the searching process. Therefore, it is necessary to revoke tabu tenure of a particular movement if the new solution generated from that movement satisfies an aspiration criterion. The criterion adopted in this paper is that if the objective value of the new solution is better than that of the current bestknown solution. 4.2.5. Diversification The process may become trapped in a space of local optimum, i.e. no improvement in the objective value for a large number of iterations. Therefore, to expand the solution space to look for the global optimum, it is necessary to diversify the search process, driving it into new regions of solutions, which is referred as the diversification procedure. This is implemented in 19

the current problem using long-term ‘frequency based’ memory, which is stored in the lower triangle of the tabu structure given in Figure 9 (each time a move is made, 2-5 in this case, the frequency memory of that move is increased by 1). Once the diversification process is activated, the frequency information is used to penalise moves by assigning a larger penalty (frequency count adjusted by a suitable penalty factor d) to swaps with greater frequency counts (Eq. (34)). The diversification process is deactivated once a better solution is found than the current best solution, and the standard tabu search takes over. Diversification strategies are particularly helpful when better solutions can be reached only by crossing barriers or “humps” in the solution space topology. (

)

(34)

4.2.6. Termination Criteria Similar to the genetic algorithm, tabu search terminates if a pre-specified number of iterations is reached. The detailed pseudocode for the tabu search is given in the Supplementary Information Table S3.

4.3.

Hybrid Algorithm

The GA-tabu hybrid algorithm used in this paper is illustrated in Figure 10. In summary, the optimal solution of the GA is used as the starting point for the tabu search, which explores the solution search space in an organised manner.

20

Using the solution from GA as the initial solution, initialize the parameters of the tabu search

START

Set input parameters

Yes

If diversification?

Initial population generation

No

Guided search based on the shortterm memory and the aspiration criteria

Calculation of the fitness value of each individual

Diversification search based on both the short-term and long-term memory as well as the aspiration criteria

Selection of fine individuals for crossover Crossover operation

Update the best solution of tabu search

Mutation operation No

Feasibility check

No

Construct the population of the next generation/store the best solution

If GA termination condition met?

Yes

If tabu search termination condition met?

Output the solution END

Figure 10: the flow chart of the GA-tabu hybrid algorithm

21

5. Case Study In this section, medium-sized sugarcane farmland in KwaZulu Natal region in South Africa ( ) is studied to demonstrate the model’s functionality. KwaZulu Natal is a typical sugarcane farm region with an annual precipitation of higher than 1200 mm (historical statistics), which offers desirable climatic conditions for non-irrigated sugarcane growth. As given in Figure 11, a polygon map of the selected farmland is presented, which was generated using ArcGIS based on a satellite image. The red dot indicates the location of mill/warehouse, which is supplied by the sugarcane paddocks nearby (light blue fields). By calculating the harvestable area of each paddock using the image and removing the ones with harvestable area that is less than 5 ha, the resulting case contains 21 sugarcane paddocks and are labelled accordingly.

5.1.

Parameters Initialisation

In Figure 11, the numbering of each field is based on its area in ascending order and the harvestable area of each field is summarised in Table S4 in Supplementary Information. In this paper, the objective is to generate a harvesting plan and an operational schedule that can optimally allocate limited harvesting resources to fulfil the harvesting need of a given sugarcane region. In this particular case, the harvesting window of the entire farm is set to be a month (0 – 720 hrs), i.e. and for all fields. Additionally, the given fields are divided into four regions (R1 – R4) based on each field’s location. Within each region, given sugarcane varieties are planted in each field at similar time, resulting in the slight deviation in the ideal harvesting window of each field, which represents the time when the sugar content is at its peak. The general range of the harvesting time window of each region is summarised in Table S5 in Supplementary Information. The detailed starting time and finishing time of the ideal harvesting window of each field are given in Tables S6 and S7 in Supplementary

Warehouse

R 2

R 1

R 3

R 4

Figure 11: Sugarcane fields that are included in the optimisation model (numbered 1- 21, 0 is the mill)

22

Information.

The travelling distances between fields were approximated using the straight-line distances between each field’s centroid and are summarised in supplementary information (SI Table S1). All other model input parameters are summarised in Supplementary Information Table S8.

5.2.

The Exact Solution

The number of available crews is set to be 4 with the total number of event points being 14. The model is formulated in Python 3.7 using Pyomo 5.6.2 and is solved on a 3.20GHz 128GB RAM computer (16-Core Intel(R) Xeon(R) CPU E5-2667 v4). The case study is solved with two different optimality gaps i.e. 5% and 3.5%. The different resulting optimal schedules are given in Figure 12 and Figure 13 respectively with their solution statistics summarised in Table 6. As demonstrated in the figures, with an increment of just 1.5% in the optimality gap, the solution schedule improves significantly where more fields are harvested within their designated harvesting window. As the optimal solutions detailed in Supplementary Information Table S9 and Table S10, more efficient harvesting schedule (less travelling and idle time) is achieved with optimality gap of 3.5% than that of 5%. However, even with only 21 fields, 4 crews and 14 event points, the CPU time already reaches almost 3 hours at an optimality gap of 5%, and a vast leap to more than 4 days with a further 1.5% improvement in the objective function. The high CPU time could be attributable to: 1) the problem’s resemblance towards a generalised assignment problem, which is proven to be NPhard and difficult to solve; 2) a swap between the schedules of any crew would still be the optimal solution, leading to a degeneracy issue which can significantly slow down the search speed of branch and bound; 3) the complex formulation of the soft harvesting window in section 3.7 generates a large number of alternative solutions with similar objective values, which increases the solution time.

23

Computational time would not hinder the decision making process if the modelled system remains unchanged over operational horizon. However, the agriculture systems change continuously in response to environmental variables. Emerging digital technologies and data advances (e.g. smart machinery, remote sensing) provide tremendous opportunities - enabling the coordinated agriculture systems and bringing real-time data and information to precision decision. Such advances also demand responsive decision-making tools, which require a significantly reduced computational time with its trade-offs with solution optimality. This is of particular interest for industries and often real-time fast near-optimal solutions are more desirable for decision-support than slow-response optimal solutions. Thereby, in this study, we explored the solutions based on meta-heuristic approaches, which are discussed in the next sections.

Figure 12: The optimal harvesting scheduling obtained using CPLEX with and optimality gap of 5%

24

Figure 13: The optimal harvesting schedule obtained using CPLEX with and optimality gap of 3.5%

Table 6: Solution statistics of the full solution

Solution Statistics Objective Value Best Bound Gap CPU Time(s) Number of Constraints Number of Variables Number of Binary Variables

5.3.

5% Case 3.5% Case 357,900 375,518 4.9% 9278 106,943 35,285 31,860

363,006 375,518 3.4% 368,336 106,943 35,285 31,860

GA Solutions

The proposed GA given in section 4.1 was applied to the case study. However, the discrepancy between the objective function of the original MILP model and the fitness function of the GA algorithm does not support a direct performance comparison between the GA algorithms and the exact method. Therefore, a secondary fitness indicator named as the ‘Objective Value’ (OV) is introduced. The OV is designed to be equivalent to the objective function value in the exact optimisation to facilitate comparison of solutions using GA and the exact algorithms. The OV is defined in Eq. (35). : ∑

∑

∑

∑

∑

(35)

25

where is the harvesting time spent by each crew at each event point, is the profit ) in the objective function of the exact model. The performance coefficient equivalent to ( of the algorithm is tested using different sized population ranging from 25 to 400 with the number of generations varying from 50 to 200 (Table 7). Table 7: Case study results and model statistics for the proposed GA algorithm using different population size and number of generations

Avg. Worst Best Gap Avg. CPU Populat Termination Avg. OVb a c d fitness OV OV (%) time (s) ion size generation -61,954 352,314 977 350,822 354,689 -1.56 5.73 25 50 -61,592 352,431 984 347,831 354,071 -1.53 11.31 50 50 -60,207 353,478 819 352,228 355,540 -1.24 21.00 100 50 -1.01 -58,999 354,284 819 352,769 356,485 43.43 200 50 -58,637 354,525 735 352,707 355,782 -0.94 88.48 400 50 -59,068 354,238 985 352,760 357,195 -1.02 39.16 100 100 -57,755 355,124 761 353,709 356,867 -0.78 78.78 200 100 -56,996 355,620 799 354,379 357,526 -0.64 167.29 400 100 -57,252 355,449 964 353,738 358,288 -0.68 76.75 100 200 -56,291 356,089 1021 354,487 358,548 -0.51 166.06 200 200 -55,193 356,821 944 355,097 359,656 -0.30 352.58 400 200 a Average of the best fitness values from 50 independent GA runs b Average objective values from 50 independent GA runs (mean standard deviation) c Worst objective value obtained from 50 independent GA runs d Best objective value obtained from 50 independent GA runs

By increasing the population size from 25 to 400, the values of the average fitness and OV improved. This can be explained by the fact that the greater the number of solutions at each iteration, the better the chance that a good solution could enter the generation pool. However, the improvement hits a plateau after the population size reaches 200 in all cases with different generation numbers. On the other hand, the average OV does improve gradually with the

Figure 14: average OV with different population sizes and numbers of generations

26

increase of the number of generations per run as shown in Figure 14. We compared the solutions using exact methods with one of the best solutions derived from GA, which is generated using a population size of 400 with the generation number of 200 and has an objective value of 359,656. The comparison suggests that the objective of the best GA solution is already better than the solution obtained from the exact method, which is 5% from optimality. Nevertheless, the GA still underperforms when its solution quality is compared to that of the exact method. The sacrifice of optimality using GA does bring computational time advantages where the solution time was significantly reduced from 9278s to 352s when applying GA instead of the exact method. To further enhance the average performance of the GA model, the proposed tabu search algorithm is incorporated.

5.4.

Genetic Algorithm with tabu Search

Section 4.3 details the hybrid approach incorporating tabu search with GA and the model statistics are summarised in Table 8. A sample solution is given in Figure 16 and Supplementary Information Table S11. The performance of the hybrid algorithm is compared with the solutions using GA with the same population size and number of generations (Figure 15). The effect of increasing population size on the performance of the GA algorithm is more significant than its impact on GA-tabu search hybrid approach. The number of generations in GA also plays a marginal role in the hybrid algorithm, which is less significant in comparison with the case where only GA algorithm is used. By comparing the highlighted results in Table 8 with the best solution obtained using GA only, it is clear that tabu search, even when applied to intermediate solutions from GA still outperforms the best solution from GA. The improvement in terms of the average OV is significant at around 1% while the CPU time drops from 352s to 204s on average. The solution comparison shown in Figure 15 suggests that the hybrid algorithms deliver good performance which is comparable to the exact solutions obtained with the optimality gap of 5%. However, the result is still far from the global optimality as it falls short when compared to the exact solution with 3.5% optimality gap (Figure 15). Table 8: Case study results and model statistics for the proposed hybrid algorithm using different population size and number of generations

Avg. fitnessa

Avg.

-52,220 -52,088 -52,071 -51,948 -51,794 -51,668 -51,655 -51,113 -51,213

358,804 358,891 358,904 358,985 359088 359,172 359,185 359,542 359,475

OVb 1411 1178 1216 1294 1173 1100 1127 1346 1159

Worst

OVc

355,510 355,238 361,643 356,244 356,550 356,567 356,577 356,325 357,487

OVd

Gap (%)

361,747 361,164 355,980 362,479 361,274 361,089 361,680 362,790 362,391

0.25 0.28 0.28 0.30 0.33 0.36 0.36 0.46 0.44

Best

Avg. CPU time (s) 132.85 152.73 212.29 157.77 204.41 301.71 201.88 291.20 481.69

Populat ion size

Termination generation

100 200 400 100 200 400 100 200 400

50 50 50 100 100 100 200 200 200

Average of the best fitness values from 50 independent runs of the hybrid algorithm Average objective values from 50 independent runs (mean standard deviation) of the hybrid algorithm c Worst objective value obtained from 50 independent runs of the hybrid algorithm d Best objective value obtained from 50 independent runs of the hybrid algorithm a

b

27

GA

Hybrid

364000

Exact solution (3.5%)

362000

Avg.OV

360000 358000

Exact solution (5%)

356000 354000 352000 350000 348000 P= 100, P= 100, P= 100, P= 200, P= 200, P= 200, P= 400, P= 400, P= 400, G = 50 G = 100 G = 200 G = 50 G = 100 G = 200 G = 50 G = 100 G = 200

P= population, G = generation Figure 15: performance comparison of the exact, the GA and the hybrid approaches

Figure 16: a sample solution of the hybrid algorithm with P=100 and G=200

28

6. Conclusion This study addresses the modelling gaps and presents an event-based continuous time MILP model to optimise short-term sugarcane harvesting operations. This represents a step forward in comparison with the harvesting optimisation models with discrete-time formulations, which by definition are approximations of actual operational problems. Our proposed model not only captures the continuously changing agricultural system but also enables data-driven optimisation (i.e. data on crop development and harvesting periods). Our research also highlights precision decision-making, which brings real-time data and information into operational optimisation. Specifically, we proposed two meta-heuristic algorithms i.e. tailored genetic algorithm and a hybrid genetic-tabu search algorithm for solving the deterministic optimisation. We presented a hypothetical case study to demonstrate the model functionality and solutions using different search algorithms. In comparison with the optimal solutions derived from the exact algorithm, the proposed meta-heuristic methods lead to near optimal solutions (less than 5% from optimality) and significantly reduced computational time by over 95% (in comparison with the solution time of almost 2.5hrs using exact method). Our results demonstrate that the tailored genetic algorithm enables rapid solution searching but the solution quality on average cannot compete with the exact method; however, the hybrid genetic-tabu search algorithm delivers a good performance with solutions comparable to those from exact methods. Although the hybrid algorithm falls short when compared with the exact solution with 3.5% optimality gap it achieves a good trade-off between computational time and solution quality. Thus the proposed hybrid metaheuristic algorithm has the potential to enhance responsive decision making in precision farming and reduce the computational time of deterministic optimisation for sugarcane harvesting. Based on the research gaps identified in this study, several research challenges emerge representing frontier research directions worth exploring in future work – 







Development of a model for the deterministic optimisation of precision farming in particular the precision planning and operations in agro-ecosystems which are interconnected by smart machinery and sensors. Improvement of the model’s solvability using exact methods by means such as the incorporation of valid inequalities, utilisation of indicator constraints and reformulation to avoid symmetry. Further development in meta-heuristic algorithms notably the hybrid approach e.g. fine tuning in genetic-tabu search algorithms to enable responsive decision making alongside rapid development in digital technologies and data advances. Development of a generic tool for scheduling problems in precision agriculture which accounts for both economic feasibility and environmental sustainability. This would benefit wider societies particularly in response to the UN food sustainable development goals.

7. Acknowledgement Authors are grateful to the UK Engineering and Physical Sciences Research Council (EPSRC) for providing financial support for the research project ‘Towards Agriculture 4.0 – Computer-Aided Decision-Making Tool To Support Precision Farming’ through Knowledge Transfer Secondment funding scheme. M.G would like to acknowledges UK EPSRC for providing financial support for 29

her research through the EPSRC Fellowship project ‘Resilient and Sustainable Biorenewable Systems Engineering Model (ReSBio)’ (grant reference: EP/N034740/1).

8. References Al Chami, Z., Manier, H., Manier, M. A., & Chebib, E. (2018). An advanced GRASP-HGA combination to solve a multi-period Pickup and Delivery Problem. Expert Systems with Applications, 105, 262-272. Belotti, P., Kirches, C., Leyffer, S., Linderoth, J., Luedtke, J., & Mahajan, A. (2013). Mixed-integer nonlinear optimization. Acta Numerica, 22, 1-131. Blum, C., Puchinger, J., Raidl, G. R., & Roli, A. (2011). Hybrid metaheuristics in combinatorial optimization: A survey. Applied Soft Computing, 11, 4135-4151. da Silva, A. F., Marins, F. A. S., & Dias, E. X. (2015). Addressing uncertainty in sugarcane harvest planning through a revised multi-choice goal programming model. Applied Mathematical Modelling, 39, 5540-5558. Devika, K., Jafarian, A., & Nourbakhsh, V. (2014). Designing a sustainable closed-loop supply chain network based on triple bottom line approach: A comparison of metaheuristics hybridization techniques. European Journal of Operational Research, 235, 594-615. Duan, Q., Liao, T. W., & Yi, H. Z. (2013). A comparative study of different local search application strategies in hybrid metaheuristics. Applied Soft Computing, 13, 1464-1477. Ebadi, Y., & Navimipour, N. J. (2019). An energy-aware method for data replication in the cloud environments using a Tabu search and particle swarm optimization algorithm. Concurrency and Computation-Practice & Experience, 31. Farzad, S., Mandegari, M. A., Guo, M., Haigh, K. F., Shah, N., & Gorgens, J. F. (2017). Multi-product biorefineries from lignocelluloses: a pathway to revitalisation of the sugar industry? Biotechnology for Biofuels, 10. Florentino, Helenice d. O., Irawan, C., Aliano, A. F., Jones, D. F., Cantane, D. R., & Nervis, J. J. (2018). A multiple objective methodology for sugarcane harvest management with varying maturation periods. Annals of Operations Research, 267, 153-177. Floudas, C. A., & Lin, X. (2004). Continuous-time versus discrete-time approaches for scheduling of chemical processes: a review. Computers & Chemical Engineering, 28, 2109-2129. Gao, H., Kwong, S., Fan, B. J., & Wang, R. (2014). A Hybrid Particle-Swarm Tabu Search Algorithm for Solving Job Shop Scheduling Problems. Ieee Transactions on Industrial Informatics, 10, 2044-2054. Glover, F. (1990). TABU SEARCH - A TUTORIAL. Interfaces, 20, 74-94. Glover, F., & Laguna, M. (1997). Tabu Search Principles. In F. Glover & M. Laguna (Eds.), Tabu Search (pp. 125-151). Boston, MA: Springer US. Higgins, A. (2006). Scheduling of road vehicles in sugarcane transport: A case study at an Australian sugar mill. European Journal of Operational Research, 170, 987-1000. Higgins, A. J. (1999). Optimizing cane supply decisions within a sugar mill region. Journal of Scheduling, 2, 229-244. Interactive, R. (2019). Order 1Crossover. In (Vol. 2019). Jena, S. D., & Poggi, M. (2013). Harvest planning in the Brazilian sugar cane industry via mixed integer programming. European Journal of Operational Research, 230, 374-384. Junqueira, R. d. Á. R., & Morabito, R. (2017). Abordagens de otimização para a programação e sequenciamento das frentes de colheita de cana-de-açúcar. Gestão & Produção, 24, 407422. Karlsson, J., Rönnqvist, M., & Bergström, J. (2003). Short-term harvest planning including scheduling of harvest crews. International Transactions in Operational Research, 10, 413-431. Kondili, E., Pantelides, C. C., & Sargent, R. W. H. (1993). A general algorithm for short-term scheduling of batch operations—I. MILP formulation. Computers & Chemical Engineering, 17, 211-227. 30

Kusumastuti, R. D., Donk, D. P. v., & Teunter, R. (2016). Crop-related harvesting and processing planning: a review. International Journal of Production Economics, 174, 76-92. Lamsal, K., Jones, P. C., & Thomas, B. W. (2016a). Continuous time scheduling for sugarcane harvest logistics in Louisiana AU. International Journal of Production Research, 54, 616627. Lamsal, K., Jones, P. C., & Thomas, B. W. (2016b). Sugarcane Harvest Logistics in Brazil. Transportation Science, 51, 771-789. Liu, S. Q., & Kozan, E. (2016). Parallel-identical-machine job-shop scheduling with different stage-dependent buffering requirements. Computers & Operations Research, 74, 31-41. Lopes Silva, M. A., de Souza, S. R., Freitas Souza, M. J., & de França Filho, M. F. (2018). Hybrid metaheuristics and multi-agent systems for solving optimization problems: A review of frameworks and a comparative analysis. Applied Soft Computing, 71, 433-459. Ma, S., Karkee, M., Scharf, P., & zhang, q. (2014). Sugarcane harvester technology: A critical overview (Vol. 30). Marimuthu, P., Arumugam, R., & Ali, J. (2018). Hybrid Metaheuristic Algorithm for Real Time Task Assignment Problem in Heterogeneous Multiprocessors. International Arab Journal of Information Technology, 15, 445-453. Nations., F. a. A. O. o. t. U. (2018). OECD‑FAO Agricultural Outlook 2018‑2027. In. Nouri, H. E., Driss, O. B., & Ghedira, K. (2016). Hybrid metaheuristics for scheduling of machines and transport robots in job shop environment. Applied Intelligence, 45, 808-828. O Ting, T., Yang, X.-S., Cheng, S., & Huang, K. (2015). Hybrid Metaheuristic Algorithms: Past, Present, and Future. In (Vol. 585). Parejo, J. A., Ruiz-Cortés, A., Lozano, S., & Fernandez, P. (2012). Metaheuristic optimization frameworks: a survey and benchmarking. Soft Computing, 16, 527-561. Parvasi, S. P., Mahmoodjanloo, M., & Setak, M. (2017). A bi-level school bus routing problem with bus stops selection and possibility of demand outsourcing. Applied Soft Computing, 61, 222-238. Ramburan, S. (2011). Sugarcane cultivar x time of harvest interactions in South Africa. South African Journal of Plant and Soil, 28, 75-84. Ramos, R. I., Paulo. (2016). An optimization model for the combined planning and harvesting of sugarcane with maturity considerations. African Journal of Agricultural Research, 11, 3950-3958. . Recommendations for implementing the strategic initiative INDUSTRIE 4.0. In. (2013): National Academy of Science and Engineering. Sahin, M., & Kellegoz, T. (2019). A new mixed-integer linear programming formulation and particle swarm optimization based hybrid heuristic for the problem of resource investment and balancing of the assembly line with multi-manned workstations. Computers & Industrial Engineering, 133, 107-120. Şahin, M., & Kellegöz, T. (2019). A new mixed-integer linear programming formulation and particle swarm optimization based hybrid heuristic for the problem of resource investment and balancing of the assembly line with multi-manned workstations. Computers & Industrial Engineering, 133, 107-120. Santoro, E., Soler, E. M., & Cherri, A. C. (2017). Route optimization in mechanized sugarcane harvesting. Computers and Electronics in Agriculture, 141, 140-146. Schneider, M., Stenger, A., & Goeke, D. (2014). The Electric Vehicle-Routing Problem with Time Windows and Recharging Stations. Transportation Science, 48, 500-520. Sethanan, K., & Neungmatcha, W. (2016). Multi-objective particle swarm optimization for mechanical harvester route planning of sugarcane field operations. European Journal of Operational Research, 252, 969-984. Solomon, S. (2009). Post-harvest deterioration of sugarcane. Sugar Tech, 11, 109-123. Statistica. (2018). World sugar cane production from 1965 to 2017 (in million metric tons). In. Stray, B. J., van Vuuren, J. H., & Bezuidenhout, C. N. (2012). An optimisation-based seasonal sugarcane harvest scheduling decision support system for commercial growers in South Africa. Computers and Electronics in Agriculture, 83, 21-31. 31

Talbi, E. G. (2002). A Taxonomy of Hybrid Metaheuristics. Journal of Heuristics, 8, 541-564. Vahdani, B., Tavakkoli-Moghaddam, R., Zandieh, M., & Razmi, J. (2012). Vehicle routing scheduling using an enhanced hybrid optimization approach. Journal of Intelligent Manufacturing, 23, 759-774. Yanik, S., Bozkaya, B., & deKervenoael, R. (2014). A new VRPPD model and a hybrid heuristic solution approach for e-tailing. European Journal of Operational Research, 236, 879-890. Zhang, Q., Manier, H., & Manier, M. A. (2012). A genetic algorithm with tabu search procedure for flexible job shop scheduling with transportation constraints and bounded processing times. Computers & Operations Research, 39, 1713-1723. Zhang, Y. P., Yang, Y., & Yang, R. (2019). DISTRIBUTION PATH OPTIMIZATION METHOD OF GAS CYLINDER BASED ON GENETIC-TABU HYBRID ALGORITHM. International Journal of Innovative Computing Information and Control, 15, 773-782. Zhou, K., Leck Jensen, A., Sørensen, C. G., Busato, P., & Bothtis, D. D. (2014). Agricultural operations planning in fields with multiple obstacle areas. Computers and Electronics in Agriculture, 109, 12-22.

32