Solve a conservation planning problem().
# S4 method for OptimizationProblem,Solver
solve(a, b, ...)
# S4 method for ConservationProblem,missing
solve(a, b, ..., run_checks = TRUE, force = FALSE)Arguments
| a |
|
|---|---|
| b |
|
| ... | arguments passed to |
| run_checks |
|
| force |
|
Value
A numeric, matrix,
RasterLayer, Spatial,
or sf::sf() object containing the solution to the problem.
Additionally, the returned object will have the following additional
attributes: "objective" containing the solution's objective,
"runtime" denoting the number of seconds that elapsed while solving
the problem, and "status" describing the status of the solution
(e.g. "OPTIMAL" indicates that the optimal solution was found).
In most cases, the first solution (e.g. "solution_001")
will contain the best solution found by the solver (note that this
may not be an optimal solution depending on the gap used to solve
the problem and noting that the default gap is 0.1).
Details
After formulating a conservation planning problem(),
it can be solved using an exact algorithm solver (see solvers
for available solvers). If no solver has been explicitly specified,
then the best available exact algorithm solver will be used by default
(see add_default_solver(). Although these exact algorithm
solvers will often display a lot of information that isn't really that
helpful (e.g. nodes, cutting planes), they do display information
about the progress they are making on solving the problem (e.g. the
performance of the best solution found at a given point in time). If
potential issues were detected during the
presolve checks (see presolve_check())
and the problem is being forcibly solved (i.e. with force = TRUE),
then it is also worth checking for any warnings displayed by the solver
to see if these potential issues are actually causing issues
(e.g. Gurobi can display warnings that include
"Warning: Model contains large matrix coefficient range" and
"Warning: Model contains large rhs").
The object returned from this function depends on the argument to
a. If the argument to a is an
OptimizationProblem object, then the
solution is returned as a logical vector showing the status
of each planning unit in each zone. However, in most cases, the argument
to a is an ConservationProblem object, and so
the type of object returned depends on the number of solutions
generated and the type data used to represent the planning units:
numericvectorcontaining the solution. Here, Each element corresponds to a different planning unit. If multiple solutions are generated, then the solution is returned as alistofnumericvectors.matrixcontaining
numericvalues for the solution. Here, rows correspond to different planning units, and fields (columns) correspond to different management zones. If multiple solutions are generated, then the solution is returned as alistofmatrixobjects.Rasterobject containing the solution in pixel values. If the argument to
xcontains a single management zone, then aRasterLayerobject will be returned. Otherwise, if the argument toxcontains multiple zones, then aRasterStackobject will be returned containing a different layer for each management zone. If multiple solutions are generated, then the solution is returned as alistofRasterobjects.Spatial,sf::sf(), ordata.framecontaining the solution in fields (columns). Here, each row corresponds to a different planning unit. If the argument to
xcontains a single zone, the fields containing solutions are named"solution_XXX"where"XXX"corresponds to the solution number. If the argument toxcontains multiple zones, the fields containing solutions are named"solution_XXX_YYY"where"XXX"corresponds to the solution and"YYY"is the name of the management zone.
After solving problems that contain multiple zones,
it may be useful to use the category_layer() or
category_vector() function to reformat the output.
See also
Examples
# set seed for reproducibility
set.seed(500)
# load data
data(sim_pu_raster, sim_pu_polygons, sim_pu_sf, sim_features,
sim_pu_zones_stack, sim_pu_zones_sf, sim_features_zones)
# build minimal conservation problem with raster data
p1 <- problem(sim_pu_raster, sim_features) %>%
add_min_set_objective() %>%
add_relative_targets(0.1) %>%
add_binary_decisions() %>%
add_default_solver(verbose = FALSE)
# \dontrun{
# solve the problem
s1 <- solve(p1)
# print solution
print(s1)
#> class : RasterLayer
#> dimensions : 10, 10, 100 (nrow, ncol, ncell)
#> resolution : 0.1, 0.1 (x, y)
#> extent : 0, 1, 0, 1 (xmin, xmax, ymin, ymax)
#> crs : NA
#> source : memory
#> names : layer
#> values : 0, 1 (min, max)
#>
# print attributes describing the optimization process and the solution
print(attr(s1, "objective"))
#> solution_1
#> 1987.399
print(attr(s1, "runtime"))
#> solution_1
#> 0.006
print(attr(s1, "status"))
#> solution_1
#> "OPTIMAL"
# calculate feature representation in the solution
r1 <- eval_feature_representation_summary(p1, s1)
print(r1)
#> # A tibble: 5 x 5
#> summary feature total_amount absolute_held relative_held
#> <chr> <chr> <dbl> <dbl> <dbl>
#> 1 overall layer.1 83.3 8.91 0.107
#> 2 overall layer.2 31.2 3.13 0.100
#> 3 overall layer.3 72.0 7.34 0.102
#> 4 overall layer.4 42.7 4.35 0.102
#> 5 overall layer.5 56.7 6.01 0.106
# plot solution
plot(s1, main = "solution", axes = FALSE, box = FALSE)
# }
# build minimal conservation problem with polygon (Spatial) data
p2 <- problem(sim_pu_polygons, sim_features, cost_column = "cost") %>%
add_min_set_objective() %>%
add_relative_targets(0.1) %>%
add_binary_decisions() %>%
add_default_solver(verbose = FALSE)
# \dontrun{
# solve the problem
s2 <- solve(p2)
# print first six rows of the attribute table
print(head(s2))
#> cost locked_in locked_out solution_1
#> 1 215.8638 FALSE FALSE 0
#> 2 212.7823 FALSE FALSE 0
#> 3 207.4962 FALSE FALSE 0
#> 4 208.9322 FALSE TRUE 0
#> 5 214.0419 FALSE FALSE 0
#> 6 213.7636 FALSE FALSE 0
# calculate feature representation in the solution
r2 <- eval_feature_representation_summary(p2, s2[, "solution_1"])
print(r2)
#> # A tibble: 5 x 5
#> summary feature total_amount absolute_held relative_held
#> <chr> <chr> <dbl> <dbl> <dbl>
#> 1 overall layer.1 74.5 8.05 0.108
#> 2 overall layer.2 28.1 2.83 0.101
#> 3 overall layer.3 64.9 6.65 0.103
#> 4 overall layer.4 38.2 3.87 0.101
#> 5 overall layer.5 50.7 5.41 0.107
# plot solution
spplot(s2, zcol = "solution_1", main = "solution", axes = FALSE, box = FALSE)
# }
# build minimal conservation problem with polygon (sf) data
p3 <- problem(sim_pu_sf, sim_features, cost_column = "cost") %>%
add_min_set_objective() %>%
add_relative_targets(0.1) %>%
add_binary_decisions() %>%
add_default_solver(verbose = FALSE)
# \dontrun{
# solve the problem
s3 <- solve(p3)
# print first six rows of the attribute table
print(head(s3))
#> Simple feature collection with 6 features and 4 fields
#> Geometry type: POLYGON
#> Dimension: XY
#> Bounding box: xmin: 0 ymin: 0.9 xmax: 0.6 ymax: 1
#> CRS: NA
#> cost locked_in locked_out solution_1 geometry
#> 1 215.8638 FALSE FALSE 0 POLYGON ((0 1, 0.1 1, 0.1 0...
#> 2 212.7823 FALSE FALSE 0 POLYGON ((0.1 1, 0.2 1, 0.2...
#> 3 207.4962 FALSE FALSE 0 POLYGON ((0.2 1, 0.3 1, 0.3...
#> 4 208.9322 FALSE TRUE 0 POLYGON ((0.3 1, 0.4 1, 0.4...
#> 5 214.0419 FALSE FALSE 0 POLYGON ((0.4 1, 0.5 1, 0.5...
#> 6 213.7636 FALSE FALSE 0 POLYGON ((0.5 1, 0.6 1, 0.6...
# calculate feature representation in the solution
r3 <- eval_feature_representation_summary(p3, s3[, "solution_1"])
print(r3)
#> # A tibble: 5 x 5
#> summary feature total_amount absolute_held relative_held
#> <chr> <chr> <dbl> <dbl> <dbl>
#> 1 overall layer.1 74.5 8.05 0.108
#> 2 overall layer.2 28.1 2.83 0.101
#> 3 overall layer.3 64.9 6.65 0.103
#> 4 overall layer.4 38.2 3.87 0.101
#> 5 overall layer.5 50.7 5.41 0.107
# plot solution
plot(s3[, "solution_1"])
# }
# build multi-zone conservation problem with raster data
p4 <- problem(sim_pu_zones_stack, sim_features_zones) %>%
add_min_set_objective() %>%
add_relative_targets(matrix(runif(15, 0.1, 0.2), nrow = 5,
ncol = 3)) %>%
add_binary_decisions() %>%
add_default_solver(verbose = FALSE)
# \dontrun{
# solve the problem
s4 <- solve(p4)
# print solution
print(s4)
#> class : RasterStack
#> dimensions : 10, 10, 100, 3 (nrow, ncol, ncell, nlayers)
#> resolution : 0.1, 0.1 (x, y)
#> extent : 0, 1, 0, 1 (xmin, xmax, ymin, ymax)
#> crs : NA
#> names : zone_1, zone_2, zone_3
#> min values : 0, 0, 0
#> max values : 1, 1, 1
#>
# calculate feature representation in the solution
r4 <- eval_feature_representation_summary(p4, s4)
print(r4)
#> # A tibble: 20 x 5
#> summary feature total_amount absolute_held relative_held
#> <chr> <chr> <dbl> <dbl> <dbl>
#> 1 overall feature_1 250. 46.2 0.185
#> 2 overall feature_2 93.6 17.5 0.187
#> 3 overall feature_3 216. 41.6 0.193
#> 4 overall feature_4 128. 22.1 0.172
#> 5 overall feature_5 170. 30.5 0.179
#> 6 zone_1 feature_1 83.3 16.3 0.195
#> 7 zone_1 feature_2 31.2 5.40 0.173
#> 8 zone_1 feature_3 72.0 14.2 0.198
#> 9 zone_1 feature_4 42.7 6.55 0.154
#> 10 zone_1 feature_5 56.7 10.7 0.189
#> 11 zone_2 feature_1 83.3 15.5 0.186
#> 12 zone_2 feature_2 31.2 6.14 0.197
#> 13 zone_2 feature_3 72.0 14.6 0.203
#> 14 zone_2 feature_4 42.7 7.83 0.184
#> 15 zone_2 feature_5 56.7 9.90 0.175
#> 16 zone_3 feature_1 83.3 14.4 0.173
#> 17 zone_3 feature_2 31.2 5.96 0.191
#> 18 zone_3 feature_3 72.0 12.7 0.177
#> 19 zone_3 feature_4 42.7 7.68 0.180
#> 20 zone_3 feature_5 56.7 9.83 0.173
# plot solution
plot(category_layer(s4), main = "solution", axes = FALSE, box = FALSE)
# }
# build multi-zone conservation problem with polygon (sf) data
p5 <- problem(sim_pu_zones_sf, sim_features_zones,
cost_column = c("cost_1", "cost_2", "cost_3")) %>%
add_min_set_objective() %>%
add_relative_targets(matrix(runif(15, 0.1, 0.2), nrow = 5,
ncol = 3)) %>%
add_binary_decisions() %>%
add_default_solver(verbose = FALSE)
# \dontrun{
# solve the problem
s5 <- solve(p5)
# print first six rows of the attribute table
print(head(s5))
#> Simple feature collection with 6 features and 9 fields
#> Geometry type: POLYGON
#> Dimension: XY
#> Bounding box: xmin: 0 ymin: 0.9 xmax: 0.6 ymax: 1
#> CRS: NA
#> cost_1 cost_2 cost_3 locked_1 locked_2 locked_3 solution_1_zone_1
#> 1 215.8638 183.3344 205.4113 FALSE FALSE FALSE 0
#> 2 212.7823 189.4978 209.6404 FALSE FALSE FALSE 0
#> 3 207.4962 193.6007 215.4212 TRUE FALSE FALSE 0
#> 4 208.9322 197.5897 218.5241 FALSE FALSE FALSE 0
#> 5 214.0419 199.8033 220.7100 FALSE FALSE FALSE 0
#> 6 213.7636 203.1867 224.6809 FALSE FALSE FALSE 0
#> solution_1_zone_2 solution_1_zone_3 geometry
#> 1 0 0 POLYGON ((0 1, 0.1 1, 0.1 0...
#> 2 1 0 POLYGON ((0.1 1, 0.2 1, 0.2...
#> 3 1 0 POLYGON ((0.2 1, 0.3 1, 0.3...
#> 4 1 0 POLYGON ((0.3 1, 0.4 1, 0.4...
#> 5 1 0 POLYGON ((0.4 1, 0.5 1, 0.5...
#> 6 1 0 POLYGON ((0.5 1, 0.6 1, 0.6...
# calculate feature representation in the solution
r5 <- eval_feature_representation_summary(
p5, s5[, c("solution_1_zone_1", "solution_1_zone_2", "solution_1_zone_3")])
print(r5)
#> # A tibble: 20 x 5
#> summary feature total_amount absolute_held relative_held
#> <chr> <chr> <dbl> <dbl> <dbl>
#> 1 overall feature_1 225. 36.7 0.163
#> 2 overall feature_2 83.9 14.2 0.169
#> 3 overall feature_3 195. 32.9 0.169
#> 4 overall feature_4 114. 18.4 0.161
#> 5 overall feature_5 154. 24.3 0.159
#> 6 zone_1 feature_1 75.1 12.4 0.166
#> 7 zone_1 feature_2 28.0 4.35 0.155
#> 8 zone_1 feature_3 65.0 10.5 0.162
#> 9 zone_1 feature_4 38.0 6.13 0.161
#> 10 zone_1 feature_5 51.2 8.32 0.163
#> 11 zone_2 feature_1 75.1 10.2 0.136
#> 12 zone_2 feature_2 28.0 5.07 0.181
#> 13 zone_2 feature_3 65.0 9.83 0.151
#> 14 zone_2 feature_4 38.0 6.33 0.167
#> 15 zone_2 feature_5 51.2 6.79 0.133
#> 16 zone_3 feature_1 75.1 14.0 0.187
#> 17 zone_3 feature_2 28.0 4.76 0.170
#> 18 zone_3 feature_3 65.0 12.6 0.193
#> 19 zone_3 feature_4 38.0 5.96 0.157
#> 20 zone_3 feature_5 51.2 9.23 0.180
# create new column representing the zone id that each planning unit
# was allocated to in the solution
s5$solution <- category_vector(s5[, c("solution_1_zone_1",
"solution_1_zone_2",
"solution_1_zone_3")])
s5$solution <- factor(s5$solution)
# plot solution
plot(s5[, "solution"])
# }