Package 'prioritizr'

Description

Add targets to a conservation planning problem expressed as the same values as the underlying feature data (ignoring any specified feature units). For example, setting a target of 10 for a feature specifies that a solution should ideally select a set of planning units that contain a total (summed) value of, at least, 10 for the feature.

Usage

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)

Arguments

Details

This function is used to set targets for each feature (separately). For problems associated with a single management zone, this function may be useful to specify individual targets for each feature. For problems associated with multiple management zones, this function can also be used to specify a target for each feature within each zone (separately). For example, this may be useful in planning exercises where it is important to ensure that some of the features are adequately represented by multiple zones. For example, in a marine spatial planning exercise, it may be important for some features (e.g., commercial important fish species) to be adequately represented by a conservation zone for ensuring their long-term persistence, and also by a fishing zone to for ensure food security. For greater flexibility in target setting (such as setting targets that can be met through the allocation of multiple zones), see the add_manual_targets() function.

Value

An updated problem() object with the targets added to it.

Targets format

The targets for a problem can be specified using the following formats.

targets as a numeric vector

Here a target value is specified for each feature. Additionally, for convenience, this format can be a single numeric value to assign the same target to each feature. Note that this format cannot be used to specify targets if x has multiple zones.

targets as a matrix object

Here a target value is specified for each feature in each zone. Each row corresponds to a different feature in x, each column corresponds to a different zone in x, and each cell contains a target value for representing a given feature in a given zone.

targets as a character vector

Here target values are specified based on the name(s) of column(s) in the feature data in x. This format can only be used when the feature data in x is a sf::st_sf() or data.frame object. If x has a single zone, then targets must contain a single character value. Otherwise, if x has multiple zones, then targets must contain a character value for each zone in x.

Target setting

Many conservation planning problems require targets. Targets are used to specify the minimum amount, or proportion, of a feature's spatial distribution that should ideally be protected. This is important so that the optimization process can weigh the merits and trade-offs between improving the representation of one feature over another feature. Although it can be challenging to set meaningful targets, this is a critical step for ensuring that prioritizations meet the stakeholder objectives that underpin a prioritization exercise (Carwardine et al. 2009). In other words, targets play an important role in ensuring that a priority setting process is properly tuned according to stakeholder requirements. For example, targets provide a mechanism for ensuring that a prioritization secures enough habitat to promote the long-term persistence of each threatened species, culturally important species, or economically important ecosystem services under consideration. Since there is often uncertainty regarding stakeholder objectives (e.g., how much habitat should be protected for a given species) or the influence of particular target on a prioritization (e.g., how would setting a 90% or 100% for a threatened species alter priorities), it is often useful to generate and compare a suite of prioritizations based on different target scenarios.

References

Carwardine J, Klein CJ, Wilson KA, Pressey RL, Possingham HP (2009) Hitting the target and missing the point: target‐based conservation planning in context. Conservation Letters, 2: 4–11.

See Also

Other functions for adding targets: add_auto_targets(), add_group_targets(), add_manual_targets(), add_relative_targets()

Examples

# 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 (cell 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 (cell values correspond to zone identifiers)
plot(category_layer(s5), main = "varying targets", axes = FALSE)

Description

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.

Usage

## 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)

Arguments

Details

This function adds penalties to a 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).

Value

An updated problem() object with the penalties added to it.

Mathematical formulation

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 penalty, $D$ represent data, and $W$ represent zones.

If data is specified as a matrix or Matrix object, then the penalties are calculated as follows.

$\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 data is specified as an array object, then the penalties are calculated as follows.

$\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$. Additionally, to linearize the problem, the expression $X_{iz} \times X_{jy}$ is modeled using a set of continuous variables (bounded between 0 and 1) based on Beyer et al. (2016).

Data format

The following formats can be used to specify data.

data as a matrix/Matrix object

Here rows and columns correspond to different planning units and cell values denote the strength of connectivity between two planning units. Cells that occur along the matrix diagonal are treated as weights which indicate that planning units are more desirable in the solution. With this format, zones can be used to control the strength of connectivity between planning units in different zones. Note that the default for zones is to treat planning units allocated to different zones as having zero connectivity.

data as a data.frame object

Here rows correspond to a pair of planning units and columns provide information about each pair of planning units. In particular, data must have 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 strength of connectivity between them (following the Marxan format). If x has 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 particular zones. Note that if the "zone1" and "zone2" columns are present, then zones must be NULL.

data as an array object

Here a four-dimension array is used to specify connectivity data, 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.

References

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.

Beyer HL, Dujardin Y, Watts ME, and Possingham HP (2016) Solving conservation planning problems with integer linear programming. Ecological Modelling, 228: 14–22.

See Also

See penalties for an overview of all functions for adding penalties. Also see add_connectivity_penalties() to account for symmetric connectivity between planning units. Also see calibrate_cohon_penalty() for assistance with selecting an appropriate penalty value.

Other functions for adding penalties: add_boundary_penalties(), add_connectivity_penalties(), add_cost_penalties(), add_feature_weights(), add_linear_penalties(), add_neighbor_penalties()

Examples

# 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 x-coordinate
      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 multi-zone 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 multi-zone 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)

Description

Add targets to a conservation planning problem based on a particular target setting method.

Usage

## S4 method for signature 'ConservationProblem,character'
add_auto_targets(x, method)

## S4 method for signature 'ConservationProblem,list'
add_auto_targets(x, method)

## S4 method for signature 'ConservationProblem,TargetMethod'
add_auto_targets(x, method)

Arguments

Value

An updated problem() object with the targets added to it.

Target methods

A variety of options are available for specifying target setting methods.

method is a character value

This option involves specifying a target setting method based on its name. It is particularly useful when all features should be assigned targets based on the same method, and the method should rely on default parameters. For example, if x was a problem(), then the following code could be used to specify that all features should have their targets calculated based on Jung et al. (2021) with default parameters.

x %>% add_auto_targets(method = "jung"))

The following character values can be used to specify target setting methods: "jung" (per Jung et al. 2021), "polak" (per Polak et al. 2016), "rodrigues" (per Rodrigues et al. 2004), "ward" (per Ward et al. 2025), and "watson" (per Watson et al. 2010). (2010). Additionally, the following values can be used to set targets based on criteria from the IUCN Red List of Threatened Species (IUCN 2025): "rl_species_VU_A1_B1" "rl_species_EN_A1_B1", "rl_species_CR_A1_B1", "rl_species_VU_A1_B2" "rl_species_EN_A1_B2", "rl_species_CR_A1_B2", "rl_species_VU_A2_B1" "rl_species_EN_A2_B1", "rl_species_CR_A2_B1", "rl_species_VU_A2_B2" "rl_species_EN_A2_B2", "rl_species_CR_A2_B2", "rl_species_VU_A3_B1" "rl_species_EN_A3_B1", "rl_species_CR_A3_B1", "rl_species_VU_A3_B2" "rl_species_EN_A3_B2", "rl_species_CR_A3_B2", "rl_species_VU_A4_B1" "rl_species_EN_A4_B1", "rl_species_CR_A4_B1", "rl_species_VU_A4_B2" "rl_species_EN_A4_B2", and "rl_species_CR_A4_B2". Furthermore, the following values can be used to set targets based on criteria from the IUCN Red List of Ecosystems (IUCN 2024): "rl_ecosystem_VU_A1_B1", "rl_ecosystem_EN_A1_B1", "rl_ecosystem_CR_A1_B1", "rl_ecosystem_VU_A1_B2", "rl_ecosystem_EN_A1_B2", "rl_ecosystem_CR_A1_B2", "rl_ecosystem_VU_a2a_B1", "rl_ecosystem_EN_a2a_B1", "rl_ecosystem_CR_a2a_B1", "rl_ecosystem_VU_a2a_B2", "rl_ecosystem_EN_a2a_B2", "rl_ecosystem_CR_a2a_B2", "rl_ecosystem_VU_a2b_B1", "rl_ecosystem_EN_a2b_B1", "rl_ecosystem_CR_a2b_B1", "rl_ecosystem_VU_a2b_B2", "rl_ecosystem_EN_a2b_B2", "rl_ecosystem_CR_a2b_B2", "rl_ecosystem_VU_A3_B1", "rl_ecosystem_EN_A3_B1", "rl_ecosystem_CR_A3_B1", "rl_ecosystem_VU_A3_B2", "rl_ecosystem_EN_A3_B2", and "rl_ecosystem_CR_A3_B2".' For convenience, these options can also be specified with lower case letters (e.g., "rl_species_VU_A1_B1" may alternatively be specified as "rl_species_vu_a1_b1"). Note that target setting methods that require additional data or parameters cannot be specified with character values.

method is a character vector

This option involves specifying a target setting method for each feature based on the name of the method. It is particularly useful when particular features should be assigned targets based on different methods, and all methods rely on default parameters. For example, if x was a problem() with three features, then the following code could be used to specify a target based on Jung et al. (2021) for the first feature, target based on Polak et al. (2015) for the second feature, and target based on Jung et al. (2021) for the third feature (note that all targets are based on default parameters).

x %>% add_auto_targets(method = c("jung", "polak", "jung"))

Note that method must specify a value for each feature in x, and the order of the method values should follow the order of the features in x (i.e, per feature_names(x)). For details on the available method values, please refer to the first option where method is a character value.

method is an object for specifying a method

This option involves specifying a target setting method based on an object. It is particularly useful when all features should be assigned targets based on the same method, and the parameters for calculating the targets should be customized. For example, if x was a problem(), then the following could could be used to specify that all features should have their targets calculated based on Jung et al. (2021) with an uplift of 5%.

x %>% add_auto_targets(method = jung_targets(prop_uplift = 0.05))

The following functions can be used to create an object for specifying a target setting method: spec_absolute_targets(), spec_relative_targets(), spec_area_targets(), spec_interp_absolute_targets(), spec_interp_area_targets(), spec_duran_targets(), spec_jung_targets(), spec_polak_targets(), spec_pop_size_targets(), spec_rl_ecosystem_targets(), spec_rl_species_targets(), spec_rodrigues_targets(), spec_rule_targets(), spec_ward_targets(), spec_watson_targets(), spec_wilson_targets(), spec_min_targets(), and spec_max_targets().

method is a list

This option involves specify a target setting method based on a list of objects. It is particularly useful when particular features should be assigned different targets based on different methods, and at least one of the methods requires customized parameters. For example, if x was a problem() with three features, then the following code could be used to specify (i) a target for the first feature based on Jung et al. (2021) with customized parameters (i.e, jung_targets(prop_uplift = 0.05)), (ii) a target for the second feature based on Polak et al. (2015) with default parameters via the name method name (i.e., "polak"), and (iii) a target based on Jung et al. (2021) with default parameters via an object (i.e., jung_targets()).

x %>% add_auto_targets(
  method = list(jung_targets(prop_uplift = 0.05), "polak", jung_targets())
)

Note that method must specify a value for each feature in x, and the order of the method values should follow the order of the features in x (i.e, per feature_names(x)). For details on the available method values, please refer to the previous options.

Data calculations

Many of the functions for specifying target setting methods involve calculating targets based on the spatial extent of the features in x (e.g., spec_jung_targets(), [spec_rodrigues_targets(), and others). Although this function for adding targets can be readily applied to problem() objects that have the feature data provided as a terra::rast() object, you will need to specify the spatial units for the features when building a problem() object if the feature data are provided in a different format. In particular, if the feature data are provided as a data.frame or character vector, then you will need to specify feature_units when using the problem() function. See the Examples section below for a demonstration of using the feature_units parameter.

Target setting

Many conservation planning problems require targets. Targets are used to specify the minimum amount, or proportion, of a feature's spatial distribution that should ideally be protected. This is important so that the optimization process can weigh the merits and trade-offs between improving the representation of one feature over another feature. Although it can be challenging to set meaningful targets, this is a critical step for ensuring that prioritizations meet the stakeholder objectives that underpin a prioritization exercise (Carwardine et al. 2009). In other words, targets play an important role in ensuring that a priority setting process is properly tuned according to stakeholder requirements. For example, targets provide a mechanism for ensuring that a prioritization secures enough habitat to promote the long-term persistence of each threatened species, culturally important species, or economically important ecosystem services under consideration. Since there is often uncertainty regarding stakeholder objectives (e.g., how much habitat should be protected for a given species) or the influence of particular target on a prioritization (e.g., how would setting a 90% or 100% for a threatened species alter priorities), it is often useful to generate and compare a suite of prioritizations based on different target scenarios.

References

Carwardine J, Klein CJ, Wilson KA, Pressey RL, Possingham HP (2009) Hitting the target and missing the point: target‐based conservation planning in context. Conservation Letters, 2: 4–11.

Jung M, Arnell A, de Lamo X, García-Rangel S, Lewis M, Mark J, Merow C, Miles L, Ondo I, Pironon S, Ravilious C, Rivers M, Schepaschenko D, Tallowin O, van Soesbergen A, Govaerts R, Boyle BL, Enquist BJ, Feng X, Gallagher R, Maitner B, Meiri S, Mulligan M, Ofer G, Roll U, Hanson JO, Jetz W, Di Marco M, McGowan J, Rinnan DS, Sachs JD, Lesiv M, Adams VM, Andrew SC, Burger JR, Hannah L, Marquet PA, McCarthy JK, Morueta-Holme N, Newman EA, Park DS, Roehrdanz PR, Svenning J-C, Violle C, Wieringa JJ, Wynne G, Fritz S, Strassburg BBN, Obersteiner M, Kapos V, Burgess N, Schmidt- Traub G, Visconti P (2021) Areas of global importance for conserving terrestrial biodiversity, carbon and water. Nature Ecology and Evolution, 5:1499–1509.

IUCN (2025) The IUCN Red List of Threatened Species. Version 2025-1. Available at https://www.iucnredlist.org. Accessed on 23 July 2025.

IUCN (2024). Guidelines for the application of IUCN Red List of Ecosystems Categories and Criteria, Version 2.0. Keith DA, Ferrer-Paris JR, Ghoraba SMM, Henriksen S, Monyeki M, Murray NJ, Nicholson E, Rowland J, Skowno A, Slingsby JA, Storeng AB, Valderrábano M, Zager I (Eds.). Gland, Switzerland: IUCN.

Polak T, Watson JEM, Fuller RA, Joseph LN, Martin TG, Possingham HP, Venter O, Carwardine J (2015) Efficient expansion of global protected areas requires simultaneous planning for species and ecosystems. Royal Society Open Science, 2: 150107.

Ward M, Possingham HP, Wintle BA, Woinarski JCZ, Marsh JR, Chapple DG, Lintermans M, Scheele BC, Whiterod NS, Hoskin CJ, Aska B, Yong C, Tulloch A, Stewart R, Watson JEM (2025) The estimated cost of preventing extinction and progressing recovery for Australia's priority threatened species. Proceedings of the National Academy of Sciences, 122: e2414985122.

Watson JEM, Evans MC, Carwardine J, Fuller RA, Joseph LN, Segan DB, Taylor MFJ, Fensham RJ, Possingham HP (2010) The capacity of Australia's protected-area system to represent threatened species. Conservation Biology,25: 324–332.

See Also

Other functions for adding targets: add_absolute_targets(), add_group_targets(), add_manual_targets(), add_relative_targets()

Examples

# set seed for reproducibility
set.seed(500)

# load data
sim_complex_pu_raster <- get_sim_complex_pu_raster()
sim_complex_features <- get_sim_complex_features()
sim_pu_polygons <- get_sim_pu_polygons()

# create base problem
p0 <-
  problem(sim_complex_pu_raster, sim_complex_features) %>%
  add_min_set_objective() %>%
  add_binary_decisions() %>%
  add_default_solver(gap = 0, verbose = FALSE)

# create problem with Polak et al. (2015) targets by using a function to
# specify the target setting method
p1 <-
  p0 %>%
  add_auto_targets(method = spec_polak_targets())

# create problem with Polak et al. (2015) targets by using a character value
# to specify the target setting method
# (note that this yields exactly the same targets as p1, and simply
# offers an alternative for specifying targets with default parameters)
p2 <-
  p0 %>%
  add_auto_targets(method = "polak")

# create problem with modified Polak et al. (2015) targets by using a
# a function to customize the parameters
# (note that this yields exactly the same targets as p1, and simply
# offers an alternative for specifying targets with default parameters)
p3 <-
  p0 %>%
  add_auto_targets(method = spec_polak_targets(common_relative_target = 0.3))

# solve problems
s <- c(solve(p1), solve(p2), solve(p3))
names(s) <- c(
  "default Polak targets", "default Polak targets", "modified Polak targets"
)

# plot solutions
plot(s, axes = FALSE)

# in the previous examples, we specified the targets for each of the features
# based on the same target setting method. we can also specify a particular
# target setting method for each feature. to achieve this, we can
# provide the targets as a list. for example, here we will specify that
# the first 10 features should have their targets based on
# Jung et al. (2021) targets (with default parameters), the following 15
# features should have their targets based on Ward et al. (2025) targets
# (with default parameters), and all the remaining features should have
# their targets based on modified Jung et al. (2021) targets.

# note that add_group_targets provides a more convenient interface for
# specifying targets for multiple features

# create a list with the target setting methods
target_list <- append(
  append(
    rep(list(spec_jung_targets()), 10),
    rep(list(spec_ward_targets()), 15)
  ),
   rep(
    list(spec_jung_targets(prop_uplift = 0.3)),
    terra::nlyr(sim_complex_features) - 25
  )
)

# create problem with the list of target setting methods
p4 <- p0 %>% add_auto_targets(method = target_list)

# solve problem
s4 <- solve(p4)

# plot solution
plot(s4, main = "solution", axes = FALSE)

# here we will show how to set the feature_units when the feature
# are not terra::rast() objects. in this example, we have planning units
# stored in an sf object (i.e., sim_pu_polygons) and the feature data will
# be stored as columns in the sf object.

# we will start by simulating feature data for the planning units.
# in particular, the simulated values will describe the amount of habitat
# for each feature expressed as acres (e.g., a value of 30 means 30 acres)
sim_pu_polygons$feature_1 <- runif(nrow(sim_pu_polygons), 0, 500)
sim_pu_polygons$feature_2 <- runif(nrow(sim_pu_polygons), 0, 600)
sim_pu_polygons$feature_3 <- runif(nrow(sim_pu_polygons), 0, 300)

# we will now build a problem with these data and specify the
# the feature units as acres because that those are the units we
# used for simulating the data. also, we will specify targets
# of 2 km^2 of habitat for each feature. although the feature units are
# acres, the function is able to able to automatically convert the units.
p5 <-
  problem(
    sim_pu_polygons,
    c("feature_1", "feature_2", "feature_3"),
    cost_column = "cost",
    feature_units = "acres"
  ) %>%
  add_min_set_objective() %>%
  add_auto_targets(
    method = spec_area_targets(targets = 2, area_units = "km^2")
  ) %>%
  add_binary_decisions() %>%
  add_default_solver(gap = 0, verbose = FALSE)

# solve problem
s5 <- solve(p5)

# plot solution
plot(s5[, "solution_1"], axes = FALSE)

Description

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 selecting planning units for protected area establishment.

Usage

add_binary_decisions(x)

Arguments

Details

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 during optimization. Also, if no decision is added to a problem(), then this decision will be used by default.

Value

An updated problem() object with the decisions added to it.

See Also

See decisions for an overview of all functions for adding decisions.

Other decisions: add_proportion_decisions(), add_semicontinuous_decisions()

Examples

# 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 multi-zone conservation problem
targs <- matrix(runif(15, 0.1, 0.2), nrow = 5, ncol = 3)

# build multi-zone 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)

Description

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).

Usage

## S4 method for signature 'ConservationProblem,ANY,ANY,ANY,ANY,data.frame'
add_boundary_penalties(x, penalty, edge_factor, formulation, zones, data)

## S4 method for signature 'ConservationProblem,ANY,ANY,ANY,ANY,matrix'
add_boundary_penalties(x, penalty, edge_factor, formulation, zones, data)

## S4 method for signature 'ConservationProblem,ANY,ANY,ANY,ANY,ANY'
add_boundary_penalties(x, penalty, edge_factor, formulation, zones, data)

Arguments

Details

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). Indeed, penalty 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.

Value

An updated problem() object with the penalties added to it.

Data format

The following formats can be used to specify data. Note that boundary data must always describe symmetric relationships between planning units.

data as a NULL value

Here boundary length data are automatically calculated using the boundary_matrix() function and then rescaled with rescale_matrix(). This is the default for data. Note that the boundary data must be supplied using one of the other formats below if x does not contain planning units that are spatially referenced. (e.g., planning unit data are a data.frame object or numeric vector).

data as a matrix/Matrix object

Here rows and columns correspond to different planning units and cell values represent the amount of boundary length shared between two planning units. Cells along the matrix diagonal denote the total boundary length associated with each planning unit. For example, boundary data in this format can be generated with the boundary_matrix() function.

data as a data.frame object

Here rows correspond to a pair of planning units and columns provide information about each pair of planning units. In particular, data must have 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 "id1" and "id2" columns contain the same values, then the value denotes the amount of exposed boundary length (not total boundary) for that particular planning unit. This format follows the the standard Marxan format for boundary data (i.e., per the "bound.dat" file).

Mathematical formulation

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 penalty, $E_z$ represent edge_factor, $B_{ij}$ represent data in matrix format (e.g., generated using boundary_matrix()), and $W_{zz}$ represent zones in matrix format.

$\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$. Additionally, to linearize the problem, the expression $X_{iz} \times X_{jy}$ is modeled using a set of continuous variables (bounded between 0 and 1) based on Beyer et al. (2016).

References

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 Also

See penalties for an overview of all functions for adding penalties. Also see add_neighbor_penalties() for a penalty that can reduce spatial fragmentation and has faster solver run times. Additionally, see calibrate_cohon_penalty() for assistance with selecting an appropriate penalty value.

Other functions for adding penalties: add_asym_connectivity_penalties(), add_connectivity_penalties(), add_cost_penalties(), add_feature_weights(), add_linear_penalties(), add_neighbor_penalties()

Examples

# 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 run-time 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)

Description

Specify that the CBC (COIN-OR branch and cut) software should be used to solve a conservation planning problem (Forrest & Lougee-Heimer 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).

Usage

add_cbc_solver(
  x,
  gap = 0.1,
  time_limit = .Machine$integer.max,
  presolve = 2,
  threads = 1,
  first_feasible = FALSE,
  start_solution = NULL,
  verbose = TRUE,
  control = list()
)

Arguments

Details

CBC is an open-source mixed integer programming solver that is part of the Computational Infrastructure for Operations Research (COIN-OR) project. This solver seems to have much better performance than the other open-source solvers (i.e., add_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.

Value

An updated problem() or multi_problem() object with the solver added to it.

Installation

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.

Start solution format

Broadly speaking, start_solution must be in the same format as the planning unit data in x. Further details on the correct format are described below.

x has numeric planning units

Here 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 x. Additionally, any planning units with missing cost (NA) values should also have missing (NA) values in the start_solution.

x has matrix planning units

Here 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 x. Additionally, any planning units with missing cost (NA) values for a particular zone should also have a missing (NA) values in start_solution.

x has terra::rast() planning units

Here start_solution be a terra::rast() object where different cells 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 x. Additionally, any planning units with missing cost (NA) values for a particular zone should also have missing (NA) values in start_solution.

x has data.frame planning units

Here 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 with missing cost (NA) values for a particular zone should also have missing (NA) values in start_solution.

x has sf::sf() planning units

Here 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, start_solution must also have the same coordinate reference system as the planning unit data. Furthermore, any planning units with missing cost (NA) values for a particular zone should also have missing (NA) values in start_solution.

References

Forrest J and Lougee-Heimer R (2005) CBC User Guide. In Emerging theory, Methods, and Applications (pp. 257–277). INFORMS, Catonsville, MD. doi:10.1287/educ.1053.0020.

See Also

Other functions for adding solvers: add_cplex_solver(), add_default_solver(), add_gurobi_solver(), add_highs_solver(), add_lsymphony_solver, add_rsymphony_solver()

Examples

# 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)

Description

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()).

Usage

## 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)

Arguments

Details

This function adds penalties to conservation planning problem to penalize solutions that have low connectivity. Specifically, it favors pair-wise connections between planning units that have high connectivity values (based on Önal and Briers 2002).

Value

An updated problem() object with the penalties added to it.

Data format

The following formats can be used to specify data.

data as a matrix/Matrix object

Here rows and columns correspond to different planning units and cell values denote the strength of connectivity between two planning units. Cells that occur along the matrix diagonal are treated as weights which indicate that planning units are more desirable in the solution. With this format, zones can be used to control the strength of connectivity between planning units in different zones. Note that the default for zones is to treat planning units allocated to different zones as having zero connectivity.

data as a data.frame object

Here rows correspond to a pair of planning units and columns provide information about each pair of planning units. In particular, data must have 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 strength of connectivity between them (following the Marxan format). If x has 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 particular zones. Note that if the "zone1" and "zone2" columns are present, then zones must be NULL.

data as an array object

Here a four-dimension array is used to specify connectivity data, 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.

Mathematical formulation

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 penalty, $D$ represent data, and $W$ represent zones.

If data is specified 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 data is specified 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$. Additionally, to linearize the problem, the expression $X_{iz} \times X_{jy}$ is modeled using a set of continuous variables (bounded between 0 and 1) based on Beyer et al. (2016).

Notes

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.

References

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.

Beyer HL, Dujardin Y, Watts ME, and Possingham HP (2016) Solving conservation planning problems with integer linear programming. Ecological Modelling, 228: 14–22.

Ö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 Also

See penalties for an overview of all functions for adding penalties. Also see add_asym_connectivity_penalties() to account for asymmetric connectivity between planning units. Additionally, see calibrate_cohon_penalty() for assistance with selecting an appropriate penalty value.

Other functions for adding penalties: add_asym_connectivity_penalties(), add_boundary_penalties(), add_cost_penalties(), add_feature_weights(), add_linear_penalties(), add_neighbor_penalties()

Examples

# 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 run-time
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 multi-zone 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 multi-zone 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 real-world 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 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)

Description

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.

Usage

## 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)

Arguments

Details

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).

Value

An updated problem() object with the constraints added to it.

Data format

The following formats can be used to specify data.

data as a NULL value

Here connection data are calculated automatically using the adjacency_matrix() function. This is the default for data. Note that the connection data must be manually specified using one of the other formats below when the planning unit data in x are not spatially referenced (e.g., in data.frame or numeric format).

data as a matrix/Matrix object

Here rows and columns correspond to different planning units and cell values indicate if two planning units are connected or not. In particular, 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 because each planning unit cannot be a connected with itself.

data as a data.frame object

Here rows correspond to a pair of planning units and columns provide information about each pair of planning units. In particular, data must have the columns: "id1", "id2", and "boundary". 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 specified (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).

Notes

In early versions, this function was named as the add_connected_constraints() function.

References

Önal H and Briers RA (2006) Optimal selection of a connected reserve network. Operations Research, 54: 379–388.

See Also

Other functions for adding constraints: add_cost_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()

Examples

# 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)

Description

Add constraints to a conservation planning problem to ensure that the cost of selected planning units meets certain criteria.

Usage

add_cost_constraints(x, budget, sense)

Arguments

Details

This function adds constraints constraints that can be used to ensure that the cost of solutions meet certain criteria (see Examples section below for details). For example, these constraints can be used to specify both minimum and maximum budget thresholds. Note that this function provided as a convenient alternative for adding linear constraints (per add_linear_constraints()) to a problem().

Value

An updated problem() object with the constraints added to it.

Mathematical formulation

The cost 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 costs associated with planning units $i \in I$ for zones $z \in Z$ (per data, if supplied as a matrix object), $\theta$ denote the constraint sense (per sense), and $B$ denote the budget threshold (per budget).

$\sum_{i}^{I} \sum_{z}^{Z} (D_{iz} \times X_{iz}) \space \theta \space t$

See Also

Other functions for adding 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(), add_neighbor_constraints()

Examples

# set seed for reproducibility
set.seed(600)

# load data
sim_complex_pu_raster <- get_sim_complex_pu_raster()
sim_complex_features <- get_sim_complex_features()

# create layer with 1s for all planning units
sim_ones_complex_raster <- (sim_complex_pu_raster * 0) + 1

# here we will formulate a multi-objective optimization problem
# that (i) minimizes the largest target shortfall for feature representation,
# (ii) minimizes the overall target shortfalls for feature representation,
# and (iii) minimizes the cost of the solution. since the
# first objective is to minimize the largest shortfall and the second
# objective is to minimize overall target shortfalls,
# this helps balance shortfalls among all features and better
# achieve complementarity. additionally, we will specify that
# (approximately) 30% of the study area should be selected (i.e., by
# specifying a budget for the upper threshold and a linear constraint for
# the lower threshold on the number of selected planning units).

# calculate budget based on 30% of the number of planning units
budget <-
  0.3 * terra::global(sim_ones_complex_raster, "sum", na.rm = TRUE)[[1]]

# build multi-objective conservation planning problem
mp <-
  multi_problem(
    obj1 =
      problem(sim_ones_complex_raster, sim_complex_features) %>%
      add_min_largest_shortfall_objective(budget = budget) %>%
      add_auto_targets("jung") %>%
      # note that this constraint only needs to be specified once
      add_cost_constraints(sense = ">=", budget = budget * 0.9) %>%
      add_binary_decisions(),
    obj2 =
      problem(sim_ones_complex_raster, sim_complex_features) %>%
      add_min_shortfall_objective(budget = budget) %>%
      # note that we use the same targets for both obj1 and obj2
      add_auto_targets("jung") %>%
      add_binary_decisions(),
    obj3 =
      problem(sim_complex_pu_raster, sim_complex_pu_raster) %>%
      add_min_penalties_objective() %>%
      # note a value of 1 is here because only the costs minimized
      add_cost_penalties(1) %>%
      add_binary_decisions()
  ) %>%
  add_default_solver(gap = 0.01, verbose = FALSE)

# to explore trade-offs between how well the feature targets
# are met and cost, we will generate a matrix of relative tolerance values
# for the hierarchical approach. note that the first column of this
# matrix will have only zeros to help promote balanced
# shortfalls across different features, and the second column
# will have non-zeros because we are interested in trade-offs between
# overall feature shortfalls and cost
rel_tol_matrix <- matrix(0, ncol = 2, nrow = 10)
rel_tol_matrix[, 2] <- seq(0, 0.5, length.out = nrow(rel_tol_matrix))

# display matrix
print(rel_tol_matrix)

# add hierarchical approach to multi-objective problem
mp <-
  mp %>%
  add_hier_approach(rel_tol = rel_tol_matrix)

# generate solutions and remove duplicates
ms <- solve(mp, remove_duplicates = TRUE)

# plot the solutions
plot(terra::rast(ms), axes = FALSE)

# extract objective values for the solutions
obj_matrix <- attributes(ms)$objective

# print the objective values
print(obj_matrix)

# plot the objectives values to visualize trade-offs
# (note that smaller values are better for both objectives)
plot(
  obj_matrix[, 2:3],
  main = "Trade-offs between objectives",
  xlab = "Species representation (overall shortfall)",
  ylab = "Solution cost"
)

Description

Add penalties to a conservation planning problem to penalize solutions that select planning units with higher cost values. These penalties assume a linear trade-off between the cost and the primary objective of the conservation planning problem (e.g., number of targets met for add_max_n_targets_met_objective().

Usage

add_cost_penalties(x, penalty)

Arguments

Details

This function penalizes solutions that have higher values according to the sum of the cost values associated with each planning unit, weighted by status of each planning unit in the solution. Note that this function provided as a convenient alternative for adding linear penalties (per add_linear_penalties()) to a problem().

Value

An updated problem() object with the penalties added to it.

Mathematical formulation

The cost 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$ (per penalty), and $D_{iz}$ represent the cost data for allocating planning unit $i \in I$ to zones $z \in Z$ (per data in matrix format).

$\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 Also

Other functions for adding penalties: add_asym_connectivity_penalties(), add_boundary_penalties(), add_connectivity_penalties(), add_feature_weights(), add_linear_penalties(), add_neighbor_penalties()

Examples

# set seed for reproducibility
set.seed(600)

# load data
sim_complex_pu_raster <- get_sim_complex_pu_raster()
sim_complex_features <- get_sim_complex_features()

# create layer with 1s for all planning units
sim_ones_complex_raster <- (sim_complex_pu_raster * 0) + 1

# here we will formulate a multi-objective optimization problem
# that (i) minimizes the largest target shortfall for feature representation,
# (ii) minimizes the overall target shortfalls for feature representation,
# and (iii) minimizes the cost of the solution. since the
# first objective is to minimize the largest shortfall and the second
# objective is to minimize overall target shortfalls,
# this helps balance shortfalls among all features and better
# achieve complementarity. additionally, we will specify that
# (approximately) 30% of the study area should be selected (i.e., by
# specifying a budget for the upper threshold and a linear constraint for
# the lower threshold on the number of selected planning units).

# calculate budget based on 30% of the number of planning units
budget <-
  0.3 * terra::global(sim_ones_complex_raster, "sum", na.rm = TRUE)[[1]]

# build multi-objective conservation planning problem
mp <-
  multi_problem(
    obj1 =
      problem(sim_ones_complex_raster, sim_complex_features) %>%
      add_min_largest_shortfall_objective(budget = budget) %>%
      add_auto_targets("jung") %>%
      # note that this constraint only needs to be specified once
      add_cost_constraints(sense = ">=", budget = budget * 0.9) %>%
      add_binary_decisions(),
    obj2 =
      problem(sim_ones_complex_raster, sim_complex_features) %>%
      add_min_shortfall_objective(budget = budget) %>%
      # note that we use the same targets for both obj1 and obj2
      add_auto_targets("jung") %>%
      add_binary_decisions(),
    obj3 =
      problem(sim_complex_pu_raster, sim_complex_pu_raster) %>%
      add_min_penalties_objective() %>%
      # note a value of 1 is here because only the costs minimized
      add_cost_penalties(1) %>%
      add_binary_decisions()
  ) %>%
  add_default_solver(gap = 0.01, verbose = FALSE)

# to explore trade-offs between how well the feature targets
# are met and cost, we will generate a matrix of relative tolerance values
# for the hierarchical approach. note that the first column of this
# matrix will have only zeros to help promote balanced
# shortfalls across different features, and the second column
# will have non-zeros because we are interested in trade-offs between
# overall feature shortfalls and cost
rel_tol_matrix <- matrix(0, ncol = 2, nrow = 10)
rel_tol_matrix[, 2] <- seq(0, 0.5, length.out = nrow(rel_tol_matrix))

# display matrix
print(rel_tol_matrix)

# add hierarchical approach to multi-objective problem
mp <-
  mp %>%
  add_hier_approach(rel_tol = rel_tol_matrix)

# generate solutions and remove duplicates
ms <- solve(mp, remove_duplicates = TRUE)

# plot the solutions
plot(terra::rast(ms), axes = FALSE)

# extract objective values for the solutions
obj_matrix <- attributes(ms)$objective

# print the objective values
print(obj_matrix)

# plot the objectives values to visualize trade-offs
# (note that smaller values are better for both objectives)
plot(
  obj_matrix[, 2:3],
  main = "Trade-offs between objectives",
  xlab = "Species representation (overall shortfall)",
  ylab = "Solution cost"
)

Description

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).

Usage

add_cplex_solver(
  x,
  gap = 0.1,
  time_limit = .Machine$integer.max,
  presolve = TRUE,
  threads = 1,
  verbose = TRUE
)

Arguments

Details

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/ilog-cplex-optimization-studio/cplex-optimizer).

Value

An updated problem() or multi_problem() object with the solver added to it.

Installation

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/ilog-cplex-optimization-studio/cplex-optimizer). 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/x86-64_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.

References

IBM (2017) IBM ILOG CPLEX Optimization Studio CPLEX User's Manual. Version 12 Release 8. IBM ILOG CPLEX Division, Incline Village, NV.

See Also

Other functions for adding solvers: add_cbc_solver(), add_default_solver(), add_gurobi_solver(), add_highs_solver(), add_lsymphony_solver, add_rsymphony_solver()

Examples

# 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)

Description

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.

Usage

add_cuts_portfolio(x, number_solutions = 10, verbose = TRUE)

Arguments

Details

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.

Value

An updated problem() object with the portfolio added to it.

Notes

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.

References

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 Also

See portfolios for an overview of all functions for adding a portfolio.

Other functions for adding portfolios: add_default_portfolio(), add_extra_portfolio(), add_gap_portfolio(), add_shuffle_portfolio(), add_single_portfolio(), add_top_portfolio()

Examples

# 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 multi-layer raster object
s1 <- terra::rast(s1)

# plot solutions in portfolio
plot(s1, axes = FALSE)

# build multi-zone 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)

Description

Generate a portfolio based on defaults.

Usage

add_default_portfolio(x)

Arguments

Details

By default, this is portfolio is added to problem() objects if no other portfolios is manually specified. In particular, this function adds the add_single_portfolio() function to x so that only a single solution is generated.

Value

An updated problem() object with the portfolio added to it.

See Also

Other functions for adding portfolios: add_cuts_portfolio(), add_extra_portfolio(), add_gap_portfolio(), add_shuffle_portfolio(), add_single_portfolio(), add_top_portfolio()

Examples

# set seed for reproducibility
set.seed(600)

# load data
sim_pu_raster <- get_sim_pu_raster()
sim_features <- get_sim_features()

# create minimal problem with default portfolio
p <-
  problem(sim_pu_raster, sim_features) %>%
  add_min_set_objective() %>%
  add_relative_targets(0.05) %>%
  add_default_portfolio() %>%
  add_default_solver(gap = 0, verbose = FALSE)

# solve problem
s <- solve(p)

# plot solution
plot(s)

Description

Specify that the best solver currently available should be used to solve a conservation planning problem.

Usage

add_default_solver(x, ...)

Arguments

Details

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).

Value

An updated problem() or multi_problem() object with the solver added to it.

References

Schuster R, Hanson JO, Strimas-Mackey M, and Bennett JR (2020). Exact integer linear programming solvers outperform simulated annealing for solving conservation planning problems. PeerJ, 8: e9258.

See Also

See solvers for an overview of all functions for adding a solver.

Other functions for adding solvers: add_cbc_solver(), add_cplex_solver(), add_gurobi_solver(), add_highs_solver(), add_lsymphony_solver, add_rsymphony_solver()

Examples

# set seed for reproducibility
set.seed(600)

# load data
sim_pu_raster <- get_sim_pu_raster()
sim_features <- get_sim_features()

# create minimal problem with default portfolio
p <-
  problem(sim_pu_raster, sim_features) %>%
  add_min_set_objective() %>%
  add_relative_targets(0.05) %>%
  add_binary_decisions() %>%
  add_default_solver(gap = 0, verbose = FALSE)

# solve problem
s <- solve(p)

# plot solution
plot(s)

Description

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.

Usage

add_extra_portfolio(x)

Arguments

Details

This strategy for generating a portfolio requires problems to be solved using the Gurobi software (i.e., using add_gurobi_solver(). Specifically, version 8.0.0 (or greater) of the gurobi package must be installed.

Value

An updated problem() object with the portfolio added to it.

See Also

Other functions for adding portfolios: add_cuts_portfolio(), add_default_portfolio(), add_gap_portfolio(), add_shuffle_portfolio(), add_single_portfolio(), add_top_portfolio()

Examples

# 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 multi-layer raster object
s1 <- terra::rast(s1)

# print number of solutions found
print(terra::nlyr(s1))

# plot solutions
plot(s1, axes = FALSE)

# create multi-zone 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)

Description

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.

Usage

## 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)

Arguments

Details

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.

Value

An updated problem() object with the constraints added to it.

Data format

The following formats can be used to specify data.

data as a NULL value

Here connection data are calculated automatically using the adjacency_matrix() function. This is the default 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 x is not spatially referenced (e.g., in data.frame or numeric format).

data as amatrix/Matrix object

Here rows and columns correspond to different planning units and cell values indicates if a pair of planning units should be treated as connected or not. In particular, 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 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 object

Here rows correspond to a pair of planning units and columns provide information about each pair of planning units. In particular, data must have 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 binary numeric values that indicate if the two planning units specified in the "id1" and "id2" columns should be treated as connected or not. Note that pairs of connected planning units are treated as being potentially dispersible for all features.

data as a list object

Here a matrix, Matrix, or data.frame object is specified for each feature to indicate which planning units should be treated as connected for that feature. In particular, 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 formats described previously.

Notes

In early versions, it was named as the add_corridor_constraints function.

References

Ö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 reserve-network 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 Also

Other functions for adding constraints: add_contiguity_constraints(), add_cost_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()

Examples

# 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)

Description

Add features weights to a conservation planning problem. Specifically, some objective functions aim to maximize (or minimize) a metric that evaluates how well a set of features are represented by a solution (e.g., maximize the number of feature targets that are met, add_max_n_targets_met_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.

Usage

## S4 method for signature 'ConservationProblem,numeric'
add_feature_weights(x, weights)

## S4 method for signature 'ConservationProblem,matrix'
add_feature_weights(x, weights)

Arguments

Details

Weights are only considered during optimization if a budget-limited objective is specified (e.g., add_max_cover_objective(), add_min_shortfall_objective()). Although weights can be added to problems that have the minimum set objective (i.e., add_min_set_objective()), they will have no effect during optimization and a warning will be thrown.

Value

An updated problem() with the weights added to it.

Weights format

The following formats can be used to specify weights. Note that weights must have values that are greater than, or equal to, zero (in other words, non-negative values).

weights as a numeric vector

Here weight values are specified for each feature. Note that this format cannot be used if x has multiple zones.

weights as a matrix object

Here weight values are specified for each feature in each zone. In particular, each row corresponds to a different feature in x, each column corresponds to a different zone in x, and cell values specify the weight value for representing a particular feature in a particular zone. Note that if the problem contains targets created using add_manual_targets(), then weights should be a matrix that has a single column containing the weight value for meeting each target.

See Also

See penalties for an overview of all functions for adding penalties.

Other functions for adding penalties: add_asym_connectivity_penalties(), add_boundary_penalties(), add_connectivity_penalties(), add_cost_penalties(), add_linear_penalties(), add_neighbor_penalties()

Examples

# 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_n_targets_met_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 pre-calculated 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., eco-tourism). 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
# trade-off 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 multi-zone problem with maximum number of targets met 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_n_targets_met_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_n_targets_met_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)

Description

Generate a portfolio of solutions for a conservation planning problem by finding a certain number of solutions that are all within a pre-specified optimality gap. This method is useful for generating multiple solutions that can be used to calculate selection frequencies for small-sized problems (similar to Marxan).

Usage

add_gap_portfolio(x, number_solutions = 10, pool_gap = 0.1)

Arguments

Details

This strategy for generating a portfolio requires problems to be solved using the Gurobi software (i.e., using add_gurobi_solver(). Specifically, version 9.0.0 (or greater) of the gurobi package must be installed. Note if the total number of solutions that meet the optimality gap is fewer than number_solutions, then the number of solutions returned may be fewer than requested. Also, note that this portfolio function only works with problems that have binary decisions (i.e., specified using add_binary_decisions()).

Value

An updated problem() object with the portfolio added to it.

See Also

Other functions for adding portfolios: add_cuts_portfolio(), add_default_portfolio(), add_extra_portfolio(), add_shuffle_portfolio(), add_single_portfolio(), add_top_portfolio()

Examples

# 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 multi-layer raster
s1 <- terra::rast(s1)

# print number of solutions found
print(terra::nlyr(s1))

# plot solutions
plot(s1, axes = FALSE)

# create multi-zone  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 multi-layer 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)

Description

Add targets to a conservation planning problem, wherein each feature is assigned to a particular group and a target setting method is specified for each feature group. This function is designed to provide a convenient alternative to add_auto_targets().

Usage

add_group_targets(x, groups, method)

Arguments

Value

An updated problem() object with the targets added to it.

Target methods

Here method is used to specify a target setting method for each group. In particular, each element should correspond to a different group, with the name of the element indicating the group and the value denoting the target setting method. One option for specifying a particular target setting method is to use a character value that denotes the name of the method (e.g., "jung" to set targets following Jung et al. 2021). This option is particularly useful if default parameters should be considered during target calculations. Alternatively, another option for specifying a particular target setting method is to use a function to define an object (e.g., spec_jung_targets() to set targets following Jung et al. 2021). This alternative option is particularly useful if customized parameters should be considered during target calculations. Note that a method can contain a mix of target setting methods defined using character values and functions.

The following character values can be used to specify target setting methods: "jung" (per Jung et al. 2021), "polak" (per Polak et al. 2016), "rodrigues" (per Rodrigues et al. 2004), "ward" (per Ward et al. 2025), and "watson" (per Watson et al. 2010). (2010). Additionally, the following values can be used to set targets based on criteria from the IUCN Red List of Threatened Species (IUCN 2025): "rl_species_VU_A1_B1" "rl_species_EN_A1_B1", "rl_species_CR_A1_B1", "rl_species_VU_A1_B2" "rl_species_EN_A1_B2", "rl_species_CR_A1_B2", "rl_species_VU_A2_B1" "rl_species_EN_A2_B1", "rl_species_CR_A2_B1", "rl_species_VU_A2_B2" "rl_species_EN_A2_B2", "rl_species_CR_A2_B2", "rl_species_VU_A3_B1" "rl_species_EN_A3_B1", "rl_species_CR_A3_B1", "rl_species_VU_A3_B2" "rl_species_EN_A3_B2", "rl_species_CR_A3_B2", "rl_species_VU_A4_B1" "rl_species_EN_A4_B1", "rl_species_CR_A4_B1", "rl_species_VU_A4_B2" "rl_species_EN_A4_B2", and "rl_species_CR_A4_B2". Furthermore, the following values can be used to set targets based on criteria from the IUCN Red List of Ecosystems (IUCN 2024): "rl_ecosystem_VU_A1_B1", "rl_ecosystem_EN_A1_B1", "rl_ecosystem_CR_A1_B1", "rl_ecosystem_VU_A1_B2", "rl_ecosystem_EN_A1_B2", "rl_ecosystem_CR_A1_B2", "rl_ecosystem_VU_a2a_B1", "rl_ecosystem_EN_a2a_B1", "rl_ecosystem_CR_a2a_B1", "rl_ecosystem_VU_a2a_B2", "rl_ecosystem_EN_a2a_B2", "rl_ecosystem_CR_a2a_B2", "rl_ecosystem_VU_a2b_B1", "rl_ecosystem_EN_a2b_B1", "rl_ecosystem_CR_a2b_B1", "rl_ecosystem_VU_a2b_B2", "rl_ecosystem_EN_a2b_B2", "rl_ecosystem_CR_a2b_B2", "rl_ecosystem_VU_A3_B1", "rl_ecosystem_EN_A3_B1", "rl_ecosystem_CR_A3_B1", "rl_ecosystem_VU_A3_B2", "rl_ecosystem_EN_A3_B2", and "rl_ecosystem_CR_A3_B2".' For convenience, these options can also be specified with lower case letters (e.g., "rl_species_VU_A1_B1" may alternatively be specified as "rl_species_vu_a1_b1"). Note that target setting methods that require additional data or parameters cannot be specified with character values.

The following functions can be used to create an object for specifying a target setting method: spec_absolute_targets(), spec_relative_targets(), spec_area_targets(), spec_interp_absolute_targets(), spec_interp_area_targets(), spec_duran_targets(), spec_jung_targets(), spec_polak_targets(), spec_pop_size_targets(), spec_rl_ecosystem_targets(), spec_rl_species_targets(), spec_rodrigues_targets(), spec_rule_targets(), spec_ward_targets(), spec_watson_targets(), spec_wilson_targets(), spec_min_targets(), and spec_max_targets().

Target setting

Many conservation planning problems require targets. Targets are used to specify the minimum amount, or proportion, of a feature's spatial distribution that should ideally be protected. This is important so that the optimization process can weigh the merits and trade-offs between improving the representation of one feature over another feature. Although it can be challenging to set meaningful targets, this is a critical step for ensuring that prioritizations meet the stakeholder objectives that underpin a prioritization exercise (Carwardine et al. 2009). In other words, targets play an important role in ensuring that a priority setting process is properly tuned according to stakeholder requirements. For example, targets provide a mechanism for ensuring that a prioritization secures enough habitat to promote the long-term persistence of each threatened species, culturally important species, or economically important ecosystem services under consideration. Since there is often uncertainty regarding stakeholder objectives (e.g., how much habitat should be protected for a given species) or the influence of particular target on a prioritization (e.g., how would setting a 90% or 100% for a threatened species alter priorities), it is often useful to generate and compare a suite of prioritizations based on different target scenarios.

Data calculations

Many of the functions for specifying target setting methods involve calculating targets based on the spatial extent of the features in x (e.g., spec_jung_targets(), [spec_rodrigues_targets(), and others). Although this function for adding targets can be readily applied to problem() objects that have the feature data provided as a terra::rast() object, you will need to specify the spatial units for the features when building a problem() object if the feature data are provided in a different format. In particular, if the feature data are provided as a data.frame or character vector, then you will need to specify feature_units when using the problem() function. See the Examples section below for a demonstration of using the feature_units parameter.

References

Carwardine J, Klein CJ, Wilson KA, Pressey RL, Possingham HP (2009) Hitting the target and missing the point: target‐based conservation planning in context. Conservation Letters, 2: 4–11.

Jung M, Arnell A, de Lamo X, García-Rangel S, Lewis M, Mark J, Merow C, Miles L, Ondo I, Pironon S, Ravilious C, Rivers M, Schepaschenko D, Tallowin O, van Soesbergen A, Govaerts R, Boyle BL, Enquist BJ, Feng X, Gallagher R, Maitner B, Meiri S, Mulligan M, Ofer G, Roll U, Hanson JO, Jetz W, Di Marco M, McGowan J, Rinnan DS, Sachs JD, Lesiv M, Adams VM, Andrew SC, Burger JR, Hannah L, Marquet PA, McCarthy JK, Morueta-Holme N, Newman EA, Park DS, Roehrdanz PR, Svenning J-C, Violle C, Wieringa JJ, Wynne G, Fritz S, Strassburg BBN, Obersteiner M, Kapos V, Burgess N, Schmidt- Traub G, Visconti P (2021) Areas of global importance for conserving terrestrial biodiversity, carbon and water. Nature Ecology and Evolution, 5:1499–1509.

Polak T, Watson JEM, Fuller RA, Joseph LN, Martin TG, Possingham HP, Venter O, Carwardine J (2015) Efficient expansion of global protected areas requires simultaneous planning for species and ecosystems. Royal Society Open Science, 2: 150107.

See Also

Other functions for adding targets: add_absolute_targets(), add_auto_targets(), add_manual_targets(), add_relative_targets()

Examples

# set seed for reproducibility
set.seed(500)

# load data
sim_complex_pu_raster <- get_sim_complex_pu_raster()
sim_complex_features <- get_sim_complex_features()
sim_pu_polygons <- get_sim_pu_polygons()

# create grouping data, wherein each feature will be randomly
# assigned to the group A, B, C, or D
sim_groups <- sample(
  c("A", "B", "C", "D"), terra::nlyr(sim_complex_features), replace = TRUE
)

# create problem where each feature in group A is assigned a 20% target,
# each feature in group B is assigned a target based on Jung et al. (2021),
# each feature in group C is assigned a target based on Ward et al. (2025),
# and each feature in group D is assigned a target based on Rodrigues
# et al. (2004)
p1 <-
  problem(sim_complex_pu_raster, sim_complex_features) %>%
  add_min_set_objective() %>%
  add_group_targets(
    groups = sim_groups,
    method = list(
      A = spec_relative_targets(0.2),
      B = spec_jung_targets(),
      C = spec_ward_targets(),
      D = spec_rodrigues_targets()
    )
  ) %>%
  add_binary_decisions() %>%
  add_default_solver(verbose = FALSE)

# solve problem
s1 <- solve(p1)

# plot solution
plot(s1, main = "solution", axes = FALSE)

# here we will show how to set the feature_units when the feature
# are not terra::rast() objects. in this example, we have planning units
# stored in an sf object (i.e., sim_pu_polygons) and the feature data will
# be stored as columns in the sf object.

# we will start by simulating feature data for the planning units.
# in particular, the simulated values will describe the amount of habitat
# for each feature expressed as acres (e.g., a value of 30 means 30 acres).
sim_pu_polygons$feature_1 <- runif(nrow(sim_pu_polygons), 0, 500)
sim_pu_polygons$feature_2 <- runif(nrow(sim_pu_polygons), 0, 600)
sim_pu_polygons$feature_3 <- runif(nrow(sim_pu_polygons), 0, 300)

# we will now build a problem with these data and specify the
# the feature units as acres because that those are the units we
# used for simulating the data. also, we will specify targets
# of 2 km^2 of habitat for the first feature, and 3 km^2 for the
# remaining two features by assigning the features to groups.
# although the feature units are acres, the function is able to able
# to automatically convert the units.
p2 <-
  problem(
    sim_pu_polygons,
    c("feature_1", "feature_2", "feature_3"),
    cost_column = "cost",
    feature_units = "acres"
  ) %>%
  add_min_set_objective() %>%
  add_group_targets(
    groups = c("A", "B", "B"),
    method = list(
      A = spec_area_targets(targets = 2, area_units = "km^2"),
      B = spec_area_targets(targets = 3, area_units = "km^2")
    )
  ) %>%
  add_binary_decisions() %>%
  add_default_solver(gap = 0, verbose = FALSE)

# solve problem
s2 <- solve(p2)

# plot solution
plot(s2[, "solution_1"], axes = FALSE)