TUNING OF THE METAHEURISTIC VARIABLE NEIGHBORHOOD SEARCH FOR A FOREST PLANNING PROBLEM

In forest science it is important evaluate new technologies from computational science. This work aimed to test a different kind of metaheuristic called Variable Neighborhood Search in a forest planning problem. The management total area has 4.210 ha distributed in 120 stands in ages between 1 and 6 years old and site index since 22 m to 31 m. The problem was modelled considering the maximization of the net present value subject to the restrictions: annual cut volume between 140.000 m3 and 160.000 m3, harvester ages equal to 5, 6 or 7 years, and the impossibility of division of the management unity at harvester time. It was evaluated different settings for the Variable Neighborhood Search, varying the quantity of neighbours, the neighbourhood structure and number or generations. 30 repetitions were performed for each setting. The results were compared to the one obtained from integer linear programming and linear programming. The integer linear programming considered the best solution obtained after 1 hour of processing. The best setting to the Variable Neighborhood Search was 100 neighbours, a neighbourhood structure with changes in 1%, 2%, 3% and 4% of prescriptions and 500 iterations. The results shown by the Variable Neighborhood Search was 2,77% worse than one obtained by the integer linear programming with 1 hours of processing, and 2,84% worse than the linear programming. It is possible to conclude that the presented metaheuristic can be used satisfactorily in a resolution of forest scheduling problem when the best parameters are chosen. v.24 n.3 2018 TUNING OF THE METAHEURISTIC VARIABLE NEIGHBORHOOD SEARCH FOR A FOREST PLANNING PROBLEM


INTRODUCTION
Forest planning problems can be considered as complex in function of the amount of decision variables, imposed constraints and difficult to obtain all necessary information to develop work plans, mainly because of intrinsic characteristics of forest activities.In this sense, the effective courses-of-actions determination that maximize the enterprise economic returns or that optimize ecological and social aspects, as mentioned by Kaya et al. (2016), Ezquerro et al. (2016), Dong et al. (2016) and Shan et al. (2009), has been done with use of mathematical programming, like Linear Programming (LP), Integer Linear Programming (IP) and Mixed Integer Linear Programming (MIP) for instance (Troncoso et al., 2016).
Although their wide application, these tools had been replaced by more flexible techniques, mainly when there are more impeditive constraints in the models, like constraints of adjacency and singularity, or when nonlinear aspects are taking in account.It occurs because, with these considerations, the problem complexity is augmented and this can avoid a quick convergence to an optimal global solution in a feasible time when we use exact algorithms like branch-and-bound (B&B).
In these cases, when the problems became as belong to NP class (nondeterministic polynomial time), algorithms that can find a good feasible solution and that cannot require the full compliance with some assumption of classical mathematical programming, like additivity and proportionality (Jin et al., 2016), are indicated.This suggests the development and application of heuristics methods (Yoshimoto et al., 2016;Jin et al., 2016;Shan et al., 2009).
A heuristic is an iterative method that employs logic and rules that guide the search of feasible good solutions to a problem, near to optimal solutions but without ensuring the optimality (Kaya et al., 2016;Hillier and Lieberman, 2013;Jin et al., 2016;Ezquerro et al., 2016).Its application for combinatorial problems has been growing in the last years and already pass the number of scientific publications that consider the development and application of exact algorithms (Ezquerro et al., 2016).This occurs because sometimes the acceptance of a too much good solution obtained in a short time can be more interesting than the optimal solution obtained from a process with high computational efforts and time processing.
Heuristics and metaheuristics have been used in forest sector since 1980 (Jin et al., 2016;Ezquerro et al., 2016).Specifically, in forest management their applications have been highlighted as computational resources have been developed.It allows that new works could be done using these methods, mainly in bigger and complex management plans (Dong et al., 2016).In terms of forest production planning, Kangas et al. (2008) mentioned that the integer nature of forest planning problems and the use of spatial criteriums are important reasons by the increase of popularity of metaheuristics in these problems.
Among the vast number of already developed metaheuristics it is possible to highlight (Boussaid et al., 2013): Tabu Search (TS), Simulated Annealing (SA), Genetic Algorithm (GA), Ant Colony (AC), GRASP and HERO, Particle Swarm Optimization (PSO), Bee Colony (BC), Multi-Agent Systems (MAS) and Artificial Neural Networks (ANN).SA is the most cited methodology (Ezquerro et al., 2016), but considering the advances in studies about artificial intelligence, new metaheuristics arise and it become possible their evaluation in a wide of operation research applications.
Thus, algorithms relatively simple, as local search ones, had been widely used.However, in too many cases it is interesting that these algorithms can be modified or adapted to a specific problem in order to improve their performance.Among these algorithms, the Variable Neighborhood Search (VNS) has been highlighted by its simplicity, efficiency, and robustness in a wide of NP-Hard problems (Doerner et al., 2007).Its basic idea is to change the neighborhood structure to search a better solution for the problem (Affi et al., 2017;Doerner et al., 2007;Glover and Kochenberger, 2003).
There are some works already published about VNS as Affi et al. (2017) for vehicle routing problem, Amous et al. (2017) for capacitated vehicle routing problem, and Brimberg et al. (2017) for capacitated grouping problem, for instance.Despite that, this algorithm was not evaluated on forest production planning problem.In that case, it is important that studies be done in order to evaluate different set of algorithm's parameters for specific kind of problem, as suggested by Jin et al. (2016) and Shan et al. (2009), since the heuristic performance is highly dependent on the parameters used (Dong et al., 2016) and of the problem considered.
Thus, we evaluated the behavior of VNS metaheuristic in function of variations in parameters values and its relative efficiency comparing with Linear Programming and Integer Linear Programming in solving a forest planning problem.

Data
The present work was done with data from a eucalyptus forest with 4,210 ha divided into 120 management units (stands) and with ages between 1 and 6 years-old: 339 ha with 1-year-old; 768 ha with 2-year-old; 1,031 ha with 3-year-old; 601 ha with 4-yearold; 958 ha with 5-year-old; and 513 ha with 6-yearold.Wood production estimates for each stand in each age was obtained from the equation (1), where V i is the volume for the stand i, I i is the age of stand i, and S i is the site index for the stand i.It was developed by the authors considering data from a continuous forest inventory.Values of site index ranged from 22 to 31 m considering an index age equals to six-year-old.
Constraint (3) guarantees that all stand areas receives a determinate prescription.Constraints (4) and (5) limit the annual cutting volume between the minimum and maximum demand.Constraint (6), that is considered only in IP model, imposes a unique management alternative for each stand along the horizon plan. [1] Prescriptions for each stand considered harvester in one of three ages 5, 6 and 7-year-old, an immediately planting activity after the cut and a time horizon of the forest plan equals to sixteen years.It resulted in 81 management prescription for each stand totaling 9,720 decision variables for each model.Mathematical formulations for the optimization models were based on Model I mentioned by Jonhson and Scheurman (1977) in order to preserve the physical identity of management unit along the horizon plan.The objective was to maximize the net present value under constraints of annual demand between 140,000 m 3 and 160,000 m³.
The annual discount rate used was equals to 8 percent, the price of wood sales equals to R$ 80.00 per cubic meter and harvester cost equals to R$ 30.00 per cubic meter.Silvicultural costs ranged according to the age of each stand and were obtained from Binoti (2010): R$ 4,059.05ha -1 on first-year; R$ 1,627.81ha -1 on second-year; R$ 757.95 ha -1 on third-year; e R$ 88.12 ha -1 since on fourth-year of forest growth.

Linear and Linear Integer Programming
LP and IP models were formulated as suggest by Rodrigues et al. (2004).A unique difference between ones is the absence of constraint 6 for the LP model, where: GNPV is the global net present value for all the forest (2), in reais; C ij is NPV for stand i when assigned the prescription j, in reais; X ij is the decision variable and represents the proportion area of stand i that will be managed with prescription j; M is the total number of stands; N is the total number of different prescriptions for each stand; V ij(k) is the total volume of wood for the stand i, when assigned the prescription j, in the period k of planning horizon; Dmin k and Dmax k are minimum and maximum wood demand for the period k of the planning horizon. [2] [3] The LP and IP models were solved in CPLEX software version 12.7.1 (IBM Corporation, 2017) considering simplex and branch-and-bound algorithms, respectively.It was done on a computer with Windows 10, 64 bits, processor Intel Core i7 with 2.0 GHz and 8Gb of RAM memory.The LP model was used to test if the optimal solution was an integer solution.If this was true, neither IP nor metaheuristic solutions were necessary for our problem.

Variable Neighborhood Search
VNS is an algorithm based on a local search with different neighborhood structures (Doerner et al., 2007).They define a set of modifications that can be applied to a solution in order to create new solutions (Meignan et al., 2012).According to Mladenovic and Hansen (1997), the basic algorithm can be described as a set of steps that starts by the definition about the amount and types of neighborhood structures that will be considered.Then a random solution is obtained (Figure 1).With that, for each iteration and neighborhood structure, the algorithm executes a local search procedure in order to find a solution better than the considered before.If the algorithm does not find no better solution after n iterations, or after evaluates all neighbors, the structure is changed and the local search is done again.This process repeats until all considered neighborhood structures are used.The algorithm finishes when a stop criterium is reached, like as a number of iterations or time of processing.
This study considered a proposed VNS version with four neighborhood structures evaluated in two sets.This way, considering set 1 of neighborhood structure (S 1 ) and the fi rst local search procedure (LS 1 ), 1% of stands from one solution, randomly chosen, had its prescription modifi ed.Local search procedures LS 2 , LS 3 , and LS 4 considered modifi cations in 2%, 3% and 4% of stands.Set 2 of neighborhood structure (S 2 ) considered local search procedures that execute modifi cations in 10%, 20%, 30% and 40% of stands at a solution.Each stand modifi cation means to change a stand prescription (adopted for all plan horizon) that are chosen randomly for other prescription, also chosen in a random way.
The VNS algorithm evaluated all individuals from a defi ned neighborhood size in each local search procedure.Thus, to identify the infl uence of the number of neighbors, it was tested different sizes for the neighborhoods.It was considered neighborhoods with 1, 10, 30, 50 and 100 neighbors.
It was considered 30 repetitions and different stop criterium (20, 50, 100, 300 and 500 iterations) in each evaluation (arrangement of a set of neighborhood structure and a number of neighbors).The results were evaluated using the Kruskal-Wallis test (K-W) with 5% of probability.Confi guration with best values for average and maximum fi tness, as a loss value for standard deviation at end of processing, was chosen.
Constraints imposed on LP and IP models also were considered for VNS model.In that case, it was considered a penalty for the objective function in function of every broken constraint.Thus, for each volumetric demand constraint that was broken, objective function was decreased in R$ 100.00 per cubic meter in excess or lack at end of each period in the planning horizon.This method is the same adopted by Rodrigues et al. (2004).
Metaheuristic processing was done using MeP (Metaheuristics for Forest Planning) software.It was developed at Operation Research and Forest Modelling Laboratory at the Federal University of Minas Gerais in java language programming.

RESULTS
The VNS has its processing interrupted after reach the stop criterium and the processing time ranged according to the parameters values used (Table 1).The smaller processing time was obtained for the confi guration with 1 neighbor, set 2 of neighborhood structure and 20 iterations.The bigger processing time was obtained for 100 neighbors, set 1 of neighborhood structure and 500 iterations.
The parametrization that reached the best results considered a neighborhood with 100 neighbors, set 1 of neighborhood structure and 500 iterations (Table 2).The best solution was 231% superior to the worse solution found (which presented 1 neighbor, set 1 of neighborhood and 20 iterations).In relation to the average and maximum values for fi tness, the worse result was presented for the parametrization with 1 neighbor, set 2 of neighborhood structure and 20 iterations.
The increase in the number of iterations provided improvements on obtained results for all cases analyzed.However, there were not signifi cative differences (p<0.05) between the results obtained for 300 or 500 iterations according to the K-W test.The same occurred when the number of neighbors was increased and, in this case, there weren't differences between the results found with 50 or 100 neighbors considering the same test (p<0.05).Considering the different sets of neighborhood structure, the increase in the number of solution modifi cations in order to create a new neighbor  made the results worse and statistically different between the two sets according to the K-W test (p<0.05).
The finding of all test observed that increasing the number of neighbors made the algorithm faster in converge to the best solution with a less number of iterations (Figure 2).The change in the neighborhood structure from type 1 to type 2 decreased the algorithm performance, and it was necessary more iterations in order to reach same results.
The IP model solution was obtained after 1 hour of processing without the algorithm had converged to the global optimal solution.The solution found by the LP model was the global optimal (Table 3).The method that found the best solution was LP followed by the solution obtained by IP model.The best VNS solution was approximately 2.77% inferior to the IP value.
The solution obtained from LP model was not binary for all decision variables.It occurred in 25 ones
(0.26%).They were divided into 12 from 120 stands and one stand received three different management prescriptions and the other received two prescriptions.
For all methods, the annual harvested volume constraint was reached at best solution found.The period with more or less wood volume production ranged among the methods (Figure 3).Also, there were differences in the total volume produced and

DISCUSSION
Metaheuristics are a field of stochastics optimization which is a general class of algorithms and techniques that employ some randomness in order to find very good solutions for complex problems (Luke, 2009).When well designed, a heuristic method can be able to find an optimal solution for the problem.However, there is not a mathematical procedure that guarantee the optimality of any solution (Kaia et al., 2017).It is necessary to define the stopping criteria of any optimization algorithm.In this case, its definition can be done empirically, but this can cause a non-necessary use of computational resources without guarantee a solution significantly better than the one found before.
It is important to highlight that the results obtained in our work demonstrate that the increase in the number of iterations for VNS did not guarantee a significative improvement in relation to the solution that was found.If time is a determinant factor, as discussed by Jin et al. (2016), it is possible to consider 300 iterations instead 500.It will decrease the processing time near to 33%, although the average fitness decreases near to 0,3%.
The neighborhood size evaluated in each VNS algorithm iteration will define how much of the solutions space is verified during its processing.According to Hansen et al. (2017), a big amount of neighbors can increase the chances in find better solutions in the neighborhood.However, in cases where an initial solution is in a low promises region, evaluate too many neighbor solutions could not cause the expected effect.In this way, it is possible to spend computational resources nonnecessarily without a big return in terms of fitness value.
Considering the analysed problem, the evaluation of a wide of solutions in each iteration cause an increment in the processing time, from 119.99 to 235.23 seconds on average.There wasn't a statistical difference (p>0,05) between the results obtained with 50 and 100 neighbors.It allows us to indicate a less amount of neighbors.It saves processing time without losses in relation to the best solution obtained.In another hand, the use of few neighbors, like 1 and 10, can make the algorithm unable to escape from local optima.Indeed, Costa et al. (2017) mentioned that VNS try to explore a solution neighborhood aiming to find a better path to reach the optimal solution for the problem.If this neighborhood is restricted to a few individuals, the search procedure could not offer an expected effect.
The best configuration in terms of processing time is that considers 50 neighbors and 300 iterations.This option resulted in a maximum fitness equals to R$ 30,249,780, which is only 0.27% smaller than the best solution found considering 100 neighbors and 500 iterations.
In terms of neighborhood structure, when it was consider modifications in prescriptions for 10%, 20%, 30% and 40% of 120 stands, the algorithm lost performance.It can be explained by the high rate of changes at each solution, that can create new solutions near to randomness.Thus, the search process can be displaced to regions few interesting of solutions space and finding solutions with lower values for fitness function.Dong et al. (2016) observed an increment in the processing time for the SA when compared a local search procedure with modification of only one prescription for a management unit with the one with changes in two management units.In that case, the second option was 4 times slower the first one.
We already not expect that a deterministic algorithm like B&B was able to solve the problem, in a feasible time, because of its combinatorial aspect.Therefore, the process was interrupted after 1 hour.This procedure also was adopted by Caro et al. (2003) because the B&B algorithm not return an optimal solution in a time smaller than 44 hours.This suggests the necessity of development of search methods that can obtain very good solutions in a feasible time, as the case of the metaheuristics.
The fact that different parametrizations of VNS presented values close to ones presented by the IP solution in a short processing time can put the algorithm as too much interesting to solve forest planning problems.In this sense, Jin et al. (2016) mentioned that, if the search time is a determinant factor, there is the necessity of identifying a heuristic procedure that can generate good solutions quickly, as the case of SA and TS algorithms.This is important when there is a demand for distinct forest planning scenarios or alternative work plans, as the ones that consider variations on annual demand of wood, on discount rates used, and costs and revenues values, for instance.

FIGURE 1
FIGURE 1Flowchart of processing routine of VNS.

TABLE 1
Processing time obtained from evaluations considering different parameters of VNS.

TABLE 2
Processing time obtained from evaluations considering different parameters of VNS.

TABLE 3
Best solutions results founded by each method and its proportion in relation to linear programming (LP) and integer linear programming (ILP).Evolution of solutions average for each iteration of VNS considering different settings.
FIGURE 3 Annual cut volume considering the best solutions founded by linear programming (LP), integer linear programming (ILP) and VNS.The dashed lines indicate the values of maximum and minimum annual demand.