Skip to main content

2021 UQ Spatial Planning Workshop using Prioritizr

6 Spatial prioritizations

6.1 Introduction

Here we will develop prioritizations to identify priority areas for protected area establishment.

6.2 Starting out simple

To start things off, let’s keep things simple. Let’s create a prioritization using the minimum set formulation of the reserve selection problem. This formulation means that we want a solution that will meet the targets for our biodiversity features for minimum cost. Here, we will set 5% targets for each vegetation class and use the data in the cost column to specify acquisition costs. Unlike Marxan, we do not have to calibrate species penalty factors (SPFs) to ensure that our target are met—prioritizr should always return solutions to minimum set problems where all the targets are met.

Although we strongly recommend using Gurobi to solve problems (with add_gurobi_solver), you can also use the lpsymphony solver in this workshop (it’s much easier to install). The Gurobi solver is much faster than the lpsymphony solver (see here for installation instructions).

# print planning unit data
print(tas_pu)
## class       : SpatialPolygonsDataFrame 
## features    : 1130 
## extent      : 1080623, 1399989, -4840595, -4497092  (xmin, xmax, ymin, ymax)
## crs         : +proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
## variables   : 5
## names       :   id,              cost, status, locked_in, locked_out 
## min values  :    1, 0.192488262910798,      0,         0,          0 
## max values  : 1130,  61.9272727272727,      2,         1,          1
# make prioritization problem
p1 <- problem(tas_pu, tas_features, cost_column = "cost") %>%
      add_min_set_objective() %>%
      add_relative_targets(0.05) %>% # 5% representation targets
      add_binary_decisions() %>%
      add_gurobi_solver(verbose = FALSE)
      # add_lpsymphony_solver(verbose = FALSE)

# print problem
print(p1)
## Conservation Problem
##   planning units: SpatialPolygonsDataFrame (1130 units)
##   cost:           min: 0.19249, max: 61.92727
##   features:       vegetation.1, vegetation.2, vegetation.3, ... (62 features)
##   objective:      Minimum set objective 
##   targets:        Relative targets [targets (min: 0.05, max: 0.05)]
##   decisions:      Binary decision 
##   constraints:    <none>
##   penalties:      <none>
##   portfolio:      default
##   solver:         Gurobi [first_feasible (0), gap (0.1), numeric_focus (0), presolve (2), threads (1), time_limit (2147483647), verbose (0)]
# solve problem
s1 <- solve(p1)

# print solution, the solution_1 column contains the solution values
# indicating if a planning unit is (1) selected or (0) not
print(s1)
## class       : SpatialPolygonsDataFrame 
## features    : 1130 
## extent      : 1080623, 1399989, -4840595, -4497092  (xmin, xmax, ymin, ymax)
## crs         : +proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
## variables   : 6
## names       :   id,              cost, status, locked_in, locked_out, solution_1 
## min values  :    1, 0.192488262910798,      0,         0,          0,          0 
## max values  : 1130,  61.9272727272727,      2,         1,          1,          1
# calculate number of planning units selected in the prioritization
sum(s1$solution_1)
## [1] 35
# calculate total cost of the prioritization
sum(s1$solution_1 * s1$cost)
## [1] 808.4013
# plot solution
plot_solution(s1)

Now let’s examine the solution.

  1. How many planing units were selected in the prioritization? What proportion of planning units were selected in the prioritization?
  2. Is there a pattern in the spatial distribution of the priority areas?
  3. Can you verify that all of the targets were met in the prioritization (hint: eval_feature_representation_summary(p1, s1[, "solution_1"]))?

6.3 Adding complexity

Our first prioritization suffers many limitations, so let’s add additional constraints to the problem to make it more useful.

First, let’s lock in planing units that are already by covered protected areas. If some vegetation communities are already secured inside existing protected areas, then we might not need to add as many new protected areas to the existing protected area system to meet their targets. Since our planning unit data (tas_pu) already contains this information in the locked_in column, we can use this column name to specify which planning units should be locked in.

# make prioritization problem
p2 <- problem(tas_pu, tas_features, cost_column = "cost") %>%
      add_min_set_objective() %>%
      add_relative_targets(0.05) %>%
      add_locked_in_constraints("locked_in") %>%
      add_binary_decisions() %>%
      add_gurobi_solver(verbose = FALSE)
      # add_lpsymphony_solver(verbose = FALSE)

# print problem
print(p2)
## Conservation Problem
##   planning units: SpatialPolygonsDataFrame (1130 units)
##   cost:           min: 0.19249, max: 61.92727
##   features:       vegetation.1, vegetation.2, vegetation.3, ... (62 features)
##   objective:      Minimum set objective 
##   targets:        Relative targets [targets (min: 0.05, max: 0.05)]
##   decisions:      Binary decision 
##   constraints:    <Locked in planning units [257 locked units]>
##   penalties:      <none>
##   portfolio:      default
##   solver:         Gurobi [first_feasible (0), gap (0.1), numeric_focus (0), presolve (2), threads (1), time_limit (2147483647), verbose (0)]
# solve problem
s2 <- solve(p2)

# plot solution
plot_solution(s2)

Let’s pretend that we talked to an expert on the vegetation communities in our study system and they recommended that a 20% target was needed for each vegetation class. So, armed with this information, let’s set the targets to 20%.

# make prioritization problem
p3 <- problem(tas_pu, tas_features, cost_column = "cost") %>%
      add_min_set_objective() %>%
      add_relative_targets(0.2) %>%
      add_locked_in_constraints("locked_in") %>%
      add_binary_decisions() %>%
      add_gurobi_solver(verbose = FALSE)
      # add_lpsymphony_solver(verbose = FALSE)

# print problem
print(p3)
## Conservation Problem
##   planning units: SpatialPolygonsDataFrame (1130 units)
##   cost:           min: 0.19249, max: 61.92727
##   features:       vegetation.1, vegetation.2, vegetation.3, ... (62 features)
##   objective:      Minimum set objective 
##   targets:        Relative targets [targets (min: 0.2, max: 0.2)]
##   decisions:      Binary decision 
##   constraints:    <Locked in planning units [257 locked units]>
##   penalties:      <none>
##   portfolio:      default
##   solver:         Gurobi [first_feasible (0), gap (0.1), numeric_focus (0), presolve (2), threads (1), time_limit (2147483647), verbose (0)]
# solve problem
s3 <- solve(p3)

# plot solution
plot_solution(s3)

Next, let’s lock out highly degraded areas. Similar to before, this data is present in our planning unit data so we can use the locked_out column name to achieve this.

# make prioritization problem
p4 <- problem(tas_pu, tas_features, cost_column = "cost") %>%
      add_min_set_objective() %>%
      add_relative_targets(0.2) %>%
      add_locked_in_constraints("locked_in") %>%
      add_locked_out_constraints("locked_out") %>%
      add_binary_decisions() %>%
      add_gurobi_solver(verbose = FALSE)
      # add_lpsymphony_solver(verbose = FALSE)
# print problem
print(p4)
## Conservation Problem
##   planning units: SpatialPolygonsDataFrame (1130 units)
##   cost:           min: 0.19249, max: 61.92727
##   features:       vegetation.1, vegetation.2, vegetation.3, ... (62 features)
##   objective:      Minimum set objective 
##   targets:        Relative targets [targets (min: 0.2, max: 0.2)]
##   decisions:      Binary decision 
##   constraints:    <Locked in planning units [257 locked units]
##                    Locked out planning units [45 locked units]>
##   penalties:      <none>
##   portfolio:      default
##   solver:         Gurobi [first_feasible (0), gap (0.1), numeric_focus (0), presolve (2), threads (1), time_limit (2147483647), verbose (0)]
# solve problem
s4 <- solve(p4)

# plot solution
spplot(s4, "solution_1", col.regions = c("white", "darkgreen"), main = "s4")
plot_solution(s4)

Now, let’s compare the solutions.

  1. What is the cost of the planning units selected in s2, s3, and s4?
  2. How many planning units are in s2, s3, and s4?
  3. Do the solutions with more planning units have a greater cost? Why or why not?
  4. Why does the first solution (s1) cost less than the second solution with protected areas locked into the solution (s2)?
  5. Why does the third solution (s3) cost less than the fourth solution solution with highly degraded areas locked out (s4)?
  6. Since planning units covered by existing protected areas have already been purchased, what is the cost for expanding the protected area system based on on the fourth prioritization (s4) (hint: total cost minus the cost of locked in planning units)?
  7. What happens if you specify targets that exceed the total amount of vegetation in the study area and try to solve the problem? You can do this by modifying the code to make p4 with add_absolute_targets(1000) instead of add_relative_targets(0.2) and generating a new solution.

6.4 Penalizing fragmentation

Plans for protected area systems should facilitate gene flow and dispersal between individual reserves in the system. However, the prioritizations we have made so far have been highly fragmented.

We can add penalties to our conservation planning problem to penalize fragmentation (i.e. total exposed boundary length) and we also need to set a useful penalty value when adding such penalties. If we set our penalty value too low, then we will end up with a solution that is identical to the solution with no added penalties. If we set our penalty value too high, then prioritizr will take a long time to solve the problem and we will end up with a solution that contains lots of extra planning units that are not needed (since the penalty value is so high that minimizing fragmentation is more important than cost). As a rule of thumb, we generally want penalty values between 0.00001 and 0.01 but finding a useful penalty value requires calibration. The “correct” penalty value depends on the size of the planning units, the main objective values (e.g. cost values), and the effect of fragmentation on biodiversity persistence.

Let’s create a new problem that is similar to our previous problem (p4)—except that it contains boundary length penalties and a slightly higher optimality gap to reduce runtime (default is 0.1)—and solve it. Since our planning unit data is in a spatial format (i.e. vector or raster data), prioritizr can automatically calculate the boundary data for us.

# make prioritization problem
p5 <- problem(tas_pu, tas_features, cost_column = "cost") %>%
      add_min_set_objective() %>%
      add_boundary_penalties(penalty = 0.0005) %>%
      add_relative_targets(0.2) %>%
      add_locked_in_constraints("locked_in") %>%
      add_locked_out_constraints("locked_out") %>%
      add_binary_decisions() %>%
      add_gurobi_solver(verbose = FALSE)
      # add_lpsymphony_solver(verbose = FALSE)

# print problem
print(p5)
## Conservation Problem
##   planning units: SpatialPolygonsDataFrame (1130 units)
##   cost:           min: 0.19249, max: 61.92727
##   features:       vegetation.1, vegetation.2, vegetation.3, ... (62 features)
##   objective:      Minimum set objective 
##   targets:        Relative targets [targets (min: 0.2, max: 0.2)]
##   decisions:      Binary decision 
##   constraints:    <Locked out planning units [45 locked units]
##                    Locked in planning units [257 locked units]>
##   penalties:      <Boundary penalties [edge factor (min: 0.5, max: 0.5), penalty (5e-04), zones]>
##   portfolio:      default
##   solver:         Gurobi [first_feasible (0), gap (0.1), numeric_focus (0), presolve (2), threads (1), time_limit (2147483647), verbose (0)]
# solve problem,
# note this will take around 30 seconds
s5 <- solve(p5)

# print solution
print(s5)
## class       : SpatialPolygonsDataFrame 
## features    : 1130 
## extent      : 1080623, 1399989, -4840595, -4497092  (xmin, xmax, ymin, ymax)
## crs         : +proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
## variables   : 6
## names       :   id,              cost, status, locked_in, locked_out, solution_1 
## min values  :    1, 0.192488262910798,      0,         0,          0,          0 
## max values  : 1130,  61.9272727272727,      2,         1,          1,          1
# plot solution
plot_solution(s5)

Now let’s compare the solutions to the problems with (s5) and without (s4) the boundary length penalties.

  1. What is the cost the fourth (s4) and fifth (s5) solutions? Why does the fifth solution (s5) cost more than the fourth (s4) solution?
  2. Try setting the penalty value to 0.000000001 (i.e. 1e-9) instead of 0.0005. What is the cost of the solution now? Is it different from the fourth solution (s4) (hint: try plotting the solutions to visualize them)? Is this is a useful penalty value? Why?
  3. Try setting the penalty value to 0.5. What is the cost of the solution now? Is it different from the fourth solution (s4) (hint: try plotting the solutions to visualize them)? Is this a useful penalty value? Why?

6.5 Budget limited prioritizations

In the real-world, the funding available for conservation is often very limited. As a consequence, decision makers often need prioritizations where the total cost of priority areas does not exceed a budget.

In our fourth prioritization (s4), we found that we would need to spend an additional $1317 million AUD to ensure that each vegetation community is adequately represented in the protected area system. But what if the funds available for establishing new protected areas were limited to $100 million AUD? In this case, we need a “budget limited prioritization”. Budget limited prioritizations aim to maximize some measure of conservation benefit subject to a budget (e.g. number of species with at least one occurrence in the protected area system, or phylogenetic diversity). Let’s create a prioritization by maximizing the number of adequately represented features whilst keeping within a pre-specified budget.

# funds for additional land acquisition (same units as cost data)
funds <- 100

# calculate the total budget for the prioritization
budget <- funds + sum(s4$cost * s4$locked_in)
print(budget)
## [1] 8575.56
# make prioritization problem
p6 <- problem(tas_pu, tas_features, cost_column = "cost") %>%
      add_max_features_objective(budget) %>%
      add_relative_targets(0.2) %>%
      add_locked_in_constraints("locked_in") %>%
      add_locked_out_constraints("locked_out") %>%
      add_binary_decisions() %>%
      add_gurobi_solver(verbose = FALSE)
      # add_lpsymphony_solver(verbose = FALSE)

# print problem
print(p6)
## Conservation Problem
##   planning units: SpatialPolygonsDataFrame (1130 units)
##   cost:           min: 0.19249, max: 61.92727
##   features:       vegetation.1, vegetation.2, vegetation.3, ... (62 features)
##   objective:      Maximum representation objective [budget (8575.56009869836)]
##   targets:        Relative targets [targets (min: 0.2, max: 0.2)]
##   decisions:      Binary decision 
##   constraints:    <Locked in planning units [257 locked units]
##                    Locked out planning units [45 locked units]>
##   penalties:      <none>
##   portfolio:      default
##   solver:         Gurobi [first_feasible (0), gap (0.1), numeric_focus (0), presolve (2), threads (1), time_limit (2147483647), verbose (0)]
# solve problem
s6 <- solve(p6)

# plot solution
plot_solution(s6)
# calculate feature representation
r6 <- eval_feature_representation_summary(p6, s6[, "solution_1"])

# calculate number of features with targets met
sum(r6$relative_held >= 0.2, na.rm = TRUE)
## [1] 28
# find out which features have their targets met when we add weights,
# note that NA is for vegetation.61
print(r6$feature[r6$relative_held >= 0.2])
##  [1] "vegetation.1"  "vegetation.2"  "vegetation.3"  "vegetation.4"  "vegetation.5"  "vegetation.6"  "vegetation.7"  "vegetation.8"  "vegetation.11" "vegetation.12" "vegetation.15"
## [12] "vegetation.17" "vegetation.22" "vegetation.23" "vegetation.25" "vegetation.28" "vegetation.29" "vegetation.30" "vegetation.32" "vegetation.33" "vegetation.34" "vegetation.35"
## [23] "vegetation.36" "vegetation.37" "vegetation.38" "vegetation.39" "vegetation.40" "vegetation.45"

We can also add weights to specify that it is more important to meet the targets for certain features and less important for other features. A common approach for weighting features is to assign a greater importance to features with smaller spatial distributions. The rationale behind this weighting method is that features with smaller spatial distributions are at greater risk of extinction.

So, let’s calculate some weights for our vegetation communities and see how weighting the features changes our prioritization.

# calculate weights as the log inverse number of grid cells that each vegetation
# class occupies, rescaled between 1 and 100
wts <- 1 / cellStats(tas_features, "sum")
wts <- rescale(wts, to = c(1, 10))

# print the name of the feature with smallest weight
names(tas_features)[which.min(wts)]
## [1] "vegetation.20"
# print the name of the feature with greatest weight
names(tas_features)[which.max(wts)]
## [1] "vegetation.52"
# plot histogram of weights
hist(wts, main = "feature weights")
# make prioritization problem with weights
p7 <- problem(tas_pu, tas_features, cost_column = "cost") %>%
      add_max_features_objective(budget) %>%
      add_relative_targets(0.2) %>%
      add_feature_weights(wts) %>%
      add_locked_in_constraints("locked_in") %>%
      add_locked_out_constraints("locked_out") %>%
      add_binary_decisions() %>%
       add_gurobi_solver(verbose = FALSE)
      # add_lpsymphony_solver(verbose = FALSE)

# print problem
print(p7)
## Conservation Problem
##   planning units: SpatialPolygonsDataFrame (1130 units)
##   cost:           min: 0.19249, max: 61.92727
##   features:       vegetation.1, vegetation.2, vegetation.3, ... (62 features)
##   objective:      Maximum representation objective [budget (8575.56009869836)]
##   targets:        Relative targets [targets (min: 0.2, max: 0.2)]
##   decisions:      Binary decision 
##   constraints:    <Locked in planning units [257 locked units]
##                    Locked out planning units [45 locked units]>
##   penalties:      <Feature weights [weights (min: 1, max: 10)]>
##   portfolio:      default
##   solver:         Gurobi [first_feasible (0), gap (0.1), numeric_focus (0), presolve (2), threads (1), time_limit (2147483647), verbose (0)]
# solve problem
s7 <- solve(p7)

# plot solution
plot_solution(s7)
# calculate feature representation
r7 <- eval_feature_representation_summary(p7, s7[, "solution_1"])

# calculate number of features with targets met
sum(r7$relative_held >= 0.2, na.rm = TRUE)
## [1] 22
# find out which features have their targets met when we add weights,
# note that NA is for vegetation.61
print(r7$feature[r7$relative_held >= 0.2])
##  [1] "vegetation.2"  "vegetation.6"  "vegetation.11" "vegetation.23" "vegetation.28" "vegetation.29" "vegetation.30" "vegetation.32" "vegetation.33" "vegetation.34" "vegetation.35"
## [12] "vegetation.36" "vegetation.37" "vegetation.38" "vegetation.39" "vegetation.40" "vegetation.45" "vegetation.49" "vegetation.50" "vegetation.52" "vegetation.53" "vegetation.61"
  1. What is the name of the feature with the smallest weight?
  2. What is the cost of the sixth (s6) and seventh (s7) solutions?
  3. Does there seem to be a big difference in which planning units were selected in the sixth (s6) and seventh (s7) solutions?
  4. Is there a difference between which features are adequately represented in the sixth (s6) and seventh (s7) solutions? If so, what is the difference?

6.6 Solution portfolios

In systematic conservation planning, only rarely do we have data on all of the stakeholder preferences and biodiversity features that we are interested in conserving. As a consequence, it is often useful to generate a portfolio of near optimal solutions to present to decision makers to guide the reserve selection process.

Generally we would want many solutions in our portfolio (e.g. 1000) to ensure that our portfolio contains a range of spatially distinct solutions, but here we will generate a portfolio containing just six near-optimal solutions so the code doesn’t take too long to run. We will also increase the optimality gap to obtain solutions that are more suboptimal than earlier (the default gap value is 0.1).

# make problem with a shuffle portfolio
p8 <- problem(tas_pu, tas_features, cost_column = "cost") %>%
      add_max_features_objective(budget) %>%
      add_relative_targets(0.2) %>%
      add_feature_weights(wts) %>%
      add_binary_decisions() %>%
      add_shuffle_portfolio(number_solutions = 6,
                            remove_duplicates = FALSE) %>%
      add_gurobi_solver(verbose = TRUE, gap = 10)
      # add_lpsymphony_solver(verbose = TRUE, gap = 10)
# print problem
print(p8)
## Conservation Problem
##   planning units: SpatialPolygonsDataFrame (1130 units)
##   cost:           min: 0.19249, max: 61.92727
##   features:       vegetation.1, vegetation.2, vegetation.3, ... (62 features)
##   objective:      Maximum representation objective [budget (8575.56009869836)]
##   targets:        Relative targets [targets (min: 0.2, max: 0.2)]
##   decisions:      Binary decision 
##   constraints:    <none>
##   penalties:      <Feature weights [weights (min: 1, max: 10)]>
##   portfolio:      Shuffle portfolio [number_solutions (6), remove_duplicates (0), threads (1)]
##   solver:         Gurobi [first_feasible (0), gap (10), numeric_focus (0), presolve (2), threads (1), time_limit (2147483647), verbose (1)]
# solve problem
# note that this will contain six solutions since we added a portfolio
s8 <- solve(p8)
## Gurobi Optimizer version 9.1.1 build v9.1.1rc0 (mac64)
## Thread count: 4 physical cores, 8 logical processors, using up to 1 threads
## Optimize a model with 63 rows, 1192 columns and 3308 nonzeros
## Model fingerprint: 0x0e2ebd61
## Variable types: 0 continuous, 1192 integer (1192 binary)
## Coefficient statistics:
##   Matrix range     [1e-05, 2e+02]
##   Objective range  [7e-08, 1e+01]
##   Bounds range     [1e+00, 1e+00]
##   RHS range        [9e+03, 9e+03]
## Found heuristic solution: objective -0.0000000
## Presolve removed 0 rows and 143 columns
## Presolve time: 0.02s
## Presolved: 63 rows, 1049 columns, 3152 nonzeros
## Variable types: 0 continuous, 1049 integer (1049 binary)
## Presolved: 63 rows, 1049 columns, 3152 nonzeros
## 
## 
## Root relaxation: objective 1.319894e+02, 63 iterations, 0.00 seconds
## 
##     Nodes    |    Current Node    |     Objective Bounds      |     Work
##  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
## 
##      0     0  131.98943    0   51   -0.00000  131.98943      -     -    0s
## H    0     0                      77.7382485  131.98943  69.8%     -    0s
## 
## Explored 1 nodes (63 simplex iterations) in 0.03 seconds
## Thread count was 1 (of 8 available processors)
## 
## Solution count 2: 77.7382 -0 
## 
## Optimal solution found (tolerance 1.00e+01)
## Best objective 7.773824845163e+01, best bound 1.319894286177e+02, gap 69.7870%
## Gurobi Optimizer version 9.1.1 build v9.1.1rc0 (mac64)
## Thread count: 4 physical cores, 8 logical processors, using up to 1 threads
## Optimize a model with 63 rows, 1192 columns and 3308 nonzeros
## Model fingerprint: 0x440aa218
## Variable types: 0 continuous, 1192 integer (1192 binary)
## Coefficient statistics:
##   Matrix range     [1e-05, 2e+02]
##   Objective range  [7e-08, 1e+01]
##   Bounds range     [1e+00, 1e+00]
##   RHS range        [9e+03, 9e+03]
## Found heuristic solution: objective -0.0000000
## Presolve removed 0 rows and 143 columns
## Presolve time: 0.02s
## Presolved: 63 rows, 1049 columns, 3152 nonzeros
## Variable types: 0 continuous, 1049 integer (1049 binary)
## Presolved: 63 rows, 1049 columns, 3152 nonzeros
## 
## 
## Root relaxation: objective 1.319894e+02, 65 iterations, 0.00 seconds
## 
##     Nodes    |    Current Node    |     Objective Bounds      |     Work
##  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
## 
##      0     0  131.98943    0   50   -0.00000  131.98943      -     -    0s
## H    0     0                      86.8777788  131.98943  51.9%     -    0s
## 
## Explored 1 nodes (65 simplex iterations) in 0.03 seconds
## Thread count was 1 (of 8 available processors)
## 
## Solution count 2: 86.8778 -0 
## 
## Optimal solution found (tolerance 1.00e+01)
## Best objective 8.687777882953e+01, best bound 1.319894280514e+02, gap 51.9254%
## Gurobi Optimizer version 9.1.1 build v9.1.1rc0 (mac64)
## Thread count: 4 physical cores, 8 logical processors, using up to 1 threads
## Optimize a model with 63 rows, 1192 columns and 3308 nonzeros
## Model fingerprint: 0xaecc1fbe
## Variable types: 0 continuous, 1192 integer (1192 binary)
## Coefficient statistics:
##   Matrix range     [1e-05, 2e+02]
##   Objective range  [7e-08, 1e+01]
##   Bounds range     [1e+00, 1e+00]
##   RHS range        [9e+03, 9e+03]
## Found heuristic solution: objective -0.0000000
## Presolve removed 0 rows and 143 columns
## Presolve time: 0.02s
## Presolved: 63 rows, 1049 columns, 3152 nonzeros
## Variable types: 0 continuous, 1049 integer (1049 binary)
## Presolved: 63 rows, 1049 columns, 3152 nonzeros
## 
## 
## Root relaxation: objective 1.319894e+02, 63 iterations, 0.00 seconds
## 
##     Nodes    |    Current Node    |     Objective Bounds      |     Work
##  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
## 
##      0     0  131.98943    0   51   -0.00000  131.98943      -     -    0s
## H    0     0                      81.7644068  131.98943  61.4%     -    0s
## 
## Explored 1 nodes (63 simplex iterations) in 0.03 seconds
## Thread count was 1 (of 8 available processors)
## 
## Solution count 2: 81.7644 -0 
## 
## Optimal solution found (tolerance 1.00e+01)
## Best objective 8.176440680698e+01, best bound 1.319894285520e+02, gap 61.4265%
## Gurobi Optimizer version 9.1.1 build v9.1.1rc0 (mac64)
## Thread count: 4 physical cores, 8 logical processors, using up to 1 threads
## Optimize a model with 63 rows, 1192 columns and 3308 nonzeros
## Model fingerprint: 0x718cd801
## Variable types: 0 continuous, 1192 integer (1192 binary)
## Coefficient statistics:
##   Matrix range     [1e-05, 2e+02]
##   Objective range  [7e-08, 1e+01]
##   Bounds range     [1e+00, 1e+00]
##   RHS range        [9e+03, 9e+03]
## Found heuristic solution: objective -0.0000000
## Presolve removed 0 rows and 143 columns
## Presolve time: 0.02s
## Presolved: 63 rows, 1049 columns, 3152 nonzeros
## Variable types: 0 continuous, 1049 integer (1049 binary)
## Presolved: 63 rows, 1049 columns, 3152 nonzeros
## 
## 
## Root relaxation: objective 1.319894e+02, 70 iterations, 0.00 seconds
## 
##     Nodes    |    Current Node    |     Objective Bounds      |     Work
##  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
## 
##      0     0  131.98943    0   53   -0.00000  131.98943      -     -    0s
## H    0     0                      76.0886867  131.98943  73.5%     -    0s
## 
## Explored 1 nodes (70 simplex iterations) in 0.03 seconds
## Thread count was 1 (of 8 available processors)
## 
## Solution count 2: 76.0887 -0 
## 
## Optimal solution found (tolerance 1.00e+01)
## Best objective 7.608868668868e+01, best bound 1.319894289957e+02, gap 73.4679%
## Gurobi Optimizer version 9.1.1 build v9.1.1rc0 (mac64)
## Thread count: 4 physical cores, 8 logical processors, using up to 1 threads
## Optimize a model with 63 rows, 1192 columns and 3308 nonzeros
## Model fingerprint: 0x91c93c4e
## Variable types: 0 continuous, 1192 integer (1192 binary)
## Coefficient statistics:
##   Matrix range     [1e-05, 2e+02]
##   Objective range  [7e-08, 1e+01]
##   Bounds range     [1e+00, 1e+00]
##   RHS range        [9e+03, 9e+03]
## Found heuristic solution: objective -0.0000000
## Presolve removed 0 rows and 143 columns
## Presolve time: 0.02s
## Presolved: 63 rows, 1049 columns, 3152 nonzeros
## Variable types: 0 continuous, 1049 integer (1049 binary)
## Presolved: 63 rows, 1049 columns, 3152 nonzeros
## 
## 
## Root relaxation: objective 1.319894e+02, 70 iterations, 0.00 seconds
## 
##     Nodes    |    Current Node    |     Objective Bounds      |     Work
##  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
## 
##      0     0  131.98943    0   53   -0.00000  131.98943      -     -    0s
## H    0     0                      76.0886932  131.98943  73.5%     -    0s
## 
## Explored 1 nodes (70 simplex iterations) in 0.03 seconds
## Thread count was 1 (of 8 available processors)
## 
## Solution count 2: 76.0887 -0 
## 
## Optimal solution found (tolerance 1.00e+01)
## Best objective 7.608869321288e+01, best bound 1.319894289957e+02, gap 73.4679%
## Gurobi Optimizer version 9.1.1 build v9.1.1rc0 (mac64)
## Thread count: 4 physical cores, 8 logical processors, using up to 1 threads
## Optimize a model with 63 rows, 1192 columns and 3308 nonzeros
## Model fingerprint: 0x793614d2
## Variable types: 0 continuous, 1192 integer (1192 binary)
## Coefficient statistics:
##   Matrix range     [1e-05, 2e+02]
##   Objective range  [7e-08, 1e+01]
##   Bounds range     [1e+00, 1e+00]
##   RHS range        [9e+03, 9e+03]
## Found heuristic solution: objective -0.0000000
## Presolve removed 0 rows and 143 columns
## Presolve time: 0.02s
## Presolved: 63 rows, 1049 columns, 3152 nonzeros
## Variable types: 0 continuous, 1049 integer (1049 binary)
## Presolved: 63 rows, 1049 columns, 3152 nonzeros
## 
## 
## Root relaxation: objective 1.319894e+02, 63 iterations, 0.00 seconds
## 
##     Nodes    |    Current Node    |     Objective Bounds      |     Work
##  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
## 
##      0     0  131.98943    0   51   -0.00000  131.98943      -     -    0s
## H    0     0                      83.4974777  131.98943  58.1%     -    0s
## 
## Explored 1 nodes (63 simplex iterations) in 0.02 seconds
## Thread count was 1 (of 8 available processors)
## 
## Solution count 2: 83.4975 -0 
## 
## Optimal solution found (tolerance 1.00e+01)
## Best objective 8.349747769044e+01, best bound 1.319894285133e+02, gap 58.0759%
# print solution
print(s8)
## class       : SpatialPolygonsDataFrame 
## features    : 1130 
## extent      : 1080623, 1399989, -4840595, -4497092  (xmin, xmax, ymin, ymax)
## crs         : +proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
## variables   : 11
## names       :   id,              cost, status, locked_in, locked_out, solution_1, solution_2, solution_3, solution_4, solution_5, solution_6 
## min values  :    1, 0.192488262910798,      0,         0,          0,          0,          0,          0,          0,          0,          0 
## max values  : 1130,  61.9272727272727,      2,         1,          1,          1,          1,          1,          1,          1,          1
# calculate the cost of the first solution
sum(s8$solution_1 * s8$cost)
## [1] 680.883
# calculate the cost of the second solution
sum(s8$solution_2 * s8$cost)
## [1] 736.733
# calculate the proportion of planning units with the same solution values
# in the first and second solutions
mean(s8$solution_1 == s8$solution_2)
## [1] 0.9867257
# plot first solution
plot_solution(s8)
# plot all solutions
# s8_plots <- lapply(paste0("solution_", seq_len(6)), function(x) {
#   spplot(s8, x, main = x, col.regions = c("white", "darkgreen"))
# })
# do.call(grid.arrange, append(s8_plots, list(ncol = 3)))

s8$solution_1 <- factor(s8$solution_1)
s8$solution_2 <- factor(s8$solution_2)
s8$solution_3 <- factor(s8$solution_3)
s8$solution_4 <- factor(s8$solution_4)
s8$solution_5 <- factor(s8$solution_5)
s8$solution_6 <- factor(s8$solution_6)

plot(st_as_sf(s8[, c("solution_1", "solution_2", "solution_3", "solution_4", "solution_5", "solution_6")]),
  pal = c("grey90", "darkgreen"))
  1. What is the cost of each of the six solutions in portfolio? Are their costs very different?
  2. Are the solutions in the portfolio very different?
  3. What could we do to obtain a portfolio with more different solutions?