This study proposes a discrete optimal control model to obtain harvest strategies that maximize the net present value (NPV) of the timber harvested from unevenaged Pinus nigra stands located in the Spanish Iberian System, between two stable positions. The model was constructed using an objective function that integrates financial data on the harvesting operations with a matrix model describing the population dynamics. The initial and final states are given by the stable diameter distribution of the stand, and the planning horizon is 70 years. The scenario analysis corresponding to the optimal solutions revealed that the stand diameter distribution does not deviate substantially from the equilibrium position over time and that the NPV for the optimal harvesting schedule was always greater than the NPV for the “sustainable/stable” harvesting strategy. The NPV increase for the different scenarios is between 5.36% and 14.43%, showing a greater increase in higher site index scenarios and higher recruitments.
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  Schütz 2006). These processes of evolution over time can be easily modeled using a matrix population model (e.g., Buongiorno & Michie 1980, Caswell 2001, Liang & Picard 2012, Vanclay 2012), 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_{0}. 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 (López et al. 2007, López et al. 2008, López et al. 2012, López et al. 2013) 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):
s = \frac{ \lambda_{o}  1 }{ \lambda_{o} }
where λ_{0} 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_{0}.
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. (2013) showed that the economically optimal harvesting strategies underlying the solutions of that model changed the stand towards distributions more characteristic of evenaged 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 unevenaged 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_{0}. 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^{2} 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 (Haight 1985, Sethi & Thompson 2000) and in discrete case (Kronrad & Huang 2002, Tahvonen 2004, López et al. 2010, López et al. 2013).
In the following sections, we describe the model and the estimates of the corresponding parameters. The general calculations, including regressions, employed the software Maple (2014) version 18, and the solutions to the optimization problems were derived using the Solver^{®} tool running under MS Excel 2013^{®}.
Materials and methodsDiscrete 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 northeastern edge of the Central Plateau, including 60% of the area occupied by Pinus nigra in Spain (Grande & García Abril 2005). 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 (Gómez Loranca 1998). 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. ( López et al. 2007, López et al. 2008) for Fagus sylvatica L., and by López et al. ( López et al. 2012, López et al. 2013) 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 BertalanffyRichards model (Von Bertalanffy 1938, Richards 1959); (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 shadetolerant species, and the gaps opened in the canopy through singletree 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: 06, 612, 1218, …, 3036, …, 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 (Caswell 2001). 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):
and I is the identity matrix; H(t) = diag[h_{1}(t), h_{2}(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 (Siegel et al. 1995).
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_{0}. Therefore, the first constraints are given by (eqn. 5):
X(0) = X(T) = W_0
We also set the following obvious constraints (eqn. 6):
m_{k} \le h_{k}(t) \le 1
for k = 1, 2, … , n and t = 0, 10, 20, 30, … , T  10, where m_{k} is the natural mortality in kth class.
On the other hand, by applying the “sustainable/stable” harvesting strategy to each scenario and setting X(0) = W_{0} 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^{2} 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 estimationTransition probabilities
We calculated the transition probabilities by applying the method introduced in López et al. (2012) 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 BertalanffyRichards function (Von Bertalanffy 1938, Richards 1959  eqn. 7):
D = f(t) = a(1e^{bt})^c
adjusted using regression analysis to data points [t, D(t)] obtained from tables by Gómez Loranca (1998). The main characteristics of these growth models are shown in Fig. 1.
Using these diametergrowth 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):
where t is the age in years and D is the dbh in cm.
The abovementioned transition probability equation is a direct application of the formula introduced by López et al. ( López et al. 2007, López et al. 2008) for a BertalanffyRichards 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 (Retana et al. 2002, Tíscar 2007), and recruitment may be very abundant even at high densities, as described by Tíscar Oliver et al. (2010) who reported more than 1875 stems ha^{1} during a regeneration event in Pinus nigra stands with G = 51.96 m^{2} ha^{1}. Nonetheless, there is a critical basal area for sufficient light transmittance for regeneration, approximately 2030 m^{2} ha^{1} (Schütz 1989, Serrada et al. 1994). Therefore, we have considered in the calculations three levels of stand basal area in that range, G = 22, 24, and 26 m^{2} 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 (Van Mantgem & Stephenson 2005), we attributed the global amount of recruitment R entirely to the last diameter class (r_{2} = r_{3} = … = r_{n1} = 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 (López et al. 2007).
Natural mortalit
We estimated natural mortalities based on previous studies. Specifically, Misir et al. (2007) 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 10year period. Trasobares & Pukkala (2004) studied Pinus nigra Arn. in northeastern Spain and suggested natural mortality rates of 2% (including minimal intermediate thinning) for the D > 30 cm diameter classes. Therefore, the following 10year individual tree constant mortality rates for each diameter class were assumed: m_{1} = 0.20; m_{2} = 0.14; m_{3} = 0.08; m_{4} = 0.05; m_{5} = 0.03; and m_{6} = m_{7} = … = 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_{0} for each time step.
By applying the “sustainable/stable” harvesting strategy to each scenario and setting X(0) = W_{0} (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^{2} ha^{1}), are shown in Tab. 2 along with the population growth rates λ_{0}, 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):
where x_{i} and w_{i} are the components of the diameter distributions X (optimal) and W_{0} (stable), scaled to 1 (i.e., Σx_{i} = Σw_{i}  Keyfitz 1968). The maximum value of Δ is 1, and the minimum of 0 occurs when the diameter distributions X and W_{0} (scaled to 1) are identical. Smaller values of Δ result when the distributions X and W_{0} (scaled to 1) become more similar.
Stumpage value model
The stumpage value model applied in this study comes from López et al. (2013). It is a generalization of the model proposed for site index 20 in López et al. (2010) for Pinus nigra in the Spanish Iberian System. In Tab. 3, we provide the stumpage values we used to obtain the economic data (Montero et al. 1992, Trasobares & Pukkala 2004). 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 (MAGRAMA 2012).
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 highquality 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_{0} 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^{2} ha^{1} and i = 3%; (b) Fig. 4 depicts the optimal harvesting strategies h_{k}(t) for G = 24 m^{2} 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_{0}) 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):
A_{T} = \prod_{t=T10}^{t=0} A[I  H(t)]
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_{0}], 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^{2} 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_{0}. 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. ( López et al. 2007, López et al. 2008, López et al. 2012, López et al. 2013), 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. (2008) also show that results obtained for λ_{0} are sensitive to the width w of the diameter classes, such that the most conservative values of λ_{0} 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 λ_{0} and the stable diameter distribution are not substantially affected by small or even moderate changes in the width of the diameter classes (López et al. 2008), 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 lowquality 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 lowquality sites, the optimal management strategy leads to stand distributions that are closer to the stable distribution over the planning horizon than in medium or highquality sites. On the other hand, it is in these highquality 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. (2013), 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 unevenaged 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_{0}], 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 highquality sites, although these NPVs are lower than the corresponding values for the economically optimal harvesting strategies.
ReferencesBuongiorno J, Michie BRA matrix model of unevenaged forest management. Forest Science 26: 609625.1980Caswell HMatrix population models: construction, analysis, and interpretation (2nd edn). Sinauer Associates Inc., Sunderland, MA, USA, pp. 710.2001Gómez Loranca JAModelo de crecimiento y producción para Pinus nigra Arn. en el Sistema Ibérico. [Growth and yield model for Pinus nigra Arn. in the Iberian System]. Montes Revista de Ámbito Forestal 54: 518.[in Spanish]1998Grande MA, García Abril ALos pinares de Pinus nigra Arn. en España: ecología, uso y gestión. [The forests of Pinus nigra Arn. in Spain: ecology, use and management]. Fundación Conde del Valle de Salazar  FUCOVASA, Madrid, Spain, pp. 699. [in Spanish]2005Haight RGA comparison of dynamic and static economic models of unevenaged stand management. Forest Science 31 (4): 957974.1985Keyfitz NIntroduction to the mathematics of population. AddisonWesley, Reading, MA, USA, pp. 450.1968Kronrad GD, Huang CHFinancially optimal thinning and final harvest schedules for loblolly pine plantations on nonindustrial private forestland in east Texas. Southern Journal of Applied Forestry 26 (1): 1317.2002Liang J, Picard NMatrix model of forest dynamics: an overview and outlook. Forest Science 59 (3): 359378.2012López I, Ortuño S, Martín AJ, Fullana CEstimating the sustainable harvesting and the stable diameter distribution of European beech with projection matrix models. Annals of Forest Science 64: 593599.2007López I, Fullana C, Ortuño SF, Martín AJChoosing Fagus sylvatica L. matrix model dimension by sensitivity analysis of the population growth rate with respect to the width of the diameter classes. Ecological Modelling 218: 307314.2008López I, Ortuño S, Martín A, Fullana CEstimating the optimal rotation age of Pinus nigra in the Spanish Iberian System applying discrete optimal control. Forest Systems 19 (3): 306314.2010López I, Ortuño S, García F, Fullana CIs De Liocourt’s distribution stable? Forest Science 58 (1): 3446.2012López I, Ortuño S, García F, Fullana CAre the economically optimal harvesting strategies of unevenaged Pinus nigra stands always sustainable and stabilizing? Forests 4: 830848.2013MAGRAMABoletín mensual de estadística 19942011. [Monthly statistical bulletin 19942011]. Technical General Secretary of the Spanish Ministry of Agriculture, Food and Environment, Madrid, Spain, pp. 48. [in Spanish]2012MapleMaplesoft version 18.0. Waterloo Maple Inc., Ontario, Canada.2014Misir M, Misir N, Yavuz HModeling individual tree mortality for Crimean pine plantations. Journal of Environmental Biology 28 (2): 167172.2007Montero G, Rojo A, Alía RDeterminación del turno de Pinus sylvestris L. en el Sistema Central. [Determination of rotation age for Pinus sylvestris L. in the Central System]. Revista Montes 29: 4247. [in Spanish]1992Retana J, Espelta JM, Habrouk A, Ordóñez JL, SolàMorales FRegeneration patterns of three Mediterranean pines and forest changes after a large wildfire in northeastern Spain. Ecoscience 9: 8997.2002Richards FJA flexible growth function for empirical use. Journal of Experimental Botany 10 (29): 290300.1959Schütz JPLe regime du jardinage [The gardening regime]. Chaire de sylviculture, ETH Zürich, Switzerland, pp. 110. [in French]1989Schütz JPModelling the demographic sustainability of pure beech plenter forests in Eastern Germany. Annals of Forest Science 63: 93100.2006Serrada R, Domínguez Lerena S, Sánchez Resco MI, Ruiz Ortiz JEl problema de la regeneración natural de Pinus nigra Arn. [The problem of natural regeneration of Pinus nigra Arn.]. Revista Montes 36: 5257. [in Spanish]1994Sethi SP, Thompson GLOptimal control theory: applications to management science and economics (2nd edn). Kluwer Academic Publishers, Boston, MA, USA, pp. 505.2000Siegel W, Hoover W, Haney H, Liu KForest owners’ guide to the federal income tax. Agriculture Handbook No. 708, USDA Forest Service, Washington, DC, USA, pp. 138.1995Tahvonen OOptimal harvesting of forest age classes: a survey of some recent results. Mathematical Population Studies 11: 205232.2004Tíscar PADinámica de regeneración de Pinus nigra subs. Salzmannii al sur de su área de distribución: etapas, procesos y factores implicados. [Regeneration dynamics of Pinus nigra subs. Salzmannii south of its distribution area: stages, processes and related factors]. Investigación Agraria: Sistemas y Recursos Forestales 16(2): 124135.2007Tíscar Oliver PA, Lucas ME, Tíscar Soria MALa alteración del suelo y la espesura como factores de regeneración de Pinus nigra subs. Salzmannii a lo largo de su área de distribución. [Soil alterations and stem density as regeneration factors for Pinus nigra subs. Salzmannii along its distribution area]. Revista Montes 103: 1015. [in Spanish]2010Trasobares A, Pukkala TOptimising the management of unevenaged Pinus sylvestris L. and Pinus nigra Arn. mixed stands in Catalonia, northeast Spain. Annals of Forest Science 61: 747758.2004Van Mantgem PJ, Stephenson NTThe accuracy of matrix population model projections for coniferous trees in the Sierra Nevada, California. Journal of Ecology 93: 737747.2005Vanclay JKModelling continuous cover forests. In: “Continuous cover forestry” (Pukkala T, Gadow KV eds). Managing Forest Ecosystems 23, Springer, Berlin, Germany, pp. 229241.2012Von Bertalanffy LA quantitative study of organic growth (Inquiries on growth laws, II). Human Biology 10: 181213.1938
Diameter growth models (dbh D in cm and time t in years). I: Site index 20, D = 51.68 (1  e^{0.015259t})^{1.255111}; II: Site index 17, D = 46.645633 (1  e^{0.014318t})^{1.337062}; III: Site index 14, D = 40.644134 (1  e^{0.013838t})^{1.456382} . The coefficient of determination was always greater than 0.999.
Stumpage value models (stumpage values v in € stem^{1} and diameter D in cm). (I) Site index 20: v = D^{3.186471} exp(7.704952  8.678687 · 10^{3}D); (II) Site index 17: v = D^{3.114196} exp(7.476506  9.903125 · 10^{3}D); (III) Site index 14: v = D^{2.987053} exp(7.110977  1.078752 · 10^{2}D). The coefficient of determination was always greater than 0.9999.
Population dynamics x_{k} to maximize the NPV function throughout a harvest cycle of 70 years for G = 24 m^{2} ha^{1} and i = 3% (left: Site index 20; middle: Site index 17; right: Site index 14; first row: R = 200 stems ha^{1}; second row: R = 520 stems ha^{1}; third row: R = 840 stems ha^{1}). The initial and final conditions are defined by the stable diameter distributions.
Optimal harvesting strategies h_{k}(t) for G = 24 m^{2} ha^{1} and i = 3% (left: Site index 20; middle: Site index 17; right: Site index 14; first row: R = 200 stems ha^{1}; second row: R = 520 stems ha^{1}; third row: R = 840 stems ha^{1}). The dashed lines depict the “sustainable/stable” harvest rates s.
Transition probabilities (p_{k}) between diameter classes for each site index.
Transition
Site index 20
Site index 17
Site index 14
(06) → (612)
p_{1} = 0.7697
p_{1} = 0.5951
p_{1} = 0.4564
(612) → (1218)
p_{2} = 0.8602
p_{2} = 0.6824
p_{2} = 0.5326
(1218) → (1824)
p_{3} = 0.7913
p_{3} = 0.6200
p_{3} = 0.4697
(1824) → (2430)
p_{4} = 0.6828
p_{4} = 0.5190
p_{4} = 0.3692
(2430) → (3036)
p_{5} = 0.5533
p_{5} = 0.3971
p_{5} = 0.2475
(3036) → (3642)
p_{6} = 0.4106
p_{6} = 0.2618
p_{6} = 0.1119
(3642) → (4248)
p_{7} = 0.2587
p_{7} = 0.1171

(4248) → (48,→)
p_{8} = 0.1000


Numerical values for λ_{0}, s, and stable distributions (G in m^{2} ha^{1}, Est. Dis. in stems ha^{1}).
Stumpage prices (€ m^{3}) of Pinus nigra in Spain according to its diameter class and industrial destination.
Products
Diameter classes (cm)
<20
2040
>40
Particle
6
1.2
0.6
Poles
0
11.85
0
Sawlog
0
9.8
11.76
Highquality sawlog
0
0
9.96
Total average price
6
22.85
22.32
NPV values (€ ha^{1}) for the “sustainable/stable” and the optimal harvesting strategies (Max. NPV) and NPV increase (%) for i = 3%, and different levels of site index, R, and G (G in m^{2} ha^{1}).
Site Index
G
Strategy
R = 200 stems ha^{1}
R = 520 stems ha^{1}
R = 840 stems ha^{1}
Value(€ ha^{1})
Increase(%)
Value(€ ha^{1})
Increase(%)
Value(€ ha^{1})
Increase(%)
20
22
Sust/Stable
5347.3646
12.48
6118.0328
14.33
6383.5522
14.43
Max. NPV
6014.7773
6994.8646
7304.4446
24
Sust/Stable
5743.3336
12.24
6611.1260
14.36
6918.9460
14.35
Max. NPV
6446.4525
7560.6334
7911.8264
26
Sust/Stable
6130.4531
11.93
7096.2876
14.35
7447.4339
14.32
Max. NPV
6861.7018
8114.2591
8514.1776
17
22
Sust/Stable
4239.2159
9.69
4986.6137
12.75
5284.1332
13.84
Max. NPV
4649.9201
5622.2967
6015.6138
24
Sust/Stable
4542.7283
9.31
5373.9477
12.61
5710.9061
13.65
Max. NPV
4965.6346
6051.5239
6490.3456
26
Sust/Stable
4838.9236
8.97
5754.1367
12.51
6131.0587
13.44
Max. NPV
5272.7622
6473.7901
6955.3522
14
22
Sust/Stable
3202.5660
6.20
3906.4753
10.12
4224.8230
11.44
Max. NPV
3401.1319
4301.8564
4708.1096
24
Sust/Stable
3421.6868
5.77
4194.9510
9.74
4548.7650
11.36
Max. NPV
3619.2075
4603.4830
5065.3056
26
Sust/Stable
3635.0537
5.36
4477.2298
9.39
4866.5857
11.27
Max. NPV
3829.9429
4897.7062
5415.2350
Set of values for Keyfitz’s Δ corresponding to the distance between the diameter distribution associated with the optimal harvesting strategy and the stable diameter distribution (given by W_{0}) for every period of the planning horizon (note that Δ = 0 for t = 0 and t = T = 70).