An Application of Optimal Control to Sugarcane Harvesting in Thailand

: The sugar industry is of great importance to the Thai economy. In general, the government sets sugarcane prices at the beginning of each harvesting season based on type (fresh or ﬁred), sweetness (sugar content) and gross weight. The main aim of the present research is to use optimal control to ﬁnd optimal sugarcane harvesting policies for fresh and ﬁred sugarcane for the four sugarcane producing regions of Thailand, namely North, Central, East and North-east, for harvesting seasons 2012/13, 2013/14, 2014/15, 2017/18 and 2018/19. The optimality problem is to determine the harvesting policy which gives maximum proﬁt to the farmers subject to constraints on the maximum amount that can be cut in each day, where a harvesting policy is deﬁned as the amount of each type of sugarcane harvested and delivered to the sugar factories during each day of a harvesting season. The results from the optimal control methods are also compared with results from three optimization methods, namely bi-objective, linear programming and quasi-Newton. The results suggest that discrete optimal control is the most effective of the ﬁve methods considered. The data used in this paper were obtained from the Ministry of Industry and the Ministry of Agriculture and Co-operatives of the Royal Thai government.


Introduction
The main aims of this paper are: (1) to use two optimal control methods to determine optimal harvesting policies for sugarcane in Thailand, where an optimal harvesting policy is defined as the amount of sugarcane harvested each day during a crop year that gives maximum profit to farmers; (2) To compare the efficiency and results of the two optimal control methods with three optimization methods. The two optimal control methods we consider are based on the discrete and continuous Pontryagin maximum principles (see, e.g., [1][2][3]). The three optimization methods we consider are the bi-objective ε-constraints method (see, e.g., [4]), a quasi-Newton optimization method (see, e.g., [5][6][7]) and linear programming (see, e.g., [8]). We use the five optimal control and optimization methods to determine optimal harvesting policies for the two types of sugarcane (fresh and fired) in the four regions (North, Central, East and North-east) of Thailand for crop years 2012/13, 2013/14, 2014/15, 2017/18 and 2018/19. We use data obtained from the Ministry of Industry and the Ministry of Agriculture and Co-operatives of the Royal Thai government [9][10][11][12][13]. The present paper is an extension of a previous paper [14] in which bi-objective and quasi-Newton methods were used to find optimal harvesting policies for the crop years 2012/13, 2013/14 and 2014/15. For these two optimization methods, it was necessary to divide a harvesting season into 15-day periods and to assume that there was zero growth of sugarcane during a season.
Although there have been few, if any, applications of optimal control to sugarcane harvesting, there have been applications of optimal control to similar types of problems in other areas which could be adapted to sugarcane farming [3]. For example, in 2011, Kiataramkul et al. [15] studied an optimal control problem for optimal nutritional intake for fetal growth in sheep. In 2014, Dine, Lenhart and Behncke [16] investigated discretetime optimal harvesting of fish populations with age structure with the two objectives of maximizing the profit of fishing by finding the optimal harvesting strategy for each age class and of finding the optimal net size. In 2015, Kiataramkul and Matkhao [17] studied the optimal control problem of food intake of swine during the post weaning period with the objective of minimizing the amount of food fed to the swine and achieving an optimal weight at time of sale. In 2016, Puengpo et al. [18] applied continuous and discrete optimal control models to find optimal methods for sheep and swine feeding problems.
In contrast to optimal control, there have been a number of applications of optimization methods to sugar production. Some examples of applications of optimization methods to sugarcane harvesting and sugar production include the following. In 2012, Gomes [19] studied a bi-objective mathematical model for choosing sugarcane varieties in Brazil which have a harvest biomass residual that could be used in electricity generation. The bi-objective optimization model aimed to minimize the cost of harvesting the residual biomass and to maximize the revenue from sale of the generated electricity. In 2016, Sungnul et al. [20] studied a bi-objective optimization model to find an optimal time of harvesting for sugarcane growers in the North-east region of Thailand. The aim of this work was to help farmers to find the optimal harvesting time in order to maximize the revenue and to minimize the cost. Sungnul et al. used the ε-constraints method [4] to find the optimal cutting pattern by using the revenue as the objective function and the costs as constraints. In 2017, Sungnul et al. [21] extended the work in [20] to find the optimal harvesting times for all of the four regions of Thailand. Quasi-Newton optimization methods are well-known methods of optimization that have been used for many years to find optimal solutions for problems in many areas of science, finance and industry (see, e.g., [5][6][7]). In 2019, Pornprakun et al. [14] used the bi-objective and quasi-Newton optimization methods (see, e.g., [5][6][7]) to determine optimal harvesting policies for sugarcane in Thailand in order to maximize the total revenue and minimize harvesting cost for crop years 2012/13, 2013/14 and 2014/15.

Mathematical Models
We used the following mathematical model as the basic model for maximizing profit from sugarcane harvesting.
where J[u] is the profit functional to be maximized by selecting the control variable u(t), which is the rate (tonnes/day) of cutting sugarcane at time t. Also, P(t, x(t), u(t)) is the profit for sugarcane cut at time t (baht/day), P(t, x(t), u(t)) is the profit for sugarcane cut at time t (baht/day), x(t) is the total amount of sugarcane (tonnes) on the farms in a region at time t, u max is the maximum rate of cutting sugarcane (tonnes/day), t 0 is the initial time at start of harvesting season (day), t f is the final time at end of harvesting season (day), a is the total amount of sugarcane on farms at the initial time t 0 (tonnes), is a logistic growth function for the rate of increase in weight of sugarcane (tonnes per day) on the farms in a region at time t, where r is a specific growth rate (1/day), and K is a constant which represents the carrying capacity of the farms in the region (tonnes).
We also considered the following discrete version of (1).
where x k = x(t k ) and u k = u(t k ) and t k , k = 0, 1, 2, . . . , N − 1 are N periods of a cutting season and t N is the end of the cutting season. For the optimal control methods, we used 1-day periods, for the bi-objective and quasi-Newton methods we found it necessary to assume 15-day periods, and for linear programming we considered both 1-day and 15-day periods.

Optimization and Optimal Control Methods
In this section, we give details of the three optimization and two optimal control methods that we used to obtain optimal policies for sugarcane harvesting in Thailand.

Bi-Objective Optimization
For the bi-objective optimization problem, we separated the profit P(t k , x k , u k ) per 15-day period into a revenue term R(t k , x k , u k ) and a cost term C(t k , x k , u k ) and assumed that r = 0, that is, no growth during the harvesting season. We then used the ε-constraints method [4,14] to solve the bi-objective optimization problem where the ε r are a set of values of costs between a minimum and a maximum value of the cost for a feasible cutting pattern.

Quasi-Newton Optimization
As in the previous paper [14], we used the constrained optimization function fmincon in Matlab with the "active-set" algorithm based on quasi-Newton method [5][6][7] to find the optimal cutting patterns to maximize the profit in (2) for 15-day periods and we assumed that there was no growth during the harvesting season.

Linear Programming
We used the Matlab program linprog with the "dual-simplex" algorithm based on linear programming [8] to solve the optimization problem in (2) for both 15-day periods and daily periods. We again assumed that there was no growth during the harvesting season because growth would make the optimization problem a nonlinear problem.
The results of three optimization methods for the years 2012/13, 2013/14 and 2014/15 have been published in [14] and the results for 2017/18, 2018/19 are included in this paper.

Pontryagin Maximum Principle for Continuous Optimal Control
For a given type, region and harvesting season, the continuous optimal control problem is of the form given in (1). We used the continuous Pontryagin maximum principle to solve this continuous optimal control problem (see, e.g., [1][2][3]). It should be noted that the optimal control problem for sugar harvesting has state variable constraints on the amount that can be cut each day and on the total amount that is available for cutting.
The first step is to define a Hamiltonian where − u(t) and the λ(t) are called costate variables.
Then, the equations for the state variables x(t) and costate variables λ(t) are given by Since the boundary conditions on the state variables in (1) are given values at fixed initial and final times, the appropriate boundary conditions on the state and costate equations are Then, the Pontryagin maximum principle states that the optimal control u * (t) is obtained by finding the maximum of H(t, x * (t), u(t), λ * (t)) with respect to u(t), where the x * (t), λ * (t) are the solutions of the state and costate equations for u * (t).
In general, the maximum of the Hamiltonian can occur at the minimum or maximum values of the cutting constraints or at an internal point. For the problem in (1) the state equation is a linear function of the control u. In this case, there is no internal optimal value of u and the control is called a "bang-bang" control and the optimal values u * (t) are either the minimum u(t) = 0 or maximum u(t) = u max . However, as noted above, the sugar harvesting problem has the extra complication that the constraints on the control u * (t) also depend on the values of the state variables x * (t) because of the conditions x * (t) ≥ 0, x * (t f ) = 0. We have included these conditions as explained in step 4 of the following algorithm.

Continuous optimal control algorithm
The algorithm that we used for the solution of the continuous optimal control problem is as follows.

2.
Solve the state and costate equations for this pattern to obtain values of x(t k ) and λ(t k ). For reasons of numerical stability, the state equations must be integrated forwards in time and the costate equations must be integrated backwards in time.

4.
In our examples, we have found that all partial derivatives are positive. Therefore, if there were no constraints on the state variables, the optimal cutting pattern would be to cut the maximum amount available for cutting each day. In order to satisfy the state variable constraints x(t) ≥ 0 and x(t f ) = 0, we chose the optimal u(t k ) values as follows.
(a) Sort the derivatives ∂H(t k , x(t k ), u(t k ), λ(t k )) ∂u(t k ) in decreasing order.
(b) Find the number of days M required to cut all available cane at the maximum rate u max and select the times t k corresponding to the top M values of the derivatives. (c) If Mu max = a, set the new cutting pattern with u(t k ) = u max at these M values of t k and u(t k ) = 0 at all other times. However, if Mu max > a, set u(t k ) = u max for the top M − 1 derivative values and Note: A slight modification is required if growth can occur (r > 0).

5.
Then, stop if the new cutting pattern is within a selected tolerance of the old cutting pattern or return to step 2 and repeat if the new pattern is not within the selected tolerance of the old.

Pontryagin Maximum Principle for Discrete Optimal Control
The discrete version of the continuous optimal control problem in (1) is as follows.
where x k = x(t k ) and u k = u(t k ).
As for the continuous case, the first step is to define a Hamiltonian where f (k, x k , u k ) = x k + rx k 1 − x k K − u k . Then, the equations for the state variables x k and costate variables λ k are given by Since the boundary conditions on the state variables in (8) are given values at fixed initial and final times, the appropriate boundary conditions on the state and costate equations are As for the continuous case, the optimal cutting pattern will be "bang-bang". The algorithm for finding the optimal cutting pattern for the continuous case can also be used for the discrete case.

The Commercial Cane Sugar System (CCS) in Thailand
As details of the CCS system have been given in [14,22], we will only give a brief review here. Thai farmers harvest two types of sugarcane, namely fresh and fired, where fired sugarcane is burnt to remove leaves before it is cut so that it can be cut manually by workers, whereas fresh sugarcane is usually cut by machines which can remove the leaves as the cane is cut. In the CCS system, the Royal Thai government sets the price of fresh and fired sugarcane for each of the four sugarcane producing regions in Thailand, namely North, Central, East and North-east. The price of the sugarcane is based on the two main factors of weight and sweetness, where sweetness or CCS is the percentage of sugar in the sugarcane. All data in this paper have been obtained from the Office of the Cane and Sugar Board (OCSB), Ministry of Industry and the Ministry of Agriculture and Co-operatives of the Royal Thai government [9][10][11][12][13]. The OCSB reports data for 15-day periods of a harvesting season which typically starts around 15th November and ends around 15th June of the following year.
In the following, we will use the notation, i = A for fresh and i = B for fired sugarcane. For the four sugarcane-growing regions of Thailand, we will use the notation j = 1 for North, j = 2 for Central, j = 3 for East, and j = 4 for North-east.

Price of Sugarcane
As stated above, the price of sugarcane is set on the basis of type, weight, and sweetness (CCS) [14,20].

1.
Price based on weight: The basic price per weight of sugarcane (baht/tonne) is fixed by the Royal Thai government for each crop year. This basic price is the same for all regions. However, the actual price paid to farmers for fired sugarcane is 30 baht/tonne less than the basic price. Then, at the end of harvesting for the year, factories in each region will share the total amount of money deducted from fired sugarcane sales in that region to farmers who sold fresh sugarcane at a rate not exceeding 70 baht/tonne of fresh sugarcane delivered. We will use the notation P w (i, j) to denote actual price based on weight (baht/tonne) for type i in region j. The actual prices per tonne based on weight for sugarcane for a given crop year are: where P w is the basic price, a(A, j) is the total amount of fresh sugarcane (tonnes) harvested in region j and a(B, j) is the total amount of fired sugarcane (tonnes) harvested in region j for a given year.

2.
Price based on sweetness (CCS): Each year, the Royal Thai government sets a basic price per tonne for sugarcane with 10 CCS, where CCS is the percentage by weight of sugar in sugarcane. This price is the same for fresh and fired sugarcane and for all regions. The actual price per tonne received by farmers is then adjusted if the CCS is different from 10. For sugarcane harvested in period k in region j the actual price per tonne for a given year is where P c is the basic price per tonne for sugarcane with 10 CCS set by the government for the year, and y(k, j) = CCS(k, j) − 10, where CCS(k, j) is the average CCS from sugarcane harvested in period k in region j and the factor 0.06 is the rate of change of price per 1 CCS from the base level of 10.
Therefore, the total prices to farmers (baht/tonne) for a given crop year are: The basic prices per weight P w and for sweetness P c set by the Royal Thai government and the total amounts delivered to the factories are shown in Table 1. It can be seen from this data that the amount of fresh sugarcane is appreciably less than the amount of fired in all regions in all crop years and that the amount of sugarcane was increasing from 2012/13 to 2017/18 and then decreased in 2018/19. The CCS values are shown in Tables 2 and 3 [14]. It can be seen that in 2017/18 the CCS value in all regions reached a maximum in April and then decreased slightly, whereas in 2018/19 the CCS value increased in all regions during the harvesting season. The total prices for the fresh and fired sugarcane are shown in Tables 4 and 5 for 15-day periods for all regions for crop years 2017/18 and 2018/19. It can be seen that the prices in all regions increase rapidly at the beginning of a crop year and then slowly at the end due to the changes in CCS values. It can also be seen that the fresh prices are appreciably higher than the fired prices due to the price adjustments discussed in Section 4.1.

Costs of Production
The costs of production (baht/tonne) can be separated into cutting costs, transport costs and maintenance costs. In this paper, we assume that the cutting and transport costs are fixed costs that depend only on the amount cut in a given period, whereas the maintenance costs (baht/tonne) are variable costs that depend on the amount of uncut sugarcane remaining on the farm after cutting in previous periods.
We assume that the total cost of sugarcane production (baht) of type i in period k in region j in a given year is where u k (i, j) is the weight of sugarcane (tonnes) cut in period k, x k (i, j) is the weight of uncut sugarcane remaining on the farm in period k, C h (i, j) is harvesting cost of sugarcane (baht/tonne), C t (i, j) is transport cost for delivering sugarcane to the mills (baht/tonne) and C m (i, j) is maintenance cost (baht/tonne) for sugarcane remaining on the farms. The harvesting, transport and maintenance costs obtained from the OCSB are shown in Table 6 for each region for each crop year. It can be seen that the harvesting costs are lower in all regions for the years 2017/18 and 2018/19 for both fresh and fired cane, whereas there is no clear pattern for the transport and maintenance costs.

Total Revenue and Total Profit
For each type i and region j, the total revenue for a given year is and the total profit is

Cutting Constraints
We assume that there are constraints u max (i, j) on the amount of sugarcane of type i that can be cut each day in a given region j. These constraints could be due, for example, to the number of machines available for cutting fresh sugarcane or the number of workers avaiable for cutting fired sugarcane. In addition, there can be constraints due to the number of trucks available for transporting the cut cane to the mill and the cutting capacities of the mills. We also add the constraints that 0 ≤ u k (i, j) ≤ x k (i, j) and that all sugarcane must be cut, that is, x N (i, j) = 0 at the end of the harvesting season.

Profits
Examples of the optimal profits computed from the three optimization and two optimal control methods are shown in Table 7 for fresh sugarcane for the four regions of Thailand for the crop year 2017/18. It can be seen that all optimization and optimal control methods give the same results. A comparison of the theoretical results with the actual profit show that the best agreement is for mcf = 1, that is, full maintenance cost with an upper limit of approximately half of the total amount available cut each 15-day period.

Optimal Daily Cutting Patterns
In this section, we show the optimal daily cutting patterns that we obtained from discrete optimal control. We found that continuous optimal control gave the same optimal cutting patterns. Also, for the special case of zero growth, we found that linear programming gave the same optimal daily cutting patterns. For comparison, optimal cutting patterns computed from the bi-objective and quasi-Newton methods have been published in [14] for 15-day cutting periods for the crop years 2012/13, 2013/14 and 2014/15.
The results in Figures 1 and 2 show the optimal cutting patterns for the crop year 2017/18 if prices, costs, and profits and an upper bound on the amount cut per day are the main factors involved in determining the optimal profits of the farmers. We have computed optimal daily cutting patterns for the crop years 2012/13, 2013/14, 2014/15, 2018/19 and obtained similar results. For each type, region and year, we examined the effect of changing the values of the upper bound on the maximum cutting per day u max (tonnes per day) and the effect of reducing the maintenance costs by a factor mcf (0 ≤ mcf ≤ 1) of the actual maintenance cost. It can be seen that the optimal cutting patterns are sensitive to changes in the maintenance cost as they range from cutting all sugarcane as early as possible for mcf = 1 to cutting all sugarcane as late as possible for mcf = 0. Some examples of actual cutting patterns obtained from OCSB data [12,13] are shown in Figures 3 and 4 for comparison with the theoretical optimal cutting patterns. It can be seen that the actual cutting patterns for both fresh and fired sugarcane in 2017/18 are similar to the cutting patterns for mcf = 0.01 in Figures 1 and 2 with a maximum cutting of between 0.15-0.2 total available sugarcane in each region. As far as the authors are aware, the actual cutting patterns are due to constraints on cutting on the farms and to the capacities of the sugar mills.

Linear Progamming Dual Variables and Hamiltonian Derivatives
In the optimal control algorithms in Section 3, we used the maximum values of the Hamiltonian derivatives ∂H(t k , x k , u k , λ k ) ∂u k to find the optimal cutting patterns. It is also well known in linear programming that nonzero values of dual variables λ correspond to active constraints and that the values of dual variables correspond to the extra profit (baht) if a cutting constraint can be increased by one unit (tonne). It is therefore expected that maximum Hamiltonian derivative values should correspond to maximum λ values. As an example, Figure 5 shows a plot of the Hamiltonian derivative values for fresh sugarcane for the North-east region in 2017/18 and Figure 6 shows a similar plot for the dual variables from linear programming.

Computation Times
As we have seen, the three optimization methods and the two optimal control methods all give similar optimal cutting patterns and profits. However, there is a big difference in computation times between the different methods. As shown in Table 8, the linear programming, discrete and continuous optimal control methods all give very fast computation times in the order discrete optimal control, linear programming, continuous optimal control. However, the optimal control methods have the advantage that the objective functions and constraints can be nonlinear. As shown in Table 9, the linear programming again gives a fast computation time whereas the bi-objective and quasi-Newton methods are much slower. However, both the bi-objective and quasi-Newton methods again have the advantage that the objective functions and constaints can be nonlinear.

Discussion and Conclusions
We have compared five different optimization and optimal control methods for optimization of sugarcane harvesting in the four sugarcane growing areas of Thailand, namely, North, Central, East and North-east regions for fresh and fired sugarcane for the crop years 2012/13, 2013/14, 2014/15, 2017/18 and 2018/19. To build the mathematical models, we have used price and cost data from the Office of the Cane and Sugar Board [9][10][11][12][13].
For the bi-objective and quasi-Newton methods, we have assumed that there is no growth in sugarcane during the crop year and have divided the crop year into 15-day periods. For the discrete and continuous optimal control methods, we have allowed for the possibility of growth during a crop year through a logistic function and have divided the crop year into 1-day periods. We have also used linear programming to find optimal cutting patterns for both 1-day and 15-day periods for a no-growth model.
As shown in Table 7, the profits obtained from the five different optimization methods have been the same to at least 4-digit accuracy for both fresh and fired for all regions for the crop year 2017/18. We have obtained similar results for both types of sugarcane for all four regions of Thailand for the crop years 2012/13, 2013/14, 2014/15, 2017/18 and 2018/19. It can also be seen from the table that the actual profits are within the range of computed profits for mcf = 1 (full maintenance cost). It is also clear from the results that changes in the maintenance cost are one of the main factors affecting the optimal cutting patterns, with high maintenance costs suggesting cutting as early as possible and low maintenance costs suggesting cutting later.
Although the bi-objective, quasi-Newton, linear programming and optimal control methods give the same optimal cutting patterns and profits for the case of zero growth, there is a big difference in computation times as shown in Tables 8 and 9 with the linear programming and optimal control methods being orders of magnitude faster than the biobjective and quasi-Newton methods. Further, for the bi-objective, quasi-Newton methods, we have considered 15-day cutting periods for a total of approximately 12 decision variables, whereas for the linear programming and optimal control methods we have considered daily cutting periods for a total of approximately 180 decision variables. However, linear programming can only be used if the growth is zero and the profit function is a linear function of the cutting patterns and the constraints. Finally, we have found that the programming for the discrete optimal control is much simpler than for the continuous optimal control and that the computation times are also shorter.
There are a number of important factors that have not been considered in this paper. For example, there are likely to be changes in the harvesting and transport costs per tonne as a function of amount cut per day and also changes in the costs of harvesting and transport during the cutting season. Including these effects would make the profit a nonlinear function of control and change the solution from the bang-bang control observed in the models in this paper. Unfortunately, we have not been able to obtain the cost data that would be required to include these effects.