feature request: specify a function for `constrained_method` parameter?

Hi,

Thank you very much for developing this package!

I'm working on solving a non-linear mixed integer programming problem, and I think this package would be very helpful. In particular, my problem has a single linear constraint and many binary variables (it's like a knapsack problem with a non-linear objective function). I'm finding that the solver generally has a difficult time generating feasible solutions, and so often ends up returning one of the initial solutions that I manually specify as the best solution. To address this, I was wondering if there was some way I could manually specify a function to tell the optimization process how to generate a feasible solution from an infeasible solution. This way I can use domain knowledge to help the optimization process generate better candidate solutions (i.e., at the step where it generates new candidate solutions in a new iteration). How does that sound?

That said, I'm only very new to using the package, so maybe I'm using the wrong settings or the wrong algorithm for my problem? For reference, I've been using the genetic algorithm (algorithm_id = "GA") solver. What do you think?

Hi,

I think there is no need for that as you can provide your own initial population via initial_population in minimize method. Please check ?MinimizerOpts.

Cheers, Davide.

Thanks so much for the quick response @dr4kan!

Yeah, I've been specifying my own initial population to ensure that the starting solutions are feasible. This definitely helps! If I don't do this, then minimize() function crashes (with an error about C stack usage) -- which I'm guessing occurs because the optimization process has trouble generating feasible solutions for the initial population. However, when I manually specify an initial population -- which contains only feasible solutions -- it seems appears that the best solution found during the optimization process is often one of the starting solutions in the initial population. This suggests -- to me, but maybe I'm wrong? -- that the process for generating new solutions in a given iteration is having trouble generating feasible solutions (because an infeasible solution would receive an additional penalty that reduces the quality of the optimization process, meaning that it is unlikely to result in an improvement).

To help explain this, I've included a reprex below. Please let me know if you have any issues running it, or questions about it?

# Initialization
## set seed for reproducibility
set.seed(500)

## load required packages
library(surveyvoi)
library(assertthat)
library(EmiR)
library(scales)

## define parameters for simulating data for optimization problem
### number of planning units
n_pu <- 300
### number of species
n_f <- 10
### sparsity of species in planning units
sparsity <- 0.7
### target level of representation for species by solution
target <- 5
### specify budget as a proportion of total costs
budget_prop <- 0.05

## define functions

#' Generate a prioritization
#'
#' Generate a prioritization for protected area establishment.
#'
#' @param rij `matrix` containing the probability of each species occurring
#' within each planning unit.
#'
#' @param pu_costs `numeric` vector containing the planning unit costs.
#' within each planning unit.
#'
#' @param pu_locked_in `logical` vector indicating which planning units
#' are locked in to the solution.
#'
#' @param pu_locked_out `logical` vector indicating which planning units
#' are locked out from the solution.
#'
#' @param target `numeric` vector indicating the representation target for
#' each species.
#'
#' @param budget `numeric` value indicating the total budget for the
#' prioritization.
#'
#' @param control `list` containing parameters for configuring the
#' optimization algorithm. These are passed to `EmiR::config_algo()`.
#' Defaults to `default_algo_config()`.
#'
#' @details
#' Briefly, this function aims to identify a set of planning units (i.e.
#' candidate  locations for protection) to protect that have a high probability
#' of safeguarding each of the species. To ensure that solutions are
#' economically feasible, the planning units selected for protection should not #' exceed a total budget. The \pkg{EmiR} package is used to perform
#' the optimization process.
#'
#' @return A `list` containing (`x`) the solution and (`objval`) the objective
#' value.
#'
#' @export
generate_prioritization <- function(rij, pu_costs, pu_locked_in, pu_locked_out,
                                    target, budget,
                                    control = default_algo_config()) {
  # assert arguments are valid
  assertthat::assert_that(
    ## rij
    is.matrix(rij),
    ## pu_costs
    is.numeric(pu_costs),
    assertthat::noNA(pu_costs),
    identical(ncol(rij), length(pu_costs)),
    ## pu_locked_in
    is.logical(pu_locked_in),
    assertthat::noNA(pu_locked_in),
    identical(ncol(rij), length(pu_locked_in)),
    ## pu_locked_out
    is.logical(pu_locked_out),
    assertthat::noNA(pu_locked_out),
    identical(ncol(rij), length(pu_locked_out)),
    ## target
    is.numeric(target),
    assertthat::noNA(target),
    identical(length(target), nrow(rij)),
    ## budget
    assertthat::is.number(budget),
    assertthat::noNA(budget),
    isTRUE(budget > 0)
  )

  # validate problem feasibility
  assertthat::assert_that(
    sum(pu_costs[pu_locked_in > 0.5]) <= budget,
    msg = "cost of locked in planning units exceeds budget"
  )

  # define parameters for problem
  n_pu <- length(pu_costs)
  n_spp <- nrow(rij)
  param_idx <- which(!pu_locked_in & !pu_locked_in)
  n_params <- length(param_idx)
  params <- EmiR::parameters(
    matrix(c(0, 1, 1), nrow = 3, ncol = n_params)
  )

  # define objective function
  obj <- function(x) {
    ## generate solution (accounting for locked in planning units
    v <- logical(n_pu)
    v[pu_locked_in] <- TRUE
    v[param_idx[x > 0.5]] <- TRUE
    ## calculate results for each species and sum togeather, and
    ## then multiply by minus because EmiR only does minimization
    surveyvoi:::rcpp_expected_value_of_action(
      solution = v,
      pij = rij,
      target = target
    )
  }

  # define constraint
  constr <- EmiR::constraint(
    func = function(x) {
      sum(pu_costs[param_idx[x > 0.5]]) + sum(pu_costs[pu_locked_in]) - budget
    },
    inequality = "<="
  )

  # define config for optimization
  config <- do.call(EmiR::config_algo, control)

  # generate initial solutions
  init <- generate_initial_solutions(
    pu_costs, pu_locked_in, pu_locked_out, n = control$population_size
  )

  # validate initial solutions
  ## calculate cost of initial solutiosn
  init_cost <-
    rowSums(init * matrix(
        pu_costs[param_idx],
        nrow = control$population_size,
        ncol = length(param_idx),
        byrow = TRUE
      )
    ) + sum(pu_costs[pu_locked_in > 0.5])
  ## validation
  assertthat::assert_that(
    all(init_cost <= budget),
    msg = "cost of initial solutions exceeds budget"
  )
  assertthat::assert_that(
    all(
      vapply(
        seq_len(control$population_size), FUN.VALUE = numeric(1), function(i) {
          constr@func(init[i, ])
        }
      ) <= 0
    ),
    msg = "failed to generate valid constraint"
  )

  # run optimization
  res <- EmiR::minimize(
    algorithm_id = control$algorithm_id,
    obj_func = obj,
    config = config,
    parameters = params,
    maximize = TRUE,
    constraints = list(constr),
    constrained_method = "PENALTY",
    initial_population = init,
    seed = 500
  )

  # return result
  res
}

#' Default algorithm configuration
#'
#' This function provides default configuration settings for
#' running `EmiR::minimize()`.
#'
#' @return A `list` containing configuration settings.
#'
#' @export
default_algo_config <- function() {
  list(algorithm_id = "GA", population_size = 200, iterations = 1000)
}

#' Generate initial solutions
#'
#' @inheritParams generate_prioritization
#'
#' @param n `numeric` value indicating the number of desired solutions.
#'
#' @details
#' This is useful for generating
#' an initial population of solutions for `EmiR::minimize()`.
#'
#' @return A `matrix` containing the initial solutions.
#'
#' @export
generate_initial_solutions <- function(pu_costs, pu_locked_in,
                                       pu_locked_out, n = 1000
) {
  # initialization
  n_pu <- length(pu_costs)
  param_idx <- which(!pu_locked_in & !pu_locked_in)
  sols <- matrix(NA_real_, nrow = n, ncol = n_pu)
  sols[, which(pu_locked_in)] <- 1
  sols[, which(pu_locked_out)] <- 0

  # generate starting solutions
  for (i in seq_len(n)) {
    ## find remaining planning units
    cand_pu <- which(is.na(sols[i, ]))
    ## calculate cost of current solution
    curr_cost <- sum(sols[i, ] * pu_costs, na.rm = TRUE)
    ## remove candidate planning units that cost too much
    cand_pu <- cand_pu[which((pu_costs[cand_pu] + curr_cost) <= budget)]
    ## loop while candidates remain
    while (length(cand_pu) > 0) {
      ## randomly pick a planning unit
      curr_idx <- sample(seq_along(cand_pu), size = 1, replace = FALSE)
      curr_pu <- cand_pu[curr_idx]
      ## add it to the solution
      sols[i, curr_pu] <- 1
      ## remove selected unit
      cand_pu <- cand_pu[-curr_idx]
      ## update cost
      curr_cost <- sum(sols[i, ] * pu_costs, na.rm = TRUE)
      ## remove candidate planning units that cost too much
      cand_pu <- cand_pu[which(pu_costs[cand_pu] + curr_cost <= budget)]
    }
  }

  ## set remaining planning units with NAs to zeros
  sols[is.na(sols)] <- 0

  # return starting values for parameters
  sols[, param_idx, drop = FALSE]
}

# Preliminary processing
### these values denote the probability of each species occurring
### in each planning unit. This is generally zero-inflated and
### and skewed towards lower probability values
rij <- matrix(
  scales::rescale(exp(runif(n_pu * n_f)), to = c(0, 1)),
  nrow = n_f,
  ncol = n_pu
)
### note that we loop over each species individually to ensure that
### each species has at least one non-zero probability value
for (i in seq_len(n_f)) {
  rij[i, runif(ncol(rij)) < sparsity] <- 0
}

### these values denote the cost of each planning unit
pu_costs <- exp(runif(n_pu, 0.1, 3))

### this value denotes the total budget for the prioritization.
### the total cost of the selected planning units cannot exceed this
### threshold
budget <- sum(pu_costs) * budget_prop

### these values control if certain planning units should always
### be selected or not selected. Although this useful for real-world
### analyses, we will just set them to FALSE since we're just looking at
### simulated data
pu_locked_in <- rep(FALSE, n_pu)
pu_locked_out <- rep(FALSE, n_pu)

# Main processing
## generate prioritization by running the optimization process
res <- generate_prioritization(
  rij = rij, pu_costs = pu_costs,
  pu_locked_in = pu_locked_in, pu_locked_out = pu_locked_out,
  target = rep(target, n_f), budget = budget,
  control = list(algorithm_id = "GA", population_size = 200, iterations = 10000)
)
#> 0%   10   20   30   40   50   60   70   80   90   100%
#> [----|----|----|----|----|----|----|----|----|----|
#> **************************************************|
#> 163.815 sec elapsed

## view summary of best solution at each iteration
print(summary(res@cost_history))
#>   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
#>  3.394   3.394   3.394   3.394   3.394   3.394

## we can see that the cost history has the same minimum and maximum value,
## indicating that the best solution found over the course of the entire
## optimization is simply one of the initial solutions

## if we try running this with a higher `budget_prop` value
## (e.g., `budget_prop = 0.7), then the res@cost_history shows an improvement
## in the best solution over subsequent iterations. Since this issue
## only occurs when `budget_prop` is low (i.e., when the criteria for
## feasibility is more strict; e.g., a random generated solution is less likely
## to be feasible), this suggests that the issue is because the optimization
## process is unable to generate new feasible solutions when a low `budget_prop`
## parameter is used.

Also, here's my session information.

R version 4.3.3 (2024-02-29)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.4 LTS

Matrix products: default
BLAS:   /opt/R/R-4.3.3/lib/R/lib/libRblas.so 
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3;  LAPACK version 3.10.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

time zone: Pacific/Auckland
tzcode source: system (glibc)

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] scales_1.3.0          PoissonBinomial_1.2.6 EmiR_1.0.4           
 [4] assertthat_0.2.1      surveyvoi_1.0.6       nloptr_2.0.3         
 [7] sf_1.0-16             Matrix_1.6-5          testthat_3.2.1       
[10] devtools_2.4.5        usethis_2.2.2        

loaded via a namespace (and not attached):
 [1] class_7.3-22       KernSmooth_2.23-22 stringi_1.8.3      lattice_0.22-5    
 [5] digest_0.6.35      magrittr_2.0.3     grid_4.3.3         pkgload_1.3.4     
 [9] fastmap_1.1.1      pkgbuild_1.4.3     sessioninfo_1.2.2  e1071_1.7-14      
[13] brio_1.1.4         DBI_1.2.2          urlchecker_1.0.1   promises_1.2.1    
[17] purrr_1.0.2        Rdpack_2.6         cli_3.6.2          shiny_1.8.0       
[21] rlang_1.1.3        rbibutils_2.2.16   units_0.8-5        munsell_0.5.0     
[25] ellipsis_0.3.2     remotes_2.4.2.1    cachem_1.0.8       tools_4.3.3       
[29] memoise_2.0.1      colorspace_2.1-0   mathjaxr_1.6-0     httpuv_1.6.13     
[33] vctrs_0.6.5        R6_2.5.1           mime_0.12          proxy_0.4-27      
[37] lifecycle_1.0.4    classInt_0.4-10    stringr_1.5.1      tictoc_1.2        
[41] fs_1.6.3           htmlwidgets_1.6.4  miniUI_0.1.1.1     later_1.3.2       
[45] glue_1.7.0         profvis_0.3.8      Rcpp_1.0.12        xtable_1.8-4      
[49] htmltools_0.5.7    compiler_4.3.3

Hi,

I see the problem now. I added the possibility of specifying a custom generator function for new solution. I will do some tests and upload to CRAN in few days.

Cheers, Davide

Brilliant - thanks so much @dr4kan!! I really appreciate your efforts. If it would be helpful, I'm happy to help with testing before you submit to CRAN? Would you mind providing (or linking to) an an example on how to use it? Sorry, I couldn't find an example in the recent changes.

Hi,

with the latest commit (release 1.1.0), everything is working and we will upload to cran tomorrow. You can take a look at how to use the new generation_functionhere:

https://github.com/dr4kan/EmiR/blob/master/demo/Test_custom_gen_func.R

Cheers, Davide.

Hi @dr4kan,

Thanks so much for implementing this functionality so quickly! Before you upload to CRAN, I was just wondering if you could please help me understand the inteded purpose of the generation_function parameter? After playing around with the example code, it seems the generation_function is used to generate an initial solutions for the optimization process. In other words, it provides an alternate means to achieve the same functionality as the initial_population parameter? Is that right? Here's a reprex showing how I arrived at this this conclusion - sorry if I'm misunderstanding something?

library(EmiR)

schaffer2 <- function(x) {
  x1 <- x[1]
  x2 <- x[2]
  fact1 <- (sin(x1^2-x2^2))^2 - 0.5
  fact2 <- (1 + 0.001*(x1^2+x2^2))^2
  y <- 0.5 + fact1/fact2
  return(y)
}

g1 <- function(x) {
  x1 <- x[1]
  x2 <- x[2]
  return(x1^2 + x2^2 - 1)
}

custom_gen <- function(x) {
  return(c(5, 6))
}

c1 <- constraint(g1, "<=")

p1 <- parameter("x1", -1000, 1000)
p2 <- parameter("x2", -1000, 1000)

conf <- config_algo(algorithm_id = "PS", population_size = 4, iterations = 3)
results1 <- minimize(
  algorithm_id = "PS",
  obj_func = schaffer2,
  config = conf,
  parameters = list(p1,p2),
  constraints = list(c1),
  save_pop_history = TRUE,
  constrained_method = "PENALTY",
  constr_init_pop = TRUE,
  oob_solutions = "RBC",
  penalty_scale = 5,
  generation_function = custom_gen,
  seed = 1
  )

results2 <- minimize(
  algorithm_id = "PS",
  obj_func = schaffer2,
  config = conf,
  parameters = list(p1,p2),
  constraints = list(c1),
  save_pop_history = TRUE,
  constrained_method = "PENALTY",
  constr_init_pop = TRUE,
  oob_solutions = "RBC",
  penalty_scale = 5,
  initial_population = t(matrix(c(5, 6), ncol = 4, nrow = 2)),
  seed = 1
)

# print results for first optimization run using the generation_function parameter
print("results1@pop_history")
print(results1@pop_history)

# print results for secondoptimization run using the initial_population parameter
print("results2@pop_history")
#> [[1]]
#> [[1]][[1]]
#> [1] 5 6
#> 
#> [[1]][[2]]
#> [1] 5 6
#> 
#> [[1]][[3]]
#> [1] 5 6
#> 
#> [[1]][[4]]
#> [1] 5 6
#> 
#> 
#> [[2]]
#> [[2]][[1]]
#> [1] 121.9557 440.6631
#> 
#> [[2]][[2]]
#> [1]  848.7635 -102.3339
#> 
#> [[2]][[3]]
#> [1] -99.87292 478.41838
#> 
#> [[2]][[4]]
#> [1] 679.86798  41.57036
#> 
#> 
#> [[3]]
#> [[3]][[1]]
#> [1]  -46.94668 -226.00352
#> 
#> [[3]][[2]]
#> [1] 182.0969 145.5347
#> 
#> [[3]][[3]]
#> [1]  196.7938 -188.2483
#> 
#> [[3]][[4]]
#> [1]   13.20132 -100.68211

print(results2@pop_history)
#> [[1]]
#> [[1]][[1]]
#> [1] 5 6
#> 
#> [[1]][[2]]
#> [1] 5 6
#> 
#> [[1]][[3]]
#> [1] 5 6
#> 
#> [[1]][[4]]
#> [1] 5 6
#> 
#> 
#> [[2]]
#> [[2]][[1]]
#> [1] -0.6656158  2.5738958
#> 
#> [[2]][[2]]
#> [1] -4.710028  1.150294
#> 
#> [[2]][[3]]
#> [1] 0.5573530 0.5407451
#> 
#> [[2]][[4]]
#> [1] -3.7734869  0.7684152
#> 
#> 
#> [[3]]
#> [[3]][[1]]
#> [1] 6.436982 3.471733
#> 
#> [[3]][[2]]
#> [1] 18.310339  2.731418
#> 
#> [[3]][[3]]
#> [1] 3.889338 4.635186
#> 
#> [[3]][[4]]
#> [1] 7.615351 3.506206

As we can see, the both results1 andresults2yield similar population histories where the specified solution (i.e.,c(5,6)`) is only being used in the first/initial iteration of the populaiton history, which suggests to me they are being implemented as alternative methods to achieve the same desired functionality?

Assuming I'm correctly understanding the intended purpose of the generation_function parameter, I'm not sure if this is exactly the functionality I need to fix my issue?

I'm sorry about the confusion, I'll try explain the feature request in more detail. A typical metaheuristical algorithm involves (at least) the following steps: (1) generate an initial population of solutions, (2) evaluate solutions in the population, (3) if any of the solutions have better performance than ever previously recorded, store the solution as the best solution, (4) generate new solutions from the population of solutions, (6) repeat steps 2 -- 4 until termination criteria are met (e.g., certain number of interations are completed), and (5) return best solution. Does that sound right to you?

My problem is not that I cannot generate a feasible initial population of solutions (i.e., not in step 1). My problem is that when generating new solutions in a subsequent iteration (i.e., step 4), I cannot seem to generate feasible solutions. One potential method to address this would be to add an additional step that can transform infeasible solutions into feasible solutions. So, the algorithm might look like this: (1) generate an initial population of solutions, (2) evaluate solutions in the population, (3) if any of the solutions have better performance than ever previously recorded, store the solution as the best solution, (4) generate new solutions from the population of solutions, (5) if any new solutions are infeasible, then modify them so that they are feasible, (6) repeat steps 2 -- 5 until termination criteria are met, and (6) return best solution.

As an example, perhaps the API for this could look something like this?

library(EmiR)

# specify objective function
obj <- function(x) {
  x1 <- x[1]
  x2 <- x[2]
  x1^2 + x2^2 - abs(x1 - 0.5) - abs(x2 - 0.5)
}

# specify a simple linear constraint,
# the sum of x1 and x2 must be less less than 100
g1 <- function(x) {
  x1 <- x[1]
  x2 <- x[2]
  return(x1 + x2 - 100)
}
c1 <- constraint(g1, "<=")

# define parameters
p1 <- parameter("x1", 0, 1000)
p2 <- parameter("x2", 0, 1000)

# specify a custom transformation function to turn
# infeasible solutions into feasible solutions
## here x is a matrix containing infeasible solutions
##  in the same format as the initial_population parameter
## (i.e. each row is different solution, column 1 = x1, column 2 = x2),
## this function will return a new matrix of exactly the same 
## dimensions/length as x (just it will have different values)
custom_transform <- function(x) {
   # if x1+x2 > 100, then we will return new parameters
   # that maintain the relative proportions of x1 and x2, but 
   # ensure they sum to 100
   x <- x / colSums(x)
   return(x * 100)
}

# configure optimization process
conf <- config_algo(algorithm_id = "PS", population_size = 4, iterations = 3)

# run optimization process
results <- minimize(
  algorithm_id = "PS",
  obj_func = schaffer2,
  config = conf,
  parameters = list(p1,p2),
  constraints = list(c1),
  save_pop_history = TRUE,
  constrained_method = "PENALTY",
  constr_init_pop = TRUE,
  oob_solutions = "RBC",
  penalty_scale = 5,
  custom_transform = custom_transform,
  seed = 1
)

Hi, the new feature should do exactly what you are asking for, but I cannot exclude there are still few things to fix. Try this.

Use the release 1.1.1 and as oob_solutions use DIS.

In principle, when a you have a solution which violates a constraint, a new one is generated and if you provided a custom generation_function that function will be used to generate a new solution.

Awesome - thanks so much for the quick response and confirming that! I really appreciate all your efforts!

I had a go at trying to solve the original problem I posted earlier. You're absolutely right - by using the generation_function, the optimization process is able to generate new solutions and avoid getting suck with just the initial population of solutions. It seems that I had to use constrained_method = "ACCREJ" in order to get minimize() to call the generation function - is that correct?

Also, I was wondering if it would be possible for the generation function to have the infeasible solution as a parameter/argument. This way, rather than generating a totally new/random solution, it might be possible to generate a solution that is a modified version of the original infeasible solution? I'll see if I can put togeather a PR.

Thanks again!

I've just opened a PR (#3) - please let me know if you have any questions, or suggestions for improving it?

dr4kan / EmiR

feature request: specify a function for `constrained_method` parameter? #2