Spatial Forest Planning for Optimized Harvest Scheduling

ABSTRACT

The aim of this study was to develop a mathematical model for the spatial forest planning of harvest activity scheduling. Thirty-eight (38) blocks were formed from stands aged 13 to 17 years considering a total area limit set to 350 hectares. The model was built in Excel ^® spreadsheet and solved by CPLEX^® optimization software using the Branch and bound method for the linear programming of problems with integer variables. The results of optimization are presented in the form of a map of spatial distribution of the selected blocks. The total harvested area was 5,989 hectares and the volume generated was 2,956,913 m³ in a planning horizon of five years. Higher results compared with those of traditional planning and efficiency to optimize forest harvesting proved to be a viable alternative to the inclusion of space and operational issues associated with forest harvesting.

Keywords:
geographic information system; tactical planning; linear programming

1. INTRODUCTION

A characteristic of spatial forest planning is the generation of problems with many constraints and integer decision variables. Several attempts have been made to represent spatial constraints within formulations of planning problems.

Linear Programming (LP) is one of the classic optimization techniques, and is one of the most widely used among the techniques developed by operational research. It is often used in forest management, especially in forest planning for timber harvest scheduling ( Ezquerro et al., 2016 ).

It is truly challenging to express problems through mathematical models involving the size, shape, and distribution of management units (e.g., stands, harvest blocks, wildlife habitat, and age class), minimum and maximum harvest limits, adjacency constraints, connectivity, and road proximity ( Baskent & Keles, 2005 ).

For Weintraub & Murray (2006) , at the tactical level, the decisions involve the planning period, during which each unit will be harvested in order to meet expected future demands, and demands a definition of access roads to harvesting units. This adds spatial characteristics to the planning process and requires linear programming models with binary variables.

Grouping of harvesting activities into operational blocks aims at the concentration and connectivity of the forest areas, so as to reduce the movement of the harvesting machines between the stands selected for clear cutting in the same year as the planning horizon, as well as efficiency in the use of forest roads. Notably, efficiency of harvesting activities involves the organization of machinery and log transport, and is significantly improved when stands are aggregated ( Smaltschinski et al., 2012 ).

Alternatives to include area connectivity requirements in forest planning models are always useful, given the scant of research addressing the aggregation of harvesting activities. In this sense, Smaltschinski et al. (2012) applied a cluster analysis to aggregate the selected stands for harvest by applying the Dijkstra minimum path algorithm to select the stands.

The present work had the objective to develop a mathematical model for spatial forestry planning aiming at the scheduling of harvesting activities. Additionally, stands were grouped into harvest blocks in order to form more concentrated and continuous areas, and to minimize the average transport distance of wood.

2. MATERIAL AND METHODS

2.1. Study area

The study was developed in a forest farm located in the municipality of Sengés, northeastern Paraná state, Brazil. Climate in the region is Cfb according to the climatic classification of Köeppen, i.e., it is a temperate climate with average temperature in the coldest month below 18 °C (mesothermic), with frost, fresh summers with average temperature in the warmer months below 22 °C, and without a defined dry season.

2.2. Data specification

Figure 1 shows the initial spatial distribution of the forest areas classified by the planting year together with the road network of the study area. This information was obtained by a GIS database.

Figure 1
Characterization of the study area located in the municipality of Sengés, Paraná state, Brazil.

The area is composed of 1,832 stands, of which nine stands are of the species Eucalyptus sp., five stands of Pinus sp., and 1,818 stands of Pinus taeda, with ​​approximately 13,510 hectares of production forest and 22,871 hectares of total area. The average area of ​​the stands is 7.25 ha and the average number of adjacencies per stand is four. Table 1 shows the distributions of the forest harvest area by age class, considering the frequency of Pinus and Eucalyptus stands.

Table 1
Initial structure of the forest for harvest planning.

Age Classes	Area				Stands
(years)	(ha)		(%)		Frequency		(%)
[0 to 5)	2,833	21.43		504		27.51
[5 to 10)	3,407	25.78		583		31.82
[10 to 15)	5,815	43.99		604		32.97
[15 to 20)	1,148	8.68		130		7.10
[20 to 25)	2	0.01		2		0.11
> 25	14	0.11		9		0.49
Total	13,219	100		1832		100

Annual volumes by assortment in each management unit were obtained from a production table considering a site index of 23m at the reference age of 15 years.

2.3. Generating harvest blocks

Aiming to prevent stands from being selected for cutting in the same period that they were very distant from each other or difficult to access, some criteria were considered in forming the blocks, such as presence of the road network, as well as of the hydrographic network, preservation areas, and management units. The inclusion of these criteria ensured a priori formation of continuous blocks and accessibility, because there is no internal displacement impediment of the cutting front or for transporting the wood.

The blocks were composed of stands aged 13 to 17 years within a planning horizon of five years. In order to ensure that the blocks contained the same data as the original forest cover (area and age), these layers were combined using the Geographic Information System (GIS).

The real road distance between all pairs of blocks was measured using the GIS. The respective bifurcations of the main road with the secondary road to where the wood is forwarded were considered as the points of origin and end of the route, thereby enabling the construction of a matrix of main road distances between blocks.

2.4. Problem description

The main focus of this work was to schedule harvest in order to minimize the average distance of timber transportation, respecting the demands of the products for each year and the management alternatives established initially. The scenario must meet the demand of three multi-products, classified as process, sawmill 1, and sawmill 2 ( Table 2 ).

Table 2
Assortment classification.

Description	Diameter (cm)	Length (m)	Total volume (%)
Process	8 to 18	2.65	0.52
Sawmill 1	18 to 25	2.65	0.25
Sawmill 2	>25	2.65	0.19

A definition of the annual volumetric production goal was proposed in relation to the demand for each assortment class obtained with the control planning. Figure 2 displays the production goals resulting from the tactical plan provided for the study area for the five annual periods considered and for the three different products. Based on the results of a tactical harvesting plan, the annual production targets for the study area were determined in m³ per product.

Figure 2
Volumetric demand in thousands of m³ according to assortment class and period.

2.5. Suggested model

For this study, the problem was represented by a model according to the integer linear programming (ILP) structure, following the model I formulation of Johnson & Scheuerman (1977) for the annual sequencing of the harvesting operations of the blocks selected for clear cutting in tactical planning. The configuration of the harvest problem addressed requires the use of binary variables for yes or no decisions. The indices used in the model are:

• i = number of the block to be harvested, i = 1, ..., 38;

• j = annual planning period, j = 1, ..., 5;

• p = multi-product or assortment number, p = 1,2,3;

• Nb = total number of blocks;

• Np = total number of years of the planning horizon;

• V_ijp = total volume (m³) of product p existing in block i in year j;

• K_pj = demand for product p in year j ;

• D_i1i2 = distance (m) between blocks i1 and i 2;

• X_ij = binary variable: 1, if block i is harvested in year j; 0, otherwise;

• Y_i1i2j = binary variable: 1, if the stretch of road between the pair of blocks i = 1 and i = 2 is used in year j ; 0, otherwise.

Objective Function:

$\text{Minimize Z}={\displaystyle \sum }_{{\text{i}}_1}^{\text{Nb}}{\displaystyle \sum }_{{\text{i}}_{2\ne }{\text{i}}_1}^{\text{Nb}}{\displaystyle \sum }_{\text{j}=1}^{\text{Np}}\left({\text{Y}}_{{\text{i}}_1{\text{i}}_2\text{j}}{\text{*D}}_{{\text{i}}_1{\text{i}}_2}\right)$ $\left(1\right)$

Subject to ${\displaystyle \sum }_{\text{i}=1}^{\text{Nb}}{X}_{\text{ij}}{\text{V}}_{\text{ijp}}\ge {\text{K}}_{\text{pj}}\forall \text{i},\forall \text{p}$ $\left(2\right)$

${\displaystyle \sum }_{\text{j}=1}^{\text{Np}}{X}_{\text{ij}}\le 1\forall \text{i}$ $\left(3\right)$

${\text{X}}_{{\text{i}}_1\text{j}}-{\text{Y}}_{{\text{i}}_1{\text{i}}_2\text{j}}\le 0$ (4)

${\text{X}}_{{\text{i}}_2\text{j}}-{\text{Y}}_{{\text{i}}_1{\text{i}}_2\text{j}}\le 0$ (5)

${\text{X}}_{\text{ij}},{\text{Y}}_{{\text{i}}_1{\text{i}}_2\text{j}}\in \left\{\mathrm{0,1}\right\}\forall \text{i},\forall \text{j}$ $\left(6\right)$

The formulated optimization model has an objective function of minimizing the total distance of main roads for timber production flow. This total distance is obtained by summing the distances between the block pairs that will be harvested in a given year (1). For each year, each of the products must meet a defined demand. A restriction that ensures the supply of products is represented by expression (2). A cut-off constraint (3) ensures that each block is harvested in its entirety and not again in the planning horizon, and defines decision variables in binary form, which force a choice of a single decision variable (management alternative). Constraints (4) and (5) allow a combination of harvesting block pairs with the distance of the road stretch between them. In addition, expression (6) guarantees that variables ${\text{X}}_{\text{ij}}$ and ${\text{Y}}_{{\text{i}}_1{\text{i}}_2\text{j}}$ only assume binary values.

2.6. Implemented hardware and software

The integer linear programming (ILP) mathematical model was built in a Microsoft Excel ^® 2013 spreadsheet and was solved using the IBM ILOG CPLEX^® Optimization Studio 12.6.1 optimization software with an IBM Academic Initiative license. Data were processed in a personal computer with an Intel Core^® i7 2Ghz processor with 12 GB of RAM and 1 TB of HD. The distances, arrays, and maps were generated by the ArcGIS 10.3 platform.

3. RESULTS AND DISCUSSION

Thirty-eight (38) blocks were created from the map with the provided geographic information, with adjacent stands and ages varying between 13 and 17 years, respecting the maximum limit of contiguous area of ​​350 ha. Figure 3 depicts the blocks and roads used to transport the forest production in the study area.

Figure 3
Spatial distribution of the blocks.

The area to be harvested was grouped into blocks with total area varying between 100 and 350 ha. Thus, if a stand has an area of ​​20 ha for example, it is added to the neighboring (adjacent) stands until it reaches an area close to the desired limits. Adopting maximum and minimum cutting ages made it possible to meet the requirements of the company, excluding young stands with low harvest stock. The concept of minimum cutting age is commonly used in spatial planning, as seen in the studies by Nelson (2001) and Binoti et al. (2014) .

The distances between each block pair obtained with the geographic information system were represented by a matrix, which is summarily presented in Table 3 . These distances were used in the objective function.

Table 3
Representation of the matrix BxB (km) elements for the first 10 blocks.

Block/Block	1	2	3	4	5	6	7	8	9	10
1	–	0.7	0.9	5.1	1.6	13.2	15.8	14.3	18.3	26.5
2	0.7	–	1.4	6.0	7.6	6.3	4.7	2.3	6.3	14.4
3	0.9	1.4	–	1.4	3.1	12.0	10.3	7.9	11.9	20.0
4	5.1	6.0	1.4	–	3.1	13.4	11.7	9.4	11.9	21.5
5	1.6	7.6	3.1	3.1	–	11.8	13.4	11.0	15.0	23.1
6	13.2	6.3	12.0	13.4	11.8	–	0.5	3.2	4.2	11.2
7	15.8	4.7	10.3	11.7	13.4	0.5	–	1.9	2.5	9.6
8	14.3	2.3	7.9	9.4	11.0	3.2	1.9	–	0.6	7.7
9	18.3	6.3	11.9	11.9	15.0	4.2	2.5	0.6	–	6.8
10	26.5	14.4	20.0	21.5	23.1	11.2	9.6	7.7	6.8	–

Formulation of the planning model resulted in 7,410 binary decision variables and 14,494 constraints, of which 38 are cut restrictions, 15 are demand constraints, and 14,440 are distance constraints. After more than nine hours of processing, the CPLEX^® software could not find an optimal solution, which led to some modifications to the model accelerate the resolution.

The model was dismembered into five annual models with the intention of minimizing the total transport distance of the objective function considering one year at a time, and beginning in year 1. Restrictions were added to the model from year 2 onwards in order to keep the optimal obtained solutions unchanged from previous years. In this way, each block selected by the model to be harvested in year 1 became a constraint for the model of year 2, and so on. The choice of optimal solution for each year is made with decreasing priority over the planning horizon.

With the dismemberment of the main model into five annual models, resolution times were significantly reduced, and varied from 0.25 to 1.5 seconds, as shown in Table 4 . Table 4 also presents the number of interactions and the GAP, which represents the percentage difference between the optimal value obtained and the best boundary.

Table 4
Summary of the results from the CPLEX^® software.

Year	Interactions	Solution Time (s)	GAP (%)
1	11,066	1.09	1.93
2	28,296	1.50	1.16
3	8,784	0.56	3.28
4	3,986	0.36	0.34
5	870	0.25	0.00

Optimization results are presented as an easy-to-understand map showing the spatial distribution of harvesting units over the five years of HP. Figure 4 illustrates the harvest scenario after optimizing the final model for each year of the planning horizon.

Figure 4
Spatial harvest scheduling for the 5-year planning horizon.

As shown in Figure 4 , the layout of the blocks selected for harvesting is considered satisfactory, both for proximity to the roads as well as between the blocks, expressing the relation of the access roads to each of the collection blocks.

Although the presented formulation does not guarantee adjacency of the harvest areas, it generated good quality responses, showing a viable alternative for the inclusion of spatial and operational issues associated with the forest harvest, and considering that this block arrangement meets the requirements of the areas that are ideal for cutting.

Figure 4 shows a solution with the blocks to be collected in each period, and which are presented in Table 5 together with the harvested area, the total volume obtained, and the total distance used to harvest the group of blocks.

Table 5
Selected blocks for each year of the planning horizon.

	Selected Blocks
	Year 1	Year 2	Year 3	Year 4	Year 5
	6	8	9	1	4
	14	10	13	5	11
	21	16	22	12	17
	25	26	24	23	19
	27	31	29	-	-
	30	-	-	-	-
Area (ha)	1,597.17	1,285.94	1,153.79	956.92	994.92
Volume (m3)	683,720.1	587,449.7	564,213.6	560,964.8	560,564.7
Road (km)	29.5	32.5	26.9	33.9	27.2

Regarding the area values ​​in Table 5 , by adding all the years of the planning horizon, a surface of 5,989 ha would be harvested, resulting in a total volume of 2,956,912.9 m³.

According to Augustynczik et al. (2016) , by aggregating the forest harvesting areas, it is possible to increase the efficiency of the operation due to the reduced movement of the machines between the harvest fronts, thereby reducing the number of unproductive hours. In addition, maintenance and road construction costs are reduced, because the production of a large number of stands can be transported on the same road. The savings in transportation time obtained by the proximity to the roads also brings a considerable reduction in the operating costs.

Results obtained from the employed methodology were compared with the control data of the same study area provided by the company, resulting from a traditional harvesting plan in which the management units are the stands, without considering the proximity between them and the range of 13 to 17 years for harvesting ( Figure 5 ).

Figure 5
Traditional annual planning.

A typical empirical planning is illustrated in Figure 5 , where very young stands were harvested aiming to reduce the displacement of the harvest front. This can result in low wood stock and thin logs. It is interesting to observe that the control data also showed good results in years 1 and 4, but it should be emphasized that this only occurred because the stand ages harvested in the control were aged 10, 11, and 12 years, which is below the age practiced by the company. Mathematical models can eliminate these stands by means of imposed constraints.

It can also be observed that in years 2, 3, and 5, the control data harvested younger trees and was still surpassed by the optimized scenario. In relation to the total harvested area in the five years of the planning horizon, the company obtained 5,496.8 hectares and a volume of 2,922,201 m³.

In addition to the area and volume harvested by blocks that are superior to traditional planning, the greatest advantage of the methodology presented herein was the proximity of the stands to be harvested in the same period grouped in blocks, as this proximity should reduce transportation costs, road maintenance, and personnel displacement. Thus, the initial requirements were met by harvesting more concentrated and contiguous areas, with access through main roads and optimal cutting ages.

According to Ohman & Eriksson (2010) , another effect of aggregating harvesting areas is the reduced application of thinning in the forest. This is due to application of clear cutting, because with aggregation of the harvest areas it is possible to obtain a larger volume from smaller harvest areas.

Problems of forest harvesting optimization involving roads are complex due to the large number of restrictions in forming the models. A good example was demonstrated by Andalaft et al. (2002) , which considered horizons of two to five years and formulated MILP models with more than 6,000 restrictions relating the roads to the forest harvest. Forming the blocks prior to constructing the model decreases the number of restrictions in the model as well as the computational effort.

4. CONCLUSION

This study proposes an effective approach to scheduling harvest considering the imposed production requirements. Flexibility in harvesting ages occurred only once during the conversion of the current block structure, considering that the blocks should last over time after the optimized scenario is implemented, and that management decisions are maintained.

The CPLEX^® software licensed by the IBM Academic Initiative to solve the problem is very efficient and easy to handle and understand. Forest spatial data can be used in tactical forest planning in block formation from adjacent stands to produce results similar, and even superior, to those normally obtained by traditional non-spatial planning scenarios.

REFERENCES

• Andalaft N, Andalaft P, Guignard M, Magendzo A, Wainer A, Weintraub A. A problem of forest harvesting and road building solved through model strengthening and lagrangean relaxation. Operations Research 2002; 51(4): 613-628.
• Augustynczik ALD, Arce JE, Silva AC. Aggregating forest harvesting activities in forest plantations through integer linear programming and goal programming. Journal of Forest Economics 2016; 24: 72-81.
• Baskent EZ, Keles S. Spatial forest planning: A review. Ecological Modelling 2005; 188: 145-173.
• Binoti DHB, Silva Binoti MLM, Leite HG, Gleriani JM, Oliveira MLR, Campos JCC et al. Regulação de florestas equiâneas com restrição de classes de idade por compartimento. Cerne 2014; 20(4): 567-573.
• Ezquerro M, Pardos M, Diaz-Balteiro L. Operational research techniques used for addressing biodiversity objectives into forest management: an overview. Forests 2016; 7(10): 1-18.
• Johnson KN, Scheuerman HL. Techniques for prescribing optimal timber harvest and investment under different objectives: discussion and synthesis. Forest Science 1977; 23: a0001-z0001.
• Nelson JD. Assessment of forest harvest blocks generated from operational polygons and forest cover polygons in tactical and strategic planning. Canadian Journal of Forest Research 2001; 31: 682-693.
• Smaltschinski T, Seeling U, Becker G. Clustering forest harvest stands on spatial networks for optimized harvest scheduling. Annals of Forest Science 2012; 69: 651-657.
• Ohman K, Eriksson LO. Aggregating harvest activities in long term forest planning by minimizing harvest area perimeters. Silva Fennica 2010; 44(1): 77-89.
• Weintraub A, Murray AT. Review of combinatorial problems induced by spatial forest harvesting planning. Discrete Applied Mathematics 2006; 154(5): 867-879.

Publication Dates

• Publication in this collection
  07 Feb 2019
• Date of issue
  2019

History

• Received
  15 Mar 2016
• Accepted
  07 Dec 2017