Introduction 

In a broad sense, a system is considered to be stable if it always returns to a reference position (equilibrium) after small disturbances. In a forest ecosystem the main factors influencing this equilibrium are recruitment, transitions between diameter classes, and mortalities (natural or induced by harvesting - [21]). These processes of evolution over time can be easily modeled using a matrix population model (e.g., [1], [2], [8], [30]), for which the population growth rate is the dominant eigenvalue λ of the transition matrix and the equilibrium position is determined by the stable diameter distribution W₀. Stand management practices aimed at conserving this equilibrium are important because they represent the stand’s ability to persist in the face of external disturbances. Managing a stand by applying the “sustainable/stable” harvesting strategy ([9], [10], [12], [13]) is a good way to reach and maintain that equilibrium with sustainability. In fact, the “sustainable/stable” harvesting strategy is aimed at reaching the stable diameter distribution of the stand for each harvest, yielding for all the diameter classes and time steps the following “sustainable/stable” harvest rate (eqn. 1):

where λ₀ is the dominant eigenvalue of the transition matrix without harvests, and s also includes the natural mortalities. This is a sustainable harvesting strategy too (in the sense of population persistence over time) because s corresponds to the case such that the right eigenvector associated with the dominant eigenvalue λ = 1 of the transition matrix with harvests, is W₀.

By contrast, many harvesting strategies do not satisfy the requirements of sustainability or stability. In fact, using an objective function that integrates financial data of harvesting operations with a matrix model to describe the population dynamics, a recent discrete optimal control model proposed by López et al. ([13]) showed that the economically optimal harvesting strategies underlying the solutions of that model changed the stand towards distributions more characteristic of even-aged stands and substantially deviated from the stable diameter distribution. Moreover, the corresponding economically optimal harvesting strategies were always unsustainable, in the sense that the population growth rate λ_T for the whole harvest cycle was always λ_T < 1, being λ_T the dominant eigenvalue of the transition matrix from t = 0 to t = T, and T the planning horizon.

This study investigates the existence of a new sustainable, stable, and economically optimal harvesting strategy for uneven-aged Pinus nigra stands in the Spanish Iberian System. We rely on an objective function that integrates financial data of harvesting operations with a matrix model that describes the population dynamics, which creates a discrete optimal control model. In addition, it is investigated whether the sustainability and stability conditions can be satisfied by assuming the stable diameter distribution of the stand as the initial and final conditions, that is X(0) = X(T) = W₀. Furthermore, we examine the effects of diameter growth, stand density, and recruitment on those strategies, considering a wide variety of typical management scenarios in the study area, obtained by combining three rates of diameter growth, three levels of stand basal area (G = 22, 24, and 26 m² ha^-1), and three levels of global recruitment (scarce: R = 200 stems ha^-1; normal: R = 520 stems ha^-1; and abundant: R = 840 stems ha^-1). Finally, an economic comparison has been performed between this economically optimal harvesting strategy between two stable positions and the “sustainable/stable” strategy for each scenario.

In general, optimal control theory enables the solution of a wide variety of dynamic problems, for which the evolution of the dynamic system (discrete or continuous) can be partially controlled by the agent’s decision. In every moment t, the system can be described by a set of state variables x_k(t) (e.g., stems ha^-1 in diameter class k at year t), so the manager chooses a specific set of control variables h_k(t) (e.g., harvest rate in diameter class k at year t). Different values for the control variables imply different paths in the phase space for the dynamic system and the manager must determine the control values that maximize the selected objective, according to the constraints of the problem.

Many authors have considered optimal control theory to model forest stands management both in continuous case ([5], [23]) and in discrete case ([7], [25], [11], [13]).

In the following sections, we describe the model and the estimates of the corresponding parameters. The general calculations, including regressions, employed the software Maple ([15]) version 18, and the solutions to the optimization problems were derived using the Solver^® tool running under MS Excel 2013^®.

  Materials and methods 

Discrete optimal control model

This model has three main features: (1) an underlying discrete time dynamic system, given by the stand population dynamics; (2) an objective function, which is in this case the net present value (NPV) of timber income over a fixed time horizon T; and (3) a set of constraints.

The combination of the objective function with the population dynamics model and the constraints produces a discrete optimal control problem, from which we have to determine control variables h_k(t) that globally maximize the function NPV subject to the population dynamics and constraints.

Stand population dynamics

The stands of the study area are located in the Spanish Iberian System, a mountain range extending about 400 km along the north-eastern edge of the Central Plateau, including 60% of the area occupied by Pinus nigra in Spain ([4]). The management of Pinus nigra stands in this area is currently based on the growth and yield tables for Pinus nigra Arn. in the Spanish Iberian System ([3]). The data in these tables are classified according to growth rates into five site qualities, of which only the first three have been selected as our data source. These three qualities were defined by the dominant heights reached at the reference age of 60 years (20, 17, and 14 m, respectively). The study area is very large, covering the whole range of P. nigra in the Spanish Iberian System. Therefore, the spectrum of values taken by environmental variables such as slope, soil type, and even climatic characteristics is also very wide. This variability is captured in the growth and yield tables used by considering different site qualities that reflect variations in topography, height, slope, aspect, weather, and vegetation.

The population dynamics model that integrates recruitment, harvesting, and stem migration throughout the diameter classes over time is based on the model proposed by López et al. ( [9], [10]) for Fagus sylvatica L., and by López et al. ( [12], [13]) for Pinus nigra Arn. stands.

The starting assumptions were as follows: (a) the forest is in a steady state; (b) the average diameter growth curves for site indices 20 to 14 are defined by the Bertalanffy-Richards model ([31], [19]); (c) within each diameter class, the probability distribution for the diameter of the trees is uniform (rectangular); and (d) harvesting operations occur at the beginning of every projection interval. Group selection cutting is the silvicultural system used because P. nigra is not a shade-tolerant species, and the gaps opened in the canopy through single-tree selection cutting would not be large enough to ensure regeneration.

Harvesting operations in the study area generally take place every 10 years, thus we adopted this period as the range of the projection intervals in the model. Considering this time period and the diameter growth functions corresponding to site indices 20 to 14, trees were grouped into n diameter classes of equal width w = 6 cm: 0-6, 6-12, 12-18, …, 30-36, …, with the last class being more than 48 cm (48,→) for site index 20, more than 42 cm (42,→) for site index 17 and more than 36 cm (36,→) for site index 14 (note that n = 9 for site index 20, n = 8 for site index 17, and n = 7 for site index 14). Therefore, an individual tree in class k can remain in class k or progress to class k + 1 during the projection interval (t, t + 10). The number of trees in each class changed for each projection interval because some were harvested, some remained in the same diameter class, and others grew past the boundary to the next diameter class. In such conditions, p_k refers to the probability that an individual tree in class k at time t will appear in class k + 1 at time t + 10. The term r_k, or recruitment coefficient, is the number of offspring (stems ha^-1) living at time t + 10, produced in the interval t, t + 10 by an average tree in class k at time t. Similar to many standard size classified matrix models, we did not consider the smallest diameter class to be fertile ([2]). The variable h_k(t) defines the proportion of harvested trees in class k at year t, natural mortality included. Finally, x_k(t) and x_k(t + 10) describe the stem densities in class k at the initial and final projection times. By analysing the dynamics of the projections, we find that the model is described accurately by the following matrix model (eqn. 2):

where (eqn. 3):

and I is the identity matrix; H(t) = diag[h₁(t), h₂(t), …, h_n(t)] is a diagonal matrix with the harvest rates h_k(t); and X(t) and X(t + 10) are column vectors indicating the stem densities at the initial and final times of projection, respectively.

Objective function

Considering that the main goal for the stands in the study area is the production of timber, the objective function of the optimization problem is the NPV of all management operations over a planning horizon of T = 70 years, discounted to the beginning of the period and from the landowner’s perspective, that is (eqn. 4):

where v_k(t) is the stumpage value that corresponds to class k at year t (€/stem) and i is the discount rate. The first summation represents the income derived from timber sale and the second is the final stocking value. It is well known that the NPV is a justifiable management objective for a single economic goal ([24]).

Harvesting costs are not considered because harvesting operations are not paid by the owner. The costs for the owner of the land are the opportunity cost of the invested capital and the opportunity cost of the land. None of them are explicit costs, but they are incorporated into the NPV obtained. This hypothesis is not restrictive since active forest management in the Mediterranean area is not frequent and the owner does not invest any effort in silvicultural actions.

The cost of harvesting operations would be paid by the buyer of the stumpage and that is the reason we did not consider it in the present study.

Set of constraints

As mentioned above, the initial and final conditions must represent the stable diameter distributions of the stand W₀. Therefore, the first constraints are given by (eqn. 5):

We also set the following obvious constraints (eqn. 6):

for k = 1, 2, … , n and t = 0, 10, 20, 30, … , T - 10, where m_k is the natural mortality in k-th class.

On the other hand, by applying the “sustainable/stable” harvesting strategy to each scenario and setting X(0) = W₀ as the initial condition, the basal area of the stand continuously oscillates between values G_min (after harvesting) and G_max (before harvesting), and the stable diameter distribution is reached at each time step. Therefore, to compare the economically optimal and the “sustainable/stable” harvesting strategies, it is essential to assume an additional constraint to maintain the stand basal area throughout the harvest cycle within the range G_min - G_max during the optimization process, where G_max was 22, 24, and 26 for the three levels G = 22, 24, and 26 m² ha^-1, respectively, and G_min is the minimum basal area reached for each scenario under the “sustainable/stable” harvesting strategy.

Optimization method

For each discount rate, the tool Solver^® of MS Excel 2013^® was used to obtain a globally optimal solution of these nonlinear models by applying multistart methods. These methods provide a “convergence in probability” to a globally optimal solution, which means that as the number of runs of the nonlinear Solver routine increases, the probability that the globally optimal solution has been found also increases towards 100%.

The control variables h_k(t) that maximize the objective function NPV (always under the population dynamics and the constraints mentioned above) are the solutions to the corresponding problem (one optimal solution for each scenario). Finally, the optimal harvesting schedules are defined by the combination of the control variables h_k(t) corresponding to the scenario analyzed.

Input estimation

Transition probabilities

We calculated the transition probabilities by applying the method introduced in López et al. ([12]) considering that, given the basal area constraints and the low basal area scenarios of this study, stand density has a small influence on diameter growth. Thus, the diameter growths were defined by the Bertalanffy-Richards function ([31], [19] - eqn. 7):

adjusted using regression analysis to data points [t, D(t)] obtained from tables by Gómez Loranca ([3]). The main characteristics of these growth models are shown in Fig. 1.

Using these diameter-growth models, the transition probabilities from class i to class (i + 1) or from interval [6(i- 1), 6i] to interval [6i, 6(i + 1)] are given by (eqn. 8):

for i = 1, 2, 3, 4, …, n-1, where f ^-1 is the inverse function of f in eqn. 7 defined as (eqn. 9):

where t is the age in years and D is the dbh in cm.

The above-mentioned transition probability equation is a direct application of the formula introduced by López et al. ( [9], [10]) for a Bertalanffy-Richards model. We provide the calculations for site indices 20, 17, and 14 in Tab. 1.

According to the model assumptions and the methodology used, we had to consider that the same transition probability p_k applies to all of the individuals in the same class k (where there might exist individuals under very different conditions), i.e., the transition probability p_k is a global property shared by all of the individuals in the same class k, independently of their local conditions. We consider that the stratification into three site index classes must result in low variability around the mean p_k values.

Recruitment and basal area

Pinus nigra regenerates well under its own canopy ([18], [26]), and recruitment may be very abundant even at high densities, as described by Tíscar Oliver et al. ([27]) who reported more than 1875 stems ha^-1 during a regeneration event in Pinus nigra stands with G = 51.96 m² ha^-1. Nonetheless, there is a critical basal area for sufficient light transmittance for regeneration, approximately 20-30 m² ha^-1 ([20], [22]). Therefore, we have considered in the calculations three levels of stand basal area in that range, G = 22, 24, and 26 m² ha^-1. With these basal area values, we assumed the same three levels of recruitment (R) as in our previous studies: (i) 200 stems ha^-1 (scarce recruitment); (ii) 520 stems ha^-1 (normal); (iii) and 840 stems ha^-1 (abundant). For these cases, similar to other matrix models for tree species ([29]), we attributed the global amount of recruitment R entirely to the last diameter class (r₂ = r₃ = … = r_n-1 = 0; r_n = R/[x_n(0)]), which had no influence on the stable diameter distribution of the stand, which is established by transition probabilities, global amount of recruitment, and stand basal area ([9]).

Natural mortalit

We estimated natural mortalities based on previous studies. Specifically, Misir et al. ([16]) developed a mortality model for Pinus nigra subsp. pallasiana in Turkey from a logistic model and the observed data. They estimated an overall annual mortality rate of 1.41%, equivalent to a 13.24% mortality rate in a 10-year period. Trasobares & Pukkala ([28]) studied Pinus nigra Arn. in north-eastern Spain and suggested natural mortality rates of 2% (including minimal intermediate thinning) for the D > 30 cm diameter classes. Therefore, the following 10-year individual tree constant mortality rates for each diameter class were assumed: m₁ = 0.20; m₂ = 0.14; m₃ = 0.08; m₄ = 0.05; m₅ = 0.03; and m₆ = m₇ = … = m_n = 0.02.

“Sustainable/stable” harvesting strategy. Stable diameter distribution

As mentioned above, the “sustainable/ stable” harvesting strategy is aimed at reaching in each harvest the proportion of stems ha^-1 in each class corresponding to the stable diameter distribution of the stand, yielding the following “sustainable/stable” harvest rate shown in eqn. 1, for all of the diameter classes and time steps. This rate corresponds to the case such that the right eigenvector associated with the dominant eigenvalue λ = 1 of the transition matrix with harvests is W₀ for each time step.

By applying the “sustainable/stable” harvesting strategy to each scenario and setting X(0) = W₀ (stable diameter distribution) as the initial condition, the basal area of the stand continuously oscillates between the values G_min (after harvesting) and G_max (before harvesting), and the stable diameter distribution is reached at each time step. The calculations for the 27 scenarios of this study, obtained by combining three rates of diameter growth with three levels of recruitment (R = 200, 520, or 840 stems ha^-1), and three levels of stand basal area (G = G_max = 22, 24, or 26 m² ha^-1), are shown in Tab. 2 along with the population growth rates λ₀, the “sustainable/stable” harvest rates s, and the minimum basal areas G_min.

On the other hand, to compare the economically optimal and the “sustainable/stable” harvesting strategies, the NPV values corresponding to the “sustainable/stable” strategies are given by (eqn. 10):

and the distance between these distributions can be evaluated with Keyfitz’s Δ (eqn. 11):

where x_i and w_i are the components of the diameter distributions X (optimal) and W₀ (stable), scaled to 1 (i.e., Σ x_i = Σ w_i - [6]). The maximum value of Δ is 1, and the minimum of 0 occurs when the diameter distributions X and W₀ (scaled to 1) are identical. Smaller values of Δ result when the distributions X and W₀ (scaled to 1) become more similar.

Stumpage value model

The stumpage value model applied in this study comes from López et al. ([13]). It is a generalization of the model proposed for site index 20 in López et al. ([11]) for Pinus nigra in the Spanish Iberian System. In Tab. 3, we provide the stumpage values we used to obtain the economic data ([17], [28]). The regression models adjusted to these data appear in Fig. 2. The stumpage prices of Pinus nigra in Spain remained fairly stable between 1994 and 2011 ([14]).

The last row in Tab. 3 shows the average price per cubic meter for trees belonging to the different diameter classes. This price is obtained by adding up the four rows above showing the contributions to the final value of a tree of the four industrial destinations according to the percentage of wood of each diameter class going into each product type. The prices of particle board, poles, sawlog, and high-quality sawlog are 6.0, 39.5, 19.6, and 33.2 € m^-3, respectively. It is necessary to point out that the demand for timber poles in Spain is very high and therefore they reach a higher price than larger size sawlogs.

  Results 

The transition probabilities in Tab. 1, the “sustainable/stable” harvest rates s, and the stable diameter distributions W₀ in Tab. 2 were replaced into the proposed discrete optimal control model. Global optimal solutions to each corresponding optimization problem were found using Solver^® by changing the h_k(t) variables for k = 1, 2, …, n, and t = 0, 10, 20, …, 70 years. Using this procedure, we obtained the optimal state x_k(t) and the control variables h_k(t); basal areas before and after each harvest; optimal NPV and the NPV corresponding to the “sustainable/stable” harvesting strategy; and Keyfitz’s Δ for all of the scenarios considered. The discount rates applied fell between 2% and 4%, as is typical of forest management studies.

The main results may be summarized as follows: (a) Fig. 3 shows the population dynamics x_k needed to maximize the NPV function over 70 years with G = 24 m² ha^-1 and i = 3%; (b) Fig. 4 depicts the optimal harvesting strategies h_k(t) for G = 24 m² ha^-1 and i = 3%; (c) the NPV values and the NPV increase between the “sustainable/stable” harvesting strategy s and the optimal harvesting strategy (maximizing NPV function) for i = 3% and different values of site index, R and G are displayed in Tab. 4; (d) Tab. 5 shows the set of values for Keyfitz’s Δ that correspond to the distance between the diameter distribution associated with the optimal harvesting strategy (for t = 70 years) and the stable diameter distribution (W₀) for every period of the planning horizon. It is worth to notice that the values obtained were always very low and yielded Δ = 0 for t = 0 and t = 70, as expected. Moreover, it should be stressed that the transition matrix A_T for the whole harvest cycle from t = 0 to t = 70 years is the following (eqn. 12):

and that λ_T = 1 for all of the scenarios, where λ_T denotes the population growth rate for the whole harvest cycle of T = 70 years. Thus, considering the whole harvest cycle, these optimal harvesting strategies are sustainable. In addition to that, the stable diameter distribution is recovered at the end of the harvest cycle.

The optimal harvesting strategy aimed at maximizing NPV subject to the initial and final stage conditions [X(0) = X(T) = W₀], results in higher NPV figures than the “sustainable/stable” conditions. For i = 3%, this increment was always over 5.36% and reached a maximum of 14.43% for site index 20, R = 840 stems ha^-1 and G = 22 m² ha^-1. The NPV increase was greater for site index 20 than for site index 17 and site index 14.

The influence of small variations in the discount rate or the natural mortality rates on the optimal harvesting strategies was minimal because they did not change the general pattern of the solutions found.

  Discussion and conclusions 

This study proposes a model that combines an objective function comprising the NPV of all the management operations in a planning horizon of T = 70 years with a matrix model for the population dynamics, constrained in the initial and final states to be the stable diameter distribution of the stand, i.e., X(0) = X(T) = W₀. The resulting discrete optimal control model provides a simple tool for heuristic analysis to answer the “what if” questions in relation to forest management plans.

Following López et al. ( [9], [10], [12], [13]), we calculated the transition probabilities considering that p_k is an average global property shared by all of the individuals in the same class k, and that, given the basal area constraints and the low basal area scenarios of this study, stand density has a small influence on diameter growth. Moreover, in the case of constrained optimization matrix models, in which the search for optimal paths throughout the phase space is mainly governed by the objective function, optimal decisions are essentially conditioned by the stumpage prices, i.e., by selective harvesting, and thus stand density has a minor influence on the objective function.

López et al. ([10]) also show that results obtained for λ₀ are sensitive to the width w of the diameter classes, such that the most conservative values of λ₀ and s result from the lowest values of w, which does not allow an individual tree to jump over more than one diameter class in a single projection. The minimum (integer) widths of the diameter classes compatible with this model assumption were w = 6 cm for site index 20, w = 5 cm for site index 17, and w = 4 cm for site index 14. However, because the dominant eigenvalue λ₀ and the stable diameter distribution are not substantially affected by small or even moderate changes in the width of the diameter classes ([10]), we assumed a width of w = 6 cm for our study to facilitate the comparison and presentation of results.

As shown in Tab. 4, the NPV increase between the optimal and the “sustainable/stable” harvesting strategies was related to site quality. The highest NPV increments in absolute terms were obtained for the scenarios corresponding to site index 20, G = 26 and high to medium recruitment (R of 840 and 520 stems ha^-1). Those increments were over 1000 € ha^-1, with NPV values for the optimal strategies over 8000 € ha^-1. However, the NPV increase in relative terms lies between 5.36% and 14.43%, showing a higher increase for higher quality sites and higher recruitments. For site index 20 and R between 520 and 840 stems ha^-1, the NPV increase was always over 14%.

For high quality sites, the highest relative advance on the NPV between the optimal harvesting strategy and the “sustainable/ stable” is reached for medium and high recruitment (R between 520 and 840 stems ha^-1), while for low-quality sites, this relative advantage is higher for high recruitment (R = 840 stems ha^-1).

By applying the optimal harvesting strategies, the maximum value obtained for Keyfitz’s Δ was 0.32, with small deviations from the stable distribution throughout the harvest cycle. The highest Δ values were obtained for site index 20, R = 840 and 520 stems ha^-1. On the other hand, the lowest values of Δ were obtained for site index 14, and R = 200 stems ha^-1 (with Δ always less than 0.125).

In low-quality sites, the optimal management strategy leads to stand distributions that are closer to the stable distribution over the planning horizon than in medium- or high-quality sites. On the other hand, it is in these high-quality sites where the relative increase in the NPV is higher.

Regarding the optimal harvesting strategy (see Fig. 3, Fig. 4), the harvest rates are low for the first four diameter classes throughout the planning horizon. As expected, harvests for the seventh and eighth diameter classes are higher at the beginning of the planning horizon and lower during the last three periods. Most harvests for the fifth and sixth diameter class take place along the planning horizon. For medium and high recruitments, the highest harvest at the end of the planning horizon (t = 60) is performed mainly on the fourth and fifth diameter classes.

It is interesting to note that the optimal harvesting strategy results in a reasonable harvesting schedule, which follows from the dynamics of the stand, initial and final conditions, and the time preference stated by the NPV function, and that the NPV values obtained for the optimal harvesting strategies from this study lie always between those of the optimal (but not sustainable) harvesting strategies shown in López et al. ([13]), where no restrictions were introduced for the final state, and the NPV values for the “sustainable/stable” harvesting strategy.

In conclusion, the proposed optimization model was used in this study to determine economically optimal harvesting strategies of uneven-aged Pinus nigra stands in the Spanish Iberian System, when the population dynamics are governed by a matrix model, the initial and final diameter distributions are constrained to be the stable diameter distribution [X(0) = X(T) = W₀], and a set of constraints is stated to maintain the stand basal area within narrow ranges commonly used in the management of this species. Under such conditions, this study shows that there are sustainable harvesting strategies that retain the stable diameter distribution of the stand at the end of the planning horizon, yielding a higher NPV than the “sustainable/stable” scenario, especially for high-quality sites, although these NPVs are lower than the corresponding values for the economically optimal harvesting strategies.

  References