Set the objective of a conservation planning problem() to maximize the phylogenetic diversity of the features represented in the solution subject to a budget. This objective is similar to add_max_features_objective() except that emphasis is placed on representing a phylogenetically diverse set of species, rather than as many features as possible (subject to weights). This function was inspired by Faith (1992) and Rodrigues et al. (2002).

add_max_phylo_div_objective(x, budget, tree)

Arguments

x problem() (i.e. ConservationProblem) object. numeric value specifying the maximum expenditure of the prioritization. For problems with multiple zones, the argument to budget can be a single numeric value to specify a budget for the entire solution or a numeric vector to specify a budget for each each management zone. phylo() object specifying a phylogenetic tree for the conservation features.

Value

Object (i.e. ConservationProblem) with the objective added to it.

Details

A problem objective is used to specify the overall goal of the conservation planning problem. Please note that all conservation planning problems formulated in the prioritizr package require the addition of objectives---failing to do so will return an error message when attempting to solve problem.

The maximum phylogenetic diversity objective finds the set of planning units that meets representation targets for a phylogenetic tree while staying within a fixed budget. If multiple solutions can meet all targets while staying within budget, the cheapest solution is chosen. Note that this objective is similar to the maximum features objective (add_max_features_objective()) in that it allows for both a budget and targets to be set for each feature. However, unlike the maximum feature objective, the aim of this objective is to maximize the total phylogenetic diversity of the targets met in the solution, so if multiple targets are provided for a single feature, the problem will only need to meet a single target for that feature for the phylogenetic benefit for that feature to be counted when calculating the phylogenetic diversity of the solution. In other words, for multi-zone problems, this objective does not aim to maximize the phylogenetic diversity in each zone, but rather this objective aims to maximize the phylogenetic diversity of targets that can be met through allocating planning units to any of the different zones in a problem. This can be useful for problems where targets pertain to the total amount held for each feature across multiple zones. For example, each feature might have a non-zero amount of suitable habitat in each planning unit when the planning units are assigned to a (i) not restored, (ii) partially restored, or (iii) completely restored management zone. Here each target corresponds to a single feature and can be met through the total amount of habitat in planning units present to the three zones.

The maximum phylogenetic diversity objective for the reserve design problem can be expressed mathematically for a set of planning units ($$I$$ indexed by $$i$$) and a set of features ($$J$$ indexed by $$j$$) as:

$$\mathit{Maximize} \space \sum_{i = 1}^{I} -s \space c_i \space x_i + \sum_{j = 1}^{J} m_b l_b \\ \mathit{subject \space to} \\ \sum_{i = 1}^{I} x_i r_{ij} \geq y_j t_j \forall j \in J \\ m_b \leq y_j \forall j \in T(b) \\ \sum_{i = 1}^{I} x_i c_i \leq B$$

Here, $$x_i$$ is the decisions variable (e.g. specifying whether planning unit $$i$$ has been selected (1) or not (0)), $$r_{ij}$$ is the amount of feature $$j$$ in planning unit $$i$$, $$t_j$$ is the representation target for feature $$j$$, $$y_j$$ indicates if the solution has meet the target $$t_j$$ for feature $$j$$. Additionally, $$T$$ represents a phylogenetic tree containing features $$j$$ and has the branches $$b$$ associated within lengths $$l_b$$. The binary variable $$m_b$$ denotes if at least one feature associated with the branch $$b$$ has met its representation as indicated by $$y_j$$. For brevity, we denote the features $$j$$ associated with branch $$b$$ using $$T(b)$$. Finally, $$B$$ is the budget allocated for the solution, $$c_i$$ is the cost of planning unit $$i$$, and $$s$$ is a scaling factor used to shrink the costs so that the problem will return a cheapest solution when there are multiple solutions that represent the same amount of all features within the budget.

Notes

In early versions, this function was named as the add_max_phylo_div_objective function.

References

Faith DP (1992) Conservation evaluation and phylogenetic diversity. Biological Conservation, 61: 1--10.

Rodrigues ASL and Gaston KJ (2002) Maximising phylogenetic diversity in the selection of networks of conservation areas. Biological Conservation, 105: 103--111.

objectives, branch_matrix().

Examples

# load ape package
require(ape)

data(sim_pu_raster, sim_features, sim_phylogeny, sim_pu_zones_stack,
sim_features_zones)

# plot the simulated phylogeny
# \dontrun{
par(mfrow = c(1, 1))
plot(sim_phylogeny, main = "phylogeny")# }
# create problem with a maximum phylogenetic diversity objective,
# where each feature needs 10 % of its distribution to be secured for
# it to be adequately conserved and a total budget of 1900
p1 <- problem(sim_pu_raster, sim_features) %>%
# \dontrun{
# solve problem
s1 <- solve(p1)#> Gurobi Optimizer version 9.0.2 build v9.0.2rc0 (linux64)
#> Optimize a model with 14 rows, 103 columns and 565 nonzeros
#> Model fingerprint: 0x19bea6ab
#> Variable types: 0 continuous, 103 integer (103 binary)
#> Coefficient statistics:
#>   Matrix range     [2e-01, 2e+02]
#>   Objective range  [4e-06, 1e+00]
#>   Bounds range     [1e+00, 1e+00]
#>   RHS range        [2e+03, 2e+03]
#> Found heuristic solution: objective -0.0000000
#> Presolve removed 5 rows and 5 columns
#> Presolve time: 0.00s
#> Presolved: 9 rows, 98 columns, 555 nonzeros
#> Variable types: 0 continuous, 98 integer (98 binary)
#> Presolved: 9 rows, 98 columns, 555 nonzeros
#>
#>
#> Root relaxation: objective 5.344516e+00, 23 iterations, 0.00 seconds
#>
#>     Nodes    |    Current Node    |     Objective Bounds      |     Work
#>  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
#>
#>      0     0    5.34452    0    5   -0.00000    5.34452      -     -    0s
#>      0     0    5.24678    0    6   -0.00000    5.24678      -     -    0s
#> H    0     0                       2.9402664    5.24678  78.4%     -    0s
#>      0     0    5.24429    0    7    2.94027    5.24429  78.4%     -    0s
#>      0     0    5.24428    0    8    2.94027    5.24428  78.4%     -    0s
#>      0     0    5.24405    0    8    2.94027    5.24405  78.4%     -    0s
#>      0     0    5.24401    0    9    2.94027    5.24401  78.4%     -    0s
#>      0     0    5.24399    0   10    2.94027    5.24399  78.4%     -    0s
#>      0     0    5.24397    0   10    2.94027    5.24397  78.4%     -    0s
#>      0     2    5.24023    0   10    2.94027    5.24023  78.2%     -    0s
#> H    6     3                       3.5363171    4.18960  18.5%   6.8    0s
#>
#> Cutting planes:
#>   Cover: 1
#>   MIR: 3
#>   RLT: 1
#>
#> Explored 9 nodes (111 simplex iterations) in 0.02 seconds
#> Thread count was 1 (of 4 available processors)
#>
#> Solution count 3: 3.53632 2.94027 -0
#>
#> Optimal solution found (tolerance 1.00e-01)
#> Best objective 3.536317127664e+00, best bound 3.833000595646e+00, gap 8.3896%
# plot solution
plot(s1, main = "solution", axes = FALSE, box = FALSE)
# find which features have their targets met
r1 <- feature_representation(p1, s1)
r1$target_met <- r1$relative_held > 0.1
print(r1)#> # A tibble: 5 x 4
#>   feature absolute_held relative_held target_met
#>   <chr>           <dbl>         <dbl> <lgl>
#> 1 layer.1          5.87        0.0705 FALSE
#> 2 layer.2          3.18        0.102  TRUE
#> 3 layer.3          4.39        0.0611 FALSE
#> 4 layer.4          4.37        0.102  TRUE
#> 5 layer.5          4.43        0.0781 FALSE
# plot the phylogeny and color the adequately represented features in red
plot(sim_phylogeny, main = "adequately represented features",
tip.color = replace(
rep("black", nlayers(sim_features)),
sim_phylogeny$tip.label %in% r1$feature[r1$target_met], "red"))# } # rename the features in the example phylogeny for use with the # multi-zone data sim_phylogeny$tip.label <- feature_names(sim_features_zones)

# create targets for a multi-zone problem. Here, each feature needs a total
# of 10 units of habitat to be conserved among the three zones to be
targets <- tibble::tibble(
feature = feature_names(sim_features_zones),
zone = list(zone_names(sim_features_zones))[rep(1,
number_of_features(sim_features_zones))],
type = rep("absolute", number_of_features(sim_features_zones)),
target = rep(10, number_of_features(sim_features_zones)))

# create a multi-zone problem with a maximum phylogenetic diversity
# objective, where the total expenditure in all zones is 5000.
p2 <- problem(sim_pu_zones_stack, sim_features_zones) %>%
# \dontrun{
# solve problem
s2 <- solve(p2)#> Gurobi Optimizer version 9.0.2 build v9.0.2rc0 (linux64)
#> Optimize a model with 104 rows, 283 columns and 1915 nonzeros
#> Model fingerprint: 0x20546489
#> Variable types: 0 continuous, 283 integer (283 binary)
#> Coefficient statistics:
#>   Matrix range     [2e-01, 2e+02]
#>   Objective range  [1e-06, 1e+00]
#>   Bounds range     [1e+00, 1e+00]
#>   RHS range        [1e+00, 5e+03]
#> Found heuristic solution: objective -0.0000000
#> Presolve removed 97 rows and 187 columns
#> Presolve time: 0.00s
#> Presolved: 7 rows, 96 columns, 460 nonzeros
#> Variable types: 0 continuous, 96 integer (96 binary)
#> Presolved: 7 rows, 96 columns, 460 nonzeros
#>
#>
#> Root relaxation: objective 5.055579e+00, 13 iterations, 0.00 seconds
#>
#>     Nodes    |    Current Node    |     Objective Bounds      |     Work
#>  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
#>
#>      0     0    5.05558    0    1   -0.00000    5.05558      -     -    0s
#> H    0     0                       3.8834554    5.05558  30.2%     -    0s
#> H    0     0                       5.0555775    5.05558  0.00%     -    0s
#>
#> Explored 1 nodes (13 simplex iterations) in 0.00 seconds
#> Thread count was 1 (of 4 available processors)
#>
#> Solution count 3: 5.05558 3.88346 -0
#>
#> Optimal solution found (tolerance 1.00e-01)
#> Best objective 5.055577499895e+00, best bound 5.055578671837e+00, gap 0.0000%
# plot solution
plot(category_layer(s2), main = "solution", axes = FALSE, box = FALSE)
# calculate total amount of habitat conserved for each feature among
# all three management zones
amount_held2 <- numeric(number_of_features(sim_features_zones))
for (z in seq_len(number_of_zones(sim_features_zones)))
amount_held2 <- amount_held2 +
cellStats(sim_features_zones[[z]] * s2[[z]], "sum")

# find which features have their targets met
targets_met2 <- amount_held2 >= targets$target print(targets_met2)#> layer.1 layer.2 layer.3 layer.4 layer.5 #> TRUE FALSE TRUE TRUE TRUE # plot the phylogeny and color the adequately represented features in red plot(sim_phylogeny, main = "adequately represented features", tip.color = replace(rep("black", nlayers(sim_features)), which(targets_met2), "red"))# } # create a multi-zone problem with a maximum phylogenetic diversity # objective, where each zone has a separate budget. p3 <- problem(sim_pu_zones_stack, sim_features_zones) %>% add_max_phylo_div_objective(c(2500, 500, 2000), sim_phylogeny) %>% add_manual_targets(targets) %>% add_binary_decisions() # \dontrun{ # solve problem s3 <- solve(p3)#> Gurobi Optimizer version 9.0.2 build v9.0.2rc0 (linux64) #> Optimize a model with 106 rows, 283 columns and 1915 nonzeros #> Model fingerprint: 0xa9bb83d0 #> Variable types: 0 continuous, 283 integer (283 binary) #> Coefficient statistics: #> Matrix range [2e-01, 2e+02] #> Objective range [1e-06, 1e+00] #> Bounds range [1e+00, 1e+00] #> RHS range [1e+00, 2e+03] #> Found heuristic solution: objective -0.0000000 #> Presolve removed 5 rows and 5 columns #> Presolve time: 0.00s #> Presolved: 101 rows, 278 columns, 1905 nonzeros #> Variable types: 0 continuous, 278 integer (278 binary) #> Presolved: 101 rows, 278 columns, 1905 nonzeros #> #> #> Root relaxation: objective 5.319634e+00, 88 iterations, 0.00 seconds #> #> Nodes | Current Node | Objective Bounds | Work #> Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time #> #> 0 0 5.31963 0 6 -0.00000 5.31963 - - 0s #> H 0 0 5.0555789 5.31963 5.22% - 0s #> #> Explored 1 nodes (88 simplex iterations) in 0.01 seconds #> Thread count was 1 (of 4 available processors) #> #> Solution count 2: 5.05558 -0 #> #> Optimal solution found (tolerance 1.00e-01) #> Best objective 5.055578907883e+00, best bound 5.319633591566e+00, gap 5.2230% # plot solution plot(category_layer(s3), main = "solution", axes = FALSE, box = FALSE) # calculate total amount of habitat conserved for each feature among # all three management zones amount_held3 <- numeric(number_of_features(sim_features_zones)) for (z in seq_len(number_of_zones(sim_features_zones))) amount_held3 <- amount_held3 + cellStats(sim_features_zones[[z]] * s3[[z]], "sum") # find which features have their targets met targets_met3 <- amount_held3 >= targets$target
print(targets_met3)#> layer.1 layer.2 layer.3 layer.4 layer.5
#>    TRUE   FALSE    TRUE    TRUE    TRUE
# plot the phylogeny and color the adequately represented features in red
plot(sim_phylogeny, main = "adequately represented features",
tip.color = replace(rep("black", nlayers(sim_features)),
which(targets_met3), "red"))# }