International Journal of Forestry Research

Volume 2015 (2015), Article ID 173042, 15 pages

http://dx.doi.org/10.1155/2015/173042

## An Evaluation of the Use of Simulated Annealing to Optimize Thinning Rates for Single Even-Aged Stands

Interdisciplinary Graduate School of Science and Technology, Shinshu University, 8304 Minamiminowa, Kami-ina-gun, Nagano 399-4511, Japan

Received 18 June 2015; Revised 27 October 2015; Accepted 10 November 2015

Academic Editor: Ignacio García-González

Copyright © 2015 Kai Moriguchi et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

We evaluated the potential of simulated annealing as a reliable method for optimizing thinning rates for single even-aged stands. Four types of yield models were used as benchmark models to examine the algorithm’s versatility. Thinning rate, which was constrained to 0–50% every 5 years at stand ages of 10–45 years, was optimized to maximize the net present value for one fixed rotation term (50 years). The best parameters for the simulated annealing were chosen from 11^{3} patterns, using the mean of the net present value from 39 runs to ensure the best performance. We compared the solutions with those from coarse full enumeration to evaluate the method’s reliability and with 39 runs of random search to evaluate its efficiency. In contrast to random search, the best run of simulated annealing for each of the four yield models resulted in a better solution than coarse full enumeration. However, variations in the objective function for two yield models obtained with simulated annealing were significantly larger than those of random search. In conclusion, simulated annealing with optimized parameters is more efficient for optimizing thinning rates than random search. However, it is necessary to execute multiple runs to obtain reliable solutions.

#### 1. Introduction

Optimizing the thinning rate for single stands is a fundamental problem in improving the economy of forest management. In Japan, the proportion of small-scale forest holders is large. For example, 65.4% of private forest owners in Nagano Prefecture own areas < 1 ha in extent, and 90.6% own areas smaller than 5 ha [1]. They therefore have little leeway in performing spatial optimization, and it is important to optimize the thinning rates applied to single stands to improve profitability. Fortunately, recent developments in computing enable rapid quantitative optimization with no experimental knowledge.

The problem may be optimized in a straightforward manner using full enumeration (FE), that is, enumerating and examining all patterns. Because calculation cost increases exponentially with the number of control variables (referred to as “the curse of dimensionality”), this however tends to be impractical. Approaches that reduce the computing costs are therefore required. Many algorithms for solving this problem have been proposed and can be grouped as follows: dynamic programming (DP) [2–6]; nonlinear programming (NP) [7, 8]; and heuristic methods such as random search (RS) [9], tabu search, or greedy algorithms [10]. Of these, DP and NP have been widely used and have been further developed based on comparison with each other.

Roise [7] reported that two NP approaches, the methods of Hooke and Jeeves [8] and Powell [11], were more efficient than the discrete DP algorithm which served as a baseline comparison. On the other hand, Paredes et al. [4] proposed a new DP algorithm named PATH (projection alternative technique) that was more efficient than traditional DP algorithms. Yoshimoto et al. [3] developed an extended PATH algorithm named MSPATH (multistage PATH), which can reflect long-term effects such as an increase in log price with diameter and showed that it could provide better solutions than PATH [6, 12]. Yoshimoto et al. [6] developed an improved PATH algorithm named RLS-PATH (region-limiting strategies PATH) to deal with yield models including many variables and reported that in many cases it provided better solutions than Hooke and Jeeves’ method. Each of these methods has been used to optimize thinning rates for single stands (e.g., NP: [13–21], DP: [22–32]). In particular, the MSPATH algorithm is more appropriate than PATH when long-term effects caused by competition-density effects or fluctuations of log price with diameter are taken into account [12]. Yoshimoto [12] developed a way to apply it to yield models based on stand density management diagram (SDMD) [33], a type of forest growth model incorporating competition-density effects that is commonly used in Japan, which has been applied in several studies [27–32]. However, it is necessary to consider carefully whether such optimized results represent a “better choice” or “the best choice.”

With respect to the spatial planning of forests, it has been argued that it is important to validate optimization algorithms [34]. This also applies to algorithms for optimizing thinning rates for single stands, but existing algorithms have as yet only been compared with each other, not against exact optimal solutions as provided by FE. Producing a solution that is validated as close to the optimal solution is potentially not always essential, if the user is only using the optimization algorithm as a way to improve the current situation; in other words, for them “optimized” means “better choice.” However, in some cases a solution close to the optimum is needed. For example, it is difficult to manage forests without subsidy in Japan. In this case, the result of a stand-level optimization analysis can determine whether a stand can be managed successfully under optimized conditions and whether it is worth subsidizing. Such an analysis can provide important information for regional zoning, and in this case a high degree of reliability, that is, confidence that the method provides solutions close to “the best choice,” is necessary to satisfy accountability requirements, especially towards forest owners who have applied for rezoning and whose forest is not subsidized. The same considerations may apply to countries or regions in situations similar to Japan, that is, little leeway in performing spatial optimization and a limited supply of funding.

In a previous study [35], RS with 10^{6} iterations and MSPATH were compared with coarse FE (CFE), by splitting variables into a lattice, comparing them using two simple SDMD-base small-scale yield models (differing only in their price functions) as benchmark, and setting net present value (NPV) for 50 years as the objective function. The NPV of the nonmonotonic pricing yield model provided by MSPATH was 60,419 yen/ha less than that derived using CFE, and the monotonic pricing yield model provided by MSPATH was 30,920 yen/ha less. In contrast, those provided by RS were 3,755 yen/ha and 5,487 yen/ha less, respectively, than that derived using CFE. On the other hand, MSPATH reduced the calculation cost by 10^{−5.2}~10^{−5.1} times relative to CFE and 10^{−2.6}~10^{−2.5} times relative to RS. MSPATH therefore confers a great advantage in terms of calculation cost, but it is not possible to rely on it providing “the best choice.” Thus we need to establish a more reliable method of optimizing thinning rates for single stands.

There are diverse thinning rate models for single stands as a result of differences in the complexity of the models (e.g., whether they consider variation in trees [26, 27] or in logging and lumbering [31]), in the objective function used (e.g., NPV [10, 12, 29, 32] or in soil expectation value (SEV) [14, 36]), or in other criteria [26, 27]. It is difficult to deal with the necessary complexity using NP algorithms, since the various trade-off relationships, for example, between yield volume at thinning and that at final cutting, may require multimodality of yield models, and the incorporation of logging processes makes them discontinuous. On the other hand, theoretical proofs are required to ensure the reliability of the DP algorithm; however it may be difficult to obtain sufficient proof to support complex yield models. We therefore suggest that simulated annealing (SA) [37, 38], a stochastic algorithm inspired by the annealing process of metals, may satisfy the requirements. It is worth noting that SA has been theoretically proven to provide the global optimum for a wide variety of problems [38], such as the traveling salesman problem [37, 38] and the knapsack problem [39], as long as the “temperature” parameter is decreased gradually enough.

In the field of forestry, SA has usually been used for forest planning under spatial constraints [40–42] or with a special objective function such as that defined based on “forest spatial value” [43, 44]. Lately the performance of SA itself has been analyzed, for example, its sensitivity to the parameters used [45] and its application in neighborhood methods [46]. Thus, SA is commonly used in the field and its performance is evaluated. However, this approach has never been used to optimize thinning rates for single stands and its performance in this area has not been evaluated. It is therefore necessary first to evaluate whether SA can be used reliably to optimize thinning rates for single stands.

In this study, we evaluate the potential performance of SA in optimizing thinning rates for single even-aged stands, using plural yield models based on SDMD. We set the thinning interval in each model at five years and the objective function as the NPV before planting for one rotation term. The age of the rotation is one of the most important variables for improving the financial viability of stand-level forest management. However, we fixed this variable in our model, because a change in rotation age leads to a change in the number of dimensions of the model (the number of variables, in this case the number of thinning ages). To optimize both rotation age and thinning rates, it would thus be necessary to define several models with different rotation ages, optimize them, and then choose the best one. The searching performance of the model with a fixed rotation age therefore determines performance in optimizing a model that treats rotation age as a variable.

Whereas the “performance” of an optimization method includes both reliability and efficiency, we must first improve its reliability and only then enhance its efficiency while maintaining reliability. The emphasis is therefore on securing reliable optima with proven, orthodox SA methods. However, the combination of parameters used has strong effects on the performance of SA approaches [45]. In addition, SA is a stochastic algorithm, and as a result the solution obtained from each searching process varies. We therefore search for “the best” combination of SA parameters using CFE, to ensure the best possible performance, using the mean NPV from 39 runs of the SA as the objective function. We also compare the 39 NPVs from the solutions obtained using optimized SA parameters with those obtained using CFE to evaluate whether SA can provide solutions close to the global optimum. Finally, we statistically compared these values with those obtained using RS with five times iterations, to evaluate the comparative merit of SA.

#### 2. Materials and Methods

##### 2.1. Benchmark Model

A hypothetical even-aged stand of Japanese larch (*Larix kaempferi*) in Nagano Prefecture with medium site quality and a density of 2500 planted trees/ha was used as a benchmark stand. Minimum thinning age was fixed in 10 years and the thinning interval in five years. The rotation age was fixed in 50 years, since this value would still allow application of CFE.

###### 2.1.1. Objective Function

We defined the objective function as NPV before planting for one rotation term:where is the income from the th cutting (yen/ha), is the annual interest (%), is the stand age at the th cut (year), and is the cost of regeneration (yen/ha). Annual interest was set at 0.9% [47]. As mentioned above, values were fixed to increments of five years over a range of 10–50 years.

###### 2.1.2. Growth Model

We used SDMD to define the growth model. The formulae for larch in Nagano Prefecture are as follows [48]:where is the stand stem volume per area (m^{3}/ha), is the mean height of dominant trees (m), is the present number of trees per area (trees/ha), is the stand form height (m), is the basal area (m^{2}/ha), is the diameter at breast height (DBH) of the mean basal area (cm), is the mean DBH (cm), is the number of trees per area on the full-density curve (trees/ha), is the stand volume per area on the full density curve (m^{3}/ha), and Ry is the relative yield index (the standard density index for the SDMD). and are parameters that change with the mean height of the dominant trees and are defined as follows:

Note that these formulae are common to the aforementioned previous study [35], but the aspects of the model described in the rest of this section and in the next one have been modified. The present model is more appropriate because it simulates self-thinning after the first thinning.

We modeled self-thinning using Tadaki’s model [49], calculated as follows [50]:where is “initial number of trees per area” (trees/ha), a parameter of the self-thinning curve. is given as follows:

Equation (11) is the self-thinning curve before the stand reaches the full-density condition. After that, self-thinning proceeds according to (12), the full-density curve, which takes the mean height of the dominant trees as a variable. is what the number of trees per area at age 0 would be if the stand had only been subject to self-thinning but not to additional thinning by humans. is therefore only ever equal to the planted number of trees per area if the forest is never thinned. We call this parameter to avoid confusion with planted number of trees per area and initial values of optimization.

SDMD is a function of the mean height of the dominant trees and the number of trees per area. It hypothesizes that only the number of trees per area, not the mean height of the dominant trees, is affected by thinning. In other words, SDMD simulates thinning by translating the self-thinning curve under a specific mean height of dominant trees. The value of corresponds to the self-thinning curve in (11). As the relationship between the number of trees per area under a specific mean height of dominant trees and is monotonic, we can simulate thinning by decreasing .

It is easy to calculate the number of trees per area based on under an arbitrary mean height of dominant trees using (11) and (12). However, if number of trees per area at an arbitrary mean height of dominant trees is given, we first need to identify with a numerical approximation with (11) to simulate self-thinning after that age. Moreover, after the stand has reached full-density condition, we cannot compute using (12) and given the number of trees per area, because the full-density curve as (12) is independent of . To address these problems, we define thinning rate as the ratio of the difference between before and after thinning as follows:where is the thinning rate (%), is before the thinning, and is after the thinning.

To simulate the growth of the stand with age using SDMD, the mean height growth curve of the dominant trees is required. The growth curve of Japanese larch at an average site in Nagano Prefecture can be predicted using the following formula [51]:where is the stand age (years). Since SDMD assumes lower thinning, the mean height of the dominant trees will not be affected by thinning.

###### 2.1.3. Volume and Diameter of Yield Trees

The total stem volume of cut trees per area is calculated as follows:where is the cut stem volume per area when trees per area are cut at time , is the standing tree volume per area calculated using (2) and (15), where the number of trees per area is at time , is the number of cut trees per area, is the number of standing trees per area, and is the total number of standing and cut trees per area. This equation can be applied to both thinning and final cutting. Mean stem volume per tree can be calculated by dividing the left side of (16) by . The mean DBH of the cut trees was calculated using the formula:where the subscripts of are the same as those of in (16). SDMD does not describe variation in stem volume, so we assumed a constant size for all cut trees for a given stand age.

###### 2.1.4. Constraints of Thinning Rates and the Number of Trees per Area

Because the parameters of SDMD were estimated using data from real forests subjected to standard silvicultural processes, unusual conditions such as extremely small number of trees per area should be avoided when using this model. Accordingly, the lower bound of was set to 200 trees/ha. In addition, it is recommended that a high thinning rate be avoided when using SDMD; however intensive thinning is often performed to reduce thinning costs. For this study, we placed emphasis on actual practice and restricted the thinning rates to the range 0–50%.

###### 2.1.5. Stem Profile

We used a stem profile curve to calculate the small end diameter of the logs in a complex yield model (defined in Section 2.1.8(2)). The Behre equation given below is a widely used relative stem profile curve:where is the height for which stem width is being calculated relative to total tree height (scaled from 0 to 1), is the diameter at relative height , and is the diameter at relative height 0.9. Parameters and are calculated as follows [52]:where is stem volume, is the total height of the tree, is the relative height at breast height, and is the DBH.

Equation (18) can be transformed as follows to provide the relationship between the diameter at any height and the DBH:

###### 2.1.6. Log Price

The relationship between the end diameter of 4 m long Japanese larch logs and the value in a real market in Nagano Prefecture [53] is shown in Table 1. Although there is a brisk market for 6–11 cm logs for log piles, the value of 12–14 cm logs is lower as there is less demand for this size. This price model adds multimodality to the NPV.