ROBUST OPTIMIZATION FOR FOREST RESOURCES DECISION-MAKING UNDER UNCERTAINTY by CRISTIAN DERECK PALMA Forest Engineer, Catholic University of Chile, 1999 M.Sc. University of Chile, 2004 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES (Forestry) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) March 2010 © Cristian Dereck Palma, 2010 Abstract There is a general consensus that management decisions concerning forest resources are made in an intrinsically uncertain environment. However, decision-making tools used in forest management assume perfect information, leaving decision-makers to explore the most likely scenarios of uncertainty and determine the most reasonable management alternative. Although techniques that explicitly consider uncertainty exist, they increase the complexity of the models hence precluding their application to large-size problems. This dissertation describes the application of robust optimization concepts that explicitly consider uncertainty in forest management problems while keeping the models computationally tractable. By introducing some simplifying assumptions about uncertainty distributions, i.e. independency and uniformity, this approach allows for including uncertainty in many coefficients of the model. The methodology modifies the constraints for which feasibility is desirable and incorporates uncertainty in the technical coefficients by introducing an additional term. This term is an optimization problem in itself that introduces new constraints into the original model and acts as a buffer that guarantees constraint satisfaction for different uncertainty realizations. By changing the value of a robustness parameter, the trade-off between cost and robustness can be analyzed. The performance of this approach is explored through three structurally different problems: (a) a non- spatial harvest scheduling problem with uncertain volume yields and demands, (b) a multi- objective problem with uncertain preferences, and (c) a spatial harvest scheduling and road building problem with uncertain volume yields. Deterministic and robust formulations of these problems are provided and the performance of their solutions is evaluated under simulated scenarios of uncertain data. In all cases, robust decisions are less sensitive to uncertain data and hence protected from the occurrence of infeasibilities, with a modest reduction in the objective function value. Moreover, deterministic and robust decisions greatly differ, suggesting that traditional solutions may require major corrections to adapt to changing future conditions with a consequent decrease in the quality of the decisions. The effect of the methodology assumptions are discussed and future work is suggested. ii Table of Contents Abstract.......................................................................................................................................... ii Table of Contents ......................................................................................................................... iii List of Tables ................................................................................................................................ vi List of Figures.............................................................................................................................. vii Acknowledgements ...................................................................................................................... ix Co-authorship Statement ............................................................................................................. x 1 Introduction............................................................................................................................ 1 1.1 Forest resources decision making and uncertainty ......................................................... 1 1.2 Incorporating uncertainty in forest management decision models ................................. 4 1.2.1 Stochastic programming............................................................................................ 4 1.2.1.1 General formulation .......................................................................................... 4 1.2.1.2 Applications ...................................................................................................... 6 1.2.2 Probabilistic programming........................................................................................ 7 1.2.2.1 General formulation .......................................................................................... 8 1.2.2.2 Applications ...................................................................................................... 9 1.2.3 Fuzzy set theory ........................................................................................................ 9 1.2.3.1 General formulation .......................................................................................... 9 1.2.3.2 Applications .................................................................................................... 11 1.2.4 Final remarks........................................................................................................... 12 1.3 Robust optimization ...................................................................................................... 13 1.3.1 Bertsimas and Sim robust optimization approach................................................... 15 1.4 Thesis objective ............................................................................................................ 19 1.5 Thesis structure ............................................................................................................. 19 1.6 References..................................................................................................................... 22 2 A Robust Optimization Approach to Protect Harvest Scheduling Decisions against Uncertainty.............................................................................................................. 28 2.1 Introduction................................................................................................................... 28 2.2 Robust optimization ...................................................................................................... 30 2.2.1 Basic concepts ......................................................................................................... 31 iii 2.2.2 Robust formulation.................................................................................................. 32 2.2.2.1 Formulating the protection function (S1) ....................................................... 33 2.2.2.2 Transforming the protection function into an optimization model (S2) ......... 34 2.2.2.3 Obtaining the dual model of the linear version of the protection function (S3) ................................................................................................... 35 2.2.2.4 Embedding the dual model in the original model (S4) ................................... 35 2.3 Methods......................................................................................................................... 36 2.3.1 Study area................................................................................................................ 37 2.3.2 Deterministic model formulation ............................................................................ 38 2.3.3 Uncertain data ......................................................................................................... 39 2.3.4 Robust model formulation....................................................................................... 40 2.3.4.1 Protection against minimum demand infeasibility ......................................... 40 2.3.4.2 Protection against infeasibility for fluctuation constraints ............................. 41 2.3.4.3 Controlling the protection level ...................................................................... 41 2.3.5 Simulation experiments........................................................................................... 42 2.4 Results........................................................................................................................... 43 2.5 Discussion..................................................................................................................... 49 2.6 Conclusions................................................................................................................... 53 2.7 References..................................................................................................................... 54 3 Bi-objective Multi-Period Planning with Uncertain Weights: A Robust Optimization Approach....................................................................................................... 58 3.1 Introduction................................................................................................................... 58 3.2 Robust optimization approach ...................................................................................... 60 3.3 Methods......................................................................................................................... 62 3.3.1 Study area................................................................................................................ 62 3.3.2 Bi-objective multi-period planning model .............................................................. 63 3.3.3 Modeling unclear preferences ................................................................................. 65 3.3.4 Uncertain weights and analysis ............................................................................... 66 3.4 Results........................................................................................................................... 68 3.5 Discussion..................................................................................................................... 72 3.6 Conclusions................................................................................................................... 75 3.7 References..................................................................................................................... 77 4 A Robust Model to Protect Road Building and Harvest Decisions from Timber Estimate Errors.................................................................................................................... 79 4.1 Introduction................................................................................................................... 79 4.2 Approaches to deal with uncertainty............................................................................. 80 4.3 The robust optimization approach ................................................................................ 81 iv 4.4 Methods......................................................................................................................... 84 4.4.1 Study area................................................................................................................ 84 4.4.2 Harvest scheduling and road building models ........................................................ 85 4.4.2.1 Deterministic model........................................................................................ 87 4.4.2.2 Robust model .................................................................................................. 89 4.4.3 Uncertain coefficients and model parameters ......................................................... 90 4.4.4 Simulation experiments........................................................................................... 91 4.5 Results........................................................................................................................... 92 4.6 Discussion..................................................................................................................... 96 4.7 Conclusions................................................................................................................... 99 4.8 References................................................................................................................... 101 5 Conclusions......................................................................................................................... 104 5.1 General conclusions .................................................................................................... 104 5.2 Strengths and weaknesses of the research .................................................................. 107 5.3 Future research and potential applications of this research ........................................ 109 5.4 References................................................................................................................... 111 Appendices................................................................................................................................. 114 Appendix A.......................................................................................................................... 114 v List of Tables Table 2.1. Differences in harvest decisions for robust models with different protection levels... 48 Table 2.2. Size and solution time for the problems (reported by the software)............................ 49 Table 3.1. Value of the objective function and scores produced by deterministic and robust solutions.. ..................................................................................................................... 70 Table 4.1. Numerical results of the optimization process............................................................. 92 Table 4.2. Summary of harvest results for the planning models (SCE10 / SCE15) ..................... 96 vi List of Figures Fig. 2.1. Initial age class distribution shows a mature forest with almost 80% of the area in stands at least 40 years-old........................................................................................... 37 Fig. 2.2. The non-restricted timber flow (a) was constrained by minimum demand and maximum negative fluctuation between period constraints (b) to define the base case for the analysis. .................................................................................................... 44 Fig. 2.3. Lower average infeasibility rates were obtained with the robust model (ROB) when the same buffer was used to modify the demand of the deterministic model (BUF). Although the objective value reduction is greater in the robust model for different protection levels, this reduction is lower when the same level of infeasibility rates is considered................................................................................................................. 45 Fig. 2.4. Increasing timber production in periods with active minimum demand constraints were obtained when more protection was required. This increase was clearly larger when volume and demands were uncertain (c and d) than when only volume was uncertain (a and b)........................................................................................................ 46 Fig. 2.5. Reductions in the average percentage of infeasible simulations were observed as higher protection levels were used when volume (a and b) and volume and demands (c and d) were assumed uncertain. The impact in the objective function was higher for higher discount rates (b and d)............................................................. 47 Fig. 2.6. Reductions in the average number of infeasible periods were also observed as higher protection levels were used when volume (a) and volume and demands (b) were uncertain. Similar results were obtained for a 2% discount rate (not shown).............. 48 Fig. 3.1. The two scenarios described two levels of uncertainty for the objective preferences over time. ..................................................................................................................... 67 Fig. 3.2. The score dispersion (gray area) was wider for the deterministic solutions (a and c) than for robust solution ROB30 (b and d) when optimal decisions (black line) were evaluated using simulated weights............................................................................... 69 Fig. 3.3. Changes in harvest decisions relative to deterministic decisions were larger as more robustness was required and more uncertainty was assumed (b)................................. 70 Fig. 3.4. In both scenarios of uncertainty, more employment (a,c) and lower proportions of old-aged forest (b,d) were produced as more robustness was required. ...................... 71 vii Fig. 3.5. Area harvested through the planning horizon in the two scenarios of uncertainty. ....... 72 Fig. 4.1. A harvestable area of 7,554 ha was used as study area. It included 431 stands and 412 potential roads. ...................................................................................................... 85 Fig. 4.2. In a tree structure (a) only one flow direction is possible so the maximum flow in an arc can be determined as a cumulative flow (Max Flowcb = Fc; Max Flowba = Fc + Fd). If cycles are present (b), the maximum flow in an arc is more difficult to determine as complex flow interactions might result. (Andalaft et al. 2003) .............. 88 Fig. 4.3. Reductions in the percentage of infeasible simulations and in the objective function were observed in both scenarios when higher protection levels were used. The buffer strategy showed higher infeasibility rates and higher losses in the objective function compared to robust decisions......................................................................... 93 Fig. 4.4. Harvest levels of robust optimization models increased in both scenarios when a higher protection level (PL) was used. Lower timber harvest was required when the volume estimate error was smaller (SCE10).......................................................... 94 Fig. 4.5. Harvest and road building decisions were different among deterministic (DET and BUF) and robust models (ROB). Buffer strategy (BUF) tended to produce a similar decision pattern as the deterministic model (DET).......................................... 95 viii Acknowledgements This work would not have been possible without the financial support of the Natural Sciences and Engineering Research Council of Canada, which has financially supported me throughout all my Ph.D. studies. I want to express my gratitude to my thesis committee members: Dr. John Nelson, Dr. Thomas Maness and Dr. Peter Marshall, for their time, suggestions and valuable comments on my work. I am especially grateful to my supervisor, Dr. John Nelson, for trusting me and giving me the opportunity of pursuing my Ph.D. under his supervision. Thanks for being there when I needed your help and for guiding me while allowing me the room to work in my own way. I want to extend my thanks to my parents, for their constant support and love from afar; to my parents-in-law and my brother-in-law, for their understanding and for being always present; and to all the friends that have made my stay in Vancouver more pleasant. Finally, I want to dedicate this work and express my heartfelt gratitude to my wife Carolina and my daughter Tiara. Without their love, support, encouragement and incredible patience I would not have made it through this program so smoothly. Part of this Ph.D. belongs to them. ix Co-authorship Statement I want to recognize the contribution of Dr. John Nelson in providing me the data used to run the models and in the preparation of all the manuscripts contained in this dissertation. x 1 Introduction 1.1 Forest resources decision making and uncertainty The forest management planning process consists of providing decision makers with the information needed to select an appropriate course of action. Key decisions to be made have traditionally included land allocation to different uses, silvicultural practices, timber harvest scheduling, road building and environment protection, among others. Most recently, forest management has gradually incorporated non-timber products and ecosystem services reflecting the society’s increasing demand for a more sustainable and diverse use of forest ecosystems. The complex structure of forests, the long-term effect of the decisions, the diverse nature of objectives, and legal constraints regarding land use and management practices make this a complicated problem in which decision-making tools play an important role (Davis et al. 2001). Due to this complexity, forest decision-making is divided into a hierarchical structure (Weintraub and Cholaky 1991, Gunn 1991) in which each of the three levels in this hierarchy focuses on specific questions that usually concern different decision-makers. Consequently, decision-making tools and data requirements are different for each hierarchical level (Nelson 2001). At the highest level in this hierarchy, strategic planning looks for general strategies to be applied to the landbase and explores which outputs should be produced over a long-term planning horizon. For instance, allowable cut and conservation targets are usually set at this level. At a tactical level, strategic goals need to be translated to specific, spatially located management units in a shorter time frame. Although there exist an increasing trend to include spatial detail in strategic decisions (Nelson 2003), it is in the tactical level in which spatial detail supports conservation goals and forest management decisions like harvesting, thinning, fuel treatment (Church 2007), road management decisions, and adjacency issues (Davis et al. 2001). Finally, the operational level deals with real-world operative actions like a detailed schedule of stand harvests, wood supply to specific facilities, transportation/routing decisions, and production planning at mills (D`Amours et al 2008). A variety of decision-making tools are used to assist decision-makers with these complex decisions. For a good overview of the main models 1 used in forest planning the reader is referred to the forestry section of Weintraub et al. (2007). Mathematical programming models, including linear, mixed-integer and dynamic programming are among the most widely used techniques. Large problems, especially those considering spatial issues, usually become hard to solve by traditional solution techniques. In this case heuristics (e.g., simulated annealing and tabu search) that quickly find good rather than optimal solutions can be used (Bettinger et al. 2002). Part of the complexity of the forest decision-making process comes from the inherent uncertainty of the information needed to make a decision. Forest management decisions often concern large areas, long time horizons and multiple stakeholders, and the sources of uncertainty include social, economic, biological, technological and catastrophic factors (Mowrer 2000, Regan et al. 2002, Kangas and Kangas 2004, Marshall 1987). In most of the cases, decision models assume all information to be deterministic (Kangas and Kangas 1999, Lohmander 2000). This means that any plan or decision implemented now will probably be sub-optimal in the ex-post analysis (Marshall 1987), since all assumptions and initial parameters will not necessarily be observed. In addition, perturbations of supposedly ‘certain’ data can make an optimal solution infeasible (Thompson and Haynes 1971, Hof et al. 1988, Pickens and Dress 1988). It is important to recognize that mathematical models in a decision-making context are abstract representations of the real, complex systems of decisions. The most relevant aspects of the system are described by using a set of (in)equalities for which the best solution is found based on certain criteria. Thus the solution to these models should only serve as a guide when a decision is made and implemented. Moreover, model infeasibilities do not necessarily mean that the system may become infeasible, but that the consequences of applying the model solution to the system may be awkward and expensive. This is particularly true in the forest management context where the decision structure allows the possibility of recourse. In a sequential decision-making framework, the first-period solution of the model is the basis for the decisions to be implemented. Once the outcomes of these decisions and some of the uncertainties are realized the model is run again and the decisions modified accordingly. This rolling planning horizon approach allows decision-makers to deal with uncertainty in a practical way. However, the more the changes in the decisions, the higher the cost of adapting to uncertainty. Less sensitive solutions to uncertain model inputs will therefore translate into a lower cost of adaptation. 2 A traditional approach to study the stability of optimal solutions to data perturbations is the post- optimal or sensitivity analysis, which explores the manner in which changes in individual data modify the optimal solution. This can be very useful to address efforts towards better estimation of key parameters, but its application in real world problems is limited (Pickens and Dress 1988). Another way to consider uncertainty in the planning process has been the use of scenario analysis or simulation. Since we know that different conditions (e.g., changing timber prices, uncertain yield) or events (e.g., fires occurrences or insect attacks) can occur, alternative scenarios of these conditions and events are evaluated. This allows managers to determine the most reasonable management option (Mowrer 2000) and to answer ‘what if’ questions relating to a particular path of a given resource (von Gadow 2000). This approach has been extensively used in forestry (von Gadow 2000, Peter and Nelson 2005, Smith and Zollner 2005, Klenner et al. 2000), although its application usually requires a large number of runs, and its usefulness is limited to relatively stable environments (Courtney et al. 1997). However, these approaches are not methodologies to explicitly deal with uncertainty, but ways to manually explore it. In this sense, they do not directly provide decision makers with good enough decisions for a range of uncertain values. Current methodologies that explicitly consider uncertainty are discussed in the next section. As new approaches to consider uncertainty remain a challenge (Martell et al. 1998, Nelson 2003), it is relevant for researchers to know the general framework, advantages and drawbacks of these methodologies in order to conduct research on this topic. Although the term ‘uncertainty’ is traditionally associated with a total ignorance about the future and ‘risk’ is referred to a quantified uncertainty (Knight 1921), I use ‘uncertainty’ here simply as a lack of certainty, measurable or not, and ‘risk’ as a state of uncertainty associated with outcomes that have an undesired effect. In the following sections, the main approaches used to explicitly consider uncertainty in decision models are presented, and the robust optimization approach (Bertsimas and Sim 2004) used in this research is described. This methodology is applied in the following chapters to some of the decision problems described above. 3 1.2 Incorporating uncertainty in forest management decision models In this section 1 describe current methodologies that consider uncertainty in the modeling process and that have been used in forest resources management (e.g., stochastic programming, chance-constrained programming, fuzzy sets, and information gap theories). Description, general formulation, and applications in forest resources management problems are provided for each methodology. 1.2.1 Stochastic programming Stochastic programming (SP) is a modeling technique to explicitly consider uncertainty in optimization models. It was introduced in the 1950s by Dantzig (1955) and Beale (1955). The basic concept in SP is the opportunity to adjust decisions to the information that is received after a random event or scenario has occurred, which is referred as ‘recourse’. In other words, a set of decisions (known as first stage decisions) are made before the realization of the uncertain parameters, and another set of decisions (second stage decisions) are made after that, in order to make corrections of previous decisions in the light of the realization of uncertainty (Sahinidis 2004). When uncertainty is revealed sequentially over time, a multi-stage approach is more appropriate; however, this approach has scarcely been used in forestry (see Gassmann (1989) for an application). 1.2.1.1 General formulation The general formulation of a discrete two-stage SP is as follows (Birge and Louveaux 1997): [ ]T bAx ≥]2.1[ 0]3.1[ ≥x T ),(]1.1[ sxhExcMin + subject to where the recourse function E[.] represents the expected value and h(x,s) equals ss yqMin]4.1[ 4 subject to xThyW −≥]5.1[ 0]6.1[ ≥y + bAx ≥]8.1[ ssss s Each possible course of action is represented by a scenario s in this formulation, and x and y represent the decision vectors of the first and second stage, respectively. Note that only variable y responds to scenario s and to the decision in the first stage, x. Matrices c, b, and A correspond to the fixed data of the first-stage decision. In the second-stage, a number of random events may occur, and for a given realization of s, the second-stage data qs, hs , Ts and Ws become known (Birge and Louveaux 1997). The objective function of SP minimizes the cost of the first stage decision (deterministic term) plus the expected cost of the second-stage, while satisfying the first stage constraints [1.2]. The expected cost depends on the random event or scenario s that will happen in the future and on the first stage decision x. In other words, h(x,s) looks for the optimal solution y for each scenario s (for each scenario, y is actually the solution of a linear programming model), subject to a set of recourse constraints that make ‘corrective actions’ for each scenario given the first stage decision x (Birge and Louveaux 1997). The greatest difficulty with SP is to evaluate the recourse function, for which a probabilistic description of future events is required. Random events may constitute a very large set of scenarios in the case of discrete random variables, or a continuum in the case of continuous random variables (Sen and Higle 1999). However, when a small set of scenarios can be defined SP can be solved by formulating the Deterministic Equivalent Problem (DEP, Walkup and Wets 1967), a deterministic model in which a problem is formulated for each possible scenario. In other words, assuming that scenario s will occur with probability ps, the stochastic problem is explicitly expressed in its full form: ∑ ∈Ss ssss yqcxpMin )(]7.1[ subject to s 5 SsxThyW ∈∀−≥]9.1[ Ssxx ∈∀=− 0]10.1[ Ssyx ∈∀≥ sssss s ss 0,]11.1[ The first-stage decision xs is defined for each scenario s, which may produce s different decisions for the first stage. However, one first-stage decision, independent of the second-stage scenario that actually occurs, must be found, which is forced through the constraint xs – x = 0, ∀s. A more compact formulation by defining a single variable x for the first stage (instead of xs) is recommended, but I used xs here for clarity. The above formulation ensures feasibility for all second-stage constraints, even for the unlikely scenarios. This produces what has been called ‘fat solutions’ (Madansky 1962), very conservative and expensive solutions. To avoid these worst-scenario solutions, penalty methods can be used in which violations of some constraints are allowed and penalized in the objective function (Birge and Louveaux 1997, Kall and Wallace 1994). In real problems, the number of scenarios may be extremely large and, as a result, models become computationally difficult to solve. In this case, efficient solution methods must exploit the special structure of stochastic programs. These methods mainly include decomposition, where second-stage problems become sub-problems, and statistical methods where the recourse function, E[h(x,s)], is approximated by statistical estimations (Birge and Louveaux 1997, Kall and Wallace 1994). For a more detailed description of SP the reader is referred to Kall and Wallace (1994), Birge and Louveaux (1997), Sen and Higle (1999) and the many resources available in the Stochastic Programming Community webpage (2007). 1.2.1.2 Applications Although the title of many publications suggests the opposite, SP has not been widely used in forestry. A two-stage stochastic model applied to the harvest scheduling problem with uncertain forest outputs was presented in Hoganson and Rose (1987). A Model I formulation (Johnson and 6 Scheurman 1977) was used to represent short-run (first-stage decisions) and long-run (second- stage decisions) prescriptions, and a decomposition-based approach was used to solve the model. The harvest scheduling problem was also formulated as a multi-stage stochastic model in Gassmann (1989). In this case, the author used different proportions of timber burned by fires as scenarios, and formulated a small version of the model presented in Reed and Errico (1986). To avoid the infeasibility produced in worst case scenarios, a penalty term that allowed violations of the flow constraints was included. A similar model was formulated by Boychuk and Martell (1996). In their model, the first few periods of the planning horizon had stochastic fire losses (two scenarios with high and low proportion burned), and for the remaining periods the expected fire loss was considered. The volume flow constraint was also encouraged by penalizing deviations from the desired volume. They solved the model using a traditional linear programming solver by formulating it as a Deterministic Equivalent Problem. Uncertainty in forest growth was considered as a multi-stage model in Eriksson (2006) by modifying a Model I formulation. Four scenarios were considered and a small sample problem was solved using the DEP formulation. Road upgrade decisions have also been examined using a two-stage stochastic model (Olsson 2007). The uncertain length of the period with poor road conditions was represented by three equiprobable scenarios and the DEP was solved using a commercial linear programming solver. 1.2.2 Probabilistic programming Although probabilistic or chance constrained programming (CCP) is sometimes referred as a kind of stochastic programming, here it is treated as a separate methodology. Different assumptions, mathematical structure of the resulting models, and solution techniques used to solve them support this distinction. In CCP, first proposed by Charnes and Cooper (1959), constraints with at least one random coefficient are modeled as probabilistic statements and are required to hold with a minimum probability. As in SP, probability distributions of uncertain coefficients are assumed to be known. Unlike SP, in this approach no recourse or ‘corrective actions’ are explicitly assumed. 7 1.2.2.1 General formulation In CCP, the set of constraints that contains random coefficients can take the form of individual or joint chance constraints. In the first case, constraints are of the form ibxa ij ijij ∀≥∑ ( ) ibxa iij ijij ∀≥≥∑ αPr]12.1[ where Pr(.) denotes a probability function and αi∈(0,1] is the minimum requirement, exogenously determined, on the probability of satisfying constraint i. Each constraint is transformed into a chance constraint individually and required to hold with a probability αi. In the second case, constraints are transformed in the single joint chance constraint ( ) α≥∀≥∑ ibxa ij ijijPr]13.1[ where 0 < α ≤ 1 is the minimum requirement on the probability of satisfying all constraints. Individual chance formulation [1.12] guarantees that the probability of violation of each constrain remains small, but the probability of violation of several constraints may remain high. The latter point is addressed using the joint formulation [1.13]. Although this formulation produces a single constraint, it is harder to solve since it is required to deal with multivariate distributions (Birge and Louveaux 1997). When only b is random in the individual chance constraint formulation, the deterministic equivalent can be obtained easily by replacing b with a new value based on its cumulative distribution, let us say F, given by F-1(α). In this case, the single constraint turns out to be a linear constraint. However, in most cases when A-matrix coefficients are uncertain or the joint chance constraint approach is used, deterministic equivalents become nonlinear models and nonlinear programming algorithms are required to solve them (Birge and Louveaux 1997, Kall and Wallace 1994). For a deeper description of CCP the reader may consult Kall and Wallace (1994), Birge and Louveaux (1997) and the Stochastic Programming Community webpage (2007). 8 1.2.2.2 Applications The first applications of CCP in forestry dealt with randomness in future timber growth. In Weintraub and Vera (1991) and Weintraub and Abramovich (1995) an individual chance constrained formulation was proposed and a cutting plane algorithm was used to solve the resulting nonlinear problem. Also an individual chance constrained approach was used in Pickens et al. (1991), and a simulation approach was applied to solve the resulting model. In Hof et al. (1992) individual and joint chance constrained models were formulated, and a reduced gradient algorithm was applied. In Hof and Pickens (1991) the chance of meeting random right-hand sides (output targets) is maximized by comparing three different objectives (i.e., maximize the minimum probability of meeting a target, maximize the joint probability that some targets are met, and maximize the total probability of meeting some targets). Models were solved using a nonlinear algorithm as well as linear piecewise approximations. 1.2.3 Fuzzy set theory Fuzzy set theory is an extension of the classical notion of crisp membership of set theory to a more vague and fuzzy one. Rather than a dichotomic membership of an element, a membership function that takes values in the interval [0,1] represents a ‘grade of membership’ of the element. It was first proposed by Zadeh (1965) and applied in the decision-making context by Bellman and Zadeh (1970). It has been extensively applied in operations research methodologies, for instance, linear, non-linear and dynamic programming, graph theory, etc. (Zimmermann 1996). 1.2.3.1 General formulation In the context of decision-making, many variations of the fuzzy set concepts are possible (Zimmermann 1996), so only the most common approach is described next. The membership function of an object x, µ(x) ∈ [0,1], in a fuzzy set represents the grade of membership, from the lowest (0) to the highest (1), of x in a given set. As in traditional set theory, operations like intersection (logical ‘AND’ operator) and union (logical ‘OR’ operator) of fuzzy sets are defined and become relevant to the decision-making context. Interception is therefore defined as the MIN operator in a way that the membership function of the intersection 9 of two fuzzy sets A and B is µA∩B(x) = Min(µA(x), µB(x)). Union, on the other hand, is defined as the MAX operator. Other operators have also been defined (Bellman and Zadeh 1970, Zimmermann 1996). In the context of fuzziness, a fuzzy goal is identified as a fuzzy set from which a good enough solution is to be found. Borrowing from Bellman and Zadeh (1970), if x is a solution from the set of all possible solutions X, a fuzzy goal might be expressed in words as ‘x should be substantially larger than Z’, where Z is an aspiration level of the objective function. The membership function of the fuzzy goal, µG(x), can be defined, for instance, in a way that obtaining less than Z represents no grade of membership and values over Z represent increasing memberships, as follows: ⎩⎨ ⎧ >−+= −− ZxZxxG ,))(1()(]14.1[ 12µ ≤ Zx,0 Similarly, a membership function of a fuzzy constraint, µ(x), must be defined as an increasing monotonously function from 0 to 1, in a way that should be 0 if it is strongly violated and 1 if it is very well satisfied. Although different mathematical expressions can be used, the simplest type of function is linearly increasing in a tolerance interval p, as follows (Zimmermann 1996): pbAxif pbAxbif bAxif p bAxx +> +≤< ≤ ⎪⎪⎩ ⎪⎪⎨ ⎧ −−= 0 1 1 )(]15.1[ µ Under this approach, objectives and constraints are treated similarly in the formulation of a decision, as they are of the same importance, unlike the conventional approach in which the objective is explicitly more important than constraints. Nevertheless, when some objectives or constraints are more important than others, a weighting coefficient can represent their relative importance (Bellman and Zadeh 1970). To make a decision, fuzzy goals and fuzzy constraints need to be aggregated in order to rank the decision alternatives, and operators like those mentioned above are used. Let µi(x) be the membership function of row i, including the objective function and constraints. The membership function of the fuzzy decision set can be described, for example, by the intersection of all the membership functions. This means that { })(min)( xx iiD , and the solution of the decision problem can be seen as maximizing the µµ = 10 { })(xdegree of membership of this decision set, or max 0 Dx µ > . This is the traditional maxmin problem, in which the minimum degree of membership of the all the rows is maximized: λMax]16.1[ ix ∀≤ )(]17.1[ subject to iµλ 0]18.1[ ≥x The interception operator is the most widely used, although it does not allow for compensation (Zimmermann 1996). That is, a low membership degree of an element cannot be compensated by a high degree of membership of another. The operator is not able to discriminate among solutions that differ with respect to the fulfillment of membership to the various constraints, except for the smallest membership degrees (Dubois et al. 1996). Different operators, criteria for selecting an appropriate operator as well as different membership functions are discussed in Zimmermann (1996). 1.2.3.2 Applications Harvest scheduling with fuzzy timber yields is among the first applications of fuzzy sets in forestry. In Hof et al. (1986) and Pickens and Hof (1991) a MAXMIN formulation of the harvest scheduling problem was formulated to produce non-declining timber flow assuming fuzzy yield. Similarly, Bare and Mendoza (1992) compared a crispy and a fuzzy formulation considering fuzziness in the non-declining flow constraints. By relaxing these constraints, as expected, increases in objective values were obtained. Fuzzy theory was also used in Mendoza and Sprouse (1989) as a methodology to generate more robust management alternatives to a multiple-use forest planning problem. In addition, Boyland et al. (2006) compared the effect of using crispy and fuzzy seral-class definitions in harvest scheduling models, dealing with the vagueness of the traditional age class boundaries. In the context of multi-objective models, a fuzzy approach was used in Ells et al. (1997), where the land use allocation problem with vagueness in the objectives definition and in the effect of management options on these objectives was addressed, and in Maness and Farrell (2004) where 11 the minimum target levels of different goals or indicators were not known with certainty. Similarly, the vagueness of the acceptance thresholds of different criteria to evaluate management alternatives was addressed in Kangas et al. (2006). Uncertainty and fuzziness of the objective function coefficients were also considered in Mendoza et al. (1993) and Stirn (2006). 1.2.4 Final remarks Uncertainty has reappeared as a relevant topic in decision making modeling. The most traditional methodologies to explicitly include uncertainty into mathematical modeling in forestry applications were broadly summarized here, and some of their applications were provided. The aim of any methodology that explicitly incorporates uncertainty should be to provide decision makers with more robust or stable decisions than ‘deterministic’ methodologies. A more robust or stable decision should be one that guarantees a desired level of output regardless of the value of the parameters that define the performance of the system. Since robust decisions are required to be good (not optimal) for multiple possible future values of uncertain parameters, their objective values will certainly be negatively affected. However, decision makers should be willing to assume this cost in exchange for a reduction in the risk associated with their decisions. The approaches discussed have advantages and drawbacks, essentially based on their assumptions and the effort needed to find solutions, which defines the appropriateness of their use in different decision problems. Stochastic programming, for example, provides an intuitive and very elegant representation of an uncertain future through the use of scenarios, for which a quantitative description of their probability of occurrence is needed. Unlike other approaches, it represents the dynamic nature of the decisions by allowing the possibility of recourse. However, as the number of scenarios increases this representation produces large models very quickly, forcing the use of specially designed solution algorithms that require additional implementation efforts. The degree of risk aversion (tradeoff between a more certain output with a lower objective value and a more uncertain one with a better objective value) is not explicitly considered in stochastic programming. Although chance-constrained models are also quite intuitive, they usually imply the use of nonlinear solution techniques except for specific problems that can be easily solved. As in stochastic programming, the probability distribution of uncertainties must be known, but in this approach risk aversion levels can be explicitly defined. 12 Fuzzy set theory, on the other hand, is designed to deal with vagueness and fuzziness rather than uncertainty. It is more appropriate when there is vagueness, for example, in the meaning of certain events, phenomena or statements, but not when representing the lack of information about the value of the parameters. When used improperly, this methodology actually results in a relaxed version of traditional decision problems, in which constraint violations are allowed and, consequently, better objective functions are obtained. The way in which uncertainty is modeled in decision problems with fuzzy set theory does not actually produce stable solutions. It should be noted here that a large group of other probabilistic modeling approaches has also been used to model uncertainty. For instance, Markov based decision models have been used to evaluate harvest decisions when the forest growth (Zhou and Buongiorno 2006), the timber price (Teeter and Caulfield 1991, Insley and Rollins 2005), or both (Lembersky and Johnson 1975, Kaya and Buongiorno 1987) are uncertainty. In these cases, Markov chains are used to represent the transition among uncertain future states, and other techniques are used to solve the resulting models (Hillier and Lieberman 2005). Research on new methodologies to include uncertainty in planning models is highly relevant. With more computer power available today, better forecast techniques to estimate some parameters, and new recognized sources of uncertainty (e.g., climate change, demand for new environmental services, technology, etc.), the decision process can benefit from research on this topic. These methodologies should consider the tradeoff between the risk reduction and the cost of this reduction (in other words, risk-aversion level), and different approaches should be applied to different kinds of uncertainties (i.e., catastrophic versus non-catastrophic events). As some kinds of uncertainty can be well described by distribution functions, and others not, approaches that can deal with both these issues are important. 1.3 Robust optimization The first efforts to deal with decisions with limited information were probably made by Wald in the late 1930’s (Wald 1939) and middle 1940’s (Wald 1945a, 1945b) in the field of statistical decision-making. In his work, Wald defined loss functions, risk functions and introduced the minmax principle, in which the course of action that minimizes the maximum loss should be preferred. In the context of hypothesis testing, he created the concept of sequential analysis, where the sample size is not fixed in advance but defined as the collected data are evaluated 13 (Wald 1945c). In its basics this concept involves the reduction of uncertainty by knowledge of past events when choosing an action. His principles were later applied to the optimal control of stochastic systems (Sworder 1965, Witsenhausen 1968, Bar-Shalom and Sivan 1969), therefore extending the already existing theory of optimal control (Sussmann and Willems 1997), in which control and state trajectories for dynamic systems are to be determined over a period of time. In the middle of the 1970’s the idea of including unknown parameters appeared in optimization problems in a more static framework. Soyster (1973) and Falk (1976) assumed that uncertain data in linear optimization models belonged to a given interval and solved linear models with a modified feasible region. Solutions to these models were guaranteed to be feasible for all possible scenarios of uncertain data. After a period of apparent inactivity, Mulvey et al. (1995) used goal programming and a scenario-based description of the problem data to generate less sensitive solutions to the realization of the model data. Since then, the concept of robust optimization (RO) has been used to represent a wide modeling framework applied to explicitly include uncertainty in optimization models in cases where uncertain data are assumed to belong to a given set. Although the re-emergence of RO in recent decades has been principally due to applications in the engineering design field (Park et al 2006), its application has been extended to a wide range of fields, including dynamic decision problems. For instance, it has been applied to multi-period portfolio management (Bertsimas and Pachamanova 2008), inventory theory (Bertsimas and Thiele 2006) and water resources planning (Chung et al. 2009). The first RO approaches (Soyster 1973, Falk 1976) suffered from providing extremely conservative solutions; new approaches have focused on identifying conditions and uncertainty sets resulting in less conservative solutions and models that can be solved efficiently (Beyer and Sendhoff 2007). In a general constraint Ax ≤ b, x ∈ X, we might assume that uncertainty affects the coefficients A (without loss of generality), and that the feasibility of the constraint needs to be guaranteed for any possible value of A within an uncertain set U. In other words, the robust counterpart of the original constraint becomes Ax ≤ b, ∀ A ∈ U, x ∈ X. The degree of conservatism and the class of the robust counterpart that needs to be solved depends on how set U is specified. For instance, Soyster’s formulation (1973) assumed the worst possible value of each uncertain data, therefore considering the simplest uncertainty representation (box uncertainty), i.e. },,ˆˆ:{ jiaaaaaAaU ijijijijijij ∀+≤≤−∈= where ija â a and are the nominal value and the deviation from the nominal value of , respectively. The original constraint is ij ij 14 then reformulated as },;,ˆ{ jyxyibyaxa jjjij jijj jij ∀≤≤−∀≤+∑∑ in its robust version, where yj≥0 is an additional variable. Although the resulting optimization model is linear, solutions are protected for all data scenarios (including the most unlikely ones) and the objective function value is seriously affected (Bertsimas and Sim 2004). Ben-Tal and Nemirovski (1997, 1998, 1999), or similarly El Ghaoui and Lebret (1997), proposed less conservative formulations assuming that the uncertainty set U is ellipsoidal, therefore removing the less likely outcomes from consideration. In this case, },)()(:{ iaaVaaAaU ijijiijijij ∀≤−−∈= θ 21T − where V is a positive definite matrix and 0 ≤ θ ≤ 1 is a “safety parameter” that controls the degree of conservatism. Note that the variability of individual coefficients is combined within a constraint i to obatain a constraint-wise variability. The robust counterpart of the original constrain is then },,;,ˆˆ{ 22 jiyzxyibzayaxa ijijjijiijj ijij jijj jij ∀≤−≤−∀≤Ω++ ∑∑∑ , where the probability that constraint i is violated is at most exp(-Ωi2/2) (Ben-Tal and Nemirovski 1999). The drawback of using ellipsoidal uncertainty sets is that it increases the complexity of the problems (e.g., linear programming problems transform into second-order cone problems1). Bertsimas and Sim (2004), on the other hand, control the level of conservatism of the solutions based on polyhedral uncertainty sets, i.e. ;,,ˆˆ:{ jiaaaaaAaU ijijijijijij ∀+≤≤−∈= },ˆ|| iaaaj ijijij ∀Γ≤−∑ , where Γ is the parameter that controls the conservatism of the solutions. This definition of uncertainty sets is in fact a case of ellipsoidal uncertainty (Ben-Tal and Nemirovski 1999) that has the advantage that it preserves the complexity of the original problem (e.g., linear programming problems remain linear). This approach will be applied to different types of forest management problems in this thesis, and it is explained in more detail in the next section. 1.3.1 Bertsimas and Sim robust optimization approach This description closely follows Bertsimas and Sim (2004). Consider the nominal linear optimization problem {Max c’x | Ax ≤ b, l ≤ x ≤ u } where data in matrix A are assumed uncertain. Without loss of generality the objective function coefficients c are not subject to uncertainty since we can use the objective maximize z, add the constraint z - c’x ≤ 0 and thus include this constraint into Ax ≤ b. For each row i of the matrix A, Ji represents the set of 1 A second-order cone problem is a class of convex optimization problems that requires a greater computational effort to be solved than a linear programming problem (Alizadeh and Goldfarb 2003) 15 coefficients in i subject to uncertainty. Each coefficient , j ∈ Jija i, is modeled as a symmetric random variable that takes values in ]ˆ,ˆ aaaa[ ijijijij +− . Associated with each uncertain is a random variable ija ijijijij )( aaa ˆ−=η which takes values in [-1, 1]. For each row, the aggregation of this random variables could, in theory, take values between –card(Ji) and card(Ji), where card represents the number of elements in a set. A new parameter [ ])(,0 Jcard ii ∈Γ , not necessarily integer, can be used to control the robustness of the solution against the level of conservatism in a way that iJj iji Γ≤∑ ∈ η . The authors consider the following (non-linear) robust formulation of the original problem: [1.19] Max c’x subject to ⎣ ⎦ ⎣ ⎦( ) ibyayaxa n j i Sj jitiijijScardSJtJStSjij i i iiiiiiiii ∀≤ ⎭⎬ ⎫ ⎩⎨ ⎧ Γ−Γ++∑ ∑ = ∈Γ=∈⊆1 )(,\,|}{ ˆˆmax]20.1[ U jyxy ∀≤≤−]21.1[ juxl ∀≤≤ jy ∀≥ 0]23.1[ ),x( jjj jjj]22.1[ j where the second term in [1.20] is the so called protection function of constrain i, ii Γβ . If is integer then the ith constraint is protected by iΓ { }∑ ∈Γ=⊆=Γ iiiiii Sj jijScardJSSii xâmax),x( })(,|{β . In the general case, a constraint will be protected against all cases that up to of the uncertain coefficients are allowed to change and one extra coefficient changes by ( ⎣ ⎦iΓ ita iΓ - ) . Note that if the method provides the most conservative solution and if then ⎣ ⎦iΓ â )(Jcard=Γ 0=Γ 0),x( =Γ it ii i iiβ and the constraints are equivalent to that the nominal problem. For the tradeoff between robustness and conservatism of the solution can be explored. )(0 Jcard<Γ< ii 16 In order to reformulate model [1.19]-[1.23] as a linear model, the authors proved (see Bertsimas and Sim 2004) that given a vector x*, the protection function of the ith constraint ⎣ ⎦ ⎣ ⎦( ) ⎭⎬ ⎫ ⎩⎨ ⎧ Γ−Γ+=Γ ∑ ∈Γ=∈⊆ i i iiiiiiiii Sj jitiijijScardSJtJStSii xaxa ** )(,\,|}{ ˆˆmax),x(]24.1[ U β equals the objective function of the following linear optimization problem: ij Jj jij wxaMax i ∑ ∈ ˆ]25.1[ * Jjw ∈∀≤≤ 10]27.1[ [ ])(,0 Jcard∈Γ subject to i Jj ij i w Γ≤∑ ∈ ]26.1[ iij By strong duality, since problem [1.25]-[1.27] is feasible and bounded for all then its dual, showed below, is also feasible and bounded and their objective functions coincide: ii ∑ ∈ +Γ iJj ijii pzMin]28.1[ subject to ijijiji Jjxapz ∈∀≥+ ˆ]29.1[ iij Jjp ∈∀≥ 0]30.1[ 0]31.1[ ≥z i where zi and pij are the dual variables associated with constraints [1.26] and [1.27], respectively. This dual of the protection function defined by [1.25]-[1.27] is finally substituted to the robust non-linear formulation in [1.19]-[1.23] to produce the following linear robust counterpart of the nominal problem: 17 [1.32] Max c’x subject to ibpzxa j i Jj ijiijij i ∀≤+Γ+∑ ∑ = ∈1 ]33.1[ n Jjiyapz ∈∀≥+ ,ˆ]34.1[ jyxy ∀≤≤−]35.1[ juxl ∀≤≤]36.1[ Jjip ∈∀≥ ,0]37.1[ jy ∀≥ 0]38.1[ iz ∀≥ ijijiji jjj jjj iij j i 0]39.1[ The new terms added in [1.33] define the appropriate buffer to guarantee that the left-hand side of the constraint is feasible for different uncertainty realizations. At the same time, the robust counterpart introduces new bounds to decision variables that prevent them from being large in directions that have considerable uncertainty in the parameters. The approach proposed by Bertsimas and Sim (2004) considers the variability of individual coefficients in the same way as other RO approaches (Soyster 1973, Ben-Tal and Nemirovski 1997, El Ghaoui and Lebret 1997) (i.e., ]ˆ,ˆ[ aaaaa +−∈ ijijijijij ), but controls the degree of conservatism by defining the uncertainty set U differently. In opposition to the more general formulation of Ben-Tal and Nemirovski (1997) and El Ghaoui and Lebret (1997) that allows interactions among uncertain data (ellipsoidal sets), the Bertsimas and Sim (2004) approach represents a particular case in which uncertainty is not correlated. By doing so, a computationally tractable robust counterpart is obtained. This approach proceeds by bounding the number of coefficients in the data subject to variability, that is, a fraction of the uncertain coefficients are allowed to change. As one of the principles behind the approach is that the errors in the estimate will compensate for each other, a large number of parameters must be uncertain in each constraint for the resulting formulation to achieve a solution that is not too conservative. Probability bounds on the constraint violations 18 can also be obtained as in other approaches (Bertsimas and Sim 2004), but, as pointed out in Chapter 2, they are loose and might produce conservative solutions. 1.4 Thesis objective As mentioned in the previous sections, current approaches to explicitly consider uncertainty in decisions models impose some difficulties (e.g., probability distribution assumptions and complexity of the resulting models) and new methodologies should be explored. The objective of the work described in this thesis is to apply a recently developed technique of robust optimization (Bertsimas and Sim 2004) to different types of forest management decision problems. This technique considerably reduces the computational complexity of the resulting models, but relies on some simplifying assumptions about the structure of the uncertainties, which appropriateness to the forest management context should be explored. In this context, this dissertation has the following specific objectives: (1) to determine how the robust optimization approach that explicitly considers uncertainty in strategic and tactical forest decision problems performs relative to their deterministic counterparts. This includes solution times, objective function values and the ability of the model to stay feasible under various assumptions of uncertainty and protection levels; and (2) to determine the potential benefits and limitations of considering uncertainty in strategic and tactical forest decision problems by analyzing how robust solutions differ from their deterministic counterparts and what implications this has for management decisions. 1.5 Thesis structure To achieve the thesis objectives three different decision problems are considered: (a) a strategic non-spatial harvest scheduling problem, (b) a strategic, bi-objective problem, and (c) a tactical spatial harvest scheduling and road building problem. These problems differ in the type of decisions, coefficients affected by uncertainty and modeling technique required for their formulation (i.e., linear and mixed-integer programming). Each problem becomes an independent manuscript in this thesis. 19 Chapter 2 describes a non-spatial harvest scheduling problem with uncertain volume yield and demand. In this problem, the feasibility of the constraints on timber production levels is highly desirable despite the uncertainty observed in volumes and demands of two products over the entire 20-period planning horizon. Deterministic and robust models were applied to a 245,090 ha forest in British Columbia and the performance of these decisions was compared through the simulation of the uncertain data. I conducted two main analyses. The first demonstrates the advantage of using this methodology over the strategy of manually imposing a buffer to the original demand levels. For simplicity, I focused only on protecting minimum demand constraints when volume was uncertain. I used the buffer determined by the robust model for each constraint to correct the minimum demand level of the same constraint in the deterministic model. I then showed the difference in decision stability between using a (more complex) robust model and a traditional model with a corrected demand level that acts as a security buffer. In the second analysis I used a complete model with protection for minimum demand and volume fluctuation constraints. I explored infeasibility rates and the change in the objective function and harvest decisions when volume and demand were uncertain under two discount rates. Chapter 3 explores a robust formulation of a bi-objective, multi-period planning problem with uncertain weights, providing an example of a formulation with uncertainty and protection requirements in the objective function. Since weights are hard to estimate and play a major role in determining the model solutions, in this chapter I explored a way of formulating a model that explicitly considers the sensitivity of the solutions to the weights when looking for “optimal” decisions. I present a robust formulation to a bi-objective multi-period harvest scheduling problem where weights were not accurately known and were allowed to change within a pre- defined range throughout the planning horizon. The model was applied to a 245,090 ha forest in British Columbia, where the amount of employment and the proportion of old forest through the planning horizon are the two management objectives. Evaluated under simulated scenarios of weights, robust solutions produced a more stable weighted sum of the objectives through the planning periods than traditional deterministic solutions. I demonstrated that robust decisions perform almost as well as the traditional fixed-value weight approach in terms of the objective function value, however, they are more consistent in a scenario of uncertain future weights. Chapter 4 presents a robust formulation of a tactical harvest scheduling and road building problem where volume yield is uncertain and minimum demand constraints need to be guaranteed. Unlike the previous chapters, in this chapter I extended the application of the robust 20 approach to a mixed-integer model and discussed its consequences. Deterministic and robust models were applied to a 7,552 ha area in British Columbia, and their solutions were compared under simulated scenarios of volume yields. Although no exact comparison is possible as mixed- integer models are unlikely to be solved to optimality, I showed that robust solutions outperformed deterministic solutions both in terms of infeasibility rates and objective function value. However, I noted that the spatial independency of the volume assumed in the robust approach becomes a critical limitation of the methodology when addressing spatial problems, and further research in this regard is needed. In the concluding chapter I summarize the main findings of the thesis and discuss the strengths and weakness of the robust optimization approach proposed. Potential applications of the methodology and future research are also addressed. 21 1.6 References Alizadeh,F., Goldfarb, D. 2003. Second-order cone programming. Math.Program., Ser. B 95(1): 3-51 Bare, B.B., and Mendoza, G.A. 1992. Timber harvest scheduling in a fuzzy decision environment. Can. J. For. Res. 22: 423-428. Bar-Shalom, Y, and Sivan, R. 1969. On the optimal control of discrete time linear systems with random parameters. IEEE Trans. Autom. Control 14(1): 3-8. Beale, E.M.L. 1955. On minimizing a convex function subject to linear inequalities. Journal of the Royal Statistical Society. Series B (Methodological) 17: 173-184. Bellman, R.E., and Zadeh, L.A. 1970. Decision-making in a fuzzy environment. Management Science Series B-Application 17: B141-B164. Ben-Tal, A. and Nemirovski, A. 1997. Robust truss topology design via semidefinite programming. SIAM J. Optim. 7: 991-1016. Ben-Tal, A. and Nemirovski, A. 1998. Robust convex optimization. Math. Oper. Res. 23(4): 769-805. Ben-Tal, A. and Nemirovski, A. 1999. Robust solutions of uncertain linear programs. Oper. Res. Lett. 25: 1-13. Bertsimas, D., and Pachamanova, D. 2008. Robust multiperiod portfolio management in the presence of transaction costs. Comput. Oper. Res. 35: 3-17. Bertsimas, D., and Sim, M. 2004. The price of robustness. Oper. Res. 52: 35-53. Bertsimas, D., and Thiele, A. 2006. A robust optimization approach to inventory theory. Oper. Res. 54: 150-168. Bettinger, P., Graetz, D., Boston, K., Sessions, J., and Chung, W.D. 2002. Eight heuristic planning techniques applied to three increasingly difficult wildlife planning problems. Silva Fennica 36: 561-584. Beyer, H.G., Sendhoff B. 2007. Robust optimization - a comprehensive survey. Comput. Methods Appl. Mech. Engrg. 196(33-34): 3190-3218. Birge, J.R., and Louveaux, F. 1997. Introduction to stochastic programming. First ed. Springer, New York. Boychuk, D., and Martell, D.L. 1996. A multistage stochastic programming model for sustainable forest-level timber supply under risk of fire. For. Sci. 42: 10-26. Boyland, M., Nelson, J., Bunnell, F.L., and D'Eon, R.G. 2006. An application of fuzzy set theory for seral-class constraints in forest planning models. For. Ecol. Manage. 223: 395-402. 22 Charnes, A., and Cooper, W.W. 1959. Chance-constrained programming. Manage. Sci. 6: 73-79. Chung, G., Lansey, K., Bayraksan, G. 2009. Reliable water supply system design under uncertainty. Environ. Modell. Softw. 24: 449-462. Church, R.L. 2007. Tactical-level forest management models. Chapter 17 in Handbook on Operations Research in Natural Resources, Eds. A. Weintraub, C. Romero, T. Bjørndal and R. Epstein. Kluwer Academic Publisher, New York. Courtney, H., Kirkland, J., and Viguerie, P. 1997. Strategy under uncertainty. Harv. Bus. Rev. 75: 67-79. D`Amours, S., Ronnqvist, M., and Weintraub, A. 2008. Using operational research for supply chain planning in the forest products industry. INFOR 46: 265-281. Dantzig, G.B. 1955. Linear programming under uncertainty. Manage. Sci. 1: 197-206. Davis, L.S., Johnson, K.N., Bettinger, P.S. and Howard, T.E. 2001. Forest Management. Fourth Edition. McGraw-Hill. New York. 804 pp. Dubois, D., Fargier, H., and Prade, H. 1996. Refinements of the maximin approach to decision- making in a fuzzy environment. Fuzzy Sets Syst. 81: 103-122. El Ghaoui, L., and Lebret, H. 1997. Robust solutions to least-squares problems with uncertain data. SIAM Journal on Matrix Analysis and Applications 18: 1035-1064. Ells, A., Bulte, E., and vanKooten, G.C. 1997. Uncertainty and forest land use allocation in British Columbia: Vague priorities and imprecise coefficients. For. Sci. 43: 509-520. Eriksson, L.O. 2006. Planning under uncertainty at the forest level: A systems approach. Scand. J. For. Res. 21: 111-117. Falk, J.E. 1976. Exact solutions of inexact linear programs. Oper. Res. 24: 783-787. Gassmann, H.I. 1989. Optimal harvest of a forest in the presence of uncertainty. Can. J. For. Res. 19: 1267-1274. Gunn, E.A. 1991. Some aspects of hierarchical planning in forest management. In proceedings of the 1991 Symposium on Systems Analysis in Forest Resources, Compiled by M.A. Buford. USDA Forest Service General Technical Report SE-74. pp 54-62. Hillier, F.S., and Lieberman, G. 2005. Introduction to operations research. 8th Edition. McGraw Hill. New York. Hof, J.G., and Pickens, J.B. 1991. Chance-constrained and chance-maximizing mathematical programs in renewable resource-management. For. Sci. 37: 308-325. Hof, J.G., Kent, B.M., and Pickens, J.B. 1992. Chance constraints and chance maximization with random yield coefficients in renewable resource optimization. For. Sci. 38: 305-323. Hof, J.G., Robinson, K.S., and Betters, D.R. 1988. Optimization with expected values of random yield coefficients in renewable resource linear-programs. For. Sci. 34: 634-646. 23 Hof, J.G., Pickens, J.B., and Bartlett, E.T. 1986. A maxmin approach to nondeclining yield timber harvest scheduling problems. For. Sci. 32: 653-666. Hoganson, H.M., and Rose, D.W. 1987. A model for recognizing forestwide risk in timber management scheduling. For. Sci. 33: 268-282. Insley, M., and Rollins, K. 2005. On solving the multirotational timber harvesting problem with stochastic prices: a linear complementarity formulation. Amer. J. Agr. Econ. 87(3): 735- 755. Johnson, K.N., and Scheurman, H.L. 1977. Techniques for prescribing optimal timber harvest and investment under different objectives – discussion and synthesis. For. Sci. 1-31. Kall, P., and Wallace, S.W. 1994. Stochastic programming. Second ed. John Wiley & Sons, Chichester. Kangas, A., Kangas, J., and Laukkanen, S. 2006. Fuzzy multicriteria approval method and its application to two forest planning problems. For. Sci. 52: 232-242. Kangas, A.S., and Kangas, J. 2004. Probability, possibility and evidence: approaches to consider risk and uncertainty in forestry decision analysis. For. Policy Econ. 6: 169-188. Kangas, A.S., and Kangas, J. 1999. Optimization bias in forest management planning solutions due to errors in forest variables. Silva Fenn. 33: 303-315. Kaya, I., and Buongiorno, J. 1987. Economic harvesting of uneven-aged Northern hardwood stands under risk: a Markovian decision model. For. Sci. 33: 889–907. Klenner, W., Kurz, W., and Beukema, S. 2000. Habitat patterns in forested landscapes: management practices and the uncertainty associated with natural disturbances. Comput. Electron. Agric. 27: 243-262. Knight, F.H. 1921. Risk, uncertainty and profit. New York, Houghton Mifflin Co, Boston. Lembersky, M.R., and Johnson, K.N. 1975. Optimal policies for managed stands: an infinite horizon Markov decision process approach. For. Sci. 21: 109–122 Lohmander, P. 2000. Optimal sequential forestry decisions under risk. Ann. Oper. Res. 95: 217- 228. Madansky, A. 1962. Methods of solution of linear programs under uncertainty. Oper. Res. 10: 463-471. Maness, T., and Farrell, R. 2004. A multi-objective scenario evaluation model for sustainable forest management using criteria and indicators. Can. J. For. Res. 34: 2004-2017. Marshall, P.L. 1987. Sources of uncertainty in timber supply projections. For. Chron. 63: 165- 168. Martell, D.L., Gunn, E.A., and Weintraub, A. 1998. Forest management challenges for operational researchers. Eur. J. Oper. Res. 104: 1-17. 24 Mendoza, G.A., and Sprouse, W. 1989. Forest planning and decision-making under fuzzy environments - An overview and illustration. For. Sci. 35: 481-502. Mendoza, G.A., Bare, B.B., and Zhou, Z.H. 1993. A fuzzy multiple objective linear- programming approach to forest planning under uncertainty. Agricultural Systems 41: 257-274. Moilanen, A., Runge, M.C., Elith, J., Tyre, A., Carmel, Y., Fegraus, E., Wintle, B.A., Burgman, M., and Ben-Haim, Y. 2006. Planning for robust reserve networks using uncertainty analysis. Ecol. Model. 199: 115-124. Mowrer, H.T. 2000. Uncertainty in natural resource decision support systems: sources, interpretation, and importance. Comput. Electron. Agric. 27: 139-154. Nelson, J. 2001. Assessment of harvest blocks generated from operational polygons and forest- cover polygons in tactical and strategic planning. Can. J. For. Res. 31: 682-693. Nelson, J. 2003. Forest-level models and challenges for their successful application. Can. J. For. Res. 33: 422-429. Olsson, L. 2007/5. Optimal upgrading of forest road networks: Scenario analysis vs. stochastic modelling. For. Policy Econ. 9: 1071-1078. Park, G.J., Lee, T.H., Lee, K., Hwang, K.H. 2006. Robust design: an overview. American Institute of Aeronautics and Astronautics Journal 44: 181-191. Peter, B., and Nelson, J. 2005. Estimating harvest schedules and profitability under the risk of fire disturbance. Can. J. For. Res. 35: 1378-1388. Pickens, J.B., and Hof, J.G. 1991. Fuzzy goal programming in forestry - An application with special solution problems. Fuzzy Sets Syst. 39: 239-246. Pickens, J.B., and Dress, P.E. 1988. Use of stochastic production coefficients in linear- programming models - objective function-distribution, feasibility, and dual activities. For. Sci. 34: 574-591. Pickens, J.B., Hof, J.G., and Kent, B.M. 1991. Use of chance-constrained programming to account for stochastic variation in the A-matrix of large-scale linear programs: A forestry application. Ann. Oper. Res. 31: 511-526. Reed, W.J., and Errico, D. 1986. Optimal harvest scheduling at the forest level in the presence of the risk of fire. Can. J. For. Res. 16: 266-278. Regan, H.M., Colyvan, M., and Burgman, M.A. 2002. A taxonomy and treatment of uncertainty for ecology and conservation biology. Ecol. Appl. 12: 618-628. Regan, H.M., Ben-Haim, Y., Langford, B., Wilson, W.G., Lundberg, P., Andelman, S.J., and Burgman, M.A. 2005. Robust decision-making under severe uncertainty for conservation management. Ecol. Appl. 15: 1471-1477. 25 Sahinidis, N.V. 2004. Optimization under uncertainty: state-of-the-art and opportunities. Comput. Chem. Eng. 28: 971-983. Sen, S., and Higle, J.L. 1999. An introductory tutorial on stochastic linear programming models. Interfaces 29: 33-61. Smith, W.P., and Zollner, P.A. 2005. Sustainable management of wildlife habitat and risk of extinction. Biol. Conserv. 125: 287-295. Soyster, A.L. 1973. Convex programming with set-inclusive constraints and applications to inexact linear-programming. Oper. Res. 21: 1154-1157. Stirn, L.Z. 2006. Integrating the fuzzy analytic hierarchy process with dynamic programming approach for determining the optimal forest management decisions. Ecol. Model. 194: 296-305. Stochastic Programming Community Home Page. 2007. http://stoprog.org. Accessed on November 02, 2007. Sussmann, H.J., and Willems, J.C. 1997. 300 years of optimal control: from the brachystochrone to the maximum principle. IEEE Contr. Syst. Mag. 17(3): 32-44. Sworder, D.D. 1965. Minimax control of discrete time stochastic systems. J. SIAM Control 2(3): 433-449. Taguchi, G. 1995. Quality engineering (Taguchi methods) for the development of electronic circuit technology". IEEE Trans. Reliab. 44: 225-229. Teeter, L.D., and Caulfield, J.P. 1991. Stand density management strategies under risk: effects of stochastic prices. Can. J. For. Res. 21: 1373–79. Thompson, E.F., and Haynes, R.W. 1971. Linear programming-probabilistic approach to decision making under uncertainty. For. Sci. 17: 224-229. von Gadow, K. 2000. Evaluating risk in forest planning models. Silva Fenn. 34: 181-191. Wald, A. 1939.Contributions to the theory of statistical estimation and testing hypotheses. Ann. of Math. Stat. 10(4): 299-326. Wald, A. 1945a. Statistical decision functions which minimize the maximum risk. Ann. of Math. 46(2): 265-280. Wald, A. 1945b. Generalization of a theorem by v. Neumann concerning zero sum person games. Ann. of Math. 46(2): 281-286. Wald, A. 1945c. Sequential tests of statistical hypotheses. Ann. of Math. Stat. 16(2): 117-186. Walkup, D.W., and Wets, R.J.B. 1967. Stochastic programs with recourse. SIAM J. Appl. Math. 15: 1299-1314. Weintraub, A., and Abramovich, A. 1995. Analysis of uncertainty of future timber yields in forest management. For. Sci. 41: 217-234. 26 Weintraub, A., Cholaky, A. 1991. A hierarchical approach to forest planning. For. Sci. 37: 439- 460. Weintraub, A., and Vera, J. 1991. A cutting plane approach for chance constrained linear- programs. Oper. Res. 39: 776-785. Weintraub, A., Romero, C., Bjorndal, T., and Epstein, R. 2007. Handbook of Operations Research in Natural Resources. International Series in Operations Research and Management Science, Vol. 99. Springer, New York. Witsenhausen, H.S. 1968. A minimax control problem for sampled linear systems. IEEE Trans. Autom. Control 13(1): 5-21. Zadeh, L.A. 1965. Fuzzy sets. Information and Control 8: 338-353. Zimmermann, H. 1996. Fuzzy set theory and its applications. Third edition ed. Kluwer academic, Boston, USA. Zhou, M. and Buongiorno, J. 2006. Forest landscape management in a stochastic environment, with an application to mixed loblolly pine-hardwood forests. For. Ecol. Manage. 223: 170-182. 27 2 A Robust Optimization Approach to Protect Harvest Scheduling Decisions against Uncertainty2 2.1 Introduction Uncertainty is an important issue in natural resources management since decisions often concern long time horizons and biological processes not totally understood. Although large areas may reduce risk by offering diversification, the decision making process must be carried out in an stochastic environment, where the sources of uncertainty include social, economic, biological, technological and catastrophic factors (Kangas and Kangas 2004, Marshall 1987, Mowrer 2000). ‘Uncertainty’ is traditionally associated with a lack of quantitative information to numerically describe the future and this distinguishes it from ‘risk’, which is referred to a quantified uncertainty (Kangas and Kangas 2004). We use ‘uncertainty’ here simply as a lack of certainty, measurable or not. Decision models have been widely used in forest management planning, but in most of the cases these models assume all information to be deterministic (Kangas and Kangas 1999). This means that decisions implemented now will probably be sub-optimal in the ex-post analysis (Marshall 1987), since all assumptions and initial parameters will not necessarily be observed. In addition, perturbations of supposedly ‘certain’ data can make an optimal solution infeasible (Pickens and Dress 1988). A traditional approach to study the stability of optimal solutions under data perturbations is the post-optimal or sensitivity analysis, which explores the manner in which changes in individual data modify the optimal solution. This can be very useful to address efforts towards a better estimation of key parameters, but its application in real world problems is limited (Pickens and Dress 1988). Another way to consider uncertainty in the planning process has been the use of 2 A version of this chapter has been published. Palma, C.D., and Nelson, J.D. 2009. A robust optimization approach protected harvest scheduling decisions against uncertainty. Canadian Journal of Forest Research 39: 342-355 28 scenario analysis or simulation. Since we know that different conditions (changing market conditions, uncertain yield) or events (fires occurrences or insect attacks) can occur, alternative scenarios of these conditions and events are evaluated. This approach has been extensively used in forestry (Klenner et al. 2000, Peter and Nelson 2005, von Gadow 2000), although its application usually requires a large number of runs, and its usefulness is limited to relatively stable environments (Courtney et al. 1997). This approach has also been used in combination with deterministic optimization techniques (Boyland et al. 2005, Hoganson and Rose 1984, Pukkala 1998), where uncertain coefficients are randomly generated and optimal solutions for a set of scenarios are found. Solutions can be progressively improved and probability distributions of optimal outputs can be obtained. A different group of techniques that explicitly consider uncertainty in the modeling process has also been applied to the harvest scheduling problem. In stochastic programming (Dantzig 1955) uncertainty is considered by scenarios of random events with a given probability of occurrence and a good solution for all of them is found. Some applications of this methodology have dealt with uncertain timber yield (Eriksson 2006, Hoganson and Rose 1987) and stochastic fire losses (Boychuk and Martell 1996, Gassmann 1989). Although this approach provides an intuitive representation of uncertainty, as the number of scenarios increases larger models are produced very quickly and the solution process becomes difficult. Moreover, no explicit representation of the risk-aversion level is considered. In probabilistic or chance constrained programming (Charnes and Cooper 1959) constraints with at least one random coefficient are modeled as probabilistic statements and are required to hold with a minimum probability. This technique has been used to deal with randomness in timber growth (Hof et al. 1996, Pickens et al. 1991, Weintraub and Abramovich 1995, Weintraub and Vera 1991) and uncertain production requirements (Hof and Pickens 1991). As in stochastic programming, probability distributions of uncertain data are assumed to be known which is rarely the case, and finding exact solutions to chance-constrained problems is typically intractable (Chen et al. 2007), although approximation approaches have been provided (see Birge and Louveaux 1997, Chapter 9.4). Finally, fuzzy set theory (Bellman and Zadeh 1970) allows the inclusion of some kind of uncertainties (i.e. vagueness) by defining a grade of membership of a parameter to a certain set of possible values. It has been applied to consider fuzziness in timber yields (Hof et al. 1986, Pickens and Hof 1991), vagueness in age classes boundaries (Boyland et al. 2006) and unclearness of objectives definition in the context of multi-objective problems (Ells et al. 1997, Maness and Farrell 2004). 29 In our view, this approach seems to be more appropriate when there is vagueness, for example, in the meaning of certain events or statements, but not to represent lack of information about the value of the parameters in a model. In this case, the methodology results in a relaxed version of traditional decision problems in which moderate constraint violations are allowed. This will only produce optimal solutions to a new deterministic problem with new ‘relaxed parameters’, and there is no guarantee that these solutions will perform well in different uncertainty scenarios. A different approach to consider uncertainty in optimization models is robust optimization, which is used in this paper and it is described in detail in the next section. As in chance constrained programming, the same idea of increasing the probability of constraint feasibility is used in robust optimization. Based on some assumptions about uncertainty distributions, robust models are able to handle many uncertain parameters and still be computationally tractable. Application of this technique can be found in engineering (Ben-Tal and Nemirovski 2002), network design (Bertsimas and Sim 2003), portfolio management (Bertsimas and Sim 2004) and inventory theory (Bertsimas and Thiele 2006). To our knowledge, this technique has not been applied to natural resources management. In this paper we describe and apply a robust optimization approach to a real size harvest scheduling problem when volume and demand are randomly uncertain. We show in detail how new variables and constraints are added to traditional linear models to produce robust models, and analyze the effect of robust solutions on both harvest decisions and the value of the objective function. The paper is structured as follows. In the second section we describe the robust optimization approach, and the detailed methodology is presented in the third section. Results are presented in the fourth section and then discussed, as well as the limitations of the approach, in the fifth section. Conclusions and future research opportunities are included in the sixth section. 2.2 Robust optimization Robust optimization (RO) is a wide modeling framework used to include uncertainty in optimization models in cases where uncertain data are assumed to belong to an uncertainty set. The term seems to be first used by Mulvey et al. (1995), who used goal programming and a scenario-based description of the problem data to generate less sensitive solutions to the realization of the model data. Previous works by Soyster (1973) and Falk (1976), however, also attempted to find immunized solutions against uncertainty using RO concepts. In these papers 30 uncertain data are assumed to belong to a uniform symmetric range, and linear problems with a modified feasible region are solved. The approach proposed by these authors produces solutions to the worst case scenario that are, as one of the author states, ‘ultraconservative strategies’ (Soyster 1973). Less conservative approaches have been suggested (Ben-Tal and Nemirovski 2000, El Ghaoui et al. 1998), but the resulting robust models become harder to solve (nonlinear) than traditional formulations. For a recent review of RO the reader is referred to Beyer and Sendhoff (2007) and Ben-Tal and Nemirovski (2008). A less conservative approach is also proposed in Bertsimas and Sim (2004), and is used in this paper. A remarkable feature of this approach is that it transforms the original linear model into another robust linear, and therefore computationally tractable, model. We next describe this approach in general and then apply it to a harvesting scheduling problem. 2.2.1 Basic concepts We first define basic concepts used throughout the paper. The protection function of a constraint is a function that determines the size of the buffer required to immunize a constraint against uncertainty. It is defined for each constraint where uncertain coefficients are present and for which feasibility is desired. As we will see later, this function is an optimization problem in itself that through a dual transformation results in additional variables and constraints to the original problem. The magnitude of the buffer for each constraint will depend on the degree of uncertainty of the parameters involved in the constraint and the degree of protection desired, which we will refer as protection level. By changing this protection level, we obtain different degrees of conservatism of the solutions, that is, different expected rates of constraints violation. However, the protection level is not equivalent to the desired probability of constraints violation. Although these two concepts are directly related (the higher the protection level, the lower the probability of constraints violation), the protection level is just a parameter with no direct meaning about the resulting probability of constraints violation. For each set of constraints for which protection is desired, the general steps required to formulate a robust model are: (S1) Formulate the protection function, (S2) Transform the protection function into an optimization model, 31 (S3) Obtain the dual model of the linear version of the protection function, and (S4) Embed the dual model in the original model. The resulting model is called the robust counterpart of the deterministic model and despite having more variables and constraints (as result of the transformations needed), it is linear and can be solved using commercial linear optimization packages. New parameters like the degree of uncertainty of the data and user-defined protection levels are also required. The above steps are described in more details below and then applied to a harvest scheduling model. The methodology is based on the following assumptions. Uncertain coefficients are to be uniform and independently distributed within a symmetrical range of values. These assumptions allow the methodology to rely on the fact that only some of the uncertain coefficients will change (from the coefficient estimate) to adversely affect the solution. We will come back on the implication of these assumptions in the Discussion section. 2.2.2 Robust formulation The following description of the RO approach is mainly based on Bertsimas and Sim (2004), with slight changes introduced. As decision variables in forest planning models are almost always positive (usually the area to be harvested) we will consider only non-negative variables in our explanation. Although the original formulation produces equivalent models, the following formulation leads to both less variables and constraints. Let us consider the general linear problem: ∑ =j jj xcMax 1 ]1.2[ n n jx ∀≥ 0]3.2[ subject to ibxa j ijij ∀≤∑ =1 ]2.2[ j 32 We will consider that only aij are uncertain, but the case in which cj and bi are also uncertain can be derived by expressing the objective function as a constraint and moving bi to the left-hand side of [2.2]. Let us assume that belongs to the symmetric range ija ]ˆ,ˆ aaaa +−[ , with ijijijij ija âthe coefficient estimate (usually used in deterministic formulations) and the precision of the estimate. Associated with the uncertain data a , a random variable of scaled deviation that takes values in [-1,1] can be defined as ij ij ijijijij aaa ˆ)( −=η . If all coefficients in constraint i are uncertain, then can theoretically take values between -n and n. If only a set J∑ =j ij1ηn [ ] i of coefficients are uncertain in constraint i, this interval will be between –card(Ji) and card(Ji), where card represents the cardinality or number of elements in a set. Having this in mind, a new parameter can be used in the mathematical expression )(,0 ii Jcard∈Γ iJj iji Γ≤∑ ∈ η to control the protection level of constraint i against uncertainty as follows: (i) if Γ , the 0=i ijη for all j in constraint i are forced to 0, so that ijij aa = for all j, and constraint i has no protection against uncertainty, (ii) if Γ constraint i is totally protected against uncertainty since all )( ii Jcard= ijη in constraint i are allowed to take a non-zero value, and (iii) if )(0 Jcard ii <Γ< constraint i is partially protected against uncertainty, in which case some of the ijη can be different from zero. 2.2.2.1 Formulating the protection function (S1) To obtain the robust counterpart of model [2.1]-[2.3], we first split the left-hand side of the constraints [2.2] into the contribution of the coefficient estimate ija and the noise produced by the precision of the estimate, . The constraint set is then reformulated as: ijâ ibxa j iiijij ∀≤Γ+∑ =1 * ),x(]4.2[ βn where ⎣ ⎦ ⎣ ⎦( ) ⎭⎬ ⎫ ⎩⎨ ⎧ Γ−Γ+=Γ ∑ ∈Γ=∈⊆ i i iiiiiiiii Sj jitiijijScardSJtJStSii xaxa ** )(,\,|}{ * ˆˆmax),x(]5.2[ U β 33 is the protection function of constraint i. This function is specific for a given solution vector x* and depends on the user-defined protection level iΓ . Note that when no protection is required , and the protection function β0=Γi i equals 0 and the constraint is the same as in the deterministic case. The protection function [2.5] looks for the appropriate buffer to keep the left-hand side of [2.2] lower than bi for different values of aij. As we will see later in our example, if the inequality is ≥ then the protection function has to be subtracted in the original constraint to keep the left-hand side term greater than bi. Since some uncertain parameters will finally take values above and below the estimate ija Γ , it is unlikely that the scenario with worst-case values occurs. Based on this intuitive idea, the protection function considers only a subset Si ⊆ Ji of uncertain coefficients, with the number of elements in Si given by the integer part of the user-defined level of protection, , that is, the higher the protection level, the more coefficients are considered in the protection function. The constraint i is then protected against all cases in which up to ⎣ ⎦i ⎣ ⎦iΓ of the coefficients change and one extra coefficient (ti in the protection function iβ ) changes by ( - ) a . Even if more than iΓ Γ ˆ⎣ ⎦i it ⎣ ⎦iΓ coefficients change, this approach provides a solution that will be feasible with high probability (Bertsimas and Sim 2004). 2.2.2.2 Transforming the protection function into an optimization model (S2) To find the value of the protection function of each constraint i, the value of [2.5] can be replaced by the objective function of the following linear optimization problem. Note that x*j is not a variable in this problem. ij Jj jij wxaMax i ∑ ∈ ˆ]6.2[ * Jjw ∈ subject to i Jj ij i w Γ≤∑ ∈ ]7.2[ iij ∀≤≤ 10]8.2[ 34 where wij are new variables that represent the random variable ijη described above. As proof of the equivalence of these problems, notice that the optimal solution value of problem [2.6]-[2.8] clearly consists of variables at 1 and one more variable at ⎣ ⎦iΓ iΓ - ⎣ ⎦iΓ . This is equivalent to select a subset {Si U {ti} | Si ⊆ Ji, |Si| = ⎣Гi⎦, ti ∈ Ji\Si} with the corresponding objective value (Bertsimas and Sim 2004). ⎣ ⎦∑ ∈ Γ−Γ+i iSj jitiijij xaxa ˆˆ ( ) ** * Jjp ∈∀≥ 0]11.2[ 0]12.2[ ≥z n 2.2.2.3 Obtaining the dual model of the linear version of the protection function (S3) Although model [2.6]-[2.8] is linear for the given solution vector x*, it is not in the general case where x is variable. However, by using its dual, equation [2.4] can be linearly expressed. Since problem [2.6]-[2.8] is feasible and bounded for all Гi, then its dual, shown below, is also feasible and bounded. ∑ ∈ +Γ iJj ijii pzMin]9.2[ subject to ijijiji Jjixapz ∈∀≥+ ,ˆ]10.2[ iij i where zi and pij are the dual variables associated with constraints [2.7] and [2.8], respectively. 2.2.2.4 Embedding the dual model in the original model (S4) We finally replace the protection function of [2.4] with the objective function of the dual model. The robust counterpart of [2.1]-[2.3] is then the following linear model: ∑ =j jj xcMax 1 ]13.2[ subject to 35 ibpzxa j i Jij ijiijij ∀≤+Γ+∑ ∑ = ∈1 ]14.2[ n Jjixapz ∈∀≥+ ,ˆ]15.2[ jx ∀≥ 0]16.2[ Jjip ∈∀≥ ,0]17.2[ iz ∀≥ 0]18.2[ ijijiji j iij i Model [2.13]-[2.18] is larger than the original deterministic formulation, but it is of the same type (linear) and is computationally tractable. In the ‘Controlling the protection level’ section we explain how the protection level can be determined. 2.3 Methods We applied the above approach to the well known harvest scheduling problem to allow readers to focus on the application of RO. The same steps described above can be applied to other problems to produce robust models. We considered uncertainty in volumes and demands of two products over the entire planning horizon. We applied deterministic and robust models to a timber-oriented forest area, and compared the performance of these decisions through the simulation of the uncertain data. We conducted two main analyses. The first (A1) demonstrates the advantage of using this methodology over the strategy of manually imposing a buffer to the original demand levels. For simplicity, we focused only in protecting minimum demand constraints when volume was uncertain. We used the buffer determined by the robust model for each constraint to correct the minimum demand level of the same constraint in the deterministic model. We then showed the difference in the decisions’ stability between using a (more complex) robust model and a traditional model with a corrected demand level that acts as a security buffer. In the second analysis (A2) we considered a complete model with protection for minimum demand and fluctuation constraints. We explored infeasibility rates and the change in the objective function and harvest decisions when volume and demand were uncertain under two discount rates. 36 All the optimization models described below were implemented in ILOG OPL Development Studio 5.2 (using CPLEX 10.2.0 as optimizer) while simulation experiments were conducted using Visual Studio 2005. All models were run in a 3.2 Ghz Pentium D with 1GB of RAM. Details of the different components of our methodology are presented next. 2.3.1 Study area The models were applied to 245,090 ha of harvestable forest corresponding to the Tree Farm License 48, located in the north-east area of British Columbia, Canada. The forest is mainly mature (Fig. 2.1), and the predominant species are white spruce (Picea glauca), lodgepole pine (Pinus contorta), subalpine fir (Abies lasiocarpa), trembling aspen (Populus tremuloides) and cottonwood (Populus sp.) (Schuetz 2002). Fig. 2.1. Initial age class distribution shows a mature forest with almost 80% of the area in stands at least 40 years-old. The area was divided into 800 strata with an average area of 306.4 ha and ranging from 92.1 to 4,255.5 ha. Strata were defined as those stands with similar age class, species composition and harvest system, regardless of their spatial location. For each stratum the following information was available: (1) volume of two products [m3] over the planning horizon for different management options, (2) management cost [$/ha] and (3) harvesting cost [$/m3]. Since not much 37 detail is required for long-term decisions, we considered conifers and deciduous species as two different products. 2.3.2 Deterministic model formulation We opted for a model I formulation (Johnson and Scheurman 1977) of the harvest scheduling problem as the base for our robust model. Since the robust formulation significantly increases the number of constraints of the original model we preferred model I as it produces fewer constraints than model II. Thus, our decision variable xij represents the number of ha of strata i to which management j is applied, where each management defines a sequence of harvests throughout the planning horizon. Based on their initial age and growing rates, we defined fifteen management options for each of the 800 strata in order to represent a range of potential rotation ages. A twenty 10-year period planning horizon was used and harvests were assumed to occur in the middle of the period. Only clearcuting was considered as management option. The objective is to maximize the net present value of the harvest plus the value of the ending inventory: ∑∑ + ji ijijT t tt xlzMax , ]19.2[ αα where α represents the discount factor, lij the value of the ending inventory if management j is applied to stratum i, and zt the net income in period t, computed as follows: ∑∑ −−= ji ijijt pji ijijptiptt xcmxvchbz ,,, )(]20.2[ where bpt is the price ($/m3) of the product p in t, vijpt is the volume (m3) in period t of product p in stratum i if managed with j, and chi and cmijt are the harvesting ($/m3) and management costs ($/ha), respectively. The following sets of constraints – area, minimum demand and flow constraints – were considered: iax j iij ∀=∑]21.2[ 38 tpdxv ptij ji ijpt ,]22.2[ , ∀≥∑ 1,)1(]23.2[ , 1 , >∀−≥ ∑∑ − tpxvxv ij ji ijptij ji ijpt δ where is the area of stratum i (ha), dia pt is the minimum requirement of product p in period t (m3), and δ is the allowable reduction in harvest volume as a proportion of the production in the previous period. We allowed 10% reduction, that is, δ = 0.1. The minimum demand was assumed the same for all periods and determined on the basis of the production of the non-restricted model, as explained in the Results section. We considered two discount rates. A low 2% discount rate that includes more information about the future into the analysis (Gassmann 1989) and the more commonly used 4%. The discount factor was therefore defined by αt = (1+r)-(t-0.5) where (t-0.5) represents the midpoint of the period. In addition, we used local timber prices for these products (B.C. Ministry of Forests and Range 2007). Although some of these constraints may not be relevant for all decisions makers (e.g. private forest owners may not be interested in satisfying minimum demand requirements), they are important to many other forest managers and allow us to show how RO can be applied to different kind of constraints. 2.3.3 Uncertain data As mentioned previously, we assumed that volume and demands in the above model were uncertain. Since independency is one of the assumptions for uncertain coefficients, we addressed the uncertainty coming from unbiased errors existing in yield forecast models. As strata are made up of stands from different locations in the forest, we also rely on the assumption that volume is not spatially correlated. Although demand can be correlated, we assumed it is independently distributed. The consequences of these assumptions are discussed in the Discussion section. In regard to the accuracy of volume estimates over long-term, studies are scarce (Yaussy 2000) and have suggested dissimilar values. Estimate accuracy of total volume can range from less than 0.5% in plantations (Trincado et al. 2003) to 38% in natural forests (Pretzsch et al. 2002) for 5- year forecasts. For 30-year periods, Yaussy (2000) reports forecast errors in total volume ranging from 3% to 98% depending on the model and the site. Moreover, it is reasonable to think that estimate errors increase for longer periods of estimation (Trincado et al. 2003). Hence, we used 39 an increasing volume estimate error, starting with a 10% accuracy for the first planning period, and an additional 2% for each future period. This implies that the estimate error was age- independent and it was only related to how far the estimation was done from the first period. For the last planning period, for instance, the volume estimates had a 48% error. Estimate error of demands followed the same shape, with 1% accuracy for the first period. In all cases, uncertain coefficients were assumed independent and uniformly distributed within the range defined for their estimate precision. 2.3.4 Robust model formulation New variables and constraints need to be defined to build the robust counterpart of the previous model. We applied the robust optimization approach previously described in order to get protection against infeasibility in constraints [2.22] and [2.23]. Detailed steps of how this formulation was obtained are presented in the Appendix A. Following the notation used, v and d represent the coefficient estimate of volume and demand, and and the respective precision ranges. v̂ d̂ 2.3.4.1 Protection against minimum demand infeasibility We have to consider two cases when protection for demand constraints is formulated: (a) when only volume is uncertain, and (b) when both volume and demands are uncertain. Since an extra coefficient (demand) which is not associated with harvest decisions is uncertain in (b), the new equations differ slightly for each case. For case (a), constraint [2.22] was replaced with the following set of constraints: tpdwuxv pt Jji ijptptpt ji ijijpt D pt ,]24.2[ ),(, ∑∑ ∈ D ∀≥−Γ− D tpu ,0]26.2[ ∀≥ D tpJjixvwu ptijijptijptpt ,,),(ˆ]25.2[ ∈∀≥+ pt tpJjiw ptijpt ,,),(0]27.2[ ∈∀≥ 40 where upt and wijpt are new variables that originated in the dual problem, and JDpt is the set of uncertain coefficients of the demand constraint pt. For case (b) constraints [2.25]-[2.27] are still valid, but constraint [2.24] was replaced with the following constraint: tpdwudxv pt Jji ijptptptpt ji ijijpt D pt ,)1(]28.2[ ),(, ∀≥−−Γ−− Dˆ ∑∑ ∈ See Appendix A for more details. 2.3.4.2 Protection against infeasibility for fluctuation constraints To get protection for timber flow between consecutive periods, constraint [2.23] was replaced with the following set of constraints: 1,)1(]29.2[ , 1 ),(, >∀−≥−Γ− ∑∑∑ − ∈ tpxvzyxv ij ji ijpt Jji ijptptptij ji ijpt F pt δF ( ) 1,,),(ˆ)1(ˆ]30.2[ )1( >∈∀−+≥+ − tpJjixvvzy ptijtijpijptijptpt δ F 1,0]31.2[ >∀≥ tpy F pt 1,,),(0]32.2[ >∈∀≥ tpJjiz ptijpt where ypt and zijpt come from the dual problem, and JFpt is the set of uncertain coefficients of the fluctuation constraint pt. 2.3.4.3 Controlling the protection level The degree of conservatism in satisfying a constraint is controlled by the level of protection Гi. This value refers to the number of uncertain coefficients allowed to change while the constraint remains feasible. The good news is that this number can be associated with probability bounds (Bertsimas and Sim 2004) to determine the value of Гi to meet the constraint i with a given probability. The bad news is that even the best of these bounds is quite conservative and represents a weak estimate of the probability of meeting the constraint. We tested (results are not presented) bounds proposed by Bertsimas and Sim (2004) and determined Гi for different 41 probabilities of constraint violation. However, when we allowed 40% of constraint violation, 100% of solutions were feasible when uncertain coefficients were simulated. Only when 45% of constraint violations were allowed, less than 0.5% of simulations were infeasible. In other words, the use of current probability bounds overestimates the level of protection Гi required to hold a constraint within a given probability. As a result, objective values are excessively affected and the real tradeoff between conservatism and solution quality is hidden. We therefore opted to gradually increase the protection levels, Γi, without using any probability bounds to determine constraints violation. Instead we tested infeasibility rates by simulating the uncertain coefficients as described later. Since for a constraint violation of at most εi we need a protection level at least equal to )()1(1 ii Jcardε−Φ+ 1− , where Φ is the cumulative distribution function of a standard normal (Bertsimas and Sim 2004), we estimated a maximum protection level (MPLi) by considering εi = 0.01. That is, for each constraint i, it is sufficient to select Γi at least equal to MPLi to get a 99% protection. We then simply determined different protection levels for constraint i as a percentage, e.g. 0%, 10%, 20% and 30%, of the MPLi value. Disadvantages of this approach are discussed in the Discussion section. 2.3.5 Simulation experiments Both deterministic and robust solutions were tested by simulating the uncertain coefficients of the model within the uncertainty sets previously described. That is, we evaluated the optimal decisions from the deterministic and robust models using randomly generated timber yields and demands. For A1 we performed 1,000 simulations of the volume coefficients and infeasibility rates and changes in the objective function were examined for different protection levels. For A2 two scenarios were considered; the first one dealt with uncertainty in timber volume and the second one with uncertainty in both volume and demand. For each scenario we ran 1,000 simulations and evaluated the performance of deterministic and robust decisions with different protection levels. The resulting objective function, change in harvest decisions and occurrence of infeasibility were examined. 42 2.4 Results For all analyses we determined the minimum demand requirement for our two base cases, with 2% and 4% discount rate (Fig. 2.2b), to be around 85% of the average production per period by running the non-restricted deterministic models (Fig. 2.2a). That is, demands were defined as 3,600,000 and 750,000 m3 per period for product 1 and 2, respectively. For product 2, however, this constraint was considered from period 2 onward, as the problem became infeasible otherwise because of the initial forest conditions. The same demand levels were used in all scenarios evaluated. Minimum demand as well as fluctuation constraints became active in many of the planning periods in our base cases (Fig. 2.2b). As expected, 4% discount rate moved production forward in time both for the unrestricted and the restricted models. 43 Fig. 2.2. The non-restricted timber flow (a) was constrained by minimum demand and maximum negative fluctuation between period constraints (b) to define the base case for the analysis. In analysis A1 protection for demand constraints and uncertain volume with a 2% discount rate were considered. Recall that the buffer (value of the protection function) determined by the robust model for each product and period was recorded and then used to modify the minimum demand of the deterministic model. Thus, robust and deterministic models only differed in the additional constraints and variables of the robust model. With the same buffer levels in both models, robust decisions (ROB) reported less infeasibility rates (a solution was considered infeasible if at least one period was infeasible for the corresponding constraint) than deterministic decisions with the buffer added to the original demand (BUF) (Fig. 2.3). Although 44 ROB reported a higher reduction in the objective function than BUF as protection level increased, this reduction was proportionally lower for the same infeasibility rate. For example, at 20% protection level, the same buffer size produced less than 10% of infeasibilities in ROB and more than 30% in BUF. On the other hand, for a 20% infeasibility rate, ROB reduced the objective value in about 0.37% (protection level 17%) while the reduction with BUF was around 0.43% (protection level 27%). Fig. 2.3. Lower average infeasibility rates were obtained with the robust model (ROB) when the same buffer was used to modify the demand of the deterministic model (BUF). Although the objective value reduction is greater in the robust model for different protection levels, this reduction is lower when the same level of infeasibility rates is considered. In A2 we applied the full robust model with protection against infeasibility in both demand and fluctuation constraints. For both discount rates, the robust model increased production levels in critical periods – where demand constraints were active in the base case – as more protection against uncertainty was required (Fig. 2.4). When only volume was uncertain (Fig. 2.4a and b) this increase was slight, but including uncertain demands imposed higher levels of production in order to ensure constraint feasibility (Fig. 2.4c and d). Smoother production was also obtained to avoid large fluctuations between consecutive periods. 45 Fig. 2.4. Increasing timber production in periods with active minimum demand constraints were obtained when more protection was required. This increase was clearly larger when volume and demands were uncertain (c and d) than when only volume was uncertain (a and b). The occurrence of infeasibilities was reduced with robust solutions, both in the amount of simulations with at least one infeasible period (Fig. 2.5) and in the average number of infeasible periods in each simulation (Fig. 2.6). The constraints in the deterministic solution (protection level 0%) were almost always infeasible in at least one period, both in the scenario with uncertain volume (Fig. 2.5a and b) and in the scenario with uncertain volume and demand (Fig. 2.5c and d). Demand of product 2 showed infeasibility rates between 70 and 80% only for scenarios with a 4% discount rate. In all cases, the rate of infeasibility was significantly reduced as higher protection levels were considered. For the highest protection level we tested (30%), infeasibilities were virtually nill for demand constraints and a bit higher for volume fluctuation constraints. The objective function was clearly affected when higher protection levels were imposed. This reduction was higher when a higher discount rate was used (Fig. 2.5b and d), both 46 when volume (Fig. 2.5b) and volume and demand were uncertain (Fig. 2.5d). The reduction in the average number of infeasible periods was even more important. Results were similar for both discount rates, so for the sake of space we only showed results for 4% discount rate (Fig. 2.6). For example, an average of around 7 infeasible periods observed by simulating the deterministic solution was reduced to less than 0.19 and 0.11 periods for the scenarios with uncertain volume (Fig. 2.6a) and volume and demands (Fig. 2.6b), respectively. Fig. 2.5. Reductions in the average percentage of infeasible simulations were observed as higher protection levels were used when volume (a and b) and volume and demands (c and d) were assumed uncertain. The impact in the objective function was higher for higher discount rates (b and d). 47 Fig. 2.6. Reductions in the average number of infeasible periods were also observed as higher protection levels were used when volume (a) and volume and demands (b) were uncertain. Similar results were obtained for a 2% discount rate (not shown). Harvest decisions were increasingly different as more protection level was considered (Table 2.1). At 30% protection level, more than 155,000 ha were harvested in a different planning period during the planning horizon when volume was assumed uncertain and more than 225,000 ha when both volume and demands were. Also at 30% protection and 4% discount rate, first- period harvest decisions, which are most likely to be implemented, changed in about 2,200 and 3,400 ha when volume and volume and demands were uncertain, respectively. These changes nearly doubled in the 2% discount rate scenario. Table 2.1. Differences in harvest decisions for robust models with different protection levels Change in harvest decisions [ha] (r = 2% / r = 4%) Uncertain volume Uncertain volume & demand Protection level [%]* First Period Total First Period Total 10 2,226 / 1,645 87,494 / 95,008 2,736 / 965 160,818 / 200,080 20 3,813 / 1,666 128,138 / 123,447 4,547 / 2,172 197,341 / 235,308 30 4,030 / 2,207 158,939 / 155,897 6,099 / 3,417 225,397 / 258,426 * Percentage of the Maximum Protection Level in a constraint As expected, the dimension and solution time of the models increased when uncertain coefficients were introduced (Table 2.2). The size of the robust models with uncertainty in 48 volume and in both volume and demands is the same since neither new variables nor constraints are required when demand becomes uncertain. Table 2.2. Size and solution time for the problems (reported by the software) Model Variables Constraints Solution time* [sec] Deterministic 12,001 879 1-2 Uncert. volume 95,355 84,154 103-404 Uncert. vol/dem 95,355 84,154 239-492 * A range of solution time is provided since models were run for two discount rates and with different protection levels. In general, the higher the protection level the more difficult to solve the problem (higher solution time) 2.5 Discussion Our results suggest that the robust optimization approach produces harvest decisions that are less sensitive to uncertain data and therefore protected from the occurrence of infeasibilities, with a moderate reduction in the objective function. In other words, for different degrees of conservatism, robust decisions performed almost as well as deterministic optimal ones, but for a range of possible volume and demand values. As expected, in order to get protection against infeasibility, production levels increased in critical periods in which demand constraints were tightly met. In periods with high negative fluctuation, production was smoothed by moving it from periods with a volume surplus to periods where only the necessary volume to meet constraints was produced. Although this ‘buffer’ strategy to protect decisions from an uncertain future is known (Boychuk and Martell 1996), the robust approach contributes a way to estimate this buffer in conjunction with the decision variables and based on the degree of uncertainty of the data. We showed that even in the case that the optimal buffer obtained by robust models is used to modify the constraint level of a deterministic model, solutions are not necessarily robust and higher infeasibility rates are obtained. Unlike deterministic models that look for optimal solutions, robust models look for optimal ‘stable’ solutions. When demand was uncertain, the need for markedly higher production levels was evident. Since demand coefficients are always present regardless of the value of the decision variables, the model needs to consider the worst case scenario of demand in order to get protection. This is also 49 clear from constraint [2.28], where deterministic demand requirements are corrected by their corresponding error estimates. Because of this, the reduction in the number of infeasible periods became more important when demand uncertainty was added to the problem (Fig. 2.6). As it has been shown in other studies (Pickens and Dress 1988), deterministic solutions were highly infeasible when tested by simulating the uncertain data. Robust decisions, however, proved to be more resistant to random changes in data and thus reduced infeasibility rates to almost zero for both discount rates. Since these decisions are required to be good (not optimal) for multiple possible future values of uncertain parameters, their objective values will certainly be negatively affected, and this represents the cost that decision makers should be willing to assume if a reduction in the risk associated with their decisions is desired. This reduction was greater for the higher discount rate (4%) because robust decisions moved harvests from initial (more valuable) to later (less valuable) periods to satisfy the increasing buffer requirements. The impact on the objective function is likely specific to the initial age class distribution, and we expect it would differ in different forests. We note here that a number of environmental requirements (i.e. spatial constraints) are being imposed on forest plans and these also reduce value as our approach does. We cannot say, however, that the reduction in the objective function produced by robust models should be added to the reduction produced by environmental requirements. In this case, an ‘optimal’ solution would probably exploit the interaction in decisions and therefore produce a value reduction lower than the simple addition of the independent reductions. Although we used the same protection levels for all constraints in each of our analyses, they can be different if we want to express different importance to each constraint. For instance, if the occurrence of infeasibility of fluctuation constraints is not a critical issue or the feasibility of one of the products does not need to be guaranteed, then lower protection levels can be used for these constraints. This will certainly produce better solutions in terms of the value of the objective function. Whether or not to take uncertainty into account in the decision process is irrelevant for management purposes if decisions remain the same in both cases. Future harvest levels are considered as reference for the current decision making process in the context of a rolling planning horizon, so they have less effect on immediate management decisions. Changes in the first-period harvest decisions, however, are relevant as they are most likely to be implemented. 50 Depending on the scenario, between 2,200 and 6,100 ha were harvested differently in the first period, that is, the harvest of some stands was moved forward or backward to meet the protection requirements. As expected, for the most uncertain scenario (uncertain volume and demand), these changes in decisions were highest. For the 2% discount rate, changes were also more important as the future and the degree of the uncertainty becomes more relevant when low discount rates are used. We gained no clear insight about what kind of stands either shortened or lengthened their rotation ages, so no general harvest policy insight could be obtained. We suspect that the magnitude of the decision changes is scenario-specific and largely influenced by the interaction between the forest state and the constraints imposed on timber production. Only a slight increase in the average rotation age (about half a planning period) was observed in different scenarios, maybe due to the need for more volume all along the planning horizon and the maturity of the initial forest (Fig. 2.1). This produced a timber shortage in the second part of the planning horizon that forced longer rotation ages in order to keep the harvest levels. Robust models are larger than traditional models. If a set of i constraints with j variables each (all multiplied for an uncertain coefficient) are to be protected, the robust model will have i+i*j constraints (not considering non-negativity) and variables. Although this can significantly enlarge the size of the models, like in our case, they remain linear. This is in our view an important advantage of this approach. Other methodologies (e.g. stochastic) would produce enormous linear problems with only a fraction (discrete scenarios) of the uncertain data we considered, or non-linear models much harder to solve (e.g. chance constrained programming). The methodology basically works as a buffering strategy. An additional term, the protection function, is added to critical constraints for which feasibility is highly desirable and uncertainty is present in the technical coefficients. This term is an optimization problem in itself. It relates decision variables to the buffer needed to meet the constraint for a given protection level defined by the modeler. Embedded in the original optimization problem, the protection function identifies decisions that provide the appropriate buffer, but at the same time, contribute to the general objective function. We believe this approach would allow decision makers interested in reducing the risk of their decisions to objectively determine the amount of buffer required to guarantee critical requirements. Without the need of solving complex problems and using sophisticated solution techniques, modified versions of traditional deterministic problems (following the steps described here) can be built and solved by a commercial linear optimization software, and then solutions tested by simulating uncertain coefficients. 51 However, our approach includes a number of simplifying assumptions that should be addressed in future research. The most important one refers to the assumption of independency of the uncertain coefficients, which limits its application to random, unbiased and uncorrelated errors in the coefficients estimates. Although this kind of error is very important when forecast models are used, catastrophic events and uncertain trends in data cannot be properly modeled. If this assumption is not met, optimistic results will be obtained. Another drawback of this methodology is the lack of a clear relationship between protection levels and the resulting infeasibility rates. Although probability bounds exist that relate constraints violation rates to protection levels, they are loose and produce more conservative solutions than desired. As explained in the Methods section, we did not use probability bounds but selected protection levels as different percentages of the square root of the number of uncertain coefficients in a constraint. Although the concept has a clear meaning, the amount of protection level to be used is not intuitive, and makes difficult for the decision maker to a priori control the conservatism of the solutions. This leads to the need for simulation experiments to test the solution performance. Further research looking for a more accurate way to relate the probability of constraints violation and the protection level would certainly improve this methodology. If probability bounds are used (as opposed to simulation experiments as used here) it would also be of interest to compare this approach to chance constrained programming. We did not address this comparison here since current probability bounds for robust optimization are quite weak, which leads us to clearly envisage more conservative solutions than a chance constrained approach. We plan, however, to do a structured comparative analysis in a future work. On the other hand, the assumption of uniform and symmetrically distributed uncertainty can be considered a limitation when the nature of the uncertainties is well described by a different distribution, or extra information other than extreme points is available. We suggest this can be handled by considering a wide enough uniform range for random values that embraces, for instance, 99% of the observations of the original distribution. This will certainly produce more conservative solutions than if a true distribution is used, but it will retain the modeling simplicity required to produce robust linear models. The use of extra information other than extreme points of uncertainty sets (e.g. asymmetry) should also be explored. Finally, the possibility to include uncertainty in coefficients arisen from the product of two uncertain coefficients would be useful. This would allow considering uncertainty in timber price and volume, and therefore get protection against the value of the objective function. 52 2.6 Conclusions We presented a harvest scheduling robust model that considered uncertainty in volume and demands and produced immunized decisions to uncertain data. By increasing protection levels against constraint infeasibility, the tradeoff between conservatism and cost of robustness was explored, showing that robust decisions had a modest effect on the value of the objective function. We also showed that this approach produces stable solutions that are not found when constraint levels are modified with a buffer amount to protect constraints from infeasibility. The robust optimization approach increased the number of variables and constraints of the traditional deterministic harvest scheduling model, but still produced a linear model that could be solved in reasonable solution time by using commercial optimization software. Although simulation experiments are required to finally estimate the expected rate of constraint infeasibility, we consider the methodology is an attractive way to find stable decisions in the face of uncertainty while keeping optimization models computationally tractable. 53 2.7 References B.C. Ministry of forests and range. 2007. B.C. Interior log market. Report for the 3 month period March 1, 2007 to May 31, 2007. Ministry of Forests and Range. Revenue Branch, Canada. 1p. http://www.for.gov.bc.ca/hva/timberp/llsp/interior/historical.htm (accessed Jan 20, 2008) Bellman, R.E., and Zadeh, L.A. 1970. Decision-making in a fuzzy environment. Management Science Series B-Application 17: B141-B164. Ben-Tal, A., and Nemirovski, A. 2008. Selected topics in robust convex optimization. Math. Program. 112: 125-158. Ben-Tal, A., and Nemirovski, A. 2002. Robust optimization - methodology and applications. Math. Program. 92: 453-480. Ben-Tal, A., and Nemirovski, A. 2000. Robust solutions of linear programming problems contaminated with uncertain data. Math. Program. 88: 411-424. Bertsimas, D., and Thiele, A. 2006. A robust optimization approach to inventory theory. Oper. Res. 54: 150-168. Bertsimas, D., and Sim, M. 2004. The price of robustness. Oper. Res. 52: 35-53. Bertsimas, D., and Sim, M. 2003. Robust discrete optimization and network flows. Math. Program. 98: 49-71. Beyer, H.G., and Sendhoff, B. 2007. Robust optimization - a comprehensive survey. Comput. Methods Appl. Mech. Eng. 196: 3190-3218. Birge, J.R., and Louveaux, F. 1997. Introduction to stochastic programming. Springer, New York. Boychuk, D., and Martell, D.L. 1996. A multistage stochastic programming model for sustainable forest-level timber supply under risk of fire. For. Sci. 42: 10-26. Boyland, M., Nelson, J., and Bunnell, F.L. 2005. A test for robustness in harvest scheduling models. For. Ecol. Manage. 207: 121-132. Boyland, M., Nelson, J., Bunnell, F.L., and D'Eon, R.G. 2006. An application of fuzzy set theory for seral-class constraints in forest planning models. For. Ecol. Manage. 223: 395-402. Charnes, A., and Cooper, W.W. 1959. Chance-constrained programming. Management Science 6: 73-79. Chen, X., Sim, M., and Sun, P. 2007. A robust optimization perspective on stochastic programming. Oper. Res. 55: 1058-1071. 54 Courtney, H., Kirkland, J., and Viguerie, P. 1997. Strategy under uncertainty. Harv. Bus. Rev. 75: 67-79. Dantzig, G.B. 1955. Linear programming under uncertainty. Management Science 1: 197-206. El Ghaoui, L., Oustry, F., and Lebret, H. 1998. Robust solutions to uncertain semidefinite programs. Siam Journal on Optimization 9: 33-52. Ells, A., Bulte, E., and vanKooten, G.C. 1997. Uncertainty and forest land use allocation in British Columbia: Vague priorities and imprecise coefficients. For. Sci. 43: 509-520. Eriksson, L.O. 2006. Planning under uncertainty at the forest level: A systems approach. Scand. J. For. Res. 21: 111-117. Falk, J.E. 1976. Exact solutions of inexact linear programs. Oper. Res. 24: 783-787. Gassmann, H.I. 1989. Optimal harvest of a forest in the presence of uncertainty. Canadian Journal of Forest Research 19: 1267-1274. Hof, J., Bevers, M., and Pickens, J. 1996. Chance-constrained optimization with spatially autocorrelated forest yields. For. Sci. 42: 118-123. Hof, J.G., and Pickens, J.B. 1991. Chance-constrained and chance-maximizing mathematical programs in renewable resource-management. For. Sci. 37: 308-325. Hof, J.G., Pickens, J.B., and Bartlett, E.T. 1986. A maxmin approach to nondeclining yield timber harvest scheduling problems. For. Sci. 32: 653-666. Hoganson, H.M., and Rose, D.W. 1987. A model for recognizing forestwide risk in timber management scheduling. For. Sci. 33: 268-282. Hoganson, H.M., and Rose, D.W. 1984. A simulation approach for optimal timber management Scheduling. For. Sci. 30: 220-238. Johnson, K.N., and Scheurman, H.L. 1977. Techniques for prescribing optimal timber harvest and investment under different objectives – discussion and synthesis. For. Sci. 1-31. Kangas, A.S., and Kangas, J. 2004. Probability, possibility and evidence: approaches to consider risk and uncertainty in forestry decision analysis. Forest Policy and Economics 6: 169- 188. Kangas, A.S., and Kangas, J. 1999. Optimization bias in forest management planning solutions due to errors in forest variables. Silva Fenn. 33: 303-315. Klenner, W., Kurz, W., and Beukema, S. 2000. Habitat patterns in forested landscapes: management practices and the uncertainty associated with natural disturbances. Comput. Electron. Agric. 27: 243-262. Maness, T., and Farrell, R. 2004. A multi-objective scenario evaluation model for sustainable forest management using criteria and indicators. Canadian Journal of Forest Research 34: 2004-2017. 55 Marshall, P.L. 1987. Sources of uncertainty in timber supply projections. For. Chron. 63: 165- 168. Mowrer, H.T. 2000. Uncertainty in natural resource decision support systems: sources, interpretation, and importance. Comput. Electron. Agric. 27: 139-154. Mulvey, J.M., Vanderbei, R.J., and Zenios, S.A. 1995. Robust optimization of large-scale systems. Oper. Res. 43: 264-281. Peter, B., and Nelson, J. 2005. Estimating harvest schedules and profitability under the risk of fire disturbance. Canadian Journal of Forest Research 35: 1378-1388. Pickens, J.B., and Hof, J.G. 1991. Fuzzy goal programming in forestry - an application with special solution problems. Fuzzy Sets Syst. 39: 239-246. Pickens, J.B., and Dress, P.E. 1988. Use of stochastic production coefficients in linear- programming models - objective function-distribution, feasibility, and dual activities. For. Sci. 34: 574-591. Pickens, J.B., Hof, J.G., and Kent, B.M. 1991. Use of chance-constrained programming to account for stochastic variation in the A-matrix of large-scale linear programs: A forestry application. Ann Oper Res 31: 511-526. Pretzsch, H., Biber, P., and Dursky, J. 2002. The single tree-based stand simulator SILVA: construction, application and evaluation. For. Ecol. Manage. 162: 3-21. Pukkala, T. 1998. Multiple risks in multi-objective forest planning: integration and importance. For. Ecol. Manage. 111: 265-284. Reed, W.J., and Errico, D. 1986. Optimal harvest scheduling at the forest level in the presence of the risk of fire. Canadian Journal of Forest Research 16: 266-278. Schuetz, RJ. 2002. Type II forest estate analysis for tree farm license 48. Canadian forest products Ltd. Available from http://www.for.gov.bc.ca/hfp/silstrat/pdffiles/prgeorge- TFL-48.pdf. Accessed on January 20, 2008] Soyster, A.L. 1973. Convex programming with set-inclusive constraints and applications to inexact linear-programming. Oper. Res. 21: 1154-1157. Trincado, G., Quezada, R., and von Gadow, K. 2003. A comparison of two stand table projection methods for young Eucalyptus nitens (Maiden) plantations in Chile. For. Ecol. Manage. 180: 443-451. von Gadow, K. 2000. Evaluating risk in forest planning models. Silva Fenn. 34: 181-191. Weintraub, A., and Abramovich, A. 1995. Analysis of uncertainty of future timber yields in forest management. For. Sci. 41: 217-234. Weintraub, A., and Vera, J. 1991. A cutting plane approach for chance constrained linear- programs. Oper. Res. 39: 776-785. 56 Yaussy, D.A. 2000. Comparison of an empirical forest growth and yield simulator and a forest gap simulator using actual 30-year growth from two even-aged forests in Kentucky. For. Ecol. Manage. 126: 385-398. 57 3 Bi-objective Multi-Period Planning with Uncertain Weights: A Robust Optimization Approach3 3.1 Introduction Forest managers usually face the challenge of conjugating the views, goals and demands of different interest groups involved in the decision-making process. This challenge arises when the forest outcomes of the different objectives have to be estimated, and when defining the objective weights and trade-offs among the planning goals. In this paper, we focus on finding decisions that perform well even when these weights are not accurately known. These decisions are therefore called robust decisions. Multiple criteria decision making (MCDM) techniques provide a formal way to incorporate multiple criteria in the planning process (Steuer 1986). In particular, multi-objective programming methods (MOP) (Cohon 1978) include techniques applied to optimization problems as distinguished from multi-attribute utility theory methods (MAUT) used when a discrete set of alternatives are evaluated (Steuer 1986). It is worth noting that MAUT explores and ranks decision alternatives previously identified by the decision-maker, that is, no optimization of decisions is performed in contrast to MOP that generates a solution that optimizes a set of different objectives. A detailed explanation of these techniques can be found in Cohon (1978) and Steuer (1986), and for recent reviews of multiple-criteria methodologies in forest management the reader is referred to Diaz-Balteiro and Romero (2008) and Ananda and Herath (2009). One of the simplest MOP techniques is the weighted sum approach in which different objectives are combined in a single objective by assigning to each of them a positive weight that represents its relative importance or preference relationship. Although this approach looks quite straightforward, its successful application relies on the appropriateness of the weights (Steuer 3 A version of this chapter has been submitted for publication. Palma, C.D., and Nelson, J.D. 2009. Bi-objective multi-period planning with uncertain weights: A robust optimization approach. 58 1986) as these values affect the solution obtained (Butler et al. 1997, Kangas 1994, Pukkala 1998). Weights estimation imposes important difficulties (Cohon 1978) and different approaches to determine them have been proposed, including rating and ranking methods (Eckenrode 1965), pairwise comparison between objectives (Eckenrode 1965), fuzzy set theory (Rao and Roy 1989) and mathematical programming (Pekelman and Sen 1974, Ma et al. 1999). In addition to this complexity, the long time frame of forest management decisions adds more uncertainty to the weights estimation as they may change in the future. This uncertainty in the value of the weights has been addressed with sensitivity analysis (Butler et al. 1997, Kangas 1994), a method called stochastic multi-criteria acceptability analysis (SMAA) (Lahdelma et al. 1998, Kangas et al. 2003, Kangas 2006, Kangas et al. 2006), and using a combination of simulation and heuristic search (Pukkala and Miina 1997). All these approaches, except the latter, have been applied to situations where a set of previously known decision alternatives are evaluated under uncertain weights. Apparently, no applications have explored uncertainty in the objective weights in optimization problems (Diaz-Balteiro and Romero 2008, Ananda and Herath 2009), those in which an optimal decision is generated rather than many decisions evaluated. To our knowledge, only uncertain future state of the forest by using multi-objective dynamic programming (Gong 1992) and vagueness in the objective definitions by using fuzzy set theory (Ells et al. 1997) have been explored in multi-objective optimization problems. Our interest in this work is to find solutions to bi-objective problems that perform well under different values of objectives weights. Since weights play a major role in determining the model solutions and extensive sensitivity analysis are suggested when using this approach (Cohon 1978), we want to explore a way of formulating a model that explicitly considers the sensitivity of the solutions to the weights when looking for “optimal” decisions. We present a robust optimization approach (Bertsimas and Sim 2004) applied to a bi-objective multi-period harvest scheduling problem when weights are not accurately known and may change within a pre- defined range throughout the planning horizon. We demonstrate that robust decisions perform almost as well as the traditional fixed-value weight approach in terms of the objective function value, however, they are more consistent in a scenario of uncertain future weights. 59 3.2 Robust optimization approach Different robust optimization approaches that seek solutions that are not overly sensitive to the realization of uncertainties have been proposed (Mulvey et al. 1995, El Ghaoui et al. 1998, Ben- Tal and Nemirovski 2000, Bertsimas and Sim 2004). By modifying the structure of the original models, robust optimal solutions produce the best objective value that simultaneously tolerates changes in the model parameters up to a given bound, known a priori. Although the first attempts to immunize solutions against uncertainty were too conservative as they considered the worst-case parameters values (Soyster 1973), newer approaches have addressed this issue and less conservative solutions can be obtained (Ben-Tal and Nemirovski 1998, El Ghaoui and Lebret 1997). The disadvantage of these approaches is that they increase the complexity of the problem. The approach we use here does not increase complexity, e.g. the robust counterpart of linear problems remains linear, allowing us to solve large-scale problems in a computationally tractable way (Bertsimas and Sim 2004). This approach has been applied, for instance, to network design (Bertsimas and Sim 2003), portfolio management (Bertsimas and Sim 2004), inventory theory (Bertsimas and Thiele 2006) and forest harvest scheduling (Palma and Nelson 2009). The latter study considered uncertainty in timber yield and demand in the context of a single objective model. The approach works by adding a new term, called the protection function, to each constraint for which protection against the realization of uncertain coefficients is desirable. Note that if uncertainty is in the objective function, as in our case, we can still express it as a constraint. To do so, we add the constraint c’x ≥ z and maximize z, where c is the vector of uncertain coefficients of the objective function and x is the vector of decisions. The robust version of this constraint adds a protection function, β(x*,Γ), so the new constraint becomes [3.1] c`x - β(x*,Γ) ≥ z The protection function acts as a buffer and its magnitude depends on the user-defined parameter Γ and the observed value of the decision vector, x*. Since this buffer is required to be as low as possible in the optimal robust solution, decision vectors that satisfy all the model constraints at the “minimum cost” for the objective are to be selected. 60 We will assume that each coefficient ci is independently distributed in the symmetric interval ii cc ˆ± where ic ĉis the coefficient estimate and the error of the estimate. Although independency may represent an important assumption as discussed in the Discussion section, it allows keeping the models solvable. For each of the n coefficients, we can define a scaled deviation i iiii ccc ˆ/)( −=η so that they belong to [-1, 1]. The sum of all ηi could theoretically take values between –n and n if all the coefficients are uncertain. However, as uncertainties are assumed independently distributed, some coefficients will exceed and others fall below their point estimates, tending to cancel each other out, therefore narrowing the observed value of the sum of scaled deviations. Using the protection level Γ, we can limit the number of uncertain coefficients allowed to change by considering Γ≤∑i iη . If Γ=0, ηi are forced to 0 for all i and so there is no protection against uncertainty. If Γ=n, all ηi are allowed to have a non-negative value and the constraint is totally protected against uncertainty. Finally, for Γ∈(0, n) partial protection against uncertainty is defined as some of the ηi can be non-zero. For a given solution x*, the protection function β(x*,Γ) is defined to allow the simultaneous change of up to Γ uncertain coefficients, that is, it corresponds to the sum of the Γ largest deviations produced by a fixed value of the decision vector. In other words, the protection function equals the objective value of the following optimization problem: ∑i iii wxcMax ˆ]2.3[ * iw ∀≤≤ 10]4.3[ Subject to Γ≤∑i iw]3.3[ i Where wi is a new variable that represents the random scaled deviation described above. Note that this model is linear when seen as a sub-problem with a fixed x*, but it turns non-linear when x is variable as in the original problem. Its dual, however, can be used to linearly include this sub-problem into the original model. If u ≥ 0 and vi ≥ 0 are the dual variables of constraints [3.3] and [3.4], respectively, the dual of model [3.2]-[3.4] is: 61 ∑+Γ i ivuMin]5.3[ Subject to ixcvu iii]6.3[ ∀≥+ *ˆ ≥−Γ− ixcvu ∀≥+ ˆ]8.3[ iv ∀≥ 0]9.3[ 0]10.3[ ≥u Since the objective value of models [3.2]-[3.4] and [3.5]-[3.6] is the same, the latter can be embedded in [3.1] to replace the protection function β(x*,Γ). Thus the robust model constraint [3.1] is replaced by the following set of constraints while z is maximized: zvuxc i ii ii ∑∑]7.3[ iii i The robust formulation takes into consideration the trade-off between the quality of the decisions and the variability of their outcomes. Solutions that maximize the original objective function c’x, but simultaneously require small buffers against uncertain values will be preferred. For the range of possible values of ci, z represents the minimum value that the objective function c’x might reach if up to Γ parameters are allowed to change. 3.3 Methods We considered two representative objectives of the community desires in the study area in a multi-period forest planning problem, i.e. the creation of employment and the increase of the proportion of old forest in the area. In this section we describe the study area, the data used and present both the deterministic and the robust models. 3.3.1 Study area The study area consists of a 245,090 ha forested area within Tree Farm License 48, located in northeastern British Columbia, Canada. The area is divided in 800 strata with an average size of 62 306.4 ha (range 92.1-4,255.5 ha). Regardless of the spatial location of the stands, those with similar age-class, species composition and harvest system were grouped into strata. For each stratum, the inventory and the employment required based on the harvest system were available. 3.3.2 Bi-objective multi-period planning model For simplicity we used a model I formulation (Johnson and Scheurman 1977) in which we have to determine the number of hectares of each stratum i to be managed with management option j, Xij. We used a twenty period planning horizon (10 year periods) and defined 15 possible management options for each of the 800 strata. For each stratum, a management option defines the complete sequence of harvests throughout the planning horizon and their corresponding outputs in each period (employment in this case). This sequence, for instance, may state that a stratum will be harvested in periods p1 and p2. Another sequence or management option may consist of harvesting the stratum on period p3 only, or even not harvest it at all. Thus, different combinations of valid rotation ages (that is, subject to a minimum harvest age) throughout the planning horizon were considered. For the purpose of explaining our methodology we considered 15 management options for each stratum, but the effect of using more options is discussed in the Discussion section. As will be explained later, the two objectives are technically addressed by representing each of them through the 20 periods of the planning horizon, leading to 40 individual objectives. If eijt and aijt represent the amount of employment produced and the age of the stratum i in period t when managed with j, respectively, we can express the objectives level of each period as txeze ji ijijtt ∀= ∑ , ]11.3[ and ( ) txSzo ijtaji ijt ∀= ∑ ≥80|, 1]12.3[ where zet and zot are the amount of employment and the proportion of old forest in the area in period t, respectively, and S is the total area of the forest. In this case 80 years old strata and older were considered old-aged strata. 63 In addition, the area assigned to management options must equal the total stratum area, that is, isx i j ij∑]13.3[ ∀= 19.0]14.3[ >∀⋅≥ tzeze where si is the area of stratum i. We also constrained the reduction in employment rates between consecutive periods to avoid abrupt differences. We included the following constraint that guarantees that the employment level in a period cannot decrease more than 10% of the level observed in the previous period. )1( −tt For the objective function we used a weighted-sums approach in which each objective is multiplied by a positive scalar weight and then the weighted objectives are summed to form a composite objective function (Steuer 1986). As we considered that weights may change over time, they were defined by objective and period, in other words, one weight for each of the 40 individual objective levels. To overcome the effect of the relative magnitudes of the different objectives function gradients, we normalized the objectives achievement based on their ideal and anti-ideal levels (Martinson 1993), and defined our objective function as ( )∑ t tttt zozeMax]15.3[ + oe '' λλ where λ is the vector of weights and z’ is the vector of normalized objectives, that is −+ − −= tt tt t zeze zezeze )(]16.3[ ' − and −+ − −= tt tt t zozo zozozo )(]17.3[ ' − where the vectors z+ and z- are the ideal and anti-ideal objective levels. These levels are obtained by solving the individual single objective problems, in this case, by maximizing each objective for each period. The anti-ideal employment and the ideal proportion of old-aged forest levels are quite simple to estimate. In the first case, no employment is produced in any of the planning 64 periods if the non-harvest decision is selected in all strata. This would be the worst case from the employment creation point of view. In the second case, the ideal proportion of old forest corresponds to the projection of the initial proportion based on the strata age, which is also the case when the non-harvest decision is implemented. Ideal employment and the anti-ideal proportion of old forest levels have to be determined by solving the corresponding linear programming models. 3.3.3 Modeling unclear preferences To deal with uncertain weights we formulated the robust version of the objective function as described before. We note that although the sum of the weights is usually set to one, what really matters is the relative weights of the objectives (Cohon 1978). Two problems made us disregard this widely used way to define weights. First, as random changes in the weights are allowed, the value of the second weight (or the last one when more than two objectives are considered) will be perfectly correlated as they have to sum up 1. This violates the independency assumption mentioned before and could be avoided if weights can freely take values within their ranges. However, doing this raises a second problem. If we scale the random weights produced, it might happen that the new scaled weights fall outside their range of possible values. For instance, let us consider without lost of generality two scaled weights 0.6 and 0.4 that can take values in the range [0.48-0.72] and [0.32-0.48], respectively. If 0.49 and 0.47 are the random, non-scaled values of each weight, respectively, then the new scaled weights would be 0.51 and 0.49, with the last weight out of its range. As shown, it becomes unclear the real values and relative importance that weights can take when many weights change at a time. We therefore used a modified weight vector λ’. We set for all periods, while allowing only the other weight, , to change. In this way the user knows exactly how the relative importance of the two weights varies. The extension to this approach to more objectives is commented in the Discussion section. By fixing the employment weight for each period, the other weight can be obtained as 1t ' =eλ 'oλt eoo '' eo ttt λλλ = where . Therefore, the robust version of the objective function has 20 uncertain coefficients given by the uncertain weights of the proportion of old- aged forest throughout the planning horizon. The final model then corresponds to equations [3.11]-[3.17] in which equation [3.15] was replaced with the following set of equations: 1=+ tt λλ 65 wMax]18.3[ Subject to ( ) wpyzoze t t t tttt ≥−Γ−+ oe ∑∑]19.3[ λλ '''' o ''ˆ negativenon,]21.3[ −py 'ˆo 'o 'oλ '' ôo 'e 'o '' ooo λλλ += oe λλ −= tzopy ttt ∀≥+]20.3[ λ t where Γ is the protection level and is the error associated with the weight estimate . In other words, the weight is assumed to belong to the symmetrical and uniformly distributed range . Recall that is fixed to one and it is not allowed to change. In what follows, w will be referred as the objective function, and the weighted sum of objectives using the original vector λ will be referred to as the solution score. As the sum of modified weights (λ’) is variable, the score is a more appropriate way to compare solutions. For a variable , original weights can be obtained, for instance, as and . tλ tλ t tt λλ ± tλ tλ )1/( ttt tt 1 The new formulation adds t+1 new constraints and variables, plus the auxiliary variable w to the original model. Models were implemented in ILOG OPL Development Studio version 5.2 (CPLEX optimizer 10.2.0) and run in a 3.2 Ghz Pentium D with 1 GB of RAM. 3.3.4 Uncertain weights and analysis Original weights (λ) were set to 0.6 and 0.4 for the employment and proportion of old-aged forest objectives, respectively, in all the periods. Consequently, the modified weights (λ’) were 1 and 0.667 for employment and old forest in each period, respectively. That means that although employment is the most important management objective in the area, the existence of an old forest is also desirable. The latter objective guarantees good quality wildlife habitat in the area and recreation spots for the community. The use of any other weight combination does not affect the methodology proposed. 66 We considered two scenarios of weight uncertainty. In both cases we assumed that weights are not totally known even in the first planning period, and that this imprecision gradually increases over time until period 10. For the first planning period, was defined as 20% (scenario SCE20) and 50% (scenario SCE50) of in the first and second scenario, respectively. That is, the value of was allowed to change in the range 0.667±0.133 and 0.667±0.333 in the two scenarios, respectively. For periods 2 through 10, the percentage of estimate error increased by two units per period, reaching an error of 38% and 68% for SCE20 and SCE50, respectively. These values were maintained for the rest of the planning horizon, as we considered them wide enough to represent long-term uncertainties. For SCE20, the first period estimate error was defined in such a way that original weights (λ) showed a relative agreement with the real values of the objectives weights (Fig. 3.1a). For SCE50, this error was the minimum required to produce changes in the first period harvest decisions as discussed later (Fig. 3.1b). 'ô 'o 'oλ tλ tλ 1 Fig. 3.1. The two scenarios described two levels of uncertainty for the objective preferences over time. The deterministic solution was compared with robust solutions generated with three levels of robustness. The protection level parameter Γ, which represents the number of uncertain coefficients allowed to change, was set to 10, 20 and 30% of the number of uncertain coefficients in the objective function. In other words, the three levels of robustness were defined by setting Γ=2, 4 and 6. We refer to these levels of robustness as ROB10, ROB20 and ROB30, respectively. These values were sufficient to provide nearly full protection against the value of the objective function as detailed later. 67 We used 1,000 sets of simulated weights (500 by scenario, each set containing weights for all the periods of the planning horizon) to evaluate the performance of deterministic and robust solutions in terms of their stability when facing uncertain future objective weights. We also compared deterministic and robust solutions in terms of the change in the optimal score, level of employment and proportion of old-aged forest, and change in the management decisions. 3.4 Results Robust solutions were less sensitive to changing objective weights than deterministic solutions. In Fig. 3.2, the black lines represent the weighted sum of objectives by period, or score, produced by the optimization models using the initial fixed-value weights, and the gray area represents the scores produced by the same optimal decision but evaluated using the simulated weights. In both scenarios of simulated weights, the score obtained in each period by robust decisions (Fig. 3.2b and d) was more stable than the one obtained by deterministic decisions (Fig. 3.2a and c). The other levels of robustness not shown in the figure (ROB10 and ROB20) showed an intermediate behavior between deterministic and robust decisions, that is, the more the level of robustness the more stable the score. 68 Fig. 3.2. The score dispersion (gray area) was wider for the deterministic solutions (a and c) than for robust solution ROB30 (b and d) when optimal decisions (black line) were evaluated using simulated weights. The value of the objective function (with modified weights λ’) decreased when more protection was required (Table 3.1) as it protected different values of the objectives weights. Under simulated scenarios, the value of the deterministic objective function (14.08) was not reached about 50% of the time, unlike the value of the robust objective functions that were achieved in almost all the simulations (ROB30). However, a reduced impact in the optimal score (solution evaluated with λ) was observed in robust solutions. In the most extreme case of uncertainty and protection level (SCE50, ROB30) this reduction was close to 1%, in all other cases it was less than 0.5%. 69 Table 3.1. Value of the objective function and scores produced by deterministic and robust solutions Scenario SCE20 Scenario SCE50 DET ROB10 ROB20 ROB30 DET ROB10 ROB20 ROB30 Objective Function (OF) 14.08 13.79 13.52 13.27 14.08 13.46 12.90 12.40 Simulations below OF (%) 49.0 17.1 3.3 0.3 51.0 19.6 3.6 0.3 Score (reduction) 8.45 (0%) 8.44 (0.12%) 8.43 (0.24%) 8.41 (0.47%) 8.45 (0%) 8.44 (0.12%) 8.41 (0.47%) 8.37 (0.95%) Deterministic and robust models scheduled harvests in a different way, that is, some of the strata were scheduled to be harvested in a different period. The change in harvest decisions relative to deterministic decisions (Fig. 3.3) was more important when more uncertainty was assumed (Fig. 3.3b) and also when more protection was required. In the scenario of more certainty in weights (Fig 3a) no different harvest decisions were scheduled in the first period, about 250 ha in the second period and 5,800 ha per period on average in the rest of the planning horizon with ROB30. In SCE50 (Fig. 3.3b), almost 1,300 ha were harvested differently in period 1 with ROB30, and up to 16,000 ha differences were observed in future periods. Fig. 3.3. Changes in harvest decisions relative to deterministic decisions were larger as more robustness was required and more uncertainty was assumed (b). 70 Deterministic solutions produced less employment and a higher proportion of old forest in both scenarios (Fig. 3.4). Consistently, as more robustness was required more employment and less proportion of old forest were produced. In the most uncertain scenario (SCE50) the change in the objective levels was slightly more pronounced. Fig. 3.4. In both scenarios of uncertainty, more employment (a,c) and lower proportions of old- aged forest (b,d) were produced as more robustness was required. Employment levels were closely related to the total area harvested. Peak and bottom points of harvests (Fig. 3.5), essentially caused by the forest structure and its dynamics, determined the peaks and bottoms of employment levels (Fig. 3.4a and c). The 10% decline constraints prevented the employment from drastically dropping in periods 6-15. The proportion of old forest increased (Fig. 3.4b and d, periods 1-3 and 11-15) in periods with low harvest levels (Fig. 3.5, same periods), and the opposite occurred when harvest levels increased. 71 Fig. 3.5. Area harvested through the planning horizon in the two scenarios of uncertainty. Robust formulations produced a minor increase in the model size. The number of variables increased from 12,041 in the deterministic formulation to 12,063 in the robust model. Constraints also increased from 859 to 880. The increase in solution time was imperceptible as all models solved in 2-2.5 seconds. 3.5 Discussion Our results suggest that the robust formulation of the bi-objective multi-period planning problem produces management decisions that are less sensitive to uncertain objective weights. Under the premise that social preferences may change in the future, and are even hard to estimate at the present time, robust models found decisions that produced more stable scores along the planning horizon. This guarantees that decisions will still produce a good and more consistent combination of the two objectives considered, independently of the relative importance that society gives to the objectives. Since robust decisions are meant to perform well in multiple scenarios of weights, the value of their score is certainly affected. As expected, the score reduction is greater as more uncertainty is assumed and more robustness is desired. However, in our case this reduction was less than 1% in the most uncertain scenario, suggesting that robust solutions do not drastically score lower than traditional deterministic solutions. This is possible because the protection function introduced in the robust formulation is an optimization problem in itself, embedded in the original model, that 72 identifies less sensitive decisions to changing weights and, at the same time contributes to the model objective. As suggested by Pukkala and Miina (1997), considering uncertainty in the objective weights affected the management decisions. Although the solution score was slightly reduced in robust models, management decisions showed important differences in both scenarios of uncertainty. Thousands of hectares were scheduled to be harvested in a different period, which is also evidence that robust models look for different attributes when selecting optimal decisions. In the most uncertain scenario, first-period harvest decisions that are to be implemented changed in about 1300 ha. That is, some strata scheduled for harvesting in the first period by the deterministic model were scheduled for a later period by the robust model, and vice versa. However, we note that first-period harvest decisions did not change in the scenario with low level of uncertainty. This suggests that uncertain objective weights might be overcome in some cases by a continuous planning process that keeps track of preference changes. The magnitude of the decision change is expected to be scenario specific in our view and influenced by the forest state and the set of management options available for each stratum. We considered 15 management options for each stratum, but if more options are provided we would expect more changes in the management decisions and less impact in the solution score. To reduce the negative impact on the score produced by increasing levels of robustness, robust decisions increased the amount of the most weighted objective, in this case the level of employment. As a trade-off between optimality and stability of the solutions is desirable, the increase in the level of employment through the planning periods guarantees that the value of the objective function will not be excessively affected. Higher levels of employment therefore compensate the natural reduction of the objective function that would be produced if a non- robust solution is generated. The robust version of a constraint guarantees that the constraint will be met for a range of realizations of the uncertain coefficients. When dealing with the objective function, robustness translates into a guarantee that the value of the objective function will be achieved for a range of uncertain coefficients. We think that this guarantee is an important concept in a multi-objective framework where usually public areas are managed to satisfy the needs of a variety of social actors. Since social preferences may not be accurately determined, decision makers can use this 73 approach to find less risky decisions and therefore reduce negative outcomes if preferences change. Our methodology finds less sensitive decisions to uncertain objective weights in bi-objective optimization problems. Current methodologies that deal with uncertain weights have focused on the comparison of previously defined decision alternatives rather than in generating robust alternatives. Stochastic multi-criteria acceptability analysis (Lahdelma et al. 1998), for instance, assigns to each alternative an acceptability index representing the variety of possible weights that support the preference of that alternative. A higher index means that the alternative is preferred in more scenarios of uncertain weights. Based on the same concept, other approaches (see for example Kangas 2006, Kurttila et al. 2009 and Liesiö et al. 2007) rank the robustness of the decision alternatives providing decision-makers with valuable information about their stability. However, the quality of the final decision is highly dependent on the set of alternatives available for comparison. Moreover, when a decision is required at a landscape level, that is, a plan consisting of selecting one alternative for each stratum in the landscape, the number of landscape-level alternatives gets large very quickly. This number might easily be much larger than the number of alternatives under comparison reported in the literature. We believe that when decisions can be represented by continuous variables at stratum level the robust formulation proposed here may provide a practical way to find good quality robust decisions without the need of enumerating decision alternatives. The approach proposed includes some assumptions that should be addressed in future research. The assumption that uncertain coefficients have to be independently distributed limits the type of uncertainty to be considered to random and unbiased errors, discarding the consideration, for instance, of uncertain trends in social preferences. In addition, the assumption that estimate errors are to be uniformly distributed denies the use of information about the probability distribution of these coefficients. Although this information is rarely available, one can define a wide enough range of uncertain values that contains a high proportion of the observations of the original distribution (Palma and Nelson, 2009). Another drawback of the approach is the lack of a clear relationship between the protection level, Γ, and the resulting level of robustness or degree of constraint satisfaction. Probability bounds to estimate the value of Γ needed to obtain a certain degree of robustness have been proposed (Bertsimas and Sim 2004), but they are loose and produce conservative solutions (Palma and Nelson, 2009). An alternative way to explore this relationship is by conducting simulation experiments. Solutions obtained using different 74 protection levels can be tested in simulated scenarios of uncertain coefficients and their robustness evaluated. Although conducting simulation experiments is quite simple, a more accurate way to determine the protection level required to obtain a desired level of robustness would improve this methodology. The possible values that a weight can take are intended to reflect what the modeler expects about the change of the relative importance of an objective relative to other objectives. Allowing weights to change within a certain range of values may produce confusing results in the final relative importance of weights as shown in the Method section. We avoided this by allowing only one of the two objectives to change while keeping the other weight fixed to one. Although this approach works well for two objectives, when considering more objectives the relative importance of the weights is still uncontrollable and special care has to be taken if this represents a problem for the modeler. We emphasize that the use of non-scaled weights does not affect the model solution, but only the value of the objective function. This value is also not expected to mean something rather than guiding the solution search process to find non-inferior solutions (Cohon 1978). As we did in our case, there always exists the chance to evaluate and compare solutions obtained with non-scaled weights by re-calculating the objective value using a set of scaled weights. Despite the limitations of the methodology, we consider that it has important benefits. Since good solutions (not optimal ones) for a range of objective weights are obtained, decision makers can reduce the risk of obtaining poor outcomes and will need less adjustment to their decisions if future conditions change. This approach would also contribute to reduce the efforts needed in conducting sensitive analyses. Decision makers would not need to evaluate different weights combinations before selecting a solution. The robust approach explicitly considers different weight combinations and based on an objective criterion provides a stable solution. Finally, we envision that these benefits can translate into a better acceptance of management decisions from stakeholders or other social actors involved in this type of decision problems. 3.6 Conclusions We showed a robust formulation of a bi-objective multi-period problem in which uncertainty in the objective weights was considered. Our method produced solutions that are less sensitive to future change in the weights than traditional deterministic formulations. Under simulated 75 scenarios of objective weights, the score or weighted sum of the objectives obtained in each period by robust decisions was less variable than deterministic decisions at the expenses of a minor reduction in the objective function score. The robust formulation presented may significantly reduce the efforts of sensitivity analyses as it explicitly considers how decisions change when weights change. We consider that decisions that involve a consensus about the importance of the planning objectives may be more easily accepted by the different actors involved in the decision-making process with this approach. 76 3.7 References Ananda, J., and Herath, G. 2009. A critical review of multi-criteria decision making methods with special reference to forest management and planning. Ecol Econ 68: 2235-2248. Ben-Tal, A., and Nemirovski, A. 2000. Robust solutions of linear programming problems contaminated with uncertain data. Math Program 88: 411-424. Ben-Tal, A., and Nemirovski, A. 1998. Robust convex optimization. Math Oper Res 23: 769- 805. Bertsimas, D., and Sim, M. 2003. Robust discrete optimization and network flows. Math Program 98: 49-71 Bertsimas, D., and Sim, M. 2004. The price of robustness. Oper Res 52: 35-53. Bertsimas, D., and Thiele, A. 2006. A robust optimization approach to inventory theory. Oper Res 54: 150-168. Butler, J., Jia, J.M., and Dyer, J. 1997. Simulation techniques for the sensitivity analysis of multi-criteria decision models. Eur J Oper Res 103: 531-546. Cohon, J.L. 1978. Multiobjective programming and planning. Academic Press, New York. Diaz-Balteiro, L., and Romero, C. 2008. Making forestry decisions with multiple criteria: A review and an assessment. For Ecol Manage 255: 3222-3241. Eckenrode, R.T. 1965. Weighting multiple criteria. Manage Sci 12: 180-192. El Ghaoui, L., Oustry, F., and Lebret, H. 1998. Robust solutions to uncertain semidefinite programs. Siam J Optim 9: 33-52. El Ghaoui, L., and Lebret, H. 1997. Robust solutions to least-squares problems with uncertain data. Siam J Matrix Anal Appl 18: 1035-1064. Ells, A., Bulte, E., and vanKooten, G.C. 1997. Uncertainty and forest land use allocation in British Columbia: vague priorities and imprecise coefficients. For Sci 43: 509-520. Gong, P. 1992. Multiobjective dynamic programming for forest resource management. For Ecol Manage 48: 43-54. Johnson, K.N., and Scheurman, H.L. 1977. Techniques for prescribing optimal timber harvest and investment under different objectives – discussion and synthesis. For Sci 1-31. Kangas, J. 1994. An approach to public-participation in strategic forest management planning. For Ecol Manage 70: 75-88. Kangas, A. 2006. The risk of decision making with incomplete criteria weight information. Can J for Res 36: 195-205. 77 Kangas, A.S, Kangas, J., Lahdelma, R., Salminen, P. 2006. Using SMAA-2 method with dependent uncertainties for strategic forest planning. For Policy Econ 9(2): 113-125. Kangas, J., Kurttila, M., Kajanus, M., and Kangas, A. 2003. Evaluating the management strategies of a forestland estate - the S-O-S approach. J Environ Manage 69: 349-358. Kurttila, M., Muinonen, E., Leskinen, P., Kilpeläinen, H. & Pykäläinen, J. 2009. An approach for examining the effects of preferential uncertainty on the contents of forest management plan at stand and holding level. Eur J For Res 128: 37-50. Liesiö J, Mild P, Salo A 2007. Preference programming for robust portfolio modeling and project selection. Eur J Oper Res 181:1488-1505. Lahdelma, R., Hokkanen, J., and Salminen, P. 1998. SMAA - Stochastic multiobjective acceptability analysis. Eur J Oper Res 106: 137-143. Ma, J., Fan, Z., and Huang, L. 1999. A subjective and objective integrated approach to determine attribute weights. Eur J Oper Res 112: 397-404. Martinson, F.K. 1993. Fuzzy v/s minmax weighted multiobjective linear-programming: illustrative comparisons. Decis Sci 24: 809-824. Mulvey, J.M., Vanderbei, R.J., and Zenios, S.A. 1995. Robust pptimization of large-scale systems. Oper Res 43: 264-281. Palma, C.D., and Nelson, J.D. 2009. A robust optimization approach protected harvest scheduling decisions against uncertainty. Can J For Res 39: 342-355. Pekelman, D., and Sen, S.K. 1974. Mathematical programming models for determination of attribute weights. Manage Sci Ser B-Application 20: 1217-1229. Pukkala, T. 1998. Multiple risks in multi-objective forest planning: integration and importance. For Ecol Manage 111: 265-284. Pukkala, T., and Miina, J. 1997. A method for stochastic multiobjective optimization of stand management. For Ecol Manage 98: 189-203. Rao, J.R., and Roy, N. 1989. Fuzzy set theoretic approach of assigning weights to objectives in multicriteria decision-making. Int J Syst Sci 20: 1381-1386. Soyster, A.L. 1973. Convex programming with set-inclusive constraints and applications to inexact linear-programming. Oper Res 21: 1154-1157. Steuer, R.E. 1986. Multiple criteria optimization: theory, computation, and application. Wiley, New York. 78 4 A Robust Model to Protect Road Building and Harvest Decisions from Timber Estimate Errors4 4.1 Introduction Medium-term decisions typically involve harvest scheduling decisions and road access construction. Although these two aspects of the planning problem can be considered separately, their integration has proved to be significantly more efficient (Jones et al 1986). Different mixed binary models to solve this problem have been proposed and solved using different approaches (Weintraub and Navon 1976, Richards and Gunn 2000, Guignard et al. 1998, Andalaft et al. 2003), but in all cases assuming the model data are perfectly known and ignoring the inherent uncertainty in the coefficients of the model. Uncertainty in natural resources management comes from a range of sources well described in the literature (Mowrer 2000, Regan et al. 2002, Marshall 1987). Biological processes that are not totally understood and changing social and economic conditions in combination with long time horizons affect the consequences of forest planning decisions made today. In addition, uncertainties result from the process of estimation of current and future resource levels. Different sampling methods will report different estimate errors of current inventory, and statistical models used to project this inventory also contribute to this inexactness. Although extensive areas might produce a diversification effect that can reduce the effect of some of these uncertainties, the decision making process still has to be carried out in an uncertain environment. As a consequence, although decisions implemented now might seem optimal, they will probably be sub-optimal once the uncertainties are realized as the information used to make these decisions will not necessarily be observed (Marshall 1987). In the context of sequential decision making that takes place in the real world, there is potential recourse when forecasted outcomes are not realized and the consequences may not necessarily be large. Since many of the real system 4 A version of this chapter will be submitted for publication. Palma, C.D., and Nelson, J.D. 2010. A robust model to protect road building and harvest decisions from timber estimate errors. 79 decisions are guided by model forecasts, we can gain understanding of uncertainty in the real system by examining feasibility of the models when information changes. Solutions of deterministic constrained models are likely to become infeasible (Pickens and Dress 1988) when evaluated with observed data (e.g. upper or lower bound constraints may become infeasible). This makes the search for good and stable solutions rather than strictly optimal ones highly relevant for implementing decisions. 4.2 Approaches to deal with uncertainty Stochastic programming (SP), chance-constraint programming (CCP) and fuzzy set theory (FS) are the best known techniques to include uncertainty in optimization models. By considering a set (usually discrete) of scenarios of random events or uncertainties, SP allows finding a good solution for all or the most likely of them. The great disadvantage of this technique is that as the number of scenarios increases, models become computationally difficult to solve and solution methods like decomposition and statistical approximation have to be used (Birge and Louveaux 1997). In addition, the probability distribution of the scenarios needs to be known, which is rarely the case. In CCP, constraints with at least one random coefficient are modeled as probabilistic statements. Probability distributions of uncertain coefficients are assumed to be known and constraints are required to hold with a minimum probability, exogenously determined. Although particular cases of CCP models are easy to solve (when only the right-hand side of a constraint is uncertain), in most cases models become nonlinear (Birge and Louveaux 1997, Kall and Wallace 1994), and finding exact solutions is typically intractable (Chen et al. 2007), thus motivating the search for approximate solution techniques (Birge and Louveaux 1997). FS provides a different approach to deal with some types of uncertainty by defining a degree of membership of a parameter to a set of possible values. The more likely the value of a parameter, the greater its degree of membership. This concept is used to model objectives and constraints for which a combination of their degree of membership can be maximized (Zimmermann 1996). This approach may be appropriate when there is vagueness, for example, in the meaning of certain events, phenomena or statements, but its appropriateness is not clear when representing the lack of information about the value of the parameters. In the latter case, FS actually results in 80 a relaxed version of the traditional deterministic problem in which constraint violations are allowed and, consequently, better objective functions are obtained. These three techniques have been mainly used in harvesting scheduling problems. For instance, SP has been used to deal with uncertainty in timber yield (Hoganson and Rose 1987, Eriksson 2006) and fire losses (Gassmann 1989, Boychuk and Martell 1996), while uncertainty in timber yield has been also considered by using CPP (Weintraub and Vera 1991, Weintraub and Abramovich 1995, Pickens et al. 1991, Hof et al. 1992). Most recently, CCP with approximation solution techniques was applied to a fire budget allocation problem with uncertain fire suppression costs, and to a habitat restoration problem where survival parameters were uncertain (Bevers 2007). Using FS, uncertainty in timber yield (Hof et al. 1986, Pickens and Hof 1991, Bare and Mendoza 1992) and vagueness in seral-class definitions (Boyland et al. 2006) have been considered in harvest scheduling problems. In the context of multi-objective models, fuzziness in the objective function coefficients (Mendoza et al. 1993, Stirn 2006) as well as in the goal definitions (Ells et al. 1997, Maness and Farrell 2004) has been addressed using FS. However, uncertainty in road building decisions has been scarcely addressed, probably because of the additional difficulty in solving integer models. Only road upgrade decisions with uncertain length of the period of poor road conditions has been examined using SP (Olsson 2007). The aim of this paper is twofold. First, we explore the effect of random volume uncertainties on harvesting and road building decisions; and second, we do this by using robust optimization (RO), a mathematical approach to immunize decisions against parametric uncertainty that has not been extensively used in forest resources management. It has been applied in engineering (Ben-Tal and Nemirovski 2002), network design (Bertsimas and Sim 2003, Ordonez and Zhao 2007), inventory theory (Bertsimas and Thiele 2006) and recently to forest harvest scheduling (Palma and Nelson 2009). 4.3 The robust optimization approach We base our development in the RO approach proposed by Bertsimas and Sim (2004). Other approaches have also been proposed (Mulvey et al. 1995, Ben-Tal and Nemirovski 2000, El Ghaoui et al. 1998), but they result in computationally more complex models. The approach used here produces models of the same type of the original model (i.e. linearity), without the need of specific solution techniques other than the one used to solve the original problem. 81 The RO approach works as a buffering strategy. An additional term, called the protection function, is added to each constraint affected by uncertainty and for which feasibility is highly desirable. Let us consider the general lower bound constraint: ibxa j ijij ∀≥∑ =1 ]1.4[ n a ]ˆ,ˆ aaaawhere are assumed uniformly distributed in the range ij [ ijijijij +− with ija â a the coefficient estimate and the accuracy of the estimate. If for each a scaled deviation is defined as ij ij ijijijij aaa ˆ)( −=η , then ]1,1[−∈ijη and for each constraint∑ =j ij1ηn can theoretically take values between –n and n. We can now limit the number of uncertain coefficients allowed to change in constraint i by considering ij ij Γ≤∑ =1ηn , where Γi is called the protection level of constraint i. Three possibilities can be identified: (i) if 0=Γi , the ijη for all j in constraint i are forced to 0, so that ijij a= na for all j, and constraint i has no protection against uncertainty; (ii) if i =Γ constraint i is totally protected against uncertainty since all ijη in constraint i are allowed to take a non-zero value; and (iii) if ni <Γ<0 constraint i is partially protected against uncertainty, in which case some of the ijη can be different from zero. The robust version of equation [4.1] is then: ibxa j iiijij ∀≥Γ−∑ =1 * ),x(]2.4[ βn * where , the protection function of constraint i, depends on the user-defined protection level ),x( ii Γβ and is specific for a given solution vector x*. RO looks for the optimal buffer to keep the left-hand side of [4.2] greater or equal than bi for different values of . For a given solution vector, the protection function determines the buffer required if up to ija iΓ parameters are allowed to change, that is, the buffer size corresponds to the sum of the iΓ largest deviations produced by a fixed value of x. In other words, the protection function of constraint i equals the objective function of the following optimization problem: 82 ∑ =j ijjij wxaMax 1 *ˆ]3.4[ n n jw Subject to i j ijw Γ≤∑ =1 ]4.4[ ij ∀≤≤ 10]5.4[ where wij are new variables that represent the random variable of scaled deviation ijη described above. Although this model is linear for the given solution vector x*, it is not the case when x is variable. Its dual, however, can be used to linearly express this function in [4.2]. If zi and pij are the dual variables of the constraints [4.4] and [4.5], respectively, then the dual of [4.3]-[4.5] is: ∑ = +Γ j ijii pzMin 1 ]6.4[ n * jp ∀≥ 0]8.4[ ≥z ],0[ n Subject to jxapz jijiji ∀≥+ ˆ]7.4[ ij 0]9.4[ i Since model [4.3]-[4.5] is feasible and bounded for all i ∈Γ , then its dual, model [4.6]- [4.9], is also feasible and bounded, and their objective functions coincide (Bertsimas and Sim 2004). Therefore, is equal to the objective function of model [4.6]-[4.9] and can be substituted into equation [4.2] to produce its robust version. ),x( ii Γβ * ibpzxa j i j ijiijij ∀≥−Γ−∑ ∑ = =1 1 ]10.4[ n n jixapz ,ˆ]11.4[ ∀≥+ jijiji 83 jip ,0]12.4[ ∀≥ iz ∀≥ 0]13.4[ ij i 4.4 Methods In this section we describe the study area, both deterministic and robust formulations, the data used and the simulation experiment performed to test the resulting infeasibility rates. 4.4.1 Study area We determined road and harvest decisions in a 11,675 ha area located on mid-Vancouver Island, British Columbia, Canada, of which 7,554 ha were available for harvesting (Fig. 4.1). The area is divided into 431 harvestable stands with an average size of 17.5 ha and ranging from 1.6 to 52.3 ha. We considered 412 potential road segments and two demand nodes for which minimum timber supply has to be guaranteed over three planning periods. To cover the planning horizon usually considered in tactical decisions while providing more detail in the first periods, the periods were defined as two, three and four years, respectively, that correspond to a nine-year planning horizon. Since a considerable portion of this area could be harvested in a relatively short period, 10% of the volume is to be retained when harvesting a stand. Inventory data and harvesting and transportation costs are available for the area. 84 Fig. 4.1. A harvestable area of 7,554 ha was used as study area. It included 431 stands and 412 potential roads. 4.4.2 Harvest scheduling and road building models The road building problem is a network problem in which arcs represent segments of potential roads and nodes represent either timber sources or the intersection of arcs. Roads have to be built and harvesting decisions, as well as timber flow, have to be determined in order to transport the harvested timber towards exit nodes that connect the forest to demand nodes. In our formulation we minimize the discounted harvesting and road building costs subject to a minimum timber demand. Adjacency constraints are not included, and only one timber product and one road standard are assumed. As integrated harvest scheduling and road building models are hard to solve to optimality, different solution approaches have been proposed to find good solutions, e.g. mixed-integer solvers (Weintraub and Navon 1976), heuristics and metaheuristics (Richards and Gunn 2000, Weintraub et al. 1995, Clark et al. 2000) and Lagrangean relaxation (Andalaft et al. 2003). In addition, strengthening and lifting techniques are also simple steps to facilitate the model solution. By adding logical inequalities (or ‘triggers’) and increasing the dimension of a constraint space (or ‘lifting’), the solution space of the LP relaxation is reduced and therefore fewer solutions are evaluated and better bounds for the integer solution techniques can be obtained (Guignard et al. 1998). We opted for a strengthened formulation of the problem and solved it using a commercial optimization software as detailed later. 85 The following nomenclature will be used in both the deterministic and robust formulations. Sets: i,j supply, intersection and exit nodes. This set does not include demand nodes. m demand nodes. r potential road. It corresponds to an undirected arc connecting two nodes. That is, each road r supports timber flow on directed arcs ij and ji. s stand. t period. s(i) set of stands that supply timber to node i. Coefficients and Parameters: mtd v R T H minimum demand in node m in period t. [m3] st total volume in stand s in period t. [m 3] rc cost of building road r. [$] tmijc )( cost of transportation between nodes i and j (or demand node m) in period t. [$/m 3] stc total cost of harvesting stand s in period t. [$] tα discount factor applied to period t. Decision Variables: rtX Y F binary variable that represents if road r is built (1) in t or not (0). st binary variable that represents if stand s is harvested (1) in t or not (0). ijt timber flow (directed) between node i and j in period t [m 3]. 86 imtF timber flow (directed) from node i to destination m in period t [m 3]. 4.4.2.1 Deterministic model We minimize the sum of road building, harvesting and transportation costs as follows: ∑∑∑∑ +++ tmi imtimtt tji ijtijtt ts ststt tr rtrt FcFcYcXcMin ,,,,,, ]14.4[ αααα TTHR Subject to the following set of constraints: sY t st ∀≤∑ 1]15.4[ tiFFFYv m imt j ijt j jit iss stst ,]16.4[ )( ∀+=+ ∑∑∑∑ ∈ tmdF mt i imt ,]17.4[ ∑ ∀≥ trjiXMFF tk rkrt tk jikijk ,),()(]18.4[ ∈∀≤+ ∑∑ ≤≤ tsXY tksrr rkst ,]19.4[ ),( ∀≤ ∑ ≤∈ Constraint [4.15] limits the number of times that each stand can be harvested during the planning horizon. Constraint [4.16] imposes the conservation of flow at each node (except demand nodes) and in each period. Inflow to a node (left-hand side) may come directly from harvesting stands or from other nodes, while the outflow (right-hand side) may go to other nodes in the network or to demand nodes. Constraint [4.17] defines the lower bound of timber production in each period and demand node. Constraints [4.18] and [4.19] correspond to the road construction trigger (lifted with respect to time) and a logical inequality, respectively. They both contribute to make a strong formulation of the harvesting and road building problem (Andalaft et al. 2003). In [4.18], timber flow in any arc and up to a given period is only possible if the corresponding road has been built in the same or in a previous period. Mrt is a big number that allows the flow to be greater than 1. Since its value has a major impact on the solution process of the model, the lowest possible value that preserves the original integer optimum should be used (Andalaft et al. 2003). When maximum demands are present, the sum of all of them can be used as an upper bound for 87 the timber flow in each period. If there are no maximum demand constraints, as in our case, the total volume of the forest in each period, ∑s stv , is another option. However, as we require a minimum demand in each period, not all the forest can be harvested in the same period because some timber has to be saved to meet demands in other periods. Therefore, a better flow upper bound for every arc and period t would be ∑∑ ≠− tlm mls st dv , . Moreover, the tightest upper bound can be determined when arcs have a treelike structure (Fig. 4.2). Unlike cycle structures, in treelike structures only one flow direction is possible and the potential flow in an arc cannot be larger than the maximum cumulative production of the source nodes from the leaves to the tree base. Fig. 4.2. In a tree structure (a) only one flow direction is possible so the maximum flow in an arc can be determined as a cumulative flow (Max Flowcb = Fc; Max Flowba = Fc + Fd). If cycles are present (b), the maximum flow in an arc is more difficult to determine as complex flow interactions might result. (Andalaft et al. 2003) For each arc and period, the flow upper bound was then determined as follows: ⎪⎩ ⎪⎨ ⎧ − += ∑∑ ∑ ≠ ∈∈ structure cycle ain isr if structure treein the isr if , ),(|)( tlm mls sst rjiiss sttr rt dav vM M where r is the preceding arc of r from the leaves to the tree base and node i is the node in r closest to the leaf. 0=M tr for terminal nodes. These calculations were automatically computed. 88 Constraint [4.19], called a project-to-road trigger (Andalaft et al. 2003), allows harvesting a stand in a period only if at least one road connecting the stand to the road network has been built in the same or in a previous period. In this equation, r(s) represents the set of roads that connects stand s to the road network. 4.4.2.2 Robust model To protect minimum demand constraints from infeasibilities, equation [4.17] was modified and new variables and constraints were added as previously described. We assumed that the real volume, , belongs to and is symmetrically distributed in the rangestv ]ˆ,ˆ vvvv +−[ , where stststst stv is the volume estimate and tstst evv ⋅=ˆ =δ tsYvpz ,ˆ]21.4[ ∀≥+ tsp ,0]22.4[ ∀≥ is its estimate error. Two special features of this problem deserve some attention. First, the network structure of the problem makes it such that uncertain coefficients (v ) are not in the equation to be buffered (eq. 4.17). Although only flow variables (F st imt) are present in this constraint, they are linked to the harvest variables Yst, which are associated with uncertain coefficients, through constraint [4.16]. The additional timber flow to be sent to each destination (eq. 4.17) forces the additional harvest (eq. 4.16) and the required road construction (eq. 4.18). Second, since uncertain coefficients do not depend on the demand destination m, the protection function βt will not provide a harvest buffer at the demand destination level, but only a total buffer at a period level. However, we can assign a portion of this total buffer to each demand by introducing the parameter δmt, with . In our case, both demand centers were assumed equally important, so δ1∑m mt mt took the value of the proportion that each destination represents from the total demand (i.e., 0.6 for node 1 and 0.4 for node 2 for all periods). Equation [4.17] of the original model was then replaced by the following set of constraints: tmdpzF mt s stttmt i imt ,]20.4[ ∀≥⎟⎠ ⎞⎜⎝ ⎛ +Γ− ∑∑ δ stststt st 89 tz ∀≥ 0]23.4[ t where zt and pst are new variables. All other equations of the original model remain the same. Models were implemented in OPL Development Studio 5.2 (CPLEX 10.2.0 as optimizer) in a 3.2 Ghz Pentium D with 1 GB of RAM. 4.4.3 Uncertain coefficients and model parameters A stable annual timber supply requirement was determined in such a way that around 70% of the total area needs be harvested during the planning horizon. Sixty percent of this supply requirement was arbitrarily assigned to demand node 1, and the rest to demand node 2. Demand levels for each period therefore resulted in 396,000, 594,000 and 792,000 m3 for demand node 1, and 264,000, 396,000 and 528,000 m3 for demand node 2. The annual discount rate was assumed 4%. The timber forecast was based on both the current volume estimate and on a growth and yield model that projected this volume to future periods. The error of the initial volume estimate, e1, was used to define two scenarios that might represent two sampling intensities. We assumed e1=10% in the first scenario (referred to as SCE10) and e1=15% in the second (referred to as SCE15), and in both cases we assumed that errors increased 1% annually. Assuming that calculations occur at the midpoint of each period, the errors for the second and third periods were 12% and 16% in SCE10, and 17.5% and 21% in SCE15, respectively. The volume of a stand s in period t was therefore assumed random and uniformly distributed in the range tstst evv ⋅± Γ as a result of unbiased errors in the forecast model or in the estimation of the stand area or the stand volume. As uncertainties have to be non-correlated in the same constraint (i.e. the volume must be independent within the same period), we also assumed spatial independence. As mentioned before, the degree of conservatism in satisfying a constraint is controlled by the user-defined protection level, . The immediate question is how big this parameter should be to get a desired feasibility rate. Although there is no exact expression to obtain this parameter, there exist bounds that relate a desired probability of constraint violation to the protection level required (e.g., Bertsimas et al. (2004) and Bertsimas and Sim (2004)). However, as noted in (Palma and Nelson 2009), these bounds represent only a weak estimate of this probability and therefore overestimate the protection level and excessively affect the objective function value. 90 We discarded the use of probability bounds in this work and determined protection levels as different percentages (e.g. 0.5%, 1.0% and 1.5%) of the number of uncertain coefficients of each constraint. These numbers provided a good description of the trade-off between robustness and optimality. For these three protection levels, we determined infeasibility rates by simulation as described next. 4.4.4 Simulation experiments Since we did not use probability bounds to determine protection levels, we had no idea about the infeasibility rates produced by the protection levels used. We therefore estimated these rates by simulation experiments. We simulated the volume for each scenario (SCE10 and SCE15) and used this volume in the deterministic and robust solutions to evaluate the performance of harvest decisions in satisfying demand constraints. That is, road decisions and the selection of stands to be harvested were fixed and the harvest levels were recomputed with the simulated volume coefficients. Since no maximum capacity on arcs was assumed, flow feasibility did not need to be checked. The total simulated timber production was compared with the total minimum demand and a simulation was considered infeasible if the production fell under the minimum requirement. As specific flow details are not required to determine whether a simulation was feasible or not, the flow variable was not recomputed. For each scenario, 1000 simulations were performed and the solution of the deterministic model was compared with the solutions of the robust models with the three protection levels used (labeled as PL0.5, PL1.0 and PL1.5). In addition, we used the buffer estimated by robust models (i.e. the value of the protection function) in a deterministic framework, that is, the deterministic model was run with modified demand levels (demand plus buffer). By doing this, we compared the effect of using a robust approach instead of a deterministic one with manually imposed buffers on the quality of harvest and road construction decisions. For each simulated scenario, we evaluated the performance of the deterministic decisions (DET), three robust decisions (ROB) and three deterministic solutions with manually imposed buffers (BUF). Simulations were performed in MS Excel in a 3.2 Ghz Pentium D with 1 GB of RAM. 91 4.5 Results Deterministic models were run for 30 min to obtain gaps around 2-3%. Robust models were harder to solve as more constraints and variables are included. In this case, models were run for 12 h to obtain similar gaps (3-4%). The possible effect of this gap difference is discussed in the Discussion section. Results of the optimization process are presented in Table 4.1. Table 4.1. Numerical results of the optimization process. Number of Problem Constraints Binary variables Total variables CPU time [min] Residual gap [%]* DET 4,025 2,532 5,017 30 2.4,2.4 ROB_PL0.5 5,318 2,532 6,313 720 3.1,3.2 ROB_PL1.0 5,318 2,532 6,313 720 3.3,3.4 ROB_PL1.5 5,318 2,532 6,313 720 3.7,4.4 BUF_PL0.5 4,025 2,532 5,017 30 2.6,2.4 BUF_PL1.0 4,025 2,532 5,017 30 2.9,2.6 BUF_PL1.5 4,025 2,532 5,017 30 2.5,2.7 * Residual gap is shown for the two scenarios of volume estimate precision (SCE10, SCE15) Simulation experiments showed that infeasibility rates of deterministic decisions were considerably reduced as higher protection levels were used to find robust decisions (Fig. 4.3). Infeasibility rates of around 49% observed in deterministic decisions (protection level 0) dropped to 1.4% and 1.5% for scenarios SCE10 and SCE15, respectively, when protection level was 1.5%. Decisions of deterministic models with modified demand levels (BUF) also showed lower infeasibility rates, although higher than those from robust models (3% in SCE10 and 2.3% in SCE15). The need for higher harvest levels than in the optimal of the deterministic solution caused a reduction in the objective function of robust and buffered solutions (Fig. 4.3). This reduction is slight and consistently lower for robust decisions than for the buffer solutions. 92 Fig. 4.3. Reductions in the percentage of infea were observed in both scenarios when higher showed higher infeasibility rates and higher l decisions. Robust models harvested more timber than th guarantee the fulfillment of minimum demand coefficients (Fig. 4.4). As expected, the highe production. Furthermore, the higher the volum production required for an equivalent protectiSCE10 SCE15 sible simulations and in the objective function protection levels were used. The buffer strategy osses in the objective function compared to robust e deterministic model as the former needed to constraints for different values of the volume r the protection level, the higher the timber e estimate error (SCE15) the higher the timber on level. 93 SCE10 SCE15 Fig. 4.4. Harvest levels of robust optimization models increased in both scenarios when a higher protection level (PL) was used. Lower timber harvest was required when the volume estimate error was smaller (SCE10). Road and harvest decisions of robust models were different from deterministic decisions (Fig. 4.5). For the sake of brevity, only decisions of scenario SCE10 are shown in Fig. 4.5. Decisions of SCE15 were similar. In the first planning period, 61 (72) stands and 98 (112) segment roads were scheduled differently in SCE10 (SCE15), respectively, when the higher protection level was used to obtain robust solutions. Decisions of the buffer strategy were also different from the traditional deterministic model because modified demand levels were used. Although the buffer size is the same as the one determined by the robust model, decisions of the buffer strategy tended to follow a similar spatial pattern as the deterministic decisions (Fig. 4.5). 94 Model: DET Model: ROB Model: BUF Fig. 4.5. Harvest and road building decisions were different among deterministic (DET and BUF) and robust models (ROB). Buffer strategy (BUF) tended to produce a similar decision pattern as the deterministic model (DET). 95 Because of the higher harvest levels required to get protection against infeasibility, both robust and buffer models increased the number of stands and the area harvested throughout the planning horizon (Table 4.2). Surprisingly, in both scenarios robust decisions generally required a smaller road network than deterministic and buffer strategy models. The buffer model required even more roads than deterministic solutions. No major changes in the costs were observed among the different models. However, ROB solutions consistently outperformed BUF solutions. Road building cost tended to be lower as more robustness was required. Table 4.2. Summary of harvest results for the planning models (SCE10 / SCE15) Model # Stands Total Area [ha] Roads [km] Average Cost [$/m3] Harvest Transport. Roads Det 218 4,804 108.7 2.85 1.77 2.55 Rob0.5 223 / 223 4,827 / 4,835 110.2 / 108.7 2.89 / 2.87 1.72 / 1.71 2.62 / 2.55 Rob1.0 224 / 226 4,911 / 4,885 106.7 / 110.0 2.86 / 2.88 1.69 / 1.68 2.50 / 2.62 Rob1.5 226 / 231 4,923 / 4,995 107.6 / 105.8 2.81 / 3.01 1.75 / 1.77 2.48 / 2.45 Buf0.5 223 / 219 4,835 / 4,821 112.9 / 110.9 2.88 / 2.85 1.82 / 1.78 2.59 / 2.58 Buf1.0 223 / 227 4,897 / 4,903 119.5 / 110.1 2.85 / 2.91 1.82 / 1.79 2.74 / 2.57 Buf1.5 224 / 228 4,868 / 4,976 113.1 / 116.0 2.90 / 2.88 1.76 / 1.77 2.59 / 2.62 4.6 Discussion Our results suggest that the robust optimization approach produces solutions that are less sensitive to random errors in the volume estimates at the expense of a slight reduction in the objective function. Even in the hypothetical case where the appropriate buffer is used to modify the minimum demand levels of deterministic models, the robust approach produced more stable solutions and lower losses in the objective function. Since meeting the minimum demand constraints is highly desirable in our case, harvest levels of robust solutions were higher than those of the deterministic solutions. As expected, the greater the need for constraint fulfillment, the higher the timber production. In addition, when volume coefficients are more uncertain (i.e. scenario SCE15) more timber harvest is needed to guarantee a certain level of constraint satisfaction. This higher harvest happens because the wider range of possible values increases the magnitude of eventual negative deviations against which we want 96 protection to satisfy the minimum demand levels. Although buffers can be used to manually modify the minimum demand levels, the robust approach estimates this buffer based on the decisions actually made as well as the degree of data uncertainty, providing optimal ‘stable’ decisions rather than the traditional optimal solutions. We showed this through the comparison between the robust approach and what we called the buffer strategy. When the same buffer obtained by robust models was used to modify demand levels of the deterministic model, robust solutions outperformed deterministic decisions both in terms of infeasibility rates and objective function value. Even though one might think that the buffer strategy performed almost as well as the robust approach, we note that finding the buffer levels by other means than the approach presented here is difficult. The amount of buffer obtained by robust models was “optimal” within the range of the solution gap, so the use of any other buffer level would likely produce inferior results. As integer models are unlikely to be solved to optimality, the presence of a gap in the solutions precludes an exact comparison of their quality. Obtaining the same gap for different models to standardize the comparison is also difficult as gaps progress by discrete steps rather than continuous moves. However, the greater gap of robust models might represent more room to improve their solutions and therefore enhances the benefits of robust formulations. The volume estimate error affected management decisions. Robust models scheduled both harvest and road construction in a different way than the deterministic model in order to take advantage of the uncertain timber yield through the planning horizon in a cost-efficient way. Although the robust models and the buffer strategy were set up with the same minimum level of production, their decisions largely differed. While the buffer strategy tended to follow a similar spatial pattern to the deterministic decisions (with the extra requirement to harvest more timber), the robust approach favored the concentration of the harvest operations. Although clustering operations may sound counter-intuitive if spatial diversification is desired (we will come back to this when discussing our assumptions), in our particular case this reduced the cost of road construction (Table 4.2) and therefore balanced the trade-off between optimality and robustness of the solutions. What should be clear is that robust models looked for different attributes when searching for optimal solutions, preferring decisions that are less sensitive to changing future conditions. However, we note that the magnitude of the decision changes should be scenario specific and influenced by the spatial distribution of uncertainties and the minimum demand levels. In addition, the discount rate will also affect the difference between robust and deterministic decisions. Higher discount rates are expected to reduce this difference as less 97 emphasis is put on the future and therefore in the magnitude of the uncertainties (Palma and Nelson 2009). For each constraint affected by uncertainty and for which feasibility is highly desirable a new term is added. This term, called the protection function, is an optimization problem in itself. For each solution vector, it looks for the buffer required to hold the constraint for a given user- defined protection level. Decisions that contribute to both higher objective value and small buffers are preferred, thereby providing high-quality solutions. Although the same protection levels were used in our analysis for all constraints, they can differ to express different degrees of importance. For instance, protection levels can emphasize the importance of the first period over the rest of the planning horizon, or the differences among products if a multi-product model is used. The importance of demand centers can be handled with the parameter δmt, that represents the fraction of the buffer to be sent to each destination point. Robust models have more variables and constraints than the original models. In our case, since the harvested volume in each period t comes from the potential harvest of any of the s stands (all stands are available and all have uncertain volume) s x t (431 x 3 = 1,293) new constraints need to be added. For each of these constraints, as well as for each period, a new variable has also to be created, that is, s x t + t (1293 + 3 = 1,296). Although the increase in the model size does not seem to be particularly important, robust models were harder to solve than the deterministic model as more protection levels were used. Despite the extended solution time used for robust models, residual gaps were greater than for deterministic models. This suggests that robust formulations are weaker than deterministic formulations and extra efforts should be done to get the same gaps of the deterministic approach (i.e., increased solution times, other strengthening and solution techniques). Our approach includes some important assumptions. For instance, uncertain coefficients need to be non-correlated in the same row, which means that volume coefficients of each constraint should be independently distributed. This forced us to assume that the stand volume is spatially independent. In our view, this is a critical limitation of the methodology to address spatial problems. As a consequence, clustering strategies are preferred and the value of the decisions is overestimated. It is important to note, however, that some degree of independency has to be assumed if we are looking for stable solutions to deal with uncertainty. If they are highly correlated then we cannot take advantage of diversification strategies and only worst-case 98 scenario solutions would become relevant. The assumption of independent coefficients also limits the type of uncertainty that we can consider to only random errors. Uncertain trends in data and catastrophic events like fires and pest attacks cannot be properly modeled with the current assumptions, hence further research is needed to address this limitation. Uncertain values are also assumed to be uniformly and symmetrically distributed, which can lead to the loss of information if uncertainties are well described by a different probability distribution. This could be handled by considering a uniform range that embraces, for instance, 99% of the observations of the original distribution of uncertain values (Palma and Nelson 2009). Although more conservative solutions will be obtained than when the true distribution is used, this simplification helps keep robust models easier to solve. Another disadvantage of the methodology is the impossibility to accurately determine protection levels required to hold specific infeasibility rates. Although probability bounds can be used to estimate them (Bertsimas and Sim 2004), they are loose and produce more conservative decisions than desired (Palma and Nelson 2009). As used here, simulation experiments can be used to evaluate infeasibility rates. Further research should address the limitations previously mentioned. Ways to include some levels of independency as well as asymmetry of the uncertain coefficients have recently been proposed (Chen et al. 2007) and their applicability should be explored. In addition, more accurate probability bounds to relate the probability of constraints violation and protection levels would improve this methodology by leaving out the need for simulation experiments. Finally, the application of other solution techniques (i.e. lagrangean relaxation, heuristics) to solve robust road building formulations would also be of interest in order to solve larger problems quickly and with reduced gaps. 4.7 Conclusions We presented a robust formulation of a road building and harvest scheduling problem with random volume coefficients when feasibility of minimum demand constraints is highly desired. Solutions of this model were both less sensitive to volume uncertainties and more efficient in terms of the objective value than manually imposed buffering strategies. The appropriate amount of buffer to satisfy a constraint is determined in conjunction with the decision variables in such a way that the trade-off between cost and robustness is explicitly considered in the formulation. Although larger models are obtained, the approach does not require sophisticated solution 99 techniques, and a commercial integer optimization software can be used. However, larger gaps were observed suggesting that more strengthening or different solution approaches may be needed. Since current probability bounds that relate protection levels and feasibility rates represent weak estimates of the probability of constraint satisfaction, simulation experiments were used to obtain infeasibility rates. The independency assumption among uncertain coefficients resulted in clustered spatial decisions, suggesting that further research is needed in this area. 100 4.8 References Andalaft, N., Andalaft, P., Guignard, M., Magendzo, A., Wainer, A., and Weintraub, A. 2003. A problem of forest harvesting and road building solved through model strengthening and Lagrangean relaxation. Oper. Res. 51: 613-628. Bare, B.B., and Mendoza, G.A. 1992. Timber harvest scheduling in a fuzzy decision environment. Canadian Journal of Forest Research 22: 423-428. Ben-Tal, A., and Nemirovski, A. 2002. Robust optimization - methodology and applications. Math. Program. 92: 453-480. Ben-Tal, A., and Nemirovski, A. 2000. Robust solutions of linear programming problems contaminated with uncertain data. Math. Program. 88: 411-424. Bertsimas, D., and Thiele, A. 2006. A robust optimization approach to inventory theory. Oper. Res. 54: 150-168. Bertsimas, D., and Sim, M. 2004. The price of robustness. Oper. Res. 52: 35-53. Bertsimas, D., and Sim, M. 2003. Robust discrete optimization and network flows. Math. Program. 98: 49-71. Bertsimas, D., Pachamanova, D., and Sim, M. 2004. Robust linear optimization under general norms. Oper. Res. Lett. 32: 510-516. Bevers, M. 2007. A chance constraint estimation approach to optimizing resource management under uncertainty. Canadian Journal of Forest Research-Revue Canadienne De Recherche Forestiere 37: 2270-2280. Birge, J.R., and Louveaux, F. 1997. Introduction to stochastic programming. First ed. Springer, New York. Boychuk, D., and Martell, D.L. 1996. A multistage stochastic programming model for sustainable forest-level timber supply under risk of fire. For. Sci. 42: 10-26. Boyland, M., Nelson, J., Bunnell, F.L., and D'Eon, R.G. 2006. An application of fuzzy set theory for seral-class constraints in forest planning models. For. Ecol. Manage. 223: 395-402. Chen, X., Sim, M., and Sun, P. 2007. A robust optimization perspective on stochastic programming. Oper. Res. 55: 1058-1071. Clark, M.M., Meller, R.D., and McDonald, T.P. 2000. A three-stage heuristic for harvest scheduling with access road network development. For. Sci. 46: 204-218. El Ghaoui, L., Oustry, F., and Lebret, H. 1998. Robust solutions to uncertain semidefinite programs. Siam Journal on Optimization 9: 33-52. 101 Ells, A., Bulte, E., and vanKooten, G.C. 1997. Uncertainty and forest land use allocation in British Columbia: Vague priorities and imprecise coefficients. For. Sci. 43: 509-520. Eriksson, L.O. 2006. Planning under uncertainty at the forest level: A systems approach. Scand. J. For. Res. 21: 111-117. Gassmann, H.I. 1989. Optimal harvest of a forest in the presence of uncertainty. Canadian Journal of Forest Research 19: 1267-1274. Guignard, M., Ryu, C., and Spielberg, K. 1998. Model tightening for integrated timber harvest and transportation planning. Eur. J. Oper. Res. 111: 448-460. Hof, J.G., Kent, B.M., and Pickens, J.B. 1992. Chance constraints and chance maximization with random yield coefficients in renewable resource optimization. For. Sci. 38: 305-323. Hof, J.G., Pickens, J.B., and Bartlett, E.T. 1986. A maxmin approach to nondeclining yield timber harvest scheduling problems. For. Sci. 32: 653-666. Hoganson, H.M., and Rose, D.W. 1987. A model for recognizing forestwide risk in timber management scheduling. For. Sci. 33: 268-282. Kall, P., and Wallace, S.W. 1994. Stochastic programming. Second ed. John Wiley & Sons, Chichester. Maness, T., and Farrell, R. 2004. A multi-objective scenario evaluation model for sustainable forest management using criteria and indicators. Canadian Journal of Forest Research 34: 2004-2017. Marshall, P.L. 1987. Sources of uncertainty in timber supply projections. For. Chron. 63: 165- 168. Mendoza, G.A., Bare, B.B., and Zhou, Z.H. 1993. A fuzzy multiple objective linear- programming approach to forest planning under uncertainty. Agricultural Systems 41: 257-274. Mowrer, H.T. 2000. Uncertainty in natural resource decision support systems: sources, interpretation, and importance. Comput. Electron. Agric. 27: 139-154. Mulvey, J.M., Vanderbei, R.J., and Zenios, S.A. 1995. Robust optimization of large-scale systems. Oper. Res. 43: 264-281. Olsson, L. 2007. Optimal upgrading of forest road networks: Scenario analysis vs. stochastic modelling. Forest Policy and Economics 9: 1071-1078. Ordonez, F., and Zhao, J.M. 2007. Robust capacity expansion of network flows. Networks 50: 136-145. Palma, C.D., and Nelson, J.D. 2009. A robust optimization approach protected harvest scheduling decisions against uncertainty. Canadian Journal of Forest Research 39: 342- 355. 102 Pickens, J.B., and Hof, J.G. 1991. Fuzzy goal programming in forestry - An application with special solution problems. Fuzzy Sets Syst. 39: 239-246. Pickens, J.B., and Dress, P.E. 1988. Use of stochastic production coefficients in linear- programming models - objective function-distribution, feasibility, and dual activities. For. Sci. 34: 574-591. Pickens, J.B., Hof, J.G., and Kent, B.M. 1991. Use of chance-constrained programming to account for stochastic variation in the A-matrix of large-scale linear programs: A forestry application. Ann Oper Res 31: 511-526. Regan, H.M., Colyvan, M., and Burgman, M.A. 2002. A taxonomy and treatment of uncertainty for ecology and conservation biology. Ecol. Appl. 12: 618-628. Richards, E.W., and Gunn, E.A. 2000. A model and tabu search method to optimize stand harvest and road construction schedules. For. Sci. 46: 188-203. Stirn, L.Z. 2006. Integrating the fuzzy analytic hierarchy process with dynamic programming approach for determining the optimal forest management decisions. Ecol. Model. 194: 296-305. Weintraub, A., and Abramovich, A. 1995. Analysis of uncertainty of future timber yields in forest management. For. Sci. 41: 217-234. Weintraub, A., and Vera, J. 1991. A cutting plane approach for chance constrained linear- programs. Oper. Res. 39: 776-785. Weintraub, A., and Navon, D. 1976. Forest Management Planning Model Integrating Silvicultural and Transportation Activities. Management Science 22: 1299-1309. Weintraub, A., Jones, G., Meacham, M., Magendzo, A., Magendzo, A., and Malchuk, D. 1995. Heuristic Procedures for Solving Mixed-Integer Harvest Scheduling - Transportation- Planning Models. Canadian Journal of Forest Research 25: 1618-1626. Zimmermann, H.-. 1996. Fuzzy set theory and its applications. Third edition ed. Kluwer academic, Boston, USA. 103 5 Conclusions 5.1 General conclusions The robust optimization (RO) approach proposed in this thesis produced solutions that were less sensitive to uncertain data and therefore more stable with respect to the changing future conditions that forest managers face. Surprisingly, robust decisions produced only minor reductions in the objective function suggesting that choosing more stable decisions in some situations might not represent a major cost for the decision maker. The methodology proved to be a useful way to explicitly consider some types of uncertainty in forest management models as discussed later. To meet the thesis objectives (i.e., to determine how the RO approach performs relative to the deterministic approach, and to determine its potential benefits, limitations and implications), three different problems were considered: (a) a strategic, non-spatial harvest scheduling problem with uncertain volume yields and demands (chapter 2), (b) a strategic, bi-objective multi-period problem with uncertain preferences (chapter 3), and (c) a tactical, spatial harvest scheduling and road building problem with uncertain volume yields (chapter 4). Deterministic and robust models to solve these problems were formulated and their solutions were compared under simulated scenarios of uncertain data. This comparison included changes in the objective function, in the rates of constraint satisfaction and in the management decisions for different degrees of robustness. In (a), minimum demand and timber flow fluctuation constraints were protected against infeasibility. The robust formulation increased timber production levels in critical periods in which demand constraints were tightly met. In periods with high negative fluctuation, timber production was smoothed by moving it from periods with a volume surplus to periods where only the necessary volume to meet constraints was produced. Although the use of buffer strategies as a response to risk is not new (e.g., Boychuk and Martell 1996), RO allows decision makers to estimate this buffer in conjunction with the decision variables and based on the degree of uncertainty of the data. It was shown that even when the same optimal buffer obtained by 104 robust models was used to modify the constraint level of a deterministic model, solutions were not necessarily robust and higher infeasibility rates were obtained. Higher discount rates diminished the effect of uncertainty on management decisions as suggested by Gassmann (1989). It is worth noting that although a Model I formulation was used to solve this problem, a Model II formulation might be more convenient as the problem becomes bigger in terms of the number of strata and management options. In this case, the number of variables of a Model I formulation will be probably larger than the number of variables observed in a Model II. This will translate into more constraints in its robust counterpart therefore making a Model II formulation more attractive as base for a robust model. In the context of multiple objectives decision-making, the difficulties in estimating the real value of the objective weights (Eckenrode 1965, Cohon 1978, Kangas 1994) together with the changing communities’ demand for outputs from the forest (Mowrer 2000, Kangas and Kangas 2004) have motivated the study of methodologies to incorporate this uncertainty into multi- objective decision models. In (b), uncertain weights of the two objectives were assumed, and unlike most of the research in the area that evaluates the robustness of sets of discrete decisions alternatives previously known (Lahdelma et al. 1998, Kangas et al. 2003, Kangas 2006), I addressed weight uncertainty in a continuous decision problem where decisions were generated rather than evaluated. The weighted sum of objectives produced by robust decisions was more stable than deterministic decisions along the planning horizon when evaluated under uncertain weights. This means that the robust decisions produced a good and consistent combination of objectives independently of the relative importance that society gives to them. It was shown that first-period management decisions did not change in scenarios with low levels of uncertainty, suggesting that uncertain weights might be overcome in some cases by a continuous planning process that keeps track of the change in social preferences. To compensate for the reduction in the objective function (weighted sum of objectives), robust decisions produced higher levels of the most important objective as more robustness was required. In (c), as in (a), the feasibility of minimum demand constraints was desirable and a slight increase in the timber production levels was observed with consequent reduction in the infeasibility rates. Since volume yield was assumed spatially independent as I discuss later, harvest decisions tended to be clustered instead of scattered as one may expect, suggesting that the methodology needs revision when spatial decisions are to be made. Unlike the modest increase in the solution efforts observed in models (a) and (b), a significant increase in the 105 difficulty to solve model (c) was observed. The strength of integer models is crucial to successfully solve integer problems (Guignard et al. 1998), and my results suggest that robust formulations of integer problems may be weaker in terms of solvability than deterministic formulations. Solving robust formulations of integer problems might therefore require extra efforts such as increased solution times and other strengthening and solution techniques if the same solution quality as the deterministic version is desired. Unlike some studies that have focused on the effect of uncertainty on the objective function and harvest levels (Boychuk and Martell 1996, Gassmann 1989, Pickens et al. 1991, Teeter and Caulfield 1991, Eriksson 2006), I also explored its effect on the specific harvest scheduling decisions as in Hoganson and Rose (1987) and Hof et al. (1995). The need for considering uncertainty in the planning process is irrelevant for management purposes if decisions remain the same in both cases. I demonstrated that explicitly considering uncertainty affected forest management decisions right at the beginning of the planning horizon, that is, even decisions to be implemented in the current period are different if uncertainty is taken into account. This suggests that not considering uncertainty may require major corrections to future management decisions through time in order to adapt to changing future conditions. The more corrections needed, the greater the reduction in the value of the management plan. By providing decisions that are less sensitive to uncertain data, robust solutions reduce the need of future modifications to these decisions and therefore reduce the cost of adapting management to realized outcomes. Solutions to robust models represent diversification strategies to reduce the effect of uncertainty as described in the 1950’s (Heady 1957, Anoff 1959). In all the evaluated cases these solutions increased the number of variables taking a non-zero value, thereby increasing the original set of selected actions. In Chapter 2 and 3 (especially Chapter 2) this translated into more management prescriptions applied to each stratum. For instance, in the harvest scheduling problem of Chapter 2 the average number of prescription applied to a single stratum increased from 1.03 to 1.9 in the most severe scenario of uncertainty with a 4% discount rate. In some cases a stratum was prescribed up to 14 management prescriptions. The increase in the number of prescriptions observed in the bi-objective problem (Chapter 3) was less significant and only visible at the most severe scenarios of uncertainty, probably because this problem was less constrained and did not include any specific output level requirement. In Chapter 4, the diversification effect led the robust model to schedule more, but slightly smaller, stands for harvesting. This effect is expected to be more evident as spatially correlated uncertainties are considered. Environmental and 106 practical implications of the change in the nature of the decisions will depend on the problem under consideration. No major implications should be observed in strategic decision problems as they basically deal with long-term trends and non-spatial aggregated information (Nelson 2003). However, environmental and practical implications might be observed in medium and short term decisions that are spatially located as access to more stands needs to be available. In this case, environmental issues should be explicitly incorporated into the mathematical formulation of the problem. 5.2 Strengths and weaknesses of the research The RO approach has important assumptions that determine its strengths and weaknesses and the appropriateness of its application in the forest management context. Uncertainties have to be independent, uniform and symmetrically distributed around the parameter estimate. The assumptions of uniformity and symmetry may represent a limitation of the methodology when the nature of the uncertainties is well described by a different distribution or when extra information other than the extreme points is available. In this case, one can use a wide enough uniform range of uncertain values that contains a high proportion of observations of the original distribution. These assumptions could also be seen as an advantage of this methodology over the most traditional methodologies used to account for uncertainty. For instance, stochastic and chance-constrained programming (Birge and Louveaux 1997) require the probability of occurrence of each scenario and the full description of the probability distributions of uncertainties, respectively. However, this information is not usually available (Regan et al. 2005). The independency assumption is an important weakness of the methodology as decisions in the forest management context are usually spatial and temporally correlated. This assumption requires that uncertain coefficients in the same constraint are independent, which implies that the outcomes from different stands in a given constraint (usually a fixed period) are non-correlated. If a stand appears in a constraint more than once (e.g. if thinning is also considered then one management alternative can harvest part of the stand and another alternative can thin the other part in the same period), then a certain degree of dependency should be expected, and optimistic results will be obtained. The independency assumption may be overcome when dealing with non-spatial decisions where stands from different locations are usually grouped into non-spatial 107 entities like in Chapters 2 and 3. In this case, these new entities tend to cancel the local correlation among stands. However, when dealing with spatial decisions such as in Chapter 4 the independency assumption becomes critical. Instead of promoting spatial diversification to offset uncertainty realizations, harvest decisions tend to be spatially clustered to reduce the cost of the higher level of timber production required by robust solutions. This suggests that the RO approach as presented here should be revised if spatial decisions need to be addressed. This is clearly a drawback of the methodology as uncertainty correlation can be more easily included in stochastic or chance-constrained programming. Another weakness of the approach presented is the lack of a clear relationship between the protection level used to represent different degrees of robustness and the resulting infeasibility rates. Current probability bounds (see for example Ben-Tal and Nemirovski 2000, Bertsimas and Sim 2004) are loose and therefore produce more conservative solutions than expected. Although these bounds can still be used to determine the protection level required to guarantee an a priori probability of constraint violation, simulation experiments might be needed to avoid the excessive impact on the objective function produced by over conservative solutions. More precise infeasibility rates may be obtained using chance-constrained programming. The simplifying assumptions described above lead to what is in my view the most important strength of this methodology over the existing ones. While considering a continuous set of scenarios, the RO approach presented in this thesis keeps the models computationally tractable, therefore only requiring common solutions techniques used for solving deterministic models. Other methodologies that account for uncertainty (i.e., stochastic and chance-constrained programming), impose serious difficulties when solving the resulting models (Birge and Louveaux 1997). In the first case, the number of scenarios may quickly become extremely large requiring specific solution methods. In the second case, models usually become non-linear also requiring specific non-linear solution algorithms. With RO, modified versions of traditional deterministic problems can be built and solved by commercial linear optimization software, and then solutions can be tested by simulating uncertain coefficients without the need to solve complex problems with sophisticated solution techniques. I believe this approach might allow decision makers interested in reducing the risk of their decisions to objectively determine the amount of buffer required to guarantee critical requirements. As good solutions for a set of model parameters are obtained, this methodology may significantly reduce the efforts invested in conducting the sensitivity analyses required by deterministic models. 108 5.3 Future research and potential applications of this research Future research should address the impacts of violating the assumptions previously described. The most important one is the assumption of independency of the uncertain coefficients, which limits the application of this approach to random, unbiased and uncorrelated errors in the coefficients estimates. Although this kind of error is important, especially when forecast models are used, catastrophic events (e.g., fire occurrence, pest attacks) and uncertain trends in data (e.g., price, systematic errors) will only be properly modeled if temporal dependence is allowed. In tactical and operational decision-making spatial correlation plays a key role, so ways to include it in a way that does not affect the computational tractability of the approach is highly relevant. It may also be important to explore approaches to consider other distributions (e.g., triangular) and asymmetry of uncertain ranges. Although a full description of the uncertainty probability distributions is unlikely to be available, an approach that considers information other than the point estimate while retaining the modeling simplicity would be an important contribution. Another disadvantage of the methodology that should be addressed in future research is the need to perform simulation experiments to determine the feasibility rates of robust models. Current probability bounds that determine feasibility rates a priori have been proposed (Bertsimas and Sim 2004) but they are loose and produce conservative solutions (Palma and Nelson 2009). Tighter probability bounds might reduce the need for simulation to estimate the resulting feasibility rates and therefore make this approach more appealing. The approach proposed here considers a static decision structure as opposed to a dynamic structure observed in most multi-period decision problems. In other words, real decisions are made sequentially as part of the uncertainty is observed, just like the stochastic programming approach is implemented (see section 1.2.1). Although the RO approach as presented here has been applied to dynamic multi-period decision problems (Bertsimas and Thiele 2006, Bertsimas and Pachamanova 2008, Chung et al. 2009), a way to formally represent the recourse structure of the decision-making process would contribute to a better conceptual justification of the approach. The application of this methodology to other important forest resources management problems would also be of interest. In particular, decisions that involve a higher degree of irreversibility, like protection of endangered species or habitat conservation, would be benefit from using this approach as they are currently made in the presence of many uncertain parameters (Moilanen et 109 al. 2006, Nicholson and Possingham 2007). In these cases, irreversible decisions might be infinitely costly to reverse (e.g. extinction of a species) and represent the loss of future management options (Arrow and Fisher 1974). Less sensitive solutions to uncertainty, like the ones produced by RO, may reduce the chance of undesirable outcomes. Another possible application of this method is to explore the effect of data quality on management decisions (e.g., the quality of current and projected volume yields or other forest outcomes). The trade-off between the cost of producing more accurate data (more intense sample methods, more complex forecast models) and the improvement in the value of the decisions can be determined by modifying the range of data uncertainty used in RO. Finally, the implementation of RO in optimization-based forest planning software currently used by forest managers, e.g. Spectrum (USDA Forest Service, 1999), MELA (Siitonen et al. 1996), FOLPI (García 1984), Woodstock (Remsoft 2006) and others, would make it easier to extend the application of this methodology to real management decision making. 110 5.4 References Ansoff, H.I. 1959. Strategies for diversification. Harvard Bus Rev. 35(5): 113-124. Arrow, K.J., and Fisher, A.C. 1974. Environmental preservation, uncertainty, and irreversibility. Quart. J. Econ. 88(2): 312-319. Ben-Tal, A., Nemirovski, A. 2000. Robust solutions of linear programming problems contaminated with uncertain data. Math. Programming 88: 411-424 Bertsimas, D., and Pachamanova, D. 2008. Robust multiperiod portfolio management in the presence of transaction costs. Comput. Oper. Res. 35: 3-17. Bertsimas, D., and Sim, M. 2004. The price of robustness. Oper. Res. 52: 35-53. Bertsimas, D., and Thiele, A. 2006. A robust optimization approach to inventory theory. Oper. Res. 54: 150-168. Birge, J.R., and Louveaux, F. 1997. Introduction to stochastic programming. First ed. Springer, New York. Boychuk, D., and Martell, D.L. 1996. A multistage stochastic programming model for sustainable forest-level timber supply under risk of fire. For. Sci. 42: 10-26. Chung, G., Lansey, K., Bayraksan, G. 2009. Reliable water supply system design under uncertainty. Environ. Modell. Softw. 24: 449-462. Eriksson, L.O. 2006. Planning under uncertainty at the forest level: A systems approach. Scand. J. For. Res. 21: 111-117. García, O. 1984. FOLPI, a forestry-oriented Linear Programming interpreter. In Proceedings IUFRO Symposium on Forest Management Planning and Managerial Economics. Japan. Edited by University of Tokyo. pp. 293-305. Gassmann, H.I. 1989. Optimal harvest of a forest in the presence of uncertainty. Canadian Journal of Forest Research 19: 1267-1274. Guignard, M., Ryu, C., and Spielberg, K. 1998. Model tightening for integrated timber harvest and transportation planning. Eur. J. Oper. Res. 111: 448-460. Heady, E.O. 1957. Diversification in resource allocation and minimization of income variability. J. Farm Econ. 34(4): 482-496 Hof, J., Bevers, M., and Pickens, J. 1995. Pragmatic approaches to optimization with random yield coefficients. For. Sci. 41: 501-512. Hoganson, H.M., and Rose, D.W. 1987. A model for recognizing forestwide risk in timber management scheduling. For. Sci. 33: 268-282. 111 Kangas, J. 1994. An approach to public-participation in strategic forest management planning. For Ecol Manage 70: 75-88. Kangas, A. 2006. The risk of decision making with incomplete criteria weight information. Can J for Res 36: 195-205. Kangas, J., Kurttila, M., Kajanus, M., and Kangas, A. 2003. Evaluating the management strategies of a forestland estate - the S-O-S approach. J Environ Manage 69: 349-358. Kangas, A.S., and Kangas, J. 2004. Probability, possibility and evidence: approaches to consider risk and uncertainty in forestry decision analysis. For. Policy Econ. 6: 169-188. Lahdelma, R., Hokkanen, J., and Salminen, P. 1998. SMAA - Stochastic multiobjective acceptability analysis. Eur J Oper Res 106: 137-143. Moilanen, A., Runge, M.C., Elith, J., Tyre, A., Carmel, Y., Fegraus, E., Wintle, B.A., Burgman, M., and Ben-Haim, Y. 2006. Planning for robust reserve networks using uncertainty analysis. Ecological Modelling 199: 115-124. Mowrer, H.T. 2000. Uncertainty in natural resource decision support systems: sources, interpretation, and importance. Comput. Electron. Agric. 27: 139-154. Nelson, J. 2003. Forest-level models and challenges for their successful application. Can. J. For. Res. 33: 422-429. Nicholson, E., and Possingham, H.P. 2007. Making conservation decisions under uncertainty for the persistence of multiple species. Ecol. Appl. 17: 251-265. Palma, C.D., and Nelson, J.D. 2009. A robust optimization approach protected harvest scheduling decisions against uncertainty. Canadian Journal of Forest Research 39: 342- 355. Pickens, J.B., Hof, J.G., and Kent, B.M. 1991. Use of chance-constrained programming to account for stochastic variation in the A-matrix of large-scale linear programs: A forestry application. Ann Oper Res 31: 511-526. Regan, H.M., Ben-Haim, Y., Langford, B., Wilson, W.G., Lundberg, P., Andelman, S.J., and Burgman, M.A. 2005. Robust decision-making under severe uncertainty for conservation management. Ecol. Appl. 15: 1471-1477. REMSOFT. 2006. Woodstock v2006.8. Modeling references. Remsoft Inc. Fredericton, NB, Canada. 104 pp. Siitonen M, Härkönen K, Hirvelä H, Jämsä J, Kilpeläinen H, Salminen O, Teuri M (1996) MELA Handbook – 1996 Edition, The Finnish Forest Research Institute, Research Papers 622, p 452. Teeter, L.D., and Caulfield, J.P. 1991. Stand density management strategies under risk - effects of stochastic prices. Canadian Journal of Forest Research 21: 1373-1379. 112 USDA Forest Service, 1999 USDA Forest Service, 1999. Spectrum Overview. Inventory and Monitoring Institute, Fort Collins, CO. Available from http://www.fs.fed.us/institute/planning_center/files/Spectrum26_Overview.pdf (accessed May 13, 2009). 113 Appendices Appendix A This appendix describes in detail the way the robust formulation of the model constraints were obtained. Protection against minimum demand infeasibility Uncertain volume: a protection function was included in the minimum demand constraint as follows: tpdxvA ptptpt ji ijijpt ,),x(]1[ , ∀≥Γ−∑ β D * where v represents the coefficient estimate of the volume. Note that (superscript comes from demand) is subtracted from the volume as this is a ‘≥’ inequality. For all p and t, is equivalent to the following problem: D D DΓ≤ D v̂ ptβ ptβ ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ ∑ ∈ DptJji ijptijijpt rxvA ),( *ˆmax]2[ subject to pt Jji ijpt D pt rA ∑ ∈),( ]3[ ptijpt JjirA ∈∀≤≤ ),(10]4[ where rijpt is the decision variables, is the precision range of the volume and JDpt is the set of uncertain coefficients for the demand constraint pt. The dual of model [A2]-[A4] is then: 114 ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ +Γ ∑ ∈ DptJji ijptpt D pt wuA ),( min]5[ subject to D* 0]7[ ≥uA D D ptijijptijptpt JjixvwuA ∈∀≥+ ),(ˆ]6[ pt ptijpt JjiwA ∈∀≥ ),(0]8[ where upt and wijpt are the dual variables of equations [A3] and [A4], respectively. The value of can be then replaced with the objective value of the model [A5]-[A8] to obtain the final robust formulation of [A1]: ptβ tpdwuxvA pt Jji ijptptpt ji ijijpt D pt ,]9[ ),(, ∀≥−Γ− ∑∑ ∈ D D tpuA ,0]11[ ∀≥ D tpJjixvwuA ptijijptijptpt ,,),(ˆ]10[ ∈∀≥+ pt tpJjiwA ptijpt ,,),(0]12[ ∈∀≥ Uncertain volume and demand: the following is the demand constraint with the protection function included: tpdxvA ptptpt ji ijijpt ,0),x(]13[ , ∀≥Γ−−∑ β D * where v and d represent the coefficient estimate of volume and demand. Considering that the demand coefficient is always present, for all p and t, is equivalent to the following problem: Dptβ ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ +∑ ∈ ptpt Jji ijptijijpt sdrxvA D pt ˆˆmax]14[ ),( * 115 subject to DΓ≤+ D 1]17[ =sA ˆ *ˆ D D D D* ≥uA D ptpt Jji ijpt srA D pt ∑ ∈),( ]15[ ptijpt JjirA ∈∀≤≤ ),(10]16[ pt where rijpt and spt are the decision variables, and and are the precision ranges of volume and demand. Because of equality [A17] the model [A14]-[A17] can be reformulated as follows: v̂ d ∑ ∈ + D ptJji ijptijijptpt rxvdA ),( ˆmax]18[ subject to 1]19[ ),( −Γ≤∑ ∈ pt Jji ijpt D pt rA ptijpt JjirA ∈∀≤≤ ),(10]20[ The dual of [A18]-[A20], which substitutes in [A13], is then: ptβ ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ +−Γ+ ∑ ∈ DptJji ijptpt D ptpt wudA ),( )1(minˆ]21[ subject to ptijijptijptpt JjixvwuA ∈∀≥+ ),(ˆ]22[ 0]23[ pt ptijpt JjiwA ∈∀≥ ),(0]24[ 116 where upt and wijpt are the dual variables of equations [A19] and [A20], respectively. The value of can be then replaced with the objective value of the model [A21]-[A24] to obtain the final robust formulation of [A13]: D ptβ tpdwudxvA pt Jji ijptptptpt ji ijijpt D pt ,)1(]25[ ),(, ∀≥−−Γ−− Dˆ ∑∑ ∈ tpJjixvwuA ptijijptijptpt ,,),(]26[ ∈∀≥+ Dˆ tpuA ,0]27[ ∀≥ D pt tpJjiwA ptijpt ,,),(0]28[ ∈∀≥ Protection against infeasibility of production fluctuation A protection function was included in the fluctuation of production constraint as follows: 1,)1(),x(]29[ , 1 , >∀−≥Γ− * ∑∑ − tpxvxvA ij ji ijptptptij ji ijpt δβ FF F Considering JFpt as the set of uncertain coefficients of the fluctuation constraint pt, for all p and t>1, the equivalent linear problem of is: ptβ ( )∑ ∈ −−+ F ptJji ijptijtijpijpt qxvvMaxA ),( )1(ˆ)1(ˆ]30[ δ * F F subject to pt Jji ijpt F pt qA Γ≤∑ ∈),( ]31[ ptijpt JjiqA ∈∀≤≤ ),(10]32[ where qijpt is the decision variable. If ypt and zijpt are the dual variables of equations [A31] and [A32], respectively, the dual problem of [A30]-[A32] for all p and t>1 is: 117 ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ +Γ ∑ ∈ FptJji ijptpt F pt zyMinA ),( ]33[ subject to ( ) F* 0]35[ ≥yA F ptijtijpijptijptpt JjixvvzyA ∈∀−+≥+ − ),(ˆ)1(ˆ]34[ )1(δ pt ptijpt JjizA ∈∀≥ ),(0]36[ Replacing the protection function in [A29] with the objective function of [A33]-[A36], the robust counterpart of [A29] is: 1,)1(]37[ , 1 ),(, >∀−≥−Γ− ∑∑∑ − ∈ tpxvzyxvA ij ji ijpt Jji ijptptptij ji ijpt F pt δF ( ) 1,,),(ˆ)1(ˆ]38[ )1( >∈∀−+≥+ − tpJjixvvzyA Fptijtijpijptijptpt δ 1,0]39[ >∀≥ tpyA F pt 1,,),(0]40[ >∈∀≥ tpJjizA ptijpt 118
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Robust optimization for forest resources decision-making...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Robust optimization for forest resources decision-making under uncertainty Palma, Cristian Dereck 2010
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Robust optimization for forest resources decision-making under uncertainty |
Creator |
Palma, Cristian Dereck |
Publisher | University of British Columbia |
Date Issued | 2010 |
Description | There is a general consensus that management decisions concerning forest resources are made in an intrinsically uncertain environment. However, decision-making tools used in forest management assume perfect information, leaving decision-makers to explore the most likely scenarios of uncertainty and determine the most reasonable management alternative. Although techniques that explicitly consider uncertainty exist, they increase the complexity of the models hence precluding their application to large-size problems. This dissertation describes the application of robust optimization concepts that explicitly consider uncertainty in forest management problems while keeping the models computationally tractable. By introducing some simplifying assumptions about uncertainty distributions, i.e. independency and uniformity, this approach allows for including uncertainty in many coefficients of the model. The methodology modifies the constraints for which feasibility is desirable and incorporates uncertainty in the technical coefficients by introducing an additional term. This term is an optimization problem in itself that introduces new constraints into the original model and acts as a buffer that guarantees constraint satisfaction for different uncertainty realizations. By changing the value of a robustness parameter, the trade-off between cost and robustness can be analyzed. The performance of this approach is explored through three structurally different problems: (a) a non-spatial harvest scheduling problem with uncertain volume yields and demands, (b) a multi-objective problem with uncertain preferences, and (c) a spatial harvest scheduling and road building problem with uncertain volume yields. Deterministic and robust formulations of these problems are provided and the performance of their solutions is evaluated under simulated scenarios of uncertain data. In all cases, robust decisions are less sensitive to uncertain data and hence protected from the occurrence of infeasibilities, with a modest reduction in the objective function value. Moreover, deterministic and robust decisions greatly differ, suggesting that traditional solutions may require major corrections to adapt to changing future conditions with a consequent decrease in the quality of the decisions. The effect of the methodology assumptions are discussed and future work is suggested. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-04-08 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivs 3.0 Unported |
DOI | 10.14288/1.0070920 |
URI | http://hdl.handle.net/2429/23329 |
Degree |
Doctor of Philosophy - PhD |
Program |
Forestry |
Affiliation |
Forestry, Faculty of |
Degree Grantor | University of British Columbia |
GraduationDate | 2010-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/3.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2010_spring_palma_cristian.pdf [ 2.48MB ]
- Metadata
- JSON: 24-1.0070920.json
- JSON-LD: 24-1.0070920-ld.json
- RDF/XML (Pretty): 24-1.0070920-rdf.xml
- RDF/JSON: 24-1.0070920-rdf.json
- Turtle: 24-1.0070920-turtle.txt
- N-Triples: 24-1.0070920-rdf-ntriples.txt
- Original Record: 24-1.0070920-source.json
- Full Text
- 24-1.0070920-fulltext.txt
- Citation
- 24-1.0070920.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0070920/manifest