Maintaining an Optimal Flow of Forest Products under a Carbon Market : Approximating a Pareto Set of Optimal Silvicultural Regimes for Eucalyptus fastigata

A competitive co-evolutionary Multi-Objective Genetic Algorithm (cc-MOGA) was used to approximate a Pareto front of efficient silvicultural regimes for Eucalyptus fastigata. The three objectives to be maximised included, sawlog, pulpwood and carbon sequestration payment. Three carbon price scenarios (3CPS), i.e. NZ $25, NZ $50 and NZ $100 for a tonne of CO2 sequestered, were used to assess the impact on silvicultural regimes, against a fourth non-carbon Pareto set of efficient regimes (nonCPS), determined from a cc-MOGA with two objectives, i.e. competing sawlog and pulpwood productions. Carbon prices included in stand valuation were found to influence the silvicultural regimes by increasing the rotation length and lowering the final crop number before clearfell. However, there were no significant changes in the frequency, timing, and intensity of thinning operations amongst all the four Pareto sets of solutions. However, the 3CPS were not significantly different from each other, which meant that these silvicultural regimes were insensitive to the price of carbon. This was because maximising carbon sequestration was directly related to the biological growth rate. As such an optimal mix of frequency, intensity, and timing of thinning maintained maximum growth rate for as long as possible for any one rotation.


Introduction
The central focus of our analysis was to approximate a set of optimal silvicultural regimes for a Eucalyptus fastigata forest stand under a carbon market.Each estimated regime was expressed as a set of values that included an initial planting stocking, frequency of thinning, timing of thinning, intensity of thinning, final crop number prior to clear-felling, and rotation length.We, therefore, crafted a three-objective optimisation problem, which simultaneously maximised, sawlog, pulpwood and carbon sequestration payment (under three different payment scenarios).This optimisation problem, described later, was based on a two-objective optimisation problem that was successfully solved by simultaneously optimising competing sawlog and pulpwood products (Chikumbo & Nicholas, 2011).The results from the two-objective and three-objective optimisation runs were statistically analysed to decipher the nuances of silvicultural strategies under a carbon market.

Assumptions and Forest Holding Value
To carry out the analysis, we assumed fixed prices for liquid fuels and fossil fuel-based fertilisers such that the carbon price would remain static over the rotation period.This assumption was based on the observation of the EU Emissions Trading Scheme where the carbon price was heavily influenced by fossil fuel prices, which tend to be volatile (White, 2007).The carbon price would in turn influence the forest holding value, and ultimately impact the stand silvicultural regimes of E. fastigata.Forest holding values enable valuation of timber as real property, where timber for immediate harvesting has a liquidation value and the immature resource has a holding value (Mayo & Straka, 2007).
Any forest has an immediate liquidation value if the existing timber is clearfelled and sold along with the land.The forest holding value is the present value of holding the forest until the optimal rotation age (maximum present value) and then selling the timber and land (Klemperer, 1996).The concept is consistent with standard forestry valuation concepts such as land expectation value (Faustmann, 1995;De Jong, Tipper & Montoya-Gómez, 2000).Note that the forest is financially immature for as long as the forest holding value exceeds the forest liquiddation value.Therefore, the rotation age should be allowed to increase until the two values are equal (Mayo & Straka, 2007).Thus, the forest holding value provides an ideal financial criterion to evaluate the impact of carbon sequestration payments on the optimal rotation length.

Contribution to Forest Literature
There is an expectation that any forester/land owner wishing to engage in forest-carbon trading, in order to take advantage of a new income stream from carbon sequestration, would want to know the ideal/optimal silvicultural regimes for his/her crop that will not only maximise carbon sequestration (for maximum pay out), but also maximise production of sawlog and/or pulpwood (De Jong, Tipper & Montoya-Gómez, 2000).However, we do know that many forest analysts have shown that increasing the rotation length would be the sensible thing to do (Appel, 2001;Asante, Armstrong & Adamowicz, 2011;Gutrich & Hoswarth, 2007).What is scarce in literature is: 1) How the frequency, intensity and timing of thinning for a silvicultural regime are affected; 2) Whether the longer rotation length is linked to a higher or relatively lower final crop number before clearfell; 3) The ability to simultaneously cater for pulpwood and sawlog products under a prevailing carbon market from a single regime instead of different regimes where each specifically caters for a unique product; and 4) How the silvicultural regimes are affected by different carbon prices.
Our paper addresses these issues.The specifics of determineing the optimal initial planting stocking, optimal rotation length, frequency, timing and intensity of thinning, tree species and site, have been modeled by forest analysts using multi-stage optimisation since the 60s' (Hool, 1965), with mixed successes.The reasons for these mixed successes boiled down to the use of inappropriate growth functions, and an inability to do an exhaustive search for all possible states in a dynamic programming formulation (Chikumbo, 1996;Chen, Rose & Leary, 1980).Bellman (1957) coined this exhaustive search problem, the "curse of dimensionality".

Historical Background of Problem Solving
A pulpwood production was characterised by short rotations and a relatively longer rotation with thinning (i.e., partial harvesting) for a sawlog production regime (Newman & Wear, 1993).To overcome the curse of dimensionality, the determination of optimal silvilcultural regimes for separate pulpwood and sawlog production was pursued using a specialised mathematical formulation, i.e. a combined optimal control and optimal parameter selection optimisation (Chikumbo & Mareels, 2003).The growth dynamics of a tree crop stand were described with difference equations in discrete-time and the complete mathematical formulation was as follows: subject to the growth dynamics, where, J 0 = the cost functional (i.e. a function of state and control variables or simply the objective function); o  = continuously differentiable function; j 0 = continuously differentiable function with respect to the state and control variables; t = 0, 1, The state vari le consisted of the mean stand height ab , which as the mean height in metres per hectare (mht(t)), the stand basal area in square metres per hectare (sba(t)), and stand volume in cubic metres per hectare (vol(t)), i.e. where, , , The problem (1)-(10) was solved using Pontrya m gin's Maxium Principle (PMP) (Chikumbo & Mareels, 2003).Only a single objective problem was solved at any one time, i.e. either a value production or a volume production cost functional.It is possible to solve a multi-objective optimal control problem using PMP (Malinowska & Torres, 2007) for a finite number of cost functionals, but this was never meant to be because of one problem.Trying to estimate the optimal rotation length as a system parameter led to ill-conditioning (i.e. an accumulation

Research Focus and Problem Solving
orithm eliminated th ando t of efficient silvicultural regimes that satisfi of round-off errors from a sequence of matrix inversions that The critical part of approximating the set of efficient thinning regimes was to find "trade-off" solutions (i.e.non-dominated solutions) where for each solution an improvement in one objective did not lead to worsening in the other (Osborne & Rubenstein, 1994).The set of solutions to the three-objective problem was determined using a competitive co-evolutionary genetic algorithm with five subpopulations of 100 individuals each, computed over 1000 generations.These sub-populations evolved independently for a certain number of generations (isolation time).After the isolation time a number of individuals were distributed between the sub-populations (a process called migration).Each sub-population exerted selective pressure on the other, thereby maintaining diversity a lot longer than each sub-population would do solitarily, thereby guarding against premature convergence.When competition was superimposed between the sub-populations, the ones with higher mean fitness values were allowed to maintain larger sub-population sizes and received more capable individuals, since they had more chances of finding the global optimum (Chikumbo, 2009b;Chikumbo, 2012).
are numerically unstable and that may result in meaningless solutions).We will call this the phase-1 ill-conditioning.The stopgap measure was to solve a series of combined optimal control and a single optimal selection problem (i.e. the initial planting stocking), with a fixed T at increasing intervals, so as to locate the optimal rotation length.Optimal T was the point at which a unit increment in T resulted in "phase-2 ill-conditioning", caused by growth dynamics that could not perform outside their range.The whole process was time-consuming, giving a clear signal that a better way had to be found.
A switch to a single objective genetic alg e phase-1 ill-conditioning problem (Chikumbo, 2009a).Chikumbo and Nicholas ( 2011) demonstrated a two-objective genetic algorithm that simultaneously optimised for a value and volume production for Eucalyptus fastigata.Figure 1 shows a summary of a generic genetic algorithm and how it works.
A genetic algorithm is initialized with a population of r The Fonseca and Fleming ranking scheme (Fonseca & Fleming, 1993) was used to determine the non-dominated solutions, also referred to as the Pareto front.A conflict in the objectives results in a trade-off set (i.e.Pareto), which means that the solutions in the set are optimal in the wider sense that no other solutions in the search space are superior to them when all the objectives are considered.mly generated individuals which is a guided process of "selection", "crossover/recombination" and "mutation".Individuals are selected on the basis of their fitness for reproduction.The parent individuals are recombined to produce offspring where only some of them are mutated with a certain probability.The fitness of the offspring is then computed, resulting in the parents being replaced, thus producing a new generation.If the criteria of the objective function(s) are not met, this cycle is performed again until the optimisation criteria are reached (Polheim, 2006).
Therefore, a se Fonseca and Fleming called their ranking scheme, Multi-Objective Genetic Algorithm (MOGA) and it involves assigning an individual's rank (in the objective function space) equal to the number of population individuals that dominated that individual.What this means is that ranking of the individuals prior to selection for recombination is done according to the degree of domination; the more members of the current population that dominate a particular individual, the lower its rank.MOGA, therefore, uses fitness sharing in the objective function space and recombination is also restricted.Reproduction probabilities are determined by means of exponential ranking.Afterwards the fitness values are averaged and shared among individuals having identical ranks (Zitler, Deb, & Thiele, 2000).Finally, stochastic universal sampling, which provides zero bias ed carbon sequestration payments, sawlog/value production, and pulpwood/volume production at different levels of "tradeoffs" was estimated using genetic algorithms, eliminating the need for determining separate regimes for sawlog and pulpwood under a carbon market.Such regimes give the forester/ land owner the ability to satisfy the market with all products whilst maximizing carbon sequestration payments, with the flexibility of meeting an increased/decreased supply of any one of the products as dictated by demand.The Model There were two main c-MOGA was used to so orithms mainta Table 1.
f the 286 data points measured from the Nelder trial (Chi-

Minimum
Mean Maximum a (i.e. the range of possible v number of offspring of an individual), is used to fill the sampling pool.The main strength of MOGA is that it is efficient and relatively easy to implement.It has also been successfully implemented in solving optimal control problems with good overall performance (Coello, 1996).
In this paper we follow the form icholas, 2011), a competitive co-evolutionary Multi-Objective Genetic Algorithm (cc-MOGA), but with an additional objective, i.e. maximisation of carbon sequestration payment.We focus on the threeobjectives and the forest holding values (held static because of the assumption of fixed prices of liquid fuel and fossil fuel-based fertilisers) on how they are formulated and discuss the results of the three-objective cc-MOGA.The same species, E. fastigata, and growth functions by (Chikumbo & Nicholas, 2011) are used in our analysis.
& Maid ost suitable eucalypt for a wide range of sites in New Zealand because it has a wide site tolerance, performing well from near sea level to altitudes up to 500 m (Chikumbo & Nicholas, 2011).It can be utilised for solid timber as well as providing short-fibre pulp for the production of fine printing paper (Miller, Hay, & Ecroyd, 2000;Haslett, 1988).
The data used for the growth dynam om a 1979 Nelder trial (Nelder, 1962) established in Kaingaroa Forest (latitude 38˚27.6'S,longitude 176˚39.9'Eand an altitude of 280 m), near Murupara, New Zealand.A Nelder trial consists of trees planted in a series of concentric circles where the growing space available to each tree is determined by the distance to the nearest eight neighbouring trees.
Further detail on the trial is found in (Chikum 11) and the summary of the data is shown in Table 1.The costs and revenue for silvicultural treatments used in determining forest holding values for the valuation of timber products and carbon sequestered were obtained from (Turner, West, Dungey, Wakelin, Maclaren, Adams & Silcock, 2008).reasons why c lve the three-objective optimisation problem: 1) A competitive co-evolutionary genetic alg ins diversity and controls the selective pressure by balancing avoids premature convergence (Chikumbo, 2009b;Menczer, Degeratu, & Street, 2000); and 2) Splitting the population into diverse sub-populations that communicate through migration may result in parallel speedups (Menczer, Degeratu, & Street The three objectives or cost functionals were as follows: for volume production; and for carbon sequestration payments.subject to the constraints,   for the rotation length, where t 0 is t m with the three cost funcir (x, tionals, J n (u) , J n (u) , and J n (u) , there exists no solution pa u) that renders a global minimum value to each of the functionals simultaneously.Rather, there exists a finite set of solutions that represent trade-offs.A key concept in determining this set of solutions is that of Pareto optimality: Assume two solutions (x, u), (x′, y′)   , where  is denoted as the solution space.Then (x, u) is said to dominate (x′, y′) (also written (x, u)  (x′, y′) where, est at time,

PV
= Discounted costs at time t costs ; and LOC = Land opportunity cost.At an point in time carbon sequestered ange in biomass and the amount of carbo is not just the age of the trees per se or standing timber volume that is important, but rather the rate of tree growth (van Kooten, Binkley & Delcourt, 1995).As trees grow they sequester carbon, but once carbon has been sequestered no further benefits are forthcoming.In other words the income generated by sequestering carbon at time, t is the value of the extra amount of carbon sequestered between t and t -1 multiplied by the price of carbon, as shown in the following equation: where,

S t  
or sequestering carbon ($•m -3 ); d in $ rice f S(t) = Tota = Carbon sequestration payment at time, t; Pc = P l carbon sequestered at time, t, for a stand; an t -1).S(t -1) = Total carbon sequestered at time, ( he present value of the above equation at the T time of plantg the crop then becomes: where, (t) = Discounted carbon sequestration payment at time, t.

qu (25)
An as a

Results and Discussion
The results o problem with th ag

PV carb
Therefore, the forest holding value that includes carbon seestration payments then takes the following form: sumption was implicit in the calculation that there was zero harvesting emissions liability.All it meant was that there was no carbon sequestration payment for the proportion of harvested timber through thinning or clearfelling.
f the two-objective optimisation e objectives ( 11) and ( 12) for value and volume products are shown in Table 2 with only 6 regimes included in the Pareto set.The other Pareto sets for the three carbon prices scenarios, NZ $25, NZ $50, and NZ $100, are shown in Appendices A (Table A1), B (Table B1), and C (Table C1), because of the larger size of their sets, at 30, 30, and 29 regimes respectively.
All these results in Tables 1, A1, B1 and C1 were not avered, as this would defeat the purpose of finding a Pareto set of solutions in the first place.Therefore, investigating the impact of carbon prices on the silvicultural regimes of E. fastigata involved comparisons of the Pareto sets using the Kruskal-Wallis test (Kruskal & Wallis, 1952).This is a non-parametric version of the classical one-way analysis of variance and an extension of the Wilcoxon rank sum test (Gibbons 1985), specifically for testing equality of population medians among groups (which in our case are the four scenarios in Tables 1, A1, B1 and C1).
Ordering the data from smallest to largest across all groups and taking the numeric index of this ordering determines the ranks.f 0.4018 suggests that we accep thesis that the initial planting stockings from the different groups have the same mean.The box plots in Figure 2 confirm this with the overlap of the medians at one standard deviation of the group means.The Wilcoxon rank sum test, plotted in Figure 3, shows more clearly the degree of overlap and no statistical difference between the initial planting stockings of the different groups.We conclude here that different carbon prices do not seem to have an impact on the initial planting stocking of E. fastigata.T T able 2. hinning r tives, sawlog (value) and pulpwood (volume) productions (Chikumbo & Nicholas, 2011).

REGIME
Age at T1   inal Crop Number suggested that there were differences be difference in the rotation length medians am The non-carbon scenario showed a higher number of retained crops before final harvesting, than all the carbon scenarios, which showed an overlap of the one standard deviation of their final crop numbers.The mean ranking in Figure 5 showed more clearly the degree of overlap among the three carbon scenarios.

Rotation Length
There was also a ongst the different groups with a p-value of 0.003.The non-carbon scenario showed no variation in the rotation length, which remained at 35 years for all the individual regimes in the Pareto set, whereas the other three carbon scenarios had higher medians, which overlapped at one standard deviation of the the mean ranks of the rotation lengths in Figure 7 confirmed the differences between the non-carbon scenario with a lower rotation length, and the three carbon scenarios with higher ranges of rotation lengths.This observation agrees with conventional wisdom that rotation lengths will be longer in order to sequester more carbon.What is interesting with our results here is that the rotation lengths are insensitive to the variations in the carbon price.
We assert that this insensitivity to the price of carbon is because our optimisation model is maxing out the sequestration of carbon as much as the equations for the growth dynamics of E. fastigata will allow.This is good news for the forester, in that for any forest signed up to an emissions trading scheme, there is only one set of optimal regimes to consider that will simultaneously satisfy a sawlog and pulpwood market, regardless of fluctuations in the carbon price.

Frequency, Timing a sity of Thinning
The frequency, timi ios.Cargy where it carbon scenario, and that the initial planting stockings of all the y the same, it therefore stands to nd Inten ng and intensity of thinnings did not four scenarios were statisticall show any statistical differences for all the four scenar bon prices do not seem to influence the thinning strate is already optimized for value and volume productions.This might be explained by the fact that though stand volume is reduced through thinning, the sudden availability of more nutrients, light and moisture to the residual trees boosts their growth.It is this growth that will guarantee more sequestration and possibly more payment.Given that the final crop numbers of all the carbon scenarios were lower than those of the non-confirm this assertion.Also with less number of trees as a final crop, it is possible to keep the trees a little longer than one would normally do in a non-carbon market environment because this may guarantee more growth until full-site occupancy is reached.This might encourage fertilisation following a lateage thinning, in order to boost growth and subsequently sequester more carbon.Implications of late-age fertilisation following a thinning may also mean a premium sawlog/veneer product at the end of the rotation.1) A dynamic stumpage model; and 2) A fossil fuel-based fertiliser price model, under a carbon trading scheme, and revisit our analysis.Both (1 lding value to a m ) and (2) will expose the forest ho alistic output, and maybe provide us with fur e impact of carbon prices on the silvicultural regimes of E. stigata.
stands at different age classes.With each stand with a Pareto set of possible silvicultural regimes, it is possible to optimise at an estate level, assigning the appropriate regime to each stand, which may ents and optimally sequestering carbon under prevailing market constraints.In other words the Pareto optimality at a stand level gives flexibility at an estate level planning.We have not touched on environmental constraints because of the focus of this paper.However, environmental issues are best dealt with at a forest estate level both temporally and spatially, given their long-term gestation period.

Conclusion
The consequences on the silvicultural regimes for Eucalyptus fastigata, when the crop is simultaneously managed for, carbon sequestration, sawlog and pulpwood were, decreased final crop numbers, and increased rotation lengths.It will be well worth it to investigate fertilisation following late-age thinning in order to boost growth.This would mean sequestering more carbon and guarantee of a premium sawlog/veneer product at the end of the rotation.The reg rbon market were also fo

A
p-value o t the null hypo egimes (showing number of trees thinned at T1, T2 and T3) for E. fastigata derived from a two-objective optimisation problem with objec-

Figure 2 .
Figure 2.f the four scenarios where, the central mark in each box is the median; the Box plots o edges of the box are the 25th and 75th percentiles; the notches in the boxes represent one standard deviation of the mean; the ends of the whiskers are the minimum and maximum; and outliers are plotted individually with the "+" sign.

Figure 3 .
Figure 3. nk sum test for the initial planting stockings of the different groups/scenarios.

Figure 6 .
Figure 6.The summary of

Figure 4 .
Figure 4. Box plots of the four scenarios.

Figure 5 .
Figure 5. Wilcoxon rank sum test for the final crop numbers of the different groups/scenarios.

Figure 6 .
Figure 6.Box plots of the four scenarios.

Figure 7 .
Figure 7. Wilcoxon rank sum test for the rotation lengths of the different groups/scenarios.Fur ussionThese results may ha ere matched with forecasted the prices of liquid fuel and fos s this would have m