1. Introduction
Forests sequester carbon from the atmosphere as part of photosynthesis, and carbon is stored within several ecosystem pools: aboveground, belowground, deadwood, litter, and soil. In the face of climate change that is triggered by increasing concentrations of atmospheric carbon, the role of forests in storing and emitting carbon is critical in the global carbon cycle. While changes in climate are projected to be more pronounced in higher latitudes [
1], the carbon cycling process in boreal forests is poorly understood, even though these ecosystems are estimated to contain more carbon than any other terrestrial ecosystem on earth [
2]. Forest biomass and soil carbon down to a depth of 1 m in boreal forests was estimated to be 559 gigatons by Watson et al. [
3]. In interior Alaska, an estimated 1.3 billion tons of carbon are sequestered within live trees alone [
4]. Soils are estimated to contain even more carbon, mostly in permafrost layers that have been frozen for many years and remained mainly undisturbed by human activities [
5]. However, forests in Alaska are subject to various disturbance processes, such as fire, insect activity, and permafrost degradation, that can significantly affect carbon storage and emissions within this boreal region. Climate forecasts predict that air temperatures will rise more rapidly at higher latitudes [
6] and rising annual temperatures across Alaska during the past 60 years confirm that this warming is already happening.
Due to the complex nature of the interactions between a changing climate, carbon cycle, fire regimes, and human disturbance, it is difficult to predict how carbon stocks in Alaskan boreal forests will change in the future. Monitoring these forests for changes in carbon stocks is crucial in understanding how, when, and at what pace climate and human changes will affect this ecosystem. In turn, these ecosystem changes will likely affect rural Alaskans and have a feedback effect on the velocity of climate change itself.
The Forest Inventory and Analysis (FIA) program, which is the national forest inventory of the United States conducted by the USDA Forest Service, aims to monitor the status and trends of all public and private forests in the nation via annual assessments [
7]. The program collects data to estimate the weight of forest carbon pools, among other forest attributes, according to a nationally standardized field sampling design and nationally consistent plot and subplot sizes [
7]. However, as of 2020, Alaska is the only state in the United States where the annual FIA system has not been fully implemented. Although it is estimated that Alaska contains 17% of all forest area in the United States [
8], most of its forest lands have not been inventoried in a systematic manner that is consistent with the national FIA program [
4,
9]. At the same time, the lack of infrastructure in this region makes the application of the standardized national field sampling design cost-prohibitive. Remote sensing techniques, such as airborne lidar (light detection and ranging) and digital imaging, may allow for more cost-effective estimation of carbon across large areas. However, lidar techniques can only be used to accurately estimate some, but not all, of the carbon pools that comprise total forest carbon. For example, aboveground tree carbon can be estimated with lidar quite well (the correlation coefficient for tree carbon based on observed pilot data in our study was 0.91), but field surveys are more appropriate for soil carbon pools. In other words, estimation of total forest carbon across the multiple pools must entail a combination of lidar methods, such as lidar, aerial or satellite imaging, and ground sampling techniques. Because different sampling designs have varying projected standard errors and costs, finding the optimal combination of ground vs. lidar protocols subject to budgetary, logistical, and quality constraints is a non-trivial task. In many remote regions, logistical issues might be the binding constraint on sampling. Internationally, developing nations with significant forest inventories but limited budgets and infrastructures wishing to adopt REDD+ (reducing emissions from deforestation and forest degradation, plus the sustainable management of forests and the conservation and enhancement of forest carbon stocks) regimes in exchange for financial compensation will also have to face the same constraints. The REDD+ program requires national sampling, reporting and certification protocols for forest carbon stocks and emissions that are transparent and verifiable [
10]. While the requirements of REDD+ programs are often a primary driver for the implementation of forest carbon monitoring systems, there are other national-scale reporting activities that depend upon forest monitoring, including the FAO global forest resources assessment [
11] and nationally determined contributions (NDC) submissions to reduce national emissions and adapt to the impacts of climate change as required by the Paris Agreement [
12].
The application of mathematical programming can both (1) inform the design of new carbon sampling strategies in remote areas with no monitoring systems (i.e., REDD+), and (2) provide a quantitative basis for evaluating tradeoffs between sampling alternatives in existing inventory programs (i.e., the FIA program in the United States). However, while previous studies have evaluated sampling alternatives for forest inventory through a straightforward comparison of standard errors for various combinations of lidar and field plots [
13], to date very little work has been done to document the tradeoffs behind the survey design choices and the associated costs using more advanced mathematical optimization techniques. One notable exception was the work of Köhl et al. [
14], who compared the projected percent standard errors of five sampling design alternatives in two different cost scenarios for remotely sensed data acquisitions in Puerto Rico. While the authors’ tradeoff analysis did provide future users with a rough strategic guide based on the expected precision of the five design alternatives, it failed to take actual logistical and budgetary constraints into account. For example, available budgets and skilled field labor might change, and some field plots might simply be inaccessible, even in developed countries such as the United States (Alaska). These problems could be even more pronounced in developing nations aspiring to be part of REDD+.
In this study, we present a multi-objective mathematical programming model that can help analysts identify the set of so-called Pareto optimal combinations of field and lidar sampling strategies for boreal forest carbon inventory. We seek those alternative combinations of field and lidar plot intensities that would lead to the smallest possible expected standard error in the estimated mean carbon tonnage for the different pools. A field-lidar plot intensity combination is “Pareto optimal” (i.e., efficient or non-dominated) if no other intensity combination would yield a smaller expected standard error for one carbon pool without increasing the standard errors for any one of the other pools. Multi-objective mathematical programs are sets of functions that represent the objectives of the analysis, such as minimizing standard errors, and inequalities that represent the constraints, such as budgetary restrictions, within which the objectives can be pursued. Since lidar and field sampling methods have different associated costs and are projected to lead to different standard errors for the different carbon pools, our objective is to minimize standard errors. The constraints are finite budgets of both money and time to be spent on ground plots, as well as time to be spent on processing remotely sensed data. Using data collected from an FIA pilot project in the Tanana Valley of interior Alaska (USA), we demonstrate how multi-objective mathematical programming can help identify the best compromise sampling intensity combinations.
2. Materials and Methods
For simplicity of illustration, we combined the five carbon pools, mentioned previously, into three that can be measured directly: tree carbon, down woody material (DWM) carbon, and soil carbon. The first pool, “tree”, included both live and dead standing trees; DWM included both fine (<3 inch (7.62 cm) diameter) and coarse (≥3 inch (7.62 cm) diameter) down woody materials; and “soil” included litter. We chose not to include belowground biomass in our analysis because it is typically predicted based on aboveground tree biomass. Our case study site for the estimation of the three carbon pools listed above was a 13.7 million ha forest area in the Tanana valley in Alaska (
Figure 1). Forty percent of this area has been previously inventoried [
4]. The region is one of the most sparsely populated in the United States and has very few roads and only one small city. As such, measuring nominal 5642 ground plots per standard FIA program specifications is not economically viable (standard FIA plot intensity is 1 plot per 6000 acres (2428 ha)). Thus, there is interest in using airborne lidar technologies to monitor the area in conjunction with a limited number of field plots.
The rest of this Materials and Methods section is structured as follows. We start by providing a mathematical description of the proposed multi-objective linear programming model. We then parameterize the model for our case study using parameter estimates derived exclusively from observed data (provided by a pilot study conducted within the area of interest). The primary objective of this exercise is to give the reader a sense of how this model works and what benefits it can offer. Third, we use a computer algorithm, implemented in R, to iterate through all possible solutions to the problem and identify those that are Pareto optimal. Last, we perform a sensitivity analysis to determine how changes in model parameter values can affect the solution set. The goal of the sensitivity analyses is to demonstrate the utility of the proposed model in situations where parameter estimates based on observed data are not available.
2.1. The Multi-Objective Linear Programming Model
To formulate this problem as a linear program, we first enumerated the feasible combinations of field and lidar plot intensities, where intensity refers to the number of plots of each type to sample in the region. We use the term “lidar plot” to describe a location where only remotely sensed metrics (i.e., no field measurements) are available. In our model, we use a combination of lidar and field plots in a so-called “double sampling” inventory design, where remotely sensed metrics (lidar-based canopy heights, etc.) are collected on a large sample of relatively inexpensive “lidar plots” within a region of interest, and then expensive field measurements (i.e., direct measurements of forest carbon) are collected on a smaller subsample of these plots (“field plots”) [
15]. For field plots, the maximum intensity was set to be approximately the maximum number of field plots that could fit within the given budget: 250 plots. For lidar plots, the maximum intensity was set to 5480, as requested by the decision makers. Next, the minimum intensity was set at a small number relative to the maximum intensity; we chose 10 for field plots and 180 for lidar plots. Since the objective is to minimize standard error, the model formulation will seek to maximize the number of plots subject to the constraints. Thus, the solutions the model formulation returns should not be sensitive to the chosen minimum intensity values. The decision variables, denoted by
, are therefore defined as binary, with one decision variable for every combination of field plot and lidar plot intensity that we considered, as represented by subscripts
i and
j, respectively. If the intensities evaluated included every integer value between the minimum and maximum intensity values, there would be more than one million decision variables. Since our goal was to illustrate the model, we chose to only examine intensities at certain intervals (10 for field plots, 100 for lidar plots) in order to reduce the number of decision variables to a more feasible number. In other words, we considered 10, 20, 30, 40 … and up to 220 field plots and 180, 280, 380 … and up to 5480 lidar plots. This resulted in 25 choices for field plots ((250 − 10)/10 + 1 = 25) and 54 for lidar plots ((5480 − 180)/100 + 1 = 54)—a total of 1350 field vs. remote sampling combinations (25 × 54 = 1350).
The linear programming formulation is as follows:
subject to:
In the above series of equations, i and j denote the discrete intensity options for the number of field and lidar plots to sample, respectively. Objective function coefficients and are the calculated standard errors for tree carbon, down woody material carbon, and soil carbon estimates (in tC/ha) at field plot intensity i and lidar intensity j. Constraint coefficient is the cost of measuring i field plots and j lidar plots (in USD), and B is the total budget in USD. Coefficient is the time needed to sample i field plots (in hours), and F is the total number of hours that can be spent in the field in two seasons, assuming three field crews who work 8 h per day. Coefficient is the time needed to collect data from j lidar plots (in hours), and L is the total number of hours allocated to this purpose during the two-year time period. Equation (1) is a composite objective function that seeks to minimize the standard error in estimates of tree carbon, down woody material carbon, and soil carbon, respectively. Inequality (2) is the budget constraint that forces the sum of the total field plot measurement cost for i field plots and the cost of lidar analyses for j lidar plots to be below budget B. Inequality (3) ensures that the field plot measurement time fits within the time allocated for fieldwork during the two-year planning horizon (F). Inequality (4) ensures that the time needed to collect remotely sensed data does not exceed what is available in a two-year time window (L). Inequality (5) allows the model to choose only one combination of field and lidar plot densities. Equation (6) defines the decision variables in the formulation as binary.
2.2. Values of Parameters and Constants
The standard error coefficients in objective function (1) were calculated as a function of field plot intensity
i and lidar plot intensity
j, based on the following formula developed for double sampling (Cochran, 1977) [
15]:
Parameter
is carbon variance for pool
k, where
k =
t,
w, and
s for tree carbon, down woody material carbon, and soil carbon variances, respectively;
is the lidar correlation coefficient for pool
k, where
k =
t,
w, and
s for tree carbon, DWM carbon, and soil carbon correlations; and
Z is the maximum number of lidar plots that can fit into the 13.7 million hectare area (assumed to be one per ha). We calculated the carbon variances and lidar correlations based on a preliminary analysis of tree carbon, down woody material carbon, and soil carbon using a pilot field and lidar dataset in the Tanana valley of interior Alaska [
16]. Our estimates for tree carbon, down woody material carbon, and soil carbon variances were 845.8, 32.4, and 2575.6 Mg/ha, respectively.
The correlation between lidar metrics and field plot-based estimations for each carbon pool is also a component of the standard error,
, which denotes the standard error for carbon pool
k using field plot intensity
i and lidar plot intensity
j. In this study, we assumed that lidar technology would be used for remotely sensed sampling. The lidar-based metric used for the lidar correlations was mean plot-level tree segment height. An individual tree segmentation algorithm was run on a 1 m lidar canopy height model, and the maximum height within each segment (above 1.37 m or breast height) was extracted as the segment height [
17]. The plot-level lidar metric was then calculated as the mean lidar segment height for all segments within the 0.06744 ha footprint of the FIA plot (this corresponds to the combined area of all four subplots on an FIA plot). Although there are a large number of possible lidar metrics, we chose this metric, representing the mean tree height within the plot, since it is easy to interpret and has been shown to be highly correlated with aboveground biomass (AGB) in boreal forests [
18]. For the purpose of this study, tree carbon was based on all trees with a breast-height diameter of 2.54 cm. Likewise, although there are many other possible approaches to estimating soil carbon and DWM carbon with satellite and/or airborne lidar, restricting our analysis to a single lidar-based estimation provides a straightforward yet realistic basis for comparing the relative costs and trade-offs of lidar and field sampling in an optimization framework. Furthermore, our intention in this paper is to demonstrate the application of a multi-objective optimization approach that could be parameterized with different standard error values, correlations, cost levels, and logistical constraints, depending on the specific inventory context within which it is applied. Based on our pilot observation data, our estimates for tree carbon, down woody material carbon, and soil carbon correlations were 0.91, 0.23, and 0.05 Mg/ha, respectively.
The rest of the model parameters were provided by USFS Forest Inventory and Analysis staff in Alaska. The total budget, B, USD 2 M. The amount of time to measure i field plots, was calculated as h, and the amount of time to analyze j lidar plots, , was calculated as h. Costs were assumed to be USD 9285 for each field plot and USD 100 for each lidar plot. The total time for fieldwork, F, and lidar analysis, L, were calculated by assuming a two-year time period with 45 field days in each year and 50 weeks of analysis time per year, assuming 40 h workweeks as the standard.
2.3. Finding Pareto Optimal Sampling Intensity Combinations
We solved the multi-objective program (1)–(6) for Pareto optimality by first enumerating all feasible solutions with respect to the budget and time constraints. Then, we checked which subset of these solutions were non-dominated or Pareto optimal. This was done by iterating through each feasible solution and testing if any of the other feasible solutions had smaller standard errors with respect to one of the three carbon pools without having a larger standard error with respect to any of the other pools. We were able to do this due to the relatively small number of discrete plot intensity combinations that we decided to analyze for illustration purposes. If a greater number of plot intensity combinations are to be evaluated, specialized multi-objective solution algorithms could be used, such as those discussed in Tóth et al. [
19] or in Tóth and McDill [
20].
2.4. Sensitivity Analyses
As briefly discussed in
Section 2, estimates of variance for forest carbon pools are not always available in remote regions. In order to test our model’s robustness in the face of this uncertainty, we re-solved the model with perturbed variance estimates of 300, 400, 500, 600, 700, 800, 900, and 1000 Mg/ha for tree carbon (in addition to our best estimate of 845.8 Mg/ha), with 10, 20, 30, 40, 50, 100, 200, 300, 400, and 500 Mg/ha for down woody material carbon (in addition to our best estimate of 32.4 Mg/ha), and with 1500, 2000, 2500, 3000, 3500, and 4000 Mg/ha for soil carbon (in addition to our best estimate of 2575.6 Mg/ha). We also examined the effect of changing the carbon correlations between remote and field estimates of soil carbon to 0, 0.1, 0.2, 0.4, 0.6, and 0.9 (in addition to our best estimate of 0.05), to 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.88, and 0.9 for down woody material carbon (in addition to our best estimate of 0.23), and to 0.4, 0.5, 0.6, 0.7, 0.88, and 0.95 for tree carbon (in addition to our best estimate of 0.91). Last, we also re-solved the model for more optimistic budget scenarios of USD 2.25 M, USD 2.5 M, and USD 2.75 M (in addition to our baseline budget of USD 2.0 M).
3. Results
We found four Pareto optimal sampling intensity combinations (
Table 1). These were 210 field plots and 480 lidar plots, 200 field plots and 1380 lidar plots, 190 field plots and 2280 lidar plots, and 180 field plots and 3280 lidar plots. At the baseline budget of USD 2 M, having more than 210 or fewer than 160 field plots is not optimal in any of the variance and correlation scenarios that we considered, even though there were feasible solutions at greater and smaller sampling intensities.
We also found that only the budget constraint affected feasibility; the time constraints were not binding in any of the variance and correlation scenarios. When we increased the budget, the number of Pareto optimal sampling intensity combinations first increased from four at USD 2 M to five at USD 2.25 M; then it went back down to four at USD 2.5 M and two at USD 2.75 M (
Table 2).
As expected, the standard errors decreased with increasing budgets from a range of 1.01–1.47 tC/ha for tree carbon at USD 2 M to 0.9–0.99 tC/ha at USD 2.5 and 0.86 tC/ha at USD 2.75. For DWM carbon, the expected standard errors decreased from 0.39–0.41 tC/ha at USD 2 M to 0.35–0.37 tC/ha at USD 2.5 M and 0.35–0.36 tC/ha at USD 2.75 M. Last, for soil carbon, we found that the standard errors went down from 3.5–3.78 tC/ha at USD 2 M to 3.21–3.42 tC/ha at USD 2.5 M and 3.21–3.27 at USD 2.75 M. Please note that standard error ranges are provided because different Pareto optima at a given budget resulted in a range of different associated standard errors (
Table 1 and
Table 2).
Changing the variances of the three carbon pools did not affect which sampling schemes were Pareto optimal. Increasing the correlation between field and remotely sensed tree carbon estimates from 0.91 to 0.95, however, changed the set of Pareto optima from four to six field-remote sampling intensity combinations: (210, 480), (200, 1380), (190, 2280), (180, 3280), (170, 4180), and (160, 5080). Increasing the correlation between field and remotely sensed soil carbon estimates to 0.4 removed one efficient point (210 field plots and 480 lidar plots), while increasing all three correlations to 0.9 led to a unique efficient solution with 180 field plots and 3280 lidar plots. While we do not necessarily consider either of these scenarios to be realistic, the exercise does show that other Pareto optimal solutions can exist. Lastly, we note that the estimates of soil carbon always have the greatest associated standard errors, since correlation of soil carbon with lidar data is small.
4. Discussion and Conclusions
In this study, we illustrated how multi-objective mathematical programming can be used to optimize forest carbon sampling with a combination of field and lidar measurements in order to minimize the standard error of estimates of three specific carbon pools: tree carbon, down woody material carbon, and soil carbon pools. These results help illustrate the utility of lidar observations in supporting data collection in the face of limited time and budget constraints.
Using the Tanana Valley in Alaska, United States, as an illustrative case study, we identified four Pareto optimal combinations of field and lidar sampling intensities (
Figure 2).
We also found, via sensitivity analyses, that these four sampling intensity combinations were rather robust to changing parameter values. Since none of these four sampling strategies consist solely of using either the field or the lidar method, we conclude that combining the two is the optimal course of action. This result matches what Stephen et al. [
21] found for New Zealand forest inventories. Given the budget constraints, field sampling alone was not feasible. We also found that spending the entire budget on field sampling produces standard errors that can be reduced for all three carbon stocks simultaneously (unless the correlation between field and lidar observations is less than zero). In decision scientific terminology, this means that sampling strategies that comprise solely field measurements are dominated. The same is true for sampling schemes with very low field plot intensities, unless field and lidar observations are nearly perfectly correlated across carbon stocks. The location of the four Pareto optimal sampling strategies in 3-dimensional standard error space (
Figure 2) suggests that choosing (210, 480) for field and lidar plot sampling intensities would likely be advantageous for the decision maker, as this option would lead to minimum standard errors for both down woody material carbon and soil carbon. Nonetheless, it is noteworthy that there is a lot of standard error space left between the (200, 1380) and the (210, 480) sampling intensity combinations. Alternatively, if the decision maker was particularly interested in minimizing standard error in tree C, given that this carbon pool is most influenced by dominant disturbance processes in interior Alaska (fires and insects) and management activities (fire suppression, bioenergy production, and timber harvesting), this could motivate them to choose the alternative with smaller standard error in tree C estimation (200, 1380) while still maintaining reasonable standard error for the DWM and soil C pools. If the decision makers were interested in exploring further options between these alternatives, then our model could be run for all nine field sampling intensities (201, 202 … 209) to see if any one of them would lead to Pareto optimality. We also note that the standard error ranges for the three carbon pools we examined are relatively narrow across the four Pareto optimal sampling intensity combinations. In other words, the tradeoffs with regard to our ability to accurately estimate carbon in these three pools are not that great. This finding speaks to the power of multi-objective optimization. No matter which of the four strategies the analyst chooses, the expected standard errors will be reasonable regardless of what the underlying carbon variances or correlations with lidar are, within extremely broad ranges.
As to which of the Pareto optimal sampling intensity combinations should be chosen, this is ultimately a decision that the inventory program (e.g., the USFS FIA) and their data user groups will have to make. While the national FIA program has established precision guidelines for area and volume estimates in the lower 48 states (Bechtold and Patterson, 2005) [
7], these precision guidelines are not applicable to other inventory attributes (e.g., soil carbon) and in cases where the standard FIA sampling design is not logistically or economically feasible (e.g., interior Alaska). In addition, the Global Forest Observations Initiative has published methods and guidance to inform the development of monitoring, reporting, and verifying (MRV) systems in the context of REDD+ and national greenhouse gas reporting, including consideration of costs, sampling design, bias, and precision (GFOI, 2020) [
22]. While these precision guidelines are useful points of reference, in the context of a multi-resource inventory (such as the FIA program) the decision maker must weigh the value of precision for a number of inventory attributes against the costs in arriving at an acceptable solution, taking into account the needs of a wide variety of stakeholders. Our contribution with the present work is to provide a mathematical framework to quantify and visualize the Pareto optimal tradeoffs that come with any one of those decisions. Unlike other multi-objective sampling optimization models that were documented in the literature (e.g., Lark 2006) [
23], ours maps the tradeoffs associated with the expected standard errors across inventory attributes such as tree carbon, down woody material carbon, and soil carbon.
We acknowledge that when weighing which sampling scheme to employ, knowing which carbon stock is the most important in minimizing standard errors would be helpful, as standard errors do change under different assumptions of carbon variability and field-lidar correlation. To this end, we created an optimization framework that is flexible and can be updated as improved estimates of carbon stock variability and correlations between field and remotely sensed data become available. In addition, the proposed multi-objective mathematical programming method can optimize sampling for other metrics as well, such as invasive species cover or mean forest stand age. The model can also be used to choose optimal intensities for any two sampling strategies with known estimates of standard error rates and correlations. Thus, it would also be applicable to other remote sensing technologies that might become available in the future, as long as information about their error and correlation profiles is provided. More broadly, we suggest that the multi-objective linear programming framework we established here is even more generalizable to other spatial sampling applications where the logistical and budgetary constraints are entirely different. This is due to the general mathematical structure of linear programs, comprising a set of linear functions that can represent various sampling objectives accompanied by a set of linear inequalities that can represent various logistical and budgetary constraints. Another example of the flexibility of this mathematical structure is that the model can easily be modified to minimize the sum total of standard errors (Equation (8)) or the sum of weighted standard errors (Equation (9)), or even to minimize the maximum standard error over the three pools (also known as the minimax formulation; see Equation (10):
where
,
, and
denote the relative weights assigned to estimating tree carbon, DWM carbon, and soil carbon, and
Whether any of these make sense will depend on the objectives and the unique features of the survey at hand. Although standard errors do propagate in an additive fashion across the pools, and thus the estimates of total carbon can be optimized by minimizing the total of standard errors, inventory analysts might still want to optimize the individual standard errors associated with each pool, thereby developing a more holistic understanding of the tradeoffs in the totality of the expected standard error space.
Lastly, we would like to emphasize that multi-objective linear programming is not the only optimization method that can be used to aid forest carbon sampling design. We would argue, however, that it is likely the best choice today for the following reasons. Ad hoc heuristics, such as Monte Carlo simulation, simulated annealing, tabu search, and genetic algorithms, cannot guarantee optimality, which is a problem when one seeks to make the best possible use of limited budgets for a specific purpose. Unlike heuristics, LPs are exact mathematical models that can find proven optimal solutions given enough computing time. Time can be a barrier to finding true optima when some of the variables in an LP are restricted to be discrete (also known as integer programs, or IPs) and the size of the problem is large in terms of the number of variables. However, even under these circumstances, IPs provide an upper bound on the gap that remains between the true optimum and the best incumbent found in the allotted time. Thus, one would know how good the solution is in terms of optimality. Heuristics cannot provide such information. What they can do, though, is find good solutions quickly for some hard problems that might otherwise be intractable by IPs (some combinatorial problems fall in this category, where worst-case IP solve times are more than a polynomial function of problem size; these problems are called NP-Hard). However, the forest carbon inventory problem we examined in this paper is computationally trivial; we were able to find the true optima in a fraction of a second under all parameter value combinations. Apart from being able to guarantee optimality or an upper bound on the objective value of the incumbent, perhaps the biggest advantage of LPs over heuristics is that encoding and solving LP models do not require custom application development. There are many off-the-shelf solvers that are easy to access and free to use (e.g., Excel Solver). Model representation, parameterization, the addition of constraints and variables, and sensitivity analyses are trivial and streamlined for the average user. With heuristics, on the other hand, not only does the developer have to build the model using custom application, but the developer must also develop and tune the solution algorithm itself. While the general algorithmic design is widely available for most popular heuristics (e.g., simulated annealing, or SA, and evolutionary algorithms), coding and executing these algorithms is not trivial. Even if the developer succeeds in finely tuning the parameters (e.g., setting the temperature parameter in SA), any addition of new constraints or variables to the model will require going back to the drawing board to adjust the code and re-tune the parameters. This is a big issue for sampling optimization, because logistical and budgetary constraints must often be added, modified, or dropped on a case-by-case basis. Unlike the situation with LP solvers, this can only be done by the developer. Last, the computational performance of software and hardware resources have been increasing at a breakneck speed for the last 30 years, rapidly making ad hoc custom heuristics less and less advantageous when compared with exact mathematics in optimization.
In sum, the principal value of our work is a novel application of multi-objective mathematical programming that not only allows forest inventory analysts to determine what combination of field and lidar plot intensities should be used for sampling design to minimize expected standard errors in estimating boreal forest carbon pools, but also to streamline what prior data must be collected or estimated, and with what accuracy, in order to find optimal sampling designs. As an example, our demonstrative case study in the Tanana Valley, Alaska, United States revealed that the accuracy in estimating variances in the amount of carbon in the various pools via pilot surveys is much less important for finding optimal sampling designs than having more accurate estimates of the budget that is available for the surveys.