Estimating the optimal rotation age of Pinus nigra in the Spanish Iberian System applying discrete optimal control

This study proposes a model based on discrete optimal control for estimating the optimal rotation age of Pinus nigra in the Spanish Iberian System. The model was constructed using the NPV (alternately LEV) of forest management as objective function and a model of stand evolution represented by a finite difference equation. The optimal rotation ages obtained corresponding to the 2-4% discount rate interval were between 60 and 100 years for the NPV strategy (50 and 70 years for the LEV strategy). With regard to the optimal management schedule, the results have shown the need to intensify the harvesting at 40-60 years in order to maximize the production of posts in view of the high price they command on the market.


Introduction
The concept of optimal rotation age of even-aged forest stands has evolved considerably since it was first introduced by Faustman (1849) and, from a rotation age maximizing the net present value of the perpetual returns from the land, it has come to include new financial aspects (Samuelson, 1976;Plantinga, 1998;Tahvonen et al., 1998;Koskella and Ollikanen, 2001;Clarke and Reed, 2009).
This last concept has been the subject of a large number of contributions, particularly the studies of Clark (1976), Anderson (1976), Hartman (1976), Bhattacharyya et al. (1989) and Ohlin (1995), which have always attempted to maximize economic benefits functions which corresponded as closely as possible to reality.
A wide range of criteria have been used to determine the rotation age for the purpose of maximizing the productivity of the forest stand, and this has given rise to several different methods for determining the optimal rotation age: maximize the physical yield of a single rotation, maximize the annualized yield of a single rotation, maximize the discounted net revenues of a single rotation or from an infinite series of like rotations, maximize the Net Present Value of the investment, etc. (Newman, 1988).This study is based on the optimization of the Net Present Value (NPV) of all the management operations in a rotation discounted to the beginning of the rotation, under constraints defined by the tree population dynamics, and from the landowner's perspective.Alternately, we have also considered the optimization strategy based on the Land Expectation Value (LEV).
Optimal control theory makes possible to solve a wide variety of dynamic problems, where the evolution of a dynamical system (discrete or continuous) through time can be partially controlled by an agent's decision.In every moment t , the system is described by a set of state variables x t (e.g.stems ha -1 at year t) and the planner chooses a set of control variables h t (e.g.harvest rate at year t).Different values of the control variables imply different paths in the phase space for the dynamical system and the manager must determine the values of the control variables that maximize the selected objective, taking into account the constraints of the problem.
Many authors have used optimal control theory for modelling forest stands (e.g., Anderson, 1976;Clark, 1976;Haight, 1985;Synder and Bhattacharyya, 1990;Sethi and Thompson, 2000;Termansen, 2007).Most of these models are based on the theory of continuous optimal control.However, when dealing with a biological system (e.g. a forest), it is often relevant not to observe the system behavior at all instants of time t but only at a sequence of instants, e.g.t = 0, 10, 20,… Commonly in such cases it is possible to characterize the system behavior by expressions involving those instants, in particular relating each instant to the previous.In fact, for managed tree population, the harvest operations generally occur in regular and predictable cycles (e.g. 10 years).Thus discrete time models are well suited to describe the corresponding dynamics.The discrete optimal control has been less commonly applied to determine optimal forest harvesting schedules and/or optimal rotation ages, due to the large number of decision variables involved (Kronrad and Huang, 2002;Tahvonen, 2004).Therefore, increasing attention has been given to computational methods and software packages for working with the complex problems of discrete optimal control.
In this context, the main purpose of this study is to estimate the optimal rotation age and the corresponding optimal harvesting schedules of even-aged Pinus nigra stands in the Spanish Iberian System.To examine these issues, we have combined a model of the stand evolution with an objective function involving the Net Present Value (NPV) of all the harvesting operations in a rotation, giving rise to a discrete optimal control problem.Alternately, we have also considered the optimization strategy based on the Land Expectation Value (LEV).The results obtained show that, under the harvesting strategy based on the NPV, the optimal rotation ages for Pinus nigra stands in the Spanish Iberian System are between 60 and 100 years, corresponding to timber devoted to posts (between 50 and 70 years for the LEV strategy).
In contrast to other current similar methodologies (Palahí and Pukkala, 2003;Trasobares and Pukkala, 2004), whose results are based on overall economic variables, this study highlights the variation in stock in the stand at each point of the harvesting schedule, linked to the iterative processes of calculation used to obtain the rotation age which maximizes the NPV or LEV.
The management of Pinus nigra stands in the study area is currently based on the growth and yield tables for Pinus nigra Arn. in the Spanish Iberian System (Gómez Loranca, 1998), which depend mainly on site index classes (Qualities I to V) and on the intermediate harvesting schedule (moderate or intense).Under these management conditions, the optimal rotation age lies somewhere between 105 years for quality I to 120 years for quality IV (Gómez Loranca, 1998).

Study area and data
The study area is located in the Spanish Iberian System, a mountain range extending about 400 km along the north-eastern edge of the central plateau, concentrating 60% of the area occupied by Pinus nigra in Spain, mainly in the provinces of Cuenca, Guadalajara and Teruel (Grande et al., 2004).
These forests support various ecological and economic functions, but the stands in our study area are predominantly devoted to timber production.Therefore, the main management objective is timber production.
Input data used to build the model have been obtained from the yield tables for Pinus nigra Arn. in the Iberian System (Gómez Loranca, 1998).The data in these tables are classified according to growth level into five site qualities, of which only the first has been selected as our data source.This Quality I level was defined by the dominant height reached at the reference age of 60 years (20 meters).
The thinnings always began at year 30, and the initial stand conditions are shown in Table 1.

The model
Considering that the main objective of the stands in the study area is the production of timber, the objective Optimal rotation Pinus nigra by discrete optimal control function of the optimization problem is the Net Present Value (NPV) of all the management operations in a rotation discounted to the beginning of the rotation, and from the landowner's perspective.Alternately, we have also considered the optimization strategy based on the Land Expectation Value (LEV).The combination of each one of these objective functions with the model of the population dynamics introduced below gives rise to a discrete optimal control problem, from which we must determine the control variables h t that maximize the function , or, alternately, the function , where x t , for t = 30, 40, 50,…, T, is the number of stems ha -1 at year t (according to the growth and yield tables of Gómez Loranca (1998), the time step between consecutive thinnings adopted in the model was 10 years); h t for t = 30, 40, 50,…, T, is the harvest rate at year t; v t is the stumpage value at year t (€ stem -1 ); i is the discount rate; and T is the rotation age in years.It is well known that both the NPV and LEV are justified management objectives if a single economic goal is selected (Siegel et al., 1995;Kangas et al., 2008).
The population dynamics are given by where t = 30, 40, 50,…, T. The initial conditions are specified in Table 1.
Natural regeneration is assumed, so costs associated with replanting are negligible in the model.On the other hand, since Gómez Loranca's (1998) tables begin at age 30 (years), from this moment the natural mortality rates are commonly stabilized, taking small (and constant) values.Therefore they do not influence the results of the optimization process.Additionally, we may consider that the natural mortalities are included in the harvest rates.
The following sections provide a description of the model inputs used and the estimates of the corresponding parameters.The general calculations, including regressions, were run using Maple version 13 (Maple, 2009) software, and the optimization problems were solved using Solver Premium Platform version 7.1 (Front Line Systems, 2008) running under Excel 2007.

Growth model
The growth model was obtained using regression analysis from data from the growth and yield tables of Gómez Loranca (1998), yielding where id10 represents the quadratic mean diameter increment (in centimeters) for each 10-year time step, D the quadratic mean diameter (in centimeters), G the stand basal area (in m 2 ha -1 ) and ln the Neperian logarithm.Two different thinning regimes moderate and strong were used in the adjustment.The main characteristics of this model are shown in Table 2, and Figure 1 shows the data and the model adjusted.

Stumpage value model
The stumpage value model applied in this study has been proposed for Qualities I, II and III in López et al. (under review) for Pinus nigra in the Spanish Iberian System.Table 3 shows the figures for the stumpage  value used to obtain the economic data (Montero et al., 1992;Trasobares and Pukkala, 2004;and own data).Gómez Loranca's (1998) volume expressions were used for the conversion between € (m 3 ) -1 of timber and € stem -1 .The price differences between high quality sawlog and sawlog are due to the following circumstances: (a) The timber devoted to high quality sawlog (finally posts) is highly demanded because its specific characteristics and unique destiny, whereas timber destined to sawlog has very different destinations, thus it has an average lower value; and (b) high quality sawlog (post) has greater demand in the market because its durability, due to its high resin contents, in contrast to sawlog, that has alternatives in other species, e.g.P. sylvestris.
The regression model adjusted to these data is where v t is the stumpage value function (€ stem -1 ), D t is the quadratic mean diameter (centimeters) at year t (with 0 ≤ D t ≤ 50).Note that D t+10 = D t + id10.The main characteristics of this model are shown in Table 4.

Quadratic mean diameter after vs. before thinning
The quadratic mean diameter of the stand increases slightly in each thinning, as the stems with the lowest diameters are harvested in each time step.The differences are very small, as shown by the following regression model for the quadratic mean diameter after thinnings (D a ) (in centimeters) in relation to the quadratic mean diameter before thinnings (D b ) (in centimeters) and harvest rate (h), adjusted for each time step of 10 years using data from the growth and yield tables of Gómez Loranca (1998): Optimal rotation Pinus nigra by discrete optimal control 309   The main characteristics of this model are shown in Table 5.

Discount rate
International Financial Reporting Standard 41 (IASB, 2004), in applying the expectation value approach to forest valuation, when referring to the process of deriving the present value of expected net cash flows for the forest stand (paragraphs 20 to 23), prescribes the use of an appropriate discount rate in each case, and specifies this as a «current market determined pre-tax rate».But there are a large number of current market determined pre-tax rates.For example, the USDA Forest Service has typically used i = 4% for long-term forest service investments (Row et al., 1981), and in the United Kingdom, HM Treasury has recommended the use of a declining long-term discount rate (3% for 31-75 years, 2.5% for 76-125 years, 2% for 126-200 years, 1.5% for 201-300 years and 1% for more than 300 years) (HM Treasury, 2007).Thus, to consider the incorporation of eventual risks, we have initially chosen in this study a wide range of discount rates, from 2% to 4%.

Optimization method
For each discount rate, the software Solver Premium Platform version 7.1 (Front Line Systems, 2008) 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 increases, the probability that the globally optimal solution has been found also increases towards 100%.
The control variables h t which maximize the objective function NPV or LEV described above, for each of the cases T = 40, 50, 60,…, 100,… (and always under the population dynamics and the constraints mentioned), are the solutions to the corresponding problems (one optimal solution for each T value), and the optimal rotation age is defined by the value of T which gives the highest value of NPV or LEV of all those obtained previously.Finally the optimal harvesting schedules are defined by the combination of the control variables h t corresponding to this value of T.

Results
Applying the procedure described in the previous section, the optimal rotation ages were obtained for each type of discount above mentioned (2-4%).The results are summarized in Table 6.In each of the previous cases, Table 7 shows the combinations of control variables h t which give the maximum value of the objective function, under the corresponding population dynamics and constraints (optimal harvesting schedules).
The optimal rotation ages obtained show a wide range of variation, between 60 and 100 years for the NPV stra-  tegy, and between 50 and 70 for the LEV strategy.When the discount rate is less than 2%, the optimal rotation age clearly exceeds the optimal rotation age of 105 years established in the study of Gómez Loranca (1998) (for Quality I).In fact, when i = 1.5%, for example, the optimal rotation age was T = 130 years for the NPV strategy.
As for the optimal management schedule, the results obtained in this study suggest a thinning model which differs from the proposals of Gómez Loranca (1998).This difference lies in the harvesting intensity of the model analyzed, which consists of a slight initial decrease from t = 40-50 years followed by a stabilization for the basal area of the stand in relation to its age, in order to maximize the production of posts in view of the high price they command on the market.And this decrease in the basal area, from 35-45 m 2 ha -1 at t = 40-50 years to 24-27 m 2 ha -1 at the end of the rotation for the NPV strategy, depending on the discount rate, is the reverse of what is contemplated in the tables of Gómez Loranca (continual growth from 32.6 m 2 ha -1 initially to 55 m 2 ha -1 at the end).

Discussion
This study proposes a model that combines an objective function comprising the NPV (alternately, the LEV) of all the management operations in a rotation with a model of the population dynamics defined by a first order finite difference equation.This model gives rise to a discrete optimal control problem that enables us to study the optimal rotation age and the optimal management schedule for even-aged Pinus nigra stands in the Spanish Iberian System.
The model of diameter growth was adjusted, using two different thinning regimes, moderate and strong, based on the tables of Gómez Loranca (1998).Although for diameters of less than 14 cm or higher than 45 cm these growth model probably do not have a good degree of fit, in all the results obtained in this study the diameters were between about 18 cm (for very high discount rates) and 37.27 cm (for i = 2%).Likewise, the basal areas were always clearly over 20 m 2 ha -1 , threshold below which natural regeneration of the stand may not always be possible (Schütz, 1997).On the other hand, since Pinus nigra is classif ied as shade intolerant, equally important to conserve this minimum coverage up to the rotation age, is to keep the stand density within reasonable bounds.
Regarding the model used for the value of the standing timber, this value is known to remain fairly stable over time, as shown in Figure 2 (MAPA, 2008).Though the possible variation in future stumpage prices has traditionally not been considered in Central Europe  (Sagl, 1995), given the obtained rotations, it is possible that timber prices may be altered (see for instance Sohngen and Sedjo (2005), for the American market).Therefore, we have verified in the context of this study that a change in timber prices, in the narrow range defined by Figure 2, does not affect the obtained optimal rotation periods.
As for the discount rates, the analysis carried out appears to suggest the interval 2-4% as a typical interval for use, since below these levels the discounted incomes are very small.This interval is similar to that proposed by other authors (e.g., Díaz Balteiro, 1997).This discount rate does not include the incorporation of possible risks by using specific indicators as in other studies (e.g., forest fires, as in González-Olabarría et al., 2008).
The results obtained regarding the optimal rotation age are between 60 and 100 years for the NPV strategy (between 50 and 70 for the LEV strategy), for discount rates between 2% and 4%.These figures for i < 2% (and/or lower quality levels), are clearly higher (e.g., T = 130 for i = 1.5% and the NPV strategy).It is well known that optimal rotation ages under 70 years, in general are not technologically viable due to the diffi-culty to commercialize timber of such dimensions.In a global economic context with very low discount rates, the possible increase in the rotation age derived from these results, which was also suggested by Gómez Loranca in his study, although based fundamentally on biological considerations, should be supported by the following considerations: a) The graph of the NPV function in relation to T, as shown in Figure 3, is relatively flat on the right side of its maximum (however, note that the graph of the LEV is very different from the graph of the NPV; and the low discount rates are more suitable for long time horizons, and therefore, for multiple rotations).
b) The price of sawn timber is greater for diameters of over 45 cm.c) Pinus nigra maintains its productive and reproductive capacity until an advanced age, around 150-170 years.Previous studies (Tiscar, 2004(Tiscar, , 2007) ) assign Pinus nigra a longevity of 500-600 years (although in Spain, trees of up to 900 years old have been inventoried), during which time there was no appreciable reduction in its reproductive capacity.
d) The diameter class of over 45 cm is obtained in most sites from 130 years of age (Gómez Loranca, 1998).One of the main consequences of the possible lengthening of the optimal rotation age, under the clearing regime suggested by these results, is ecological, as this extension might favour a greater conservation of the stand (Grande et al., 2004;Schütz, 1997).This extension of the optimal rotation age has also been proposed by other authors, although based on forest management criteria (Grande et al., 2004;Schütz, 1997) or biological considerations (Gómez Loranca, 1998), and not generally on criteria of economic optimization.
In conclusion, these results appear to suggest that, in a global economic context with very low discount rates, the optimal rotation age of even-aged Pinus nigra stands in the Spanish Iberian System should be extended, and the thinning schedules should be re-examined considering in particular the evolution of the Pinus nigra Arn.timber prices.

Figure 1 .
Figure 1.Quadratic mean diameter increment (id10) in relation to the quadratic mean diameter (D) and the stand basal area (G) (id10 and D in centimeters, G in m 2 ha -1 ).

Table 1 .
Initial conditions

Table 4 .
Stumpage value model v t = D t a 8 exp(a 9 + a 10 D t ) R 2 : coef. of determination.R 2 adj : adjusted coef. of determination.D-W: Durbin-Watson statistic.D in centimeters and v in euro stem -1 .

Table 3 .
Stumpage Prices [€ (m 3 ) -1 ] of Pinus nigra in Spain according to its diameter class and industrial destination.Note that timber devoted to high quality sawlog (post) is obtained from rotation ages considerably lower than the corresponding to sawlog.This explains the difference between their present values, in favour of high quality sawlog

Table 5 .
Characteristics of the model for the quadratic mean diameter after thinnings (D a ) (centimeters) in relation to the quadratic mean diameter before thinnings (D b ) (centimeters) and harvest rate (h), D a = a 11 + a 12 D b + a 13 h + + a 14 D b 2 + a 15 h 2 + a 16 D b h, for D b < 30 and h ≠ 0, D a = D b , otherwise Durbin-Watson statistic.D in centimeters and v in euro stem -1 .

Table 6 .
Optimal rotation age T (years) in relation to the discount rate i.Optimal rotation ages under 70 years, in general are not technologically viable due to the difficulty to commercialize timber of such dimensions

Table 7 .
Optimal rotation Pinus nigra by discrete optimal control 311 Optimal harvesting schedule: harvest rates h t which maximize the objective function (NPV or LEV) described in section 2.2., for each discount rate.
Figure 2. Evolution of the stumpage prices of Pinus nigra in Spain between 1994 and 2008.