# A Parallelized Variable Fixing Process for Solving Multistage Stochastic Programs with Progressive Hedging

Authors: Martin B. Bagaram; Sándor F. Tóth; Weikko S. Jaross; Andrés Weintraub
Journal: Advances in Operations Research
 (2020)

---

## Abstract

Long time horizons, typical of forest management, make planning more difficult due to added exposure to climate uncertainty. Current methods for stochastic programming limit the incorporation of climate uncertainty in forest management planning. To account for climate uncertainty in forest harvest scheduling, we discretize the potential distribution of forest growth under different climate scenarios and solve the resulting stochastic mixed integer program. Increasing the number of scenarios allows for a better approximation of the entire probability space of future forest growth but at a computational expense. To address this shortcoming, we propose a new heuristic algorithm designed to work well with multistage stochastic harvest-scheduling problems. Starting from the root-node of the scenario tree that represents the discretized probability space, our progressive hedging algorithm sequentially fixes the values of decision variables associated with scenarios that share the same path up to a given node. Once all variables from a node are fixed, the problem can be decomposed into subproblems that can be solved independently. We tested the algorithm performance on six forests considering different numbers of scenarios. The results showed that our algorithm performed well when the number of scenarios was large.


---

## Body

## 1. Introduction

Uncertainty is common in all disciplines that involve decision making. In forestry, planners need to prescribe, several decades in advance, actions that should be taken in a forest in order to achieve a given management objective. Those actions include the segments of roads that should be built in a given period to allow hauling, the forest units or stands that should be treated to reduce the risk of fire, the stands that should be logged in each time period to produce timber and secure employment to the local community, etc. The most common objective is the maximization of the net present value subject to environmental, budgetary, and logistic restrictions [39]. The prescriptive models allowing to assign forest units to different actions in different time periods are known as harvest-scheduling models.To build harvest-scheduling models, forest planners have traditionally used expected growth and yield coefficients to predict future merchantable timber volumes. However, uncertainties in long-term temperature and precipitation coupled with increasing wildfire, windstorms, or landslides due to climate change may affect forest development. The traditional approach fails to account for uncertainty in forest growth and leaves forest planners unprepared in case of occurrence of any of these uncertainties. The uncertainty in forest planning can be addressed by a particular class of mathematical programming known as stochastic programming. One can distinguish between two-stage stochastic programming and multistage stochastic programming. In two-stage stochastic programming, a decision is made, and then the uncertainty is revealed and a recourse action that depends on the revealed uncertainty is taken. However, in the case of multistage stochastic programming the sequence of decision and uncertainty revealing itself occurs more than once giving therefore more flexibility to the decision maker to take recourse actions as uncertainty unfolds progressively.Multistage stochastic programs are mainly composed of scenarios and stages. Scenarios are the set of possible future outcomes of the uncertain parameter, and stages represent the level at which decisions can be made and/or recourse actions can be taken. We call “period,” the time that separates two consecutive stages. Uncertainty unfolds progressively in each period. Since at the beginning of the planning, we cannot anticipate which scenarios would unfold, we require that the decision made at the first stage must be the same for all scenarios. This requirement is known as nonanticipativity constraints. Nonanticipativity constraints (NACs) impose that if two scenarios cannot be distinguished up to any given stage, then the decision made in the two scenarios up to that stage must be the same. Multistage stochastic programming has been commonly used to model uncertainty in forest harvest scheduling because of the long planning horizon that characterizes forest harvest scheduling and the presence of interdependent relationships between periods such as the maximum contiguous area harvested from one period to another should not exceed a given limit. For instance, [41] solved a multistage stochastic harvest-scheduling model with uncertainty in wood price, wood demand, and productivity while Alonso-Ayuso et al. [1] focused on the uncertainty in the price and risk aversion. Finally, Álvarez-Miranda et al. [2] assessed how key ecosystem services change in forest management if there is growth uncertainty stemming from climate change.Although climate change might be one of the biggest challenges of forest management, it has received little attention in forest harvest planning in part because stochastic programs, especially multistage stochastic programs, are considered one of the most challenging classes of optimization problems to solve [5, 23, 40, 45]. For instance, the number of scenarios in [2] was limited to 32 although climate scientists forecast at least four climate pathways [31] which may translate into hundreds of possible forest growths.Real-life applications of multistage stochastic problems lead to large models that are hard to solve directly [37]. The size of the models is directly related to the number of scenarios used to represent the uncertainty. Although uncertain parameters may be continuous, to make the problem suitable for stochastic programming, the discrete realization of the uncertainty is cast in a structure known as a scenario tree. In the scenario tree, each node represents a state decision and the branches of the tree represent the realization of the uncertainty. As more scenarios are included in the tree, the model will better approximate the entire probability space of the uncertain parameter [26]. However, increasing the number of scenarios makes the resulting stochastic programming models hard to solve. It is therefore necessary to develop some decomposition algorithms. The two most common decomposition methods are Benders decomposition [7] also known as stagewise decomposition or vertical decomposition and progressive hedging [36] also called scenario-wise decomposition or horizontal decomposition. The study [38] provided an overview of the algorithms for stochastic programming decomposition. Additional summary of different solution methods for stochastic programming is available in [16].Benders decomposition (BD) is a delayed constraint generation approach for solving mixed integer programs. It has been mainly applied to two-stage stochastic programming problems with some assumptions on the nature of the first-stage and second-stage variables. The authors in [27] give a summary of conditions for application of BD to two-stage stochastic programs. According to [15], BD is well suitable when there is only a small set of constraints that prevent the decomposition of the problem into blocks. This is not the case in forest harvest scheduling characterized by long planning horizons and constraints such as even flow of wood and ending inventory linking variables from one decision stage to another of multistage stochastic programs. For instance, Egging [15] had limited success when extended their BD algorithm to multistage stochastic programs. We acknowledge, however, that there are multistage applications where BD had been used [17, 43, 44].Progressive hedging (PH), on the other hand, was developed by [36] for convex two-stage stochastic programs and is proven to produce a global optimal solution for continuous problems. For nonconvex problems such as stochastic mixed integer programs (SMIP), the algorithm is not proven to converge. It relaxes nonanticipativity constraints and iteratively solves the stochastic program by independently solving its scenarios and penalizing the violation of nonanticipativity constraints. PH is appealing because it is a method based on scenario-wise decomposition and thus at each iteration the stochastic program solved is equivalent to a risk neutral problem that is the same as the deterministic problem which ignores uncertainty. Consequently, the algorithm can deal with a large number of scenarios.Despite its benefits, PH performs poorly for nonconvex multistage stochastic programs. Several researchers have proposed PH-heuristics for solving SMIP in different applications. The most promising PH-based heuristic explored is fixing variables that participate in defining nonanticipativity constraints as they meet consensus (see Section2.2 for definition of nonanticipativity variables). For instance, Veliz et al. [41] fixed nonanticipativity variables during PH iterations as they meet consensus. After a given percentage of variables, 80%, for instance, consent on the value they should take, the reduced problem, which is the SMIP with some NAC variables fixed (named here and after reduced extensive form, REF), is then solved directly. The inconvenience of this approach is that the reduced extensive form may be infeasible because too many variables have been fixed (for proof, see Appendix). The infeasibility occurs when there is uncertainty in the yield of the forest such as in the case of climate change. In this case, the algorithm wasted considerable time iterating. On the other hand, solving the reduced extensive form problem can still be difficult if the number of variables fixed is too low. The third limitation stems from the fact that the reduced form problem is not separable and thus solving it cannot be parallelized. Similarly, for a two-stage problem in resource management, the authors in [42] proposed a scheme for fixing variables. Their algorithm applicability was limited to a special class of resource management where constraints are one sided. They proposed “slamming” technique which forces nonanticipativity variables to converge, although they have not met consensus yet. They showed that this technique accelerated PH convergence. Recently, the authors in [30] proposed fixing variables as well. In their case, they had both binary and continuous decision variables and proposed therefore to fix the binary NAC variables, and the resulting linear stochastic program, since convex, can easily be solved using progressive hedging.Aside from fixing variables, researchers have investigated other strategies for improving PH application to nonconvex problems. Some of the strategies are only applicable to two-stage stochastic programs. For instance, Atakan and Sen [4] proposed a branch-and-bound algorithm for stochastic mixed integer programs and tested the algorithm on stochastic server location problem. Similarly, Barnett et al. [6] proposed a combination of branch-and-bound and PH using the former as a wrapper. In the same context, Gade et al. [19] proposed an algorithm for computing the lower bound of PH for two-stage SMIP. Other researchers invested into reducing the duality gap that arises from relaxation of the nonanticipativity constraints [9, 19]. In the same spirit, Boland et al. [10] looked into the minimum number of nonanticipativity constraints that need to be reinforced to ensure the quality of SMIP. The authors in [8, 33] looked at reducing the number of scenarios by clustering them into bundles that can be solved independently. All these algorithms remain dependent on some assumptions of the nature of the first-stage variables or the second-stage variables, or both. In addition, their applicability is limited to two-stage SMIP and specific to certain disciplines. However, the main challenge that remains is that PH does not converge for nonconvex problems which are common in forest planning. In forest planning, there are some problems that possess only binary decision variables (there could be some accounting continuous variables which do not participate in decision making). Some of those problems are found in spatial explicit forest management planning where the decision variables are the forest units that should receive a specific treatment. This class of problems with pure binary decision variables is the focus of this research.In this work, we aim at overcoming some limitations that are posed by using PH for multistage stochastic harvest-scheduling problems with uncertainty in the yield. The objective is to have a heuristic that efficiently solves multistage stochastic harvest scheduling with little loss of optimality. (i) We propose decomposing SMIP into scenarios and solving them by a special heuristic form of progressive hedging that is completely parallelizable. We capitalize on the idea of fixing variables as they meet consensus and extend it so that the reduced extensive form problem is parallelizable. Our form of PH heuristic exploits the structure of the scenario tree by fixing variables starting from the root-node. (ii) We assess the impact of fixing variables on the value of the optimal solution. (iii) We investigate as well some acceleration strategies that allow to accelerate the algorithm by extending the slamming technique proposed for two-stage stochastic programs by [42] to multistage stochastic programming. (iv) We use climate change to show one source of uncertainty in the yield although there can be other sources of yield uncertainty such as the errors in measurement of forest growth or uncertainty of future yield due to prediction models. The algorithm is suitable as well to price uncertainty. Our aim is not to model perfectly the uncertainty in forest growth due to climate change but to provide an algorithm that can be used for real-life application in stochastic forest harvest scheduling.The remainder of this paper is organized as follows. In Section2, we provide a background on stochastic programming problem formulation and formally introduce the progressive hedging algorithm. In Section 3, the special form of PH for solving forest harvest planning is presented. We provide in Section 4 the computational experiment and the interpretation of the results. Finally, Section 5 concludes the paper and presents the limitations of the algorithm and future work.


## 2. Background

To illustrate variable fixation variant of the progressive hedging algorithm, it is necessary to represent the scenario tree that is the abstraction of the realization of uncertainty in growth and yield and formerly introduce progressive hedging (PH) algorithm.


### 2.1. Scenario Representation

The stochastic program representation can be visualized as a tree, the so-called “scenario tree.” It can be represented as follows. LetT denote the set of periods in the planning horizon with T=T being the number of periods. In the tree, a node represents the realization of the uncertain parameter and variable at a given time period. Let n and N describe the node and the lexicographically numbered set of nodes 1,…,N in the tree, respectively. From each node n, for t∈T∖T, there is at least one branch leading to another node m with a conditional probability pm with m∈ℐn, that is, the set of nodes, immediate successors of the node n. ∑m∈ℐnpm=1. We denote by Nt the subset of nodes belonging to period t such that N=∪t∈TNt and Nt∩Nt+1=Φ for t∈T∖T. Let Ω represent the finite set of representative scenarios in the tree. A scenario ω∈Ω is a particular realization of the uncertain parameter represented as a path from the root-node to a leaf-node. Each scenario ω has an associated probability or weight denoted by wω. Note that ∑ω∈Ωwω=1. Similarly, let Nω represent the set of nodes forming the scenario ω. In other words, Nω is the set of nodes in a path from the root-node to a leaf-node. There is a number of scenarios that traverse each node. Let Ωn denote the set of scenarios that traverse the node n. We have Ω1=Ω, Ωn∩Ωn′=Φ∀n≠n′ and n,n′∈Nt. Finally, let S˜n denote the set of successor nodes to the node n, for n∈N. We have N1 is a singleton and S˜n=Φ for n∈NT (leaf-nodes). Furthermore, let αn denote the immediate ancestor node of node n for n∈N∖1, since n=1 is the root-node of scenario tree. To introduce nonanticipativity, we need to define η as the set of nodes having more than one leaf-node as successor S˜n∩NT>1⟹n∈ηA. For instance, in Figure 1, nodes 2 and 3 have as leaf-node successors the set of nodes {5, 6} and {7}, respectively. Therefore, node 2∈η but node 3 does not. Finally, Let Xω represent the matrix of variables for all t∈T for the scenario ω∈Ω. As the result, the stochastic program can be stated as follows:(1)maxX∑ω∈ΩwωfωX,(2)s.t.Xω∈Cω∀ω∈Ω,(3)Xω=Xω′,∀ω,ω′∈Ωn,ω≠ω′,n∈η.Figure 1
Scenario tree representation of stochastic programming.Equation (1) maximizes the expectation from all scenarios. Equation (2) says that all the solutions should be feasible with respect to the constraints of the scenario. Cω is the feasible set for scenario ω. Finally, equation (3) imposes the nonanticipativity constraints (NACs) which require that the solution up to a period t from two scenarios should be the same if the two scenarios are indistinguishable up to period t. We call this model the extensive form (EF).Notice that we could write constraint (3) in a different form. Let Xnω:=x1nω,x2nω,…,xsnω be the vector of variables at node n∈η under scenario ω, and xsnω is the state variable of the forest unit or “stand” s. Let Zn represent a vector of binary variables at node n. Imposing constraint (4) is equivalent to reinforcing NAC. Progressive hedging exploits the following formulation:(4)Xnω−Zn=0,∀ω∈Ωn,n∈η.


### 2.2. Progressive Hedging

Progressive hedging (PH) is a scenario-based decomposition algorithm proposed by [36] for stochastic programming models. The idea of progressive hedging is to relax the nonanticipativity constraints (NACs, equations (4)) in an augmented Lagrangian manner [18, 22, 35, 36] so that each scenario (subproblem) can be solved independently. This assumes that solving subproblems independently is much easier and faster. Nonanticipativity constraints require that values of variables that share the same ancestor nodes should be equal across scenarios up to that node. In other words, if two scenarios ω and ω′ are indistinguishable up to period t, then the solutions of the two scenarios should be the same up to that period. For instance, in Figure 2, each variable at time t=1 should be the same across all the five scenarios and values of variables in period t=2 should be the same for scenarios 1 and 2. However, these constraints are not reinforced when scenarios are solved independently. Through progressive penalization of violations of NAC, the algorithm is proven to ultimately converge to the optimal solution for convex stochastic programs. In a nutshell, progressive hedging follows the following steps:(1)
Solve each scenario without penalization(2)
Compute the averagez for each variable(3)
If solutions have sufficiently converged, then stop(4)
Update penalization termsλ=ρx−z+λ(5)
Solve scenarios with penalization terms(6)
Go to step 2Figure 2
Example of scenario representation in a compact form (a) and in a detached form (b). The red-dotted ellipses show the nonanticipativity constraints that need to be imposed.
(a)(b)Algorithm1 describes progressive hedging for multistage stochastic programming. The inputs of the algorithm are the penalty factor ρ, the maximum number of iterations kmax, and the termination criterion ϵ which indicates the level of consensus of nonanticipativity constraints that is acceptable. ε=0 means that the algorithm stops if all the nonanticipativity constraints are satisfied. In the algorithm, line 2 initializes a Lagrangian multiplier λ for each NAC constraint. A Lagrangian multiplier is associated with each equation (4). To compute the average, of each variable, we need to have the conditional probability associated to each branch of the scenario tree in the detached form (Figure 2). Equation (5) computes the conditional probability from node n to node m which is an immediate successor. That probability qnmω is proportional to the number of leaf-nodes associated with node m because the number of leaf-nodes informs on the number of scenarios passing by node m. For instance, from Figure 2, q121=q122=p2/2, q251=p5, and q133=p3. Lines 6 and 17 compute the average of each variable at node n∈η. Lines 7 and 18 update the Lagrangian multiplier associated with each nonanticipativity constraint. Line 14 computes the NAC convergence euclidean distance. This distance is zero if all the NACs are satisfied:(5)qnmω=pmΩm,∀m∈ℐn,m∈Nω.Algorithm 1: Progressive hedging multistage.
(1)
function ph ρ,kmax, ε(2)
Initializeλn0,ω=0,∀n∈η,∀ω∈Ωn(3)
forω∈Ωdo(4)
x0,ω=argmaxxfωx:x∈Cω(5)
end for(6)
zn0=∑ω∈Ωnqnmωxn0,ω,∀n∈η,∀m∈ℐn(7)
λn1,ω=λn0,ω+ρxn0,ω−zn0,∀n∈η,∀ω∈Ωn(8)
fork=1 to kmaxdo(9)
forω=1 to Ωdo(10)
xk,ω=argmaxxgωx=fωx−λk,ωx−zk−1−ρ/2x−zk−122:x∈Cω(11)
ϕωk=maxxgωx(12)
end for(13)
ϕk=∑ω∈Ωwωϕωk(14)
if∑ω∈Ωwωxk,ω−zk−122<εthen(15)
returnxk,λk,ϕk(16)
end if(17)
znk=∑ω∈Ωnqnmωxnk,ω,∀n∈η,∀m∈ℐn(18)
λnk+1,ω=λnk,ω+ρxnk,ω−znk,∀n∈η,∀ω∈Ωn(19)
end if(20)
returnxkmax,λkmax,ϕkmax(21)
end function


## 2.1. Scenario Representation

The stochastic program representation can be visualized as a tree, the so-called “scenario tree.” It can be represented as follows. LetT denote the set of periods in the planning horizon with T=T being the number of periods. In the tree, a node represents the realization of the uncertain parameter and variable at a given time period. Let n and N describe the node and the lexicographically numbered set of nodes 1,…,N in the tree, respectively. From each node n, for t∈T∖T, there is at least one branch leading to another node m with a conditional probability pm with m∈ℐn, that is, the set of nodes, immediate successors of the node n. ∑m∈ℐnpm=1. We denote by Nt the subset of nodes belonging to period t such that N=∪t∈TNt and Nt∩Nt+1=Φ for t∈T∖T. Let Ω represent the finite set of representative scenarios in the tree. A scenario ω∈Ω is a particular realization of the uncertain parameter represented as a path from the root-node to a leaf-node. Each scenario ω has an associated probability or weight denoted by wω. Note that ∑ω∈Ωwω=1. Similarly, let Nω represent the set of nodes forming the scenario ω. In other words, Nω is the set of nodes in a path from the root-node to a leaf-node. There is a number of scenarios that traverse each node. Let Ωn denote the set of scenarios that traverse the node n. We have Ω1=Ω, Ωn∩Ωn′=Φ∀n≠n′ and n,n′∈Nt. Finally, let S˜n denote the set of successor nodes to the node n, for n∈N. We have N1 is a singleton and S˜n=Φ for n∈NT (leaf-nodes). Furthermore, let αn denote the immediate ancestor node of node n for n∈N∖1, since n=1 is the root-node of scenario tree. To introduce nonanticipativity, we need to define η as the set of nodes having more than one leaf-node as successor S˜n∩NT>1⟹n∈ηA. For instance, in Figure 1, nodes 2 and 3 have as leaf-node successors the set of nodes {5, 6} and {7}, respectively. Therefore, node 2∈η but node 3 does not. Finally, Let Xω represent the matrix of variables for all t∈T for the scenario ω∈Ω. As the result, the stochastic program can be stated as follows:(1)maxX∑ω∈ΩwωfωX,(2)s.t.Xω∈Cω∀ω∈Ω,(3)Xω=Xω′,∀ω,ω′∈Ωn,ω≠ω′,n∈η.Figure 1
Scenario tree representation of stochastic programming.Equation (1) maximizes the expectation from all scenarios. Equation (2) says that all the solutions should be feasible with respect to the constraints of the scenario. Cω is the feasible set for scenario ω. Finally, equation (3) imposes the nonanticipativity constraints (NACs) which require that the solution up to a period t from two scenarios should be the same if the two scenarios are indistinguishable up to period t. We call this model the extensive form (EF).Notice that we could write constraint (3) in a different form. Let Xnω:=x1nω,x2nω,…,xsnω be the vector of variables at node n∈η under scenario ω, and xsnω is the state variable of the forest unit or “stand” s. Let Zn represent a vector of binary variables at node n. Imposing constraint (4) is equivalent to reinforcing NAC. Progressive hedging exploits the following formulation:(4)Xnω−Zn=0,∀ω∈Ωn,n∈η.


## 2.2. Progressive Hedging

Progressive hedging (PH) is a scenario-based decomposition algorithm proposed by [36] for stochastic programming models. The idea of progressive hedging is to relax the nonanticipativity constraints (NACs, equations (4)) in an augmented Lagrangian manner [18, 22, 35, 36] so that each scenario (subproblem) can be solved independently. This assumes that solving subproblems independently is much easier and faster. Nonanticipativity constraints require that values of variables that share the same ancestor nodes should be equal across scenarios up to that node. In other words, if two scenarios ω and ω′ are indistinguishable up to period t, then the solutions of the two scenarios should be the same up to that period. For instance, in Figure 2, each variable at time t=1 should be the same across all the five scenarios and values of variables in period t=2 should be the same for scenarios 1 and 2. However, these constraints are not reinforced when scenarios are solved independently. Through progressive penalization of violations of NAC, the algorithm is proven to ultimately converge to the optimal solution for convex stochastic programs. In a nutshell, progressive hedging follows the following steps:(1)
Solve each scenario without penalization(2)
Compute the averagez for each variable(3)
If solutions have sufficiently converged, then stop(4)
Update penalization termsλ=ρx−z+λ(5)
Solve scenarios with penalization terms(6)
Go to step 2Figure 2
Example of scenario representation in a compact form (a) and in a detached form (b). The red-dotted ellipses show the nonanticipativity constraints that need to be imposed.
(a)(b)Algorithm1 describes progressive hedging for multistage stochastic programming. The inputs of the algorithm are the penalty factor ρ, the maximum number of iterations kmax, and the termination criterion ϵ which indicates the level of consensus of nonanticipativity constraints that is acceptable. ε=0 means that the algorithm stops if all the nonanticipativity constraints are satisfied. In the algorithm, line 2 initializes a Lagrangian multiplier λ for each NAC constraint. A Lagrangian multiplier is associated with each equation (4). To compute the average, of each variable, we need to have the conditional probability associated to each branch of the scenario tree in the detached form (Figure 2). Equation (5) computes the conditional probability from node n to node m which is an immediate successor. That probability qnmω is proportional to the number of leaf-nodes associated with node m because the number of leaf-nodes informs on the number of scenarios passing by node m. For instance, from Figure 2, q121=q122=p2/2, q251=p5, and q133=p3. Lines 6 and 17 compute the average of each variable at node n∈η. Lines 7 and 18 update the Lagrangian multiplier associated with each nonanticipativity constraint. Line 14 computes the NAC convergence euclidean distance. This distance is zero if all the NACs are satisfied:(5)qnmω=pmΩm,∀m∈ℐn,m∈Nω.Algorithm 1: Progressive hedging multistage.
(1)
function ph ρ,kmax, ε(2)
Initializeλn0,ω=0,∀n∈η,∀ω∈Ωn(3)
forω∈Ωdo(4)
x0,ω=argmaxxfωx:x∈Cω(5)
end for(6)
zn0=∑ω∈Ωnqnmωxn0,ω,∀n∈η,∀m∈ℐn(7)
λn1,ω=λn0,ω+ρxn0,ω−zn0,∀n∈η,∀ω∈Ωn(8)
fork=1 to kmaxdo(9)
forω=1 to Ωdo(10)
xk,ω=argmaxxgωx=fωx−λk,ωx−zk−1−ρ/2x−zk−122:x∈Cω(11)
ϕωk=maxxgωx(12)
end for(13)
ϕk=∑ω∈Ωwωϕωk(14)
if∑ω∈Ωwωxk,ω−zk−122<εthen(15)
returnxk,λk,ϕk(16)
end if(17)
znk=∑ω∈Ωnqnmωxnk,ω,∀n∈η,∀m∈ℐn(18)
λnk+1,ω=λnk,ω+ρxnk,ω−znk,∀n∈η,∀ω∈Ωn(19)
end if(20)
returnxkmax,λkmax,ϕkmax(21)
end function


## 3. Variable Fixation

Progressive hedging variable fixing (PHVF) algorithm is identical to the classic progressive hedging algorithm (Algorithm1). However, instead of letting variables converge progressively, PHVF fixes variables as they converge. Line 10 in the progressive hedging algorithm (Algorithm 1) is replaced by Algorithm 2. The algorithm starts by fixing variables in the root-node as they converge. For instance, for a given node n∈η, if a variable has a value of 1 across all scenarios ω∈Ωn, that variable will be fixed to 1 across all those scenarios. This process is better described in Algorithm 2. However, unless a variable belongs to the root-node, it can only be fixed if all variables in the immediate ancestor αn have been entirely fixed. The first thing to notice is that once the scenario tree is decomposed, there are more scenarios traversing the root-node (Figure 2). In addition, as we move away from the root-node toward the leaf-nodes, there is a fewer number of scenarios passing by any given node. Hence, if for instance, there are three branches originating from the root-node (as displayed in Figure 2), ℐ1=3, then after fixing the root-node, we get three distinct subextensive forms (SEFs) (Ω2,Ω3, and Ω4) that can be solved independently (Figure 3). Furthermore, the three SEFs are much easier to solve compared to the original problem Ω1. The algorithm description exploits that property of the scenario tree.Algorithm 2: Variable fixation module.
(1)
function fixVariables θ,ω,k(2)
stop = False;(3)
fort=1 to Tdo(4)
n:n∈Nt∩η∩Nω(5)
ifαn is fixed OR n is root-node then(6)
fori=1 to number of stands do(7)
ifzn,i≥θtthen(8)
FixXi,tω to 1(9)
else ifzn,i≤1−θtthen(10)
FixXi,tω to 0(11)
else(12)
stop = True(13)
end if(14)
end for(15)
end if(16)
if stop is True then(17)
break(18)
else ifn not marked(19)
Markn as fixed(20)
end if(21)
end for(22)
xk,ω=argmaxxgωx=fωx−λk,ωx−zk−1−ρ/2x−zk−122:x∈C∗ω(23)
return fixed nodes(24)
end functionFigure 3
Structure of stochastic problem after fixing the root-node (node 1). The three subextensive forms (SEFs) denoted byΩ2, Ω3, and Ω4 are independent.


### 3.1. Algorithm Description

During PHVF, at each iteration and for each scenario, the problem online 10 of Algorithm1 needs to be solved. We refer to Cω as the feasible set created by constraint (2) for each scenario. Note that if the uncertainty is only in the objective function (for instance, price uncertainty), then Cω=Cω′,∀ω≠ω′. In this case each solution of ω is feasible for Cω′ω≠ω′. Since z⋅ is the average of the values of the variables across scenarios (see line 6 of Algorithm 1), if its value is 1 (or 0), then it means all the scenarios met consensus that the variable should take a value of 1 (or 0), respectively. This realization is exploited in Algorithm 2 if it is supposed that the slamming factor θt=1,∀t. The algorithm moves to the node at the next stage and performs variable fixation if that node belongs to the set of nonanticipativity nodes η and its predecessor is entirely fixed. The new feasible set created by constraint (2) and fixing some variables is denoted by C∗ω. Note that this subproblem is much smaller compared to the original subproblem. In Algorithm 2, the inputs are θ, a function defining the value of the slamming factor depending on the stage t0≤θt≤1, see Section 3.3 for more information on this parameter; ω which is the scenario; and k which indicates the iteration.Fixing variables impact the optimality since some variables may be fixed to values that they may not have taken in the optimal solution. Hence, the algorithm has an impact on the optimality. The objective of fixing variables is not to completely solve the model through the procedure but to obtain a reduced extensive form (REF) that is tractable. Letτ denote the percentage of NAC variables that are fixed before switching to solving the REF. τ is a proxy to the number of nodes fixed. The question is what is the appropriate τ that will make the REF tractable while not severely impacting the optimality? The answer to this question is investigated in Section 4.4.1.


### 3.2. Parallel Implementation

Algorithm3 is designed for parallel implementation of PHVF. Hence, line 10 from Algorithm 1 has to be replaced by Algorithm 3. At each iteration of progressive hedging, the subproblems can be solved in parallel because they are independent. However, fixing τ variables and switching to solve the reduced extensive form (the extensive form with some variables fixed, REF) is not parallelized. Nonetheless, since our algorithm (PHVF) fixes variables from the root-node to the leaf-nodes, after fixing each node, the problem can be separated into subextensive forms (SEFs) that can be solved independently. For example, as illustrated in Figure 3, after fixing the root-node, we get three separate SEF models. Recursively, each subextensive form model can be solved through variable fixation leading to reduced subextensive form (RSEF). Each RSEF can then be solved directly if its size is deemed tractable. At the end, the solutions from all the REFs can be combined to get the solution of the stochastic program.Algorithm 3: PH variable fixing.
(1)
function PHVF θ,Ωj,k(2)
fixed nodes = fixVaribles(θ,ω,k,∀ω∈Ωj)(3)
if problem small enough then(4)
Solve reduced extensive form.(5)
else(6)
forn in fixed nodes which successor are not fixed do(7)
form∈ℐndo(8)
PHVFθ,Ωm,k+1(9)
end for(10)
end for(11)
end if(12)
returnxω(13)
end function


### 3.3. Acceleration Methods

Note from Algorithm2 that if θt=1, then we fix a variable to 1 (or alternatively to 0) only if all scenarios of interest agree that the value of the variable should be 1 (or 0 alternatively), respectively. This requirement is hard to meet when dealing with hundreds of scenarios, especially for the root-node variables which are replicated in all the scenarios. A variable may have the same value in all but one scenario for many iterations. Hence, that scenario slows the convergence. To avoid such a situation, “slamming” [42] forces variables to converge if the percentage of concordant scenarios for that variable reaches a threshold θ≤1. However, instead of defining a scalar θ as in [42], we defined θt=fθ0,Qt as a function that depends on the stage because if a low value of θ is acceptable for the root-node, such a value is not acceptable when closer to the leaf-nodes (to avoid infeasibility). θ0 is the initial value of θ, whereas Qt can be a linear or exponential function term that increases θ value as a function of the t. Notice that if the uncertainty is in the objective function and we apply slamming, then the problem will always be feasible. In our case, preliminary tests allowed us to find the range of acceptable values of θ0. Low values of θ0 led to infeasibility.In the same line of thought, even with slamming, there is a chance that the algorithm is locked in a situation where there is no improvement of NAC convergence for many consecutive iterations. We define a “cascading” effect which lowers the value ofθ to θmin for one iteration if the algorithm does not improve (convergence of some nonanticipativity) for α consecutive iterations. This behavior allows to avoid getting stuck, since lowering θ allows to fix some variables for which almost all the scenario already reached consensus. For those variables, the percentage of variables that agree on the value the variable should take is close to θt and the consensus could eventually be reached after several iterations. We have noticed that the cascading effect may lead to infeasibility because it may be forcing some variables to a value they would not take otherwise. When infeasibility arises, we roll back and eliminate the cascading effect. In general, there is a tradeoff between infeasibility and the value of θ. Low values of θ lead to a risk of infeasibility while raising that value may make the acceleration methods less efficient.


### 3.4. Stopping Criteria

In the case of classic progressive hedging, the algorithm stops because we have reached an acceptable level of consensus for NAC or because all the NACs are satisfied. The algorithm terminates as well if the maximum number of iterationskmax is reached. In our algorithm, we keep these two stopping criteria, although we know that they may never be reached. We instead rely on the fact that at some points, many nodes will be fixed (NAC consensus is reached for some nodes) and the reduced subextensive forms can be solved directly. The number of NAC consensus reached is checked through the parameter τ which is the percentage of variables fixed. When the percentage of variables fixed is greater than or equal to τ, the algorithm switches to solve the REF or RSEF.


## 3.1. Algorithm Description

During PHVF, at each iteration and for each scenario, the problem online 10 of Algorithm1 needs to be solved. We refer to Cω as the feasible set created by constraint (2) for each scenario. Note that if the uncertainty is only in the objective function (for instance, price uncertainty), then Cω=Cω′,∀ω≠ω′. In this case each solution of ω is feasible for Cω′ω≠ω′. Since z⋅ is the average of the values of the variables across scenarios (see line 6 of Algorithm 1), if its value is 1 (or 0), then it means all the scenarios met consensus that the variable should take a value of 1 (or 0), respectively. This realization is exploited in Algorithm 2 if it is supposed that the slamming factor θt=1,∀t. The algorithm moves to the node at the next stage and performs variable fixation if that node belongs to the set of nonanticipativity nodes η and its predecessor is entirely fixed. The new feasible set created by constraint (2) and fixing some variables is denoted by C∗ω. Note that this subproblem is much smaller compared to the original subproblem. In Algorithm 2, the inputs are θ, a function defining the value of the slamming factor depending on the stage t0≤θt≤1, see Section 3.3 for more information on this parameter; ω which is the scenario; and k which indicates the iteration.Fixing variables impact the optimality since some variables may be fixed to values that they may not have taken in the optimal solution. Hence, the algorithm has an impact on the optimality. The objective of fixing variables is not to completely solve the model through the procedure but to obtain a reduced extensive form (REF) that is tractable. Letτ denote the percentage of NAC variables that are fixed before switching to solving the REF. τ is a proxy to the number of nodes fixed. The question is what is the appropriate τ that will make the REF tractable while not severely impacting the optimality? The answer to this question is investigated in Section 4.4.1.


## 3.2. Parallel Implementation

Algorithm3 is designed for parallel implementation of PHVF. Hence, line 10 from Algorithm 1 has to be replaced by Algorithm 3. At each iteration of progressive hedging, the subproblems can be solved in parallel because they are independent. However, fixing τ variables and switching to solve the reduced extensive form (the extensive form with some variables fixed, REF) is not parallelized. Nonetheless, since our algorithm (PHVF) fixes variables from the root-node to the leaf-nodes, after fixing each node, the problem can be separated into subextensive forms (SEFs) that can be solved independently. For example, as illustrated in Figure 3, after fixing the root-node, we get three separate SEF models. Recursively, each subextensive form model can be solved through variable fixation leading to reduced subextensive form (RSEF). Each RSEF can then be solved directly if its size is deemed tractable. At the end, the solutions from all the REFs can be combined to get the solution of the stochastic program.Algorithm 3: PH variable fixing.
(1)
function PHVF θ,Ωj,k(2)
fixed nodes = fixVaribles(θ,ω,k,∀ω∈Ωj)(3)
if problem small enough then(4)
Solve reduced extensive form.(5)
else(6)
forn in fixed nodes which successor are not fixed do(7)
form∈ℐndo(8)
PHVFθ,Ωm,k+1(9)
end for(10)
end for(11)
end if(12)
returnxω(13)
end function


## 3.3. Acceleration Methods

Note from Algorithm2 that if θt=1, then we fix a variable to 1 (or alternatively to 0) only if all scenarios of interest agree that the value of the variable should be 1 (or 0 alternatively), respectively. This requirement is hard to meet when dealing with hundreds of scenarios, especially for the root-node variables which are replicated in all the scenarios. A variable may have the same value in all but one scenario for many iterations. Hence, that scenario slows the convergence. To avoid such a situation, “slamming” [42] forces variables to converge if the percentage of concordant scenarios for that variable reaches a threshold θ≤1. However, instead of defining a scalar θ as in [42], we defined θt=fθ0,Qt as a function that depends on the stage because if a low value of θ is acceptable for the root-node, such a value is not acceptable when closer to the leaf-nodes (to avoid infeasibility). θ0 is the initial value of θ, whereas Qt can be a linear or exponential function term that increases θ value as a function of the t. Notice that if the uncertainty is in the objective function and we apply slamming, then the problem will always be feasible. In our case, preliminary tests allowed us to find the range of acceptable values of θ0. Low values of θ0 led to infeasibility.In the same line of thought, even with slamming, there is a chance that the algorithm is locked in a situation where there is no improvement of NAC convergence for many consecutive iterations. We define a “cascading” effect which lowers the value ofθ to θmin for one iteration if the algorithm does not improve (convergence of some nonanticipativity) for α consecutive iterations. This behavior allows to avoid getting stuck, since lowering θ allows to fix some variables for which almost all the scenario already reached consensus. For those variables, the percentage of variables that agree on the value the variable should take is close to θt and the consensus could eventually be reached after several iterations. We have noticed that the cascading effect may lead to infeasibility because it may be forcing some variables to a value they would not take otherwise. When infeasibility arises, we roll back and eliminate the cascading effect. In general, there is a tradeoff between infeasibility and the value of θ. Low values of θ lead to a risk of infeasibility while raising that value may make the acceleration methods less efficient.


## 3.4. Stopping Criteria

In the case of classic progressive hedging, the algorithm stops because we have reached an acceptable level of consensus for NAC or because all the NACs are satisfied. The algorithm terminates as well if the maximum number of iterationskmax is reached. In our algorithm, we keep these two stopping criteria, although we know that they may never be reached. We instead rely on the fact that at some points, many nodes will be fixed (NAC consensus is reached for some nodes) and the reduced subextensive forms can be solved directly. The number of NAC consensus reached is checked through the parameter τ which is the percentage of variables fixed. When the percentage of variables fixed is greater than or equal to τ, the algorithm switches to solve the REF or RSEF.


## 4. Numerical Experiment

We describe here an empirical performance analysis of our proposed algorithm. For easiness to follow our experiment and replicate the results, we formulate the stochastic version of the so-called “Model I” of forest harvest scheduling.


### 4.1. Problem Definition and Formulation

Indicesω,ω′: scenarios: standt: time periodSetsS: set of standsΩ: set of scenariosΩn: set of scenarios passing by the node nη: set of nodes on which NAC should be reinforcedNt: set of nodes at stage tT: set of time in the planning horizonParametersrt: profit from selling wood in period t ($/mbf). It is the discounted profit that includes selling costcst: cost of harvesting and hauling wood from stand s in period t ($). It is the discounted costystω: volume of wood harvestable per area from stand s in period t according to scenario ω (mbf/ac). This volume depends on the climate that materializes; therefore, it is a parameter that depends on the scenario ωas: area of stand s (acres)agest: age of stand s at the end of the planning horizon if harvested in year tages·: current age of stand sages0: age of stand s if not harvested during the planning horizonfmin: allowable percentage of decrease of volume harvested from one period to anotherfmax: allowable percentage of increase of volume harvested from one period to anotherwω: probability or weight of scenario ωVariablesxstω: binary variable taking a value of 1 if stand s should be harvested in period t according to scenario ω, and 0 otherwisensω: binary variable: 1 if management unit s should not be harvested during the planning horizon under scenario ω, and 0 otherwiseHtω: accounting variable storing the volume of wood harvested in period t according to scenario ω (mbf):(6)max∑ω∈Ωwω∑t∈TrtHtω−∑s∈Scstxstω,which is subject to(7)nsω+∑t∈Txstω=1,∀s∈S,∀ω∈Ω,(8)∑s∈Sasystωxstω=Htω,∀t∈T,∀ω∈Ω,(9)1−fminHtω≥Ht+1ω,∀ω∈Ω,∀t=1,…,T−1,(10)1+fmaxHtω≤Ht+1ω,∀ω∈Ω,∀t=1,…,T−1,(11)∑s∈Sas∑t∈Tagestxstω+ages0nsω≥∑sasages.,∀ω∈Ω,(12)xstω=xstω′,∀ω≠ω′;ω,ω′∈Ωn;n∈η;t=t|n∈Nt;∀s∈S,(13)xstω,nsω=0,1,Htω∈R+.Expression (6) maximizes the profit from timber harvest from all scenarios weighted by their respective probabilities. Constraints (7) require that a stand is at most harvested once during the planning horizon. We use the variable ns to capture stands that are not prescribed to be harvested during the whole planning horizon. We need that variable to compute the average ending age of the forest as shown in constraints (11). Constraints (8) compute the volume harvested in each period during the planning horizon. It uses the parameter ystω which values depend on the forest growth scenario ω of interest. Note that the set of constraints (8) is not necessary. We could have written the same model without using that set of constraints. However, doing so would require rewriting constraints (9) and (10), and finally, it would negatively affect the readability of the model. Constraints (9) and (10) impose the even flow constraints so that the volume harvested from one period to another remain within an allowed fluctuation range. Constraints (11) require that on average, the forest at the end of the planning horizon is at least as old as the forest at the beginning of the planning horizon. Finally, constraints (12) impose the nonanticipativity constraints. It requires that for each stand s, at time t, if two scenarios are indistinguishable, then the decision should be the same for the two scenarios at that time. If Nω∩Nω′≠Φ, then there exists t such that the two scenarios are indistinguishable at t. Constraints (13) enforce that the decision variables are binary, and the accounting variables are continuous. We remind the reader that variables Htω are not required for this model.


### 4.2. Climate Change Data

The potential mean annual increment, which is an indicator of forest growth in a year might change in the Pacific North West because of climate change. The change depends on temperature, precipitation, and air moisture content, all driven by human activities and economic development [28]. In addition, the change is not geographically uniform. Hence, the change tends to be negative in Oregon compared to Washington State. Similarly, the change tends to be negative in low altitudes in Oregon. We used the data from [28] which forecast the potential mean annual increment change (pMAI) by the year 2100. The pMAI m3/ha/year is the potential change in forest growth that will be observed in a given year; hence, it is a volume given as a function of time and area. We assumed linear change of the growth from now until that year. For instance, the growth change for next two decades is the double of the growth change for the next decade. The climate paths defined in the Pacific Northwest are A2, A1B, B1, and Commit (or C). Tables 1 and 2 present the values of potential mean annual increment change for each one of the climate paths. In Table 2, climate paths D1 to D4 are artificial climate paths that suppose higher pMAI. We built forest growth scenarios by assuming it is possible to transit from one climate path to another because of mitigation or intensification of climate change due to human actions.Table 1
Value of pMAIm3/ha/year for small and medium instances.
ForestsP1P34P36P75P83P100A20.53.1−0.31.80.42.1A1B0.42.5−0.31.60.31.6B10.41.30.01.00.31.0C0.20.20.10.20.10.3Table 2
Values of pMAIm3/ha/year of all the six forests for large instances.
A2A1BB1CD1D2D3D43.12.51.30.24.05.06.07.0


### 4.3. Experimental Design

The experiments were conducted under a DELL desktop computer running on Windows with Intel(R) Core(TM) 2 Quad CPU @ 3.70 GHz and 8 GB of memory. During PHVF iterations, each scenario is solved at optimality gap of 1%. This is a premature stop. Nevertheless, this criterion proved to accelerate solution time for each scenario. Furthermore, the last solutions of a mixed integer program are the most difficult ones with no to little improvement of the objective function value. All the models were solved using IBM ILOG CPLEX 12.6 (CPLEX, [12]). The memory allocated for storing the nodes was 3,000 MB, and the nodes were set to be stored in a compressed format on the hard drive. All other parameters were left to their default values. The code was implemented in Java 10 using Concert Technology of CPLEX.For parallel computation, we used the paradigm of master-workers. The master is in charge of coordinating the PHVF algorithm while distributing the task of solving individual submodels or REF to the workers. Each worker sends back its solution upon completion. The workers compete for access to the memory. Therefore, the choice of the number of workers must be judicious to avoid the overhead which is the amount of time required to coordinate parallel tasks, as opposed to doing useful work. Preliminary experiments showed that two workers was the optimal number given the configurations of the computer. The general framework of the model is presented in Figure4. The input data contain the information on the forest, the climate change (growth change), and the regulations that ought to be met by the harvest planning. The input is fed to the PHVF module that is the master governing the optimization process. The workers are the ones interacting with the optimizer (CPLEX in this case). They send the model to the optimizer which returns the solution. The advantage of this framework is that we could change the optimizer without having to readapt our algorithm. At the end, the output module collects the harvest planning and informs the decision making.Figure 4
Framework of the PHVF algorithm for harvest scheduling under climate change.Forest growth scenarios were generated using the information on the growth change reported in [28] (see Tables 1 and 2). We assumed that due to climate change and climate change mitigation efforts, it is possible to transition from one climate path to another in two consecutive periods. Hence, for instance, it is possible to transition from climate path A2 in year 2020 to climate path B1 in year 2030. We further assume that the probability of transiting from one climate path to the other is the same. As result, the scenarios are considered equally probable. Note that this assumption does not have much incidence on the performance of the algorithm. Based on the number of scenarios generated, we defined small, medium, and big instances. The small instances have 64 scenarios (1 × 4 × 4 × 4 × 1) which is four branching for the first three periods. The medium instances have 256 scenarios which correspond to four branching for the first four periods (1 × 4 × 4 × 4 × 4). For small and medium instances, pMAI corresponding to each climate path is reported in Table 1. The big instances have 512 scenarios made of eight branching for the first three periods (1 × 8 × 8 × 8 × 1). The in-depth description of the scenario tree generation methods and its quality are beyond the scope of this research. The interested reader could refer to [11, 14, 21, 24, 25, 29, 33] and more importantly [34] who describe scenario generation in the forest management framework. We tested our algorithm on six different forests with three being real forests, and the data of which are publicly available1 (P1, P34, P36). The other three were computer-generated forests (P75, P83, P100). The six forests are set to be located in three different altitude classes. The forests are supposed to be located in Oregon and in Washington State. The values of pMAI for the six forests for large instances are reported in Table 2. Table 3 shows the characteristics of the forests and the scenario instances. In the tables, values show the potential mean annual increment change by 2100. Hence, negative values mean there will be a decrease in forest growth, whereas positive values mean that there will be an increase in forest growth.Table 3
Description of problem size.
ForestΩBinary colsTotal colsRowsNonzerosStandsP1646,552,0156,552,3356,349,76061,375,4881,36325662,209,42362,210,70361,747,96866,199,80851264,286,63564,289,19562,798,592614,320,640P3464612,960613,28069,024633,47232256651,872653,152644,2886,150,2725126,100,6406,103,200672,7046,322,560P3664636,045636,365623,616696,576892566,144,2696,145,5496,117,2486,431,8725126,279,9056,282,4656,189,4406,920,576P75646,121,5006,121,820677,6326,311,8083002566,486,3006,487,5806,387,32861,400,8325126,943,5006,946,0606,621,56863,288,576P83646,121,5006,121,820677,6326,311,1043002566,486,3006,487,5806,387,32861,398,0165126,943,5006,946,0606,621,56863,277,312P100646,121,5006,121,820677,6326,310,4003002566,486,3006,487,5806,387,32861,395,2005126,943,5006,946,0606,621,56863,266,048The medium instances were used as reference to investigate how many variables should be fixed in order to have a problem that is computationally tractable. For that purpose, we solved the stochastic models, using the serial version of our algorithm (Algorithms1 and 2), setting τ to just 20%, 40%, 60%, and 70% before switching to solving the reduced extensive form. τ=20% means some variables are fixed but the root-node is not completely fixed, and τ=40% means that all the root-node variables are completely fixed and some variables from the second period nodes are fixed (as stated before, τ is proxy to the number of nodes fixed. This scheme is possible because there are 25% of NAC variables at each stage and thus at the root-node as well).We made the experiment by allocating 15, 30, 45, and 60 minutes for each problem. The time is chosen to allow the effect ofτ to manifest itself. The tests were run with five repetitions, and the average of the objective function value was reported. Notice that if the time is too long, then eventually all problems can be solved, even the extensive form. Conversely, if the time is too short, then PHVF may not have enough time to finish fixing variables before the time limit and hence low values of τ will be favored. During the whole experiment, we defined θt=min0.999,1.05t−1θ0. As results, for the first stage t=1, θ1=θ0. Similarly, we defined θmin=max0.75,θt−0.05t. These values were defined from empirical experimentation.For assessment of our algorithm performance, each problem was solved using our algorithm (PHVF) and comparing that solution to the one obtained by CPLEX solving directly the EF for the equivalent wall clock time. In addition, we report the optimality gap produced by CPLEX and the gain which is calculated according togain=PHVF−EF/EF where EF is the objective function value obtained by solving directly the extensive form model and PHVF is the solution obtained by PHVF algorithm for the same runtime. For medium and large instances, we report as well the EF value after solving the EF problems for 86,000 seconds (24 hours).


### 4.4. Results of the Experiment

In addition to the number of scenarios, SMIP size grows as well with the number of units (stands) because we need to define the nonanticipativity constraint for each stand at the time periods on which constraints (12) should be imposed. The problem size increases almost by tenfold when going from small instances to the medium ones. The smallest problem in this experiment has 600k + binary variables and over 69k constraints (Table 3).


#### 4.4.1. Impact ofτ on the Optimality and Solution Time

Tables4–9 report the effect of τ on the value of the objective function with respect to the runtime for forests P1, P34, P36, P75, P83, and P100, respectively. We report as well, improvement% which is the percentage of increase (if positive) or decrease (if negative) of the objective function value from the objective function value for τ=20% for the same runtime.Table 4
Effects ofτ and allocated time for the forest P1.
ForestρTime (min)τIterationsαθ0ValueImprovement (%)P110000001520160.95371,524,817—40960.95386,068,4333.9501860.95TNA602060.95TNA702160.95TNA3020160.95373,450,417—40960.95386,156,7793.4501860.95385,676,1183.3602060.95385,394,6863.2702160.95385,258,7713.24520160.95375,427,413—40960.95386,156,7792.9501860.95385,676,1182.7602060.95385,394,6862.7702160.95385,258,7712.66020160.95375,804,390—40960.95386,156,7792.8501860.95385,676,1182.6602060.95385,394,6862.6702160.95385,258,7712.5T: termination because PHVF ran out of time while iterating; NA: not applicable;ρ: penalty factor; τ: percentage of variables fixed; α: cascading factor; θ0: initial slamming factor.Table 5
Effects ofτ and allocated time for the forest P34.
ForestρTime (min)τIterationsαθ0ValueImprovement (%)P341000015201100.951,333,129—401100.951,332,692−0.0506100.95NNA609100.95NNA7018100.951,347,1521.130201100.951,344,498—401100.951,344,260−0.0506100.95NNA609100.95NNA7018100.951,347,4910.245201100.951,345,620—401100.951,345,7610.0506100.95NNA609100.95NNA7018100.951,347,5550.160201100.951,346,799—401100.951,346,489−0.0506100.951,343,802−0.2609100.95NNA7018100.951,347,7590.1N: terminated because could not find feasible solution; NA: not applicable;ρ: penalty factor; τ: percentage of variables fixed; α: cascading factor; θ0: initial slamming factor.Table 6
Effects ofτ and allocated time for the forest P36.
ForestρTime (min)τIterationsαθ0ValueImprovement (%)P36500015201100.94,468,454—4022100.94,513,7171.05041100.94,532,8191.46051100.94,522,0131.27060100.94,525,6691.330201100.94,570,040—4022100.94,534,327−0.85041100.94,534,156−0.86051100.94,523,000−1.07060100.94,530,669−0.945201100.94,616,899—4022100.94,534,470−1.85041100.94,534,161−1.86051100.94,532,251−1.87060100.94,530,669−1.960201100.94,662,056–4022100.94,534,574−2.75041100.94,534,219−2.76051100.94,532,251−2.87060100.94,530,669−2.8ρ: penalty factor; τ: percentage of variables fixed; α: cascading factor; θ0: initial slamming factor.Table 7
Effects ofτ and allocated time for the forest P75.
ForestρTime (min)τIterationsαθ0ValueImprovement (%)P75200015201100.958,202,734—4019100.959,007,4031.45055100.9TNA6057100.9TNA7070100.9TNA30201100.958,272,824—4019100.959,106,2771.45055100.959,387,8841.96057100.959,324,5441.87070100.959,337,7821.845201100.958,297,317—4019100.959,131,6011.45055100.959,394,9421.96057100.959,336,5921.87070100.959,354,2491.860201100.958,310,384—4019100.959,248,4441.65055100.959,395,8161.96057100.959,339,1451.87070100.959,354,5251.8T: termination because run out of time iterating; NA: not applicable;ρ: penalty factor; τ: percentage of variables fixed; α: cascading factor; θ0: initial slamming factor.Table 8
Effects ofτ and allocated time for the forest P83.
ForestρTime (min)τIterationsαθ0ValueImprovement (%)P83200015201100.957,970,162—4014100.958,558,1701.05030100.958,516,9760.96032100.958,559,0691.07034100.958,592,5621.130201100.958,305,853—4014100.958,891,3811.05030100.958,833,3170.96032100.958,733,5630.77034100.958,661,3470.645201100.959,052,205—4014100.958,898,935−0.35030100.958,834,898−0.46032100.958,733,746−0.57034100.958,661,964−0.760201100.959,052,215—4014100.958,904,481−0.35030100.958,835,134−0.46032100.958,733,734−0.57034100.958,661,851−0.7ρ: penalty factor; τ: percentage of variables fixed; α: cascading factor; θ0: initial slamming factor.Table 9
Effects ofτ and allocated time for the forest P100.
ForestρTime (min)τIterationsαθ0ValueImprovement (%)P1005000015201100.9542,331,212—4027100.95TNA5063100.95TNA6069100.95TNA7079100.95TNA30201100.9542,494,130—4027100.9543,307,7441.95063100.9543,298,4391.96069100.95NNA7079100.9543,288,8411.945201100.9542,819,969—4027100.9543,309,4551.15063100.9543,299,1891.16069100.95NNA7079100.9543,290,1271.160201100.9542,931,057—4027100.9543,309,4020.95063100.9543,299,7590.96069100.95NNA7079100.9543,290,1270.8T: termination because run out of time iterating; N: terminated because could not find feasible solution; NA: not applicable,ρ: penalty factor; τ: percentage of variables fixed; α: cascading factor; θ0: initial slamming factor.In summary, as expected, everything else being equal, longer runtimes allow a higher objective function. As we can see forτ>40% which corresponds to fixing entirely the root-node and fixing some variables from the second stage, we have numerical issues. The numerical issues have two sources. On one hand, the algorithm may run out of time while fixing the variables or finished fixing the variables but does not have enough time to find a feasible solution to the REF. On the other hand, the basis of the REF may be disturbed in a way that the number of constraints and variables is not balanced so it is hard for branch and cut algorithm used by CPLEX to find an integral solution for the REF in the allocated time. This is the case for forest P34 with no solution when τ=60% even when the time is 1 hour. However, for the same forest, there is an integral solution when τ is raised to 70%. This behavior occurred mainly for forest P34 and P100 and is not observed when τ=40%.The second aspect is the impact ofτ on the optimality. The advantage that higher values of τ have over the small ones tends to vanish when the runtime is long. Similarly, the improvement in optimality from τ=40% to higher values of τ is low. In fact, for some forest such as P83, over τ=40%, the objective function value diminishes as τ increases which can be interpreted as during PHVF, many variables are fixed to values that are suboptimal. It clearly appears that fixing just the root-node is the best choice for these problems. It allows sufficient time for REF to find an integral solution while avoiding to impact the optimality, and limiting the disturbance on the structure of the basis matrix.


#### 4.4.2. PHVF versus EF Solved Directly

Table10 presents the results from PHVF using τ=40% which is equivalent to fixing the root-node and then solving in parallel the SEF. As we can see for the small instances (64 scenarios), solving directly the extensive form outperformed using PHVF. This means the overhead of decomposing the problem into scenarios and solving each one overweighs the benefit. For the medium instances, the results are mixed. Out of the six forests, PHVF outperformed the EF in two cases. Even when PHVF under performed, it had results that were less than 1% away from the optimal solution. The benefit of using PHVF over the state of art commercial solver such as CPLEX is highlighted when dealing with big instances (512 scenarios). PHVF outperformed the EF for the six forests. In the case of P1 which is the biggest model with 1,363 units (Table 3), EF completely failed to solve the model because it could not fit it in the memory. For other cases for the same runtime, PHVF reached gains ranging from 0.84% to over 767% corresponding to forests P34 and P100, respectively. Even after leaving EF for 24 h, PHVF run for 15 min still outperformed in some cases. Comparing the gain to the gap from EF suggests that PHVF almost reached the optimal solution or its solution has a relative optimal gap less than 1%.Table 10
Comparison of solutions from PHVF and the extensive form.
Same duration with PHVF24 hΩTime (s)ForestPHVFEFGap (%)Gain (%)EFGap (%)Gain (%)64120P12.99E + 083.08E + 080.14−2.69———64120P341.04E + 061.05E + 060.14−1.64———64120P363.48E + 063.73E + 060.32−6.77———64120P754.71E + 074.75E + 071.59−0.81———64120P834.70E + 074.72E + 070.43−0.46———64120P1006.53E + 076.54E + 070.51−0.05———256600P13.86E + 083.87E + 080.05−0.173.87E + 080.02−0.28256600P341.35E + 061.33E + 062.381.561.35E + 060.73−0.05256600P364.73E + 064.64E + 062.051.754.70E + 060.820.52256600P755.94E + 075.94E + 070.69−0.065.97E + 070.19−0.55256600P835.90E + 075.93E + 070.28−0.475.93E + 070.22−0.50256600P1004.34E + 074.34E + 070.34−0.064.35E + 070.14−0.22512900P12.91E + 08NANA∞NANA∞512900P342.07E + 062.05E + 061.20.842.07E + 060.27-0.08512900P366.85E + 066.78E + 061.8761.136.79E + 061.640.89512900P758.48E + 074.50E + 0789.5788.638.50E + 070.22−0.28512900P838.47E + 074.49E + 0789.4488.768.45E + 070.630.27512900P1004.74E + 075.49E + 06767.9763.35.49E + 06767.9763.3NA: no results because of memory limitation.


## 4.1. Problem Definition and Formulation

Indicesω,ω′: scenarios: standt: time periodSetsS: set of standsΩ: set of scenariosΩn: set of scenarios passing by the node nη: set of nodes on which NAC should be reinforcedNt: set of nodes at stage tT: set of time in the planning horizonParametersrt: profit from selling wood in period t ($/mbf). It is the discounted profit that includes selling costcst: cost of harvesting and hauling wood from stand s in period t ($). It is the discounted costystω: volume of wood harvestable per area from stand s in period t according to scenario ω (mbf/ac). This volume depends on the climate that materializes; therefore, it is a parameter that depends on the scenario ωas: area of stand s (acres)agest: age of stand s at the end of the planning horizon if harvested in year tages·: current age of stand sages0: age of stand s if not harvested during the planning horizonfmin: allowable percentage of decrease of volume harvested from one period to anotherfmax: allowable percentage of increase of volume harvested from one period to anotherwω: probability or weight of scenario ωVariablesxstω: binary variable taking a value of 1 if stand s should be harvested in period t according to scenario ω, and 0 otherwisensω: binary variable: 1 if management unit s should not be harvested during the planning horizon under scenario ω, and 0 otherwiseHtω: accounting variable storing the volume of wood harvested in period t according to scenario ω (mbf):(6)max∑ω∈Ωwω∑t∈TrtHtω−∑s∈Scstxstω,which is subject to(7)nsω+∑t∈Txstω=1,∀s∈S,∀ω∈Ω,(8)∑s∈Sasystωxstω=Htω,∀t∈T,∀ω∈Ω,(9)1−fminHtω≥Ht+1ω,∀ω∈Ω,∀t=1,…,T−1,(10)1+fmaxHtω≤Ht+1ω,∀ω∈Ω,∀t=1,…,T−1,(11)∑s∈Sas∑t∈Tagestxstω+ages0nsω≥∑sasages.,∀ω∈Ω,(12)xstω=xstω′,∀ω≠ω′;ω,ω′∈Ωn;n∈η;t=t|n∈Nt;∀s∈S,(13)xstω,nsω=0,1,Htω∈R+.Expression (6) maximizes the profit from timber harvest from all scenarios weighted by their respective probabilities. Constraints (7) require that a stand is at most harvested once during the planning horizon. We use the variable ns to capture stands that are not prescribed to be harvested during the whole planning horizon. We need that variable to compute the average ending age of the forest as shown in constraints (11). Constraints (8) compute the volume harvested in each period during the planning horizon. It uses the parameter ystω which values depend on the forest growth scenario ω of interest. Note that the set of constraints (8) is not necessary. We could have written the same model without using that set of constraints. However, doing so would require rewriting constraints (9) and (10), and finally, it would negatively affect the readability of the model. Constraints (9) and (10) impose the even flow constraints so that the volume harvested from one period to another remain within an allowed fluctuation range. Constraints (11) require that on average, the forest at the end of the planning horizon is at least as old as the forest at the beginning of the planning horizon. Finally, constraints (12) impose the nonanticipativity constraints. It requires that for each stand s, at time t, if two scenarios are indistinguishable, then the decision should be the same for the two scenarios at that time. If Nω∩Nω′≠Φ, then there exists t such that the two scenarios are indistinguishable at t. Constraints (13) enforce that the decision variables are binary, and the accounting variables are continuous. We remind the reader that variables Htω are not required for this model.


## 4.2. Climate Change Data

The potential mean annual increment, which is an indicator of forest growth in a year might change in the Pacific North West because of climate change. The change depends on temperature, precipitation, and air moisture content, all driven by human activities and economic development [28]. In addition, the change is not geographically uniform. Hence, the change tends to be negative in Oregon compared to Washington State. Similarly, the change tends to be negative in low altitudes in Oregon. We used the data from [28] which forecast the potential mean annual increment change (pMAI) by the year 2100. The pMAI m3/ha/year is the potential change in forest growth that will be observed in a given year; hence, it is a volume given as a function of time and area. We assumed linear change of the growth from now until that year. For instance, the growth change for next two decades is the double of the growth change for the next decade. The climate paths defined in the Pacific Northwest are A2, A1B, B1, and Commit (or C). Tables 1 and 2 present the values of potential mean annual increment change for each one of the climate paths. In Table 2, climate paths D1 to D4 are artificial climate paths that suppose higher pMAI. We built forest growth scenarios by assuming it is possible to transit from one climate path to another because of mitigation or intensification of climate change due to human actions.Table 1
Value of pMAIm3/ha/year for small and medium instances.
ForestsP1P34P36P75P83P100A20.53.1−0.31.80.42.1A1B0.42.5−0.31.60.31.6B10.41.30.01.00.31.0C0.20.20.10.20.10.3Table 2
Values of pMAIm3/ha/year of all the six forests for large instances.
A2A1BB1CD1D2D3D43.12.51.30.24.05.06.07.0


## 4.3. Experimental Design

The experiments were conducted under a DELL desktop computer running on Windows with Intel(R) Core(TM) 2 Quad CPU @ 3.70 GHz and 8 GB of memory. During PHVF iterations, each scenario is solved at optimality gap of 1%. This is a premature stop. Nevertheless, this criterion proved to accelerate solution time for each scenario. Furthermore, the last solutions of a mixed integer program are the most difficult ones with no to little improvement of the objective function value. All the models were solved using IBM ILOG CPLEX 12.6 (CPLEX, [12]). The memory allocated for storing the nodes was 3,000 MB, and the nodes were set to be stored in a compressed format on the hard drive. All other parameters were left to their default values. The code was implemented in Java 10 using Concert Technology of CPLEX.For parallel computation, we used the paradigm of master-workers. The master is in charge of coordinating the PHVF algorithm while distributing the task of solving individual submodels or REF to the workers. Each worker sends back its solution upon completion. The workers compete for access to the memory. Therefore, the choice of the number of workers must be judicious to avoid the overhead which is the amount of time required to coordinate parallel tasks, as opposed to doing useful work. Preliminary experiments showed that two workers was the optimal number given the configurations of the computer. The general framework of the model is presented in Figure4. The input data contain the information on the forest, the climate change (growth change), and the regulations that ought to be met by the harvest planning. The input is fed to the PHVF module that is the master governing the optimization process. The workers are the ones interacting with the optimizer (CPLEX in this case). They send the model to the optimizer which returns the solution. The advantage of this framework is that we could change the optimizer without having to readapt our algorithm. At the end, the output module collects the harvest planning and informs the decision making.Figure 4
Framework of the PHVF algorithm for harvest scheduling under climate change.Forest growth scenarios were generated using the information on the growth change reported in [28] (see Tables 1 and 2). We assumed that due to climate change and climate change mitigation efforts, it is possible to transition from one climate path to another in two consecutive periods. Hence, for instance, it is possible to transition from climate path A2 in year 2020 to climate path B1 in year 2030. We further assume that the probability of transiting from one climate path to the other is the same. As result, the scenarios are considered equally probable. Note that this assumption does not have much incidence on the performance of the algorithm. Based on the number of scenarios generated, we defined small, medium, and big instances. The small instances have 64 scenarios (1 × 4 × 4 × 4 × 1) which is four branching for the first three periods. The medium instances have 256 scenarios which correspond to four branching for the first four periods (1 × 4 × 4 × 4 × 4). For small and medium instances, pMAI corresponding to each climate path is reported in Table 1. The big instances have 512 scenarios made of eight branching for the first three periods (1 × 8 × 8 × 8 × 1). The in-depth description of the scenario tree generation methods and its quality are beyond the scope of this research. The interested reader could refer to [11, 14, 21, 24, 25, 29, 33] and more importantly [34] who describe scenario generation in the forest management framework. We tested our algorithm on six different forests with three being real forests, and the data of which are publicly available1 (P1, P34, P36). The other three were computer-generated forests (P75, P83, P100). The six forests are set to be located in three different altitude classes. The forests are supposed to be located in Oregon and in Washington State. The values of pMAI for the six forests for large instances are reported in Table 2. Table 3 shows the characteristics of the forests and the scenario instances. In the tables, values show the potential mean annual increment change by 2100. Hence, negative values mean there will be a decrease in forest growth, whereas positive values mean that there will be an increase in forest growth.Table 3
Description of problem size.
ForestΩBinary colsTotal colsRowsNonzerosStandsP1646,552,0156,552,3356,349,76061,375,4881,36325662,209,42362,210,70361,747,96866,199,80851264,286,63564,289,19562,798,592614,320,640P3464612,960613,28069,024633,47232256651,872653,152644,2886,150,2725126,100,6406,103,200672,7046,322,560P3664636,045636,365623,616696,576892566,144,2696,145,5496,117,2486,431,8725126,279,9056,282,4656,189,4406,920,576P75646,121,5006,121,820677,6326,311,8083002566,486,3006,487,5806,387,32861,400,8325126,943,5006,946,0606,621,56863,288,576P83646,121,5006,121,820677,6326,311,1043002566,486,3006,487,5806,387,32861,398,0165126,943,5006,946,0606,621,56863,277,312P100646,121,5006,121,820677,6326,310,4003002566,486,3006,487,5806,387,32861,395,2005126,943,5006,946,0606,621,56863,266,048The medium instances were used as reference to investigate how many variables should be fixed in order to have a problem that is computationally tractable. For that purpose, we solved the stochastic models, using the serial version of our algorithm (Algorithms1 and 2), setting τ to just 20%, 40%, 60%, and 70% before switching to solving the reduced extensive form. τ=20% means some variables are fixed but the root-node is not completely fixed, and τ=40% means that all the root-node variables are completely fixed and some variables from the second period nodes are fixed (as stated before, τ is proxy to the number of nodes fixed. This scheme is possible because there are 25% of NAC variables at each stage and thus at the root-node as well).We made the experiment by allocating 15, 30, 45, and 60 minutes for each problem. The time is chosen to allow the effect ofτ to manifest itself. The tests were run with five repetitions, and the average of the objective function value was reported. Notice that if the time is too long, then eventually all problems can be solved, even the extensive form. Conversely, if the time is too short, then PHVF may not have enough time to finish fixing variables before the time limit and hence low values of τ will be favored. During the whole experiment, we defined θt=min0.999,1.05t−1θ0. As results, for the first stage t=1, θ1=θ0. Similarly, we defined θmin=max0.75,θt−0.05t. These values were defined from empirical experimentation.For assessment of our algorithm performance, each problem was solved using our algorithm (PHVF) and comparing that solution to the one obtained by CPLEX solving directly the EF for the equivalent wall clock time. In addition, we report the optimality gap produced by CPLEX and the gain which is calculated according togain=PHVF−EF/EF where EF is the objective function value obtained by solving directly the extensive form model and PHVF is the solution obtained by PHVF algorithm for the same runtime. For medium and large instances, we report as well the EF value after solving the EF problems for 86,000 seconds (24 hours).


## 4.4. Results of the Experiment

In addition to the number of scenarios, SMIP size grows as well with the number of units (stands) because we need to define the nonanticipativity constraint for each stand at the time periods on which constraints (12) should be imposed. The problem size increases almost by tenfold when going from small instances to the medium ones. The smallest problem in this experiment has 600k + binary variables and over 69k constraints (Table 3).


### 4.4.1. Impact ofτ on the Optimality and Solution Time

Tables4–9 report the effect of τ on the value of the objective function with respect to the runtime for forests P1, P34, P36, P75, P83, and P100, respectively. We report as well, improvement% which is the percentage of increase (if positive) or decrease (if negative) of the objective function value from the objective function value for τ=20% for the same runtime.Table 4
Effects ofτ and allocated time for the forest P1.
ForestρTime (min)τIterationsαθ0ValueImprovement (%)P110000001520160.95371,524,817—40960.95386,068,4333.9501860.95TNA602060.95TNA702160.95TNA3020160.95373,450,417—40960.95386,156,7793.4501860.95385,676,1183.3602060.95385,394,6863.2702160.95385,258,7713.24520160.95375,427,413—40960.95386,156,7792.9501860.95385,676,1182.7602060.95385,394,6862.7702160.95385,258,7712.66020160.95375,804,390—40960.95386,156,7792.8501860.95385,676,1182.6602060.95385,394,6862.6702160.95385,258,7712.5T: termination because PHVF ran out of time while iterating; NA: not applicable;ρ: penalty factor; τ: percentage of variables fixed; α: cascading factor; θ0: initial slamming factor.Table 5
Effects ofτ and allocated time for the forest P34.
ForestρTime (min)τIterationsαθ0ValueImprovement (%)P341000015201100.951,333,129—401100.951,332,692−0.0506100.95NNA609100.95NNA7018100.951,347,1521.130201100.951,344,498—401100.951,344,260−0.0506100.95NNA609100.95NNA7018100.951,347,4910.245201100.951,345,620—401100.951,345,7610.0506100.95NNA609100.95NNA7018100.951,347,5550.160201100.951,346,799—401100.951,346,489−0.0506100.951,343,802−0.2609100.95NNA7018100.951,347,7590.1N: terminated because could not find feasible solution; NA: not applicable;ρ: penalty factor; τ: percentage of variables fixed; α: cascading factor; θ0: initial slamming factor.Table 6
Effects ofτ and allocated time for the forest P36.
ForestρTime (min)τIterationsαθ0ValueImprovement (%)P36500015201100.94,468,454—4022100.94,513,7171.05041100.94,532,8191.46051100.94,522,0131.27060100.94,525,6691.330201100.94,570,040—4022100.94,534,327−0.85041100.94,534,156−0.86051100.94,523,000−1.07060100.94,530,669−0.945201100.94,616,899—4022100.94,534,470−1.85041100.94,534,161−1.86051100.94,532,251−1.87060100.94,530,669−1.960201100.94,662,056–4022100.94,534,574−2.75041100.94,534,219−2.76051100.94,532,251−2.87060100.94,530,669−2.8ρ: penalty factor; τ: percentage of variables fixed; α: cascading factor; θ0: initial slamming factor.Table 7
Effects ofτ and allocated time for the forest P75.
ForestρTime (min)τIterationsαθ0ValueImprovement (%)P75200015201100.958,202,734—4019100.959,007,4031.45055100.9TNA6057100.9TNA7070100.9TNA30201100.958,272,824—4019100.959,106,2771.45055100.959,387,8841.96057100.959,324,5441.87070100.959,337,7821.845201100.958,297,317—4019100.959,131,6011.45055100.959,394,9421.96057100.959,336,5921.87070100.959,354,2491.860201100.958,310,384—4019100.959,248,4441.65055100.959,395,8161.96057100.959,339,1451.87070100.959,354,5251.8T: termination because run out of time iterating; NA: not applicable;ρ: penalty factor; τ: percentage of variables fixed; α: cascading factor; θ0: initial slamming factor.Table 8
Effects ofτ and allocated time for the forest P83.
ForestρTime (min)τIterationsαθ0ValueImprovement (%)P83200015201100.957,970,162—4014100.958,558,1701.05030100.958,516,9760.96032100.958,559,0691.07034100.958,592,5621.130201100.958,305,853—4014100.958,891,3811.05030100.958,833,3170.96032100.958,733,5630.77034100.958,661,3470.645201100.959,052,205—4014100.958,898,935−0.35030100.958,834,898−0.46032100.958,733,746−0.57034100.958,661,964−0.760201100.959,052,215—4014100.958,904,481−0.35030100.958,835,134−0.46032100.958,733,734−0.57034100.958,661,851−0.7ρ: penalty factor; τ: percentage of variables fixed; α: cascading factor; θ0: initial slamming factor.Table 9
Effects ofτ and allocated time for the forest P100.
ForestρTime (min)τIterationsαθ0ValueImprovement (%)P1005000015201100.9542,331,212—4027100.95TNA5063100.95TNA6069100.95TNA7079100.95TNA30201100.9542,494,130—4027100.9543,307,7441.95063100.9543,298,4391.96069100.95NNA7079100.9543,288,8411.945201100.9542,819,969—4027100.9543,309,4551.15063100.9543,299,1891.16069100.95NNA7079100.9543,290,1271.160201100.9542,931,057—4027100.9543,309,4020.95063100.9543,299,7590.96069100.95NNA7079100.9543,290,1270.8T: termination because run out of time iterating; N: terminated because could not find feasible solution; NA: not applicable,ρ: penalty factor; τ: percentage of variables fixed; α: cascading factor; θ0: initial slamming factor.In summary, as expected, everything else being equal, longer runtimes allow a higher objective function. As we can see forτ>40% which corresponds to fixing entirely the root-node and fixing some variables from the second stage, we have numerical issues. The numerical issues have two sources. On one hand, the algorithm may run out of time while fixing the variables or finished fixing the variables but does not have enough time to find a feasible solution to the REF. On the other hand, the basis of the REF may be disturbed in a way that the number of constraints and variables is not balanced so it is hard for branch and cut algorithm used by CPLEX to find an integral solution for the REF in the allocated time. This is the case for forest P34 with no solution when τ=60% even when the time is 1 hour. However, for the same forest, there is an integral solution when τ is raised to 70%. This behavior occurred mainly for forest P34 and P100 and is not observed when τ=40%.The second aspect is the impact ofτ on the optimality. The advantage that higher values of τ have over the small ones tends to vanish when the runtime is long. Similarly, the improvement in optimality from τ=40% to higher values of τ is low. In fact, for some forest such as P83, over τ=40%, the objective function value diminishes as τ increases which can be interpreted as during PHVF, many variables are fixed to values that are suboptimal. It clearly appears that fixing just the root-node is the best choice for these problems. It allows sufficient time for REF to find an integral solution while avoiding to impact the optimality, and limiting the disturbance on the structure of the basis matrix.


### 4.4.2. PHVF versus EF Solved Directly

Table10 presents the results from PHVF using τ=40% which is equivalent to fixing the root-node and then solving in parallel the SEF. As we can see for the small instances (64 scenarios), solving directly the extensive form outperformed using PHVF. This means the overhead of decomposing the problem into scenarios and solving each one overweighs the benefit. For the medium instances, the results are mixed. Out of the six forests, PHVF outperformed the EF in two cases. Even when PHVF under performed, it had results that were less than 1% away from the optimal solution. The benefit of using PHVF over the state of art commercial solver such as CPLEX is highlighted when dealing with big instances (512 scenarios). PHVF outperformed the EF for the six forests. In the case of P1 which is the biggest model with 1,363 units (Table 3), EF completely failed to solve the model because it could not fit it in the memory. For other cases for the same runtime, PHVF reached gains ranging from 0.84% to over 767% corresponding to forests P34 and P100, respectively. Even after leaving EF for 24 h, PHVF run for 15 min still outperformed in some cases. Comparing the gain to the gap from EF suggests that PHVF almost reached the optimal solution or its solution has a relative optimal gap less than 1%.Table 10
Comparison of solutions from PHVF and the extensive form.
Same duration with PHVF24 hΩTime (s)ForestPHVFEFGap (%)Gain (%)EFGap (%)Gain (%)64120P12.99E + 083.08E + 080.14−2.69———64120P341.04E + 061.05E + 060.14−1.64———64120P363.48E + 063.73E + 060.32−6.77———64120P754.71E + 074.75E + 071.59−0.81———64120P834.70E + 074.72E + 070.43−0.46———64120P1006.53E + 076.54E + 070.51−0.05———256600P13.86E + 083.87E + 080.05−0.173.87E + 080.02−0.28256600P341.35E + 061.33E + 062.381.561.35E + 060.73−0.05256600P364.73E + 064.64E + 062.051.754.70E + 060.820.52256600P755.94E + 075.94E + 070.69−0.065.97E + 070.19−0.55256600P835.90E + 075.93E + 070.28−0.475.93E + 070.22−0.50256600P1004.34E + 074.34E + 070.34−0.064.35E + 070.14−0.22512900P12.91E + 08NANA∞NANA∞512900P342.07E + 062.05E + 061.20.842.07E + 060.27-0.08512900P366.85E + 066.78E + 061.8761.136.79E + 061.640.89512900P758.48E + 074.50E + 0789.5788.638.50E + 070.22−0.28512900P838.47E + 074.49E + 0789.4488.768.45E + 070.630.27512900P1004.74E + 075.49E + 06767.9763.35.49E + 06767.9763.3NA: no results because of memory limitation.


## 4.4.1. Impact ofτ on the Optimality and Solution Time

Tables4–9 report the effect of τ on the value of the objective function with respect to the runtime for forests P1, P34, P36, P75, P83, and P100, respectively. We report as well, improvement% which is the percentage of increase (if positive) or decrease (if negative) of the objective function value from the objective function value for τ=20% for the same runtime.Table 4
Effects ofτ and allocated time for the forest P1.
ForestρTime (min)τIterationsαθ0ValueImprovement (%)P110000001520160.95371,524,817—40960.95386,068,4333.9501860.95TNA602060.95TNA702160.95TNA3020160.95373,450,417—40960.95386,156,7793.4501860.95385,676,1183.3602060.95385,394,6863.2702160.95385,258,7713.24520160.95375,427,413—40960.95386,156,7792.9501860.95385,676,1182.7602060.95385,394,6862.7702160.95385,258,7712.66020160.95375,804,390—40960.95386,156,7792.8501860.95385,676,1182.6602060.95385,394,6862.6702160.95385,258,7712.5T: termination because PHVF ran out of time while iterating; NA: not applicable;ρ: penalty factor; τ: percentage of variables fixed; α: cascading factor; θ0: initial slamming factor.Table 5
Effects ofτ and allocated time for the forest P34.
ForestρTime (min)τIterationsαθ0ValueImprovement (%)P341000015201100.951,333,129—401100.951,332,692−0.0506100.95NNA609100.95NNA7018100.951,347,1521.130201100.951,344,498—401100.951,344,260−0.0506100.95NNA609100.95NNA7018100.951,347,4910.245201100.951,345,620—401100.951,345,7610.0506100.95NNA609100.95NNA7018100.951,347,5550.160201100.951,346,799—401100.951,346,489−0.0506100.951,343,802−0.2609100.95NNA7018100.951,347,7590.1N: terminated because could not find feasible solution; NA: not applicable;ρ: penalty factor; τ: percentage of variables fixed; α: cascading factor; θ0: initial slamming factor.Table 6
Effects ofτ and allocated time for the forest P36.
ForestρTime (min)τIterationsαθ0ValueImprovement (%)P36500015201100.94,468,454—4022100.94,513,7171.05041100.94,532,8191.46051100.94,522,0131.27060100.94,525,6691.330201100.94,570,040—4022100.94,534,327−0.85041100.94,534,156−0.86051100.94,523,000−1.07060100.94,530,669−0.945201100.94,616,899—4022100.94,534,470−1.85041100.94,534,161−1.86051100.94,532,251−1.87060100.94,530,669−1.960201100.94,662,056–4022100.94,534,574−2.75041100.94,534,219−2.76051100.94,532,251−2.87060100.94,530,669−2.8ρ: penalty factor; τ: percentage of variables fixed; α: cascading factor; θ0: initial slamming factor.Table 7
Effects ofτ and allocated time for the forest P75.
ForestρTime (min)τIterationsαθ0ValueImprovement (%)P75200015201100.958,202,734—4019100.959,007,4031.45055100.9TNA6057100.9TNA7070100.9TNA30201100.958,272,824—4019100.959,106,2771.45055100.959,387,8841.96057100.959,324,5441.87070100.959,337,7821.845201100.958,297,317—4019100.959,131,6011.45055100.959,394,9421.96057100.959,336,5921.87070100.959,354,2491.860201100.958,310,384—4019100.959,248,4441.65055100.959,395,8161.96057100.959,339,1451.87070100.959,354,5251.8T: termination because run out of time iterating; NA: not applicable;ρ: penalty factor; τ: percentage of variables fixed; α: cascading factor; θ0: initial slamming factor.Table 8
Effects ofτ and allocated time for the forest P83.
ForestρTime (min)τIterationsαθ0ValueImprovement (%)P83200015201100.957,970,162—4014100.958,558,1701.05030100.958,516,9760.96032100.958,559,0691.07034100.958,592,5621.130201100.958,305,853—4014100.958,891,3811.05030100.958,833,3170.96032100.958,733,5630.77034100.958,661,3470.645201100.959,052,205—4014100.958,898,935−0.35030100.958,834,898−0.46032100.958,733,746−0.57034100.958,661,964−0.760201100.959,052,215—4014100.958,904,481−0.35030100.958,835,134−0.46032100.958,733,734−0.57034100.958,661,851−0.7ρ: penalty factor; τ: percentage of variables fixed; α: cascading factor; θ0: initial slamming factor.Table 9
Effects ofτ and allocated time for the forest P100.
ForestρTime (min)τIterationsαθ0ValueImprovement (%)P1005000015201100.9542,331,212—4027100.95TNA5063100.95TNA6069100.95TNA7079100.95TNA30201100.9542,494,130—4027100.9543,307,7441.95063100.9543,298,4391.96069100.95NNA7079100.9543,288,8411.945201100.9542,819,969—4027100.9543,309,4551.15063100.9543,299,1891.16069100.95NNA7079100.9543,290,1271.160201100.9542,931,057—4027100.9543,309,4020.95063100.9543,299,7590.96069100.95NNA7079100.9543,290,1270.8T: termination because run out of time iterating; N: terminated because could not find feasible solution; NA: not applicable,ρ: penalty factor; τ: percentage of variables fixed; α: cascading factor; θ0: initial slamming factor.In summary, as expected, everything else being equal, longer runtimes allow a higher objective function. As we can see forτ>40% which corresponds to fixing entirely the root-node and fixing some variables from the second stage, we have numerical issues. The numerical issues have two sources. On one hand, the algorithm may run out of time while fixing the variables or finished fixing the variables but does not have enough time to find a feasible solution to the REF. On the other hand, the basis of the REF may be disturbed in a way that the number of constraints and variables is not balanced so it is hard for branch and cut algorithm used by CPLEX to find an integral solution for the REF in the allocated time. This is the case for forest P34 with no solution when τ=60% even when the time is 1 hour. However, for the same forest, there is an integral solution when τ is raised to 70%. This behavior occurred mainly for forest P34 and P100 and is not observed when τ=40%.The second aspect is the impact ofτ on the optimality. The advantage that higher values of τ have over the small ones tends to vanish when the runtime is long. Similarly, the improvement in optimality from τ=40% to higher values of τ is low. In fact, for some forest such as P83, over τ=40%, the objective function value diminishes as τ increases which can be interpreted as during PHVF, many variables are fixed to values that are suboptimal. It clearly appears that fixing just the root-node is the best choice for these problems. It allows sufficient time for REF to find an integral solution while avoiding to impact the optimality, and limiting the disturbance on the structure of the basis matrix.


## 4.4.2. PHVF versus EF Solved Directly

Table10 presents the results from PHVF using τ=40% which is equivalent to fixing the root-node and then solving in parallel the SEF. As we can see for the small instances (64 scenarios), solving directly the extensive form outperformed using PHVF. This means the overhead of decomposing the problem into scenarios and solving each one overweighs the benefit. For the medium instances, the results are mixed. Out of the six forests, PHVF outperformed the EF in two cases. Even when PHVF under performed, it had results that were less than 1% away from the optimal solution. The benefit of using PHVF over the state of art commercial solver such as CPLEX is highlighted when dealing with big instances (512 scenarios). PHVF outperformed the EF for the six forests. In the case of P1 which is the biggest model with 1,363 units (Table 3), EF completely failed to solve the model because it could not fit it in the memory. For other cases for the same runtime, PHVF reached gains ranging from 0.84% to over 767% corresponding to forests P34 and P100, respectively. Even after leaving EF for 24 h, PHVF run for 15 min still outperformed in some cases. Comparing the gain to the gap from EF suggests that PHVF almost reached the optimal solution or its solution has a relative optimal gap less than 1%.Table 10
Comparison of solutions from PHVF and the extensive form.
Same duration with PHVF24 hΩTime (s)ForestPHVFEFGap (%)Gain (%)EFGap (%)Gain (%)64120P12.99E + 083.08E + 080.14−2.69———64120P341.04E + 061.05E + 060.14−1.64———64120P363.48E + 063.73E + 060.32−6.77———64120P754.71E + 074.75E + 071.59−0.81———64120P834.70E + 074.72E + 070.43−0.46———64120P1006.53E + 076.54E + 070.51−0.05———256600P13.86E + 083.87E + 080.05−0.173.87E + 080.02−0.28256600P341.35E + 061.33E + 062.381.561.35E + 060.73−0.05256600P364.73E + 064.64E + 062.051.754.70E + 060.820.52256600P755.94E + 075.94E + 070.69−0.065.97E + 070.19−0.55256600P835.90E + 075.93E + 070.28−0.475.93E + 070.22−0.50256600P1004.34E + 074.34E + 070.34−0.064.35E + 070.14−0.22512900P12.91E + 08NANA∞NANA∞512900P342.07E + 062.05E + 061.20.842.07E + 060.27-0.08512900P366.85E + 066.78E + 061.8761.136.79E + 061.640.89512900P758.48E + 074.50E + 0789.5788.638.50E + 070.22−0.28512900P838.47E + 074.49E + 0789.4488.768.45E + 070.630.27512900P1004.74E + 075.49E + 06767.9763.35.49E + 06767.9763.3NA: no results because of memory limitation.


## 5. Conclusions

In this paper, we have developed a method for solving multistage stochastic mixed integer programs that arise in natural resource management such as forest harvest planning. Although tested for growth uncertainty due to climate change, the method is valid as well for other sources of uncertainty such as errors in the model predicting forest growth, and uncertainty in the price of wood. This algorithm is also applicable to other disciplines where the decision variables are binary.In the forest industry, managers solve multiple instances of the same forest model. It is therefore necessary to have an algorithm that allows one to quickly solve the SMIP. Furthermore, having an algorithm that is faster allows the practitioners to explore different management options. However, problems tested in this experiment are quite small compared to the models in the industrial standard. Nevertheless, this requirement may be compensated by the availability of more powerful computation resources in the industry.The PHVF presented here overcomes some limitations listed in the literature regarding the choice of the penalty termρ. It is documented that large values of ρ lead to the phenomenon of cycles [9]. However, because PHVF fixes variables, such a problem is avoided. Notwithstanding all these benefits, PHVF has more parameters that need to be set compared to the classic PH. Fortunately, most of those parameters are easy to set and the limitation remains on the choice of the penalty term ρ. In the preliminary explorations, too low values of ρ led to some numerical issues while too high ones, although accelerated the convergence, negatively impacted the objective function value.The future extension of this work is to investigate the impact of different parameters on the solution quality. Similarly, having an algorithm that can dynamically determine the value of the penalty term would be an improvement for practitioners. On practical aspects, the question of the optimal number of scenarios necessary to approximate the growth change under climate change remains open. This can be done by computing the value of the stochastic solution, the expected value of perfect information for different scenario trees [3, 13, 20, 32].



---

*Source: 8965679-2020-12-12.xml*