The `prioriactions`

package allows to address to planning
goals: **recovery** and **conservation**. In
order to understand the difference between each of them, let us consider
the following figures:

According to the figure, for a conservation feature there are
**two** planning units where it co-occurs with the threat
it is sensitive to (and where, therefore, it can be impacted by the
threat) and in **four** other units where the feature does
not co-occur with the threat. Therefore, **recovery
targets** could be achieved by prescribing management actions
against the threat present in the two planning units where the feature
co-occurs with that threat, while **conservation targets**
could be achieved by selecting for conservation any of the planning
units where the species does not co-occur with the threat. As we will
show in this document, selecting planning units only for conservation
purposes plays a crucial role in the conservation plan.

If we aim at minimizing the conservation costs
(`minimize costs`

), the `prioriactions`

package
allows to explicitly impose a target for both, the
**recovery** and the **conservation** levels
that we aim at achieving. Please note that, by the default, only the
**recovery targets** is requested, given the main focus of
`prioriactions`

is on the prioritization of management
actions against threats, while conservation targets are optional. There
are other software, like marxan, and packages, like prioritizr, that already
solve conservation only problems. On the contrary, if we aim at
maximizing the benefits (`maximize benefits`

), we assume that
the benefits are obtained only from **recovery**
actions.

The remainder of this document is organized as follows. In Section 1
we present the results obtained when only imposing **recovery
targets**. In Section 2 we present the results obtained when
imposing **recovery targets** and **conservation
targets**. In Section 3 we present the results obtained when
considering **recovery targets** along with a connectivity
optimization. Finally, in Section 4 we present the results obtained when
considering **recovery targets** and **conservation
targets** along with a connectivity optimization.

In order to analyze how planning decision behave when no
**conservation target** is considered, we present some
results using the example available in the Get
started vignette.

```
# load package
library(prioriactions)
data(sim_boundary_data, sim_dist_features_data, sim_dist_threats_data,
sim_features_data, sim_pu_data, sim_sensitivity_data, sim_threats_data)
#create raster
library(raster)
r <- raster(ncol=10, nrow=10, xmn=0, xmx=10, ymn=0, ymx=10)
values(r) <- 0
# create data class using input data from prioriactions
b <- inputData(pu = sim_pu_data,
features = sim_features_data,
dist_features = sim_dist_features_data,
threats = sim_threats_data,
dist_threats = sim_dist_threats_data,
sensitivity = sim_sensitivity_data,
boundary = sim_boundary_data)
# print problem
print(b)
```

As explained the package description, we can use the
`getPotentialBenefit()`

function to know what is the maximum
benefit that can be achieved for both types of planning objectives:

```
# get potential benefit
getPotentialBenefit(b)
```

```
## feature dist dist_threatened maximum.conservation.benefit maximum.recovery.benefit maximum.benefit
## 1 1 47 47 0 47 47
## 2 2 30 28 2 28 30
## 3 3 66 56 10 56 66
## 4 4 33 33 0 33 33
```

As we can observe, features 2 and 3 are present in 2 and 10 units,
respectively, where no threat is present. Therefore, if select these
units just for monitoring, they could contribute to achieving the
**conservation** target. In order to to show the impact of
including recovery requirements, we now present the results obtained
when setting the value of `target_recovery`

to 10. After
setting these values, we solve the corresponding model. These steps are
shown below.

```
# setting the targets
features_data.base <- sim_features_data
features_data.base$target_recovery <- 10
# create conservation problem using data from prioriactions
b.base <- inputData(pu = sim_pu_data,
features = features_data.base,
dist_features = sim_dist_features_data,
threats = sim_threats_data,
dist_threats = sim_dist_threats_data,
sensitivity = sim_sensitivity_data,
boundary = sim_boundary_data)
# create optimization problem
c.base <- problem(b.base, model_type = "minimizeCosts")
```

`## Warning: The blm argument was set to 0, so the boundary data has no effect`

`## Warning: Some blm_actions argument were set to 0, so the boundary data has no effect for these cases`

```
# solve optimization problem
d.base <- solve(c.base, solver = "gurobi", verbose = TRUE, output_file = FALSE, cores = 2)
```

```
## Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (linux64)
## Thread count: 2 physical cores, 4 logical processors, using up to 2 threads
## Optimize a model with 284 rows, 396 columns and 785 nonzeros
## Model fingerprint: 0xe4778898
## Variable types: 176 continuous, 220 integer (220 binary)
## Coefficient statistics:
## Matrix range [5e-01, 2e+00]
## Objective range [1e+00, 1e+01]
## Bounds range [1e+00, 1e+00]
## RHS range [1e+01, 1e+01]
## Found heuristic solution: objective 964.0000000
## Found heuristic solution: objective 371.0000000
## Presolve removed 250 rows and 277 columns
## Presolve time: 0.00s
## Presolved: 34 rows, 119 columns, 237 nonzeros
## Variable types: 0 continuous, 119 integer (101 binary)
##
## Root relaxation: objective 1.310000e+02, 30 iterations, 0.00 seconds
##
## Nodes | Current Node | Objective Bounds | Work
## Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
##
## 0 0 131.00000 0 6 371.00000 131.00000 64.7% - 0s
## H 0 0 143.0000000 131.00000 8.39% - 0s
## H 0 0 142.0000000 131.00000 7.75% - 0s
## 0 0 134.00000 0 4 142.00000 134.00000 5.63% - 0s
## H 0 0 138.0000000 134.00000 2.90% - 0s
## H 0 0 134.0000000 134.00000 0.00% - 0s
## 0 0 134.00000 0 4 134.00000 134.00000 0.00% - 0s
##
## Cutting planes:
## Gomory: 1
## Cover: 5
##
## Explored 1 nodes (39 simplex iterations) in 0.00 seconds
## Thread count was 2 (of 4 available processors)
##
## Solution count 6: 134 138 142 ... 964
##
## Optimal solution found (tolerance 0.00e+00)
## Best objective 1.340000000000e+02, best bound 1.340000000000e+02, gap 0.0000%
```

Once the corresponding optimization problem is solved, we can
retrieve the spatial distribution of the recovery actions using the
`getActions()`

function. The results obtained for our test
instance are shown below:

```
# get actions from solution
actions.base <- getActions(d.base, format = "wide")
# plot actions
group_rasters <- raster::stack(r, r)
values(group_rasters[[1]]) <- actions.base$`1`
values(group_rasters[[2]]) <- actions.base$`2`
names(group_rasters) <- c("action 1", "action 2")
plot(group_rasters)
```

In this section, we present the results obtained when only
considering **conservation targets**. Note that the value
of these targets are bounded by the maximum possible benefits that can
be obtained according to the structure of the planning area. The
corresponding value of this maximum possible benefits can be retrieved
using the `getPotentialBenefit()`

function. As we have
already seen in the previous section, these values correspond to 0, 2,
10 and 0 respectively.Once the **conservation targets**
values are defined, the corresponding optimization model is solved.
These steps are shown below.

```
# setting the targets
features_data.cons <- sim_features_data
features_data.cons$target_recovery <- 10
features_data.cons$target_conservation <- c(0, 2, 10, 0)
# create conservation problem using data from prioriactions
b.cons <- inputData(pu = sim_pu_data,
features = features_data.cons,
dist_features = sim_dist_features_data,
threats = sim_threats_data,
dist_threats = sim_dist_threats_data,
sensitivity = sim_sensitivity_data,
boundary = sim_boundary_data)
# create optimization problem
c.cons <- problem(b.cons, model_type = "minimizeCosts")
```

`## Warning: The blm argument was set to 0, so the boundary data has no effect`

`## Warning: Some blm_actions argument were set to 0, so the boundary data has no effect for these cases`

```
# solve optimization problem
d.cons <- solve(c.cons, solver = "gurobi", verbose = TRUE, output_file = FALSE, cores = 2)
```

```
## Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (linux64)
## Thread count: 2 physical cores, 4 logical processors, using up to 2 threads
## Optimize a model with 284 rows, 396 columns and 785 nonzeros
## Model fingerprint: 0xf6bb92ed
## Variable types: 176 continuous, 220 integer (220 binary)
## Coefficient statistics:
## Matrix range [5e-01, 2e+00]
## Objective range [1e+00, 1e+01]
## Bounds range [1e+00, 1e+00]
## RHS range [2e+00, 1e+01]
## Found heuristic solution: objective 964.0000000
## Found heuristic solution: objective 381.0000000
## Presolve removed 250 rows and 277 columns
## Presolve time: 0.00s
## Presolved: 34 rows, 119 columns, 237 nonzeros
## Variable types: 0 continuous, 119 integer (101 binary)
##
## Root relaxation: objective 1.540000e+02, 30 iterations, 0.00 seconds
##
## Nodes | Current Node | Objective Bounds | Work
## Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
##
## 0 0 154.00000 0 6 381.00000 154.00000 59.6% - 0s
## H 0 0 166.0000000 154.00000 7.23% - 0s
## H 0 0 165.0000000 154.00000 6.67% - 0s
## 0 0 157.00000 0 4 165.00000 157.00000 4.85% - 0s
## H 0 0 161.0000000 157.00000 2.48% - 0s
## H 0 0 157.0000000 157.00000 0.00% - 0s
## 0 0 157.00000 0 4 157.00000 157.00000 0.00% - 0s
##
## Cutting planes:
## Gomory: 1
## Cover: 5
##
## Explored 1 nodes (39 simplex iterations) in 0.01 seconds
## Thread count was 2 (of 4 available processors)
##
## Solution count 6: 157 161 165 ... 964
##
## Optimal solution found (tolerance 0.00e+00)
## Best objective 1.570000000000e+02, best bound 1.570000000000e+02, gap 0.0000%
```

As can be seen from the log, the objective function value of the attained solution is 157. This value is 23 units larger than the value obtained for the same instance but without imposing conservation targets (the values was 134). This difference is explained by the fact that we now impose the implementation of monitoring actions.

```
# get actions from solution
actions.cons <- getActions(d.cons, format = "wide")
# plot actions
group_rasters <- raster::stack(r, r, r)
values(group_rasters[[1]]) <- actions.cons$`1`
values(group_rasters[[2]]) <- actions.cons$`2`
values(group_rasters[[3]]) <- actions.cons$`conservation`
names(group_rasters) <- c("action 1", "action 2", "conservation")
plot(group_rasters)
```

We can observe that, in this case, the units selected for
**conservation** do not overlap with those selected for
applying action 1 and action 2. When comparing this solution with the
one shown in the model I, we observe that in both cases the spatial
distribution of the **recovery** actions is the same.

To demonstrate the influence of including connectivity constrains we
will use the same inputs as above, but we will set the connectivity
penalty factor (*blm*) to 100 in the `problem()`

function. Next we will modify this parameter in the model that uses only
**recovery targets**.

```
# create optimization problem
c.base_conn <- problem(b.base,
model_type = "minimizeCosts",
blm = 100)
```

`## Warning: Some blm_actions argument were set to 0, so the boundary data has no effect for these cases`

```
# solve optimization problem
d.base_conn <- solve(c.base_conn, solver = "gurobi", verbose = TRUE, output_file = FALSE, cores = 2)
```

```
## Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (linux64)
## Thread count: 2 physical cores, 4 logical processors, using up to 2 threads
## Optimize a model with 29984 rows, 10296 columns and 70085 nonzeros
## Model fingerprint: 0x0d6e31a0
## Variable types: 176 continuous, 10120 integer (10120 binary)
## Coefficient statistics:
## Matrix range [5e-01, 2e+00]
## Objective range [1e+00, 7e+04]
## Bounds range [1e+00, 1e+00]
## RHS range [1e+00, 1e+01]
## Found heuristic solution: objective 964.0000000
## Presolve removed 190 rows and 181 columns
## Presolve time: 0.22s
## Presolved: 29794 rows, 10115 columns, 69709 nonzeros
## Variable types: 0 continuous, 10115 integer (10110 binary)
##
## Root relaxation: objective 2.573214e+02, 253 iterations, 0.09 seconds
##
## Nodes | Current Node | Objective Bounds | Work
## Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
##
## * 0 0 0 547.0000000 547.00000 0.00% - 0s
##
## Explored 0 nodes (306 simplex iterations) in 0.43 seconds
## Thread count was 2 (of 4 available processors)
##
## Solution count 2: 547 964
##
## Optimal solution found (tolerance 0.00e+00)
## Best objective 5.470000000163e+02, best bound 5.469999999816e+02, gap 0.0000%
```

```
# get actions from solution
actions.base_conn <- getActions(d.base_conn, format = "wide")
# plot actions
group_rasters <- raster::stack(r, r)
values(group_rasters[[1]]) <- actions.base_conn$`1`
values(group_rasters[[2]]) <- actions.base_conn$`2`
names(group_rasters) <- c("action 1", "action 2")
plot(group_rasters)
```

We can notice that the selected actions are more connected compared to the solutions shown for the previous models.

Same as above, we set the *blm* parameter to 100, but now
using **conservation targets**.

```
# create optimization problem
c.cons_conn <- problem(b.cons,
model_type = "minimizeCosts",
blm = 100)
```

`## Warning: Some blm_actions argument were set to 0, so the boundary data has no effect for these cases`

```
# solve optimization problem
d.cons_conn <- solve(c.cons_conn, solver = "gurobi", verbose = TRUE, output_file = FALSE, cores = 2)
```

```
## Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (linux64)
## Thread count: 2 physical cores, 4 logical processors, using up to 2 threads
## Optimize a model with 29984 rows, 10296 columns and 70085 nonzeros
## Model fingerprint: 0x8ba6f640
## Variable types: 176 continuous, 10120 integer (10120 binary)
## Coefficient statistics:
## Matrix range [5e-01, 2e+00]
## Objective range [1e+00, 7e+04]
## Bounds range [1e+00, 1e+00]
## RHS range [1e+00, 1e+01]
## Found heuristic solution: objective 964.0000000
## Presolve removed 5860 rows and 2081 columns
## Presolve time: 0.21s
## Presolved: 24124 rows, 8215 columns, 56479 nonzeros
## Variable types: 0 continuous, 8215 integer (8210 binary)
##
## Root relaxation: objective 5.470000e+02, 9 iterations, 0.01 seconds
##
## Nodes | Current Node | Objective Bounds | Work
## Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
##
## * 0 0 0 547.0000000 547.00000 0.00% - 0s
##
## Explored 0 nodes (9 simplex iterations) in 0.28 seconds
## Thread count was 2 (of 4 available processors)
##
## Solution count 2: 547 964
##
## Optimal solution found (tolerance 0.00e+00)
## Best objective 5.470000000163e+02, best bound 5.469999999524e+02, gap 0.0000%
```

```
# get actions from solution
actions.cons_conn <- getActions(d.cons_conn, format = "wide")
# plot actions
group_rasters <- raster::stack(r, r, r)
values(group_rasters[[1]]) <- actions.cons_conn$`1`
values(group_rasters[[2]]) <- actions.cons_conn$`2`
values(group_rasters[[3]]) <- actions.cons_conn$`conservation`
names(group_rasters) <- c("action 1", "action 2", "conservation")
plot(group_rasters)
```

Note that in the last two models (those that incorporate connectivity
requirements) the solutions differ. Because in both we have used the
same parameters, it is inferred that the sites selected for the
conservation of species affect in a certain way the selection of units
for **recovery targets**, which is why it is very important
to carefully think about the planning objectives, as they might lead to
different solutions.