Summary
The prioritizr R package uses mixed integer linear programming (MILP) techniques to provide a flexible interface for building and solving conservation planning problems (Rodrigues et al. 2000; Billionnet 2013). It supports a broad range of objectives, constraints, and penalties that can be used to custom-tailor conservation planning problems to the specific needs of a conservation planning exercise. Once built, conservation planning problems can be solved using a variety of commercial and open-source exact algorithm solvers. In contrast to the algorithms conventionally used to solve conservation problems, such as heuristics or simulated annealing (Ball et al. 2009), the exact algorithms used here are guaranteed to find optimal solutions. Furthermore, conservation problems can be constructed to optimize the spatial allocation of different management zones (or actions), meaning that conservation practitioners can identify solutions that benefit multiple stakeholders. Finally, this package has the functionality to read input data formatted for the Marxan conservation planning program (Ball et al. 2009), and find much cheaper solutions in a much shorter period of time than Marxan (Beyer et al. 2016).
Introduction
Systematic conservation planning is a rigorous, repeatable, and structured approach to designing new protected areas that efficiently meet conservation objectives (Margules & Pressey 2000). Historically, conservation decision-making has often evaluated parcels opportunistically as they became available for purchase, donation, or under threat. Although purchasing such areas may improve the status quo, such decisions may not substantially enhance the long-term persistence of target species or communities. Faced with this realization, conservation planners began using decision support tools to help simulate alternative reserve designs over a range of different biodiversity and management goals and, ultimately, guide protected area acquisitions and management actions. Due to the systematic, evidence-based nature of these tools, conservation prioritization can help contribute to a transparent, inclusive, and more defensible decision making process.
A conservation planning exercise typically starts by defining a study area. This study area should encompass all the areas relevant to the decision maker or the hypothesis being tested. The extent of a study area could encompass a few important localities (e.g., Stigner et al. 2016), a single state (e.g., Kirkpatrick 1983), an entire country (Fuller et al. 2010), or the entire planet (Butchart et al. 2015). Next, the study area is divided into a set of discrete areas termed planning units. Each planning unit represents a discrete locality in the study area that can be managed independently of other areas. The general idea is that some combination of the planning units can be selected for conservation actions (e.g., protected area establishment, habitat restoration). Planning units are often created as square or hexagon cells that are sized according to the scale of the conservation actions, and the resolution of the data that underpin the planning exercise (but see Klein et al. 2009).
Cost data (or a surrogate thereof) are needed to inform the prioritization process. Specifically, these cost data describe the relative expenditure associated with managing each planning unit for conservation. For example, if the goal of the conservation planning exercise is to identify priority areas for expanding a local protected area system, then the cost data could represent the physical cost of purchasing the land. Alternatively, if such data are not available, then surrogate data are often used instead (e.g., human population density, opportunity cost of foregone commercial activities, or planning unit size).
Conservation planning exercises also require data on the biodiversity elements that are of conservation interest (termed conservation features). These features could be species (e.g., Neofelis nebulosa, the Clouded Leopard), populations, or habitats (e.g., mangroves or cloud forest). After identifying the set of relevant conservation features for a conservation planning exercise, spatially explicit data need to be obtained for each and every feature to describe their spatial distribution (e.g., habitat suitability data, probability of occurrence data, presence/absence data). This is important to ensure that conservation features are adequately covered (represented) by prioritizations. After assembling all the data, the next step is to define the conservation planning problem.
The prioritizr R package is designed to help you build and solve conservation planning problems. Specifically, prioritizations are generated by formulating a mathematical optimization problem and then solving it to generate a solution. These mathematical optimization problems are formulated using the planning unit data, cost data, and feature data, and with information related to the overarching aim of the prioritization process. In general, the goal of an optimization problem is to minimize (or maximize) an objective function that is calculated using a set of decision variables, subject to a series of constraints to ensure that solution exhibits specific properties. The objective function describes the quantity which we are trying to minimize (e.g., cost of the solution) or maximize (e.g., number of features conserved). The decision variables describe the entities that we can control, and indicate which planning units are selected for conservation management and which are not. The constraints can be thought of as rules that the decision variables need to follow. For example, a commonly used constraint is specifying that the solution must not exceed a certain budget.
A wide variety of approaches have been developed for solving optimization problems. Reserve design problems are frequently solved using simulated annealing (Kirkpatrick et al. 1983) or heuristics (Nicholls & Margules 1993; Moilanen 2007). These methods are conceptually simple and can be applied to a wide variety of optimization problems. However, they do not scale well for large or complex problems (Beyer et al. 2016). Additionally, these methods cannot tell you how close any given solution is to the optimal solution. The prioritizr package uses exact algorithms to efficiently solve conservation planning problems to within a pre-specified optimality gap. In other words, you can specify that you need the optimal solution (i.e., a gap of 0%) and the algorithms will, given enough time, find a solution that meets this criteria. In the past, exact algorithms have been too slow for conservation planning exercises (Pressey et al. 1996). However, improvements over the last decade mean that they are now much faster (Achterberg & Wunderling 2013; Beyer et al. 2016).
In this package, optimization problems are expressed using integer linear programming (ILP) so that they can be solved using (linear) exact algorithm solvers. The general form of an integer programming problem can be expressed in matrix notation using the following equation.
\[\text{Minimize} \space \boldsymbol{c}^\text{T} \boldsymbol{x} \space \space \text{subject to} \space A\boldsymbol{x} \space \Box \space \boldsymbol{b}\]
Where \(x\) is a vector of decision variables, \(c\) and \(b\) are vectors of known coefficients, and \(A\) is the constraint matrix. The final term specifies a series of structural constants and the \(\Box\) symbol is used to indicate that the relational operators for the constraints can be either \(\geq\), \(=\), or \(\leq\). In the context of a conservation planning problem, \(c\) could be used to represent the planning unit costs, \(A\) could be used to store the data showing the presence / absence (or amount) of each feature in each planning unit, \(b\) could be used to represent minimum amount of habitat required for each species in the solution, the \(\Box\) could be set to \(\geq\) symbols to indicate that the total amount of each feature in the solution must exceed the quantities in \(b\). But there are many other ways of formulating the reserve selection problem (Rodrigues et al. 2000).
A grammar for conservation planning
The prioritizr R package uses a grammar to describe elements
of conservation planning. This means that functions are organized into
verbs that relate to specific concepts. For example, all of the
functions used to specify the primary objective for optimization end
with the _objective suffix (e.g.,
add_min_set_objective() and
add_min_shortfall_objective()). By combining multiple
functions together, they can be used to formulate a complete
conservation planning problem. Specifically, the verbs for formulating
problems are described below.
- Create a new conservation planning problem by specifying the planning units, features, and management zones of conservation interest.
- Add a primary objective to a conservation planning problem (e.g., minimize overall cost).
- Add penalties to a conservation planning problem to penalize solutions according to specific metric (e.g., connectivity).
- Add targets to a conservation planning problem to specify how much of each feature should ideally be represented by solutions.
- Add constraints to a conservation planning problem to ensure that solutions exhibit specific properties (e.g., select specific planning units for protection).
- Add decisions to a conservation planning problem to specify the nature of the decisions in the problem (e.g., binary decisions indicate the planning units should be selected or not selected for management).
- Add a portfolio to a conservation planning problem to specify a methodological approach for generating multiple solutions (e.g., generate multiple solutions by finding 100 solutions within 10% of optimality).
- Add a solver to a conservation problem to specify the optimization software (e.g., Gurobi) for generating solutions and customize the optimization settings (e.g., generate an optimal solution).
After building a conservation planning problem, it can be solved to
generate a prioritization (using the solve() function).
There are also verbs available to help evaluate and interpret solutions.
These verbs are described below.
- Evaluate performance by computing summary statistics (e.g., overall cost, feature representation, or total boundary length).
- Evaluate the relative importance of planning units selected by a solution (e.g., based on irreplaceability).
Workflow
The general workflow when using the prioritizr R package
starts with creating a new conservation planning problem()
object using data. Specifically, the problem() object
should be constructed using data that specify the planning units,
biodiversity features, management zones (if applicable), and costs.
After creating a new problem() object, it can be
customized—by adding objectives, penalties, constraints, and other
information—to build a precise representation of the conservation
planning problem required, and then solved to obtain a solution.
All conservation planning problems require an objective. An objective
specifies the property which is used to compare different feasible
solutions. Simply put, the objective is the property of the solution
which should be maximized or minimized during the optimization process.
For instance, with the minimum set objective (specified using
add_min_set_objective()), we are seeking to minimize the
cost of the solution (similar to Marxan). On the other hand,
with the minimum shortfall objective (specified using
add_min_shortfall_objective()), we are seeking to minimize
the average target shortfall for all features represented in the
solution, subject to a budget.
Many objectives require representation targets (e.g., the minimum set
objective). These targets are a specialized set of constraints that
relate to the total quantity of each feature secured in the solution
(e.g., amount of suitable habitat or number of individuals). In the case
of the minimum set objective ( add_min_set_objective()),
they are used to ensure that solutions secure a sufficient quantity of
each feature, and in other objectives, such as the maximum features
objective ( add_max_features_objective()) they are used to
assess whether a feature has been adequately conserved by a candidate
solution. Targets can be expressed numerically as the total amount
required for a given feature (using
add_absolute_targets()), or as a proportion of the total
amount found in the planning units (using
add_relative_targets()). Note that not all objectives
require targets, and a warning will be thrown if an attempt is made to
add targets to a problem with an objective that does not use them.
Constraints and penalties can be added to a conservation planning
problem to ensure that solutions exhibit a specific property or penalize
solutions which don’t exhibit a specific property (respectively). The
difference between constraints and penalties, strictly speaking, is
constraints are used to rule out potential solutions that don’t exhibit
a specific property. For instance, constraints can be used to ensure
that specific planning units are selected in the solution for
prioritization (using add_locked_in_constraints()) or not
selected in the solution for prioritization (using
add_locked_out_constraints()). On the other hand, penalties
are combined with the objective of a problem, with a penalty factor, and
the overall objective of the problem then becomes to minimize (or
maximize) the primary objective function and the penalty function. For
example, penalties can be added to a problem to penalize solutions that
are excessively fragmented (using
add_boundary_penalties()). These penalties have a
penalty argument that specifies the relative importance of
having spatially clustered solutions. When the argument to
penalty is high, then solutions which are less fragmented
are valued more highly – even if they cost more – and when the argument
to penalty is low, then the solutions which are more
fragmented are valued less highly.
After building a conservation problem, it can then be solved to obtain a solution (or portfolio of solutions if desired). The solution is returned in the same format as the planning unit data used to construct the problem. For instance, this means that if raster data was used to initialize the problem, then the solution will also be output in raster format. This can be very helpful when it comes to interpreting and visualizing solutions because it means that the solution data does not first have to be merged with spatial data before they can be plotted on a map.
Usage
Here we will provide an introduction to using the prioritizr R package to build and solve a conservation planning problem. Please note that we will not discuss conservation planning with multiple zones in this vignette, for more information on working with multiple management zones please see the Management zones tutorial.
First, we will load the prioritizr package.
# load package
library(prioritizr)Data
Now we will load some built-in data sets that are distributed with the prioritizr R package. This package contains several different planning unit data sets. To provide a comprehensive overview of the different ways that we can initialize a conservation planning problem, we will load each of them.
First, we will load the raster planning unit data
(sim_pu_raster). Here, the planning units are represented
as a single-layer raster object (i.e., a terra::rast()
object) and each pixel corresponds to the spatial extent of each panning
unit. The pixel values correspond to the acquisition costs of each
planning unit.
# load raster planning unit data
sim_pu_raster <- get_sim_pu_raster()
# print data
print(sim_pu_raster)## class : SpatRaster
## dimensions : 10, 10, 1 (nrow, ncol, nlyr)
## resolution : 0.1, 0.1 (x, y)
## extent : 0, 1, 0, 1 (xmin, xmax, ymin, ymax)
## coord. ref. : Undefined Cartesian SRS
## source : sim_pu_raster.tif
## name : layer
## min value : 190.1328
## max value : 215.8638
# plot data
plot(
sim_pu_raster, main = "Planning unit costs",
xlim = c(-0.1, 1.1), ylim = c(-0.1, 1.1), axes = FALSE
)
Secondly, we will load one of the spatial vector planning unit data
sets (sim_pu_polygons). Here, each polygon (i.e., feature
using ArcGIS terminology) corresponds to a different planning unit. This
data set has an attribute table that contains additional information
about each polygon. Namely, the cost field (column) in the
attribute table contains the acquisition cost for each planning
unit.
# load polygon planning unit data
sim_pu_polygons <- get_sim_pu_polygons()
# print data
print(sim_pu_polygons)## Simple feature collection with 90 features and 3 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 0 ymin: 0 xmax: 1 ymax: 1
## Projected CRS: Undefined Cartesian SRS
## # A tibble: 90 × 4
## cost locked_in locked_out geom
## * <dbl> <lgl> <lgl> <POLYGON [m]>
## 1 216. FALSE FALSE ((0 1, 0.1 1, 0.1 0.9, 0 0.9, 0 1))
## 2 213. FALSE FALSE ((0.1 1, 0.2 1, 0.2 0.9, 0.1 0.9, 0.1 1))
## 3 207. FALSE FALSE ((0.2 1, 0.3 1, 0.3 0.9, 0.2 0.9, 0.2 1))
## 4 209. FALSE TRUE ((0.3 1, 0.4 1, 0.4 0.9, 0.3 0.9, 0.3 1))
## 5 214. FALSE FALSE ((0.4 1, 0.5 1, 0.5 0.9, 0.4 0.9, 0.4 1))
## 6 214. FALSE FALSE ((0.5 1, 0.6 1, 0.6 0.9, 0.5 0.9, 0.5 1))
## 7 210. FALSE FALSE ((0.6 1, 0.7 1, 0.7 0.9, 0.6 0.9, 0.6 1))
## 8 211. FALSE TRUE ((0.7 1, 0.8 1, 0.8 0.9, 0.7 0.9, 0.7 1))
## 9 210. FALSE FALSE ((0.8 1, 0.9 1, 0.9 0.9, 0.8 0.9, 0.8 1))
## 10 204. FALSE FALSE ((0.9 1, 1 1, 1 0.9, 0.9 0.9, 0.9 1))
## # ℹ 80 more rows
# plot the planning units, and color them according to cost values
plot(sim_pu_polygons[, "cost"])
Thirdly, we will load some planning unit data stored in tabular
format (i.e., data.frame format). For those familiar with
Marxan or dealing with very large conservation planning
problems (> 10 million planning units), it may be useful to work with
data in this format because it does not contain any spatial information
which will reduce computational burden. When using tabular data to
initialize conservation planning problems, the data must follow the
conventions used by Marxan. Specifically, each row in the
planning unit table must correspond to a different planning unit. The
table must also have an “id” column to provide a unique integer
identifier for each planning unit, and it must also have a column that
indicates the cost of each planning unit. For more information, please
see the official Marxan
documentation.
# specify file path for planning unit data
pu_path <- system.file("extdata/marxan/input/pu.dat", package = "prioritizr")
# load in the tabular planning unit data
# note that we use the data.table::fread function, as opposed to the read.csv
# function, because it is much faster
pu_dat <- data.table::fread(pu_path, data.table = FALSE)
# preview first six rows of the tabular planning unit data
# note that it has some extra columns other than id and cost as per the
# Marxan format
head(pu_dat)## id cost status xloc yloc
## 1 3 0.000 0 1116623 -4493479
## 2 30 7527.275 3 1110623 -4496943
## 3 56 37349.075 0 1092623 -4500408
## 4 58 16959.021 0 1116623 -4500408
## 5 84 34220.256 0 1098623 -4503872
## 6 85 178907.584 0 1110623 -4503872Finally, we will load data showing the spatial distribution of the
conservation features. Our conservation features
(sim_features) are represented as a multi-layer raster
object (i.e., a terra::rast() object), where each layer
corresponds to a different feature. The pixel values in each layer
correspond to the amount of suitable habitat available in a given
planning unit. Note that our planning unit and our feature data have
exactly the same spatial properties (i.e., resolution, extent,
coordinate reference system) so their pixels line up perfectly.
# load feature data
sim_features <- get_sim_features()
# plot the distribution of suitable habitat for each feature
plot(
sim_features, nr = 2,
axes = FALSE, xlim = c(-0.1, 1.1), ylim = c(-0.1, 1.1)
)
Initialize a problem
After having loaded our planning unit and feature data, we will now
try initializing some conservation planning problems. There are a lot of
different ways to initialize a conservation planning problem, so here we
will just showcase a few of the more commonly used methods. For an
exhaustive description of all the ways you can initialize a conservation
problem, see the help file for the problem() function.
First off, we will initialize a conservation planning problem using the
raster data.
## A conservation problem (<ConservationProblem>)
## ├•data
## │├•features: "feature_1", "feature_2", "feature_3", "feature_4", and "feature_5" (5 total)
## │└•planning units:
## │ ├•data: <SpatRaster> (90 total)
## │ ├•costs: continuous values (between 190.1328 and 215.8638)
## │ ├•extent: 0, 0, 1, 1 (xmin, ymin, xmax, ymax)
## │ └•CRS: Undefined Cartesian SRS (projected)
## ├•formulation
## │├•objective: none specified
## │├•penalties: none specified
## │├•targets: none specified
## │├•constraints: none specified
## │└•decisions: binary decision
## └•optimization
## ├•portfolio: default portfolio
## └•solver: gurobi solver (`gap` = 0.1, `time_limit` = 2147483647, `first_feasible` = FALSE, …)
## # ℹ Use `summary(...)` to see complete formulation.
# print number of planning units
number_of_planning_units(p1)## [1] 90
# print number of features
number_of_features(p1)## [1] 5Generally, we recommend initializing problems using raster data where
possible. This is because the problem() function needs to
calculate the amount of each feature in each planning unit, and by
providing both the planning unit and feature data in raster format with
the same spatial resolution, extents, and coordinate systems, this means
that the problem() function does not need to do any
geoprocessing behind the scenes. But sometimes we can’t use raster
planning unit data, because our planning units aren’t equal-sized grid
cells. So, below is an example showing how we can initialize a
conservation planning problem using planning units that are formatted as
spatial vector data. Note that we could reduce run-time by pre-computing
the amount of each feature in each planning unit and storing the data in
the attribute table (e.g., by performing zonal statistics with
R or ESRI ArcGIS), and then passing in the names of
the columns as an argument to the problem() function (see
Examples section for problem() for details).
# create problem with spatial vector data
# note that we have to specify which column in the attribute table contains
# the cost data
p2 <- problem(sim_pu_polygons, sim_features, cost_column = "cost")
# print problem
print(p2)## A conservation problem (<ConservationProblem>)
## ├•data
## │├•features: "feature_1", "feature_2", "feature_3", "feature_4", and "feature_5" (5 total)
## │└•planning units:
## │ ├•data: <sftbl_dftbldata.frame> (90 total)
## │ ├•costs: continuous values (between 190.1328 and 215.8638)
## │ ├•extent: 0, 0, 1, 1 (xmin, ymin, xmax, ymax)
## │ └•CRS: Undefined Cartesian SRS (projected)
## ├•formulation
## │├•objective: none specified
## │├•penalties: none specified
## │├•targets: none specified
## │├•constraints: none specified
## │└•decisions: binary decision
## └•optimization
## ├•portfolio: default portfolio
## └•solver: gurobi solver (`gap` = 0.1, `time_limit` = 2147483647, `first_feasible` = FALSE, …)
## # ℹ Use `summary(...)` to see complete formulation.We can also initialize a conservation planning problem using tabular
planning unit data (i.e., data.frame format). Since the
tabular planning unit data does not contain any spatial information, we
also have to provide the feature data in tabular format (i.e.,
data.frame format) and data showing the amount of each
feature in each planning unit in tabular format (i.e.,
data.frame format). The feature data must have an “id”
column containing a unique integer identifier for each feature, and the
planning unit by feature data must contain the following three columns:
pu corresponding to the planning unit identifiers,
species corresponding to the feature identifiers, and
amount showing the amount of a given feature in a given
planning unit.
# set file path for feature data
spec_path <- system.file(
"extdata/marxan/input/spec.dat", package = "prioritizr"
)
# load in feature data
spec_dat <- data.table::fread(spec_path, data.table = FALSE)
# print first six rows of the data
# note that it contains extra columns
head(spec_dat)## id prop spf name
## 1 10 0.3 1 bird1
## 2 11 0.3 1 nvis2
## 3 12 0.3 1 nvis8
## 4 13 0.3 1 nvis9
## 5 14 0.3 1 nvis14
## 6 15 0.3 1 nvis20
# set file path for planning unit vs. feature data
puvspr_path <- system.file(
"extdata/marxan/input/puvspr.dat", package = "prioritizr"
)
# load in planning unit vs feature data
puvspr_dat <- data.table::fread(puvspr_path, data.table = FALSE)
# print first six rows of the data
head(puvspr_dat)## species pu amount
## 1 26 56 120.344884
## 2 26 58 45.167010
## 3 26 84 68.047375
## 4 26 85 9.735624
## 5 26 86 7.803476
## 6 26 111 478.327417
# create problem
p3 <- problem(pu_dat, spec_dat, cost_column = "cost", rij = puvspr_dat)
# print problem
print(p3)## A conservation problem (<ConservationProblem>)
## ├•data
## │├•features: "bird1", "nvis2", "nvis8", "nvis9", "nvis14", "nvis20", … (17 total)
## │└•planning units:
## │ ├•data: <data.frame> (1751 total)
## │ ├•costs: continuous values (between 0 and 415692.1938)
## │ ├•extent: NA
## │ └•CRS: NA
## ├•formulation
## │├•objective: none specified
## │├•penalties: none specified
## │├•targets: none specified
## │├•constraints: none specified
## │└•decisions: binary decision
## └•optimization
## ├•portfolio: default portfolio
## └•solver: gurobi solver (`gap` = 0.1, `time_limit` = 2147483647, `first_feasible` = FALSE, …)
## # ℹ Use `summary(...)` to see complete formulation.For more information on initializing problems, please see the help
page for the problem() function (which you can open by
entering the code: ?problem). Now that we have initialized
a conservation planning problem, we will show you how you can customize
it to suit the exact needs of your conservation planning scenario.
Although we initialized the conservation planning problems using several
different methods, moving forward, we will only use raster-based
planning unit data to keep things simple.
Add an objective
The next step is to add a primary objective to the problem. A problem objective is used to specify the primary goal of the problem (i.e., the quantity that is to be maximized or minimized). All conservation planning problems involve minimizing or maximizing some kind of objective. For instance, we might require a solution that conserves enough habitat for each species while minimizing the overall cost of the reserve network. Alternatively, we might require a solution that maximizes the number of conserved species while ensuring that the cost of the reserve network does not exceed the budget. Please note that objectives are added in the same way regardless of the type of data used to initialize the problem. The following objectives are available.
- Minimum set objective: Minimize the cost of the solution whilst ensuring that all targets are met (Rodrigues et al. 2000). This objective is similar to that used in Marxan (Ball et al. 2009). For example, we can add a minimum set objective to a problem using the following code.
# create a new problem that has the minimum set objective
p3 <-
problem(sim_pu_raster, sim_features) %>%
add_min_set_objective()
# print problem
print(p3)## A conservation problem (<ConservationProblem>)
## ├•data
## │├•features: "feature_1", "feature_2", "feature_3", "feature_4", and "feature_5" (5 total)
## │└•planning units:
## │ ├•data: <SpatRaster> (90 total)
## │ ├•costs: continuous values (between 190.1328 and 215.8638)
## │ ├•extent: 0, 0, 1, 1 (xmin, ymin, xmax, ymax)
## │ └•CRS: Undefined Cartesian SRS (projected)
## ├•formulation
## │├•objective: minimum set objective
## │├•penalties: none specified
## │├•targets: none specified
## │├•constraints: none specified
## │└•decisions: binary decision
## └•optimization
## ├•portfolio: default portfolio
## └•solver: gurobi solver (`gap` = 0.1, `time_limit` = 2147483647, `first_feasible` = FALSE, …)
## # ℹ Use `summary(...)` to see complete formulation.- Maximum cover objective: Represent at least one instance of as many features as possible within a given budget (Church et al. 1996).
# create a new problem that has the maximum coverage objective and a budget
# of 5000
p4 <-
problem(sim_pu_raster, sim_features) %>%
add_max_cover_objective(5000)
# print problem
print(p4)## A conservation problem (<ConservationProblem>)
## ├•data
## │├•features: "feature_1", "feature_2", "feature_3", "feature_4", and "feature_5" (5 total)
## │└•planning units:
## │ ├•data: <SpatRaster> (90 total)
## │ ├•costs: continuous values (between 190.1328 and 215.8638)
## │ ├•extent: 0, 0, 1, 1 (xmin, ymin, xmax, ymax)
## │ └•CRS: Undefined Cartesian SRS (projected)
## ├•formulation
## │├•objective: maximum coverage objective (`budget` = 5000)
## │├•penalties: none specified
## │├•targets: none specified
## │├•constraints: none specified
## │└•decisions: binary decision
## └•optimization
## ├•portfolio: default portfolio
## └•solver: gurobi solver (`gap` = 0.1, `time_limit` = 2147483647, `first_feasible` = FALSE, …)
## # ℹ Use `summary(...)` to see complete formulation.- Maximum features objective: Fulfill as many targets as possible while ensuring that the cost of the solution does not exceed a budget (inspired by Cabeza & Moilanen 2001). This object is similar to the maximum cover objective except that we have the option of later specifying targets for each feature. In practice, this objective is more useful than the maximum cover objective because features often require a certain amount of area for them to persist and simply capturing a single instance of habitat for each feature is generally unlikely to enhance their long-term persistence.
# create a new problem that has the maximum features objective and a budget
# of 5000
p5 <-
problem(sim_pu_raster, sim_features) %>%
add_max_features_objective(budget = 5000)
# print problem
print(p5)## A conservation problem (<ConservationProblem>)
## ├•data
## │├•features: "feature_1", "feature_2", "feature_3", "feature_4", and "feature_5" (5 total)
## │└•planning units:
## │ ├•data: <SpatRaster> (90 total)
## │ ├•costs: continuous values (between 190.1328 and 215.8638)
## │ ├•extent: 0, 0, 1, 1 (xmin, ymin, xmax, ymax)
## │ └•CRS: Undefined Cartesian SRS (projected)
## ├•formulation
## │├•objective: maximum representation objective (`budget` = 5000)
## │├•penalties: none specified
## │├•targets: none specified
## │├•constraints: none specified
## │└•decisions: binary decision
## └•optimization
## ├•portfolio: default portfolio
## └•solver: gurobi solver (`gap` = 0.1, `time_limit` = 2147483647, `first_feasible` = FALSE, …)
## # ℹ Use `summary(...)` to see complete formulation.- Minimum shortfall objective: Minimize the overall (weighted sum) shortfall for as many targets as possible while ensuring that the cost of the solution does not exceed a budget. In practice, this objective is useful when there is a large amount of left-over budget when using the maximum feature representation objective and the remaining funds need to be allocated to places that will enhance the representation of features with unmet targets.
# create a new problem that has the minimum shortfall objective and a budget
# of 5000
p6 <-
problem(sim_pu_raster, sim_features) %>%
add_min_shortfall_objective(budget = 5000)
# print problem
print(p6)## A conservation problem (<ConservationProblem>)
## ├•data
## │├•features: "feature_1", "feature_2", "feature_3", "feature_4", and "feature_5" (5 total)
## │└•planning units:
## │ ├•data: <SpatRaster> (90 total)
## │ ├•costs: continuous values (between 190.1328 and 215.8638)
## │ ├•extent: 0, 0, 1, 1 (xmin, ymin, xmax, ymax)
## │ └•CRS: Undefined Cartesian SRS (projected)
## ├•formulation
## │├•objective: minimum shortfall objective (`budget` = 5000)
## │├•penalties: none specified
## │├•targets: none specified
## │├•constraints: none specified
## │└•decisions: binary decision
## └•optimization
## ├•portfolio: default portfolio
## └•solver: gurobi solver (`gap` = 0.1, `time_limit` = 2147483647, `first_feasible` = FALSE, …)
## # ℹ Use `summary(...)` to see complete formulation.- Minimum largest shortfall objective: Minimize the largest (maximum) shortfall while ensuring that the cost of the solution does not exceed a budget. In practice, this objective useful when the minimum shortfall objective returns solutions that focus too much on representing a small number of features (e.g., because they occur in much cheaper planning units), and solutions are needed to spread conservation effort out more evenly among all features—even if it means that all features will have (relatively) poor representation.
# create a new problem that has the minimum largest shortfall objective and a
# budget of 5000
p7 <-
problem(sim_pu_raster, sim_features) %>%
add_min_largest_shortfall_objective(budget = 5000)
# print problem
print(p7)## A conservation problem (<ConservationProblem>)
## ├•data
## │├•features: "feature_1", "feature_2", "feature_3", "feature_4", and "feature_5" (5 total)
## │└•planning units:
## │ ├•data: <SpatRaster> (90 total)
## │ ├•costs: continuous values (between 190.1328 and 215.8638)
## │ ├•extent: 0, 0, 1, 1 (xmin, ymin, xmax, ymax)
## │ └•CRS: Undefined Cartesian SRS (projected)
## ├•formulation
## │├•objective: minimum largest shortfall objective (`budget` = 5000)
## │├•penalties: none specified
## │├•targets: none specified
## │├•constraints: none specified
## │└•decisions: binary decision
## └•optimization
## ├•portfolio: default portfolio
## └•solver: gurobi solver (`gap` = 0.1, `time_limit` = 2147483647, `first_feasible` = FALSE, …)
## # ℹ Use `summary(...)` to see complete formulation.-
Maximum phylogenetic diversity objective: Maximize
the phylogenetic diversity of the features represented in the solution
subject to a budget (inspired by Faith 1992;
Rodrigues & Gaston 2002). This objective is similar to the
maximum features objective except that emphasis is placed on protecting
features which are associated with a diverse range of evolutionary
histories. The package contains a simulated phylogeny that can be used
with the simulated feature data (
sim_phylogny).
# load simulated phylogeny data
sim_phylogeny <- get_sim_phylogeny()
# create a new problem that has the maximum phylogenetic diversity
# objective and a budget of 5000
p8 <-
problem(sim_pu_raster, sim_features) %>%
add_max_phylo_div_objective(budget = 5000, tree = sim_phylogeny)
# print problem
print(p8)## A conservation problem (<ConservationProblem>)
## ├•data
## │├•features: "feature_1", "feature_2", "feature_3", "feature_4", and "feature_5" (5 total)
## │└•planning units:
## │ ├•data: <SpatRaster> (90 total)
## │ ├•costs: continuous values (between 190.1328 and 215.8638)
## │ ├•extent: 0, 0, 1, 1 (xmin, ymin, xmax, ymax)
## │ └•CRS: Undefined Cartesian SRS (projected)
## ├•formulation
## │├•objective: phylogenetic diversity objective (`budget` = 5000, `tree` = phylogenetic tree (branch lengths between 0.041 and 1.4211))
## │├•penalties: none specified
## │├•targets: none specified
## │├•constraints: none specified
## │└•decisions: binary decision
## └•optimization
## ├•portfolio: default portfolio
## └•solver: gurobi solver (`gap` = 0.1, `time_limit` = 2147483647, `first_feasible` = FALSE, …)
## # ℹ Use `summary(...)` to see complete formulation.- Maximum phylogenetic endemism objective: Maximize the phylogenetic endemism of the features represented in the solution subject to a budget (inspired by Faith 1992; Rodrigues & Gaston 2002; Rosauer et al. 2009). This objective is similar to the maximum phylogenetic diversity except that emphasis is placed on conserving features that are associated with geographically restricted periods of evolutionary history rather than a diverse range of evolutionary histories.
# load simulated phylogeny data
sim_phylogeny <- get_sim_phylogeny()
# create a new problem that has the maximum phylogenetic diversity
# objective and a budget of 5000
p9 <-
problem(sim_pu_raster, sim_features) %>%
add_max_phylo_end_objective(budget = 5000, tree = sim_phylogeny)
# print problem
print(p9)## A conservation problem (<ConservationProblem>)
## ├•data
## │├•features: "feature_1", "feature_2", "feature_3", "feature_4", and "feature_5" (5 total)
## │└•planning units:
## │ ├•data: <SpatRaster> (90 total)
## │ ├•costs: continuous values (between 190.1328 and 215.8638)
## │ ├•extent: 0, 0, 1, 1 (xmin, ymin, xmax, ymax)
## │ └•CRS: Undefined Cartesian SRS (projected)
## ├•formulation
## │├•objective: phylogenetic endemism objective (`budget` = 5000, `tree` = phylogenetic tree (branch lengths between 0.041 and 1.4211))
## │├•penalties: none specified
## │├•targets: none specified
## │├•constraints: none specified
## │└•decisions: binary decision
## └•optimization
## ├•portfolio: default portfolio
## └•solver: gurobi solver (`gap` = 0.1, `time_limit` = 2147483647, `first_feasible` = FALSE, …)
## # ℹ Use `summary(...)` to see complete formulation.- Maximum utility objective: Secure as much of the features as possible without exceeding a budget. This objective is functionally equivalent to selecting the planning units with the greatest amounts of each feature (e.g., species richness). Generally, we don’t encourage the use of this objective because it will only rarely identify complementary solutions–solutions which adequately conserve a range of different features—except perhaps to explore trade-offs or provide a baseline solution with which to compare other solutions.
# create a new problem that has the maximum utility objective and a budget
# of 5000
p10 <-
problem(sim_pu_raster, sim_features) %>%
add_max_utility_objective(budget = 5000)
# print problem
print(p10)## A conservation problem (<ConservationProblem>)
## ├•data
## │├•features: "feature_1", "feature_2", "feature_3", "feature_4", and "feature_5" (5 total)
## │└•planning units:
## │ ├•data: <SpatRaster> (90 total)
## │ ├•costs: continuous values (between 190.1328 and 215.8638)
## │ ├•extent: 0, 0, 1, 1 (xmin, ymin, xmax, ymax)
## │ └•CRS: Undefined Cartesian SRS (projected)
## ├•formulation
## │├•objective: maximum utility objective (`budget` = 5000)
## │├•penalties: none specified
## │├•targets: none specified
## │├•constraints: none specified
## │└•decisions: binary decision
## └•optimization
## ├•portfolio: default portfolio
## └•solver: gurobi solver (`gap` = 0.1, `time_limit` = 2147483647, `first_feasible` = FALSE, …)
## # ℹ Use `summary(...)` to see complete formulation.Add targets
Most conservation planning problems require targets. Targets are used to specify the minimum amount or proportion of a feature’s distribution that needs to be protected in the solution. For example, we may want to develop a reserve network that will secure 20% of the distribution for each feature for minimal cost. The following methods are available for specifying targets.
- Absolute targets: Targets are expressed as the total amount of each feature in the study area that need to be secured. For example, if we had binary feature data that showed the absence or presence of suitable habitat across the study area, we could set an absolute target of 5 to mean that we require 5 planning units with suitable habitat in the solution.
# create a problem with targets which specify that the solution must conserve
# a sum total of 3 units of suitable habitat for each feature
p11 <-
problem(sim_pu_raster, sim_features) %>%
add_min_set_objective() %>%
add_absolute_targets(3)
# print problem
print(p11)## A conservation problem (<ConservationProblem>)
## ├•data
## │├•features: "feature_1", "feature_2", "feature_3", "feature_4", and "feature_5" (5 total)
## │└•planning units:
## │ ├•data: <SpatRaster> (90 total)
## │ ├•costs: continuous values (between 190.1328 and 215.8638)
## │ ├•extent: 0, 0, 1, 1 (xmin, ymin, xmax, ymax)
## │ └•CRS: Undefined Cartesian SRS (projected)
## ├•formulation
## │├•objective: minimum set objective
## │├•penalties: none specified
## │├•targets: absolute targets (between 3 and 3)
## │├•constraints: none specified
## │└•decisions: binary decision
## └•optimization
## ├•portfolio: default portfolio
## └•solver: gurobi solver (`gap` = 0.1, `time_limit` = 2147483647, `first_feasible` = FALSE, …)
## # ℹ Use `summary(...)` to see complete formulation.- Relative targets: Targets are set as a proportion (between 0 and 1) of the total amount of each feature in the study area. For example, if we had binary feature data and the feature occupied a total of 20 planning units in the study area, we could set a relative target of 0.5 (i.e., 50%) to specify that the solution must secure 10 planning units for the feature. We could alternatively specify an absolute target of 10 to achieve the same result, but sometimes proportions are easier to work with.
# create a problem with the minimum set objective and relative targets of 10 %
# for each feature
p12 <-
problem(sim_pu_raster, sim_features) %>%
add_min_set_objective() %>%
add_relative_targets(0.1)
# print problem
print(p12)## A conservation problem (<ConservationProblem>)
## ├•data
## │├•features: "feature_1", "feature_2", "feature_3", "feature_4", and "feature_5" (5 total)
## │└•planning units:
## │ ├•data: <SpatRaster> (90 total)
## │ ├•costs: continuous values (between 190.1328 and 215.8638)
## │ ├•extent: 0, 0, 1, 1 (xmin, ymin, xmax, ymax)
## │ └•CRS: Undefined Cartesian SRS (projected)
## ├•formulation
## │├•objective: minimum set objective
## │├•penalties: none specified
## │├•targets: relative targets (between 0.1 and 0.1)
## │├•constraints: none specified
## │└•decisions: binary decision
## └•optimization
## ├•portfolio: default portfolio
## └•solver: gurobi solver (`gap` = 0.1, `time_limit` = 2147483647, `first_feasible` = FALSE, …)
## # ℹ Use `summary(...)` to see complete formulation.
# create a problem with targets which specify that we need 10 % of the habitat
# for the first feature, 15 % for the second feature, 20 % for the third feature
# 25 % for the fourth feature and 30 % of the habitat for the fifth feature
targets <- c(0.1, 0.15, 0.2, 0.25, 0.3)
p13 <-
problem(sim_pu_raster, sim_features) %>%
add_min_set_objective() %>%
add_relative_targets(targets)
# print problem
print(p13)## A conservation problem (<ConservationProblem>)
## ├•data
## │├•features: "feature_1", "feature_2", "feature_3", "feature_4", and "feature_5" (5 total)
## │└•planning units:
## │ ├•data: <SpatRaster> (90 total)
## │ ├•costs: continuous values (between 190.1328 and 215.8638)
## │ ├•extent: 0, 0, 1, 1 (xmin, ymin, xmax, ymax)
## │ └•CRS: Undefined Cartesian SRS (projected)
## ├•formulation
## │├•objective: minimum set objective
## │├•penalties: none specified
## │├•targets: relative targets (between 0.1 and 0.3)
## │├•constraints: none specified
## │└•decisions: binary decision
## └•optimization
## ├•portfolio: default portfolio
## └•solver: gurobi solver (`gap` = 0.1, `time_limit` = 2147483647, `first_feasible` = FALSE, …)
## # ℹ Use `summary(...)` to see complete formulation.- Log-linear targets: Targets are expressed using scaling factors and log-linear interpolation. This method for specifying targets is commonly used for global prioritization analyses (Rodrigues et al. 2004).
# create problem with added log-linear targets
p14 <-
problem(sim_pu_raster, sim_features) %>%
add_min_set_objective() %>%
add_loglinear_targets(10, 0.9, 100, 0.2)
# print problem
print(p14)## A conservation problem (<ConservationProblem>)
## ├•data
## │├•features: "feature_1", "feature_2", "feature_3", "feature_4", and "feature_5" (5 total)
## │└•planning units:
## │ ├•data: <SpatRaster> (90 total)
## │ ├•costs: continuous values (between 190.1328 and 215.8638)
## │ ├•extent: 0, 0, 1, 1 (xmin, ymin, xmax, ymax)
## │ └•CRS: Undefined Cartesian SRS (projected)
## ├•formulation
## │├•objective: minimum set objective
## │├•penalties: none specified
## │├•targets: absolute targets (between 17.2905 and 21.5906)
## │├•constraints: none specified
## │└•decisions: binary decision
## └•optimization
## ├•portfolio: default portfolio
## └•solver: gurobi solver (`gap` = 0.1, `time_limit` = 2147483647, `first_feasible` = FALSE, …)
## # ℹ Use `summary(...)` to see complete formulation.- Manual targets: Targets are manually specified. This is only really recommended for advanced users or problems that involve multiple management zones. See the Management zones tutorial for more information on these targets.
As with the functions for specifying the objective of a problem, if we try adding multiple targets to a problem, only the most recently added set of targets are used.
Add constraints
A constraint can be added to a conservation planning problem to ensure that all solutions exhibit a specific property. For example, they can be used to make sure that all solutions select a specific planning unit or that all selected planning units in the solution follow a certain configuration. The following constraints are available.
- Locked in constraints: Add constraints to ensure that certain planning units are prioritized in the solution. For example, it may be desirable to lock in planning units that are inside existing protected areas so that the solution fills in the gaps in the existing reserve network.
# create problem with constraints which specify that the first planning unit
# must be selected in the solution
p15 <-
problem(sim_pu_raster, sim_features) %>%
add_min_set_objective() %>%
add_relative_targets(0.1) %>%
add_locked_in_constraints(1)
# print problem
print(p15)## A conservation problem (<ConservationProblem>)
## ├•data
## │├•features: "feature_1", "feature_2", "feature_3", "feature_4", and "feature_5" (5 total)
## │└•planning units:
## │ ├•data: <SpatRaster> (90 total)
## │ ├•costs: continuous values (between 190.1328 and 215.8638)
## │ ├•extent: 0, 0, 1, 1 (xmin, ymin, xmax, ymax)
## │ └•CRS: Undefined Cartesian SRS (projected)
## ├•formulation
## │├•objective: minimum set objective
## │├•penalties: none specified
## │├•targets: relative targets (between 0.1 and 0.1)
## │├•constraints:
## ││└•1: locked in constraints (1 planning units)
## │└•decisions: binary decision
## └•optimization
## ├•portfolio: default portfolio
## └•solver: gurobi solver (`gap` = 0.1, `time_limit` = 2147483647, `first_feasible` = FALSE, …)
## # ℹ Use `summary(...)` to see complete formulation.- Locked out constraints: Add constraints to ensure that certain planning units are not prioritized in the solution. For example, it may be useful to lock out planning units that have been degraded and are not suitable for conserving species.
# create problem with constraints which specify that the second planning unit
# must not be selected in the solution
p16 <-
problem(sim_pu_raster, sim_features) %>%
add_min_set_objective() %>%
add_relative_targets(0.1) %>%
add_locked_out_constraints(2)
# print problem
print(p16)## A conservation problem (<ConservationProblem>)
## ├•data
## │├•features: "feature_1", "feature_2", "feature_3", "feature_4", and "feature_5" (5 total)
## │└•planning units:
## │ ├•data: <SpatRaster> (90 total)
## │ ├•costs: continuous values (between 190.1328 and 215.8638)
## │ ├•extent: 0, 0, 1, 1 (xmin, ymin, xmax, ymax)
## │ └•CRS: Undefined Cartesian SRS (projected)
## ├•formulation
## │├•objective: minimum set objective
## │├•penalties: none specified
## │├•targets: relative targets (between 0.1 and 0.1)
## │├•constraints:
## ││└•1: locked out constraints (1 planning units)
## │└•decisions: binary decision
## └•optimization
## ├•portfolio: default portfolio
## └•solver: gurobi solver (`gap` = 0.1, `time_limit` = 2147483647, `first_feasible` = FALSE, …)
## # ℹ Use `summary(...)` to see complete formulation.- Neighbor constraints: Add constraints to a conservation problem to ensure that all selected planning units have at least a certain number of neighbors.
# create problem with constraints which specify that all selected planning units
# in the solution must have at least 1 neighbor
p17 <-
problem(sim_pu_raster, sim_features) %>%
add_min_set_objective() %>%
add_relative_targets(0.1) %>%
add_neighbor_constraints(1)
# print problem
print(p17)## A conservation problem (<ConservationProblem>)
## ├•data
## │├•features: "feature_1", "feature_2", "feature_3", "feature_4", and "feature_5" (5 total)
## │└•planning units:
## │ ├•data: <SpatRaster> (90 total)
## │ ├•costs: continuous values (between 190.1328 and 215.8638)
## │ ├•extent: 0, 0, 1, 1 (xmin, ymin, xmax, ymax)
## │ └•CRS: Undefined Cartesian SRS (projected)
## ├•formulation
## │├•objective: minimum set objective
## │├•penalties: none specified
## │├•targets: relative targets (between 0.1 and 0.1)
## │├•constraints:
## ││└•1: neighbor constraints (`k` = 1, `clamp` = TRUE, …)
## │└•decisions: binary decision
## └•optimization
## ├•portfolio: default portfolio
## └•solver: gurobi solver (`gap` = 0.1, `time_limit` = 2147483647, `first_feasible` = FALSE, …)
## # ℹ Use `summary(...)` to see complete formulation.- Contiguity constraints: Add constraints to a conservation problem to ensure that all selected planning units are spatially connected to each other and form spatially contiguous unit.
# create problem with constraints which specify that all selected planning units
# in the solution must form a single contiguous unit
p18 <-
problem(sim_pu_raster, sim_features) %>%
add_min_set_objective() %>%
add_relative_targets(0.1) %>%
add_contiguity_constraints()
# print problem
print(p18)## A conservation problem (<ConservationProblem>)
## ├•data
## │├•features: "feature_1", "feature_2", "feature_3", "feature_4", and "feature_5" (5 total)
## │└•planning units:
## │ ├•data: <SpatRaster> (90 total)
## │ ├•costs: continuous values (between 190.1328 and 215.8638)
## │ ├•extent: 0, 0, 1, 1 (xmin, ymin, xmax, ymax)
## │ └•CRS: Undefined Cartesian SRS (projected)
## ├•formulation
## │├•objective: minimum set objective
## │├•penalties: none specified
## │├•targets: relative targets (between 0.1 and 0.1)
## │├•constraints:
## ││└•1: contiguity constraints
## │└•decisions: binary decision
## └•optimization
## ├•portfolio: default portfolio
## └•solver: gurobi solver (`gap` = 0.1, `time_limit` = 2147483647, `first_feasible` = FALSE, …)
## # ℹ Use `summary(...)` to see complete formulation.-
Feature contiguity constraints: Add constraints to
ensure that each feature is represented in a contiguous unit of
dispersible habitat. These constraints are a more advanced version of
those implemented in the
add_contiguity_constraintsfunction, because they ensure that each feature is represented in a contiguous unit and not that the entire solution should form a contiguous unit.
# create problem with constraints which specify that the planning units used
# to conserve each feature must form a contiguous unit
p19 <-
problem(sim_pu_raster, sim_features) %>%
add_min_set_objective() %>%
add_relative_targets(0.1) %>%
add_feature_contiguity_constraints()
# print problem
print(p19)## A conservation problem (<ConservationProblem>)
## ├•data
## │├•features: "feature_1", "feature_2", "feature_3", "feature_4", and "feature_5" (5 total)
## │└•planning units:
## │ ├•data: <SpatRaster> (90 total)
## │ ├•costs: continuous values (between 190.1328 and 215.8638)
## │ ├•extent: 0, 0, 1, 1 (xmin, ymin, xmax, ymax)
## │ └•CRS: Undefined Cartesian SRS (projected)
## ├•formulation
## │├•objective: minimum set objective
## │├•penalties: none specified
## │├•targets: relative targets (between 0.1 and 0.1)
## │├•constraints:
## ││└•1: feature contiguity constraints
## │└•decisions: binary decision
## └•optimization
## ├•portfolio: default portfolio
## └•solver: gurobi solver (`gap` = 0.1, `time_limit` = 2147483647, `first_feasible` = FALSE, …)
## # ℹ Use `summary(...)` to see complete formulation.- Linear constraints: Add constraints to ensure that all selected planning units meet certain criteria. For example, they can be used to add multiple budgets, or limit the number of planning units selected in different administrative areas within a study region (e.g., different countries).
# create problem with constraints which specify that the sum of
# values in sim_features[[1]] among selected planning units must not exceed a
# threshold value of 190.
p20 <-
problem(sim_pu_raster, sim_features) %>%
add_min_shortfall_objective(budget = 1800) %>%
add_relative_targets(0.1) %>%
add_linear_constraints(190, "<=", sim_features[[1]])
# print problem
print(p20)## A conservation problem (<ConservationProblem>)
## ├•data
## │├•features: "feature_1", "feature_2", "feature_3", "feature_4", and "feature_5" (5 total)
## │└•planning units:
## │ ├•data: <SpatRaster> (90 total)
## │ ├•costs: continuous values (between 190.1328 and 215.8638)
## │ ├•extent: 0, 0, 1, 1 (xmin, ymin, xmax, ymax)
## │ └•CRS: Undefined Cartesian SRS (projected)
## ├•formulation
## │├•objective: minimum shortfall objective (`budget` = 1800)
## │├•penalties: none specified
## │├•targets: relative targets (between 0.1 and 0.1)
## │├•constraints:
## ││└•1: linear constraints (`threshold` = 190, `sense` = "<=", …)
## │└•decisions: binary decision
## └•optimization
## ├•portfolio: default portfolio
## └•solver: gurobi solver (`gap` = 0.1, `time_limit` = 2147483647, `first_feasible` = FALSE, …)
## # ℹ Use `summary(...)` to see complete formulation.- Mandatory allocation constraints: Add constraints to ensure that every planning unit is allocated to a management zone in the solution. Please note that this function can only be used with problems that contain multiple zones. For more information on problems with multiple zones and an example using this function, see the Management zones tutorial.
In particular, The add_locked_in_constraints and
add_locked_out_constraints functions are incredibly useful
for real-world conservation planning exercises, so it’s worth pointing
out that there are several ways we can specify which planning units
should be locked in or out of the solutions. If we use raster planning
unit data, we can also use raster data to specify which planning units
should be locked in or locked out.
# load data to lock in or lock out planning units
sim_locked_in_raster <- get_sim_locked_in_raster()
sim_locked_out_raster <- get_sim_locked_out_raster()
# plot the locked data
plot(
c(sim_locked_in_raster, sim_locked_out_raster),
xlim = c(-0.1, 1.1), ylim = c(-0.1, 1.1),
axes = FALSE, main = c("Locked in", "Locked out")
)
# create a problem using raster planning unit data and use the locked raster
# data to lock in some planning units and lock out some other planning units
p21 <-
problem(sim_pu_raster, sim_features) %>%
add_min_set_objective() %>%
add_relative_targets(0.1) %>%
add_locked_in_constraints(sim_locked_in_raster) %>%
add_locked_out_constraints(sim_locked_out_raster)
# print problem
print(p21)## A conservation problem (<ConservationProblem>)
## ├•data
## │├•features: "feature_1", "feature_2", "feature_3", "feature_4", and "feature_5" (5 total)
## │└•planning units:
## │ ├•data: <SpatRaster> (90 total)
## │ ├•costs: continuous values (between 190.1328 and 215.8638)
## │ ├•extent: 0, 0, 1, 1 (xmin, ymin, xmax, ymax)
## │ └•CRS: Undefined Cartesian SRS (projected)
## ├•formulation
## │├•objective: minimum set objective
## │├•penalties: none specified
## │├•targets: relative targets (between 0.1 and 0.1)
## │├•constraints:
## ││├•1: locked in constraints (10 planning units)
## ││└•2: locked out constraints (10 planning units)
## │└•decisions: binary decision
## └•optimization
## ├•portfolio: default portfolio
## └•solver: gurobi solver (`gap` = 0.1, `time_limit` = 2147483647, `first_feasible` = FALSE, …)
## # ℹ Use `summary(...)` to see complete formulation.If our planning unit data are in a spatial vector format (similar to
the sim_pu_polygons data) or a tabular format (similar to
pu_dat), we can use the field names in the data to refer to
which planning units should be locked in and / or out. For example, the
sim_pu_polygons object has TRUE /
FALSE values in the “locked_in” field which indicate which
planning units should be selected in the solution. We could use the data
in this field to specify that those planning units with
TRUE values should be locked in using the following
methods.
# preview sim_pu_polygons object
print(sim_pu_polygons)## Simple feature collection with 90 features and 3 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 0 ymin: 0 xmax: 1 ymax: 1
## Projected CRS: Undefined Cartesian SRS
## # A tibble: 90 × 4
## cost locked_in locked_out geom
## * <dbl> <lgl> <lgl> <POLYGON [m]>
## 1 216. FALSE FALSE ((0 1, 0.1 1, 0.1 0.9, 0 0.9, 0 1))
## 2 213. FALSE FALSE ((0.1 1, 0.2 1, 0.2 0.9, 0.1 0.9, 0.1 1))
## 3 207. FALSE FALSE ((0.2 1, 0.3 1, 0.3 0.9, 0.2 0.9, 0.2 1))
## 4 209. FALSE TRUE ((0.3 1, 0.4 1, 0.4 0.9, 0.3 0.9, 0.3 1))
## 5 214. FALSE FALSE ((0.4 1, 0.5 1, 0.5 0.9, 0.4 0.9, 0.4 1))
## 6 214. FALSE FALSE ((0.5 1, 0.6 1, 0.6 0.9, 0.5 0.9, 0.5 1))
## 7 210. FALSE FALSE ((0.6 1, 0.7 1, 0.7 0.9, 0.6 0.9, 0.6 1))
## 8 211. FALSE TRUE ((0.7 1, 0.8 1, 0.8 0.9, 0.7 0.9, 0.7 1))
## 9 210. FALSE FALSE ((0.8 1, 0.9 1, 0.9 0.9, 0.8 0.9, 0.8 1))
## 10 204. FALSE FALSE ((0.9 1, 1 1, 1 0.9, 0.9 0.9, 0.9 1))
## # ℹ 80 more rows
# specify locked in data using the field name
p22 <-
problem(sim_pu_polygons, sim_features, cost_column = "cost") %>%
add_min_set_objective() %>%
add_relative_targets(0.1) %>%
add_locked_in_constraints("locked_in")
# print problem
print(p22)## A conservation problem (<ConservationProblem>)
## ├•data
## │├•features: "feature_1", "feature_2", "feature_3", "feature_4", and "feature_5" (5 total)
## │└•planning units:
## │ ├•data: <sftbl_dftbldata.frame> (90 total)
## │ ├•costs: continuous values (between 190.1328 and 215.8638)
## │ ├•extent: 0, 0, 1, 1 (xmin, ymin, xmax, ymax)
## │ └•CRS: Undefined Cartesian SRS (projected)
## ├•formulation
## │├•objective: minimum set objective
## │├•penalties: none specified
## │├•targets: relative targets (between 0.1 and 0.1)
## │├•constraints:
## ││└•1: locked in constraints (10 planning units)
## │└•decisions: binary decision
## └•optimization
## ├•portfolio: default portfolio
## └•solver: gurobi solver (`gap` = 0.1, `time_limit` = 2147483647, `first_feasible` = FALSE, …)
## # ℹ Use `summary(...)` to see complete formulation.
# specify locked in data using the values in the field
p23 <-
problem(sim_pu_polygons, sim_features, cost_column = "cost") %>%
add_min_set_objective() %>%
add_relative_targets(0.1) %>%
add_locked_in_constraints(which(sim_pu_polygons$locked_in))
# print problem
print(p23)## A conservation problem (<ConservationProblem>)
## ├•data
## │├•features: "feature_1", "feature_2", "feature_3", "feature_4", and "feature_5" (5 total)
## │└•planning units:
## │ ├•data: <sftbl_dftbldata.frame> (90 total)
## │ ├•costs: continuous values (between 190.1328 and 215.8638)
## │ ├•extent: 0, 0, 1, 1 (xmin, ymin, xmax, ymax)
## │ └•CRS: Undefined Cartesian SRS (projected)
## ├•formulation
## │├•objective: minimum set objective
## │├•penalties: none specified
## │├•targets: relative targets (between 0.1 and 0.1)
## │├•constraints:
## ││└•1: locked in constraints (10 planning units)
## │└•decisions: binary decision
## └•optimization
## ├•portfolio: default portfolio
## └•solver: gurobi solver (`gap` = 0.1, `time_limit` = 2147483647, `first_feasible` = FALSE, …)
## # ℹ Use `summary(...)` to see complete formulation.Add penalties
We can also add penalties to a problem to favor or penalize solutions
according to a secondary objective. Unlike the constraint functions,
these functions will add extra information to the objective function of
the optimization function to penalize solutions that do not exhibit
specific characteristics. For example, penalties can be added to a
problem to avoid highly fragmented solutions at the expense of accepting
slightly more expensive solutions. All penalty functions have a
penalty argument that controls the relative importance of
the secondary penalty function compared to the primary objective
function. It is worth noting that incredibly low or incredibly high
penalty values – relative to the main objective function –
can cause problems to take a very long time to solve, so when trying out
a range of different penalty values it can be helpful to limit the
solver to run for a set period of time. The following penalties are
available.
- Boundary penalties: Add penalties to penalize solutions that are excessively fragmented. These penalties are similar to those used in Marxan (Ball et al. 2009; Beyer et al. 2016).
# create problem with penalties that penalize fragmented solutions with a
# penalty factor of 0.01
p24 <-
problem(sim_pu_raster, sim_features) %>%
add_min_set_objective() %>%
add_relative_targets(0.1) %>%
add_boundary_penalties(penalty = 0.01)
# print problem
print(p24)## A conservation problem (<ConservationProblem>)
## ├•data
## │├•features: "feature_1", "feature_2", "feature_3", "feature_4", and "feature_5" (5 total)
## │└•planning units:
## │ ├•data: <SpatRaster> (90 total)
## │ ├•costs: continuous values (between 190.1328 and 215.8638)
## │ ├•extent: 0, 0, 1, 1 (xmin, ymin, xmax, ymax)
## │ └•CRS: Undefined Cartesian SRS (projected)
## ├•formulation
## │├•objective: minimum set objective
## │├•penalties:
## ││└•1: boundary penalties (`penalty` = 0.01, `edge_factor` = 0.5, …)
## │├•targets: relative targets (between 0.1 and 0.1)
## │├•constraints: none specified
## │└•decisions: binary decision
## └•optimization
## ├•portfolio: default portfolio
## └•solver: gurobi solver (`gap` = 0.1, `time_limit` = 2147483647, `first_feasible` = FALSE, …)
## # ℹ Use `summary(...)` to see complete formulation.- Connectivity penalties: Add penalties to favor solutions that select combinations of planning units with high connectivity between them. This function is used for symmetric connectivities among planning units (based on Önal & Briers 2002).
# create problem with penalties for symmetric connectivity,
# here we will use only the first four layers in
# sim_features for the features and we will use the fifth layer in sim_features
# to represent the connectivity data, where the connectivity_matrix function
# will create a matrix showing the average strength of connectivity between
# adjacent planning units using the data in the fifth layer of sim_features
p25 <-
problem(sim_pu_raster, sim_features[[1:4]]) %>%
add_min_set_objective() %>%
add_relative_targets(0.1) %>%
add_connectivity_penalties(
penalty = 5,
data = connectivity_matrix(sim_pu_raster, sim_features[[5]])
)
# print problem
print(p25)## A conservation problem (<ConservationProblem>)
## ├•data
## │├•features: "feature_1", "feature_2", "feature_3", and "feature_4" (4 total)
## │└•planning units:
## │ ├•data: <SpatRaster> (90 total)
## │ ├•costs: continuous values (between 190.1328 and 215.8638)
## │ ├•extent: 0, 0, 1, 1 (xmin, ymin, xmax, ymax)
## │ └•CRS: Undefined Cartesian SRS (projected)
## ├•formulation
## │├•objective: minimum set objective
## │├•penalties:
## ││└•1: connectivity penalties (`penalty` = 5, …)
## │├•targets: relative targets (between 0.1 and 0.1)
## │├•constraints: none specified
## │└•decisions: binary decision
## └•optimization
## ├•portfolio: default portfolio
## └•solver: gurobi solver (`gap` = 0.1, `time_limit` = 2147483647, `first_feasible` = FALSE, …)
## # ℹ Use `summary(...)` to see complete formulation.- Asymmetric connectivity penalties: Add penalties to penalize solutions that select planning units that share high connectivity values with other planning units that are not selected by the solution. This function is used for asymmetric connectivities among planning units (based on Beger et al. 2010).
# create a matrix containing asymmetric connectivity data,
asym_con_matrix <- matrix(
runif(ncell(sim_pu_raster)^2),
nrow = ncell(sim_pu_raster),
ncol = ncell(sim_pu_raster)
)
# create problem with penalties for asymmetric connectivity
p26 <-
problem(sim_pu_raster, sim_features[[1:4]]) %>%
add_min_set_objective() %>%
add_relative_targets(0.1) %>%
add_asym_connectivity_penalties(penalty = 5, data = asym_con_matrix)
# print problem
print(p26)## A conservation problem (<ConservationProblem>)
## ├•data
## │├•features: "feature_1", "feature_2", "feature_3", and "feature_4" (4 total)
## │└•planning units:
## │ ├•data: <SpatRaster> (90 total)
## │ ├•costs: continuous values (between 190.1328 and 215.8638)
## │ ├•extent: 0, 0, 1, 1 (xmin, ymin, xmax, ymax)
## │ └•CRS: Undefined Cartesian SRS (projected)
## ├•formulation
## │├•objective: minimum set objective
## │├•penalties:
## ││└•1: asymmetric connectivity penalties (`penalty` = 5, …)
## │├•targets: relative targets (between 0.1 and 0.1)
## │├•constraints: none specified
## │└•decisions: binary decision
## └•optimization
## ├•portfolio: default portfolio
## └•solver: gurobi solver (`gap` = 0.1, `time_limit` = 2147483647, `first_feasible` = FALSE, …)
## # ℹ Use `summary(...)` to see complete formulation.- Linear penalties: Add penalties to penalize solutions that select planning units according to a certain variable (e.g., anthropogenic pressure).
# create data for penalizing planning units
pen_raster <- simulate_cost(sim_pu_raster)
# create problem with penalties that penalize solutions that select
# planning units with high values in the pen_raster object,
# here we will use a penalty value of 5 to indicate the trade-off (scaling)
# between the penalty values (in the sim_pu_raster) and the main objective
# (i.e., the cost of the solution)
p27 <-
problem(sim_pu_raster, sim_features) %>%
add_min_set_objective() %>%
add_relative_targets(0.1) %>%
add_linear_penalties(penalty = 5, data = pen_raster)
# print problem
print(p27)## A conservation problem (<ConservationProblem>)
## ├•data
## │├•features: "feature_1", "feature_2", "feature_3", "feature_4", and "feature_5" (5 total)
## │└•planning units:
## │ ├•data: <SpatRaster> (90 total)
## │ ├•costs: continuous values (between 190.1328 and 215.8638)
## │ ├•extent: 0, 0, 1, 1 (xmin, ymin, xmax, ymax)
## │ └•CRS: Undefined Cartesian SRS (projected)
## ├•formulation
## │├•objective: minimum set objective
## │├•penalties:
## ││└•1: linear penalties (`penalty` = 5, …)
## │├•targets: relative targets (between 0.1 and 0.1)
## │├•constraints: none specified
## │└•decisions: binary decision
## └•optimization
## ├•portfolio: default portfolio
## └•solver: gurobi solver (`gap` = 0.1, `time_limit` = 2147483647, `first_feasible` = FALSE, …)
## # ℹ Use `summary(...)` to see complete formulation.Add the decision types
Conservation planning problems involve making decisions on how planning units will be managed. These decisions are then associated with management actions (e.g., turning a planning unit into a protected area). The type of decision describes how the action is applied to planning units. For instance, the default decision-type is a binary decision type, meaning that we are either selecting or not selecting planning units for management. The following decision types are available.
- Binary decisions: Add a binary decision to a conservation planning problem. This is the classic decision of either prioritizing or not prioritizing a planning unit. Typically, this decision has the assumed action of buying the planning unit to include in a protected area network. If no decision is added to a problem object, then this decision class will be used by default.
# add binary decisions to a problem
p28 <-
problem(sim_pu_raster, sim_features) %>%
add_min_set_objective() %>%
add_relative_targets(0.1) %>%
add_binary_decisions()
# print problem
print(p28)## A conservation problem (<ConservationProblem>)
## ├•data
## │├•features: "feature_1", "feature_2", "feature_3", "feature_4", and "feature_5" (5 total)
## │└•planning units:
## │ ├•data: <SpatRaster> (90 total)
## │ ├•costs: continuous values (between 190.1328 and 215.8638)
## │ ├•extent: 0, 0, 1, 1 (xmin, ymin, xmax, ymax)
## │ └•CRS: Undefined Cartesian SRS (projected)
## ├•formulation
## │├•objective: minimum set objective
## │├•penalties: none specified
## │├•targets: relative targets (between 0.1 and 0.1)
## │├•constraints: none specified
## │└•decisions: binary decision
## └•optimization
## ├•portfolio: default portfolio
## └•solver: gurobi solver (`gap` = 0.1, `time_limit` = 2147483647, `first_feasible` = FALSE, …)
## # ℹ Use `summary(...)` to see complete formulation.- Proportion decisions: Add a proportion decision to a problem. This is a relaxed decision where a part of a planning unit can be prioritized, as opposed to the default of the entire planning unit. Typically, this decision has the assumed action of buying a fraction of a planning unit to include in a protected area network. Generally, problems can be solved much faster with proportion-type decisions than binary-type decisions, so they can be very useful when commercial solvers are not available.
# add proportion decisions to a problem
p29 <-
problem(sim_pu_raster, sim_features) %>%
add_min_set_objective() %>%
add_relative_targets(0.1) %>%
add_proportion_decisions()
# print problem
print(p29)## A conservation problem (<ConservationProblem>)
## ├•data
## │├•features: "feature_1", "feature_2", "feature_3", "feature_4", and "feature_5" (5 total)
## │└•planning units:
## │ ├•data: <SpatRaster> (90 total)
## │ ├•costs: continuous values (between 190.1328 and 215.8638)
## │ ├•extent: 0, 0, 1, 1 (xmin, ymin, xmax, ymax)
## │ └•CRS: Undefined Cartesian SRS (projected)
## ├•formulation
## │├•objective: minimum set objective
## │├•penalties: none specified
## │├•targets: relative targets (between 0.1 and 0.1)
## │├•constraints: none specified
## │└•decisions: proportion decision
## └•optimization
## ├•portfolio: default portfolio
## └•solver: gurobi solver (`gap` = 0.1, `time_limit` = 2147483647, `first_feasible` = FALSE, …)
## # ℹ Use `summary(...)` to see complete formulation.- Semi-continuous decisions: Add a semi-continuous decision to a problem. This decision is similar to proportion decisions except that it has an upper bound parameter. By default, the decision can range from prioritizing none (0%) to all (100%) of a planning unit. However, a upper bound can be specified to ensure that at most only a fraction (e.g., 80%) of a planning unit can be purchased. This type of decision may be useful when it is not practical to conserve the entire area indicated by a planning unit.
# add semi-continuous decisions to a problem, where we can only manage at most
# 50 % of the area encompassed by a planning unit
p30 <-
problem(sim_pu_raster, sim_features) %>%
add_min_set_objective() %>%
add_relative_targets(0.1) %>%
add_semicontinuous_decisions(0.5)
# print problem
print(p30)## A conservation problem (<ConservationProblem>)
## ├•data
## │├•features: "feature_1", "feature_2", "feature_3", "feature_4", and "feature_5" (5 total)
## │└•planning units:
## │ ├•data: <SpatRaster> (90 total)
## │ ├•costs: continuous values (between 190.1328 and 215.8638)
## │ ├•extent: 0, 0, 1, 1 (xmin, ymin, xmax, ymax)
## │ └•CRS: Undefined Cartesian SRS (projected)
## ├•formulation
## │├•objective: minimum set objective
## │├•penalties: none specified
## │├•targets: relative targets (between 0.1 and 0.1)
## │├•constraints: none specified
## │└•decisions: semicontinuous decision (`upper_limit` = 0.5)
## └•optimization
## ├•portfolio: default portfolio
## └•solver: gurobi solver (`gap` = 0.1, `time_limit` = 2147483647, `first_feasible` = FALSE, …)
## # ℹ Use `summary(...)` to see complete formulation.Add a solver
Next, after specifying the mathematical formulation that underpins your conservation planning problem, you can specify how the problem should be solved. If you do not specify this information, the prioritizr R package will automatically use the best solver currently installed on your system with some reasonable defaults. We strongly recommend installing the Gurobi software suite and the gurobi R package to solve problems, and for more information on this topic please refer to the Gurobi Installation Guide. The following solvers are available.
- Gurobi solver: Gurobi is a state of the art commercial optimization software. It is by far the fastest of the solvers that can be used to solve conservation problems. However, it is not freely available. That said, special licenses are available to academics at no cost.
# create a problem and specify that Gurobi should be used to solve the problem
# and specify an optimality gap of zero to obtain the optimal solution
p31 <-
problem(sim_pu_raster, sim_features) %>%
add_min_set_objective() %>%
add_relative_targets(0.1) %>%
add_binary_decisions() %>%
add_gurobi_solver(gap = 0)
# print problem
print(p31)## A conservation problem (<ConservationProblem>)
## ├•data
## │├•features: "feature_1", "feature_2", "feature_3", "feature_4", and "feature_5" (5 total)
## │└•planning units:
## │ ├•data: <SpatRaster> (90 total)
## │ ├•costs: continuous values (between 190.1328 and 215.8638)
## │ ├•extent: 0, 0, 1, 1 (xmin, ymin, xmax, ymax)
## │ └•CRS: Undefined Cartesian SRS (projected)
## ├•formulation
## │├•objective: minimum set objective
## │├•penalties: none specified
## │├•targets: relative targets (between 0.1 and 0.1)
## │├•constraints: none specified
## │└•decisions: binary decision
## └•optimization
## ├•portfolio: default portfolio
## └•solver: gurobi solver (`gap` = 0, `time_limit` = 2147483647, `first_feasible` = FALSE, …)
## # ℹ Use `summary(...)` to see complete formulation.- IBM CPLEX solver: IBM CPLEX is a commercial optimization software. It is faster than the open source solvers available for generating prioritizations (see below), however, it is not freely available. Similar to the Gurobi software, special licenses are available to academics at no cost.
# create a problem and specify that IBM CPLEX should be used to solve the
# problem and specify an optimality gap of zero to obtain the optimal solution
p32 <-
problem(sim_pu_raster, sim_features) %>%
add_min_set_objective() %>%
add_relative_targets(0.1) %>%
add_binary_decisions() %>%
add_cplex_solver(gap = 0)
# print problem
print(p32)## A conservation problem (<ConservationProblem>)
## ├•data
## │├•features: "feature_1", "feature_2", "feature_3", "feature_4", and "feature_5" (5 total)
## │└•planning units:
## │ ├•data: <SpatRaster> (90 total)
## │ ├•costs: continuous values (between 190.1328 and 215.8638)
## │ ├•extent: 0, 0, 1, 1 (xmin, ymin, xmax, ymax)
## │ └•CRS: Undefined Cartesian SRS (projected)
## ├•formulation
## │├•objective: minimum set objective
## │├•penalties: none specified
## │├•targets: relative targets (between 0.1 and 0.1)
## │├•constraints: none specified
## │└•decisions: binary decision
## └•optimization
## ├•portfolio: default portfolio
## └•solver: cplex solver (`gap` = 0, `time_limit` = 2147483647, …)
## # ℹ Use `summary(...)` to see complete formulation.-
CBC solver: CBC is an open-source
mixed integer programming solver that is part of the Computational
Infrastructure for Operations Research (COIN-OR) project. This solver
seems to have much better performance than the other open-source solvers
(i.e.,
add_rsymphony_solver(),add_lpsymphony_solver()) (see the Solver benchmarks vignette for details). As such, it is strongly recommended to use this solver if the Gurobi and IBM CPLEX solvers are not available. This solver requires the rcbc R package, which is currently only available on GitHub.
# create a problem and specify that CBC should be used to solve the
# problem and specify an optimality gap of zero to obtain the optimal solution
p33 <-
problem(sim_pu_raster, sim_features) %>%
add_min_set_objective() %>%
add_relative_targets(0.1) %>%
add_binary_decisions() %>%
add_cbc_solver(gap = 0)
# print problem
print(p33)## A conservation problem (<ConservationProblem>)
## ├•data
## │├•features: "feature_1", "feature_2", "feature_3", "feature_4", and "feature_5" (5 total)
## │└•planning units:
## │ ├•data: <SpatRaster> (90 total)
## │ ├•costs: continuous values (between 190.1328 and 215.8638)
## │ ├•extent: 0, 0, 1, 1 (xmin, ymin, xmax, ymax)
## │ └•CRS: Undefined Cartesian SRS (projected)
## ├•formulation
## │├•objective: minimum set objective
## │├•penalties: none specified
## │├•targets: relative targets (between 0.1 and 0.1)
## │├•constraints: none specified
## │└•decisions: binary decision
## └•optimization
## ├•portfolio: default portfolio
## └•solver: cbc solver (`gap` = 0, `time_limit` = 2147483647, `first_feasible` = FALSE, …)
## # ℹ Use `summary(...)` to see complete formulation.-
HiGHS solver: HiGHS is an open source
optimization software. Although this solver can have comparable
performance to the CBC solver (i.e.,
add_cbc_solver()) for particular problems and is generally faster than the SYMPHONY based solvers (i.e.,add_rsymphony_solver(),add_lpsymphony_solver()), it can sometimes take much longer than the CBC solver for particular problems. This solver is recommended if theadd_gurobi_solver(),add_cplex_solver(),add_cbc_solver()cannot be used.
# create a problem and specify that HiGHS should be used to solve the
# problem and specify an optimality gap of zero to obtain the optimal solution
p34 <-
problem(sim_pu_raster, sim_features) %>%
add_min_set_objective() %>%
add_relative_targets(0.1) %>%
add_binary_decisions() %>%
add_highs_solver(gap = 0)
# print problem
print(p34)## A conservation problem (<ConservationProblem>)
## ├•data
## │├•features: "feature_1", "feature_2", "feature_3", "feature_4", and "feature_5" (5 total)
## │└•planning units:
## │ ├•data: <SpatRaster> (90 total)
## │ ├•costs: continuous values (between 190.1328 and 215.8638)
## │ ├•extent: 0, 0, 1, 1 (xmin, ymin, xmax, ymax)
## │ └•CRS: Undefined Cartesian SRS (projected)
## ├•formulation
## │├•objective: minimum set objective
## │├•penalties: none specified
## │├•targets: relative targets (between 0.1 and 0.1)
## │├•constraints: none specified
## │└•decisions: binary decision
## └•optimization
## ├•portfolio: default portfolio
## └•solver: highs solver (`gap` = 0, `time_limit` = 2147483647, …)
## # ℹ Use `summary(...)` to see complete formulation.- lpsymphony solver: SYMPHONY is an open-source integer programming solver that is also part of the COIN-OR project. This solver uses the lpsymphony R package (available on Bioconductor) to interface with the SYMPHONY software.
# create a problem and specify that lpsymphony should be used to solve the
# problem and specify an optimality gap of zero to obtain the optimal solution
p35 <-
problem(sim_pu_raster, sim_features) %>%
add_min_set_objective() %>%
add_relative_targets(0.1) %>%
add_binary_decisions() %>%
add_lpsymphony_solver(gap = 0)
# print problem
print(p35)## A conservation problem (<ConservationProblem>)
## ├•data
## │├•features: "feature_1", "feature_2", "feature_3", "feature_4", and "feature_5" (5 total)
## │└•planning units:
## │ ├•data: <SpatRaster> (90 total)
## │ ├•costs: continuous values (between 190.1328 and 215.8638)
## │ ├•extent: 0, 0, 1, 1 (xmin, ymin, xmax, ymax)
## │ └•CRS: Undefined Cartesian SRS (projected)
## ├•formulation
## │├•objective: minimum set objective
## │├•penalties: none specified
## │├•targets: relative targets (between 0.1 and 0.1)
## │├•constraints: none specified
## │└•decisions: binary decision
## └•optimization
## ├•portfolio: default portfolio
## └•solver: lpsymphony (`gap` = 0, `time_limit` = 2147483647, `first_feasible` = FALSE, …)
## # ℹ Use `summary(...)` to see complete formulation.- Rsymphony solver: This solver provides a different interface to the SYMPHONY software. It uses the Rsymphony R package which is available on The Comprehensive R Archive Network (CRAN). This solver is generally slower than the other solvers, because it cannot use parallel processing.
# create a problem and specify that Rsymphony should be used to solve the
# problem and specify an optimality gap of zero to obtain the optimal solution
p36 <-
problem(sim_pu_raster, sim_features) %>%
add_min_set_objective() %>%
add_relative_targets(0.1) %>%
add_binary_decisions() %>%
add_rsymphony_solver(gap = 0)
# print problem
print(p36)## A conservation problem (<ConservationProblem>)
## ├•data
## │├•features: "feature_1", "feature_2", "feature_3", "feature_4", and "feature_5" (5 total)
## │└•planning units:
## │ ├•data: <SpatRaster> (90 total)
## │ ├•costs: continuous values (between 190.1328 and 215.8638)
## │ ├•extent: 0, 0, 1, 1 (xmin, ymin, xmax, ymax)
## │ └•CRS: Undefined Cartesian SRS (projected)
## ├•formulation
## │├•objective: minimum set objective
## │├•penalties: none specified
## │├•targets: relative targets (between 0.1 and 0.1)
## │├•constraints: none specified
## │└•decisions: binary decision
## └•optimization
## ├•portfolio: default portfolio
## └•solver: rsymphony solver (`gap` = 0, `time_limit` = 2147483647, `first_feasible` = FALSE, …)
## # ℹ Use `summary(...)` to see complete formulation.Add a portfolio
Many conservation planning exercises require a portfolio of solutions. For example, real-world exercises can involve presenting decision makers with a range of near-optimal decisions. Additionally, the number of times that different planning units are selected in different solutions can provide insight into their relative importance. The following portfolio methods are available.
- Gap portfolio: Generate a portfolio of solutions by finding a certain number of solutions that are all within a pre-specified optimality gap. This method is especially useful for generating multiple solutions that can used be to calculate selection frequencies (similar to Marxan). Note that this method requires that the Gurobi optimization software is used to generate solutions.
# create a problem and specify that a portfolio should be created by
# finding five solutions within 10% of optimality
p37 <-
problem(sim_pu_raster, sim_features) %>%
add_min_set_objective() %>%
add_relative_targets(0.1) %>%
add_binary_decisions() %>%
add_gap_portfolio(number_solutions = 5, pool_gap = 0.2)
# print problem
print(p37)## A conservation problem (<ConservationProblem>)
## ├•data
## │├•features: "feature_1", "feature_2", "feature_3", "feature_4", and "feature_5" (5 total)
## │└•planning units:
## │ ├•data: <SpatRaster> (90 total)
## │ ├•costs: continuous values (between 190.1328 and 215.8638)
## │ ├•extent: 0, 0, 1, 1 (xmin, ymin, xmax, ymax)
## │ └•CRS: Undefined Cartesian SRS (projected)
## ├•formulation
## │├•objective: minimum set objective
## │├•penalties: none specified
## │├•targets: relative targets (between 0.1 and 0.1)
## │├•constraints: none specified
## │└•decisions: binary decision
## └•optimization
## ├•portfolio: gap portfolio (`number_solutions` = 5, `pool_gap` = 0.2)
## └•solver: gurobi solver (`gap` = 0.1, `time_limit` = 2147483647, `first_feasible` = FALSE, …)
## # ℹ Use `summary(...)` to see complete formulation.- Top portfolio: Generate a portfolio of solutions by finding a pre-specified number of solutions that are closest to optimality (i.e the top solutions). Note that this method requires that the Gurobi optimization software is used to generate solutions.
# create a problem and specify that a portfolio should be created using
# the top five solutions
p38 <-
problem(sim_pu_raster, sim_features) %>%
add_min_set_objective() %>%
add_relative_targets(0.1) %>%
add_binary_decisions() %>%
add_top_portfolio(number_solutions = 5)
# print problem
print(p38)## A conservation problem (<ConservationProblem>)
## ├•data
## │├•features: "feature_1", "feature_2", "feature_3", "feature_4", and "feature_5" (5 total)
## │└•planning units:
## │ ├•data: <SpatRaster> (90 total)
## │ ├•costs: continuous values (between 190.1328 and 215.8638)
## │ ├•extent: 0, 0, 1, 1 (xmin, ymin, xmax, ymax)
## │ └•CRS: Undefined Cartesian SRS (projected)
## ├•formulation
## │├•objective: minimum set objective
## │├•penalties: none specified
## │├•targets: relative targets (between 0.1 and 0.1)
## │├•constraints: none specified
## │└•decisions: binary decision
## └•optimization
## ├•portfolio: top portfolio (`number_solutions` = 5)
## └•solver: gurobi solver (`gap` = 0.1, `time_limit` = 2147483647, `first_feasible` = FALSE, …)
## # ℹ Use `summary(...)` to see complete formulation.- Extra portfolio: Generate a portfolio of solutions by storing feasible solutions found during the optimization process. Note that this method requires that the Gurobi optimization software is used to generate solutions.
# create a problem and specify that a portfolio should be created using
# extra solutions found while solving the problem
p39 <-
problem(sim_pu_raster, sim_features) %>%
add_min_set_objective() %>%
add_relative_targets(0.1) %>%
add_binary_decisions() %>%
add_extra_portfolio()
# print problem
print(p39)## A conservation problem (<ConservationProblem>)
## ├•data
## │├•features: "feature_1", "feature_2", "feature_3", "feature_4", and "feature_5" (5 total)
## │└•planning units:
## │ ├•data: <SpatRaster> (90 total)
## │ ├•costs: continuous values (between 190.1328 and 215.8638)
## │ ├•extent: 0, 0, 1, 1 (xmin, ymin, xmax, ymax)
## │ └•CRS: Undefined Cartesian SRS (projected)
## ├•formulation
## │├•objective: minimum set objective
## │├•penalties: none specified
## │├•targets: relative targets (between 0.1 and 0.1)
## │├•constraints: none specified
## │└•decisions: binary decision
## └•optimization
## ├•portfolio: extra portfolio
## └•solver: gurobi solver (`gap` = 0.1, `time_limit` = 2147483647, `first_feasible` = FALSE, …)
## # ℹ Use `summary(...)` to see complete formulation.- Cuts portfolio: Generate a portfolio of distinct solutions within a pre-specified optimality gap. This method is only recommended if the Gurobi optimization solver is not available.
# create a problem and specify that a portfolio containing 10 solutions
# should be created using using Bender's cuts
p40 <-
problem(sim_pu_raster, sim_features) %>%
add_min_set_objective() %>%
add_relative_targets(0.1) %>%
add_binary_decisions() %>%
add_cuts_portfolio(number_solutions = 10)
# print problem
print(p40)## A conservation problem (<ConservationProblem>)
## ├•data
## │├•features: "feature_1", "feature_2", "feature_3", "feature_4", and "feature_5" (5 total)
## │└•planning units:
## │ ├•data: <SpatRaster> (90 total)
## │ ├•costs: continuous values (between 190.1328 and 215.8638)
## │ ├•extent: 0, 0, 1, 1 (xmin, ymin, xmax, ymax)
## │ └•CRS: Undefined Cartesian SRS (projected)
## ├•formulation
## │├•objective: minimum set objective
## │├•penalties: none specified
## │├•targets: relative targets (between 0.1 and 0.1)
## │├•constraints: none specified
## │└•decisions: binary decision
## └•optimization
## ├•portfolio: cuts portfolio (`number_solutions` = 10)
## └•solver: gurobi solver (`gap` = 0.1, `time_limit` = 2147483647, `first_feasible` = FALSE, …)
## # ℹ Use `summary(...)` to see complete formulation.- Shuffle portfolio: Generate a portfolio of solutions by randomly reordering the data prior to attempting to solve the problem. If the Gurobi optimization solver is not available, this method is the fastest method for generating a set number of solutions within a specified distance from optimality.
# create a problem and specify a portfolio should be created that contains
# 10 solutions and that any duplicate solutions should not be removed
p41 <-
problem(sim_pu_raster, sim_features) %>%
add_min_set_objective() %>%
add_relative_targets(0.1) %>%
add_binary_decisions() %>%
add_shuffle_portfolio(number_solutions = 10, remove_duplicates = FALSE)
# print problem
print(p41)## A conservation problem (<ConservationProblem>)
## ├•data
## │├•features: "feature_1", "feature_2", "feature_3", "feature_4", and "feature_5" (5 total)
## │└•planning units:
## │ ├•data: <SpatRaster> (90 total)
## │ ├•costs: continuous values (between 190.1328 and 215.8638)
## │ ├•extent: 0, 0, 1, 1 (xmin, ymin, xmax, ymax)
## │ └•CRS: Undefined Cartesian SRS (projected)
## ├•formulation
## │├•objective: minimum set objective
## │├•penalties: none specified
## │├•targets: relative targets (between 0.1 and 0.1)
## │├•constraints: none specified
## │└•decisions: binary decision
## └•optimization
## ├•portfolio: shuffle portfolio (`number_solutions` = 10, …)
## └•solver: gurobi solver (`gap` = 0.1, `time_limit` = 2147483647, `first_feasible` = FALSE, …)
## # ℹ Use `summary(...)` to see complete formulation.Solve the problem
After formulating our conservation planning problem and specifying
how the problem should be solved, we can use the solve()
function to obtain a solution. Note that the solver will typically print
out some information describing the size of the problem and report its
progress when searching for a suitable solution.
# formulate the problem
p42 <-
problem(sim_pu_raster, sim_features) %>%
add_min_set_objective() %>%
add_relative_targets(0.1) %>%
add_boundary_penalties(penalty = 500, edge_factor = 0.5) %>%
add_binary_decisions()
# print problem
print(p42)## A conservation problem (<ConservationProblem>)
## ├•data
## │├•features: "feature_1", "feature_2", "feature_3", "feature_4", and "feature_5" (5 total)
## │└•planning units:
## │ ├•data: <SpatRaster> (90 total)
## │ ├•costs: continuous values (between 190.1328 and 215.8638)
## │ ├•extent: 0, 0, 1, 1 (xmin, ymin, xmax, ymax)
## │ └•CRS: Undefined Cartesian SRS (projected)
## ├•formulation
## │├•objective: minimum set objective
## │├•penalties:
## ││└•1: boundary penalties (`penalty` = 500, `edge_factor` = 0.5, …)
## │├•targets: relative targets (between 0.1 and 0.1)
## │├•constraints: none specified
## │└•decisions: binary decision
## └•optimization
## ├•portfolio: default portfolio
## └•solver: gurobi solver (`gap` = 0.1, `time_limit` = 2147483647, `first_feasible` = FALSE, …)
## # ℹ Use `summary(...)` to see complete formulation.
# solve the problem (using the default solver)
s42 <- solve(p42)## Set parameter Username
## Set parameter TimeLimit to value 2147483647
## Set parameter MIPGap to value 0.1
## Set parameter Presolve to value 2
## Set parameter Threads to value 1
## Academic license - for non-commercial use only - expires 2025-04-21
## Gurobi Optimizer version 11.0.3 build v11.0.3rc0 (linux64 - "Ubuntu 22.04.4 LTS")
##
## CPU model: 11th Gen Intel(R) Core(TM) i7-1185G7 @ 3.00GHz, instruction set [SSE2|AVX|AVX2|AVX512]
## Thread count: 4 physical cores, 8 logical processors, using up to 1 threads
##
## Optimize a model with 293 rows, 234 columns and 1026 nonzeros
## Model fingerprint: 0xb4b90756
## Variable types: 0 continuous, 234 integer (234 binary)
## Coefficient statistics:
## Matrix range [2e-01, 1e+00]
## Objective range [1e+02, 4e+02]
## Bounds range [1e+00, 1e+00]
## RHS range [3e+00, 8e+00]
## Found heuristic solution: objective 20287.197006
## Found heuristic solution: objective 3087.9617767
## Presolve time: 0.00s
## Presolved: 293 rows, 234 columns, 1026 nonzeros
## Variable types: 0 continuous, 234 integer (234 binary)
## Root relaxation presolved: 293 rows, 234 columns, 1026 nonzeros
##
##
## Root relaxation: objective 2.265862e+03, 228 iterations, 0.00 seconds (0.00 work units)
##
## Nodes | Current Node | Objective Bounds | Work
## Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
##
## 0 0 2265.86242 0 234 3087.96178 2265.86242 26.6% - 0s
## 0 0 2328.27423 0 232 3087.96178 2328.27423 24.6% - 0s
## 0 0 2357.91894 0 229 3087.96178 2357.91894 23.6% - 0s
## H 0 0 3034.0204010 2357.91894 22.3% - 0s
## 0 0 2379.43275 0 197 3034.02040 2379.43275 21.6% - 0s
## 0 0 2379.43275 0 197 3034.02040 2379.43275 21.6% - 0s
## H 0 0 2890.8754120 2379.43275 17.7% - 0s
## H 0 0 2814.4777374 2379.43275 15.5% - 0s
## 0 2 2379.78730 0 197 2814.47774 2379.78730 15.4% - 0s
## H 20 20 2648.2596588 2382.48313 10.0% 23.4 0s
##
## Cutting planes:
## Gomory: 3
##
## Explored 34 nodes (1117 simplex iterations) in 0.09 seconds (0.12 work units)
## Thread count was 1 (of 8 available processors)
##
## Solution count 6: 2648.26 2814.48 2890.88 ... 20287.2
##
## Optimal solution found (tolerance 1.00e-01)
## Best objective 2.648259658813e+03, best bound 2.391338223393e+03, gap 9.7015%
# plot solution
plot(
s42, col = c("grey90", "darkgreen"), main = "Solution",
xlim = c(-0.1, 1.1), ylim = c(-0.1, 1.1), axes = FALSE
)
We can plot this solution because the planning unit input data are
spatially referenced in a raster format. The output format will always
match the planning unit data used to initialize the problem. For
example, the solution to a problem with planning units in a spatial
vector (shapefile) format would also be in a spatial vector format.
Similarly, if the planning units were in a tabular format (i.e.,
data.frame), the solution would also be returned in a
tabular format.
We can also extract attributes from the solution that describe the quality of the solution and the optimization process.
# extract the objective (numerical value being minimized or maximized)
print(attr(s42, "objective"))## solution_1
## 2648.26## solution_1
## 0.111
# extract state message from the solver that describes why this specific
# solution was returned
print(attr(s42, "status"))## solution_1
## "OPTIMAL"Evaluate the solution
Conservation planning involves making trade-offs between different criteria (e.g., overall cost, feature representation, connectivity). After obtaining a solution to a conservation planning problem, it is important to evaluate it to help understand the trade-offs made by the prioritization. This is also useful to compare different solutions with each other.
Evaluating performance
Summary statistics can be computed to evaluate the overall performance of a solution based on certain criteria. The following summaries can be computed.
The following functions are available to summarize a solution:
- Number summary: Calculate the number of planning units selected by a solution.
# calculate statistic
eval_n_summary(p42, s42)## # A tibble: 1 × 2
## summary n
## <chr> <dbl>
## 1 overall 10- Cost summary: Calculate the total cost of a solution.
# calculate statistic
eval_cost_summary(p42, s42)## # A tibble: 1 × 2
## summary cost
## <chr> <dbl>
## 1 overall 2048.- Feature representation summary: Calculate how well features are represented by a solution. This function can be used for any type of problem.
# calculate statistics
eval_feature_representation_summary(p42, s42)## # A tibble: 5 × 5
## summary feature total_amount absolute_held relative_held
## <chr> <chr> <dbl> <dbl> <dbl>
## 1 overall feature_1 83.3 8.57 0.103
## 2 overall feature_2 31.2 3.18 0.102
## 3 overall feature_3 72.0 7.28 0.101
## 4 overall feature_4 42.7 5.01 0.118
## 5 overall feature_5 56.7 5.71 0.101- Target coverage summary: Calculate how well representation targets are met by a solution. This function can only be used with problems contain targets.
# calculate statistics
eval_target_coverage_summary(p42, s42)## # A tibble: 5 × 9
## feature met total_amount absolute_target absolute_held absolute_shortfall
## <chr> <lgl> <dbl> <dbl> <dbl> <dbl>
## 1 feature_1 TRUE 83.3 8.33 8.57 0
## 2 feature_2 TRUE 31.2 3.12 3.18 0
## 3 feature_3 TRUE 72.0 7.20 7.28 0
## 4 feature_4 TRUE 42.7 4.27 5.01 0
## 5 feature_5 TRUE 56.7 5.67 5.71 0
## # ℹ 3 more variables: relative_target <dbl>, relative_held <dbl>,
## # relative_shortfall <dbl>- Boundary summary: Calculate the total exposed boundary length (perimeter) associated with a solution.
# calculate statistic
eval_boundary_summary(p42, s42)## # A tibble: 1 × 2
## summary boundary
## <chr> <dbl>
## 1 overall 1.2- Connectivity summary: Calculate the connectivity of a solution using symmetric connectivity data.
# create symmetric connectivity data
# here we parametrize connectivity based on which planning units are
# adjacent to each other
cm <- adjacency_matrix(sim_pu_raster)
# calculate statistic
eval_connectivity_summary(p42, s42, data = cm)## # A tibble: 1 × 2
## summary connectivity
## <chr> <dbl>
## 1 overall 24- Asymmetric connectivity summary: Calculate the connectivity of a solution using asymmetric connectivity data.
# create asymmetric connectivity data
# here we parametrize connectivity by randomly simulating values
acm <- matrix(runif(ncell(sim_pu_raster) ^ 2), ncol = ncell(sim_pu_raster))
# calculate statistic
eval_asym_connectivity_summary(p42, s42, data = acm)## # A tibble: 1 × 2
## summary asym_connectivity
## <chr> <dbl>
## 1 overall 50.1Evaluating relative importance
Conservation plans can take a long time to implement. Since funding availability and habitat quality can decline over time, it is critical that the most important places in a prioritization are scheduled for management as early as possible. For instance, some planning units in a solution might contain many rare species which do not occur in any other planning units. Alternatively, some planning units might offer an especially high return on investment that reduces costs considerably. As a consequence, conservation planners often need information on which planning units selected in a prioritization are most important to the overall success of the prioritization. To achieve this, conservation planners can use importance scores for each planning unit selected by a solution.
Let’s generate a prioritization so that can compare the different importance methods.
# formulate the problem
p43 <-
problem(sim_pu_raster, sim_features) %>%
add_min_set_objective() %>%
add_relative_targets(0.1) %>%
add_binary_decisions()
# solve the problem
s43 <- solve(p43)## Set parameter Username
## Set parameter TimeLimit to value 2147483647
## Set parameter MIPGap to value 0.1
## Set parameter Presolve to value 2
## Set parameter Threads to value 1
## Academic license - for non-commercial use only - expires 2025-04-21
## Gurobi Optimizer version 11.0.3 build v11.0.3rc0 (linux64 - "Ubuntu 22.04.4 LTS")
##
## CPU model: 11th Gen Intel(R) Core(TM) i7-1185G7 @ 3.00GHz, instruction set [SSE2|AVX|AVX2|AVX512]
## Thread count: 4 physical cores, 8 logical processors, using up to 1 threads
##
## Optimize a model with 5 rows, 90 columns and 450 nonzeros
## Model fingerprint: 0x4bb5d283
## Variable types: 0 continuous, 90 integer (90 binary)
## Coefficient statistics:
## Matrix range [2e-01, 9e-01]
## Objective range [2e+02, 2e+02]
## Bounds range [1e+00, 1e+00]
## RHS range [3e+00, 8e+00]
## Found heuristic solution: objective 2337.9617767
## Presolve time: 0.00s
## Presolved: 5 rows, 90 columns, 450 nonzeros
## Variable types: 0 continuous, 90 integer (90 binary)
## Found heuristic solution: objective 2332.1004028
## Root relaxation presolved: 5 rows, 90 columns, 450 nonzeros
##
##
## Root relaxation: objective 1.931582e+03, 12 iterations, 0.00 seconds (0.00 work units)
##
## Nodes | Current Node | Objective Bounds | Work
## Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
##
## 0 0 1931.58191 0 4 2332.10040 1931.58191 17.2% - 0s
## H 0 0 1987.3985291 1931.58191 2.81% - 0s
##
## Explored 1 nodes (12 simplex iterations) in 0.00 seconds (0.00 work units)
## Thread count was 1 (of 8 available processors)
##
## Solution count 3: 1987.4 2332.1 2337.96
##
## Optimal solution found (tolerance 1.00e-01)
## Best objective 1.987398529053e+03, best bound 1.931581907658e+03, gap 2.8085%
# plot solution
plot(
s43, col = c("grey90", "darkgreen"), main = "Solution",
xlim = c(-0.1, 1.1), ylim = c(-0.1, 1.1), axes = FALSE
)
The following methods are available for computing importance scores.
- Replacement cost: Evaluate importance using the replacement cost method (Cabeza & Moilanen 2006). This method quantifies the importance of a given planning unit as the decrease in the performance of the solution (based on the objective function) if the planning unit cannot be acquired (e.g., in terms of the additional costs required to meet feature targets). The advantages of this method are that it (i) accounts for the costs of the planning units, (ii) can be applied to multiple management zones, (iii) can be applied to any objective function, and (iv) can identify truly irreplaceable planning units (denoted with infinite values). However, the key disadvantage of this method, is that can take an infeasible amount of time to complete for large and complex problems.
# calculate replacement cost scores and make the solver quiet
rc43 <-
p43 %>%
add_default_solver(gap = 0, verbose = FALSE) %>%
eval_replacement_importance(s43)
# plot replacement cost scores
plot(
rc43, main = "Replacement cost scores",
xlim = c(-0.1, 1.1), ylim = c(-0.1, 1.1), axes = FALSE
)
* Ferrier method: Evaluate importance by computing
irreplaceability scores following Ferrier et al. (2000). The advantages of this method are that
it (i) can be computed relatively quickly for moderate and large-sized
problems, and (ii) calculates a score for each feature within each
planning unit to provide insight into why certain planning units are
more important than others. The disadvantage with this method is that it
can only be applied to conservation problems that use targets and have a
single zone (i.e., similar to Marxan-type problems).
# calculate Ferrier scores and extract total score
fs43 <- eval_ferrier_importance(p43, s43)[["total"]]
# plot Ferrier scores
plot(
fs43, main = "Ferrier scores",
xlim = c(-0.1, 1.1), ylim = c(-0.1, 1.1), axes = FALSE
)
* Rarity weighted richness: Evaluate importance by
computing rarity weighted richness scores (Williams et al. 1996). The only
advantage with this method is that it can be computed very quickly for
very large problems. The key disadvantage with this approach is that it
merely describes the spatial patterns of biodiversity, and does not
consider any of the goals that underpin conservation planning exercise.
For instance, it does not account for planning costs, management zones,
objective functions, or feature representation targets.
# calculate rarity weighted richness scores
rwr43 <- eval_rare_richness_importance(p43, s43)
# plot rarity weighted richness scores
plot(
rwr43, main = "Rarity weighted richness scores",
xlim = c(-0.1, 1.1), ylim = c(-0.1, 1.1), axes = FALSE
)
In general, we recommend using replacement cost scores for small and moderate sized problems (e.g., less than 30,000 planning units) when it is feasible to do so. It can take a very long time to compute replacement cost scores, and so it is simply not feasible to compute these scores for particularly large problems. For moderate and large sized problems (e.g., more than 30,000 planning units), we recommend using the Ferrier method. This is because it explicitly accounts for representation targets, unlike the rarity weighted richness scores. We almost never recommend using the rarity weighted richness scores. This is because they do not consider criteria needed to inform conservation decision making (Brown et al. 2015).
Marxan problems
Although we encourage users to build and tailor conservation planning
problems to suit their own needs, sometimes it just simply easier to use
something you’re already familiar with. The
marxan_problem() function is provided as a convenient
wrapper for building and solving Marxan-style conservation
problems. If users already have their conservation planning data
formatted for use with Marxan, this function can also read
Marxan data files and solve the Marxan-style problems
using exact algorithm solvers. Please note that problems built
using the marxan_problem() function are still solved the
same way as a problem initialized using the problem()
function, and therefore still require the installation of one of the
solver packages.
Here is a short example showing how the marxan_problem()
function can be used to read Marxan input files and the
solve function can be used to solve the problem.
# set file path for Marxan input file
minput <- system.file("extdata/marxan/input.dat", package = "prioritizr")
# read Marxan input file
mp <- marxan_problem(minput)
# print problem
print(mp)## A conservation problem (<ConservationProblem>)
## ├•data
## │├•features: "bird1", "nvis2", "nvis8", "nvis9", "nvis14", "nvis20", … (17 total)
## │└•planning units:
## │ ├•data: <data.frame> (1751 total)
## │ ├•costs: continuous values (between 0 and 415692.1938)
## │ ├•extent: NA
## │ └•CRS: NA
## ├•formulation
## │├•objective: minimum set objective
## │├•penalties:
## ││└•1: boundary penalties (`penalty` = 1, `edge_factor` = 1, …)
## │├•targets: relative targets (between 0.3 and 0.3)
## │├•constraints:
## ││├•1: locked in constraints (317 planning units)
## ││└•2: locked out constraints (1 planning units)
## │└•decisions: binary decision
## └•optimization
## ├•portfolio: default portfolio
## └•solver: gurobi solver (`gap` = 0.1, `time_limit` = 2147483647, `first_feasible` = FALSE, …)
## # ℹ Use `summary(...)` to see complete formulation.
# solve the problem
ms <- solve(mp)## Set parameter Username
## Set parameter TimeLimit to value 2147483647
## Set parameter MIPGap to value 0.1
## Set parameter Presolve to value 2
## Set parameter Threads to value 1
## Academic license - for non-commercial use only - expires 2025-04-21
## Gurobi Optimizer version 11.0.3 build v11.0.3rc0 (linux64 - "Ubuntu 22.04.4 LTS")
##
## CPU model: 11th Gen Intel(R) Core(TM) i7-1185G7 @ 3.00GHz, instruction set [SSE2|AVX|AVX2|AVX512]
## Thread count: 4 physical cores, 8 logical processors, using up to 1 threads
##
## Optimize a model with 10075 rows, 6780 columns and 24778 nonzeros
## Model fingerprint: 0x7cea1a5d
## Variable types: 0 continuous, 6780 integer (6780 binary)
## Coefficient statistics:
## Matrix range [5e-05, 4e+03]
## Objective range [8e+03, 4e+05]
## Bounds range [1e+00, 1e+00]
## RHS range [5e+03, 3e+05]
## Found heuristic solution: objective 1.255825e+08
## Presolve removed 4707 rows and 3103 columns
## Presolve time: 0.05s
## Presolved: 5368 rows, 3677 columns, 12704 nonzeros
## Variable types: 0 continuous, 3677 integer (3677 binary)
## Found heuristic solution: objective 1.065946e+08
## Root relaxation presolved: 5368 rows, 3677 columns, 12704 nonzeros
##
##
## Root relaxation: interrupted, 11 iterations, 0.01 seconds (0.01 work units)
##
## Nodes | Current Node | Objective Bounds | Work
## Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
##
## 0 0 - 0 1.0659e+08 9.5941e+07 10.0% - 0s
##
## Explored 1 nodes (11 simplex iterations) in 0.06 seconds (0.09 work units)
## Thread count was 1 (of 8 available processors)
##
## Solution count 2: 1.06595e+08 1.25582e+08
##
## Optimal solution found (tolerance 1.00e-01)
## Best objective 1.065946314895e+08, best bound 9.594065554616e+07, gap 9.9949%
# since the Marxan data was in a tabular format, the solution is also returned
# in a tabular format, so we will print the first six rows of the table
# containing the solution
head(ms)## id cost status xloc yloc locked_in locked_out solution_1
## 1 3 0.000 0 1116623 -4493479 FALSE FALSE 0
## 2 30 7527.275 3 1110623 -4496943 FALSE TRUE 0
## 3 56 37349.075 0 1092623 -4500408 FALSE FALSE 0
## 4 58 16959.021 0 1116623 -4500408 FALSE FALSE 0
## 5 84 34220.256 0 1098623 -4503872 FALSE FALSE 0
## 6 85 178907.584 0 1110623 -4503872 FALSE FALSE 0Alternatively, rather then using a Marxan input file to
construct the problem, we can manually read in the Marxan data
files and input these to the marxan_problem() function.
# load data
pu <-
system.file("extdata/marxan/input/pu.dat", package = "prioritizr") %>%
read.table(sep = ",", header = TRUE)
features <-
system.file("extdata/marxan/input/spec.dat", package = "prioritizr") %>%
read.table(sep = ",", header = TRUE)
bound <-
system.file("extdata/marxan/input/bound.dat", package = "prioritizr") %>%
read.table(sep = "\t", header = TRUE)
rij <-
system.file("extdata/marxan/input/puvspr.dat", package = "prioritizr") %>%
read.table(sep = ",", header = TRUE)
# build Marxan problem using data.frame objects
mp2 <- marxan_problem(
x = pu, spec = features, puvspr = rij, bound = bound, blm = 0
)
# print problem
print(mp2)## A conservation problem (<ConservationProblem>)
## ├•data
## │├•features: "bird1", "nvis2", "nvis8", "nvis9", "nvis14", "nvis20", … (17 total)
## │└•planning units:
## │ ├•data: <data.frame> (1751 total)
## │ ├•costs: continuous values (between 0 and 415692.1938)
## │ ├•extent: NA
## │ └•CRS: NA
## ├•formulation
## │├•objective: minimum set objective
## │├•penalties:
## ││└•1: boundary penalties (`penalty` = 0, `edge_factor` = 1, …)
## │├•targets: relative targets (between 0.3 and 0.3)
## │├•constraints:
## ││├•1: locked in constraints (317 planning units)
## ││└•2: locked out constraints (1 planning units)
## │└•decisions: binary decision
## └•optimization
## ├•portfolio: default portfolio
## └•solver: gurobi solver (`gap` = 0.1, `time_limit` = 2147483647, `first_feasible` = FALSE, …)
## # ℹ Use `summary(...)` to see complete formulation.Conservation planning problems that are built using the
marxan_problem() function can also be customized. For
example, we could change the decision type for mp2 to
involve selecting a proportion of each planing unit (using the
add_proportion_decisions() function).
Conclusion
Hopefully, this vignette has provided an informative overview of the prioritizr R package. For more examples using the package, please see the other vignettes. Perhaps, one of the best ways to learn a new piece of software is to just try it out. Test it, try breaking it, make mistakes, and learn from them. We would recommend trying to build conservation planning problems that resemble those you face in your own work—but using the built-in example data sets. This way you can quickly verify that the problems you build actually mean what you think they mean. For instance, you can try playing around with the targets and see what effect they have on the solutions, or try playing around with penalties and see what effect they have on the solutions. Finally, if you have any questions about using the package or suggestions for it, please post an issue on the package’s online coding repository.