CAN LINEAR PROGRAMMING ASSIST METAHEURISTICS IN FOREST PRODUCTION PLANNING PROBLEM?

The planning of forest production requires the adoption of mathematical models to optimize the utilization of available resources. Hence, studies involving the improvement of decision-making processes must be performed. Herein, we evaluate an alternative method for improving the performance of metaheuristics when they are applied for identifying solutions to problems in forest production planning. The inclusion of a solution obtained by rounding the optimal solution of linear programming to a relaxed problem is investigated. Such a solution is included in the initial population of the clonal selection algorithm, genetic algorithm, simulated annealing, and variable neighborhood search metaheuristics when it is used to generate harvest and planting plans in an area measuring 4,210 ha comprising 120 management units with ages varying between 1 and 6 years. The same algorithms are executed without including the solutions mentioned in the initial population. Results show that the performance of the clonal selection algorithm, genetic algorithm, and variable neighborhood search algorithms improved significantly. Positive effects on the performance of the simulated annealing metaheuristic are not indicated. Hence, it is concluded that rounding off the solution to a relaxed problem is a good alternative for generating an initial solution for metaheuristics.


INTRODUCTION
Timber production planning instigates research pertaining to mathematical modeling. Since the early studies by Curtis (1962), Loucks (1964), Kidd et al. (1966), Nautiyal and Pearse (1967), and Ware and Clutter (1971), mathematical planning models have been proposed and applied to identify planting, harvesting, and forestry activities that must be applied to each forest compartment.
The models proposed by Johnson and Scheurman (1977), i.e., Model Types I and II, serve as basic models for forest planning and references for several relevant studies , Silva et al. (2003), and ). Areas of the management unit subject to a certain silvicultural prescription can be fractionated in these models. In several cases, the best management alternative must be defined for each stand without permission to fractionate its area, which results in binary responses (harvesting or not harvesting a stand in a specified year). This requires the use of whole mixed programming (Silva et al., 2003) or metaheuristics , Matos et al. (2019), and Ferreira et al. (2018)).
Previously, the abovementioned problem was solved by rounding the mathematical solution. However, this alternative is not applicable when the management unit is a stand (Silva et al., 2003). Integer linear FLORESTA,Curitiba,PR,v. 51,n. 3,jul programming and its variations (primarily mixed and binary whole) can be employed depending on the complexity of the modeled situation. However, it can result in a combinatorial nondeterministic polynomial (NP) problem, with no guarantee of an algorithm that returns an answer to the problem in a timely manner.
Researchers of artificial intelligence have proposed several algorithms for solving combinatorial problems. Examples include heuristics and metaheuristics (Jin et al. (2016), Dong et al. (2016), Caro et al. (2003), Meignan et al. (2012), Rodrigues et al. (2004a), Rodrigues et al. (2004b), andAraújo Júnior et al. (2017)).An initial solution must be determined to ensure the efficiency of metaheuristics (Gaspar-Cunha et al., 2016). A good algorithm allows an algorithm to converge quickly to regions of interest in the solution space. Several authors indicated that analysts must consider all aspects that render the generated solutions viable when modeling a complex problem.
A method for generating good initial solutions to a problem with one or more binary or integer variables is to apply linear programming to a relaxed problem (Akbulut et al., 2017). Subsequently, deterministic algorithms (i.e., simplex) can be used to solve the problem as they are not restricted by the fact that the variables are integers.
Herein, we present an evaluation complementary to that developed by Akbulut et al. (2017). The objective of this study was to test the methodology proposed by previous authors with some modifications to solve a forest production planning problem. The test was performed by employing four metaheuristics, i.e., clonal selection algorithm (CSA), genetic algorithm (GA), simulated annealing (SA), and variable neighborhood search (VNS), and then their results were compared.

Case Study
In the present study, data from a eucalyptus forest measuring 4,210 ha were used. The data contained 120 management units (stands) with ages ranging from 1 to 6 years. The distributions of area per age were as follows: 339, 768, 1031, 601, 958, and 513 ha for 1, 2, 3, 4, 5 and 6 years of age, respectively. The estimates of wood production in each stand and age were obtained using the equation Vi = 6,09 -117.55*Ii-1*Si-1, where Vi is the production of stand i (m3 ha-1) at age Ii (years), and Si is the site index of stand i (m). The site index values ranged from 22 to 31 m for the index age of 6 years.
The possibility of harvesting wood at the ages of 5, 6, and 7 years was considered, with immediate planting after cutting was performed and a planning horizon of 16 years. This resulted in 81 management prescriptions for each stand and 9.720 decision variables. The mathematical formulation was devised based on Model I to maintain the physical identity of the management units along the planning horizon (Johnson and Scheurman, 1977). The objective was to maximize the net present value while considering the minimum and maximum annual demands of 140,000 and 160,000 m³, respectively.
The discount rate used was 8% per year, the sale price of the wood was R$ 80.00/m³, and the cost of harvesting was R$ 30.00/m³. The forestry costs varied based on the stand forest age, and they were obtained from Araújo Júnior et al. (2017), i.e., R$ 4,059.05, R$ 1,627.81, R$ 757.95, and R$ 88.12 per hectare in the first, second, third, and fourth years of plantation development, respectively.

Linear Programming and Integer Linear Programming
The mathematical LP and ILP models were based on the study by . The only difference between the models is the absence of restrictions 2 and 5 for the LP.
where VPLG is the overall net present value for the entire forest, in reais; Cij is the net present value (NPV) for the i-field when the prescription j is marked, in reais; Xij is the decision variable of the model, represented by the FLORESTA,Curitiba,PR,v. 51,n. 3,jul proportion of the area of the stand i that will receive the alternative management j; M is the total number of stands; N is the total number of management alternatives; Vij(k) is the total volume of wood (m3) of stand i under prescription j in period k of the planning horizon; Dmínk and Dmaxk are the minimum and maximum wood demand values in the k period of the planning horizon, respectively. Constraint (2) ensures that the entire area of each management unit receives a unique prescription. Restrictions (3) and (4) limit the annual harvest volume between the minimum and maximum demand values. Finally, constraint (5), considered only for the ILP model, ensures that the decision variable is of the binary type.
Additionally, the metaheuristics accounted for the restrictions imposed on the LP and ILP models. In these cases, the objective function was penalized because of the non-observation of any of the constraints. Hence, the fitness function increased by R$ 100.00 for each cubic meter, which caused excess or inadequate wood in a specified period of the planning horizon. This method is the same as that adopted by .
The decision variable values obtained using the simplex algorithm were rounded to the nearest integer (0 or 1). The new solution created, herein referred to as rounded linear programming (RLP), was included as an initial solution for the processing of metaheuristics. Furthermore, the processing was performed without including RLPs to evaluate whether the results obtained were different. The "branch and bound" algorithm was used to obtain the solution of the ILP model.

Clonal Selection Algorithm
The CSA is based on the natural process of clonal selection of antibodies, in which only cells that recognize a specified antigen proliferate (Castro and von Zuben (2002); Wang et al. (2016)). During optimization, only the best solutions of the current generation are considered in the search process evolution to obtain better values. Araújo Júnior et al. (2017) demonstrated its use in a forest planning problem. Other relevant studies include those of Boussaid et al. (2013) and Qiu and Lau (2014).
The values considered for hypermutation, cloning, selection, and substitution rates were 0.20, 0.50, 0.80, and 0.20, respectively. The initial population comprised 80 individuals, and processing was interrupted after 100 generations. These values were defined as described by Araújo et al. (2017).

Genetic Algorithm
Metaheuristic GA is considered efficient for the exploration of various parts of a region comprising viable solutions, where it evolves gradually toward the best solution. It is inspired by the process of natural selection described by Charles Darwin and the study of genetics initiated by Gregor Mendel. The general concept is that the best individuals, i.e., the most adapted ones, tend to survive and propagate their genes through crosses. Eventually, mutation occurs during the formation of a new individual.
In the current study, the GA configuration used accounted for an initial population of 50 individuals and the selection of parents via the tournament method at a rate of 0.80; for the crossover with one breakpoint, a geneto-gene mutation rate of 0.01 with elitism was applied, and the final solution was obtained after 100 generations. These values were obtained from the study of Matos et al. (2019).

Simulated Annealing (SA)
Metaheuristic SA is analogous to physical annealing. It involves initially blowing a metal or glass at high temperatures and then slowly cooling the substance until it reaches a low-energy state and exhibits the desired physical properties. The energy levels of the atoms fluctuate with a tendency to decrease during the process. An increase or decrease in energy level occurs based on a probability function.
The current energy level corresponds to the objective function value for a viable solution in an optimization problem. The objective of achieving a stable state in the substance using minimal energy is to obtain a viable solution for the problem using the best possible value of the objective function.
For the SA algorithm, we considered an initial temperature of 108, a final temperature of 10-5, an acceptance probability of lower solutions of 0.05, and 80 neighbors. More details regarding the SA are available in the study by Ferreira et al. (2018).

Variable Neighborhood Search (VNS)
The VNS algorithm is based on a local search with different neighborhood structures (Doerner et al., 2007). These structures define a set of modifications that can be applied to generate new solutions (Meignan et al., 2012).
The proposed version of the VNS considers four structures to modify the neighborhood of the current solution. This implies that the prescriptions of a certain percentage of stands are randomly modified (Araújo et al., 2018) for each neighborhood structure, i.e., structures 1, 2, 3, and 4 allow changes in 1%, 2%, 3%, and 4% of stands, respectively. In this study, the structures were considered in an orderly manner during each iteration. The best solution obtained after the evaluation of the fourth structure was maintained. The algorithm evaluated all individuals in a neighborhood of size equivalent to 100 individuals in each local search procedure. The evaluation was terminated after 100 generations.

Evaluation
The best solution for each repetition was analyzed per metaheuristic. The minimum, maximum, and average values, as well as the standard deviation of the fitness function were obtained. The efficiency (Ef) of each algorithm was calculated based on the equation by  as follows: where fM is the value of the best solution obtained by metaheuristics (R$), and f0 is the value of the best solution generated by the branch-and-bound algorithm (R$). Boxplot graphs were constructed to compare the results of each metaheuristic with and without RLP. The Kruskal-Wallis (K-W) test was performed to test the hypothesis of equality between processing. The evolution of the maximum, average, and minimum solutions of each algorithm in each generation was evaluated.

Processing
LP and ILP models were solved using the CPLEX software (version 12.7.1; IBM Corporation, 2017) using the simplex and branch-and-bound algorithms. For the metaheuristics, the Metaheuristics for Forest Planning software, developed at the Laboratory of Operational Research and Forest Modeling of the Federal University of Minas Gerais, was used. All algorithms were executed on a 64-bit Windows 10 operating system computer with a 2.0 GHz Intel Core i7 processor and 8 GB of RAM.

RESULTS
The solution obtained using the LP model presented a global NPV of R$ 32,191,790. However, this solution did not contain any stands with a single management prescription. Therefore, the rounding of the solution caused did not satisfy the wood demand in five periods of the planning horizon, which totaled 15,608 m³ of wood. The net present value considering the penalty due to the non-attendance of the volumetric restriction was R$ 24,312,036.
The results obtained with the introduction of the optimal solution of LP with subsequent rounding were superior to those obtained without such consideration (Table 1), except when the SA algorithm was used. The solutions were similar to those of SA.
In terms of the maximum fitness in each algorithm, the greatest difference in percentage between using and not using RLP was presented by metaheuristic GA (4.60%). For the mean value of 25 repetitions, the greatest difference was presented by the CSA (3.93%). For the lowest fitness value, the greatest difference was observed in the CSA (2.05%). Such differences were reflected in the maximum, average, and minimum efficiencies of each algorithm. The efficiencies were calculated based on the value of R$ 32,170,883 for the best solution, which was obtained using the branch-and-bound algorithm. However, such a value is not an optimal solution to this problem. The inclusion of a relatively good initial solution benefits the processing (Figure 1). For the CSA, all 25 repetitions with RLP inclusion were superior to the solutions obtained for processing without inclusion. For the VNS algorithm, only one repetition without the inclusion of the RLP was higher than the worst repetition with inclusion. For the GA, a metaheuristic that presented greater variability among the repetitions with the inclusion of RLP, the best solution without inclusion was superior to only eight replicates with inclusion. For SA, solutions with and without inclusions alternated, with no emphasis on one or the other methodology. The K-W test indicated differences (p < 0.05) between the use and non-use of RLP for the AG, CSA, and VNS algorithms (Table 2). FLORESTA,Curitiba,PR,v. 51,n. 3,jul Considering the evolution of the fitness values for each algorithm (Figure 3), the inclusion of RLP improved the initial behavior of the AG, SA, and VNS algorithms; however, after some iterations, the solutions evolved to those obtained without the inclusion of RLP. To achieve maximum fitness values based on less than 20 iterations (or generations), differentiation occurred throughout processing, causing the algorithms that considered RLP to yield better final solutions. Metaheuristic SA with the inclusion or absence of RLP demonstrated identical behaviors.

Maximum Fitness
Average

DISCUSSION
The application of LP to a problem that can be strictly regarded as typical of ILP does not generate integer values for all decision variables. This implies that not all stands have only one prescription, i.e., certain sites should be partitioned such that more than one management alternative is implemented. This implication is extremely inconvenient when one wishes to manage a forest, since the stand is assumed as the smallest management unit. Partitioning a stand causes serious operational and management problems, such as difficulty in limiting the harvesting area and in applying specific forestry tracts, as well as difficulty in mapping forest inventory attributes to each substand.
Hence, ILP or one of its variations, such as binary programming, rather than LP, may be considered. This is because the problem pertains to whether the harvest of a particular stand should be defined within 1 h of the planning horizon.
Nonetheless, ILP presents some difficulty when the problem belongs to the NP class. For the problem investigated, the optimal solution of ILP was not achieved even within 1 h of processing using the branch-andbound algorithm. The best solution obtained presented viability, i.e., it fulfilled all the restrictions imposed by the formulated ILP model.
When rounding the optimal solution provided by LP, an unfeasible solution is obtained, i.e., a deteriorated or degenerated solution. This occurs when at least one of the restrictions imposed on the solutions is not fulfilled, as in the case of non-compliance with wood demand in a specified period of the planning horizon. This is consistent to the discussion by Silva et al. (2003), thereby reinforcing the fact that rounding the optimal solution of LP can generate non-viable solutions to a forest production planning problem.
Although some authors perceived that the solutions obtained thus far are mediocre initial solutions for metaheuristics, the solution derived from the rounding of the solution via LP is still the best possible solution to the problem. Other methods for generating initial viable solutions can be applied, as discussed by Gaspar-Cunha et al. (2013); however, they are not applicable to this study.It is noteworthy that the inclusion of such a solution in the initial set of metaheuristic solutions significantly facilitated the identification of new viable solutions that were better than the initial ones. For the AG, CSA, and VNS, the best solutions were obtained without the inclusion of RLP, thereby justifying their use. The highest gain was yielded by metaheuristic GA. Meanwhile, the SA presented no contribution, and this was related exclusively to the algorithm used.
The SA considers the possibility of selecting a lower solution over an extremely good initial solution. This is due to the probability of acceptance of lower solutions, which is high at the beginning of processing and decreases gradually with the system temperature. Hence, such a result was anticipated, and it was assumed that the routines implemented in each algorithm were intrinsically associated with the problems being solved, i.e., for each class of problems, an algorithm exists that can yield a good solution in a relatively short time.

CONCLUSIONS
• The use of a solution derived from the rounding of a solution generated via an LP algorithm to solve the problem of forest production planning contributes significantly to the improvement in the results of the GA, CSA, and VNS. • Using an extremely good initial solution did not affect the results obtained by metaheuristic SA. • Rounding the solution to the relaxed problem is a good alternative for generating an initial solution for metaheuristics.