Harvest scheduling is a resource planning process aimed at deciding what areas of the forest (i.e.; stands) to harvest and when. On a small scale, Forest Managers may simply rely on experience, memory, and/or field notes when planning harvests. However, on a larger scale, this is inefficient and highly frowned upon.

At its core, the Forester’s harvest schedule is a plan about what resources to allocate and where. Linear Programming (LP) is a mathematical method used to “... allocate limited resources to competing activities in an optimal manner.” (J. Buongiorno & J.K. Gilless, 2003, p.9). Thus, LP models provide a good solution for the harvest scheduling problem.

An article in the Northern Journal of Applied Forestry, written by Macmillian and Fairweather (1988), provides an example of a short-term harvest schedule for an Industrial forest ownership in Pennsylvania using an LP model. In this article, I would like to expand on the example they provided, with a slightly modified version using a Python library called PuLP. For details about Macmillian and Fairweather’s model I refer you to their article, which can be accessed here.

Problem Description

The example provided by Macmillian and Fairweather was structured using 15 harvest units of northern hardwood stands. In their words, “The objective was to develop a harvest schedule which maximized the net present value (NPV) of timber harvested from stands over a 5-year period.” (Macmillian & Fairweather, p.145). They were also attempting to satisfy the annual timber demand for a sawmill.

The project setup described in their article involved calculating NPV for each stand over a 5-year period using local stumpage prices, an assumed 3% volume growth rate, and a 4% discount rate. This process resulted in 75 possible NPV combinations, one for each combination of stand and harvest year (15 stands x 5 years = 75). The final structure of their model included the 75 decision variables, 15 harvest area constraints, and 10 constraints satisfying the mill requirement for harvesting 3800-4200 MBF each year from company land. Their problem was solved using LINDO, a commercial LP software package.

For my example, I plan to adopt the same objective of maximizing NPV of timber harvested from stands over a 5-year period. My example evaluates 17 stands for a total of 85 stand-harvest year combinations (17 stands x 5 years = 85). The NPV calculations are based on current pine and hardwood inventory data and current stumpage prices from a privately managed forest located in the Pineywoods region of Texas. However, I am using the same growth and discount assumptions that Macmillian and Fairweather used, 3% and 4% respectively.

Finally, I chose to setup my problem as an Integer Programming problem, meaning the decision variables (acres harvested from a stand in a given year) are discrete (a.k.a integers). This just means that the model will not consider a fractional acre (e.g.; 99 instead of 99.12). This is a personal preference, because in operational planning I don’t like dealing with fractional acres.

Python PuLP Package

PuLP is a modeling framework for solving Linear and Integer Programming problems. PuLP is capable of interfacing with several known and recognized solvers including CPLEX, COIN and Gurobi. The strength of using a package like PuLP is that, once the problem has been defined, the code is very minimal (assuming you understand Python data structures and syntax).

PuLP can be installed with Pip or Conda, or similar Python package installers. For additional information about PuLP, please check out the documentation here.

Basic Harvest Schedule using PuLP

As indicated above, once the problem is formulated, the coding is very minimal. In this case I was able to run the model in as few as 19 lines of code (not counting formatting and comments – which make the code easier for us humans to read).

Here’s the code line by line.

# module imports
from pulp import LpMaximize, LpProblem, LpStatus, lpSum, LpVariable

Next, I initialize the model and the data structures needed for the problem. Note that the “npv” variable below is a list that contains each stand-year combination of the calculated NPV. The first element of the list is the NPV value of harvesting an acre from Stand A in Year 1 (A1), and the last element is the NPV value for harvesting Stand Q in Year 5 (Q5). These naming conventions could be anything you choose, but I’ve found that using alpha characters for stands and numeric for years, makes the most sense to me.

# initialize the model
model = LpProblem(name='harvest_schedule', sense=LpMaximize)

# initialize the model variables
stands = ['A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'J', 'K', 
     'L', 'M', 'N', 'O', 'P', 'Q']
years = [1, 2, 3, 4, 5]
var_list = [s + str(y) for s in stands for y in years]
acres = [107, 113, 9, 105, 43, 69, 102, 66, 227, 96, 288,
     20, 119, 160, 15, 57, 9]
npv = [
    3523, 3258, 2896, 2476, 2035,
    2824, 2611, 2321, 1984, 1631,
    2824, 2611, 2321, 1984, 1631,
    2721, 2516, 2236, 1912, 1571,
    2824, 2611, 2321, 1984, 1631,
    2126, 1965, 1747, 1494, 1228,
    3246, 3001, 2668, 2281, 1875,
    2824, 2611, 2321, 1984, 1631,
    3581, 3310, 2943, 2515, 2068,
    2796, 2585, 2298, 1964, 1615,
    2972, 2748, 2443, 2088, 1716,
    810, 749, 665, 569, 468,
    3434, 3175, 2823, 2413, 1983,
    3281, 3033, 2697, 2305, 1895,
    2824, 2611, 2321, 1984, 1631,
    2824, 2611, 2321, 1984, 1631,
    2824, 2611, 2321, 1984, 1631
]

npv_by_harvest = dict(zip(var_list, npv))
harvest_acres = dict(zip(stands, acres))

Next, I create the decision variables. These will be integer variables with a lower bound of 0, satisfying the model requirements that harvested acreage cannot be negative, and must be a whole number.

# create decision variables
x = LpVariable.dicts('hvst', [x for x in npv_by_harvest.keys()], 
     lowBound=0, cat='Integer')

Then I add Stand Constraints. These could be anything really, but below are the stand constraints that I chose to use for this example.

Stand Constraints

1. stand harvests cannot exceed stand acres
2. a harvest must be greater than the smallest stand size to be operable
3. annual harvest cannot exceed 425 acres

# add stand constraints: must be less than total acres in a stand
# but greater than the minimum stand acreage
for s in stands:
    model += lpSum([x[std] for std in var_list if std.startswith(s)]) 
      == harvest_acres[s]

min_acres = min(acres)
for s in stands:
    model += lpSum([x[std] for std in var_list if std.startswith(s)])
      >= min_acres

# add annual harvest constraint: cannot exceed 425 acres per 
# harvest year
for y in years:
    model += lpSum([x[yr] for yr in var_list if yr.endswith(str(y))])
      <= 425

Now, I add the Objective Function. This is the function that will seek to maximize NPV for our model. This is maximizing the sum of the acres harvested in each stand by the NPV for a particular stand-year combination. In other words it does something like this:

MAXIMIZE SUM(A1 * NPV of A1 + A2 * NPV of A2 + … Q5 * NPV of Q5)

# objective function
model += lpSum([x[z] * npv_by_harvest[z] for z in var_list])

The final steps are to solve the model and output the results to standard output.

# solve the model
status = model.solve()

# print the results
print(f'status: {model.status}, {LpStatus[model.status]}\n')
print(f'objective: {model.objective.value()}\n')
print('variables:')
for i in model.variables():    
    print(f'\t{i.name}: {i.value()}') 
print('\nconstraints:')
for name, constraint in model.constraints.items():
    print(f'\t{name}: {constraint.value()}')

Instead, it’s possible to output the results and model to a Text file as shown in the alternative version below.

# ouptut the results
with open('output.txt', 'w') as writer:
    writer.write(f'status: {model.status}, {LpStatus[model.status]}\n')
    writer.write('\n')
    writer.write(f'objective: {model.objective.value()}\n')
    writer.write('\n')
    writer.write('variables:\n')
    for i in model.variables():
        writer.write(f'    {i.name}: {i.value()}\n')
    writer.write('\n')
    writer.write('constraints:\n')
    for name, constraint in model.constraints.items():
        writer.write(f'    {name}: {constraint.value()}\n')
    writer.write('\n' * 3)
    writer.write(str(model))

Summary

The resulting output, as printed, can take some time getting used to. With the solution stored in the status variable, users have the option to output the information as they see fit. Following is a sample of my output of the model results. Note that I am only showing the non-zero decision variables here for the sake of brevity.

These results show that the model found an optimal solution, and maximized the objective function resulting in $4,341,825 of Net Present Value. The variables indicate a harvest of 107 acres in stand A in year 1, 113 acres of stand B in year 3, and so on.

As you can see, Linear Programming using the PuLP package is a great tool for resource scheduling problems, including Harvest Scheduling. The PuLP software package is free and available to anyone, but requires users to have some basic understanding of the Python programming language to interact with a model.

References:

Buongiorno, J., Gilless, J.K. (2003). Decision methods for Forest Resource Management. Academic Press.

Macmillian, D. C., Fairweather, S.E. (1998). An application of linear programming for short-term harvest scheduling. Northern Journal of Applied Forestry, 5(2), 148-152.

Categories: Data Analytics Forest Inventory Python

Tagged as: Data Analysis Harvest Schedule Linear Programming Python

JwL