Heuristic Optimization of Thinning Individual Douglas-Fir

Abstract

Research Highlights: (1) Optimizing mid-rotation thinning increased modeled land expectation values by as much as 5.1–10.1% over a representative reference prescription on plots planted at 2.7 and 3.7 m square spacings. (2) Eight heuristics, five of which were newly applied to selecting individual trees for thinning, produced thinning prescriptions of near identical quality. (3) Based on heuristic sampling properties, we introduced a variant of the hero heuristic with a 5.3–20% greater computational efficiency. Background and Objectives: Thinning, which is arguably the most subjective human intervention in the life of a stand, is commonly executed with limited decision support in tree selection. This study evaluated heuristics’ ability to support tree selection in a factorial experiment that considered the thinning method, tree density, thinning age, and rotation length. Materials and Methods: The Organon growth model was used for the financial optimization of even age Douglas-fir (Pseudotsuga menziesii (Mirb.) Franco) harvest rotations consisting of a single thinning followed by clearcutting on a high-productivity site. We evaluated two versions of the hero heuristic, four Monte Carlo heuristics (simulated annealing, record-to-record travel, threshold accepting, and great deluge), a genetic algorithm, and tabu search for their efficiency in maximizing land expectation value. Results: With 50–75 years rotations and a 4% discount rate, heuristic tree selection always increased land expectation values over other thinning methods. The two hero heuristics were the most computationally efficient methods. The four Monte Carlo heuristics required 2.8–3.4 times more computation than hero. The genetic algorithm and the tabu search required 4.2–8.4 and 21–52 times, respectively, more computation than hero. Conclusions: The accuracy of the resulting thinning prescriptions was limited by the quality of stand measurement, and the accuracy of the growth and yield models was linked to the heuristics rather than to the choice of heuristic. However, heuristic performance may be sensitive to the chosen models.

1. Introduction

While tree harvest scheduling has been studied since the late 1600s [1], it remains a topic of current interest due to its complexity. Given a planning unit containing thousands of trees, one form of forest scheduling seeks to determine which, if any, silvicultural interventions on the unit produce the most desirable outcomes. Most desirable has often been defined as producing the greatest amount of revenue from harvesting the trees [2]. Previous studies have shown revenue varies with the species composition of the trees, their current density and sizes, the future growth rates anticipated on the site, the present and anticipated pricing of the wood products made from harvested trees, the number and frequency of harvests, and how many trees of which types are removed in each harvest (e.g., [3,4,5,6,7,8,9,10]). This broad range of parameters has motivated the ongoing demand for techniques that can identify the most desirable management choices, ideally at the level of individual trees.

Forest planning abilities are closely associated with tree measurement capabilities, forest models predicting the effects of silvicultural actions, mathematical optimization techniques, and available computing power. When a unit is thinned by removing some of its trees, optimization seeks to identify the best set of trees to cut and the best harvest timing under a given set of planning objectives. These decisions also determine the trees retained for later harvest and other purposes. Recent works on thinning (Table 1) have therefore included laser scanning to identify individual trees, the projection of individual trees’ growth, the estimation of timber assortments provided by harvest, and the use of 16 different optimization methods.

This study expanded the range of heuristic optimization techniques that have been evaluated for developing thinning prescriptions based on selecting individual trees for harvest. While Pascual [3 and Vaukhonen et al. 8] selected individual trees and Xue et al. [7] used heuristics, Fransson et al. [5] conducted the only other study we are aware of which combines the two approaches. More commonly, thinning prescriptions are optimized by dividing trees into 2–10 size classes and determining the proportion of trees that should be cut within each size class [11]. While size class-based approaches may produce good prescriptions, their accuracy is restricted by the granularity of the used size classes. The use of size classes also reduces the amount of decision support provided because the optimization method typically does not identify specific trees to cut. Given a harvest unit containing N trees, harvesting staff must then identify the best implementation among the 2^N possible tree selections. This is a challenging task [12] but it can be supported by heuristic optimization of individual trees’ selections using a combination of onsite computing, mobile computing, and greater offsite computing power. Unfortunately, heuristics are not guaranteed to find the best possible solution and may instead converge to a good solution that is a local, rather than global, optimum.

Because it is desirable to optimize each unit on a managed landscape, the amount of time a heuristic requires to determine a good tree selection is important. For a given amount of computing power and planning time, a more computationally efficient heuristic can be run more times. These additional runs can be allocated to increasing the range of growth, harvest, and pricing combinations available to inform management decisions. Runs can also be used to replicate the optimization of scenarios of interest, increasing the probability of finding the global optimum and mitigating the risk of reaching only a relatively undesirable local optimum. While a range of techniques have been developed to increase the efficiency of subsequent heuristic runs by making of use of information from prior runs [13], this study was devoted to properties of initial heuristic runs.

A second concern is the characterization of optimal thinning prescriptions for Douglas-fir. Because tree selections depend upon species- and site-specific responses, prescriptions do not necessarily transfer between studies [4]. Douglas-fir is a major commercial species in the south-central portion of the Pacific Temperate Rainforest along the west coast of North America and a worldwide plantation species [14]. Within our study area of western Oregon and Washington, investigations of Douglas-fir thinning in the past 30 years have emphasized ecosystem services [15]. Comparatively little research attention has been given to linking developments in growth models, numerical optimization methods, and computing power in the support of selecting individual Douglas-fir trees.

Our objectives were therefore to (1) determine how the choice of heuristic affected tree selection quality, (2) quantify differences in computational efficiency among heuristics, (3) re-examine thinning prescriptions for Douglas-fir in the context of selecting individual trees, and (4) identify areas for further investigation when integrating heuristic tree selection into the planning and execution of thinning sales. The following sections summarize how heuristics were used to perform tree selection, present corresponding results for thinning three Douglas-fir densities, discuss the limitations of the methods used in this study, and conclude with recommendations for improved modeling.

2. Materials and Methods

2.1. The Financially Optimal Thinning Problem

A common objective when optimizing thinning prescriptions is to maximize the financial return from a timber harvest unit. Often, the maximization is done under the assumption of an infinite sequence of identical harvests, in which case the financial return is the land expectation value (LEV, also referred to as the bare land value or the soil expectation value) defined by the Faustmann formula (e.g., [4,5,7,9]). In the case of an infinite sequence of tree planting, growth followed by one thinning, continued growth, and then clearcutting at a rotation length of R years, Faustmann’s formula is:

(1)$LE{V}_{prescription,i}=\frac{{\left(1+r\right)}^R}{{\left(1+r\right)}^R-1}\left(NP{V}_{reforestation}+NP{V}_{thin,i}+NP{V}_{final,i}\right)+\frac{C_{annual}}{r}$

where r is the annual discount rate; C_annual is the annual management cost for holding the land; and NPV_{reforestation}, NPV_thin,i, and NPV_final,i are the net present values (NPVs) of the planting, thinning, and final harvest at the end of the rotation, respectively (Table 2). Here, the NPVs associated with thinning prescription i are calculated from the real mean appreciated prices $\overline{P}$(•) for the wood volume harvested, the volumes (V_j) of individual trees, and the timber sale administration, harvest, and haul costs C_variable and C_fixed.

(2)$NP{V}_{thin,i}=\frac{1}{{\left(1+r\right)}^T}\left(\left(\overline{P}\left(T\right){\left(1+a\right)}^T-C_{thin,variable}\right)\left({\displaystyle \sum }_{j=1}^N{u}_i\left(j\right)V_{thin,j}\right)-C_{thin,fixed}\right)$

(3)$NP{V}_{final,i}=\frac{1}{{\left(1+r\right)}^R}\left(\left(\overline{P}\left(R\right){\left(1+a\right)}^R-C_{final,variable}\right)\left({\displaystyle \sum }_{j=1}^N\neg u_i\left(j\right)V_{final,i,j}\right)-C_{final,fixed}\right)$

where T is the age of thinning, a is the real appreciation rate for Douglas-fir sawtimber, and N is again the number of trees representing the stand. The tree selection made by a heuristic at thinning is u_i, with u_i(j) taking the value 1 if tree j is thinned and 0 if it is retained until final harvest. At final harvest, the logical complement of u_i(j), ¬u_i(j), takes the value 0 if the tree was previously removed at thinning and is 1 if the tree is removed at final harvest. Unlike its volume when retained to final harvest, V_final,i,j, a tree’s volume at the time of thinning, V_thin,j, is independent of prescription i because the stand has not yet been silviculturally modified.

We obtained the individual tree volumes, V_thin,j and V_final,i,j, required by Equations (2) and (3) from the individual tree growth model Organon [16] and the Scribner board-foot volume regressions of Brackett [17]. Since this study emphasized Douglas-fir in western Oregon and Washington, we present results from Organon’s northwest Oregon variant as it applies to the middle of the study region (Figure 1). We initialized Organon from height and diameter at breast height measurements of individual trees on three Douglas-fir plots of differing densities (Figure 1 and Table 3; University of British Columbia [18] with additional information in Reukema [19]).

After initializing Organon, LEV maximization proceeded by simulating the plot’s growth to a specified thinning age, removing the trees indicated by u_i, and continuing simulation to the final harvest age. The resulting LEV was then calculated, the heuristic generated a new u_i, and the process repeated from the tree removal step until the heuristic ceased generating new u_i values to investigate.

2.2. Adaptation of Heuristics to Select Individual Trees

Because growth and yield calculations are computationally expensive, we sought heuristics that could quickly find a desirable choice of u. Since most heuristics incorporate randomness into their search process, a computationally efficient heuristic is one that produces high LEVs at low variance while requiring the smallest number of growth simulations. For this study, we elected to compare the efficiency of eight different heuristics (Table 4). Seven of these heuristics have been demonstrated in other forest scheduling problems (e.g. [22,23]). Simulated annealing and threshold accepting have been used for selecting individual trees with non-financial objectives [24,25] and Fransson et al. [5] used a genetic algorithm for the financial optimization of individual tree selection.

To our knowledge, this study was the first to investigate use of hero, the four Monte Carlo heuristics (simulated annealing, record-to-record travel, threshold accepting, and great deluge), and tabu search for financial optimization by selecting individual trees for thinning. As originally described, hero sequentially evaluates the elements of u. We introduce a variation of hero which randomizes the order of evaluation. For clarity, we refer to this new variant as hero stochastic and the original version of hero as hero sequential. Because of limited prior use and because we made adjustments for computational efficiency in individual tree selection, we describe implementation details and provide pseudocode for all eight heuristics in the Supplementary Material.

For the six heuristics with parameters, we sought combinations of parameters that resulted in good performance across all 69 available combinations of the three plots, four thinning ages, and six rotation lengths. Since LEV depends on these three variables, we made many heuristic parameters relative to a heuristic’s internal working values of LEV. Most other parameters were made relative to the number of trees on each plot. This approach allowed for a single parameterization of each heuristic to run across all combinations, avoiding the tuning of parameters for different factor levels. We identified parameter values by running a heuristic 100 or more times at each of several parameter levels, evaluating the resulting LEV distributions, and selecting the most effective combination of observed parameters. The resulting default parameters are provided in the Supplementary Material (Listings S1–S8). After selecting a heuristic’s parameterization, we evaluated its empirical LEV distribution by running it 1000 times on each combination of plot and harvest timing. The one exception was tabu search which, due to its computational cost, was run only 100 times per combination of thinning age and rotation length on the Nelder plot.

As controls on heuristic performance, we included a reference prescription representative of common silvicultural practices used in commercial thinning of Douglas-fir [32], as well as an optimized conventional prescription. The reference prescription primarily thinned from below (preferentially removing the smallest diameter trees) to a basal area retention of 27.5 m² ha⁻¹ and was intended to provide mid-rotation cashflow from units planted at spacings around 3.2 × 3.2 m square (1000 trees per hectare). It did not apply to the Nelder plot due to the plot’s variable planting density, and it is typically used in circumstances most similar to thinning the 2.7 and 3.7 m plots in this study by age 30. We therefore did not consider use of the reference prescription on the plots at later thinning ages. The reference prescription is also not intended for long intervals between thinning and final harvest, so we excluded its use with rotations beyond 60 years. While the prescription’s intent is to thin from below, other trees will likely be removed due to defects or to maintain spacing. We approximated defect or spacing constraints as a proportional removal (equal thinning of trees across the diameter range) of 10% of the trees in each plot. This created reference thinning intensities of 37% from below plus 10% proportional intensities, thus equaling 47% of trees on the 2.7 m plot and 44% plus 10% equaling 54% on the 3.7 m plot.

To produce optimized conventional prescriptions, we considered 1% steps in the intensity of thinning from below, proportional thinning, and thinning from above (preferential removal of the largest diameter trees) for a total of 8.61 million prescriptions. Each combination of plot, thinning age, and rotation length resulted in 124,807 prescriptions that removed 30–90% of trees from a plot. We selected the prescription with the highest LEV among these as the optimized conventional prescription.

3. Results

For all combinations of plot, thinning age, and rotation length, every run of all eight heuristics resulted in a thinning prescription with a higher LEV than available from conventional prescription optimization (Figure 2 and Figure 3). Conventional prescription optimization always resulted in higher LEVs than the reference prescription (Figure 2). Heuristic prescriptions increased LEV by 0.5–2.8% over optimized conventional prescriptions. In most cases, heuristics increased LEV by increasing total volume growth, but the relationship between volume and LEV lacked a clear pattern across plots, thinning ages, and rotation lengths (Figure 4). LEV increases on the 2.7 and 3.7 m plots were 5.1–10.1% relative to the reference prescription (Table 5).

Because the choice of heuristic had little effect on LEV (Figure 3), the primary criteria for selecting a heuristic is likely to be its computational cost (Figure 5). While hero did not guarantee an O(N² log N) runtime and exhibited median LEVs of 0.004–0.024% below the overall median, its N² log N behavior made hero one-to-two orders more computationally efficient than tabu search’s O(N³/log N). Tabu search’s O(N³) behavior resulted from evaluating N growth simulations in order to choose the next u_i, a computationally intensive process that increased the median LEV over by at most 0.003% and decreased the median LEV by 0.003% on the 3.7 m plot. The genetic algorithm’s expected runtimes were 2.5–12 times shorter than tabu search, with the greatest reduction occurring on the Nelder plot.

The four Monte Carlo heuristics all performed interchangeably (Figure 3, Figure 4 and Figure 5). During parameterization, we found their implementations converged to functional equivalence and that varying computational effort by ±20% had little effect on solution quality. However, the Monte Carlo heuristics required at least 1.5–2 times more computation to reach solutions of quality comparable to the hero variants (Figures S1–S3). Because the heuristics were configured to operate in an otherwise identical fashion during their initial LEV maximization phase (Supplementary Material, S1.3), hero’s greater efficiency was attributable to its more systematic assessment of each tree’s cut decision rather than performing random sampling with replacement as the Monte Carlo heuristics did.

Additionally, hero stochastic often converged to a solution more quickly than hero sequential, decreasing the expected runtime by 5.0–17% (Figure 5). Because the only difference between the two variants was hero’s stochastic iterative evaluation of all of a plot’s trees using sampling without replacement rather than in a fixed order, this increase in efficiency was attributable to hero’s stochastic randomization. Hero sequential lacked this feature, repeatedly evaluating trees in the order in which they appeared in the plot data (Supplementary Material, S1.2). The Monte Carlo heuristics and genetic algorithm also randomized against tree order, and tabu search’s steepest ascent approach rendered it independent of ordering as well.

The LEV increases of heuristic prescriptions over the reference prescription were due to increasing use of thinning from above, with later thinning ages and longer intervals between thinning and final harvest lengthened (Figure 6 and Figure 7). Fransson et al. reported optimal thinning ratios of 1.25–1.33 for selecting individual Norway spruce trees (Picea abies L.) with thinning ages of 40–48 years, a final harvest at 65 or 75 years, and plots with age 30 densities in the vicinity of 2000–2200 TPH (trees per hectare) ([5]; see Table 1 and Table 3). Matching Fransson et al.’s use of arithmetic means, we found thinning ratios of 1.20–1.50 for ages of 40 and 45 years and the same rotation lengths on the 3.7 m plot. Since the 3.7 m plot had an age 30 density of 596 TPH (Table 3), we interpreted these higher Douglas-fir thinning ratios at lower density as supporting the non-transferability of the results between species and sites asserted by Halbritter [4].

While Organon neglected the Nelder plot’s spatial density variation, it did consider a tree’s social position. The resulting prescriptions removed some of the small trees in dense areas of the plot and favored higher thinning ratios than on the 2.7 or 3.7 m plots, increasing removal of large trees in less dense areas of the plot (Figure 1, Figure 6, and Figure 7). This pattern was consistent with the prescriptions obtained for the more homogenous 2.7 and 3.7 m plots, and the Nelder plot’s comparatively smaller ΔLEV range (Figure 3) provided no evidence of heuristic optimization difficulty.

4. Discussion

We found the maximum LEV difference among heuristic runs to be less than 1%, an order of magnitude smaller than uncertainties in stand measurement [33,34], and LEV differences could be reduced below 0.1% by selecting the best result from multiple heuristic runs. It is also unlikely growth, yield, and financial models maintain a better-than-1% accuracy, both when preparing a thinning sale in the near term and when predicting a final harvest 20–40 years in the future. Therefore, the accuracy of a thinning prescription is most likely limited by the stand models used and not by choice of heuristic. However, realizing the LEV increase offered by heuristic tree selection is likely to be difficult due to the complexity of translating a prescription based on a few hundred trees, as was the case here, to a production harvest unit containing perhaps 10,000–40,000 trees.

Given sufficient computational resources, scanning a harvest unit with LiDAR (light detection and ranging) or SfM (structure from motion photogrammetry) and applying a spatially explicit growth model would allow for a heuristic to automate tree selection across the entire unit. Because scanning a harvest unit georeferences the trees within it and because a heuristic selects specific trees, thinning could then be performed by georeferenced harvest equipment [35], thus implementing the heuristic’s cut list. While Fransson et al. [5] implement elements of this approach in Norway spruce, a spatial model calibrated for Douglas-fir is needed. Extrapolation from the runtimes observed here also suggests the optimization of an entire harvest unit would require heuristic runs lasting several hours to a few days. Because scanning errors will likely be encountered during thinning and cutting errors may additionally occur, a need for an investigation into the rapid re-optimization of a thinning prescription is apparent, as is the desirability of performing initial optimizations more quickly.

Certain silvicultural concerns specific to our study area were not modeled in this study due to the use of an LEV formulation and Organon being a regression-based model without a carbon cycle. A common management objective is the restoration of large, structurally complex Douglas-firs and their eventual conversion to snags and down wood [15]. Because such trees continue to grow over multiple harvest rotations, the Faustmann assumption of an infinite series of identical rotations no longer holds. A constant climate is also implicit in the Faustmann assumption and, while Organon accepts site productivity in the form a site index, it is unclear whether Organon remains accurate if the site index changes under future climate scenarios. This suggests combining a physiological process-based model, the consideration of stand structural objectives over periods of 100–500 years or longer, and more flexible financial accounting [10] could provide both more desirable and more accurate thinning prescriptions.

Because heuristic performance can be problem-specific [36], changes in the used models may alter heuristics’ optimization success. This study’s finding of near equal LEVs across heuristics does not appear to be unique, as Xue et al. [7] obtained LEVs matched within 0.82% on five of seven Scots pine (Pinus sylvestris L.) stands across five heuristics and 10 diameter classes. However, we are not aware of any thinning optimization study that has evaluated a given heuristic across multiple growth, yield, or economic models. Our results suggest regression-based growth and yield models provide a smooth enough optimization surface the cost of carrying a population of solutions, as a genetic algorithm does, and that this not offset by a population’s greater resilience to local minima [37]. It follows single-solution, hill-climbing heuristics such as hero and Monte Carlo methods are more computationally efficient than the genetic algorithm in the current case. This property may not hold when models include weather, episodic disturbance, timber price fluctuations, or other effects creating more complexity in how trees are valued over time.

The thinning prescriptions obtained in this study may also be affected by four sources of bias or error. First, Organon relies on stand-level multipliers to increase diameter growth and reduce height growth after thinning, thus leading to overestimation of post-thinning growth in small diameter classes ([38] and personal communication [39]). Second, thinning may damage retained trees and increase their exposure. Greater retention may therefore be required to mitigate injury or windthrow risks on some sites. Third, in our study area, logs with scaling diameters near 55 cm or larger often incur higher costs and lower mill-delivered prices, which may favor thinning from above to remove trees before they exceed this size. Fourth, while the simple linear regressions used for $\overline{P}$(•) captured nearly all variation in log assortment prices (Table 2; R² > 0.997 and p < 0.001), they did so under the idealized assumptions of perfectly predictable pricing and tree taper. Providing heuristics with a refined growth model, information about disturbance and fluctuating price risks, or log values obtained from individual trees’ taper when a harvest unit is scanned would very likely alter the optimum tree selection.

However, despite these uncertainties, it appears use of proportional thinning or thinning from above could increase financial returns compared to thinning primarily from below. Emmingham et al. [32] evaluated thinning from above versus below in Douglas-fir experimental plots within our study area. Their 10-year results appear broadly consistent with the prescriptions obtained here but do suggest heuristic tree selection could be improved by considering crown ratios or other metrics of tree vigor available when stands are scanned.

5. Conclusions

We found selection of individual trees using heuristics consistently increased the modeled LEV over prescriptions employing combinations of thinning from above, proportional thinning, and thinning from below. While the LEV of tree selections varied between heuristic runs, differences among heuristic prescriptions were an order of magnitude smaller than typical volume prediction and valuation uncertainties. The eight heuristics also differed by up to two orders of magnitude in the computation required to find the prescriptions. The choice of heuristic was therefore of limited importance to the quality of thinning prescriptions but had large effects on the speed at which solutions could be obtained. We introduced a stochastic variant of hero that required less computation per run than hero’s sequential form or the other six heuristics considered.

It appears this study’s primary limitation was the growth, yield, and economic models used. Repetition of this study or conducting similar studies with a process-based spatial growth model, more detailed yield estimation, and a more flexible economic formulation is therefore desirable to verify the observed LEV increases. While the heuristic optimization methods employed in this study will likely improve in the future, changing the models used may also require changing the heuristics to maintain performance. Further work is also recommended to incorporate non-monetary values such as large tree recruitment and retention, snag creation, understory development, wildfire fuel loads and geometry, and resilience to other ecological risks. Over the past few decades, forest policy in our study area has sought to balance financial returns with other ecosystem services and we caution against applying this study’s optimization techniques without considering them within this broader context.

Supplementary Materials

The following are available online at https://www.mdpi.com/1999-4907/12/3/280/s1, Listings S1–S9: heuristic and enumeration pseudocode; Figure S1: heuristic convergence on the 2.7 m plot; Figure S2: heuristic convergence on the 3.7 m plot; Figure S3: heuristic convergence on the Nelder plot; Figure S4: LEV optimization of the 2.7 and 3.7 m plots combined; Section S3: soft-ware and regression imputation used. References [40,41,42,43,44,45,46,47,48] are cited in the supplementary materi-als.

Author Contributions

Conceptualization, B.M.S.; methodology, J.S.; software, T.W.; validation, T.W.; formal analysis, T.W.; investigation, T.W.; resources, B.M.S.; data curation, T.W.; writing—original draft preparation, T.W.; writing—review and editing, T.W., J.S. and B.M.S.; visualization, T.W.; supervision, B.M.S.; project administration, B.M.S.; funding acquisition, J.S. and B.M.S. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Oregon State University Graduate School and U.S. Department of Agriculture, grant number 2019-67019-29462.

Data Availability Statement

Data was obtained from [18] and is available from https://www.mkrf.forestry.ubc.ca/research/research-projects-database/ with the permission of the University of British Columbia.

Acknowledgments

The tree measurements used to initialize plot simulations were provided by the University of British Columbia’s Malcolm Knapp Research Forest and the British Columbia’s Ministry of Forests under BC Research Forests Project 2020-02. We thank Doug Maguire for his comments on Organon and suggestion of [32].

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. (a) Location of study area [20,21] and (b–d) stem maps of live trees on the 2.7 m, 3.7 m, and Nelder spacing trial plots, respectively, used for model initialization as of their measurement at age 30 years. Plots are shown to scale, but their placement relative to each other is not shown. Spacing in the Nelder plot varies from 0.91 to 4.94 m.

Figure 2. Land expectation values (LEVs) obtained from thinning the three Douglas-fir plots in this study as a function of the method used to create a thinning prescription, a plot’s stand age at the time of thinning, and length of the even age rotation. For heuristics, the median LEV was taken across all eight heuristics to integrate with Figure 3.

Figure 3. Empirical distributions of the differences (ΔLEV) between a plot’s median LEV and the LEVs obtained from 531,300 individual heuristic runs. Each heuristic was run 23,000 times on each plot, except for tabu search on the Nelder plot, which was limited to 2300 runs due to its computational cost (light grey shading); 98.8% of the runs converged within ±0.2% of the median LEV. The largest differences of −0.65%, −0.72%, and −0.69% were outliers on the 2.7 m, 3.7 m, and Nelder plots, respectively.

Figure 4. Changes in total merchantable volume harvested (thinned plus final harvest) and increases in LEV over optimized conventional prescriptions due to use of heuristic tree selection by plot, rotation length, and thinning age. Correlations between harvest volume and LEV were 0.08, 0.71, and 0.08 for the 2.7 m, 3.7 m, and Nelder plots, respectively.

Figure 5. Boxplots of the elapsed time required to run each of the eight heuristics on each plot using an Intel i7-3300 quad-core processor with a 3.3 GHz base clock. One core-second represents the use of one processor core for one second. Since all heuristics obtained similar LEVs, heuristics with shorter runtimes would typically be preferred.

Figure 6. Retention of trees on each plot, by diameter class, resulting from heuristic optimization across all combinations of thinning age and rotation length. While the reference prescription thinned from below, heuristic tree selection typically identified complex prescriptions thinning predominantly from above.

Figure 7. Thinning ratios obtained on each plot as a function of thinning method, a plot’s stand age at time of thinning, and the length of the even age rotation (the harvest interval was the number of years between thinning and end of rotation). Thinning ratios greater than 1.0 (dashed line) indicate thinning from above, and ratios below 1.0 indicate thinning from below.

Table 1

Comparison of this study to other recent investigations of thinning optimization. The authors of Çağlayan et al. and Yoshimoto et al. [2,11] non-exhaustively surveyed an additional 182 studies of optimal management from 2016 and earlier.

Study	Year	Tree Species	Objectives	Optimization Methods
this study	2020	Douglas-fir	LEV	Eight heuristics
Pascual [3]	2020	Stone pine	Value increment and spacing	Mixed integer programming
Halbritter [4]	2020	Any	LEV	Analytic
Fransson et al. [5]	2019	Norway spruce	LEV	Genetic algorithm (heuristic)
Jin et al. [6]	2019	Changbai larch	NPV	Particle swarm (heuristic)
Xue et al. [7]	2019	Scots pine	LEV	Five population heuristics
Vaukhonen and Pukkala [8]	2016	≥95% Norway spruce	Value increment	Iterative

LEV: land expectation value; NPV: net present value; value increment: tree growth in the first model time step after thinning. Tree species are Douglas-fir (Pseudotsuga menziesii (Mirb.) Franco), stone pine (Pinus pinea L.), Norway spruce (Picea abies L.), Changbai larch (Larix olgensis A. Henry) and Scots pine (Pinus sylvestris L.).

Table 2

Economic parameters used in calculating per hectare land expectation value (LEV). Timber volumes in thousands of board feet (MBF) are measured using the Scribner Decimal C long log scale.

Parameter	Value
Reforestation cost, NPVreforestation	US$ 813 ha−1 + US$0.50 seedling ha−1
Discount rate, r	4% year−1
Thinning age, T	30–45 years in 5-year increments
Rotation length, R	50–75 years in 5-year increments
Selection vector of individual trees, ui	controlled by heuristic
Volume of Douglas-fir j at thinning, Vthin,j	MBF ha−1 from growth and yield models
Mean Douglas-fir pond value at thinning, $\overline{P}$(T)	US$453.67 + 2.776T MBF−1, 7.3 m logs
Douglas-fir real appreciation rate, a	1% year−1 within rotation
Fixed thinning cost, Cthin,fixed	US$148 ha−1
Variable thinning cost, Cthin,variable	US$275 MBF−1
Volume of Douglas-fir j at final harvest, Vfinal,i,j	MBF ha−1 from growth and yield models
Mean Douglas-fir pond value at final harvest, $\overline{P}$(R)	US$567.56 + 0.365R MBF−1, 12.2 m logs
Fixed final harvest cost, Cfinal,fixed	US$247 ha−1
Variable final harvest cost, Cfinal,variable	US$250 MBF−1
Annual management cost, Cannual	US$18.50 ha−1

Table 3

Summary of the three Douglas-fir spacing trial plots used for growth model initialization, as measured at age 30. All three plots are located in the Malcolm Knapp Research Forest in southwestern British Columbia, Canada (Figure 1) and have 50 years Douglas-fir site indices of 39.5 ± 0.5 m.

PlotSpacing	Area(ha)	LiveTrees	Mean Density(TPH)	QMD (cm)	H100 (m)	Mean Planting Density (TPH)	Tree HeightsImputed
2.7 m square	0.25	222	996	29.5	27.0	1278	0
3.7 m square	0.22	147	596	23.3	28.3	730	0
Nelder radial	0.75	571	757	24.1	28.6	1035	4

TPH: trees per hectare; QMD: quadratic mean diameter; H₁₀₀: average height of the tallest 100 trees ha⁻¹ measured. The regression for imputing Nelder plot heights is provided in Section S3 of the Supplementary Materials.

Table 4

Comparison of heuristics evaluated in this study in order of increasing computational cost. Details of heuristic implementations and optimization mechanisms are provided in the Supplementary Material.

Heuristic	Cost IncludingGrowth Model	Optimization Mechanisms	Parameters
hero stochastic [this study]	~O(N2 log N)	sampling without replacement	0
hero sequential [26]		sequential iteration	0
Simulated annealing [27]	~O(19N2)	sampling with replacement, reheating	6
Record-to-record travel [28]			4
Threshold accepting [29]			1 or 4+
Great deluge [28]			5
Steady state genetic algorithm [30]	~O(60N2.6)	parent selection, uniform crossover, mutation, replacement	4
Tabu search [31]	~O(N3/log N)	steepest ascent, tenure	1–2

O(•): empirically observed cost expressed in Landau notation; N: number of trees used to describe stand.

Table 5

Changes in merchantable volume harvested and increases of median LEVs relative to the reference prescription from use of heuristic tree selection on the 2.7 and 3.7 m plots. The reference prescription was not valid for the Nelder plot, other rotation lengths, or other thinning ages.

Rotation Length(Years)	2.7 m Plot Increase in		3.7 m Plot Increase in
	Volume (%)	LEV (%)	Volume (%)	LEV (%)
50	2.37	5.05	3.26	4.55
55	1.02	7.00	2.66	5.27
60	0.59	10.1	2.20	6.42

Minimum and maximum values across both plots are indicated in bold.