Combining spatial and other objectives in forest design

Combining spatial and other objectives in forest design The paper presents a forest planning model for Northern Germany with an objective function which includes an economic, even flow and spatial component. The spatial modeling approach followed in this analysis relies on specific geographic relations among neighboring stands, including the distance between stand centroids and the length of the common boundary between adjacent compartments. To our knowledge, this is a new and more realistic approach in spatial optimization when compared with earlier approaches. The geographic data are derived from a GIS system and stored in a database. The database also includes details of all management regimes and the corresponding objective function coefficients. The method is first applied in a relatively small forest with only 41 compartments. The Simulated Annealing algorithm is used to identify the optimum combination of management regimes. Solutions are presented for different objective function components. When the even flow and spatial components were added in the objective function, the economic objective was only moderately reduced. The method is then also applied using a bigger example, a forest with 591 compartments. The results, obtained within a reasonable processing time, were similar and, because of the large number of compartments, the even flow was perfectly balanced.


Introduction
Most of the world's forests are utilized by man and this implies that the dynamics of a forest ecosystem is not so much an ecological, but predominantly a cultural problem. Traditional forest planning is based on principles of constancy and long term stability. In reality, however, periodic reorientation and frequent changes of forest policy are quite common. Examples of this vicissitude are changes in the preferred silvicultural systems (plantations; near-natural forest management), the preferred tree species (beech, spruce) and the preferred forest structures (even-aged monocultures; uneven-aged multi-species forest). The history of forestry shows that policies change, often several times within the life of a tree (Gadow et al., 2007).
The reality of forest management is not constancy, but constant change. Social, economic and environmental conditions are continuously changing. This implies that the traditional practice of standardizing silviculture is ineffective. We propose that one of the basic requirements of forest planning is the ability to generate and evaluate silvicultural regimes. A silvicultural regime, a "management path" in our terminology, is a specific succession of silvicultural activities during a given period of time. A forest typically consists of a number of compartments. If we define several management paths for each compartment, then each combination of management paths involving all compartments represents one possible future design for the forest as a whole. The obvious task is to identify the best forest design, i.e. the optimum combination of management paths. It has been shown that this concept, known as the Multiple Path theory of forest design (Gadow & Pukkala, 2008), provides a suitable basis for linking the levels of the spatial hierarchy in a forested landscape, for balancing a desired mix of forest services and for integrating varied forms of scientific expertise in the planning process. During the past three decades, the Multiple Path concept has emerged as a universal theory of forest design and successful applications have been described involving a variety of forest types and geographical regions (Ware & Clutter, 1971;Hoganson & Rose, 1984;Gadow, 1991;Lappi, 1992;Rodriguez, 1996;Gadow & Pukkala, 2008). The practical implementation, however, is not a trivial task.
The objective of this contribution is to present a specific modeling approach based on the concept proposed by Pukkala & Kangas (1993). We will develop spatial relations among neighboring compartments, including the distance between compartment centroids and the length of the common boundary between adjacent compartments. The Simulated Annealing algorithm will be used to identify the optimum combination of management paths for two forests in North-Western Germany, and solutions will be presented for different objective function components. We will evaluate the processing times and test the sensitivity of the solution by adding different components in the objective function.

The Problem Formulation
In this section we present the specific problem formulation, including details of the objective function, with particular emphasis on the spatial component, and the growth and thinning models used in generating alternative management paths. Our approach relies on the basic structure which has been used by Pukkala & Kangas (1993).

Objective Function
An important part of the Multiple Path model is the objective function. In this study we consider an optimisation problem with three objective function components:

,
( Where a, b and c are constant factors which are used to balance the relative importance of the three components. EC is an economic criterion which evaluates the cash flow within the planning period. In our specific applications near-natural silviculture and selective harvesting is being practiced and clearfellings are not allowed. Furthermore, public ownership of the forest guarantees that the land may not be sold or used for other purposes except forestry. We therefore use the present net worth criterion proposed by Hille et al. (1999), which ignores the value of the land. EC is thus equal to the discounted cash flow derived from selective harvests during the planning period plus the terminal stumpage value at the end of the planning period.
EF is an even flow penalty function which quantifies the degree of output fluctuation during the planning period. It is defined as follows: , ( Where M is the planning horizon (the number of years in the planning period); V j is the timber yield in each year and is the average yearly timber yield during the entire planning period.
FSV is a forest spatial criterion which can be used to achieve spatial aggregation or dispersion of forestry operations during the planning period. Aggregation of operations is desirable to reduce the cost of moving harvesting machines and equipment during a given point in time. Dispersion may be desirable to achieve conservation objectives (Öhman, 2000). The FSV-value of a real forest is rather more complicated than that of a theoretical forest where the land parcels may be arranged in some regular fashion (Chen & Gadow, 2002). We define the variable FSV for a forest with known geographical details of the compartments as follows: , Where N is number of land parcels (commonly known as stands); NV i is the neighbor value for stand i; R ik is a value which quantifies the relationship between two treatment options of two neighbouring stands i and k; L ik is the common boundary length (CBL) of stands i and k and D ik is the distance between the two centroids of stands i and k. The variable R ik is used to encourage or discourage pairs of neighboring stands to have the same, similar, or different management paths. R ik defines a value for management path i and management path k, whenever i and k refer to a pair of neighboring stands. If a spatial concentration of management activities at a given point in time is desired, the pair of adjacent stands will contribute a positive value to the objective function (R ik >0), otherwise R ik is assigned a negative value. The magnitude of the positive or negative contribution is determined by the value of R ik . The entire matrix of R ik among all management paths is determined beforehand. A detailed description of the variable R ik and the method of aggregation or dispersion of management activities in the landscape may be found in Chen & Gadow (2002).
When considering the spatial arrangement of stands with alternative management paths it is useful to calculate the common boundary length (CBL) of every pair of polygons. The CBL will provide not only the adjacent relation between two neighboring stands, i.e. connected (CBL=0) or not connected (CBL>0), but also a quantitative value of adjacency. This quantitative value can be used to evaluate the treatment options in a group of neighboring stands.
We illustrate the approach using a small example called Winnefeld south. Winnefeld south is part of the Winnefeld forest district in the Solling region in Northern Germany, comprising 41 stands of different area, age, site index, dominant height, stem number and basal area. A typical forest map in ArcView GIS can provide of the required geographic information about the forest stands, such as the area and the centroid coordinates. However, it is not easy to directly obtain the CBL. The first step is to generate a unique ID for each parcel in the forest. One way to do this is to generate a unique link string with a format which guarantees that there are no duplicate polygon labels (this requirement is not as trivial as it may appear).
It seems that there is no simple method in ArcView for calculating the CBL of every two neighboring polygons. The Arc/INFO system is better suited to such cal-culations. There are several methods that can be used to convert a shape file into an ArcInfo coverage. The forest map in Figure 1 with 41 polygons originally specified in the ArcView shape format has been converted into an Arc/INFO coverage 1 . A coverage of the Arc/INFO topology usually has two attributes, the polygon and the polyline. The latter one includes the required CBL information. The polyline properties of an Arc/INFO coverage are stored in it's a so-called AAT 2 file, the polygon properties in the corresponding PAT 3 file.
The distance between two polygons can be defined as the distance between their centroids. The Centroid of a polygon can be directly calculated in the ArcView system using a standard function. Then the Euclidian distance can be calculated between two centroids. An example of part of the required data file is presented in  Polyline number (small digit) and centroid position (#). There are 108 polylines in this coverage. The CBL is derived from the polyline length. The distance between every two centroids is calculated using the Euclidian distance formula.

Growth and Thinning Models
The models used for estimating Spruce growth and yield were developed on the basis of a comprehensive data set from permanent growth trials maintained by the Nordwestdeutsche Forstliche Versuchsanstalt (Schübeler, 1997;Sánchez et al., 2001;Vilčko, 2002). The system includes equations for dominant height growth, tree mortality (maximum density), basal area growth, and a stand volume equation.
The thinning options were defined by a specific weight and a thinning type. The thinning weight specifies how much of the available basal area is removed. Based on proposals made by Spellmann et al. (1999), rules referring to the thinning weight were developed for three thinning phases. The thinning phases are defined by certain top height intervals. The rules defining the thinning weights in each top height intervals are based on the assumption that G max the maximum possible basal area is known. The thinning rules can be summarized as follows: A thinning is possible when the stand top height > 14 m and G maxent > 0; basal area removed in each of three top height intervals (14-20 m; 20-26 m; > 26 m); We define W to represent the percentage of the maximum allowable thinning in terms of basal area removed (G maxent ). Thus, "W=100" means that 100% of the allowed basal area is removed during a thinning. The planning period has been set to 20 years, the rate of interest to 4%, the timber price is assumed to be 100 €/m 3 . 21 alternative management paths were specified for all stands using the approach presented by Chen & Gadow (2002). The 1 st path is always a "do nothing" option. Other option types include: 1 thinning, 2 thinnings, 3 thinnings and 4 thinnings during the planning period. The same option ID for different stands represents the same specific thinning regime and timing of thinnings. This is essential to ensure aggregation of operations in neighboring stands in the same year.

Optimisation
In this section we first present the results of the optimization for the Winnefeld South forest and discuss the effects of introducing different objective function components. This analysis is followed by a similar approach involving a bigger forest area.

A forest with 41 compartments
With 21 alternative management paths for each stand, and 41 stand polygons, the total number of possible path combinations is equal to 21 41 . Simulated annealing (SA) can be used to search for the best combination in terms of the objective function. A specific path number refers to a specific thinning regime which specifies the timing of harvesting operations. For example, no thinnings will take place for all management paths labeled as "1". The label "2" specifies a specific sequence of thinnings scheduled for year 1, 6, and 11 during the planning period. The operations labeled as path number "3" are scheduled to take place in year 2, 7 and 12. This method of assigning management activities facilitates the spatial aggregation of harvesting operations. Depending on the objective function, neighboring stands may be encouraged (aggregation) or discouraged (dispersion) to carry the same path. In our current model formulation, all relationships between two management paths of two neighboring stands are specified as follows: R ik =1, if stand i and stand k are assigned the same management path R ik =0, if stand i and stand k are assigned different management paths Figure 2 presents an example of a simulated annealing process. All objective function components converged after less than 100 iterations. The processing time for this problem was less than 1 minute using a P3-500 /Window98/VB6 environment.
The result of the optimization may be expressed in terms of annual operational plans. The operational schedules for the first 5 years are presented in Figure 3. The maps in Figure 3 indicate the (shaded) compartments where thinnings are scheduled to take place in the different years. The maps in time thus facilitate the interpretation of the spatial effect of adding specific components to the objective function. Thinnings were not scheduled in period two when only EC featured in the objective function (Figure 3, first row, second map). There are no timber harvests during the second year in the planning period and the timber output is rather unbalanced. Adding the even flow component results in a more balanced output of timber over time. Thus, a simplified objective function may cause a high fluctuation of timber yields and thinning areas during the years, which is often undesirable.  Figure 4. If the economic criterion is the only component in the objective function, the value of EC is the highest, as expected, while the other components are ignored. Adding the even-flow component improves the corresponding desirable property of the harvest schedule while decreasing the economic component value to some extent. When compared with the results of the one-component objective function, the yearly operational areas of the two-component objective function show some similarities (Figure 3: 2nd row, Figure 4).
When the spatial component is added into the objective function, the spatial aggregation of thinnings improves notably (Figure 3: maps in 3rd row). The economic component decreased somewhat while the yield fluctuation showed only a slight increase (Figure 4).
The optimization for multiple components is a process where different objectives compete with each other. The choice of the weight coefficients is obviously important to balance the relative importance among the different objective function components. Different objectives produce different optimization directions. The solution of a multiple objective problem will be decided by the balance among the responding coefficients (in our case by the constants a, b and c). In our example, the final values of a, b, and c where established using sensitivity analysis. We were not able to determine a priori values of the weighting constants. To realize spatial aggregation, it is necessary to assign a high relation value when the paths of neighboring stands are the same. However, specifying the same option for neighboring stands is not enough. We must also make sure that neighboring stands will be harvested (thinned) in the same time period. To realize periodic spatial clustering of operations in neighboring stands, it is essential that the same option identifier has the same meaning in different stands, i.e. the same treatment in a single year or the same succession of operations in successive periods (e.g. thinning in periods 1, 6 and 11).

A bigger case
A similar approach was followed using a bigger forest with more compartments. The forest, called "Winnefeld 2", includes 739 compartments. The total area is 3182.2 ha and the compartment map is presented in Figure 5.
The same procedure concerning definition of management paths was followed. the timber price, interest rate and growth and thinning models were the same as in the previous example. The Simulated Annealing method was also used to optimise the three-component objective function. The map in Figure 6 shows stands with same color representing the same treatment option. The fact that some groups of neighboring stands are assigned the same option indicates that the desired aggregation of harvesting operations is realized in the same year. White polygons with dots are assigned the "Do nothing" option. Some of these stands are too young for thinning during the next five years.
The processing time for this problem was less than 30 minutes in a P3-500 / Window98/VB6 environment. As expected, due to the large number of stands, the fluctuation of periodic timber yields can be reduced to almost zero ( Figure 6).
The result of changing the number of objective function components is similar to that shown in Figure 4. As the economic, even flow and spatial components are added into the objective function, the operational schedule approaches more and more the desirable one. The bigger example shows that the optimization process is stable, producing satisfactory schedules for different size problems.

Discussion and Conclusions
The spatial modeling approach followed in this analysis is an important part of the Multiple Path theory of forest design. Our specific application utilizes geographic relations among neighboring stands, including the distance between stand centroids and the length of the common boundary between adjoining stands. This approach appears to be more realistic, than previous attempts in spatial optimization.
In previous studies, the shape of compartments has often been simplified, representing rectangles or other regular shapes. Yoshimoto & Brodie (1993) generated spatial constraints using simplified forests consisting of stands of a rectangular and ladder-like shape (see also Jones et al., 1991). Hof & Joyce (1992; and Hof et al. (1994) simulated the spatial allocation of forest stands as a cellular grid. Murray (1999) applied spatial restrictions in harvest scheduling with a map of irregular spatial units. Carter et al. (1997) and Yoshimoto (2001) used a real map in their approach of modeling wildlife habitat and timber production. Van Deusen (2001) optimized the spatial problem by defining forests and lakes as rectangles. All the above authors consider only the quality of adjacency relations (connected or not) with a simulated map of squares. Irregular polygon shapes were used by Barrett & Gilless (2000) who applied the "4-color" theorem with a restriction of the maximum size of even-aged harvest units. The distances between the units are not included in the calculations. Chen and Gadow (2002) also used a grid map in their approach, but considered the distances among neighboring stands. Irregular polygon shapes were accepted by Öhman & Erikson (1998) who considered the edge widths for calculating the core area using accurate geographic information. They used 1 and 2 pixels to represent 32 and 64 m to respectively with a map which closely resembles a real forest map.
In these applications, the interpretation of spatial modeling is quite different to ours. We refer to aggregation or dispersion of harvest operations in time, while some authors consider the road system or wildlife habitat and forage area in their spatial Figure 5. Map of the "Winnefeld 2" example. Stands with the same color represent the same treatmvent path. Compartments with dots were assigned the "do nothing" option during the 5 year planning period.
models. Bettinger et al. (1997; use an objective function to evaluate the NPV of road maintenance costs.
The main focus of this paper is the quantitative analysis of neighbourhood relations for a real forest using geographical information about individual stands with irregular boundaries, derived from a GIS system. As far as we know, this is a new approach. Economics, even flow and specific temporal-spatial requirements are all important forest management goals. The examples show that it is possible, and quite revealing, to aggregate silvicultural treatments in space and time and to evaluate different combinations of objective function components. The method of Simulated Annealing, can be a useful tool when searching for a satisfactory operational schedule.
The greater the number of management paths specified for each individual stand, and the greater the number of stands, the greater is the chance that a fine-tuned global optimum will be obtained. However, a higher number of possible combinations may also considerably increase the processing time (Hinrichs, 2006;Seo et al., 2005). We found that about 21 options per stand represents a reasonable compromise and that processing times even for large forests with up to 1500 compartments are satisfactory using an ordinary personal computer.