An Extended Unit Restriction Model with Environmental Considerations for Forest Harvesting

Abstract

This paper addresses a forest harvesting problem with adjacency constraints, including additional environmental constraints to protect wildlife habitats and minimize infrastructure deployment costs. To this end, we propose an integer programming model to include those considerations during the optimization of the harvest regime of a Mexican forest. The model considered was based on the Unit Restriction Model, a benchmark approach that merges the management units before the optimization process. The resulting model, namely the Green Unit Restriction Model (GURM) and the benchmark model (URM) from the literature were tested with the forest Las Bayas, using information obtained from the SiPlaFor project from Universidad Juárez. The proposed model was solvable in all tested instances. Furthermore, a sensitivity analysis study over a core data set of test instances was carried out on the different parameters of the GURM model to determine optimal configurations for the specific case study. Several environmental measures were assessed in our experimental work. The parameters evaluated were the distance value between pairs of units harvested in the same period, the distance value between those considered natural reserve units, the timber volume to be harvested, the green-up period, and the minimum forest reserve area. An interesting observation from the experiments was that the maximum area inversely affected the URM and GURM models; larger regions resulted in a reduced number of management units in the URM model, thus reducing the computational time to solve the instance of the problem, but in this case, at the expense of a reduced profit. One of the interesting findings was that, in all experiments under all different factors, harvesting every 5 or 6 years yields better profits than harvesting every 10 or 12 years. The current standard in the Mexican system is to harvest every 10 years.

1. Introduction

The collaboration of the operations research and forest science communities has been very successful over the past few decades. The use of operations research analytical tools has positively impacted many areas of forestry management problems, such as strategic, tactical, and operational planning, wildfire management, conservation, and the use of OR to address environmental concerns [1].

One particular area that has received considerable attention is harvest scheduling problems, which in essence involve strategic and tactical decisions on how to harvest a forest, seen as a collection of stands, so that the maximum timber volume or profits are obtained. These decisions are subject to requirements that assure sustainable forest growth. Many of these planning requirements are set forth by government authorities and are designed to preserve wild habitats.

Among these problems, a particular class considers adjacency constraints, which are requirements that prevent the harvesting of adjacent stands or forest management units. The original motivation for these types of constraints was to prevent soil erosion. However, recently there have been other environmental benefits, such as preserving the wildlife habitat or creating a barrier around the recently harvested area until it recovers some of its biomass. These problems consider a regime of clear-cuts distributed along a planning horizon that is divided into periods. An important requirement limits the maximum adjacent area allowed to be cut at any given period.

Several approaches have been studied in the literature to address the adjacency issue. Two well-known approaches are the Unit Restriction Model (URM) and the Area Restriction Model (ARM) [2]. Both approaches differ in how the forest management units, or simply units, were defined out of single forest stands. The URM assumes that the management units are the product of a priori and manual fusion of sub-stands. Contrariwise, the ARM defines the management units from sub-stands, but during the optimization process, while considering maximum area constraints. This paper considers a harvest scheduling problem with URM-type adjacency requirements and additional environmental considerations to help preserve the forest itself and wildlife. Our model includes the preservation of old-growth forest reserves [3], green-up period requirements [4,5], tree biological maturity preservation [6], and distance-based requirements [7,8]. These are further explained in Section 3. These requirements have been proposed independently for other forest management problems; however, to the best of our knowledge, ours is the first model to consider all these simultaneously under the URM approach. In addition, the model is studied under a case study from the state of Durango, Mexico. Despite being the eleventh most abundant country [9] with timber resources globally, Mexico has not systematically used these tools to manage its forests at the level of the nations mentioned above.

To assess our model, we used a case study from Las Bayas. Las Bayas is a piece of land owned by Universidad Juárez, located in Durango, Mexico. The information on geographic and forest measurements is available through the SiPlaFor project [10]. We used real-world data estimates of timber pricing. A sensitivity analysis study over a core data set of test instances was carried out on the different parameters of the GURM model to determine optimal configurations for the specific case study.

The rest of the paper is organized as follows: A literature review is presented in Section 2. Section 3 describes the problem, its requirements, and modeling assumptions. It also includes the integer programming model. Section 4 describes the experimental work in detail, including how the test instances were constructed and a full assessment of the case study. Finally, Section 5 summarizes the main findings and discusses future work.

2. Literature Review

2.1. Review of Basic Adjacency-Based Models (ARM and URM)

The use of operations research tools and methodologies to address many decision-making problems in forestry management has been active over the past 50 years [7]. Many problems involving operational, tactical, and planning decisions have been addressed from the optimization perspective. Some recent surveys can be found in Ezquerro et al. [11], Belavenutti et al. [12], and Rönnqvist et al. [1]. In particular, this section is focused on a unique set of harvesting problems involving adjacency constraints. Adjacency-constrained models based on spatial units were first studied by O’Hara et al. [13], Nelson and Brodie [14], Jones et al. [15], Murray and Church [16], and Barret [17]. Murray [2] coined the term Unit Restriction Model (URM) for these types of models, and also introduced the Area Restriction Model (ARM) in that paper.

The ARM approach is based on limiting, at a particular period, the total area of adjacent harvesting units to a given upper bound and on selecting, at most, a single harvesting set per maximal clique (Figure 1). In contrast, the URM handles the adjacency constraints by establishing that two adjacent units cannot be simultaneously harvested in the same period (Figure 2). There is a third class of adjacency-constrained models that attempts to improve on the limitations of the ARM. These are the so-called hyper unit models [18,19,20], where one hyper unit is generated a priori as a candidate cluster for aggregated forest areas from each forest unit by a predefined rule to meet lower and upper size requirements.

Although the ARM perspective increases the number of potential solutions, the addition of clustering processes within the computational optimization stage causes the problem to become more complex. Many solutions to ARM-based models have been studied in the past, including heuristic approaches [21,22,23,24,25] and exact optimization schemes such as the path formulation [26,27,28], the clique-cluster packing formulation [29,30], and the bucket formulation [31]. Stochastic adjacency models are out of the scope of this paper, but the reader is referred to the work of Wei and Murray [32].

In terms of addressing URMs, which is the focus of this work, in the early years, significant effort was devoted to developing exact and heuristic algorithms. For instance, Weintraub et al. [33] present an excellent survey on exact and heuristic algorithms devoted to solving URMs. These efforts have continued over the past twenty years as we have seen further development of exact [27,28,34], heuristic, and metaheuristic approaches for versions of URMs [25,35].

2.2. Environmental Requirements and Related Works

Most of the URM models from the literature do not consider environmental considerations in their respective formulations beyond adjacency restrictions; among those few, we can cite Bettinger et al. [36], who consider a harvest scheduling problem subject to environmental requirements and adjacency constraints based on URM. Within ARM models, we have seen some studies that consider environmental requirements such as wildlife habitat [37,38,39]. This work proposes a new URM-based model with the following additional considerations. Unlike the ARM, our URM-derived formulation does not use area-related unit aggregation constraints (see Figure 3).

Primary forest reserve: A traditional measure to protect wildlife species is the spatial assignment of natural reserves, known in the literature as an old-growth or primary forest reserve. Primary forest reserve constraints have become more relevant in recent years because these areas may shelter valuable ancient trees and endangered animal species. Carvajal et al. [3] discuss important considerations for integrating old-growth concerns into forest planning with adjacency constraints. They found that forest management practices constantly seek to maintain some link between reserve stands to protect the habitat of vulnerable natural species. They also discuss establishing a requirement that assures a minimum primary forest area, typically between 10 and 30% of the total area. In our model, the requirements for primary forests are considered simultaneously for the maximum distance between reserve stands, where a single group of units must lie within a preset mutual distance.

Green-up period: Green-up constraints induce all units adjacent to a recently harvested unit to remain standing for many consecutive periods [4]. This way, a protective barrier is assigned until the harvested unit restores its forest mass, while we avoid deforesting large areas at any given time of the planning horizon. In addition, this measure reinforces a better harvesting distribution by promoting continuous areas of mature woodland at seeding age adjacent to growing stands. In Murray [5], green-up constraints are defined in the URM formulation by extending the periods in which adjacent clear-cuts to a recently harvested unit are proscribed. We follow this approach in our paper.

Biological maturity: The forest age distribution is one of the most important environmental aspects in harvesting models. A well-distributed age structure at the landscape level ensures a sustained yield of forest products and provides a diverse habitat for wildlife [6]. A way to achieve this goal is to protect the trees of a unit before reaching their age of biological maturity. This parameter can be seen as a threshold from which it is feasible to harvest a basic unit. To the best of our knowledge, the previous optimization work on forest management uses the age threshold mainly to establish the optimal rotation age for carbon sequestration [40], or as a part of heuristic methods for penalizing cuts in immature stands [36,37,39]. We incorporate such a perspective in this work.

Distance Requirements: The consideration of distance in forest management can be raised from the environmental point of view, which aims to maintain a reasonable distance between primary forest reserves and allow wild species, such as elk and deer, to have a safe habitat, or from the economic perspective, given that (i) each harvest implies the deployment of infrastructure, such as new routes and collection areas, and (ii) the reduction of the dispersion of harvested management sets reduces operating costs. Our framework uses the distance notion to allocate reserve patches within a single perimeter to favor wild species regarding survival and reproduction. It is assumed that beyond a certain specific threshold, the functional habitat for such species starts decreasing [8]. We also establish a diameter to create a maximum dispersal area per period between units to be harvested to reduce transportation costs [7]. Since we use a single set of reserve units for all periods and one set of harvest units per period, our approach is to manage small and medium-sized forests, i.e., no larger than 1000 stands.

2.3. Other Related Works

As for larger forest areas, Cova and Church [41] recommend a non-exact method to select a single set of contiguous cells into large territories. The scheme searches through several pre-located seeds until a cell neighborhood is established around one seed. However, by only covering small regions, the approach must be run multiple times to assess the whole territory. Shirabe [42] creates the first exact formulation to create a single patch of contiguous units for a medium-sized territory with 179 land parcels in Maryland, US. Although the problem is solved in a reasonable time, the author acknowledges that combining his model with heuristics would be necessary in more extensive case studies. Later, Ligmann-Zielinska et al. [43] developed a hybrid multiobjective spatial framework named SMOLA, which was applied to a forested area in Washington State, USA. The model, enhanced by a heuristic, uses density-based constraints to assign multiple separated sets of stands with different land use and the distance concept to minimize the dispersion of linked stands. Duque et al. [44] propose three exact models to solve the p-regions problem for territorial aggregation. The models allow the joining of several basic units into p-regions, that is, within multiple sets. As the problem is NP-hard, the exact methods could solve problems up to 50 units and 10 regions in a considerable time (over 3 h), despite commercial solvers and parallel computing advances. Carvajal et al. [3] present an MILP model that imposes connectivity constraints between reserve units dispersed into a single-static set and within forests sized up to 1363 stands.

More recently, Ager et al. have [45] proposed a large-scale bi-objective approach for 79 US national forests, covering more than 59 million ha and 1,443,208 stands. The authors prioritize two conflicting goals: maximizing the timber volume, and promoting fire-reducing activities such as thinning and pile burning. While the approach does not use distance or connectivity constraints, it generated unique assignment solutions for 2350 planning areas. Yet, according to the actual workforce distribution, they found unit assignments infeasible. In Cabral et al. [46], the authors present a framework based on the ARM for extensive forests with several common landowners and management bodies. However, their model does not involve reserve areas nor space-related sets, and it must be pre-adapted by hand depending on the common land at issue. Lastly, Murray et al. [47] extend the work of Shirabe [42] to propose an exact solution for forest management based on MILP to assign p to multiple sets of spatial units. Instead of distance constraints, the authors use connectivity requirements to warranty the units’ linkage within sets. Although the unit allotment use may vary, e.g., creating reserve patches, the method was tested in the allocation of antifire activities for the USDA Forest Service. The authors report minimum run times for a single set (p = 1) and around one hour for multiple sets (p = 10, 50, 100) in forests up to 10,642 spatial units.

2.4. Statement of Contribution

The only work, to the best of our knowledge, that studies a forest harvesting problem using URM adjacency constraints and addressing environmental concerns is due to Bettinger et al. [36]. In their work, the authors study a forest planning process, the development of structural conditions consistent with mid- to late-successional forests, the development of nesting, roofing, and foraging (NRF) habitat within a certain distance of an owl nest location, and the production of revenue. Thus, in their model, among the various constraints, they use constraints to impose minimum and maximum harvest age, to ensure a minimum NRF habitat level for the owl, and to impose green-up delays for adjacent clearcuts. For handling adjacency, they use a URM type of constraint. The difference with our work is that we propose a model that handles primary forest reserve, biological maturity, and distance requirements within a URM adjacency modeling context.

3. Materials and Methods

3.1. Modeling Assumptions

For our model, we assume that the problem is deterministic; that is, all parameters are known with certainty. Clearcut per management unit is the harvesting method to be used. The management policy must respect the maximum clearcut area, adopted by law or voluntary guidelines. The spatial management units come from human-made aggregations of smaller units called sub-stands. All other model assumptions, related to area constraints and economic definitions, are explained in the model below.

3.2. Green Unit Restriction Model (GURM)

The following notation is used for the mathematical formulation:

Indices and Sets

I: Set of harvesting units; $i\in I$.

T: Set of time periods in the planning horizon; $t\in T$.

$N_i$: Set of harvesting units adjacent to harvesting unit $i\in I$.

Parameters

$v_{it}:$ Volume of timber collected in harvesting unit i in period t; $i\in I$, $t\in T$.

${\beta }_{it}:$ Profit obtained by harvesting unit i in period t; $i\in I$, $t\in T$.

$L_t:$ Lower bound for the timber volume to be harvested in period $t\in T$.

$U_t:$ Upper bound for the timber volume to be harvested in period $t\in T$.

r: Green-up period in years.

$m_i:$ Maturity age of trees in harvest unit $i\in I$.

${ϵ}_{it}:$ Age of trees in unit i in period t; $i\in I$ I, $t\in T$.

$a_i:$ Area (in hectares) of unit $i\in I$.

$A_{\min }:$ Minimum area of native forest reservation (in hectares).

${\delta }_{\mathrm{PH}}:$ Planning horizon length in years.

${ϵ}_{\mathrm{nat} }$: Age of trees at which the unit can be considered a primary forest reserve.

${\phi }_{ij}:$ Euclidean distance between units i and j;. $i,j\in I$.

${\mathsf{\Delta }}_{\max }:$ Maximum distance between each pair of units harvested in the same period.

${\mathsf{\Gamma }}_{\max }:$ Maximum distance between each pair of units considered as primary forest reserve.

M: A user-defined large enough constant. This should be large enough to make Constraints (10) and (11) redundant when $x_{jt}=0$ and $z_j=0$, respectively. In our case, we use $M=\max \left\{{\phi }_{ij}:i,j\in I\right\}$.

Decision Variables

$x_{it}$: Binary variable equal to 1 if unit I is harvested in period t, and 0 otherwise; $i\in I$, $t\in T$.

$z_i$: Binary variable equal to 1 if unit i is considered as a primary forest reserve, and 0 otherwise; $i\in I$.

Integer Programming Model

$$\mathrm{Maximize}\hspace{1em}z={{\displaystyle \sum }}_{i\in I}{{\displaystyle \sum }}_{t\in T}{\beta }_{it}{x}_{it}$$

(1)

$$\mathrm{subject} \mathrm{to}\hspace{1em}{{\displaystyle \sum }}_{t^{\prime }=t-r}^{t+r}\left(x_{it\prime }+x_{jt\prime }\right)\le 1\hspace{1em}\hspace{1em}t\in T,i\in I,j\in N_i$$

(2)

$${{\displaystyle \sum }}_{t\in T}{x}_{it}\le 1\hspace{1em}\hspace{1em}\hspace{1em}i\in I$$

(3)

$${{\displaystyle \sum }}_{i\in I}{v}_{it}{x}_{it}\ge L_t\hspace{1em}\hspace{1em}\hspace{1em}t\in T$$

(4)

$${{\displaystyle \sum }}_{i\in I}{v}_{it}{x}_{it}\le U_t\hspace{1em}\hspace{1em}\hspace{1em}t\in T$$

(5)

$$m_i{x}_{it}\le {ϵ}_{it}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}t\in T,i\in I$$

(6)

$$z_i+{{\displaystyle \sum }}_{t\in T}{x}_{it}\le 1\hspace{1em}\hspace{1em}\hspace{1em}i\in I$$

(7)

$${{\displaystyle \sum }}_{i\in I}{a}_i{z}_i\ge A_{\min }$$

(8)

$$z_i=0\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}i\in I : {ϵ}_{i1}+{\delta }_{\mathrm{PH}}<{ϵ}_{\mathrm{nat}}$$

(9)

$${\phi }_{ij}{x}_{it}-M\left(1-x_{jt}\right)\le {\mathsf{\Delta }}_{\max } i,j\in I,t\in T$$

(10)

$${\phi }_{ij}{z}_i-M\left(1-z_j\right)\le {\mathsf{\Gamma }}_{\max }\hspace{1em}\hspace{1em}i,j\in I$$

(11)

$$x_{it}\in \left\{0,1\right\}\hspace{1em}\hspace{1em}\hspace{1em}i\in I,t\in T$$

(12)

$$z_j\in \left\{0,1\right\}\hspace{1em}\hspace{1em}\hspace{1em}j\in I$$

(13)

The objective function (1) seeks to maximize the harvesting profit throughout the planning horizon. It is assumed that the wood volume’s growth pace depends on the average annual increase in the number of standing trees. We also assume that the harvesting profit in a specific period is a function of the wood potential, prices, and production costs per cubic meter of wood of the tree species in the unit, and that selling prices per cubic meter of wood are static over the planning horizon. Constraints (2) prevent any neighboring unit pair from being harvested simultaneously before the green-up period r. It is assumed that the green-up or regeneration barrier length should allow the natural regeneration of most of the trees of a unit after its harvesting. Strong and weak adjacency types (see Section 1) are used to avoid harvests between neighboring units, such that simultaneous harvests have no point of contact to reinforce the continuity of the forest. For simplicity, we assume that all management units have a surface between 51% and 100% of the maximum clearcut area. Namely, the simultaneous harvest of any pair of adjacent units violates the maximum area. Constraints (3) assure that a unit is harvested at most once during the planning horizon. Constraints (4) and (5) guarantee the lower and upper limits of the timber volume harvested in each period. Constraints (6) guarantee that only the units that have reached the maturity age in the current period are considered for harvesting. It is assumed that age of (manually constructed) units corresponds to the age of the youngest trees into it, so that no clear-cut interferes in the maturing vegetation process. The age of trees in units is assumed homogeneous. Constraints (7)–(9) help us model the primary forest reserve constraints. Constraints (7) assure that no unit can be harvested (at any time) and simultaneously be considered a primary forest. Constraint (8) guarantees a minimum area exclusively dedicated to the primary forest reserve. As young or immature trees do not satisfy primary forest requirements, minimum age constraints must be imposed to prevent impractical assignments. To this end, Constraints (9) ensure that units with relatively immature trees are not assigned to the primary forest reserve. Note that this is handled in fact not as a constraint, but simply fixing those variables to zero and thus eliminating them from the model. Constraints (10) establish the maximum distance allowed between any pair of units harvested in the same period to decrease the operational costs. Constraints (11) set the distance limit between any pair of primary forest reserve management units to guarantee the habitat of the wildlife. Finally, Constraints (12) and (13) express the nature of the binary variables. Note that Constraints (6) can be eliminated from the model by defining only those units which are mature for harvesting, that is, setting $x_{it}=0$ for those i and t, such that $m_i>{ϵ}_{it}$. This allows us to reduce the number of constraints and variables in the model.

The URM, as proposed by Murray [2], considers maximizing (1) subject to (3)–(5) and constraints

$$x_{it}+x_{jt}\le 1, t\in T,i\in I,j\in N_i$$

instead of Constraints (2). Therefore, our model can also be seen as an extension of the URM.

The URM is NP-hard [29,48]. Despite its inherent computational complexity, the URM has been solved relatively well in instances with up to 80,000 stands over 10 time periods by branch and bound [28]. Typically, integer programming models can be solved by the branch-and-bound algorithm [49]. In our case, we use a commercial branch-and-bound method (from the CPLEX library) for solving this integer programming model.

3.3. Preprocessing

The implementation of the solution framework can be divided into two stages: (a) the construction of instances from the geographical, economic, and measurement information of the forest, and (b) the solution of the corresponding integer programming model.

The information about the individual stands, such as the projected volume of wood that a stand can yield at each period of the planning horizon, the profit that can be obtained from the timber mentioned above, and the average age of each stand, was calculated first using the available databases. In this case study in Mexico, the geographical data were obtained from the SiPlaFor project. SiPlaFor is an open-access system designed to support the decision-making process in preparing and executing sustainable forest management planning on temperate forests of Mexico [10]. SiPlaFor also contains information on the species, age, and samples of trees from each stand, among other factors, to calculate the potential volume of timber harvested at each period during the planning horizon, considering an even-aged whole-stand model. For this project, we obtained the profit from the previously calculated volume of wood per stand during each harvesting period and the proportion of species in the stand. We used a selling price of $1163 Mexican pesos/m³ [50].

3.4. Case Study

In this paper, we present a case study of Las Bayas Forest, a wooded terrain located in the state of Durango in northwest Mexico, as depicted in Figure 4. As seen in Figure 4, Las Bayas is located in an area between −104.9 and −104.8 of longitude, and 23.4 and 23.5 of latitude. Additionally, detailed UTM coordinates are also available in Corral Rivas et al. [51]. Las Bayas is a forest comprised mainly of pines and oaks. Its selection was motivated by the following reasons:

• Las Bayas belongs to the Universidad Juárez del Estado de Durango (UJED), under the management of the Faculty of Forestry. Therefore, its database contains high-quality information from the stands, freely available to the public [10].
• Las Bayas information is available in the SiPlaFor project, also managed by the UJED. The project contains information about other forests in Mexico, which makes it easier for our approach to be applied to other forests if the information in the databases meets the minimum requirements for the model.

For this case study, we selected 118 stands, as shown in the shaded areas in Figure 5. Those are some of the stands authorized to carry out timber extraction activities. Surrounding stands, instead, are considered primary forest reserves for agriculture or tourism.

3.5. Experimental Settings and Test Instance Generation

All models were implemented in the C++14 programming language, compiled in g++ version 7.4.0, and solved with the ILOG CPLEX 12.8 solver through its API. For test instance generation and data recollection, the programming language used was Python 3.7.4. The election of this language was mainly because of the availability of powerful GIS (Geographic Information System) libraries, such as GeoPandas and pyproj. The experimentation was conducted in a server with an Intel Xeon E3-1240 v3 3.5 at GHz with 8 cores, with 16 GB of DDR3 RAM at 1600 MHz, in the Ubuntu 14.04.6 LTS operating system.

As mentioned in Section 1 and Section 3, the GURM requires harvesting units to be hand-built by management staff, merging smaller adjacent units called sub-stands. With this in mind, we formed forest management units by aggregating sub-stands whose total area did not exceed a specified maximum size. Some average statistics per hectare about Las Bayas are: 592 trees, 19.24 m² of basal area, 151.78 m³ of volume, and 3.83 m³ of current annual increment. The main tree species are Pinus durangensis, P. engelmannii, P. cooperi, P. teocote, Quercus sideroxyla, and Q. crassifolia, with a suggested rotation of 60 years for all species. We studied three maximum area values for this: 30, 40, and 50 ha. The length of the planning horizon was established at 60 years. This length of the horizon period corresponds to the usual rotation in the study area. The current practice establishes a harvesting period every ten years [52]; however, for this study, we have an interest in investigating the effect that different values of rotation periods could have. Therefore, for this planning horizon of 60 years, different harvesting periods were evaluated: (i) 5 periods (that is, a rotation of 12 years), (ii) 6 periods (a rotation of 10 years), (iii) 10 periods (a rotation of 6 years), and (iv) 12 periods (a rotation of 5 years). Since (number of rotation periods) x (rotation length) = 60, and to avoid confusion, for the purpose of displaying results and discussion, in the remainder of the paper we will use the number of rotation periods (identified as P).

With these values, we formed the core instances, from which we built the instances for each experiment by varying the values described next and their default values.

The core instances are described in

Table 1. Description of the core instances.

Instance	Max Area (ha)	P	MU	Vars	Cons
bayas-urm-30-1	30	5	56	336	38,579
bayas-urm-30-2	30	6	56	392	44,909
bayas-urm-30-3	30	10	56	616	70,637
bayas-urm-30-4	30	12	56	728	83,433
bayas-urm-40-1	40	5	53	318	34,586
bayas-urm-40-2	40	6	53	371	40,259
bayas-urm-40-3	40	10	53	583	63,323
bayas-urm-40-4	40	12	53	689	74,793
bayas-urm-50-1	50	5	54	324	35,885
bayas-urm-50-2	50	6	54	378	41,773
bayas-urm-50-3	50	10	54	594	65,703
bayas-urm-50-4	50	12	54	702	77,605

Note: Max Area = Maximum area of management units; P = Number of rotation periods in the 60 year planning horizon; MU = Number of management units; Vars = Number of decision variables; Cons = Number of constraints in the model.

, where the first column is the name of the core instance, and the second, third, and fourth columns show the maximum area of the management units of the instance, the number of periods, and the number of management units that it contains. Finally, the fifth and sixth columns show the number of decision variables and constraints of the resulting model for that instance.

The upper and lower bounds of the wood volume harvested at each period were calculated by considering the average potential yield of the forest at each period, divided by the number of periods, giving the average potential per period (ap). We also considered a volume variation percentage (vvp) with respect to vp, that is, a value to be added or subtracted to the ap value to obtain the upper and lower bounds, e.g., upper = ap(1 + vvp), lower = ap(1 − vvp). Note that vpp is a value in the [0,1] range, where a 20% deviation, for instance, corresponds to vpp = 0.20. In this work, the default vvp value was established at 0.10 (or 10%).

The set of distances (${\mathsf{\Delta }}_{\max }$, ${\mathsf{\Gamma }}_{\max }$) that represents the maximum distances between any two pair of management units and native forest, respectively, at any period, was established with the default value (8000 m, 5500 m). These parameters appear in Constraints (10) and (11). The minimum required fraction of timberland total area, to be assigned as the native forest parameter $A_{\min }$ in Constraint (8), was set at $A_{\min }$ = 0.10 or 10%. The green-up years, that is, the number of years required before harvesting a management unit that is adjacent to a recently treated management unit, was set at r = 20 years. This parameter r plays a role in Constraints (2).

The time threshold was set at 3600 s. If this time is higher at the start of a new iteration, the optimization process is stopped. A second stopping criterion considered was a relative optimality gap of 0.01%.

4. Results and Discussion

4.1. Experiments and Results

This section describes the experiments conducted to determine the impact that different parameters of the instances, complementary to the core parameters, had on the objective value of the Green URM (GURM). In particular, we investigate the effect of the following factors or parameters on net revenue given by objective function: (a) the maximum area of management units, (b) the maximum distance between each pair of units harvested in the same period (parameter ${\mathsf{\Delta }}_{\max }$ in Constraints (10)), (c) the lower and upper variation of limits of timber volume per period (Constraints (4) and (5)), (d) the green-up period length (parameter r in Constraints (2)), € the minimum forest reserve area (parameter $A_{\min }$ in Constraint (8)), and (f) the rotation length.

To achieve this, different experiments were carried out, by varying a different parameter for each experiment while keeping the others fixed with the default values mentioned in the previous section. This is explained in detailed below.

4.2. Core Parameters with Default Values

For this experiment, we evaluated the core instances shown in Table 1, with the default values of the complementary parameters described earlier in Section 4.3. Such a configuration resulted in a total of 12 test instances.

The core parameters include three different maximum areas: A1 = 30 ha, A2 = 40 ha, and A3 = 50 ha; and four different numbers of rotation periods: P1 = 5, P2 = 6, P3 = 10, P4 = 12, which resulted in 12 test instances.

This experiment aims at having a baseline on which to compare the results obtained by varying the parameters individually.

Table 2. Results for the core instances experiment.

Max Area (ha)	P	OV (mill. MXN)	T (sec.)
30	5	$90.82	22.58
30	6	$90.98	31.07
30	10	$92.47	1455.84
30	12	$92.49	2857.95
40	5	$90.80	22.19
40	6	$90.82	39.01
40	10	$92.34	595.38
40	12	$92.48	129.92
50	5	$90.86	31.71
50	6	$90.82	57.96
50	10	$92.42	607.46
50	12	$92.37	2446.49

Note: Max Area = Maximum area of management units; P = Number of rotation periods in the 60 year planning horizon; OV = Objective function value; T = Time employed to solve the instance (CPU seconds).

shows the results obtained from the 12 instances evaluated in this experiment. The default configuration was feasible. The first and second columns show the maximum area and the number of periods, respectively. The third column shows the objective value in millions of Mexican pesos (OV (mill. MXN)), and the fourth column shows the time in seconds of CPU (T). For each instance tested, a global optimal solution was found (within a 0.01% precision), and the runtimes were reasonably short.

4.3. Assessing the Effect of Distance Parameters

As mentioned before, distance is an essential factor to consider in forest management, both from an economic and environmental point of view. From the economic perspective, our model addresses, in Constraints (10), the cost of deploying the necessary infrastructure to carry out the harvest in a certain period, as the more dispersed the management areas to be harvested, the higher the cost of building roads, making trips between them or installing collection points.

From the environmental point of view, our model considers, in Constraints (11), that native forest areas must maintain a certain maximum distance to fulfill their purpose of safeguarding the habitat of forest wildlife.

For this experiment, the pairs of distances of the 12 core instances were varied with four different pairs of values {${\mathsf{\Delta }}_{\max }$, ${\mathsf{\Gamma }}_{\max }$}, namely: {8500 m, 5750 m}, {9000 m, 6000 m}, {9500 m, 6250 m}, and {10,000 m, 6500 m}, resulting in a total of 48 test instances.

Figure 6 shows the results obtained. Each row corresponds to a group of instances with the same maximum area. The different pairs of ${\mathsf{\Delta }}_{\max }$ and ${\mathsf{\Gamma }}_{\max }$ are represented as follows: D1 (the default, shown in Table 2) is {8000, 5500}, D2 is {8500, 5750}, D3 is {9000, 6000}, D4 is {9500, 6250}, and D5 is {10,000, 6500}. All solutions were feasible and optimal (within 0.01% precision).

The first four subplots correspond to the objective function value, and the last four subplots to the running time within each row. Each subplot displays the results for a given value of P (5, 6, 10, 12). Within each subplot, the set of distances defined in the experiment was evaluated.

As shown in Figure 6, the variation among the different distance values is negligible for a given fixed area and harvesting period. However, for a fixed distance value, it can be observed that a p value of 12 or 10 (that is, harvesting every 5 or 6 years) yields better results than a p value of 6 or 5 (that is, harvesting every 10 and 12 years).

4.4. Assessing the Effect of Volume Percentage Threshold

Maintaining a constant harvest volume during each period is an important consideration so that the economic benefit also remains constant. However, it is not always possible for the timber volume to be the same every time.

Therefore, in Constraints (4) and (5), our model includes lower and upper limits of timber volume per period, which, as they get close to each other, the harvest between each period will be more constant. However, this has the disadvantage that, if they were too strict, it could result in unfeasible solutions for that particular instance. On the other hand, the more separated they are, the easier it will be to find feasible solutions, but the volume will likely vary too much between periods.

In this experiment, the volume variation percentage of the 12 core instances was varied with four different values: vvp2 = 12%, vvp3 = 16%, vvp4 = 20%, and vvp5 = 24%, to determine how tight the lower and upper volume limits per period can be for Las Bayas. This configuration resulted in a total of 48 test instances.

Figure 7 shows the results obtained. Each row corresponds to a group of instances with the same maximum area. The first four subplots correspond to the objective function value, and the last four subplots to the running time within each row. All solutions were feasible and optimal (within 0.01% precision). Each subplot displays the results for a given value of P (5, 6, 10, 12). Within each subplot, the set of different values of vp previously defined in the experiment was evaluated.

As we can see from Figure 7, for a given fixed area and harvesting period, the best results are observed for the largest value of vp = 24. The decrease between this value and the smallest is less than 2%. About the variation among periods for a fixed value of volume, it can be observed that a p value of 12 or 10 (that is, harvesting every 5 or 6 years) yields better results than a p value of 6 or 5 (that is, harvesting every 10 and 12 years). This result is consistent with the previous experiment.

4.5. Assessing the Effect of Green-Up Period

To avoid the formation of large deforested areas, even when adjacency restrictions are respected during the same period, the recovery period or green-up is included alongside the adjacency constraints of our model in Constraints (2).

The green-up period consists of avoiding, for several consecutive periods, the harvesting of management units adjacent to another one that has been recently treated, in order to leave a protective barrier of standing trees that allows the unit to recover some of its forest mass.

The longer the green-up period, the larger the forest mass is recovered. However, it also limits many of the options for harvesting management units in future periods, resulting in unfeasible solutions.

For this experiment, the green-up lengths of the 12 core instances were varied with four different values for parameter r in Constraints (2), namely: r = 10, 15, 25, 30. This resulted in a total of 48 test instances.

The results of this experiment are shown in

Table 3. Results for the green-up experiment.

		r = 10		r = 15		r = 25		r = 30
MA	P	OV	T	OV	T	OV	T	OV	T
30	5	$90.99	17.07	$90.81	22.12	$89.66	6.81	$89.66	5.87
30	6	$91.62	21.12	$91.62	24.44	$90.97	22.39	$88.26	8.45
30	10	$93.06	1844.31	$92.94	2423.81	$91.24	46.89	-	-
30	12	$93.34	759.48	$93.08	1755.67	$90.29	830.10	-	-
40	5	$90.98	19.74	$90.79	21.98	$89.35	10.48	$89.35	13.59
40	6	$91.60	24.36	$91.60	21.44	$90.81	35.99	$87.44	9.97
40	10	$93.02	889.15	$92.85	958.92	$91.22	52.44	-	-
40	12	$93.27	2889.52	$92.92	3601.41	$90.24	3228.85	-	-
50	5	$91.02	28.67	$90.85	12.86	$89.58	3.97	$89.59	3.97
50	6	$91.65	25.72	$91.65	20.68	$90.81	12.99	$88.11	5.64
50	10	$93.05	2152.11	$92.89	663.24	$91.36	143.24	-	-
50	12	$93.32	959.59	$93.00	1043.45	$91.27	45.94	-	-

Note: MA = Maximum area of management units (ha); P = Number of rotation periods in the 60 year planning horizon; OV = Objective function value (millions of Mexican pesos); T = Time employed to solve the instance (CPU seconds).

. In this case, 42 configurations were feasible. The first and second columns show the maximum area and the number of periods, respectively. The following columns indicate the results for the different values of r tested. Each two-column set shows the objective value in millions of Mexican pesos (OV) and the time in seconds (T). For each instance tested but one, an optimal solution was found (within 0.01% precision). The instance that was stopped by the time limit had a 0.23% optimality gap.

Figure 8 visually shows the results obtained. Each row corresponds to a group of instances with the same maximum area. The first column corresponds to the objective value obtained, and the second is the time it took for the method to find the optimal solution. The groups are divided into four subgroups that correspond to the periods considered for the core instances. For each period, the different green-up lengths of r were evaluated.

As we can see from Table 3 and Figure 8, for a given fixed area and harvesting period, the best results are observed for the smallest value of r = 10. As r gets larger, the objective function decreases by about less than 2%. About the variation among periods for a fixed value of the green-up period, it can be observed that a p value of 12 or 10 (that is, harvesting every 5 or 6 years) yields better results than a p value of 6 or 5 (that is, harvesting every 10 and 12 years). This result again is consistent with the previous experiments.

4.6. Assessing the Effect of Minimum Forest Reserve Area

As already mentioned, maintaining a minimum forest reserve area percentage (Constraints (7)–(9) of our model) is important as a measure to protect plant and animal species. In recent years, policies in which a minimum percentage of the forest must be preserved for this purpose have been established.

It is intuitive to deduce that the larger the minimum area of native forest, the lower the economic benefit, such that the purpose of this experiment was to evaluate how relevant the impact is by testing different percentages of native forest.

To this end, the minimum forest reserve area of the 12 core instances was varied with four different values of $A_{\min }$, namely: 2%, 6%, 14%, and 16%, resulting in a total of 48 test instances.

Figure 9 shows the results obtained. Each row corresponds to a group of instances with the same maximum area. The first column corresponds to the objective value obtained, and the second is the time it took for the method to find the optimal solution. The groups are divided into four subgroups that correspond to the periods considered for the core instances. For each period, the set of four $A_{\min }$ percentages previously defined in this experiment was evaluated. For each instance tested, an optimal solution was found (within 0.01% precision).

As we can see from Figure 9, for a given fixed area and harvesting period, the best results are observed for the smallest value of $A_{\min }$ = 2%. As $A_{\min }$ gets larger, the objective function decreases by about 5–9%. About the variation among periods for a fixed value of minimum forest reserve area, it can be observed that a p value of 12 or 10 (that is, harvesting every 5 or 6 years) yields better results than a p value of 6 or 5 (that is, harvesting every 10 and 12 years). This result again is consistent with the previous experiments.

4.7. Example of a Harvesting Planning

Figure 10 shows an optimal plan obtained, for instance in bayas-urm-30-1 (see Table 1), with a max. area = 30, periods = 6, management units = 56, and default complementary parameters of the instance required by the GURM (defined in Section 4.4), which correspond to the values of volume variation percentage, the maximum distance between native forest units, the maximum distance between management units, the minimum required area percentage of native forest, and years required before harvesting a management unit that is adjacent to a recently treated management unit.

We also solved the same instance but using the traditional URM model instead; that is, ignoring the additional constraints stated in Constraints (2) and Constraints (6)–(11), with default values. This solution is shown in Figure 10a. Each color represents a period of harvesting, and in the case of the GURM, the white color represents a management unit turned into a primary forest reserve, and the gray color represents a sub-stand that was not considered in the harvesting.

By contrasting the harvesting plans in both figures, the GURM enables a minimum area of native forest within a limited perimeter. The GURM also hinders the harvesting of too-distant units in the same period, and thanks to the green-up barrier, it prevents the harvesting of neighboring units in consecutive periods. For instance, in Figure 10a, we can find harvesting prescriptions for adjacent units in consecutive periods (1 and 2, 4 and 5, 5 and 6), while this is not observed in its counterpart in Figure 10b, thanks to the green-up constraints incorporated in our model.

4.8. Discussion

As for the GURM, the change in the values of the different constraints caused decreases in the OV of −3.23% at most for the default configuration. However, in most cases, the OV was hardly affected, and there were instances in which this value increased when some constraint thresholds were relaxed. Regarding the management and native forest distance constraints, the OV decreased by 2.17% for the benchmark when the corresponding thresholds varied. Regarding the volume threshold constraints by period, the OV decreased by 1.79% and increased by 0.28% compared to the benchmark instance, as the values were restricted and relaxed, respectively. Green-up constraints were among those that most negatively affected the objective value. In such a case, the economic benefit decreased from 1.93 to 3.11% when different configurations were compared, versus the default setup. It was also observed that the native forest constraints significantly impacted the objective value. When suppressed, the OV increased by 2.02%. When strengthened, the OV decreased up to −3.23%.

An extra finding was that the m”ximum area inversely affected the URM and GURM returns, as larger areas resulted in fewer management units. Nevertheless, as the maximum area grew, the runtime to solve the problem instances decreased. Similarly, for extensive forests, say 1000 plus stands, GRUM Constraints (11), which is aimed at establishing a max. distance between native forest cells, is not designed for multiple sets but for a single static set, which could be limiting at some point. For example, if we increase the max. distance value, a single old-growth set might still be functional for medium-sized forests. However, such an adjustment could be insufficient in forests with thousands of cells. Moreover, a single but broad set of reserve patches could include dispersed native forest stands or far from regions with the same foliage variety.

Another interesting observation is as follows. Note that as the p value increases, the model has more integer variables and constraints (Table 1), and it is in general harder to solve, meaning that the integer programming method takes longer to find an optimal solution. This behavior can be observed, for instance, in Figure 6. However, if looking at Figure 7, this behavior is not quite observed in some cases. Take the results corresponding to vvp = 12%, for instance. We can observe that solving an instance with a p value of 10 takes longer than solving an instance with a p value of 12. This behavior, although not too common, is sometimes observed when solving integer programming problems. The explanation for this has to do with how branch-and-bound solvers work. Sometimes, instead of making the problem more difficult, having a larger number of constraints benefits in making a better adjustment to the convex hull of integer solutions, providing tighter bounds that make the branch-and-bound enumeration tree prune and eliminate more nodes in its enumeration tree. Thus, in this particular case, for the data corresponding to those instances, this behavior is observed.

Finally, in all experiments under all different factors, it was observed that a p value of 12 or 10 (that is, harvesting every 5 or 6 years) yielded better results than a p value of 6 or 5 (that is, harvesting every 10 and 12 years). The current standard in the Mexican system is to harvest every 10 years. The results indicate that the main factor observed to be significantly affecting the revenue was the length of the harvest rotation. The other tested factors did not cause a significant variation in the revenue. However, it would be interesting to assess the model and methodology in other forests under different conditions. In that regard, the work developed here can be valuable for similar studies.

5. Conclusions

In this paper, we proposed extending a model from the literature with environmental constraints to solve the forest harvesting problem with adjacency, to optimize the harvesting of a Mexican forest during a planning horizon. The model (GURM) was based on the Unit Restriction Model (URM), which considers the management of units beforehand. The resulting model GURM does not intend to compete with the ones from the literature regarding maximizing harvesting profit, but to maximize the profit while considering environmental constraints that protect the wildlife and the forest. In this regard, our model contributes to the literature on URM-based adjacency models considering some environmental constraints.

The current standard in the Mexican system is to harvest every 10 years. The results from this study indicate that the main factor observed to be having a significant effect on the revenue was the length of the harvest rotation, and that better results could be obtained if we change the rotation to harvest every 5 or 6 years instead of the current standard. This is a very interesting finding. The other tested factors did not cause a significant variation in the revenue. However, assessing the model and methodology in other forests under different conditions would be interesting. In that regard, the work developed here can be valuable for similar studies.

Future work in this area will include implementing and reformulating other models with spatial constraints, specifically the Area Restriction Model (ARM), which sets up a more extensive set of management units that could potentially increase the profit obtained during the harvesting planning. Moreover, although our distance constraints can cover many instances, adding multiple sets will also be considered for extensive forest cases.