Solve

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	`ConservationProblem-class` or an `OptimizationProblem-class` object.
b	`Solver-class` object. Not used if `a` is an `ConservationProblem-class` object.
...	arguments passed to `compile`.
run_checks	`logical` flag indicating whether presolve checks should be run prior solving the problem. These checks are performed using the `presolve_check` function. Defaults to `TRUE`. Skipping these checks may reduce run time for large problems.
force	`logical` flag indicating if an attempt to should be made to solve the problem even if potential issues were detected during the presolve checks. Defaults to `FALSE`.

Value

A numeric, matrix, RasterLayer-class, Spatial-class, or 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).

Details

The object returned from this function depends on the argument to a. If the argument to a is an OptimizationProblem-class object, then the solution is returned as a logical vector showing the status of each planning unit in each zone. On the other hand, if the argument to a is an ConservationProblem-class object, then the type of object returned depends on the number of solutions generated and the type data used to represent planning unit costs in the argument to a.

numeric

vector containing the solution. Here, Each element corresponds to a different planning unit. If multiple solutions are generated, then the solution is returned as a list of numeric vectors.

matrix

containing numeric values 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 a list of matrix objects.

Raster-class

object containing the solution in pixel values. If the argument to x contains a single management zone, then a RasterLayer object will be returned. Otherwise, if the argument to x contains multiple zones, then a RasterStack-class object will be returned containing a different layer for each management zone. If multiple solutions are generated, then the solution is returned as a list of Raster objects.

Spatial-class, sf, or data.frame

containing the solution in fields (columns). Here, each row corresponds to a different planning unit. If the argument to x contains a single zone, the fields containing solutions are named "solution_XXX" where "XXX" corresponds to the solution number. If the argument to x contains 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.

Since this function returns an object that specifies how much of each planning unit is allocated to each management zone, it may be useful to use the category_layer function to reformat the output for problems containing multiple zones.

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()
# \donttest{
# solve the problem
s1 <- solve(p1)
#> Gurobi Optimizer version 9.0.1 build v9.0.1rc0 (linux64)
#> Optimize a model with 5 rows, 90 columns and 450 nonzeros
#> Model fingerprint: 0x6442bf6e
#> Variable types: 0 continuous, 90 integer (90 binary)
#> Coefficient statistics:
#>   Matrix range     [2e-01, 9e-01]
#>   Objective range  [2e+02, 2e+02]
#>   Bounds range     [1e+00, 1e+00]
#>   RHS range        [3e+00, 8e+00]
#> Found heuristic solution: objective 2337.9617505
#> Presolve time: 0.00s
#> Presolved: 5 rows, 90 columns, 450 nonzeros
#> Variable types: 0 continuous, 90 integer (90 binary)
#> Presolved: 5 rows, 90 columns, 450 nonzeros
#> 
#> 
#> Root relaxation: objective 1.931582e+03, 12 iterations, 0.00 seconds
#> 
#>     Nodes    |    Current Node    |     Objective Bounds      |     Work
#>  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
#> 
#>      0     0 1931.58191    0    4 2337.96175 1931.58191  17.4%     -    0s
#> H    0     0                    1987.3985265 1931.58191  2.81%     -    0s
#> 
#> Explored 1 nodes (12 simplex iterations) in 0.00 seconds
#> Thread count was 1 (of 4 available processors)
#> 
#> Solution count 2: 1987.4 2337.96 
#> 
#> Optimal solution found (tolerance 1.00e-01)
#> Best objective 1.987398526526e+03, best bound 1.931581908865e+03, gap 2.8085%

# 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.001976967 
print(attr(s1, "status"))
#> solution_1 
#>  "OPTIMAL" 

# calculate feature representation in the solution
r1 <- feature_representation(p1, s1)
print(r1)
#> # A tibble: 5 x 3
#>   feature absolute_held relative_held
#>   <chr>           <dbl>         <dbl>
#> 1 layer.1          8.91         0.107
#> 2 layer.2          3.13         0.100
#> 3 layer.3          7.34         0.102
#> 4 layer.4          4.35         0.102
#> 5 layer.5          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()
# \donttest{
# solve the problem
s2 <- solve(p2)
#> Gurobi Optimizer version 9.0.1 build v9.0.1rc0 (linux64)
#> Optimize a model with 5 rows, 90 columns and 450 nonzeros
#> Model fingerprint: 0xa41351c6
#> Variable types: 0 continuous, 90 integer (90 binary)
#> Coefficient statistics:
#>   Matrix range     [2e-01, 9e-01]
#>   Objective range  [2e+02, 2e+02]
#>   Bounds range     [1e+00, 1e+00]
#>   RHS range        [3e+00, 7e+00]
#> Found heuristic solution: objective 2145.2678910
#> Presolve time: 0.00s
#> Presolved: 5 rows, 90 columns, 450 nonzeros
#> Variable types: 0 continuous, 90 integer (90 binary)
#> Presolved: 5 rows, 90 columns, 450 nonzeros
#> 
#> 
#> Root relaxation: objective 1.732564e+03, 12 iterations, 0.00 seconds
#> 
#>     Nodes    |    Current Node    |     Objective Bounds      |     Work
#>  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
#> 
#>      0     0 1732.56396    0    4 2145.26789 1732.56396  19.2%     -    0s
#> H    0     0                    1792.5345202 1732.56396  3.35%     -    0s
#> 
#> Explored 1 nodes (12 simplex iterations) in 0.00 seconds
#> Thread count was 1 (of 4 available processors)
#> 
#> Solution count 2: 1792.53 2145.27 
#> 
#> Optimal solution found (tolerance 1.00e-01)
#> Best objective 1.792534520218e+03, best bound 1.732563961206e+03, gap 3.3456%

# 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 <- feature_representation(p2, s2[, "solution_1"])
print(r2)
#> # A tibble: 5 x 3
#>   feature absolute_held relative_held
#>   <chr>           <dbl>         <dbl>
#> 1 layer.1          8.05         0.108
#> 2 layer.2          2.83         0.101
#> 3 layer.3          6.65         0.103
#> 4 layer.4          3.87         0.101
#> 5 layer.5          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()
# \donttest{
# solve the problem
s3 <- solve(p3)
#> Gurobi Optimizer version 9.0.1 build v9.0.1rc0 (linux64)
#> Optimize a model with 5 rows, 90 columns and 450 nonzeros
#> Model fingerprint: 0xa41351c6
#> Variable types: 0 continuous, 90 integer (90 binary)
#> Coefficient statistics:
#>   Matrix range     [2e-01, 9e-01]
#>   Objective range  [2e+02, 2e+02]
#>   Bounds range     [1e+00, 1e+00]
#>   RHS range        [3e+00, 7e+00]
#> Found heuristic solution: objective 2145.2678910
#> Presolve time: 0.00s
#> Presolved: 5 rows, 90 columns, 450 nonzeros
#> Variable types: 0 continuous, 90 integer (90 binary)
#> Presolved: 5 rows, 90 columns, 450 nonzeros
#> 
#> 
#> Root relaxation: objective 1.732564e+03, 12 iterations, 0.00 seconds
#> 
#>     Nodes    |    Current Node    |     Objective Bounds      |     Work
#>  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
#> 
#>      0     0 1732.56396    0    4 2145.26789 1732.56396  19.2%     -    0s
#> H    0     0                    1792.5345202 1732.56396  3.35%     -    0s
#> 
#> Explored 1 nodes (12 simplex iterations) in 0.00 seconds
#> Thread count was 1 (of 4 available processors)
#> 
#> Solution count 2: 1792.53 2145.27 
#> 
#> Optimal solution found (tolerance 1.00e-01)
#> Best objective 1.792534520218e+03, best bound 1.732563961206e+03, gap 3.3456%

# 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
#> bbox:           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 <- feature_representation(p3, s3[, "solution_1"])
print(r3)
#> # A tibble: 5 x 3
#>   feature absolute_held relative_held
#>   <chr>           <dbl>         <dbl>
#> 1 layer.1          8.05         0.108
#> 2 layer.2          2.83         0.101
#> 3 layer.3          6.65         0.103
#> 4 layer.4          3.87         0.101
#> 5 layer.5          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()
# \donttest{
# solve the problem
s4 <- solve(p4)
#> Gurobi Optimizer version 9.0.1 build v9.0.1rc0 (linux64)
#> Optimize a model with 105 rows, 270 columns and 1620 nonzeros
#> Model fingerprint: 0xa478f36b
#> Variable types: 0 continuous, 270 integer (270 binary)
#> Coefficient statistics:
#>   Matrix range     [2e-01, 1e+00]
#>   Objective range  [2e+02, 2e+02]
#>   Bounds range     [1e+00, 1e+00]
#>   RHS range        [1e+00, 2e+01]
#> Found heuristic solution: objective 12614.662551
#> Presolve removed 3 rows and 0 columns
#> Presolve time: 0.00s
#> Presolved: 102 rows, 270 columns, 1350 nonzeros
#> Variable types: 0 continuous, 270 integer (270 binary)
#> Presolved: 102 rows, 270 columns, 1350 nonzeros
#> 
#> 
#> Root relaxation: objective 1.096484e+04, 94 iterations, 0.00 seconds
#> 
#>     Nodes    |    Current Node    |     Objective Bounds      |     Work
#>  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
#> 
#>      0     0 10964.8406    0   14 12614.6626 10964.8406  13.1%     -    0s
#> H    0     0                    11209.793128 10964.8406  2.19%     -    0s
#> 
#> Explored 1 nodes (94 simplex iterations) in 0.01 seconds
#> Thread count was 1 (of 4 available processors)
#> 
#> Solution count 2: 11209.8 12614.7 
#> 
#> Optimal solution found (tolerance 1.00e-01)
#> Best objective 1.120979312806e+04, best bound 1.096484061295e+04, gap 2.1852%

# 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      : layer.1.1, layer.1.2, layer.1.3 
#> min values :         0,         0,         0 
#> max values :         1,         1,         1 
#> 

# calculate feature representation in the solution
r4 <- feature_representation(p4, s4)
print(r4)
#> # A tibble: 15 x 4
#>    feature   zone   absolute_held relative_held
#>    <chr>     <chr>          <dbl>         <dbl>
#>  1 feature_1 zone_1         16.3          0.195
#>  2 feature_2 zone_1          5.40         0.173
#>  3 feature_3 zone_1         14.2          0.198
#>  4 feature_4 zone_1          6.55         0.154
#>  5 feature_5 zone_1         10.7          0.189
#>  6 feature_1 zone_2         15.5          0.186
#>  7 feature_2 zone_2          6.14         0.197
#>  8 feature_3 zone_2         14.6          0.203
#>  9 feature_4 zone_2          7.83         0.184
#> 10 feature_5 zone_2          9.90         0.175
#> 11 feature_1 zone_3         14.4          0.173
#> 12 feature_2 zone_3          5.96         0.191
#> 13 feature_3 zone_3         12.7          0.177
#> 14 feature_4 zone_3          7.68         0.180
#> 15 feature_5 zone_3          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()
# \donttest{
# solve the problem
s5 <- solve(p5)
#> Gurobi Optimizer version 9.0.1 build v9.0.1rc0 (linux64)
#> Optimize a model with 105 rows, 270 columns and 1620 nonzeros
#> Model fingerprint: 0x00214866
#> Variable types: 0 continuous, 270 integer (270 binary)
#> Coefficient statistics:
#>   Matrix range     [2e-01, 1e+00]
#>   Objective range  [2e+02, 2e+02]
#>   Bounds range     [1e+00, 1e+00]
#>   RHS range        [1e+00, 1e+01]
#> Found heuristic solution: objective 10229.621561
#> Presolve removed 6 rows and 0 columns
#> Presolve time: 0.00s
#> Presolved: 99 rows, 270 columns, 1080 nonzeros
#> Variable types: 0 continuous, 270 integer (270 binary)
#> Presolved: 99 rows, 270 columns, 1080 nonzeros
#> 
#> 
#> Root relaxation: objective 8.838392e+03, 43 iterations, 0.00 seconds
#> 
#>     Nodes    |    Current Node    |     Objective Bounds      |     Work
#>  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
#> 
#>      0     0 8838.39193    0    8 10229.6216 8838.39193  13.6%     -    0s
#> H    0     0                    9015.0191263 8838.39193  1.96%     -    0s
#> 
#> Explored 1 nodes (43 simplex iterations) in 0.01 seconds
#> Thread count was 1 (of 4 available processors)
#> 
#> Solution count 2: 9015.02 10229.6 
#> 
#> Optimal solution found (tolerance 1.00e-01)
#> Best objective 9.015019126339e+03, best bound 8.838391930390e+03, gap 1.9593%

# 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
#> bbox:           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 <- feature_representation(p5, s5[, c("solution_1_zone_1",
                                        "solution_1_zone_2",
                                        "solution_1_zone_3")])
print(r5)
#> # A tibble: 15 x 4
#>    feature   zone   absolute_held relative_held
#>    <chr>     <chr>          <dbl>         <dbl>
#>  1 feature_1 zone_1         12.4          0.166
#>  2 feature_2 zone_1          4.35         0.155
#>  3 feature_3 zone_1         10.5          0.162
#>  4 feature_4 zone_1          6.13         0.161
#>  5 feature_5 zone_1          8.32         0.163
#>  6 feature_1 zone_2         10.2          0.136
#>  7 feature_2 zone_2          5.07         0.181
#>  8 feature_3 zone_2          9.83         0.151
#>  9 feature_4 zone_2          6.33         0.167
#> 10 feature_5 zone_2          6.79         0.133
#> 11 feature_1 zone_3         14.0          0.187
#> 12 feature_2 zone_3          4.76         0.170
#> 13 feature_3 zone_3         12.6          0.193
#> 14 feature_4 zone_3          5.96         0.157
#> 15 feature_5 zone_3          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"])
# }

Arguments

Value

Details

See also

Examples