Title:  Systematic Conservation Prioritization in R 

Description:  Systematic conservation prioritization using mixed integer linear programming (MILP). It provides a flexible interface for building and solving conservation planning problems. Once built, conservation planning problems can be solved using a variety of commercial and opensource exact algorithm solvers. By using exact algorithm solvers, solutions can be generated that are guaranteed to be optimal (or within a prespecified optimality gap). Furthermore, conservation problems can be constructed to optimize the spatial allocation of different management actions or zones, meaning that conservation practitioners can identify solutions that benefit multiple stakeholders. To solve largescale or complex conservation planning problems, users should install the Gurobi optimization software (available from <https://www.gurobi.com/>) and the 'gurobi' R package (see Gurobi Installation Guide vignette for details). Users can also install the IBM CPLEX software (<https://www.ibm.com/products/ilogcplexoptimizationstudio/cplexoptimizer>) and the 'cplexAPI' R package (available at <https://github.com/cran/cplexAPI>). Additionally, the 'rcbc' R package (available at <https://github.com/dirkschumacher/rcbc>) can be used to generate solutions using the CBC optimization software (<https://github.com/coinor/Cbc>). For further details, see Hanson et al. (2024) <doi:10.1111/cobi.14376>. 
Authors:  Jeffrey O Hanson [aut] , Richard Schuster [aut, cre] , Nina Morrell [aut], Matthew StrimasMackey [aut] , Brandon P M Edwards [aut] , Matthew E Watts [aut], Peter Arcese [aut] , Joseph R Bennett [aut] , Hugh P Possingham [aut] 
Maintainer:  Richard Schuster <[email protected]> 
License:  GPL3 
Version:  8.0.4.4 
Built:  20240930 23:35:47 UTC 
Source:  https://github.com/prioritizr/prioritizr 
Set targets expressed as the actual value of features in the study area that need to be represented in the prioritization. For instance, setting a target of 10 requires that the solution secure a set of planning units for which their summed feature values are equal to or greater than 10.
add_absolute_targets(x, targets) ## S4 method for signature 'ConservationProblem,numeric' add_absolute_targets(x, targets) ## S4 method for signature 'ConservationProblem,matrix' add_absolute_targets(x, targets) ## S4 method for signature 'ConservationProblem,character' add_absolute_targets(x, targets)
add_absolute_targets(x, targets) ## S4 method for signature 'ConservationProblem,numeric' add_absolute_targets(x, targets) ## S4 method for signature 'ConservationProblem,matrix' add_absolute_targets(x, targets) ## S4 method for signature 'ConservationProblem,character' add_absolute_targets(x, targets)
x 

targets 
object that specifies the targets for each feature. See the Targets format section for more information. 
Targets are used to specify the minimum amount or proportion of a
feature's distribution that needs to be protected. Most conservation
planning problems require targets with the exception of the maximum cover
(see add_max_cover_objective()
) and maximum utility
(see add_max_utility_objective()
) problems. Attempting to solve
problems with objectives that require targets without specifying targets
will throw an error.
For problems associated with multiple management zones,
add_absolute_targets()
can
be used to set targets that each pertain to a single feature and a single
zone. To set targets that can be met through allocating different
planning units to multiple zones, see the add_manual_targets()
function. An example of a target that could be met through allocations
to multiple zones might be where each management zone is expected to
result in a different amount of a feature and the target requires that
the total amount of the feature in all zones must exceed a certain
threshold. In other words, the target does not require that any single
zone secure a specific amount of the feature, but the total amount held
in all zones must secure a specific amount. Thus the target could,
potentially, be met through allocating all planning units to any specific
management zone, or through allocating the planning units to different
combinations of management zones.
An updated problem()
object with the targets added to it.
The targets
for a problem can be specified using the following formats.
targets
as a numeric
vectorcontaining target values for each feature. Additionally, for convenience, this format can be a single value to assign the same target to each feature. Note that this format cannot be used to specify targets for problems with multiple zones.
targets
as a matrix
objectcontaining a target for each feature
in each zone.
Here, each row corresponds to a different feature in argument to
x
, each column corresponds to a different zone in argument to
x
, and each cell contains the target value for a given feature
that the solution needs to secure in a given zone.
targets
as a character
vectorcontaining the column name(s) in the
feature data associated with the argument to x
that
contain targets. This format can only be used when the
feature data associated with x
is a sf::st_sf()
or data.frame
.
For problems that contain a single zone, the argument to targets
must
contain a single column name. Otherwise, for problems that
contain multiple zones, the argument to targets
must
contain a column name for each zone.
See targets for an overview of all functions for adding targets.
Other targets:
add_loglinear_targets()
,
add_manual_targets()
,
add_relative_targets()
## Not run: # set seed for reproducibility set.seed(500) # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_features < get_sim_zones_features() # create minimal problem with no targets p0 < problem(sim_pu_raster, sim_features) %>% add_min_set_objective() %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # create problem with targets to secure 3 amounts for each feature p1 < p0 %>% add_absolute_targets(3) # create problem with varying targets for each feature targets < c(1, 2, 3, 2, 1) p2 < p0 %>% add_absolute_targets(targets) # solve problem s1 < c(solve(p1), solve(p2)) names(s1) < c("equal targets", "varying targets") # plot solution plot(s1, axes = FALSE) # create a problem with multiple management zones p3 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_min_set_objective() %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # create a problem with targets that specify an equal amount of each feature # to be represented in each zone p4_targets < matrix( 2, nrow = number_of_features(sim_zones_features), ncol = number_of_zones(sim_zones_features), dimnames = list( feature_names(sim_zones_features), zone_names(sim_zones_features) ) ) print(p4_targets) p4 < p3 %>% add_absolute_targets(p4_targets) # solve problem s4 < solve(p4) # plot solution (pixel values correspond to zone identifiers) plot(category_layer(s4), main = "equal targets", axes = FALSE) # create a problem with targets that require a varying amount of each # feature to be represented in each zone p5_targets < matrix( rpois(15, 1), nrow = number_of_features(sim_zones_features), ncol = number_of_zones(sim_zones_features), dimnames = list( feature_names(sim_zones_features), zone_names(sim_zones_features) ) ) print(p5_targets) p5 < p3 %>% add_absolute_targets(p5_targets) # solve problem s5 < solve(p5) # plot solution (pixel values correspond to zone identifiers) plot(category_layer(s5), main = "varying targets", axes = FALSE) ## End(Not run)
## Not run: # set seed for reproducibility set.seed(500) # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_features < get_sim_zones_features() # create minimal problem with no targets p0 < problem(sim_pu_raster, sim_features) %>% add_min_set_objective() %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # create problem with targets to secure 3 amounts for each feature p1 < p0 %>% add_absolute_targets(3) # create problem with varying targets for each feature targets < c(1, 2, 3, 2, 1) p2 < p0 %>% add_absolute_targets(targets) # solve problem s1 < c(solve(p1), solve(p2)) names(s1) < c("equal targets", "varying targets") # plot solution plot(s1, axes = FALSE) # create a problem with multiple management zones p3 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_min_set_objective() %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # create a problem with targets that specify an equal amount of each feature # to be represented in each zone p4_targets < matrix( 2, nrow = number_of_features(sim_zones_features), ncol = number_of_zones(sim_zones_features), dimnames = list( feature_names(sim_zones_features), zone_names(sim_zones_features) ) ) print(p4_targets) p4 < p3 %>% add_absolute_targets(p4_targets) # solve problem s4 < solve(p4) # plot solution (pixel values correspond to zone identifiers) plot(category_layer(s4), main = "equal targets", axes = FALSE) # create a problem with targets that require a varying amount of each # feature to be represented in each zone p5_targets < matrix( rpois(15, 1), nrow = number_of_features(sim_zones_features), ncol = number_of_zones(sim_zones_features), dimnames = list( feature_names(sim_zones_features), zone_names(sim_zones_features) ) ) print(p5_targets) p5 < p3 %>% add_absolute_targets(p5_targets) # solve problem s5 < solve(p5) # plot solution (pixel values correspond to zone identifiers) plot(category_layer(s5), main = "varying targets", axes = FALSE) ## End(Not run)
Add penalties to a conservation planning problem to account for asymmetric connectivity between planning units. Asymmetric connectivity data describe connectivity information that is directional. For example, asymmetric connectivity data could describe the strength of rivers flowing between different planning units. Since river flow is directional, the level of connectivity from an upstream planning unit to a downstream planning unit would be higher than that from a downstream planning unit to an upstream planning unit.
## S4 method for signature 'ConservationProblem,ANY,ANY,matrix' add_asym_connectivity_penalties(x, penalty, zones, data) ## S4 method for signature 'ConservationProblem,ANY,ANY,Matrix' add_asym_connectivity_penalties(x, penalty, zones, data) ## S4 method for signature 'ConservationProblem,ANY,ANY,data.frame' add_asym_connectivity_penalties(x, penalty, zones, data) ## S4 method for signature 'ConservationProblem,ANY,ANY,dgCMatrix' add_asym_connectivity_penalties(x, penalty, zones, data) ## S4 method for signature 'ConservationProblem,ANY,ANY,array' add_asym_connectivity_penalties(x, penalty, zones, data)
## S4 method for signature 'ConservationProblem,ANY,ANY,matrix' add_asym_connectivity_penalties(x, penalty, zones, data) ## S4 method for signature 'ConservationProblem,ANY,ANY,Matrix' add_asym_connectivity_penalties(x, penalty, zones, data) ## S4 method for signature 'ConservationProblem,ANY,ANY,data.frame' add_asym_connectivity_penalties(x, penalty, zones, data) ## S4 method for signature 'ConservationProblem,ANY,ANY,dgCMatrix' add_asym_connectivity_penalties(x, penalty, zones, data) ## S4 method for signature 'ConservationProblem,ANY,ANY,array' add_asym_connectivity_penalties(x, penalty, zones, data)
x 

penalty 

zones 

data 

This function adds penalties to conservation planning problem to penalize solutions that have low connectivity. Specifically, it penalizes solutions that select planning units that share high connectivity values with other planning units that are not selected by the solution (based on Beger et al. 2010).
An updated problem()
object with the penalties added to it.
The connectivity penalties are implemented using the following equations.
Let $I$
represent the set of planning units
(indexed by $i$
or $j$
), $Z$
represent the set
of management zones (indexed by $z$
or $y$
), and $X_{iz}$
represent the decision variable for planning unit $i$
for in zone
$z$
(e.g., with binary
values one indicating if planning unit is allocated or not). Also, let
$p$
represent the argument to penalty
, $D$
represent the
argument to data
, and $W$
represent the argument
to zones
.
If the argument to data
is supplied as a matrix
or
Matrix
object, then the penalties are calculated as:
$\sum_{i}^{I} \sum_{j}^{I} \sum_{z}^{Z} \sum_{y}^{Z}
(p \times X_{iz} \times D_{ij} \times W_{zy}) 
\sum_{i}^{I} \sum_{j}^{I} \sum_{z}^{Z} \sum_{y}^{Z}
(p \times X_{iz} \times X_{jy} \times D_{ij} \times W_{zy})$
Otherwise, if the argument to data
is supplied as an
array
object, then the penalties are
calculated as:
$\sum_{i}^{I} \sum_{j}^{I} \sum_{z}^{Z} \sum_{y}^{Z}
(p \times X_{iz} \times D_{ijzy}) 
\sum_{i}^{I} \sum_{j}^{I} \sum_{z}^{Z} \sum_{y}^{Z}
(p \times X_{iz} \times X_{jy} \times D_{ijzy})$
Note that when the problem objective is to maximize some measure of
benefit and not minimize some measure of cost, the term $p$
is
replaced with $p$
.
The argument to data
can be specified using several different formats.
data
as a matrix
/Matrix
objectwhere rows and columns represent
different planning units and the value of each cell represents the
strength of connectivity between two different planning units. Cells
that occur along the matrix diagonal are treated as weights which
indicate that planning units are more desirable in the solution.
The argument to zones
can be used to control
the strength of connectivity between planning units in different zones.
The default argument for zones
is to treat planning units
allocated to different zones as having zero connectivity.
data
as a data.frame
objectcontaining columns that are named
"id1"
, "id2"
, and "boundary"
. Here, each row
denotes the connectivity between a pair of planning units
(per values in the "id1"
and "id2"
columns) following the
Marxan format.
If the argument to x
contains multiple zones, then the
"zone1"
and "zone2"
columns can optionally be provided to manually
specify the connectivity values between planning units when they are
allocated to specific zones. If the "zone1"
and
"zone2"
columns are present, then the argument to zones
must be
NULL
.
data
as an array
objectcontaining fourdimensions where cell values
indicate the strength of connectivity between planning units
when they are assigned to specific management zones. The first two
dimensions (i.e., rows and columns) indicate the strength of
connectivity between different planning units and the second two
dimensions indicate the different management zones. Thus
the data[1, 2, 3, 4]
indicates the strength of
connectivity between planning unit 1 and planning unit 2 when planning
unit 1 is assigned to zone 3 and planning unit 2 is assigned to zone 4.
Beger M, Linke S, Watts M, Game E, Treml E, Ball I, and Possingham, HP (2010) Incorporating asymmetric connectivity into spatial decision making for conservation, Conservation Letters, 3: 359–368.
See penalties for an overview of all functions for adding penalties.
Other penalties:
add_boundary_penalties()
,
add_connectivity_penalties()
,
add_feature_weights()
,
add_linear_penalties()
## Not run: # set seed for reproducibility set.seed(600) # load data sim_pu_polygons < get_sim_pu_polygons() sim_features < get_sim_features() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_features < get_sim_zones_features() # create basic problem p1 < problem(sim_pu_polygons, sim_features, "cost") %>% add_min_set_objective() %>% add_relative_targets(0.2) %>% add_default_solver(verbose = FALSE) # create an asymmetric connectivity matrix. Here, connectivity occurs between # adjacent planning units and, due to rivers flowing southwards # through the study area, connectivity from northern planning units to # southern planning units is ten times stronger than the reverse. acm1 < matrix(0, nrow(sim_pu_polygons), nrow(sim_pu_polygons)) acm1 < as(acm1, "Matrix") centroids < sf::st_coordinates( suppressWarnings(sf::st_centroid(sim_pu_polygons)) ) adjacent_units < sf::st_intersects(sim_pu_polygons, sparse = FALSE) for (i in seq_len(nrow(sim_pu_polygons))) { for (j in seq_len(nrow(sim_pu_polygons))) { # find if planning units are adjacent if (adjacent_units[i, j]) { # find if planning units lay north and south of each other # i.e., they have the same xcoordinate if (centroids[i, 1] == centroids[j, 1]) { if (centroids[i, 2] > centroids[j, 2]) { # if i is north of j add 10 units of connectivity acm1[i, j] < acm1[i, j] + 10 } else if (centroids[i, 2] < centroids[j, 2]) { # if i is south of j add 1 unit of connectivity acm1[i, j] < acm1[i, j] + 1 } } } } } # rescale matrix values to have a maximum value of 1 acm1 < rescale_matrix(acm1, max = 1) # visualize asymmetric connectivity matrix Matrix::image(acm1) # create penalties penalties < c(1, 50) # create problems using the different penalties p2 < list( p1, p1 %>% add_asym_connectivity_penalties(penalties[1], data = acm1), p1 %>% add_asym_connectivity_penalties(penalties[2], data = acm1) ) # solve problems s2 < lapply(p2, solve) # create object with all solutions s2 < sf::st_sf( tibble::tibble( p2_1 = s2[[1]]$solution_1, p2_2 = s2[[2]]$solution_1, p2_3 = s2[[3]]$solution_1 ), geometry = sf::st_geometry(s2[[1]]) ) names(s2)[1:3] < c("basic problem", paste0("acm1 (", penalties,")")) # plot solutions based on different penalty values plot(s2, cex = 1.5) # create minimal multizone problem and limit solver to one minute # to obtain solutions in a short period of time p3 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_min_set_objective() %>% add_relative_targets(matrix(0.15, nrow = 5, ncol = 3)) %>% add_binary_decisions() %>% add_default_solver(time_limit = 60, verbose = FALSE) # crate asymmetric connectivity data by randomly simulating values acm2 < matrix( runif(ncell(sim_zones_pu_raster) ^ 2), nrow = ncell(sim_zones_pu_raster) ) # create multizone problems using the penalties p4 < list( p3, p3 %>% add_asym_connectivity_penalties(penalties[1], data = acm2), p3 %>% add_asym_connectivity_penalties(penalties[2], data = acm2) ) # solve problems s4 < lapply(p4, solve) s4 < lapply(s4, category_layer) s4 < terra::rast(s4) names(s4) < c("basic problem", paste0("acm2 (", penalties,")")) # plot solutions plot(s4, axes = FALSE) ## End(Not run)
## Not run: # set seed for reproducibility set.seed(600) # load data sim_pu_polygons < get_sim_pu_polygons() sim_features < get_sim_features() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_features < get_sim_zones_features() # create basic problem p1 < problem(sim_pu_polygons, sim_features, "cost") %>% add_min_set_objective() %>% add_relative_targets(0.2) %>% add_default_solver(verbose = FALSE) # create an asymmetric connectivity matrix. Here, connectivity occurs between # adjacent planning units and, due to rivers flowing southwards # through the study area, connectivity from northern planning units to # southern planning units is ten times stronger than the reverse. acm1 < matrix(0, nrow(sim_pu_polygons), nrow(sim_pu_polygons)) acm1 < as(acm1, "Matrix") centroids < sf::st_coordinates( suppressWarnings(sf::st_centroid(sim_pu_polygons)) ) adjacent_units < sf::st_intersects(sim_pu_polygons, sparse = FALSE) for (i in seq_len(nrow(sim_pu_polygons))) { for (j in seq_len(nrow(sim_pu_polygons))) { # find if planning units are adjacent if (adjacent_units[i, j]) { # find if planning units lay north and south of each other # i.e., they have the same xcoordinate if (centroids[i, 1] == centroids[j, 1]) { if (centroids[i, 2] > centroids[j, 2]) { # if i is north of j add 10 units of connectivity acm1[i, j] < acm1[i, j] + 10 } else if (centroids[i, 2] < centroids[j, 2]) { # if i is south of j add 1 unit of connectivity acm1[i, j] < acm1[i, j] + 1 } } } } } # rescale matrix values to have a maximum value of 1 acm1 < rescale_matrix(acm1, max = 1) # visualize asymmetric connectivity matrix Matrix::image(acm1) # create penalties penalties < c(1, 50) # create problems using the different penalties p2 < list( p1, p1 %>% add_asym_connectivity_penalties(penalties[1], data = acm1), p1 %>% add_asym_connectivity_penalties(penalties[2], data = acm1) ) # solve problems s2 < lapply(p2, solve) # create object with all solutions s2 < sf::st_sf( tibble::tibble( p2_1 = s2[[1]]$solution_1, p2_2 = s2[[2]]$solution_1, p2_3 = s2[[3]]$solution_1 ), geometry = sf::st_geometry(s2[[1]]) ) names(s2)[1:3] < c("basic problem", paste0("acm1 (", penalties,")")) # plot solutions based on different penalty values plot(s2, cex = 1.5) # create minimal multizone problem and limit solver to one minute # to obtain solutions in a short period of time p3 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_min_set_objective() %>% add_relative_targets(matrix(0.15, nrow = 5, ncol = 3)) %>% add_binary_decisions() %>% add_default_solver(time_limit = 60, verbose = FALSE) # crate asymmetric connectivity data by randomly simulating values acm2 < matrix( runif(ncell(sim_zones_pu_raster) ^ 2), nrow = ncell(sim_zones_pu_raster) ) # create multizone problems using the penalties p4 < list( p3, p3 %>% add_asym_connectivity_penalties(penalties[1], data = acm2), p3 %>% add_asym_connectivity_penalties(penalties[2], data = acm2) ) # solve problems s4 < lapply(p4, solve) s4 < lapply(s4, category_layer) s4 < terra::rast(s4) names(s4) < c("basic problem", paste0("acm2 (", penalties,")")) # plot solutions plot(s4, axes = FALSE) ## End(Not run)
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 then this decision class will be used by default.
add_binary_decisions(x)
add_binary_decisions(x)
x 

Conservation planning problems involve making decisions on planning
units. These decisions are then associated with actions (e.g., turning a
planning unit into a protected area). Only a
single decision should be added to a problem()
object.
Note that if multiple decisions are added to an object, then the
last one to be added will be used.
An updated problem()
object with the decisions added to it.
See decisions for an overview of all functions for adding decisions.
Other decisions:
add_proportion_decisions()
,
add_semicontinuous_decisions()
## Not run: # set seed for reproducibility set.seed(500) # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_features < get_sim_zones_features() # create minimal problem with binary decisions p1 < problem(sim_pu_raster, sim_features) %>% add_min_set_objective() %>% add_relative_targets(0.1) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s1 < solve(p1) # plot solution plot(s1, main = "solution", axes = FALSE) # create a matrix with targets for a multizone conservation problem targs < matrix(runif(15, 0.1, 0.2), nrow = 5, ncol = 3) # build multizone conservation problem with binary decisions p2 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_min_set_objective() %>% add_relative_targets(targs) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve the problem s2 < solve(p2) # print solution print(s2) # plot solution plot(category_layer(s2), main = "solution", axes = FALSE) ## End(Not run)
## Not run: # set seed for reproducibility set.seed(500) # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_features < get_sim_zones_features() # create minimal problem with binary decisions p1 < problem(sim_pu_raster, sim_features) %>% add_min_set_objective() %>% add_relative_targets(0.1) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s1 < solve(p1) # plot solution plot(s1, main = "solution", axes = FALSE) # create a matrix with targets for a multizone conservation problem targs < matrix(runif(15, 0.1, 0.2), nrow = 5, ncol = 3) # build multizone conservation problem with binary decisions p2 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_min_set_objective() %>% add_relative_targets(targs) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve the problem s2 < solve(p2) # print solution print(s2) # plot solution plot(category_layer(s2), main = "solution", axes = FALSE) ## End(Not run)
Add penalties to a conservation planning problem to favor solutions that spatially clump planning units together based on the overall boundary length (i.e., total perimeter).
## S4 method for signature 'ConservationProblem,ANY,ANY,ANY,data.frame' add_boundary_penalties(x, penalty, edge_factor, zones, data) ## S4 method for signature 'ConservationProblem,ANY,ANY,ANY,matrix' add_boundary_penalties(x, penalty, edge_factor, zones, data) ## S4 method for signature 'ConservationProblem,ANY,ANY,ANY,ANY' add_boundary_penalties(x, penalty, edge_factor, zones, data)
## S4 method for signature 'ConservationProblem,ANY,ANY,ANY,data.frame' add_boundary_penalties(x, penalty, edge_factor, zones, data) ## S4 method for signature 'ConservationProblem,ANY,ANY,ANY,matrix' add_boundary_penalties(x, penalty, edge_factor, zones, data) ## S4 method for signature 'ConservationProblem,ANY,ANY,ANY,ANY' add_boundary_penalties(x, penalty, edge_factor, zones, data)
x 

penalty 

edge_factor 

zones 

data 

This function adds penalties to a conservation planning problem
to penalize fragmented solutions. It was is inspired by Ball et al.
(2009) and Beyer et al. (2016). The penalty
argument is
equivalent to the boundary length modifier (BLM
) used in
Marxan.
Note that this function can only
be used to represent symmetric relationships between planning units. If
asymmetric relationships are required, use the
add_connectivity_penalties()
function.
An updated problem()
object with the penalties added to it.
The argument to data
can be specified using the following formats.
Note that boundary data must always describe symmetric relationships
between planning units.
data
as a NULL
valueindicating that the data should be
automatically calculated using the boundary_matrix()
function.
This argument is the default.
Note that the boundary data must be supplied
using one of the other formats below if the planning unit data
in the argument to x
do not explicitly contain spatial information
(e.g., planning unit data are a data.frame
or numeric
class).
data
as a matrix
/Matrix
objectwhere rows and columns represent different planning units and the value of each cell represents the amount of shared boundary length between two different planning units. Cells that occur along the matrix diagonal denote the total boundary length associated with each planning unit.
data
as a data.frame
objectwith the columns "id1"
,
"id2"
, and "boundary"
. The "id1"
and "id2"
columns contain
identifiers (indices) for a pair of planning units, and the "boundary"
column contains the amount of shared boundary length between these
two planning units.
Additionally, if the values in the "id1"
and "id2"
columns
contain the same values, then the value denotes the
amount of exposed boundary length (not total boundary).
This format follows the the standard Marxan format for boundary
data (i.e., per the "bound.dat" file).
The boundary penalties are implemented using the following equations. Let
$I$
represent the set of planning units
(indexed by $i$
or $j$
), $Z$
represent
the set of management zones (indexed by $z$
or $y$
), and
$X_{iz}$
represent the decision
variable for planning unit $i$
for in zone $z$
(e.g., with binary
values one indicating if planning unit is allocated or not). Also, let
$p$
represent the argument to penalty
, $E_z$
represent the
argument to edge_factor
, $B_{ij}$
represent the matrix argument
to data
(e.g., generated using boundary_matrix()
), and
$W_{zz}$
represent the matrix argument to zones
.
$\sum_{i}^{I} \sum_{z}^{Z} (p \times W_{zz} B_{ii}) +
\sum_{i}^{I} \sum_{j}^{I} \sum_{z}^{Z}
\sum_{y}^{Z} (2 \times p \times X_{iz} \times X_{jy} \times W_{zy} \times
B_{ij})$
Note that when the problem objective is to maximize some measure of
benefit and not minimize some measure of cost, the term $p$
is
replaced with $p$
.
Ball IR, Possingham HP, and Watts M (2009) Marxan and relatives: Software for spatial conservation prioritisation in Spatial conservation prioritisation: Quantitative methods and computational tools. Eds Moilanen A, Wilson KA, and Possingham HP. Oxford University Press, Oxford, UK.
Beyer HL, Dujardin Y, Watts ME, and Possingham HP (2016) Solving conservation planning problems with integer linear programming. Ecological Modelling, 228: 14–22.
See penalties for an overview of all functions for adding penalties.
Other penalties:
add_asym_connectivity_penalties()
,
add_connectivity_penalties()
,
add_feature_weights()
,
add_linear_penalties()
## Not run: # set seed for reproducibility set.seed(500) # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_features < get_sim_zones_features() # create minimal problem p1 < problem(sim_pu_raster, sim_features) %>% add_min_set_objective() %>% add_relative_targets(0.2) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # create problem with low boundary penalties p2 < p1 %>% add_boundary_penalties(50, 1) # create problem with high boundary penalties but outer edges receive # half the penalty as inner edges p3 < p1 %>% add_boundary_penalties(500, 0.5) # create a problem using precomputed boundary data bmat < boundary_matrix(sim_pu_raster) p4 < p1 %>% add_boundary_penalties(50, 1, data = bmat) # solve problems s1 < c(solve(p1), solve(p2), solve(p3), solve(p4)) names(s1) < c("basic solution", "small penalties", "high penalties", "precomputed data" ) # plot solutions plot(s1, axes = FALSE) # create minimal problem with multiple zones and limit the runtime for # solver to 10 seconds so this example doesn't take too long p5 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_min_set_objective() %>% add_relative_targets(matrix(0.2, nrow = 5, ncol = 3)) %>% add_binary_decisions() %>% add_default_solver(time_limit = 10, verbose = FALSE) # create zone matrix which favors clumping planning units that are # allocated to the same zone together  note that this is the default zm6 < diag(3) print(zm6) # create problem with the zone matrix and low penalties p6 < p5 %>% add_boundary_penalties(50, zone = zm6) # create another problem with the same zone matrix and higher penalties p7 < p5 %>% add_boundary_penalties(500, zone = zm6) # create zone matrix which favors clumping units that are allocated to # different zones together zm8 < matrix(1, ncol = 3, nrow = 3) diag(zm8) < 0 print(zm8) # create problem with the zone matrix p8 < p5 %>% add_boundary_penalties(500, zone = zm8) # create zone matrix which strongly favors clumping units # that are allocated to the same zone together. It will also prefer # clumping planning units in zones 1 and 2 together over having # these planning units with no neighbors in the solution zm9 < diag(3) zm9[upper.tri(zm9)] < c(0.3, 0, 0) zm9[lower.tri(zm9)] < zm9[upper.tri(zm9)] print(zm9) # create problem with the zone matrix p9 < p5 %>% add_boundary_penalties(500, zone = zm9) # create zone matrix which favors clumping planning units in zones 1 and 2 # together, and favors planning units in zone 3 being spread out # (i.e., negative clumping) zm10 < diag(3) zm10[3, 3] < 1 print(zm10) # create problem with the zone matrix p10 < p5 %>% add_boundary_penalties(500, zone = zm10) # solve problems s2 < list(solve(p5), solve(p6), solve(p7), solve(p8), solve(p9), solve(p10)) #convert to category layers for visualization s2 < terra::rast(lapply(s2, category_layer)) names(s2) < c( "basic solution", "within zone clumping (low)", "within zone clumping (high)", "between zone clumping", "within + between clumping", "negative clumping" ) # plot solutions plot(s2, axes = FALSE) ## End(Not run)
## Not run: # set seed for reproducibility set.seed(500) # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_features < get_sim_zones_features() # create minimal problem p1 < problem(sim_pu_raster, sim_features) %>% add_min_set_objective() %>% add_relative_targets(0.2) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # create problem with low boundary penalties p2 < p1 %>% add_boundary_penalties(50, 1) # create problem with high boundary penalties but outer edges receive # half the penalty as inner edges p3 < p1 %>% add_boundary_penalties(500, 0.5) # create a problem using precomputed boundary data bmat < boundary_matrix(sim_pu_raster) p4 < p1 %>% add_boundary_penalties(50, 1, data = bmat) # solve problems s1 < c(solve(p1), solve(p2), solve(p3), solve(p4)) names(s1) < c("basic solution", "small penalties", "high penalties", "precomputed data" ) # plot solutions plot(s1, axes = FALSE) # create minimal problem with multiple zones and limit the runtime for # solver to 10 seconds so this example doesn't take too long p5 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_min_set_objective() %>% add_relative_targets(matrix(0.2, nrow = 5, ncol = 3)) %>% add_binary_decisions() %>% add_default_solver(time_limit = 10, verbose = FALSE) # create zone matrix which favors clumping planning units that are # allocated to the same zone together  note that this is the default zm6 < diag(3) print(zm6) # create problem with the zone matrix and low penalties p6 < p5 %>% add_boundary_penalties(50, zone = zm6) # create another problem with the same zone matrix and higher penalties p7 < p5 %>% add_boundary_penalties(500, zone = zm6) # create zone matrix which favors clumping units that are allocated to # different zones together zm8 < matrix(1, ncol = 3, nrow = 3) diag(zm8) < 0 print(zm8) # create problem with the zone matrix p8 < p5 %>% add_boundary_penalties(500, zone = zm8) # create zone matrix which strongly favors clumping units # that are allocated to the same zone together. It will also prefer # clumping planning units in zones 1 and 2 together over having # these planning units with no neighbors in the solution zm9 < diag(3) zm9[upper.tri(zm9)] < c(0.3, 0, 0) zm9[lower.tri(zm9)] < zm9[upper.tri(zm9)] print(zm9) # create problem with the zone matrix p9 < p5 %>% add_boundary_penalties(500, zone = zm9) # create zone matrix which favors clumping planning units in zones 1 and 2 # together, and favors planning units in zone 3 being spread out # (i.e., negative clumping) zm10 < diag(3) zm10[3, 3] < 1 print(zm10) # create problem with the zone matrix p10 < p5 %>% add_boundary_penalties(500, zone = zm10) # solve problems s2 < list(solve(p5), solve(p6), solve(p7), solve(p8), solve(p9), solve(p10)) #convert to category layers for visualization s2 < terra::rast(lapply(s2, category_layer)) names(s2) < c( "basic solution", "within zone clumping (low)", "within zone clumping (high)", "between zone clumping", "within + between clumping", "negative clumping" ) # plot solutions plot(s2, axes = FALSE) ## End(Not run)
Specify that the CBC (COINOR branch and cut) software should be used to solve a conservation planning problem (Forrest & LougeeHeimer 2005). This function can also be used to customize the behavior of the solver. It requires the rcbc package to be installed (only available on GitHub, see below for installation instructions).
add_cbc_solver( x, gap = 0.1, time_limit = .Machine$integer.max, presolve = TRUE, threads = 1, first_feasible = FALSE, start_solution = NULL, verbose = TRUE )
add_cbc_solver( x, gap = 0.1, time_limit = .Machine$integer.max, presolve = TRUE, threads = 1, first_feasible = FALSE, start_solution = NULL, verbose = TRUE )
x 

gap 

time_limit 

presolve 

threads 

first_feasible 

start_solution 

verbose 

CBC is an
opensource mixed integer programming solver that is part of the
Computational Infrastructure for Operations Research (COINOR) project.
This solver seems to have much better performance than the other opensource
solvers (i.e., add_highs_solver()
, 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.
An updated problem()
object with the solver added to it.
The rcbc package is required to use this solver. Since the rcbc package is not available on the the Comprehensive R Archive Network (CRAN), it must be installed from its GitHub repository. To install the rcbc package, please use the following code:
if (!require(remotes)) install.packages("remotes") remotes::install_github("dirkschumacher/rcbc")
Note that you may also need to install several dependencies – such as the Rtools software or system libraries – prior to installing the rcbc package. For further details on installing this package, please consult the online package documentation.
Broadly speaking, the argument to start_solution
must be in the same
format as the planning unit data in the argument to x
.
Further details on the correct format are listed separately
for each of the different planning unit data formats:
x
has numeric
planning unitsThe argument to start_solution
must be a
numeric
vector with each element corresponding to a different planning
unit. It should have the same number of planning units as those
in the argument to x
. Additionally, any planning units missing
cost (NA
) values should also have missing (NA
) values in the
argument to start_solution
.
x
has matrix
planning unitsThe argument to start_solution
must be a
matrix
vector with each row corresponding to a different planning
unit, and each column correspond to a different management zone.
It should have the same number of planning units and zones
as those in the argument to x
. Additionally, any planning units
missing cost (NA
) values for a particular zone should also have a
missing (NA
) values in the argument to start_solution
.
x
has terra::rast()
planning unitsThe argument to start_solution
be a terra::rast()
object where different grid cells (pixels) correspond
to different planning units and layers correspond to
a different management zones. It should have the same dimensionality
(rows, columns, layers), resolution, extent, and coordinate reference
system as the planning units in the argument to x
. Additionally,
any planning units missing cost (NA
) values for a particular zone
should also have missing (NA
) values in the argument to start_solution
.
x
has data.frame
planning unitsThe argument to start_solution
must
be a data.frame
with each column corresponding to a different zone,
each row corresponding to a different planning unit, and cell values
corresponding to the solution value. This means that if a data.frame
object containing the solution also contains additional columns, then
these columns will need to be subsetted prior to using this function
(see below for example with sf::sf()
data).
Additionally, any planning units missing cost
(NA
) values for a particular zone should also have missing (NA
)
values in the argument to start_solution
.
x
has sf::sf()
planning unitsThe argument to start_solution
must be
a sf::sf()
object with each column corresponding to a different
zone, each row corresponding to a different planning unit, and cell values
corresponding to the solution value. This means that if the
sf::sf()
object containing the solution also contains additional
columns, then these columns will need to be subsetted prior to using this
function (see below for example).
Additionally, the argument to start_solution
must also have the same
coordinate reference system as the planning unit data.
Furthermore, any planning units missing cost
(NA
) values for a particular zone should also have missing (NA
)
values in the argument to start_solution
.
Forrest J and LougeeHeimer R (2005) CBC User Guide. In Emerging theory, Methods, and Applications (pp. 257–277). INFORMS, Catonsville, MD. doi:10.1287/educ.1053.0020.
Other solvers:
add_cplex_solver()
,
add_default_solver()
,
add_gurobi_solver()
,
add_highs_solver()
,
add_lsymphony_solver
,
add_rsymphony_solver()
## Not run: # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() # create problem p1 < problem(sim_pu_raster, sim_features) %>% add_min_set_objective() %>% add_relative_targets(0.1) %>% add_binary_decisions() %>% add_cbc_solver(gap = 0, verbose = FALSE) # generate solution %>% s1 < solve(p1) # plot solution plot(s1, main = "solution", axes = FALSE) # create a similar problem with boundary length penalties and # specify the solution from the previous run as a starting solution p2 < problem(sim_pu_raster, sim_features) %>% add_min_set_objective() %>% add_relative_targets(0.1) %>% add_boundary_penalties(10) %>% add_binary_decisions() %>% add_cbc_solver(gap = 0, start_solution = s1, verbose = FALSE) # generate solution s2 < solve(p2) # plot solution plot(s2, main = "solution with boundary penalties", axes = FALSE) ## End(Not run)
## Not run: # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() # create problem p1 < problem(sim_pu_raster, sim_features) %>% add_min_set_objective() %>% add_relative_targets(0.1) %>% add_binary_decisions() %>% add_cbc_solver(gap = 0, verbose = FALSE) # generate solution %>% s1 < solve(p1) # plot solution plot(s1, main = "solution", axes = FALSE) # create a similar problem with boundary length penalties and # specify the solution from the previous run as a starting solution p2 < problem(sim_pu_raster, sim_features) %>% add_min_set_objective() %>% add_relative_targets(0.1) %>% add_boundary_penalties(10) %>% add_binary_decisions() %>% add_cbc_solver(gap = 0, start_solution = s1, verbose = FALSE) # generate solution s2 < solve(p2) # plot solution plot(s2, main = "solution with boundary penalties", axes = FALSE) ## End(Not run)
Add penalties to a conservation planning problem to account for
symmetric connectivity between planning units.
Symmetric connectivity data describe connectivity information that is not
directional. For example, symmetric connectivity data could describe which
planning units are adjacent to each other (see adjacency_matrix()
),
or which planning units are within threshold distance of each other (see
proximity_matrix()
).
## S4 method for signature 'ConservationProblem,ANY,ANY,matrix' add_connectivity_penalties(x, penalty, zones, data) ## S4 method for signature 'ConservationProblem,ANY,ANY,Matrix' add_connectivity_penalties(x, penalty, zones, data) ## S4 method for signature 'ConservationProblem,ANY,ANY,data.frame' add_connectivity_penalties(x, penalty, zones, data) ## S4 method for signature 'ConservationProblem,ANY,ANY,dgCMatrix' add_connectivity_penalties(x, penalty, zones, data) ## S4 method for signature 'ConservationProblem,ANY,ANY,array' add_connectivity_penalties(x, penalty, zones, data)
## S4 method for signature 'ConservationProblem,ANY,ANY,matrix' add_connectivity_penalties(x, penalty, zones, data) ## S4 method for signature 'ConservationProblem,ANY,ANY,Matrix' add_connectivity_penalties(x, penalty, zones, data) ## S4 method for signature 'ConservationProblem,ANY,ANY,data.frame' add_connectivity_penalties(x, penalty, zones, data) ## S4 method for signature 'ConservationProblem,ANY,ANY,dgCMatrix' add_connectivity_penalties(x, penalty, zones, data) ## S4 method for signature 'ConservationProblem,ANY,ANY,array' add_connectivity_penalties(x, penalty, zones, data)
x 

penalty 

zones 

data 

This function adds penalties to conservation planning problem to penalize solutions that have low connectivity. Specifically, it favors pairwise connections between planning units that have high connectivity values (based on Önal and Briers 2002).
An updated problem()
object with the penalties added to it.
The argument to data
can be specified using several different formats.
data
as a matrix
/Matrix
objectwhere rows and columns represent
different planning units and the value of each cell represents the
strength of connectivity between two different planning units. Cells
that occur along the matrix diagonal are treated as weights which
indicate that planning units are more desirable in the solution.
The argument to zones
can be used to control
the strength of connectivity between planning units in different zones.
The default argument for zones
is to treat planning units
allocated to different zones as having zero connectivity.
data
as a data.frame
objectcontaining columns that are named
"id1"
, "id2"
, and "boundary"
. Here, each row
denotes the connectivity between a pair of planning units
(per values in the "id1"
and "id2"
columns) following the
Marxan format.
If the argument to x
contains multiple zones, then the
"zone1"
and "zone2"
columns can optionally be provided to manually
specify the connectivity values between planning units when they are
allocated to specific zones. If the "zone1"
and
"zone2"
columns are present, then the argument to zones
must be
NULL
.
data
as an array
objectcontaining fourdimensions where cell values
indicate the strength of connectivity between planning units
when they are assigned to specific management zones. The first two
dimensions (i.e., rows and columns) indicate the strength of
connectivity between different planning units and the second two
dimensions indicate the different management zones. Thus
the data[1, 2, 3, 4]
indicates the strength of
connectivity between planning unit 1 and planning unit 2 when planning
unit 1 is assigned to zone 3 and planning unit 2 is assigned to zone 4.
The connectivity penalties are implemented using the following equations.
Let $I$
represent the set of planning units
(indexed by $i$
or $j$
), $Z$
represent the set
of management zones (indexed by $z$
or $y$
), and $X_{iz}$
represent the decision variable for planning unit $i$
for in zone
$z$
(e.g., with binary
values one indicating if planning unit is allocated or not). Also, let
$p$
represent the argument to penalty
, $D$
represent the
argument to data
, and $W$
represent the argument
to zones
.
If the argument to data
is supplied as a matrix
or
Matrix
object, then the penalties are calculated as:
$\sum_{i}^{I} \sum_{j}^{I} \sum_{z}^{Z} \sum_{y}^{Z} (p \times X_{iz}
\times X_{jy} \times D_{ij} \times W_{zy})$
Otherwise, if the argument to data
is supplied as a
data.frame
or array
object, then the penalties are
calculated as:
$\sum_{i}^{I} \sum_{j}^{I} \sum_{z}^{Z} \sum_{y}^{Z} (p \times X_{iz}
\times X_{jy} \times D_{ijzy})$
Note that when the problem objective is to maximize some measure of
benefit and not minimize some measure of cost, the term $p$
is
replaced with $p$
.
In previous versions, this function aimed to handle both symmetric and
asymmetric connectivity data. This meant that the mathematical
formulation used to account for asymmetric connectivity was different
to that implemented by the Marxan software
(see Beger et al. for details). To ensure that asymmetric connectivity is
handled in a similar manner to the Marxan software, the
add_asym_connectivity_penalties()
function should now be used for
asymmetric connectivity data.
Beger M, Linke S, Watts M, Game E, Treml E, Ball I, and Possingham, HP (2010) Incorporating asymmetric connectivity into spatial decision making for conservation, Conservation Letters, 3: 359–368.
Önal H, and Briers RA (2002) Incorporating spatial criteria in optimum reserve network selection. Proceedings of the Royal Society of London. Series B: Biological Sciences, 269: 2437–2441.
See penalties for an overview of all functions for adding penalties.
Additionally, see add_asym_connectivity_penalties()
to account for
asymmetric connectivity between planning units.
Other penalties:
add_asym_connectivity_penalties()
,
add_boundary_penalties()
,
add_feature_weights()
,
add_linear_penalties()
## Not run: # set seed for reproducibility set.seed(600) # load data sim_pu_polygons < get_sim_pu_polygons() sim_features < get_sim_features() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_features < get_sim_zones_features() # create basic problem p1 < problem(sim_pu_polygons, sim_features, "cost") %>% add_min_set_objective() %>% add_relative_targets(0.2) %>% add_default_solver(verbose = FALSE) # create a symmetric connectivity matrix where the connectivity between # two planning units corresponds to their shared boundary length b_matrix < boundary_matrix(sim_pu_polygons) # rescale matrix values to have a maximum value of 1 b_matrix < rescale_matrix(b_matrix, max = 1) # visualize connectivity matrix Matrix::image(b_matrix) # create a symmetric connectivity matrix where the connectivity between # two planning units corresponds to their spatial proximity # i.e., planning units that are further apart share less connectivity centroids < sf::st_coordinates( suppressWarnings(sf::st_centroid(sim_pu_polygons)) ) d_matrix < (1 / (Matrix::Matrix(as.matrix(dist(centroids))) + 1)) # rescale matrix values to have a maximum value of 1 d_matrix < rescale_matrix(d_matrix, max = 1) # remove connections between planning units with values below a threshold to # reduce runtime d_matrix[d_matrix < 0.8] < 0 # visualize connectivity matrix Matrix::image(d_matrix) # create a symmetric connectivity matrix where the connectivity # between adjacent two planning units corresponds to their combined # value in a column of the planning unit data # for example, this column could describe the extent of native vegetation in # each planning unit and we could use connectivity penalties to identify # solutions that cluster planning units together that both contain large # amounts of native vegetation c_matrix < connectivity_matrix(sim_pu_polygons, "cost") # rescale matrix values to have a maximum value of 1 c_matrix < rescale_matrix(c_matrix, max = 1) # visualize connectivity matrix Matrix::image(c_matrix) # create penalties penalties < c(10, 25) # create problems using the different connectivity matrices and penalties p2 < list( p1, p1 %>% add_connectivity_penalties(penalties[1], data = b_matrix), p1 %>% add_connectivity_penalties(penalties[2], data = b_matrix), p1 %>% add_connectivity_penalties(penalties[1], data = d_matrix), p1 %>% add_connectivity_penalties(penalties[2], data = d_matrix), p1 %>% add_connectivity_penalties(penalties[1], data = c_matrix), p1 %>% add_connectivity_penalties(penalties[2], data = c_matrix) ) # solve problems s2 < lapply(p2, solve) # create single object with all solutions s2 < sf::st_sf( tibble::tibble( p2_1 = s2[[1]]$solution_1, p2_2 = s2[[2]]$solution_1, p2_3 = s2[[3]]$solution_1, p2_4 = s2[[4]]$solution_1, p2_5 = s2[[5]]$solution_1, p2_6 = s2[[6]]$solution_1, p2_7 = s2[[7]]$solution_1 ), geometry = sf::st_geometry(s2[[1]]) ) names(s2)[1:7] < c( "basic problem", paste0("b_matrix (", penalties,")"), paste0("d_matrix (", penalties,")"), paste0("c_matrix (", penalties,")") ) # plot solutions plot(s2) # create minimal multizone problem and limit solver to one minute # to obtain solutions in a short period of time p3 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_min_set_objective() %>% add_relative_targets(matrix(0.15, nrow = 5, ncol = 3)) %>% add_binary_decisions() %>% add_default_solver(time_limit = 60, verbose = FALSE) # create matrix showing which planning units are adjacent to other units a_matrix < adjacency_matrix(sim_zones_pu_raster) # visualize matrix Matrix::image(a_matrix) # create a zone matrix where connectivities are only present between # planning units that are allocated to the same zone zm1 < as(diag(3), "Matrix") # print zone matrix print(zm1) # create a zone matrix where connectivities are strongest between # planning units allocated to different zones zm2 < matrix(1, ncol = 3, nrow = 3) diag(zm2) < 0 zm2 < as(zm2, "Matrix") # print zone matrix print(zm2) # create a zone matrix that indicates that connectivities between planning # units assigned to the same zone are much higher than connectivities # assigned to different zones zm3 < matrix(0.1, ncol = 3, nrow = 3) diag(zm3) < 1 zm3 < as(zm3, "Matrix") # print zone matrix print(zm3) # create a zone matrix that indicates that connectivities between planning # units allocated to zone 1 are very high, connectivities between planning # units allocated to zones 1 and 2 are moderately high, and connectivities # planning units allocated to other zones are low zm4 < matrix(0.1, ncol = 3, nrow = 3) zm4[1, 1] < 1 zm4[1, 2] < 0.5 zm4[2, 1] < 0.5 zm4 < as(zm4, "Matrix") # print zone matrix print(zm4) # create a zone matrix with strong connectivities between planning units # allocated to the same zone, moderate connectivities between planning # unit allocated to zone 1 and zone 2, and negative connectivities between # planning units allocated to zone 3 and the other two zones zm5 < matrix(1, ncol = 3, nrow = 3) zm5[1, 2] < 0.5 zm5[2, 1] < 0.5 diag(zm5) < 1 zm5 < as(zm5, "Matrix") # print zone matrix print(zm5) # create vector of penalties to use creating problems penalties2 < c(5, 15) # create multizone problems using the adjacent connectivity matrix and # different zone matrices p4 < list( p3, p3 %>% add_connectivity_penalties(penalties2[1], zm1, a_matrix), p3 %>% add_connectivity_penalties(penalties2[2], zm1, a_matrix), p3 %>% add_connectivity_penalties(penalties2[1], zm2, a_matrix), p3 %>% add_connectivity_penalties(penalties2[2], zm2, a_matrix), p3 %>% add_connectivity_penalties(penalties2[1], zm3, a_matrix), p3 %>% add_connectivity_penalties(penalties2[2], zm3, a_matrix), p3 %>% add_connectivity_penalties(penalties2[1], zm4, a_matrix), p3 %>% add_connectivity_penalties(penalties2[2], zm4, a_matrix), p3 %>% add_connectivity_penalties(penalties2[1], zm5, a_matrix), p3 %>% add_connectivity_penalties(penalties2[2], zm5, a_matrix) ) # solve problems s4 < lapply(p4, solve) s4 < lapply(s4, category_layer) s4 < terra::rast(s4) names(s4) < c( "basic problem", paste0("zm", rep(seq_len(5), each = 2), " (", rep(penalties2, 2), ")") ) # plot solutions plot(s4, axes = FALSE) # create an array to manually specify the connectivities between # each planning unit when they are allocated to each different zone # for realworld problems, these connectivities would be generated using # data  but here these connectivity values are assigned as random # ones or zeros c_array < array(0, c(rep(ncell(sim_zones_pu_raster[[1]]), 2), 3, 3)) for (z1 in seq_len(3)) for (z2 in seq_len(3)) c_array[, , z1, z2] < round( runif(ncell(sim_zones_pu_raster[[1]]) ^ 2, 0, 0.505) ) # create a problem with the manually specified connectivity array # note that the zones argument is set to NULL because the connectivity # data is an array p5 < list( p3, p3 %>% add_connectivity_penalties(15, zones = NULL, c_array) ) # solve problems s5 < lapply(p5, solve) s5 < lapply(s5, category_layer) s5 < terra::rast(s5) names(s5) < c("basic problem", "connectivity array") # plot solutions plot(s5, axes = FALSE) ## End(Not run)
## Not run: # set seed for reproducibility set.seed(600) # load data sim_pu_polygons < get_sim_pu_polygons() sim_features < get_sim_features() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_features < get_sim_zones_features() # create basic problem p1 < problem(sim_pu_polygons, sim_features, "cost") %>% add_min_set_objective() %>% add_relative_targets(0.2) %>% add_default_solver(verbose = FALSE) # create a symmetric connectivity matrix where the connectivity between # two planning units corresponds to their shared boundary length b_matrix < boundary_matrix(sim_pu_polygons) # rescale matrix values to have a maximum value of 1 b_matrix < rescale_matrix(b_matrix, max = 1) # visualize connectivity matrix Matrix::image(b_matrix) # create a symmetric connectivity matrix where the connectivity between # two planning units corresponds to their spatial proximity # i.e., planning units that are further apart share less connectivity centroids < sf::st_coordinates( suppressWarnings(sf::st_centroid(sim_pu_polygons)) ) d_matrix < (1 / (Matrix::Matrix(as.matrix(dist(centroids))) + 1)) # rescale matrix values to have a maximum value of 1 d_matrix < rescale_matrix(d_matrix, max = 1) # remove connections between planning units with values below a threshold to # reduce runtime d_matrix[d_matrix < 0.8] < 0 # visualize connectivity matrix Matrix::image(d_matrix) # create a symmetric connectivity matrix where the connectivity # between adjacent two planning units corresponds to their combined # value in a column of the planning unit data # for example, this column could describe the extent of native vegetation in # each planning unit and we could use connectivity penalties to identify # solutions that cluster planning units together that both contain large # amounts of native vegetation c_matrix < connectivity_matrix(sim_pu_polygons, "cost") # rescale matrix values to have a maximum value of 1 c_matrix < rescale_matrix(c_matrix, max = 1) # visualize connectivity matrix Matrix::image(c_matrix) # create penalties penalties < c(10, 25) # create problems using the different connectivity matrices and penalties p2 < list( p1, p1 %>% add_connectivity_penalties(penalties[1], data = b_matrix), p1 %>% add_connectivity_penalties(penalties[2], data = b_matrix), p1 %>% add_connectivity_penalties(penalties[1], data = d_matrix), p1 %>% add_connectivity_penalties(penalties[2], data = d_matrix), p1 %>% add_connectivity_penalties(penalties[1], data = c_matrix), p1 %>% add_connectivity_penalties(penalties[2], data = c_matrix) ) # solve problems s2 < lapply(p2, solve) # create single object with all solutions s2 < sf::st_sf( tibble::tibble( p2_1 = s2[[1]]$solution_1, p2_2 = s2[[2]]$solution_1, p2_3 = s2[[3]]$solution_1, p2_4 = s2[[4]]$solution_1, p2_5 = s2[[5]]$solution_1, p2_6 = s2[[6]]$solution_1, p2_7 = s2[[7]]$solution_1 ), geometry = sf::st_geometry(s2[[1]]) ) names(s2)[1:7] < c( "basic problem", paste0("b_matrix (", penalties,")"), paste0("d_matrix (", penalties,")"), paste0("c_matrix (", penalties,")") ) # plot solutions plot(s2) # create minimal multizone problem and limit solver to one minute # to obtain solutions in a short period of time p3 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_min_set_objective() %>% add_relative_targets(matrix(0.15, nrow = 5, ncol = 3)) %>% add_binary_decisions() %>% add_default_solver(time_limit = 60, verbose = FALSE) # create matrix showing which planning units are adjacent to other units a_matrix < adjacency_matrix(sim_zones_pu_raster) # visualize matrix Matrix::image(a_matrix) # create a zone matrix where connectivities are only present between # planning units that are allocated to the same zone zm1 < as(diag(3), "Matrix") # print zone matrix print(zm1) # create a zone matrix where connectivities are strongest between # planning units allocated to different zones zm2 < matrix(1, ncol = 3, nrow = 3) diag(zm2) < 0 zm2 < as(zm2, "Matrix") # print zone matrix print(zm2) # create a zone matrix that indicates that connectivities between planning # units assigned to the same zone are much higher than connectivities # assigned to different zones zm3 < matrix(0.1, ncol = 3, nrow = 3) diag(zm3) < 1 zm3 < as(zm3, "Matrix") # print zone matrix print(zm3) # create a zone matrix that indicates that connectivities between planning # units allocated to zone 1 are very high, connectivities between planning # units allocated to zones 1 and 2 are moderately high, and connectivities # planning units allocated to other zones are low zm4 < matrix(0.1, ncol = 3, nrow = 3) zm4[1, 1] < 1 zm4[1, 2] < 0.5 zm4[2, 1] < 0.5 zm4 < as(zm4, "Matrix") # print zone matrix print(zm4) # create a zone matrix with strong connectivities between planning units # allocated to the same zone, moderate connectivities between planning # unit allocated to zone 1 and zone 2, and negative connectivities between # planning units allocated to zone 3 and the other two zones zm5 < matrix(1, ncol = 3, nrow = 3) zm5[1, 2] < 0.5 zm5[2, 1] < 0.5 diag(zm5) < 1 zm5 < as(zm5, "Matrix") # print zone matrix print(zm5) # create vector of penalties to use creating problems penalties2 < c(5, 15) # create multizone problems using the adjacent connectivity matrix and # different zone matrices p4 < list( p3, p3 %>% add_connectivity_penalties(penalties2[1], zm1, a_matrix), p3 %>% add_connectivity_penalties(penalties2[2], zm1, a_matrix), p3 %>% add_connectivity_penalties(penalties2[1], zm2, a_matrix), p3 %>% add_connectivity_penalties(penalties2[2], zm2, a_matrix), p3 %>% add_connectivity_penalties(penalties2[1], zm3, a_matrix), p3 %>% add_connectivity_penalties(penalties2[2], zm3, a_matrix), p3 %>% add_connectivity_penalties(penalties2[1], zm4, a_matrix), p3 %>% add_connectivity_penalties(penalties2[2], zm4, a_matrix), p3 %>% add_connectivity_penalties(penalties2[1], zm5, a_matrix), p3 %>% add_connectivity_penalties(penalties2[2], zm5, a_matrix) ) # solve problems s4 < lapply(p4, solve) s4 < lapply(s4, category_layer) s4 < terra::rast(s4) names(s4) < c( "basic problem", paste0("zm", rep(seq_len(5), each = 2), " (", rep(penalties2, 2), ")") ) # plot solutions plot(s4, axes = FALSE) # create an array to manually specify the connectivities between # each planning unit when they are allocated to each different zone # for realworld problems, these connectivities would be generated using # data  but here these connectivity values are assigned as random # ones or zeros c_array < array(0, c(rep(ncell(sim_zones_pu_raster[[1]]), 2), 3, 3)) for (z1 in seq_len(3)) for (z2 in seq_len(3)) c_array[, , z1, z2] < round( runif(ncell(sim_zones_pu_raster[[1]]) ^ 2, 0, 0.505) ) # create a problem with the manually specified connectivity array # note that the zones argument is set to NULL because the connectivity # data is an array p5 < list( p3, p3 %>% add_connectivity_penalties(15, zones = NULL, c_array) ) # solve problems s5 < lapply(p5, solve) s5 < lapply(s5, category_layer) s5 < terra::rast(s5) names(s5) < c("basic problem", "connectivity array") # plot solutions plot(s5, axes = FALSE) ## End(Not run)
Add constraints to a conservation planning problem to ensure that all selected planning units are spatially connected with each other and form a single contiguous unit.
## S4 method for signature 'ConservationProblem,ANY,ANY' add_contiguity_constraints(x, zones, data) ## S4 method for signature 'ConservationProblem,ANY,data.frame' add_contiguity_constraints(x, zones, data) ## S4 method for signature 'ConservationProblem,ANY,matrix' add_contiguity_constraints(x, zones, data)
## S4 method for signature 'ConservationProblem,ANY,ANY' add_contiguity_constraints(x, zones, data) ## S4 method for signature 'ConservationProblem,ANY,data.frame' add_contiguity_constraints(x, zones, data) ## S4 method for signature 'ConservationProblem,ANY,matrix' add_contiguity_constraints(x, zones, data)
x 

zones 

data 

This function uses connection data to identify solutions that form a single contiguous unit. It was inspired by the mathematical formulations detailed in Önal and Briers (2006).
An updated problem()
object with the constraints added to it.
The argument to data
can be specified using the following formats.
data
as a NULL
valueindicating that connection data should be
calculated automatically using the adjacency_matrix()
function.
This is the default argument.
Note that the connection data must be manually defined
using one of the other formats below when the planning unit data
in the argument to x
is not spatially referenced (e.g.,
in data.frame
or numeric
format).
data
as a matrix
/Matrix
objectwhere rows and columns represent
different planning units and the value of each cell indicates if the
two planning units are connected or not. Cell values should be binary
numeric
values (i.e., one or zero). Cells that occur along the
matrix diagonal have no effect on the solution at all because each
planning unit cannot be a connected with itself.
data
as a data.frame
objectcontaining columns that are named
"id1"
, "id2"
, and "boundary"
. Here, each row
denotes the connectivity between two planning units following the
Marxan format. The "boundary"
column should contain
binary numeric
values that indicate if the two planning units
specified in the "id1"
and "id2"
columns are connected
or not. This data can be used to describe symmetric or
asymmetric relationships between planning units. By default,
input data is assumed to be symmetric unless asymmetric data is
also included (e.g., if data is present for planning units 2 and 3, then
the same amount of connectivity is expected for planning units 3 and 2,
unless connectivity data is also provided for planning units 3 and 2).
In early versions, this function was named as the
add_connected_constraints()
function.
Önal H and Briers RA (2006) Optimal selection of a connected reserve network. Operations Research, 54: 379–388.
See constraints for an overview of all functions for adding constraints.
Other constraints:
add_feature_contiguity_constraints()
,
add_linear_constraints()
,
add_locked_in_constraints()
,
add_locked_out_constraints()
,
add_mandatory_allocation_constraints()
,
add_manual_bounded_constraints()
,
add_manual_locked_constraints()
,
add_neighbor_constraints()
## Not run: # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_features < get_sim_zones_features() # create minimal problem p1 < problem(sim_pu_raster, sim_features) %>% add_min_set_objective() %>% add_relative_targets(0.2) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # create problem with added connected constraints p2 < p1 %>% add_contiguity_constraints() # solve problems s1 < c(solve(p1), solve(p2)) names(s1) < c("basic solution", "connected solution") # plot solutions plot(s1, axes = FALSE) # create minimal problem with multiple zones, and limit the solver to # 30 seconds to obtain solutions in a feasible period of time p3 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_min_set_objective() %>% add_relative_targets(matrix(0.2, ncol = 3, nrow = 5)) %>% add_binary_decisions() %>% add_default_solver(time_limit = 30, verbose = FALSE) # create problem with added constraints to ensure that the planning units # allocated to each zone form a separate contiguous unit z4 < diag(3) print(z4) p4 < p3 %>% add_contiguity_constraints(z4) # create problem with added constraints to ensure that the planning # units allocated to each zone form a separate contiguous unit, # except for planning units allocated to zone 3 which do not need # form a single contiguous unit z5 < diag(3) z5[3, 3] < 0 print(z5) p5 < p3 %>% add_contiguity_constraints(z5) # create problem with added constraints that ensure that the planning # units allocated to zones 1 and 2 form a contiguous unit z6 < diag(3) z6[1, 2] < 1 z6[2, 1] < 1 print(z6) p6 < p3 %>% add_contiguity_constraints(z6) # solve problems s2 < lapply(list(p3, p4, p5, p6), solve) s2 < lapply(s2, category_layer) s2 < terra::rast(s2) names(s2) < c("basic solution", "p4", "p5", "p6") # plot solutions plot(s2, axes = FALSE) # create a problem that has a main "reserve zone" and a secondary # "corridor zone" to connect up import areas. Here, each feature has a # target of 50% of its distribution. If a planning unit is allocated to the # "reserve zone", then the prioritization accrues 100% of the amount of # each feature in the planning unit. If a planning unit is allocated to the # "corridor zone" then the prioritization accrues 40% of the amount of each # feature in the planning unit. Also, the cost of managing a planning unit # in the "corridor zone" is 30% of that when it is managed as the # "reserve zone". Finally, the problem has constraints which # ensure that all of the selected planning units form a single contiguous # unit, so that the planning units allocated to the "corridor zone" can # link up the planning units allocated to the "reserve zone" # create planning unit data pus < sim_zones_pu_raster[[c(1, 1)]] pus[[2]] < pus[[2]] * 0.3 print(pus) # create biodiversity data fts < zones( sim_features, sim_features * 0.4, feature_names = names(sim_features), zone_names = c("reserve zone", "corridor zone") ) print(fts) # create targets targets < tibble::tibble( feature = names(sim_features), zone = list(zone_names(fts))[rep(1, 5)], target = terra::global(sim_features, "sum", na.rm = TRUE)[[1]] * 0.5, type = rep("absolute", 5) ) print(targets) # create zones matrix z7 < matrix(1, ncol = 2, nrow = 2) print(z7) # create problem p7 < problem(pus, fts) %>% add_min_set_objective() %>% add_manual_targets(targets) %>% add_contiguity_constraints(z7) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problems s7 < category_layer(solve(p7)) # plot solutions plot(s7, main = "solution", axes = FALSE) ## End(Not run)
## Not run: # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_features < get_sim_zones_features() # create minimal problem p1 < problem(sim_pu_raster, sim_features) %>% add_min_set_objective() %>% add_relative_targets(0.2) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # create problem with added connected constraints p2 < p1 %>% add_contiguity_constraints() # solve problems s1 < c(solve(p1), solve(p2)) names(s1) < c("basic solution", "connected solution") # plot solutions plot(s1, axes = FALSE) # create minimal problem with multiple zones, and limit the solver to # 30 seconds to obtain solutions in a feasible period of time p3 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_min_set_objective() %>% add_relative_targets(matrix(0.2, ncol = 3, nrow = 5)) %>% add_binary_decisions() %>% add_default_solver(time_limit = 30, verbose = FALSE) # create problem with added constraints to ensure that the planning units # allocated to each zone form a separate contiguous unit z4 < diag(3) print(z4) p4 < p3 %>% add_contiguity_constraints(z4) # create problem with added constraints to ensure that the planning # units allocated to each zone form a separate contiguous unit, # except for planning units allocated to zone 3 which do not need # form a single contiguous unit z5 < diag(3) z5[3, 3] < 0 print(z5) p5 < p3 %>% add_contiguity_constraints(z5) # create problem with added constraints that ensure that the planning # units allocated to zones 1 and 2 form a contiguous unit z6 < diag(3) z6[1, 2] < 1 z6[2, 1] < 1 print(z6) p6 < p3 %>% add_contiguity_constraints(z6) # solve problems s2 < lapply(list(p3, p4, p5, p6), solve) s2 < lapply(s2, category_layer) s2 < terra::rast(s2) names(s2) < c("basic solution", "p4", "p5", "p6") # plot solutions plot(s2, axes = FALSE) # create a problem that has a main "reserve zone" and a secondary # "corridor zone" to connect up import areas. Here, each feature has a # target of 50% of its distribution. If a planning unit is allocated to the # "reserve zone", then the prioritization accrues 100% of the amount of # each feature in the planning unit. If a planning unit is allocated to the # "corridor zone" then the prioritization accrues 40% of the amount of each # feature in the planning unit. Also, the cost of managing a planning unit # in the "corridor zone" is 30% of that when it is managed as the # "reserve zone". Finally, the problem has constraints which # ensure that all of the selected planning units form a single contiguous # unit, so that the planning units allocated to the "corridor zone" can # link up the planning units allocated to the "reserve zone" # create planning unit data pus < sim_zones_pu_raster[[c(1, 1)]] pus[[2]] < pus[[2]] * 0.3 print(pus) # create biodiversity data fts < zones( sim_features, sim_features * 0.4, feature_names = names(sim_features), zone_names = c("reserve zone", "corridor zone") ) print(fts) # create targets targets < tibble::tibble( feature = names(sim_features), zone = list(zone_names(fts))[rep(1, 5)], target = terra::global(sim_features, "sum", na.rm = TRUE)[[1]] * 0.5, type = rep("absolute", 5) ) print(targets) # create zones matrix z7 < matrix(1, ncol = 2, nrow = 2) print(z7) # create problem p7 < problem(pus, fts) %>% add_min_set_objective() %>% add_manual_targets(targets) %>% add_contiguity_constraints(z7) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problems s7 < category_layer(solve(p7)) # plot solutions plot(s7, main = "solution", axes = FALSE) ## End(Not run)
Specify that the IBM CPLEX software should be used to solve a conservation planning problem (IBM 2017) . This function can also be used to customize the behavior of the solver. It requires the cplexAPI package to be installed (see below for installation instructions).
add_cplex_solver( x, gap = 0.1, time_limit = .Machine$integer.max, presolve = TRUE, threads = 1, verbose = TRUE )
add_cplex_solver( x, gap = 0.1, time_limit = .Machine$integer.max, presolve = TRUE, threads = 1, verbose = TRUE )
x 

gap 

time_limit 

presolve 

threads 

verbose 

IBM CPLEX is a
commercial optimization software. It is faster than
the available open source solvers (e.g., add_lpsymphony_solver()
and
add_rsymphony_solver()
.
Although formal benchmarks examining the performance of this solver for
conservation planning problems have yet to be completed, preliminary
analyses suggest that it performs slightly slower than the Gurobi
solver (i.e., add_gurobi_solver()
).
We recommend using this solver if the Gurobi solver is not available.
Licenses are available for the IBM CPLEX software to academics at no cost
(see < https://www.ibm.com/products/ilogcplexoptimizationstudio/cplexoptimizer>).
An updated problem()
object with the solver added to it.
The cplexAPI package is used to interface with IBM CPLEX software.
To install the package, the IBM CPLEX software must be installed
(see https://www.ibm.com/products/ilogcplexoptimizationstudio/cplexoptimizer). Next, the CPLEX_BIN
environmental variable must be set to specify the file path for the
IBM CPLEX software. For example, on a Linux system,
this variable can be specified by adding the following text to the
~/.bashrc
file:
export CPLEX_BIN="/opt/ibm/ILOG/CPLEX_Studio128/cplex/bin/x8664_linux/cplex"
Please Note that you may need to change the version number in the file path
(i.e., "CPLEX_Studio128"
). After specifying the CPLEX_BIN
environmental variable, the cplexAPI package can be installed.
Since the cplexAPI package is not available on the
the Comprehensive R Archive Network (CRAN), it must be installed from
its GitHub repository. To
install the cplexAPI package, please use the following code:
if (!require(remotes)) install.packages("remotes") remotes::install_github("cran/cplexAPI")
For further details on installing this package, please consult the installation instructions.
IBM (2017) IBM ILOG CPLEX Optimization Studio CPLEX User's Manual. Version 12 Release 8. IBM ILOG CPLEX Division, Incline Village, NV.
Other solvers:
add_cbc_solver()
,
add_default_solver()
,
add_gurobi_solver()
,
add_highs_solver()
,
add_lsymphony_solver
,
add_rsymphony_solver()
## Not run: # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() # create problem p < problem(sim_pu_raster, sim_features) %>% add_min_set_objective() %>% add_relative_targets(0.1) %>% add_binary_decisions() %>% add_cplex_solver(gap = 0.1, time_limit = 5, verbose = FALSE) # generate solution s < solve(p) # plot solution plot(s, main = "solution", axes = FALSE) ## End(Not run)
## Not run: # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() # create problem p < problem(sim_pu_raster, sim_features) %>% add_min_set_objective() %>% add_relative_targets(0.1) %>% add_binary_decisions() %>% add_cplex_solver(gap = 0.1, time_limit = 5, verbose = FALSE) # generate solution s < solve(p) # plot solution plot(s, main = "solution", axes = FALSE) ## End(Not run)
Generate a portfolio of solutions for a conservation planning
problem using Bender's cuts (discussed in Rodrigues
et al. 2000). This is recommended as a replacement for
add_gap_portfolio()
when the Gurobi software is not
available.
add_cuts_portfolio(x, number_solutions = 10)
add_cuts_portfolio(x, number_solutions = 10)
x 

number_solutions 

This strategy for generating a portfolio of solutions involves solving the problem multiple times and adding additional constraints to forbid previously obtained solutions. In general, this strategy is most useful when problems take a long time to solve and benefit from having multiple threads allocated for solving an individual problem.
An updated problem()
object with the portfolio added to it.
In early versions (< 4.0.1), this function was only compatible with
Gurobi (i.e., add_gurobi_solver()
). To provide functionality with
exact algorithm solvers, this function now adds constraints to the
problem formulation to generate multiple solutions.
Rodrigues AS, Cerdeira OJ, and Gaston KJ (2000) Flexibility, efficiency, and accountability: adapting reserve selection algorithms to more complex conservation problems. Ecography, 23: 565–574.
See portfolios for an overview of all functions for adding a portfolio.
Other portfolios:
add_default_portfolio()
,
add_extra_portfolio()
,
add_gap_portfolio()
,
add_shuffle_portfolio()
,
add_top_portfolio()
## Not run: # set seed for reproducibility set.seed(500) # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_features < get_sim_zones_features() # create minimal problem with cuts portfolio p1 < problem(sim_pu_raster, sim_features) %>% add_min_set_objective() %>% add_relative_targets(0.2) %>% add_cuts_portfolio(10) %>% add_default_solver(gap = 0.2, verbose = FALSE) # solve problem and generate 10 solutions within 20% of optimality s1 < solve(p1) # convert portfolio into a multilayer raster object s1 < terra::rast(s1) # plot solutions in portfolio plot(s1, axes = FALSE) # build multizone conservation problem with cuts portfolio p2 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_min_set_objective() %>% add_relative_targets(matrix(runif(15, 0.1, 0.2), nrow = 5, ncol = 3)) %>% add_binary_decisions() %>% add_cuts_portfolio(10) %>% add_default_solver(gap = 0.2, verbose = FALSE) # solve the problem s2 < solve(p2) # print solution str(s2, max.level = 1) # convert each solution in the portfolio into a single category layer s2 < terra::rast(lapply(s2, category_layer)) # plot solutions in portfolio plot(s2, main = "solution", axes = FALSE) ## End(Not run)
## Not run: # set seed for reproducibility set.seed(500) # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_features < get_sim_zones_features() # create minimal problem with cuts portfolio p1 < problem(sim_pu_raster, sim_features) %>% add_min_set_objective() %>% add_relative_targets(0.2) %>% add_cuts_portfolio(10) %>% add_default_solver(gap = 0.2, verbose = FALSE) # solve problem and generate 10 solutions within 20% of optimality s1 < solve(p1) # convert portfolio into a multilayer raster object s1 < terra::rast(s1) # plot solutions in portfolio plot(s1, axes = FALSE) # build multizone conservation problem with cuts portfolio p2 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_min_set_objective() %>% add_relative_targets(matrix(runif(15, 0.1, 0.2), nrow = 5, ncol = 3)) %>% add_binary_decisions() %>% add_cuts_portfolio(10) %>% add_default_solver(gap = 0.2, verbose = FALSE) # solve the problem s2 < solve(p2) # print solution str(s2, max.level = 1) # convert each solution in the portfolio into a single category layer s2 < terra::rast(lapply(s2, category_layer)) # plot solutions in portfolio plot(s2, main = "solution", axes = FALSE) ## End(Not run)
Generate a portfolio containing a single solution.
add_default_portfolio(x)
add_default_portfolio(x)
x 

By default, this is portfolio is added to problem()
objects if no
other portfolios is manually specified.
An updated problem()
object with the portfolio added to it.
See portfolios for an overview of all functions for adding a portfolio.
Other portfolios:
add_cuts_portfolio()
,
add_extra_portfolio()
,
add_gap_portfolio()
,
add_shuffle_portfolio()
,
add_top_portfolio()
Specify that the best solver currently available should be used to solve a conservation planning problem.
add_default_solver(x, ...)
add_default_solver(x, ...)
x 

... 
arguments passed to the solver. 
Ranked from best to worst, the available solvers that can be used are:
add_gurobi_solver()
, add_cplex_solver()
, add_cbc_solver()
,
add_highs_solver()
, add_lpsymphony_solver()
, and finally
add_rsymphony_solver()
.
For information on the performance of different solvers,
please see Schuster et al. (2020).
An updated problem()
object with the solver added to it.
Schuster R, Hanson JO, StrimasMackey M, and Bennett JR (2020). Exact integer linear programming solvers outperform simulated annealing for solving conservation planning problems. PeerJ, 8: e9258.
See solvers for an overview of all functions for adding a solver.
Other solvers:
add_cbc_solver()
,
add_cplex_solver()
,
add_gurobi_solver()
,
add_highs_solver()
,
add_lsymphony_solver
,
add_rsymphony_solver()
Generate a portfolio of solutions for a conservation planning problem by storing feasible solutions discovered during the optimization process. This method is useful for quickly obtaining multiple solutions, but does not provide any guarantees on the number of solutions, or the quality of solutions.
add_extra_portfolio(x)
add_extra_portfolio(x)
x 

This strategy for generating a portfolio requires problems to
be solved using the Gurobi software suite (i.e., using
add_gurobi_solver()
. Specifically, version 8.0.0 (or greater)
of the gurobi package must be installed.
An updated problem()
object with the portfolio added to it.
See portfolios for an overview of all functions for adding a portfolio.
Other portfolios:
add_cuts_portfolio()
,
add_default_portfolio()
,
add_gap_portfolio()
,
add_shuffle_portfolio()
,
add_top_portfolio()
## Not run: # set seed for reproducibility set.seed(600) # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_features < get_sim_zones_features() # create minimal problem with a portfolio for extra solutions p1 < problem(sim_pu_raster, sim_features) %>% add_min_set_objective() %>% add_relative_targets(0.05) %>% add_extra_portfolio() %>% add_default_solver(gap = 0, verbose = FALSE) # solve problem and generate portfolio s1 < solve(p1) # convert portfolio into a multilayer raster object s1 < terra::rast(s1) # print number of solutions found print(terra::nlyr(s1)) # plot solutions plot(s1, axes = FALSE) # create multizone problem with a portfolio for extra solutions p2 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_min_set_objective() %>% add_relative_targets(matrix(runif(15, 0.1, 0.2), nrow = 5, ncol = 3)) %>% add_extra_portfolio() %>% add_default_solver(gap = 0, verbose = FALSE) # solve problem and generate portfolio s2 < solve(p2) # convert each solution in the portfolio into a single category layer s2 < terra::rast(lapply(s2, category_layer)) # print number of solutions found print(terra::nlyr(s2)) # plot solutions in portfolio plot(s2, axes = FALSE) ## End(Not run)
## Not run: # set seed for reproducibility set.seed(600) # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_features < get_sim_zones_features() # create minimal problem with a portfolio for extra solutions p1 < problem(sim_pu_raster, sim_features) %>% add_min_set_objective() %>% add_relative_targets(0.05) %>% add_extra_portfolio() %>% add_default_solver(gap = 0, verbose = FALSE) # solve problem and generate portfolio s1 < solve(p1) # convert portfolio into a multilayer raster object s1 < terra::rast(s1) # print number of solutions found print(terra::nlyr(s1)) # plot solutions plot(s1, axes = FALSE) # create multizone problem with a portfolio for extra solutions p2 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_min_set_objective() %>% add_relative_targets(matrix(runif(15, 0.1, 0.2), nrow = 5, ncol = 3)) %>% add_extra_portfolio() %>% add_default_solver(gap = 0, verbose = FALSE) # solve problem and generate portfolio s2 < solve(p2) # convert each solution in the portfolio into a single category layer s2 < terra::rast(lapply(s2, category_layer)) # print number of solutions found print(terra::nlyr(s2)) # plot solutions in portfolio plot(s2, axes = FALSE) ## End(Not run)
Add constraints to a problem 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_constraints()
function, because they ensure that
each feature is represented in a contiguous unit and not that the entire
solution should form a contiguous unit. Additionally, this function
can use data showing the distribution of dispersible habitat for each
feature to ensure that all features can disperse throughout the areas
designated for their conservation.
## S4 method for signature 'ConservationProblem,ANY,data.frame' add_feature_contiguity_constraints(x, zones, data) ## S4 method for signature 'ConservationProblem,ANY,matrix' add_feature_contiguity_constraints(x, zones, data) ## S4 method for signature 'ConservationProblem,ANY,ANY' add_feature_contiguity_constraints(x, zones, data)
## S4 method for signature 'ConservationProblem,ANY,data.frame' add_feature_contiguity_constraints(x, zones, data) ## S4 method for signature 'ConservationProblem,ANY,matrix' add_feature_contiguity_constraints(x, zones, data) ## S4 method for signature 'ConservationProblem,ANY,ANY' add_feature_contiguity_constraints(x, zones, data)
x 

zones 

data 

This function uses connection data to identify solutions that represent features in contiguous units of dispersible habitat. It was inspired by the mathematical formulations detailed in Önal and Briers (2006) and Cardeira et al. 2010. For an example that has used these constraints, see Hanson et al. (2019). Please note that these constraints require the expanded formulation and therefore cannot be used with feature data that have negative vales. Please note that adding these constraints to a problem will drastically increase the amount of time required to solve it.
An updated problem()
object with the constraints added to it.
The argument to data
can be specified using the following formats.
data
as a NULL
valueconnection
data should be calculated automatically
using the adjacency_matrix()
function. This is the default
argument and means that all adjacent planning units are treated
as potentially dispersible for all features.
Note that the connection data must be manually defined
using one of the other formats below when the planning unit data
in the argument to x
is not spatially referenced (e.g.,
in data.frame
or numeric
format).
data
as amatrix
/Matrix
objectwhere rows and columns represent
different planning units and the value of each cell indicates if the
two planning units are connected or not. Cell values should be binary
numeric
values (i.e., one or zero). Cells that occur along the
matrix diagonal have no effect on the solution at all because each
planning unit cannot be a connected with itself. Note that pairs
of connected planning units are treated as being potentially dispersible
for all features.
data
as a data.frame
objectcontaining columns that are named
"id1"
, "id2"
, and "boundary"
. Here, each row
denotes the connectivity between two planning units following the
Marxan format. The "boundary"
column should contain
binary numeric
values that indicate if the two planning units
specified in the "id1"
and "id2"
columns are connected
or not. This data can be used to describe symmetric or
asymmetric relationships between planning units. By default,
input data is assumed to be symmetric unless asymmetric data is
also included (e.g., if data is present for planning units 2 and 3, then
the same amount of connectivity is expected for planning units 3 and 2,
unless connectivity data is also provided for planning units 3 and 2).
Note that pairs of connected planning units are treated as being
potentially dispersible for all features.
data
as a list
objectcontaining matrix
, Matrix
, or
data.frame
objects showing which planning units
should be treated as connected for each feature. Each element in the
list
should correspond to a different feature (specifically,
a different target in the problem), and should contain a matrix
,
Matrix
, or data.frame
object that follows the conventions
detailed above.
In early versions, it was named as the add_corridor_constraints
function.
Önal H and Briers RA (2006) Optimal selection of a connected reserve network. Operations Research, 54: 379–388.
Cardeira JO, Pinto LS, Cabeza M and Gaston KJ (2010) Species specific connectivity in reservenetwork design using graphs. Biological Conservation, 2: 408–415.
Hanson JO, Fuller RA, & Rhodes JR (2019) Conventional methods for enhancing connectivity in conservation planning do not always maintain gene flow. Journal of Applied Ecology, 56: 913–922.
See constraints for an overview of all functions for adding constraints.
Other constraints:
add_contiguity_constraints()
,
add_linear_constraints()
,
add_locked_in_constraints()
,
add_locked_out_constraints()
,
add_mandatory_allocation_constraints()
,
add_manual_bounded_constraints()
,
add_manual_locked_constraints()
,
add_neighbor_constraints()
## Not run: # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_features < get_sim_zones_features() # create minimal problem p1 < problem(sim_pu_raster, sim_features) %>% add_min_set_objective() %>% add_relative_targets(0.3) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # create problem with contiguity constraints p2 < p1 %>% add_contiguity_constraints() # create problem with constraints to represent features in contiguous # units p3 < p1 %>% add_feature_contiguity_constraints() # create problem with constraints to represent features in contiguous # units that contain highly suitable habitat values # (specifically in the top 5th percentile) cm4 < lapply(seq_len(terra::nlyr(sim_features)), function(i) { # create connectivity matrix using the i'th feature's habitat data m < connectivity_matrix(sim_pu_raster, sim_features[[i]]) # convert matrix to 0/1 values denoting values in top 5th percentile m < round(m > quantile(as.vector(m), 1  0.05, names = FALSE)) # remove 0s from the sparse matrix m < Matrix::drop0(m) # return matrix m }) p4 < p1 %>% add_feature_contiguity_constraints(data = cm4) # solve problems s1 < c(solve(p1), solve(p2), solve(p3), solve(p4)) names(s1) < c( "basic solution", "contiguity constraints", "feature contiguity constraints", "feature contiguity constraints with data" ) # plot solutions plot(s1, axes = FALSE) # create minimal problem with multiple zones, and limit the solver to # 30 seconds to obtain solutions in a feasible period of time p5 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_min_set_objective() %>% add_relative_targets(matrix(0.1, ncol = 3, nrow = 5)) %>% add_binary_decisions() %>% add_default_solver(time_limit = 30, verbose = FALSE) # create problem with contiguity constraints that specify that the # planning units used to conserve each feature in different management # zones must form separate contiguous units p6 < p5 %>% add_feature_contiguity_constraints(diag(3)) # create problem with contiguity constraints that specify that the # planning units used to conserve each feature must form a single # contiguous unit if the planning units are allocated to zones 1 and 2 # and do not need to form a single contiguous unit if they are allocated # to zone 3 zm7 < matrix(0, ncol = 3, nrow = 3) zm7[seq_len(2), seq_len(2)] < 1 print(zm7) p7 < p5 %>% add_feature_contiguity_constraints(zm7) # create problem with contiguity constraints that specify that all of # the planning units in all three of the zones must conserve first feature # in a single contiguous unit but the planning units used to conserve the # remaining features do not need to be contiguous in any way zm8 < lapply( seq_len(number_of_features(sim_zones_features)), function(i) matrix(ifelse(i == 1, 1, 0), ncol = 3, nrow = 3) ) print(zm8) p8 < p5 %>% add_feature_contiguity_constraints(zm8) # solve problems s2 < lapply(list(p5, p6, p7, p8), solve) s2 < terra::rast(lapply(s2, category_layer)) names(s2) < c("p5", "p6", "p7", "p8") # plot solutions plot(s2, axes = FALSE) ## End(Not run)
## Not run: # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_features < get_sim_zones_features() # create minimal problem p1 < problem(sim_pu_raster, sim_features) %>% add_min_set_objective() %>% add_relative_targets(0.3) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # create problem with contiguity constraints p2 < p1 %>% add_contiguity_constraints() # create problem with constraints to represent features in contiguous # units p3 < p1 %>% add_feature_contiguity_constraints() # create problem with constraints to represent features in contiguous # units that contain highly suitable habitat values # (specifically in the top 5th percentile) cm4 < lapply(seq_len(terra::nlyr(sim_features)), function(i) { # create connectivity matrix using the i'th feature's habitat data m < connectivity_matrix(sim_pu_raster, sim_features[[i]]) # convert matrix to 0/1 values denoting values in top 5th percentile m < round(m > quantile(as.vector(m), 1  0.05, names = FALSE)) # remove 0s from the sparse matrix m < Matrix::drop0(m) # return matrix m }) p4 < p1 %>% add_feature_contiguity_constraints(data = cm4) # solve problems s1 < c(solve(p1), solve(p2), solve(p3), solve(p4)) names(s1) < c( "basic solution", "contiguity constraints", "feature contiguity constraints", "feature contiguity constraints with data" ) # plot solutions plot(s1, axes = FALSE) # create minimal problem with multiple zones, and limit the solver to # 30 seconds to obtain solutions in a feasible period of time p5 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_min_set_objective() %>% add_relative_targets(matrix(0.1, ncol = 3, nrow = 5)) %>% add_binary_decisions() %>% add_default_solver(time_limit = 30, verbose = FALSE) # create problem with contiguity constraints that specify that the # planning units used to conserve each feature in different management # zones must form separate contiguous units p6 < p5 %>% add_feature_contiguity_constraints(diag(3)) # create problem with contiguity constraints that specify that the # planning units used to conserve each feature must form a single # contiguous unit if the planning units are allocated to zones 1 and 2 # and do not need to form a single contiguous unit if they are allocated # to zone 3 zm7 < matrix(0, ncol = 3, nrow = 3) zm7[seq_len(2), seq_len(2)] < 1 print(zm7) p7 < p5 %>% add_feature_contiguity_constraints(zm7) # create problem with contiguity constraints that specify that all of # the planning units in all three of the zones must conserve first feature # in a single contiguous unit but the planning units used to conserve the # remaining features do not need to be contiguous in any way zm8 < lapply( seq_len(number_of_features(sim_zones_features)), function(i) matrix(ifelse(i == 1, 1, 0), ncol = 3, nrow = 3) ) print(zm8) p8 < p5 %>% add_feature_contiguity_constraints(zm8) # solve problems s2 < lapply(list(p5, p6, p7, p8), solve) s2 < terra::rast(lapply(s2, category_layer)) names(s2) < c("p5", "p6", "p7", "p8") # plot solutions plot(s2, axes = FALSE) ## End(Not run)
Add features weights to a conservation planning problem. Specifically,
some objective functions aim to maximize (or minimize) a metric that
measures how well a set of features are represented by a solution
(e.g., maximize the number of features that are adequately represented,
add_max_features_objective()
). In such cases,
it may be desirable to prefer the representation of some features
over other features (e.g., features that have higher extinction risk
might be considered more important than those with lower extinction risk).
To achieve this, weights can be used to specify how much more important
it is for a solution to represent particular features compared with other
features.
## S4 method for signature 'ConservationProblem,numeric' add_feature_weights(x, weights) ## S4 method for signature 'ConservationProblem,matrix' add_feature_weights(x, weights)
## S4 method for signature 'ConservationProblem,numeric' add_feature_weights(x, weights) ## S4 method for signature 'ConservationProblem,matrix' add_feature_weights(x, weights)
x 

weights 

Weights can only be applied to problems that have an objective
that is budget limited (e.g., add_max_cover_objective()
,
add_min_shortfall_objective()
).
They can also be applied to problems that aim to maximize phylogenetic
representation (add_max_phylo_div_objective()
) to favor the
representation of specific features over the representation of
some phylogenetic branches. Weights cannot be negative values
and must have values that are equal to or larger than zero.
Note that planning unit costs are scaled to 0.01 to identify
the cheapest solution among multiple optimal solutions. This means
that the optimization process will favor cheaper solutions over solutions
that meet feature targets (or occurrences) when feature weights are
lower than 0.01.
An updated problem()
with the weights added to it.
The argument to weights
can be specified using the following formats.
weights
as a numeric
vectorcontaining weights for each feature. Note that this format cannot be used to specify weights for problems with multiple zones.
weights
as a matrix
objectcontaining weights
for each feature in each zone.
Here, each row corresponds to a different feature in argument to
x
, each column corresponds to a different zone in argument to
x
, and each cell contains the weight value for a given feature
that the solution can to secure in a given zone. Note that
if the problem contains targets created using
add_manual_targets()
then a matrix
should be
supplied containing a single column that indicates that weight for
fulfilling each target.
See penalties for an overview of all functions for adding penalties.
Other penalties:
add_asym_connectivity_penalties()
,
add_boundary_penalties()
,
add_connectivity_penalties()
,
add_linear_penalties()
## Not run: # load package require(ape) # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() sim_phylogeny < get_sim_phylogeny() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_features < get_sim_zones_features() # create minimal problem that aims to maximize the number of features # adequately conserved given a total budget of 3800. Here, each feature # needs 20% of its habitat for it to be considered adequately conserved p1 < problem(sim_pu_raster, sim_features) %>% add_max_features_objective(budget = 3800) %>% add_relative_targets(0.2) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # create weights that assign higher importance to features with less # suitable habitat in the study area w2 < exp((1 / terra::global(sim_features, "sum", na.rm = TRUE)[[1]]) * 200) # create problem using rarity weights p2 < p1 %>% add_feature_weights(w2) # create manually specified weights that assign higher importance to # certain features. These weights could be based on a precalculated index # (e.g., an index measuring extinction risk where higher values # denote higher extinction risk) w3 < c(0, 0, 0, 100, 200) p3 < p1 %>% add_feature_weights(w3) # solve problems s1 < c(solve(p1), solve(p2), solve(p3)) names(s1) < c("equal weights", "rarity weights", "manual weights") # plot solutions plot(s1, axes = FALSE) # plot the example phylogeny par(mfrow = c(1, 1)) plot(sim_phylogeny, main = "simulated phylogeny") # create problem with a maximum phylogenetic diversity objective, # where each feature needs 10% of its distribution to be secured for # it to be adequately conserved and a total budget of 1900 p4 < problem(sim_pu_raster, sim_features) %>% add_max_phylo_div_objective(1900, sim_phylogeny) %>% add_relative_targets(0.1) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s4 < solve(p4) # plot solution plot(s4, main = "solution", axes = FALSE) # find out which features have their targets met r4 < eval_target_coverage_summary(p4, s4) print(r4, width = Inf) # plot the example phylogeny and color the represented features in red plot( sim_phylogeny, main = "represented features", tip.color = replace( rep("black", terra::nlyr(sim_features)), which(r4$met), "red" ) ) # we can see here that the third feature ("layer.3", i.e., # sim_features[[3]]) is not represented in the solution. Let us pretend # that it is absolutely critical this feature is adequately conserved # in the solution. For example, this feature could represent a species # that plays important role in the ecosystem, or a species that is # important commercial activities (e.g., ecotourism). So, to generate # a solution that conserves the third feature whilst also aiming to # maximize phylogenetic diversity, we will create a set of weights that # assign a particularly high weighting to the third feature w5 < c(0, 0, 10000, 0, 0) # we can see that this weighting (i.e., w5[3]) has a much higher value than # the branch lengths in the phylogeny so solutions that represent this # feature be much closer to optimality print(sim_phylogeny$edge.length) # create problem with high weighting for the third feature and solve it s5 < p4 %>% add_feature_weights(w5) %>% solve() # plot solution plot(s5, main = "solution", axes = FALSE) # find which features have their targets met r5 < eval_target_coverage_summary(p4, s5) print(r5, width = Inf) # plot the example phylogeny and color the represented features in red # here we can see that this solution only adequately conserves the # third feature. This means that, given the budget, we are faced with the # tradeoff of conserving either the third feature, or a phylogenetically # diverse set of three different features. plot( sim_phylogeny, main = "represented features", tip.color = replace( rep("black", terra::nlyr(sim_features)), which(r5$met), "red" ) ) # create multizone problem with maximum features objective, # with 10% representation targets for each feature, and set # a budget such that the total maximum expenditure in all zones # cannot exceed 3000 p6 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_max_features_objective(3000) %>% add_relative_targets(matrix(0.1, ncol = 3, nrow = 5)) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # create weights that assign equal weighting for the representation # of each feature in each zone except that it does not matter if # feature 1 is represented in zone 1 and it really important # that feature 3 is really in zone 1 w7 < matrix(1, ncol = 3, nrow = 5) w7[1, 1] < 0 w7[3, 1] < 100 # create problem with weights p7 < p6 %>% add_feature_weights(w7) # solve problems s6 < solve(p6) s7 < solve(p7) # convert solutions to category layers c6 < category_layer(s6) c7 < category_layer(s7) # plot solutions plot(c(c6, c7), main = c("equal weights", "manual weights"), axes = FALSE) # create minimal problem to show the correct method for setting # weights for problems with manual targets p8 < problem(sim_pu_raster, sim_features) %>% add_max_features_objective(budget = 3000) %>% add_manual_targets( data.frame( feature = c("feature_1", "feature_4"), type = "relative", target = 0.1) ) %>% add_feature_weights(matrix(c(1, 200), ncol = 1)) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s8 < solve(p8) # plot solution plot(s8, main = "solution", axes = FALSE) ## End(Not run)
## Not run: # load package require(ape) # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() sim_phylogeny < get_sim_phylogeny() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_features < get_sim_zones_features() # create minimal problem that aims to maximize the number of features # adequately conserved given a total budget of 3800. Here, each feature # needs 20% of its habitat for it to be considered adequately conserved p1 < problem(sim_pu_raster, sim_features) %>% add_max_features_objective(budget = 3800) %>% add_relative_targets(0.2) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # create weights that assign higher importance to features with less # suitable habitat in the study area w2 < exp((1 / terra::global(sim_features, "sum", na.rm = TRUE)[[1]]) * 200) # create problem using rarity weights p2 < p1 %>% add_feature_weights(w2) # create manually specified weights that assign higher importance to # certain features. These weights could be based on a precalculated index # (e.g., an index measuring extinction risk where higher values # denote higher extinction risk) w3 < c(0, 0, 0, 100, 200) p3 < p1 %>% add_feature_weights(w3) # solve problems s1 < c(solve(p1), solve(p2), solve(p3)) names(s1) < c("equal weights", "rarity weights", "manual weights") # plot solutions plot(s1, axes = FALSE) # plot the example phylogeny par(mfrow = c(1, 1)) plot(sim_phylogeny, main = "simulated phylogeny") # create problem with a maximum phylogenetic diversity objective, # where each feature needs 10% of its distribution to be secured for # it to be adequately conserved and a total budget of 1900 p4 < problem(sim_pu_raster, sim_features) %>% add_max_phylo_div_objective(1900, sim_phylogeny) %>% add_relative_targets(0.1) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s4 < solve(p4) # plot solution plot(s4, main = "solution", axes = FALSE) # find out which features have their targets met r4 < eval_target_coverage_summary(p4, s4) print(r4, width = Inf) # plot the example phylogeny and color the represented features in red plot( sim_phylogeny, main = "represented features", tip.color = replace( rep("black", terra::nlyr(sim_features)), which(r4$met), "red" ) ) # we can see here that the third feature ("layer.3", i.e., # sim_features[[3]]) is not represented in the solution. Let us pretend # that it is absolutely critical this feature is adequately conserved # in the solution. For example, this feature could represent a species # that plays important role in the ecosystem, or a species that is # important commercial activities (e.g., ecotourism). So, to generate # a solution that conserves the third feature whilst also aiming to # maximize phylogenetic diversity, we will create a set of weights that # assign a particularly high weighting to the third feature w5 < c(0, 0, 10000, 0, 0) # we can see that this weighting (i.e., w5[3]) has a much higher value than # the branch lengths in the phylogeny so solutions that represent this # feature be much closer to optimality print(sim_phylogeny$edge.length) # create problem with high weighting for the third feature and solve it s5 < p4 %>% add_feature_weights(w5) %>% solve() # plot solution plot(s5, main = "solution", axes = FALSE) # find which features have their targets met r5 < eval_target_coverage_summary(p4, s5) print(r5, width = Inf) # plot the example phylogeny and color the represented features in red # here we can see that this solution only adequately conserves the # third feature. This means that, given the budget, we are faced with the # tradeoff of conserving either the third feature, or a phylogenetically # diverse set of three different features. plot( sim_phylogeny, main = "represented features", tip.color = replace( rep("black", terra::nlyr(sim_features)), which(r5$met), "red" ) ) # create multizone problem with maximum features objective, # with 10% representation targets for each feature, and set # a budget such that the total maximum expenditure in all zones # cannot exceed 3000 p6 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_max_features_objective(3000) %>% add_relative_targets(matrix(0.1, ncol = 3, nrow = 5)) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # create weights that assign equal weighting for the representation # of each feature in each zone except that it does not matter if # feature 1 is represented in zone 1 and it really important # that feature 3 is really in zone 1 w7 < matrix(1, ncol = 3, nrow = 5) w7[1, 1] < 0 w7[3, 1] < 100 # create problem with weights p7 < p6 %>% add_feature_weights(w7) # solve problems s6 < solve(p6) s7 < solve(p7) # convert solutions to category layers c6 < category_layer(s6) c7 < category_layer(s7) # plot solutions plot(c(c6, c7), main = c("equal weights", "manual weights"), axes = FALSE) # create minimal problem to show the correct method for setting # weights for problems with manual targets p8 < problem(sim_pu_raster, sim_features) %>% add_max_features_objective(budget = 3000) %>% add_manual_targets( data.frame( feature = c("feature_1", "feature_4"), type = "relative", target = 0.1) ) %>% add_feature_weights(matrix(c(1, 200), ncol = 1)) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s8 < solve(p8) # plot solution plot(s8, main = "solution", axes = FALSE) ## End(Not run)
Generate a portfolio of solutions for a conservation planning problem by finding a certain number of solutions that are all within a prespecified optimality gap. This method is useful for generating multiple solutions that can be used to calculate selection frequencies for moderate and largesized problems (similar to Marxan).
add_gap_portfolio(x, number_solutions = 10, pool_gap = 0.1)
add_gap_portfolio(x, number_solutions = 10, pool_gap = 0.1)
x 

number_solutions 

pool_gap 

This strategy for generating a portfolio requires problems to
be solved using the Gurobi software suite (i.e., using
add_gurobi_solver()
. Specifically, version 9.0.0 (or greater)
of the gurobi package must be installed.
Note that the number of solutions returned may be less than the argument to
number_solutions
, if the total number of solutions that
meet the optimality gap is less than the number of solutions requested.
Also, note that this portfolio function only works with problems
that have binary decisions (i.e., specified using
add_binary_decisions()
).
An updated problem()
object with the portfolio added to it.
See portfolios for an overview of all functions for adding a portfolio.
Other portfolios:
add_cuts_portfolio()
,
add_default_portfolio()
,
add_extra_portfolio()
,
add_shuffle_portfolio()
,
add_top_portfolio()
## Not run: # set seed for reproducibility set.seed(600) # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_features < get_sim_zones_features() # create minimal problem with a portfolio containing 10 solutions within 20% # of optimality p1 < problem(sim_pu_raster, sim_features) %>% add_min_set_objective() %>% add_relative_targets(0.05) %>% add_gap_portfolio(number_solutions = 5, pool_gap = 0.2) %>% add_default_solver(gap = 0, verbose = FALSE) # solve problem and generate portfolio s1 < solve(p1) # convert portfolio into a multilayer raster s1 < terra::rast(s1) # print number of solutions found print(terra::nlyr(s1)) # plot solutions plot(s1, axes = FALSE) # create multizone problem with a portfolio containing 10 solutions within # 20% of optimality p2 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_min_set_objective() %>% add_relative_targets(matrix(runif(15, 0.1, 0.2), nrow = 5, ncol = 3)) %>% add_gap_portfolio(number_solutions = 5, pool_gap = 0.2) %>% add_default_solver(gap = 0, verbose = FALSE) # solve problem and generate portfolio s2 < solve(p2) # convert portfolio into a multilayer raster of category layers s2 < terra::rast(lapply(s2, category_layer)) # print number of solutions found print(terra::nlyr(s2)) # plot solutions in portfolio plot(s2, axes = FALSE) ## End(Not run)
## Not run: # set seed for reproducibility set.seed(600) # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_features < get_sim_zones_features() # create minimal problem with a portfolio containing 10 solutions within 20% # of optimality p1 < problem(sim_pu_raster, sim_features) %>% add_min_set_objective() %>% add_relative_targets(0.05) %>% add_gap_portfolio(number_solutions = 5, pool_gap = 0.2) %>% add_default_solver(gap = 0, verbose = FALSE) # solve problem and generate portfolio s1 < solve(p1) # convert portfolio into a multilayer raster s1 < terra::rast(s1) # print number of solutions found print(terra::nlyr(s1)) # plot solutions plot(s1, axes = FALSE) # create multizone problem with a portfolio containing 10 solutions within # 20% of optimality p2 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_min_set_objective() %>% add_relative_targets(matrix(runif(15, 0.1, 0.2), nrow = 5, ncol = 3)) %>% add_gap_portfolio(number_solutions = 5, pool_gap = 0.2) %>% add_default_solver(gap = 0, verbose = FALSE) # solve problem and generate portfolio s2 < solve(p2) # convert portfolio into a multilayer raster of category layers s2 < terra::rast(lapply(s2, category_layer)) # print number of solutions found print(terra::nlyr(s2)) # plot solutions in portfolio plot(s2, axes = FALSE) ## End(Not run)
Specify that the Gurobi software should be used to solve a conservation planning problem (Gurobi Optimization LLC 2021). This function can also be used to customize the behavior of the solver. It requires the gurobi package to be installed (see below for installation instructions).
add_gurobi_solver( x, gap = 0.1, time_limit = .Machine$integer.max, presolve = 2, threads = 1, first_feasible = FALSE, numeric_focus = FALSE, node_file_start = Inf, start_solution = NULL, verbose = TRUE )
add_gurobi_solver( x, gap = 0.1, time_limit = .Machine$integer.max, presolve = 2, threads = 1, first_feasible = FALSE, numeric_focus = FALSE, node_file_start = Inf, start_solution = NULL, verbose = TRUE )
x 

gap 

time_limit 

presolve 

threads 

first_feasible 

numeric_focus 

node_file_start 

start_solution 

verbose 

Gurobi is a stateoftheart commercial optimization software with an R package interface. It is by far the fastest of the solvers available for generating prioritizations, however, it is not freely available. That said, licenses are available to academics at no cost. The gurobi package is distributed with the Gurobi software suite. This solver uses the gurobi package to solve problems. For information on the performance of different solvers, please see Schuster et al. (2020) for benchmarks comparing the run time and solution quality of different solvers when applied to different sized datasets.
An updated problem()
object with the solver added to it.
Please see the Gurobi Installation Guide vignette for details on installing the Gurobi software and the gurobi package. You can access this vignette online or using the following code:
vignette("gurobi_installation_guide", package = "prioritizr")
Broadly speaking, the argument to start_solution
must be in the same
format as the planning unit data in the argument to x
.
Further details on the correct format are listed separately
for each of the different planning unit data formats:
x
has numeric
planning unitsThe argument to start_solution
must be a
numeric
vector with each element corresponding to a different planning
unit. It should have the same number of planning units as those
in the argument to x
. Additionally, any planning units missing
cost (NA
) values should also have missing (NA
) values in the
argument to start_solution
.
x
has matrix
planning unitsThe argument to start_solution
must be a
matrix
vector with each row corresponding to a different planning
unit, and each column correspond to a different management zone.
It should have the same number of planning units and zones
as those in the argument to x
. Additionally, any planning units
missing cost (NA
) values for a particular zone should also have a
missing (NA
) values in the argument to start_solution
.
x
has terra::rast()
planning unitsThe argument to start_solution
be a terra::rast()
object where different grid cells (pixels) correspond
to different planning units and layers correspond to
a different management zones. It should have the same dimensionality
(rows, columns, layers), resolution, extent, and coordinate reference
system as the planning units in the argument to x
. Additionally,
any planning units missing cost (NA
) values for a particular zone
should also have missing (NA
) values in the argument to start_solution
.
x
has data.frame
planning unitsThe argument to start_solution
must
be a data.frame
with each column corresponding to a different zone,
each row corresponding to a different planning unit, and cell values
corresponding to the solution value. This means that if a data.frame
object containing the solution also contains additional columns, then
these columns will need to be subsetted prior to using this function
(see below for example with sf::sf()
data).
Additionally, any planning units missing cost
(NA
) values for a particular zone should also have missing (NA
)
values in the argument to start_solution
.
x
has sf::sf()
planning unitsThe argument to start_solution
must be
a sf::sf()
object with each column corresponding to a different
zone, each row corresponding to a different planning unit, and cell values
corresponding to the solution value. This means that if the
sf::sf()
object containing the solution also contains additional
columns, then these columns will need to be subsetted prior to using this
function (see below for example).
Additionally, the argument to start_solution
must also have the same
coordinate reference system as the planning unit data.
Furthermore, any planning units missing cost
(NA
) values for a particular zone should also have missing (NA
)
values in the argument to start_solution
.
Gurobi Optimization LLC (2021) Gurobi Optimizer Reference Manual. https://www.gurobi.com.
Schuster R, Hanson JO, StrimasMackey M, and Bennett JR (2020). Exact integer linear programming solvers outperform simulated annealing for solving conservation planning problems. PeerJ, 8: e9258.
See solvers for an overview of all functions for adding a solver.
Other solvers:
add_cbc_solver()
,
add_cplex_solver()
,
add_default_solver()
,
add_highs_solver()
,
add_lsymphony_solver
,
add_rsymphony_solver()
## Not run: # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() # create problem p1 < problem(sim_pu_raster, sim_features) %>% add_min_set_objective() %>% add_relative_targets(0.1) %>% add_binary_decisions() %>% add_gurobi_solver(gap = 0, verbose = FALSE) # generate solution s1 < solve(p1) # plot solution plot(s1, main = "solution", axes = FALSE) # create a similar problem with boundary length penalties and # specify the solution from the previous run as a starting solution p2 < problem(sim_pu_raster, sim_features) %>% add_min_set_objective() %>% add_relative_targets(0.1) %>% add_boundary_penalties(10) %>% add_binary_decisions() %>% add_gurobi_solver(gap = 0, start_solution = s1, verbose = FALSE) # generate solution s2 < solve(p2) # plot solution plot(s2, main = "solution with boundary penalties", axes = FALSE) ## End(Not run)
## Not run: # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() # create problem p1 < problem(sim_pu_raster, sim_features) %>% add_min_set_objective() %>% add_relative_targets(0.1) %>% add_binary_decisions() %>% add_gurobi_solver(gap = 0, verbose = FALSE) # generate solution s1 < solve(p1) # plot solution plot(s1, main = "solution", axes = FALSE) # create a similar problem with boundary length penalties and # specify the solution from the previous run as a starting solution p2 < problem(sim_pu_raster, sim_features) %>% add_min_set_objective() %>% add_relative_targets(0.1) %>% add_boundary_penalties(10) %>% add_binary_decisions() %>% add_gurobi_solver(gap = 0, start_solution = s1, verbose = FALSE) # generate solution s2 < solve(p2) # plot solution plot(s2, main = "solution with boundary penalties", axes = FALSE) ## End(Not run)
Specify that the HiGHS software should be used to solve a conservation planning problem (Huangfu and Hall 2018). This function can also be used to customize the behavior of the solver. It requires the highs package to be installed.
add_highs_solver( x, gap = 0.1, time_limit = .Machine$integer.max, presolve = TRUE, threads = 1, verbose = TRUE )
add_highs_solver( x, gap = 0.1, time_limit = .Machine$integer.max, presolve = TRUE, threads = 1, verbose = TRUE )
x 

gap 

time_limit 

presolve 

threads 

verbose 

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
the add_gurobi_solver()
, add_cplex_solver()
, add_cbc_solver()
cannot
be used.
An updated problem()
object with the solver added to it.
Huangfu Q and Hall JAJ (2018). Parallelizing the dual revised simplex method. Mathematical Programming Computation, 10: 119142.
Other solvers:
add_cbc_solver()
,
add_cplex_solver()
,
add_default_solver()
,
add_gurobi_solver()
,
add_lsymphony_solver
,
add_rsymphony_solver()
## Not run: # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() # create problem p < problem(sim_pu_raster, sim_features) %>% add_min_set_objective() %>% add_relative_targets(0.1) %>% add_binary_decisions() %>% add_highs_solver(gap = 0, verbose = FALSE) # generate solution s < solve(p) # plot solution plot(s, main = "solution", axes = FALSE) ## End(Not run)
## Not run: # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() # create problem p < problem(sim_pu_raster, sim_features) %>% add_min_set_objective() %>% add_relative_targets(0.1) %>% add_binary_decisions() %>% add_highs_solver(gap = 0, verbose = FALSE) # generate solution s < solve(p) # plot solution plot(s, main = "solution", axes = FALSE) ## End(Not run)
Add constraints to a conservation planning problem to ensure that all selected planning units meet certain criteria.
## S4 method for signature 'ConservationProblem,ANY,ANY,character' add_linear_constraints(x, threshold, sense, data) ## S4 method for signature 'ConservationProblem,ANY,ANY,numeric' add_linear_constraints(x, threshold, sense, data) ## S4 method for signature 'ConservationProblem,ANY,ANY,matrix' add_linear_constraints(x, threshold, sense, data) ## S4 method for signature 'ConservationProblem,ANY,ANY,Matrix' add_linear_constraints(x, threshold, sense, data) ## S4 method for signature 'ConservationProblem,ANY,ANY,Raster' add_linear_constraints(x, threshold, sense, data) ## S4 method for signature 'ConservationProblem,ANY,ANY,SpatRaster' add_linear_constraints(x, threshold, sense, data) ## S4 method for signature 'ConservationProblem,ANY,ANY,dgCMatrix' add_linear_constraints(x, threshold, sense, data)
## S4 method for signature 'ConservationProblem,ANY,ANY,character' add_linear_constraints(x, threshold, sense, data) ## S4 method for signature 'ConservationProblem,ANY,ANY,numeric' add_linear_constraints(x, threshold, sense, data) ## S4 method for signature 'ConservationProblem,ANY,ANY,matrix' add_linear_constraints(x, threshold, sense, data) ## S4 method for signature 'ConservationProblem,ANY,ANY,Matrix' add_linear_constraints(x, threshold, sense, data) ## S4 method for signature 'ConservationProblem,ANY,ANY,Raster' add_linear_constraints(x, threshold, sense, data) ## S4 method for signature 'ConservationProblem,ANY,ANY,SpatRaster' add_linear_constraints(x, threshold, sense, data) ## S4 method for signature 'ConservationProblem,ANY,ANY,dgCMatrix' add_linear_constraints(x, threshold, sense, data)
x 

threshold 

sense 

data 

This function adds general purpose constraints that can be used to
ensure that solutions meet certain criteria
(see Examples section below for details).
For example, these constraints can be used to add multiple budgets.
They can also be used to ensure that the total number of planning units
allocated to a certain administrative area (e.g., country) does not exceed
a certain threshold (e.g., 30% of its total area). Furthermore,
they can also be used to ensure that features have a minimal level
of representation (e.g., 30%) when using an objective
function that aims to enhance feature representation given a budget
(e.g., add_min_shortfall_objective()
).
An updated problem()
object with the constraints added to it.
The linear constraints are implemented using the following
equation.
Let $I$
denote the set of planning units
(indexed by $i$
), $Z$
the set of management zones (indexed by
$z$
), and $X_{iz}$
the decision variable for allocating
planning unit $i$
to zone $z$
(e.g., with binary
values indicating if each planning unit is allocated or not). Also, let
$D_{iz}$
denote the constraint data associated with
planning units $i \in I$
for zones $z \in Z$
(argument to data
, if supplied as a matrix
object),
$\theta$
denote the constraint sense
(argument to sense
, e.g., $<=$
), and $t$
denote the constraint
threshold (argument to threshold
).
$\sum_{i}^{I} \sum_{z}^{Z} (D_{iz} \times X_{iz}) \space \theta \space t$
The argument to data
can be specified using the following formats.
data
as character
vectorcontaining column name(s) that
contain penalty values for planning units. This format is only
compatible if the planning units in the argument to x
are a
sf::sf()
or data.frame
object. The column(s) must have numeric
values, and must not contain any missing (NA
) values.
For problems that contain a single zone, the argument to data
must
contain a single column name. Otherwise, for problems that
contain multiple zones, the argument to data
must
contain a column name for each zone.
data
as a numeric
vectorcontaining values for
planning units. These values must not contain any missing
(NA
) values. Note that this format is only available
for planning units that contain a single zone.
data
as a matrix
/Matrix
objectcontaining numeric
values
that specify data for each planning unit.
Each row corresponds to a planning unit, each column corresponds to a
zone, and each cell indicates the data for penalizing a planning unit
when it is allocated to a given zone.
data
as a terra::rast()
objectcontaining values for planning
units. This format is only
compatible if the planning units in the argument to x
are
sf::sf()
, or terra::rast()
objects.
If the planning unit data are a sf::sf()
object,
then the values are calculated by overlaying the
planning units with the argument to data
and calculating the sum of the
values associated with each planning unit.
If the planning unit data are a terra::rast()
object, then the values
are calculated by extracting the cell
values (note that the planning unit data and the argument to data
must
have exactly the same dimensionality, extent, and missingness).
For problems involving multiple zones, the argument to data
must
contain a layer for each zone.
See constraints for an overview of all functions for adding constraints.
Other constraints:
add_contiguity_constraints()
,
add_feature_contiguity_constraints()
,
add_locked_in_constraints()
,
add_locked_out_constraints()
,
add_mandatory_allocation_constraints()
,
add_manual_bounded_constraints()
,
add_manual_locked_constraints()
,
add_neighbor_constraints()
## Not run: # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() # create a baseline problem with minimum shortfall objective p0 < problem(sim_pu_raster, sim_features) %>% add_min_shortfall_objective(1800) %>% add_relative_targets(0.2) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s0 < solve(p0) # plot solution plot(s0, main = "solution", axes = FALSE) # now let's create some modified versions of this baseline problem by # adding additional criteria using linear constraints # first, let's create a modified version of p0 that contains # an additional budget for a secondary cost dataset # create a secondary cost dataset by simulating values sim_pu_raster2 < simulate_cost(sim_pu_raster) # plot the primary cost dataset (sim_pu_raster) and # the secondary cost dataset (sim_pu_raster2) plot( c(sim_pu_raster, sim_pu_raster2), main = c("sim_pu_raster", "sim_pu_raster2"), axes = FALSE ) # create a modified version of p0 with linear constraints that # specify that the planning units in the solution must not have # values in sim_pu_raster2 that sum to a total greater than 500 p1 < p0 %>% add_linear_constraints( threshold = 500, sense = "<=", data = sim_pu_raster2 ) # solve problem s1 < solve(p1) # plot solutions s1 and s2 to compare them plot(c(s0, s1), main = c("s0", "s1"), axes = FALSE) # second, let's create a modified version of p0 that contains # additional constraints to ensure that each feature definitely has # at least 8% of its overall distribution represented by the solution # (in addition to the 20% targets which specify how much we would # ideally want to conserve for each feature) # to achieve this, we need to calculate the total amount of each feature # within the planning units so we can, in turn, set the constraint thresholds feat_abund < feature_abundances(p0)$absolute_abundance # create a modified version of p0 with additional constraints for each # feature to specify that the planning units in the solution must # secure at least 8% of the total abundance for each feature p2 < p0 for (i in seq_len(terra::nlyr(sim_features))) { p2 < p2 %>% add_linear_constraints( threshold = feat_abund[i] * 0.08, sense = ">=", data = sim_features[[i]] ) } # overall, p2 could be described as an optimization problem # that maximizes feature representation as much as possible # towards securing 20% of the total amount of each feature, # whilst ensuring that (i) the total cost of the solution does # not exceed 1800 (per cost values in sim_pu_raster) and (ii) # the solution secures at least 8% of the total amount of each feature # (if 20% is not possible due to the budget) # solve problem s2 < solve(p2) # plot solutions s0 and s2 to compare them plot(c(s0, s2), main = c("s1", "s2"), axes = FALSE) # third, let's create a modified version of p0 that contains # additional constraints to ensure that the solution equitably # distributes conservation effort across different administrative areas # (e.g., countries) within the study region # to begin with, we will simulate a dataset describing the spatial extent of # four different administrative areas across the study region sim_admin < sim_pu_raster sim_admin < terra::aggregate(sim_admin, fact = 5) sim_admin < terra::setValues(sim_admin, seq_len(terra::ncell(sim_admin))) sim_admin < terra::resample(sim_admin, sim_pu_raster, method = "near") sim_admin < terra::mask(sim_admin, sim_pu_raster) # plot administrative areas layer, # we can see that the administrative areas subdivide # the study region into four quadrants, and the sim_admin object is a # SpatRaster with integer values denoting ids for the administrative areas plot(sim_admin, axes = FALSE) # next we will convert the sim_admin SpatRaster object into a SpatRaster # object (with a layer for each administrative area) indicating which # planning units belong to each administrative area using binary # (presence/absence) values sim_admin2 < binary_stack(sim_admin) # plot binary stack of administrative areas plot(sim_admin2, axes = FALSE) # we will now calculate the total amount of planning units associated # with each administrative area, so that we can set the constraint threshold # since we are using raster data, we won't bother explicitly # accounting for the total area of each planning unit (because all # planning units have the same area in raster formats)  but if we were # using vector data then we would need to account for the area of each unit admin_total < Matrix::rowSums(rij_matrix(sim_pu_raster, sim_admin2)) # create a modified version of p0 with additional constraints for each # administrative area to specify that the planning units in the solution must # not encompass more than 10% of the total extent of the administrative # area p3 < p0 for (i in seq_len(terra::nlyr(sim_admin2))) { p3 < p3 %>% add_linear_constraints( threshold = admin_total[i] * 0.1, sense = "<=", data = sim_admin2[[i]] ) } # solve problem s3 < solve(p3) # plot solutions s0 and s3 to compare them plot(c(s0, s3), main = c("s0", "s3"), axes = FALSE) ## End(Not run)
## Not run: # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() # create a baseline problem with minimum shortfall objective p0 < problem(sim_pu_raster, sim_features) %>% add_min_shortfall_objective(1800) %>% add_relative_targets(0.2) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s0 < solve(p0) # plot solution plot(s0, main = "solution", axes = FALSE) # now let's create some modified versions of this baseline problem by # adding additional criteria using linear constraints # first, let's create a modified version of p0 that contains # an additional budget for a secondary cost dataset # create a secondary cost dataset by simulating values sim_pu_raster2 < simulate_cost(sim_pu_raster) # plot the primary cost dataset (sim_pu_raster) and # the secondary cost dataset (sim_pu_raster2) plot( c(sim_pu_raster, sim_pu_raster2), main = c("sim_pu_raster", "sim_pu_raster2"), axes = FALSE ) # create a modified version of p0 with linear constraints that # specify that the planning units in the solution must not have # values in sim_pu_raster2 that sum to a total greater than 500 p1 < p0 %>% add_linear_constraints( threshold = 500, sense = "<=", data = sim_pu_raster2 ) # solve problem s1 < solve(p1) # plot solutions s1 and s2 to compare them plot(c(s0, s1), main = c("s0", "s1"), axes = FALSE) # second, let's create a modified version of p0 that contains # additional constraints to ensure that each feature definitely has # at least 8% of its overall distribution represented by the solution # (in addition to the 20% targets which specify how much we would # ideally want to conserve for each feature) # to achieve this, we need to calculate the total amount of each feature # within the planning units so we can, in turn, set the constraint thresholds feat_abund < feature_abundances(p0)$absolute_abundance # create a modified version of p0 with additional constraints for each # feature to specify that the planning units in the solution must # secure at least 8% of the total abundance for each feature p2 < p0 for (i in seq_len(terra::nlyr(sim_features))) { p2 < p2 %>% add_linear_constraints( threshold = feat_abund[i] * 0.08, sense = ">=", data = sim_features[[i]] ) } # overall, p2 could be described as an optimization problem # that maximizes feature representation as much as possible # towards securing 20% of the total amount of each feature, # whilst ensuring that (i) the total cost of the solution does # not exceed 1800 (per cost values in sim_pu_raster) and (ii) # the solution secures at least 8% of the total amount of each feature # (if 20% is not possible due to the budget) # solve problem s2 < solve(p2) # plot solutions s0 and s2 to compare them plot(c(s0, s2), main = c("s1", "s2"), axes = FALSE) # third, let's create a modified version of p0 that contains # additional constraints to ensure that the solution equitably # distributes conservation effort across different administrative areas # (e.g., countries) within the study region # to begin with, we will simulate a dataset describing the spatial extent of # four different administrative areas across the study region sim_admin < sim_pu_raster sim_admin < terra::aggregate(sim_admin, fact = 5) sim_admin < terra::setValues(sim_admin, seq_len(terra::ncell(sim_admin))) sim_admin < terra::resample(sim_admin, sim_pu_raster, method = "near") sim_admin < terra::mask(sim_admin, sim_pu_raster) # plot administrative areas layer, # we can see that the administrative areas subdivide # the study region into four quadrants, and the sim_admin object is a # SpatRaster with integer values denoting ids for the administrative areas plot(sim_admin, axes = FALSE) # next we will convert the sim_admin SpatRaster object into a SpatRaster # object (with a layer for each administrative area) indicating which # planning units belong to each administrative area using binary # (presence/absence) values sim_admin2 < binary_stack(sim_admin) # plot binary stack of administrative areas plot(sim_admin2, axes = FALSE) # we will now calculate the total amount of planning units associated # with each administrative area, so that we can set the constraint threshold # since we are using raster data, we won't bother explicitly # accounting for the total area of each planning unit (because all # planning units have the same area in raster formats)  but if we were # using vector data then we would need to account for the area of each unit admin_total < Matrix::rowSums(rij_matrix(sim_pu_raster, sim_admin2)) # create a modified version of p0 with additional constraints for each # administrative area to specify that the planning units in the solution must # not encompass more than 10% of the total extent of the administrative # area p3 < p0 for (i in seq_len(terra::nlyr(sim_admin2))) { p3 < p3 %>% add_linear_constraints( threshold = admin_total[i] * 0.1, sense = "<=", data = sim_admin2[[i]] ) } # solve problem s3 < solve(p3) # plot solutions s0 and s3 to compare them plot(c(s0, s3), main = c("s0", "s3"), axes = FALSE) ## End(Not run)
Add penalties to a conservation planning problem to penalize
solutions that select planning units with higher values from a specific
data source (e.g., anthropogenic impact). These penalties assume
a linear tradeoff between the penalty values and the primary
objective of the conservation planning problem (e.g.,
solution cost for minimum set problems; add_min_set_objective()
.
## S4 method for signature 'ConservationProblem,ANY,character' add_linear_penalties(x, penalty, data) ## S4 method for signature 'ConservationProblem,ANY,numeric' add_linear_penalties(x, penalty, data) ## S4 method for signature 'ConservationProblem,ANY,matrix' add_linear_penalties(x, penalty, data) ## S4 method for signature 'ConservationProblem,ANY,Matrix' add_linear_penalties(x, penalty, data) ## S4 method for signature 'ConservationProblem,ANY,Raster' add_linear_penalties(x, penalty, data) ## S4 method for signature 'ConservationProblem,ANY,SpatRaster' add_linear_penalties(x, penalty, data) ## S4 method for signature 'ConservationProblem,ANY,dgCMatrix' add_linear_penalties(x, penalty, data)
## S4 method for signature 'ConservationProblem,ANY,character' add_linear_penalties(x, penalty, data) ## S4 method for signature 'ConservationProblem,ANY,numeric' add_linear_penalties(x, penalty, data) ## S4 method for signature 'ConservationProblem,ANY,matrix' add_linear_penalties(x, penalty, data) ## S4 method for signature 'ConservationProblem,ANY,Matrix' add_linear_penalties(x, penalty, data) ## S4 method for signature 'ConservationProblem,ANY,Raster' add_linear_penalties(x, penalty, data) ## S4 method for signature 'ConservationProblem,ANY,SpatRaster' add_linear_penalties(x, penalty, data) ## S4 method for signature 'ConservationProblem,ANY,dgCMatrix' add_linear_penalties(x, penalty, data)
x 

penalty 

data 

This function penalizes solutions that have higher values according to the sum of the penalty values associated with each planning unit, weighted by status of each planning unit in the solution.
An updated problem()
object with the penalties added to it.
The argument to data
can be specified using the following formats.
data
as character
vectorcontaining column name(s) that
contain penalty values for planning units. This format is only
compatible if the planning units in the argument to x
are a
sf::sf()
or data.frame
object. The column(s) must have numeric
values, and must not contain any missing (NA
) values.
For problems that contain a single zone, the argument to data
must
contain a single column name. Otherwise, for problems that
contain multiple zones, the argument to data
must
contain a column name for each zone.
data
as a numeric
vectorcontaining values for
planning units. These values must not contain any missing
(NA
) values. Note that this format is only available
for planning units that contain a single zone.
data
as a matrix
/Matrix
objectcontaining numeric
values
that specify data for each planning unit.
Each row corresponds to a planning unit, each column corresponds to a
zone, and each cell indicates the data for penalizing a planning unit
when it is allocated to a given zone.
data
as a terra::rast()
objectcontaining values for planning
units. This format is only
compatible if the planning units in the argument to x
are
sf::sf()
, or terra::rast()
objects.
If the planning unit data are a sf::sf()
object,
then the values are calculated by overlaying the
planning units with the argument to data
and calculating the sum of the
values associated with each planning unit.
If the planning unit data are a terra::rast()
object, then the values
are calculated by extracting the cell
values (note that the planning unit data and the argument to data
must
have exactly the same dimensionality, extent, and missingness).
For problems involving multiple zones, the argument to data
must
contain a layer for each zone.
The linear penalties are implemented using the following
equations.
Let $I$
denote the set of planning units
(indexed by $i$
), $Z$
the set of management zones (indexed by
$z$
), and $X_{iz}$
the decision variable for allocating
planning unit $i$
to zone $z$
(e.g., with binary
values indicating if each planning unit is allocated or not). Also, let
$P_z$
represent the penalty scaling value for zones
$z \in Z$
(argument to penalty
), and
$D_{iz}$
the penalty data for allocating planning unit
$i \in I$
to zones $z \in Z$
(argument to
data
, if supplied as a matrix
object).
$\sum_{i}^{I} \sum_{z}^{Z} P_z \times D_{iz} \times X_{iz}$
Note that when the problem objective is to maximize some measure of
benefit and not minimize some measure of cost, the term $P_z$
is
replaced with $P_z$
.
See penalties for an overview of all functions for adding penalties.
Other penalties:
add_asym_connectivity_penalties()
,
add_boundary_penalties()
,
add_connectivity_penalties()
,
add_feature_weights()
## Not run: # set seed for reproducibility set.seed(600) # load data sim_pu_polygons < get_sim_pu_polygons() sim_features < get_sim_features() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_features < get_sim_zones_features() # add a column to contain the penalty data for each planning unit # e.g., these values could indicate the level of habitat sim_pu_polygons$penalty_data < runif(nrow(sim_pu_polygons)) # plot the penalty data to visualise its spatial distribution plot(sim_pu_polygons[, "penalty_data"], axes = FALSE) # create minimal problem with minimum set objective, # this does not use the penalty data p1 < problem(sim_pu_polygons, sim_features, cost_column = "cost") %>% add_min_set_objective() %>% add_relative_targets(0.1) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # print problem print(p1) # create an updated version of the previous problem, # with the penalties added to it p2 < p1 %>% add_linear_penalties(100, data = "penalty_data") # print problem print(p2) # solve the two problems s1 < solve(p1) s2 < solve(p2) # create a new object with both solutions s3 < sf::st_sf( tibble::tibble( s1 = s1$solution_1, s2 = s2$solution_1 ), geometry = sf::st_geometry(s1) ) # plot the solutions and compare them, # since we supplied a very high penalty value (i.e., 100), relative # to the range of values in the penalty data and the objective function, # the solution in s2 is very sensitive to values in the penalty data plot(s3, axes = FALSE) # for real conservation planning exercises, # it would be worth exploring a range of penalty values (e.g., ranging # from 1 to 100 increments of 5) to explore the tradeoffs # now, let's examine a conservation planning exercise involving multiple # management zones # create targets for each feature within each zone, # these targets indicate that each zone needs to represent 10% of the # spatial distribution of each feature targ < matrix( 0.1, ncol = number_of_zones(sim_zones_features), nrow = number_of_features(sim_zones_features) ) # create penalty data for allocating each planning unit to each zone, # these data will be generated by simulating values penalty_raster < simulate_cost( sim_zones_pu_raster[[1]], n = number_of_zones(sim_zones_features) ) # plot the penalty data, each layer corresponds to a different zone plot(penalty_raster, main = "penalty data", axes = FALSE) # create a multizone problem with the minimum set objective # and penalties for allocating planning units to each zone, # with a penalty scaling factor of 1 for each zone p4 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_min_set_objective() %>% add_relative_targets(targ) %>% add_linear_penalties(c(1, 1, 1), penalty_raster) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # print problem print(p4) # solve problem s4 < solve(p4) # plot solution plot(category_layer(s4), main = "multizone solution", axes = FALSE) ## End(Not run)
## Not run: # set seed for reproducibility set.seed(600) # load data sim_pu_polygons < get_sim_pu_polygons() sim_features < get_sim_features() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_features < get_sim_zones_features() # add a column to contain the penalty data for each planning unit # e.g., these values could indicate the level of habitat sim_pu_polygons$penalty_data < runif(nrow(sim_pu_polygons)) # plot the penalty data to visualise its spatial distribution plot(sim_pu_polygons[, "penalty_data"], axes = FALSE) # create minimal problem with minimum set objective, # this does not use the penalty data p1 < problem(sim_pu_polygons, sim_features, cost_column = "cost") %>% add_min_set_objective() %>% add_relative_targets(0.1) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # print problem print(p1) # create an updated version of the previous problem, # with the penalties added to it p2 < p1 %>% add_linear_penalties(100, data = "penalty_data") # print problem print(p2) # solve the two problems s1 < solve(p1) s2 < solve(p2) # create a new object with both solutions s3 < sf::st_sf( tibble::tibble( s1 = s1$solution_1, s2 = s2$solution_1 ), geometry = sf::st_geometry(s1) ) # plot the solutions and compare them, # since we supplied a very high penalty value (i.e., 100), relative # to the range of values in the penalty data and the objective function, # the solution in s2 is very sensitive to values in the penalty data plot(s3, axes = FALSE) # for real conservation planning exercises, # it would be worth exploring a range of penalty values (e.g., ranging # from 1 to 100 increments of 5) to explore the tradeoffs # now, let's examine a conservation planning exercise involving multiple # management zones # create targets for each feature within each zone, # these targets indicate that each zone needs to represent 10% of the # spatial distribution of each feature targ < matrix( 0.1, ncol = number_of_zones(sim_zones_features), nrow = number_of_features(sim_zones_features) ) # create penalty data for allocating each planning unit to each zone, # these data will be generated by simulating values penalty_raster < simulate_cost( sim_zones_pu_raster[[1]], n = number_of_zones(sim_zones_features) ) # plot the penalty data, each layer corresponds to a different zone plot(penalty_raster, main = "penalty data", axes = FALSE) # create a multizone problem with the minimum set objective # and penalties for allocating planning units to each zone, # with a penalty scaling factor of 1 for each zone p4 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_min_set_objective() %>% add_relative_targets(targ) %>% add_linear_penalties(c(1, 1, 1), penalty_raster) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # print problem print(p4) # solve problem s4 < solve(p4) # plot solution plot(category_layer(s4), main = "multizone solution", axes = FALSE) ## End(Not run)
Add constraints to a conservation planning problem to ensure
that specific planning units are selected (or allocated
to a specific zone) 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. If specific
planning units should be locked out of a solution, use
add_locked_out_constraints()
. For problems with nonbinary
planning unit allocations (e.g., proportions), the
add_manual_locked_constraints()
function can be used to lock
planning unit allocations to a specific value.
add_locked_in_constraints(x, locked_in) ## S4 method for signature 'ConservationProblem,numeric' add_locked_in_constraints(x, locked_in) ## S4 method for signature 'ConservationProblem,logical' add_locked_in_constraints(x, locked_in) ## S4 method for signature 'ConservationProblem,matrix' add_locked_in_constraints(x, locked_in) ## S4 method for signature 'ConservationProblem,character' add_locked_in_constraints(x, locked_in) ## S4 method for signature 'ConservationProblem,Spatial' add_locked_in_constraints(x, locked_in) ## S4 method for signature 'ConservationProblem,sf' add_locked_in_constraints(x, locked_in) ## S4 method for signature 'ConservationProblem,Raster' add_locked_in_constraints(x, locked_in) ## S4 method for signature 'ConservationProblem,SpatRaster' add_locked_in_constraints(x, locked_in)
add_locked_in_constraints(x, locked_in) ## S4 method for signature 'ConservationProblem,numeric' add_locked_in_constraints(x, locked_in) ## S4 method for signature 'ConservationProblem,logical' add_locked_in_constraints(x, locked_in) ## S4 method for signature 'ConservationProblem,matrix' add_locked_in_constraints(x, locked_in) ## S4 method for signature 'ConservationProblem,character' add_locked_in_constraints(x, locked_in) ## S4 method for signature 'ConservationProblem,Spatial' add_locked_in_constraints(x, locked_in) ## S4 method for signature 'ConservationProblem,sf' add_locked_in_constraints(x, locked_in) ## S4 method for signature 'ConservationProblem,Raster' add_locked_in_constraints(x, locked_in) ## S4 method for signature 'ConservationProblem,SpatRaster' add_locked_in_constraints(x, locked_in)
x 

locked_in 
Object that determines which planning units should be locked in. See the Data format section for more information. 
An updated problem()
object with the constraints added to it.
The locked planning units can be specified using the following formats.
Generally, the locked data should correspond to the planning units
in the argument to x.
To help make working with
terra::rast()
planning unit data easier,
the locked data should correspond to cell indices in the
terra::rast()
data. For example, integer
arguments
should correspond to cell indices and logical
arguments should have
a value for each cell—regardless of which planning unit cells contain
NA
values.
data
as an integer
vectorcontaining indices that indicate which planning units should be locked for the solution. This argument is only compatible with problems that contain a single zone.
data
as a logical
vectorcontaining TRUE
and/or
FALSE
values that indicate which planning units should be locked
in the solution. This argument is only compatible with problems that
contain a single zone.
data
as a matrix
objectcontaining logical
TRUE
and/or
FALSE
values which indicate if certain planning units are
should be locked to a specific zone in the solution. Each row
corresponds to a planning unit, each column corresponds to a zone, and
each cell indicates if the planning unit should be locked to a given
zone. Thus each row should only contain at most a single TRUE
value.
data
as a character
vectorcontaining column name(s)
that indicates if planning units should be locked for the solution.
This format is only
compatible if the planning units in the argument to x
are a
sf::sf()
or data.frame
object. The columns
must have logical
(i.e., TRUE
or FALSE
)
values indicating if the planning unit is to be locked for the solution.
For problems that contain a single zone, the argument to data
must
contain a single column name. Otherwise, for problems that
contain multiple zones, the argument to data
must
contain a column name for each zone.
data
as a sf::sf()
objectcontaining geometries that will be used to lock planning units for
the solution. Specifically, planning units in x
that spatially
intersect with y
will be locked (per intersecting_units()
).
Note that this option is only available
for problems that contain a single management zone.
data
as a terra::rast()
objectcontaining cells used to lock planning units for the solution.
Specifically, planning units in x
that intersect with cells that have nonzero and nonNA
values are
locked.
For problems that contain multiple zones, the
data
object must contain a layer
for each zone. Note that for multiband arguments, each pixel must
only contain a nonzero value in a single band. Additionally, if the
cost data in x
is a terra::rast()
object, we
recommend standardizing NA
values in this dataset with the cost
data. In other words, the pixels in x
that have NA
values
should also have NA
values in the locked data.
See constraints for an overview of all functions for adding constraints.
Other constraints:
add_contiguity_constraints()
,
add_feature_contiguity_constraints()
,
add_linear_constraints()
,
add_locked_out_constraints()
,
add_mandatory_allocation_constraints()
,
add_manual_bounded_constraints()
,
add_manual_locked_constraints()
,
add_neighbor_constraints()
## Not run: # set seed for reproducibility set.seed(500) # load data sim_pu_polygons < get_sim_pu_polygons() sim_features < get_sim_features() sim_locked_in_raster < get_sim_locked_in_raster() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_pu_polygons < get_sim_zones_pu_polygons() sim_zones_features < get_sim_zones_features() # create minimal problem p1 < problem(sim_pu_polygons, sim_features, "cost") %>% add_min_set_objective() %>% add_relative_targets(0.2) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # create problem with added locked in constraints using integers p2 < p1 %>% add_locked_in_constraints(which(sim_pu_polygons$locked_in)) # create problem with added locked in constraints using a column name p3 < p1 %>% add_locked_in_constraints("locked_in") # create problem with added locked in constraints using raster data p4 < p1 %>% add_locked_in_constraints(sim_locked_in_raster) # create problem with added locked in constraints using spatial polygon data locked_in < sim_pu_polygons[sim_pu_polygons$locked_in == 1, ] p5 < p1 %>% add_locked_in_constraints(locked_in) # solve problems s1 < solve(p1) s2 < solve(p2) s3 < solve(p3) s4 < solve(p4) s5 < solve(p5) # create single object with all solutions s6 < sf::st_sf( tibble::tibble( s1 = s1$solution_1, s2 = s2$solution_1, s3 = s3$solution_1, s4 = s4$solution_1, s5 = s5$solution_1 ), geometry = sf::st_geometry(s1) ) # plot solutions plot( s6, main = c( "none locked in", "locked in (integer input)", "locked in (character input)", "locked in (raster input)", "locked in (polygon input)" ) ) # create minimal multizone problem with spatial data p7 < problem( sim_zones_pu_polygons, sim_zones_features, cost_column = c("cost_1", "cost_2", "cost_3") ) %>% add_min_set_objective() %>% add_absolute_targets(matrix(rpois(15, 1), nrow = 5, ncol = 3)) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # create multizone problem with locked in constraints using matrix data locked_matrix < as.matrix(sf::st_drop_geometry( sim_zones_pu_polygons[, c("locked_1", "locked_2", "locked_3")] )) p8 < p7 %>% add_locked_in_constraints(locked_matrix) # solve problem s8 < solve(p8) # create new column representing the zone id that each planning unit # was allocated to in the solution s8$solution < category_vector(sf::st_drop_geometry( s8[, c("solution_1_zone_1", "solution_1_zone_2", "solution_1_zone_3")] )) s8$solution < factor(s8$solution) # plot solution plot(s8[ "solution"], axes = FALSE) # create multizone problem with locked in constraints using column names p9 < p7 %>% add_locked_in_constraints(c("locked_1", "locked_2", "locked_3")) # solve problem s9 < solve(p9) # create new column representing the zone id that each planning unit # was allocated to in the solution s9$solution < category_vector(sf::st_drop_geometry( s9[, c("solution_1_zone_1", "solution_1_zone_2", "solution_1_zone_3")] )) s9$solution[s9$solution == 1 & s9$solution_1_zone_1 == 0] < 0 s9$solution < factor(s9$solution) # plot solution plot(s9[, "solution"], axes = FALSE) # create multizone problem with raster planning units p10 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_min_set_objective() %>% add_absolute_targets(matrix(rpois(15, 1), nrow = 5, ncol = 3)) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # create multilayer raster with locked in units locked_in_raster < sim_zones_pu_raster[[1]] locked_in_raster[!is.na(locked_in_raster)] < 0 locked_in_raster < locked_in_raster[[c(1, 1, 1)]] names(locked_in_raster) < c("zone_1", "zone_2", "zone_3") locked_in_raster[[1]][1] < 1 locked_in_raster[[2]][2] < 1 locked_in_raster[[3]][3] < 1 # plot locked in raster plot(locked_in_raster) # add locked in raster units to problem p10 < p10 %>% add_locked_in_constraints(locked_in_raster) # solve problem s10 < solve(p10) # plot solution plot(category_layer(s10), main = "solution", axes = FALSE) ## End(Not run)
## Not run: # set seed for reproducibility set.seed(500) # load data sim_pu_polygons < get_sim_pu_polygons() sim_features < get_sim_features() sim_locked_in_raster < get_sim_locked_in_raster() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_pu_polygons < get_sim_zones_pu_polygons() sim_zones_features < get_sim_zones_features() # create minimal problem p1 < problem(sim_pu_polygons, sim_features, "cost") %>% add_min_set_objective() %>% add_relative_targets(0.2) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # create problem with added locked in constraints using integers p2 < p1 %>% add_locked_in_constraints(which(sim_pu_polygons$locked_in)) # create problem with added locked in constraints using a column name p3 < p1 %>% add_locked_in_constraints("locked_in") # create problem with added locked in constraints using raster data p4 < p1 %>% add_locked_in_constraints(sim_locked_in_raster) # create problem with added locked in constraints using spatial polygon data locked_in < sim_pu_polygons[sim_pu_polygons$locked_in == 1, ] p5 < p1 %>% add_locked_in_constraints(locked_in) # solve problems s1 < solve(p1) s2 < solve(p2) s3 < solve(p3) s4 < solve(p4) s5 < solve(p5) # create single object with all solutions s6 < sf::st_sf( tibble::tibble( s1 = s1$solution_1, s2 = s2$solution_1, s3 = s3$solution_1, s4 = s4$solution_1, s5 = s5$solution_1 ), geometry = sf::st_geometry(s1) ) # plot solutions plot( s6, main = c( "none locked in", "locked in (integer input)", "locked in (character input)", "locked in (raster input)", "locked in (polygon input)" ) ) # create minimal multizone problem with spatial data p7 < problem( sim_zones_pu_polygons, sim_zones_features, cost_column = c("cost_1", "cost_2", "cost_3") ) %>% add_min_set_objective() %>% add_absolute_targets(matrix(rpois(15, 1), nrow = 5, ncol = 3)) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # create multizone problem with locked in constraints using matrix data locked_matrix < as.matrix(sf::st_drop_geometry( sim_zones_pu_polygons[, c("locked_1", "locked_2", "locked_3")] )) p8 < p7 %>% add_locked_in_constraints(locked_matrix) # solve problem s8 < solve(p8) # create new column representing the zone id that each planning unit # was allocated to in the solution s8$solution < category_vector(sf::st_drop_geometry( s8[, c("solution_1_zone_1", "solution_1_zone_2", "solution_1_zone_3")] )) s8$solution < factor(s8$solution) # plot solution plot(s8[ "solution"], axes = FALSE) # create multizone problem with locked in constraints using column names p9 < p7 %>% add_locked_in_constraints(c("locked_1", "locked_2", "locked_3")) # solve problem s9 < solve(p9) # create new column representing the zone id that each planning unit # was allocated to in the solution s9$solution < category_vector(sf::st_drop_geometry( s9[, c("solution_1_zone_1", "solution_1_zone_2", "solution_1_zone_3")] )) s9$solution[s9$solution == 1 & s9$solution_1_zone_1 == 0] < 0 s9$solution < factor(s9$solution) # plot solution plot(s9[, "solution"], axes = FALSE) # create multizone problem with raster planning units p10 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_min_set_objective() %>% add_absolute_targets(matrix(rpois(15, 1), nrow = 5, ncol = 3)) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # create multilayer raster with locked in units locked_in_raster < sim_zones_pu_raster[[1]] locked_in_raster[!is.na(locked_in_raster)] < 0 locked_in_raster < locked_in_raster[[c(1, 1, 1)]] names(locked_in_raster) < c("zone_1", "zone_2", "zone_3") locked_in_raster[[1]][1] < 1 locked_in_raster[[2]][2] < 1 locked_in_raster[[3]][3] < 1 # plot locked in raster plot(locked_in_raster) # add locked in raster units to problem p10 < p10 %>% add_locked_in_constraints(locked_in_raster) # solve problem s10 < solve(p10) # plot solution plot(category_layer(s10), main = "solution", axes = FALSE) ## End(Not run)
Add constraints to a conservation planning problem to ensure
that specific planning units are not selected
(or allocated to a specific zone) 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. If specific planning units should be locked
in to the solution, use add_locked_in_constraints()
. For
problems with nonbinary planning unit allocations (e.g., proportions), the
add_manual_locked_constraints()
function can be used to lock
planning unit allocations to a specific value.
add_locked_out_constraints(x, locked_out) ## S4 method for signature 'ConservationProblem,numeric' add_locked_out_constraints(x, locked_out) ## S4 method for signature 'ConservationProblem,logical' add_locked_out_constraints(x, locked_out) ## S4 method for signature 'ConservationProblem,matrix' add_locked_out_constraints(x, locked_out) ## S4 method for signature 'ConservationProblem,character' add_locked_out_constraints(x, locked_out) ## S4 method for signature 'ConservationProblem,Spatial' add_locked_out_constraints(x, locked_out) ## S4 method for signature 'ConservationProblem,sf' add_locked_out_constraints(x, locked_out) ## S4 method for signature 'ConservationProblem,Raster' add_locked_out_constraints(x, locked_out) ## S4 method for signature 'ConservationProblem,SpatRaster' add_locked_out_constraints(x, locked_out)
add_locked_out_constraints(x, locked_out) ## S4 method for signature 'ConservationProblem,numeric' add_locked_out_constraints(x, locked_out) ## S4 method for signature 'ConservationProblem,logical' add_locked_out_constraints(x, locked_out) ## S4 method for signature 'ConservationProblem,matrix' add_locked_out_constraints(x, locked_out) ## S4 method for signature 'ConservationProblem,character' add_locked_out_constraints(x, locked_out) ## S4 method for signature 'ConservationProblem,Spatial' add_locked_out_constraints(x, locked_out) ## S4 method for signature 'ConservationProblem,sf' add_locked_out_constraints(x, locked_out) ## S4 method for signature 'ConservationProblem,Raster' add_locked_out_constraints(x, locked_out) ## S4 method for signature 'ConservationProblem,SpatRaster' add_locked_out_constraints(x, locked_out)
x 

locked_out 
Object that determines which planning units that should be locked out. See the Data format section for more information. 
An updated problem()
object with the constraints added to it.
The locked planning units can be specified using the following formats.
Generally, the locked data should correspond to the planning units
in the argument to x.
To help make working with
terra::rast()
planning unit data easier,
the locked data should correspond to cell indices in the
terra::rast()
data. For example, integer
arguments
should correspond to cell indices and logical
arguments should have
a value for each cell—regardless of which planning unit cells contain
NA
values.
data
as an integer
vectorcontaining indices that indicate which planning units should be locked for the solution. This argument is only compatible with problems that contain a single zone.
data
as a logical
vectorcontaining TRUE
and/or
FALSE
values that indicate which planning units should be locked
in the solution. This argument is only compatible with problems that
contain a single zone.
data
as a matrix
objectcontaining logical
TRUE
and/or
FALSE
values which indicate if certain planning units are
should be locked to a specific zone in the solution. Each row
corresponds to a planning unit, each column corresponds to a zone, and
each cell indicates if the planning unit should be locked to a given
zone. Thus each row should only contain at most a single TRUE
value.
data
as a character
vectorcontaining column name(s)
that indicates if planning units should be locked for the solution.
This format is only
compatible if the planning units in the argument to x
are a
sf::sf()
or data.frame
object. The columns
must have logical
(i.e., TRUE
or FALSE
)
values indicating if the planning unit is to be locked for the solution.
For problems that contain a single zone, the argument to data
must
contain a single column name. Otherwise, for problems that
contain multiple zones, the argument to data
must
contain a column name for each zone.
data
as a sf::sf()
objectcontaining geometries that will be used to lock planning units for
the solution. Specifically, planning units in x
that spatially
intersect with y
will be locked (per intersecting_units()
).
Note that this option is only available
for problems that contain a single management zone.
data
as a terra::rast()
objectcontaining cells used to lock planning units for the solution.
Specifically, planning units in x
that intersect with cells that have nonzero and nonNA
values are
locked.
For problems that contain multiple zones, the
data
object must contain a layer
for each zone. Note that for multiband arguments, each pixel must
only contain a nonzero value in a single band. Additionally, if the
cost data in x
is a terra::rast()
object, we
recommend standardizing NA
values in this dataset with the cost
data. In other words, the pixels in x
that have NA
values
should also have NA
values in the locked data.
See constraints for an overview of all functions for adding constraints.
Other constraints:
add_contiguity_constraints()
,
add_feature_contiguity_constraints()
,
add_linear_constraints()
,
add_locked_in_constraints()
,
add_mandatory_allocation_constraints()
,
add_manual_bounded_constraints()
,
add_manual_locked_constraints()
,
add_neighbor_constraints()
## Not run: # set seed for reproducibility set.seed(500) # load data sim_pu_polygons < get_sim_pu_polygons() sim_features < get_sim_features() sim_locked_out_raster < get_sim_locked_out_raster() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_pu_polygons < get_sim_zones_pu_polygons() sim_zones_features < get_sim_zones_features() # create minimal problem p1 < problem(sim_pu_polygons, sim_features, "cost") %>% add_min_set_objective() %>% add_relative_targets(0.2) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # create problem with added locked out constraints using integers p2 < p1 %>% add_locked_out_constraints(which(sim_pu_polygons$locked_out)) # create problem with added locked out constraints using a column name p3 < p1 %>% add_locked_out_constraints("locked_out") # create problem with added locked out constraints using raster data p4 < p1 %>% add_locked_out_constraints(sim_locked_out_raster) # create problem with added locked out constraints using spatial polygon data locked_out < sim_pu_polygons[sim_pu_polygons$locked_out == 1, ] p5 < p1 %>% add_locked_out_constraints(locked_out) # solve problems s1 < solve(p1) s2 < solve(p2) s3 < solve(p3) s4 < solve(p4) s5 < solve(p5) # create single object with all solutions s6 < sf::st_sf( tibble::tibble( s1 = s1$solution_1, s2 = s2$solution_1, s3 = s3$solution_1, s4 = s4$solution_1, s5 = s5$solution_1 ), geometry = sf::st_geometry(s1) ) # plot solutions plot( s6, main = c( "none locked out", "locked out (integer input)", "locked out (character input)", "locked out (raster input)", "locked out (polygon input)" ) ) # reset plot par(mfrow = c(1, 1)) # create minimal multizone problem with spatial data p7 < problem( sim_zones_pu_polygons, sim_zones_features, cost_column = c("cost_1", "cost_2", "cost_3") ) %>% add_min_set_objective() %>% add_absolute_targets(matrix(rpois(15, 1), nrow = 5, ncol = 3)) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # create multizone problem with locked out constraints using matrix data locked_matrix < as.matrix(sf::st_drop_geometry( sim_zones_pu_polygons[, c("locked_1", "locked_2", "locked_3")] )) p8 < p7 %>% add_locked_out_constraints(locked_matrix) # solve problem s8 < solve(p8) # create new column representing the zone id that each planning unit # was allocated to in the solution s8$solution < category_vector(sf::st_drop_geometry( s8[, c("solution_1_zone_1", "solution_1_zone_2", "solution_1_zone_3")] )) s8$solution < factor(s8$solution) # plot solution plot(s8[, "solution"], main = "solution", axes = FALSE) # create multizone problem with locked out constraints using column names p9 < p7 %>% add_locked_out_constraints(c("locked_1", "locked_2", "locked_3")) # solve problem s9 < solve(p9) # create new column in s8 representing the zone id that each planning unit # was allocated to in the solution s9$solution < category_vector(sf::st_drop_geometry( s9[, c("solution_1_zone_1", "solution_1_zone_2", "solution_1_zone_3")] )) s9$solution[s9$solution == 1 & s9$solution_1_zone_1 == 0] < 0 s9$solution < factor(s9$solution) # plot solution plot(s9[, "solution"], main = "solution", axes = FALSE) # create multizone problem with raster planning units p10 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_min_set_objective() %>% add_absolute_targets(matrix(rpois(15, 1), nrow = 5, ncol = 3)) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # create multilayer raster with locked out units locked_out_raster < sim_zones_pu_raster[[1]] locked_out_raster[!is.na(locked_out_raster)] < 0 locked_out_raster < locked_out_raster[[c(1, 1, 1)]] names(locked_out_raster) < c("zones_1", "zones_2", "zones_3") locked_out_raster[[1]][1] < 1 locked_out_raster[[2]][2] < 1 locked_out_raster[[3]][3] < 1 # plot locked out raster plot(locked_out_raster) # add locked out raster units to problem p10 < p10 %>% add_locked_out_constraints(locked_out_raster) # solve problem s10 < solve(p10) # plot solution plot(category_layer(s10), main = "solution", axes = FALSE) ## End(Not run)
## Not run: # set seed for reproducibility set.seed(500) # load data sim_pu_polygons < get_sim_pu_polygons() sim_features < get_sim_features() sim_locked_out_raster < get_sim_locked_out_raster() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_pu_polygons < get_sim_zones_pu_polygons() sim_zones_features < get_sim_zones_features() # create minimal problem p1 < problem(sim_pu_polygons, sim_features, "cost") %>% add_min_set_objective() %>% add_relative_targets(0.2) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # create problem with added locked out constraints using integers p2 < p1 %>% add_locked_out_constraints(which(sim_pu_polygons$locked_out)) # create problem with added locked out constraints using a column name p3 < p1 %>% add_locked_out_constraints("locked_out") # create problem with added locked out constraints using raster data p4 < p1 %>% add_locked_out_constraints(sim_locked_out_raster) # create problem with added locked out constraints using spatial polygon data locked_out < sim_pu_polygons[sim_pu_polygons$locked_out == 1, ] p5 < p1 %>% add_locked_out_constraints(locked_out) # solve problems s1 < solve(p1) s2 < solve(p2) s3 < solve(p3) s4 < solve(p4) s5 < solve(p5) # create single object with all solutions s6 < sf::st_sf( tibble::tibble( s1 = s1$solution_1, s2 = s2$solution_1, s3 = s3$solution_1, s4 = s4$solution_1, s5 = s5$solution_1 ), geometry = sf::st_geometry(s1) ) # plot solutions plot( s6, main = c( "none locked out", "locked out (integer input)", "locked out (character input)", "locked out (raster input)", "locked out (polygon input)" ) ) # reset plot par(mfrow = c(1, 1)) # create minimal multizone problem with spatial data p7 < problem( sim_zones_pu_polygons, sim_zones_features, cost_column = c("cost_1", "cost_2", "cost_3") ) %>% add_min_set_objective() %>% add_absolute_targets(matrix(rpois(15, 1), nrow = 5, ncol = 3)) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # create multizone problem with locked out constraints using matrix data locked_matrix < as.matrix(sf::st_drop_geometry( sim_zones_pu_polygons[, c("locked_1", "locked_2", "locked_3")] )) p8 < p7 %>% add_locked_out_constraints(locked_matrix) # solve problem s8 < solve(p8) # create new column representing the zone id that each planning unit # was allocated to in the solution s8$solution < category_vector(sf::st_drop_geometry( s8[, c("solution_1_zone_1", "solution_1_zone_2", "solution_1_zone_3")] )) s8$solution < factor(s8$solution) # plot solution plot(s8[, "solution"], main = "solution", axes = FALSE) # create multizone problem with locked out constraints using column names p9 < p7 %>% add_locked_out_constraints(c("locked_1", "locked_2", "locked_3")) # solve problem s9 < solve(p9) # create new column in s8 representing the zone id that each planning unit # was allocated to in the solution s9$solution < category_vector(sf::st_drop_geometry( s9[, c("solution_1_zone_1", "solution_1_zone_2", "solution_1_zone_3")] )) s9$solution[s9$solution == 1 & s9$solution_1_zone_1 == 0] < 0 s9$solution < factor(s9$solution) # plot solution plot(s9[, "solution"], main = "solution", axes = FALSE) # create multizone problem with raster planning units p10 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_min_set_objective() %>% add_absolute_targets(matrix(rpois(15, 1), nrow = 5, ncol = 3)) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # create multilayer raster with locked out units locked_out_raster < sim_zones_pu_raster[[1]] locked_out_raster[!is.na(locked_out_raster)] < 0 locked_out_raster < locked_out_raster[[c(1, 1, 1)]] names(locked_out_raster) < c("zones_1", "zones_2", "zones_3") locked_out_raster[[1]][1] < 1 locked_out_raster[[2]][2] < 1 locked_out_raster[[3]][3] < 1 # plot locked out raster plot(locked_out_raster) # add locked out raster units to problem p10 < p10 %>% add_locked_out_constraints(locked_out_raster) # solve problem s10 < solve(p10) # plot solution plot(category_layer(s10), main = "solution", axes = FALSE) ## End(Not run)
Add targets to a conservation planning problem by loglinearly interpolating the targets between thresholds based on the total amount of each feature in the study area (Rodrigues et al. 2004). Additionally, caps can be applied to targets to prevent features with massive distributions from being overrepresented in solutions (Butchart et al. 2015).
add_loglinear_targets( x, lower_bound_amount, lower_bound_target, upper_bound_amount, upper_bound_target, cap_amount = NULL, cap_target = NULL, abundances = feature_abundances(x, na.rm = FALSE)$absolute_abundance )
add_loglinear_targets( x, lower_bound_amount, lower_bound_target, upper_bound_amount, upper_bound_target, cap_amount = NULL, cap_target = NULL, abundances = feature_abundances(x, na.rm = FALSE)$absolute_abundance )
x 

lower_bound_amount 

lower_bound_target 

upper_bound_amount 

upper_bound_target 

cap_amount 

cap_target 

abundances 

Targets are used to specify the minimum amount or proportion of a
feature's distribution that needs to be protected. All conservation
planning problems require adding targets with the exception of the maximum
cover problem (see add_max_cover_objective()
), which maximizes
all features in the solution and therefore does not require targets.
Seven parameters are used to calculate the targets:
lower_bound_amount
specifies the first range size threshold,
lower_bound_target
specifies the relative target required for
species with a range size equal to or less than the first threshold,
upper_bound_amount
specifies the second range size threshold,
upper_bound_target
specifies the relative target required for
species with a range size equal to or greater than the second threshold,
cap_amount
specifies the third range size threshold,
cap_target
specifies the absolute target that is uniformly applied
to species with a range size larger than that third threshold, and finally
abundances
specifies the range size for each feature
that should be used when calculating the targets.
The target calculations do not account for the
size of each planning unit. Therefore, the feature data should account for
the size of each planning unit if this is important (e.g., pixel values in
the argument to features
in the function problem()
could
correspond to amount of land occupied by the feature in $km^2$
units).
Additionally, the function can only be applied to
problem()
objects that are associated with a
single zone.
An updated problem()
object with the targets added to it.
Early versions (< 5.0.2.4) used different equations for calculating targets.
Rodrigues ASL, Akcakaya HR, Andelman SJ, Bakarr MI, Boitani L, Brooks TM, Chanson JS, Fishpool LDC, da Fonseca GAB, Gaston KJ, and others (2004) Global gap analysis: priority regions for expanding the global protectedarea network. BioScience, 54: 1092–1100.
Butchart SHM, Clarke M, Smith RJ, Sykes RE, Scharlemann JPW, Harfoot M, Buchanan, GM, Angulo A, Balmford A, Bertzky B, and others (2015) Shortfalls and solutions for meeting national and global conservation area targets. Conservation Letters, 8: 329–337.
See targets for an overview of all functions for adding targets.
Other targets:
add_absolute_targets()
,
add_manual_targets()
,
add_relative_targets()
## Not run: # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() # create problem using loglinear targets p < problem(sim_pu_raster, sim_features) %>% add_min_set_objective() %>% add_loglinear_targets(10, 0.9, 100, 0.2) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s < solve(p) # plot solution plot(s, main = "solution", axes = FALSE) ## End(Not run)
## Not run: # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() # create problem using loglinear targets p < problem(sim_pu_raster, sim_features) %>% add_min_set_objective() %>% add_loglinear_targets(10, 0.9, 100, 0.2) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s < solve(p) # plot solution plot(s, main = "solution", axes = FALSE) ## End(Not run)
Specify that the SYMPHONY software – using the lpsymphony package – should be used to solve a conservation planning problem (Ralphs & Güzelsoy 2005). This function can also be used to customize the behavior of the solver. It requires the lpsymphony package to be installed (see below for installation instructions).
add_lpsymphony_solver( x, gap = 0.1, time_limit = .Machine$integer.max, first_feasible = FALSE, verbose = TRUE )
add_lpsymphony_solver( x, gap = 0.1, time_limit = .Machine$integer.max, first_feasible = FALSE, verbose = TRUE )
x 

gap 

time_limit 

first_feasible 

verbose 

SYMPHONY is an
opensource mixed integer programming solver that is part of the
Computational Infrastructure for Operations Research (COINOR) project.
This solver is provided because it may be easier to install
on some systems than the Rsymphony package. Additionally –
although the lpsymphony package doesn't provide the functionality
to specify the number of threads for solving a problem – the
lpsymphony package will solve problems using parallel processing
(unlike the Rsymphony package). As a consequence, this
solver will likely generate solutions much faster than the
add_rsymphony_solver()
.
Although formal benchmarks examining the performance of this solver
have yet to be completed,
please see Schuster et al. (2020) for benchmarks comparing the
run time and solution quality of the Rsymphony solver.
An updated problem()
object with the solver added to it.
The lpsymphony package is distributed through Bioconductor. To install the lpsymphony package, please use the following code:
if (!require(remotes)) install.packages("remotes") remotes::install_bioc("lpsymphony")
Ralphs TK and Güzelsoy M (2005) The SYMPHONY callable library for mixed integer programming. In The Next Wave in Computing, Optimization, and Decision Technologies (pp. 61–76). Springer, Boston, MA.
Schuster R, Hanson JO, StrimasMackey M, and Bennett JR (2020). Exact integer linear programming solvers outperform simulated annealing for solving conservation planning problems. PeerJ, 8: e9258.
Other solvers:
add_cbc_solver()
,
add_cplex_solver()
,
add_default_solver()
,
add_gurobi_solver()
,
add_highs_solver()
,
add_rsymphony_solver()
## Not run: # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() # create problem p < problem(sim_pu_raster, sim_features) %>% add_min_set_objective() %>% add_relative_targets(0.05) %>% add_proportion_decisions() %>% add_lpsymphony_solver(time_limit = 5, verbose = FALSE) # generate solution s < solve(p) # plot solution plot(s, main = "solution", axes = FALSE) ## End(Not run)
## Not run: # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() # create problem p < problem(sim_pu_raster, sim_features) %>% add_min_set_objective() %>% add_relative_targets(0.05) %>% add_proportion_decisions() %>% add_lpsymphony_solver(time_limit = 5, verbose = FALSE) # generate solution s < solve(p) # plot solution plot(s, main = "solution", axes = FALSE) ## End(Not run)
Add constraints to a conservation planning problem to ensure that every planning unit is allocated to a management zone in the solution. Note that this function can only be used with problems that contain multiple zones.
add_mandatory_allocation_constraints(x)
add_mandatory_allocation_constraints(x)
x 

For a conservation planning problem()
with multiple
management zones, it may sometimes be desirable to obtain a solution that
assigns each and every planning unit to a zone. For example, when
developing landuse plans, some decision makers may require that
every parcel of land is allocated a specific landuse type.
In other words are no "left over" areas. Although it might seem tempting
to simply solve the problem and manually assign "left over" planning units
to a default zone afterwards (e.g., an "other", "urban", or "grazing"
landuse), this could result in highly suboptimal solutions if there are
penalties for siting the default landuse adjacent to other zones.
Instead, this function can be used to specify that all planning units in a
problem with multiple zones must be allocated to a management zone (i.e.,
zone allocation is mandatory).
An updated problem()
object with the constraints added to it.
See constraints for an overview of all functions for adding constraints.
Other constraints:
add_contiguity_constraints()
,
add_feature_contiguity_constraints()
,
add_linear_constraints()
,
add_locked_in_constraints()
,
add_locked_out_constraints()
,
add_manual_bounded_constraints()
,
add_manual_locked_constraints()
,
add_neighbor_constraints()
## Not run: # set seed for reproducibility set.seed(500) # load data sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_features < get_sim_zones_features() # create multizone problem with minimum set objective targets_matrix < matrix(rpois(15, 1), nrow = 5, ncol = 3) # create minimal problem with minimum set objective p1 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_min_set_objective() %>% add_absolute_targets(targets_matrix) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # create another problem that is the same as p1, but has constraints # to mandate that every planning unit in the solution is assigned to # zone p2 < p1 %>% add_mandatory_allocation_constraints() # solve problems s1 < solve(p1) s2 < solve(p2) # convert solutions into category layers, where each pixel is assigned # value indicating which zone it was assigned to in the zone c1 < category_layer(s1) c2 < category_layer(s2) # plot solution category layers plot(c(c1, c2), main = c("default", "mandatory allocation"), axes = FALSE) ## End(Not run)
## Not run: # set seed for reproducibility set.seed(500) # load data sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_features < get_sim_zones_features() # create multizone problem with minimum set objective targets_matrix < matrix(rpois(15, 1), nrow = 5, ncol = 3) # create minimal problem with minimum set objective p1 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_min_set_objective() %>% add_absolute_targets(targets_matrix) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # create another problem that is the same as p1, but has constraints # to mandate that every planning unit in the solution is assigned to # zone p2 < p1 %>% add_mandatory_allocation_constraints() # solve problems s1 < solve(p1) s2 < solve(p2) # convert solutions into category layers, where each pixel is assigned # value indicating which zone it was assigned to in the zone c1 < category_layer(s1) c2 < category_layer(s2) # plot solution category layers plot(c(c1, c2), main = c("default", "mandatory allocation"), axes = FALSE) ## End(Not run)
Add constraints to a conservation planning problem to ensure
that the planning unit values (e.g., proportion, binary) in a solution
range between specific lower and upper bounds. This function offers more
finegrained control than the add_manual_locked_constraints()
function and is is most useful for problems involving proportiontype
or semicontinuous decisions.
add_manual_bounded_constraints(x, data) ## S4 method for signature 'ConservationProblem,data.frame' add_manual_bounded_constraints(x, data) ## S4 method for signature 'ConservationProblem,tbl_df' add_manual_bounded_constraints(x, data)
add_manual_bounded_constraints(x, data) ## S4 method for signature 'ConservationProblem,data.frame' add_manual_bounded_constraints(x, data) ## S4 method for signature 'ConservationProblem,tbl_df' add_manual_bounded_constraints(x, data)
x 

data 

An updated problem()
object with the constraints added to it.
The argument to data
should be a data.frame
with the following columns:
integer
planning unit identifier.
character
names of zones. Note that this
argument is optional for arguments to x
that contain a single
zone.
numeric
values indicating the minimum
value that each planning unit can be allocated to in each zone
in the solution.
numeric
values indicating the maximum
value that each planning unit can be allocated to in each zone
in the solution.
See constraints for an overview of all functions for adding constraints.
Other constraints:
add_contiguity_constraints()
,
add_feature_contiguity_constraints()
,
add_linear_constraints()
,
add_locked_in_constraints()
,
add_locked_out_constraints()
,
add_mandatory_allocation_constraints()
,
add_manual_locked_constraints()
,
add_neighbor_constraints()
## Not run: # set seed for reproducibility set.seed(500) # load data sim_pu_polygons < get_sim_pu_polygons() sim_features < get_sim_features() sim_zones_pu_polygons < get_sim_zones_pu_polygons() sim_zones_features < get_sim_zones_features() # create minimal problem p1 < problem(sim_pu_polygons, sim_features, "cost") %>% add_min_set_objective() %>% add_relative_targets(0.2) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # create problem with locked in constraints using add_locked_constraints p2 < p1 %>% add_locked_in_constraints("locked_in") # create identical problem using add_manual_bounded_constraints bounds_data < data.frame( pu = which(sim_pu_polygons$locked_in), lower = 1, upper = 1 ) p3 < p1 %>% add_manual_bounded_constraints(bounds_data) # solve problems s1 < solve(p1) s2 < solve(p2) s3 < solve(p3) # create object with all solutions s4 < sf::st_sf( tibble::tibble( s1 = s1$solution_1, s2 = s2$solution_1, s3 = s3$solution_1 ), geometry = sf::st_geometry(s1) ) # plot solutions ## s1 = none locked in ## s2 = locked in constraints ## s3 = manual bounds constraints plot(s4) # create minimal problem with multiple zones p5 < problem( sim_zones_pu_polygons, sim_zones_features, c("cost_1", "cost_2", "cost_3") ) %>% add_min_set_objective() %>% add_relative_targets(matrix(runif(15, 0.1, 0.2), nrow = 5, ncol = 3)) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # create data.frame with the following constraints: # planning units 1, 2, and 3 must be allocated to zone 1 in the solution # planning units 4, and 5 must be allocated to zone 2 in the solution # planning units 8 and 9 must not be allocated to zone 3 in the solution bounds_data2 < data.frame( pu = c(1, 2, 3, 4, 5, 8, 9), zone = c(rep("zone_1", 3), rep("zone_2", 2), rep("zone_3", 2)), lower = c(rep(1, 5), rep(0, 2)), upper = c(rep(1, 5), rep(0, 2)) ) # print bounds data print(bounds_data2) # create problem with added constraints p6 < p5 %>% add_manual_bounded_constraints(bounds_data2) # solve problem s5 < solve(p5) s6 < solve(p6) # create two new columns representing the zone id that each planning unit # was allocated to in the two solutions s5$solution < category_vector(sf::st_drop_geometry( s5[, c("solution_1_zone_1","solution_1_zone_2", "solution_1_zone_3")] )) s5$solution < factor(s5$solution) s5$solution_bounded < category_vector(sf::st_drop_geometry( s6[, c("solution_1_zone_1", "solution_1_zone_2", "solution_1_zone_3")] )) s5$solution_bounded < factor(s5$solution_bounded) # plot solutions plot(s5[, c("solution", "solution_bounded")], axes = FALSE) ## End(Not run)
## Not run: # set seed for reproducibility set.seed(500) # load data sim_pu_polygons < get_sim_pu_polygons() sim_features < get_sim_features() sim_zones_pu_polygons < get_sim_zones_pu_polygons() sim_zones_features < get_sim_zones_features() # create minimal problem p1 < problem(sim_pu_polygons, sim_features, "cost") %>% add_min_set_objective() %>% add_relative_targets(0.2) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # create problem with locked in constraints using add_locked_constraints p2 < p1 %>% add_locked_in_constraints("locked_in") # create identical problem using add_manual_bounded_constraints bounds_data < data.frame( pu = which(sim_pu_polygons$locked_in), lower = 1, upper = 1 ) p3 < p1 %>% add_manual_bounded_constraints(bounds_data) # solve problems s1 < solve(p1) s2 < solve(p2) s3 < solve(p3) # create object with all solutions s4 < sf::st_sf( tibble::tibble( s1 = s1$solution_1, s2 = s2$solution_1, s3 = s3$solution_1 ), geometry = sf::st_geometry(s1) ) # plot solutions ## s1 = none locked in ## s2 = locked in constraints ## s3 = manual bounds constraints plot(s4) # create minimal problem with multiple zones p5 < problem( sim_zones_pu_polygons, sim_zones_features, c("cost_1", "cost_2", "cost_3") ) %>% add_min_set_objective() %>% add_relative_targets(matrix(runif(15, 0.1, 0.2), nrow = 5, ncol = 3)) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # create data.frame with the following constraints: # planning units 1, 2, and 3 must be allocated to zone 1 in the solution # planning units 4, and 5 must be allocated to zone 2 in the solution # planning units 8 and 9 must not be allocated to zone 3 in the solution bounds_data2 < data.frame( pu = c(1, 2, 3, 4, 5, 8, 9), zone = c(rep("zone_1", 3), rep("zone_2", 2), rep("zone_3", 2)), lower = c(rep(1, 5), rep(0, 2)), upper = c(rep(1, 5), rep(0, 2)) ) # print bounds data print(bounds_data2) # create problem with added constraints p6 < p5 %>% add_manual_bounded_constraints(bounds_data2) # solve problem s5 < solve(p5) s6 < solve(p6) # create two new columns representing the zone id that each planning unit # was allocated to in the two solutions s5$solution < category_vector(sf::st_drop_geometry( s5[, c("solution_1_zone_1","solution_1_zone_2", "solution_1_zone_3")] )) s5$solution < factor(s5$solution) s5$solution_bounded < category_vector(sf::st_drop_geometry( s6[, c("solution_1_zone_1", "solution_1_zone_2", "solution_1_zone_3")] )) s5$solution_bounded < factor(s5$solution_bounded) # plot solutions plot(s5[, c("solution", "solution_bounded")], axes = FALSE) ## End(Not run)
Add constraints to a conservation planning problem to ensure
that solutions allocate (or do not allocate) specific planning units to
specific management zones. This function offers more finegrained control
than the add_locked_in_constraints()
and
add_locked_out_constraints()
functions.
add_manual_locked_constraints(x, data) ## S4 method for signature 'ConservationProblem,data.frame' add_manual_locked_constraints(x, data) ## S4 method for signature 'ConservationProblem,tbl_df' add_manual_locked_constraints(x, data)
add_manual_locked_constraints(x, data) ## S4 method for signature 'ConservationProblem,data.frame' add_manual_locked_constraints(x, data) ## S4 method for signature 'ConservationProblem,tbl_df' add_manual_locked_constraints(x, data)
x 

data 

An updated problem()
object with the constraints added to it.
The argument to data
should be a data.frame
with the following columns:
integer
planning unit identifier.
character
names of zones. Note that this
argument is optional for arguments to x
that contain a single
zone.
numeric
values indicating how much
of each planning unit should be allocated to each zone in the solution.
For example, the numeric
values could be binary values (i.e., zero
or one) for problems containing binarytype decision variables
(using the add_binary_decisions()
function). Alternatively,
the numeric
values could be proportions (e.g., 0.5) for problems
containing proportiontype decision variables (using the
add_proportion_decisions()
).
See constraints for an overview of all functions for adding constraints.
Other constraints:
add_contiguity_constraints()
,
add_feature_contiguity_constraints()
,
add_linear_constraints()
,
add_locked_in_constraints()
,
add_locked_out_constraints()
,
add_mandatory_allocation_constraints()
,
add_manual_bounded_constraints()
,
add_neighbor_constraints()
## Not run: # set seed for reproducibility set.seed(500) # load data sim_pu_polygons < get_sim_pu_polygons() sim_features < get_sim_features() sim_zones_pu_polygons < get_sim_zones_pu_polygons() sim_zones_features < get_sim_zones_features() # create minimal problem p1 < problem(sim_pu_polygons, sim_features, "cost") %>% add_min_set_objective() %>% add_relative_targets(0.2) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # create problem with locked in constraints using add_locked_constraints p2 < p1 %>% add_locked_in_constraints("locked_in") # create identical problem using add_manual_locked_constraints locked_data < data.frame( pu = which(sim_pu_polygons$locked_in), status = 1 ) p3 < p1 %>% add_manual_locked_constraints(locked_data) # solve problems s1 < solve(p1) s2 < solve(p2) s3 < solve(p3) # create object with all solutions s4 < sf::st_sf( tibble::tibble( s1 = s1$solution_1, s2 = s2$solution_1, s3 = s3$solution_1 ), geometry = sf::st_geometry(s1) ) # plot solutions ## s1 = none locked in ## s2 = locked in constraints ## s3 = manual locked constraints plot(s4) # create minimal problem with multiple zones p5 < problem( sim_zones_pu_polygons, sim_zones_features, c("cost_1", "cost_2", "cost_3") ) %>% add_min_set_objective() %>% add_relative_targets(matrix(runif(15, 0.1, 0.2), nrow = 5, ncol = 3)) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # create data.frame with the following constraints: # planning units 1, 2, and 3 must be allocated to zone 1 in the solution # planning units 4, and 5 must be allocated to zone 2 in the solution # planning units 8 and 9 must not be allocated to zone 3 in the solution locked_data2 < data.frame( pu = c(1, 2, 3, 4, 5, 8, 9), zone = c(rep("zone_1", 3), rep("zone_2", 2),rep("zone_3", 2)), status = c(rep(1, 5), rep(0, 2)) ) # print locked constraint data print(locked_data2) # create problem with added constraints p6 < p5 %>% add_manual_locked_constraints(locked_data2) # solve problem s5 < solve(p5) s6 < solve(p6) # create two new columns representing the zone id that each planning unit # was allocated to in the two solutions s5$solution < category_vector(sf::st_drop_geometry( s5[, c("solution_1_zone_1", "solution_1_zone_2", "solution_1_zone_3")] )) s5$solution < factor(s5$solution) s5$solution_locked < category_vector(sf::st_drop_geometry( s6[, c("solution_1_zone_1", "solution_1_zone_2", "solution_1_zone_3")] )) s5$solution_locked < factor(s5$solution_locked) # plot solutions plot(s5[, c("solution", "solution_locked")], axes = FALSE) ## End(Not run)
## Not run: # set seed for reproducibility set.seed(500) # load data sim_pu_polygons < get_sim_pu_polygons() sim_features < get_sim_features() sim_zones_pu_polygons < get_sim_zones_pu_polygons() sim_zones_features < get_sim_zones_features() # create minimal problem p1 < problem(sim_pu_polygons, sim_features, "cost") %>% add_min_set_objective() %>% add_relative_targets(0.2) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # create problem with locked in constraints using add_locked_constraints p2 < p1 %>% add_locked_in_constraints("locked_in") # create identical problem using add_manual_locked_constraints locked_data < data.frame( pu = which(sim_pu_polygons$locked_in), status = 1 ) p3 < p1 %>% add_manual_locked_constraints(locked_data) # solve problems s1 < solve(p1) s2 < solve(p2) s3 < solve(p3) # create object with all solutions s4 < sf::st_sf( tibble::tibble( s1 = s1$solution_1, s2 = s2$solution_1, s3 = s3$solution_1 ), geometry = sf::st_geometry(s1) ) # plot solutions ## s1 = none locked in ## s2 = locked in constraints ## s3 = manual locked constraints plot(s4) # create minimal problem with multiple zones p5 < problem( sim_zones_pu_polygons, sim_zones_features, c("cost_1", "cost_2", "cost_3") ) %>% add_min_set_objective() %>% add_relative_targets(matrix(runif(15, 0.1, 0.2), nrow = 5, ncol = 3)) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # create data.frame with the following constraints: # planning units 1, 2, and 3 must be allocated to zone 1 in the solution # planning units 4, and 5 must be allocated to zone 2 in the solution # planning units 8 and 9 must not be allocated to zone 3 in the solution locked_data2 < data.frame( pu = c(1, 2, 3, 4, 5, 8, 9), zone = c(rep("zone_1", 3), rep("zone_2", 2),rep("zone_3", 2)), status = c(rep(1, 5), rep(0, 2)) ) # print locked constraint data print(locked_data2) # create problem with added constraints p6 < p5 %>% add_manual_locked_constraints(locked_data2) # solve problem s5 < solve(p5) s6 < solve(p6) # create two new columns representing the zone id that each planning unit # was allocated to in the two solutions s5$solution < category_vector(sf::st_drop_geometry( s5[, c("solution_1_zone_1", "solution_1_zone_2", "solution_1_zone_3")] )) s5$solution < factor(s5$solution) s5$solution_locked < category_vector(sf::st_drop_geometry( s6[, c("solution_1_zone_1", "solution_1_zone_2", "solution_1_zone_3")] )) s5$solution_locked < factor(s5$solution_locked) # plot solutions plot(s5[, c("solution", "solution_locked")], axes = FALSE) ## End(Not run)
Set targets for a conservation planning problem by manually
specifying all the required information for each target. This function
is useful because it can be used to customize all aspects of a target. For
most cases, targets can be specified using the
add_absolute_targets()
and add_relative_targets()
functions. However, this function can be used to (i) mix absolute and
relative targets for different features and zones, (ii) set targets that
pertain to the allocations of planning units in multiple zones, and (iii)
set targets that require different senses (e.g., targets which specify the
solution should not exceed a certain quantity using "<="
values).
add_manual_targets(x, targets) ## S4 method for signature 'ConservationProblem,data.frame' add_manual_targets(x, targets) ## S4 method for signature 'ConservationProblem,tbl_df' add_manual_targets(x, targets)
add_manual_targets(x, targets) ## S4 method for signature 'ConservationProblem,data.frame' add_manual_targets(x, targets) ## S4 method for signature 'ConservationProblem,tbl_df' add_manual_targets(x, targets)
x 

targets 

Targets are used to specify the minimum amount or proportion of a
feature's distribution that needs to be protected. Most conservation
planning problems require targets with the exception of the maximum cover
(see add_max_cover_objective()
) and maximum utility
(see add_max_utility_objective()
) problems. Attempting to solve
problems with objectives that require targets without specifying targets
will throw an error.
For problems associated with multiple management zones,
add_absolute_targets()
can
be used to set targets that each pertain to a single feature and a single
zone. To set targets that can be met through allocating different
planning units to multiple zones, see the add_manual_targets()
function. An example of a target that could be met through allocations
to multiple zones might be where each management zone is expected to
result in a different amount of a feature and the target requires that
the total amount of the feature in all zones must exceed a certain
threshold. In other words, the target does not require that any single
zone secure a specific amount of the feature, but the total amount held
in all zones must secure a specific amount. Thus the target could,
potentially, be met through allocating all planning units to any specific
management zone, or through allocating the planning units to different
combinations of management zones.
An updated problem()
object with the targets added to it.
The targets
argument should be a data.frame
with the following
columns:
character
name of features in argument
to x
.
character
name of zones in the argument
x
. It can also be a list
of character
vectors if
targets should correspond to multiple zones (see Examples section below).
This column is optional for arguments to x
that do not contain multiple zones.
character
describing the type of target.
Acceptable values include "absolute"
and "relative"
.
These values correspond to add_absolute_targets()
,
and add_relative_targets()
respectively.
character
sense of the target. Acceptable
values include: ">="
, "<="
, and "="
. This
column is optional and if it is missing then target senses will
default to ">="
values.
numeric
target threshold.
See targets for an overview of all functions for adding targets.
Other targets:
add_absolute_targets()
,
add_loglinear_targets()
,
add_relative_targets()
## Not run: # set seed for reproducibility set.seed(500) # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_features < get_sim_zones_features() # create problem with 10% relative targets p1 < problem(sim_pu_raster, sim_features) %>% add_min_set_objective() %>% add_relative_targets(0.1) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s1 < solve(p1) # plot solution plot(s1, main = "solution", axes = FALSE) # create equivalent problem using add_manual_targets p2 < problem(sim_pu_raster, sim_features) %>% add_min_set_objective() %>% add_manual_targets( data.frame( feature = names(sim_features), type = "relative", sense = ">=", target = 0.1 ) ) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s2 < solve(p2) # plot solution plot(s2, main = "solution", axes = FALSE) # create problem with targets set for only a few features p3 < problem(sim_pu_raster, sim_features) %>% add_min_set_objective() %>% add_manual_targets( data.frame( feature = names(sim_features)[1:3], type = "relative", sense = ">=", target = 0.1 ) ) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s3 < solve(p3) # plot solution plot(s3, main = "solution", axes = FALSE) # create problem that aims to secure at least 10% of the habitat for one # feature whilst ensuring that the solution does not capture more than # 20 units habitat for different feature # create problem with targets set for only a few features p4 < problem(sim_pu_raster, sim_features[[1:2]]) %>% add_min_set_objective() %>% add_manual_targets( data.frame( feature = names(sim_features)[1:2], type = "relative", sense = c(">=", "<="), target = c(0.1, 0.2) ) ) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s4 < solve(p4) # plot solution plot(s4, main = "solution", axes = FALSE) # create a multizone problem that requires a specific amount of each # feature in each zone targets_matrix < matrix(rpois(15, 1), nrow = 5, ncol = 3) p5 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_min_set_objective() %>% add_absolute_targets(targets_matrix) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s5 < solve(p5) # plot solution plot(category_layer(s5), main = "solution", axes = FALSE) # create equivalent problem using add_manual_targets targets_dataframe < expand.grid( feature = feature_names(sim_zones_features), zone = zone_names(sim_zones_features), sense = ">=", type = "absolute" ) targets_dataframe$target < c(targets_matrix) p6 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_min_set_objective() %>% add_manual_targets(targets_dataframe) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s6 < solve(p6) # plot solution plot(category_layer(s6), main = "solution", axes = FALSE) # create a problem that requires a total of 20 units of habitat to be # captured for two species. This can be achieved through representing # habitat in two zones. The first zone represents a full restoration of the # habitat and a second zone represents a partial restoration of the habitat # Thus only half of the benefit that would have been gained from the full # restoration is obtained when planning units are allocated a partial # restoration # create data spp_zone1 < as.list(sim_zones_features)[[1]][[1:2]] spp_zone2 < spp_zone1 * 0.5 costs < sim_zones_pu_raster[[1:2]] # create targets targets_dataframe2 < tibble::tibble( feature = names(spp_zone1), zone = list(c("z1", "z2"), c("z1", "z2")), sense = c(">=", ">="), type = c("absolute", "absolute"), target = c(20, 20) ) # create problem p7 < problem( costs, zones( spp_zone1, spp_zone2, feature_names = names(spp_zone1), zone_names = c("z1", "z2") ) ) %>% add_min_set_objective() %>% add_manual_targets(targets_dataframe2) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s7 < solve(p7) # plot solution plot(category_layer(s7), main = "solution", axes = FALSE) ## End(Not run)
## Not run: # set seed for reproducibility set.seed(500) # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_features < get_sim_zones_features() # create problem with 10% relative targets p1 < problem(sim_pu_raster, sim_features) %>% add_min_set_objective() %>% add_relative_targets(0.1) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s1 < solve(p1) # plot solution plot(s1, main = "solution", axes = FALSE) # create equivalent problem using add_manual_targets p2 < problem(sim_pu_raster, sim_features) %>% add_min_set_objective() %>% add_manual_targets( data.frame( feature = names(sim_features), type = "relative", sense = ">=", target = 0.1 ) ) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s2 < solve(p2) # plot solution plot(s2, main = "solution", axes = FALSE) # create problem with targets set for only a few features p3 < problem(sim_pu_raster, sim_features) %>% add_min_set_objective() %>% add_manual_targets( data.frame( feature = names(sim_features)[1:3], type = "relative", sense = ">=", target = 0.1 ) ) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s3 < solve(p3) # plot solution plot(s3, main = "solution", axes = FALSE) # create problem that aims to secure at least 10% of the habitat for one # feature whilst ensuring that the solution does not capture more than # 20 units habitat for different feature # create problem with targets set for only a few features p4 < problem(sim_pu_raster, sim_features[[1:2]]) %>% add_min_set_objective() %>% add_manual_targets( data.frame( feature = names(sim_features)[1:2], type = "relative", sense = c(">=", "<="), target = c(0.1, 0.2) ) ) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s4 < solve(p4) # plot solution plot(s4, main = "solution", axes = FALSE) # create a multizone problem that requires a specific amount of each # feature in each zone targets_matrix < matrix(rpois(15, 1), nrow = 5, ncol = 3) p5 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_min_set_objective() %>% add_absolute_targets(targets_matrix) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s5 < solve(p5) # plot solution plot(category_layer(s5), main = "solution", axes = FALSE) # create equivalent problem using add_manual_targets targets_dataframe < expand.grid( feature = feature_names(sim_zones_features), zone = zone_names(sim_zones_features), sense = ">=", type = "absolute" ) targets_dataframe$target < c(targets_matrix) p6 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_min_set_objective() %>% add_manual_targets(targets_dataframe) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s6 < solve(p6) # plot solution plot(category_layer(s6), main = "solution", axes = FALSE) # create a problem that requires a total of 20 units of habitat to be # captured for two species. This can be achieved through representing # habitat in two zones. The first zone represents a full restoration of the # habitat and a second zone represents a partial restoration of the habitat # Thus only half of the benefit that would have been gained from the full # restoration is obtained when planning units are allocated a partial # restoration # create data spp_zone1 < as.list(sim_zones_features)[[1]][[1:2]] spp_zone2 < spp_zone1 * 0.5 costs < sim_zones_pu_raster[[1:2]] # create targets targets_dataframe2 < tibble::tibble( feature = names(spp_zone1), zone = list(c("z1", "z2"), c("z1", "z2")), sense = c(">=", ">="), type = c("absolute", "absolute"), target = c(20, 20) ) # create problem p7 < problem( costs, zones( spp_zone1, spp_zone2, feature_names = names(spp_zone1), zone_names = c("z1", "z2") ) ) %>% add_min_set_objective() %>% add_manual_targets(targets_dataframe2) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s7 < solve(p7) # plot solution plot(category_layer(s7), main = "solution", axes = FALSE) ## End(Not run)
Set the objective of a conservation planning problem to represent at least one instance of as many features as possible within a given budget. This objective does not use targets, and feature weights should be used instead to increase the representation of certain features by a solution.
add_max_cover_objective(x, budget)
add_max_cover_objective(x, budget)
x 

budget 

The maximum coverage objective seeks to find the set of planning units that
maximizes the number of represented features, while keeping cost within a
fixed budget. Here, features are treated as being represented if
the reserve system contains at least a single instance of a feature
(i.e., an amount greater than 1). This formulation has often been
used in conservation planning problems dealing with binary biodiversity
data that indicate the presence/absence of suitable habitat
(e.g., Church & Velle 1974). Additionally, weights can be used to favor the
representation of certain features over other features (see
add_feature_weights()
). Check out the
add_max_features_objective()
for a more
generalized formulation which can accommodate userspecified representation
targets.
An updated problem()
object with the objective added to it.
This objective is based on the maximum coverage reserve
selection problem (Church & Velle 1974; Church et al. 1996).
The maximum coverage objective for the reserve design problem can be
expressed mathematically for a set of planning units ($I$
indexed by
$i$
) and a set of features ($J$
indexed by $j$
) as:
$\mathit{Maximize} \space \sum_{i = 1}^{I} s \space c_i \space x_i +
\sum_{j = 1}^{J} y_j w_j \\
\mathit{subject \space to} \\
\sum_{i = 1}^{I} x_i r_{ij} \geq y_j \times 1 \forall j \in J \\
\sum_{i = 1}^{I} x_i c_i \leq B$
Here, $x_i$
is the decisions variable (e.g.,
specifying whether planning unit $i$
has been selected (1) or not
(0)), $r_{ij}$
is the amount of feature $j$
in planning
unit $i$
, $y_j$
indicates if the solution has meet
the target $t_j$
for feature $j$
, and $w_j$
is the
weight for feature $j$
(defaults to 1 for all features; see
add_feature_weights()
to specify weights). Additionally,
$B$
is the budget allocated for the solution, $c_i$
is the
cost of planning unit $i$
, and $s$
is a scaling factor used
to shrink the costs so that the problem will return a cheapest solution
when there are multiple solutions that represent the same amount of all
features within the budget.
In early versions (< 3.0.0.0), the mathematical formulation
underpinning this function was very different. Specifically,
as described above, the function now follows the formulations outlined in
Church et al. (1996). The old formulation is now provided by the
add_max_utility_objective()
function.
Church RL and Velle CR (1974) The maximum covering location problem. Regional Science, 32: 101–118.
Church RL, Stoms DM, and Davis FW (1996) Reserve selection as a maximum covering location problem. Biological Conservation, 76: 105–112.
See objectives for an overview of all functions for adding objectives.
Also, see add_feature_weights()
to specify weights for different features.
Other objectives:
add_max_features_objective()
,
add_max_phylo_div_objective()
,
add_max_phylo_end_objective()
,
add_max_utility_objective()
,
add_min_largest_shortfall_objective()
,
add_min_set_objective()
,
add_min_shortfall_objective()
## Not run: # load data sim_pu_raster < get_sim_pu_raster() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_features < get_sim_features() sim_zones_features < get_sim_zones_features() # threshold the feature data to generate binary biodiversity data sim_binary_features < sim_features thresholds < terra::global( sim_features, fun = quantile, probs = 0.5, na.rm = TRUE ) for (i in seq_len(terra::nlyr(sim_features))) { sim_binary_features[[i]] < terra::as.int( sim_features[[i]] > thresholds[[1]][[i]] ) } # create problem with maximum cover objective p1 < problem(sim_pu_raster, sim_binary_features) %>% add_max_cover_objective(500) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s1 < solve(p1) # plot solution plot(s1, main = "solution", axes = FALSE) # threshold the multizone feature data to generate binary biodiversity data sim_binary_features_zones < sim_zones_features for (z in seq_len(number_of_zones(sim_zones_features))) { thresholds < terra::global( sim_zones_features[[z]], fun = quantile, probs = 0.5, na.rm = TRUE ) for (i in seq_len(number_of_features(sim_zones_features))) { sim_binary_features_zones[[z]][[i]] < terra::as.int( sim_zones_features[[z]][[i]] > thresholds[[1]][[i]] ) } } # create multizone problem with maximum cover objective that # has a single budget for all zones p2 < problem(sim_zones_pu_raster, sim_binary_features_zones) %>% add_max_cover_objective(800) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s2 < solve(p2) # plot solution plot(category_layer(s2), main = "solution", axes = FALSE) # create multizone problem with maximum cover objective that # has separate budgets for each zone p3 < problem(sim_zones_pu_raster, sim_binary_features_zones) %>% add_max_cover_objective(c(400, 400, 400)) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s3 < solve(p3) # plot solution plot(category_layer(s3), main = "solution", axes = FALSE) ## End(Not run)
## Not run: # load data sim_pu_raster < get_sim_pu_raster() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_features < get_sim_features() sim_zones_features < get_sim_zones_features() # threshold the feature data to generate binary biodiversity data sim_binary_features < sim_features thresholds < terra::global( sim_features, fun = quantile, probs = 0.5, na.rm = TRUE ) for (i in seq_len(terra::nlyr(sim_features))) { sim_binary_features[[i]] < terra::as.int( sim_features[[i]] > thresholds[[1]][[i]] ) } # create problem with maximum cover objective p1 < problem(sim_pu_raster, sim_binary_features) %>% add_max_cover_objective(500) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s1 < solve(p1) # plot solution plot(s1, main = "solution", axes = FALSE) # threshold the multizone feature data to generate binary biodiversity data sim_binary_features_zones < sim_zones_features for (z in seq_len(number_of_zones(sim_zones_features))) { thresholds < terra::global( sim_zones_features[[z]], fun = quantile, probs = 0.5, na.rm = TRUE ) for (i in seq_len(number_of_features(sim_zones_features))) { sim_binary_features_zones[[z]][[i]] < terra::as.int( sim_zones_features[[z]][[i]] > thresholds[[1]][[i]] ) } } # create multizone problem with maximum cover objective that # has a single budget for all zones p2 < problem(sim_zones_pu_raster, sim_binary_features_zones) %>% add_max_cover_objective(800) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s2 < solve(p2) # plot solution plot(category_layer(s2), main = "solution", axes = FALSE) # create multizone problem with maximum cover objective that # has separate budgets for each zone p3 < problem(sim_zones_pu_raster, sim_binary_features_zones) %>% add_max_cover_objective(c(400, 400, 400)) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s3 < solve(p3) # plot solution plot(category_layer(s3), main = "solution", axes = FALSE) ## End(Not run)
Set the objective of a conservation planning problem to fulfill as many targets as possible, whilst ensuring that the cost of the solution does not exceed a budget.
add_max_features_objective(x, budget)
add_max_features_objective(x, budget)
x 

budget 

The maximum feature representation objective is an enhanced version of the
maximum coverage objective add_max_cover_objective()
because
targets can be used to ensure that a certain amount of each feature is
required in order for them to be adequately represented (similar to the
minimum set objective (see add_min_set_objective()
). This
objective finds the set of planning units that meets representation targets
for as many features as possible while staying within a fixed budget
(inspired by Cabeza and Moilanen 2001). Additionally, weights can be used
add_feature_weights()
). If multiple solutions can meet the same
number of weighted targets while staying within budget, the cheapest
solution is returned.
An updated problem()
object with the objective added to it.
This objective can be expressed mathematically for a set of planning units
($I$
indexed by
$i$
) and a set of features ($J$
indexed by $j$
) as:
$\mathit{Maximize} \space \sum_{i = 1}^{I} s \space c_i \space x_i +
\sum_{j = 1}^{J} y_j w_j \\
\mathit{subject \space to} \\
\sum_{i = 1}^{I} x_i r_{ij} \geq y_j t_j \forall j \in J \\
\sum_{i = 1}^{I} x_i c_i \leq B$
Here, $x_i$
is the decisions variable (e.g.,
specifying whether planning unit $i$
has been selected (1) or not
(0)), $r_{ij}$
is the amount of feature $j$
in planning
unit $i$
, $t_j$
is the representation target for feature
$j$
, $y_j$
indicates if the solution has meet
the target $t_j$
for feature $j$
, and $w_j$
is the
weight for feature $j$
(defaults to 1 for all features; see
add_feature_weights()
to specify weights). Additionally,
$B$
is the budget allocated for the solution, $c_i$
is the
cost of planning unit $i$
, and $s$
is a scaling factor used
to shrink the costs so that the problem will return a cheapest solution
when there are multiple solutions that represent the same amount of all
features within the budget.
Cabeza M and Moilanen A (2001) Design of reserve networks and the persistence of biodiversity. Trends in Ecology & Evolution, 16: 242–248.
See objectives for an overview of all functions for adding objectives.
Also, see targets for an overview of all functions for adding targets, and
add_feature_weights()
to specify weights for different features.
Other objectives:
add_max_cover_objective()
,
add_max_phylo_div_objective()
,
add_max_phylo_end_objective()
,
add_max_utility_objective()
,
add_min_largest_shortfall_objective()
,
add_min_set_objective()
,
add_min_shortfall_objective()
## Not run: # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_features < get_sim_zones_features() # create problem with maximum features objective p1 < problem(sim_pu_raster, sim_features) %>% add_max_features_objective(1800) %>% add_relative_targets(0.1) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s1 < solve(p1) # plot solution plot(s1, main = "solution", axes = FALSE) # create multizone problem with maximum features objective, # with 10% representation targets for each feature, and set # a budget such that the total maximum expenditure in all zones # cannot exceed 3000 p2 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_max_features_objective(3000) %>% add_relative_targets(matrix(0.1, ncol = 3, nrow = 5)) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s2 < solve(p2) # plot solution plot(category_layer(s2), main = "solution", axes = FALSE) # create multizone problem with maximum features objective, # with 10% representation targets for each feature, and set # separate budgets for each management zone p3 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_max_features_objective(c(3000, 3000, 3000)) %>% add_relative_targets(matrix(0.1, ncol = 3, nrow = 5)) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s3 < solve(p3) # plot solution plot(category_layer(s3), main = "solution", axes = FALSE) ## End(Not run)
## Not run: # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_features < get_sim_zones_features() # create problem with maximum features objective p1 < problem(sim_pu_raster, sim_features) %>% add_max_features_objective(1800) %>% add_relative_targets(0.1) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s1 < solve(p1) # plot solution plot(s1, main = "solution", axes = FALSE) # create multizone problem with maximum features objective, # with 10% representation targets for each feature, and set # a budget such that the total maximum expenditure in all zones # cannot exceed 3000 p2 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_max_features_objective(3000) %>% add_relative_targets(matrix(0.1, ncol = 3, nrow = 5)) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s2 < solve(p2) # plot solution plot(category_layer(s2), main = "solution", axes = FALSE) # create multizone problem with maximum features objective, # with 10% representation targets for each feature, and set # separate budgets for each management zone p3 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_max_features_objective(c(3000, 3000, 3000)) %>% add_relative_targets(matrix(0.1, ncol = 3, nrow = 5)) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s3 < solve(p3) # plot solution plot(category_layer(s3), main = "solution", axes = FALSE) ## End(Not run)
Set the objective of a conservation planning problem to
maximize the phylogenetic diversity of the features represented in the
solution subject to a budget. This objective is similar to
add_max_features_objective()
except
that emphasis is placed on representing a phylogenetically diverse set of
species, rather than as many features as possible (subject to weights).
This function was inspired by Faith (1992) and Rodrigues et al.
(2002).
add_max_phylo_div_objective(x, budget, tree)
add_max_phylo_div_objective(x, budget, tree)
x 

budget 

tree 

The maximum phylogenetic diversity objective finds the set of
planning units that meets representation targets for a phylogenetic tree
while staying within a fixed budget. If multiple solutions can meet all
targets while staying within budget, the cheapest solution is chosen.
Note that this objective is similar to the maximum
features objective (add_max_features_objective()
) in that it
allows for both a budget and targets to be set for each feature. However,
unlike the maximum feature objective, the aim of this objective is to
maximize the total phylogenetic diversity of the targets met in the
solution, so if multiple targets are provided for a single feature, the
problem will only need to meet a single target for that feature
for the phylogenetic benefit for that feature to be counted when
calculating the phylogenetic diversity of the solution. In other words,
for multizone problems, this objective does not aim to maximize the
phylogenetic diversity in each zone, but rather this objective
aims to maximize the phylogenetic diversity of targets that can be met
through allocating planning units to any of the different zones in a
problem. This can be useful for problems where targets pertain to the total
amount held for each feature across multiple zones. For example,
each feature might have a nonzero amount of suitable habitat in each
planning unit when the planning units are assigned to a (i) not restored,
(ii) partially restored, or (iii) completely restored management zone.
Here each target corresponds to a single feature and can be met through
the total amount of habitat in planning units present to the three
zones.
An updated problem()
object with the objective added to it.
This objective can be expressed mathematically for a set of planning units
($I$
indexed by $i$
) and a set of features ($J$
indexed by $j$
) as:
$\mathit{Maximize} \space \sum_{i = 1}^{I} s \space c_i \space x_i +
\sum_{j = 1}^{J} m_b l_b \\
\mathit{subject \space to} \\
\sum_{i = 1}^{I} x_i r_{ij} \geq y_j t_j \forall j \in J \\
m_b \leq y_j \forall j \in T(b) \\
\sum_{i = 1}^{I} x_i c_i \leq B$
Here, $x_i$
is the decisions variable (e.g.,
specifying whether planning unit $i$
has been selected (1) or not
(0)), $r_{ij}$
is the amount of feature $j$
in planning
unit $i$
, $t_j$
is the representation target for feature
$j$
, $y_j$
indicates if the solution has meet
the target $t_j$
for feature $j$
. Additionally,
$T$
represents a phylogenetic tree containing features $j$
and has the branches $b$
associated within lengths $l_b$
.
The binary variable $m_b$
denotes if
at least one feature associated with the branch $b$
has met its
representation as indicated by $y_j$
. For brevity, we denote
the features $j$
associated with branch $b$
using
$T(b)$
. Finally, $B$
is the budget allocated for the
solution, $c_i$
is the cost of planning unit $i$
, and
$s$
is a scaling factor used to shrink the costs so that the problem
will return a cheapest solution when there are multiple solutions that
represent the same amount of all features within the budget.
In early versions, this function was named as the
add_max_phylo_div_objective
function.
Faith DP (1992) Conservation evaluation and phylogenetic diversity. Biological Conservation, 61: 1–10.
Rodrigues ASL and Gaston KJ (2002) Maximising phylogenetic diversity in the selection of networks of conservation areas. Biological Conservation, 105: 103–111.
See objectives for an overview of all functions for adding objectives.
Also, see targets for an overview of all functions for adding targets, and
add_feature_weights()
to specify weights for different features.
Other objectives:
add_max_cover_objective()
,
add_max_features_objective()
,
add_max_phylo_end_objective()
,
add_max_utility_objective()
,
add_min_largest_shortfall_objective()
,
add_min_set_objective()
,
add_min_shortfall_objective()
## Not run: # load ape package require(ape) # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() sim_phylogeny < get_sim_phylogeny() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_features < get_sim_zones_features() # plot the simulated phylogeny par(mfrow = c(1, 1)) plot(sim_phylogeny, main = "phylogeny") # create problem with a maximum phylogenetic diversity objective, # where each feature needs 10% of its distribution to be secured for # it to be adequately conserved and a total budget of 1900 p1 < problem(sim_pu_raster, sim_features) %>% add_max_phylo_div_objective(1900, sim_phylogeny) %>% add_relative_targets(0.1) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s1 < solve(p1) # plot solution plot(s1, main = "solution", axes = FALSE) # find out which features have their targets met r1 < eval_target_coverage_summary(p1, s1) print(r1, width = Inf) # plot the phylogeny and color the adequately represented features in red plot( sim_phylogeny, main = "adequately represented features", tip.color = replace( rep("black", terra::nlyr(sim_features)), sim_phylogeny$tip.label %in% r1$feature[r1$met], "red" ) ) # rename the features in the example phylogeny for use with the # multizone data sim_phylogeny$tip.label < feature_names(sim_zones_features) # create targets for a multizone problem. Here, each feature needs a total # of 10 units of habitat to be conserved among the three zones to be # considered adequately conserved targets < tibble::tibble( feature = feature_names(sim_zones_features), zone = list(zone_names(sim_zones_features))[ rep(1, number_of_features(sim_zones_features))], type = rep("absolute", number_of_features(sim_zones_features)), target = rep(10, number_of_features(sim_zones_features)) ) # create a multizone problem with a maximum phylogenetic diversity # objective, where the total expenditure in all zones is 5000. p2 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_max_phylo_div_objective(5000, sim_phylogeny) %>% add_manual_targets(targets) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s2 < solve(p2) # plot solution plot(category_layer(s2), main = "solution", axes = FALSE) # find out which features have their targets met r2 < eval_target_coverage_summary(p2, s2) print(r2, width = Inf) # plot the phylogeny and color the adequately represented features in red plot( sim_phylogeny, main = "adequately represented features", tip.color = replace( rep("black", terra::nlyr(sim_features)), which(r2$met), "red" ) ) # create a multizone problem with a maximum phylogenetic diversity # objective, where each zone has a separate budget. p3 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_max_phylo_div_objective(c(2500, 500, 2000), sim_phylogeny) %>% add_manual_targets(targets) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s3 < solve(p3) # plot solution plot(category_layer(s3), main = "solution", axes = FALSE) # find out which features have their targets met r3 < eval_target_coverage_summary(p3, s3) print(r3, width = Inf) # plot the phylogeny and color the adequately represented features in red plot( sim_phylogeny, main = "adequately represented features", tip.color = replace( rep("black", terra::nlyr(sim_features)), which(r3$met), "red" ) ) ## End(Not run)
## Not run: # load ape package require(ape) # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() sim_phylogeny < get_sim_phylogeny() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_features < get_sim_zones_features() # plot the simulated phylogeny par(mfrow = c(1, 1)) plot(sim_phylogeny, main = "phylogeny") # create problem with a maximum phylogenetic diversity objective, # where each feature needs 10% of its distribution to be secured for # it to be adequately conserved and a total budget of 1900 p1 < problem(sim_pu_raster, sim_features) %>% add_max_phylo_div_objective(1900, sim_phylogeny) %>% add_relative_targets(0.1) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s1 < solve(p1) # plot solution plot(s1, main = "solution", axes = FALSE) # find out which features have their targets met r1 < eval_target_coverage_summary(p1, s1) print(r1, width = Inf) # plot the phylogeny and color the adequately represented features in red plot( sim_phylogeny, main = "adequately represented features", tip.color = replace( rep("black", terra::nlyr(sim_features)), sim_phylogeny$tip.label %in% r1$feature[r1$met], "red" ) ) # rename the features in the example phylogeny for use with the # multizone data sim_phylogeny$tip.label < feature_names(sim_zones_features) # create targets for a multizone problem. Here, each feature needs a total # of 10 units of habitat to be conserved among the three zones to be # considered adequately conserved targets < tibble::tibble( feature = feature_names(sim_zones_features), zone = list(zone_names(sim_zones_features))[ rep(1, number_of_features(sim_zones_features))], type = rep("absolute", number_of_features(sim_zones_features)), target = rep(10, number_of_features(sim_zones_features)) ) # create a multizone problem with a maximum phylogenetic diversity # objective, where the total expenditure in all zones is 5000. p2 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_max_phylo_div_objective(5000, sim_phylogeny) %>% add_manual_targets(targets) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s2 < solve(p2) # plot solution plot(category_layer(s2), main = "solution", axes = FALSE) # find out which features have their targets met r2 < eval_target_coverage_summary(p2, s2) print(r2, width = Inf) # plot the phylogeny and color the adequately represented features in red plot( sim_phylogeny, main = "adequately represented features", tip.color = replace( rep("black", terra::nlyr(sim_features)), which(r2$met), "red" ) ) # create a multizone problem with a maximum phylogenetic diversity # objective, where each zone has a separate budget. p3 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_max_phylo_div_objective(c(2500, 500, 2000), sim_phylogeny) %>% add_manual_targets(targets) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s3 < solve(p3) # plot solution plot(category_layer(s3), main = "solution", axes = FALSE) # find out which features have their targets met r3 < eval_target_coverage_summary(p3, s3) print(r3, width = Inf) # plot the phylogeny and color the adequately represented features in red plot( sim_phylogeny, main = "adequately represented features", tip.color = replace( rep("black", terra::nlyr(sim_features)), which(r3$met), "red" ) ) ## End(Not run)
Set the objective of a conservation planning problem to
maximize the phylogenetic endemism of the features represented in the
solution subject to a budget. This objective is similar to
add_max_phylo_div_objective()
except
that emphasis is placed on representing species with geographically
restricted evolutionary histories, instead representing as much evolutionary
history as possible. This function was inspired by Faith (1992),
Rodrigues et al. (2002), and Rosauer et al. (2009).
add_max_phylo_end_objective(x, budget, tree)
add_max_phylo_end_objective(x, budget, tree)
x 

budget 

tree 

The maximum phylogenetic endemism objective finds the set of
planning units that meets representation targets for a phylogenetic tree
while staying within a fixed budget. If multiple solutions can meet all
targets while staying within budget, the cheapest solution is chosen.
Note that this objective is similar to the maximum
features objective (add_max_features_objective()
) in that it
allows for both a budget and targets to be set for each feature. However,
unlike the maximum feature objective, the aim of this objective is to
maximize the total phylogenetic endemism of the targets met in the
solution, so if multiple targets are provided for a single feature, the
problem will only need to meet a single target for that feature
for the phylogenetic benefit for that feature to be counted when
calculating the phylogenetic endemism of the solution. In other words,
for multizone problems, this objective does not aim to maximize the
phylogenetic endemism in each zone, but rather this objective
aims to maximize the phylogenetic endemism of targets that can be met
through allocating planning units to any of the different zones in a
problem. This can be useful for problems where targets pertain to the total
amount held for each feature across multiple zones. For example,
each feature might have a nonzero amount of suitable habitat in each
planning unit when the planning units are assigned to a (i) not restored,
(ii) partially restored, or (iii) completely restored management zone.
Here each target corresponds to a single feature and can be met through
the total amount of habitat in planning units present to the three
zones.
An updated problem()
object with the objective added to it.
This objective can be expressed mathematically for a set of planning units
($I$
indexed by $i$
) and a set of features ($J$
indexed by $j$
) as:
$\mathit{Maximize} \space \sum_{i = 1}^{I} s \space c_i \space x_i +
\sum_{j = 1}^{J} m_b l_b \frac{1}{a_b} \\
\mathit{subject \space to} \\
\sum_{i = 1}^{I} x_i r_{ij} \geq y_j t_j \forall j \in J \\
m_b \leq y_j \forall j \in T(b) \\
\sum_{i = 1}^{I} x_i c_i \leq B$
Here, $x_i$
is the decisions variable (e.g.,
specifying whether planning unit $i$
has been selected (1) or not
(0)), $r_{ij}$
is the amount of feature $j$
in planning
unit $i$
, $t_j$
is the representation target for feature
$j$
, $y_j$
indicates if the solution has meet
the target $t_j$
for feature $j$
. Additionally,
$T$
represents a phylogenetic tree containing features $j$
and has the branches $b$
associated within lengths $l_b$
.
Each branch $b \in B$
is associated with a total amount
$a_b$
indicating the total geographic extent or amount of habitat.
The $a_b$
variable for a given branch is calculated by summing the
$r_{ij}$
data for all features $j \in J$
that are
associated with the branch. The binary variable $m_b$
denotes if
at least one feature associated with the branch $b$
has met its
representation as indicated by $y_j$
. For brevity, we denote
the features $j$
associated with branch $b$
using
$T(b)$
. Finally, $B$
is the budget allocated for the
solution, $c_i$
is the cost of planning unit $i$
, and
$s$
is a scaling factor used to shrink the costs so that the problem
will return a cheapest solution when there are multiple solutions that
represent the same amount of all features within the budget.
Faith DP (1992) Conservation evaluation and phylogenetic diversity. Biological Conservation, 61: 1–10.
Rodrigues ASL and Gaston KJ (2002) Maximising phylogenetic diversity in the selection of networks of conservation areas. Biological Conservation, 105: 103–111.
Rosauer D, Laffan SW, Crisp, MD, Donnellan SC and Cook LG (2009) Phylogenetic endemism: a new approach for identifying geographical concentrations of evolutionary history. Molecular Ecology, 18: 4061–4072.
See objectives for an overview of all functions for adding objectives.
Also, see targets for an overview of all functions for adding targets, and
add_feature_weights()
to specify weights for different features.
Other objectives:
add_max_cover_objective()
,
add_max_features_objective()
,
add_max_phylo_div_objective()
,
add_max_utility_objective()
,
add_min_largest_shortfall_objective()
,
add_min_set_objective()
,
add_min_shortfall_objective()
## Not run: # load ape package require(ape) # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() sim_phylogeny < get_sim_phylogeny() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_features < get_sim_zones_features() # plot the simulated phylogeny par(mfrow = c(1, 1)) plot(sim_phylogeny, main = "phylogeny") # create problem with a maximum phylogenetic endemism objective, # where each feature needs 10% of its distribution to be secured for # it to be adequately conserved and a total budget of 1900 p1 < problem(sim_pu_raster, sim_features) %>% add_max_phylo_end_objective(1900, sim_phylogeny) %>% add_relative_targets(0.1) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s1 < solve(p1) # plot solution plot(s1, main = "solution", axes = FALSE) # find out which features have their targets met r1 < eval_target_coverage_summary(p1, s1) print(r1, width = Inf) # plot the phylogeny and color the adequately represented features in red plot( sim_phylogeny, main = "adequately represented features", tip.color = replace( rep("black", terra::nlyr(sim_features)), sim_phylogeny$tip.label %in% r1$feature[r1$met], "red" ) ) # rename the features in the example phylogeny for use with the # multizone data sim_phylogeny$tip.label < feature_names(sim_zones_features) # create targets for a multizone problem. Here, each feature needs a total # of 10 units of habitat to be conserved among the three zones to be # considered adequately conserved targets < tibble::tibble( feature = feature_names(sim_zones_features), zone = list(zone_names(sim_zones_features))[ rep(1, number_of_features(sim_zones_features))], type = rep("absolute", number_of_features(sim_zones_features)), target = rep(10, number_of_features(sim_zones_features)) ) # create a multizone problem with a maximum phylogenetic endemism # objective, where the total expenditure in all zones is 5000. p2 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_max_phylo_end_objective(5000, sim_phylogeny) %>% add_manual_targets(targets) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s2 < solve(p2) # plot solution plot(category_layer(s2), main = "solution", axes = FALSE) # find out which features have their targets met r2 < eval_target_coverage_summary(p2, s2) print(r2, width = Inf) # plot the phylogeny and color the adequately represented features in red plot( sim_phylogeny, main = "adequately represented features", tip.color = replace( rep("black", terra::nlyr(sim_features)), which(r2$met), "red" ) ) # create a multizone problem with a maximum phylogenetic endemism # objective, where each zone has a separate budget. p3 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_max_phylo_end_objective(c(2500, 500, 2000), sim_phylogeny) %>% add_manual_targets(targets) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s3 < solve(p3) # plot solution plot(category_layer(s3), main = "solution", axes = FALSE) # find out which features have their targets met r3 < eval_target_coverage_summary(p3, s3) print(r3, width = Inf) # plot the phylogeny and color the adequately represented features in red plot( sim_phylogeny, main = "adequately represented features", tip.color = replace( rep("black", terra::nlyr(sim_features)), which(r3$met), "red" ) ) ## End(Not run)
## Not run: # load ape package require(ape) # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() sim_phylogeny < get_sim_phylogeny() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_features < get_sim_zones_features() # plot the simulated phylogeny par(mfrow = c(1, 1)) plot(sim_phylogeny, main = "phylogeny") # create problem with a maximum phylogenetic endemism objective, # where each feature needs 10% of its distribution to be secured for # it to be adequately conserved and a total budget of 1900 p1 < problem(sim_pu_raster, sim_features) %>% add_max_phylo_end_objective(1900, sim_phylogeny) %>% add_relative_targets(0.1) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s1 < solve(p1) # plot solution plot(s1, main = "solution", axes = FALSE) # find out which features have their targets met r1 < eval_target_coverage_summary(p1, s1) print(r1, width = Inf) # plot the phylogeny and color the adequately represented features in red plot( sim_phylogeny, main = "adequately represented features", tip.color = replace( rep("black", terra::nlyr(sim_features)), sim_phylogeny$tip.label %in% r1$feature[r1$met], "red" ) ) # rename the features in the example phylogeny for use with the # multizone data sim_phylogeny$tip.label < feature_names(sim_zones_features) # create targets for a multizone problem. Here, each feature needs a total # of 10 units of habitat to be conserved among the three zones to be # considered adequately conserved targets < tibble::tibble( feature = feature_names(sim_zones_features), zone = list(zone_names(sim_zones_features))[ rep(1, number_of_features(sim_zones_features))], type = rep("absolute", number_of_features(sim_zones_features)), target = rep(10, number_of_features(sim_zones_features)) ) # create a multizone problem with a maximum phylogenetic endemism # objective, where the total expenditure in all zones is 5000. p2 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_max_phylo_end_objective(5000, sim_phylogeny) %>% add_manual_targets(targets) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s2 < solve(p2) # plot solution plot(category_layer(s2), main = "solution", axes = FALSE) # find out which features have their targets met r2 < eval_target_coverage_summary(p2, s2) print(r2, width = Inf) # plot the phylogeny and color the adequately represented features in red plot( sim_phylogeny, main = "adequately represented features", tip.color = replace( rep("black", terra::nlyr(sim_features)), which(r2$met), "red" ) ) # create a multizone problem with a maximum phylogenetic endemism # objective, where each zone has a separate budget. p3 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_max_phylo_end_objective(c(2500, 500, 2000), sim_phylogeny) %>% add_manual_targets(targets) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s3 < solve(p3) # plot solution plot(category_layer(s3), main = "solution", axes = FALSE) # find out which features have their targets met r3 < eval_target_coverage_summary(p3, s3) print(r3, width = Inf) # plot the phylogeny and color the adequately represented features in red plot( sim_phylogeny, main = "adequately represented features", tip.color = replace( rep("black", terra::nlyr(sim_features)), which(r3$met), "red" ) ) ## End(Not run)
Set the objective of a conservation planning problem to secure as much of the features as possible without exceeding a budget. This objective does not use targets, and feature weights should be used instead to increase the representation of certain features by a solution. Note that this objective does not aim to maximize as much of each feature as possible, and so often results in solutions that are heavily biased towards just a few features.
add_max_utility_objective(x, budget)
add_max_utility_objective(x, budget)
x 

budget 

The maximum utility objective seeks to maximize the overall level of
representation across a suite of conservation features, while keeping cost
within a fixed budget.
Additionally, weights can be used to favor the
representation of certain features over other features (see
add_feature_weights()
). It is essentially calculated as a weighted
sum of the feature data inside the selected planning units.
An updated problem()
object with the objective added to it.
This objective can be expressed mathematically for a set of planning units
($I$
indexed by $i$
) and a set of features ($J$
indexed
by $j$
) as:
$\mathit{Maximize} \space \sum_{i = 1}^{I} s \space c_i \space x_i +
\sum_{j = 1}^{J} a_j w_j \\
\mathit{subject \space to} \\ a_j = \sum_{i = 1}^{I} x_i r_{ij} \space
\forall j \in J \\ \sum_{i = 1}^{I} x_i c_i \leq B$
Here, $x_i$
is the decisions variable (e.g.,
specifying whether planning unit $i$
has been selected (1) or not
(0)), $r_{ij}$
is the amount of feature $j$
in planning
unit $i$
, $A_j$
is the amount of feature $j$
represented in in the solution, and $w_j$
is the weight for
feature $j$
(defaults to 1 for all features; see
add_feature_weights()
to specify weights). Additionally, $B$
is the budget allocated for
the solution, $c_i$
is the cost of planning unit $i$
, and
$s$
is a scaling factor used to shrink the costs so that the problem
will return a cheapest solution when there are multiple solutions that
represent the same amount of all features within the budget.
In early versions (< 3.0.0.0), this function was named as
the add_max_cover_objective
function. It was renamed to avoid
confusion with existing terminology.
See objectives for an overview of all functions for adding objectives.
Also, see add_feature_weights()
to specify weights for different features.
Other objectives:
add_max_cover_objective()
,
add_max_features_objective()
,
add_max_phylo_div_objective()
,
add_max_phylo_end_objective()
,
add_min_largest_shortfall_objective()
,
add_min_set_objective()
,
add_min_shortfall_objective()
## Not run: # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_features < get_sim_zones_features() # create problem with maximum utility objective p1 < problem(sim_pu_raster, sim_features) %>% add_max_utility_objective(5000) %>% add_binary_decisions() %>% add_default_solver(gap = 0, verbose = FALSE) # solve problem s1 < solve(p1) # plot solution plot(s1, main = "solution", axes = FALSE) # create multizone problem with maximum utility objective that # has a single budget for all zones p2 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_max_utility_objective(5000) %>% add_binary_decisions() %>% add_default_solver(gap = 0, verbose = FALSE) # solve problem s2 < solve(p2) # plot solution plot(category_layer(s2), main = "solution", axes = FALSE) # create multizone problem with maximum utility objective that # has separate budgets for each zone p3 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_max_utility_objective(c(1000, 2000, 3000)) %>% add_binary_decisions() %>% add_default_solver(gap = 0, verbose = FALSE) # solve problem s3 < solve(p3) # plot solution plot(category_layer(s3), main = "solution", axes = FALSE) ## End(Not run)
## Not run: # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_features < get_sim_zones_features() # create problem with maximum utility objective p1 < problem(sim_pu_raster, sim_features) %>% add_max_utility_objective(5000) %>% add_binary_decisions() %>% add_default_solver(gap = 0, verbose = FALSE) # solve problem s1 < solve(p1) # plot solution plot(s1, main = "solution", axes = FALSE) # create multizone problem with maximum utility objective that # has a single budget for all zones p2 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_max_utility_objective(5000) %>% add_binary_decisions() %>% add_default_solver(gap = 0, verbose = FALSE) # solve problem s2 < solve(p2) # plot solution plot(category_layer(s2), main = "solution", axes = FALSE) # create multizone problem with maximum utility objective that # has separate budgets for each zone p3 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_max_utility_objective(c(1000, 2000, 3000)) %>% add_binary_decisions() %>% add_default_solver(gap = 0, verbose = FALSE) # solve problem s3 < solve(p3) # plot solution plot(category_layer(s3), main = "solution", axes = FALSE) ## End(Not run)
Set the objective of a conservation planning problem to minimize the largest target shortfall while ensuring that the cost of the solution does not exceed a budget. Note that if the target shortfall for a single feature cannot be decreased beyond a certain point (e.g., because all remaining planning units occupied by that feature are too costly or are locked out), then solutions may only use a small proportion of the specified budget.
add_min_largest_shortfall_objective(x, budget)
add_min_largest_shortfall_objective(x, budget)
x 

budget 

The minimum largest shortfall objective aims to
find the set of planning units that minimize the largest
shortfall for any of the representation targets—that is, the fraction of
each target that remains unmet—for as many features as possible while
staying within a fixed budget. This objective is different from the
minimum shortfall objective (add_min_shortfall_objective()
) because this
objective minimizes the largest (maximum) target shortfall,
whereas the minimum shortfall objective
minimizes the total (weighted sum) of the target shortfalls.
Note that this objective function is not compatible with feature weights
(add_feature_weights()
).
An updated problem()
object with the objective added to it.
This objective can be expressed mathematically for a set of planning units
($I$
indexed by $i$
) and a set of features ($J$
indexed by $j$
) as:
$\mathit{Minimize} \space l \\
\mathit{subject \space to} \\
\sum_{i = 1}^{I} x_i r_{ij} + y_j \geq t_j \forall j \in J \\
l \geq \frac{y_j}{t_j} \forall j \in J \\
\sum_{i = 1}^{I} x_i c_i \leq B$
Here, $x_i$
is the decisions variable (e.g.,
specifying whether planning unit $i$
has been selected (1) or not
(0)), $r_{ij}$
is the amount of feature $j$
in planning
unit $i$
, and $t_j$
is the representation target for feature
$j$
.
Additionally, $y_j$
denotes the target shortfall for
the target $t_j$
for feature $j$
, and
$l$
denotes the largest target shortfall.
Furthermore, $B$
is the budget allocated for the solution,
$c_i$
is the cost of planning unit $i$
. Note that
$y_j$
and $s$
are continuous variables bounded between zero
and infinity.
See objectives for an overview of all functions for adding objectives.
Also, see targets for an overview of all functions for adding targets, and
add_feature_weights()
to specify weights for different features.
Other objectives:
add_max_cover_objective()
,
add_max_features_objective()
,
add_max_phylo_div_objective()
,
add_max_phylo_end_objective()
,
add_max_utility_objective()
,
add_min_set_objective()
,
add_min_shortfall_objective()
## Not run: # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_features < get_sim_zones_features() # create problem with minimum largest shortfall objective p1 < problem(sim_pu_raster, sim_features) %>% add_min_largest_shortfall_objective(1800) %>% add_relative_targets(0.1) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s1 < solve(p1) # plot solution plot(s1, main = "solution", axes = FALSE) # create multizone problem with minimum largest shortfall objective, # with 10% representation targets for each feature, and set # a budget such that the total maximum expenditure in all zones # cannot exceed 1800 p2 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_min_largest_shortfall_objective(1800) %>% add_relative_targets(matrix(0.1, ncol = 3, nrow = 5)) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s2 < solve(p2) # plot solution plot(category_layer(s2), main = "solution", axes = FALSE) # create multizone problem with minimum largest shortfall objective, # with 10% representation targets for each feature, and set # separate budgets of 1800 for each management zone p3 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_min_largest_shortfall_objective(c(1800, 1800, 1800)) %>% add_relative_targets(matrix(0.1, ncol = 3, nrow = 5)) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s3 < solve(p3) # plot solution plot(category_layer(s3), main = "solution", axes = FALSE) ## End(Not run)
## Not run: # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_features < get_sim_zones_features() # create problem with minimum largest shortfall objective p1 < problem(sim_pu_raster, sim_features) %>% add_min_largest_shortfall_objective(1800) %>% add_relative_targets(0.1) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s1 < solve(p1) # plot solution plot(s1, main = "solution", axes = FALSE) # create multizone problem with minimum largest shortfall objective, # with 10% representation targets for each feature, and set # a budget such that the total maximum expenditure in all zones # cannot exceed 1800 p2 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_min_largest_shortfall_objective(1800) %>% add_relative_targets(matrix(0.1, ncol = 3, nrow = 5)) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s2 < solve(p2) # plot solution plot(category_layer(s2), main = "solution", axes = FALSE) # create multizone problem with minimum largest shortfall objective, # with 10% representation targets for each feature, and set # separate budgets of 1800 for each management zone p3 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_min_largest_shortfall_objective(c(1800, 1800, 1800)) %>% add_relative_targets(matrix(0.1, ncol = 3, nrow = 5)) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s3 < solve(p3) # plot solution plot(category_layer(s3), main = "solution", axes = FALSE) ## End(Not run)
Set the objective of a conservation planning problem to minimize the cost of the solution whilst ensuring that all targets are met. This objective is similar to that used in Marxan and is detailed in Rodrigues et al. (2000).
add_min_set_objective(x)
add_min_set_objective(x)
x 

The minimum set objective – in the the context of systematic reserve design – seeks to find the set of planning units that minimizes the overall cost of a reserve network, while meeting a set of representation targets for the conservation features. This objective is equivalent to a simplified Marxan reserve design problem with the Boundary Length Modifier (BLM) set to zero. The difference between this objective and the Marxan software is that the targets for the features will always be met (and as such it does not use Species Penalty Factors).
An updated problem()
object with the objective added to it.
This objective can be expressed
mathematically for a set of planning units ($I$
indexed by
$i$
) and a set of features ($J$
indexed by $j$
) as:
$\mathit{Minimize} \space \sum_{i = 1}^{I} x_i c_i \\
\mathit{subject \space to} \\
\sum_{i = 1}^{I} x_i r_{ij} \geq T_j \space \forall \space j \in J$
Here, $x_i$
is the decisions variable (e.g.,
specifying whether planning unit $i$
has been selected (1) or not
(0)), $c_i$
is the cost of planning unit $i$
,
$r_{ij}$
is the amount of feature $j$
in planning unit
$i$
, and $T_j$
is the target for feature $j$
. The
first term is the objective function and the second is the set of
constraints. In words this says find the set of planning units that meets
all the representation targets while minimizing the overall cost.
Rodrigues AS, Cerdeira OJ, and Gaston KJ (2000) Flexibility, efficiency, and accountability: adapting reserve selection algorithms to more complex conservation problems. Ecography, 23: 565–574.
See objectives for an overview of all functions for adding objectives. Also see targets for an overview of all functions for adding targets.
Other objectives:
add_max_cover_objective()
,
add_max_features_objective()
,
add_max_phylo_div_objective()
,
add_max_phylo_end_objective()
,
add_max_utility_objective()
,
add_min_largest_shortfall_objective()
,
add_min_shortfall_objective()
## Not run: # set seed for reproducibility set.seed(500) # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_features < get_sim_zones_features() # create minimal problem with minimum set objective p1 < problem(sim_pu_raster, sim_features) %>% add_min_set_objective() %>% add_relative_targets(0.1) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s1 < solve(p1) # plot solution plot(s1, main = "solution", axes = FALSE) # create multizone problem with minimum set objective targets_matrix < matrix(rpois(15, 1), nrow = 5, ncol = 3) p2 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_min_set_objective() %>% add_absolute_targets(targets_matrix) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s2 < solve(p2) # plot solution plot(category_layer(s2), main = "solution", axes = FALSE) ## End(Not run)
## Not run: # set seed for reproducibility set.seed(500) # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_features < get_sim_zones_features() # create minimal problem with minimum set objective p1 < problem(sim_pu_raster, sim_features) %>% add_min_set_objective() %>% add_relative_targets(0.1) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s1 < solve(p1) # plot solution plot(s1, main = "solution", axes = FALSE) # create multizone problem with minimum set objective targets_matrix < matrix(rpois(15, 1), nrow = 5, ncol = 3) p2 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_min_set_objective() %>% add_absolute_targets(targets_matrix) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s2 < solve(p2) # plot solution plot(category_layer(s2), main = "solution", axes = FALSE) ## End(Not run)
Set the objective of a conservation planning problem to minimize the overall shortfall for as many targets as possible while ensuring that the cost of the solution does not exceed a budget.
add_min_shortfall_objective(x, budget)
add_min_shortfall_objective(x, budget)
x 

budget 

The minimum shortfall objective aims to
find the set of planning units that minimize the overall
(weighted sum) shortfall for the
representation targets—that is, the fraction of each target that
remains unmet—for as many features as possible while staying within a
fixed budget (inspired by Table 1, equation IV, Arponen et al.
2005). Additionally, weights can be used
to favor the representation of certain features over other features (see
add_feature_weights()
.
An updated problem()
object with the objective added to it.
This objective can be expressed mathematically for a set of planning units
($I$
indexed by $i$
) and a set of features ($J$
indexed
by $j$
) as:
$\mathit{Minimize} \space \sum_{j = 1}^{J} w_j \frac{y_j}{t_j} \\
\mathit{subject \space to} \\
\sum_{i = 1}^{I} x_i r_{ij} + y_j \geq t_j \forall j \in J \\
\sum_{i = 1}^{I} x_i c_i \leq B$
Here, $x_i$
is the decisions variable (e.g.,
specifying whether planning unit $i$
has been selected (1) or not
(0)), $r_{ij}$
is the amount of feature $j$
in planning
unit $i$
, $t_j$
is the representation target for feature
$j$
, $y_j$
denotes the representation shortfall for
the target $t_j$
for feature $j$
, and $w_j$
is the
weight for feature $j$
(defaults to 1 for all features; see
add_feature_weights()
to specify weights). Additionally,
$B$
is the budget allocated for the solution, $c_i$
is the
cost of planning unit $i$
. Note that $y_j$
is a continuous
variable bounded between zero and infinity, and denotes the shortfall
for target $j$
.
Arponen A, Heikkinen RK, Thomas CD, and Moilanen A (2005) The value of biodiversity in reserve selection: representation, species weighting, and benefit functions. Conservation Biology, 19: 2009–2014.
See objectives for an overview of all functions for adding objectives.
Also, see targets for an overview of all functions for adding targets, and
add_feature_weights()
to specify weights for different features.
Other objectives:
add_max_cover_objective()
,
add_max_features_objective()
,
add_max_phylo_div_objective()
,
add_max_phylo_end_objective()
,
add_max_utility_objective()
,
add_min_largest_shortfall_objective()
,
add_min_set_objective()
## Not run: # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_features < get_sim_zones_features() # create problem with minimum shortfall objective p1 < problem(sim_pu_raster, sim_features) %>% add_min_shortfall_objective(1800) %>% add_relative_targets(0.1) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s1 < solve(p1) # plot solution plot(s1, main = "solution", axes = FALSE) # create multizone problem with minimum shortfall objective, # with 10% representation targets for each feature, and set # a budget such that the total maximum expenditure in all zones # cannot exceed 3000 p2 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_min_shortfall_objective(3000) %>% add_relative_targets(matrix(0.1, ncol = 3, nrow = 5)) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s2 < solve(p2) # plot solution plot(category_layer(s2), main = "solution", axes = FALSE) # create multizone problem with minimum shortfall objective, # with 10% representation targets for each feature, and set # separate budgets for each management zone p3 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_min_shortfall_objective(c(3000, 3000, 3000)) %>% add_relative_targets(matrix(0.1, ncol = 3, nrow = 5)) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s3 < solve(p3) # plot solution plot(category_layer(s3), main = "solution", axes = FALSE) ## End(Not run)
## Not run: # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_features < get_sim_zones_features() # create problem with minimum shortfall objective p1 < problem(sim_pu_raster, sim_features) %>% add_min_shortfall_objective(1800) %>% add_relative_targets(0.1) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s1 < solve(p1) # plot solution plot(s1, main = "solution", axes = FALSE) # create multizone problem with minimum shortfall objective, # with 10% representation targets for each feature, and set # a budget such that the total maximum expenditure in all zones # cannot exceed 3000 p2 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_min_shortfall_objective(3000) %>% add_relative_targets(matrix(0.1, ncol = 3, nrow = 5)) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s2 < solve(p2) # plot solution plot(category_layer(s2), main = "solution", axes = FALSE) # create multizone problem with minimum shortfall objective, # with 10% representation targets for each feature, and set # separate budgets for each management zone p3 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_min_shortfall_objective(c(3000, 3000, 3000)) %>% add_relative_targets(matrix(0.1, ncol = 3, nrow = 5)) %>% add_binary_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s3 < solve(p3) # plot solution plot(category_layer(s3), main = "solution", axes = FALSE) ## End(Not run)
Add constraints to a conservation planning problem to ensure that all selected planning units in the solution have at least a certain number of neighbors that are also selected in the solution.
## S4 method for signature 'ConservationProblem,ANY,ANY,ANY,ANY' add_neighbor_constraints(x, k, clamp, zones, data) ## S4 method for signature 'ConservationProblem,ANY,ANY,ANY,data.frame' add_neighbor_constraints(x, k, clamp, zones, data) ## S4 method for signature 'ConservationProblem,ANY,ANY,ANY,matrix' add_neighbor_constraints(x, k, clamp, zones, data) ## S4 method for signature 'ConservationProblem,ANY,ANY,ANY,array' add_neighbor_constraints(x, k, clamp, zones, data)
## S4 method for signature 'ConservationProblem,ANY,ANY,ANY,ANY' add_neighbor_constraints(x, k, clamp, zones, data) ## S4 method for signature 'ConservationProblem,ANY,ANY,ANY,data.frame' add_neighbor_constraints(x, k, clamp, zones, data) ## S4 method for signature 'ConservationProblem,ANY,ANY,ANY,matrix' add_neighbor_constraints(x, k, clamp, zones, data) ## S4 method for signature 'ConservationProblem,ANY,ANY,ANY,array' add_neighbor_constraints(x, k, clamp, zones, data)
x 

k 

clamp 

zones 

data 

This function uses neighborhood data to identify solutions that surround planning units with a minimum number of neighbors. It was inspired by the mathematical formulations detailed in Billionnet (2013) and Beyer et al. (2016).
An updated problem()
object with the constraints added to it.
The argument to data
can be specified using the following formats:
data
as a NULL
valueneighborhood data should be calculated
automatically
using the adjacency_matrix()
function. This is the default
argument. Note that the neighborhood data must be manually defined
using one of the other formats below when the planning unit data
in the argument to x
is not spatially referenced (e.g.,
in data.frame
or numeric
format).
data
as a matrix
/Matrix
objectwhere rows and columns represent
different planning units and the value of each cell indicates if the
two planning units are neighbors or not. Cell values should be binary
numeric
values (i.e., one or zero). Cells that occur along the
matrix diagonal have no effect on the solution at all because each
planning unit cannot be a neighbor with itself.
data
as a data.frame
objectcontaining columns that are named
"id1"
, "id2"
, and "boundary"
. Here, each row
denotes the connectivity between two planning units following the
Marxan format. The "boundary"
column should contain
binary numeric
values that indicate if the two planning units
specified in the "id1"
and "id2"
columns are neighbors
or not. This data can be used to describe symmetric or
asymmetric relationships between planning units. By default,
input data is assumed to be symmetric unless asymmetric data is
also included (e.g., if data is present for planning units 2 and 3, then
the same amount of connectivity is expected for planning units 3 and 2,
unless connectivity data is also provided for planning units 3 and 2).
If the argument to x
contains multiple zones, then the
"zone1"
and "zone2"
columns can optionally be provided to manually
specify if the neighborhood data pertain to specific zones. The
"zone1"
and "zone2"
columns should contain the character
names of the zones. If the columns "zone1"
and "zone2"
are present, then the argument to zones
must be NULL
.
data
as an array
objectcontaining fourdimensions where binary
numeric
values indicate if planning unit should be treated
as being neighbors with every other planning unit when they
are allocated to every combination of management zone. The first two
dimensions (i.e., rows and columns) correspond to the planning units,
and second two dimensions correspond to the management zones. For
example, if the argument to data
had a value of 1 at the index
data[1, 2, 3, 4]
this would indicate that planning unit 1 and
planning unit 2 should be treated as neighbors when they are
allocated to zones 3 and 4 respectively.
Beyer HL, Dujardin Y, Watts ME, and Possingham HP (2016) Solving conservation planning problems with integer linear programming. Ecological Modelling, 228: 14–22.
Billionnet A (2013) Mathematical optimization ideas for biodiversity conservation. European Journal of Operational Research, 231: 514–534.
Other constraints:
add_contiguity_constraints()
,
add_feature_contiguity_constraints()
,
add_linear_constraints()
,
add_locked_in_constraints()
,
add_locked_out_constraints()
,
add_mandatory_allocation_constraints()
,
add_manual_bounded_constraints()
,
add_manual_locked_constraints()
## Not run: # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_features < get_sim_zones_features() # create minimal problem p1 < problem(sim_pu_raster, sim_features) %>% add_min_set_objective() %>% add_relative_targets(0.1) %>% add_default_solver(verbose = FALSE) # create problem with constraints that require 1 neighbor # and neighbors are defined using a rookstyle neighborhood p2 < p1 %>% add_neighbor_constraints(1) # create problem with constraints that require 2 neighbor # and neighbors are defined using a rookstyle neighborhood p3 < p1 %>% add_neighbor_constraints(2) # create problem with constraints that require 3 neighbor # and neighbors are defined using a queenstyle neighborhood p4 < p1 %>% add_neighbor_constraints( 3, data = adjacency_matrix(sim_pu_raster, directions = 8) ) # solve problems s1 < terra::rast(list(solve(p1), solve(p2), solve(p3), solve(p4))) names(s1) < c("basic solution", "1 neighbor", "2 neighbors", "3 neighbors") # plot solutions plot(s1, axes = FALSE) # create minimal problem with multiple zones p5 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_min_set_objective() %>% add_relative_targets(matrix(0.1, ncol = 3, nrow = 5)) %>% add_default_solver(verbose = FALSE) # create problem where selected planning units require at least 2 neighbors # for each zone and planning units are only considered neighbors if they # are allocated to the same zone z6 < diag(3) print(z6) p6 < p5 %>% add_neighbor_constraints(rep(2, 3), zones = z6) # create problem where the planning units in zone 1 don't explicitly require # any neighbors, planning units in zone 2 require at least 1 neighbors, and # planning units in zone 3 require at least 2 neighbors. As before, planning # units are still only considered neighbors if they are allocated to the # same zone p7 < p5 %>% add_neighbor_constraints(c(0, 1, 2), zones = z6) # create problem given the same constraints as outlined above, except # that when determining which selected planning units are neighbors, # planning units that are allocated to zone 1 and zone 2 can also treated # as being neighbors with each other z8 < diag(3) z8[1, 2] < 1 z8[2, 1] < 1 print(z8) p8 < p5 %>% add_neighbor_constraints(c(0, 1, 2), zones = z8) # solve problems s2 < list(p5, p6, p7, p8) s2 < lapply(s2, solve) s2 < lapply(s2, category_layer) s2 < terra::rast(s2) names(s2) < c("basic problem", "p6", "p7", "p8") # plot solutions plot(s2, main = names(s2), axes = FALSE) ## End(Not run)
## Not run: # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_features < get_sim_zones_features() # create minimal problem p1 < problem(sim_pu_raster, sim_features) %>% add_min_set_objective() %>% add_relative_targets(0.1) %>% add_default_solver(verbose = FALSE) # create problem with constraints that require 1 neighbor # and neighbors are defined using a rookstyle neighborhood p2 < p1 %>% add_neighbor_constraints(1) # create problem with constraints that require 2 neighbor # and neighbors are defined using a rookstyle neighborhood p3 < p1 %>% add_neighbor_constraints(2) # create problem with constraints that require 3 neighbor # and neighbors are defined using a queenstyle neighborhood p4 < p1 %>% add_neighbor_constraints( 3, data = adjacency_matrix(sim_pu_raster, directions = 8) ) # solve problems s1 < terra::rast(list(solve(p1), solve(p2), solve(p3), solve(p4))) names(s1) < c("basic solution", "1 neighbor", "2 neighbors", "3 neighbors") # plot solutions plot(s1, axes = FALSE) # create minimal problem with multiple zones p5 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_min_set_objective() %>% add_relative_targets(matrix(0.1, ncol = 3, nrow = 5)) %>% add_default_solver(verbose = FALSE) # create problem where selected planning units require at least 2 neighbors # for each zone and planning units are only considered neighbors if they # are allocated to the same zone z6 < diag(3) print(z6) p6 < p5 %>% add_neighbor_constraints(rep(2, 3), zones = z6) # create problem where the planning units in zone 1 don't explicitly require # any neighbors, planning units in zone 2 require at least 1 neighbors, and # planning units in zone 3 require at least 2 neighbors. As before, planning # units are still only considered neighbors if they are allocated to the # same zone p7 < p5 %>% add_neighbor_constraints(c(0, 1, 2), zones = z6) # create problem given the same constraints as outlined above, except # that when determining which selected planning units are neighbors, # planning units that are allocated to zone 1 and zone 2 can also treated # as being neighbors with each other z8 < diag(3) z8[1, 2] < 1 z8[2, 1] < 1 print(z8) p8 < p5 %>% add_neighbor_constraints(c(0, 1, 2), zones = z8) # solve problems s2 < list(p5, p6, p7, p8) s2 < lapply(s2, solve) s2 < lapply(s2, category_layer) s2 < terra::rast(s2) names(s2) < c("basic problem", "p6", "p7", "p8") # plot solutions plot(s2, main = names(s2), axes = FALSE) ## End(Not run)
Add a proportion decision to a conservation planning problem. This is a relaxed decision where a part of a planning unit can be prioritized, as opposed to the entire planning unit. Typically, this decision has the assumed action of buying a fraction of a planning unit to include in decisions will solve much faster than problems that use binarytype decisions.
add_proportion_decisions(x)
add_proportion_decisions(x)
x 

Conservation planning problems involve making decisions on planning
units. These decisions are then associated with actions (e.g., turning a
planning unit into a protected area). Only a
single decision should be added to a problem()
object.
Note that if multiple decisions are added to an object, then the
last one to be added will be used.
An updated problem()
object with the decisions added to it.
See decisions for an overview of all functions for adding decisions.
Other decisions:
add_binary_decisions()
,
add_semicontinuous_decisions()
## Not run: # set seed for reproducibility set.seed(500) # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_features < get_sim_zones_features() # create minimal problem with proportion decisions p1 < problem(sim_pu_raster, sim_features) %>% add_min_set_objective() %>% add_relative_targets(0.1) %>% add_proportion_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s1 < solve(p1) # plot solutions plot(s1, main = "solution", axes = FALSE) # build multizone conservation problem with proportion decisions p2 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_min_set_objective() %>% add_relative_targets(matrix(runif(15, 0.1, 0.2), nrow = 5, ncol = 3)) %>% add_proportion_decisions() %>% add_default_solver(verbose = FALSE) # solve the problem s2 < solve(p2) # print solution print(s2) # plot solution # panels show the proportion of each planning unit allocated to each zone plot(s2, axes = FALSE) ## End(Not run)
## Not run: # set seed for reproducibility set.seed(500) # load data sim_pu_raster < get_sim_pu_raster() sim_features < get_sim_features() sim_zones_pu_raster < get_sim_zones_pu_raster() sim_zones_features < get_sim_zones_features() # create minimal problem with proportion decisions p1 < problem(sim_pu_raster, sim_features) %>% add_min_set_objective() %>% add_relative_targets(0.1) %>% add_proportion_decisions() %>% add_default_solver(verbose = FALSE) # solve problem s1 < solve(p1) # plot solutions plot(s1, main = "solution", axes = FALSE) # build multizone conservation problem with proportion decisions p2 < problem(sim_zones_pu_raster, sim_zones_features) %>% add_min_set_objective() %>% add_relative_targets(matrix(runif(15, 0.1, 0.2), nrow = 5, ncol = 3)) %>% add_proportion_decisions() %>% add_default_solver(verbose = FALSE) # solve the problem s2 < solve(p2