A Progressive Hedging Approach to Solve Harvest Scheduling Problem under Climate Change

Abstract

Due to the long time horizon typically characterizing forest planning, uncertainty plays an important role when developing forest management plans. Especially important is the uncertainty related to recently human-induced global warming since it has a clear impact on forest capacity to contribute to biogenic and anthropogenic ecosystem services. If the forest manager ignores uncertainty, the resulting forest management plan may be sub-optimal, in the best case. This paper presents a methodology to incorporate uncertainty due to climate change into forest management planning. Specifically, this paper addresses the problem of harvest planning, i.e., defining which stands are to be cut in each planning period in order to maximize expected net revenues, considering several climate change scenarios. This study develops a solution approach for a planning problem for a eucalyptus forest with 1000 stands located in central Portugal where expected future conditions are anticipated by considering a set of climate scenarios. The model including all the constraints that link all the scenarios and spatial adjacency constraints leads to a very large problem that can only be solved by decomposing it into scenarios. For this purpose, we solve the problem using Progressive Hedging (PH) algorithm, which decomposes the problem into scenario sub-problems easier to solve. To analyze the performance of PH versus the use of the extensive form (EF), we solve several instances of the original problem using both approaches. Results show that PH outperforms the EF in both solving time and final optimality gap. In addition, the use of PH allows to solve the most difficult problems while the commercial solvers are not able to solve the EF. The approach presented allows the planner to develop more robust management plans that incorporate the uncertainty due to climate change in their plans.

1. Introduction

Forest management is one of the fields characterized by future uncertainty. Typical uncertainties may be related to market conditions such as timber prices and demand bounds [1] as well as natural factors such as climate, forest fires, droughts, insects, and diseases. Human-induced global warming is another source of uncertainty since climate directly affects forest growth and productivity through the photosynthesis process. In addition, forest growth is also affected by water availability and temperature, which are uncertain factors due to climate change [2]. Christensen et al. [3] point to global warming, which implies increases in temperature both in winter and summer as well as longer and warmer growing seasons and the occurrence of droughts in some regions. The weather conditions influence forest yields, and thus, the timber supply may vary depending on climatic scenarios. Therefore, consideration of climate change should lead to adjustments in forest plans. Under the stochastic nature of future conditions, decision makers are thus challenged to find robust solutions.

One of the most typical planning problems is the so-called harvest scheduling problem. The problem is related to which stands are going to be harvested in each period. This is a well-known problem, and a deterministic model is presented in Davis et al. [4]. Considering explicitly climate change and, consequently, variations in tree growth in the forest planning process (e.g., using stochastic programming) enables the planner to make more robust decisions based on the range of climate scenarios rather than using a single average scenario. This calls for the use of techniques that allow to address risk and uncertainty.

There are many examples of studies focusing on the management of uncertainty in forest planning (see [5,6,7], for thorough reviews). Some of them have focused on the use of chance-constrained to satisfy timber demand levels [8]. Other authors used robust optimization to ensure satisfying timber demands constraints subject to a certain protection level [9]. Stochastic goal programming has been also applied to develop forest management plans that minimize the expected deviations from the goals/targets set for different criteria defined by the decision maker [10]. Kangas et al. [11] used a two-stage stochastic programming method to define the optimal measurement strategy to maximize revenues of the derived management plan at forest level. An advantage of stochastic programming is that it deals with the uncertainty of the problem, and therefore, it allows managing risk [12]. For example, approaches as conditional value at risk [13] could be utilized for risk management.

Different approaches have been used to consider explicitly climate change in forest planning [5,6]. One of the most widely used techniques to explicitly consider risk and uncertainty is the use of scenarios [14]. There are many examples of the use of scenario analysis to deal with uncertainty in climate change [15,16,17,18]. In addition, decision support systems that incorporated optimization routines to generate optimized management plans at forest level considering different climate scenarios have been recently developed [19,20]. Nevertheless, these two DSSs used deterministic approaches (i.e., optimizing forest management for each scenario), which means that each climate change scenario was assumed to occur with certainty, and its stochasticity was not addressed. Eriksson [21] presented how stochastic programming may be used to address climate change uncertainty. In stochastic optimization, uncertainty in a dynamic on time a problem can be modeled using a scenario tree where each scenario represents a realization of the stochastic parameters. Stochastic programming takes all the scenarios into account without being subordinated to any one of them allowing developing policies that lead to solutions with a better expected objective value and that can be feasible for all scenarios. Few studies have focused on the use of stochastic programming to address climate change uncertainty [22,23]. Garcia-Gonzalo et al. [22] developed a stochastic linear programming model in a harvest scheduling problem to demonstrate the advantage of using a stochastic versus a deterministic model based on expected values. Álvarez-Miranda et al. [23] presented a stochastic goal programming model to deal with climate change uncertainty and included the use of conditional value at risk.

In this paper, a spatial large-scale multi-stage stochastic harvest scheduling problem considering uncertainty due to climate change was solved. In order to address climate change uncertainty when developing the harvest scheduling model studied in this paper, the 32 climate scenarios presented by Garcia-Gonzalo et al. [22] were considered for tree growth. Even though Garcia-Gonzalo et al. 2016 [4] used the same case study area and the same number of scenarios as in the present study; they simplified the problem (i.e., reduced the size of the mode by aggregating 1000 stands into 24 macro-stands (i.e., strata). The aggregation implied that the model did not take into account spatial considerations (e.g., exact location of individual stands), nor adjacency considerations (i.e., total adjacent area harvested within the same period or between periods must be smaller than a maximum size), which are important to reduce impacts of forest operations on landscape visual quality as well as to limit environmental impacts. The model aggregated stands into strata having the same rotation cycle, tree age, and soil quality. In the present paper, the original 1000 stands were considered, where the spatial definition is restored. Adjacency restriction constraints, limiting the maximum size of a contiguous harvested area, were considered. Adjacency constraints lead to much more difficult to solve mixed integer problems [24].

The presented harvest scheduling problem considered a planning horizon of 15 harvesting years and 32 scenarios. The model involved ensuring a feasible solution under each scenario and in addition the so-called non-anticipativity constraints [25]. These non-anticipativity constraints rule state that if any scenarios i and j are identical up to period t, then decisions taken under those scenarios up to period t must be equal as well, as they share the same information (for more details, see Birge and Louveaux [12] and Rockafeller and Wets [14].

The inclusion of both non-anticipativity and adjacency constraints lead to a large and difficult to solve mixed integer linear (MIP) problem. The extended formulation of this model may not be solvable. In this case, algorithms that use some form of decomposition have been proposed. One option is (PH) [14], which has been successfully applied in forest management problems including road construction [7]. Another option may be coordinating branching [26,27]. The use of PH was chosen since it has been proved to have a good performance for MIP problems (even when used as a heuristic) [7,28,29,30,31]. In addition, there is a free software, Pyomo Stochastic Programming, which implements an efficient version of PH with parallel optimization, which has also been successfully used in forestry problems [7]. The PH algorithm is based on horizontal decomposition [32], which relies on decomposing the problem into scenarios. This is done by relaxing the non-anticipativity constraints, which allow solving each scenario problem independently. The performance of the PH approach was compared to the use of the extended formulation to analyze its performance. For this purpose, several experimental instances of the original problem (in terms of forest sizes and timber demand) were solved. Different applications of PH algorithm in linear and mixed integer problems (where the algorithm is used as a heuristic) can be seen in Watson and Woodruff [32] and Mulvey and Vladimirou [28].

The main contributions of this work are: (i) developing a framework to consider explicitly climate change uncertainty in a spatially explicit stochastic multi-stage tactical forestry planning problem in Portugal by using a realistic wide range of climate scenarios, which are transformed into growth scenarios and (ii) demonstrating the usefulness of PH to solve this complex problem (due to number of constraints and variables). The rest of the paper is organized as follows. First, a problem description (Section 2.1) is provided; then, the scenario generation is presented in Section 2.2. The model is defined in Section 2.3. The solution approach is presented in Section 2.4. Results are placed in Section 3, and conclusions are in Section 4.

2. Materials and Methods

2.1. Problem Description

In this study, a tactical planning model for a eucalyptus forest located in central Portugal was considered. The forest managed by the industry (1 decision maker) encompasses 1000 stands and extends over 11,873 ha (ranging from 5.1 to 29.5 ha). In this region, the mean annual rainfall is 826 mm. Soils are mostly sandy and are characterized by low fertility, with low organic carbon content (0.23% to 0.28%) and an average of 395 mm (range between 242 and 737 mm) of water holding capacity.

The planner had to decide which stands to cut in each time period (i.e., develop a harvest scheduling management plan) while considering uncertainty due to climate change which affected tree growth. The importance of considering climate change when making decisions in this area was determined in interviews with the forest industry in the frame of the consultation process described by [33]. The industry showed their interest in considering explicitly the uncertainty due to climate change in their plans. This was to avoid the constant need of adapting their management plans due to the deviation between the real timber volumes harvested and the projections considered when developing their plans. As an output of these interviews, the objective of the management plan was set (i.e., maximization of economic returns).

The planning horizon considered 15 one-year periods. The goal was to fully harvest the forest during the horizon time where each stand had to be harvested once. The plan included the volume flow constraints to ensure smooth variations in total timber harvested from one period to the next one, to ensure a steady supply of timber and long-range sustainability of the forest. These constraints defined that the difference in the amount of volume harvested in two consecutive periods could not be larger than 15%. The plan also included timber demand constraints and adjacency constraints [34] to avoid harvesting contiguous areas exceeding a maximum size. The presence of these constraints made the problem considerably more difficult to solve. The adjacency constraints are related to environmental considerations (e.g., to avoid erosion, avoid reduction of visual quality) and state that if a stand k is harvested in period t, stands adjacent to stand k can only be harvested in the same period t if the harvest patch size including stand k is under a certain patch size. Green up period between harvests was not considered since, according with the industry, in eucalypt forests one year is enough to have a green area covered by new trees. In the present study, the maximum patch size was defined as the maximum allowable area of contiguous stands equal to 30 ha, as suggested by the pulp and paper company managing the area (see Section 2.3 for more details about how these constraints were constructed).

In Portugal, most eucalyptus stands are even-aged stands intensively managed under a coppice silvicultural system. In a coppice system, after the clear-cut, the regeneration is obtained from the shoots arising from the stump of felled trees. In general, a typical eucalyptus rotation may include up to two or three coppice cuts, each coppice cut being followed by a stool thinning that may leave an average number of shoots ranging from one to two. After the coppice cuts, the root system may be removed, and the stand is planted again. In order to calculate discounted net revenues, the average eucalypt pulpwood stumpage prices from the statistical data provided by the Portuguese Forest service were used. In this study, an average stumpage price of 36 € per ha and a discount rate of 3% were used.

Eucalyptus represents about 23% of the forest cover in Portugal (see the National Forest Inventory http://www2.icnf.pt/portal/florestas/ifn/ifn6). Raw eucalyptus material is used by the pulp and paper industry, a leading Portuguese export-driven industry. This explains why forestry growth and yield research in Portugal has targeted mostly eucalyptus forest ecosystems [35,36].

2.2. Scenario Generation

This paper focused on uncertainty due to climate change. Since climate affects forest productivity the paper deals with uncertainty in future growth and yield. The emphasis was on assessing the volume of timber [m³ ha⁻¹] harvested on each management unit. To address this problem several climate change scenarios were used, which were then transformed into growth and yield scenarios using a growth and yield model that take into account climatic data as predictor variables (i.e., a process-based model). These kind of models simulate detailed physiological processes that describe system behavior, while the empirical approaches rely on correlative relationships, which do not explicitly describe system behaviors and interactions [37].

Empirical models mainly describe the statistical relationships among data (i.e., tree measurements over time) with limited regard to an object’s internal structure, rules, or behavior [38]. In the empirical models, the projected results may be incorrect if the models are applied beyond the range of conditions represented in the calibration data [39]. That is why empirical models have not been considered suitable to predict growth under varying climate. On the contrary, process-based models predict growth and yield based on calculations of several physiological processes (e.g., canopy photosynthesis, respiratory losses, allocation of carbohydrates from photosynthesis). Recently, different options to improve the capability of empirical models to project growth under climate changing conditions have been presented (e.g., [40,41,42]. In our case, an innovative Decision Support tool DSS—SADfLOR v ecc 1.0 developed by [19] was used. This DSS uses a projection module that incorporates a process-based model [36] which is sensitive to environmental changes and therefore can be used to predict forest growth under climate change.

The climate dataset used was based on the EMSEMBLES project (http://ensembles-eu.metoffice.com/) that provided an ensemble prediction system for climate change. The climate change scenarios of the ensembles project were considered the most appropriate for Portuguese conditions [43]. The scenarios were obtained by assessing weather conditions and were characterized by maximum, average, minimum temperatures, precipitation, number of rainy days, and solar radiation. Those parameters were the input for the growth simulator needed to project the forest growth.

Based on the scenario tree presented in [22] that considered 32 scenarios and in order to cover and represent a wide range of climate variability for the study case area, other 9 different scenario trees were generated. The scenario trees differed in the data structure (i.e., number of nodes or children per period) while maintaining the probability distribution of the original tree [28]. Children are generated up to period seven. The different tree structures were generated to avoid an unfair comparison of the solving method (i.e., EF vs PH) since some particular scenario tree structures might beneficiate one solving method against the other. Each scenario tree included all stochastic parameters. A scenario was a particular realization of uncertainty through the whole time horizon [1], starting in root-node and leading up to a leaf. A stage represented a time period, where a stochastic parameter took a given value. For those scenarios that have been realized identically up to a certain stage, non-anticipativity constraints were built (see Section 2.3 for more details). Climate scenarios were transformed in growth scenarios providing values of timber volume available in each period through the horizon for each stand.

2.3. The Model

The stochastic optimization model presented each decision variable indexed by scenarios s ∈ S. The deterministic model used only the expected scenario (one scenario). The decision variables considered when each stand was going to be harvested. The objective function was to maximize the net present value, which corresponded to the revenues from timber sales. The constraints included the non-anticipativity constraints and the adjacency constraints.

Here, a brief description on how the non-anticipativity constraints were constructed is presented. The scenarios introduced in Section 2.2 may be presented as a tree diagram (Figure 1). Each scenario is the path from the first node (root node) to the end and represents the growth of the stand over the planning horizon. Each stage in the tree is a time period (in this study, one year) in which the stochastic parameters take a given value and correspond to a certain growth. A scenario group in a certain stage is formed by the set of scenarios where the growth up to that point has been equal.

The example in Figure 1 has four periods and four scenarios. the variables $x_{i,s}^{t=1}$, $x_{i,s}^{t=2}$, $x_{i,s}^{t=3}$, and $x_{i,s}^{t=4}$ for periods one, two, three, and four represent the harvest decisions. Thus, $x_{i,s}^t$ is the decision variable x in growth scenario s for period t and stand i. In this Figure, all scenarios share identical decisions in period one, whereas up to period two, two different groups are shown, on the one hand, scenarios one and two, and on the other hand, scenarios three and four. Thus, considering a single stand i, the non-anticipativity constraints are as follows: first stage, $x_{i,s=1}^{t=1}=\dots =x_{i,s=4}^{t=1}$; second stage, $x_{i,s=1}^{t=2}=x_{i,s=2}^{t=2}$ and $x_{i,s=3}^{t=2}=x_{i,s=4}^{t=2}$; third stage, $x_{i,s=3}^{t=3}$ = $x_{i,s=4}^{t=3}$; in fourth stage, no scenarios share information, and for this reason, there is no need to add any constraint.

The previous formulation known as the explicit/traditional formulation uses a constraint associated to each pair of variables that belongs to the same information set and thus, the number of rows added to the formulation is in the order of the cardinality of the information sets (which increases substantially the complexity of the model). In this paper, the compact formulation its used. This formulation allows reducing the size of the model—the number of variables and constraints —and thus improves its performance. In this formulation, each decision variable is associated to the nodes of the scenario tree (following its structure), and then, non-anticipativity constraints are implicitly satisfied (see [1,7,44] for more details). Thanks to this implementation, redundant constraints are eliminated obtaining a more efficient formulation when solving the problem using traditional solvers like CPLEX. In the model section below (Equation (8)), the explicit formulation for simplicity was included.

The main approaches to model adjacency constraints in forest management were presented by [45]. The two main methodologies to incorporate those constraints into a forest management model are the Unit Restriction Model (URM) and the Area Restriction Model (ARM). The first approach, the URM, disallows harvesting two management units that are adjacent (regardless of size). The second approach, the ARM, allows harvesting two adjacent management units only if the total size of the harvested area does not exceed the maximum opening size that is assumed. The ARM approach is difficult to formulate and solve. Nevertheless, research has been conducted to increase the efficiency of the solution of the more general case of ARM. One option called the Path Algorithm, consists of clustering contiguous stands, which if cut altogether would exceed the maximum opening size, and then adding a constraint that does not allow to harvest at those clusters [46]. Another option presented consists in enumerating all possible and feasible harvesting clusters and adding constraints that prevent harvesting adjacent pairs of clusters ([46,47]). Readers are referred to Goycoolea et al. [34] for a deeper discussion of these formulations.

The maximum patch size restriction was formulated as Path inequalities as described in Goycoolea et al. [34] which consists of listing all the minimally infeasible clusters. For each stand i and each time period t, a binary variable $x_{i,s}^t$ was defined. This variable took value one if stand s was to be harvested in period t and value zero otherwise. If a group of stands (e.g., stands one, two, and three) were all harvested in time period t = 1, and their total contiguous area was higher or equal to a certain limit (e.g., 30 ha), this cluster would be infeasible, and a constraint prohibiting to harvest all these stands simultaneously in period t = 1 would be generated.

$$x_{1,s}^{t=1}+x_{2,s}^{t=1}+x_{3,s}^{t=1}\le 2$$

Instead of defining all possible infeasible clusters (Λ+), which would be exponentially large, the Path algorithm ensured that any minimum infeasible cluster does not contain any other infeasible cluster. Then for all these clusters, an inequality was built forcing to harvest at least one stand less than the total number of stands included in the cluster.

Sets:

T: planning horizon, $t\in$ T;

I: management units, $i\in I$;

S: scenarios, $s\in S$;

Λ⁺: minimally infeasible clusters, $\mathrm{C}\in {{\mathsf{\Lambda }}^{+}}^{+}$

Nt: information set in period t.

Parameters:

$r_{i,s}^t$: discounted net revenue per ha if stand i is harvested in period t € per ha.;

$a_{i,s}^t$: volume harvested m³ in stand i in period t per ha., under scenario S;

A_i: area of stand i (ha.);

α: maximum allowed volume fluctuation (%) between two consecutive periods;

$\mathrm{Pr}\left(s\right)$: probability of scenario s;

Variables:

$x_{i,s}^t\in \left\{0,1\right\}$: equals 1 if stand i is harvested in period t under scenario S, 0 otherwise;

$V_{s,t}\ge 0$: total volume harvested in period t under scenario S (m³);

Objective function:

$$\mathbf{max}{{\displaystyle \sum }}_{s\in S}\mathbf{Pr}\left(\mathit{s}\right){{\displaystyle \sum }}_{\mathit{t}\in \mathit{T}}{{\displaystyle \sum }}_{\mathit{i}\in \mathit{I}}{\mathit{r}}_{\mathit{i},\mathit{s}}^{\mathit{t}}{\mathit{A}}_{\mathit{i}}{\mathit{x}}_{\mathit{i},\mathit{s}}^{\mathit{t}}$$

(1)

Constraints:

• Stand harvest 0–1

  $${{\displaystyle \sum }}_{tϵT}{x}_{i,s}^t=1\forall i\in I,\forall s\in S;$$

  (2)
• Volume harvested in period t

  $${{\displaystyle \sum }}_{i\in I}{a}_{i,s}^t{A}_i{x}_{i,s}^t=V_{s,t},\forall t\in T,\forall s\in S;$$

  (3)
• Even-flow of harvest

  $$V_{s,t}-\left(1+\alpha \right)V_{s,t-1}\le 0,\forall tϵ\left\{2,\dots T\right\},\forall s\in S$$

  (4)

  $$\left(1-\alpha \right)V_{s,t-1}-V_{s,t}-\le 0,\forall tϵ\left\{2,\dots T\right\},\forall s\in S;$$

  (5)
• Demand constraints

  $$D_t\le V_{s,t},\forall t\in T,\forall s\in S;$$

  (6)
• Adjacency constraints

  $${{\displaystyle \sum }}_{iϵC}{x}_{i,s}^t\le |C|-1,\forall C\in {\Lambda }^{+},\forall t\in T,\forall s\in S;$$

  (7)
• Non-anticipativity rule

  $$x_{i,s}^t=x_{i,s\prime }^t\forall s,s^{\prime }\in N_t,\forall t\in T,\forall i\in I,s\ne s\prime ;$$

  (8)
• Binary requirements

  $$x_{i,s}^t\in \left\{0,1\right\},\forall i\in I,\forall t\in T,\forall s\in S.$$

  (9)

Equation (2) states that each stand must be harvested only once. Equation (3) calculates the total volume harvested in period t and scenario s. Equations (4) and (5) are volume flow restrictions and do not allow harvest in period t more than 115% or less than 85% of volume harvested in period t−1. Constraints (6) ensure that minimum demand of volume is satisfied in each period. Constraints (7) are the maximum patch size restrictions, expressed by path formulation, fully described in [34], and preserve harvest of contiguous stands to not exceed 30 ha. This equation indicates that all the stands belonging to a certain unfeasible cluster cannot be harvested at the same time. Equation (8) are the non-anticipativity constraints. Expression (9) states the nature of the decision variables.

2.4. Progressive Hedging Approach

2.4.1. Methodology

Due to the high number of stands and scenarios, the stochastic optimization problem stated in Section 2.3 was not easy to solve. Adding scenarios to the extended model increased exponentially the number of non-anticipativity constraints, which made the problem even more difficult to solve using classic optimization techniques. To solve the extensive formulation for the original forest planning problem, CPLEX was tried, and the results were not satisfying in terms of solving time and gap (see Section 3). For some instances, CPLEX did not find any feasible solution during the time limit set (18,000 s for the problem without adjacency constraints and 36,000 s when using adjacency constraints). In some other cases, final gaps above 20% were achieved. In this context, PH was tested as it has proved to be a good scenario decomposition algorithm to employ in stochastic optimization problems [7]. By relaxing the non-anticipativity constraints, the scenarios could be solved independently. Since solutions for each scenario may not satisfy the non-anticipativity rule, a penalty term was introduced into the objective function, which attempted to minimize the infeasibilities for each variable set. Different heuristic rules were applied to gain a fast convergence to a close to optimal solution.

The PH algorithm that solved each problem separately was as follows:

$$\max {{\displaystyle \sum }}_{tϵT}{{\displaystyle \sum }}_{iϵI}{r}_{i,s}^t{A}_i{x}_{i,s}^t\mathrm{subject}\mathrm{to}:{\overrightarrow{x}}_s\in Q_s$$

(10)

The intuition behind the decomposition is the possibility of solving the original stochastic problem from the independent solution of each scenario and then calculate the expected value of the obtained solution.

The expression (10) states the problem constraints for each scenario s where $Q_s$ represents the feasible region of the problem for a particular scenario s. In each iteration, the non-anticipativity constraints were included in the penalty term of the objective function. However, convergence to the sub-problem (i.e., problem formulated for each scenario s) optimal solution under this formulation was not guaranteed. To avoid this situation, which could lead to unbounded sub-problems, a proximal quadratic term which forces variables to converge toward the right direction was added, penalizing the difference between the current solutions with the average solution from the previous iteration. (Algorithm 1).

Once each scenario was solved independently, the average solution in each node was computed. To satisfy the non-anticipativity rule, the average value of a variable should be equal to each solution found. If this was satisfied, then the problem was solved and the algorithm finished. Otherwise the next iteration of PH was penalized accordingly to the deviation of the solution and average value (${\overrightarrow{w}}_s$ penalty terms). The new problem, with an updated penalization for each scenario, was solved. The algorithm could iterate until it found a satisfying solution. Let N denote the set of nodes of scenario tree,$N_t$ the information sets in year t, ρ the penalty factor, $\mathrm{Pr}\left(\mathrm{s}\right)$ the probability of scenario s, and k denotes the iteration number. A pseudo code for the PH algorithm is given below.

In order to analyze the performance of the PH approach compared to the use of the extended formulation, several experimental instances of the original problem were solved (in terms of forest sizes and timber demand), see Section 2.5 for more details.

Algorithm 1. Progressive Hedging

begin

STEP 0

$k:=0$

$w^{0,s}(\forall s\in S);$

$\overline{x}=0$

STEP 1

Solve each scenario by max:

$ma{x}_{{\overrightarrow{x}}_s\in Q_s}{\displaystyle {\displaystyle \sum }_{tϵT}}{\displaystyle {\displaystyle \sum }_{hϵH}}{r}_{i,s}^t{A}_i{x}_{i,s}^t;$

STEP 2

Compute the average solution in each node:

$\overline{x}={\displaystyle {\displaystyle \sum }_{{}_{s\in S}}}\mathrm{Pr}\left(\text{s}\right)\overrightarrow{x_s}$

STEP 3

If the solutions are equal according to the criterion:

$||\overrightarrow{x}-\overline{x}||<\epsilon ;$

then terminate.

STEP 4

Update $\rho$

Update the penalty factor:

$w_{k,s}=\rho \left(\overrightarrow{x}-\overline{x}\right)+w_{k-1,s}\left(\forall s\in S\right);$

STEP 5

Solve each scenario with the penalty term:

$ma{x}_{{\overrightarrow{x}}_s\in Q_s}{\displaystyle {\displaystyle \sum }_{tϵT}}{\displaystyle {\displaystyle \sum }_{hϵH}}{r}_{i,s}^t{A}_i{x}_{i,s}^t+w_{k,s}{\overrightarrow{x}}_s+\frac{\rho }{2}||\overrightarrow{x}-\overline{x}{||}^2$

STEP 6

Update iteration number

k = k+1

Return to STEP 2

End

2.4.2. PH Heuristic

The PH algorithm converges to the global optimal solution for convex problems [14]. No such convergence proof exists when integer variables are present. Still, as a heuristic algorithm PH performs well for MIP problems. In this study, a series of heuristic techniques that have been successfully implemented in previous works [7,28] were tested.

Fixing Variables and Compact Model

This technique consists of fixing those variables that satisfy the non-anticipativity constraints during the iterations of the algorithm. The key is to reduce the number of free variables of the original model and solve a smaller and easier compact MIP model. In order to obtain better results, a relatively high percentage of 0–1 variables need to be fixed; then, the original Extended Formulation becomes significantly smaller and can be solved using a commercial MIP solver [32,48].

There is a trade-off between solving time and quality of the solution, depending on the number of convergent iterations before fixing variables. A convergent iteration consists of an iteration where all variables have the same value inside their information nodes, satisfying the non-anticipativity constraints. Less convergent iterations can lead to lower quality solutions since the variables can be fixed too early without reaching a good global convergence (e.g., convergence only in some nodes of the tree), but waiting for more iterations can lead to larger solving times, obtaining similar results to the ones obtained by the extended formulation. Then, the number of iterations can be determined through experimentation, analyzing this trade-off between solving time and quality of the solution.

Penalty Factor

As Watson and Woodruff [32] suggest, selecting a penalty factor ρ value depends on a given problem. Different values and methods to select ρ were tested. The PH algorithm has a particular sensitivity to the ρ parameter, which is associated with the quadratic term (Step 5 in Algorithm 1) used to force the non-anticipativity in the variables of the scenario sub-problems. Because of this condition, it is very important to set its value or define how to calculate it, since the algorithm’s performance is significantly affected by this parameter.

Finding a good value for this parameter allows a more stable and better convergence in terms of speed and quality (smoothness) of the solution. In the present study, different strategies to treat this parameter, as described by Badilla et al. [7] were tested, in particular, FX (Fixed, same penalty term for all iterations), CP (Cost proportional, updating ρ in every iteration based on decision costs), and SEP (Selecting per element, where the penalty term depends on the scenario).

Hot Starts

In every iteration, each scenario sub-problem was solved until some gap or time was reached. The solutions obtained in iteration k of each sub-problem were used as a starting point for the optimization of each scenario in iteration k + 1 (i.e., hot start) and so on. Note that a hot start for each scenario would be always feasible since all the modifications were made in the objective function, not in the sub-problems structure. The use of this technique sped up considerably the solving time (up to 50%) [28].

Dynamic Gap

As the solution for the original model obtained in the first iteration was not satisfactory—no convergence—there was no need to solve the sub-problems to optimality (or very low gap). The algorithm’s performance was improved by setting a high initial gap value and then decreasing this value in each iteration. Using this technique, better convergence was obtained and therefore, faster and better solutions.

2.5. Instances of the Original Problem

In order to analyze the performance of the PH approach proposed when dealing with different problem sizes, several experimental instances of the original problem (i.e., forests with different sizes) were developed. These instances were created by selecting 100, 200, 400, and 800 connected stands from the original forest (1000 stands). These stands were connected since we wanted to apply adjacency constraints. For each instance, a stand was randomly selected as an initial point (seed), and then stands in the “neighborhood” (adjacent stands) were selected until the total number of stands of each instance was achieved. Thus, the final area was forced to be contiguous (i.e., no isolated stands are allowed). Fifteen different forests (15 different random seeds) for each size were generated which results in a total of 61 forests (60 new sub-forests and the original one). All the stands selected were characterized with their original forest data. The objective of the planning problem for each forest size was the same as presented in previous sections (maximization of net returns subject to even-flow of harvests using different timber demand levels). Different demands, defined as 0%, 25%, 50%, 75% 90%, 95%, 100% of the potential forest production, were tested in each forest. In order to calculate the potential forest production, an auxiliary MIP model was created for each climate scenario, where the objective was the maximization of the total volume harvested (per instance) for the complete planning horizon. The combination of forest sizes and demand bounds resulted in 427 instances. In addition, each problem was run with and without adjacency constraints for all the different instances and demand bounds, which finally resulted in 854 instances. These adjacency constraints were added in the original model using the minimum infeasible cluster formulation [34] as explained in Section 2.3.

2.6. Hardware and Software

The basic experiments were carried out in a computer with the following features:

Core: Intel I7 4700 RAM Memory: 16 GB

OS: WIN 8.0 64 Bits

CPLEX: Version 12.6 64 Bits

PYTHON: Version 2.7 64 Bits

Pyomo Sandia Optimization Package Version 4.1.10527

All the PH algorithm and extensive form scripts were written in Python using Pyomo modeling language, as part of the Pyomo Sandia Optimization Package. These scripts performed all the operations needed on each iteration, like updating the penalty term, checking and computing the convergence metric, etc., and for each scenario sub-problem, calling and running CPLEX in order to solve them.

3. Results

3.1. Configurations

In order to achieve the best performance of the PH algorithm, it was necessary to select the parameters and techniques presented in Section 2.4.2 (e.g., strategies to select the ρ parameter or penalty factor). After a series of experiments using all the instances, the CP strategy was chosen due to the high performance and convergence speed achieved by the algorithm for the model under study. In the CP strategy, the ρ parameter was calculated according to costs/benefits associated with each of the decisions from the model, different penalties were produced proportionally with an increase in the penalty term when there was no improvement in a particular solution, which drives it to a feasible solution.

Every solution obtained in each iteration was used as a starting point for the next one (hot starts), saving more than 50% of the solving time per iteration. The solution obtained after running the PH algorithm was used as a starting point for the extended form. In order to accelerate PH algorithm, the “dynamic gap strategy” was used, which consisted of decreasing the value of gap while iterating. Based on the performed experiments, the strategy beginning with a 20% gap in the first iteration and decreasing in a linear way until a one per cent in the last iteration obtains the best performance, in terms of quality (value) and solving time.

As a heuristic, binary variables that satisfy non-anticipativity constraints for a certain node for k consecutive iterations were fixed, obtaining a reduced model (i.e., an extensive formulation including the fixed variables) by the end of the algorithm. Note that this could lead to a non-optimality. Based on the experiments, the best performance was achieved when binary variables having value one for five consecutive iterations were fixed. Variables with zero values were not fixed since they were the complement of the fixed ones. Using this configuration, over 90% of the binary variables were fixed for all the instances, obtaining an easy-to-solve reduced model. After any of the stopping criteria was reached, the reduced model was solved.

This technique could not be applied in the most difficult instances (model with adjacency constraints and demand bounds of 95% and 100%) due to the bad convergence of the method, obtaining poor solutions (low values for the objective function). Based on the experimental results, stopping criteria were established by fixing (at least) 80% of the binary variables or reaching a convergence metric value ɛ < 0.001 within a maximum solving time.

Due to the size and complexity of the extended model, commercial solvers could not handle the problem. In some cases, the computer ran out of memory (e.g., cases of 95% and 100% of demand with and without adjacency), and in others, it could not reach a low gap solution (less than seven per cent). In order to compare the use of PH and the EF with and without spatial constraints, time and/or gap tolerances were changed. For a model without adjacency constraints, the stop criteria were a time limit of 18,000 s and a gap of one per cent. For instances with adjacency constraints, the new time limit was set to 36,000 s and gap one per cent since spatial constraints add an important complexity to the model, thus, solvers require more time to obtain satisfactory solutions. For the most difficult instances (larger forests with higher demand bounds) a new gap threshold of two per cent was set up.

3.2. Performance Comparison

Based on the results obtained, the problem became harder to solve when the forest size and the minimum demand constraints increased, especially when the minimum demand was close to the maximum forest’s production capacity. However, the global effect was more important in the case of EF where performance was highly correlated to the size of the instance, something that was not as strong when using PH. Solving times increased exponentially for the extended model while for the PH approach were quite similar for all instances except for the instance of 1000 stands and demand levels higher than 90%. In the latter case, no solution below two per cent of gap could be reached within the maximum solving time (

Table 1. Comparison of performance of PH algorithm (PH) versus the extended formulation with CPLEX (EF) for all the instances including adjacency restriction and demand levels D. PH gap is computed using the best bound given by CPLEX when solving EF. The results shown in this table are the averages found for each of the 15 forests generated for each instance. NO SOL means that no solution was found in more than 50% of the 15 forest repetitions.

	EF		PH		Statistics
Instance	Gap [%]	Time [s]	Gap [%]	Time [s]	Saved Time [%]
100SD0%	NO SOL	NO SOL	1.2	354.5	-
100SD25%	23.2	36,000	7.5	2458.1	91.7
100SD50%	5.2	36,000	4.8	920.2	98.3
100SD75%	4.3	36,000	3.9	469.3	98.6
100SD90%	3.6	36,000	3.5	765.2	98.7
100SD95%	NO SOL	NO SOL	0.9	226.00	-
100SD100%	NO SOL	NO SOL	1.9	5974	-
200SD0%	NO SOL	NO SOL	3.8	608.4 *	-
200SD25%	49.4	36,000	5.4	3963.5	88.9
200SD50%	27.9	36,000	6.2	827.1	97.7
200SD75%	6.4	36,000	6.1	261.7	98.5
200SD90%	2.9	36,000	2.7	224.6	98.7
200SD95%	NO SOL	NO SOL	1.6	5410.6	-
200SD100%	NO SOL	NO SOL	1.3	5853.8	-
400SD0%	NO SOL	NO SOL	3.4	3006.7	-
400SD25%	6.1	36,000	5.8	3996.8	69.6
400SD50%	5.2	36,000	4.9	1082.2	95
400SD75%	7.2	36,000	6.8	514.4	97.1
400SD90%	5.0	36,000	4.7	429.0	97.5
400SD95%	NO SOL	NO SOL	2.6	6887.6 *	-
400SD100%	NO SOL	NO SOL	1.7	6068.4	-
800SD0%	NO SOL	NO SOL	2.3	3109.4	-
800SD25%	1.0	11,092	1.1	3268.6	66.6
800SD50%	1.5	9029	1.7	2275.1	76.3
800SD75%	7.8	36,000	7.5	1406.8	92.2
800SD90%	8.8	36,000	8.1	1390.9	93.8
800SD95%	NO SOL	NO SOL	4.7	9519.8 *	-
800SD100%	NO SOL	NO SOL	2.7	11,603.2 *	-
1000SD0%	NO SOL	NO SOL	1.8	4844	-
1000SD25%	0.9	6061	1.0	1574	91.4
1000SD50%	1.0	8279	1.3	5283	91.6
1000SD75%	7.9	36,000	7.5	2404	92.9
1000SD90%	8.8	36,000	8.6	1787	93.8
1000SD95%	NO SOL	NO SOL	13.7	36,000	-
1000SD100%	NO SOL	NO SOL	13.1	36,000	-

* means that memory issues where experienced during the optimization, stopping the algorithm.

).

Results showed that PH approach was more suitable than the EF model, both in terms of solving times and gap of the solution achieved. PH heuristic performed better than the EF model in most of the cases with savings on solving time between 70% and 100% (Table 1). Only in few unconstrained (in terms of minimum demand) problems, the EF performed similarly to PH in terms of final gap but performed much worse in solving times. In addition, when the demand was not constrained (i.e., 0% demand level) using the EF, CPLEX did not find a solution due to the high number of possible combinations (i.e., the solution space was too big and CPLEX did not find a feasible solution in a reasonable time). For high demand levels (95% to 100%) it could not find any feasible solution—in any instance—after running for 18,000 s. Nevertheless, the PH algorithm using fixing heuristic was able to solve all the instances in a reasonable time, with a very low gap.

In general, the instances using 200 stands were the most difficult to solve using the EF, reaching very high average gaps (up to 49.5%). PH solving times were far better in comparison with all the EF solutions, with an average (for all instances) smaller than 2000 s, reaching very similar or better gaps than the EF. For PH, most of the final gaps were between 2% and 10%.

When using PH with fixing strategy, forests with demands of 0% and 25% were the most difficult instances, due to the higher number of feasible combinations. Good convergence was obtained in instances from 50% to 90% of demand, where the algorithm fixed more than 90% of the binary variables leading to an easy to solve reduced model.

PH convergence was not good for instances of 95% and 100% of demand where no more than 40% of the variables could be fixed, and their fixed values were not part of a reasonable solution, leading to bad quality results (i.e., infeasible or with a very high final gap). Then, for these instances, no fixing strategy was followed, using PH solution as a hot start for the EF (see Section 2.4.2). This approach allowed solving all the high-demand instances, reaching an average of 2% of final gap (in comparison with the LP bound), for the forest instances with 100, 200, 400, and 800 stands. Average solving times were 6000 s (for instances of 100–400 stands) and 12000 s (for 800 stands). For the instance with 1000 stands the gap was 12.4% reaching the maximum solving time (18,000 s). The 1000 stands instances became very difficult to solve after 95% of demand threshold, due to its size and these strong lower bound constraints. The final gap showed in Table 1 is the gap computed between the final PH solution and the EF best bound or LP solution for the cases where EF could not be solved.

Comparing solving times of both approaches, PH saved at least 70% in all instances, with an average around 94%. In addition, PH final gaps were better (lower) or equal than the ones obtained using EF, thus, results demonstrated that PH was better than EF. Especially, the PH algorithm outperformed EF, in all forest sizes and high demand levels (i.e., 95% and 100% of demand) where EF could not find any solution. Then, the implementation of PH successfully showed a great performance in obtaining good solutions with better solving times than EF for all the instances.

4. Discussion and Conclusions

In this paper, a difficult multi-stage stochastic forest planning optimization problem is successfully solved under climate uncertainty for a real forest in Portugal including 1000 harvest units. The problem is solved with and without spatial constraints known as adjacency constraints for different forest sizes: 100, 200, 400, 800 (sub-forests), and 1000 stands (original forest). Uncertainty is represented by a scenario tree structure, including 32 climate scenarios per instance. This uncertainty affects the forest growth, leading to scenarios with different timber production. The planning horizon of the problem covers 15 years, where each stand must be harvested by the end of this horizon, satisfying even flow constraints that aim to have smooth flow deviations within consecutive years. This is a similar problem than presented by [22]. The main difference is that in the present study spatial integrity of the planning units (i.e., stands) was considered, which made the problem much more difficult to solve. The value of explicitly introducing climate uncertainty through the use of stochastic programming in harvest scheduling models was already presented by [22]. They showed the superiority of the plans using stochastic programming versus the use of deterministic models, superiority in terms of objective values, and feasibility of the solutions.

The methodology shown in the present study, using mixed integer stochastic model, provides valuable support for forest managers in making more robust harvesting decisions. The superiority of the use of multiple stochastic scenarios compared to the use of one scenario with average values has been also shown by [1,7,22]. Another clear advantage is the possibility to adapt forest management plans to climate dynamics. Adaptive management is a structured process of learning about the problem and adapting the management decisions based on that learning [49]. In a stochastic process, learning could be interpreted as the resolution of some uncertainty. In this context, in our study, the planner can adapt easily the decisions during the planning horizon since he can observe which one of the scenarios is closer to the real climate observed in the period and adapt the decisions accordingly.

The results of this study highlight the conclusion that stochastic programming models can be implemented into current forest management tools to be used in forest planning by forest owners and planners. Nevertheless, as stated by Eyvindson and Kangas [50], the value of using stochastic programming needs to be presented to both the forest owners and the forest planners. This may also be expanded to the need of moving from optimizing forest plans in a stochastic setting and accounting for the forest owners’ preferences towards risk [23].

Regarding the solution algorithm, the usefulness of the PH algorithm is demonstrated in solving stochastic multi-stage problems of large size (in terms of scenarios, variables, and constraints) using several climate scenarios in a spatially explicit planning harvest scheduling model. Since the problems include non-anticipativity constraints, it was not easy to solve it under the extended formulation approach using commercial solvers due to the memory limits, computational time, or both. The advantage of PH lies in decomposing the problem by scenarios and solving them separately, using non-anticipativity expressions as a penalty term in the objective function instead of explicit constraints.

In general terms, when the problem becomes harder to solve, the PH algorithm implementation performs significantly better than the extended formulation even without using distributive or parallel programming. A clear pattern can be seen in all instance sizes, both small and large (in terms of the total number of stands in the forest) with low and high demand, where EF was outperformed by PH consistently. This difference in performance is increased in the case of larger and more constrained instances (full forest with higher demand lower bounds). The average saving times are above 90%, and average gaps are up to a 4% better than the EF solutions (when a solution is found) showing the advantage of using the decomposition and heuristic approach.

As expected, the instances become more difficult and complex when the model becomes harder to solve (i.e., including a larger number of stands and higher demand’s lower bounds). However, the effects and impacts of both, solving time and final gaps were significantly smaller for the PH approach in comparison to the EF implementation.

The heuristics and enhancements implemented when using PH algorithm as a heuristic approach were crucial in order to obtain the best performance of the algorithm. In particular, techniques such as hot starts, how to select/update the quadratic penalty term, dynamic gap, and fixing variables strategies are necessary when dealing with complex non-convex problems, otherwise, the performance of the algorithm can be very poor, leading to very bad quality results.

Overall, PH showed a great performance, obtaining good solutions with better solving times and gap than EF for all the instances. Thus, it is demonstrated that PH approach is the most appropriate for the problems where a large number of scenarios and constraints are used.

This research may be expanded in the future by considering also uncertainty in future timber markets, in terms of future prices and/or demand levels or even considering fire risk. In addition, this research may be expanded by developing and implementing a suitable parallel programming scheme—inherent to the algorithm PH—which will potentially allow for direct execution in regards to the problem addressed as well as similar problems, thus generating a considerable saving of time in solving each of the sub-problems during the iterations of the method, with consequent efficient use of resources.

Investigate metrics and/or methodology for establishing “sufficient” number of scenarios is to be included in the model, in order to be able to represent in the best way real world scenarios based on the current climatic reality and to compare the quality of the solution found for many instances, according to the scenario tree used, allowing to achieve a balance between quality of solution and solving times involved.