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

Moriguchi, Kai; Ueki, Tatsuhito; Saito, Masashi

doi:https://doi.org/10.1155/2015/173042

International Journal of Forestry Research

On this page

Abstract Introduction Materials and Methods Results Discussion Conclusion References Copyright Related Articles

Research Article | Open Access

Volume 2015 | Article ID 173042 | https://doi.org/10.1155/2015/173042

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

Kai Moriguchi,¹Tatsuhito Ueki,¹and Masashi Saito¹

Academic Editor: Ignacio García-González

Received18 Jun 2015

Revised27 Oct 2015

Accepted10 Nov 2015

Published30 Nov 2015

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³ 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⁶ 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³/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²/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³/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.

2.1.7. Cost of Regeneration and Yield

The thinning cost was set to 4588 yen/m³, based on the mean total cost of logging and transportation of Japanese larch [54], assuming a 50% subsidy is provided, and the cost of final cutting was set to 5987 yen/m³, assuming no subsidies. Usually no subsidy for final cutting is provided; however for the purpose of examining a range of patterns, we also tested the model on a scenario in which the final cost was similarly reduced to 2994 yen/m³. Hereafter, we refer to the two scenarios as F (full) and H (half), respectively. The difference between scenarios affects the flatness of the NPV. The present costs for regeneration are calculated as 1,078,590 yen/ha according to standard government sources [55], which includes ground clearance of a typical meadow at the beginning of the planting year, planting 2500 trees/ha, prevention of mammal damage by spreading Ziram solution in three years, bush cutting every year from one to five years, tree trimming and crosscutting in 10 and 15 years, 8% consumption tax, 0.9% annual interest, and 50% subsidy. This total was applied to in (1).

2.1.8. Yield Model

We defined two types of yield models that differ in the degree of difficulty involved in optimizing them.

(1) Simple Yield Model. This model calculates the yield volume by simply multiplying the yield rate by the stem volume of the harvested trees. We calculated income based on the price in Table 1 corresponding to the DBH of the harvested trees. The thinning and final cutting yield rates were set to 58% and 65%, respectively. Gaps in the classes in Table 1 were linearly interpolated. The income from a cut was calculated as follows:where is the income from a cutting per area (yen/ha), is the yield rate (%), is the total volume of cut trees per area (m³/ha), is the price of a cut tree per volume (yen/m³), and is the cutting cost per volume (yen/m³). Hereafter, we refer to this model as S and in combination with the half and full scenarios described in Section 2.1.7, S-H and S-F model are defined. The models are identical to the nonmonotonic pricing yield models used in a previous study [35], except for the self-thinning model and volume calculation used.

(2) Logging Yield Model. This model takes the logging process into consideration. As many 4 m logs (with 0.1 m margin) as possible were bucked, from stump height (0.5 m) to the top of the harvested tree. Diameters were calculated at every height from 4.6 m to the top of the tree in 4.1 m increments, using the stem profile curve, and were assumed to be end diameters. The Japanese Agricultural Standard for logs [56] was applied to the calculations of diameter and volume. The thinning and final cutting yield rates were set to 80% and 90%, respectively. Income from a cut was calculated as follows:where is the index number of a log from its stump, is the price of the th log (yen/m³), is the volume of the th log (m³) calculated aswhere is the end diameter of the th log (cm) (omitted for multiples of two for diameters > 14 cm and omitted for natural numbers otherwise), and is the number of bucked logs from a harvested tree, calculated as follows:where “floor” rounds down to the next lower integer.

Because the number of bucked logs from a tree and the end diameters were restricted to integers, the NPV was a discontinuous, multimodal function. Hereafter we refer to this model as L and in combination with the half and full scenarios described in Section 2.1.7, L-H and L-F model are defined.

2.2. Optimization Method

This section details the implementation of SA, its application to the thinning rate optimization problem, and the method we used to evaluate its performance.

2.2.1. Simulated Annealing

SA was developed based on the metal annealing process, which finds the global minimum “energy” required by gradually decreasing a parameter named “temperature.” At the beginning, the temperature is high, so frequent transitions occur for both low-energy and high-energy states. This makes it possible to find the global optimum by searching a wide solution space. As the temperature cools, the frequency of transitions to a higher-energy state decreases, and the system tends to transition to a lower-energy state more frequently. Using these features, approximate global minima can be found heuristically. The flowchart is shown in Figure 1. The procedure can be described as a Metropolis algorithm [57] with changes in temperature. SA is an adaptable method that requires the following four components: a cooling function, which controls the rate of decrease in temperature; a proposal density function, which is a probability distribution function (PDF) that generates candidate variables, that is, thinning rates; an energy function, which is the metaobjective function to be minimized; and an acceptance probability, a temperature-dependent PDF that decides whether or not the proposed state is accepted.

The cooling function is defined as an exponential decay function; that is, a decreasing rate determined by the initial and final temperatures is multiplied with temperature at time to obtain temperature at time . The proposal density is defined by a normal distribution with a mean of the present value of the target variable and a standard deviation that is optimized as described later. The energy function is defined as follows:where is the energy function, is the tentative minimum NPV, and is the tentative maximum NPV. Since our aim was to maximize the NPV, the problem was to minimize the negative NPV. This definition limits the differential of the energy function to the range 0-1. We set the acceptance probability as the Boltzmann distribution, as follows:where is the acceptance probability, is the energy of the candidate state, and is the energy of the present state. It is difficult to set appropriate temperature bounds to control the acceptance probability, since this depends on the units of the objective function when it is set directly as the energy function. The limits inherent in (25) facilitate control of the acceptance probability.

The following SA parameters are required: number of iterations per temperature; number of total iterations; initial temperature; final temperature; standard deviation of the proposal density. Because of the practical constraints of computing costs, we fixed parameters 1 and 2 to 5 × 10³ and 2 × 10⁵, respectively. We then obtained the best combinations of parameters 3–5 using CFE. Taking into consideration that log-transformed parameters are more suited for optimization and that final temperature must be lower than initial temperature, we chose the best combinations from 11³ (1331) patterns as shown below, using the mean NPV of 39 runs with random initial thinning as the objective function, to ensure average performance:where is initial temperature, is final temperature, is the standard deviation of the proposal density (cut rate, %), and , , and are the variables to be selected.

2.2.2. Control of Variables

The variables to be optimized are for each thinning age. Each of these is restricted to be lower than or equal to that at the previous thinning age. Controlling them directly requires changing the bounds according to the value of at the previous and next ages. However, this restriction can always be satisfied by constraining the thinning rate at each age to the range 0–50%. No thinning is simulated by setting the thinning rate to 0%. We implemented this method because it simplifies controlling .

2.2.3. Controlling for Infeasible Solutions

If high thinning rates are used at multiple cutting ages, at the final cutting may be <200 trees/ha. Penalty functions are often defined for infeasible solutions but may be difficult to define appropriately [42, 46]. The quality of definition has a significant influence on performance. Because we placed emphasis on the potential performance of SA and on fair comparison with RS (by random sampling only; see next section), the generation of infeasible solutions was not prevented; rather, candidate variables from the proposal density were sampled repeatedly until the generated solution was feasible. This differs from defining the energy function as infinity for infeasible solutions in that the latter method does not increase iterations if generated solutions were infeasible.

2.3. Evaluation

We derived feasible approximate global optima using CFE (Figure 2). The “combinatorial state” in this case is all patterns of repeated permutations of candidate thinning rates ranging from 0 to 50% in 5% increments. Since this solution was just an approximation based on a coarse thinning rate, it is possible to obtain better solutions (or at least solutions sufficiently close to the optimum) using other algorithms, given satisfactory performance.

We also applied RS to the same yield models, for the purposes of assessing the effectiveness of SA (Figure 3). Independent combinations of thinning rates with a uniform distribution over the range 0–50% were generated at each cutting age. Since this method was implemented as random sampling rather than a random walk, all solutions for each sample were independent of each other. A total of 10⁶ feasible solutions were simulated. We ran each model 39 times, and the optimal NPVs were sampled for comparison with those generated by SA.

Because the thinning rates are initiated randomly at the beginning of each SA process, the results can be tested statistically [44, 58, 59]. We tested our results to confirm independence between the initial and optimal solutions, differences between solutions generated by SA and RS, and differences between yield models.

3. Results

3.1. Examination of the Basis of Data

The best values of , , and were 1.2, 6.0, and 1.9, respectively, for the S-H model; 0.8, 1.5, and 1.6 for L-H; 0.0, 5.5, and 1.0 for S-F; and 0.0, 2.0, and 1.0 for L-F. There was no significant correlation between the initial and the optimal NPV for any model (Spearman’s correlation test; S-H: ; L-H: ; S-F: ; L-F: ). This indicates that the optimal solutions provided by SA can be regarded as independent of their initial solutions.

3.2. Comparing SA with CFE and RS

The NPV derived using CFE was 674,765 yen/ha for the S-H model, 769,711 yen/ha for L-H, 177,065 yen/ha for S-F, and 185,302 yen/ha for L-F. at final cutting was 2000, 819, 205, and 819 trees/ha for S-H, L-H, S-F, and L-F, respectively. The numerical differences between the optimal NPVs produced by SA and RS and those derived using CFE are presented in Table 2. The median and mean optimal NPVs provided by SA were better than those produced by the RS. In addition, some of the SA runs for each model provided better solutions than the CFE. Comparison of each model’s NPVs provided by SA with those by RS reveals distinct differences in the level of variability (Figure 4).

With respect to S-H, the standard deviation of the optimal NPVs was less than 1 yen/ha, much smaller than that for the other models (Table 2). SA also provided some low-variability outcomes for the S-F model (Figure 4); however, some of the runs might have converged to local optima (Figure 4) that caused the relatively large standard deviation for this model (Table 2). SA provided better solutions than RS for the L-H and L-F models, particularly with respect to the median, maximum, and mean optimal NPVs, although the standard deviations were larger (Table 2). Each of these values followed the Gumbel distribution (Kolmogorov-Smirnov test; L-H using SA: and using RS: ; L-F using SA: and using RS: ). A Gumbel distribution is defined by two parameters: location (equal to the mode) and scale. There are four possible SA/RS parameter combinations to describe the density of the data: common scale and location; independent scale and location; common scale and independent location; and independent scale and common location. A likelihood ratio test (using the test) between all patterns for each model revealed significant differences () between SA and RS for both the scale and location parameters. Using maximum likelihood estimation, the estimated L-H location parameters were −5128 for SA and −20,390 for RS. For the L-F model they were −13,663 for SA and −19,216 for RS. The estimated L-H scale parameters were 6362 for SA and 3352 for RS, and for the L-F model they were 8013 for SA and 3796 for RS.

3.3. Number of Updates

In a comparison of the models, the order of performance of the models in terms of the median or mean number of updates (Table 3) is the same as the order of performance in terms of the standard deviation of optimal NPVs (Table 2). The same pattern was apparent for the number of iterations at last update; however, the order of performance is reversed with respect to the medians of the L-H and L-F model results (Table 3). There were significant correlations between the optimal NPVs and the number of iterations at last update or the number of updates for each model except S-H (Spearman’s test; iterations at last update: S-H: ; S-F: ; L-H: ; L-F: ; number of updates: S-H: ; S-F: ; L-H: ; L-F: ).

Fluctuations in NPVs of the worst, best, and median runs of SA for each model are shown in Figure 5. Those of the S-H model included many updates and the trajectories were smooth. This indicates that for this model SA achieved its aim at each stage, that is, rough global sampling at the beginning and local optimization at the end. However, the updates of the S-F model were concentrated in the final stages. The L-F model runs contained the fewest updates over all stages.

3.4. Trajectory

Thinning rates at each age provided by the best, median, and worst runs of SA for each model are shown in Table 4. There were some differences between the thinning rates reached for each age by CFE and the best SA runs, since the thinning rates of CFE were split into 5% increments. The NPV of the best runs was higher than that of CFE, so the solutions produced by the best runs are closer to the global optima. However, there was some variation between the SA runs of each model except S-H. Even if the thinning rate is constant, the numbers of cut trees differ if numbers of trees per area are different. Figure 6 presents the trajectories for each model, which indicate the changes in the number of trees per area as the stand ages. Those of the worst, best, and median S-H model runs were very similar (Figure 6), reflecting almost identical thinning rates (Table 4). The other models showed some discrepancy between the trajectories produced by the worst runs and those produced by CFE or the best runs. With respect to the S-F model, trajectory of the worst runs differed from the others in 35 years, but converged again in 40 years (Figure 6). In contrast, all trajectories of the L-H model differed at almost every age. trajectory of the worst runs for the L-F model varied greatly from the others and that of median also differed from those of the best and CFE at most ages. However, the trajectory of the best run for each model matched that of the CFE closely.

4. Discussion

There were no significant correlations between the NPVs calculated with the initial solution and that calculated with the optimized solution for any of the models. This indicates that an SA process using optimized parameters could provide solutions, regardless of the initial solutions.

The best runs for all models provided better solutions than the respective CFE, whereas the RS never provided such solutions. The trajectory of the best runs of each yield model was almost identical to that of the respective CFE. These results suggest that SA can provide an approximate global optimum for similar yield models to those used in this study, if the parameters are optimized and the best solution is chosen from multiple runs. The S-H and S-F models are identical to the nonmonotonic pricing yield models used in a previous study [35] (see Introduction), except for the self-thinning model and volume calculation used. Our results therefore suggest that SA may be a better algorithm than MSPATH for optimizing thinning rates for single even-aged stands, for the purposes of obtaining reliable solutions.

Compared with RS, SA provided highly reliable solutions for the S-H model, as indicated by the negligible difference between the minimum and maximum NPVs provided by the SA runs (<1 yen/ha). With respect to the S-F model, a few of SA runs may have converged on local optima, although even those NPVs were higher than most of those produced by the RS runs. In contrast, the scale and location parameters of the Gumbel distribution of NPVs produced by SA for the L-H and L-F models were both significantly larger than those produced by the RS. The difference in location parameters means that NPVs produced by SA runs were higher overall than those produced by RS. However, the scale parameter determines the variance of the Gumbel distribution, so this indicates that the variation in NPVs produced by SA was larger than those produced by RS.

The only difference between the H and F model definitions was the final cutting cost, but this resulted in substantial differences in the variance of the optimal NPVs and the amounts of variation in the trajectories. It is possible that the high degree of sensitivity of the yield model caused large changes in the form of the objective function. In fact, in the S-F model, optimal at final cutting was 205 trees/ha, based on CFE, which is close to the lower bound of values, whereas in the S-H model it was 2000 trees/ha. This difference indicates that the change in final cutting cost resulted in large changes in the optimal solution. This suggests that the degree of variation in NPVs provided by SA was strongly influenced by the yield model.

Bettinger et al. [60] examined the performance of eight optimizing algorithms for spatial forest planning, using three benchmark models. They showed that SA, the parameters of which were decided based on “several trials,” could provide values close to the global optima or assumed “best” values for two of the models with a high probability, but not for the other model or the other algorithms. Our study shows that this is a common problem for the optimization of thinning rates for single stands; hence the variation in the solutions provided by SA depends on the model. In cases where the performance of SA is influenced by the parameters [45], this can be reduced by optimizing the parameters of SA themselves, as in the present study. However, our results suggest that this problem cannot be solved entirely by optimizing the parameters.

Considering these results, we suggest that using SA as demonstrated here has the potential to achieve the approximate global optima for thinning rates for single stands, if the controlling parameters are optimized. However, it is necessary to execute multiple runs and choose the best solutions from among them.

Despite the above, the reliability of SA is weakened by the large variation between the L-H and L-F models, and it is important to decrease this variation. The results of updating processes (in Section 3.3) indicate that few updates of L-H and L-F are the main cause of the large variation. Whereas the controlling parameters have been optimized, there is little room to improve them, and a simple increase in the number of iterations leads to increased computing time. Therefore other approaches may be advisable. In this study, we used orthodox SA because it is widely used. However, different approaches to SA implementation may be more efficient, for example, fast annealing [61] or adaptive annealing [62], or other types of SA such as the great deluge algorithm [63] used by Bettinger et al. [60]. At the same time, Borges et al. [46] showed that improving the neighborhood search method can enhance the performance of SA for spatial planning. Attempting to decrease the variation in optimal NPVs of SA and enhancing its reliability using all of these methods may be worthwhile.

5. Conclusion

With the aim of establishing a reliable method for optimizing thinning rates for single even-aged stands, we evaluated the performance of SA by comparing it with RS and CFE, using optimized parameters and four benchmark models. A comparison of 39 runs of the SA (2 × 10⁵ iterations of feasible solutions) and RS (10⁶ iterations) indicated that SA generally provided better solutions than RS. Further, the best SA run of each model provided better solutions than CFE did for that model. Therefore, if the controlling parameters for the SA are optimized, as in this study, the best of multiple runs has a high probability of approximating the global optimum. Because stand-level thinning optimization methods have rarely been compared with exact solutions such as those provided by CFE, their reliability has not previously been tested. At present, SA is the only method validated to provide solutions superior to those provided by CFE for the same problem.

However, the variance of the NPVs provided by SA differed noticeably between yield models: for the S-H model, the standard deviation of the optimal NPVs was <1 yen/ha, but for the L-H and L-F models it was substantially larger than that provided by RS. The dependence of the variance on the yield model decreases the reliability of SA and makes it necessary to select the best from among multiple runs. To increase the reliability of the SA method, the variation of optimal solutions within and between yield models should therefore be decreased.

Conflict of Interests

The authors declared that there is no conflict of interests.

Acknowledgment

This work was supported by JSPS (Japan Society for the Promotion of Science) KAKENHI Grant no. 15J04607.

References

Nagano Prefectural Government, Present Conditions of Forest in Nagano Prefecture at Fiscal 2014, 2014 (Japanese), http://www.pref.nagano.lg.jp/rinsei/sangyo/ringyo/toukei/minyurin/h26.html.
T. Arimizu, “Regulation of the cut by dynamic programming,” Journal of the Operations Research Society of Japan, vol. 1, no. 4, pp. 175–182, 1958.
View at: Google Scholar
A. Yoshimoto, V. G. L. Paredes, and J. D. Brodie, “Efficient optimization of an individual tree growth model. In the 1988 symposium on systems analysis in forest resources,” General Technical Report RM, vol. 161, pp. 154–162, 1988.
View at: Google Scholar
V. Paredes, L. Gonzalo, and J. D. Brodie, “Efficient specification and solution of the even-aged rotation and thinning problem,” Forest Science, vol. 33, no. 1, pp. 14–29, 1987.
View at: Google Scholar
R. G. Haight, J. D. Brodie, and W. G. Dahms, “A dynamic programming algorithm for optimization of lodgepole pine management,” Forest Science, vol. 31, no. 2, pp. 321–330, 1985.
View at: Google Scholar
A. Yoshimoto, R. G. Haight, and J. D. Brodie, “A comparison of the pattern search algorithm and the modified PATH algorithm for optimizing an individual tree model,” Forest Science, vol. 36, no. 2, pp. 394–412, 1990.
View at: Google Scholar
J. P. Roise, “A nonlinear programming approach to stand optimization,” Forest Science, vol. 32, no. 3, pp. 735–748, 1986.
View at: Google Scholar
R. Hooke and T. A. Jeeves, “‘Direct search’ solution of numerical and statistical problems,” Journal of the ACM, vol. 8, no. 2, pp. 212–229, 1961.
View at: Publisher Site | Google Scholar
S. H. Bullard, H. D. Sherali, and W. D. Klemperer, “Estimating optimal thinning and rotation for mixed-species timber stands using a random search algorithm,” Forest Science, vol. 31, no. 2, pp. 303–315, 1985.
View at: Google Scholar
P. Wikström and L. O. Eriksson, “Solving the stand management problem under biodiversity-related considerations,” Forest Ecology and Management, vol. 126, no. 3, pp. 361–376, 2000.
View at: Publisher Site | Google Scholar
M. J. D. Powell, “An efficient method for finding the minimum of a function of several variables without calculating derivatives,” The Computer Journal, vol. 7, no. 2, pp. 155–162, 1964.
View at: Publisher Site | Google Scholar | MathSciNet
A. Yoshimoto, “A dynamic programming model for forest stand management using MSPATH algorithm,” Proceedings of the Institute of Statistical Mathematics, vol. 51, no. 1, pp. 73–94, 2003.
View at: Google Scholar | MathSciNet
D. M. Adams and A. R. Ek, “Optimizing the management of uneven-aged forest stands,” Canadian Journal of Forest Research, vol. 4, no. 3, pp. 274–287, 1974.
View at: Publisher Site | Google Scholar
J. R. González-Olabarria, M. Palahí, T. Pukkala, and A. Trasobares, “Optimising the management of Pinus nigra Arn. stands under endogenous risk of fire in Catalonia,” Investigación Agraria: Sistemas y Recursos Forestales, vol. 17, no. 1, pp. 10–17, 2008.
View at: Publisher Site | Google Scholar
J. H. Gove and M. J. Ducey, “Optimal uneven-aged stocking guides: an application to spruce-fir stands in New England,” Forestry, vol. 87, no. 1, pp. 61–70, 2014.
View at: Publisher Site | Google Scholar
J. H. Gove and S. E. Fairweather, “Optimizing the management of uneven-aged forest stands: a stochastic approach,” Forest Science, vol. 38, no. 3, pp. 623–640, 1992.
View at: Google Scholar
T. Marutani, “The effect of site quality on economically optimal stand management,” Journal of Forest Economics, vol. 16, no. 1, pp. 35–46, 2010.
View at: Publisher Site | Google Scholar
J. Miina, “Optimizing thinning and rotation in a stand of Pinus sylvestris on a drained peatland site,” Scandinavian Journal of Forest Research, vol. 11, no. 1–4, pp. 182–192, 1996.
View at: Publisher Site | Google Scholar
M. Palahí and T. Pukkala, “Optimising the management of Scots pine (Pinus sylvestris L.) stands in Spain based on individual-tree models,” Annals of Forest Science, vol. 60, no. 2, pp. 105–114, 2003.
View at: Publisher Site | Google Scholar
T. Pukkala, J. Miina, M. Kurttila, and T. Kolström, “A spatial yield model for optimizing the thinning regime of mixed stands of Pinus sylvestris and Picea abies,” Scandinavian Journal of Forest Research, vol. 13, no. 1–4, pp. 31–42, 1998.
View at: Publisher Site | Google Scholar
O. Tahvonen, T. Pukkala, O. Laiho, E. Lähde, and S. Niinimäki, “Optimal management of uneven-aged Norway spruce stands,” Forest Ecology and Management, vol. 260, no. 1, pp. 106–115, 2010.
View at: Publisher Site | Google Scholar
G. J. Arthaud and W. D. Klemperer, “Optimizing high and low thinnings in loblolly pine with dynamic programming,” Canadian Journal of Forest Research, vol. 18, no. 9, pp. 1118–1122, 1988.
View at: Publisher Site | Google Scholar
D. H. Graetz, J. Sessions, and S. L. Garman, “Using stand-level optimization to reduce crown fire hazard,” Landscape and Urban Planning, vol. 80, no. 3, pp. 312–319, 2007.
View at: Publisher Site | Google Scholar
T. Heikkinen, “On optimal forest management: a dynamic programming approach,” International Journal of Society Systems Science, vol. 3, no. 3, pp. 217–235, 2011.
View at: Publisher Site | Google Scholar
A. Manabe, “Studies on the optimal thinning schedule for Todomatsu (Abies sachalinensis) plantations,” Bulletin of Foresty and Forest Products Research Institute, vol. 328, pp. 43–106, 1984.
View at: Google Scholar
P. Bettinger, D. Graetz, and J. Sessions, “A density-dependent stand-level optimization approach for deriving management prescriptions for interior northwest (USA) landscapes,” Forest Ecology and Management, vol. 217, no. 2-3, pp. 171–186, 2005.
View at: Publisher Site | Google Scholar
N. A. Ribeiro, P. Surový, and A. Yoshimoto, “Optimal regeneration regime under continuous crown cover requirements in cork oak woodlands,” FORMATH, vol. 11, pp. 83–102, 2012.
View at: Publisher Site | Google Scholar
K. Nakama, T. Ota, N. Mizoue, and S. Yoshida, “Effects of recovering logging residues on strategy and benefits of forest stand management,” Journal of the Japanese Forestry Society, vol. 93, no. 5, pp. 226–234, 2011.
View at: Google Scholar
T. Ota, S. Takahira, K. Nakama, S. Yoshida, and N. Mizoue, “Effectiveness of low-density planting in terms of planting cost reduction and felling income reduction,” Journal of the Japanese Forestry Society, vol. 95, no. 2, pp. 126–133, 2013.
View at: Google Scholar
S. Takahira, T. Murakami, N. Mizoue, S. Yoshida, and H. Kaga, “Evaluation of yield regulation with mathematical programming—comparison with the management by Kyushu Electric Power Co., Inc,” Kyushu Journal of Forest Research, vol. 59, pp. 154–157, 2006.
View at: Google Scholar
S. Takahira, T. Murakami, N. Mizoue, S. Yoshida, and H. Kaga, “Optimization of harvest scheduling to produce wood-based materials with dynamic programming: a case study in the forest of Kyushu Electric Power Co. Inc,” FORMATH, vol. 6, pp. 57–75, 2007.
View at: Google Scholar
A. Yoshimoto and R. Marušák, “Evaluation of carbon sequestration and thinning regimes within the optimization framework for forest stand management,” European Journal of Forest Research, vol. 126, no. 2, pp. 315–329, 2007.
View at: Publisher Site | Google Scholar
T. Ando, “Ecological studies on the stand density control in even-aged pure stand,” Bulletin of the Forestry and Forest Products Research Institute, vol. 210, pp. 1–153, 1968.
View at: Google Scholar
P. Bettinger, J. Sessions, and K. Boston, “A review of the status and use of validation procedures for heuristics used in forest planning,” Mathematical and Computational Forestry & Natural-Resource Sciences, vol. 1, no. 1, pp. 26–37, 2009.
View at: Google Scholar
K. Moriguchi, “Comparison of solutions of three operations research methods for forest stand management,” Journal of the Japanese Forest Society, vol. 95, no. 4, pp. 199–205, 2013.
View at: Google Scholar
J. R. González, T. Pukkala, and M. Palahí, “Optimising the management of Pinus sylvestris L. stand under risk of fire in Catalonia (north-east of Spain),” Annals of Forest Science, vol. 62, no. 6, pp. 493–501, 2005.
View at: Publisher Site | Google Scholar
V. Černý, “Thermodynamical approach to the traveling salesman problem: an efficient simulation algorithm,” Journal of Optimization Theory and Applications, vol. 45, no. 1, pp. 41–51, 1985.
View at: Publisher Site | Google Scholar | MathSciNet
S. Kirkpatrick, C. D. Gelatt Jr., and M. P. Vecchi, “Optimization by simulated annealing,” Science, vol. 220, no. 4598, pp. 671–680, 1983.
View at: Publisher Site | Google Scholar | MathSciNet
A. Drexl, “A simulated annealing approach to the multiconstraint zero-one knapsack problem,” Computing. Archives for Scientific Computing, vol. 40, no. 1, pp. 1–8, 1988.
View at: Publisher Site | Google Scholar | MathSciNet
K. Öhman and L. O. Eriksson, “Allowing for spatial consideration in long-term forest planning by linking linear programming with simulated annealing,” Forest Ecology and Management, vol. 161, no. 1–3, pp. 221–230, 2002.
View at: Publisher Site | Google Scholar
K. Ohman and L. O. Eriksson, “The core area concept in forming contiguous areas for long-term forest planning,” Canadian Journal of Forest Research, vol. 28, no. 7, pp. 1032–1039, 1998.
View at: Publisher Site | Google Scholar
C. Lockwood and T. Moore, “Harvest scheduling with spatial constraints: a simulated annealing approach,” Canadian Journal of Forest Research, vol. 23, no. 3, pp. 468–678, 1993.
View at: Publisher Site | Google Scholar
C. BoWang and K. V. Gadow, “Timber harvest planning with spatial objectives, using the method of simulated annealing,” Forstwissenschaftliches Centralblatt vereinigt mit Tharandter forstliches Jahrbuch, vol. 121, no. 1, pp. 25–34, 2002.
View at: Publisher Site | Google Scholar
L. Dong, P. Bettinger, Z. Liu, and H. Qin, “Spatial forest harvest scheduling for areas involving carbon and timber management goals,” Forests, vol. 6, no. 4, pp. 1362–1379, 2015.
View at: Publisher Site | Google Scholar
B. M. Strimbu and M. Paun, “Sensitivity of forest plan value to parameters of simulated annealing,” Canadian Journal of Forest Research, vol. 43, no. 1, pp. 28–38, 2013.
View at: Publisher Site | Google Scholar
P. Borges, T. Eid, and E. Bergseng, “Applying simulated annealing using different methods for the neighborhood search in forest planning problems,” European Journal of Operational Research, vol. 233, no. 3, pp. 700–710, 2014.
View at: Publisher Site | Google Scholar | MathSciNet
Japan Finance Corporation, The Interest Catalogue for Agriculture, Forestry and Fisheries at 18/06/2014, Japan Finance Corporation, 2014 (Japanese).
Japan Forestry Agency, Stand Density Management Diagram, Japan Forest Technology Association, 1981 (Japanese).
Y. Tadaki, “The pre-estimating of stem yield based on the competition-density effect,” Bulletin of the Forestry and Forest Products Research Institute, vol. 154, pp. 1–19, 1963.
View at: Google Scholar
Y. Ejiri, “Decisions on optimal thinning by a recurring decision-making model (I): a maximization of total gains of stem volume,” Journal of the Japanese Forestry Society, vol. 72, no. 4, pp. 304–315, 1990.
View at: Google Scholar
Forestry Division of Nagano Prefectural Government, Guide for Long Rotation Forest Management of Private Japanese Larch Stand in Nagano Prefecture, Forestry Division of Nagano Prefectural Government, 1991 (Japanese).
A. Inoue and Y. Kurokawa, “A new method for estimating relative stem profile equations: application to system yield tables,” Journal of the Japanese Forestry Society, vol. 83, no. 1, pp. 1–4, 2001.
View at: Google Scholar
Hokushin Log Market, The Log Prices at the 932nd Market, 2012 (Japanese).
Policy Planning Division of Japan Forestry Agency, Report of Cost for Log Production at Fiscal the 2011, Policy Planning Division of Japan Forestry Agency, 2013 (Japanese).
Nagano Prefectural Government, “Standard unit cost table for ‘Shinshu no mori zukuri’,” service at fiscal 2014 (for project at fiscal 2013) (Japanese).
View at: Google Scholar
Forestry and Fisheries Ministry of Agriculture, Japanese Agricultural Standard for Log, Forestry and Fisheries Ministry of Agriculture, 2007 (Japanese).
N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, “Equation of state calculations by fast computing machines,” The Journal of Chemical Physics, vol. 21, no. 6, pp. 1087–1092, 1953.
View at: Publisher Site | Google Scholar
P. Bettinger, J. Sessions, and K. N. Johnson, “Ensuring the compatibility of aquatic habitat and commodity production goals in eastern Oregon with a tabu search procedure,” Forest Science, vol. 44, no. 1, pp. 96–112, 1998.
View at: Google Scholar
B. L. Golden and F. B. Alt, “Interval estimation of a global optimum for large combinatorial problems,” Naval Research Logistics Quarterly, vol. 26, no. 1, pp. 69–77, 1979.
View at: Publisher Site | Google Scholar
P. Bettinger, D. Graetz, K. Boston, J. Sessions, and W. Chung, “Eight heuristic planning techniques applied to three increasingly difficult wildlife planning problems,” Silva Fennica, vol. 36, no. 2, pp. 561–584, 2002.
View at: Google Scholar
H. Szu and R. Hartley, “Fast simulated annealing,” Physics Letters A, vol. 122, no. 3-4, pp. 157–162, 1987.
View at: Publisher Site | Google Scholar
L. Ingber, “Very fast simulated re-annealing,” Mathematical and Computer Modelling, vol. 12, no. 8, pp. 967–973, 1989.
View at: Publisher Site | Google Scholar | MathSciNet
G. Dueck, “New optimization heuristics: the great deluge algorithm and the record-to-record travel,” Journal of Computational Physics, vol. 104, no. 1, pp. 86–92, 1993.
View at: Publisher Site | Google Scholar

Copyright

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.

PDF Download Citation

Download other formats

Order printed copies

Views

2537

Downloads

1075

Citations