IMPLEMENTING MINIMUM AREA HARVESTING BLOCKS IN AN OPTIMIZED FOREST PLANNING MODEL 1

1 Received on 16.02.2016 accepted for publication on 16.03.2017. 2 Universidade Federal do Paraná, Programa de Pós Gradução em Engenharia Florestal, Curitiba, Paraná Brasil.. E-mail: <andreylessa@gmail.com>. 3 Universidade Federal do Paraná, Departamento de Ciências Florestais, Curitiba, Paraná Brasil. E-mail: <jarce@ufpr.br>. 4 Universidade Federal do Paraná, Departamento de Matemática, Curitiba, Paraná Brasil. E-mail: <arineicls@gmail.com>. *Corresponding author.


INTRODUCTION
The inclusion of spatial relationships has become an essential aspect of forest planning during the past decades, especially the issues related to landscape patterns resulting from forest harvesting activities (Constantino et al., 2008).Optimization models aiming to include these spatial relationships into forest planning have been proposed in the literature, notably through adjacency constraints.The focus of spatial planning problems has been the inclusion of maximum clearcut area constraints and developing efficient solution methods for these problems (Öhman and Lämas, 2003).
Problems considering the connectivity of forest areas, however, have been gaining importance in forest planning, such as the creation of ecological corridors and continuous core areas, aiming to ensure the maintenance of biodiversity and other ecosystem services (e.g.water protection and scenic beauty).Numerous studies addressed problems related to the creation of ecological corridors and core areas for biodiversity ( Rebain and Mcdill, 2003;Tóth et al., 2006;Wei and Hoganson, 2007;Carvajal et al., 2013).
While the creation of core areas aims to set contiguous stands aside from production, other important class of connectivity problems targets at clustering harvesting areas.These models also aim at protecting biodiversity, once spreading harvesting activities may reduce the habitat availability for interior forest species and enhance negative edge effects (Gustafson, 1996).Different approaches were proposed in the literature for dealing with clustering of harvesting activities, such as dynamic zoning of harvesting areas and multiobjective models for aggregating stands scheduled for harvesting (Gustafson, 1998;Öhman and Lämas, 2003;Öhman and Eriksson, 2010;Smaltschinski et al., 2012).
Problems involving clustering of harvesting areas and creation of core areas are closely related, allowing the use of similar approaches for its solution.Nevertheless, in the case of clustering harvesting activities the problem can become more complex, as different contiguous harvesting blocks must be formed at each period.On the other hand, core areas are usually maintained during the whole planning horizon (PH).Therewith, the clustering of harvesting areas may demand more processing time.
The objective of this study was to include adjacency constraints into an optimized forest planning model, aiming to maintain interior forest habitat and reduce edge effects.To this end, we applied two approaches, proposed initially for the creation of core areas, generating Integer Linear Programming models.The first approach was based on the dynamic connectivity model proposed by Carvajal et al. (2013) and the second was based in the old-growth forest formulation from Rebain and McDill (2003).We evaluated the efficiency of both formulations for creating harvesting blocks and their impacts on the objective function.

Study area
Our research area was a planted forest with Pinus taeda, Pinus elliottii and Eucalyptus dunnii, located in the municipalities of Bituruna and General Carneiro,Brazil,between the coordinates 26°13'58.31" and 26°22'5.634" S and 51°34'14.6 and 51°30'26.14"W. The climate according to Köppen classification is Cfb: mesotermic humid subtropical, with mild summer and occurrence of frequent frosts during winter.
The area was composed by 236 stands, with 18 Eucalyptus dunnii stands, 21 Pinus elliottii stands and 197 Pinus taeda stands, with a total area of 2365.8 ha, average stand area of 10.3 ha and the average number of adjacencies per stand was 3 (Figure 1).

Optimization model
The first formulation applied for ensuring the minimum contiguous harvesting area was the approach proposed by Carvajal et al. (2013), so-called ring inequalities, described in Eq. (1): Where: x ij : binary values that take value 1 case stand i is harvested in year j or 0 otherwise; n c : number of stands in set C minus one; : block of stands that violate the minimum area requirement; C (A min ): all sets of stand blocks with area inferior to the established limit; : set of stands in the neighborhood of C.
Eq. ( 1) requires that in case a block of stands is selected for harvesting in a defined year of the planning horizon, and the sum of their areas is inferior to the (1) ( min), Implementing minimum area harvesting... established minimum area, at least one stand adjacent to the considered block and not part of the block is selected for harvesting in the same year.This constraint is applied for all sets of stands that violate the minimum area limit.Therewith, the optimized forest planning model considering minimum harvesting block areas can be defined as follows: Where: N: number of stands; NP number of years in the planning horizon; c ij : NPV generated by stand i in case it is harvested in year j; x ij : binary variable that takes value 1 case stand is harvested in year j or 0 otherwise; v ij : volume generated by stand i in year j; A min : minimum area limit; ; C: block of stands that violate the minimum area requirement; C(A min ): all sets of stand blocks with area inferior to the established limit; : set of stands in the neighborhood of C.
The objective function of the optimization model is represented in Eq. ( 2), maximizing the total forest total NPV, from the sum of the NPVs generated by each stand under a certain management regime.Constraint Eq. ( 3) imposes that stands are harvested only once during the planning horizon.Constraints Eq. ( 4) and Eq. ( 5) are the wood flow constraints, ensuring that the harvested volume during each period of the planning horizon respects the established limits.In our study, the allowed fluctuation for the harvested volume at each period of the PH was 10% in comparison to the volume harvested in the first year.Constraint ( 6) is the minimum area set of constraints, ensuring that the area of stands scheduled for harvesting are greater or equal the established limit.We tested different minimum harvesting area limits, creating scenarios with minimum area constraints of 30, 40, 50 and 60 ha.Constraint Eq. ( 7) requires that the decision variables x ij take binary values.
The second approach, based on the model from Rebain and McDill (2003), creates harvesting blocks a priori applying a minimum and maximum block area threshold and ensures that all stands selected for harvesting in a defined year of the PH are part of at least one of these blocks: (5) Where x ij : bynary variable that takes value 1 case stand is harvested in year or 0 otherwise; n k : number of stands in block ; : binary variable that takes value 1 case block k is selected for harvesting in year j; C k : set of stands in block k; S: set of all blocks of stands formed applying the minimum and maximum area limits; P i : set of harvesting blocks containing stand i.
Constraint Eq. ( 8) guarantees that in case block k is selected for harvesting in a defined year of the planning horizon, all stands in this block are selected for harvesting at the same year.Constraint (9) aims to ensure that case stand i is selected for harvesting in a defined year of the planning horizon, a harvesting block containing stand i is also selected for harvesting.For the formulation based on Rebain and McDill (2003) the set of minimum area constraints Eq. ( 6) was replaced by the set of constraints Eq. ( 8) and Eq. ( 9), with j varying from 1 to 5, the objective function and other constraints remained unchanged.We highlight that this set of constraints allows harvesting blocks to overlap, thus generating contiguous harvesting areas greater than the limits applied for creating the blocks.The adopted planning horizon was 16 years, applying management regimes without thinning.The wood flow constraints allowed a fluctuation of 10% for more or less, compared to the wood production in the first year of the PH.For producing harvesting blocks according to the formulation from Rebain and McDill (2003), we applied ranges of 10 and 20 ha for the harvesting block area.Hence, we obtained 13 scenarios (Table 1).The planning scenarios were generated through the software Optimber-LP and the optimization models were solved through Gurobi 5.5 in a PC with Intel® Core™ Duo CPU 2.93 GHz processor and 4Mb of memory.The optimization time limit was established as 90 minutes.

RESULTS
Figure 2 shows the forest harvest scheduling for selected scenarios (1, 5 and 13).The stands selected for harvesting in years in which the minimum area constraints were applied appear highlighted.Implementing minimum area harvesting...
Considering scenario 1 in Figure 2, without inclusion of minimum area constraints, we perceived a dispersion of harvesting areas, especially in year 5 of the planning horizon.In the remaining years isolated harvesting stands occurred, but harvesting blocks were generated as well.The inclusion of the minimum area constraints leads to a reorganization in the harvest scheduling for all scenarios.Notably, it is possible to see the creation of harvesting blocks with area greater than 60 ha in scenarios 5 and 13, especially in year 5, in the southern part of the study area.
In all scenarios, the minimum area limit was respected and the largest changes occurred in 5 th year.The stands scheduled for harvesting during the 1 st and 2 nd years showed alteration in the block configuration, but most of harvesting areas remained unchanged.
Regarding the harvesting block formation, all models showed satisfactory results, being capable of creating harvesting blocks respecting the established thresholds, as exemplified by scenarios 5 and 13 in Figure 2. Therewith, there was no major differences between the two formulations regarding the compliance with a minimum area limit.Therefore, aiming at evaluating the formulation performance, additional criteria was assessed and the results are displayed in Figure 3.
The results shown in Figure 3 explicit the effect of the increase in the area limit on the optimization model.With the increase in the area limit, there was an exponential increase in the number of constraints, due to the growing number of possible combination of stands to form blocks respecting the established threshold, resulting in a greater complexity for solving the optimization model.The formulation based on Rebain and McDill (2003) with a 10 ha range and the formulation based on Carvajal et al. (2013) showed similar constraint number for all minimum area limits, with the formulation based on Rebain and McDill (2003) presenting constraint number slightly smaller for minimum area limits of 30 and 40ha.For all tested area limits with a 20ha block area range the number of constraints was superior.
We observed that the bounds of scenarios based on Rebain e McDill (2003) formulation with a 10ha block area range were inferior to the bounds obtained in scenarios with a 20ha block range.With smaller range, the block formation may not include all possible stand combinations that respect the minimum area limit.Additionally, stands may not be included in any block.Therefore, the solution space is reduced, generating a limitation in the model, reflected in the scenariosb ounds.The formulation based on Carvajal et al. (2013) presented the best results, with smaller number of constraints and higher bounds.This behavior can also be recognized in the optimization results, as presented in Table 2.
Based on the results displayed in Table 2, we found that the greatest impact on the NPV occurred in scenario 6, with a 5.1% reduction.In general, the minimum area constraints caused reasonable impacts on the NPV.The best outcomes were achieved applying the ring inequalities (Carvejal et al., 2013), with reductions ranging from 0.6 to 1.8%.The constraints based on Rebain and McDill (2003) with a 20 ha block range showed reductions ranging from 1.3 to 2.6% and the scenarios with a 10ha block range displayed the worst outcomes in terms of NPV.
With the increase in the area limit, the number of constraints in the model increased substantially (Table 2), at an exponential rate, and consequently the solution process demanded a higher processing time.Despite the binary variables and the high number of constraints generated for the harvesting blocks formation, it was possible to achieve suitable responses in reasonable processing times for all scenarios.We   obtained gaps between the relaxed optimum and the current solution below 2% after one hour of processing time.

DISCUSSION
All tested formulations could create harvesting blocks respecting the proposed minimum area limits.The ring inequalities formulation, proposed by Carvajal et al. (2013), presented the best results, with inferior number of constraints and suitable responses in terms of the NPV.Nevertheless, in problems with simultaneous consideration of minimum and maximum clear-cut limits, e.g.resulting from certification policies, the formulation based on Rebain and McDill (2003) may be suitable, sufficing to ensure that each stand is part of only one harvesting block.This formulation can be particularly superior when the difference between the minimum and maximum harvesting area limits is small.The range of the block area created through the Rebain and McDill (2003) formulation had a significant impact in the number of constraints and obtained solution.In this sense, with small area range, the number of constraints is reduced.However, with small area range, the probability of stands not being included in any harvesting block increases and the number of possible solutions reduces, which may lead to inferior NPV.
The increase in the minimum area limit led to a reduction in forest NPV.Considering that stands are not uniformly distributed in a forest, in terms of the standing stock, the occurrence of neighbor stands with different areas and standing volume becomes challenging when the wood production flow is required.Thus, it is expected that the impacts on the objective function are enhanced when minimum area constraints are applied in combination with wood flow constraints.For this reason, the inclusion of minimum area constraints might lead to infeasibilities in the optimization model.In this sense, it may be necessary a relaxation of limits of wood flow constraints and management regimes, meaning harvesting stands in ages outside the usual management range, so as to increase the fluctuation limits of wood production.
The impacts of the minimum area constraints were acceptable for all scenarios tested and are compatible with values reported in the literature.For example, Öhman and Eriksson (2010) obtained a 2.6% NPV reduction for achieving a 51.8% reduction in the perimeter of harvesting areas.Our results are also compatible with related adjacency problems, such as maximum harvesting area constraints (Murray and Weintraub, 2002;Zhu and Bettinger, 2007;Gomide et al., 2010;Binoti et al., 2012;Augustynczik et al., 2015).Moreover, with the aggregation of harvesting areas and increase in operational efficiency related to road maintenance and forest harvesting, a cost reduction resulting from blocking harvesting activities is expected.
The inclusion of minimum area constraints increased substantially the complexity of the optimization model, with the addition of a considerable number of constraints Implementing minimum area harvesting... with binary variables.In this sense, the application of heuristic solutions may be suitable for problems involving forests with large number of stands.A series of heuristics were proposed for solving harvest scheduling problems (Pukkala and Kurttila, 2005;Bettinger and Zhu, 2006;Gomide et al., 2009;Bachmatiuk et al., 2015), displaying suitable outcomes.

CONCLUSION
The formulations based on Carvajal et al. (2013) and Rebain and McDill (2003) were capable to successfully incorporate minimum area requirements, aiming to reduce negative edge effects and promote biodiversity maintenance, with acceptable impacts on the economic performance of the forest.Moreover, the formulation based on Carvajal et al. (2013) is superior when the unique manager objective is to incorporate minimum area constraints, due to the reduced number of constraints and a higher efficiency in the solution process.

Table 1 -
Tested scenarios.The table shows all scenarios tested in the optimization, with its respective constraints.

Table 2 -
Augustynczik ALD, et al.Optimization results.The table shows the values obtained for the total volume, Net Present Value, the difference related to the free scenario and the number of constraints for each scenario after the optimization process.Tabela 2 -Resultados da otimização.A tabela mostra os valores obtidos para o volume total, Valor Presente Líquido, a redução em relação ao cenário livre e o número de restrições para cada cenário após a otimização.