Simulated examples
Data
To investigate the behavior of the problem, we will generate prioritizations for three simulated species. We will use the unreliable formulation of the problem to understand the basics, and later move onto the reliable formulation. The first species (termed ‘uniform’) will represent a hyper-generalist. This species will inhabit all areas with equal probability. The second species (termed ‘normal’) will represent a species with a single range core. The third species (termed ‘bimodal’) will represent a species with two distinct ecotypes, each with their own range core. To reduce computational time for this example, we will use a 10 \(\times\) 10 grid of square planning units.
# make planning units
sim_pus <- sim.pus(100L)
# simulate species distributions
sim_spp <- lapply(c("uniform", "normal", "bimodal"), sim.species, n = 1,
x = sim_pus, res = 1)Let’s see what these species’ distributions look like.
# plot species
plot(stack(sim_spp),
main = c("Uniform species", "Normal species", "Bimodal species"),
addfun = function() lines(sim_pus), nc = 3)
Distribution of three simulated species. Each square represents a planning unit. The color of each square denotes the probability that individuals from each species occupy it.
Next, we will generate a set of demand points. To understand the effects of probabilities and weights on the demand points, we will generate the demand points in geographic space. These demand points will be the centroids of the planning units. Additionally, we will use the same set of demand points for each species and only vary the weights of the demand points between species. Note that we are only using the same distribution of demand points for different species for teaching purposes. It is strongly recommended to use different demand points for different species in real-world conservation planning exercises. See the case-study section of this tutorial for examples on how to generate suitable demand points.
# generate coordinates for pus/demand points
pu_coords <- gCentroid(sim_pus, byid = TRUE)
# calculate weights
sim_dps <- lapply(sim_spp, function(x) extract(x, pu_coords))
# create demand point objects
sim_dps <- lapply(sim_dps, function(x) DemandPoints(pu_coords@coords,
weights = c(x)))Now, we will construct a RapUnsolved object to store our input data and parameters. This contains all the information to generate prioritizations.
## create RapUnreliableOpts object
# this stores parameters for the unreliable formulation problem (ie. BLM)
sim_ro <- RapUnreliableOpts()
## create RapData object
# create data.frame with species info
species <- data.frame(name = c("uniform", "normal", "bimodal"))
## create data.frame with species and space targets
# amount targets at 20% (denoted with target=0)
# space targets at 20% (denoted with target=1)
targets <- expand.grid(species = 1:3, target = 0:1, proportion = 0.2)
# calculate probability of each species in each pu
pu_probabilities <- calcSpeciesAverageInPus(sim_pus, stack(sim_spp))
## create AttributeSpace object
# this stores the coordinates of the planning units in an attribute space
# and the coordinates and weights of demand points in the space
pu_points <- PlanningUnitPoints(coords = pu_coords@coords,
ids = seq_len(nrow(sim_pus@data)))
attr_spaces <- AttributeSpaces(
list(AttributeSpace(planning.unit.points = pu_points,
demand.points = sim_dps[[1]],
species = 1L),
AttributeSpace(planning.unit.points = pu_points,
demand.points = sim_dps[[2]],
species = 2L),
AttributeSpace(planning.unit.points = pu_points,
demand.points = sim_dps[[3]],
species = 3L)),
name = "geographic")
# generate boundary data information
boundary <- calcBoundaryData(sim_pus)
## create RapData object
# this stores all the input data for the prioritization
sim_rd <- RapData(sim_pus@data, species, targets, pu_probabilities,
list(attr_spaces), boundary, SpatialPolygons2PolySet(sim_pus))
## create RapUnsolved object
# this stores all the input data and parameters needed to generate
# prioritizations
sim_ru <- RapUnsolved(sim_ro, sim_rd)Single-species prioritizations
Amount-based targets
To investigate the effects of space-based targets, we will generate a prioritization for each species using only amount-based targets and compare them to prioritizations generated using amount- and space-based targets. To start off, we will generate a prioritization for the uniform species using amount-based targets. To do this, we will generate a new sim_ru object by extracting out the data for the uniform species from the sim_ru object. Then, we will update the targets in the new object. Finally, we will solve the object to generate a prioritization that fulfills the targets for minimal cost.
# create new object with just the uniform species
sim_ru_s1 <- spp.subset(sim_ru, "uniform")
# update amount targets to 20% and space targets to 0%
sim_ru_s1 <- update(sim_ru_s1, amount.target = 0.2, space.target = NA,
solve = FALSE)options(error = function() {
traceback()
q("no", 1, FALSE)
})
# solve problem to identify prioritization
sim_rs_s1_amount <- solve(sim_ru_s1, Threads = threads)## Optimize a model with 1 rows, 100 columns and 100 nonzeros
## Variable types: 0 continuous, 100 integer (100 binary)
## Coefficient statistics:
## Matrix range [5e-01, 5e-01]
## Objective range [1e+00, 1e+00]
## Bounds range [1e+00, 1e+00]
## RHS range [1e+01, 1e+01]
## Found heuristic solution: objective 20.0000000
## Presolve removed 1 rows and 100 columns
## Presolve time: 0.00s
## Presolve: All rows and columns removed
##
## Explored 0 nodes (0 simplex iterations) in 0.00 seconds
## Thread count was 1 (of 4 available processors)
##
## Solution count 1: 20
##
## Optimal solution found (tolerance 1.00e-01)
## Best objective 2.000000000000e+01, best bound 2.000000000000e+01, gap 0.0000%## Warning in validityMethod(object): some species have space.held values less
## than 0, and thus are poorly represented## show summary
# note the format for this is similar to that used by Marxan
# see ?raptr::summary for details on this table
summary(sim_rs_s1_amount)## Run_Number Status Score Cost Planning_Units Connectivity_Total
## 1 1 OPTIMAL 20 20 20 220
## Connectivity_In Connectivity_Edge Connectivity_Out
## 1 42 168 10
## Connectivity_In_Fraction
## 1 0.1909091# show amount held
amount.held(sim_rs_s1_amount)## uniform
## 1 0.2# show space held
space.held(sim_rs_s1_amount)## uniform (Space 1)
## 1 -0.2363636Now that we have generated a prioritization, let’s see what it looks like. We can use the spp.plot method to see how the prioritization overlaps with the uniform species’ distribution. Note that since all planning units have equal probabilities for this species, all planning units have the same fill color.
# plot the prioritization and the uniform species' distribution
spp.plot(sim_rs_s1_amount, 1, main = "Uniform species")
A prioritization for the uniformly distributed species generated using amount-based targets (20%). Sqaures represent planning units. Planning units with a green border are selected for prioritization, and their colour denotes the probability they are inhabited by the species.
The prioritization for the uniform species appears to be just a random selection of planning units. This behavior is due to the fact that any prioritization with 20 planning units is optimal. By relying on just amount targets, this solution may preserve a section of the species’ range core, or just focus on the range margin, or some random part of its range–no emphasis is directed towards preserving different parts of the species’ range. This behavior highlights a fundamental limitation of just using amount-based targets. In the absence of additional criteria, conventional reserve selection problems do not contain any additional information to identify the most effective prioritization.
Now, we will generate a prioritization for the normally distributed species using amount-based targets. We will use a similar process to what we used for the uniformly distributed species, but for brevity, we will use code to generate solutions immediately after updating the object.
# create new object with just the normal species
sim_ru_s2 <- spp.subset(sim_ru, "normal")# update amount targets to 20% and space targets to 0% and solve it
sim_rs_s2_amount <- update(sim_ru_s2, amount.target = 0.2, space.target = NA,
solve = TRUE, Threads = threads)## Optimize a model with 1 rows, 100 columns and 100 nonzeros
## Variable types: 0 continuous, 100 integer (100 binary)
## Coefficient statistics:
## Matrix range [7e-02, 7e-01]
## Objective range [1e+00, 1e+00]
## Bounds range [1e+00, 1e+00]
## RHS range [7e+00, 7e+00]
## Found heuristic solution: objective 27.0000000
## Presolve removed 0 rows and 86 columns
## Presolve time: 0.00s
## Presolved: 1 rows, 14 columns, 14 nonzeros
## Variable types: 0 continuous, 14 integer (0 binary)
## Presolved: 1 rows, 14 columns, 14 nonzeros
##
##
## Root relaxation: objective 9.864476e+00, 6 iterations, 0.00 seconds
##
## Nodes | Current Node | Objective Bounds | Work
## Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
##
## 0 0 9.86448 0 1 27.00000 9.86448 63.5% - 0s
## H 0 0 10.0000000 9.86448 1.36% - 0s
##
## Explored 1 nodes (6 simplex iterations) in 0.00 seconds
## Thread count was 2 (of 4 available processors)
##
## Solution count 2: 10 27
##
## Optimal solution found (tolerance 1.00e-01)
## Best objective 1.000000000000e+01, best bound 1.000000000000e+01, gap 0.0000%# show summary
summary(sim_rs_s2_amount)## Run_Number Status Score Cost Planning_Units Connectivity_Total
## 1 1 OPTIMAL 10 10 10 220
## Connectivity_In Connectivity_Edge Connectivity_Out
## 1 12 192 16
## Connectivity_In_Fraction
## 1 0.05454545# show amount held
amount.held(sim_rs_s2_amount)## normal
## 1 0.2026153# show space held
space.held(sim_rs_s2_amount)## normal (Space 1)
## 1 0.5909494Now let’s visualize the prioritization we made for the normal species.
# plot the prioritization and the normal species' distribution
spp.plot(sim_rs_s2_amount, 1, main = "Normal species")
A prioritization for the normally distributed species generated using amount-based targets (20%). See Figure 3 caption for conventions.
The amount-based prioritization for the normal species focuses only on the species’ range core. This prioritization fails to secure any peripheral parts of the species’ distribution. As a consequence, it may miss out on populations with novel adaptations to environmental conditions along the species’ range margin.
Now, let’s generate an amount-based target for the bimodally distributed species view it.
# create new object with just the bimodal species
sim_ru_s3 <- spp.subset(sim_ru, "bimodal")# update amount targets to 20% and space targets to 0% and solve it
sim_rs_s3_amount <- update(sim_ru_s3, amount.target = 0.2,
space.target = NA, Threads = threads)## Optimize a model with 1 rows, 100 columns and 100 nonzeros
## Variable types: 0 continuous, 100 integer (100 binary)
## Coefficient statistics:
## Matrix range [7e-03, 9e-01]
## Objective range [1e+00, 1e+00]
## Bounds range [1e+00, 1e+00]
## RHS range [7e+00, 7e+00]
## Found heuristic solution: objective 21.0000000
## Presolve removed 0 rows and 75 columns
## Presolve time: 0.00s
## Presolved: 1 rows, 25 columns, 25 nonzeros
## Variable types: 0 continuous, 25 integer (0 binary)
## Presolved: 1 rows, 25 columns, 25 nonzeros
##
##
## Root relaxation: objective 7.919039e+00, 17 iterations, 0.00 seconds
##
## Nodes | Current Node | Objective Bounds | Work
## Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
##
## 0 0 7.91904 0 1 21.00000 7.91904 62.3% - 0s
## H 0 0 8.0000000 7.91904 1.01% - 0s
##
## Explored 1 nodes (17 simplex iterations) in 0.00 seconds
## Thread count was 2 (of 4 available processors)
##
## Solution count 2: 8 21
##
## Optimal solution found (tolerance 1.00e-01)
## Best objective 8.000000000000e+00, best bound 8.000000000000e+00, gap 0.0000%# plot the prioritization and the bimodal species' distribution
spp.plot(sim_rs_s3_amount, 1, main = "Bimodal species")
A prioritization for the bimodally distributed species generated using amount-based targets (20%). See Figure 3 caption for conventions.
# show summary
summary(sim_rs_s3_amount)## Run_Number Status Score Cost Planning_Units Connectivity_Total
## 1 1 OPTIMAL 8 8 8 220
## Connectivity_In Connectivity_Edge Connectivity_Out
## 1 9 197 14
## Connectivity_In_Fraction
## 1 0.04090909# show amount held
amount.held(sim_rs_s3_amount)## bimodal
## 1 0.2018391# show space held
space.held(sim_rs_s3_amount)## bimodal (Space 1)
## 1 0.1200278The amount-based prioritization for the bimodally distributed species only selects planning units in the bottom left corner of the study area. This prioritization only preserves individuals belonging to one of the two ecotypes. As a consequence, this prioritization may fail to preserve a representative sample of the genetic variation found inside this species.
Amount-based and space-based targets
Now that we have generated a prioritization for each species using only amount-based targets, we will generate a prioritizations using both amount-based and space-targets. To do this we will update the space targets in our amount-based prioritizations to 85%, and store the new prioritizations in new objects.
First, let’s do this for the uniform species.
# make new prioritization
sim_rs_s1_space <- update(sim_rs_s1_amount, amount.target = 0.2,
space.target = 0.85, Threads = threads)## Optimize a model with 10102 rows, 10100 columns and 40000 nonzeros
## Variable types: 0 continuous, 10100 integer (10100 binary)
## Coefficient statistics:
## Matrix range [5e-01, 8e+01]
## Objective range [1e+00, 1e+00]
## Bounds range [1e+00, 1e+00]
## RHS range [1e+00, 1e+02]
## Found heuristic solution: objective 95.0000000
## Presolve removed 36 rows and 0 columns
## Presolve time: 1.82s
## Presolved: 10066 rows, 10100 columns, 40104 nonzeros
## Variable types: 0 continuous, 10100 integer (10100 binary)
## Presolved: 10066 rows, 10100 columns, 40104 nonzeros
##
##
## Root relaxation: objective 2.000000e+01, 621 iterations, 0.06 seconds
##
## Nodes | Current Node | Objective Bounds | Work
## Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
##
## * 0 0 0 20.0000000 20.00000 0.00% - 1s
##
## Explored 0 nodes (1154 simplex iterations) in 1.96 seconds
## Thread count was 2 (of 4 available processors)
##
## Solution count 2: 20 95
##
## Optimal solution found (tolerance 1.00e-01)
## Best objective 2.000000000000e+01, best bound 2.000000000000e+01, gap 0.0000%# show summary
summary(sim_rs_s1_space)## Run_Number Status Score Cost Planning_Units Connectivity_Total
## 1 1 OPTIMAL 20 20 20 220
## Connectivity_In Connectivity_Edge Connectivity_Out
## 1 13 147 60
## Connectivity_In_Fraction
## 1 0.05909091# show amount held
amount.held(sim_rs_s1_space)## uniform
## 1 0.2# show space held
space.held(sim_rs_s1_space)## uniform (Space 1)
## 1 0.909697Let’s take a look at the prioritization for the uniform species with amount-based and space-based targets. Then, let’s compare the solutions for the amount-based prioritization with the new prioritization using both amount and space targets.
# plot the prioritization and the uniform species' distribution
spp.plot(sim_rs_s1_space, "uniform", main = "Uniform species")
A prioritization for the uniformly distributed species generated using amount-based targets (20%) and space-based targets (85%). See Figure 3 caption for conventions.
# plot the difference between old and new prioritizations
plot(sim_rs_s1_amount, sim_rs_s1_space, 1, 1,
main = "Difference between solutions")
Difference between two prioritizations for the uniformly distributed species. Prioritisation \(X\) was generated using just amount-based targets (20%), and prioritization \(Y\) was generated using an additional space-based target (85%).
Here we can see that by including a space-target, the prioritization is spread out evenly across the species’ distribution. Unlike the amount-based prioritization, this prioritization samples all the different parts of the species’ distribution.
Now, let’s generate a prioritization for the normally distributed species that considers amount-based and space-based targets. Then, let’s visualize the new prioritization and compare it to the old amount-based prioritization.
# make new prioritization
sim_rs_s2_space <- update(sim_rs_s2_amount, amount.target = 0.2,
space.target = 0.85, Threads = threads)## Optimize a model with 10102 rows, 10100 columns and 40000 nonzeros
## Variable types: 0 continuous, 10100 integer (10100 binary)
## Coefficient statistics:
## Matrix range [7e-02, 4e+01]
## Objective range [1e+00, 1e+00]
## Bounds range [1e+00, 1e+00]
## RHS range [1e+00, 6e+01]
## Found heuristic solution: objective 91.0000000
## Presolve removed 220 rows and 0 columns
## Presolve time: 1.28s
## Presolved: 9882 rows, 10100 columns, 41664 nonzeros
## Variable types: 0 continuous, 10100 integer (10100 binary)
## Presolved: 9882 rows, 10100 columns, 41664 nonzeros
##
##
## Root relaxation: objective 1.232570e+01, 4165 iterations, 0.19 seconds
##
## Nodes | Current Node | Objective Bounds | Work
## Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
##
## 0 0 12.32570 0 294 91.00000 12.32570 86.5% - 1s
## H 0 0 13.0000000 12.32570 5.19% - 1s
##
## Explored 1 nodes (4562 simplex iterations) in 1.58 seconds
## Thread count was 2 (of 4 available processors)
##
## Solution count 2: 13 91
##
## Optimal solution found (tolerance 1.00e-01)
## Best objective 1.300000000000e+01, best bound 1.300000000000e+01, gap 0.0000%# show summary
summary(sim_rs_s2_space)## Run_Number Status Score Cost Planning_Units Connectivity_Total
## 1 1 OPTIMAL 13 13 13 220
## Connectivity_In Connectivity_Edge Connectivity_Out
## 1 5 173 42
## Connectivity_In_Fraction
## 1 0.02272727# show amount held
amount.held(sim_rs_s2_space)## normal
## 1 0.2099983# show space held
space.held(sim_rs_s2_space)## normal (Space 1)
## 1 0.8533832# plot the prioritization and the normal species' distribution
spp.plot(sim_rs_s2_space, "normal", main = "Normal species")
A prioritization for the normally distributed species generated using amount-based targets (20%) and space-based targets (85%). See Figure 3 caption for conventions.
# plot the difference between old and new prioritizations
plot(sim_rs_s2_amount, sim_rs_s2_space, 1, 1,
main = "Difference between solutions")
Difference between two prioritizations for the normally distributed species. See Figure 7 caption for conventions.
We can see by using both amount-based and space-based targets we can obtain a prioritization that secures both the species’ range core and parts of its range margin. As a consequence, it may capture any novel adaptations found along the species’ range margin–unlike the amount-based prioritization.
Finally, let’s generate a prioritization for the bimodal species using amount-based and space-based targets.
# make new prioritization
sim_rs_s3_space <- update(sim_rs_s3_amount, amount.target = 0.2,
space.target = 0.85, Threads = threads)## Optimize a model with 10102 rows, 10100 columns and 40000 nonzeros
## Variable types: 0 continuous, 10100 integer (10100 binary)
## Coefficient statistics:
## Matrix range [7e-03, 9e+01]
## Objective range [1e+00, 1e+00]
## Bounds range [1e+00, 1e+00]
## RHS range [1e+00, 7e+01]
## Found heuristic solution: objective 87.0000000
## Presolve removed 589 rows and 15 columns
## Presolve time: 2.37s
## Presolved: 9513 rows, 10085 columns, 51461 nonzeros
## Variable types: 0 continuous, 10085 integer (10085 binary)
## Presolved: 9513 rows, 10085 columns, 51461 nonzeros
##
##
## Root relaxation: objective 8.980792e+00, 5690 iterations, 0.24 seconds
##
## Nodes | Current Node | Objective Bounds | Work
## Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
##
## 0 0 8.98079 0 77 87.00000 8.98079 89.7% - 2s
## H 0 0 10.0000000 8.98079 10.2% - 2s
## 0 0 cutoff 0 10.00000 10.00000 0.00% - 2s
##
## Cutting planes:
## Gomory: 1
## MIR: 1
## Zero half: 1
##
## Explored 1 nodes (5911 simplex iterations) in 2.76 seconds
## Thread count was 2 (of 4 available processors)
##
## Solution count 2: 10 87
##
## Optimal solution found (tolerance 1.00e-01)
## Best objective 1.000000000000e+01, best bound 1.000000000000e+01, gap 0.0000%# show summary
summary(sim_rs_s3_space)## Run_Number Status Score Cost Planning_Units Connectivity_Total
## 1 1 OPTIMAL 10 10 10 220
## Connectivity_In Connectivity_Edge Connectivity_Out
## 1 7 187 26
## Connectivity_In_Fraction
## 1 0.03181818# show amount held
amount.held(sim_rs_s3_space)## bimodal
## 1 0.2226414# show space held
space.held(sim_rs_s3_space)## bimodal (Space 1)
## 1 0.8562891# plot the prioritization and the bimodal species' distribution
spp.plot(sim_rs_s3_space, 'bimodal', main='Bimodal species')
A prioritization for the normally distributed species generated using amount-based targets (20%) and space-based targets (85%). See Figure 3 caption for conventions.
# plot the difference between old and new prioritizations
plot(sim_rs_s3_amount, sim_rs_s3_space, 1, 1,
main = "Difference between solutions")
Difference between two prioritizations for the bimodally distributed species. See Figure 7 caption for conventions.
Earlier we found that the amount-based prioritization only preserved individuals from a single ecotype, and would have failed to adequately preserve the intra-specific variation for this species. However, here we can see that by including space-based targets, we can develop prioritizations that secure individuals belonging to both ecotypes. This new prioritization is much more effective at sampling the intra-specific variation for this species.
Overall, these results demonstrate that under the simplest of conditions, the use of space-based targets can improve prioritizations. However, these prioritizations were each generated for a single species. Prioritizations generated using multiple species may do a better job at preserving the intra-specific variation for individuals species by preserving them in different parts of their range. We will investigate this in the next section.
Multi-species prioritizations
Effects of including space-based targets
So far we have generated prioritizations using only a single species at a time. However, real world prioritizations are often generated using multiple species to ensure that they preserve a comprehensive set of biodiversity. Here, we will generate multi-species prioritizations that preserve all three of the simulated species. First, we will generate a prioritization using amount-based targets that only aims to preserve 20% of the area they occupy. Then, we will generate a prioritization that also incorporate space-based targets to also preserve 85% of their geographic distribution. We will then compare the two prioritizations.
# make prioritizations
sim_mrs_amount <- update(sim_ru, amount.target = c(0.2, 0.2, 0.2),
space.target = c(0, 0, 0), Threads = threads)## Optimize a model with 30306 rows, 30100 columns and 120000 nonzeros
## Variable types: 0 continuous, 30100 integer (30100 binary)
## Coefficient statistics:
## Matrix range [7e-03, 9e+01]
## Objective range [1e+00, 1e+00]
## Bounds range [1e+00, 1e+00]
## RHS range [1e+00, 8e+02]
## Found heuristic solution: objective 98.0000000
## Presolve time: 1.31s
## Presolved: 30306 rows, 30100 columns, 120000 nonzeros
## Variable types: 0 continuous, 30100 integer (30100 binary)
## Presolved: 30306 rows, 30100 columns, 120000 nonzeros
##
##
## Root relaxation: objective 2.000000e+01, 2839 iterations, 1.44 seconds
##
## Nodes | Current Node | Objective Bounds | Work
## Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
##
## H 0 0 20.0000000 0.00000 100% - 2s
## 0 0 - 0 20.00000 20.00000 0.00% - 3s
##
## Explored 0 nodes (4267 simplex iterations) in 3.35 seconds
## Thread count was 2 (of 4 available processors)
##
## Solution count 2: 20 98
##
## Optimal solution found (tolerance 1.00e-01)
## Best objective 2.000000000000e+01, best bound 2.000000000000e+01, gap 0.0000%sim_mrs_space <- update(sim_ru, amount.target = c(0.2, 0.2, 0.2),
space.target = c(0.85, 0.85, 0.85),
Threads = threads)## Optimize a model with 30306 rows, 30100 columns and 120000 nonzeros
## Variable types: 0 continuous, 30100 integer (30100 binary)
## Coefficient statistics:
## Matrix range [7e-03, 9e+01]
## Objective range [1e+00, 1e+00]
## Bounds range [1e+00, 1e+00]
## RHS range [1e+00, 1e+02]
## Found heuristic solution: objective 99.0000000
## Presolve removed 360 rows and 15 columns (presolve time = 5s) ...
## Presolve removed 360 rows and 15 columns
## Presolve time: 5.82s
## Presolved: 29946 rows, 30085 columns, 131504 nonzeros
## Variable types: 0 continuous, 30085 integer (30085 binary)
## Presolve removed 67 rows and 1 columns
## Presolved: 29879 rows, 30084 columns, 127194 nonzeros
##
##
## Root simplex log...
##
## Iteration Objective Primal Inf. Dual Inf. Time
## 0 1.0000000e+02 0.000000e+00 1.000000e+02 6s
## 5887 2.0000000e+01 0.000000e+00 0.000000e+00 8s
## 5887 2.0000000e+01 0.000000e+00 0.000000e+00 8s
##
## Root relaxation: objective 2.000000e+01, 5887 iterations, 2.32 seconds
##
## Nodes | Current Node | Objective Bounds | Work
## Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
##
## H 0 0 20.0000000 0.00000 100% - 9s
## 0 0 - 0 20.00000 20.00000 0.00% - 9s
##
## Explored 0 nodes (11017 simplex iterations) in 9.81 seconds
## Thread count was 2 (of 4 available processors)
##
## Solution count 2: 20 99
##
## Optimal solution found (tolerance 1.00e-01)
## Best objective 2.000000000000e+01, best bound 2.000000000000e+01, gap 0.0000%# show summaries
summary(sim_mrs_amount)## Run_Number Status Score Cost Planning_Units Connectivity_Total
## 1 1 OPTIMAL 20 20 20 220
## Connectivity_In Connectivity_Edge Connectivity_Out
## 1 12 144 64
## Connectivity_In_Fraction
## 1 0.05454545summary(sim_mrs_space)## Run_Number Status Score Cost Planning_Units Connectivity_Total
## 1 1 OPTIMAL 20 20 20 220
## Connectivity_In Connectivity_Edge Connectivity_Out
## 1 8 143 69
## Connectivity_In_Fraction
## 1 0.03636364# show amount held for each prioritization
amount.held(sim_mrs_amount)## uniform normal bimodal
## 1 0.2 0.2134189 0.2270925amount.held(sim_mrs_space)## uniform normal bimodal
## 1 0.2 0.2065024 0.2302924# show space held for each prioritization
space.held(sim_mrs_amount)## uniform (Space 1) normal (Space 1) bimodal (Space 1)
## 1 0.92 0.9010801 0.918714space.held(sim_mrs_space)## uniform (Space 1) normal (Space 1) bimodal (Space 1)
## 1 0.9284848 0.9021208 0.9313511# plot multi-species prioritization with amount-based targets
plot(sim_mrs_amount, 1, main = "Amount-based targets")
A multi-species prioritization for the uniformly, normally, and bimodally distributed species generated using just amount-based targets (20%). Squares represent planning units. Dark green planning units are selected for preservation.
# plot multi-species prioritization with amount- and space-based targets
plot(sim_mrs_space, 1, main = "Amount and space-based targets")
A multi-species prioritization for the uniformly, normally, and bimodally distributed species generated using amount-based targets (20%) and space-based targets (85%). See Figure 12 caption for conventions.
# difference between the two prioritizations
plot(sim_mrs_amount, sim_mrs_space, 1, 1, main = "Difference between solutions")
_Difference between two multi-species prioritizations. See Figure 7 caption for conventions.)
Here we can see that the inclusion of space-based targets changes which planning units are selected for prioritization, but also the number of planning units that are selected. The amount-based prioritization is comprised of 20 units, and the space-based prioritization is comprised of 20 units. This result suggests that an adequate and representative prioritization can be achieved for only a minor increase in cost.
Uncertainty in species’ distributions
The unreliable formulation does not consider the probability that the planning units are occupied by features when calculating how well a given solution secures a representative sample of an attribute space. Thus solutions identified using the unreliable formulation may select regions of an attribute space for a species using planning units that only have a small chance of being inhabited. As a consequence, if the prioritization is implemented, it may fail to secure regions of an attribute space if individuals do not inhabit these planning units, and ultimately fail to fulfill the space-based targets.
A simple solution to this issue would be to ensure that planning units cannot be assigned to demand points if they have a low probability of occupancy. This can be achieved by setting a probability threshold for planning units, such that planning units with a probability of occupancy below the threshold are effectively set to zero.
# make new prioritization with probability threshold of 0.1 for each species
sim_mrs_space2 <- solve(prob.subset(sim_mrs_space, species = 1:3,
threshold = c(0.1, 0.1, 0.1)),
Threads = threads)## Optimize a model with 30306 rows, 30100 columns and 119974 nonzeros
## Variable types: 0 continuous, 30100 integer (30100 binary)
## Coefficient statistics:
## Matrix range [7e-03, 9e+01]
## Objective range [1e+00, 1e+00]
## Bounds range [1e+00, 1e+00]
## RHS range [1e+00, 1e+02]
## Found heuristic solution: objective 99.0000000
## Presolve removed 360 rows and 15 columns (presolve time = 5s) ...
## Presolve removed 360 rows and 15 columns
## Presolve time: 5.45s
## Presolved: 29946 rows, 30085 columns, 131478 nonzeros
## Variable types: 0 continuous, 30085 integer (30085 binary)
## Presolve removed 67 rows and 1 columns
## Presolved: 29879 rows, 30084 columns, 127168 nonzeros
##
##
## Root simplex log...
##
## Iteration Objective Primal Inf. Dual Inf. Time
## 0 1.0000000e+02 0.000000e+00 1.000000e+02 6s
## 5510 2.0000000e+01 0.000000e+00 0.000000e+00 7s
## 5510 2.0000000e+01 0.000000e+00 0.000000e+00 7s
##
## Root relaxation: objective 2.000000e+01, 5510 iterations, 1.44 seconds
##
## Nodes | Current Node | Objective Bounds | Work
## Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
##
## H 0 0 20.0000000 0.00000 100% - 7s
## 0 0 - 0 20.00000 20.00000 0.00% - 7s
##
## Explored 0 nodes (8707 simplex iterations) in 7.43 seconds
## Thread count was 2 (of 4 available processors)
##
## Solution count 2: 20 99
##
## Optimal solution found (tolerance 1.00e-01)
## Best objective 2.000000000000e+01, best bound 2.000000000000e+01, gap 0.0000%# show summary
summary(sim_mrs_space2)## Run_Number Status Score Cost Planning_Units Connectivity_Total
## 1 1 OPTIMAL 20 20 20 220
## Connectivity_In Connectivity_Edge Connectivity_Out
## 1 6 141 73
## Connectivity_In_Fraction
## 1 0.02727273# plot prioritization
plot(sim_mrs_space2, 1)
A multi-species prioritization for the uniformly, normally, and bimodally distributed species generated using amount-based targets (20%) and space-based targets (85%). This prioritization was generated to be robust against low occupancy probabilities, by preventing planning units with low probabilities from being used to represent demand points. See Figure 12 caption for conventions.
# difference between prioritizations that use and do not use thresholds
plot(sim_mrs_space2, sim_mrs_space, 1, 1, main = "Difference between solutions")
Difference between two multi-species prioritizations. See Figure 7 caption for conventions.
But this method requires setting somewhat arbitrary thresholds. A more robust solution to this issue is to actually use the probability that species occupy planning units to generate the prioritizations. This is what the reliable formulation does. First we will try and generate a solution using existing targets and the reliable formulation. To reduce computational time, we will set the maximum backup \(R\)-level to 1.
# make new prioritization using reliable formulation
sim_mrs_space3 <- try(update(sim_mrs_space, formulation = "reliable",
max.r.level = 1L, Threads = threads))## Optimize a model with 364206 rows, 181900 columns and 3847200 nonzeros
## Variable types: 0 continuous, 60700 integer (60700 binary)
## Semi-Variable types: 121200 continuous, 0 integer
## Coefficient statistics:
## Matrix range [8e-04, 1e+02]
## Objective range [1e+00, 1e+00]
## Bounds range [1e+00, 1e+00]
## RHS range [7e-03, 1e+02]
## Presolve removed 338506 rows and 156703 columns
## Presolve time: 2.95s
##
## Explored 0 nodes (0 simplex iterations) in 3.30 seconds
## Thread count was 1 (of 4 available processors)
##
## Solution count 0
##
## Model is infeasible
## Best objective -, best bound -, gap -However, this fails. The reason why we cannot generate a prioritization that fulfills these targets is because even the solution that contains all the planning units is still insufficient when we consider probabilities. The negative maximum targets printed in the error message indicate that planning units have low probabilities of occupancy. To fulfill the targets, we must obtain more planning units with higher probabilities of occupancy. We also could attempt resolving the problem using a higher \(R\)-level. Instead, we will set lower targets and generate solution.
# make new prioritization using reliable formulation and reduced targets
sim_mrs_space3 <- update(sim_mrs_space, formulation = "reliable",
max.r.level = 1L, space.target = -1000,
Threads = threads)## Optimize a model with 364206 rows, 181900 columns and 3847200 nonzeros
## Variable types: 0 continuous, 60700 integer (60700 binary)
## Semi-Variable types: 121200 continuous, 0 integer
## Coefficient statistics:
## Matrix range [8e-04, 1e+02]
## Objective range [1e+00, 1e+00]
## Bounds range [1e+00, 1e+00]
## RHS range [7e-03, 8e+05]
## Presolve removed 333903 rows and 151800 columns (presolve time = 5s) ...
## Presolve removed 333903 rows and 151800 columns
## Presolve time: 5.69s
## Presolved: 30303 rows, 30100 columns, 90300 nonzeros
## Variable types: 0 continuous, 30100 integer (30100 binary)
## Found heuristic solution: objective 27.0000000
## Presolved: 30303 rows, 30100 columns, 90300 nonzeros
##
##
## Root simplex log...
##
## Iteration Objective Primal Inf. Dual Inf. Time
## 0 1.0000000e+02 0.000000e+00 1.000000e+02 6s
## 129 2.0000000e+01 0.000000e+00 0.000000e+00 6s
## 129 2.0000000e+01 0.000000e+00 0.000000e+00 6s
##
## Root relaxation: objective 2.000000e+01, 129 iterations, 0.28 seconds
##
## Nodes | Current Node | Objective Bounds | Work
## Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
##
## 0 0 20.00000 0 3 27.00000 20.00000 25.9% - 6s
## H 0 0 20.0000000 20.00000 0.00% - 6s
##
## Explored 1 nodes (129 simplex iterations) in 6.39 seconds
## Thread count was 2 (of 4 available processors)
##
## Solution count 2: 20 27
##
## Optimal solution found (tolerance 1.00e-01)
## Best objective 2.000000000000e+01, best bound 2.000000000000e+01, gap 0.0000%## Warning in validityMethod(object): some species have space.held values less
## than 0, and thus are poorly represented# show summary
summary(sim_mrs_space3)## Run_Number Status Score Cost Planning_Units Connectivity_Total
## 1 1 OPTIMAL 20 20 20 220
## Connectivity_In Connectivity_Edge Connectivity_Out
## 1 24 160 36
## Connectivity_In_Fraction
## 1 0.1090909# plot prioritization
plot(sim_mrs_space3, 1)
A multi-species prioritization for the uniformly, normally, and bimodally distributed species generated using amount-based targets (20%) and space-based targets (85%). This prioritization was generated to be robust against low occupancy probabilities, by explicitly using the probability of occupancy data when deriving a solution. See Figure 12 caption for conventions.
# difference between prioritizations based on unreliable and reliable formulation
plot(sim_mrs_space3, sim_mrs_space, 1, 1, main = "Difference between solutions")
Difference between two multi-species prioritizations. See Figure 7 caption for conventions.
An additional planning unit was selected using the reliable formulation. The prioritization based on the unreliable formulation had 20 planning units, but the prioritization based on the reliable formulation has 20 planning units. This difference occurs because the reliable formulation needs to ensure that all selected planning units with a low chance of being occupied have a suitable backup planning unit. While the reliable formulation can deliver more robust prioritizations, it takes much longer to solve conservation planning problems expressed using this formulation than the unreliable formulation. As a consequence, the reliable formulation is only feasible for particularly small problems, such as those involving few features and less than several hundred planning units.
Fragmentation
Fragmentation is an important consideration in real-world planning situations. Up until now, we haven’t considered the effects of fragmentation on the viability of the prioritization. As a consequence, our prioritizations have tended to contain planning units without any neighbors. We can use the BLM parameter to penalize fragmented solutions.
Let’s generate a new prioritization that heavily penalises fragmentation. Here, we will update the sim_mrs_amount object with BLM of 100.
# update prioritization
sim_mrs_amount_blm <- update(sim_mrs_amount, BLM = 100, Threads = threads)## Optimize a model with 30666 rows, 30280 columns and 120720 nonzeros
## Variable types: 0 continuous, 30280 integer (30280 binary)
## Coefficient statistics:
## Matrix range [7e-03, 9e+01]
## Objective range [1e+02, 4e+02]
## Bounds range [1e+00, 1e+00]
## RHS range [1e+00, 8e+02]
## Found heuristic solution: objective 4498.0000000
## Presolve time: 1.09s
## Presolved: 30666 rows, 30280 columns, 120720 nonzeros
## Variable types: 0 continuous, 30280 integer (30280 binary)
## Presolved: 30666 rows, 30280 columns, 120720 nonzeros
##
##
## Root relaxation: objective 4.644444e+02, 4432 iterations, 1.69 seconds
##
## Nodes | Current Node | Objective Bounds | Work
## Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
##
## 0 0 464.44444 0 1725 4498.00000 464.44444 89.7% - 3s
## H 0 0 2130.0000000 464.44444 78.2% - 10s
## 0 0 537.57397 0 1665 2130.00000 537.57397 74.8% - 11s
## 0 0 537.57397 0 1664 2130.00000 537.57397 74.8% - 11s
## 0 0 571.03030 0 1531 2130.00000 571.03030 73.2% - 12s
## 0 0 574.69174 0 1761 2130.00000 574.69174 73.0% - 12s
## 0 0 574.69174 0 1765 2130.00000 574.69174 73.0% - 12s
## 0 0 574.69174 0 1825 2130.00000 574.69174 73.0% - 13s
## 0 0 574.69174 0 1828 2130.00000 574.69174 73.0% - 13s
## 0 0 574.69174 0 1818 2130.00000 574.69174 73.0% - 13s
## 0 0 574.69174 0 1818 2130.00000 574.69174 73.0% - 13s
## 0 0 574.69174 0 1478 2130.00000 574.69174 73.0% - 16s
## H 0 0 1335.0000000 574.69174 57.0% - 18s
## 0 0 574.69174 0 1477 1335.00000 574.69174 57.0% - 18s
## 0 0 574.69174 0 1488 1335.00000 574.69174 57.0% - 21s
## H 0 0 1333.0000000 574.69174 56.9% - 23s
## 0 0 574.69174 0 1487 1333.00000 574.69174 56.9% - 23s
## 0 0 574.69174 0 1487 1333.00000 574.69174 56.9% - 25s
## H 0 0 1332.0000000 574.69174 56.9% - 27s
## 0 0 574.69174 0 1487 1332.00000 574.69174 56.9% - 27s
## 0 0 574.69174 0 1480 1332.00000 574.69174 56.9% - 29s
## 0 0 574.69174 0 1481 1332.00000 574.69174 56.9% - 31s
## 0 0 574.69174 0 1480 1332.00000 574.69174 56.9% - 33s
## 0 0 574.69174 0 1441 1332.00000 574.69174 56.9% - 34s
## H 0 2 1020.0000000 574.69174 43.7% - 39s
## 0 2 574.69174 0 1363 1020.00000 574.69174 43.7% - 39s
## 1 4 608.32371 1 1477 1020.00000 581.13743 43.0% 8064 40s
## 75 55 934.28571 8 589 1020.00000 648.92237 36.4% 939 46s
## 134 75 968.81234 9 1236 1020.00000 697.28516 31.6% 955 50s
## 278 107 1016.82857 14 881 1020.00000 752.93755 26.2% 758 55s
## 430 145 cutoff 15 1020.00000 783.90493 23.1% 686 60s
## 595 202 949.97821 10 576 1020.00000 804.33214 21.1% 621 65s
## H 712 124 920.0000000 810.24108 11.9% 596 68s
## 804 118 cutoff 7 920.00000 824.09505 10.4% 570 71s
##
## Cutting planes:
## Gomory: 3
## Zero half: 33
##
## Explored 833 nodes (582370 simplex iterations) in 71.15 seconds
## Thread count was 2 (of 4 available processors)
##
## Solution count 7: 920 1020 1332 ... 4498
##
## Optimal solution found (tolerance 1.00e-01)
## Best objective 9.200000000000e+02, best bound 8.310000000000e+02, gap 9.6739%# show summary of prioritization
summary(sim_mrs_amount_blm)## Run_Number Status Score Cost Planning_Units Connectivity_Total
## 1 1 OPTIMAL 920 20 20 220
## Connectivity_In Connectivity_Edge Connectivity_Out
## 1 31 171 18
## Connectivity_In_Fraction
## 1 0.1409091# show amount held for each prioritization
amount.held(sim_mrs_amount_blm)## uniform normal bimodal
## 1 0.2 0.2736046 0.2009364# show space held for each prioritization
space.held(sim_mrs_amount_blm)## uniform (Space 1) normal (Space 1) bimodal (Space 1)
## 1 0.4727273 0.4694515 0.3551407# plot prioritization
plot(sim_mrs_amount_blm, 1)
A multi-species prioritization for the uniformly, normally, and bimodally distributed species generated using only amount-based targets (20%). Additionally, this prioritization was specified to have high connectivity, by using a high \(BLM\) parameter. See Figure 12 caption for conventions.
# difference between the two prioritizations
plot(sim_mrs_amount_blm, sim_mrs_amount, 1, 1,
main = "Difference between solutions")
Difference between two multi-species prioritizations. See Figure 7 caption for conventions.
Here we can see that the prioritization generated using a BLM parameter of 100 is much more clustered than the prioritization generated using a BLM of 0. In practice, conservation planners will need to try a variety of BLM parameters to find a suitable prioritization.
