In: Bevers, Michael; Barrett, Tara M., comps. 2005. Systems Analysis in Forest Resources: Proceedings of the 2003 Symposium. General Technical Report 

PNW-GTR-656. Portland, OR: U.S. Department of Agriculture, Forest Service, Pacific Northwest Research Station.

1 Y.-H. Kim is a Graduate Research Assistant, Department of Forest Resources, Oregon State University, Corvallis, OR 97331.2 P. Bettinger is an Associate 

Professor, Warnell School of Forest Resources, University of Georgia, Athens, GA 30602.

Systems Analysis in Forest Resources: Proceedings of the 2003 Symposium

SPATIAL OPTIMIZATION OF FUEL

MANAGEMENT ACTIVITIES

Young-Hwan Kim1 and Pete Bettinger2

INTRODUCTION

The use of fuel management treatments in western

forests can reduce undesirable consequences of wildfires

(cost, sizes, ecological damage, threats to developed areas)

only if wildfire behavior is modified at the landscape-scale

(~103-106 acres). Individual management activities can be

expected to modify fire behavior on a local scale, but alone,

may have negligible impact on the overall growth and

behavior of larger fires. The cumulative effects of individual fuel management activities on fire behavior may depend

heavily upon their spatial pattern, as the behavior of fires

may only be sensitive to management activities when they

are scaled and arranged to disrupt the progress of the fire.

To this point in time, research on landscape effects of

fuel management activities has been mostly theoretical.

Observations of fuel patterns in the forests of California

(van Wagtendonk 1995, Parsons and van Wagtendonk 1996)

ABSTRACT

This paper presents preliminary work associated with scheduling fuel management treatments to achieve timber harvest

and landscape pattern goals. Four landscape patterns of management activities are modeled (dispersed, clumped, random,

and regular). The overall intent of this research is to examine the effects of spatial and temporal placement of fuel management activities on resulting wildfire behavior although we do not provide results on wildfire behavior here. However, in this

paper, we describe the forest planning scheduling processes that provide solutions that both place activities in a pattern on

the landscape, and achieve commodity production goals. Results indicate that the scheduling methodology developed with

this preliminary research is adequate in allowing one to optimize the spatial pattern of management units across a landscape

while attempting to also achieve a commodity production goal. The scheduling methodologies developed in this study are

expected to facilitate an understanding of the impact of the spatial variation of management activities on wildfire behavior

and management goals. 

205

support the idea that spatial fragmentation of fuels can affect

fire size and behavior. Isolated management activities, however, may have no effect on the growth and progress of

large fires (Dunn 1989). Several basic spatial patterns of

management activities can be distinguished based on their

overlap, or the ability of nearby treated units to disrupt fire

growth. Random patterns of fuels (Finney 2003) have no

requirement for overlap and therefore reduce spread rate in

a sigmoidal fashion, meaning that relatively large proportions of the landscape must be treated to substantially

reduce fire sizes. Parallel strips (Fujioka 1985, Martin

1988, Catchpole and others 1989) completely overlap in

one direction, and are the most efficient at reducing fire

spread rates, as only small fractions of the landscape need

to be treated. However, this strategy unrealistically requires

fires to always move perpendicular to the strips. Regular

patterns of dispersed treatments (Finney 2001) partially

overlap and can reduce spread rate, and are more flexible

in accommodating spatial constraints because they are not

necessarily spatially connected.
Several landscape simulation approaches are currently

used for spatially modeling fire and forest development

(i.e., Keane and others 1997, Jones and Chew 1999,

Mladenoff and He 1999). Some of these have been proposed for modeling effects of treatments and for optimizing

the scheduling of fuel treatment. None, however, account

for the topological effects of fuel management activities

with respect to landscape fire behavior. This requirement

suggests the use of integer or heuristic scheduling models

(rather than linear programming) due to the use of spatial

information, and the fine-scale conditions of landscape

units. 

The overall objective of this study is to understand how

the spatial pattern of fuel management activities influences

characteristics of fire behavior. In this paper, however, we

describe preliminary work aimed at the development of

methodologies for arranging fuel management activities 

in desired patterns across a landscape while meeting other

management goals. These scheduling processes utilize

heuristic techniques, and are expected to facilitate a variety

of studies that center on the spatial variation of management

activities, and the associated effects on management goals.

STUDY SITE AND DATA PREPARATION

In a recent study (Hirte 2002), vegetation changes and

timber harvest under several management prescriptions had

been simulated in private lands (approximately 11,000

acres) located in the Upper Grand Ronde River basin in

northeastern Oregon (fig. 1). Most of this area is surrounded

by U.S. Forest Service land (Wallowa-Whitman National

Forest). GIS databases and forest structure data developed

for the previous study were used for scheduling fuel management activities. When scheduling activities across the

landscape, we use centroids of management units to represent their location. The x, y coordinates of the centroids

were generated and saved using ArcView software and its

extensions. In addition, scheduling of fuel management

activities requires some attribute data accompanied with

GIS databases describing specific vegetation structure of

each management unit. All required attribute data (forest

structure data) were exported in the ASCII format, and made

suitable for the various scheduling procedures (described

below). The yield data represents changes in stand condition under 5 management prescriptions during 5 ten-year

management periods (50 years). These prescriptions were

not optimized for economic nor ecologic goals. They were

developed to represent private forest management activities

consistent with a survey of forest landowners conducted by

the Oregon Department of Forestry in the late 1990's. Hirte

(2002) described the results of the survey, and used the same

prescriptions we use, yet on a smaller portion of our landscape, and for optimizing different goals. In sum, the survey

of private landowners was simply used to help us determine

the types of prescriptions that would be reasonable for 

this area. Readers are encouraged to contact the Oregon

Department of Forestry (Salem, OR) for more detail on 

the survey of private landowners.

SCHEDULING OF SPATIAL PATTERNS 

OF FUEL MANAGEMENT ACTIVITIES

Four spatial patterns of fuel management activities were

examined in this research, which included three basic landscape patterns (dispersed pattern, clustered pattern, and

random pattern) and an artificial pattern (regular pattern).

These spatial patterns of fuel management activities were

scheduled with two heuristic modeling techniques: the Great

Deluge Algorithm (GDA) and Tabu Search (TS). GDA was

introduced by Dueck (1993) and applied to forest planning

problems in Bettinger and others (2002). Tabu Search was

initially developed by Glover (1989) and applied to forest

planning problems in Bettinger and others (1997, 1998, 

and 2002). 

Scheduling procedures were repeated 30 times for each

spatial pattern, to find the best solution that spatially optimizes the pattern across a landscape. Each repetition started with a random schedule of management activities to

make the 30 resulting solutions independent. For quantifying the effects of solutions more accurately, a control solution with no management activities scheduled was generated.

206

Figure 1—Study site: Upper Grand Ronde river basin in eastern Oregon. 
Dispersed Pattern of Fuel Management Activities

A dispersed pattern is a pattern in which management

units are widely spread across the landscape, minimizing

clustering. Here, ideal dispersed patterns are assumed to

maximize total distance between management units, and

minimize deviations between actual harvest volume and 

a harvest volume target. In order to generate a pattern as

close to the ideal pattern, the following objective function

was developed. 

Minimize

P Nk P Nk-1 Nk

WH ∑ ∑Hik –T – WD ∑ ∑ ∑ Dij [1]

k=1 i=1 k=1 i=1 j=i+1

Where: 

WH : Weight corresponding to the even harvest to the 

target

WD : Weight corresponding to the dispersion 

(WH + WD = 1)

Hik : Harvest volume from unit i in time period k (i = 1, 2,

..., Nk, k = 1, 2, ..., P)

T : Target volume of timber harvesting (T = 1,000

MBF (thousand board feet))

Dij : Distance between centroids of unit i and j (i = 1, 2, 

..., Nk-1, j = 2, 3, ..., Nk)

i, j : Index of management units scheduled for harvest

k : Index of management periods 

P : Total number of time periods (P = 5)

Nk : The set of management units scheduled for harvest 

in time period k

A scheduling procedure based on the above function

seeks a solution that minimizes the difference between

actual harvest volume and a harvest volume target, and

maximizes the total distance between centroids of management units. Original GDA seeks a solution with a higher

peak (higher objective function value) as increasing waterlevel (threshold value), and finally produces an optimized

solution which is expected to have highest peak (maximum

objective function value). Since the optimized solution in

this research was expected to have the minimum objective

function value, the algorithm was modified to seek a solution with a lower bottom as discharging water from a lake

(fig. 2). In addition, a spatial constraint was applied for

limiting the size of clearcuts to a maximum of 120 acres.

This constraint is based on the unit restriction model introduced by Murray (1999):

207

Figure 2—Flowchart of scheduling processes

for dispersed, clumped, and random landscape

pattern.

 

   

 








xi,k + xj,k ≤ 1 i,k,j ∈ Si [2]

xi,k , xj,k : 1, if unit i or j is scheduled for harvest in

period k. Otherwise, 0.

Si : set of harvest units adjacent to unit i

Thus, if a solution violated this spatial constraint, it

would be rejected no matter what objective function value

resulted. 3 stopping criteria were used in the modified version of GDA: total iterations, non-improved iterations, and

water-level. Parameters related to these stopping criteria

were set as follows:

• Maximum total iterations: 5,000,000

• Maximum non-improved iterations: 1,000,000

• Initial water-level: 100,000

• Discharging speed: 1

• Minimum water-level: -100,000

Objective function values might obviously vary according to the weights assumed, so 9 weight combinations (0.9,

0.8, 0.7, ..., and 0.1) were tested to determine the most

appropriate weights for both patterning and even flow

objective. The entire process described in the figure 2 was

repeated 10 times for each weight combination with the

above parameters. From these test trials, a weight values

(WH = 0.4 and WD = 0.6) were chosen for further processing. With the selected weight combination and the same

parameters, the entire process was repeated 20 more times.

The weights we have assumed are arbitrary, and obviously

one can envision larger sets of weights using more detailed

categories, but the added value of assessing this level of

detail in weighting the objectives is speculative. The choice

of weights was made by evaluating the point where dramatic differences in the objective values occurred (i.e., 

the threshold where a change in weights caused dramatic

declines in the objective function value). Therefore, the

choice of weights made here and with the other scheduling

techniques described in this paper was made using our

judgement in an assessment of trial runs of the scheduling

model (i.e., by determining at what point along the range

of weights better solutions are found).

Clumped Pattern

A clumped pattern is assumed to be a pattern in which

management units are clustered on the landscape. Here, the

ideal clumped pattern is assumed to minimize the total distance between management units and minimize the deviation between actual harvest volume and a harvest volume

target. While the dispersed pattern is expected to maximize

total distance between management units, the clumped pattern is expected to minimize it. Therefore, equation 1 was

modified to accept this distinction by adding the two portions of the objective function as follows:

Minimize

P Nk P Nk-1 Nk

WH ∑ ∑Hik –T + WD ∑ ∑ ∑ Dij [3]

k=1 i=1 k=1 i=1 j=i+1

The scheduling procedure seeks a solution that minimizes the difference between actual harvest volume and

harvest volume target and also minimizes the total distance

between centroids of management units. In addition, the

spatial constraint limiting the size of clearcuts to a maximum of 120 acres was applied. The scheduling process for

the simulation of clumped pattern was the same as that of

dispersed pattern (GDA). However, some of the parameters

related to the stopping criteria – initial water level and minimum water level – were altered based on trial runs of the

GDA model. Listed below are the parameters used in the

GDA scheduling process.

• Maximum total iterations: 5,000,000

• Maximum non-improved iterations: 1,000,000

• Initial water-level: 500,000

• Discharging speed: 1

• Minimum water-level: 0

Nine weight combinations were also tested and the most

appropriate weight values (WH = 0.1 and WD = 0.9) were

chosen from the test trials. 

Random Pattern

A random pattern is a pattern in which management

units are randomly allocated across landscape. Within the

GDA, management units are randomly chosen and random

prescriptions are assigned. Solutions generated with this

process are assumed to be random across the landscape

(although the pattern may be influenced by the actual pattern of vegetation in the study area). Therefore, dispersion

of management units is not a factor in scheduling a random

pattern. Deviation between actual harvest volume and the

harvest volume target is the only criterion for evaluating

the acceptability of a solution. Thus, the latter portion of

equations 1 and 2, corresponding to the dispersion of management units, was not included in the objective function

for the random pattern. Therefore, a solution that provides

the smallest deviation between actual harvest volume and

the harvest volume target is generated from the scheduling

process.

Minimize

P Nk

∑ ∑Hik –T [4]

k=1 i=1 

208

A

 

   

 

 

   

 








Most of the GDA parameters are the same as those used

in the development of dispersed and clumped patterns, but

a few have been altered based on the trial runs of the GDA

model. Although composed of randomly selected units,

solutions violating the spatial constraint for a large clearcut

patch were rejected.

• Maximum total iterations: 5,000,000

• Maximum non-improved iterations: 1,000,000

• Initial water-level: 100,000

• Discharging speed: 1

• Minimum water-level: 0

Regular Pattern

A regular pattern is a pattern in which management

units are systematically allocated across landscape using a

constant spatial interval. Ideally, units scheduled for harvest in the regular pattern are expected to have the same

distance to four neighbor units (northern, southern, eastern

& western). The “interval”, therefore, could be defined as

desired distance between centroids of management units

that produce an ideal regular pattern. To enable one to generate a regular pattern, a different approach was developed

and utilized for dispersing management units. It is based on

the following idea:

• Select one initial management unit

• Acquire the x, y coordinate of centroid of the unit

• Generate “systematic points” by adding or subtracting

a given interval to x, y coordinate of the centroid 

(fig. 3)

• Calculate distance between each systematic point and

centroids of all units

• Find the nearest centroid for each systematic point

• Check whether each systematic point is located within

the boundary of the study site

• Exclude systematic points located outside of the study

site

• Save the management unit nearest to each systematic

point

• Generate a solution by assigning a feasible prescription to the saved units

• Calculate the objective function value and evaluate the

solution

One of the issues related to the above idea is how to

exclude systematic points located outside of the study site.

In order to automate this process, inspecting whether a systematic point is out of the study site was possible by testing

whether a vector connecting a systematic point and its

nearest unit centroid is intersected by any boundary vector

surrounding the study site. That is to say, if a systematic

point were located outside of the boundary, the vector connecting the systematic point and its nearest unit centroid

should be intersected by at least one boundary vector (fig.

4). To inspect whether two vectors intersect, a two-step

process introduced by Loudon (1999) was used in the program: a quick rejection test and a straddle test. If both tests

succeed, two vectors intersect and thereby, the systematic

point is out of the study site.

The quick rejection test is initiated by constructing a

rectangle called a bounding box that surrounds each vector.

A vector between a systematic point and its nearest unit

centroid has two end nodes, n1 = (x1, y1) and n2 = (x2, y2).

The bounding box of the vector is a rectangle with lower

left point (min(x1, x2) (min (y1, y2)) and upper right point

(max(x1, x2), max(y1, y2)). Also, a boundary vector has two

end nodes, n3 = (x3, y3) and n4 = (x4, y4) , and a bounding

box with lower left point (min(x3, x4), min(y3, y4)) and

upper right point (max(x3, x4), max(y3, y4)). If bounding

boxes of the two vectors intersect, all of the following tests

must be true (fig. 5):

max(x1, x2) ≥ min(x3, x4) max(x3, x4) ≥ min(x1, x2)

max(y1, y2) ≥ min(y3, y4) max(y3, y4) ≥ min(y1, y2)

A straddle test follows only when the quick rejection

test succeeds. To examine whether a vector straddles another,

209

Figure 3—Systematic points generated to facilitate modeling the regular

landscape pattern.

Figure 4—Vectors between systematic points and their nearest unit c

entroids.
the orientation of n3 relative to n2 is compared with that of

n4 relative to n2. Orientation of n3 and n4 convey whether

the nodes are clockwise or counterclockwise from n2 with

respect to n1. The orientation of n3 and n4 are determined

by following equations:

z1 = (x3 – x1)(y2 – y1) – (x2 – x1)(y3 – y1)

z2 = (x4 – x1)(y2 – y1) – (x2 – x1)(y4 – y1)

If the sign of z1 and z2 are different, or either one is 0,

the vectors straddle each other, and the two vectors intersect. Figure 5 describes results of the quick rejection test

and the straddle test based on 4 different cases. 

Unlike the scheduling of dispersed, clumped, or random

patterns, Tabu Search (TS), instead of GDA, was used in

the scheduling of the regular pattern (fig. 6). TS generally

consumes much processing time for its neighborhood search.

However, limited numbers of neighbors would have been

searched in this case because prescriptions were assigned

only to the selected units. Therefore, TS might be more

efficient in scheduling than GDA when developing a regular

pattern. Since management units have already been chosen

before prescriptions are assigned to them, assigning a prescription to management units has no influence on the dispersion of management units in the scheduling of regular

pattern. As the result, dispersion of management units was

not counted in the objective function. In addition, according to the current prescriptions, a unit scheduled for harvest

in the first time period might be scheduled to be harvested

again in one of the following time periods. This means a

set of prescriptions assigned to management units for one

time period could affect scheduling of other following time

period. By this reason, a solution that guarantees the regular pattern over the entire planning horizon is rarely

obtained. Therefore, the scheduling process seeks a solution that optimizes the harvest in the first time period.

N

∑Hi1 –T [5]

i=1 

Two limitations of this scheduling process are: 1) even

though units have been identified to be contained within the

regular pattern, a treatment is not guaranteed to be assigned

to each, and 2) only those units and prescriptions available

are used to attempt to achieve even-flow harvest.

210

Figure 5—Examples of the quick rejection test and

the straddle test for use in the generation of the

regular landscape pattern.



  

 
Since limited amount of information is available to

specify the most efficient interval between management

units for reducing the fire damage, a 3.0 mile interval was

assumed, although other intervals could also be used. The

entire scheduling process described in the figure 6 was

repeated 30 times. In each run of the process, 50 management units were randomly selected and examined as the

initial unit for generating systematic points. Parameters

related to the TS scheduling process include:

• Maximum non-improved iterations: 1,000

• Size of tabu list: 10

RESULTS

The scheduling methodology developed in this research

acceptably allocated fuel management activities in desired

spatial patterns while meeting a commodity production goal.

Two criteria were used to assess how well the scheduling

process worked: a visual assessment of the pattern of

211

Figure 6—Flowchart of scheduling process for the regular landscape pattern.
212

Figure 7—Spatial patterns of management units developed by the scheduling processes.

Table 1—Harvest volume for the best solutions generated by the scheduling models.

Harvest volume (MBF)

Spatial pattern Period 1 Period 2 Period 3 Period 4 Period 5 Total 

Dispersed 1,001,607 1,001,827 1,003,161 1,007,219 1,001,348 5,015,172

Clumped 997,463 995,970 999,370 997,787 995,731 4,986,321

Random 1,000,091 999,675 1,000,175 999,948 1,000,177 5,000,066

Regular 994,765 — — — — —

Table 2—Total number of and area of management units scheduled 

for harvest.

Area (Acres)

Spatial pattern # of Units Total Thinning Clearcut 

Dispersed 80 353 34 (9.6%) 319 (90.4%)

Clumped 34 209 7 (3.3%) 202 (96.7%)

Random 61 262 21 (8.0%) 241 (92.0%)

Regular 34 246 33 (13.4%) 213 (86.6%)
scheduled activities (i.e., do the scheduled activities look 

as if they have been scheduled across the landscape in the

pattern desired?), and an assessment of the objective (i.e.,

was an even-flow of timber harvest volume achieved, and

was it high?). While no statistical tests were utilized to test

the hypotheses that these goals were met, we feel confident

that the model is adequate for optimizing the spatial pattern

of management activities. More reasonable criteria will be

necessary for evaluating the patterns produced, such as

using FRAGSTATS, and one of the landscape metrics provided, to compare the patterns of various schedules. We

leave this for further study. 

Figure 7 represents management units which were

scheduled for harvest by best solutions retained from the

scheduling process. Figure 7 clearly verifies the distinction

between the spatial patterns. Also, as shown in the table 1,

the best solutions produced an acceptable even-flow harvest level, as compared to the target volume (1,000 MBF).

Harvest volumes of each spatial pattern are quite close to

the target volume during each time period. These results

enable us to draw a conclusion that the methodology developed in this research was adequate to optimize the spatial

pattern of management units across a real landscape while

seeking to achieve a commodity production goal.

Total area of management units scheduled for thinning

or clearcut varied by the pattern modeled (table 2). The 

dispersed pattern and clumped pattern were particularly

distinctive in the total harvest area. For example, much

more area was scheduled for harvest in the dispersed pattern compared to the others. Also, the area of thinning in

the clumped pattern is much lower than that in other patterns. These observations suggest that management units

with lower stand volumes tend to be selected in the scheduling process of the dispersed pattern. On the contrary, in

scheduling the clumped pattern, management units with

higher stand volumes tend to be selected and scheduled 

for clearcut treatments.

DISCUSSION

We have described heuristic techniques for distributing

management activities across a landscape in specific patterns, while seeking to achieve a commodity production

goal. Based on preliminary observations of the results of

these scheduling processes, we feel we have an adequate

system for achieving multiple goals across large landscapes

and long time frames. In future work, we will use FARSITE

(Finney 1999) to spatially model fires on the landscape

originating from a number of ignition points. Characteristics

213

of resulting fires such as fire size, fire line intensity, and

flame length will be measured, and then statistically compared to quantify the effect of each spatial pattern of management activity on the risk of intense wildfire. In addition,

after preliminary tests of the entire modeling methodology,

we will expand this research to a larger landscape. Several

issues remain to be resolved in the further study. For example, in scheduling a regular pattern, the efficient interval

between systematic points for reducing the risk of wildfire

is not known and the pattern was not optimized through the

entire time horizon. Thus, more attention to these issues

will be applied in the further study.

LITERATURE CITED:

Bettinger, P.; Graetz, D.; Boston, K.; Sessions, J.; 

Chung, W. 2002. Eight heuristic planning techniques

applied to three increasingly difficult wildlife planning

problems. Silva Fennica 36: 561-584.

Bettinger, P.; Sessions, J.; Boston, K. 1997. Using Tabu 

search to schedule timber harvests subject to spatial

wildlife goals for big game. Ecological Modelling

94: 111-123.

Bettinger, P.; Sessions, J.; Johnson, K.N. 1998. Ensuring 

the compatibility of aquatic habitat and commodity production goals in eastern Oregon with a Tabu search procedure. Forest Science 44: 96-112.

Catchpole, E.A.; Hatton, T.J.; Catchpole, W.R. 1989.

Fire spread through nonhomogeneous fuel modeled as 

a Markov process. Ecological Modeling 48: 101-112.

Dueck, G. 1993. New optimization heuristics: The great 

deluge algorithm and the record-to-record travel.

Journal of Computational Physics 104: 86-92.

Dunn, A.T. 1989. The effects of prescribed burning on 

fire hazard in the chaparral: toward a new conceptual

synthesis. In: Berg, N.H. (Tech. Coord). Proc. of the

Symp. on Fire and Watershed Management. USDA

Forest Service, Pacific Southwest Research Station,

Berkeley, CA. General Technical Report PSW-109. 

23-29 p.

Finney, M.A. 2001. Design of regular landscape fuel treatment patterns for modifying fire growth and behavior.

Forest Science 47: 219-228.
Finney, M.A. 2003. Calculation of fire spread rates across 

random landscapes. International Journal of Wildland

Fire 12: 167-174.

Fujioka, F.M. 1985. Estimating wildland fire rate of 

spread in a spatially non-uniform environment. Forest

Science 31: 21-29.

Glover, F. 1989. Tabu search - Part I. ORSA Journal on 

Computing 1: 190-206.

Hirte, S.R. 2002. Comparison of the value of forest plans 

developed under three levels of detail regarding management intentions. Masters Thesis, Oregon State

University, Corvallis, OR.

Jones, J.G.; Chew, J.D. 1999. Applying simulation and 

optimization to evaluate the effectiveness of fuel treatments for different fuel conditions at landscape scales.

In: Neuenschwander, L.F.; Ryan, K.C. (tech. eds). Proc.

joint fire science conference and workshop. p. 89-95.

Keane, R.E.; Morgan, P.; Running, S.W. 1997. FIREBGC- A mechanistic ecological process model for simulation fire succession on coniferous forest landscapes of

the northern Rocky Mountains. USDA Forest Service,

Intermountain Research Station, Ogden, UT. General

Technical Report INT-484.

Loudon, K. 1999. Mastering Algorithms with C (1st

edition). Oram, A. (ed.). O’Reilly & Associates, Inc.,

Sebstopol, CA. 540 p.

Martin, R.E. 1988. Rate of spread calculation for two 

fuels. Western Journal of Applied Forestry 3: 54-55.

Mladenoff, D.; He, H.S. 1999. Design, behavior and 

application of LANDIS, an object-oriented model 

of forest landscape disturbance and succession. In:

Mladenoff, D.; Baker, W. (eds.). Spatial Modeling of

Forest Landscape Change: approaches and applications.

Cambridge University Press, UK. 125-162 p. 

Murray, A.T. 1999. Spatial restrictions in harvest scheduling. Forest Science 45: 45-52.

Parsons, D.J.; van Wagtendonk, J.W. 1996. Fire research 

and management in the Sierra Nevada. In: Halvorson,

W.L.; Davis, G.E. (eds). Science and ecosystem management in the National Parks. Univ. Ariz. Press,

Tucson.

van Wagtendonk, J.W. 1995. Large fires in wilderness 

areas. In: Brown, J.K.; Mutch, R.W.; Spoon, C.W.;

Wakimoto, R.H. (Tech. Coords.). Proceedings of a symposium on fire in wilderness and park management.

USDA Forest Service, Intermountain Research Station,

Ogden, UT. General Technical Report INT-GTR-320. 

p. 113-116.

214
In: Bevers, Michael; Barrett, Tara M., comps. 2005. Systems Analysis in Forest Resources: Proceedings of the 2003 Symposium. General Technical Report 

PNW-GTR-656. Portland, OR: U.S. Department of Agriculture, Forest Service, Pacific Northwest Research Station.

1 Alan A. Ager is an Operations Research Analyst, Umatilla National Forest, Pendleton, OR. 2 James Barbour is Program Manager, Focused Science 

Delivery Program, Pacific Northwest Station. 3 Jane L. Hayes is a Research Biological Scientist, Managing Disturbance Regimes Program, Pacific

Northwest Research Station, LaGrande, OR.

Systems Analysis in Forest Resources: Proceedings of the 2003 Symposium

SIMULATING FUEL REDUCTION SCENARIOS 

ON A WILDLAND-URBAN INTERFACE 

 IN NORTHEASTERN OREGON

Alan A. Ager1, R. James Barbour2, and Jane L. Hayes3

INTRODUCTION

Fuels reduction programs have been initiated on many

areas throughout the western US, even though the long-term

effects of these treatments for reducing severe wildfires on

large landscape has not received much attention. While case

studies show that stand-level treatments can reduce fire

severity (Fule and others 2001; Kalabodkidis and Omi 1998;

Pollet and Omi 2002; Stephens 1998), the effect of treatments

on large forested landscapes over time and the consideration of other resource constraints remains an experimental

topic (Finney 2001; Johnson and others 1998; Barbour and

others 2004). Also problematic are the mechanics of using

existing stand-based fire and vegetation simulation models

(Wykoff and others 1982; Reinhardt and Crookston 2003)

on large landscapes and extracting landscape-scale measures of fire behavior that are meaningful to the general

public. Methods for the latter are in wide demand for use

in project-level planning for fuels reduction treatments on

Federal and other lands.

ABSTRACT

We analyzed the long-term effects of fuels reduction treatments around a wildland-urban interface located in the Blue

Mountains near La Grande, Oregon. The study area is targeted for fuels reduction treatments on both private and federal

lands to reduce the risk of severe wildfire and associated damage to property and homes. We modeled a number of hypothetical fuel treatment scenarios and examined the resulting changes in fuel characteristics, fire potential, and stand structure

over time. Aggressive thinning alternatives showed significant reductions in stand characteristics that contribute to severe

crown fires, such as height to live crown and crown bulk density for the landscape as a whole. However, simulations with

extensive thinning showed larger overall flame lengths and torching compared to a no management scenario. Significant

changes in stand structure and other characteristics were noted for the thinning versus no management scenarios. Work is

ongoing to refine the simulation methods and test a wider range of treatment alternatives. The study motivated a discussion

of the long-term problem of managing forest fuels in areas like the Blue Mountains.

215

As part of a larger project to examine landscape management issues in the Blue Mountains of Eastern Oregon

(Hayes and others 2004) we simulated a number of fuel

management scenarios on a key wildland-urban interface

near La Grande, Oregon. The Mt. Emily area and surrounding community was identified in the National Fire Plan as

“high risk” due to the intermingling of homes and vegetation, potential fire behavior, and existing fire protection

capabilities (Wallace 2003). The forests have accumulated

high loadings of live and dead fuels after decades of fire

exclusion and a large spruce budworm epidemic in the 1980s

that resulted in extensive mortality (Quigley and others

2001). A number of state, federal, and local agencies and

organizations are coordinating efforts to reduce stand density and ladder and surface fuels with the goal of reducing

flame heights, spotting, and crown fire potential, as well as

provide defensible space for fire-fighting crews to safely

approach future wildfires (Wallace 2003). 
The Mt. Emily project provided an opportunity to

explore a number of questions related to landscape analysis

of fuel treatment initiatives. We analyzed a broad set of

hypothetical management alternatives that are beyond those

being proposed for the Mt. Emily area in terms of treatment

extent and intensity, with the goal of examining a number of strategic issues including: 1) How often are treatments needed to maintain desired vegetation conditions;

2) How effective are different types of treatments; and 3)

What is the long-term effect of fuels treatments on fire

behavior, fuels, and other resources. We report preliminary results 

of our work to address these and related questions.

MATERIALS AND METHODS

Study Area

The Mt Emily wildland-urban interface (WUI) is a 16-

mile long area immediately north of La Grande, Oregon

where the forested slopes of Mt. Emily and adjacent ridges

descend to the agricultural lands in the Grande Ronde

Valley (fig 1). For analysis purposes, a boundary around

the area was established, containing 40,368 acres including

federal, state, and privately owned lands. About 30,348

acres within the study area are classified as forested lands

based on potential vegetation data. Approximately 25,139

acres within the larger analysis area are Forest Service

lands, administered by the Wallowa-Whitman (9,346

acres) and Umatilla (14,230 acres) National Forests.

Forest Service lands are managed for a number of objectives including big game winter and summer range, dedicated old growth and roadless. Surface fuel loading ranges

from 15 to 80 tons per acre with excessive dead ladder

fuels in a large number of the stands. Fuel accumulations

accelerated after the 1980s spruce budworm epidemic that

caused extensive mortality within the grand fir (Abies

grandis) and Douglas-fir (Pseudotsuga menziesii) stands

in the project area (Wallace 2003). The area has experienced 129 documented fire starts since 1970, including

both lightning and human caused fires. Ninety-nine percent

of the fires have been contained within the first twenty-four

hours, and only one fire has grown to any significant size

(Frizell fire, 250 acres) in recorded fire history. 

Treatment Scenarios

Fuels specialists on the Wallowa - Whitman National

Forest developed a range of different fuel treatments for

about 3,000 acres that prescribe: 1) mechanical thinning

216

Figure 1—Vicinity map of the Mt Emily Study area.
trees under 21 inches to reduce crown bulk density and 

ladder fuels; 2) Site removal of fuels when they exceed 25

tons per acre; 3) Underburning after treatment to reduce

surface fuels (Wallace 2003); and 4) Piling and burning of

surface fuels. The choice of treatment(s) for a given stand

was dependent on a number of factors related to stand conditions, management goals, and management restrictions,

and was the subject of considerable debate among specialists and managers. The treatments were strategically located

to reduce fire spread, create defensible space for suppression activities, and protect utilities like the Mt. Emily electronic site. 

For the purpose of examining long-term effects of 

fuel treatments, we used the general treatment methods

described above and simulated a broader, simplified set 

of seven hypothetical scenarios as follows: 

1. No Management or Natural Disturbance [NOMAN]. 

2. Long Term Aggressive Thinning [LTAT]. Thin from 

below all stands that exceed 55 percent of the maximum

stand density index (SDI) to 35 percent of the maximum SDI; site removal of fuels followed by underburning. Management is applied to all ownerships over time

as needed. 

3. Long Term Moderate Thinning [LTMT]. Same as 

LTAT except that thinning is triggered at 65 percent of

maximum SDI and trees are removed until 45 percent

of SDI is reached. 

4. Thinning with Maintenance Burns [LTATMBURN]. 

Same as LTAT with the addition of decadal underburning of stands that are thinned for the entire period of the

simulation.

5. Long Term Aggressive Thinning with Maintenance 

Burns on Forest Service lands only [LTATMBURNFS].

Same as LTATMBURN except that treatments are only

applied on Forest Service lands administered by the

Wallowa - Whitman National Forest. The remaining

Forest Service lands in the project area are in a roadless

area administered by the Umatilla National Forest and

were not considered for treatment in any of the scenarios. 

6. Long Term Aggressive Thinning with Maintenance 

Burns on Forest Service lands only, 21-inch diameter

limit [LTATMBURNFS21]. Same as LTATMBURNFS21 except that the Region 6 old growth screens that

prevent harvesting of trees greater than 21 inches is

implemented. 

7. Short Term Aggressive Thinning [STAT]. Same as 

LTAT except that stands are treated in the first decade

only of the simulation. Management is applied to all

ownerships.

We emphasize that these scenarios are hypothetical, and

were formulated to better understand how existing models

and typical Blue Mountain landscapes respond to long term

management scenarios. 

Vegetation and Fuels Data

We developed an initial vegetation database consisting

of data on trees per acre by species and 1-inch diameter

class for each stand using data from stand exams, photointerpretation, and field reconnaissance. The resolution of

the data was coarse in some areas, but adequately represented the project area for purposes of the present study.

Where stand exam data did not exist we generated approximate tree lists from photo-interpretation of 1:12000 color

resource photos from 1998. Photo-interpreted data consisted

of canopy closure estimates by size class and species. These

data were converted to trees per acre by diameter class

using relationships between crown diameter and tree diameter developed from stand exam data at the La Grande

Ranger District. Extensive field checks of the data were

completed over the course of the study. 

Detailed fuels data were obtained from stand exams 

and supplemented with a sub-sample of stands that were

chosen to represent average fuel conditions for different

fire regimes and commonly observed fuel models (Wallace

2003, Anderson 1982). Fuels were sampled with line transects according to Hilbruner and Wordell (1992) on a

sample of 6 stands chosen to represent dominant fuel conditions. These data were used to adjust the initial fuels and

default fuel loadings for the dense conifer stands (fuel

models 8 and 10, Anderson 1982). Default values (Reinhardt

and Crookston 2003; Anderson 1982) were used for the

remaining fuel models found in the project area. 

Models

For each of the 1060 stands in the project area we projected vegetation with the Blue Mountains Variant of the

Forest Vegetation Simulator (FVS) (Wyckoff and others

1982), a distance-independent individual tree growth model.

We simulated decadal time steps and report simulations for

60 years. We used a forest regeneration model (Wilson and

Maguire in prep.) that was developed from extensive seedling

survey data obtained from the La Grande Ranger District.

Mechanical thinning was simulated within FVS using 

prescriptions developed for thinning studies in the Blue

Mountains (Rainville 2002). Specifically, stands were

thinned once they exceeded 55 percent of maximum SDI,

thinning from below until the stand SDI was 35 percent of

maximum. Maximum SDI values were specific to each plant

217
association and target species (Cochran and others 1994).

The thinning prescriptions targeted removal of late-seral

species like grand fir in mixed-species stands, favoring

early seral species like western larch (Larix occidentalis)

and ponderosa pine (Pinus ponderosa). Although the prescriptions used were relatively generic for the diversity of

ecological settings and management goals in the project

area, they were adequate for the present study. 

Potential fire effects, dead and down fuel dynamics, 

and underburns were simulated with the Fire and Fuels

extension (FFE) to FVS (Reinhardt and Crookston 2003)

and FlamMap (Finney in prep). Due to space limitations

the FlamMap results are not presented here. The FFE

FUELMOVE keyword was used to simulate removal of the

larger surface fuels, under the assumption that 90 percent

of the 3 to 6 inch and 40 percent of the 1 to 3 inch materials

would be removed from the site. Underburning followed

thinning and fuel removal treatments, and was simulated

with the FFE SIMFIRE keyword using typical weather for

burning treatments on the La Grande Ranger District

(table 1). Potential wildfire effects on each stand were

simulated using the FFE POTFIRE keyword. Weather conditions for simulating potential wildfire effects were derived

from the J Ridge (Station 351414), Black Mountain (Station

351314), and Black Mountain 2 (Station 351317) remote

automated weather stations. The Black Mountain Stations

are located 10 miles east, and the J-Ridge stations are

located 25 miles south of the project area on national

forest land. Weather data for June to September from the

years 1986 to 2002 were analyzed in Fire - Family Plus

(Bradshaw and McCormick 2000) to generate 90th and 97th

percentile temperature, wind speed, and fuel moisture values (table 1). Outputs were analyzed for crown fire activity,

flame length, crowning and torching index, and changes in

fuel models over time. The 10-minute average wind speeds

generated in Fire- Family Plus were converted to maximum

velocity gusts using the tables developed by NOAA

(http://www.seawfo.noaa.gov/fire/olm/fire/10togust.htm). 

We measured forest structure (O’Hara and others 

1996) for the scenarios using the cover extension to FVS

(Crookston and Stage 1999). 

RESULTS

Stand Development 

The scenario that simulated long-term thinning of all

lands in the study area (LTAT) resulted in the treatment 

of about 23 percent (7,032 acres) of the forested area per

decade, averaged over the 60 year simulation (table 2). The

same prescription with decadal maintenance burns thinned

stands in about 15 percent of the study area (4,568 acres)

per decade (LTATMBURN). Thus the addition of the maintenance burns reduced the area thinned by 8 percent per the

decade. The LTAT scenario removed an average 6,394 MCF

of merchantable material per decade (table 2), with the

other scenarios resulting in lesser values commensurate

with the acres treated. Thinning only Forest Service lands

within the study area (LTATMBURNFS scenario) called

for thinning only about 1,212 acres, and removed 1,154

MCF per decade (table 2). Thinning to the upper SDI management zone (LTMT) called for thinning about 6,291

acres per decade, or 21 percent of the forested area (table

2). For the landscape as a whole, the overall SDI was

maintained at 135 for the LTAT scenario, while SDI for the

no management scenario rapidly increased and eventually

tapered off at 215 over the 60-year simulation period (table

2). The thinning schedule over time was most irregular for

the LTAT scenario, and it was apparent that heavy thinning

in the beginning of the simulation produced pulses of

regeneration that required thinning in the future. The pulse

of regeneration in the first decade is consistent with field

observations. A large number of stands in the cool moist

ecological settings have experienced severe mortality from

spruce budworm, which has opened up the stands and set

the stage for regeneration. 

The LTAT and LTATMBURN thinning management

scenarios resulted in large changes in crown bulk density

and crown base height (fig. 2). By 2020 the average crown

bulk density under both scenarios was reduced to 0.04

kg/m3 compared to almost 0.1 kg/m3 for the no management scenario. Crown base height was higher for the thinning scenarios, the difference being almost 15 feet by the

year 2020. Trends in both crown bulk density and live

crown height were strongly affected by regeneration in 

the early decades. 

Stand structural characteristics showed a number of

changes over time and differences among the scenarios

were apparent (fig. 3). In the no management scenario,

simulations indicated large increases in the stem exclusion

(SE) structural stage, and by 2030, over half of the forested

lands were in this structural stage. The effect of thinning

within the LTAT or LTATMBURN scenarios was to convert

the stem exclusion stage to single-stratum old forest. Thinning also increased the proportion of stands in the understory re-initiation stage as well (UR, fig 3). 

The STAT scenarios that called for only a single thinning at the beginning of the simulation had little effect on

the landscape through time (table 2). Likewise, recurrent

218
219

Table 1—Weather and fuel moisture parameters used in potential fire simulations and prescribed fires. The 

90th and 97th percentile conditions were determined from analysis of local weather stations. See methods for

additional details on weather station data and calculation of percentile values and velocity of wind gusts. 

Variable 90th Percentile Values 97th Percentile Values Maintenance Burns

Temperature (degrees F) 85o 90o 70o

1 hour fuel moisture (percent) 4 3 6

10 hour fuel moisture (percent) 6 4 12

100 hour fuel moisture (percent) 7 6 16

1,000 hour fuel moisture (percent) 9 7 26

Live fuel moisture (percent) 74 64 100

Duff moisture (percent) 42 20 —

10-minute average wind speed 

(mph) at 20 feet above ground 8 9 —

Maximum gusts 23 24 —

thinning and underburning of only Forest Service lands

(LTATMBURNFS) showed relatively minor changes compared to the no action alternative (table 2) which could be

expected given that only about 25% of the area is available

for treatment in this scenario.

Fuels

Changes in fuel models among the scenarios generally

followed expected trends (fig. 4). For the no management

scenario, there was a rapid decrease in fuel model 5 acres

(brush type), an increase in model 8 acres (closed canopy

conifer), and, after 50 years, a dramatic increase in fuel

model 10 acres (closed canopy with high surface fuels)

(fig. 4). The change in fuel models reflected a landscape

with active regeneration in the early part of the simulation,

with stands building up down fuels and increasing mortality through time. In contrast, the LTAT alternative showed 

a large reduction in fuel model 8 as stands were thinned

and converted to fuel models 2 and 5 (fig. 4). There was

also an increase in fuel model 11 (closed canopy with high

activity fuels) in the first 20 years, from thinning residue.

The thinning treatments did not result in an overall decrease

in surface fuels (table 2) even though we assumed that a

large portion of the residual fuels (90 percent of the 3 to 6

inch and 40 percent of the 1 to 3 inch materials would be

removed from the site as part of the treatments. Scenarios

that called for lighter thinning (LTMT) and 10-year maintenance burns did not significantly reduce downed fuels

either when measured as an average for the landscape.

Potential fire analysis

Acres of potential active crown fire increased in the first

30 years for all scenarios, with the no management scenario exceeding all other scenarios (fig. 5). For the LTAT

scenario, acres increased from about 500 to 3000 acres

between years 2000 to 2020, declining thereafter, and at

2030, acres of potential active crown fires were reduced to

nearly 0. Acres of potential active crown fire also showed 

a decreasing trend for the no management scenario after

about 30 years, although the acres with a crowning index

less than 24 (meaning that they will carry a crown fire in

24 mph wind) remained high throughout the simulation

while rapidly decreasing for the thinning scenarios (table

2). Although the LTAT and LTATMBURN simulations

showed reduced acres of potential active crown fire, they

resulted in a significant increase in passive crown fires

(torching) throughout the simulation (fig. 5). The increase

in torching was not simply the net result of decreased active

crown fire (fig 5.) since the increases were much larger than

the decrease in active crown fires. Also, the future trend in

passive crown fire acres continued to increase in the LTAT

scenario after the acreages of active crown were negligible.

Passive crown fire for the no management scenario initially

decreased until 2030 and stayed relatively constant thereafter. Increased torching in the aggressive thinning alternatives was caused by the increased fuel loadings and acreages

of the grass and brush fuel models, which have a higher

inherent flame length (Anderson 1982). 

DISCUSSION

Our simulations oversimplified the design of fuel treatment scenarios on the Mt. Emily landscape by not considering the spatial arrangement of treatments over time

(Finney 2001, Finney and Cohen 2002). We also relied primarily on landscape averages to present the results, which

can mask important trends at the stand level. In addition,

there are many aspects of the FVS Fire and Fuels Extension
220Table 2—Results of simulations for selected management scenarios examined in the study. See text for definition of scenario labels. Variable CI24 is the acres that had a crowning index below 24, meaning that the stand can carry a crown fire in the presence of 24 miles per hour

wind (97th percentile gust wind speed, Table 1).

NOMAN LTAT LTATMBURN LTATMBURNFS

Year Year Year Year

2000 2020 2040 2060 2000 2020 2040 2060 2000 2020 2040 2060 2000 2020 2040 2060

Acres

thinnedper decade 0 0 0 0 3699 6566 4825 6773 3699 5614 3356 4461 1442 1377 403 775

Stand

densityindex 95 199 256 281 95 136 149 123 95 131 123 134 95 183 217 222

Decadal

thinningvolume(merch.MCF) 0 0 0 0 3973 3389 3234 8475 3973 2759 1904 5077 1630 923 277 736

Flame

lengthmoderatefire (ft) 3.2 1.8 1.9 2.3 3.4 3.1 3.1 3.5 3.4 3.1 3.3 3.5 3.2 2.2 2.3 2.6

Flame

length

severefire (ft) 11.7 11.4 10.2 9.8 12.1 11.7 11.8 12.0 12.1 11.6 11.4 11.8 11.8 11.4 10.8 9.9

CI24

(acres) 10663 20148 14958 14566 9761 12615 3652 833 9761 12600 2905 588 10355 18760 12330 11632
that have not undergone rigorous examination, and thus

results we present pertaining to fuels and potential fire

need to be viewed with some caution. One area in particular that needs attention is the assignment of fuel model 5 to

stands that are dominated by conifer regeneration. These

stands may or may not behave like the fuel model 5 in a

wildfire depending on species composition and other factors

not considered in the model. Nevertheless, the effect of

thinning treatments had marked landscape-level effects on

parameters like crown bulk density and live crown height,

221

Figure 2—Results of simulations showing changes over time in crown bulk density (top) and live crown height (bottom) for selected

scenarios examined in the study.

Year

Year
both of which are strong determinants of crown fire activity.

For instance, the average crown bulk density after 30 years

of treatment was 0.01 kg/m3, which is below the 0.050

kg/m3 level required to sustain active crown fires (Keyes

and O’Hara 2002). These changes were reflected in the

lower proportion of the stands that exhibited active crown

fire behavior in the potential fire simulations, and later on,

a decrease in the crowning index. It was interesting to find

222

Figure 3—Results of simulations showing changes over time in forest structure for the no management (NOMAN, top) and longterm aggressive thinning (LTAT, bottom) scenarios. Structure classifications follow O’Hara et al. (1996) and were determined with

the FVS cover extension (Crookston and Stage 1999).
that acreage of active crown fire sharply decreased after 40

years for the no management scenario. Examination of the

acreages that could sustain a crown fire under 97th percentile weather conditions for the no management scenario

(table 2) showed that while crown fire could not be initiated,

a large portion of the project area could sustain a crown

fire spreading from an adjacent stand. If we had simulated

small disturbances like windthrow and root disease pockets,

regeneration would probably provide the ladder fuels needed

to initiate a crown fire and alter the results. 

223

Figure 4—Results of simulations showing changes over time in fuel models for the no management (NOMAN, top) and long-term

aggressive thinning (LTAT, bottom) scenarios.

Year

Year
It was also interesting to find that thinning resulted in a

large increase in passive crown fires (torching) and a small

increase in average flame length for the landscape as a

whole. Thinning treatments resulted in a larger proportion

of the area assigned to fuel models 2 and 5 (grass, brush),

which by default, have higher flame lengths and spread rates,

contributing to a higher potential for torching (Anderson

1982). In addition, activity fuels most likely contributed

224

Figure 5—Results of simulations showing changes over time in acres of potential passive (top) and active (bottom) crown fire for

the no management (NOMAN) and long-term aggressive thinning (LTAT) scenarios.

Year

Year
to the torching potential. Although the fuel models 2 and 5

have relatively high flame lengths and spread rates, they

burn at low intensities relative to crown fires, and fire suppression in these fuel types is more effective (Finney

2001). A number of empirical observations have shown

that treatments like those simulated here lower the intensity, spread rates, and tree mortality from wildfires. 

The scenarios that simulated treatment of the Forest

Service lands only showed a relatively minor response for

the study area as a whole (table 2), although in part this

result is strongly influenced by the size of the study area

relative to the Forest Service land base. In any case, the

LTATMBURNFS scenario reinforces the importance of

cooperation among landowners to bring about fuels treatments on all areas. It is important to recognize that we did

not attempt to model an array of resource constraints on

Forest Service lands that would further limit treatments. 

The mechanical thinning scenarios provide some indication of future treatment schedules for landscapes like the

Mt. Emily WUI, which is useful information for strategic

planning efforts that are underway on the National Forests

in the Blue Mountains province. Treatment schedules for

different scenarios can be used to forecast fiber supply and

net revenues of thinning treatments over the long run. Under

a scenario that calls for thinning based on a relative density

measure like SDI, a shorthand method to calculate treatment rates can be approximated as

Treated acres/year = (Watershed Acres)/ (F’[SDImax] -

F’[SDImin])

Where SDImax and SDImin are the values of SDI at 

the upper and lower management threshold and F’(t) is the

inverse of the function that predicts SDI change over time.

The latter can be derived from stand growth simulations

and averaged for all stands in a landscape. 

The current and future density problems in areas like Mt.

Emily demand additional ecological discussions (Tiedeman

and others 2000) and perhaps new and innovative solutions.

Thinning, fuels treatments, and underburning have welldocumented beneficial effects on overstocked stands in

terms of general forest health and wildfire risk. However,

thinning overstocked stands induces waves of regeneration

and sets the stage for future development of ladder fuels

and high crown bulk density. If recurring treatments are not

applied, the initial thinning will be counter to long-term

landscape goals (Keyes and O’Hara 2002). The long - term

treatment of fuels is made difficult by the finding that most

restoration treatments in the Blue Mountains have negative

net values under current economic conditions (Rainville

2000) and a variety of resource values must be considered

on the Federal land base. 

While fire suppression and removal of fire tolerant

species has been a major factor in the development of fireprone conditions, there are many processes besides wildfire

that regulate stand density through time in coniferous forests.

For instance, thinning treatments can enhance seed production in residual trees (Reukema 1961), leave a favorable

substrate for seed germination (Schopmeyer 1974), and

degrade habitat for small mammals (Tiedemann and others

2000) that consume conifer seed (Sylvester and others 1994).

Thus treatments may set the stage for prolific conifer regeneration and, eventually, dense stands with high crown fire

potential. New research is needed to find ways to manage

seed production, seed survival, germination and seedling

establishment after underburning and mechanical treatments. 

There are many issues that need to be explored with the

simulation approach we have used here, such as the spatiotemporal patterns of treatments (Finney 2001), interactions

between wildfire and other disturbances (Quigley and 

others 2001), and secondary effects of fuels treatments on

other resources including wildlife and forest productivity

(Tiedemann and others 2000; Johnson and Miyanishi 1995).

For instance, optimal spatial patterns of shaded fuel breaks

(Finney 2001) could fragment stands of late-old structure

considered important to keystone wildlife species like

American marten (Martes americana) (Hargis and others

1999). On the other hand, these treatment patterns might

create desired foraging habitat for Canada lynx (Lynx

canadensis) (Aubry and others 2000) and Rocky Mountain

elk (Cervus elaphus). Frequent burning can reduce foraging substrate (down wood) for a number of important

avian and mammal species that feed on ants and other

insects (Torgersen and Bull 1995). Economic questions

about the kinds of investments that will be needed in the

long term to finance forest restoration treatments (Christiansen

and others 2002) have not been addressed in a way that

reflect spatially explicit harvesting and transportation costs.

All of these issues contribute to a long-term landscape

design questions of how to achieve restoration goals in the

Blue Mountains (Quigley and others 2001) within a constraint matrix that includes riparian buffers, roadless area,

visual objectives, and habitat for listed species. 

Many improvements are needed to the simulations we

presented here, and work is in progress towards that end.

We are analyzing the different scenarios in terms of wood

utilization, habitat for selected wildlife species, and insect

225
and disease risks for the Mt Emily WUI. We are also developing mortality models for insects and disease using a multiagent approach (Roberts and Weatherby 1997). Mechanical

thinning prescriptions are being refined for the different

ecological settings and management goals within the project area. These refinements will be used for additional simulations on the Mt Emily area to gain further insights into

the long-term problem of managing for multiple resource

objectives in the Blue Mountains. 

ACKNOWLEDGEMENTS

We are grateful to many people who have provided 

help with this study, including Willie Crippen, Jim Barrett,

Tom Burry, Nick Crookston, Rob Seli, and Chuck McHugh.

This work was funded by the Focused Science Delivery

Program, Pacific Northwest Research Station, the Blue

Mountain Demonstration Program, and the Joint Fire program. 

LITERATURE CITED

Anderson, H.E. 1982. Aids to determining fuel models for 

estimating fire behavior. Gen. Tech. Rep. INT-122.

Ogden, UT: U.S. Department of Agriculture Forest

Service, Intermountain Forest and Range Experiment

Station, 22p. 

Aubry, K.B.; Koehler, G.M.; Squires, J.R. 2000. Ecology 

of Canada lynx in southern boreal forests. Pages 373-

396.In Ruggerio, L.R.; Aubry, K.B.; Buskirk, S.W.;

Koehler, G.M.; Krebs, C.J.; Kelvey, K.S.; Squires, J.R.

(editors). Ecology and conservation of lynx in the

United States. University Press of Colorado, Boulder.

Bradshaw, L; McCormick, E. 2000. FireFamily Plus 

user’s guide, version 2.0. Gen. Tech. Rep. RMRS-GTR67WWW. Ogden, UT: U.S. Department of Agriculture,

Forest Service, Rocky Mountain Research Station.

Barbour, R.J.; Ager, A.A.; Hayes, J.L. 2004. A framwork 

for the development and application of INLAS: the

Interior Northwest Landscape Analsis System. Chapter

1 in Gen. Tech. Rep. GTR-PNW-610, Portland, OR: 

US Department of Agriculture, Forest Service, Pacific

Northwest Research Station. 218 p.

Christiansen, G.A.; Fight, R.; Barbour, R.J. 2002.

Simulating fire hazard reduction, wood flows, and economics of fuel treatments with FVS, FEEMA, and FIA

data. In: Crookston, N.L., Havis, R. N., (comps.) Second

Forest Vegetation Simulator Conference. Proceeding

RMRS-P-25., Ft. Collins CO: US Department of

Agriculture Forest Service, Rocky Mountain Research

Station. Pp. 91-96.

Cochran, P.H.; Geist, J.M.; Clemens, D.L.; Clausnitzer, 

R.R.; Powell, D.C. 1994. Suggested stocking levels for

Forest Stands of Northeastern Oregon and Southeastern

Washington. Gen. Tech. Rep. PNW-RN-513. Portland,

OR: US Department of Agriculture, Forest Service

Pacific Northwest Research Station.

Crookston, N.L.; Stage, A.R. 1999. Percent canopy cover 

and stand structure statistics from the forest vegetation

simulator. Gen. Tech. Rep. RMRS-GTR-24. Ft. Collins,

CO: USDA Forest Service Rocky Mountain Research

Station.

Finney, M.A. 2001. Design of regular landscape fuel 

treatment patterns for modifying fire growth and 

behavior. Forest Science 47:219-228

Finney, M.A. (in prep.). FlamMap. Rocky Mountain 

Research Station, Fire Sciences Lab, Missoula MT.

Finney, M.A.; Cohen, J.D. 2002. Expectation and evaluation of fuel management objectives. In: Fire, Fuel

Treatments, and Ecological Restoration: Conference

Proceedings. US Department of Agriculture, Forest

Service, Rocky Mountain Research Station. General

Technical Report RMRS-P-29. p. 353 – 366.

Fule, P.Z.; McHugh, C.; Heinlein T.A.; Covington, 

W.W. 2001. Potential fire behavior is reduced following

forest restoration treatments. In. Vance, R.K.: Edminister,

C.B.: Covington, W.W.; Blake J.A. comps. Ponderosa

pine ecosystems restoration and conservation: steps

towards stewardship. RMRS-P-22. Ft. Collins, CO: 

US Department of Agriculture, Forest Service, Rocky

Mountain Research Station.

Hargis, C.D.; Bissonette, J.A.; Turner, D.A. 1999. The 

influence of forest fragmentation and landscape pattern

on American martens. Journal of Applied Ecology 

36: 157-172.

226
Hilbruner, M.W.; Wordell, T.A. 1992. Planar Intercept 

Fuel Inventory Field Guide. Adapted from the Handbook

for Inventorying Downed Woody Material and R-6

Stand Examination Program Field Procedures Guide.

USDA Forest Service Region 6, Portland OR.

Johnson E.A.; Miyanishi, K. 1995. The need for consideration of fire behavior and effects in prescribed burning.

Restoration Ecology 3:271-278.

Johnson, K.N.; Sessions, J.; Franklin, J.; Gabriel, J. 

1998. Integrating wildfire into strategic planning for

Sierra Nevada forests. Journal of Forestry 96(1): 42-49.

Kalabodkidis, K.D.; Omi, P.N. 1998. Reductions in fire 

hazard through thinning/residue disposal in the urban

interface. International Journal of Wildland Fire 8: 

29-35. 

Keyes C.R.; O’Hara, K.L. 2002. Quantifying stand targets

for silvicultural prevention of crown fires. Western

Journal of Applied Forestry 17: 101-109. 

O’Hara, K.L.; Latham, P.A.; Hessburg, P.; Smith, B.G. 

1996. A structural classification for northwest forest

vegetation. Western Journal of Applied Forestry 11: 

97-102.

Pollet, J.; Omi, P.N. 2002. Effect of thinning and prescribed burning on crown fire severity in ponderosa

pine forests. International Journal of Wildland Fire 11:

1-10. 

Quigley, T.M.; Hayes, J.L.; Starr, G.L. [and others]. 

2001. Improving forest health and productivity in eastern Oregon and Washington. Northwest Science 75:

234-251.

Rainville, R. Restoration opportunities in the Blue 

Mountains. Blue Mountains Demonstration Project.

USDA Forest Service, Pacific Northwest Forest and

Range Experiment Station, La Grande OR. Report 

on file. 66 p.

Reinhardt, E.; Crookston, N.L. 2003. The fire and fuels 

extension to the forest vegetation simulator. Gen. 

Tech. Rep. RMRS-GTR-116. Ft. Collins, CO: U.S.

Department of Agriculture, Forest Service, Rocky

Mountain Research Station. 209 p. 

Reukema, D.L. 1961. Seed production of Douglas-fir 

increased after thinning. Portland, OR: U.S. Department

of Agriculture, Forest Service, PNW Forest and Range

Experiment Station, Research Note 210. 5 p. 

Roberts, J.C.; Weatherby, J.C. 1997. Successional simulation using the forest vegetation simulator and multipest scenarios for various levels of planning. In: Teck,

R.; Moeur, M.; Adams, J., comps. Proceedings—Forest

vegetation simulator conference. Gen. Tech. Rep. INTGTR-373. Ogden, UT: U.S. Department of Agriculture,

Forest Service, Intermountain Research Station: 111-

115.

Schopmeyer, C.S. 1974. Seeds of Woody Plants in the 

United States. USDA Forest Service Agriculture

Handbook No. 450. 883 p. 

Stephens, S.L. 1998. Evaluation of the effects of silviculture and fuels treatment on potential fire behavior in

Sierra Nevada mixed conifer forests. Forest ecology 

and management 105: 21-35. 

Sylvester, A.W.; Gaud, W.S.; States, J.S. 1994. Effects of 

herbivory by Abert squirrels (Sciurus aberti) on cone

crops of ponderosa pine. Journal of Mammology 75:

700-710.

Tiedemann, A.R.; Klemmedson, J.O.; Bull, E.L. 2000. 

Solution of forest health problems with prescribed fire:

Are forest productivity and wildlife at risk? Forest

Ecology and Management 127: 1-18. 

Torgersen, T.R.; Bull, E.L. 1995. Down logs as habitat 

for forest-dwelling ants - the primary prey of pileated

woodpeckers in Northeast Oregon. Northwest Science

69: 294-303. 

Wallace, P.L. 2003. Existing fuels conditions on Mt. Emily 

project area. Unpublished report on file, LaGrande

Ranger District, LaGrande, OR.

Wilson, D.; Maguire, D. (in prep.). A regeneration model 

for the Upper Grande Ronde Watershed. Oregon State

University, Department of Forest Sciences, Corvallis,

OR.

Wykoff, W.R.; Crookston, N.L.; Stage, A.R. 1982. User’s 

guide to the Stand Prognosis Model. Gen. Tech. Rep.

INT-133, Ogden, UT: U.S. Department of Agriculture,

Forest Service, Intermountain Forest and Range

Experiment Station. 112 p.

227
This page is intentionally left blank.
STAND 

OPTIMIZATION

AND HARVEST

OPERATION
This page is intentionally left blank.
In: Bevers, Michael; Barrett, Tara M., comps. 2005. Systems Analysis in Forest Resources: Proceedings of the 2003 Symposium. General Technical Report 

PNW-GTR-656. Portland, OR: U.S. Department of Agriculture, Forest Service, Pacific Northwest Research Station.

1 P. Rauch is Assistant of Industrial Management and Logistics, Institute of Forest Sector Policy and Economics, BOKU - University of Natural Resources 

and Applied Life Sciences, Vienna, Feistmantelstrasse 4,1180 Wien, Austria.

Systems Analysis in Forest Resources: Proceedings of the 2003 Symposium

ESTIMATING THE HARVESTING POTENTIAL FOR HIGHLY

MECHANIZED HARVESTING SYSTEMS IN SMALL- SCALE

FORESTRY BY USING A GIS MODEL

Peter Rauch1

INTRODUCTION

There is an increasing demand of the wood industry for

small-diameter timber in Austria. On the one hand, there is

a large green stock, but on the other hand, the continuous

supply – especially out of the numerous small forest plots -

is a problem. The Austrian Forest Inventory shows tending

arrears with a harvesting potential of 40.5 million m3 for

Austria’s small-scale forests alone. Demand for smalldiameter timber could be satisfied if in small forest plots

more thinnings were done by mechanised harvesting systems.

Analyses of current data of the forest owners’ cooperative Leoben (Country of Styria, Austria) show that only

for a small part of the cooperative’s forest plots there are

forest stand data available and in most cases these data are

not up-to-date any more. Furthermore, basic cadastral data

of the members’ forest plots is also lacking. A similar situation holds true for most of the other Austrian forest owners’

cooperatives.

ABSTRACT

There is an increasing demand by the wood industry for small-diameter timber in Austria. On the one hand, there is a

large green stock, but on the other hand, the continuous supply – especially out of the numerous small forest plots - is a

problem. In Austria, forest owners’ co-operatives bundle the supply of round wood coming out of the small forest plots. The

only information the manager of a forest owners’ co-operative actually has are some personal data from the members of his

co-operative and the size of their forests. The main reason for this lack of data is that there are no forest management data

available.

For a forest owners’ co-operative with 240 members representing 9,000 ha (22,230 acre) of forest land a basic information system has been set up. The prototyping of the model has been done for a tracked harvester. The model estimates areas

where machinability of the stand and trafficability of the terrain both are fitting to the harvester configuration. Validations

of the model results show that it matches the reality of the mountainous research area quite well.

231

Aims—The project has pursued the following aims by

implementing a GIS model for estimating the harvesting

potential for highly mechanised harvesting systems in

small-scale forestry:

• Increase the amount of small-diameter timber 

harvested in the forest owners’ cooperative (FOC)

Leoben,

• reduce tending arrears,

• develop a useful basic forest map.

THE GIS MODEL

According to Reinberg (1998) the basic requirements 

of a model that is a simplified picture of reality are (i) its

similarity to the real system, (ii) being as simple as possible, (iii) having the possibility to adapt, and (iv) having 

a variable resolution. The demands for simplicity and for

similarity to the real system contradict each other to a
certain degree. But even complex real situations should not

lead to exorbitant complex models.

Spors and others (1992) name the following hard decision factors, i.e., factors that cannot be compensated, for

estimating potential areas for harvesters with GIS: danger

level of the site (site classification on the basis of substratum and water regime), inclination, d.b.h., and tree species

characteristics. Based on those data harvester areas can be

estimated on a regional scale.

Lan (2001) names the following parameters which

determine whether harvesters can be applied: (i) topographical conditions, such as inclination and the bearing capacity

of the soil, (ii) stand parameters, (iii) soil type, and (iv)

other environmental constraints. More models for intensive

timber production will be introduced in the near future

(Johnsen and others 2001) and the question where to use

which harvesting technique is one of the central tasks in

forest harvesting operations (Eichrodt 2003).

Synthesising the work of Spors and Lan one can say

that the application of a harvester depends on both trafficability and machinability of a forest area.

Adaptation of the model for the special conditions 

in small forest plots

For the FOC Leoben there is none of the geodata available named by Spors and Lan. The only exception is a digital terrain model. To estimate harvester areas under small

forest plot conditions some adaptations of the described

models have to be made. First of all, one should try to use

as much already available data as possible in the model.

This principle ensures that the model is adaptable to other

regions with small forest plot conditions. In their study

dealing with process models in forestry Johnsen and others

(2001) stress differences in data quality between complex

and data-intensive scientific models and applications for

forest managers which need data that is simpler and more

easily available. The model developed in this paper can be

assigned to the second type, the forest management applications, for two reasons: first, because of the (poor) quality

of the input data and, second, because of the intention to

support timber harvest in practise.

The model is a prototype adapted to the machine configuration of the FOC’s tracked harvester (Robin) and to

regional characteristics of the district of Leoben. The main

parameters of the machine configuration are gradeability

that mainly limits trafficable area and minimum and maximum d.b.h. suitable to the harvesting head.

For estimating harvester trafficability of the area of the

district of Leoben the digital terrain model of the Federal

Office of Metrology and Surveying (BEV) is the only one

available at the moment. Despite its wide raster (50 meters

between two measurement points) it is still good enough

to deduce topographical and morphological terrain data

(Schreier and Lavkulich 1979). Decision criteria inclination

(deduced from the digital terrain model) is as well a socalled hard factor (Spors and others 1992) as the main topographical condition for estimating harvester trafficability

(Lan 2001).

In this project I had to do without these data because

site classification maps were not available for the study

area and because mapping the scattered forest plots of the

FOC Leoben would have been too labour intensive. The

additional value of site classification data would not make

up for the additional costs (see Meyer and others 2001).

Caused by the weather, water content of soil and

depending soil parameters like Californian bearing rate

change permanently. Therefore the model works under the

assumption of gradeability under dry soil conditions.

Since for the FOC Leoben there were no forest stand

data available, the basic stand data needed for the model

had to be collected. In order to get up-to-date stand data in

an acceptable amount of time, the stages of development of

the forest have been determined by stereoscopic interpretation of aerial images. This technique has been used frequently and is known as field tested (Schmidtke and others

2000). The harvesting area was defined on the basis of the

maximum d.b.h. suitable to the harvesting head (Logmax

3000) of the tracked harvester used. With that stands with

d.b.h. of less than 20 cm (pole stage) as well as stands with

d.b.h. of up to 34 cm (small-diameter timber tree) were

selected as harvesting areas. Property borders have been

determined on the basis of the digital land register map from

BEV which is available for the whole district of Leoben.

Process Model—To create the GIS model for estimating harvester suitable stands the following process model

was used:

1. build up a forest owners’ database,

2. join the forest owners’ database with the digital land

register map,

3. estimate borders of stands and their stage of development as well as further stand data (tree species, tree

heights) by stereoscopic interpretation of aerial

images,

4. implement GIS model,

232
5. validate GIS model using terrestrially mapped stands

as reference.

Build up a forest owners’ database—For building up

a forest owners’ database (which should include the whole

forest area owned by members of the FOC Leoben) the 

following steps had to be taken: 

1. collect the cadastral basic data from each member

(lot codes),

2. determine land register data via internet on the basis

of lot codes,

3. put the results of the inquire in a forest owners’ database using hand made software and

4. check the results of the inquiry (e.g., cancel entries of

non-members, check completeness of data).

After having collected the cadastral basic data from

each member, land register data could be determined via

internet. Because of the quite complex land register system

it took several iterations to get all the required data for the

FOC’s members. Lot codes and owner data were put in a

database. 

Joining the forest owners’ database with the digital

land register map—In Austria each parcel of land has a

unique lot code. On the basis of this lot code the forest

owners’ database was joined with the digital land register

map. A layer of the wood lots of all members and the connection with the forest owners’ database are the results of

this operation. 

Stereoscopic interpretation of aerial images—In order

to get the stand data needed for the GIS model (stand borders and stage of development) in an acceptable amount of

time and for an acceptable amount of money, stereoscopic

interpretations of aerial images were carried out. The results

of the stereoscopic interpretations of aerial images are

based on expert’s opinion and because of subjective influences flawed with inaccuracy (Schmidtke and others 2000;

Kätsch 2001). Nevertheless the data quality expected should

be good enough for the modelling purpose at hand.

Digital, panchromatic orthophotos from BEV and real

colour aerial images from a private enterprise have been

used. The real colour photo flight includes 396 images and

has an average image scale of 1:19,000.

With a stereoscopic device (AvioPret from Wild) using

the stereo pairs of real colour images the borders of the

forest stands can be seen in three dimensions. As the real

colour images are not to scale the exact geometric delimitation was supported by orthophoto that was shown on a

computer monitor in ArcView. A digital land register map

was used to show the borders of the forest plots which have

to be edited. Table 1 shows the stages of development that

were used for describing the age phase of a stand.

After the delimitation further stand data (such as tree

height, canopy density and stem number) are estimated and

put in a database. Only stands with a minimum area of

1000 m2 were utilized because smaller areas could not

have been located exactly when they are drawn onto a

layer overlaying the orthophoto. 

Implementation of the GIS model—In order to get an

assessment of the harvester suitability of the stands the different layers of themes were joined with the Simple Additive

Weighting Method. This method is the most common for

multi-attribute decision support (Malczewski 1999). Application of this method in a GIS takes the following steps:

• define evaluation criteria and standardise layers,

• define evaluation weights for criteria and 

• generate resulting layer as score of weights per grid.

Modelbuilder, an extension to ArcView, supports this

method and is used to implement the model. All layers

(inclination, stage of development, storm damages and

other damages) are transformed into grid format. This,

finally, enables spatial analyses of harvester suitability of

the stands of the FOC Leoben. Assessment of input parameters with a solving algorithm is possible if clear relations

between assessment criteria are given (Czeranka 1999).

Following Spors and Lan the suitability of stands for harvesters is mainly limited by trafficability (here estimated

only by inclination) and stand parameters (here estimated

only by stage of development). Data of damaged plots are

digitally available for the whole country of Styria and are

used because these damages are another restriction for harvester suitability. Criteria weights and rules for connecting

input parameters are documented in table 2.

233

Table 1—Stages of development estimated by stereoscopic interpretation of aerial images

Stage of development Characteristic

Regeneration, culture Before canopy is closed 

Thicket stage d.b.h. less than 10 cm

Pole stage d.b.h. between 10 and 19 cm

Small-diameter timber tree d.b.h. from 20 to 34 cm

Medium timber tree d.b.h. from 35 to 50 cm

Big timber tree d.b.h. more than 50 cm
The implemented algorithm sums up products of percentage weights of a single parameter given in an input

theme and scale values for estimating the Gridcode for each

5x5m cell. Areas with attributes that disable a harvester (for

example stands “thinned” by storm) are excluded from the

possible solution space. Their scale value is set as restricted. A cell is assumed as harvester-suitable if the d.b.h. estimated by stereoscopic interpretation of aerial images as

well as trafficability assessed by inclination are fitting the

machine configuration of Robin and if at the same time the

stand is free of damages. 

For example if a 5x5m Grid is part of a stand in pole

stage with an inclination of this part of 53% and if there

are no storm or other biotic damages the Gridcode is estimated as follows: Gridcode for that area = (0.49x(-1)+

0.49x1+0.01x0+0.01x0).

Scale values of the input theme ‘stand grid’ are set so

that Gridcodes differ if logging can be done by a forwarder

(scale value -1) or only by a cable system (scale value 1).

So the models computes for harvester-suitable cells if logging can be done with a forwarder or a cable system (see

figure 1).

Finally, the output layer of the model has to be connected with the databases to provide – besides the information

of harvester suitability – all available data on forest stand

and owner. This can be done by using functionality of geometric analyses (for example: merge, intersect, union, spatial join). As a result, for each estimated harvester-suitable

thinning area stand data as well as land register data can be

shown (figure 2 ).

Validation—To validate the quality of the results parameters used in the model were measured in some reference

stands. The terrestrially mapped stands were digitalised and

parameters were put in a database. Stands were selected

only when they could be easily identified on orthophotos

and if corresponding colour images were of good quality.

234

Table 2—Model algorithm

Input theme % Influence Input field Input label Scale value

Stand Grid 49 Value

10 Forest roads, etc. Restricted

20 Regeneration Restricted

30 Thicket stage Restricted

40 Pole stage -1

51 Small-diameter timber tree -1

52 Medium timber tree Restricted

53 Big timber tree Restricted

60 Non forest area Restricted

No data No data Restricted

Trafficability 49 Value

1 Harvester+forwarder -1

2 Harvester+cable 1

3 Motor manuell+cable Restricted

No data No data Restricted

Storm Damage 1 Value

0 Peripheral area 0

1 Main area Restricted

No data No data 0

Other non biotic Damages 1 Value

3 hail 0

5 Forest fire Restricted

No data No data 0

Value trafficability is inclination:1 means inclination from 0 to 50 percent, 2 from 51 to 60 percent, 3 more than 60 percent. The column

“% Influence” gives the weight in percent that an input theme has in respect to the overall result. Each 5x5m grid with an input value

that eliminates a harvester to work is dropped from the solution area. This is done by setting their scale values as restricted. The solution

area is showing therefore harvester suitable grids only.
That was done in order to minimize errors of inaccuracy in

position when comparing terrestrial mapping with stereoscopic interpretation or with model results. A total of 77

stands representing 36 ha could be used to validate model

results for an area of 5,000 ha with about 2,500 different

stands. Furthermore all areas where the harvester had been

already working could be used for validation (in total 40 ha

extra). Due to the high effort to draw forest maps and to

measure parameters on site the sample size is restricted to

the area indicated.

Validation is split in (i) validation of model-estimated trafficability, (ii) validation of suitability to harvesting head to

235

Figure 1—Harvester suitable areas (Grid, section):

light gray areas are also forwarder suitable, dark

gray areas for cable extraction only.

Figure 2—Harvester suitable areas (Shape, section)

with stand and owner data shown for the light gray

stand area.
estimate accuracy of input data for these parameters and

(iii) validation of the model result.

Validation of suitability to harvesting head—In order

to reduce the influence of inaccuracy in position the comparison of results of stereoscopic interpretation of stage of

development (that were used for estimating average d.b.h.

of stand) with on-site measurement of average d.b.h. was

restricted to areas larger than 1,000 m2.

Stereoscopic interpretation estimated 45 percent of the

reference area exactly according to the stage of development as classified on site (figure 3). Further 38 percent

were estimated within a range of plus / minus one stage of

development. The most errors occurred in the differentiation between pole stage and small-diameter timber tree. But

since both of these stages of development fit to the harvesting head this error can be tolerated for the model’s purpose.

So if we combine pole stage and small diameter tree to one

category stereoscopic interpretation estimated 83 percent 

of the reference area correctly. This leads to the conclusion

that the method used is appropriate to get useful stand data

cheaply and quickly.

An overestimation for two stages of development

occurred for 10 percent of reference area. This happened

because of light and open structures of stands that leads to

a canopy surface which is normally characteristic for older

stands. Stands overestimated by three or four stages of

development resulted from inaccuracy in position between

stand borders digitised on orthophoto and seen on stereoscopic images.

Validation of harvester trafficability—The model

reveals trafficability for 77 percent of the reference area

according to terrestrial measurements. In addition, the

model points out all areas where the harvester has already

been working as harvester trafficable. On the one hand,

errors occur because of inaccuracies of the digital terrain

model (mostly underestimation of inclination) and on the

other hand because some areas are not trafficable even if

the inclination is suitable. Observed underestimation of

inclination derived from the digital terrain model is in

accordance with results of other GIS models using inclination as derived parameter (Wilson and others 2000).

Validation of harvester and forwarder trafficability

—Trying to estimate if a forwarder can be used or not fails

in 52 percent. The limits of trafficability of a forwarder are

easier reached than those of a tracked harvester. Slight

cross slope and barriers as well as steep acclivities stop a

fully loaded forwarder. All parameters mentioned are not

included in the actual model and therefore the model is not

able to estimate a suitable logging system.

Figure 4 shows the comparison of the forwarder trafficability estimated on one hand by the GIS model (first digit)

and on the other hand by terrestrial judgement (second

digit). Gridcode 1 means that the area is suitable for harvester and forwarder, Gridcode 2 that it is only harvester

suitable and Gridcode 3 that is neither suitable for harvester nor for forwarder. So for example a Gridcode combination of 12 means that the GIS model estimates that

area as harvester and forwarder trafficable and the terrestrial judgement says it is only harvester trafficable.

Validation of the GIS model results—The described

GIS model estimates harvester suitability correctly for 71

percent of the reference stands. This rate of accuracy identifies the model as suitable for use in practise (at least for

the area of investigation). Errors are to equal parts errors in

estimating trafficability or suitability for the harvesting

head or a combination of both.

236

Figure 3—Difference between stereoscopic interpretation and terrestrial

measurement of stage of development.

Figure 4—Comparison of harvester and forwarder trafficability estimated

by the GIS model with terrestrial judgement.
CONCLUSIONS

The FOC Leoben can use the model results as a marketing instrument to find owners of large harvester suitable

thinning areas. To the owners of these areas the possibilities and advantages of mechanised harvesting systems can

be pointed out. For a contracted thinning area some more

nearby potential harvester areas can be found and their

owners contacted. The model has a further surplus value

because it provides digital maps of forest stands of the

FOC Leoben for the first time as a basis for various planning tasks.

More accurate basic data (like a better digital terrain

model, soil parameters and actual stand data) would be

needed to estimate forwarder suitability. But collecting

these data would cause significant additional costs. More

exact terrain models generated by laser scanning with standard deviation in height inaccuracy of less than 15 cm in

comparison to terrestrial measurements (Wever 1999) or

plus /minus 25 cm in flat wooden areas (Pfeiffer and others

1999) could help to estimate forwarder suitability.

Even automated identification of stage of development

and stands could be done by comparing crown surface

model gained from the first pulse mode of laser scanning

with last pulse mode information that usually describes terrain surface (Ziegler and others 2000). Another possibility

to automate stand identification is to use spectral analyses

from Landsat TM EWDI (Enhanced Wetness Difference

Imagery) as is done for monitoring forest harvesting and

thinning activities (Franklin and others 2001).

In the future, these technologies could be useful to 

estimate suitability of mechanised harvesting systems in a

more automated way.

LITERATURE CITED

Czeranka, M. 1999. SDSS Spatial Data Descision Sytems. 

http://www.gis-tutor.de/theorie/grundlag/sdss/kap1.htm

ff., cited 16.10.2002.

Eichrodt, A. 2003. Development of a spatial trafficability 

evaluation system. Dissertation, Technische

Wissenschaften ETH Zürich, Nr. 15009; URL: http://ecollection.ethbib.ethz.ch/show?type=diss&nr=15009

Franklin, S.E.; Lavigne, M.B.; Moskal, L.M., Wulder

M.A.; McCaffrey, T.M. 2001. Interpretation of forest

harvest conditions in New Brunswick using Landsat

TM Enhanced Wetness Difference Imagery (EWDI).

Canadian Journal of Remote Sensing 27: 118-128.

Johnsen, K.; Samuelson, L.; Teskey, R.; McNulty, S.; 

Fox, T. 2001. Process Models as Tools in Forestry

Research and Management. Forest Science 47(1): 2-8.

Kätsch, C. 2001. Stand und Entwicklungsmöglichkeiten 

der automatischen Objekterkennung auf

Fernerkundungsaufzeichnungen für Zwecke der Waldund Landschaftsinventur. Forstarchiv 72: 244-250.

Lan, Z. 2001. Heuristic harvesting System Selection 

through the KBS Approach. University of Helsinki,

Department of Forest Resource, Management

Publications 27, 112.

Malczewski, J. 1999. GIS and multicriteria desicion analysis. John Wiley & Sons, New York, Chichester, 392.

Meyer, T.; Stückelberger, J.; Hollenstein, K.; Attenberg, 

M.; Hänggli, T. 2001. effor2: Holzerntekonzept für die

Untersuchungsbetriebe im Kanton Schwyz. ETHZ,

Professur Forstliches Ingenieurwesen, Interne Berichte

Nr. 11. 98.

Pfeiffer, N.; Kraus, K.; Köstli, A. 1999. Restitution of 

Airborne Laser Scanner Data in Wooded Areas. GIS.

Geoinformationssysteme. Zeitschrift für raumbezogene

Informationen und Entscheidungen. 12,2: 18-21.

Reinberg, S. 1998. Räumliche 

Entscheidungsfindungsprozesse und Modellbildung. In:

Strobl et Dollinger (Hrsg.): Angewandte Geographische

Informationsverarbeitung. Beiträge zur AGITSymposium Salzburg 98, 296-278.

Schmidtke, H.; Budmiger, P.; Tomisawa, S. 2000. 

Interpretation von Stichprobendaten aus digitalen Ortho

photos. AFZ/Der Wald 55, 4: 212-214.

Schreier, H.; Lavkulich, L.M. 1979. A Numerical 

Approach to Terrain Analysis for Off-Road

Trafficability. Photogrammetric Engineering and

Remote Sensing, Journal of American Society of

Photogrammetry, ASP, 45, 5: 635-642.

237
Spors, H.J.; Stock R.; Sloboda, B. 1992. Räumliche 

Informationssysteme als Entscheidungsbasis für die

forstliche Praxis. Forstarchiv 63: 33-39.

Wever, C. 1999. Laserscanning – ein Verfahren setzt sich 

durch. GIS. Geoinformationssysteme. Zeitschrift für

raumbezogene Informationen und Entscheidungen.

12,2: 12-17.

Wilson, J.P.; Repetto, P.L.; Snyder, R.D. 2000. Effect 

of Data Source, Grid Resolution and Flow-Routing

Method on Computed Topographic Attributes. In: J.P.

Wilson and J.C. Gallant (eds.): Terrain Analyses:

Principles and Applications. John Wilson & Sons, 

New York: 133-161.

Ziegler, M.; Schardt, M.; Konrad, H.; Griess, O. 2000. 

Laserscanning in der Forsteinrichtung. Österreichische

Forstzeitung. 111, 10: 12-13.

238
In: Bevers, Michael; Barrett, Tara M., comps. 2005. Systems Analysis in Forest Resources: Proceedings of the 2003 Symposium. General Technical Report 

PNW-GTR-656. Portland, OR: U.S. Department of Agriculture, Forest Service, Pacific Northwest Research Station.

1 Glen Murphy is a Professor, Department of Forest Engineering, Oregon State University, Corvallis, OR 97331. 2 Hamish Marshall and Barbara Hock are 

Scientists, Ensis Forests, Private Bag 3020, Rotorua, New Zealand.

Systems Analysis in Forest Resources: Proceedings of the 2003 Symposium

MEETING ORDER-BOOK CONSTRAINTS BY ADAPTIVE

CONTROL OF BUCKING ON HARVESTERS

Glen Murphy1, Hamish Marshall2 and Barbara Hock2

INTRODUCTION

Background

There are a number of approaches to optimally bucking

trees. If the market is supply constrained, in other words,

the market will take as much volume of each log type as

each forest owner can produce at the stated market prices,

then the forest owner’s objective should be to maximize the

value of each and every individual tree. This is known as

bucking-to-value. A number of mathematical formulations

and computer models using dynamic programming (DP)

and shortest path (SPN)/longest path (LPN) network algorithms have been developed to optimize the value of individual stems (for example, Pnevmaticos and Mann 1972,

239

Briggs 1980, Geerts and Twaddle 1985, Nasberg 1985,

Sessions and others 1988). 

Unfortunately, from a forest owner’s point of view, 

supply constrained markets are the exception rather than the

rule. Markets are more likely to be demand constrained. They

are also likely to include agreements between customers

and suppliers that relate to specifications but are not reflected in market prices. In addition, the owner is likely to be

faced with operational constraints related to the minimum

economically viable volume of each log-type that can be

produced at a given site in an acceptable timeframe. 

ABSTRACT

Modern mechanized harvesters are often fitted with sensors that measure stem dimensions and with computers that optimally buck each stem based on stem dimensions, qualities, log prices, and desired specifications. Optimally bucking individual

stems, based on market prices, is unlikely to provide yields that meet order book constraints at the harvest unit or forest level. 

An adaptive control heuristic was developed by embedding an individual stem optimal bucking dynamic programming procedure in a threshold accepting algorithm which adjusts relative prices and minimum small end diameter specifications to meet

order book constraints.

The heuristic was tested on four stands where the location and detailed stem description of every tree was known. 

The adaptive control heuristic is described, results from a series of tests on its performance are presented, and the use of

stem information collected during pre-harvest inventory and by the harvester as it works its way through a harvest unit are

compared.

KEYWORDS: Threshold accepting, heuristic, optimal bucking, inventory.
Demand and operational constraints mean that each forest owner is likely to have an order book which contains a

list of customers along with the volumes and proportions

they require of each log-type at specific points in time.

Experience has shown that optimally bucking individual

stems based on market prices is unlikely to yield the products with the attributes and quantities that will meet order

book constraints at the stand or forest level (Sessions and

others 1989). Bucking-to-order therefore requires a different solution procedure than bucking-to-value.

Bucking-to-order formulations usually use a two-level

hierarchical approach with the lower level model optimizing

the bucking of individual stems and the upper level model

meeting constraints through the selection of the best bucking patterns (prices and specifications). Eng and Daellenbach

(1985) and Mendoza and Bare (1986) used a DP algorithm

at the lower level and a linear programming (LP) formulation at the upper level to solve this problem. Nasberg (1985)

and Laroze and Greber (1997) used network algorithms at

the lower level and LP algorithms at the upper level. These

four sets of authors maximized net value at the upper level

and provided one or more bucking patterns that were appropriate to each of hundreds of stem classes per stand (a class

contained stems with similar size and quality features).

Sessions and others (1989) noted, however, that there

are two problems with the practical implementation of the

type of formulation described in the previous paragraph:

(1) identifying which stem class each tree belongs to, and

(2) determining what bucking pattern should be applied to

each tree within the stem class. They proposed a single,

weighted-price set which is applicable to all trees within 

a stand. This approach is also used by Pickens and others

(1997) and Kivenen and Uusitalo (2002) in their model 

formulations.

Sessions and others (1989) used a SPN algorithm at the

lower level and a binary search algorithm at the upper level

to adjust the relative price of a single preferred log type.

Laroze and Greber (1997) used a rule-based algorithm at

the lower level to evaluate bucking patterns and a tabu

search algorithm at the upper level to adjust log specifications (minimum small-end diameters and acceptable qualities) where multiple log types were constrained. Pickens

and others (1997) used a DP algorithm at the lower level

and LP algorithm at the upper level to adjust the relative

prices of multiple log types. The objective at the upper

level for Sessions and others (1989), Laroze and Greber

(1997) and Pickens and others (1997) was to maximize net

value. 

Kivenen and Uusitalo (2002) used a DP algorithm at the

lower level and a fuzzy logic algorithm at the upper level

to adjust the relative prices of multiple log types. Instead of

maximizing net value at the upper level, however, they

maximized the apportionment degree – a measure of how

closely the actual distribution of log-types meets the target

distribution. The apportionment degree concept was originally proposed as a measure of success for optimal bucking

patterns by Bergstrand (1989). Kivinen and others (2003)

compared it with three other “goodness-of-fit” measures

and commented that none of the three alternatives outperformed apportionment degree. 

Sources of stem data and within-stand variability

If perfect information on stem location, stem dimensions,

and qualities in each stem was available on every tree in a

stand it would be theoretically possible, using a two-level

hierarchical approach, to determine the set of relative prices

and specifications that would exactly meet order book constraints at the stand and block (sub-stand felled by a harvester) level prior to harvesting. Perfect information is

rarely available, however. 

Some pre-harvest inventory systems provide detailed

descriptions of a small sample of stems which can be used

as the basis for determining relative prices and specifications prior to the harvester moving into a stand (Deadman

and Goulding 1979). How well the order book constraints

are met at the stand-level depends on how representative

the trees selected for the inventory are. 

Duffner (1980) describes an adaptive price control system which was used in a centralized log processing yard in

Europe in the 1970's to adjust prices based on whether too

many or too few logs of a certain type were produced.

Similar approaches are applicable to harvesting and bucking trees within a stand. 

Once the harvester begins to harvest the stand, more

stem information becomes available. Prices and specifications can then be adjusted based on “better” information.

Using time-of-harvest information for better decision-making is already a common feature on modern harvesters. For

example, Ponsse harvesters store information on the previous 80 stems harvested and use the “closest” eight trees as

the basis for predicting stem taper (M. Laurila, Ponsse Oy.

pers. comm. 2003). 

Natural variability due to site and genetics, along with

variability introduced by man’s uneven treatment of a stand,

mean that even if order book constraints could be met

240
exactly at the stand-level they may not be met at the block

level. Controlling the amount of variation may be of just as

much importance to the forest owner and customer as the

degree to which the order book constraints are met overall.

Malinen and others (2001) used relative root mean squared

error as a measure of variation to evaluate estimations of

sawtimber volume and sawtimber/pulpwood ratio in their

study of harvester and pre-harvest inventory generated

databases.

Both pre-harvest inventory information and time-of-harvest information could be used for adaptively controlling

prices and specifications to better meet order book constraints and minimize variation at the block level.

Study Objective

The objective of our research was to develop and test 

an adaptive control heuristic that uses stem data collected

in the past to predict appropriate prices and log specifications to meet market and operational constraints in the

future.

METHODS

Test Stands

Four stands were used to evaluate the adaptive control

heuristic. Stand details are provided in Table 1. All stands

had been pruned and were of similar mean diameter at

breast height (d.b.h.); about 45 cm. The location of every

tree in each stand was known. One of the stands (WHAKA)

was a real-world Pinus radiata plantation stand in the central North Island of New Zealand. The other three stands

were virtual stands generated to represent a variety of forest

conditions and based on radiata pine tree characteristics. 

The lower limbs had been removed (pruned) from all

trees to a height of approximately 6 m in the EVEN stand.

Selection for pruning was uneven in the UNEVEN stand;

100 percent of trees were pruned in the middle of the stand

decreasing to 70 percent at the edges of the stand. This

mimicked situations where pruning contract supervision or

funds were inadequate to ensure all final crop trees in the

stand were pruned.

The FROST stand mimicked a situation where there was

a frost effect in the center of the stand; tree size was small

in the center and increased towards the edges of the stand.

All trees were pruned in the FROST stand.

Stem Data Sources and Harvesting Pattern

Fifteen circular pre-harvest inventory plots were systematically located in each of the EVEN, UNEVEN, and

FROST stands and five square plots were located in the

WHAKA stand. The inventory plots occupied 3 percent of

total area in each stand.

To evaluate the effects of within stand variability on

adaptive control, the EVEN, UNEVEN, and FROST stands

were broken up into 25 blocks of equal area and the WHAKA

stand was broken up into 10 similar-sized blocks. A serpentine harvesting pattern where the strip width was approximately 20 m wide was used for all four stands. 

We evaluated the effects of using two types of stem data

for determining the best price and log-specification set to

use when harvesting the next block: (1) pre-harvest inventory data and (2) data from the previously harvested block.

Market Requirements for Test Stands

The same market requirements were applied to all four

test stands (Table 2). Five log-types (Pruned Domestic

Sawlog, Unpruned Export Sawlog, Domestic Sawlog #1,

Domestic Sawlog #2, and Pulp) plus waste were included.

All log-types included multiple lengths; some in multiples

of 0.3 m, others in multiples of 0.1 m. A total of 51 lengths

were included in the analyses. Each log-type, and some

sub-groups of log-types, had target proportions of volume

that were required. For example, the order book target proportion for Pruned Domestic Sawlogs was 15 percent of

the total volume harvested.

241

Table 1—Characteristics of Test Stands

EVEN UNEVEN FROST WHAKA

Total Area (ha) 10.0 10.0 10.0 1.9

Density (stems per ha) 375 375 375 249

Mean d.b.h. (cm) 45.0 44.8 45.2 48.5

Total Volume (m3) 5990 5960 6374 1062

Number of blocks 25 25 25 10
In addition to these constraints, we included two others;

a minimum average SED for Domestic Sawlog #1 logs of

250 mm and no less than 75 percent of the combined volume of the 8 m and 12 m Unpruned Export Sawlogs could

be in 12 m lengths. Penalties to the objective function were

applied if these two constraints were not met.

Adaptive Control Heuristic

A two-level hierarchical adaptive control heuristic,

which we called FASTBUCK, was developed and implemented in Microsoft Visual BasicTM 6.0. A DP algorithm

was used to control optimal bucking of single stems at the

lower level based on relative prices and log specifications.

The DP was similar to that described by Deadman and

Goulding (1979) and its aim was to maximize the total

value of each stem. The DP was imbedded in a threshold

accepting (TA) algorithm (Dueck and Scheuer 1990) (fig.

1) which aimed to maximize the apportionment degree

(AD%) from the total number of stems included in the

sample by adjusting SED’s or relative prices3. At the end

of the search procedure, the set of prices and log specifications that yielded the highest AD% at any time throughout

the search was considered to be the solution.

Apportionment degree is defined as:

n

AD % = 100*(1–∑Daj – Drj/ 2)

j=1

where:

n = number of log-types

j = log-type number

Daj = Actual decimal portion of total volume that is logtype j

Drj = Required decimal portion of total volume that was

wanted as log-type j

Penalties were applied to the AD% objective function if

the minimum average SED constraint for Domestic Sawlog

#1 or the minimum proportion constraint for Unpruned

Export Sawlog 12 m were not met. These penalties were

arbitrarily set at 1 percent per mm below the minimum

average SED and 0.5 percent per 1 percent below the target

proportion of Unpruned Export Sawlog 12 m volume.

242

Table 2—Market requirements and constraints for the four test stands

Minimum Relative

small end market Target 

diameters prices proportions 

Log-types Lengths (m) (mm) ($/m3) (percent)

Pruned Domestic Sawlog 3.7-6.1 350 100 15

Unpruned Export Sawlog 12.2m 12.2 260 67 25

Unpruned Export Sawlog 8.2m 8.2 260 61 6

Domestic Sawlog #1 3.7-6.1 200 47 12

Domestic Sawlog #2 Longs 4.9-6.1 200 33 8

Domestic Sawlog #2 Shorts 3.7-4.6 350 32 7

Pulp 3.7-6.1 100 17 24

Waste 0.1-6.0 0 0 8

Additional Constraints:

• Minimum average SED for Domestic Sawlog #1 is 250 mm

• 12.2 m Unpruned Export Sawlogs should be at least 75 percent of the combined total volume for all Unpruned 

Export Sawlogs

3 Some authors only adjust relative price to control product yield (Sessions and others 1989) while others adjust log specifications (Laroze and

Greber 1997). We have noted that logging supervisors sometimes tighten log specifications to reduce product yield if too much of a product is being

produced. We have chosen to include both price and SED adjustments in our formulation. Some limits were put on SED and relative price changes.

SED could not be dropped below the real market limit and price could not be increased above the real market price for Pruned Domestic Sawlogs.
Tuning Search Parameters

Fine tuning of search parameters is itself a significant

combinatorial problem; possible combinations of such parameters as number of threshold steps, threshold step size,

number of iterations per run, number of runs per evaluation,

number of logs-types changed per iteration, etc. can be

very large. Initial tuning was carried out on a stand with

similar characteristics to the EVEN stand. Comments from

the literature on the application of threshold accepting and

other heuristics were also taken into account in selecting

the search parameter values. Based on the tuning of our

heuristic, we carried out all further analyses using:

• five threshold steps, decreasing from 0.050 to 0.001,

• variations in both price and SED, 

• a maximum of two log-types changed per iteration,

• a maximum expansion factor of two for price or SED

increments,

• a single run per evaluation, and

• 1,000 iterations per run (a conservative decision since

500 iterations may have been adequate).

Assessing Within-Stand Variation in AD%

Within-stand variation between blocks was calculated

using the following formula:

243

Figure 1—Adaptive control heuristic

implemented in FASTBUCK.
RMSE% = 100*ADˆ%-1* √ ∑ ˆ (AD%k – AD%k)2 /(m–1)

k=1

where:

RMSE% = relative root mean square error for AD% 

between blocks within a stand

ADˆ% = mean predicted AD% for all blocks within 

a stand

AD% = actual AD% for block k

ADˆ% = predicted AD% for block k

m = number of blocks in the stand

RESULTS

Adaptive Control Heuristic Effectiveness

The adaptive control heuristic was able to provide sets

of relative prices and log specifications that could be used

on a harvester to better meet order book constraints than

would be obtained from using market prices and specifications. An example of the improvement can be seen in

Figure 2 for the pre-harvest inventory data set from the

WHAKA stand. The difference is also evident in the AD%

values for this data-set which were 70.3 percent and 96.2

percent for the market prices and adaptively controlled

prices respectively; an improvement of about 26 percent.

Similar levels of improvement were found for pre-harvest

inventory datasets for the EVEN (~27 percent), UNEVEN

(~22 percent) and FROST stands (~21 percent).

The highest AD% found was 99.6 percent in the FROST

stand for a single block. It must be noted, however, that the

maximum AD% found in a search is partly a function of

the heuristic and partly a function of what is available in

the set of stems being evaluated. For example, if none of

the trees in a stem data set had been pruned, the maximum

AD% that could be found for the market conditions included in our study would have been 85 percent.

Solution times depended on the number of stems in a

block and the number of iterations each block of trees was

optimally bucked per run. Approximately 11 stems per

second were optimally bucked using a Pentium IV processor. Optimally bucking a block of 150 stems 1000 times

(iterations) would, therefore, take about 3.8 hours per run.

Effect of Sources of Stem Data

Table 3 summarizes the results for the four stands harvested using a serpentine pattern. Since detailed descriptions of every stem in the four test stands are known with

perfect knowledge, we can compare the effect of using different sources of stem data in the adaptive control heuristic

on predicted total volume, AD%, and RMSE%.

Total predicted volume, based on the 3 percent inventory, was 5 percent and 19 percent too high for the UNEVEN

and WHAKA stands, the same as actual total volume for

the FROST stand, and 8 percent too low for the EVEN

stand.

When no attempts were made to adaptively control

prices and specifications, in other words, market prices 

and specifications were used to optimally buck all stems 

in each stand, AD% values ranging from 64 to 71 percent

were obtained. These values included penalties for not

meeting the 12 m Unpruned Export Sawlog proportion

constraint for the EVEN, UNEVEN, and FROST stands.

This constraint was met for the WHAKA stand. All stands

met the minimum average SED constraint for Domestic

Sawlog #1.

When stem data from the pre-harvest inventories were

used as the basis for deriving adaptively controlled prices

and log specifications, it was predicted that AD% values

would improve by about 22 percent; In other words, AD%

ranging from 87 to 96 percent depending on the stand.

When these inventory-based, adaptively controlled prices

and log specifications were then applied to the blocks to be

harvested in each of the four stands, all stands had actual

AD% values (AD%_Inv) that were lower than the predicted

AD% values. The EVEN and UNEVEN Stands were about

1 percent lower, the FROST stand was about 5 percent lower,

and the WHAKA stand was about 7 percent lower than

predicted. The variation between blocks (RMSE%) was

relatively low for the EVEN and UNEVEN stands (~3 percent), relatively high for the WHAKA stand (~13 percent),

and very high for the FROST stand (~30 percent).

244

Figure 2—Proportion of total volume by log-type for the WHAKA stand

pre-harvest inventory data-set as obtained through using the adaptive control heuristic (FASTBUCK), as required (Target), and as obtained by optimal bucking based solely on market prices (Market Prices Solution).

m
When stem data from the previously harvested block

were used as the basis for deriving adaptively controlled

prices and log specifications, actual AD% values (AD%_

Block) were better for all four stands than actual AD%_Inv

values, and better than predicted AD%_Inv values for all

but the FROST stand. The RMSE% was slightly higher for

the EVEN and UNEVEN stands, and substantially lower

for the FROST and WHAKA stands than was found for the

inventory-derived prices and log specifications. 

DISCUSSION AND CONCLUSIONS

The adaptive control heuristic described in this paper

was able to provide a single list of relative prices and log

specifications which could be applied to all stems within 

a stand, or sub-stand, and which could better meet order

book constraints than when market prices and log specifications were used to optimally buck individual stems. The

heuristic used AD%, with penalties for not meeting some

market constraints, as its performance measure for success.

We were able to find solutions with AD% values of over

99.5 percent for some blocks of trees within a stand but this

was partly a function of the heuristic and partly a function

of what was available in the set of stems being evaluated.

The heuristic starts by using market prices and log 

specifications to optimally buck the stems and evaluate the

AD% for that solution. Alternative start points could be to

randomly allocate relative prices to each log-type or to set

all log-types to the same relative price. Siedentopf (1995),

however, found no advantage in selecting a random start

point for job shop scheduling problems solved with a threshold accepting algorithm. We believe there is some advantage in watching the adaptive control heuristic move from

the unconstrained solution, based on market prices and log

specifications, towards those which meet order book constraints. 

Some heuristics explore a systematic grid around the

current solution and select the move which maximizes the

objective function. With five log-types and two variables

(price and SED) a full exploration of the closest neighbors

on the systematic grid would require over one thousand

lots of optimal bucking of stems per iteration at the lower

level using the DP. Since most of the solution time for the

heuristic is taken up by the DP bucking of hundreds of

stems per iteration, a full search of the systematic grid is

not feasible in practical terms. Bettinger and others (2002)

comment that 2-opt moves, where two variables were

changed at the same time, were better than 1-opt moves

when using tabu search heuristics with harvest scheduling

problems. We chose to use multiple changes in a single

move per iteration. Our limited tests showed that multiple

changes were better than a single change per iteration; in

some cases leading to a 7% improvement in AD%.

245

Table 3—Effect of different sources of data on how well the adaptive control heuristic meets order book constraints

at the stand and block level based on apportionment degree (AD%) and relative root mean squared error (RMSE%

- shown in brackets) performance measures.

Stand Name

Performance Measures EVEN UNEVEN FROST WHAKA

Predicted total volume based on 3 percent pre-harvest inventory (m3) 5,533 6,291 6,359 1,262

Actual total volume (m3) 5,990 5,970 6,374 1,063

AD% resulting from optimally bucking all stems using market 

prices and log specifications 68.8 63.9 69.7 71.0

Predicted AD% based only on the 3 percent Pre-harvest inventory, 

adaptively controlled prices, and log specifications 87.1 86.8 91.6 96.2

Actual AD% and RMSE% based on bucking all stems using 

adaptively controlled prices and log specifications derived from 

the 3 percent pre-harvest inventory (AD%_Inv) 86.5 86.0 87.1 89.4

[2.4] [3.4] [29.5] [13.2]

Actual AD% and RMSE% based on bucking all stems using 

adaptively controlled prices and log specifications derived from

the previous Block’s stems data (AD%_Block) 87.7 87.6 90.9 97.0

[4.4] [3.7] [12.9] [9.8]
Achieving 100 percent for the AD% from a stand will

not ensure that all market constraints are met, particularly

where specific quantities, as opposed to relative quantities,

are required. For example, if the total volume in a stand

was only 5,000 m3 it would be impossible to meet order

book constraints for log-types totaling 6,000 m3 even if the

log-types were available in the correct proportions. Despite

these problems, there is need for improved adaptive control

heuristics that can be used on harvesters to help meet volume as well as proportion requirements.

The evaluation of different sources of stem data showed

that predicted AD% values, based on pre-harvest inventory

information, could be obtained at the stand level provided

that relative prices and log specifications were adaptively

controlled on a block by block basis. It was found, however, that pre-harvest inventory stem data were not as good 

as stem data from the previously harvested block if the

pre-harvest inventory data were used to fine-tune prices

and specifications prior to harvest, and no further changes

were made throughout the harvest of the stand. This result

is similar to the findings of Kivinen and Uusitalo (2002)

who commented that “no gain can be achieved by pre-harvest inventorying of stands and fine-tuning the price lists

with this imperfect knowledge”. Stand inventories, of course,

have a role to play in tactical and strategic management of

forest resources. Using stem data from the previously harvested block to adjust prices and specifications for the current block gave the highest AD% values for all of the test

stands. 

Waiting for the whole stand to be completed before

finally meeting order book constraints is not an option for

many wood suppliers. Managing within stand variation is,

therefore, of considerable interest. We found that stand

conditions have a large effect on the between block variation in AD%, as measured by RMSE%. The FROST and

WHAKA stands had substantially larger RMSE% than the

EVEN and UNEVEN stands. 

Using a sub-sample of the trees harvested in the previous block, or from even earlier harvested blocks, may allow

for faster solution times if the selected sub-sample closely

reflects the trees in the current block to be harvested. As

noted earlier in this paper, Ponsse harvesters use a system

where the eight most similar trees to the current tree being

harvested are used to predict stem taper. Malinen and others (2001) provide evidence that such a system has application at the forest level – between stands. The challenge at

the sub-stand level will be to find parameters which reflect

both size and quality attributes of trees so that the trees with

the “most similar” characteristics in the harvester stem

database can be easily found.

In a report on five merchandising computers that adaptively control bucking on a tree-by-tree basis, effectively a

block size of one tree, Sondell and others (2002) noted that

AD% rose quickly to about 85 percent for three computers,

more slowly to about 80 percent for one computer, and to a

maximum of about 60 percent for the fifth computer. The

five computers were tested in a Grade 3 spruce stand in

Sweden

FURTHER WORK

While the results from this study have shown that the

application of an adaptive control procedure would better

meet order book constraints than solely using market prices

and log specifications, much work still remains to be done.

Further research in the area of adaptive control would seek

to:

• evaluate the adaptive control heuristic presented in 

this paper in a wider range of stands,

• investigate new objective functions other than maximizing AD% and evaluate the effects of different constraint penalty functions (for example, exponential

functions to heavily penalize solutions which are farremoved from the constraints),

• determine how many stems from the previously 

harvested block are required to obtain acceptable

AD% values; how small a block can be used,

• determine how short the search-length can be and still 

yield acceptable AD% values,

• investigate which data sources are important; ones that 

are spatially related, ones that are temporally related,

or ones that are most similar,

• evaluate the effects of different weighting factors on 

data sources, and

• determine the impact of harvester work methods and 

measurement accuracy on AD% and RMSE%.

LITERATURE CITED

Bergstrand, K-G. 1989. Fordelningsaptering med naroptimal. Skogsarbeten. Unpublished internal report. 11 p.

Bettinger, P.; Graetz, D.; Boston, K.; Sessions, J.; 

Chung, W. 2002. Eight heuristic planning techniques

applied to three increasingly difficult wildlife planning

problems. Silva Fennica 36(2): 561-584.

246
Briggs, D.G. 1980. A dynamic programming approach to 

optimizing stem conversion. Ph.D. thesis. Univ. of

Washington, Seattle. 393 p.

Deadman, M.W.; Goulding, C.J. 1979. A method for 

assessment of recoverable volume by log-type. N.Z. J.

For. Sci. 9: 225-239.

Dueck, G.; Scheuer, T. 1990. Threshold accepting: a general

purpose optimization algorithm appearing superior to

simulated annealing. J. Computational Physics 90: 

161-175.

Duffner, W.W. 1980. Decision making from market to 

stump. P.81-95 in Proceedings of Weyerhaeuser Science

Symposium, Tacoma, Washington, USA.

Eng, G.; Daellenbach, H.G. 1985. Forest outturn optimization by Dantzig-Wolfe decomposition and dynamic

programming column generation. Op. Res. 33: 459-64.

Geerts, J.M.P.; Twaddle, A.A. 1985. A method to assess 

log value loss caused by cross-cutting practice on the

skidsite. N.Z. J. For. 29(2): 173-184.

Kivinen, V.-P.; Uusitalo, J. 2002. Applying fuzzy logic to 

tree bucking control. For. Sci. 48(4): 673-684.

Kivinen, V.-P.; Uusitalo, J.; Nummi, T. 2003. Measuring 

the goodness of the bucking outcome. P. 87-95 in

Proceedings of the 2nd Forest Engineering Conference:

Raw Materials Utilization (ed. M. Iwarsson Wide and

B. Baryd). Arbetsrapport Nr. 537. Skogforsk, Sweden. 

Laroze, A.J.; Greber, B.J. 1997. Using tabu search to 

generate stand-level, rule-based bucking patterns. For.

Sci. 43(2): 157-169. 

Malinen, J.; Maltamo, M.; Harstela, P. 2001. 

Applications of most similar neighbor inference for 

estimating marked stand characteristics using harvester

and inventory generated stem databases. Int. J. For.

Eng. 12(2): 33-41.

Mendoza, G.A.; Bare, B.B. 1986. A two-stage decision 

model for log bucking and allocation. For. Prod. J.

36(10): 70-74.

Nasberg, M. 1985. Mathematical programming models for 

optimal log bucking. Linkoping Studies in Science and

Technology. Diss. No. 132. Department of Mathematics,

Linkoping University, Sweden. 200 p.

Pickens, J.B.; Throop, S.A.; Frendewey, J.O. 1997.

Choosing prices to optimally buck hardwood logs with

multiple log-length demand restrictions. For. Sci. 43(3):

403-413.

Pnevmaticos, S.M.; Mann, S.H. 1972. Dynamic programming in tree bucking. For. Prod. J. 22(2): 26-30.

Sessions, J.; Layton, R.; Guangda, L. 1988. Improving 

tree bucking decisions: a network approach. The

Compiler 6(1): 5-9.

Sessions, J.; Olsen, E.; Garland, J. 1989. Tree bucking 

for optimal stand value with log allocation constraints.

For. Sci. 35(1): 271-276.

Siedentopf, J. 1995. An efficient scheduling algorithm 

based upon threshold accepting. Institute of Production

Management and Industrial Information Management,

University of Liepzig, Germany. Report No. 16. 22 p.

Sondell, J.; Moller, J.J.; Arlinger, J. 2002. Third-generation merchandising computers. Skogforsk Results No.

2. 6 p.

247
This page is intentionally left blank.
In: Bevers, Michael; Barrett, Tara M., comps. 2005. Systems Analysis in Forest Resources: Proceedings of the 2003 Symposium. General Technical Report 

PNW-GTR-656. Portland, OR: U.S. Department of Agriculture, Forest Service, Pacific Northwest Research Station.

1 M.H. Pelkki is an Associate Professor at the Arkansas Forest Resources Center and the University of Arkansas-Monticello School of Forest Resources, 

Monticello, AR 71656.

Systems Analysis in Forest Resources: Proceedings of the 2003 Symposium

EFFECTS OF STOCHASTIC STUMPAGE PRICES 

ON OPTIMAL TIMBER ROTATIONS DETERMINED 

BY DYNAMIC PROGRAMMING

Matthew H. Pelkki1

INTRODUCTION

Forward recursive dynamic programming is a longaccepted method in stand-level optimization research to

solve optimal rotation-length and density problems (Amidon

and Akin 1968, Schreuder 1971). During the last thirty years,

advances in the application of dynamic programming have

included the use of state neighborhoods to permit the use

of continuous variables in discrete-state dynamic programming formulations (Brodie and Kao 1979, Pelkki 1997b),

the use of individual-tree growth and yield models to permit optimization under various stand-level constraints or

multiple objectives (Riiters et al. 1982, Kirillova and Pelkki

2003), heuristics (Arthaud and Pelkki 1996), and regionlimiting methods (Paredes and Brodie 1987). These and

other studies (Arthaud and Klemperer 1988, Valsta 1990)

have found optimal results using deterministic prices in

their return or objective functions. Stochastic price functions can be applied to both backward- and forward-recursive dynamic programming formulations.

ABSTRACT

A three-state forward-recursive dynamic programming simulation was used to investigate the effects of stochastic

stumpage prices on optimal rotation length, thinning strategies, timing, and intensity of intermediate harvests in shortleaf

pine plantations in the south-central United States. Historical timber stumpage price data shows that 2-year changes in

stumpage prices can range from -49% to +96%. Under stochastic stumpage prices, optimal rotation length and soil expectation values varied greatly. Rotation lengths ranged from 40 to 118 years and soil expectation values ranged from $173/ac 

to $2818/ac. Thinning from above was the strategy chosen in 96% of all observed intermediate harvest, with typical entries

occurring late in the second, third, and fourth decades of the plantation, with progressive thinnings removing an increasing

percentage of basal area. Perfect knowledge of future stumpage prices was found to increase land expectation value by an

average of 48%, indicating the possible returns from good “market timing” of harvests.

249

Stochastic price functions will simulate realistic price

changes according to the parameters used in the underlying

frequency distribution. These parameters reflect our assumptions about future market conditions, particularly about the

rate of future price increases or decreases, and the amount

of variation or uncertainty in future prices over short periods of time. Both overall price trends and price stability are

important factors when making intermediate or final harvest decisions. This study reports the impacts of stochastic

prices on aspects of optimal financial timber rotations (rotation length, timing, intensity, and method of thinnings) for

shortleaf pine (Pinus echinata, L.) plantations in the southcentral United States. The optimal timber rotations are

determined through the use of repeated simulations where

forward-recursive dynamic programming is used to determine the optimal solution.

METHODS

A forward-recursive dynamic programming method was

used in this study. The objective function can be defined

mathematically in equation 1:
N

f N (YN) = ∑rn(Tn) (1) n=0

Where the variables and functions can be defined as:

Tn – a management action with resulting physical outputs at stage n

rn – is the net present value of the action Tn

YN – is the ending state for the problem, in this case,

the final harvest is a clearfelling operation.

The objective function (1) is subject to several constraints, the first of which links states in a particular stage

(Yn) to states in a subsequent stage (Yn+1) through the

growth and yield function (Gn+1) and the applied management action (Tn). This relationship is shown in equation 2:

Yn + Gn+1(Yn) – Tn = Yn+1 for n = 0, 1, 2, 3, .... N-1 (2)

Two additional bookkeeping constraints, equations three

and four, define a state in a stage prior to any management

activity (Xn) and also define that the final action (TN) is a

clearfelling harvest. 

Xn – Tn = Yn (3)

XN – TN = YN = 0 (4)

Finally, a recursive relationship explains how the policies (management actions) in each stage are linked together

sequentially to form an optimal policy or set of management actions over an entire rotation:

f n(Yn) = max [r(Xn , Tn ) + fn–1(Yn–1)] (5) (Yn-1,

Tn )

Implementing this dynamic programming formulation

in a computer simulation requires several elements. First,

state and stage variables must be defined. For a forwardrecursive formulation, stand age is the most straightforward

choice for the stage variable. In this study, the stage (growth)

interval was set at two years. Many possibilities exist for

choosing state variables. Descriptors for density include basal

area, number of trees, total volume, and average diameter.

Other state variables could include number of previous

thinnings, previous thinning strategy or type, species composition, or average tree quality. Previous research (Pelkki

1997b) has shown that many combinations of state variables

can lead to an optimal solution as long as the state neighborhood is small enough that any path entering the state has

the same future growth and value potential. In this way, the

principle of optimality (Dykstra 1984) is not violated. In

this study, number of trees per acre (TPA) and cubic foot

volume per acre (ft.3/ac) were chosen as state variables.

The neighborhoods or intervals for these variables were set

at + 10 TPA and + 10 ft.3/ac. A third state descriptor, number of previous thinnings (NTh) was added to the formulation. This prevented state with different numbers of thinnings

to be compared to each other and also permitted the total

number of thinnings in any policy set for a rotation to be

restricted to three or less.

The initial stage and state condition was a 10-year old

shortleaf pine plantation tree list taken from Smalley and

Bailey (1974) for a plantation of 500 trees per acre with a

site index of 82 feet at 50 years. While multiple starting

conditions could have been considered, the purpose of this

research was to study the impacts of stochastic price functions on dynamic programming results, and previous

research (Pelkki 1997a) had indicated that this planting

density was optimal. Stand establishment costs were $131

per acre and based on South wide average costs for planting and site preparation (Dubois et al. 2003).

An individual-tree growth model based on the TWIGS

growth and yield system (Miner et al. 1987) was used to

project stand states forward in time, representing the function Gn+1() from equation 2. The model used a parameter

set fit for the Southeast United States.

Various thinning regimes were programmed to represent

possible management actions (Tn) imposed upon unthinned

but projected states (Xn) in any stage. These actions included thinning from below (TfB), thinning from above (TfA),

thinning from above and below (TAB), and mechanical

row thinning (TM). These four thinning strategies could be

applied at levels removing between 10% and 50% of the

initial state’s (Tn–1) basal area, in 10% increments. Finally,

management actions of clearcutting (CC) and “do nothing”

(DN) were also applied to each initial management state

(Tn).

The return function, rn() was based on a net present

value equation applied to the physical outputs from each

action Tn applied to initial state Xn. Thus, every completed

state Yn included a residual tree list, a harvested tree list,

state parameters, and present value for action Tn, and

cumulative present value for any previous thinning actions

in the policy leading from the starting state in stage 0 to

that particular state. Volume equations from TWIGS

(Miner et al. 1987) converted the individual trees to sawtimber and pulpwood volumes, to which stumpage prices

were applied and then a present value was calculated as

shown in equation 6.

250
Pn × V(Tn ) – F – PCT rn(Tn ) = (6) (1 + i)nw

The variable w represents the width of the stage interval, which in this case was 2 years, and the variable i represents the cost of capital or the hurdle rate, which for all

simulations was set at a nominal rate of 8%. The variable F

represents a fixed entry cost of $40 per acre on all harvests,

and PCT represents the cost of removing premerchantable

stems. Furthermore, the value Pn is reduced by 15% in all

thinning operations to reflect the reduced efficiency and

higher costs of partial removals.

The values Pn were derived from a stochastic price

function. Historic two-year stumpage price changes were

determined from timber prices reported by Timber MartSouth (Timber Mart-South 1981-2001) and from the State

of Louisiana (LDA&F, 2002). Table 1 shows the minimum,

maximum and mean 2-year nominal price changes observed

in both data sets for pine sawtimber and pulpwood stumpage

prices. Due to the length of time available for the Louisiana

dataset, it was selected to model the frequency distribution

of prices for this study with a skewed triangular distribution based on minimum, mean, and maximum two-year

price changes.

It is worthwhile to note that timber stumpage prices 

are typically correlated over time (Haight and Holmes 1991,

Yin and Newman 1995) and can be modeled with time

series analyses. The purpose of this study is not to predict

future prices but rather to observe optimal stand rotations

under conditions of stochastic prices based on a frequency

distribution that is based on historical data patterns.

As mentioned previously, the dynamic programming

formulation was constrained to three thinnings. This was

necessary as initial, unrestrained simulations selected thinning every time there was a rise in prices. Reducing the

number of allowable intermediate harvests to a reasonable

number (3) solved this problem.

With the dynamic programming model in place, thirty

stochastic runs were completed. Analysis of each run

included the economic returns under a particular price projection, the rotation length, thinning strategies, intensities,

and timing of thinnings used.

RESULTS AND DISCUSSION

General results for the 30 dynamic programming simulations are provided in table 2. An average rotation age of

71 years resulted in a soil expectation value (SEV) of

$966.50 per acre. Table 3 shows the yields for a “typical”

70-year rotation with the three thinnings from above that

were most frequently observed in the 30 simulations. Typical

rotation lengths for shortleaf pine in the central United States

range from 40 to 60 years under conditions of price uncertainly or under the assumption of constant prices. The rotation lengths and economic returns may appear high for

shortleaf pine, but the reader should remember that the stochastic price model had an average undiscounted 2-year

price nominal increase of 14% for sawtimber and 9% for

pulpwood stumpage. The two-year discount rate was 16.64%

(nominal annual rate of 8%), which accounts for the relatively high values. Table 4 provides some statistics relative

to the prices generated in the simulations. Based on the triangular distribution parameters presented in table 1 from

Louisiana, we can see that over a one hundred year time

frame that sawtimber stumpages prices increase to an average of $116 per board foot. While this seems incredible, it

is in fact an average increase of 6.4% per year or 13.3%

every two years, which is less than the average historical

251

Table 1—Two-year nominal changes in southern pine

stumpage prices.

Source: Timber Mart-South (1981 to 2001)

Product Minimum Maximum Mean

Pine sawtimber -32% +96% 9%

Pine pulpwood -51% +73% 5%

Source: Louisiana Quarterly Report of Timber

Products (1955 to 2001)1

Product Minimum Maximum Mean

Pine sawtimber -43% +84% 14%

Pine pulpwood -49% +51% 9%

1 Frequency distribution for this study based on parameters from LA

dataset.

Table 2—Summary of observations on soil expectation

value and rotation age from a shortleaf pine plantation

under repeated simulations using a stochastic price

function.

Soil Expectation Value ($/ac) Rotation Age

Mean $ 966.5 71

Min $ 173.3 40

Max $ 2818.5 118

SE $ 113.6 3

N = 30
price increases in Louisiana from 1955 to 2001. Rotation

lengths and values for bare land would undoubtedly be

lower under less optimistic future price assumptions.

Furthermore, the forward-recursive nature of dynamic

programming allowed the simulation to have “perfect

information” regarding price changes. Since a “do-nothing”

alternative was simulated from all states, the solution network preserved states that had not been harvested at every

stage, permitting the optimal policy to choose the timing 

of all harvests with perfect knowledge of prices throughout

a rotation. While this is an unrealistic assumption, a more

realistic calculation of economic returns will be discussed

below, as well as an estimate of the value of “perfect market information.”

Of the four thinning strategies permitted, thinning from

above (TfA) was observed in 96% of the cases (n = 86/90).

Table 5 shows frequency and median stand age and intensity

for all four thinning strategies. All simulations used three

thinnings, with intensities generally increasing in later thinning operations. Thinning from above and below (TAB)

was infrequently observed (n = 4/90), and with so few

observations, generalizations cannot be made. Mechanical

or row thinning and thinning from below were never

observed in the 30 simulations.

Because the dynamic programming formulation had

perfect knowledge of the future, it is understandable that

65% (n=58) of the thinning actions came at the end of a

run of stumpage price increases. However, a considerable

number of thinnings were observed either in the middle

(30%, n = 27) or during a decline (6%, n = 5) in stumpage

prices. These operations were all first or second thinnings,

indicating that that the timing of final thinnings and harvests to high market prices is more important and that midrotation thinnings play an important role in maintaining

growth rates.

In order to calculate a more realistic SEV for the results

of this study, the median observed times for thinning and

252

Table 3—Simulated per acre yields from the most frequently observed thinning schedule and rotation length.

Stage Stand Age Harvest activity Sawtimber yield (bf/ac) Pulpwood yield (ft.3/ac)

0 10 Start — —

9 28 Thin from above 20% of basal area 2,370 —

13 36 Thin from above, 50% of basal area 5,009 215

19 48 Thin from above, 50% of basal area 4,481 —

30 70 Final harvest 6,362 34

Total per acre production over rotation 18,222 249

Table 4—Simulation-based pulpwood and sawtimber price statistics for selected stages.

Rotation Pulpwood price $/CCF Sawtimber price $/MBF

Stage Age Mean Min Max St. Dev. Mean Min Max St. Dev.

0 10 $24 $24 $24 $0 $227 $227 $227 $0

5 20 $28 $10 $47 $9 $404 $198 $868 $175

10 30 $32 $8 $60 $14 $680 $213 $1380 $321

15 40 $36 $5 $81 $19 $1074 $242 $2469 $554

20 50 $44 $4 $172 $31 $1687 $263 $3062 $857

25 60 $55 $6 $359 $62 $2525 $238 $7496 $1720

30 70 $72 $5 $434 $90 $5601 $374 $22361 $5656

35 80 $89 $5 $594 $121 $11175 $679 $56807 $14162

40 90 $107 $13 $687 $153 $24421 $2119 $66429 $21461

45 100 $96 $23 $403 $112 $71514 $2883 $459913 $111,549

50 110 $75 $25 $216 $67 $116133 $23760 $260646 $73662
final harvest, as reported in table 4, were resimulated under

the recorded prices for each of the 30 simulations. By

adhering to a fixed schedule of thinnings and a pre-determined rotation age, the SEV was reduced 10% to 65%,

with an average SEV of $641 per acre, a 34% reduction

from a situation of “perfect knowledge.” Looking at this

from the basis of imperfect knowledge, soil expectation

values could increase from 12% to 287% with perfect

knowledge of future prices, with an average increase of

48%.

The practical application of these results depends on

your assumptions of the behavior of stumpage markets into

the future. In 1955, sawtimber stumpage sold for $31/MBF

in Louisiana, and in 1997, the selling price was $412/MBF.

Had landowners known this would occur (perfect knowledge), then selling 35-year old pine stands in 1990, when

stumpage prices were $184 per MBF would have been irrational. Perhaps the most critical decision for forest landowners is their “belief” in the direction of future stumpage

prices. If future softwood stumpage prices are as much as

$5.60 per board foot in the year 2064 (as suggested by

table 3), rotation lengths of 70 years will be economically

optimal. Larger questions dealing with future supply and

demand interactions, building material substitutions, and

consumer tastes and preferences will certainly effect the

rise in future stumpage prices. The U.S. Forest Service estimates that timber stumpage markets will be much weaker

through the year 2050, with prices increasing but at a much

slower rate than from the period 1952-2000 (Haynes, 2003).

If this is indeed the situation, then optimal rotation lengths

will shorten and thinning will probably be from below or

mechanical to hasten the production of sawtimber. For 

forest managers and landowners, thinning from above can

prolong rotation lengths and preserve future harvesting

options under situations where the landowner expects higher

future prices.

Future research will address some of the limitations of

this work. The number of observations using this type of

empirical approach is rather small. Increasing the number

of sample runs to 100 or 200 runs will allow for a more

thorough analysis of conditions under which thinning

strategies other than thinning from above are employed.

Altering the parameters of the triangular stochastic price

function, and examining sensitivity to site quality and cost

of capital should also yield useful results. Also, testing sensitivity to constraints on the number of intermediate harvests

and thinning strategies permitted under stochastic prices

might elucidate more information about system behavior

and optimal harvesting patterns.

LITERATURE CITED

Amidon, E.L., and G.S. Akin. 1968. Dynamic programming to determine optimal levels of growing stock. For.

Sci. 14: 287-291.

Arthaud, G.J., and W.D. Klemperer. 1988. Optimizing 

high and low thinnings in loblolly pine with dynamic

programming. Can. J. For. Res. 18: 1118-1122.

Arthaud, G.J. and M.H. Pelkki. 1996. A Comparison of 

Dynamic Programming and A* in Optimal Forest Stand

Management. Forest Science, Vol. 42(4): 498-503.

Brodie, J.D. and C. Kao. 1979. Optimizing thinning in 

Douglas-fir with three descriptor dynamic programming

to account for accelerated diameter growth. For. Sci. 25:

665-672.

Dubois, M.R., T.J. Straka, S.D. Crim, and L.J. 

Robinson. 2003. Cost and Cost Trends for Forestry

Practices in the South. Forest Landowner 62(2): 3-9.

253

Table 5—Observed thinning strategy, timing, and intensity for a shortleaf pine plantation under repeated simulations 

using stochastic prices.

Median Median Median

age age age

Thinning Frequency First Second Third median median median

Strategy of use thinning thinning thinning intensity intensity intensity

Thinning from above 96% 28 20% 36 50% 48 50%

Thinning from above

and below 4% 26 20% 30 30% — —

Mechanical thinning 0% — — — — — —

Thinning from below 0% — — — — — —
Dykstra, D.P. 1984. Mathematical programming for 

natural resource management. McGraw-Hill. 318 p.

Haight, R.G. and T.P. Holmes. 1991. Harvesting loblolly 

pine plantations with hardwood competition and stochastic prices. For. Sci. 37: 1266-1282.

Haynes, R.W. 2003. An Analysis of the Timber Situation 

in the United States, 1952-2050. U.S.D.A. Forest

Service, Pacific Northwest Research Station, Portland,

OR. PNW-GTR-560. 254 p.

Kirillova, N.V., and M.H. Pelkki. 2003. Optimal 

Financial Management of the Yellow-Poplar/Oak

Covertype in the Central Appalachians. In G.J. Arthaud

and T.M. Barrett (eds.), Spatial Analysis in Forest

Resources, proceedings of the Eighth Symposium, Sept.

27-30, 2000. Snowmass Village, Colorado, USA.

Kluwer Academic Publishers, Dordrecht, Netherlands.

Pp. 261-268.

LDA&F. 2002. Louisiana Quarterly Report of Forest 

Products. Louisiana Department of Agriculture and

Forestry. Data for 1955 to 2001 accessed on web page

www.ldaf.state.la.us/divisions/forestry/reports/quarterreport/original-report.asp on August 20, 2003.

Miner, C.L., N.R. Walters, M.L. Belli. 1988. A guide to 

the TWIGS program for the North Central United

States. USDA For. Serv. Gen. Tech. Rep. NC-125. 

105 p.

Paredes, G.L. and J.D. Brodie. 1987. Efficient 

Specification and Solution of the Even-aged Rotation

and Thinning Problem. For. Sci., Vol. 33, No. 1, pp. 

14-29.

Pelkki, M.H. 1997a. Optimal shortleaf pine management 

in the Central United States. Northern Journal of

Applied Forestry 14(2)67-71.

Pelkki, M. H. 1997b. The Effects of Neighborhood 

Storage Size on Dynamic Programming Solutions.

Forest Science 43(3): 387-395.

Riitters, K., J.D. Brodie, and D.W. Hann. 1982. Dynamic 

programming for optimization of timber production and

grazing in ponderosa pine. For. Sci. 28: 517-526.

Schreuder, G.F. 1971. The simultaneous determination of 

optimal thinning schedule and rotation for an even-aged

forest. For. Sci. 17: 333-338.

Smalley, G.W. and R.L. Bailey. 1974. Yield tables and 

stand structure for shortleaf pine plantations in

Tennessee, Alabama, and Georgia highlands. USDA

For. Serv. Res. Pap. SO-97.

Timber Mart-South. 1981- 2001. Annual stumpage price 

summaries from fourth quarter reports for all southern

states in years 1981 to 2001. Timber Mart-South, D.B.

Warnell School of Forest Resources, The University of

Georgia, Athens, GA 30602.

Valsta, L.T. 1990. A comparison of numerical methods for 

optimizing even aged stand management. Can. J. For.

Res. 20: 961-969.

Yin, R.S. and D.H. Newman. 1995. A note on the treecutting problem in a stochastic environment. J. For.

Econ. 1: 181-190.

254
In: Bevers, Michael; Barrett, Tara M., comps. 2005. Systems Analysis in Forest Resources: Proceedings of the 2003 Symposium. General Technical Report 

PNW-GTR-656. Portland, OR: U.S. Department of Agriculture, Forest Service, Pacific Northwest Research Station.

1 D. Graetz is a Graduate Research Assistant, Department of Forest Resources, Oregon State University, Corvallis, OR 97331; 2 P. Bettinger is an Associate 

Professor, Warnell School of Forest Resources, University of Georgia, Athens, GA 30602.

Systems Analysis in Forest Resources: Proceedings of the 2003 Symposium

DETERMINING THINNING REGIMES TO REACH STAND 

DENSITY TARGETS FOR ANY-AGED STAND MANAGEMENT

IN THE BLUE MOUNTAINS OF EASTERN OREGON

David Graetz1 and Pete Bettinger2

INTRODUCTION

Stand-level management planning assumes that each

stand, in isolation of other stands in a forest, is to be managed in the best manner possible. Decisions regarding the

management of forested stands have traditionally been based

on expert knowledge and field experience with a number of

different treatment regimes (Palahí and Pukkala 2003). The

optimum combination of the number of entries, the timing

of entries, and the type of treatments applied to a stand is

difficult to determine through these non-computerized processes, since generally only a small portion of the potential

combinations of treatments have been evaluated. Operations

research techniques applied to the stand-level management

started appearing in the literature in the 1960s (Chappelle

and Nelson 1964; Hool 1966; Amidon and Akin 1968).

Since that time a plethora of research has been done to

solve problems for even-aged, uneven-aged, and any-aged

stand management using a variety of optimization methods.

ABSTRACT

A stand-level optimization procedure, using the RLS-PATH algorithm, was developed to allow the development of management prescriptions that address ecological or biological goals for eastern Oregon forests. The procedure requires that

stand-level goals be defined by basic tree-level data (e.g., basal area, trees per acre, species, etc.). It then seeks to minimize

the deviations of stand conditions relative to these goals over a 100-year time frame. A distance-independent, single-tree

growth model is embedded within a user-friendly program called SLOMO that also executes and solves the stand-level 

optimization problem. A variety of constraints can be turned on and off and thus can be used to determine surrogate shadow

prices in meeting those constraints. An example stand from the Upper Grande Ronde Watershed will be used to demonstrate

the capabilities of SLOMO.

255

Stand-level optimization methods have evolved with 

the changing demands placed on forests. Initially the goals

were to maximize economic or commodity production values, yet more recently have placed emphasis on non-commodity values. The approaches that can be used to develop

optimal stand-level management prescriptions include the

Hooke and Jeeves (1961) method (Roise 1986a, 1986b;

Haight and Monserud 1990a, 1990b; Yoshimoto and others

1990; Valsta 1990; Haight and others 1992; Thorsen and

Helles 1998), other non-linear programming methods

(Adams and Ek 1974; Kao and Brodie 1980; Roise 1986a;

Haight and others 1992; Bare and Opalach 1987; Gove and

Fairweather 1992 ), dynamic programming (Hool 1966;

Amidon and Akin 1968; Brodie and others 1978; Brodie

and Kao 1979; Kao and Brodie 1979; Chen and others

1980a, 1980b; Haight and others 1985a; Arthaud and

Klemperer 1988; Yoshimoto and others 1990; Brukas and

Brodie 1999), and other heuristics or amalgamations of

various methods (Buongiorno and Michie 1980; Bullard

and others 1985; Valsta 1990; Gong 1992; Wikstrom and

Eriksson 2000; Wikstrom 2001). 
Most stand-level optimization methodology reported in

the literature focuses on meeting forest economic or commodity production goals and not so much on biological or

ecological goals. There are, however, some exceptions.

Haight and others (1985a) tracked biological indicators in

the development of stand prescriptions, although they were

not influential in developing the management prescriptions.

It wasn’t until the early 1990s when literature appears showing stand-level optimization problems that incorporate 

non-economic or non-commodity production goals and

constraints. Haight and others (1992), for example, incorporated stand density targets into the development of optimal prescriptions, using penalty functions to ensure the

attainment of goals. Buongiorno and others (1994) measured the cost of maintaining tree size diversity by solving

the stand-level problem with a measure of diversity in and

out of the objective function as well as in and out of the

constraints. The Applegate Project tackled the any-aged

stand management problem by developing a stand-alone

prescription generator that solved for multiple biological

and ecological goals in the objective function (Wedin 1999).

More recently, Wikstrom and Eriksson (2000) incorporated

biological measures of green tree retention, broad-leaved

tree importance, coarse woody debris, and tree-size diversity

into the objective function of an even-aged stand management problem.

Dynamic programming (DP) approaches have the ability to avoid local optima in the final solution to a standlevel problem (Yoshimoto and others 1990). As well, DP

has characteristics and capabilities that provide a useful

solution framework for the timing and intensity of typical

silvicultural actions such as thinning, fertilization, pest control and rotation age decisions (Hann and Brodie 1980;

Paredes and Brodie 1987). However, when DP is employed,

the number of computations required, as well as the storage

space required, increases exponentially with the number of

state variables assumed (Hann and Brodie 1980; Martin and

Ek 1981; Haight and others 1985; Yoshimoto and others

1990). One disadvantage of DP is the absence of shadow

prices that incorporate explicit conditions of first-order

optimality, data that might be provided by other optimization processes (Kao and Brodie 1980). Brodie and others

(1978) indicate that the forward procedure associated with

DP is more amenable to developing prescriptions that

include thinning than the backward procedure.

While DP has been applied in conjunction with a number of growth and yield models, the thinning approaches

have usually consisted of thinning from above, from below,

or even with the distribution of diameters. A relaxation of

these controls in stand-level optimization has been achieved

with non-linear programming. Non-linear programming

approaches attempt to solve problems with non-concave

production surfaces and have trouble locating the optimal

solutions to problems (Yoshimoto and others 1990). The

PATH (Projection Alternative TecHnique) algorithm (Paredes

and Brodie 1987) was developed to solve for the optimal

solution at every stage in a DP program, reducing the computation and storage space requirements because the problem is reduced to a one-state, one-stage DP problem.

Yoshimoto and others (1988) introduced the look-ahead

process to establish that the optimum solution has been

located at each stage. Then Yoshimoto and others (1990)

introduced a region-limited strategy combined with the

PATH algorithm (RLS-PATH) that has two parts. In part

one, a restricted domain of the state space is established. 

In part two, DP (via the PATH algorithm) is applied to the

restricted domain. Once the optimal solution has been

located, a new, tighter restricted domain of the state space

is established. The process continues until convergence on

a local optima has been established (Yoshimoto and others

1990). 

The difficulty with attempting to optimize management

prescriptions that include multiple thinnings is that the level

and timing of each thinning affect the potential intensity and

timing of subsequent thinnings (Amidon and Akin 1968).

In problems with three of more thinnings, the RLS-PATH

algorithm has been shown to perform better than the Hooke

and Jeeves method, and with less computation time

(Yoshimoto and others 1990).

This paper will look at the use of the RLS-PATH algorithm within an automated prescription generator that uses

a single-tree, distant-independent growth model to project

stand growth. The term “prescription” is used in this paper

to refer to the overall schedule and level of activities that

will occur in a stand over some planning horizon. Inherent

to the term prescription is the fact that individual tree data

is explicitly tracked over time and thus a prescription also

implies knowledge about what a stand looks like over the

planning horizon (some people use the term “regime” to

describe a set of individual silvicultural “prescriptions” to

occur in a stand over time – in this case the terms regime

and our use of prescription are the same).

The example prescription we model is aimed at generating “general” guidelines for reducing stand density index

values for stands in the Blue Mountains of eastern Oregon.

The ideal of a “general guideline” is meant to address the

following question – if a stand has certain characteristics

(which may be shared by other stands in some landscape)

and the desire is to reduce the stand density to levels that

256
are within some specified range, what are the general thinning (or not thin) timings and levels that could be applied

to any such stand in that landscape? One solution to this

question is to simply solve the problem for each individual

stand. An approach we chose to take is to solve the problem for some subset of stands that have characteristics in

common, and then evaluate whether timing of entry and

level of activity occur with some regularity in the individual

solutions. If so, that may be suggestive that certain prescriptions are highly correlated with certain stand characteristics,

and thus can be considered a general prescription for similar stands sharing those characteristics (as long as the desired

goal for the stands are the same). There is an analogy in

this method to that of “expert opinion.” Silviculturists often

determine what to do to achieve certain stand characteristics based on their own experience with similar stands. In

essence, this latter method could be considered an “expert

opinion” method as well. A rigorous statistical analysis of

this approach would be prudent for future work in this

area.

PROBLEM FORMULATION

The example problem addressed in this paper is best

described as an any-aged stand management problem. Haight

(1987) and Haight and Monserud (1990a) describe the anyage stand management problem as determining the best

time sequence of harvest and planting levels without constraints on the stand age or size structure. Without constraints on the stand age or size structure the optimal

management solution may yield stands that have irregular

structures or the solution may yield stands with conventional even- or uneven-aged structure – thus the term anyaged is used. This paper will not be determining the optimal

planting levels as regeneration is a stochastic element within our growth and yield model.

For an individual-tree growth model, a stand is represented by a tree list which contains a record for each tree

found in the stand. These records usually include an expansion factor indicating how many trees per unit area are

actually present. Additionally, a suite of attributes could be

associated with each tree record, including: species, diameter at breast height (dbh), height, and crown ratio. 

The problem formulation is couched in terms of harvesting decisions and growth predictions that occur at 

discrete time intervals (every 10 years). The problem is to

determine, for each stand, the number of trees to cut in

each period over the planning horizon, while minimizing

the penalty scalar function associated with the stand's ability

to maintain an acceptable stand density index measurement.

This paper is forgoing the traditional use of an economic

criterion as the objective function. Future work by the

authors will analyze the cost and benefits of any-aged stand

management by including economic criteria. 

The state and control vectors and variables are as 

follows:

s(t) = stand state in period t (vector) [defined by dominate species and basal area]

si(t) = the attributes of tree record i at the beginning of

period t (element)

x(t) = cutting levels in period t (vector)

xi(t) = percentage of record i cut in period t (element)

h(t) = the harvest controls in period t (vector)

hj(t) = percentage of trees cut in diameter group j and

hj

(t) = [0,1] (element)

Cj = a unique diameter class group

J = total number of diameter class groups = 4 (scalar)

P = number of look-ahead periods = 5 (scalar)

T = total number of simulation periods = 10 (scalar)

N = total number of tree records in stand

A mapping occurs within the program such that h → x by

the equation:

xi(t) = hj(t) i ∈ Cj, J [1]

The stand density index penalty is defined by comparing the stand’s calculated stand density index (SDI – using

Reineke’s summation method) value for the evaluation

period and the next P periods against the calculated minimum and maximum SDI value for those periods.

d[s(t)] = a stand’s SDI value with state s(t) at the end of

period t

bu[s(t)] = upper SDI bound for a stand with state s(t)

bl

[s(t)] = lower SDI bound for a stand with state s(t)

Therefore, the penalty scalar function ( c[s(t), h(t)] ) can be

defined as:

t+P

c =

{∑ (d[s(t)] – bu,l[s(t)])2 if d[s(t)] > bu[s(t)] t or d[s(t)] < bu[s(t)] 

0 if bl

[s(t)] ≤ d[s(t)] ≤ bu[s(t)]

[2]

Equation [2] is then used in the objective function as:

T

Min Z = ∑c[s(t), h(t)] [3] t=0

When equation [3] is solved, the thinning levels hj(t) will 

be known for each j and t.

257

A A
SOLUTION METHOD

The optimization problem of equation [3] is solved

using an RLS-PATH process described by Yoshimoto and

others (1990). This algorithm allows us to partition the

optimization problem into partial optimization problems,

one for each stage. It should be noted that the RLS-PATH

was developed out of the need to reduce computational and

storage space from that required for a traditional dynamic

programming formulation. As such, RLS-PATH works

within a restricted subspace of traditional DP and provides

“good” but not necessarily optimal solutions, but with much

less computational and storage size requirements. A more

thorough examination of the PATH and RLS-PATH algorithm’s is given in Yoshimoto and others (1990).

A program called Stand Level Optimization with

Multiple Objectives (SLOMO) was written in the C and

C++ languages to solve stand optimization problems of the

type given in this study. SLOMO is designed with a graphical user interface (GUI) that allows a user to enter multiple

stands and solve for multiple objectives. The user is required

to go through a simple 9-step process to set up and specify

all the necessary parameters, inputs, and outputs (figure 1).

Inside of SLOMO the growth and yield functions from the

Blue Mountain variant of the Forest Vegetation Simulator

(FVS) were incorporated (Dixon 2003). The original FVS

Fortran computer code was literally translated into C and

C++ code. This allowed for faster computational times as it

removes the overhead associated with running FVS in its

stand-alone format.

The harvest control vector, h(t), represents the fraction 

of trees harvested by diameter class for a specified period.

Four meaningful diameter class breaks were defined for this

study (0-4, 4-12, 12-21, 21+ in.). Meaningful in the sense

that these four classes strike a balance between having a

useful management utility (for either wildlife, ecological,

or economic uses) and the desire to have breaks that allow

enough silvicultural flexibility to manage the stand. The

potential hj

(t)’s are a function of which diameter class is

being evaluated. For the bottom three diameter classes (0-

21 in.) the elements are 0, 20, 40, 60, 80, and 98. These

values can be interpreted as follows: if hj(t) = 20 for a particular iteration during the optimization process, then 20%

of the trees per unit area from each record in class j will be

harvested. This could be considered a proportional method

of harvesting because the value of the element is applied

proportionally to each record in the class. The top diameter

class (21+ in.) has a slightly different setup. The elements

for this class are 0, 25, 50, 75, and 98. These values can be

258

Figure 1—Main GUI interface for

the SLOMO program.
interpreted as follows: if hj

(t) = 25, then 25% of the total

trees per unit area in that diameter class are harvested, either

starting from the top or bottom record within that class (i.e.

that record with the largest or smallest diameter). This could

be considered an un-proportional method of harvesting

because the element is applied un-proportionally to each

record in the class. 

A second purpose of this paper is to illustrate the flexibility of SLOMO in evaluating optimal prescriptions with

various constraints. We have already outlined the basic formulation for a problem that attempts to minimize the penalty

scalar function associated with a goal of reducing a stands

SDI value over time. The earlier formulation was for an

unconstrained problem. SLOMO has some built-in constraints that can be toggled on or off for particular stand

goals. Figure 2 shows Step 6 of the SLOMO process, whereby the user has the ability to determine the constraints, the

levels of those constraints, to a pre-defined set of stand

goals. In figure 2, the section titled “Goal 1” is the area

applicable for this study. In this section there are two constraints that can be used – an upper and a lower size limit.

In this study, the lower size limit was 7 inches (i.e. harvesting was not allowed for any tree less than 7 inches). The

upper size limit was 21 inches (i.e. no harvesting of any

tree greater than 21 inches) for one example in this study,

and was unlimited for a second example. These constraints,

when applied, now need to be included and linked in our

model formulation with the following scalar function

L[s(t), x(t)] :

T

Min Z = ∑c[s(t), h(t)] [3] t=0

subject to

L = {xi

(t) = 0 if dbh of sj(t)<7

xi(t) = 0 if dbh of sj(t)>21 [4]

The SDI calculations are embedded functions with the

FVS system and were simply transferred into the SLOMO

code. The current SDI for a stand is calculated by the summation method suggested by Long and Daniel (1990) and

formally proved to be applicable to irregular structured

stands by Shaw (2000).

dbhi ∑[TPAi

( )1.605] for i = 1,....N [5] 10

A maximum SDI is then calculated for each stand based

on a pre-calculated theoretical maximum SDI (a required

input into the SLOMO model) that is then proportionately

weighted by the basal area of the species actually present.

The study in this paper used an upper and lower bound of

55 and 35% of the calculated weighted maximum SDI as

the bu[s(t)] and bl[s(t)] values as seen in equation [2]. These

values were picked because there is evidence to suggest

that they have significant silvicultural implications. Note

from figure 2 that SLOMO allows the user the ability to

change these values for any run requesting the application

of Goal 1.

259

Figure 2—Step 6 of SLOMO –

Specify Goal Parameters.
STUDY DESIGN

Stand 21323623 (hereafter just called the “stand”) is

located in the La Grande Ranger district of the WallowaWhitman National Forest and was selected as the example

stand for this study. More specifically, the stand is located

in the Upper Grande Ronde Watershed and has an aspect of

10°, a slope of 54% and elevation of approximately 4,600

feet. Tables 1 and 2 show a detailed breakdown of the stand

components. The species present include western larch

(Larix occidentalis Nutt.), Douglas-fir [Pseudotsuga menziesii (Mirb.) Franco], grand fir [Abies grandis (Dougl. ex

D. Don) Lindl.], lodgepole pine (Pinus contorta Dougl. ex

Loud. ), and ponderosa pine (Pinus ponderosa Dougl. ex

Laws.) (hereafter referred to as WL, DF, GF, LP, and PP,

respectively).

Recall that the main objective of this paper is to demonstrate a process by which general prescription guidelines

could be created for reducing and maintaining SDI values

in stands that have similar characteristics. Additionally, we

wish to demonstrate the ease of using SLOMO for analysis

of this sort in which constraints can be turned on and off to

provide the user a surrogate measure of marginal cost for

applying those constraints (or different levels of constraints).

As can be seen in tables 1 and 2, the example stand is initially dominated by DF and falls in our basal area category

of 100 to 150 ft2 ac-1. Thus we identify this stand as having

the characteristic DF-100. Rephrasing our main objective

with this characteristic in mind we are asking, “For stands

in the Blue Mountains of eastern Oregon that are dominated

by DF and have a basal area of 100 to 150 ft2 – what are

general prescriptions guidelines that we could follow that

would reduce and maintain SDI values to levels we find

acceptable?” Furthermore, we are also asking, “What are

the tradeoffs of allowing harvest to occur in trees over 21

inches versus not allowing harvest of those trees, when

searching for an answer to the previous question?”

Ideally we would have a number of stands with the

characteristic DF-100 to run through SLOMO, and we

would analyze the results to obtain general prescription

guidelines. However, time constraints prevented us from

accomplishing that for this particular study. In replacement,

we ran one stand through SLOMO ten independent times.

Although the initial conditions are exactly the same we

suspect that the stochastic elements of the embedded growth

and yield model create enough “displacement” of stand

conditions as soon as the first growth period is encountered. In other words, even if the stand was run through the

growth and yield model with no thinning or optimization

occurring (a grow-only simulation), the stochastic elements

of the growth and yield model will create different stands

even after a single growth period. This displacement is

accentuated by our optimization process because even

slight differences in stand characteristics may cause slight

differences in the solution hj

(t)’s that are found for independent runs. 

260

Table 1—Initial stand conditions.

TPAa BAb SDIc CCFd QMDe CF/ACREf

1009 139 233 150 5.0 4261

a Trees per acre

b Basal area in ft2 per acre

c Stand density index

d Crown competition factor

e Quadratic mean diameter

f Cubic feet per acre

Table 2—Species breakdown of some initial stand attributes.

Stand Attribute WL DF GF LP PP

BA 6 107 24 0 0

QMD 14.8 13.9 29.3 0.1 0.1

TPA 5 103 5 797 100
RESULTS AND DISCUSSION

Stand prescriptions were first optimized with SLOMO

using no upper size limit constraint, and a lower size limit

of 7 inches. The current view of Step 6 as seen in figure 2

shows the exact configuration of the SLOMO interface to

accomplish this task. Table 3 shows a summary of the

stand data for one of the ten independent runs over a 100

year planning horizon. Stand prescriptions were then optimized with SLOMO using an upper size limit of 21 inches

and a lower size limit of 7 inches. Enabling the upper size

limit was simply a matter of selecting the “yes” button as

seen in the Goal 1 section of figure 2. Table 4 shows a

summary of the stand data for one of those runs.

The main objective of this study was to determine if

there were any patterns in the prescriptions that appear

with regularity for stands with similar characteristics, and

figures 3 and 4 enable us help us to locate some of these

patterns. Figure 3, for example, represents the ten runs in

which there was no upper size limit constraint imposed

upon the objective function. Here we can see that in time

period 1, 80% of the solutions scheduled harvests in the

21”+ diameter class (i.e. eight of the ten solutions). Additionally, 20% of the solutions scheduled harvests in the

12-21” class, and 60% scheduled harvests in the 4-12”

class. Note that harvesting can occur in multiple diameter

classes during any single time period. Recall also that a

lower size limit of 7” was used, so we expected harvests to

261

Table 3—Data for one independent run with no upper size limit 

constraint.

Year TPA BA SDI CCF QMD CF/ACRE

2000 1010 139 233 150 5.0 4261

2010 1437 139 240 153 4.2 4437

2020 1367 158 258 172 4.6 5352

2030 1312 177 276 191 5.0 6152

2040 1264 195 293 208 5.3 6988

2050 1201 208 304 220 5.6 7716

2060 1099 214 304 223 6.0 7917

2070 1009 218 302 225 6.3 8219

2080 929 223 303 229 6.6 8544

2090 858 226 302 230 7.0 8637

2100 801 230 301 232 7.3 8857

Table 4—Data for one independent run with a 21-inch upper size

limit.

Year TPA BA SDI CCF QMD CF/ACRE

2000 1010 139 233 150 5.0 4261

2010 928 118 189 119 4.8 3962

2020 882 134 204 158 5.3 4596

2030 829 132 201 156 5.4 4573

2040 798 146 220 168 5.8 5155

2050 762 145 217 165 5.9 5251

2060 717 152 224 169 6.2 5453

2070 666 156 227 169 6.6 5477

2080 619 158 227 169 6.9 5522

2090 578 160 228 169 7.1 5563

2100 544 163 229 170 7.4 5693
be absent from 0-4” class. What is not known from viewing this figure is the intensity of harvests scheduled for

each diameter class (the actual hj(t)’s – note: they are available as output from SLOMO, they are just not displayed

here), but only whether or not harvesting was scheduled 

in that class during that time period. Future work by the

authors will examine the intensity of harvest question.

In examining period 1 of figure 3, one might recommend that for stands with characteristics similar to DF-100,

there is a good likelihood that some level of thinning

should occur in the 21”+ diameter class, and a moderate

likelihood of thinning in the 4-12” class, during the first

decade and if the stand goal is as stated in equation [3].

Interestingly, for this example we can see that for the

remaining planning horizon there is only a very small likelihood of thinning in any later time period when prescriptions are optimized. This can be attributed to the fact that

the thinning in period 1 removed a significant number of

large trees. The stand was initially composed (table 2) of a

DF overstory, yet with a significant component of LP and

PP regeneration. After removing the larger trees, the smaller

trees require a long time to grow before they are available

for harvest (beyond the 7” dbh threshold). Here, we might

say that to achieve the goals and adhere to the constraints

of this problem, a thinning from above (to some level of

density) in the first decade is required.

In contrast, we find a very different pattern in figure 4.

Here the 21” dbh upper size limit constraint is imposed on

the objective function. Period 1 in figure 4 reveals that

every solution required thinning in the 4-12” diameter class.

The optimization process is thus trying to reduce the SDI

by thinning the smaller diameter material. Again, only those

trees with a dbh > 7” are even eligible to be thinned so the

trees being thinned in this diameter class are really those

with dbh’s between 7 and 12 inches. Periods 3 and 5 also

indicate that every solution required thinning in the 12-21”

dbh class. For this set of goals and constraints, we can confidently say that the immediate need (in the first decade)

is to thin in the 4-12” dbh class (to some level of density).

Then return at year 30 and thin in the 12-21” dbh class, and

at year 50 to again thin in the 12-21” dbh class. 

CONCLUSIONS

We have described here a process for optimizing standlevel goals with an RLS-PATH algorithm embedded in a

computer program called SLOMO. The SLOMO program

was developed with a flexible interface to allow users to

develop stand-level goals related to biological or ecological

concerns. We feel that by applying this method to a number

of stands of similar characteristics, that some general thinning guidelines can be developed for broad areas. By assessing figures such as those presented as figures 3 and 4 one

can quickly assess whether or not general guidelines can be

developed for stands that have similar structural characteristics. It is certainly recognized that individual site conditions may preclude the development of general guidelines.

In this case, SLOMO can still be viewed as a tool that can

be used to develop a stand-specific prescription.

262

Figure 3—Potential thinning (or no thinning) timings for a DF-100 stand

with no upper size limit.

Figure 4—Potential thinning (or no thinning) timings for a DF-100 stand

with a 21” upper size limit.
LITERATURE CITED

Adams, Darius M. and Ek, A.R. 1974. Optimizing the 

management of uneven-aged forest stands. Can J. For.

Res. 4: 274-287.

Amidon, Elliot L. and Akin, G.S. 1968. Dynamic programming to determine optimum levels of growing

stock. Forest Science. 14(3): 287-291.

Arthaud, Greg J. and Klemperer, W.D. 1988.

Optimizing high and low thinnings in loblolly pine with

dynamic programming. Can J. For. Res. 18: 1118-1122.

Bare, B. Bruce and Opalach, D. 1987. Optimizing 

species composition in uneven-aged forest stands.

Forest Science. 33(4): 958-970.

Brodie, J. Douglas, Adams, D.M., and Kao, C. 1978.

Analysis of economic impacts on thinning and rotation

for Douglas-fir using dynamic programming. Forest

Science. 24(4): 513-522.

Brodie, J. Douglas and Kao, C. 1979. Optimizing thinning in Douglas-fir with three-descriptor dynamic programming to account for accelerated diameter growth.

Forest Science. 25(4): 665-672.

Brukas, Vilis and Brodie, J.D. 1999. Economic optimisation of silvicultural regimes for Scots pine using

dynamic programming. Baltic Forestry. 5(1): 28-34.

Bullard, Steven H., Sherali, H.D., and Klemperer, W.D. 

1985. Estimating optimal thinning and rotation for

mixed-species timber stands using a random search

algorithm. Forest Science. 31(2): 303-315.

Buongiorno, Joseph and Michie, B.R. 1980. A matrix 

model of uneven-aged forest management. Forest

Science. 26(4): 609-625.

Buongiorno, Joseph, Dahir, S., Lu, H.-C., and Lin, C.-

R. 1994. Tree size diversity and economic returns in

uneven-aged forest stands. Forest Science. 40(1): 

83-103.

Chappelle, Daniel E. and Nelson, T.C. 1964. Estimation 

of optimal stocking levels and rotation ages of loblolly

pine. Forest Science. 10(4): 471-502.

Chen, Chung M., Rose, D.W., and Leary, R.A. 1980a. 

How to formulate and solve “optimal stand density over

time” problems for even-aged stands using dynamic

programming. U.S. Dept. of Ag., Forest Service, North

Central Forest Experiment Station. St. Paul, MN; GTRNC-56. 17 p.

Chen, Chung M., Rose, D.W., and Leary, R.A. 1980b. 

Derivation of optimal stand density over time - A discrete stage, continuous state dynamic programming

solution. Forest Science. 26(2): 217-227.

Dixon, Gary E., (Comp). 2003. Essential FVS: A user’s 

guide to the Forest Vegetation Simulator. U.S. Dept. of

Agriculture, Forest Service, Forest Management Service

Center. Fort Collins, CO. 193 p.

Gong, Peichen. 1992. Multiobjective dynamic programming for forest resource management. Forest Ecology

and Management; 48: 43-54.

Gove, Jeffrey H. and Fairweather, S.E. 1992. 

Optimizing the management of uneven-aged forest

stands: a stochastic approach. Forest Science. 38(3):

623-640.

Haight, Robert G. 1987. Evaluating the efficiency of 

even-aged and uneven-aged stand management. Forest

Science. 33(1): 116-134.

Haight, Robert G., Brodie, J.D., and Dahms, W.G. 

1985. A dynamic programming algorithm for optimization of lodgepole pine management. Forest Science.

31(2): 321-330.

Haight, Robert G. and Monserud, R.A. 1990a. 

Optimizing any-aged management of mixed-species

stands. I. Performance of a coordinate-search process.

Can J. For. Res. 20: 15-25.

Haight, Robert G. and Monserud, R.A. 1990b.

Optimizing any-aged management of mixed-species

stands: II. Effects of decision criteria. Forest Science.

36(1): 125-144.

Haight, Robert G., Monserud, R.A., and Chew, J.D. 

1992. Optimal harvesting with stand density targets:

Managing Rocky Mountain conifer stands for multiple

forest outputs. Forest Science. 38(3): 554-574.

263
Hann, David W. and Brodie, J.D. 1980. Even-aged management: Basic managerial questions and available or

potential techniques for answering them. U.S. Dept. of

Ag., Forest Service, Intermountain Forest and Range

Experiment Station. Ogden, UT; Gen. Tech. Rep. INT83. 29 p.

Hooke, Robert and Jeeves, T.A. 1961. “Direct Search” 

solution of numerical and statistical problems. J. Assoc.

Comput. Mach. 8: 212-229.

Hool, James N. A. dynamic programming-Markov chain 

approach to forest production control. Monograph 12,

Forest Science. 66.

Kao, Chiang and Brodie, J.D. 1980. Simultaneous optimization of thinnings and rotation with continuous

stocking and entry intervals. Forest Science. 26(3): 

338-346.

Long, James N. and Daniel, T.W. 1990. Assessment of 

growing stock in uneven-aged stands. West. J. Appl.

For. 5(3): 93-96.

Martin, George L. and Ek, A.R. 1981. A dynamic programming analysis of silvicultural alternatives for red

pine plantations in Wisconsin. Can J. For. Res. 11: 

370-379.

Palahí, M., and Pukkala, T. 2003. Optimising the management of Scots pine (Pinus sylvestris L.) stands in

Spain based on individual-tree models. Annals of Forest

Science. 60: 105-114.

Paredes V., Gonzalo L. and Brodie, J. D. 1987. Efficient 

specification and solution of the even-aged rotation and

thinning problem. Forest Science. 33(1): 14-29.

Roise, Joseph P. 1986a. A nonlinear programming 

approach to stand optimization. Forest Science. 32(3):

735-748.

Roise, Joseph P. 1986b. An approach for optimizing 

residual diameter class distributions when thinning

even-aged stands. Forest Science. 32(4): 871-881.

Shaw, John D. 2000. Application of stand density index to 

irregularly structured stands;. West. J. Appl. For. 15(1):

40-42.

Thorsen, Bo Jellesmark and Helles, F. 1998. Optimal 

stand management with endogenous risk of sudden

destruction. Forest Ecology and Management. 108: 

287-299.

Valsta, Lauri T. 1990. A comparison of numerical methods 

for optimizing even-aged stand management. Can J.

For. Res. 20: 961-969.

Wedin, Heidi L. 1999. Stand level prescription generation 

under multiple objectives. Oregon State University;

Corvallis, OR. M.S. Thesis. 77 p.

Wikstrom, Peder. 2001. Effect of decision variable definition and data aggregation on a search process applied to

a single-tree simulator. Can J. For. Res. 31: 1057-1066.

Wikstrom, Peder and Eriksson, L.O. 2000. Solving the 

stand management problem under biodiversity-related

considerations. Forest Ecology and Management. 126:

361-376.

Yoshimoto, Atsushi, Haight, R.G., and Brodie, J.D. 

1990. A comparison of the pattern search algorithm and

the modified PATH algorithm for optimizing an individual tree model. Forest Science. 36(2): 394-412.

Yoshimoto, Atsushi, Paredes V., G.L. and Brodie, J.D. 

1988. Efficient optimization of an individual tree

growth model. U.S. Dept. of Ag., Forest Service, 

Rocky Mountain Research Station. Ogden, UT;

in Gen. Tech. Rep. RM-161. pp. 154-162.

264
In: Bevers, Michael; Barrett, Tara M., comps. 2005. Systems Analysis in Forest Resources: Proceedings of the 2003 Symposium. General Technical Report 

PNW-GTR-656. Portland, OR: U.S. Department of Agriculture, Forest Service, Pacific Northwest Research Station.

1 Eldon Gunn is a Professor of Industrial Engineering, Dalhousie University, Halifax, NS, Canada, B3J 2X4.

Systems Analysis in Forest Resources: Proceedings of the 2003 Symposium

A NEURO-DYNAMIC PROGRAMMING APPROACH 

TO THE OPTIMAL STAND MANAGEMENT PROBLEM 

Eldon Gunn1

INTRODUCTION

Forestry is full of problems that involve the control of

dynamic, stochastic systems with the intent to either optimize some objective functional or to attempt to steer the

system to some desired state. Often these involve large

state spaces that make traditional dynamic programming

prohibitive. Recent work on neuro-dynamic programming

shows considerable promise for solving large dimension

problems with complex dynamics.

Even-age stand harvesting has been represented as a

control problem (Naslund, 1969) and the classic Faustman

(1849) formula can be seen as a control policy developed

from a deterministic view of the system. Some have been

interested in how this control policy should be modified in

the face of uncertain growth dynamics, market prices and

the threat of disaster such as fire or budworm attack. Other

researchers have noted that to properly model growth and

response to various silvicultural interventions, a single state

variable, such as volume used in the Faustman formula, is

inadequate. Most growth models require a representation of

something like stand average diameter, basal area, height as

ABSTRACT

Some ideas of neuro-dynamic programming are illustrated by considering the problem of optimally managing a forest

stand. Because reasonable growth models require state information such as height (or age), basal area, and stand diameter,

as well as an indicator variable for treatments that have been performed on the stand, they can easily lead to very large state

spaces. Researchers have developed sophisticated forward dynamic programming and related algorithms However, an interest in problems that are stochastic in their basic growth dynamics, in market prices and in disturbances, ranging from insects

to fire to hurricanes, requires a different approach. Neuro-dynamic programming algorithms are promising for problems

with large dimension that may lack a simple model of dynamics and stochastics. This paper looks at applying these ideas in

the context of a simple one species model with a view to eventual application to more complex problems.

265

well as an indicator variable of stand nature (natural, plantation, pre-commercially thinned, commercially thinned, etc.);

thus four state variables, three of which are continuous. In

representing multiple species stands, the state dimensions

grow considerably. Even-aged stands are not the only source

of large-scale dynamic stochastic control problems in forestry.

For example, Islam and Martell (1997) describe a multidimensional queueing model of forest fire attack aircraft. If

we concentrate on forest harvesting, Buongiorno (2001) has

modeled uneven-age management as a Markov decision

process. 

The main ideas of this paper are closely based on

Bertsekas and Tsitsiklis (1995) who come from a dynamic

programming tradition. They point out that many of the

ideas of neuro-dynamic programming have their roots in

the artificial intelligence field under the name of reinforcement learning theory (Sutton and Barto, 1998)

DYNAMIC PROGRAMMING

There are many technicalities in properly detailing the

nature of the problem to be solved. Because of lack of
space, we will gloss over many of the details. Those interested in filling these in should consult Bertsekas and

Tsitsiklis (1995) or Sutton and Barto (1998).

The basic dynamic programming (DP) problem involves

a situation where the decision maker is aware of a current

state i and is required to make a decision u. The decision u

puts the system in a new state j at the beginning of the next

stage with probability pij (u). In a finite horizon problem,

we are trying to estimate the value of being in state i with 

k stages remaining, assuming that we will make optimal

decisions in this and all subsequent stages. That is we need

to calculate 

J*

k(i) = Max ∑ pij (u)(g(i,u, j) + α J*k-1( j)) (1) u ∈U(i) j

The explanation of this equation is that in state i at stage k

we want to make a decision that maximizes the expected

value of all present rewards g(i,u,j) and all future rewards

J*

k( j) with the future rewards discounted back to the beginning of the current stage k using the discount factor α. The

functions J*

k( j) are referred to the cost-to-go functions.

We are particularly interested in the idea of policies. By

a policy π= {μ0,μ1, ...}, we have at each stage k, a policy

function μk(i) that maps each state i to a particular decision

μk(i)∈U(i). Given a policy, we can develop the idea of

expected value of that policy

J π

k (i) = ∑ pij (μk(i))(g(i,μk(i),j) + αJπk-1( j)) (2)

j

Another way of viewing the DP equation above is that of

searching for the optimal policy

J*

k (i) = Max Jπk (i) (3) π

As we go to an infinite horizon framework, we are really

letting k→ ∞ so that there is no distinction between having

k periods left to go and k-1 periods left to go. In that case

the equation (1) becomes Bellman’s equation:

J*(i) = Max ∑ pij (u)(g(i,u, j) + αJ*( j)) . (4) u ∈U(i) j

In an infinite horizon framework, it makes sense to think of

stationary policies, policies that depend only on the state

and not the stage, π= {μ,μ,μ,μ, ...}. 

With stationary policies, we can also think of the

expected value of a policy μ being given by: 

Jμ(i) = ∑ pij (μ(i))(g(i,μ(i),j) + α Jμ( j)) (5)

j

Once again we can think of the DP problem as the search

for the optimal policy among all stationary policies: 

J*(i) = Max Jμ(i) (6)

μ

The fundamental idea of DP is the intimate connection

between policies and values. If we know an optimal policy

μ*(i), then we need only solve equations (5) to calculate 

the optimal expected values J*(i) of being in state i. J*(i) is

referred to both as the cost-to-go function or optimal value

function. Alternatively, if we have the cost-to-go function

J*(i), then we solve for the optimal policy μ*(i) by carrying

out the maximization operation in equation (4). 

The problem of finding an optimal policy or finding the

cost-to-go function is usually approached by way of recursive approximation algorithms. There are two main classes

of algorithms that are usually applied, value iteration and

policy iteration. There is also a third class of algorithms,

linear programming, that we will not discuss here (Bertsekas

and Tsitsiklis, 1995)

Value iteration: These algorithms begin with an initial

estimate J0(i) for J*(i). They then use the maximization transformation in equation (1) for some or all of the states i to

build up a new estimate Jk+1(i) from Jk (i). As k→ ∞, if 

we update the estimates for each i often enough, then the 

Jk (i)→ J*(i). Conventionally one starts with a given Jk (i)

and then calculates Jk+1(i) for all states i. As Bertsekas and

Tsitsiklis (1995) point out, this is not necessary and various asynchronous approaches, such as Gauss- Seidel, can

be useful. 

Policy iteration: These algorithms start with a given

policy μk(i) and then calculate a cost-to-go function J μk

that corresponds to the policy μk(i). This requires solving

the linear equations (5) for J μk (i). Given the J μk (i), we

then update the policy to μk(i) by solving the max operation in (4); that is 

μk+1(i) = arg max ∑ pij (u)(g(i,μ,j) + α Jμk ( j)) u ∈U(i) j

for some or all of the states i. We then repeat. As k→∞,

μk(i)→μ*(i). Again, conventional algorithms update the

policy μk(i) for all i in the state space but it is possible to

only update some states at each iteration.

266





 









 









 




Thus, conventional DP consists of successively building

up approximations to the cost-to-go function or the policy

function. We will see that neuro-dynamic programming

involves more extensive approximations. In the case of

policy iteration, solving the equations (5) can be difficult 

to do directly if the state space is large. Some algorithms

iterate using equation (5) to build up an approximation of

J μk (i) and then use this approximation to calculate μk+1(i).

In general, value iteration algorithms are conceptually quite

simple, but often converge slowly. Policy iteration algorithms tend to converge more quickly. 

Stand Management Dynamic Programming Problems

Dynamic programming has been a favourite tool for

examining stand management problems. Although early

work used a single state variable, it is impossible to accurately model growth and the effect of interventions without

taking into account more complex states such as diameter,

height (age) and stand basal area (for example Brodie, Adams

Kao 1978; Brodie and Kao 1979). Recently there has been

a tendency to use individual tree growth models (see

Pelkki, 1994).

An interesting dichotomy in DP models can be observed.

On the one hand, we find models with complex state variables but deterministic dynamics and prices (see Pelkki,

1994, and Pelkki, 2002). On the other, we find models with

stochastic prices and/or dynamics but with very simple state

representation (for example Haight and Holmes, 1991). We

generally don’t find both. 

For the complex state vector typically associated with

detailed growth models, advantages can be gained by using

a forward DP recursion. Here we start with stage 0 corresponding to the initial state and compute

F*

k ( j) = Max{(F*k–1(i) + αk–1g(i,u(i, j), j))}

i

where F*

k ( j) is the maximum net present value of a policy

that begins stage k in state j and u(i,j) is the policy decision

that transforms state i to state j in one stage. (If it is impossible to transform state i to state j in one step, we adopt the

convention that g(i,u(i, j), j)= – ∞). The advantage of a forward recursion is that it is possible to prune states that

cannot be reached at a given stage. If we use a finite set of

possible policy decisions, then only a finite number of

states are possible at each stage. 

Forward recursion requires us to know the state resulting from a decision u in state i. Hence, it cannot be used in

a stochastic setting where the outcome of a decision depends

on a probability distribution. As a result, stochastic DP

problems must use backward recursion. This prevents the

state pruning that is at the heart of forward methods. In

order to keep the size of the state space manageable, stochastic forestry DP models often use a very simplified state

space and with a resulting requirement for overly simplified

growth dynamics. To deal with realistic growth models in a

stochastic setting, we need a way of dealing with equations

(1)-(4) with high dimensional state space. Typically we will

want to use infinite horizon methods in order to pick up the

regeneration process typified by the Faustman formula.

There are many ways in which stochastic issues can

enter into the stand management problem. A few of these

include:

i) The price(s) of the forest products are random. Something like a random walk model has been proposed by

Haight and Holmes, 1991. In this case the DP problem

has some of the nature of optimal stopping problems.

ii) Natural disasters such as budworm, fire and hurricanes 

can force a transition to the initial regeneration state.

iii) Natural regeneration can result in a random number 

of periods from the time that a stand has been regeneration harvested until the regeneration is established.

iv) Weather and natural variation can result in variability 

about the assumed mean growth. Any growth model is

the result of regression with the variability implicit in

this. 

v) Individual tree growth models are inherently stochastic. 

To grow a given forest state (diameter, basal area,

height), it is necessary to randomly generate a forest

stand with a certain number of trees and a distribution

over diameter classes. Growing such a model involves

a probabilistic choice for each tree as to whether it survives for a period or not. Trees that survive are grown

according to a formula developed through regression.

APPROXIMATION ISSUES – 

VALUE ITERATION

Forward dynamic programming typically uses the idea

of a neighbourhood, in which a given state is chosen to be

representative of the continuum of states within a certain

radius. The larger the neighbourhood radius, the smaller

the number of states to be evaluated and the easier the

computations. Pelkki has shown that there are significant

issues with accuracy as the neighbourhoods become larger.

The idea of representative states surrounded by neighbourhoods amounts to a piecewise constant approximation

to the cost-to-go function. The first idea of neuro-dynamic

programming is a continuous approximation. That is the

267
functions J*

(i) (resp. Jμ (i), and Jk (i)) are approximated 

by continuous functions J *(i,r) (resp. Jμ(i,r), and Jk(i,r)),

where the parameter vector r has been chosen to minimize

the approximation error at a finite number of states i1, i2,

..., iN.

There are many ways to do the approximation. Simple

linear models (example: quadratic fit) are one possible

approach. If the state variables are of the form (BA, D, H,

tt) with continuous variables for basal area (BA), diameter

(D) and height (H) and the discrete variable treatment type,

it may be reasonable to use a separate continuous approximation on (BA, D, H) for each of the discrete values tt.

The fitting approach that leads to the name neuro-dynamic

programming is to let r be the parameters of a neural network architecture used to approximate J*

(i) (resp. Jμ (i), and

Jπ

k (i)). 

One idea discussed in Bertselas and Tsitsiklis (1995) is

to use features, indicator variables for particular knowledge

of the system under control. For example, it is easy to characterize certain states (BA, D, H, tt) as being ineligible for

certain treatments. In a red spruce forest, there are certain

feasible diameter ranges for pre-commercial thinning, commercial thinning and final felling. A feature can either be

used to create separate linear models for each feature value

or as an additional input to the neural network architecture. 

This then suggests an approximate value iteration algorithm along the following line. 

i) Begin with estimates Jk(ij

), calculated at Nk points 

i1,i2,...,iNk.

ii) Fit the continuous function Jk(i,r) to minimize 

the approximation error over i1,i2,...,iNk.

iii) Let k=k+1 and identify states i1,i2,...,iNk.

iv) Solve the recursion 

J*

k (i) = Max ∑ pij (u)(g(i,u, j) + α J*

k-1( j,r)) u ∈U(i) j

at the points i1,i2,...,iNk.

v) Either stop or go to step ii) and repeat.

Similar approximate value iteration has been successfully used in the engineering community. Johnson and 

others (1993) report on using cubic spline approximations

for hydroelectric power systems with five continuous state

variables, evaluated at as few as five points along each

state dimension. Shoemaker (2003) indicates an extension

of these ideas to nine-dimensional stochastic DP problems. 

This approach is a straightforward extension to standard

DP ideas. When compared to the idea of neighbourhoods,

there is now no longer the requirement to map the single

state representing one neighbourhood at a given stage to

another single state representing a subsequent neighbourhood. The question of accuracy is now that a sufficient

number of states be examined so that the best fit approximating function is a good representation of the true costto-go function. However, this introduces a number of complications. The first is whether the approximating architecture Jk (i,r) is capable of modeling the actual Jk (i),

the second is whether the fitting of the parameter vector r

has been successful in minimizing the fitting error (except

in the case of linear models where this problem is trivial)

and the third is whether the i1, i2, ..., iNk are numerous

enough and chosen appropriately. 

Trying to develop this into an approach for the stand

management problem raises additional issues. We are normally more interested in policy than cost-to-go functions

(although cost-to go functions do have an important interpretation as economic rent). It is unlikely that the i1, i2, ...,

iNk are sufficiently dense to reflect all forest states at which

we would like to know a thinning/harvesting policy. There

are two solutions here. One is to compute policy at a given

state i′ on the fly by explicitly solving the equation 

μ*

k (i′ ) = arg max ∑ pi′ j (u)(g(i′ ,μ,j) + α J*

(j,r)) . u ∈U(i′ ) j

The other is to compute the policy 

μ*

k (i′ ) = arg max ∑ pi′ j (u)(g(i′ ,μ,j) + αJ*

(j,r)) u ∈U(i′ ) j

for a set of points i′ ∈ {i′1, i′2, ...} and then develop another

approximation architecture μ*(i,r). In forestry, although the

components of the state vector such as BA, D and H are

continuous variables, it is realistic to think of a reasonably

small finite set of decisions. In pre-commercial thinning

this might consist of spacing at 6, 8, 10 or 12 ft (2, 2.5, 3,

3.5, 4 meters) and in the case of commercial thinning, the

removal of 40%, 50%, 60%, 100% of basal area. In both

cases, it is unrealistic to be more precise in giving instructions to thinning crews. Thus the fitting problem μ*(i,r)

really becomes a classification problem, a task to which

neural networks are well suited. The classification problem

is to classify the states by the appropriate policy action for

that state.

268





 









 









 




POLICY ITERATION

There are several reasons to consider approximate policy

iteration algorithms. First, the main focus of forest management is on policy. Most jurisdictions will have a current

policy that governs its silviculture. A very natural question

is can this policy be improved. Second, conventional wisdom is that policy iteration algorithms often converge faster

than value iteration. There may well be a third reason. In

many cases, the probabilities that govern the state transitions may not be known analytically. This is certainly the

case when an individual tree growth model is used to model

stand dynamics. The evaluation of a cost-to-go function

(equation (5)) can be carried out using simulation when an

analytic expression for pij (u) is unavailable. To see the

possibilities, we will describe approximate policy iteration

assuming analytical probabilities are not available. 

Policy iteration (and indeed value iteration) works best

if we start with a reasonable initial policy (cost-to-go function of that policy). An example might be a simple diameter

limit policy and simple rules for pre-commercial and commercial thinning. Once we have that policy, we can simulate it to estimate a cost-to-go function. There are a variety

of issues in simulating a policy. Some of these are easier 

if we assume that the dynamic programming problem has

been formulated as a stochastic shortest path problem

(Bertsekas and Tsitsiklis, 1995, show this can always be

done). A stochastic shortest path is a DP problem in which

one of the states, say state 0, is absorbing with p00(u) =1

and g(0,u,0)=0 for all u. Starting from some initial state i,

and using the established policy, we simulate the current

returns and next state until we eventually encounter the

absorbing state. Thus we construct a simulation path (i0, i1,

i2, ..., 0). Note this simulation run also constitutes a costto-go sample for every unique state encountered on this

path. There are a number of issues about what to do if a

state is encountered more than once on a simulation path.

See Bertsekas and Tsitsiklis (1995) for a more extensive

discussion on simulation issues. Through sufficient repeated

simulations and judicious sampling, we can build up costto-go estimates at a number of states. We then use fitting

procedures as discussed earlier to produce an approximation ∼

Jμ(i,r). 

We then need to compute an improved policy that 

satisfies 

μ(i) = arg max ∑pij(u)(g(i,u,j) + α∼J(j,r)) . u ∈U(i) j

There are again many ways to do this. If probabilities pij(u)

are analytically available, then for a finite number of choices

in U(i), we can compute the so-called Q-factors directly by 

Q(i,u,r) = ∑ pij(u)(g(i,u,j) + α ∼J(j,r)). j

If the pij(u)are not readily available, this computation can

be carried out by simulation and this simulation is easier

than that discussed above because it is only a one stage

simulation. Then the improved policy is 

μ(i) = arg maxQ(i,u,r). u ∈U(i)

Given this improved policy we can repeat this process

of i) developing an estimating cost-to go function ∼J(i,r), ii)

calculating the Q-factors Q(i,u,r), and iii) calculating the

next improved policy μ(i)=arg maxQ(i,u,r). u ∈U(i)

At present, we are working on problems where the pij(u)

are available and easily computed. However, if we use the

types of individual tree models used by others, such as

Pelkki (1994), this will no longer be the case. 

EXAMPLE APPLICATION

In this section, we give a brief outline of a DP model ,

where the objective is to find an optimal policy of harvest

and silviculture for even-aged management of the stand. At

harvest (clearcut), the stand is subject to an uncertain period of regeneration. This can be eliminated through planting

the stand. Unplanted stands, once regenerated, grow as

natural stands. Within certain diameter limits, they can be

pre-commercially thinned (PCT). Plantations can either have

either early competition control or not. Those without early

competition control have a certain probability of reverting

to a well stocked natural stand. Plantations or PCT stands

can be commercially thinned (CT) once they reach certain

diameter limits. In both the case of PCT and CT, there are

a range of density reduction choices available. The model

maintains a four dimensional stand state; HT, height, D –

diameter, BA basal area and tt – treatment condition. The

DP problem is to develop a policy of what management

decision to make at the beginning of a period for any state

of the stand. The decisions include those indicated above 

as well as the decision to do nothing and let grow for one

period. We use five-year periods in this work. The current

model is based on a Red Spruce conifer stand. 

The growth models are based on the Nova Scotia

Softwood Growth and Yield Model (GNY) (NSDNR, 1993).

This model is based on extensive field analysis using the

269





 




series of research permanent sample plots established and

maintained by the N.S. Department of Natural Resources.

The main parameter is site capability expressed as Site

Index height at age 50. The model maintains height, diameter and basal area from period to period. Merchantable and

sawlog volumes are computed using Honer’s equations

(NSDNR,1996). GNY has essentially three growth models

corresponding to:

Natural Unmanaged Stands—These stands are grown

using equations developed for the provincial revised normal

yield tables (NSDNR, 1990). Age implies height (and vice

versa), height implies diameter and site density, which in

turn implies basal area. The relationship between height

and age offers the modeler the choice between using height

or breast height age as the state variable. However, height

is the easiest variable to measure and because of the stochastic issues involving the time to reach breast height, we

choose to use it.

Natural Stands That Have Been Subject To Precommercial And Commercial Thinning—These stands are

grown using a variable density model. The age-height relationship characteristic of the site is maintained. The stand

average diameter is modeled as a function of site, spacing

and free age. Free age is the mechanism that we use to deal

with density modifications. This function is invertible so

that free age can be calculated as a function of site, spacing

and diameter. Given a certain site, spacing and diameter,

we use the inverse function to compute the free age. If we

choose to leave the spacing as is, then growth amounts to

increasing the free age by one period (five years) and calculating the resulting diameter. If we do a thinning this alters

both the spacing and diameter. The nature of the diameter

change depends on the nature of the thinning (thinning from

below, or above). Given the new diameter and spacing, we

can recomputed the free age and continue to grow the stand

by incrementing the free age. The maximum density is a

function of the diameter. The current diameter and basal

area define both the current density and the spacing. Trees

are grown by advancing the age (height age) and free age

by five years. The diameter and height are then computed.

The site density index (SDI) is computed for this diameter

and, if the current density exceeds SDI, the current density

is reduced to SDI, spacing is recalculated and free age

recomputed to correspond to the new diameter, spacing 

and diameter.

Plantations That May Also Be Thinned—Plantations

are also grown using the variable density growth model but

the maximum density for plantations is computed using a

different equation.

The dynamic programming state variable is thus

(HT,BA, D, tt) where HT is height (m), BA is basal area

(m2/ha) and D is diameter (cm). The tt is a treatment indicator with the following values: 1) natural, unmanaged, 2)

natural, pre-commercially thinned, 3) natural, commercially

thinned, 4) plantation, spaced, 5) plantation, commercially

thinned

For growth purposes, there is no necessity to distinguish

between tt=2 and tt=3 or to distinguish between tt=4 and

tt=5. However, the nature of revenue generated and the

actions that can be taken are different in these circumstances. This is what Bertsekas and Tsitsiklis (1995) refer

to as a feature.

In our calculations, for each treatment tt, we have an

evaluation of the cost-to-go function at a number of discrete states i1

tt

,... , iNtt

tt where each of the states i

k

tt is of the

form i

k

tt=(HTktt,BAktt,Dktt,tt).We fit the functions Jk(i,r) as

five separate linear models, one for each tt state. Currently

we are working with the vector r being the parameters of 

a second order model; that is J is a quadratic on the three

variables HT, BA and D (ten parameters for the quadratic).

We perform a separate fit for each of the five different values of tt. This means that there are a total of 50 parameters

(10x5) that make up the r. This can be reduced since for the

states with tt=1, the BA and D are determined from HT so

that a ten term model is not necessary. In this case, J is fit

as a quadratic depending on HT only.

The main issue of probabilities in the current model

involves regeneration. When a stand is harvested this period

and successfully regenerates, i) it can exist as a five year

old natural stand in the next period. ii) if it has advanced

regeneration it may exist as a ten year old (or older) stand,

and iii) if it has a regeneration failure it may exist as a stand

with state (0,0,0,1). Stands in this state can either be planted

and converted to plantation or can be left to natural regeneration which may experience continued regeneration 

failure. The state dependent probabilities of successful

regeneration, advanced regeneration and initial and continued regeneration failure define a stochastic model with

computable pij(u). 

We have some very preliminary results from a code

written in Visual Basic for Applications using Microsoft

Excel as the interface. The code is not written for computational efficiency. All growth and volume calculations are

being done on the fly. Although we believe our implementation to be correct, it has not undergone extensive validation. The model has been implemented as a value iteration

270
algorithm with a total of 943 states used to directly evaluate

the cost-to-go function and a quadratic model fit to these

states. The states are broken down into the five treatment

type tt classes as follows: 32 (tt=1), 196 (tt=2), 260 (tt=3),

195 (tt=4), 260 (tt=5). The states are chosen to span the

attainable heights, diameters and basal areas. Note that for

three continuous state variables and five treatment types,

this is a quite modest number of states. The computations

reported below are only for the purposes of indicating that

this approach is practical and not meant to be a real study.

The growth is based on Site Index 50 = 17 m (55 ft) which

corresponds to Land Capability=6m3/ha/year. The computations were done with stumpage values of $15 per m3 for

pulp and $50 per m3 for studwood and sawlogs. Planting,

pre-commercial thinning and commercial thinning treatments

were assumed to be costless. Starting from a cost-to-go

function that is identically zero, figure 1 shows the evolution of the cost-to-go of the bare land state (0,0,0, natural)

as the DP iterations proceed. 

The cost-to go functions essentially stabilize after 60

iterations. The computed policy does not change from

iteration 21 onward.

What this result indicates is that it is possible to do these

computations with fairly modest effort and that convergence is reasonably rapid. A more extensive computational

study will be necessary before we can report on accuracy

and optimality of the computed results. More realistic costs

and revenues will be required before it makes sense to comment on the nature of the policies calculated. 

SUMMARY

Forestry will continue to pose important problems that

involve the solution of large scale stochastic dynamic programming problems. Neuro-dynamic programming is an

approach that shows considerable promise. The ability to

deal with large dimension state spaces, and models that

may exist only as a simulation, are important features of

the approach. 

ACKNOWLEDGEMENTS

This work is partially supported through a Discovery

Grant from the Natural Sciences and Engineering Research

Council of Canada. Mr. Jules Comeau carried out the computations, supported in part by a scholarship from the Baxter

and Alma Ricard Foundation.

271

Figure 1—Cost-to Go Function

Evolution.
(NSDNR) Nova Scotia Department of Natural 

Resources. 1993. Nova Scotia Forest Research Report

No 43, Nova Scotia Softwood Growth and Yield Model

Version 1.0 User Manual. 12 p.

(NSDNR) Nova Scotia Department of Natural 

Resources. 1996. Nova Scotia Forest Research Report

No. 61, Honers Standard Volume Table Estimates

Compared to Nova Scotia Stem Analyses. 16 p.

Pelkki, M., 2002. Dynamic Programming (web page). 

Available at: http://www.afrc.uamont.edu/pelkkim/ 

forest_economics.htm [Viewed October 5, 2003].

Pelkki, M.H. 1994. Exploring the effects of aggregation 

in dynamic programming. In: J. Sessions and J. Brodie

(Eds.) Management Systems for a Global Economy

with Global Resource Concerns, Pacific Grove,

California, USA, pp. 181-190.

Shoemaker, C.A. 2003. Optimization of Uncertain 

Dynamic Systems: Algorithms and Applications to

Hydropower Reservoirs, Inventory Problems. Abstract

for a talk at the Institute for Mathematics and

Applications. Available at http://www.ima.umn.edu/

decision-making/abstracts/shoemaker2.html. [Viewed

September 9, 2003].

Sutton, R.S.; A.G. Barto. 1998. Reinforcement Learning: 

An Introduction. MIT Press, Cambridge, MA. 432 p

272

LITERATURE CITED

Bertsekas, D.P.; J.N. Tsitsiklis. 1995. Neuro-Dynamic 

Programming. Athena Scientific, Belmont, MA. 512 p.

Brodie, J.D.; D.M. Adams; C. Kao. 1978. Analysis of 

economic impacts on thinning and rotation for Douglasfir using dynamic programming. For. Sci. 24: 513-522.

Brodie, J.D.; C. Kao. 1979. Optimizing thinning for 

Douglas-fir with three-descriptor dynamic programming

to account for accelerated diameter growth. For. Sci. 

25: 655-672.

Buongiorno, J. 2001. Generalisation of Faustmann’s 

formula for stochastic forest growth and prices with

Markov Decision Process models. Forest Science. 

47(4): 466-474.

Faustman, M. 1849. Calculation of the value which 

forestry land and immature stands possess for forestry.

In: Martin Faustman and the Evolution of Discounted

Cash Flow, Institute Paper No. 42, Commonwealth

Forestry Institute, University of Oxford, 1968, 27-55.

Haight, R.G.; T.P. Holmes, 1991. Stochastic price models 

and optimal tree cutting: results for loblolly pine.

Natural Resource Modeling, 5(4): 423-444.

Islam, K.M.S.; D.L. Martell. 1998. Performance of initial 

attack airtanker systems with interacting bases and variable initial attack ranges. Can. J. Forest Res. 28(10):

1448-1455.

Johnson, S.A.; J.R. Stedinger; C.A. Shoemaker; Y. Li; 

J.A. Tejada-Guibert. 1993. Numerical solution of 

continuous-state dynamic programs using linear and

spline interpolation. Operations Research, 41(3): 

484-500.

Naslund, B. 1969. Optimal rotation and thinning. Forest 

Science 15: 446-451. 

(NSDNR) Nova Scotia Department of Natural 

Resources. 1990. Nova Scotia Forest Research Report

No. 22, Revised Normal Yield Tables for Nova Scotia

Softwoods. Available at: http://www.gov.ns.ca/natr/

publications/forpubs.htm#ffr. 45 p.

Continue