{'fraction_non_alphanumeric': 0.028525668773266386, 'fraction_numerical': 0.016254003982339192, 'mean_word_length': 4.388966643340331, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Optimal Management of Naturally Regenerating Uneven-aged Forests

A shift from even-aged forest management to uneven-aged management practices leads to a problem rather different from the existing straightforward practice that follows a rotation cycle of artificial regeneration, thinning of inferior trees and a clearcut. A lack of realistic models and methods suggesting how to manage uneven-aged stands in a way that is economically viable and ecologically sustainable creates difficulties in adopting this new management practice. To tackle this problem, we make a two-fold contribution in this paper. The first contribution is the proposal of an algorithm that is able to handle a realistic uneven-aged stand management model that is otherwise computationally tedious and intractable. The model considered in this paper is an empirically estimated size-structured ecological model for uneven-aged spruce forests. The second contribution is on the sensitivity analysis of the forest model with respect to a number of important parameters. The analysis provides us an insight into the behavior of the uneven-aged forest model.


Introduction
Optimizing the use of forest resources has hundreds of years of history. The infinite horizon model specified in Faustmann (1849), reintroduced by Samuelson (1976) and extended in numerous papers like Kao and Brodie (1979) and Chen et al. (1980) served as a cornerstone both in research and in practical forestry applications. In its generic form the model determines optimal forest rotation, i.e., the length of optimal interval between clearcuts. However, it is seldom noticed that since this model is most suitable for plantations (Yoshimoto and Shoji, 1998), it has directed research to forests that actually cover only 7% of the total world forest land area. Our study presents major progress in the line of research that serves developing the management of more natural forest stands that have great potential in solving several pressing problems related to forest environment.
The alternative to plantations is to rely on native tree species, natural regeneration and continuous forest cover, i.e., to manage forests as heterogeneous uneven-aged systems. The rationale of this model depends on tree species, but for shade-tolerant trees the economic outcome may become fully competitive because of natural regeneration and more accurate targeting of cuttings to those trees that are financially mature. Additionally, managing forest resources in more natural and heterogamous state has high potential in coping with problems such as climate change (Field et al., 2014), loss of biodiversity and landscape esthetics (Thompson et al., 2009). Multi-criteria decision making approaches have also been used in the past to meet multiple objectives in forest management problems (Steuer and Schuler, 1978;Nhantumbo et al., 2001). Interest in continuous cover forestry is increasing in Nordic countries and UK, for e.g., in Finland it has been released from a 70 year of legislation ban from the beginning of 2014. According to surveys a major problem among forestry professionals is the suspense of the alternative system's economic viability (Valkonen and Cheng, 2014).
While the Faustmann approach describes a chain of exactly similar even-aged cohorts, the model for more natural forests includes the internal structure of heterogeneous trees. As shown in the seminal paper by Adams and Ek (1974) this expands model dimensions and the development of the research has been a struggle against limitations in computing capacity. This has led researchers to develop various simplifications with the cost of losing economically and mathematically sound theoretical structure as already surveyed by Getz and Haight (1989). Most studies still circumvent the problems by studying the fundamentally dynamic problem in a static setup with limited scientific progress and low practical credibility. One problem is in solving a multiple state variable infinite horizon model from any initial stand state. A straightforward solution for this problem was given already in Haight and Monserud (1990): lengthen the planning horizon until the approach path towards a stationary state (or cycle) becomes invariant from further lengthening and consider it as an approximation of the infinite horizon solution. Given a single tree species cases this leads to solvable problems with e.g., 200 periods and 24 optimized variables per period; albeit non-linearities and non-convexities require special attention. Besides dimensionality the other problem is that the optimal solution becomes cutting stand every period which does not make sense in the presence of fixed harvesting cost and the fact that too small yield is commercially invaluable. Fixed harvesting cost is taken into account in even-aged models, like Tahvonen et al. (2013), with the implication that optimal number of intermediate cuttings is between zero and five periods (25 years) depending on factors like site fertility and interest rate. However, in the even-aged problem all rotations are similar implying that the time horizon in computation is relatively short (40-150 years) and the number of combinational variables are usually six or lower. In uneven-aged models Haight and Monserud (1990) take this into account by allowing cuttings every 20 years only. Wikström (2000) includes fixed harvesting cost and computes solutions using tabu search under the simplification that regeneration is fixed to 50 trees per 5 years period and that stand volume is not allowed to decrease below a level determined by Swedish forest legislation. He does not interpret his results on harvesting interval but it seems to vary between 5 and 20 years without any systematic pattern. In Tahvonen (2011) the model includes fixed harvesting cost which leads to optimal harvesting period of 15-20 years under the constraint that the interval is constant over time.
Given these studies, the proper solution method and most general solutions for the uneven-aged management problem are still open. This is pressing in the practically most important cases where the initial forest state is a consequence of even-aged management and the problem is to solve optimal path or transition to uneven-aged management. This question has been studied in numerous works with specifications without full generality. In this paper we make a two-fold contribution. As the first contribution we develop a computational method for solving uneven-aged stand management problems that is a large scale mixed integer non-linear program; and as the second contribution we provide an analysis for the uneven-aged stand management model. The solution method for handling the problem is based on the following: 1. A two-level approach with genetic algorithm at the upper level and continuous nonlinear programing at the lower level: The approach is faster by more than an order of magnitude in terms of computation time as compared to branch-and-bound method. This supports handling of large scale uneven-aged management problems. 2. Modeling the infinite time horizon uneven-aged stand management problem into a tractable problem by assuming transition and steady states: The assumption causes no loss of generality as the transition and steady state lengths are assumed to be endogenous subject to optimization. For an earlier study on forest management practices where the time horizon is divided into transition and steady states the readers may refer to Salo and Tahvonen (2003).
A faster algorithm allowed us to perform a number of computational studies by varying the parameters in the uneven-aged stand management problem. This provided us an insight into the behavior of the uneven-aged model. These insights may play a significant role in directing future research on uneven-aged management.
The later part of the paper is structured as follows. In Section 2 we discuss the sizestructured stand model and introduce the net present value maximization problem. This is followed by the description of the algorithm in Section 3 that is used for solving the optimization problem. Thereafter, in Section 4 we present the results and provide comparisons against the standard approaches that are used to solve uneven-aged stand management problems. Finally, the conclusions are provided in Section 5, where we also highlight the future research directions on uneven-aged stand management. The paper includes appendices that provides additional computational results.

Size-structured Forestry model
The model being considered in this paper is a discrete infinite time horizon model that involves two kinds of variables that are listed below:  A forest management strategy is shown in Figure 1 that we want to optimize for maximum net present value (N P V ) over an infinite horizon in a discrete time framework. The time stages are represented as t = 0, 1, 2, . . . on an infinite time horizon. Harvesting stages are represented by δ t that takes values 0 or 1 with 1 denoting that harvesting is done and 0 denoting that no harvesting is done at a given time stage. The forest states and the extent of harvests are represented with vectors x t and h t respectively. We discuss forest land of one hectare. Larger forest areas require minor and straightforward modifications which we omit. The size-structured forestry model defined in this section utilizes a number of symbols that are described in the discussions. For ease of reference we have also provided these symbols in Table 1.
Trees in the forest are subdivided into a finite number of size classes s for s = 1, 2, . . . , n in increasing order. Let x st be the number of trees in size class s at stage t and define vector x t = (x st ). For t = 0, x 0 is the given initial state of the forest. Let vector h t = (h st ) denote the level of harvesting at stage t. Component h st is the number of trees harvested in size class s at time stage t. For all t, let δ t be a binary variable indicating whether harvesting takes place at stage t (δ t = 1) or not (δ t = 0). Then a logical requirement for harvesting levels is Before stating the optimization problem we introduce a number of endogenous auxiliary variables concerning forest dynamics and cash flow. Martin Bollandsås et al. (2008) use Norwegian National Forest Inventory data and estimate uneven-aged models for most common Nordic tree species. In our study we use their data for Norway spruce.

Other symbols
Ingrowth related: S 1 , S 2 , γ, ν, B 0 Mortality related: m Growth related: A 1 , A 2 , S, L Discounting related: r Revenues related: P 1 , P 2 , C 1 , C 2 Given basal area b s of a tree in size class s, the total basal area (per hectare) at stage t is and the total basal area (per hectare) of trees in size classes larger than s is Ingrowth φ t of trees in step ∆ into the smallest size class 1 as a function of basal area B t is where S 1 , S 2 , B 0 , γ and ν are positive parameters. Mortality µ st is the share of trees dying in size class s in one step ∆. As a function of basal area B t , it is given by where m is a positive parameter, and M s = exp(2.492 + 0.02d s − 3.2 × 10 −5 d 2 s ) depends on tree diameter d s in each size class s. For size class s, given basal area B t and the basal area B st of trees in size classes larger than s, the share of trees growing in one step ∆ from size class s to s + 1 is where G s depends on site index S and latitude L as follows G s (S, L) = 0.02(17.839 + 0.0476d s − 11.585 × 10 −5 d 2 s + 0.906S − 0.268L), while A 1 and A 2 are positive parameters. In this notation, the share of trees remaining in size class s is 1 − µ st − α st .
The gross revenue R t at any time step t is given as where v 1s is the tree volume of small diameter logs in size class s, p 1 is the unit price of small diameter logs, v 2s is the volume for saw timber in size class s and p 2 is the unit price for saw timber. Cutting and hauling cost C t depend on the overall volumes in each size class v s = v 1s + v 2s . Let v = (v s ) be the vector of volumes, then C t is given as follows: where C c is the cost function for logging cost and C h is the cost function for hauling cost, and fixed cost is given by C f . If β denotes the annual discount factor per year, then the problem is to find x t ≥ 0, h t ≥ 0, binary variables δ t as well as auxiliary variables B t , B st , φ t , µ st , α st , R t and C t , for all t and s, to subject to (1)-(8) and forest dynamics state equations x s+1,t+1 = α s,t x s,t + (1 − µ s+1,t − α s+1,t )x s+1,t − h s+1,t for all t and s < n (11) The functions C c (h t , v) and C h (h t , v) used in the model have been defined below. The harvesting cost is specified following the estimation results in Nurminen et al. (2006) assuming that in uneven-aged management cutting costs are 15% higher than in even-aged clearcut operations (Surakka and Sirén, 2007), while hauling cost is determined as in even-aged thinning operations (Tahvonen, 2011).
where C 1 and C 2 represent the cutting cost per minute and hauling cost per minute. It is noteworthy that the hauling cost has an inherent fixed term. For ease of discussions in a later section we denote hauling cost as , where the first term is the fixed term and the second term is the variable hauling cost.
The parameter values employed in our numerical illustrations are given in Tables 2-3. Table 2 provides the size class dependent parameters in the model, i.e basal area (b s ), diameter (d s ), and tree volumes (Heinonen, 1994) of small diameter logs (v 1s ) and saw timber (v 1s ). It also provides three different initial states of the forest that have been studied in this paper. The three initial states x 0 = x 1 , x 2 , and x 3 represent an old even-aged stand, unevenaged stand and young even-aged stand respectively. Table 3 provides other parameters of the model that are independent of the size class. All prices and costs are given at the level of the year 2011. Estimated costs correspond to average costs by hectare of large enough stands. Note also that economies of scale may be taken into account by varying the fixed cost.
Given that dimension of cash flow R t − C t is e, basal areas B t and B st are m 2 , volume v s is m 3 , and step size ∆ is years, parameter dimensions in Tables 2-3 are implied by (1)-(9).  All the parameters and functions have been kept fixed as suggested above. However, we have performed certain sensitivity studies by varying the parameters r, C f and S.

Proposed Algorithm
Finding an optimal forest management strategy requires the optimization of binary as well as continuous variables, which makes the problem in the previous section a mixed integer non-linear programing problem. Mixed integer programing commonly arises in optimization of forest harvesting operations and researchers have used both heuristics (Weintraub et al., 1994) and exact methods (Goycoolea et al., 2005;Constantino et al., 2008;Carvajal et al., 2013) to handle the problems in different contexts. Both the approaches have its own advantages and disadvantages. An exact method, in such cases, would guarantee an optimal solution but might be computationally intractable for large scale problems. On the other hand a heuristic might be computationally tractable, but does not guarantee optimality.
In this section, we describe a two level approach that is customized to solve unevenaged forest management problem. The harvesting decisions are determined (δ t ) at level 1 of our solution procedure that fixes the values of the binary variables at each time stage to 0 or 1. Once the binary variables are fixed, at level 2 we determine the corresponding optimal state of the stand (x s = (x st )) and the extent of harvests (h s = (h st )) at each time stage. The harvesting strategies at level 1 are generated using an evolutionary algorithm and their corresponding optimal continuous variables are determined at level 2 by solving a non-convex optimization problem. The proposed method performs this process iteratively with an intelligent update of the harvesting strategies and leads to a near optimal solution. This section is divided into two parts: the first part provides a detailed description of a customized evolutionary procedure to optimize harvesting strategy and also describes how we convert the infinite time horizon problem into a finite time horizon problem; the second part involves a discussion about the non-convex optimization problem being solved for each harvesting strategy.

Level 1: Evolutionary Optimization
In this section we describe the evolutionary optimization algorithm that uses principles from biological evolution to move towards the optimum by generating improved harvesting strategies. It is a population based approach where each member represents a harvesting strategy. The quality of a harvesting strategy is measured using level 2 that returns the maximum net present value for the given strategy. We refer the net present value corresponding to each harvesting strategy as the fitness of the strategy. The technique emphasizes better harvesting strategies in the population leading to a rise in average fitness over generations (iterations). We use a genetic representation to code harvesting strategies in our algorithm. The algorithm begins with random initialization of a population representing different harvesting strategies. The strategies are then improved over generations by repetitively applying Selection, Crossover, Mutation and Replacement that is described later in this section. For few earlier studies where evolutionary computation techniques have been used in the context of forest management, the readers may refer to Bayat et al. (2013), for instance.

Genetic Representation
It is commonly observed in discrete-time dynamic systems that the variables go through a transition phase and eventually stabilize into a steady state where a fixed pattern gets repeated. Taking insights from such behavior of discrete-time dynamic systems and our prior experience in solving forest management problems (Salo and Tahvonen, 2003), we start with an assumption that the optimal solution to the forest management problem on an infinite horizon consists of transition period and steady state period. During the transition period the state and the control variables keep changing without following any regular pattern. However, during the steady state period both the state and control variables change in a cyclic manner. Therefore, we model an infinite time horizon with a finite set of variables as shown in Figure 2. In the figure, the variables vary over time during the transition phase and once the steady state begins, the harvesting decisions, quantity harvested, as well as the states of the forest follow a repeated cycle. Such a construct allows us to explore a limited time horizon leading to a significant reduction in the search region. Using this representation, any harvesting strategy with transition and steady periods can be coded with finite binary variables. Within the algorithm the solution shown in Figure 2 will have a genetic representation as {(01000100), (10001000)}. The number of bits in the transition period denotes the transition period length, and the number of bits in the steady state period denotes the cycle length. The length of the transition period and steady state cycle is automatically adapted by our algorithm.

Population Initialization
The initial population consisting of N members is generated randomly by the algorithm. The minimum and the maximum transition length (t min , t max ) and steady state cycle length (s min , s max ) are provided as input by the user. For a given member, the algorithm decides the transition length (t len ) by generating a random integer between t min and t max . Thereafter, it generates a genetic string of length t len with 0 or 1 appearing at each location with equal probabilities. A genetic string of length s len for steady state cycle is generated in a similar manner, and then the two strings are combined to obtain a harvesting strategy for the member. This operation is repeated for every member in the population leading to a wide variety of harvesting strategies.

Fitness Assignment
Each harvesting strategy or population member generated in the algorithm has to be evaluated in terms of the maximum net present value that can be obtained with the strategy over an infinite period of time. Identifying this maximal value itself is a non-convex optimization problem that we solve at level 2. The non-convex optimization problem is solved with respect to variables x t and h t leading to their optimal values corresponding to the given harvesting strategy (δ t ). It is noteworthy that we need to solve a non-convex optimization task for every new harvesting strategy that we generate making the overall task computationally demanding. The non-convex optimization is discussed in detail in Subsection 3.2.

Genetic Operators
The genetic operators are used in an evolutionary algorithm to generate new members using the existing population members. The existing members that are used to generate new members are referred to as parents while the newly generated members are referred to as offspring. The genetic operators consist of the steps of crossover and mutation. Crossover involves information sharing between two or more parents while mutation is performed on a single member and is helpful in maintaining diversity. The crossover mechanism used in this study involves two parents that lead to two offspring. The crossover operator is applied on the parents with a probability of p c . In case no crossover is performed between the parents then the offspring are considered to be identical to the parents. Figure 3 shows the crossover operation between two parents with unequal transition and steady state cycle lengths. To perform a crossover a random crossover point is chosen on the member with a smaller transition length. We fix this crossover point on both the members and swap the genetic string that follows the crossover point to generate two offspring. Next, the offspring generated using this operation undergoes mutation.
There are two kinds of mutation that we perform in this study. In type 1 mutation, we go through each genetic bit and flip it from 0 to 1 or 1 to 0 with a probability of p m . This introduces diversity in harvesting strategies and helps in exploring the search space that is not reachable only by sharing of information between the parents. Type 2 mutation is performed to explore an appropriate transition length and cycle length. Type 2 mutation is applied first on the transition string with probability p m . Thereafter, it is applied on the steady state cycle string with probability p m . When this operation is applied, the length of the transition period or steady state period increases or decreases by 1 with equal probabilities. In case the length increases then a bit (0 or 1) is inserted at a random location. If the length decreases then a bit is removed from a randomly chosen location. The type 1 and type 2 mutation operations have been shown in Figure 4 and 5 respectively.

Algorithm Description
Next, we provide a step-by-step description of the algorithm that is used to explore the harvesting decisions and progress towards the optimum. S 1: Initialize N population members (Refer Section 3.1.2) representing different harvesting strategies. S 2: Assign fitness (Refer Section 3.1.3) to each member by solving a non-convex optimization problem at level 2. S 3: Initialize a generation counter: g ← 0. S 4: Increment the generation counter by 1: g ← g + 1. S 5: Choose two parents from the population using tournament selection (Miller and Goldberg, 1995). S 6: Perform genetic operations (Refer Section 3.1.4) to produce the final offspring. S 7: Evaluate each offspring by solving a non-convex optimization problem at level 2. S 8: Choose λ random members from the population, and pool them with the offspring. S 9: Find the best λ members from the pool based on the fitness values and update the population by replacing the chosen members in the previous step with the best λ members. S 10: If the generation counter (g) is less than the maximum number of allowed generations (g max ) then go to Step 4; otherwise terminate.

Parameters
The optimization task at level 1 requires two kinds of parameters. The first set of parameters are the evolutionary parameters, and the second set of parameters are the problemspecific parameters. The parameter setting used in the experiments performed in this paper are given in Table 4.

Level 2: Optimization for fitness evaluation
Our fitness evaluation tasks are non-convex optimization problems which need to be solved fast. In this section we define the problem, discuss the sources of non-convexities and finally propose an initialization procedure for optimization with good prospects for the solver to find a global optimum.

The problem
Let t 0 = t len denote the end of transition (beginning of first cycle) and t 1 = t len + s len the end of first cycle (beginning of second cycle). Given t 0 , t 1 and δ t , for all t < t 1 , the fitness evaluation problem is to find x t ≥ 0, for all t ≤ t 1 , h t ≥ 0 and the auxiliary variables, for all t < t 1 , to subject to (1)-(8), (10)- (14), for t < t 1 , and Here (16) is the steady state condition for state variables. In the objective function (15), the first term accounts for the transition phase and the second term accounts for the steady phase. The multiplier in front of the summation in the second term accounts for a geometric series of discounted steady state cash flows.

Non-convexities
If harvesting takes place at stage t, then the variable cutting and hauling cost (e) in (8) where vector v is given in Table 3. The graph of such concave cost function in Figure 6 (left inset) with respect to h t is almost linear. In (4), ingrowth of trees into size class 1 as a function of basal area B t (m 2 /ha) is φ t = (S 1 (B t + B 0 ) −ν )/[1 + S 2 exp(γB t )] where S 1 , S 2 , γ and ν are given in Table 3. Figure 6 (right inset) illustrates the ingrowth function for B t ∈ [0, 20] m 2 ; for B t > 10 m 2 the function is mildly non-linear.
In (5), mortality fraction of trees dying in size class s is given by µ st = 1/[1+M s exp(−mB t )] where m is given in Table 3. Mortality µ st is a convex function for relevant domain of B t . In the forest state equations (10)-(11) the reduction in the number of trees is µ st x st , where x st is the number of trees in size class s. Figure 7 shows such reduction for three size classes s, s = 1, 5, 10, as a function of B t and x st . Note that possible cases are those for which B t ≥ b s x st , where b s is the basal area of a single tree in size class s. For such region, the mortality functions are almost linear.
For size class s < 12, given basal area B t and the basal area B st of trees in size classes above s, the fraction of trees that grows in five years from size class s to s + 1 is α st = G s −A 1 B st −A 2 B t where G s , A 1 and A 2 are given in Table 3. In (10)-(11) of the optimization problem the transition in terms of number of trees is α st x st , which involves a bi-linear function.  region again the functions are almost linear. Given small values of ρ we may term the nonconvexities as relatively mild. We further confirm this (refer to Section 3.2.3 and Appendix A) by solving the non-convex optimization procedure with different random initialization schemes i.e., choosing the starting point for the solver.

Random initialization for solver
In Appendix A we test seven alternative versions of random initialization (starting point selection) for the solver Knitro (Byrd et al., 2006) to find a local optimum for the problem (15). All seven procedures performed quite well with a high chance for Knitro to end up with a global optimal solution. Thereby, we chose a randomized starting point for the solver as described in the following. Let s denote a size class such that harvesting only applies for size classes i > s. With s = 5 exogenously given, let η it denote the share of trees that are not harvested in size class i at stage t. Then η it = 1 for i ≤ s; otherwise η i+1,t = it η it where it is drawn from uniform distribution U (0, 1). Thus the share of trees harvested in size class i at stage t is (1 − η it ) and it increases in random proportions with i, for i > s. Shares it are drawn independently for each time stage t and size class i > s. The number of trees harvested is h it = δ t (1 − η it )x it for which we need the forest state vector x t unless δ t = 0. Initially, x 0 = x 0 is given and we obtain h i0 = δ 0 (1 − η i0 )x i0 as well as cash flow C t from (7)-(8). Thereafter, auxiliary variables in (2)-(6) and state equations (10)-(11) yield x 1 . Similarly, in forward recursion we obtain values x t , h t , C t and the auxiliary variables for all t. In this random initialization procedure we neglect the steady state condition x t 0 = x t 1 .
Let p denote the probability that a random initialization procedure for the solver ends up with a suboptimal solution. If the problem is initiated independently k times, then the probability of not finding the global optimum is 1 − p k . In our tests reported in Appendix A, p appears to be small (of the order of 0.1 %). Hence, 1 − p k is very close to 1 even for k = 1. This observation may be explained by relatively mild non-convexities in our valuation problems; see nonlinear function illustrations in Figures 6-8.

Results
In this section, the proposed algorithm has been applied on an empirically estimated size-structured ecological model for uneven-aged Norway spruce forests (Martin Bollandsås et al., 2008). We have considered three different initial states of the forest, as stated in Table 2. The first initial state (x 0 = x 1 ) is a young even-aged stand, the second (x 0 = x 2 ) is an uneven-aged stand, and the third (x 0 = x 3 ) is an old even-aged stand. All the forest parameters and functions in our study are fixed as given in Tables 2 and 3. Sensitivity studies have been performed by varying the parameters: interest rate (r), fixed cost (C f ) and site index (S). The different values of these parameters considered are given in Table 4. It was observed from our runs that the optimal steady state harvesting interval and size distribution is independent of the initial stand state of the forest. Steady state results corresponding to the three initial states and different values of parameters r and C f are shown in Table 6. Site index, S, was kept fixed at 15 in this study. Some of the important observations drawn from these runs have been outlined below: 1. The optimal solution converges to the steady state harvesting interval within 120 years (24 intervals) or sooner (Table 6). 2. The transition towards the steady state tree size distribution is shortest when the initial stand state is already heterogeneous and longest for the young even-aged stand. 3. The length of steady state interval varies between 10 and 25 years and increases or remains the same for higher levels of fixed cost. 4. Within the Faustmann optimal rotation framework higher interest rate implies shorter rotation suggesting that the length of the steady state interval decreases in the interest rate. However, this is not the case: for example, when fixed cost is equal to e100 and interest rate is increased from 3% to 4% harvesting interval lengthens from 10 to 15 years. 5. When interest rate increases it is optimal to allocate a larger fraction of capital from forestry to alternative sources. This is reflected in the fact that when harvesting interval lengthens with interest rate, it becomes optimal to cut smaller size class (i.e., size class with average diameter equal to 225mm) and the stand volume both before and after harvest decreases. 6. Given any fixed cost level the stand volume before and after harvest is lower or the same with higher interest rate.

Dependence of Optimal Solution on Initial State
The optimal solution depends strongly on the initial stand state (Figure 9). Given an initially dense stand (solid line) it is optimal to almost clearcut after 10 years and then wait 50 years until the next harvest. In contrast when the initial stand is already heterogeneous ( Table 2) the steady state harvesting interval (15 years) is optimal immediately from the beginning even if it takes about 140 years to reach the exact steady state tree size distribution and harvest level. Given an initial stand containing trees only in the smallest size class, it is optimal to wait 25 years until the first harvest and it takes 80 years to reach the steady state harvesting interval.

Dependence of Optimal Solution on Interest Rate and Fixed Costs
In Figure 10 the initial stand contains only young trees and the first cutting is postponed to 30, 25 and 20 years when the interest rate is increased from 1 to 4%. Varying the fixed cost level above e100 has strong effects on optimal harvest timing ( Figure 11). Given an initially young stand, interest rate equal to 3% and e100 fixed cost, there is a 45 year period with harvest every 5 years, while with fixed cost equal to e500 it is always optimal to wait at least 15 years before the next cut. Fixed cost must be decreased to e20 until the steady state solution is to harvest the stand every period.

Dependence of Optimal Solution on Site Index
Site index has as strong effect on the the optimal transitioning and steady state that can be observed in Figure 12. Given that initial stand contains trees only in the smallest size class, the time of first harvest is at year 20, regardless of site index. The stand density at the first harvest; however, differs between the sites; with more productive sites the overall density of the stand is higher compared to less productive sites. Similarly with other initial stand states, the timing of the first harvest is the same regardless of site index. However, the timing of the second and subsequent harvests during the transition as well as in the steady state differ between the site types. At less productive sites the harvests occur less frequently and target smaller size classes compared to more productive sites. During transition, the more productive sites are kept at higher densities, but in the steady state both the harvest timing and size of harvested trees are the same with site indices 15 and 17. However, annual profit and harvested volume are higher at the more productive site due to the higher growth rate (refer to Table 7).  Figure 13 shows the results for an initially uneven-aged stand (x 0 = x 2 ). We observe that the stand develops close to the steady state values in 95 years. When steady state harvesting interval is three periods, the steady state harvest is targeted to three oldest size classes. However, during the transition phase with varying harvesting interval, trees are harvested from size classes 3-5. Economic optimization tends to yield solutions where harvesting applies to larger trees and no smaller trees are harvested when natural regeneration is insufficient. The economic gain from optimizing the harvesting interval against a fixed harvesting interval was also compared for different cases. The gain was found to be remarkable when the initial state is far from the steady state. For cases shown in Figure 10, an optimal harvesting schedule increases the objective function value by about 10 − 12%.

Other Results
All the results presented in this section so far are based on the computations performed using our proposed algorithm. The computational efficiency of the approach allowed us to perform large number of runs for a variety of cases, that would have been difficult otherwise. In order to gain confidence in the results achieved by our algorithm, we evaluate our approach against the standard branch-and-bound algorithm that is commonly used to handle mixed integer programs.

Comparison of Evolutionary Optimization with Branch-and-Bound
In Appendix B we demonstrate how the problem can be solved exactly using mixed integer programing with branch-and-bound. In Table 8 we provide a comparative evaluation of our heuristic evolutionary optimization algorithm (EO) against the exact branch-andbound approach (BB) for the base case of 3% interest rate (r) and e 300 fixed cost C f . Comparisons have been drawn in terms of the best objective value achieved (1000e), run time (hours) and number of non-convex optimization calls made. The evolutionary optimization algorithm was terminated when a maximum number of non-convex optimization calls exceeded 8000, while the branch-and-bound approach was terminated either when the maximum number of active nodes in the method exceeded 100,000 or the method exited normally. In Table 8, termination of EO due to maximum number of non-convex optimization calls is denoted as NCO, termination of BB due to maximum number of active nodes is denoted as NAN, and termination of BB at confirmed optimum is denoted as NOR. For further information about the mixed integer programing solved using branch-and-bound, the readers may refer to Appendix B.
The results in Table 8 suggest that evolutionary optimization approach offers an enormous computational advantage over the branch-and-bound approach. In terms of the best objective value achieved, evolutionary optimization algorithm was found to outperform branch-and-bound marginally for initial state x 2 . In this case BB was terminated by the upper limit on the number of active nodes. For initial states, x 1 and x 3 , the difference in the best NPV found by EO and BB is truly insignificant. The run time results are not surprising given that branch-and-bound performs an exhaustive search in the decision space.

Conclusions and Future Work
In this paper, we make a two-fold contribution to the uneven-aged forest management literature. The first contribution is on the method side where we propose a customized algorithm to handle mixed integer non-linear optimization problem on optimal management of naturally regenerating uneven-aged forests. For comparison with our heuristic approach we use branch-and-bound, which is a method for which optimality can be proven. Our  approach is able to computationally surpass branch-and-bound by more than an order of magnitude, which would facilitate solving large scale forest management problems. A fast algorithm enabled us to make a more significant second contribution that is building a better insight for the uneven-aged forest model. A thorough analysis of an empirical model for uneven-aged spruce forest stands provided the following important results contributing towards uneven-aged management literature: • Optimizing the harvesting interval yields up to 12% increase in net present value income compared to a solution where the harvesting interval is fixed. The efficiency gain is highest when interest rate and fixed cost is high and initial stand is young.
• Fixed cost was found to have major implications on the optimal solution. For a given initial size class distribution, and interest rate, varying the fixed costs can lead to very different optimal solutions.
• The stationary state of the stand was found to be independent of the initial stand state.
• Stationary harvesting interval tends to be longer with higher fixed cost but may lengthen or shorten with interest rate.
• Stationary harvesting interval is reached much sooner than stationary size class structure.
• A parametric study for the site index that denotes land fertility suggests that with higher yields the harvesting interval is shorter.
The result on efficiency gain when compared with fixed interval harvesting could be an important factor for making a shift towards uneven-aged forest management practices, given that uneven-aged management also offers the advantage of biodiversity preservation. Other results on interest rate study, fixed cost study and site index study provided in the paper would be useful in directing the future research on uneven-aged forest management. The proposed method would be directly applicable for multi-species forests as well, and we intend to study optimal management of such forests in our future work.

Appendices
Computational tests for fitness evaluation and formulation for mixed integer programing In the appendices we discuss tests with alternative random initialization procedures for the fitness evaluation problem (Equation 15) and provide mixed integer programing formulation for the problem that has been solved using branch-and-bound algorithm (Section 4.4.1).

A. Testing random initialization procedures
For tests we considered seven versions of random initialization procedures. Each one has some charter for solving the problem such as harvesting applies primarily to mature trees (high size classes). For versions other than the initialization explained in Section 3.2, the smallest size class subject to harvesting is varied, for instance, it may be randomly drawn. Also the strategy for the shares of trees harvested varies: e.g., harvesting applies to largest sizes classes until a target level of harvested volume is met. Such target level may be tied to the volume growth of the forest adjusted by a random perturbation. In our tests all such initialization procedures performed well in the sense that suboptimal solutions were infrequently produced by the solver.
Initial forest states and harvesting schedules used for tests are illustrated in Tables 9-10  Table 9 defines two sets (A and B) of three initial states of the forest and Table 10 shows eight alternative cases of harvesting schedules.  Each of the seven methods was run with 1000 independent repetitions for each of the 2 × 3 × 8 = 48 test problems. Table 11 shows the best objective function value found which Table 10: Eight cases of harvesting schedules. Interpretation of the three character case code xyz is as follows: x is l (s) for long (short) transition period, y is l (s) for long (short) steady state cycle, and z is d (s) for dense (sparse) harvesting intervals. The length of transition period is t 0 (steps) and s is the steady state cycle (steps). For time stages in five year intervals, 1 indicates harvesting and 0 no harvesting. we now call the optimum. It is interesting to note that for set B of 24 problems each of the seven methods found the optimal solution. The same is true for set A with one exception: for the case e-ssd the optimal value is 14.684 but three of the methods ended up with best value of 14.680. Because the difference is only 0.03 % and given data precision, we may regard the value 14.680 optimal as well. In 34 out of the 48 test problems all seven methods never produced a suboptimal solution in 1000 trials. Here suboptimal means a local optimum for which the objective function value is less than the best found by the method in 1000 trials. Among 9 problems, frequency of suboptimal solutions was less than 0.5 % in 1000 trials. The general observation is that the likelihood for suboptimal solution is small. There are a few exceptions concerning set A. First, in the case e-ssd all methods most frequently find the second best value 14.680 (which lags 0.03 % from the optimum). Second, in cases e-lld and e-lsd two of the seven versions produce the second best value 12.537 in more than 10 % of runs; however, the optimum 12.544 is only 0.05 % higher. For the random initialization step in Section 3.2 no suboptimal solution was found in 45 out of 48 test cases. In two cases suboptimal frequency was less than 0.5 % and in the case e-ssd only the second best value 14.680 (which lags 0.03 % from the optimum) is found in 1000 trials.
Finally, we note that for any particular problem initialization by any of the seven methods does not necessarily lead to a local optimum. Due to non-convexities, the solver code may stop before optimality conditions are met. In such a case we try again and repeat until a successful completion. We observed for the seven versions that in 95 % of the case problems the average number of trials needed for success is close to one (the first trial being a success in most cases). For the initialization in Section 3.2 the average number of trials was more than 1.4 in three cases, the maximum average being 2.9 for n-sls.

B. MIP approach for fitness evaluation
In this section we consider mixed integer programing (M IP ) for the problem for which the evolutionary algorithm is proposed in Sections 3.1-3.2. Let t 0 and t 1 denote the beginning of the first and second cycle, respectively. A complicating feature for M IP is that both t 0 and t 1 are endogenous in the model. Consequently, for the steady state requirement x t 0 = x t 1 subscripts of decision vectors are endogenous. Furthermore, cash flows F t = R t −C t (refer Equation 9) are treated differently depending on whether t < t 0 , t 0 ≤ t < t 1 or t ≥ t 1 . We begin this section by reformulating the problem as a M IP problem and discuss the implementation of branch-and-bound algorithm thereafter. Computational results are presented in Section 4.4.1, including comparison with the evolutionary algorithm.
Let T = t max + s max be an exogenous upper limit for t 1 . For all t = 0, 1, . . . , T let δ t , u t and r t be binary variables such that δ t = 1 if harvesting takes place at stage t and δ t = 0 otherwise, u t = 1 if t = t 0 and u t = 0 otherwise, and r t = 1 if t = t 1 and r t = 0 otherwise. Hence, T t=0 u t = 1 and T t=0 r t = 1, and the cycle length is where s min and s max are exogenous bounds for cycle length, and t min and t max are exogenous bounds for transient period.
The level of harvest h t (number of trees at stage t by size class) is restricted by wherex is an exogenous upper bound for the forest state vectors x t (number of trees at stage t by size class). Harvests are only allowed during the period [0, t 1 ). Hence we require At the beginning of each cycle we require the state to be an endogenous vector x c . Hence for all t ≥ t min , we require For all t, the cash flow F t is split into