How to create binary constraints for optimization in R?-CodePudding

I have a function f(x) which I intend to minimize. "x" is a vector containing 50 parameters. This function has several constraints: first is that all parameters in x should be binary, so that x = (1,1,0,1,...); second is that the sum of "x" should be exactly 25, so that sum(x) = 25. The question can be illustrated as:

min f(x)

s.t. sum(x) = 25, x = 0 or 1

However when I try to solve this problem in R, I met some problems. Prevalent packages such as "optim","constrOptim" from "stats" can only input coefficients of the target function (in my case, the function is bit complex and cannot be simply illustrated using coefficient matrix), "donlp2" from "Rdonlp" does not support setting parameters to be binary. I'm wondering whether anyone has any idea of how to set binary constraints for this case?

CodePudding user response：

You can try the amazing package rgenoud. Below is an example.

I take 20 binary variables instead of your 50 for easier reading. I take f(x) = sum(1:20 * x), this is a weighted sum with increasing weights so clearly the best solution (restricted to sum(x)=10) is 1, 1, ..., 1, 0, 0, ..., 0. And rgenoud brilliantly finds it.

library(rgenoud)

f <- function(x) { # the function to be minimized
  sum(1:20 * x)
}

g <- function(x){
  c(
    ifelse(sum(x) == 10, 0, 1), # set the constraint (here sum(x)=10) in this way
    f(x) # the objective function (to minimize/maximize)
  )
}

solution <- genoud(
  g, 
  pop.size = 3000,
  lexical = 2, # see ?genoud for explanations
  nvars = 20,  # number of x_i's
  starting.values = c(rep(0, 10), rep(1, 10)), 
  Domains = cbind(rep(0, 20), rep(1, 20)), # lower and upper bounds
  data.type.int = TRUE # x_i's are integer
)

solution$par # the values of x
## [1] 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0
solution$value
## [1] 0 55 ; 0 is the value of ifelse(sum(x)=10,0,1) and 55 is the value of f(x)

CodePudding user response：

Expanding my comment, here is an example of a Local Search, as implemented in package NMOF. (I borrow Stéphane's objective function).

library("NMOF")
library("neighbours")

## Stéphane's objective function
f <- function(x)
    sum(1:20 * x)

nb <- neighbourfun(type = "logical", kmin = 10, kmax = 10)

x0 <- c(rep(FALSE, 10), rep(TRUE, 10))
sol <- LSopt(f, list(x0 = x0, neighbour = nb, nI = 1000))

## initial solution
as.numeric(x0)
## [1] 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1

## final solution
as.numeric(sol$xbest)
## [1] 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0

(Disclosure: I am the maintainer of packages NMOF and neighbours.)