Integrating risk management tools for regional forest planning: an interactive multiobjective 5 value at risk approach

Abstract


Introduction
Regional level forest planning involves balancing the extraction of resources over time and over large spatial scales.How the resources are extracted depends on the preferences of the decision maker, such as the income needs or desired forest qualities (i.e.old growth forests, mushroom and berry provision) expected from the forest.The provision of ecosystem services and biodiversity conservation inevitably involves trade-offs (Triviño et al. 2017).For instance, the provision of biodiversity competes directly with the provision of timber resources.As the management of forest can be done in a wide variety of manners, forest planners and decision makers need to consider how we can best balance the resource needs with conservation considerations.
Selecting the optimal treatment schedule for each stand is a combinatorial optimization problem, where the size of the problem depends on the number of separate stands, the number of possible treatment schedules for each stand and the number of constraints included to the problem.In forest management, linear programming has been utilized since the 1960's (Kangas et al. 2015).The majority of forest applications include multiple objectives; however these applications have generally been solved using an ε-constraint method (see, e.g., Miettinen 1999) which converts all objectives but one into constraints.As a means to compensate for computational complexity (especially with spatially specific problems), heuristic optimization has also been used.Within a heuristic framework, multiobjective optimization has typically been based on a weighted additive utility function (Kangas et al. 2015).The main reason for not utilizing multiobjective optimization more has been the lacking integration of computational tools into the forest planning systems, a lack of knowledge of the potential benefits these tools could provide and the various conceptual challenges for utilizing these tools.
In multiobjective optimization, rather than a single optimal solution, there are multiple so-called Pareto optimal solutions which form the Pareto optimal frontier (one example subset of this frontier is the Production Possibility Frontier).This frontier of solutions represents the case when none of the objectives can be improved without impairing at least one of the others.To determine which solution is preferable to the others, additional preference information is required from a decision maker (DM), to decide which of the Pareto optimal solutions is the most preferred one.Different methods incorporate DM's preferences before (a priori) or after (a posteriori) having generated Pareto optimal solutions (see e.g.Miettinen 1999).In interactive multiobjective optimization methods (e.g.Steuer 1986, Miettinen 1999, Miettinen et al. 2008), solutions are iteratively improved towards a solution satisfying the DM by alternating the steps of re-defining the preferences and producing a plan fulfilling these preferences.
The changing preferences can be expressed, for instance, by giving a satisfactory value or desirable target value, or minimum or maximum bounds for objective values or weights for the criteria (e.g.Liu and Davis 1995, Krcmar-Nazic et al. 1998, Kazana et al. 2003).
Management of forest resources requires the use of models to predict the future availability of forest resources and the expected state of the forest through time.Decisions regarding the management actions can be adjusted, and this will alter the provision of forest resources, and will influence the development of the forest.The ability to predict the future state of forest resources depends upon growth models and the quality of the initial forest information used to run these models.If poor models are used to predict the development of the forest, and if poor quality information is used to run the models, the anticipated uncertainty of the predictions can be expected to be larger than if higher quality data and models are used (Diebold 2001).
A vast majority of the applications within forest management planning assume a decision situation under certainty, even though incorporating elements of uncertainty can be accomplished through a wide variety of methods (Pasalodos-Tato et al. 2013).In forestry specific applications, research with an aim to mitigate the impact of natural hazards (e.g.fire, insects and wind throw events) has been done through stochastic programming (Gassmann 1989, Boychuk andMartel 1996) and chance-constrained programming (Bevers 2007).Robust programming (Bertsimas and Sim 2004) has been applied to develop an operational forest management plan which maximizes net present value, while ensuring that harvest constraints are met with a high level of certainty (Palma and Nelson 2009).In a multi-criteria context, stochastic goal programming has been applied to a small forest holding (Eyvindson and Kangas 2014).Additionally, through Markov Chain models, Buongournio and Zhou (2017) have managed multi-criteria risk neutral decision problems while Buongournio et al. (2017) have incorporated risk preference into the decision making problem.The study highlights the potential for improvement to the plan simply by incorporating the estimates of uncertainty which are available.In this study our use of the terms uncertainty and risk differs slightly from the Knightian definitions (Knight 1921) (where risk has a known probability, and uncertainty has an unknown probability).We use uncertainty to reflect attributes of forest data (where the exact probability is not known, however can be estimated), and where risk refers to perceptions these uncertainties have on the DM.
The measurement of risk is an important aspect which has been thoroughly applied in the realms of both engineering and business.The ISO has defined risk as the "effect of uncertainty on objectives" (ISO Guide 73:2009).This is an open definition, allowing risk to be both positive and negative.Often DMs are more interested in the negative side of risk (in a financial context this could be interpreted as a loss).From this perspective, there are several downside risk measures which are frequently used, such as the mean semi-deviation (Krzemienowski and Ogryczak 2005), the Value-at-Risk (VaR; Duffie and Pan 1997) and Conditional Value-at-Risk (CVaR;Rockafellar and Uryasev 2000).
Forest applications utilizing risk measures are increasingly often found in the forest science literature (see Eyvindson and Kangas 2016a, Eyvindson and Cheng 2016, Piazza and Pagnoncelli 2015, Roessiger et al. 2011).To our knowledge, applications including two or more different value at risk concepts have not been published in any field.In forest planning, it is quite possible that the DM is willing to accept high risk for timber cash flow but only a low risk for losing biodiversity or vice versa (this can be understood as an asymmetrical aversion to gains and losses (Kahneman and Tversky 1979)).The main reason why the uncertainties have been ignored so far is that the problem with including uncertainty requires quite heavy calculations, which may not be accomplished in a reasonable timeframe for large areas with a lot of stands.Additionally, there is sometimes a challenge to objectively (or even subjectively) quantify the probability distributions associated with the uncertainties.
The risk measure used in this application is the VaR.The VaR measurement is well known, and very commonly used in finance and business to identify the unit's exposure to risk.Quite simply, the VaR identifies a threshold for losses, so that losses larger than this threshold can only occur with a given probability.Mathematically it can be defined as: where X is the random loss variable, and ߙ is the critical value for the confidence interval (commonly 0.9, 0.95 or 0.99).This risk measure can be criticized (Krokhmal et al. 2011): it does not utilize information regarding the tail of the distribution (so that very large losses can occur), and that the ܸܴܽ ఈ ሺܺሻ is a nonconvex function of X, meaning that the VaR of the total may exceed the VaR of the sum of the parts.However, according to Durbach and Stewart (2012), while these critiques are valid, they may not be as important in a multi-criteria context, and the value of having the DM being able to understand the risk measure may supersede the value obtained by using more complicated risk measures.
The aim of this paper is to integrate the VaR as a specific objective function into a regional forest management planning problem.As mentioned earlier, in a multi-criteria context, a DM is required to provide some form of preference information in order to get more preferred solutions.We extend the research of Hartikainen et al. (2016) with more objectives and we employ an interactive multiobjective optimization method (a subset of NIMBUS, see Miettinen and Mäkelä 2006) where the DM provides preference information and then new alternative forest management plans are generated using the achievement scalarizing function approach of Wierzbicki (1982).To illustrate the impact of integrating VaR as an objective function into the decision making process, a large regional forest holding in central Finland is used as an example.

Methods
Multiobjective optimization problems can in general be formulated as maxሼ݂ ଵ ሺ‫ݔ‬ሻ, … , ݂ ሺ‫ݔ‬ሻሽ Page 7 of 30 Can.J. For.Res.Downloaded from www.nrcresearchpress.com by JYVASKYLAN YLIOPISTO on 04/10/18 For personal use only.This Just-IN manuscript is the accepted manuscript prior to copy editing and page composition.It may differ from the final official version of record.
In the above problem, ݂ are objective functions to be maximized, ݃ and ℎ are inequality and equality constraints, respectively, bounds ܽ ≤ ‫ݔ‬ ≤ ܾ are called box constraints and, finally, the decision variable vector ‫ݔ‬ = ൫‫ݔ‬ , ‫ݔ‬ ൯ consists of integer valued variables ‫ݔ‬ and real valued variables ‫ݔ‬ .We denote the feasible region determined by all the constraints by ܺ ⊂ ܴ and objective vectors by ‫ݖ‬ = ݂ሺ‫ݔ‬ሻ = ൫݂ ଵ ሺ‫ݔ‬ሻ, … , ݂ ሺ‫ݔ‬ሻ൯ ∈ ܴ .
As mentioned earlier, we focus on applying interactive methods since they allow the DM to learn about both the interdependencies involved among the objectives and the feasibility of one's preferences.
Eliciting preferences from DMs can be done through a wide variety of methods (see e.g.Luque et al, 2011, Miettinen et al. 2016, Steuer 1986).One of these approaches is to elicit the preferences as aspiration levels, the value at which the DM desires the objectives should achieve, and together the aspiration levels constitute a reference point ‫ݖ‬ ∈ ܴ .Aspiration levels are employed in so-called achievement scalarizing functions (Wierzbicki 1982(Wierzbicki , 1986)).In a theoretically justifiable manner, achievement scalarizing functions measure the preferability of a solution provided a specific reference point z ref .Wierzbicki (1986) defined this as order-consistency.
For this paper, we use the following achievement scalarizing function (asf) to be maximized objectives within the set of Pareto optimal solutions (this vector is typically estimated, see e.g.Miettinen 1999).We utilize the nadir point rather than an anti-ideal point, which contains the minimum values of the individual objectives overall, in order to focus on Pareto optimal solutions.The summation term at the end is called an augmentation term and the constant ߩ is a small positive constant, e.g.0.0001.The augmentation term guarantees that the solutions are indeed Pareto optimal (and the trade-offs among the objectives are bounded by ߩ and 1/ ߩ).The difference of the components of the ideal and nadir points normalizes objective functions to have similar magnitudes.
This leads us to solving the optimization problem max ௫∈ ‫ݏ‬ ௦ ሺ݂ሺ‫ݔ‬ሻ, ‫ݖ‬ ሻ for reference points ‫ݖ‬ given by the DM.Given a reference point, the optimal solution to the above problem is a Pareto optimal solution to the original multiobjective optimization problem (Wierzbicki, 1982(Wierzbicki, , 1986)).We apply a variant of the reference point method of Wierzbicki where the DM iterates as long as (s)he wishes by providing a reference point and seeing the corresponding Pareto optimal solution.In this way, (s)he learns about the attainability of the desired objective function values and the conflicting nature of the objectives.The method was chosen because the preference information and the solutions obtained are all objective function values and there is no need of any cognitive mapping in the mind of the DM (as e.g. when applying weighting coefficients).
We assume that we have S stands, T treatment schedules for each stand, R scenarios and P periods to consider.We have simulated values ‫ܫ‬ ௧,௦,, for the timber cash flow and ‫ܤ‬ ௧,௦,, for the biodiversity with indices denoting treatment schedules t (including one or more timed treatments), stands s, scenarios of the future development of the forest stand r and 5-year periods p.The problem of choosing the best treatment options for each of the stands can be formulated as a multiobjective optimization problem In the above problem, we have six objectives to be maximized: 1. Minimum timber cash flow in the scenarios that are in the set of scenarios I R in euros.
According to the first constraint, the set I R is a subset of the complete set of scenarios R and the number of scenarios in this subset is greater or equal to I δ times the number of scenarios in the complete set of scenarios.This means that this is the timber cash flow at risk for the risk level1 I δ − .The objective is denoted by I VaR .
2. Minimum change between subsequent periods for the biodiversity index in the scenarios that are in the set of scenarios B R .As an index, this variable is unitless, and relative differences can be more informative to decision makers.According to the first constraint, the set B R is a subset of the complete set of scenarios R and the number of scenarios in this subset is greater or equal to B δ times the number of scenarios in the complete set of scenarios.This means that this is the biodiversity at risk for the risk level1 B δ − .The objective is denoted by B VaR .
Page 10 of 30 Can.J. For.Res.Downloaded from www.nrcresearchpress.com by JYVASKYLAN YLIOPISTO on 04/10/18 For personal use only.This Just-IN manuscript is the accepted manuscript prior to copy editing and page composition.It may differ from the final official version of record.
3. Expected minimum timber cash flow for all of the periods p in the complete set of scenarios R in euros.This objective will be denoted by I E .
4. Expected minimum change between subsequent periods of the biodiversity index for all of the periods in the complete set of scenarios R .This objective will be denoted by B E .
5. Probability I δ of the set I R .The risk level for the value at risk for the timber cash flow is, thus, 1 I δ − .The objective is denoted by I δ .
6. Probability B δ of the set B R .The risk level for the value at risk for the biodiversity is, thus, 1 B δ − .The objective is denoted as B δ .
At regional levels, DMs may desire a relatively smooth flow of timber resources.This can be interpreted as an even-flow solution (see, e.g.Bettinger et al. 2009, Kangas et al. 2015).In an evenflow solution, the timber cash flow is equal over all periods.From an ecological perspective, a similar importance can be given on biodiversity stability.For this case, rather than ensuring an even-flow, we promote the idea of maximizing the minimum timber cash flow and biodiversity increase over periods.
The solution that maximizes the minimum timber cash flow over the periods could sometimes equate with the even-flow solution, if all the periodic timber cash flows are the same as the minimum timber cash flow.Alternatively, it could be the case that this solution is better than the best available evenflow solution, if the timber cash flow in one of the periods is higher than the minimum timber cash flow.With respect to biodiversity, a minimum of a no change between all periods would provide a stable biodiversity.As biodiversity is generally considered to be threatened in Finnish forests, a trend of improvements is likely to be desirable.In this way, we believe that it is more justified to maximize the minimum over the periods than seek for the even-flow solution.
Page 11 of 30 Can.J. For.Res.Downloaded from www.nrcresearchpress.com by JYVASKYLAN YLIOPISTO on 04/10/18 For personal use only.This Just-IN manuscript is the accepted manuscript prior to copy editing and page composition.It may differ from the final official version of record.
In the problem formulation, the decision variables are the sets ܴ ூ and ܴ , and the treatment decisions ‫ݔ‬ ௧,௦ for all treatment schedules t and stands s.The treatment decision ‫ݔ‬ ௧,௦ is a binary variable, assigned one, if the treatment schedule t is chosen for the stand s and 0 otherwise.As only one treatment schedule can be chosen for each stand, the third constraint requires for each stand that only one of the ‫ݔ‬ ௧,௦ values can be one, while the remaining treatment decisions must be zero. In where: ൱.
In the above problem, the constant M is a big number (here 1,000,000) that allows for the minimum over the scenarios R in both the biodiversity increase and timber cash flow to be the minimum over the scenarios, for which the variable ‫ݐ‬ ூ or ‫ݐ‬ has the value one.Because of the two first new constraints, the new variables ‫ݐ‬ ூ and ‫ݐ‬ must be one for the ratio of scenarios given by their respective probability variables ߜ ூ and ߜ .
As a mixed integer linear problem, this can be extremely computationally expensive to solve with a high number of stands.To ease this computational burden, we replace the problem above with a problem where the decision variables for the treatment schedules ‫ݔ‬ ௧,௦ are allowed to take any real value between 0 and 1, instead of being binary variables.This is a common approximation, when optimization is used in forest management.The interpretation of a treatment schedule with non-binary values is that a part of the stand is treated to a different schedule than the rest of the stand.For this reason, we can do this approximation and still get solutions which can be implemented as treatment schedules for the stands.This leaves us with a multiobjective optimization problem, where there are both integer ‫ݐ(‬ ூ and ‫ݐ‬ ) and real-valued ‫ݔ(‬ ௧,௦ ) variables.

Materials
A large region of forested land north of Jyväskylä, Finland was used to demonstrate the application of this method.The region is composed of nearly equal compositions of Scots pine (Pinus sylvestris L.) and Norway spruce (Picea abies (L.) Karst.), and only a minor component of deciduous trees (mainly Birch, Betula pendula and Betula pubescens).The entire forest area has been segmented into a set of 1,080 stands (18,815 ha) where each stand represents a relatively homogenous forest area.The segmentation is based on a relatively cloud free Landsat image (LC81880162013155LGN00), utilizing an algorithm developed by the Technical Research Centre of Finland (VTT).
For this study, the stochastic features under consideration were the initial inventory errors and the growth model errors.For this case, we have excluded the stochasticity of timber prices through time, simply to ease the understandability of the results, inclusion of additional uncertainty would also require an increase in the number of scenarios used in the calculations.A Monte Carlo simulation was used to generate a set of 25 scenarios for each stand, resulting in a problem with a size similar to a 27,000 stand problem.The selection of 25 scenarios was based on previous research (Eyvindson and Kangas 2016b), and the problem solvability.Each scenario represents an equally probable initial inventory for the forest holding (i.e.‫‬ = ଵ ଶହ , ∀݅ ∈ ‫.)ܫ‬These data are estimates of the expected value of the forest resources, and the actual forest resources may differ from the expectation.The sources of uncertainty included into this analysis were initial inventory errors of basal area, height and forest stand age.We assumed that the inventory errors for the height, the basal area and the forest stand age were normally distributed, with a mean of 0 and a standard deviation equal to 10% of the mean of variables considered.These estimates of errors reflect the current state-of-the-art inventory methods for Finnish forest conditions (Mäkinen et al. 2010).Growth model errors were modelled using an AR(1) process, using the same approach as Pietilä et al. (2010).
The simulation of future forest resources was conducted using an open-source forest simulation and optimization software (SIMO, Rasinmäki et al. 2009).The forest simulator creates a wide range of possible treatment options through the use of a decision tree (Siitonen 1993).Depending on the initial conditions of the stand, a variety of management alternatives can be created by conducting optional silvicultural treatments, delaying or lengthening harvesting actions.Silvicultural activities include planting, tending and fertilizing while harvesting activities include final felling or thinning.For this case, the maximum number of possible sets of treatment schedules for all stands was 34.
In this study, in order to estimate the biodiversity of the forest, a weighted combination of habitat suitability indices of red-listed saproxylic species for boreal forests was used (Tikkanen et al. 2007).
These models require an estimate of dead-wood in the forest, which is especially difficult to measure using remote sensing techniques.The amount of dead-wood was estimated based on average quantities of dead-wood in Finland (Tomppo et al. 1999) and the age of the forest stand.Two functions were used to estimate the quantity of deadwood, one to represent the increase of dead-wood from a middle aged forest to an old-aged forest, while the other represented the decrease of dead-wood from a recently harvest stand to a middle aged one.While not applied in this study, an alternative method to evaluating deadwood would be to monitor simulated mortality and provide assumptions on how natural mortality, harvest residues and the decay rates may be reflected in the quantity of dead-wood at the stand.

Case study
The decision support environment has been done using a Jupyter Notebook, which has been made freely available at https://github.com/eyvindson/VaR.The data used in decision making is also available at the same repository.The data of the stands is an extraction from the forest simulation database, and represents those objectives of interest.For this case, the DM was a forest management planning expert.The experiments were run on an eight core Intel Core i7 4910MQ processor with 32 GB of RAM.The computer was running Windows 7.
The multiobjective optimization problem of forest management was modelled using Pyomo.The entire process is contained within a Python framework, calling an outside solver to solve the achievement scalarizing problem.While we used a commercial solver for this application, open-source options can be used very easily through the Pyomo framework.
Prior to introducing the DM to the decision problem, the components of the estimated ideal and nadir points were calculated.To aid in the evaluation of the results, all values were provided as per hectare values.These values were estimated using the pay-off table approach (Miettinen 1999), where one first optimizes each objective individually and forms the k-dimensional ideal point from the best values of each objective and estimates the components of the nadir point from the worst objective function values obtained at the decision variable vectors found.In the pay-off table ( At each iteration of the iterative method, the reference point specified by the DM was projected to be Pareto optimal by solving the achievement scalarizing function introduced earlier.As a way to test the proof of concept with a case study, we used a forestry expert as a DM.
To be more specific, the interactive decision process can be broken down into a few steps.First, the DM was shown the potential opportunities for the objectives under consideration.In this case, we provided the DM with the ideal and nadir points to get an understanding of the range of the objectives.
Once the DM had a frame of reference, we then asked the DM to provide a reference point which reflects the desired values of the different objectives.With the specific reference point, a single solution was found and the corresponding objective vector shown to the DM.From this, an iterative process continued, with the DM providing updated reference points until she was satisfied with the solution or did not wish to continue the process.If needed, a series of Pareto optimal solutions can be generated between any two solutions (as in the NIMBUS methods, see Miettinen and Mäkelä 2006).
This could be useful if the DM sees promising values for different objectives in different solutions and is interested in seeing how objective function values change between them.
The solution process was the following: After reviewing the ideal and nadir points, the DM was asked to provide the first reference point.With the first reference point, the DM was seeking a solution with moderate values for both timber cash flow and biodiversity change, and a slightly higher importance placed on the risk associated with timber cash flow as compared to the change in biodiversity: 200, 0.0015, 200, 0.0015, 0.95, 0.8).
It turned out that the DM was rather pessimistic since the Pareto optimal solution corresponding to this reference point had better values for nearly all objectives, with the single exception of the δ parameter on the VaR of the biodiversity change.The solution was: S1 = (276, 0.0021, 263, 0.0021, 1, 0.8), which the DM thought was rather satisfactory.However, she was interested in considering an improvement in the positive change in biodiversity.To see such an alternative, the DM was willing to sacrifice in some of the objectives and a new reference point was given: z ref2 = (250, 0.0022, 250, 0.0021, 0.95, 0.95).
It turned out that a Pareto optimal solution fully corresponding to the desires in the reference point was not available and the closest one was: S2 = (240, 0.00216, 236, 0.00216, 0.92, 0.92).
Thus, improving the last objective implied sacrifices in several others.After seeing this solution, the DM thought that the trade-off between timber cash flow and biodiversity was rather high, as the deterioration in the timber cash flow objectives did not seem to result in a large enough increase in the change of biodiversity values.The DM then provided a new reference point, which focused on the other side of the trade-off: z ref3 = (350, 0.0, 350, 0.0, 0.95, 0.95).
Page 18 of 30 Can.J. For.Res.Downloaded from www.nrcresearchpress.com by JYVASKYLAN YLIOPISTO on 04/10/18 For personal use only.This Just-IN manuscript is the accepted manuscript prior to copy editing and page composition.It may differ from the final official version of record.
After this solution, the DM highlighted the general acceptability of the solution produced from this solution (S3) and the first solution (S1).
While the DM indicated a general satisfaction with these solutions, to clarify the trade-off between these solutions, a set of 12 additional Pareto optimal solutions was generated between the acceptable solutions.This was done to explicitly demonstrate the range of options available between the acceptable solutions, and to perhaps find a solution more preferred than either of the two acceptable solutions.The set of reference points used for this analysis were vectors calculated between the two interesting Pareto optimal solutions, indicated as ݂ሺ‫ݔ‬ ଵ ሻ and ݂ሺ‫ݔ‬ ଶ ሻ.The reference points were generated as ݂ሺ‫ݔ‬ ଶ + ‫ݐ‬ሺ݅ሻ݀ሻ where d is a vector of the differences between the solutions (݀ = ‫ݔ‬ ଵ − ‫ݔ‬ ଶ ), and i =1,…, P where P is the number of intermediate solutions and where ‫ݐ‬ሺ݅ሻ = ݅/ሺܲ + 1ሻ.Once each of the 12 reference points were solved, these solutions were paired down to 3 of the most different solutions.These solutions and the two satisfactory solutions generated during the interactive process are illustrated as petal diagrams in Figure 1.

FIGURE 1
From the resulting set of five Pareto optimal solutions, the DM was able to select what she considered to be the most preferred solution M3 (from Figure 1): M3 = (280, 0.00208, 272, 0.00207, 0.88, 0.8).
The interpretation of this final preferred solution is that • the expected per hectare values of the minimum timber cash flow and minimum biodiversity change (i.e., minima over the periods) are 280€/ha and 0.00208; • there is only a 12 percent risk that the minimum timber cash flow is smaller than 272€/ha; and • the minimum biodiversity is guaranteed to be over 0.00207 with a 80% probability.
This was very much what the DM was hoping for, because the expected values are at satisfactory levels and the risk levels in the value at risk are acceptable (i.e., 12% and 20%).

Discussion
Through an interactive process, an approach to handling conflicts between ensuring a minimum level of timber cash flow and ensuring a positive increase in biodiversity change was examined.An interactive multiobjective optimization process allowed the DM to examine a range of solution options.
This method incorporated a scenario-based simulation approach to incorporate uncertainty.The decision problem was converted to a six-objective optimization problem, utilizing a relatively small set of scenarios (25) to represent the uncertainty.The final problem utilized an interactive reference point based method to find preferable treatment schedules for the forests.
The DM considered the interactive method of discovering a preferred solution to be essential.One aspect the DM appreciated was the ability to understand the trade-offs between the various objectives of the problem.Thanks to the interactive nature of the method, the DM was able to learn about the interdependencies between the objectives and the attainability of her preferences.The DM felt that this process increased her perceived confidence in selecting the final solution.Additionally, as in this case the biodiversity change objective is a unitless index, the DM found it difficult to assign preferences directly.The approach the DM used to formalize her preferences was through an interpretation of the results in a relative manner.This was done by relating the obtained values of the most recent solution to the previous solutions, and through this manner the DM was able to construct a set of preferences for the unitless index.
When reflecting on the interactive process, we noticed that some of the components of the reference points were not within the range of the nadir / ideal values for the specific objective.This highlights the importance of using a graphical interface to help guide DMs to provide reflective reference points.
Although the reference points provided by the DM were not entirely rational, the method still guided the DM to a satisfactory solution.As our focus was on a proof of concept for the method, we did not include visualization techniques.However our future research focus includes consideration of appropriate graphical user interfaces.
In this application, both the risk level (i.e. the δ parameter) and the value-at-risk with the specific risk level were used as objectives for both timber cash flow and biodiversity change.This highlighted the need for an interactive approach.As the value-at-risk depends strongly on the risk level selected, the conceptual interpretation between the objectives was clarified through the interactive process.Due to using adjustable risk levels, the preferences for these objectives could be difficult to provide in an apriori fashion.While the selection of a risk level is often set prior to the conduct of the analysis, this limits the flexibility of such approaches.
For this application, computational processing times were sufficiently quick (within one to two minutes) to allow for an active interactive process.For larger problems (either through the use of a larger holding, increased stand resolution, increased number of available treatment options, longer planning horizon or increased number of scenarios for approximating the uncertainty), this kind of a problem may become computationally taxing.For interactive methods this can cause issues as the DM may not feel comfortable waiting for the computation of new solutions (Korhonen and Wallenius 1996).For these more computationally expensive tasks, alternative methods of evaluating solutions to the problems may be needed.One option could be through the use of hierarchical planning.For instance, Pantuso et al. (2015) have developed a method of solving hierarchical stochastic programs for a maritime fleet renewal problem.They were able to solve larger problems in a hierarchical framework, which could not be solved using CPLEX alone.Alternatively, one could employ surrogate based methods (see, e.g.Tabatabaei et al. 2015).
Both the timber cash flow and biodiversity change objectives were analysed from a specific perspective of sustainable management.Both objectives included the assumption that the minimum periodic values should be maximized, resulting in a form of an even-flow of timber from the timber cash flow objective and a constant improvement in biodiversity as the planning period progresses.From a purely economic perspective, an alternative for sustainable management of income could be the maximization of net present value, as it provides the highest sustainable income for specific discount rates.For this case, the components of the nadir point for both objectives were non-negative, highlighting the potential for biodiversity improvement while maximizing the sustainable harvest of timber.These results may be partially due to how uncertainty has been modelled in the biodiversity index.However, the requirement for an even-flow of timber implies a harvest level lower in some periods and higher in others than when maximizing purely economic objectives (i.e.net present value).
As the biodiversity change objective used here is an index, i.e. a proxy to the real biodiversity, the assessment of the uncertainty is especially difficult.In this study, we essentially assessed the uncertainty of the proxy based on the uncertainty in the forest variables.However, it is much more difficult to include the uncertainty of biodiversity, i.e. to assess how well the proxy actually describes the true biodiversity.If the index is based on expert judgment, the variation among the experts used may be utilized.If primary data is available for predicting the biodiversity as a function of forest characteristics, uncertainty of the predictive models can be utilized.However, good quality data on biodiversity may be hard to find.
It is thus inevitable that in the studied case the uncertainty involved is an underestimate for the biodiversity, while it may describe the uncertainty in the proxy very well.This aspect obviously requires some future research.For instance, in addition to the 25 scenarios for forest characteristics, additional scenarios could be included to describe the different experts.Especially problematic is that an expert judgment may be biased, i.e. systematically overestimating or underestimating an effect of any given forest characteristics, which is difficult to account for in the stochastic analysis.The same concerns also any non-random or fuzzy uncertainty.With sensitivity analysis it may be analysed how much this kind of uncertainty would affect the results, and see if it is possible to account for such uncertainties by adjusting the VaR and δ preferences.

Conclusions
In this paper, we have presented a model for explicitly accounting for the tradeoffs between timber cash flows and biodiversity, and, in addition, the risks related to both of these objectives.The benefit of the used interactive reference point based method utilizing an achievement scalarizing function with an augmentation term, enables using desirable objective function values (rather than more abstract weighting) in evaluating the solutions, but ensures a Pareto optimal solution rather than inside optimum (often met when using goal programming).Moreover, the VaR makes a relatively easy way to assess the trade-offs between the expected values and the risk associated.It is important, that the DMs understand the risks in the decision they are about to carry out, and can evaluate the consequences.The probability (δ) with which at least a given minimum target value (VaR) is obtained should be comprehensible also to non-expert DMs.However, as the risk is a difficult concept, an interactive approach is recommended to improve the DMs abilities to give his/her preferences.This model is applicable also to other decision problems and other sources of uncertainty, except for non-random bias-type uncertainty.
Page 24 of 30 Can.J. For.Res.Downloaded from www.nrcresearchpress.com by JYVASKYLAN YLIOPISTO on 04/10/18 For personal use only.This Just-IN manuscript is the accepted manuscript prior to copy editing and page composition.It may differ from the final official version of record.All objectives are normalized between 0 and 1 using the ideal and nadir points.
Page 30 of 30 Can.J. For.Res.Downloaded from www.nrcresearchpress.com by JYVASKYLAN YLIOPISTO on 04/10/18 For personal use only.This Just-IN manuscript is the accepted manuscript prior to copy editing and page composition.It may differ from the final official version of record.

ୀଵ
In the above problem, ‫ݖ‬ ௗ is the ideal vector of the problem containing the maximum values of the individual objectives and ‫ݖ‬ ௗ is the nadir vector containing the minimum values of the individual Page 8 of 30 Can.J. For.Res.Downloaded from www.nrcresearchpress.com by JYVASKYLAN YLIOPISTO on 04/10/18 For personal use only.This Just-IN manuscript is the accepted manuscript prior to copy editing and page composition.It may differ from the final official version of record.

Page 23 of 30
Can. J. For.Res.Downloaded from www.nrcresearchpress.com by JYVASKYLAN YLIOPISTO on 04/10/18For personal use only.This Just-IN manuscript is the accepted manuscript prior to copy editing and page composition.It may differ from the final official version of record.

Figure 1 .
Figure 1.Petal diagrams of the two most promising solutions generated through the interactive process order to solve the above problem efficiently, it can be converted into a mixed integer linear optimization problem (MILP).Once we can formulate the problem as a MILP, we can use the efficient MILP solvers available e.g., IBM ILOG CPLEX Optimization Studio (see e.g., http://www-01.ibm.com/support/knowledgecenter/SSSA5P_12.6.3/ilog.odms.studio.help/Optimization_Studio/topics/COS_home.html)or Gurobi Optimization (see e.g., http://www.gurobi.com/documentation/).The above problem can be represented as the following multiobjective MILP max ൝‫ܫ‬ ,௦௨௦௧ , ‫ܤ‬ ,௦௨௦௧ , ‫ݔ‬ ௧,௦ min ∈ Page 12 of 30 Can.J. For.Res.Downloaded from www.nrcresearchpress.com by JYVASKYLAN YLIOPISTO on 04/10/18For personal use only.This Just-IN manuscript is the accepted manuscript prior to copy editing and page composition.It may differ from the final official version of record. ൱, Can. J. For.Res.Downloaded from www.nrcresearchpress.com by JYVASKYLAN YLIOPISTO on 04/10/18For personal use only.This Just-IN manuscript is the accepted manuscript prior to copy editing and page composition.It may differ from the final official version of record.

Table 1
), the rows represent objective values calculated at the decision variable vectors where an objective function obtained the best value (ideal on the diagonal) and the components of the nadir vector are the worst values of each column.The ideal values are bolded and the nadir values are underlined.The ideal point is thus (E I ,E B ,VaR I ,VaR B ,δ I , δ B ) = (

Table 1 .
Pay-off table for the multiobjective optimization problem per hectare.Ideal values are bolded and the nadir values are underlined.Can.J. For.Res.Downloaded from www.nrcresearchpress.com by JYVASKYLAN YLIOPISTO on 04/10/18For personal use only.This Just-IN manuscript is the accepted manuscript prior to copy editing and page composition.It may differ from the final official version of record.
Can. J. For.Res.Downloaded from www.nrcresearchpress.com by JYVASKYLAN YLIOPISTO on 04/10/18For personal use only.This Just-IN manuscript is the accepted manuscript prior to copy editing and page composition.It may differ from the final official version of record.
Can. J. For.Res.Downloaded from www.nrcresearchpress.com by JYVASKYLAN YLIOPISTO on 04/10/18For personal use only.This Just-IN manuscript is the accepted manuscript prior to copy editing and page composition.It may differ from the final official version of record.