OPTIMIZATION MODEL OF SELECTIVE CUTTING FOR TIMBER HARVEST PLANNING BASED ON A SIMULATED ANNEALING APPROACH

An optimization model for Timber Harvest Planning ( THP) is used to form harvest areas with the objecti v of maximizing the tree harvest volume subject to ha rvest regulations. Most optimization models for THP currently in use are based on clear cutting and res tricted to adjacency constraints to prevent the cre ation of large clear-cut openings in the forest. However, an optimization model based on selective cuttingthe cutting technique utilized in tropical countrie s-has rarely been described in the literature. The aim of this study was to propose an optimization model for THP based on selective cutting and subject to a maximum number of trees to be harvested and a minim um number of trees to be damaged during each planning period. The model was solved using three o ptimization techniques to identify the most suitable technique for use in the process of harves t area formation: Monte Carlo Programming (MCP), Simulated Annealing (SA) and Threshold Acceptance ( TA). The obtained results indicate that the SA method provides better solutions than the MCP and T A methods, regardless of problem size. As a conclusion, the proposed model provides a tool to g enerate timber harvest plans with the improved monitoring and control techniques used in harvestin g operations in tropical countries.


INTRODUCTION
One crucial problem of timber harvest planning involves selecting and grouping a number of harvest blocks to form harvest areas. This problem requires managing a large number of harvest blocks and addressing conflicts in management objectives while simultaneously addressing several inherent constraints (Walters, 1997). Optimization models have been used to determine the optimum combination of blocks to form the harvest areas to be cut in different time periods. The objective of the model is to maximize the harvest volume subject to spatial constraints that protect the non-timber value of a forest, such as biodiversity conservation and wildlife habitat protection. These spatial constraints include the following: Patch shape and size, old forest issues, site sensitivity and adjacency constraints (Walters, 1997;Vielma et al., 2007).

Optimization Models for THP
Optimization models for THP formulate harvest areas with the objective of maximizing harvest volume subject to spatial constraints that value the non-timber interests of a forest. Most of the previous models are based on a clear-cutting technique, restricted to adjacency constraints and solved using optimization techniques. However, the development of an optimization model for THP based on the selective cutting technique utilized in tropical countries has rarely Science Publications JCS been reported in the literature. Thus, the aims of this study were to propose an optimization model based on selective cutting that reflects the cutting technique utilized in tropical countries and to test three optimization techniques and identify the technique most suitable for solving this problem.

Challenges
As reported in the literature, earlier THP models that were developed based on a LP approach, such as SPECTRUM by Grear and Meneghin (1997), focused on maximizing timber production without spatial constraints. However, because maintaining other forest values has become increasingly important, various IP-and MIPbased models incorporating spatial constraints have been reported in the literature (Goycoolea and Murray, 2005;Vielma et al., 2007;Constantino et al., 2008;Ohman and Wikstrom, 2008). Nevertheless, the performance of these solution techniques is restricted by the complexity of the problem. The complexity of the THP problem depends on the size of the problem and the number of constraints imposed in the model (Lockwood and Moore, 1993). Thus, many researchers have proposed the used of metaheuristic techniques to solve large-sized problems incorporating several spatial constraints.
Meta-heuristic optimization techniques are used to provide feasible solutions in large-size and multi-objective THP problems. However, the limitation of these solution techniques is that they do not guarantee optimal solutions (Baskent and Jordan, 2002). For example, MCP implements a basic random search known as a "blind search". SA is based on the heating process of crystalline materials at extreme temperatures and the process of slowly cooling the materials to minimise damage (Kirkpatrick et al., 1983). Lockwood and Moore (1993) first proposed the use of SA to solve the THP problem. TA, on the other hand, seeks to improve the SA technique by utilising an acceptance condition of poor solutions based on a pre-defined threshold value (Dueck and Scheuer, 1990), rather than calculating the probability based on the current temperature, as is the case with SA.
TS uses memory structures called tabu lists to store potential solutions to avoid repeating the same processes, distinguishing TS from other meta-heuristic techniques. TS is essentially deterministic because in tabu lists, previous search information is used to control the process of improving the current solution to avoid becoming trapped in local optimization problems (Glover and Laguna, 1997). GA, which was introduced by Holland (1975), is based on an analogy to the process of biological reproduction. The GA approach can be succinctly described as an initial population generated at random and used as the basis of the formation of the next population, which is created by changing the combination of information in the base population.

Related Work
Various optimization models of THP have been reported in the literature. Most of the models are based on clear cutting and impose adjacency constraints, such as maximum opening size and green-up delay. The maximum opening size constraint is used to avoid creating large open areas in a forest with the clearcutting approach. Thus, adjacent harvest areas are protected from being selected for harvesting during the planning period. The green-up delay constraint, on the other hand, mandates the duration required for trees in the selected harvest areas to grow and reach a minimum height before neighboring areas are selected for the next harvest Ohman and Lamas, 2003;Crowe et al., 2003;Goycoolea and Murray, 2005;Vielma et al., 2007;Bettinger and Jianping 2008;Konoshima et al., 2011;Nora and Toth, 2013). The implementation of adjacency constraints aims to protect living species and their habitats from the various impacts of clear cutting on forest landscapes, such as erosion and deterioration in water quality (Weintraub et al., 2000).
Many studies have proposed the used of optimization techniques based on mathematical programming and meta-heuristic approaches to solve the models. Mathematical programming techniques, such as Linear Programming (LP), are currently used to generate optimal solutions in large forest areas without spatial constraints, while Integer Programming (IP) and Mixed Integer Programming (MIP) have been proposed to solve THP models in smaller forest areas that incorporate spatial constraints.

Model Formulation
The model proposed here is based on a selective cutting technique in which only trees that meet predetermined criteria are felled. These criteria include a Minimum Tree Diameter (MinTreeD) to ensure that only mature trees are selected for harvesting, a minimum number

JCS
of residual trees for future stock in each block (MinFStock) to ensure the continuous production of timber, a minimum number of trees to be cut in each block to ensure economic harvesting (EcoH) and a maximum number of trees permitted to be harvested per planning period (AwHVol) to avoid excessive logging in the specified period. The model can be divided into two phases: The pre-processing phase and the harvest area formation phase.

Phase i. Pre-Processing Phase
The objectives of the pre-processing phase include estimating the number of Damaged trees (DmgT), calculating the number of residual trees with a diameter of at least 30 cm diameter at breast height (dbh) for Future Stock (FStock) and calculating the amount of the Harvest Volume (HVol) based on the total number of trees that can be cut (T2Cut) in each harvest block.
In this study, the pre-processing phase can be divided into two sub-processes: (i) estimating the number of trees that can be cut and calculating the number of damaged trees and (ii) identifying trees to be excluded from the harvest if the number of residual trees for future stock is lower than the minimum requirement (the tree-saving approach). The processes involves in the pre-processing phase are illustrated in Fig. 1.
In sub-process (i), each tree with a diameter of at least 45 cm dbh is selected as a potential tree to be cut (PT2Cut). The analysis performed on each of these trees involves estimating the affected area based on the tree felling direction, identifying a list of trees that are located in the affected area, estimating the number of damaged trees based on the distance from the potential tree to be cut and estimating the number of residual trees with a diameter at least 30 cm dbh for future stock. If the number for future stock is sufficient, all the trees that have been selected for harvesting will be felled and the number of trees in the harvest volume can be estimated. Otherwise, we proceed to sub-process (ii).
In sub-process (ii), the analysis identifies any potential tree to be cut that causes the highest number of damaged trees to be excluded from the harvest (T2Save). However, this process decreases the harvest volume and increases the number of residual trees. The same process is repeated until the number of trees for future stock is sufficient. If all the potential trees to be cut have been excluded from the harvest, but the number of trees for future stock is still lower than the minimum requirement, then the block would be excluded from the harvest. The outputs generated from this process include the number of trees in the harvest volume and the number of damaged trees in each harvest block. These data become inputs in the harvest area formation process. The model formulations for each task may be described as follows.

Damaged Trees Estimation
The number of damaged trees is estimated using the following steps: • Determine the affected area based on the cutting degree (CutDg): Affected area=1 if CutDg ≥ 0 and CutDg<90 (1) Affected area = 2 if CutDg≥90 and CutDg<180 (2) Affected area = 3 if CutDg≥180 and CutDg<270 (3) Affected area = 4 if CutDg ≥ 270 and CutDg <360 (4) The affected area is divided into four areas. The affected area is assigned to 1 if the cut degree is in between 0 to 89 degree (Equation 1). The affected area is assigned to 2 if the cut degree is in between 90 to 179 degree (Equation 2). The affected area is assigned to 3 if the cut degree is in between 180 to 269 degree (Equation 3). The affected area is assigned to 4 if the cut degree is in between 270 to 359 degree (Equation 4).
• Determine the tree felling direction: The tree felling direction is determined based on the affected area. Equation (5) represents the formulation to estimate the tree felling direction in affected areas 1 and 3. While Equation (6) represents the formulation to estimate the tree felling direction in affected areas 2 and 4. • Calculate the distance (Dist) of the neighboring trees from the potential trees to be cut: Where: (Dx0, Dy0) = The location of the potential tree to be cut (Dx1, Dy1) = The location of neighboring trees in the affected area  (7) represents the formulation to calculate the distance of the neighboring trees in order to identify the potential damaged trees.
• Calculate the number of damaged trees: where, stem height refers to the height of the potential tree to be cut. Equation (8) represents the formulation to calculate the damaged tree. The neighboring tree will be marked as a damaged tree, if the distance is less than the stem height of the potential tree to be cut.

Future Stock Estimation
Future stock refers to the residual trees with diameters in the range of 30 to 44 cm dbh. The formulation used to calculate the number of trees for future stock in each harvest block may be depicted as follows: where, i refers to the block number and y refers to the residual trees with diameters between 30-44 cm. The future stock in each harvest block is equal to the total number of the residual trees with diameter in between 30-44 cm (Equation 9).

Harvest Volume Estimation
Harvest volume refers to the total number of trees to be cut in each block. The harvest volume in each harvest block is estimated based on the number of trees for future stock, as follows: • If the number of FStock is greater than the minimum requirement, the formulation to calculate the harvest volume is: where, |t2cut| refers to the total number of trees allowed to be cut in each block. Thus, the harvest volume is equal to the total number of trees to be cut (Equation 10).

JCS
• If the number of FStock is smaller than the minimum requirement, the formulation to calculate the harvest volume is: where, |pt2cut| refers to the total number of potential trees to be cut and |t2save| refers to the total number of trees to be excluded from the harvest. Thus, the harvest volume is equal to the remaining number of potential trees to be cut (Equation 11). • If the number of FStock is smaller than the minimum requirement and all the PT2Cut have been excluded from the harvest, the formulation to calculate the harvest volume is: The harvest is equal to zero. This means that the block is excluded from the harvest (Equation 12).

Phase II. Harvest area Formation Phase
The objective of the harvest area formation phase is to combine a number of blocks to form the harvest area with the greatest harvest volume and the lowest number of damaged trees that adheres to the restrictions of the minimum and maximum number of trees that can be cut in each block and in a single planning period. The process involved in forming the harvest area is illustrated in Fig. 2.
For example, if the variable EcoH is set to seven and the variable MaxHVol is set to 25, the harvest areas may be formed by selecting and combining any block in which the HVol is equal to or greater than seven. Thus, three harvest areas can be generated, with area_1 consisting of blocks

JCS
The symbol Q in Equation (13) represents the objective function that attempts to maximize the total number of trees in the harvest volume subject to constraints (14-16). Constraint in Equation (14) represents the minimum number of trees to be cut in each harvest block. The harvest volume in each block must be greater than the minimum requirement to ensure for economic harvesting. Constraint in Equation (15) represents the maximum number of trees that can be cut in one planning period. The total number of harvest volume in one planning period is at most equal to the maximum number of trees permitted to be harvested per planning period to avoid excessive logging in the specified period. Constraint in Equation (16) represents the lowest number of damaged trees. The selected harvest area for harvest has the lowest number of damaged trees. The variable X in Equation (17) is used as a decision variable. The value of X is equal to 1 if the block is allowed for harvesting; otherwise, the value of X is equal to 0.

Data
The proposed model was tested using actual and hypothetical data. The actual data were obtained from the reserved forest of Bintang Hijau in Kuala Kangsar. The forest area was divided into smaller harvest blocks. Seven harvest blocks were used, comprising a total of 636 trees with an average of 90 trees in each block. A total of 165 trees were categorized as potential trees to be cut with an average of 22 trees in each block. A total of 139 trees were categorized as protected trees, with an average of 20 trees in each block. The remaining trees which had diameters below 30 cm dbh, were categorized as normal trees. The average number of normal trees in each block was 47 trees. Meanwhile, 30 blocks of hypothetical data were randomly generated to represent a larger problem.
In this study, two experiments were conducted with the following objectives: (i) to form a harvest area based on maximum HVol and the smallest number of damaged trees and (ii) to identify a suitable optimization technique to be used in the harvest area formation process and tested using actual and hypothetical data types. Three optimization techniques were used to solve the harvest area formation problem: MCP, SA and TA.
The best solution is gained through the iteration process. In the MCP solution approach, the number of iterations is controlled by a variable t, whereas in the SA approach, the number of iterations is controlled using temperature and the variables involved are starttemperature, reduction-rate and end_temperature. The implementation of the TA method in this study was conducted based on the SA technique and the number of iterations is, therefore, controlled using temperature, with an additional variable to control the threshold value.
The value for each variable was determined based on the results from several trial-and-error experiments. Therefore, the variable t was set to 145 and the variables start_temperature, reduction_rate and end_temperature were set to 10000, 0.09 and 0.001, respectively, for both the SA and TA techniques. The threshold value was set at 10% of the AwHVol and the threshold value for the small-sized problem was set to-10. Meanwhile, the threshold value was set to -200 for the large-size problem.
The implementation of the proposed solutions required seven additional variables, including curr_HVol, curr_DmgT, new_HVol, new_DmgT, MxHVol, MnDmgT and X. The variables curr_HVol and curr_DmgT were used for the current values of HVol and DmgT, respectively. The variables new_HVol and new_DmgT were used for the new values of HVol and DmgT. Meanwhile, the variables MxHVol and MnDmgT were used for the largest value of HVol and the smallest value of DmgT.

Proposed Solution Approaches
A harvest area is formed by combining a number of randomly generated blocks. The process includes calculating the number of trees to be cut (HVol) in each block. The process is repeated until the HVol is at most equal to the maximum number of trees allowed to be cut in one planning period (AwHVol). The proposed solutions were used to generate a harvest area with the highest HVol subject to the predefined AwHVol and with the lowest number of damaged trees (DmgT).
In the proposed approaches, an initial harvest area is formed by selecting blocks in which HVol is at least equal to the predefined value of EcoH. The total HVol and DmgT are determined by calculating the number of trees to be cut and the number of damaged trees. The process will be repeated until the total HVol is at most equal to AwHVol. The algorithm used to generate the initial solution is shown below. Generate initial solution: Set the value of variables EcoH, AwHVol, X Randomly form an initial harvest area Meanwhile, the current solution is obtained by comparing the new solution and the initial solution. If the new solution is better than the initial solution, the new solution will be accepted as the current solution to replace the initial solution. Otherwise, the initial solution remains the current solution. In the SA approach, however, a poor solution with a certain probability will replace the initial solution.
In the TA approach, the acceptance of a poor solution is based on the predefined threshold value. In the proposed approaches, a new solution will be accepted to replace the current solution if and only if, the difference between new_DmgT and curr_DmgT is less than zero and the difference between new_HVol and curr_HVol is equal or greater than zero.
In addition, the best solution is obtained by comparing the values of curr_HVol and curr_DmgT with the values of MxHVol and MnDmgT, respectively. The current solution is accepted as a best solution if the value of curr_HVol is greater than or equal to the value of MxHVol and the value of curr_DmgT is less than the value of MnDmgT. In the MCP approach, however, this comparison is not required because the current solution represents the best solution and the current solution will, therefore, be reported as the best solution. The algorithms for the MCP, SA and TA techniques are as follows:

RESULTS
This section shows the results obtained from the experiment that attempted to identify a suitable optimization technique to solve the harvest area formation problem.

Harvest Area
The experiment was repeated 15 times. Table 1 presents the sample results obtained from the harvest area formation process using the MCP algorithm. The results show that two harvest areas are generated. The harvest areas are shown in Fig. 3. Table 2 shows the performance of each algorithm with the data involving seven (7) blocks and 30 blocks, in term of the percentages of the solutions generated in each category. The algorithms performance comparisons based on harvest volume and damaged tree for both categories are presented in Fig. 4 and 5 respectively.  The percentage of damaged trees in the category of greater than 10% above average for the MCP, SA and TA techniques are 8, 12 and 13%, respectively. Thus, the SA technique is a better method than the MCP or TA techniques for use in solving the harvest area formation problem using the proposed model, regardless of data size.

DISCUSSION
The combination of blocks that form the harvest area is generated based on the maximum HVol and the smallest number of damaged trees. The harvest area is formed by combining a number of blocks from a sequence of blocks that is randomly generated in every process cycle. The results show that two harvest areas are generated: Area_1, consisting of a combination of blocks 1, 2, 3 and 4; and area_2, consisting of a combination of blocks 1, 3, 4 and 6. The maximum numbers of HVol in the respective harvest areas are 100 and 99 and the number of damaged trees is 121 in area_1 and 129 in area_2.
In this study, the performance of the algorithms is evaluated based on the following two factors: (i) the ability to meet the maximum number of trees allowed for harvesting in one planning period and (ii) the ability to minimise the damage to residual trees. The solutions obtained are classified into the following three categories: (a) an "Optimal" solution, in which the total number of HVol is equal to the MxHVol or the total number of damaged trees is below average; (b) a "Reasonable" solution, in which the total number of HVol is within 5% below the optimal solution or the total number of damaged trees is within 10% above average; and (c) a "Bad" solution, in which the total number of HVol is greater than 5% below the optimal solution or the total number of damaged trees is greater than 10% above average.
After 100 runs using a small-sized problem, both the MCP and SA techniques provide 100% optimal solutions compared to 99% with the TA technique. Thus, the MCP and SA techniques are better at forming the harvest area in a small-sized problem. With data involving 30 blocks, the solutions obtained from the SA technique better maximize the total number of trees in the harvest volume than the solutions obtained with the MCP and TA techniques. For example, 16% of the solutions generated using the SA technique is optimal, compared to 10 with the TA technique and 9% with the MCP technique.
However, in terms of minimizing the number of damaged trees, it appears that the MCP technique generated better results. The percentage of damaged trees that is below average generated from the MCP technique is higher than in the results obtained using the SA and TA techniques. The values are 61, 52 and 57%, respectively. However, as shown in column Science Publications

JCS
DmgT under the category of >10% above average, the SA technique obtains the lowest percentage of damaged trees compared to the MCP and TA techniques.

CONCLUSION
In this study, we proposed an optimization model for THP based on selective cutting to reflect the cutting technique utilized in tropical countries. The proposed model is used to form the harvest areas with the objective of maximizing the harvest volume, subject to the maximum number of trees allowed for harvesting and the minimum number of damaged trees per planning period. In addition, the model implements an approach to exclude one or more potential trees to be cut from the harvest to fulfil the minimum number of future stock required for timber sustainability. The minimum number of trees to be cut is also imposed in each block to ensure economical harvesting; thus, the harvest areas are formed by combining and selecting a number of blocks in which the harvest volume is at most equal with the maximum number of trees allowed for harvesting in one planning period.
Three optimization techniques have been developed and tested to identify the most suitable technique for use in solving the problem of harvest area formation. The obtained results indicate that the SA technique provides better solutions than the MCP and TA techniques, regardless of the size of the problem. For future research, a component to predict tree growth may be incorporated into the model to form the harvest areas and generate harvest plans for several planning periods simultaneously.