Introducing a non-stationary matrix model for stand-level optimization, an even-aged pine (Pinus sylvestris L.) stand in Finland

Article Introducing a Non-Stationary Matrix Model for Stand-Level Optimization, an Even-Aged Pine (Pinus Sylvestris L.) Stand in Finland Johanna Pyy 1,*, Anssi Ahtikoski 2, Erkki Laitinen 1and Jouni Siipilehto 3 1Faculty of Sciences, University of Oulu, P.O. Box 8000, FI-90014 Oulu, Finland; [email protected] 2Natural Resources Institute Finland (Luke) Oulu, Paavo Havaksen tie 3, FI-90014 Oulu, Finland; [email protected] 3Natural Resources Institute Finland (Luke), Latokartanonkaari 9, FI-00790 Helsinki, Finland; [email protected] *Correspondence: [email protected]; Tel.: +35-844-545-0982 Academic Editors: Maarten Nieuwenhuis and Timothy A. Martin Received: 10 April 2017; Accepted: 9 May 2017; Published: 11 May 2017 Abstract: In general, matrix models are commonly applied to predict tree growth for size-structured tree populations, whereas empirical–statistical models are designed to predict tree growth based on a vast amount of field observations. From the theoretical point of view, matrix models can be considered to be more generic since their dependency on ad hoc growth conditions is far less prevalent than that of empirical–statistical models. On the other hand, the main pitfall of matrix models is their inability to include variation among the individuals within a size class, occasionally resulting in less accurate predictions of tree growth compared to empirical–statistical models. Thus, the relevant question is whether a matrix model can capture essential tree-growth dynamics/characteristics so that the model produces accurate stand projections which can further be applied in practical decision-making. Such a dynamic characteristic in our model is the basal area of trees, which causes nonlinearity in time. Therefore, our matrix model is a nonlinear model. The empirical data for models was based on 20 sample plots representing 8360 tree records. Further, according to the model, stand projections were produced for three Scots pine (Pinus sylvestris L.) sapling stands (age of 25 years, stand density fluctuating from 850 to 1400 stems ha −1 ). Then, (even-aged) stand management was optimized by applying sequential quadratic programming (SQP) among those growth predictions. The objective function of the optimization task was to maximize the net present value (NPV) of the ongoing rotation. The stands were located in Northern Ostrobothnia, Finland, on nutrient-poor soil type. The results indicated that initial stand density had an effect on optimal solutions—optimal stand management varied with respect to thinnings (timing and intensity) as well as to optimal rotation. Further, an increasing discount rate shortened considerably the optimal rotation period, and relaxing the minimum thinning removal to 30 m 3 ha −1 resulted in an increase both in number of thinnings and in the maximum net present value. Keywords: forest management; nonlinear matrix model; optimization; discrete optimal control 1. Introduction The matrix population models (interchangeably transition matrix models or Usher matrix models; see [ 1 ] and [ 2 ], respectively) are widely used to study the dynamics of forest types around the world [1,3] . In general, matrix population models rely on a division of the diameter distributions into ordered classes [ 4 , 5 ]. With respect to classifying forest dynamics models (depending on their level of description of the forest), matrix population models fall between stand models and individual tree Forests 2017,8, 163; doi:10.3390/f8050163 www.mdpi.com/journal/forests Forests 2017,8, 163 2 of 13 models (e.g., [ 6 , 7 ]). Matrix population models are grounded in distribution-based population models where individual attributes are summarized by their population-level distribution (e.g., [8]). In matrix population models, tree growth is modelled as a transition from a class to upper classes (upgrowth rate), survival as the cumulated transitions from one class to another (mortality rate), and recruitment as a transition into the first class (recruitment rate) [ 2 ]. In principal, matrix models are based on four assumptions: Markov property, Usher property, stationarity and geospatial independence [ 1 ]. Throughout the history of developing matrix models, they have been criticized e.g., for their inability to include variation among the individuals within a size class [ 9 , 10 ], for the Usher property (an individual cannot move up by more than one class or move backwards), and for the arbitrariness of the class division [ 1 , 9 , 11 , 12 ]). Several solutions for these above-mentioned drawbacks have been developed and applied (e.g., [ 4 , 13 , 14 ]), and a few more approaches are under way [ 11 , 15 ]. Despite the few shortcomings associated with matrix population models, they have been applied to almost all the subject areas of forestry [ 1 ]. The advantages of population matrix models, compared to individual-based models (interchangeably empirical–statistical individual tree models, see [ 16 ]) are abundant, but dependent on the application; the best approach for a particular case should be the one that is the most consistent with modelling purposes while making the fewest assumptions according to the law of parsimony, i.e., Occam’s razor [ 1 ]. Stated differently, when the two approaches (individual-based models and population matrix models) make predictions of similar quality, Occam’s razor favours parsimonious population matrix models (see, e.g., [ 17 ]). Given the complexity of individual-based models, the large amount of information (and field measurements) they require and the long processing time still make them difficult and laborious to apply for forest management, thus simpler and more compact models dealing with e.g., size classes are more efficient and practical for the majority of purposes [ 1 , 18 ]. Especially when applying growth models within an economic analysis of forest management, population matrix models have been shown to be advantageous and efficient to apply since various optimization algorithms can straight-forwardly be combined, resulting in a sound and solid assessment framework (e.g., [ 19 , 20 ]). In such a framework, the population matrix model produces growth predictions using the optimization algorithm, which maximizes the objective function (e.g., net revenues; see [20]) by dynamic computing. In Usher matrix models, it is assumed that the growth speed of a tree is stationary, i.e., it depends only on the size of the tree [ 1 ]. In reality, the growth speed depends also implicitly on time. Therefore, we relaxed the stationarity assumption and used a non-stationary matrix model where the growth speed of the tree depends not only on the size class but also on the basal area of the forest stand (basal area evolves in time). The objective of this study was first to introduce a nonlinear matrix model for Scots pine tree-growth in an even-aged stand, and then to maximize the net present value (NPV) of the ongoing rotation by optimizing stand management with sequential quadratic programming (SQP). We applied three advanced seedling stands with varying intensity (i.e., number of stems per hectare) as an initial point for stand projections. Timing and intensity of thinnings and optimal rotation period associated with the optimal stand-level management were analyzed in detail to reveal potential patterns in results with respect to initial stand characteristics. 2. Materials and Methods 2.1. The Growth Model We suppose that the forest stand consists of a single species of trees, which has different diameters. Moreover, we assume that the diameter distribution of trees depends on the basal area. Let us assume that the time is divided into subintervals [tk , tk+1] , k= 0, . . . , M , and the diameter is divided into N discrete and non-overlapping size classes, which are denoted by subindices i= 1, . . . , N . We denote, yk= [yk 1 , yk 2 , . . . , yk N]T and hk= [hk 1 , hk 2 , . . . , hk N]T are the discrete diameter distribution and harvest vector of the forest stand at time event tk , respectively. Moreover, their elements yk i and hk i Forests 2017,8, 163 3 of 13 are the number of trees per unit area and the number of removed trees per unit area of size class i , respectively. Let us denote by bk i and ak i the probability that a tree from diameter class i grows to the next diameter class i+ 1 and the probability that a tree from diameter class i remains in the same diameter class between tk and tk+1 , respectively. Denote also by mk i the probability that a tree from diameter class i dies between tk and tk+1 . The rate ak i is called stasis rate, bk i upgrowth transition rate and mk i mortality rate [ 1 ]. We assume that the upgrowth transition rate bk i and mortality rate mk i depend linearly on the basal area of the forest stand [21], i.e., bk i=c0i+c1iΠ(yk),i=1, . . . , N−1, k=0, . . . , M−1, (1) mk i=d0i+d1iΠ(yk),i=1, . . . , N,k=0, . . . , M−1, (2) where Π(yk) = N ∑ j=1 yk jπ(xj/ 2 )2 is the basal area of the stand and xj is the centre of the diameter class j . The stasis rate ak ican be calculated from the upgrowth transition rate bk iand mortality rate mk i ak i=1−bk i−mk i,i=1, . . . , N,k=0, . . . , M−1, (3) where bk N=0 for all k=1, . . . , M. We assume the following explicit matrix equation for the size class distribution at time tk+1: yk+1=U(yk)yk−hk(4) where U(yk)is the forest growth matrix, which has the following structure at time event tk U(yk) =          ak 10 . . . 0 0 bk 1ak 2. . . 0 0 0bk 2 ...0 0 . . .. . ........ . . 0 0 . . . bk N−1ak N          . (5) The time dependency introduced in this paper contributes to the current literature on matrix growth models (e.g., [1,11,18,22]). Model Parameter Estimation The data used to estimate the model parameters of the size-structured transition matrix were derived from two long-term experiments (HARKAS series; see, e.g., [ 23 ]). The two experiments were established in 1978 and 1984 in even-aged, pure commercial Scots pine (Pinus sylvestris L.) stands located on mineral soil in Ostrobothnia region, Finland. The biological ages of the experimental stands at the time of establishment were 43 and 58 years, respectively. The sites were classified as the Vaccinium forest site type, a sub-xeric forest, a nutrient-poor soil type [ 24 ], which presents app. 25% of forest land area on mineral soils according to the 11th national forest inventory (NFI) in Finland [25]. The stands were established by sowing with seed of local origin (i.e., unimproved seed material). The experiments located in Muhos municipality, 26 ° 06 0 05 00 E and 64 ° 46 0 02 00 N, asl 60–70 m. (Figure 1). Average growing season falls into a range of 100–140 days (threshold +5 ° C, see [ 26 ]). The data for upgrowth and mortality in the matrix model were based on 8360 tree records from 20 sample plots. The stand management among the sample plots fluctuated considerably: from control plots (no thinnings) to very intensive thinnings (60% of the basal area removed). The sample plots of one experiment (altogether 12 sample plots) were measured five times during 1978 and 2014. Accordingly, in another experiment, the sample plots (eight sample plots) were measured four times. The age range Forests 2017,8, 163 4 of 13 of measurements in experiments covers together a time span from 43 to 88 years which practically presents all commercial thinnings of Scots pine during a rotation in the Ostrobothnia region. Figure 1. Locations of experiments. According to the data described above, we estimated the coefficients c0i , c1i , i= 1, . . . , N− 1 and d0i , d1i , i= 1, . . . , N in Equations (1) and (2), respectively, by using the least squares method. Technically, parameter estimations were calculated by using Matlab (R2012a 7.14.0.739, MathWorks Inc, Natick, MA, USA) . 2.2. The Optimization Problem The aim of the optimization was to maximize the revenues from the thinnings. The optimization problem was formulated as follows: max h∈RMN M−1 ∑ k=0 ∑N i=1(cpvp i+csvs i)hk i (1+r)tk−p(hk)!, (6) subject to Equation (4) y0=y0, (7) y≥0, h≥0 (8) where r is the interest rate; cp and cs are the price of pulpwood and saw log, respectively; vp i and vs i are the volume of pulpwood and saw log of a tree in diameter class i , respectively; and y0 is the initial value for size class distribution. The function p(hk) is a penalty term, which ensures that if any thinning is done at time tk, at least Bcubic meters per hectare have to be removed. It is in the form p(hk) =        N ∑ j=1 (vp j+vs j)hk j(B− N ∑ i=1 (vp i+vs i)hk i), if 0 < N ∑ i=1 (vp i+vs i)hk i<B, 0, else. (9) Forests 2017,8, 163 5 of 13 The pulpwood and sawlog volumes of a Scots pine were tabulated at different diameters at intervals of 5 cm starting from 7.5 cm (cf. [ 27 ]). By using table (original values according to pine, H100 = 20 m ) and cubic spline with not-a-knot end condition [ 28 ], we calculated the pulpwood and sawlog volumes of a tree in diameter classes i , i= 1, . . . , N . The sawlog and pulpwood volumes of a tree are shown in Figure 2. 6 8 10 12 14 16 18 20 22 24 26 28 30 32 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 diameter of a tree sawn timber and pulpwood volumes, m3 pulpwood saw log Figure 2. Pulpwood and sawlog volumes of a tree as a function of diameter at breast height, cm. We applied FMINCON solver (Matlab Optimization Toolbox) for solving the optimization problem (6) – (9) . Technical details and formal subproblems associated with SQP are presented in Appendix A. Initial Data for Stand-Level Optimization Data for stand-level optimization were sought from a permanent stand plot database for young stands (TINKA). We searched for Scots pine-dominated advanced seedling stands (dominant height of approximately 8 m) from Ostrobothnia. We calculated the rounded average stand characteristics from the artificially regenerated Scots pine stand on a dryish site. Stand characteristics were needed as input variables for creating the individual tree dimensions, stem diameters at breast height (dbh) and tree heights (h). These stand characteristics represented the average stand density having 1000 stems ha−1 at the age of 25 years. The stand basal area was 9 m ha −1 , basal area-weighted mean diameter ( ∑dbh3/∑dbh2 ) 12 cm and Lorey’s height ( ∑h·dbh2/∑dbh2 ) 7.6 m. We created variation in the number of stems while fixing the other stand characteristics in order to develop the options for one more sparse stand and also for one higher stand density (Table 1). The two-parameter Weibull distribution was solved from the stand characteristics using parameter recovery [ 29 ]. The samples of 20 randomly selected diameters were generated for each option. We sampled trees from the cumulative probability distribution by randomizing the percentile (P) from the uniform 0–1 distribution. The cumulative Weibull distribution function is of form F(dbh) = 1 −e−(dbh/b)c and tree diameter is solved as dbh =b{−ln( 1 −P)}( 1 /c) where b and c are the scale and shape parameter of the Weibull distribution [ 30 ]. From the solved parameters (see Table 1), we noticed that the average density resulted in almost normal distribution ( c= 3.61) while the highest density (1400 ha −1 ) with parameter c of 2.01 is strongly skewed to the right and finally the sparse stand (850 ha −1 ) resulted in a peaked distribution with parameter cof 7.53 (Figure 3). Forests 2017,8, 163 6 of 13 The parameters b0= 1.319 and b1= 0.296 of the Näslund’s height curve [ 31 ] were predicted from stand basal area, basal area-weighted mean dbh and Lorey’s height together with the effective temperature sum (993 ° C d) using models by [ 32 ]. Models also included prediction for the standard deviation of the residual error as a variance function (see [ 32 ]). Random deviation was added to the predicted tree height in order to result in realistic variation in stem dimensions. Table 1. The standard stand characteristics and the recovered Weibull parameters band c. Stand Age, Years Basal Area, m2ha−1 Number of Stems ha−1 Weighted Mean dbh, cm Lorey’s Height, m Dominant Height, m Scale Parameter b Shape Parameter c 25 9 850 12 7.6 8.5 12.216 7.535 25 9 1000 12 7.6 8.5 11.352 3.612 25 9 1400 12 7.6 8.5 9.057 2.011 Forests 2017,8, 163 6 of 13 The parameters b0= 1.319 and b1= 0.296 of the Näslund’s height curve [ 31 ] were predicted from stand basal area, basal area-weighted mean dbh and Lorey’s height together with the effective temperature sum (993 ° C d) using models by [ 32 ]. Models also included prediction for the standard deviation of the residual error as a variance function (see [ 32 ]). Random deviation was added to the predicted tree height in order to result in realistic variation in stem dimensions. Table 1. The standard stand characteristics and the recovered Weibull parameters band c. Stand Age, Years Basal Area, m2ha−1 Number of Stems ha−1 Weighted Mean dbh, cm Lorey’s Height, m Dominant Height, m Scale Parameter b Shape Parameter c 25 9 850 12 7.6 8.5 12.216 7.535 25 9 1000 12 7.6 8.5 11.352 3.612 25 9 1400 12 7.6 8.5 9.057 2.011 Figure 3. The shapes of the Weibull dbh distributions for the average stand (solid line) which followed approximately normal distribution ( c= 3.61) and for sparse (broken line) and dense stand (dotted line). 3. Results In our calculations, we used the following values for model parameters: interest rate r= 3%, price of pulpwood cp= 16.56 em−3 and sawlog cs= 58.44 em−3 , and minimum thinning removal B=50 m3ha−1. The results show that two or three intermediate thinnings took place during rotation (Table 2). Further, it seems that the thinning pattern (with respect to harvested and remaining trees in different size classes) was more or less similar, regardless of the initial stand density: large trees were always removed (Figure 4). With regard to the thinning intensity of the first thinning, app. 43–45% of basal area was removed. In the second thinning, from 45% (normal stand) up to 68% (dense stand) of basal area was removed. In the third thinning (sparse and normal stands), the removal percentage (expressed as relative to basal area) was on average 64%. The removal percentages fall into the original range of thinning intensities in the modelling data. The maximized revenues fluctuated between 3713 and 4198 eha−1, indicating that initial stand density has an effect on the maximum net present value. On the other hand, with regard to MAI, initial stand density has only a minor role (Table 2). Figure 3. The shapes of the Weibull dbh distributions for the average stand (solid line) which followed approximately normal distribution ( c= 3.61) and for sparse (broken line) and dense stand (dotted line). 3. Results In our calculations, we used the following values for model parameters: interest rate r= 3%, price of pulpwood cp= 16.56 em−3 and sawlog cs= 58.44 em−3 , and minimum thinning removal B=50 m3ha−1. The results show that two or three intermediate thinnings took place during rotation (Table 2). Further, it seems that the thinning pattern (with respect to harvested and remaining trees in different size classes) was more or less similar, regardless of the initial stand density: large trees were always removed (Figure 4). With regard to the thinning intensity of the first thinning, app. 43–45% of basal area was removed. In the second thinning, from 45% (normal stand) up to 68% (dense stand) of basal area was removed. In the third thinning (sparse and normal stands), the removal percentage (expressed as relative to basal area) was on average 64%. The removal percentages fall into the original range of thinning intensities in the modelling data. The maximized revenues fluctuated between 3713 and 4198 eha−1, indicating that initial stand density has an effect on the maximum net present value. On the other hand, with regard to MAI, initial stand density has only a minor role (Table 2). Forests 2017,8, 163 7 of 13 Table 2. Optimal stand-level managements generated by the matrix model. Discount rate 3%. Thinning Age of the Stand (a) Volume of Removed Trees (m3ha−1) Proportion of Saw log (%) NPV (eha−1) MAI (m3ha−1year−1 ) Sparse (Number of stems 850) First thinning 45 63.7 75 Second thinning 55 50.0 80 Third thinning 65 52.3 81 Final thinning 95 73.7 90 Total 239.6 82 4042 2.52 Normal (Number of stems 1000) First thinning 45 66.3 78 Second thinning 55 50.0 76 Third thinning 65 56.3 75 Final thinning 90 73.7 84 Total 246.3 79 4198 2.74 Dense (Number of stems 1400) First thinning 45 54.9 78 Second thinning 55 71.4 72 Final thinning 85 95.1 85 Total 221.4 79 3713 2.60 5−8 11−1417−2023−26 29− 0 100 200 300 First thinning 850 5−8 11−1417−2023−26 29− 0 100 200 300 Second thinning 5−8 11−1417−2023−26 29− 0 100 200 300 Third thinning 5−8 11−1417−2023−26 29− 0 100 200 300 Number of trees / ha 1000 5−8 11−1417−2023−26 29− 0 100 200 300 5−8 11−1417−2023−26 29− 0 100 200 300 5−8 11−1417−2023−26 29− 0 100 200 300 1400 5−8 11−1417−2023−26 29− 0 100 200 300 Size class Trees after thinning Removed Figure 4. Diameter distributions associated with optimal thinnings. Numbers (e.g., 5–8, 11–14) represent diameter in centimetres. We tested the sensitivity of the results with respect to two critical aspects. First, the effect of discount rate on optimum stand management was analysed by changing the original 3% into 4% and 5%. Then, we loosened the penalty term (Equation (9) ) so that the minimum thinning removal would be 30 m3ha−1 , which can be considered to be the “decisive limit” for contractors to execute a thinning in Finland (e.g., [ 33 ]). For simplicity, we conducted both sensitivity analyses only for the normal stand density option (number of stems is 1000 ha −1 ). The results of the sensitivity analysis are shown in Table 3. In the sensitivity analyses, the thinning removals ranged from 31% (third thinning related to 3% discounting and 30 m 3 minimum removal criterion) up to 66% (sixth thinning related to 3% discounting and 30 m3minimum removal criterion), expressed as a percentage of basal area. Forests 2017,8, 163 8 of 13 Optimal rotation shortened with increasing discount rate—for instance, with 5% discounting, the optimal rotation was 65 years whereas with 3% discounting it was as long as 90 years (Tables 2and 3). Another interesting result with discount rates was that with 4% and 5% discounting, the first thinning occurred earlier than with 3%, and the number of intermediate thinnings dropped to two with 4% and 5% discounting (Table 3). The differences in MAI between 3%, 4% and 5% discounting were, however, minor. Relaxing the minimum thinning removal criterion to 30 m 3 ha −1 resulted in a slightly higher maximum net present value as well as more intermediate thinnings to be conducted—compared to the baseline optimal solution (Table 3). Thus, from the forest owner’s point of view, milder thinnings (with respect to thinning removals) are favourable: applying more frequent thinnings with relatively small thinning removals increases the maximum net present value (Table 3). Table 3. Optimal stand-level managements generated by the matrix model with different interest rates and minimum thinning removals. “Baseline optimal solution” indicates stand management with 3% interest rate and normal initial density (number of stems 1000 ha−1). Thinning Age of the Stand (a) Volume of Removed Trees (m3ha−1) Proportion of Saw log (%) NPV (eha−1) MAI (m3ha−1year−1 ) Interest rate 4% and minimum thinning removal 50 m3ha−1 First thinning 40 61.4 70 Second thinning 50 50.2 67 Final thinning 75 109.3 78 Total 220.9 73 3168 2.95 Interest rate 5% and minimum thinning removal 50 m3ha−1 First thinning 35 50.2 58 Second thinning 45 50.0 57 Final thinning 65 77.0 68 Total 177.2 62 2513 2.73 Interest rate 3% and minimum thinning removal 30 m3ha−1 First thinning 35 36.9 66 Second thinning 45 38.6 68 Third thinning 55 30.0 80 Fourth thinning 65 36.5 86 Fifth thinning 75 30.2 88 Sixth thinning 85 43.1 88 Final thinning 105 41.2 88 Total 256.6 81 4355 2.44 Baseline optimal solution First thinning 45 66.3 78 Second thinning 55 50.0 76 Third thinning 65 56.3 75 Final thinning 90 73.7 84 Total 246.3 79 4198 2.74 4. Discussion When dealing with modelling, one should bear in mind that it is the end users, not the modelers themselves, who finally determine the value of a model (e.g., [ 1 ]). Thus, from the end users’ point of view, simplicity and accuracy (in prediction) are the key words to emphasize. From the computational and analytical viewpoint, matrix models are actually simpler by structure than empirical–statistical individual tree models [ 1 ], and they essentially require handling of less information [ 8 ]. With respect to accuracy, matrix models have been shown to predict tree growth accurately—both in the short term [ 34 ] and long term [ 18 , 22 ]. This paper introduces a new time-dependent transition matrix model for predicting pine growth in northern Finland, in which the time-dependency originates from modelling the basal area evolving in time. To our knowledge, such a model presents a novel approach in the literature of matrix population models (cf. [ 1 ]). Prior to concluding, we need to compare our results to existing literature on similar growth conditions and tree species in order to discover whether our model is applicable for end users (who are responsible for actual decision-making). Forests 2017,8, 163 9 of 13 For stand-level optimization, there are various different algorithms to choose from [ 35 , 36 ], e.g., the derivative-free direct search method such as the Hooke and Jeeves method, differential evolution, particle swarm optimization [ 37 ], hybrid optimization strategies which combine separate algorithms (e.g., [ 16 ]) or depth-first search algorithms which apply a search tree consisting of a backtracking mechanism [ 37 ]. In this study, we applied a new algorithm which has recently been introduced to forest applications: sequential quadratic programming (SQP) [ 38 ]. Tentatively, the SQP (as representative of gradient-type methods) has been proven to be robust and much more efficient than the derivative-free methods [38]. With respect to financial performance, our results are in line with existing literature on the same tree species (pine) and growth conditions in Finland (e.g., [ 36 , 39 , 40 ])—given the fact that we assessed the maximum net present value (NPV) for the ongoing rotation, whereas the existing literature focuses on assessing the maximum bare land value (BLV). However, these two measures (MaxNPV and MaxBLV) can be technically commensurated for comparison (e.g., [ 41 ]). Optimal number of thinnings varied in this study, depending on the penalty term (Equation (9) ) and particularly on the discount rate. This result is coherent with earlier studies [ 39 , 40 ]) suggesting, e.g., that higher interest rates decrease the optimal number of thinnings. Mean annual increments (MAIs) associated with the optimal stand management were here slightly lower than presented in existing literature [ 40 , 42 ], varying from app. 2.52–3.0 m 3 ha −1 year −1 to 3.6 m 3 ha −1 year −1 . The reasons for this minor discrepancy in volume output are easy to depict: first, the exact locations (in terms of temperature sum and micro-climatic conditions) of stands are slightly different between this study and the others [ 40 , 42 ]). The other reason is related to study frameworks, more precisely to growth models which are, of course, different, and thus result in slightly different outcomes. However, one can argue that the MAIs underlying optimal management produced by alternative growth models can be considered to be similar. In addition, different optimization algorithms were applied in the studies, which also has an effect on outcomes (see, e.g., [35]). Finally, the underlying rationale (and motive) for constructing a time-dependent matrix model is to later on be able to incorporate genetic gains into that particular model. The idea of incorporating genetic gains into the time-dependent matrix model stems from the fact that genetic gain in growth evolves with time (e.g., [ 43 ]), indicating that genetic gain once assessed at juvenile stage could change towards maturity [ 44 , 45 ]. Thus, in the near future, we need growth models which take into account this phenomenon. Further, we could also include variation in stem quality (due to tree breeding; see [ 46 , 47 ]), and further incorporate it into the time-dependent matrix model. For instance, this can technically be done with a subdivision of the categories, i.e., size classes (see [ 4 , 48 ]). Having incorporated the effect of time-evolving genetic gain and variation in stem quality into the time-dependent matrix model, this creates a new set of analysis tools for assessments. In this new assessment framework, it would finally be possible in forestry to value the relevant traits (such as enhanced growth and stem quality) in monetary terms, and to construct genuine trade-offs between breedable traits. This further enables efficient deployment of improved material in different climatic conditions (see [49]). Acknowledgments: We want to acknowledge the Jenny and Antti Wihuri foundation for financial support. Author Contributions: Jouni Siipilehto provided the data for optimizations; Johanna Pyy developed the matrix model and optimization problem and conducted the optimizations under supervision of Erkki Laitinen; Anssi Ahtikoski analyzed the results; all authors contributed to the writing of the manuscript. Conflicts of Interest: The authors declare no conflict of interest. The founding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.