Add max area threshold for zone and dealing with fragmented zones

[Tegan151]

Hi Jeff,

I have a multizone problem that I am trying to run. The conservation zone tends to dominate the solution, which is great, in my view, but is not practical for implementation purposes. Is there a way to set a max area threshold for the zone? So for example, the government says that they are only willing to implement a 50 hectare protected area, can I build this into the prioritizr problem?

Secondly, most of my zones are nicely congregated together but I do sometimes have a random planning unit assigned to a zone in the middle of another zone. I have added a boundary and connectivity penalties to try solve this but I still have the same problem. Is there a way around this? Also, what is the connectivity penalty value in this function? Is there any way to meaning select this value? What is the different, for example, between a connectivity penalty value of 1 and 50?

Warmest,
Tegan

[jeffreyhanson]

Hi @Tegan151,

Awesome! Yeah, to set a threshold for the maximum area selected, you can use the linear constraints (see https://prioritizr.net/reference/add_linear_constraints.html). Note that if you're using a min set problem formulation (menaing that the targets must be met), then you might need to switch to a budget limited objective if you get errors about problem feasibility (e.g. min shortfall objective, https://prioritizr.net/reference/add_min_shortfall_objective.html). These problem feasbility errors (in this particular case) mean that there's no way to meet the target whilst also meeting the maximum area constraint. How does that sound?

Yeah, may I ask if you rescaled the boundary data for the prioritization? Although we recommended using the scales::rescale() function for previous versions of prioritizr, the scales::rescale() function will produce this strange behavior for the more recent versions of prioritizr (due to recent change in how prioritizr handles the boundary data internally). Could you try using the rescale_matrix() function to rescale the boundary data prior to optimization and see if that fixes it (see https://prioritizr.net/reference/rescale_matrix.html for details)?

The penalty values basically serve as trade-off parameters between the primary objective (e.g. total cost when using the min set objective) and a secondary objective (i.e. the boundary or connectivity penalites). When you use a higher penalty value, then you're saying that it's more important to focus in the secondary objective than the first objective. Unfortunately, the penalty value don't really have a meaning (they're unitless) - so this can make it hard to find a suitable penalty value. To help with selecting a "good" penalty value, we've written up a tutorial showing some different strategies (see https://prioritizr.net/articles/calibrating_trade-offs_tutorial.html). Does that help?

Let me know if you have any further questions?

[Tegan151]

Hi Jeff,

Thanks for the quick reply.

I understand the linear constraints and have used it to make sure thresholds of certain features are met within different zones but I am still unsure how to tell the problem that I would only like the conservation zone to cover a max of 50 planning units, for example, so that our new suggested MPA does not take up too much space in my solution. I have pasted a figure of my most recent solution to help you visualize.

I did not rescale the boundary data. I will give this a go. Thank you.

I did have a look through your calibration tutorial at one stage but I found it difficult to implement for a multi-zone problem. I will give another go and let you know if I have any problems.

Thanks again for your help reply and and assistance.

[jeffreyhanson]

No worries! Yeah, I can put togeather an example showing how to use the linear constraints - what data format are your planning units in (raster, shapefile)? I can write the example so it uses the same data format to make adapting it easier.

Ah, yeah, the tutorial only covers single zones, so let me know if you run into any issues. Broadly speaking, I would first create a zones matrix with the right zone relationships, and then then use a calibration procedure to find a good penalty value (with zones matrix staying the same for every run in the calibration procedure). To assess the total boundary length or connectivity for a given run in the calibration procedure, I would recommend using the "overall" statistic (e.g. computed via eval_boundary_summary() or eval_connectivity_summary()). After finding a "good" penalty value, I would then see if the spatial clustering of zones looks good. If it doesn't I would then try updating the zones matrix, and then redoing the calibration procedure again. This might take some iterations to find a good zones matrix and penalty values. Another option could be to try generating multiple different prioritizations based on different zones matrix that represent different scenarios. How does that sound?

[jeffreyhanson]

Also, just out of curiosity, is that the Bay of Plenty in New Zealand?

[Tegan151]

Thanks, Jeff! You have been so great.

No, this is Algoa Bay in South Africa. The 'bay of plenty' is what we have named the scenario that is 'fair' to humans and biodiversity by setting targets and cost layers the same way for each feature/use and zone. We will also have scenario called 'eerie bay' (human use focused) and 'deserted bay' (biodiversity focused).

I am using shapefiles. I have paste my problem below to show you what I have done.

p0 <-
problem(
polygon,
zones(
names(p1),
names(p2),
names(p3),
names(p4),
names(p5),
names(p6),
names(p7),

  zone_names = c("ConsRes", "ConsCon",'CmmNnEx','CmmExtr','CmmFshr','Trnsprt','ADZ')
),
cost_column = c("ConsRes", "ConsCon",'CmmNnEx','CmmExtr','CmmFshr','Trnsprt','ADZ') ) %>%

add_locked_in_constraints(c("lock_in_conres",'lock_in_concon','lock_in_CmmNnEx','lock_in_CmmExtr','lock_in_CmmFshr','lock_in_Trnsprt','lock_in_ADZ')) %>%
add_locked_out_constraints(c('lock_out_conres','lock_out_concon','lock_out_CmmNnEx','lock_out_CmmExtr','lock_out_CmmFshr','lock_out_Trnsprt','lock_out_ADZ')) %>%
add_min_set_objective()%>%
add_relative_targets(targets) %>%
add_connectivity_penalties(1, zm, d_matrix)%>%
add_mandatory_allocation_constraints()%>%
add_boundary_penalties(penalty = 0.01, zones=zm6,edge_factor= rep(0.5,7))%>%
add_gurobi_solver(numeric_focus = T,gap=0.05,threads = parallel::detectCores(TRUE) - 1)

[jeffreyhanson]

Ok, here's a reproducible example.

Let's start by creating a simple multi-zone problem. Here we will use the minimum shortfall objective. This means our objective is to get as close to possible as meeting the targets, whilst ensuring that a particular budget is not exceeded. Note that we are using the minimum shortfall objective because it is a type of budget-limited objective -- meaning that the optimization problem will have feasible solution even if all of the targets cannot be met -- which is useful when we later add in linear constraints. Also, we will set a budget equal to 15% of the total cost across all zones.

# load packages
library(prioritizr)
library(sf)

# import data
pu <- get_sim_zones_pu_polygons()
fts <- get_sim_zones_features()

# define targets
targets <- matrix(
  runif(number_of_features(fts) * number_of_zones(fts), 0.2, 0.3),
  nrow = number_of_features(fts),
  ncol = number_of_zones(fts),
  dimnames = list(
    feature_names(fts), 
    zone_names(fts)
  )
)

# define budgets
budget <- 0.15 * sum(pu$cost_1, pu$cost_2, pu$cost_3)

# create minimal multi-zone problem
p1 <-
  problem(
    pu, fts,
    cost_column = c("cost_1", "cost_2", "cost_3")
  ) %>%
  add_min_shortfall_objective(budget = budget) %>%
  add_relative_targets(targets) %>%
  add_binary_decisions() %>%
  add_default_solver(verbose = FALSE)

# solve problem
s1 <- solve(p1)

# plot solution
plot(s1[, c("solution_1_zone_1", "solution_1_zone_2", "solution_1_zone_3")])

We can see that many planning units were selected in the first zone. Let's pretend that we have stakeholders who say that no more than 10 planning units should be allocated to the first zone. We can achieve this using linear constraints. To specify the constraints, we will add some columns to the planning unit data (one column per zone) that have a value of 1 for the first zone (because that's the zone we want to control) and a value of 0 for the remaining zones (because we don't want this constraint to affect the other zones). Note that if we wanted to set maximum limit for each zone, then we would need to add a seperate constraint for each zone. Alternatively, if we wanted to set maximum limits for combinations of zones, then we would set a value of 1 in multiple columns (e.g., if we wanted to ensure that the total number of planning units selected in zones 1 and zone 3 should not exceed 12 planning units, then we would have 1s in the columns for zones 1 and 3, and 0s in the column for zone 2).

# add columns to planning unit data for the linear constraint
pu$lc_zone_1 <- 1
pu$lc_zone_2 <- 0
pu$lc_zone_3 <- 0

# create multi-zone problem with added linear constraint
p2 <-
  problem(
    pu, fts,
    cost_column = c("cost_1", "cost_2", "cost_3")
  ) %>%
  add_min_shortfall_objective(budget = budget) %>%
  add_relative_targets(targets) %>%
  add_linear_constraints(
    data = c("lc_zone_1", "lc_zone_2", "lc_zone_3"),
    sense = "<=",
    threshold = 10
  ) %>%
  add_binary_decisions() %>%
  add_default_solver(verbose = FALSE)

# solve problem
s2 <- solve(p2)

# plot solution
plot(s2[, c("solution_1_zone_1", "solution_1_zone_2", "solution_1_zone_3")])

How does that sound? Let me know if you have any questions?

[jeffreyhanson]

Re scenarios: Oh I see! Thanks for explaining that. Yeah, those pretty good names for scenarios.

[Tegan151]

Thanks Jeff! I will run over the weekend and let you know if I have any questions on Monday.

[Tegan151]

Hi Jeff,

Your code worked like a dream! Thank you.

I am still getting random block in the middle of zones so I have added a connectivity penalty in addition to the boundary penalty using the following code:
centroids <- sf::st_coordinates(
suppressWarnings(sf::st_centroid(polygon))
)
d_matrix <- (1 / (Matrix::Matrix(as.matrix(dist(centroids))) + 1))

rescale matrix values to have a maximum value of 1

d_matrix <- rescale_matrix(d_matrix, max = 1)

And the problem:
p0 <-
problem(
polygon,
zones(
names(p1),
names(p2),
names(p3),
names(p4),
names(p5),
names(p6),
names(p7),

  zone_names = c("ConsRes", "ConsCon",'CmmNnEx','CmmExtr','CmmFshr','Trnsprt','ADZ')
),
cost_column = c("ConsRes", "ConsCon",'CmmNnEx','CmmExtr','CmmFshr','Trnsprt','ADZ') ) %>%

add_locked_in_constraints(c("lock_in_conres",'lock_in_concon','lock_in_CmmNnEx','lock_in_CmmExtr','lock_in_CmmFshr','lock_in_Trnsprt','lock_in_ADZ')) %>%
add_locked_out_constraints(c('lock_out_conres','lock_out_concon','lock_out_CmmNnEx','lock_out_CmmExtr','lock_out_CmmFshr','lock_out_Trnsprt','lock_out_ADZ')) %>%
add_min_set_objective()%>%
add_relative_targets(targets) %>%
add_connectivity_penalties(1, zm, d_matrix)%>%
add_linear_constraints(
data = c("lc_ConsRes", "lc_ConsCon",'lc_CmmNnEx','lc_CmmExtr','lc_CmmFshr','lc_Trnsprt','lc_ADZ'),
sense = "<=",
threshold = 3000
) %>%
add_mandatory_allocation_constraints()%>%
add_boundary_penalties(penalty = 0.01, zones=zm6,edge_factor= rep(0.5,7),data =bm_ply)%>%
add_gurobi_solver(numeric_focus = T,threads = parallel::detectCores(TRUE) - 1)

Is that the best way to add this penalty?

Thank you again for all your help! Is has been incredibly useful.

[jeffreyhanson]

Yeah that looks good to me! Note that having both connectivity and boundary penalties will slow things down a fair bit though.

Also, here's some other ideas for the random planning unit issues:

1. It looks like the d_matrix object has non-zero values for all pairs of planning units? If so, I'd suggest applying a threshold (e.g using a quantile approach) to basically set all "connectivity" values below a certain threshold to zero. This could potentially help with the random planning unit issues too, because currently the connectivity data is saying that even really distant planning units have some non-zero connectivity between them. So, e.g., you could apply a 50% quantile threshold using something like this:

d_matrix <- Matrix::drop0(d_matrix)   # ensure sparse matrix format has non-zeros
d_threshold <- quantile(d_matrix, probs = 0.5, names = FALSE)  # calculate threshold
d_matrix@x[d_matrix@x < d_threshold] <- 0 # apply threshold
d_matrix <- Matrix::drop0(d_matrix) # ensure sparse matrix format has non-zeros

2. Another strategy could be playing with the zones matrix argument for add_connectivity_penalties(). By setting off-diagonal cell values to -1, this will discourage the solution from allocating nearby planning units to different zones. Note that you'll still want cell values on the diagonal to be 1 to promote connectivity within the same zone. You could create such a matrix using something like this:

zm <- matrix(0, ncol = 7, nrow = 7)  # initialize with zeros
zm[] <- -1 # set all values to -1
diag(zm) <- 1 # set diagonal values to 1

How does that sound?

[Tegan151]

Thanks, Jeff! I will definitely give it a go and let you know how I get on.

I am leaving today for a 3 week work trip so I will be slow to response.

[jeffreyhanson]

@Tegan151, I just wanted to follow up and ask if those suggestions were able to resolve the issue? Also, did you have any further questions?

[Tegan151]

Hi Jeff,

Sorry I have been so slow to get back to you. No - the issue is not quite solved yet. I went back and reformatted and simplified all of my data so that the problem would run faster. This took a bit of time.

I am currently trying a few things to get the zones to cluster better. I will let you know how I get on.

[jeffreyhanson]

No worries! Sorry to hear it's not solved yet. Sounds good, let me know if I can help with any further questions.

[Tegan151]

Hi Jeff,
Sorry for taking so long to get back to you. I am still not able to get rid of those random red blocks in the middle of the solutions. See some examples below. I have changed to using the adjacency matrix as the problems solve much faster than the 'd_matrix'. I have tried both suggestion you gave on 8 Nov. Do you have any other suggestions? Otherwise, I can justify the random blocks as this is a decision support tools and not for final zonation. I can also show what is driving the selection of those blocks..

[jeffreyhanson]

Hi @Tegan151,

No worries! I'm sorry those suggestions weren't able to fix this. Could you remind me which objective function you're using?

If you're using min set, then I wonder if it might be useful to see if there's any features in those random red blocks that are causing them to be selected? E.g., if you have some features that only occur in the middle part of the study area and you have targets for the commercial fisheries zone saying that you need a certain percent of each feature, then the prioritization will have to prioritize some planning units in the middle of the study area for fisheries to meet the targets (no matter how you change the boundary or connectivity penalties)? After identifying which features are causing the planning units in the middle to be selected, you could then reduce the targets.

If you're using min shortfall, then I would expect increasing the connectivity penalty parameter to get rid of them. Unless, there's a mistake with any locked in constraints or linear constraints that say the commerical fisheries zone must select a planning unit in the middle of the study area? Maybe it would be worth double checking the input data for those constraints?

Also, how close is the total area of planning units allocated to the community non-extractive zone to its maximum area constraint? Is it possible that the prioritization is allocating those planning units to the commerical fisheries zone, because the maximum area constraint is preventing the prioritization from allocating all those planning units to the community non-extractive zone?

How does that sound?

[Tegan151]

Hi Jeff,

Thanks again for quick response!

I have checked my feature layers, and it is one of the commercial fisheries layers that pretty much only occurs in those blocks. Well, in a very high intensity in those blocks and much lower intensities in all other planning units. I will take this to the rest of the team and see what they think.

I am using the min_set_objective.

Thanks so much for all off your patience. I really have learnt so much in these interactions with you.

One last question, for publication, to show how I selected the connectivity and boundary penalty factors., could I just show the cost of solution and number of planning units selected per zone? Something like below? I know you have a calibrating tradeoffs tutorial on your website but it hasn't been written for zones.

Warmest regards,
Tegan

[jeffreyhanson]

No worries!

Ah - yeah that would definitely explain why the penalties aren't affecting the red blocks. In addition to changing the fisheries target for this feature, another option could be to try using the locked out constraints to specify that none of the planning units in the middle of the planning region should be assigned to the fisheries zone. However, if you're using percentage-based targets, this would probably lead to problem infeasibility (where there is no valid solution that meets all the targets).

Yeah, a figure like that is perfectly reasonable. Many studies don't bother with a proper calibration exercise and simply say something like "Preliminary analyses were used to identify a suitable boundary length penalty factor." - so a figure like that is already doing better than some published papers. Also, since you're using both boundary and connectiivity penalties, running a full calibration process could be a bit tricky.

If you're interested, the code for a boundary length calibration for multiple zones is pretty similar to that described in the tutorial. The main difference is that when you calculate the boundary length, you need to extract the overall boundary length. E.g., using something like this:

# notes on input data
## p0 is the problem
## x is the solution
## z is the zones argument used for the penalties 
## d is the boundary data used for the penalties

# here is the code for a single zone problem
total_boundary_length <- eval_boundary_summary(p, x, data = d)$boundary

# here is the equivalent code for multi-zone problem
total_boundary_length <- 
  eval_boundary_summary(p, x, zones = z, data = d) %>% 
  dplyr::filter(summary == "overall") %>% 
  dplyr::pull(boundary)

How does that sound?

[Tegan151]

Hi Jeff,

Thanks again for your reply.

I have written some code following the hierarchical approach tutorial. It is taking forever to run, so no results yet. Would you please double check it for me? It was a bit fiddley.

#calibrating boundary penalty
rm(list=ls())
polygon<- readRDS('polygon.rds')

pp1 <- readRDS('p1.rds')
pp2 <- readRDS('p2.rds')
pp3 <- readRDS('p3.rds')
pp4 <- readRDS('p4.rds')
pp5 <- readRDS('p5.rds')
pp6 <- readRDS('p6.rds')

targets <- readRDS('targets.rds')
feature <- gsub('_ConsRes','',names(pp1))

zone_features <- zones(names(pp1),names(pp2),names(pp3),names(pp4),names(pp5),names(pp6), zone_names = c("ConsRes", "ConsCon",'CmmNnEx','CmmExtr','CmmFshr','Trnsprt'),feature_names = feature)

#Set up penalty factor
bm_ply <- boundary_matrix(polygon)
bm_ply <- rescale_matrix(bm_ply)
zm6 <- diag(6)
colnames(zm6) <- c("ConsRes", "ConsCon",'CmmNnEx','CmmExtr','CmmFshr','Trnsprt')
rownames(zm6) <- colnames(zm6)

#Step one - set cost of lockin PUs to 0

polygon$ConsRes[polygon$lock_in_ConsRes ==T] <- F
polygon$ConsCon[polygon$lock_in_ConsRes ==T] <- F
polygon$CmmNnEx[polygon$lock_in_ConsRes ==T] <- F
polygon$CmmExtr[polygon$lock_in_ConsRes ==T] <- F
polygon$CmmFshr[polygon$lock_in_ConsRes ==T] <- F
polygon$Trnsprt[polygon$lock_in_ConsRes ==T] <- F

polygon$ConsRes[polygon$lock_in_ConsCon ==T] <- F
polygon$ConsCon[polygon$lock_in_ConsCon ==T] <- F
polygon$CmmNnEx[polygon$lock_in_ConsCon ==T] <- F
polygon$CmmExtr[polygon$lock_in_ConsCon ==T] <- F
polygon$CmmFshr[polygon$lock_in_ConsCon ==T] <- F
polygon$Trnsprt[polygon$lock_in_ConsCon ==T] <- F

polygon$ConsRes[polygon$lock_in_CmmNnEx ==T] <- F
polygon$ConsCon[polygon$lock_in_CmmNnEx ==T] <- F
polygon$CmmNnEx[polygon$lock_in_CmmNnEx ==T] <- F
polygon$CmmExtr[polygon$lock_in_CmmNnEx ==T] <- F
polygon$CmmFshr[polygon$lock_in_CmmNnEx ==T] <- F
polygon$Trnsprt[polygon$lock_in_CmmNnEx ==T] <- F

polygon$ConsRes[polygon$lock_in_Trnsprt ==T] <- F
polygon$ConsCon[polygon$lock_in_Trnsprt ==T] <- F
polygon$CmmNnEx[polygon$lock_in_Trnsprt ==T] <- F
polygon$CmmExtr[polygon$lock_in_Trnsprt ==T] <- F
polygon$CmmFshr[polygon$lock_in_Trnsprt ==T] <- F
polygon$Trnsprt[polygon$lock_in_Trnsprt ==T] <- F

Hierarchical approach

define a problem without boundary penalties to optimality (gap = 0)

p1 <-
problem(
polygon,
zones(names(pp1),names(pp2), names(pp3),names(pp4),names(pp5),names(pp6),
zone_names = c("ConsRes", "ConsCon",'CmmNnEx','CmmExtr','CmmFshr','Trnsprt')),
cost_column = c("ConsRes", "ConsCon",'CmmNnEx','CmmExtr','CmmFshr','Trnsprt') ) %>%
add_locked_in_constraints(c("lock_in_ConsRes",'lock_in_ConsCon','lock_in_CmmNnEx','lock_in_CmmExtr','lock_in_CmmFshr','lock_in_Trnsprt')) %>%
add_min_set_objective()%>%
add_relative_targets(targets) %>%
add_mandatory_allocation_constraints()%>%
add_gurobi_solver(numeric_focus = T,threads = parallel::detectCores(TRUE)-1,gap=0)
s1 <- solve(p1)

calculate cost

s1_cost <- eval_cost_summary(p1, s1[, c("solution_1_ConsRes","solution_1_ConsCon","solution_1_CmmNnEx", "solution_1_CmmExtr","solution_1_CmmFshr",'solution_1_Trnsprt')])$cost

calculate cost threshold values for each zone

threshold_consres <- s1_cost[1] + (s1_cost[1] * seq(1e-5, 4, length.out = 4))
threshold_consres <- ceiling(threshold_consres)

threshold_concon <- s1_cost[2] + (s1_cost[2] * seq(1e-5, 4, length.out = 4))
threshold_concon <- ceiling(threshold_concon)

threshold_cmmN<- s1_cost[3] + (s1_cost[3] * seq(1e-5, 4, length.out = 4))
threshold_cmmN <- ceiling(threshold_cmmN)

threshold_CmmE <- s1_cost[4] + (s1_cost[4] * seq(1e-5, 4, length.out = 4))
threshold_CmmE <- ceiling(threshold_CmmE)

threshold_Fish <- s1_cost[5] + (s1_cost[5] * seq(1e-5, 4, length.out = 4))
threshold_Fish <- ceiling(threshold_Fish)

threshold_Trpt <- s1_cost[6] + (s1_cost[6] * seq(1e-5, 4, length.out = 4))
threshold_Trpt <- ceiling(threshold_Trpt)

add a column with zeros

polygon$zerosCR <- 0
polygon$zerosCC <- 0
polygon$zerosCN <- 0
polygon$zerosCE <- 0
polygon$zerosFS <- 0
polygon$zerosTr <- 0

define a problem with zero cost values and boundary penalties

p2 <-
problem(
polygon,
zones(names(pp1),names(pp2), names(pp3),names(pp4),names(pp5),names(pp6),
zone_names = c("ConsRes", "ConsCon",'CmmNnEx','CmmExtr','CmmFshr','Trnsprt')),
cost_column = c("zerosCR", "zerosCC",'zerosCN','zerosCE','zerosFS','zerosTr') ) %>%
add_locked_in_constraints(c("lock_in_ConsRes",'lock_in_ConsCon','lock_in_CmmNnEx','lock_in_CmmExtr','lock_in_CmmFshr','lock_in_Trnsprt')) %>%
add_min_set_objective()%>%
add_boundary_penalties(penalty = 0.01, data = bm_ply, zones=zm6,edge_factor= rep(0.5,6)) %>%
add_relative_targets(targets) %>%
add_mandatory_allocation_constraints()%>%
add_gurobi_solver(numeric_focus = T,threads = parallel::detectCores(TRUE)-2)

generate prioritizations based on each cost threshold

list <- llist()
for(i in 1:4){
s <-
p2 %>%
add_linear_constraints(threshold_consres[i], sense = "<=", data = c("ConsRes",'zerosCC','zerosCN','zerosCE','zerosFS','zerosTr'))
s <-
s %>%
add_linear_constraints(threshold_consres[i], sense = "<=", data = c("ConsRes",'zerosCC','zerosCN','zerosCE','zerosFS','zerosTr'))
s <-
s %>%
add_linear_constraints(threshold_concon[i], sense = "<=", data = c("zerosCR",'ConsCon','zerosCN','zerosCE','zerosFS','zerosTr'))
s <-
s %>%
add_linear_constraints(threshold_cmmN[i], sense = "<=", data = c("ConsRes",'zerosCC','CmmNnEx','zerosCE','zerosFS','zerosTr'))
s <-
s %>%
add_linear_constraints(threshold_CmmE[i], sense = "<=", data = c("ConsRes",'zerosCC','zerosCN','CmmExtr','zerosFS','zerosTr'))
s <-
s %>%
add_linear_constraints(threshold_Fish[i], sense = "<=", data = c("ConsRes",'zerosCC','zerosCN','zerosCE','CmmFshr','zerosTr'))
solve(s)
s <- data.frame(ConConsRes = s1$solution_1_ConsRes, ConsCon=s1$solution_1_ConsCon, CmmNnEx=s1$solution_1_CmmNnEx,CmmExtr=s1$solution_1_CmmExtr,CmmFshr=s1$solution_1_CmmFshr,Trnsprt=s1$solution_1_Trnsprt)
names(s) <- c(paste0("ConRes_threshold_", threshold_consres[i]),paste0("ConCon_threshold_", threshold_concon[i]),paste0("CmmNnEx_threshold_", threshold_cmmN[i]),paste0("CmmExtr_threshold_", threshold_CmmE[i]),paste0("CmmFshr_threshold_", threshold_Fish[i]),paste0("ConRes_Trnsprt_", threshold_Trpt[i]))
list[[i]] <- s

}

[jeffreyhanson]

Hmm, if I understand correctly, you're setting linear constraints to specify a budget for each zone? I think there might be an issue with this block:

  s <-
    p2 %>%
    add_linear_constraints(threshold_consres[i], sense = "<=", data = c("ConsRes",'zerosCC','zerosCN','zerosCE','zerosFS','zerosTr')) 
  s <-
    s %>%
    add_linear_constraints(threshold_consres[i], sense = "<=", data = c("ConsRes",'zerosCC','zerosCN','zerosCE','zerosFS','zerosTr'))
  s <-
    s %>%
    add_linear_constraints(threshold_concon[i], sense = "<=", data = c("zerosCR",'ConsCon','zerosCN','zerosCE','zerosFS','zerosTr')) 
  s <-
    s %>%
    add_linear_constraints(threshold_cmmN[i], sense = "<=", data = c("ConsRes",'zerosCC','CmmNnEx','zerosCE','zerosFS','zerosTr')) 
  s <-
    s %>%
    add_linear_constraints(threshold_CmmE[i], sense = "<=", data = c("ConsRes",'zerosCC','zerosCN','CmmExtr','zerosFS','zerosTr')) 
  s <-
    s %>%
    add_linear_constraints(threshold_Fish[i], sense = "<=", data = c("ConsRes",'zerosCC','zerosCN','zerosCE','CmmFshr','zerosTr')) 

I would have expected these constraints to use the original cost values for each of the zones. E.g., something like this:

  s <-
    p2 %>%
    add_linear_constraints(
      threshold_consres[i], sense = "<=", 
      data = c("ConsRes",'zerosCC','zerosCN','zerosCE','zerosFS','zerosTr')
    )  %>%
    add_linear_constraints(
      threshold_concon[i], sense = "<=", 
      data = c("zerosCR",'ConsCon','zerosCN','zerosCE','zerosFS','zerosTr')
    )  %>%
    add_linear_constraints(
      threshold_cmmN[i], sense = "<=", 
      data = c("zerosCR",'zerosCC','CmmNnEx','zerosCE','zerosFS','zerosTr')
    )  %>%
    add_linear_constraints(
      threshold_CmmE[i], sense = "<=", 
      data = c("zerosCR",'zerosCC','zerosCN','CmmExtr','zerosFS','zerosTr')
    )  %>%
    add_linear_constraints(
      threshold_Fish[i], sense = "<=", 
      data = c("zerosCR",'zerosCC','zerosCN','zerosCE','CmmFshr','zerosTr')
    )  %>%
    add_linear_constraints(
      threshold_Trpt[i], sense = "<=", 
      data = c("zerosCR",'zerosCC','zerosCN','zerosCE','zerosFS','Trnsprt')
    ) %>%
   solve()

How does that sound? Sorry if I'm misunderstanding something.

[Tegan151]

Hi Jeff,

You are 100% correct except that the 'Trnsprt' zone is just a dummy zone where I lock in certain features like anchorage areas and mariculture.

I have one more question, the script above has been running forever and is just not finishing. Is there a way to add a linear constraint for the total cost of the solution rather than for each zone separately?

[jeffreyhanson]

Sorry for the slow response, I've been on holiday with limited internet access for the last couple of weeks.

Hmm, to clarify, did the script get to the step where it starts solving the problem (e.g., you can see messages from the solver), or does it hang before this?

Yeah, to add constraints for the overall cost, you could do something like this (where total_cost_budget is the maxiumum total cost across all zones):

  s <-
    p2 %>%
    add_linear_constraints(
      total_cost_budget, sense = "<=", 
      data = c("ConsRes",'ConsCon','CmmNnEx','CmmExtr','CmmFshr','Trnsprt')
    ) 

How does that sound?

[Tegan151]

Hi Jeff,

Thank you, I managed to figure out the piece of code in your response. It has been running since last Wednesday. It quickly solved the first problem but the second problem with a threshold of total_cost_budget + ( total_cost_budget * 0.5) is not solving.

The problems are starting to be solved but the gap never falls below 16% so it is never solved.

I found a paper that ran problems with a boundary penalty at different increment and then selected the appropriate boundary penalty as the inflection point between cost and sum of the perimeters of the solution zones. What do you think of this approach?

[jeffreyhanson]

Ah - I see - thanks for clarifying that!

Hmm, I can't think of any immediate reason for why it's only getting to 16% and then not making much progress. Could you please share (or copy and paste) the Gurobi solver output text shown on the R console? There might be some information that might help with diagnosing this behavior.

I wonder if you could try using the solution for the first problem as a starting solution in subsequent problems to help speed things up? E.g., replace this line:

add_gurobi_solver(numeric_focus = T,threads = parallel::detectCores(TRUE)-2)

with this line:

add_gurobi_solver(
  numeric_focus = T,
  threads = parallel::detectCores(TRUE)-2,
  start_solution = s1[, c(
    "solution_1_ConsRes",
    "solution_1_ConsCon",
    "solution_1_CmmNnEx", 
    "solution_1_CmmExtr",
    "solution_1_CmmFshr",
    "solution_1_Trnsprt"
  )]
)

Yeah, I would expect the current approach -- which involves using hierachical optimization -- to solve faster than that approach (which involves using blended optimization). You could certainly give this a go though.

[jeffreyhanson]

Just had an idea: rather than boundary penalties, you could try using connecitivity penalties with an adjacency matrix. This would be a bit simpler mathematically, and would give very similar results when you have planning units that conform to a regular grid.

E.g., you can create the adjacency matrix like this:

adj_ply <- adjacency_matrix(polygon)
zm6 <- diag(6)
colnames(zm6) <- c("ConsRes", "ConsCon",'CmmNnEx','CmmExtr','CmmFshr','Trnsprt')
rownames(zm6) <- colnames(zm6)

and then replace the lines where you use add_boundary_penalties() with:

add_connectivity_penalties(penalty = 1, data = adj_ply, zones=zm6) %>%

Broadly speaking, this approach would say that we want to promote connectivity by allocating pairs of planning units that are next to each other (i.e., adjacent, defined in adj_ply), weighted by the values in the zone matrix (i.e., zm6). Although this would that the algorithm would be biased towards "clumping"/"selecting" planning units that have (relatively) lots of neighbors (e.g., if some planning units have 20 neighbors and others only have 1 neighbor), this isn't an issue for your data since you have a regular grid (i.e., all planning units have roughly the same number of neighbors).

How does that sound?