Add penalties to a conservation planning problem to favor solutions that select planning units with high connectivity between them.

# S4 method for ConservationProblem,ANY,ANY,matrix

# S4 method for ConservationProblem,ANY,ANY,Matrix

# S4 method for ConservationProblem,ANY,ANY,dgCMatrix

# S4 method for ConservationProblem,ANY,ANY,data.frame

# S4 method for ConservationProblem,ANY,ANY,array
add_connectivity_penalties(x, penalty, zones, data)

## Arguments

x ConservationProblem-class object. numeric penalty that is used to scale the importance of selecting planning units with strong connectivity between them compared to the main problem objective (e.g. solution cost when the argument to x has a minimum set objective set using add_min_set_objective). Higher penalty values can be used to obtain solutions with a high degree of connectivity, and smaller penalty values can be used to obtain solutions with a #' small degree of connectivity. Note that negative penalty values can be used to obtain solutions that have very little connectivity. matrix or Matrix object describing the level of connectivity between different zones. Each row and column corresponds to a different zone in the argument to x, and cell values indicate the level of connectivity between each combination of zones. Cell values along the diagonal of the matrix represent the level of connectivity between planning units allocated to the same zone. Cell values must lay between 1 and -1, where negative values favor solutions with weak connectivity. The default argument to zones is an identity matrix (i.e. a matrix with ones along the matrix diagonal and zeros elsewhere), so that planning units are only considered to be connected when they are allocated to the same zone. This argument is required when the argument to data is a matrix or Matrix object. If the argument to data is an array or data.frame with zone data, this argument must explicitly be set to NULL otherwise an error will be thrown. matrix, Matrix, data.frame, or array object containing connectivity data. The connectivity values correspond to the strength of connectivity between different planning units. Thus connections between planning units that are associated with higher values are more favorable in the solution. See the Details section for more information. not used.

## Value

ConservationProblem-class object with the penalties added to it.

## Details

This function uses connectivity data to penalize solutions that have low connectivity. It can accommodate symmetric and asymmetric relationships between planning units. Although Marxan penalizes connections between planning units with high connectivity values, it is important to note that this function favors connections between planning units with high connectivity values. This function was inspired by Beger et al. (2010).

The argument to data can be specified in several different ways:

matrix, Matrix

where rows and columns represent different planning units and the value of each cell represents the strength of connectivity between two different planning units. Cells that occur along the matrix diagonal are treated as weights which indicate that planning units are more desirable in the solution. The argument to zones can be used to control the strength of connectivity between planning units in different zones. The default argument for zones is to treat planning units allocated to different zones as having zero connectivity.

data.frame

containing the fields (columns) "id1", "id2", and "boundary". Here, each row denotes the connectivity between two planning units following the Marxan format. The data can be used to denote symmetric or asymmetric relationships between planning units. By default, input data is assumed to be symmetric unless asymmetric data is also included (e.g. if data is present for planning units 2 and 3, then the same amount of connectivity is expected for planning units 3 and 2, unless connectivity data is also provided for planning units 3 and 2). If the argument to x contains multiple zones, then the columns "zone1" and "zone2" can optionally be provided to manually specify the connectivity values between planning units when they are allocated to specific zones. If the columns "zone1" and "zone2" are present, then the argument to zones must be NULL.

array

containing four-dimensions where cell values indicate the strength of connectivity between planning units when they are assigned to specific management zones. The first two dimensions (i.e. rows and columns) indicate the strength of connectivity between different planning units and the second two dimensions indicate the different management zones. Thus the data[1, 2, 3, 4] indicates the strength of connectivity between planning unit 1 and planning unit 2 when planning unit 1 is assigned to zone 3 and planning unit 2 is assigned to zone 4.

The connectivity penalties are calculated using the following equations. Let $$I$$ represent the set of planning units, $$Z$$ represent the set of management zones, and $$X_iz$$ represent the decision variable for planning unit $$i$$ for in zone $$z$$ (e.g. with binary values one indicating if planning unit is allocated or not). Also, let $$p$$ represent the argument to penalty, $$D$$ represent the argument to data, and $$W$$ represent the argument to zones.

If the argument to data is supplied as a matrix or Matrix, then the penalties are calculated as:

$$\sum_{i}^{I} \sum_{j}^{I} \sum_{z}^{Z} \sum_{y}^{Z} (-p \times X_{iz} \times X_{jy} \times D_{ij} \times W_{zy})$$

Otherwise, if the argument to data is supplied as a data.frame or array, then the penalties are calculated as:

$$\sum_{i}^{I} \sum_{j}^{I} \sum_{z}^{Z} \sum_{y}^{Z} (-p \times X_{iz} \times X_{jy} \times D_{ijzy})$$

Note that when the problem objective is to maximize some measure of benefit and not minimize some measure of cost, the term $$-p$$ is replaced with $$p$$.

## References

Beger M, Linke S, Watts M, Game E, Treml E, Ball I, and Possingham, HP (2010) Incorporating asymmetric connectivity into spatial decision making for conservation, Conservation Letters, 3: 359--368.

penalties.

## Examples

# load Matrix package for visualizing matrices
data(sim_pu_polygons, sim_pu_zones_stack, sim_features, sim_features_zones)

# define function to rescale values between zero and one so that we
# can compare solutions from different connectivity matrices
rescale <- function(x, to = c(0, 1), from = range(x, na.rm = TRUE)) {
(x - from[1]) / diff(from) * diff(to) + to[1]
}

# create basic problem
p1 <- problem(sim_pu_polygons, sim_features, "cost") %>%

# create a symmetric connectivity matrix where the connectivity between
# two planning units corresponds to their shared boundary length
b_matrix <- boundary_matrix(sim_pu_polygons)

# standardize matrix values to lay between zero and one
b_matrix[] <- rescale(b_matrix[])

# visualize connectivity matrix
image(b_matrix)
# create a symmetric connectivity matrix where the connectivity between
# two planning units corresponds to their spatial proximity
# i.e. planning units that are further apart share less connectivity
centroids <- rgeos::gCentroid(sim_pu_polygons, byid = TRUE)
d_matrix <- (1 / (as(dist(centroids@coords), "Matrix") + 1))

# standardize matrix values to lay between zero and one
d_matrix[] <- rescale(d_matrix[])

# remove connections between planning units without connectivity to
# reduce run-time
d_matrix[d_matrix < 0.7] <- 0

# visualize connectivity matrix
image(d_matrix)
# create a symmetric connectivity matrix where the connectivity
# between adjacent two planning units corresponds to their combined
# value in a field in the planning unit attribute data
# for example, this field could describe the extent of native vegetation in
# each planning unit and we could use connectivity penalties to identify
# solutions that cluster planning units together that both contain large
# amounts of native vegetation
c_matrix <- connectivity_matrix(sim_pu_polygons, "cost")

# standardize matrix values to lay between zero and one
c_matrix[] <- rescale(c_matrix[])

# visualize connectivity matrix
image(c_matrix)
# create an asymmetric connectivity matrix. Here, connectivity occurs between
# adjacent planning units and, due to rivers flowing southwards
# through the study area, connectivity from northern planning units to
# southern planning units is ten times stronger than the reverse.
ac_matrix <- matrix(0, length(sim_pu_polygons), length(sim_pu_polygons))
ac_matrix <- as(ac_matrix, "Matrix")
adjacent_units <- rgeos::gIntersects(sim_pu_polygons, byid = TRUE)
for (i in seq_len(length(sim_pu_polygons))) {
for (j in seq_len(length(sim_pu_polygons))) {
# find if planning units are adjacent
# find if planning units lay north and south of each other
# i.e. they have the same x-coordinate
if (centroids@coords[i, 1] == centroids@coords[j, 1]) {
if (centroids@coords[i, 2] > centroids@coords[j, 2]) {
# if i is north of j add 10 units of connectivity
ac_matrix[i, j] <- ac_matrix[i, j] + 10
} else if (centroids@coords[i, 2] < centroids@coords[j, 2]) {
# if i is south of j add 1 unit of connectivity
ac_matrix[i, j] <- ac_matrix[i, j] + 1
}
}
}
}
}

# standardize matrix values to lay between zero and one
ac_matrix[] <- rescale(ac_matrix[])

# visualize asymmetric connectivity matrix
image(ac_matrix)
# create penalties
penalties <- c(10, 25)

# create problems using the different connectivity matrices and penalties
p2 <- list(p1,
p1 %>% add_connectivity_penalties(penalties[1], data = b_matrix),
p1 %>% add_connectivity_penalties(penalties[2], data = b_matrix),
p1 %>% add_connectivity_penalties(penalties[1], data = d_matrix),
p1 %>% add_connectivity_penalties(penalties[2], data = d_matrix),
p1 %>% add_connectivity_penalties(penalties[1], data = c_matrix),
p1 %>% add_connectivity_penalties(penalties[2], data = c_matrix),
p1 %>% add_connectivity_penalties(penalties[1], data = ac_matrix),
p1 %>% add_connectivity_penalties(penalties[2], data = ac_matrix))

# assign names to the problems
names(p2) <- c("basic problem",
paste0("b_matrix (", penalties,")"),
paste0("d_matrix (", penalties,")"),
paste0("c_matrix (", penalties,")"),
paste0("ac_matrix (", penalties,")"))# solve problems
s2 <- lapply(p2, solve)#> Optimize a model with 5 rows, 90 columns and 450 nonzeros
#> 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        [6e+00, 1e+01]
#> Found heuristic solution: objective 3934.6218396
#> 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 3.496032e+03, 16 iterations, 0.00 seconds
#>
#>     Nodes    |    Current Node    |     Objective Bounds      |     Work
#>  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
#>
#>      0     0 3496.03193    0    4 3934.62184 3496.03193  11.1%     -    0s
#> H    0     0                    3585.9601335 3496.03193  2.51%     -    0s
#>
#> Explored 1 nodes (16 simplex iterations) in 0.00 seconds
#> Thread count was 1 (of 4 available processors)
#>
#> Solution count 2: 3585.96 3934.62
#>
#> Optimal solution found (tolerance 1.00e-01)
#> Best objective 3.585960133519e+03, best bound 3.496031931890e+03, gap 2.5078%
#> Optimize a model with 581 rows, 378 columns and 1602 nonzeros
#> Variable types: 0 continuous, 378 integer (378 binary)
#> Coefficient statistics:
#>   Matrix range     [2e-01, 1e+00]
#>   Objective range  [3e+00, 2e+02]
#>   Bounds range     [1e+00, 1e+00]
#>   RHS range        [6e+00, 1e+01]
#> Found heuristic solution: objective 17287.196992
#> Found heuristic solution: objective 3691.2885062
#> Presolve time: 0.01s
#> Presolved: 581 rows, 378 columns, 1602 nonzeros
#> Variable types: 0 continuous, 378 integer (378 binary)
#> Presolved: 581 rows, 378 columns, 1602 nonzeros
#>
#>
#> Root relaxation: objective 3.307193e+03, 159 iterations, 0.01 seconds
#>
#>     Nodes    |    Current Node    |     Objective Bounds      |     Work
#>  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
#>
#>      0     0 3307.19310    0   86 3691.28851 3307.19310  10.4%     -    0s
#> H    0     0                    3433.8940898 3307.19310  3.69%     -    0s
#>
#> Explored 1 nodes (159 simplex iterations) in 0.02 seconds
#> Thread count was 1 (of 4 available processors)
#>
#> Solution count 3: 3433.89 3691.29 17287.2
#>
#> Optimal solution found (tolerance 1.00e-01)
#> Best objective 3.433894089811e+03, best bound 3.307193102718e+03, gap 3.6897%
#> Optimize a model with 581 rows, 378 columns and 1602 nonzeros
#> Variable types: 0 continuous, 378 integer (378 binary)
#> Coefficient statistics:
#>   Matrix range     [2e-01, 1e+00]
#>   Objective range  [8e+00, 2e+02]
#>   Bounds range     [1e+00, 1e+00]
#>   RHS range        [6e+00, 1e+01]
#> Found heuristic solution: objective 15487.196992
#> Found heuristic solution: objective 3326.2885062
#> Presolve time: 0.01s
#> Presolved: 581 rows, 378 columns, 1602 nonzeros
#> Variable types: 0 continuous, 378 integer (378 binary)
#> Presolved: 581 rows, 378 columns, 1602 nonzeros
#>
#>
#> Root relaxation: objective 3.012030e+03, 174 iterations, 0.01 seconds
#>
#>     Nodes    |    Current Node    |     Objective Bounds      |     Work
#>  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
#>
#>      0     0 3012.03046    0  138 3326.28851 3012.03046  9.45%     -    0s
#>
#> Explored 1 nodes (174 simplex iterations) in 0.02 seconds
#> Thread count was 1 (of 4 available processors)
#>
#> Solution count 2: 3326.29 15487.2
#>
#> Optimal solution found (tolerance 1.00e-01)
#> Best objective 3.326288506226e+03, best bound 3.012030457072e+03, gap 9.4477%
#> Optimize a model with 1601 rows, 888 columns and 3642 nonzeros
#> Variable types: 0 continuous, 888 integer (888 binary)
#> Coefficient statistics:
#>   Matrix range     [2e-01, 1e+00]
#>   Objective range  [7e+00, 2e+02]
#>   Bounds range     [1e+00, 1e+00]
#>   RHS range        [6e+00, 1e+01]
#> Found heuristic solution: objective 11397.038660
#> Found heuristic solution: objective 2915.9453686
#> Presolve time: 0.01s
#> Presolved: 1601 rows, 888 columns, 3642 nonzeros
#> Variable types: 0 continuous, 888 integer (888 binary)
#> Presolved: 1601 rows, 888 columns, 3642 nonzeros
#>
#>
#> Root relaxation: objective 2.277325e+03, 601 iterations, 0.05 seconds
#>
#>     Nodes    |    Current Node    |     Objective Bounds      |     Work
#>  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
#>
#>      0     0 2277.32517    0  877 2915.94537 2277.32517  21.9%     -    0s
#>      0     0 2282.91603    0  874 2915.94537 2282.91603  21.7%     -    3s
#>      0     0 2288.19139    0  871 2915.94537 2288.19139  21.5%     -    4s
#>      0     0 2288.19139    0  483 2915.94537 2288.19139  21.5%     -    4s
#>      0     0 2295.81836    0  478 2915.94537 2295.81836  21.3%     -    5s
#>      0     0 2299.18104    0  478 2915.94537 2299.18104  21.2%     -    5s
#>      0     0 2310.18554    0  450 2915.94537 2310.18554  20.8%     -    5s
#>      0     0 2314.83207    0  473 2915.94537 2314.83207  20.6%     -    5s
#>      0     0 2316.71918    0  467 2915.94537 2316.71918  20.5%     -    5s
#>      0     0 2316.71918    0  467 2915.94537 2316.71918  20.5%     -    5s
#> H    0     0                    2879.5923214 2316.71918  19.5%     -    6s
#>      0     2 2316.75143    0  467 2879.59232 2316.75143  19.5%     -    6s
#> H   53    53                    2669.2652017 2367.46532  11.3%  45.8    6s
#> H  273   209                    2665.1537029 2389.53644  10.3%  56.5    7s
#>
#> Cutting planes:
#>   Gomory: 6
#>
#> Explored 330 nodes (18882 simplex iterations) in 7.68 seconds
#> Thread count was 1 (of 4 available processors)
#>
#> Solution count 5: 2665.15 2669.27 2879.59 ... 11397
#>
#> Optimal solution found (tolerance 1.00e-01)
#> Best objective 2.665153702866e+03, best bound 2.411089064269e+03, gap 9.5328%
#> Optimize a model with 1601 rows, 888 columns and 3642 nonzeros
#> Variable types: 0 continuous, 888 integer (888 binary)
#> Coefficient statistics:
#>   Matrix range     [2e-01, 1e+00]
#>   Objective range  [2e+01, 2e+02]
#>   Bounds range     [1e+00, 1e+00]
#>   RHS range        [6e+00, 1e+01]
#> Found heuristic solution: objective 761.8011610
#> Presolve time: 0.01s
#> Presolved: 1601 rows, 888 columns, 3642 nonzeros
#> Variable types: 0 continuous, 888 integer (888 binary)
#> Presolved: 1601 rows, 888 columns, 3642 nonzeros
#>
#>
#> Root relaxation: objective 1.267850e+02, 243 iterations, 0.02 seconds
#>
#>     Nodes    |    Current Node    |     Objective Bounds      |     Work
#>  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
#>
#>      0     0  126.78496    0  695  761.80116  126.78496  83.4%     -    0s
#> H    0     0                     468.0158310  126.78496  72.9%     -    0s
#>      0     0  139.58775    0  693  468.01583  139.58775  70.2%     -    3s
#>      0     0  148.83710    0  691  468.01583  148.83710  68.2%     -    3s
#>      0     0  148.83710    0  382  468.01583  148.83710  68.2%     -    3s
#>      0     0  161.29257    0  380  468.01583  161.29257  65.5%     -    3s
#>      0     0  170.44305    0  348  468.01583  170.44305  63.6%     -    3s
#>      0     0  182.28500    0  352  468.01583  182.28500  61.1%     -    4s
#> H    0     0                     454.2214760  182.28500  59.9%     -    4s
#>      0     0  194.23505    0  378  454.22148  194.23505  57.2%     -    4s
#>      0     0  194.23505    0  378  454.22148  194.23505  57.2%     -    4s
#>      0     2  195.12431    0  378  454.22148  195.12431  57.0%     -    4s
#> H    1     1                     423.4053090  195.12431  53.9%   0.0    4s
#> *   49     3               2     421.3259805  329.15212  21.9%  41.5    4s
#>
#> Cutting planes:
#>   Gomory: 5
#>
#> Explored 52 nodes (2800 simplex iterations) in 4.37 seconds
#> Thread count was 1 (of 4 available processors)
#>
#> Solution count 5: 421.326 423.405 454.221 ... 761.801
#>
#> Optimal solution found (tolerance 1.00e-01)
#> Best objective 4.213259804528e+02, best bound 3.823629791434e+02, gap 9.2477%
#> Optimize a model with 581 rows, 378 columns and 1602 nonzeros
#> Variable types: 0 continuous, 378 integer (378 binary)
#> Coefficient statistics:
#>   Matrix range     [2e-01, 1e+00]
#>   Objective range  [9e+00, 2e+02]
#>   Bounds range     [1e+00, 1e+00]
#>   RHS range        [6e+00, 1e+01]
#> Found heuristic solution: objective 15721.108025
#> Found heuristic solution: objective 3513.2283384
#> Presolve time: 0.01s
#> Presolved: 581 rows, 378 columns, 1602 nonzeros
#> Variable types: 0 continuous, 378 integer (378 binary)
#> Presolved: 581 rows, 378 columns, 1602 nonzeros
#>
#>
#> Root relaxation: objective 3.099126e+03, 223 iterations, 0.01 seconds
#>
#>     Nodes    |    Current Node    |     Objective Bounds      |     Work
#>  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
#>
#>      0     0 3099.12568    0  190 3513.22834 3099.12568  11.8%     -    0s
#>      0     0 3106.44940    0  164 3513.22834 3106.44940  11.6%     -    0s
#> H    0     0                    3244.5496513 3106.44940  4.26%     -    0s
#>
#> Cutting planes:
#>   Gomory: 1
#>
#> Explored 1 nodes (267 simplex iterations) in 0.06 seconds
#> Thread count was 1 (of 4 available processors)
#>
#> Solution count 3: 3244.55 3513.23 15721.1
#>
#> Optimal solution found (tolerance 1.00e-01)
#> Best objective 3.244549651264e+03, best bound 3.106449395648e+03, gap 4.2564%
#> Optimize a model with 581 rows, 378 columns and 1602 nonzeros
#> Variable types: 0 continuous, 378 integer (378 binary)
#> Coefficient statistics:
#>   Matrix range     [2e-01, 1e+00]
#>   Objective range  [2e+01, 2e+02]
#>   Bounds range     [1e+00, 1e+00]
#>   RHS range        [6e+00, 1e+01]
#> Found heuristic solution: objective 11571.974576
#> Found heuristic solution: objective 2881.1380866
#> Presolve time: 0.01s
#> Presolved: 581 rows, 378 columns, 1602 nonzeros
#> Variable types: 0 continuous, 378 integer (378 binary)
#> Presolved: 581 rows, 378 columns, 1602 nonzeros
#>
#>
#> Root relaxation: objective 2.309063e+03, 214 iterations, 0.01 seconds
#>
#>     Nodes    |    Current Node    |     Objective Bounds      |     Work
#>  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
#>
#>      0     0 2309.06312    0  362 2881.13809 2309.06312  19.9%     -    0s
#>      0     0 2328.00732    0  359 2881.13809 2328.00732  19.2%     -    0s
#>      0     0 2339.73282    0  365 2881.13809 2339.73282  18.8%     -    0s
#>      0     0 2339.73282    0  224 2881.13809 2339.73282  18.8%     -    0s
#>      0     0 2346.29691    0  224 2881.13809 2346.29691  18.6%     -    0s
#>      0     0 2357.20173    0  222 2881.13809 2357.20173  18.2%     -    0s
#>      0     0 2360.84727    0  226 2881.13809 2360.84727  18.1%     -    0s
#>      0     0 2376.63837    0  223 2881.13809 2376.63837  17.5%     -    0s
#>      0     0 2376.63837    0  223 2881.13809 2376.63837  17.5%     -    0s
#> H    0     0                    2633.0995862 2376.63837  9.74%     -    0s
#>
#> Cutting planes:
#>   Gomory: 4
#>
#> Explored 1 nodes (717 simplex iterations) in 0.49 seconds
#> Thread count was 1 (of 4 available processors)
#>
#> Solution count 3: 2633.1 2881.14 11572
#>
#> Optimal solution found (tolerance 1.00e-01)
#> Best objective 2.633099586223e+03, best bound 2.376638371787e+03, gap 9.7399%
#> Optimize a model with 293 rows, 234 columns and 1026 nonzeros
#> Variable types: 0 continuous, 234 integer (234 binary)
#> Coefficient statistics:
#>   Matrix range     [2e-01, 1e+00]
#>   Objective range  [1e+00, 2e+02]
#>   Bounds range     [1e+00, 1e+00]
#>   RHS range        [6e+00, 1e+01]
#> Found heuristic solution: objective 17098.178722
#> Found heuristic solution: objective 3824.6218396
#> Presolve time: 0.00s
#> Presolved: 293 rows, 234 columns, 1026 nonzeros
#> Variable types: 0 continuous, 234 integer (234 binary)
#> Presolved: 293 rows, 234 columns, 1026 nonzeros
#>
#>
#> Root relaxation: objective 3.385781e+03, 49 iterations, 0.00 seconds
#>
#>     Nodes    |    Current Node    |     Objective Bounds      |     Work
#>  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
#>
#>      0     0 3385.78126    0   24 3824.62184 3385.78126  11.5%     -    0s
#> H    0     0                    3506.9306268 3385.78126  3.45%     -    0s
#>
#> Explored 1 nodes (49 simplex iterations) in 0.01 seconds
#> Thread count was 1 (of 4 available processors)
#>
#> Solution count 3: 3506.93 3824.62 17098.2
#>
#> Optimal solution found (tolerance 1.00e-01)
#> Best objective 3.506930626812e+03, best bound 3.385781260518e+03, gap 3.4546%
#> Optimize a model with 293 rows, 234 columns and 1026 nonzeros
#> Variable types: 0 continuous, 234 integer (234 binary)
#> Coefficient statistics:
#>   Matrix range     [2e-01, 1e+00]
#>   Objective range  [2e+00, 2e+02]
#>   Bounds range     [1e+00, 1e+00]
#>   RHS range        [6e+00, 1e+01]
#> Found heuristic solution: objective 15910.178722
#> Found heuristic solution: objective 3659.6218396
#> Presolve time: 0.00s
#> Presolved: 293 rows, 234 columns, 1026 nonzeros
#> Variable types: 0 continuous, 234 integer (234 binary)
#> Presolved: 293 rows, 234 columns, 1026 nonzeros
#>
#>
#> Root relaxation: objective 3.172486e+03, 93 iterations, 0.00 seconds
#>
#>     Nodes    |    Current Node    |     Objective Bounds      |     Work
#>  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
#>
#>      0     0 3172.48586    0   64 3659.62184 3172.48586  13.3%     -    0s
#> H    0     0                    3386.1100237 3172.48586  6.31%     -    0s
#>
#> Explored 1 nodes (93 simplex iterations) in 0.01 seconds
#> Thread count was 1 (of 4 available processors)
#>
#> Solution count 3: 3386.11 3659.62 15910.2
#>
#> Optimal solution found (tolerance 1.00e-01)
#> Best objective 3.386110023737e+03, best bound 3.172485864952e+03, gap 6.3088%
# plot solutions
par(mfrow = c(3, 3))
for (i in seq_along(s2)) {
plot(s2[[i]], main = names(p2)[i], cex = 1.5, col = "white")
plot(s2[[i]][s2[[i]]\$solution_1 == 1, ], col = "darkgreen", add = TRUE)
}
# create minimal multi-zone problem and limit solver to one minute
# to obtain solutions in a short period of time
p3 <- problem(sim_pu_zones_stack, sim_features_zones) %>%
add_relative_targets(matrix(0.15, nrow = 5, ncol = 3)) %>%

# create matrix showing which planning units are adjacent to other units
a_matrix <- connected_matrix(sim_pu_zones_stack)

# visualize matrix
image(a_matrix)
# create a zone matrix where connectivities are only present between
# planning units that are allocated to the same zone
zm1 <- as(diag(3), "Matrix")

# print zone matrix
print(zm1)#>      [,1] [,2] [,3]
#> [1,]    1    .    .
#> [2,]    .    1    .
#> [3,]    .    .    1
# create a zone matrix where connectivities are strongest between
# planning units allocated to different zones
zm2 <- matrix(1, ncol = 3, nrow = 3)
diag(zm2) <- 0
zm2 <- as(zm2, "Matrix")

# print zone matrix
print(zm2)#> 3 x 3 Matrix of class "dsyMatrix"
#>      [,1] [,2] [,3]
#> [1,]    0    1    1
#> [2,]    1    0    1
#> [3,]    1    1    0
# create a zone matrix that indicates that connectivities between planning
# units assigned to the same zone are much higher than connectivities
# assigned to different zones
zm3 <- matrix(0.1, ncol = 3, nrow = 3)
diag(zm3) <- 1
zm3 <- as(zm3, "Matrix")

# print zone matrix
print(zm3)#> 3 x 3 Matrix of class "dsyMatrix"
#>      [,1] [,2] [,3]
#> [1,]  1.0  0.1  0.1
#> [2,]  0.1  1.0  0.1
#> [3,]  0.1  0.1  1.0
# create a zone matrix that indicates that connectivities between planning
# units allocated to zone 1 are very high, connectivities between planning
# units allocated to zones 1 and 2 are moderately high, and connectivities
# planning units allocated to other zones are low
zm4 <- matrix(0.1, ncol = 3, nrow = 3)
zm4[1, 1] <- 1
zm4[1, 2] <- 0.5
zm4[2, 1] <- 0.5
zm4 <- as(zm4, "Matrix")

# print zone matrix
print(zm4)#> 3 x 3 Matrix of class "dsyMatrix"
#>      [,1] [,2] [,3]
#> [1,]  1.0  0.5  0.1
#> [2,]  0.5  0.1  0.1
#> [3,]  0.1  0.1  0.1
# create a zone matrix with strong connectivities between planning units
# allocated to the same zone, moderate connectivities between planning
# unit allocated to zone 1 and zone 2, and negative connectivities between
# planning units allocated to zone 3 and the other two zones
zm5 <- matrix(-1, ncol = 3, nrow = 3)
zm5[1, 2] <- 0.5
zm5[2, 1] <- 0.5
diag(zm5) <- 1
zm5 <- as(zm5, "Matrix")

# print zone matrix
print(zm5)#> 3 x 3 Matrix of class "dsyMatrix"
#>      [,1] [,2] [,3]
#> [1,]  1.0  0.5   -1
#> [2,]  0.5  1.0   -1
#> [3,] -1.0 -1.0    1
# create vector of penalties to use creating problems
penalties2 <- c(5, 30)

# create multi-zone problems using the adjacent connectivity matrix and
# different zone matrices
p4 <- list(
p3,

# assign names to the problems
names(p4) <- c("basic problem",
paste0("zm", rep(seq_len(5), each = 2), " (",
rep(penalties2, 2), ")"))# solve problems
s4 <- lapply(p4, solve)#> Optimize a model with 105 rows, 270 columns and 1620 nonzeros
#> 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 10275.924144
#> Presolve time: 0.00s
#> Presolved: 105 rows, 270 columns, 1620 nonzeros
#> Variable types: 0 continuous, 270 integer (270 binary)
#> Presolved: 105 rows, 270 columns, 1620 nonzeros
#>
#>
#> Root relaxation: objective 8.943504e+03, 123 iterations, 0.00 seconds
#>
#>     Nodes    |    Current Node    |     Objective Bounds      |     Work
#>  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
#>
#>      0     0 8943.50352    0   15 10275.9241 8943.50352  13.0%     -    0s
#> H    0     0                    9111.3117236 8943.50352  1.84%     -    0s
#>
#> Explored 1 nodes (123 simplex iterations) in 0.01 seconds
#> Thread count was 1 (of 4 available processors)
#>
#> Solution count 2: 9111.31 10275.9
#>
#> Optimal solution found (tolerance 1.00e-01)
#> Best objective 9.111311723633e+03, best bound 8.943503517801e+03, gap 1.8418%
#> Optimize a model with 1845 rows, 1140 columns and 5100 nonzeros
#> Variable types: 0 continuous, 1140 integer (1140 binary)
#> Coefficient statistics:
#>   Matrix range     [2e-01, 1e+00]
#>   Objective range  [5e+00, 2e+02]
#>   Bounds range     [1e+00, 1e+00]
#>   RHS range        [1e+00, 1e+01]
#> Found heuristic solution: objective 10075.924144
#> Presolve time: 0.03s
#> Presolved: 1845 rows, 1140 columns, 5100 nonzeros
#> Variable types: 0 continuous, 1140 integer (1140 binary)
#> Presolved: 1845 rows, 1140 columns, 5100 nonzeros
#>
#>
#> Root relaxation: objective 8.367447e+03, 811 iterations, 0.06 seconds
#>
#>     Nodes    |    Current Node    |     Objective Bounds      |     Work
#>  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
#>
#>      0     0 8367.44688    0  465 10075.9241 8367.44688  17.0%     -    0s
#> H    0     0                    9108.8765843 8367.44688  8.14%     -    0s
#>
#> Explored 1 nodes (811 simplex iterations) in 0.11 seconds
#> Thread count was 1 (of 4 available processors)
#>
#> Solution count 2: 9108.88 10075.9
#>
#> Optimal solution found (tolerance 1.00e-01)
#> Best objective 9.108876584257e+03, best bound 8.367446876836e+03, gap 8.1396%
#> Optimize a model with 1845 rows, 1140 columns and 5100 nonzeros
#> Variable types: 0 continuous, 1140 integer (1140 binary)
#> Coefficient statistics:
#>   Matrix range     [2e-01, 1e+00]
#>   Objective range  [3e+01, 2e+02]
#>   Bounds range     [1e+00, 1e+00]
#>   RHS range        [1e+00, 1e+01]
#> Found heuristic solution: objective 9075.9241439
#> Presolve time: 0.02s
#> Presolved: 1845 rows, 1140 columns, 5100 nonzeros
#> Variable types: 0 continuous, 1140 integer (1140 binary)
#> Presolved: 1845 rows, 1140 columns, 5100 nonzeros
#>
#>
#> Root relaxation: objective 4.897305e+03, 713 iterations, 0.04 seconds
#>
#>     Nodes    |    Current Node    |     Objective Bounds      |     Work
#>  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
#>
#>      0     0 4897.30484    0  786 9075.92414 4897.30484  46.0%     -    0s
#> H    0     0                    7701.5492676 4897.30484  36.4%     -    0s
#>      0     0 4956.79223    0  788 7701.54927 4956.79223  35.6%     -    0s
#> H    0     0                    7443.0424128 4956.79223  33.4%     -    0s
#>      0     0 4973.38838    0  752 7443.04241 4973.38838  33.2%     -    0s
#> H    0     0                    7375.9473376 4973.38838  32.6%     -    0s
#>      0     0 4973.38838    0  487 7375.94734 4973.38838  32.6%     -    0s
#> H    0     0                    7136.8817068 4973.38838  30.3%     -    0s
#>      0     0 4980.10607    0  517 7136.88171 4980.10607  30.2%     -    0s
#> H    0     0                    7021.4892757 4980.10607  29.1%     -    0s
#>      0     0 4988.31454    0  524 7021.48928 4988.31454  29.0%     -    0s
#>      0     0 5000.57640    0  546 7021.48928 5000.57640  28.8%     -    0s
#>      0     0 5024.94548    0  545 7021.48928 5024.94548  28.4%     -    1s
#> H    0     0                    6909.9319607 5024.94548  27.3%     -    1s
#>      0     0 5028.64284    0  529 6909.93196 5028.64284  27.2%     -    1s
#>      0     0 5031.68205    0  541 6909.93196 5031.68205  27.2%     -    1s
#>      0     0 5031.92221    0  539 6909.93196 5031.92221  27.2%     -    1s
#>      0     0 5032.17738    0  539 6909.93196 5032.17738  27.2%     -    1s
#>      0     0 5032.17738    0  539 6909.93196 5032.17738  27.2%     -    1s
#> H    0     0                    6422.0997997 5032.17738  21.6%     -    1s
#> H    0     0                    6370.2361980 5032.17738  21.0%     -    1s
#> H    0     0                    6367.2149957 5032.17738  21.0%     -    1s
#>      0     2 5032.69826    0  539 6367.21500 5032.69826  21.0%     -    1s
#> H   27    27                    6364.5405799 5056.80746  20.5%  63.9    1s
#> H   85    85                    6060.3446757 5056.80746  16.6%  39.9    1s
#> H  108   108                    6004.2097530 5056.80746  15.8%  34.8    1s
#> H  718   579                    6000.2969364 5104.91957  14.9%  26.7    4s
#> H  771   550                    5997.6477267 5104.91957  14.9%  25.7    4s
#> H  803   532                    5985.9485617 5104.91957  14.7%  27.2    4s
#>    843   560 5627.97393   41  338 5985.94856 5104.91957  14.7%  28.6    5s
#> H  929   526                    5849.2645301 5108.97327  12.7%  29.8    5s
#> H 1037   419                    5666.5547040 5112.30253  9.78%  30.1    5s
#>
#> Cutting planes:
#>   Gomory: 21
#>
#> Explored 1037 nodes (33062 simplex iterations) in 5.86 seconds
#> Thread count was 1 (of 4 available processors)
#>
#> Solution count 10: 5666.55 5849.26 5985.95 ... 6370.24
#>
#> Optimal solution found (tolerance 1.00e-01)
#> Best objective 5.666554703976e+03, best bound 5.112302534715e+03, gap 9.7811%
#> Optimize a model with 3585 rows, 2010 columns and 8580 nonzeros
#> Variable types: 0 continuous, 2010 integer (2010 binary)
#> Coefficient statistics:
#>   Matrix range     [2e-01, 1e+00]
#>   Objective range  [5e+00, 2e+02]
#>   Bounds range     [1e+00, 1e+00]
#>   RHS range        [1e+00, 1e+01]
#> Found heuristic solution: objective 9985.9241439
#> Presolve time: 0.08s
#> Presolved: 3585 rows, 2010 columns, 12060 nonzeros
#> Variable types: 0 continuous, 2010 integer (2010 binary)
#> Presolve removed 870 rows and 0 columns
#> Presolved: 2715 rows, 2010 columns, 9450 nonzeros
#>
#>
#> Root relaxation: objective 8.346175e+03, 2630 iterations, 0.20 seconds
#>
#>     Nodes    |    Current Node    |     Objective Bounds      |     Work
#>  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
#>
#>      0     0 8346.17500    0  329 9985.92414 8346.17500  16.4%     -    0s
#> H    0     0                    8807.8710351 8346.17500  5.24%     -    0s
#>
#> Explored 1 nodes (2996 simplex iterations) in 0.32 seconds
#> Thread count was 1 (of 4 available processors)
#>
#> Solution count 2: 8807.87 9985.92
#>
#> Optimal solution found (tolerance 1.00e-01)
#> Best objective 8.807871035125e+03, best bound 8.346175000834e+03, gap 5.2419%
#> Optimize a model with 3585 rows, 2010 columns and 8580 nonzeros
#> Variable types: 0 continuous, 2010 integer (2010 binary)
#> Coefficient statistics:
#>   Matrix range     [2e-01, 1e+00]
#>   Objective range  [3e+01, 2e+02]
#>   Bounds range     [1e+00, 1e+00]
#>   RHS range        [1e+00, 1e+01]
#> Found heuristic solution: objective 8535.9241439
#> Presolve time: 0.08s
#> Presolved: 3585 rows, 2010 columns, 12060 nonzeros
#> Variable types: 0 continuous, 2010 integer (2010 binary)
#> Presolve removed 870 rows and 0 columns
#> Presolved: 2715 rows, 2010 columns, 9450 nonzeros
#>
#>
#> Root relaxation: objective 4.854963e+03, 3234 iterations, 0.30 seconds
#>
#>     Nodes    |    Current Node    |     Objective Bounds      |     Work
#>  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
#>
#>      0     0 4854.96349    0  525 8535.92414 4854.96349  43.1%     -    0s
#> H    0     0                    6916.7028058 4854.96349  29.8%     -    0s
#> H    0     0                    5326.8802741 4854.96349  8.86%     -    0s
#>
#> Explored 1 nodes (3727 simplex iterations) in 0.57 seconds
#> Thread count was 1 (of 4 available processors)
#>
#> Solution count 3: 5326.88 6916.7 8535.92
#>
#> Optimal solution found (tolerance 1.00e-01)
#> Best objective 5.326880274073e+03, best bound 4.854963486880e+03, gap 8.8592%
#> Optimize a model with 5325 rows, 2880 columns and 12060 nonzeros
#> Variable types: 0 continuous, 2880 integer (2880 binary)
#> Coefficient statistics:
#>   Matrix range     [2e-01, 1e+00]
#>   Objective range  [5e-01, 2e+02]
#>   Bounds range     [1e+00, 1e+00]
#>   RHS range        [1e+00, 1e+01]
#> Found heuristic solution: objective 10046.924144
#> Presolve time: 0.18s
#> Presolved: 5325 rows, 2880 columns, 17280 nonzeros
#> Variable types: 0 continuous, 2880 integer (2880 binary)
#> Presolve removed 870 rows and 0 columns
#> Presolved: 4455 rows, 2880 columns, 14670 nonzeros
#>
#>
#> Root relaxation: objective 8.343967e+03, 2550 iterations, 0.26 seconds
#>
#>     Nodes    |    Current Node    |     Objective Bounds      |     Work
#>  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
#>
#>      0     0 8343.96746    0  835 10046.9241 8343.96746  17.0%     -    0s
#> H    0     0                    9212.6916242 8343.96746  9.43%     -    0s
#>
#> Explored 1 nodes (2792 simplex iterations) in 0.56 seconds
#> Thread count was 1 (of 4 available processors)
#>
#> Solution count 2: 9212.69 10046.9
#>
#> Optimal solution found (tolerance 1.00e-01)
#> Best objective 9.212691624211e+03, best bound 8.343967460484e+03, gap 9.4296%
#> Optimize a model with 5325 rows, 2880 columns and 12060 nonzeros
#> Variable types: 0 continuous, 2880 integer (2880 binary)
#> Coefficient statistics:
#>   Matrix range     [2e-01, 1e+00]
#>   Objective range  [3e+00, 2e+02]
#>   Bounds range     [1e+00, 1e+00]
#>   RHS range        [1e+00, 1e+01]
#> Found heuristic solution: objective 8901.9241439
#> Presolve time: 0.14s
#> Presolved: 5325 rows, 2880 columns, 17280 nonzeros
#> Variable types: 0 continuous, 2880 integer (2880 binary)
#> Presolve removed 870 rows and 0 columns
#> Presolved: 4455 rows, 2880 columns, 14670 nonzeros
#>
#>
#> Root relaxation: objective 4.633285e+03, 2617 iterations, 0.23 seconds
#>
#>     Nodes    |    Current Node    |     Objective Bounds      |     Work
#>  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
#>
#>      0     0 4633.28464    0 1347 8901.92414 4633.28464  48.0%     -    0s
#> H    0     0                    8471.3546786 4633.28464  45.3%     -    0s
#> H    0     0                    6981.3189782 4633.28464  33.6%     -    0s
#>      0     0 4656.86431    0 1313 6981.31898 4656.86431  33.3%     -    1s
#>      0     0 4658.09598    0 1312 6981.31898 4658.09598  33.3%     -    1s
#>      0     0 4658.55762    0 1298 6981.31898 4658.55762  33.3%     -    1s
#>      0     0 4658.89951    0 1298 6981.31898 4658.89951  33.3%     -    1s
#>      0     0 4659.06805    0 1298 6981.31898 4659.06805  33.3%     -    1s
#>      0     0 4668.49591    0 1313 6981.31898 4668.49591  33.1%     -    1s
#>      0     0 4669.27758    0 1310 6981.31898 4669.27758  33.1%     -    1s
#>      0     0 4669.27758    0 1310 6981.31898 4669.27758  33.1%     -    1s
#>      0     0 4670.12210    0 1309 6981.31898 4670.12210  33.1%     -    1s
#>      0     0 4674.25742    0 1305 6981.31898 4674.25742  33.0%     -    1s
#>      0     0 4674.70511    0 1307 6981.31898 4674.70511  33.0%     -    1s
#>      0     0 4675.43856    0 1306 6981.31898 4675.43856  33.0%     -    1s
#>      0     0 4675.43856    0 1306 6981.31898 4675.43856  33.0%     -    5s
#> H    0     0                    6144.1824440 4675.43856  23.9%     -    6s
#> H    0     0                    6137.5371003 4675.43856  23.8%     -    6s
#>      0     0 4877.80932    0  773 6137.53710 4877.80932  20.5%     -    6s
#> H    0     0                    5731.2082580 4877.80932  14.9%     -    7s
#>      0     0 4928.81742    0  844 5731.20826 4928.81742  14.0%     -    7s
#>      0     0 4932.02113    0  857 5731.20826 4932.02113  13.9%     -    7s
#>      0     0 4932.05034    0  856 5731.20826 4932.05034  13.9%     -    7s
#>      0     0 4932.05034    0  853 5731.20826 4932.05034  13.9%     -    8s
#>      0     2 4932.24916    0  853 5731.20826 4932.24916  13.9%     -    9s
#>      4     6 5036.12183    2  871 5731.20826 4969.46126  13.3%   450   10s
#> H   27    27                    5667.6160091 5026.81145  11.3%   293   12s
#>    193   166 5213.24245    6  827 5667.61601 5036.27841  11.1%   108   15s
#>    378   328     cutoff   39      5667.61601 5037.64301  11.1%   115   20s
#>    513   441 5356.42776   14  710 5667.61601 5047.35478  10.9%   105   26s
#>    520   446 5560.82950   23  766 5667.61601 5047.35478  10.9%   104   30s
#>    535   458 5221.31926   22  679 5667.61601 5047.35478  10.9%   112   35s
#>    676   530 5378.58243   30  607 5667.61601 5076.55144  10.4%   123   40s
#>    842   616 5566.49748   37  740 5667.61601 5082.67839  10.3%   133   45s
#>
#> Cutting planes:
#>   Gomory: 9
#>
#> Explored 1022 nodes (141332 simplex iterations) in 49.73 seconds
#> Thread count was 1 (of 4 available processors)
#>
#> Solution count 7: 5667.62 5731.21 6137.54 ... 8901.92
#>
#> Optimal solution found (tolerance 1.00e-01)
#> Best objective 5.667616009136e+03, best bound 5.102159870824e+03, gap 9.9770%
#> Optimize a model with 5325 rows, 2880 columns and 12060 nonzeros
#> Variable types: 0 continuous, 2880 integer (2880 binary)
#> Coefficient statistics:
#>   Matrix range     [2e-01, 1e+00]
#>   Objective range  [5e-01, 2e+02]
#>   Bounds range     [1e+00, 1e+00]
#>   RHS range        [1e+00, 1e+01]
#> Found heuristic solution: objective 10118.924144
#> Presolve time: 0.13s
#> Presolved: 5325 rows, 2880 columns, 17280 nonzeros
#> Variable types: 0 continuous, 2880 integer (2880 binary)
#> Presolve removed 870 rows and 0 columns
#> Presolved: 4455 rows, 2880 columns, 14670 nonzeros
#>
#>
#> Root relaxation: objective 8.569664e+03, 2956 iterations, 0.28 seconds
#>
#>     Nodes    |    Current Node    |     Objective Bounds      |     Work
#>  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
#>
#>      0     0 8569.66432    0  836 10118.9241 8569.66432  15.3%     -    0s
#> H    0     0                    9381.9905191 8569.66432  8.66%     -    0s
#>
#> Explored 1 nodes (3447 simplex iterations) in 0.50 seconds
#> Thread count was 1 (of 4 available processors)
#>
#> Solution count 2: 9381.99 10118.9
#>
#> Optimal solution found (tolerance 1.00e-01)
#> Best objective 9.381990519089e+03, best bound 8.569664324288e+03, gap 8.6584%
#> Optimize a model with 5325 rows, 2880 columns and 12060 nonzeros
#> Variable types: 0 continuous, 2880 integer (2880 binary)
#> Coefficient statistics:
#>   Matrix range     [2e-01, 1e+00]
#>   Objective range  [3e+00, 2e+02]
#>   Bounds range     [1e+00, 1e+00]
#>   RHS range        [1e+00, 1e+01]
#> Found heuristic solution: objective 9333.9241439
#> Presolve time: 0.13s
#> Presolved: 5325 rows, 2880 columns, 17280 nonzeros
#> Variable types: 0 continuous, 2880 integer (2880 binary)
#> Presolve removed 870 rows and 0 columns
#> Presolved: 4455 rows, 2880 columns, 14670 nonzeros
#>
#>
#> Root relaxation: objective 6.027838e+03, 3416 iterations, 0.37 seconds
#>
#>     Nodes    |    Current Node    |     Objective Bounds      |     Work
#>  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
#>
#>      0     0 6027.83792    0 1239 9333.92414 6027.83792  35.4%     -    0s
#> H    0     0                    8316.1074209 6027.83792  27.5%     -    0s
#> H    0     0                    7943.5078917 6027.83792  24.1%     -    1s
#>      0     0 6085.01523    0 1291 7943.50789 6085.01523  23.4%     -    1s
#>      0     0 6085.89684    0 1288 7943.50789 6085.89684  23.4%     -    1s
#>      0     0 6115.17402    0 1368 7943.50789 6115.17402  23.0%     -    1s
#>      0     0 6117.59053    0 1391 7943.50789 6117.59053  23.0%     -    2s
#>      0     0 6117.59053    0 1392 7943.50789 6117.59053  23.0%     -    2s
#>      0     0 6125.28918    0 1396 7943.50789 6125.28918  22.9%     -    2s
#>      0     0 6125.30907    0 1425 7943.50789 6125.30907  22.9%     -    2s
#>      0     0 6125.32346    0 1425 7943.50789 6125.32346  22.9%     -    2s
#>      0     0 6125.32346    0 1423 7943.50789 6125.32346  22.9%     -    2s
#>      0     0 6125.37767    0 1420 7943.50789 6125.37767  22.9%     -    3s
#>      0     0 6125.37767    0 1403 7943.50789 6125.37767  22.9%     -    3s
#> H    0     0                    7601.8938452 6125.37767  19.4%     -    4s
#>      0     0 6173.38320    0  842 7601.89385 6173.38320  18.8%     -    4s
#> H    0     0                    7589.6868875 6173.38320  18.7%     -    5s
#>      0     0 6287.27243    0  928 7589.68689 6287.27243  17.2%     -    5s
#>      0     0 6287.71486    0  942 7589.68689 6287.71486  17.2%     -    5s
#>      0     0 6302.24509    0  909 7589.68689 6302.24509  17.0%     -    5s
#>      0     0 6302.24509    0  907 7589.68689 6302.24509  17.0%     -    8s
#>      0     2 6302.28396    0  906 7589.68689 6302.28396  17.0%     -    9s
#>      1     3 6361.29595    1 1001 7589.68689 6302.28983  17.0%   756   10s
#> H   27    27                    7574.0270291 6344.63600  16.2%   312   14s
#>     31    33 6530.42577   18  926 7574.02703 6344.63600  16.2%   319   15s
#> H   54    54                    7556.4677132 6344.63600  16.0%   285   16s
#> H   81    81                    7435.9445339 6344.63600  14.7%   252   16s
#>    251   237 7329.21229   55  527 7435.94453 6344.64453  14.7%   170   20s
#>    509   465 7114.79320   29  822 7435.94453 6360.19390  14.5%   154   25s
#>
#> Explored 512 nodes (92924 simplex iterations) in 27.65 seconds
#> Thread count was 1 (of 4 available processors)
#>
#> Solution count 8: 7435.94 7556.47 7574.03 ... 9333.92
#>
#> Optimal solution found (tolerance 1.00e-01)
#> Best objective 7.435944533885e+03, best bound 7.179092652434e+03, gap 3.4542%
#> Optimize a model with 6485 rows, 2880 columns and 15540 nonzeros
#> Variable types: 0 continuous, 2880 integer (2880 binary)
#> Coefficient statistics:
#>   Matrix range     [2e-01, 1e+00]
#>   Objective range  [2e+00, 2e+02]
#>   Bounds range     [1e+00, 1e+00]
#>   RHS range        [1e+00, 1e+01]
#> Found heuristic solution: objective 10084.071885
#> Presolve time: 0.14s
#> Presolved: 6485 rows, 2880 columns, 21920 nonzeros
#> Variable types: 0 continuous, 2880 integer (2880 binary)
#> Presolve removed 870 rows and 0 columns
#> Presolved: 5615 rows, 2880 columns, 19310 nonzeros
#>
#>
#> Root relaxation: objective 8.178938e+03, 1514 iterations, 0.22 seconds
#>
#>     Nodes    |    Current Node    |     Objective Bounds      |     Work
#>  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
#>
#>      0     0 8178.93759    0  960 10084.0719 8178.93759  18.9%     -    0s
#> H    0     0                    9010.5140391 8178.93759  9.23%     -    0s
#>
#> Explored 1 nodes (1514 simplex iterations) in 0.46 seconds
#> Thread count was 1 (of 4 available processors)
#>
#> Solution count 2: 9010.51 10084.1
#>
#> Optimal solution found (tolerance 1.00e-01)
#> Best objective 9.010514039144e+03, best bound 8.178937587478e+03, gap 9.2290%
#> Optimize a model with 6485 rows, 2880 columns and 15540 nonzeros
#> Variable types: 0 continuous, 2880 integer (2880 binary)
#> Coefficient statistics:
#>   Matrix range     [2e-01, 1e+00]
#>   Objective range  [2e+01, 2e+02]
#>   Bounds range     [1e+00, 1e+00]
#>   RHS range        [1e+00, 1e+01]
#> Found heuristic solution: objective 8859.0718845
#> Presolve time: 0.14s
#> Presolved: 6485 rows, 2880 columns, 21920 nonzeros
#> Variable types: 0 continuous, 2880 integer (2880 binary)
#> Presolve removed 870 rows and 0 columns
#> Presolved: 5615 rows, 2880 columns, 19310 nonzeros
#>
#>
#> Root relaxation: objective 3.483401e+03, 1096 iterations, 0.12 seconds
#>
#>     Nodes    |    Current Node    |     Objective Bounds      |     Work
#>  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
#>
#>      0     0 3483.40090    0 1372 8859.07188 3483.40090  60.7%     -    0s
#> H    0     0                    7810.0906779 3483.40090  55.4%     -    0s
#> H    0     0                    7064.2663330 3483.40090  50.7%     -    0s
#>      0     0 3503.81769    0 1355 7064.26633 3503.81769  50.4%     -    0s
#>      0     0 3523.48949    0 1357 7064.26633 3523.48949  50.1%     -    3s
#>      0     0 3523.48949    0 1357 7064.26633 3523.48949  50.1%     -    3s
#> H    0     0                    6356.2198076 3523.48949  44.6%     -    4s
#>      0     0 4662.41137    0  963 6356.21981 4662.41137  26.6%     -    5s
#> H    0     0                    5923.3227102 4662.41137  21.3%     -    5s
#>      0     0 4689.06065    0  978 5923.32271 4689.06065  20.8%     -    5s
#>      0     0 4699.05361    0  983 5923.32271 4699.05361  20.7%     -    5s
#>      0     0 4703.72536    0  975 5923.32271 4703.72536  20.6%     -    5s
#>      0     0 4704.30769    0  952 5923.32271 4704.30769  20.6%     -    6s
#>      0     0 4705.83010    0  951 5923.32271 4705.83010  20.6%     -    6s
#>      0     0 4705.83010    0  954 5923.32271 4705.83010  20.6%     -    6s
#>      0     0 4705.87951    0  942 5923.32271 4705.87951  20.6%     -    6s
#>      0     0 4707.56395    0  940 5923.32271 4707.56395  20.5%     -    6s
#>      0     0 4707.56395    0  944 5923.32271 4707.56395  20.5%     -    6s
#>      0     0 4707.93248    0  933 5923.32271 4707.93248  20.5%     -    6s
#>      0     0 4707.93248    0  940 5923.32271 4707.93248  20.5%     -    7s
#>      0     0 4708.05202    0  920 5923.32271 4708.05202  20.5%     -    7s
#>      0     0 4708.08292    0  949 5923.32271 4708.08292  20.5%     -    7s
#>      0     0 4708.08292    0  952 5923.32271 4708.08292  20.5%     -    7s
#>      0     0 4708.08292    0  952 5923.32271 4708.08292  20.5%     -    7s
#>      0     0 4708.08292    0  952 5923.32271 4708.08292  20.5%     -    7s
#>      0     2 4708.22379    0  952 5923.32271 4708.22379  20.5%     -    8s
#>     38    40 5016.93254   23  833 5923.32271 4808.95076  18.8%   269   10s
#> H   78    78                    5824.7673625 4808.95076  17.4%   240   11s
#>    248   227 5720.68039   57  377 5824.76736 4809.11728  17.4%   159   15s
#> H  349   306                    5810.0856611 4809.77546  17.2%   139   16s
#>    470   409 5142.64596   15  888 5810.08566 4811.51971  17.2%   153   20s
#> H  481   356                    5744.3666214 4811.51971  16.2%   154   20s
#> H  512   368                    5702.3026197 4897.30484  14.1%   155   24s
#>    513   368 5103.93582   13  710 5702.30262 4897.30484  14.1%   155   25s
#>    532   381 5676.42444   47  747 5702.30262 4979.52135  12.7%   149   30s
#> H  562   382                    5657.4025993 5080.45682  10.2%   161   33s
#> H  618   376                    5624.6902852 5081.02753  9.67%   155   35s
#>
#> Cutting planes:
#>   Gomory: 11
#>   Zero half: 2
#>
#> Explored 618 nodes (100451 simplex iterations) in 35.19 seconds
#> Thread count was 1 (of 4 available processors)
#>
#> Solution count 10: 5624.69 5657.4 5702.3 ... 7810.09
#>
#> Optimal solution found (tolerance 1.00e-01)
#> Best objective 5.624690285211e+03, best bound 5.081027527176e+03, gap 9.6656%s4 <- lapply(s4, category_layer)
s4 <- stack(s4)

# plot solutions
plot(s4, main = names(p4), axes = FALSE, box = FALSE)
# create an array to manually specify the connectivities between
# each planning unit when they are allocated to each different zone
# for real-world problems, these connectivities would be generated using
# data - but here these connectivity values are assigned as random
# ones or zeros
c_array <- array(0, c(rep(ncell(sim_pu_zones_stack[[1]]), 2), 3, 3))
for (z1 in seq_len(3))
for (z2 in seq_len(3))
c_array[, , z1, z2] <- round(runif(ncell(sim_pu_zones_stack[[1]]) ^ 2,
0, 0.505))

# create a problem with the manually specified connectivity array
# note that the zones argument is set to NULL because the connectivity
# data is an array
p5 <- list(p3,
p3 %>% add_connectivity_penalties(30, zones = NULL, c_array))

# assign names to the problems
names(p5) <- c("basic problem", "connectivity array")# solve problems
s5 <- lapply(p5, solve)#> Optimize a model with 105 rows, 270 columns and 1620 nonzeros
#> 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 10275.924144
#> Presolve time: 0.01s
#> Presolved: 105 rows, 270 columns, 1620 nonzeros
#> Variable types: 0 continuous, 270 integer (270 binary)
#> Presolved: 105 rows, 270 columns, 1620 nonzeros
#>
#>
#> Root relaxation: objective 8.943504e+03, 123 iterations, 0.00 seconds
#>
#>     Nodes    |    Current Node    |     Objective Bounds      |     Work
#>  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
#>
#>      0     0 8943.50352    0   15 10275.9241 8943.50352  13.0%     -    0s
#> H    0     0                    9111.3117236 8943.50352  1.84%     -    0s
#>
#> Explored 1 nodes (123 simplex iterations) in 0.01 seconds
#> Thread count was 1 (of 4 available processors)
#>
#> Solution count 2: 9111.31 10275.9
#>
#> Optimal solution found (tolerance 1.00e-01)
#> Best objective 9.111311723633e+03, best bound 8.943503517801e+03, gap 1.8418%
#> Optimize a model with 1479 rows, 957 columns and 4368 nonzeros
#> Variable types: 0 continuous, 957 integer (957 binary)
#> Coefficient statistics:
#>   Matrix range     [2e-01, 1e+00]
#>   Objective range  [3e+01, 2e+02]
#>   Bounds range     [1e+00, 1e+00]
#>   RHS range        [1e+00, 1e+01]
#> Found heuristic solution: objective 9765.9241439
#> Presolve time: 0.03s
#> Presolved: 1479 rows, 957 columns, 4416 nonzeros
#> Variable types: 0 continuous, 957 integer (957 binary)
#> Presolve removed 24 rows and 0 columns
#> Presolved: 1455 rows, 957 columns, 4344 nonzeros
#>
#>
#> Root relaxation: objective 5.800862e+03, 856 iterations, 0.07 seconds
#>
#>     Nodes    |    Current Node    |     Objective Bounds      |     Work
#>  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
#>
#>      0     0 5800.86247    0  718 9765.92414 5800.86247  40.6%     -    0s
#> H    0     0                    8273.8383046 5800.86247  29.9%     -    0s
#>      0     0 5811.11011    0  710 8273.83830 5811.11011  29.8%     -    1s
#>      0     0 5849.13104    0  760 8273.83830 5849.13104  29.3%     -    2s
#>      0     0 5849.13104    0  760 8273.83830 5849.13104  29.3%     -    2s
#> H    0     0                    7450.1953151 5849.13104  21.5%     -    2s
#> H    0     0                    7254.3989724 5849.13104  19.4%     -    3s
#>      0     2 5849.14013    0  760 7254.39897 5849.14013  19.4%     -    3s
#>    201   179 7085.32194   56  396 7254.39897 5903.38490  18.6%   104    5s
#> H  516   427                    7247.7170296 5979.83470  17.5%  89.6    9s
#>    586   476 6860.44564   44  326 7247.71703 5979.83470  17.5%  98.1   10s
#> H  979   539                    7065.5819604 5979.83470  15.4%  90.2   13s
#>   1065   585 6205.14691   19  560 7065.58196 5979.83470  15.4%  93.1   15s
#>   1368   761 6494.45306   25  589 7065.58196 5979.83470  15.4%   105   20s
#> H 1659   684                    6975.8820893 6021.95071  13.7%   103   23s
#>   1830   825 6739.73594   37  495 6975.88209 6022.72638  13.7%   107   25s
#>   2098  1018 6057.21237   13  628 6975.88209 6053.31604  13.2%   106   30s
#>   2283  1167 6291.70924   20  652 6975.88209 6056.64380  13.2%   107   35s
#>   2478  1329 6351.43886   22  609 6975.88209 6071.92695  13.0%   109   40s
#> H 2810  1572                    6975.8818960 6075.25248  12.9%   106   42s
#> H 2919  1663                    6975.8818689 6077.66578  12.9%   107   43s
#> H 2958  1387                    6862.9838396 6077.66578  11.4%   107   45s
#> H 2960  1319                    6832.0786795 6077.66578  11.0%   107   45s
#> H 3258  1446                    6798.3350632 6087.63386  10.5%   108   48s
#> H 3316  1470                    6794.8824564 6106.12515  10.1%   109   49s
#>   3376  1530     cutoff   38      6794.88246 6106.12515  10.1%   109   50s
#> H 3514  1594                    6787.0663425 6110.18031  10.0%   110   51s
#>
#> Cutting planes:
#>   Gomory: 2
#>   Zero half: 12
#>
#> Explored 3514 nodes (386522 simplex iterations) in 51.60 seconds
#> Thread count was 1 (of 4 available processors)
#>
#> Solution count 10: 6787.07 6794.88 6798.34 ... 7247.72
#>
#> Optimal solution found (tolerance 1.00e-01)
#> Best objective 6.787066342525e+03, best bound 6.110180313460e+03, gap 9.9732%s5 <- lapply(s5, category_layer)
s5 <- stack(s5)

# plot solutions
plot(s5, main = names(p5), axes = FALSE, box = FALSE)