if (FALSE) {
# 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)
}
Run the code above in your browser using DataLab