Combinatorial optimization with discrete options in R

Question

I have a function with five variables that I want to maximize using only an specific set of parameters for each variable.

Are there any methods in R that can do this, other than by brutal force? (e.g. Particle Swarm Optimization, Genetic Algorithm, Greedy, etc.). I have read a few packages but they seem to create their own set of parameters from within a given range. I am only interested in optimizing the set of options provided.

Here is a simplified version of the problem:

#Example of 5 variable function to optimize
Fn<-function(x){
  a=x[1]
  b=x[2]
  c=x[3]
  d=x[4]
  e=x[5]
  SUM=a+b+c+d+e
 return(SUM)
}

#Parameters for variables to optimize
Vars=list(
  As=c(seq(1.5,3, by = 0.3)),             #float
  Bs=c(1,2),                              #Binary
  Cs=c(seq(1,60, by=10)),                  #Integer
  Ds=c(seq(60,-60, length.out=5)),        #Negtive
  Es=c(1,2,3)
)

#Full combination
FullCombn= expand.grid(Vars)

Results=data.frame(I=as.numeric(),   Sum=as.numeric())
for (i in 1:nrow(FullCombn)){
  ParsI=FullCombn[i,]
  ResultI=Fn(ParsI)
  Results=rbind(Results,c(I=i,Sum=ResultI))
}

#Best iteration (Largest result)
Best=Results[Results[, 2] == max(Results[, 2]),]
#Best parameters
FullCombn[Best$I,]

Answer 1

Here is a genetic algorithm solution with package GA.
The key is to write a function decode enforcing the constraints, see the package vignette.

library(GA)
#> Loading required package: foreach
#> Loading required package: iterators
#> Package 'GA' version 3.2.2
#> Type 'citation("GA")' for citing this R package in publications.
#> 
#> Attaching package: 'GA'
#> The following object is masked from 'package:utils':
#> 
#>     de

decode <- function(x) {
  As <- Vars$As
  Bs <- Vars$Bs
  Cs <- Vars$Cs
  Ds <- rev(Vars$Ds)
  # fix real variable As
  i <- findInterval(x[1], As)
  if(x[1L] - As[i] < As[i + 1L] - x[1L])
    x[1L] <- As[i]
  else x[1L] <- As[i + 1L]
  # fix binary variable Bs
  if(x[2L] - Bs[1L] < Bs[2L] - x[2L])
    x[2L] <- Bs[1L]
  else x[2L] <- Bs[2L]
  # fix integer variable Cs
  i <- findInterval(x[3L], Cs)
  if(x[3L] - Cs[i] < Cs[i + 1L] - x[3L])
    x[3L] <- Cs[i]
  else x[3L] <- Cs[i + 1L]
  # fix integer variable Ds
  i <- findInterval(x[4L], Ds)
  if(x[4L] - Ds[i] < Ds[i + 1L] - x[4L])
    x[4L] <- Ds[i]
  else x[4L] <- Ds[i + 1L]
  # fix the other, integer variable
  x[5L] <- round(x[5L])
  setNames(x  , c("As", "Bs", "Cs", "Ds", "Es"))
}
Fn <- function(x){
  x <- decode(x)
  # a <- x[1]
  # b <- x[2]
  # c <- x[3]
  # d <- x[4]
  # e <- x[5]
  # SUM <- a + b + c + d + e
  SUM <- sum(x, na.rm = TRUE)
  return(SUM)
}

#Parameters for variables to optimize
Vars <- list(
  As = seq(1.5, 3, by = 0.3),              # Float
  Bs = c(1, 2),                            # Binary
  Cs = seq(1, 60, by = 10),                # Integer
  Ds = seq(60, -60, length.out = 5),       # Negative
  Es = c(1, 2, 3)
)

res <- ga(type = "real-valued", 
          fitness = Fn, 
          lower = c(1.5, 1, 1, -60, 1),
          upper = c(3, 2, 51, 60, 3),
          popSize = 1000,
          seed = 123)
summary(res)
#> ── Genetic Algorithm ─────────────────── 
#> 
#> GA settings: 
#> Type                  =  real-valued 
#> Population size       =  1000 
#> Number of generations =  100 
#> Elitism               =  50 
#> Crossover probability =  0.8 
#> Mutation probability  =  0.1 
#> Search domain = 
#>        x1 x2 x3  x4 x5
#> lower 1.5  1  1 -60  1
#> upper 3.0  2 51  60  3
#> 
#> GA results: 
#> Iterations             = 100 
#> Fitness function value = 119 
#> Solutions = 
#>             x1       x2       x3       x4       x5
#> [1,]  2.854089 1.556080 46.11389 49.31045 2.532682
#> [2,]  2.869408 1.638266 46.12966 48.71106 2.559620
#> [3,]  2.865254 1.665405 46.21684 49.04667 2.528606
#> [4,]  2.866494 1.630416 46.12736 48.78017 2.530454
#> [5,]  2.860940 1.650015 46.31773 48.92642 2.521276
#> [6,]  2.851644 1.660358 46.09504 48.81425 2.525504
#> [7,]  2.855078 1.611837 46.13855 48.62022 2.575492
#> [8,]  2.857066 1.588893 46.15918 48.60505 2.588992
#> [9,]  2.862644 1.637806 46.20663 48.92781 2.579260
#> [10,] 2.861573 1.630762 46.23494 48.90927 2.555612
#>  ...                                              
#> [59,] 2.853788 1.640810 46.35649 48.87381 2.536682
#> [60,] 2.859090 1.658127 46.15508 48.85404 2.590679
apply(res@solution, 1, decode) |> t() |> unique()
#>      As Bs Cs Ds Es
#> [1,]  3  2 51 60  3

^{Created on 2022-10-24 with reprex v2.0.2}

Answer 2

Two more possibilities. Both minimize by default, so I flip the sign in your objective function (i.e. return -SUM).

#Example of 5 variable function to optimize
Fn<-function(x, ...){
  a=x[1]
  b=x[2]
  c=x[3]
  d=x[4]
  e=x[5]
  SUM=a+b+c+d+e
 return(-SUM)
}

#Parameters for variables to optimize
Vars=list(
  As=c(seq(1.5,3, by = 0.3)),             #float
  Bs=c(1,2),                              #Binary
  Cs=c(seq(1,60, by=10)),                 #Integer
  Ds=c(seq(60,-60, length.out=5)),        #Negtive
  Es=c(1,2,3)
)

First, a grid search. Exactly what you did, just convenient. And the implementation allows you to distribute the evaluations of the objective function.

library("NMOF")
gridSearch(fun = Fn,
           levels = Vars)[c("minfun", "minlevels")]


## 5 variables with 6, 2, 6, 5, ... levels: 1080 function evaluations required.
## $minfun
## [1] -119
## 
## $minlevels
## [1]  3  2 51 60  3

An alternative: a simple Local Search. You start with a valid initial guess, and then move randomly through possible feasible solutions. The key ingredient is the neighbourhood function. It picks one element randomly and then, again randomly, sets this element to one allowed value.

nb <- function(x, levels, ...) {
    i <- sample(length(levels), 1)
    x[i] <- sample(levels[[i]], 1)
    x
}

(There would be better algorithms for neighbourhood functions; but this one is simple and so demonstrates the idea well.)

LSopt(Fn, list(x0 = c(1.8, 2, 11, 30, 2),  ## a feasible initial solution
               neighbour = nb,             
               nI = 200                    ## iterations
               ),
      levels = Vars)$xbest
## Local Search.
## ##...
## Best solution overall: -119
## [1]  3  2 51 60  3

(Disclosure: I am the maintainer of package NMOF, which provides functions gridSearch and LSopt.)

---

In response to the comment, a few remarks on Local Search and the neighbourhood function above (nb). Local Search, as implemented in LSopt, will start with an arbitrary solution, and then change that solution slightly. This new solution, called a neighbour, will be compared (by its objective-function value) to the old solution. If the new solution is better, it becomes the current solution; otherwise it is rejected and the old solution remains the current one. Then the algorithm repeats, for a number of iterations. So, in short, Local Search is not random sampling, but a guided random-walk through the search space. It's guided because only better solutions get accepted, worse one's get rejected. In this sense, LSopt will narrow down on good parameter values.

The implementation of the neighbourhood is not ideal for two reasons. The first is that a solution may not be changed at all, since I sample from feasible values. But for a small set of possible values as here, it might often happen that the same element is selected again. However, for larger search spaces, this inefficiency is typically negligible, since the probability of sampling the same value becomes smaller. Often so small, that the additional code for testing if the solution has changed becomes more expensive that the occasionally-wasted iteration.

A second thing could be improved, albeit through a more complicated function. And again, for this small problem it does not matter. In the current neighbourhood, an element is picked and then set to any feasible value. But that means that changes from one solution to the next might be large. Instead of picking any feasible values of the As, in realistic problems it will often be better to pick a value close to the current value. For example, when you are at 2.1, either move to 1.8 or 2.4, but not to 3.0. (This reasoning is only relevant, of course, if the variable in question is on a numeric or at least ordinal scale.)

Ultimately, what implementation works well can be tested only empirically. Many more details are in this tutorial.

Here is one alternative implementation. A solution is now a vector of positions for the original values, e.g. if x[1] is 2, it "points" to 1.8, if x[2] is 2, it points to 1, and so on.

## precompute lengths of vectors in Vars
lens <- lengths(Vars)

nb2 <- function(x, lens, ...) {
    i <- sample(length(lens), 1)
    if (x[i] == 1L) {
        x[i] <- 2
    } else if (x[i] == lens[i]) {
        x[i] <- lens[i] - 1
    } else
        x[i] <- x[i] + sample(c(1, -1), 1)
    x
}


## the objective function now needs to map the
## indices in x back to the levels in Vars
Fn2 <- function(x, levels, ...){

    y <- mapply(`[`, levels, x)
    ## => same as
    ##      y <- numeric(length(x))
    ##      y[1] <- Vars[[1]][x[1]]
    ##      y[2] <- Vars[[2]][x[2]]
    ##      ....

    SUM <- sum(y)
    return(-SUM)
}

xbest <- LSopt(Fn2,
               list(x0 = c(1, 1, 1, 1, 1),    ## an initial solution
                    neighbour = nb2,
                    nI = 200                  ## iterations
                    ),
               levels = Vars,
               lens = lens)$xbest
## Local Search.
## ....
## Best solution overall: -119

## map the solution back to the values
mapply(`[`, Vars, xbest)
## As Bs Cs Ds Es 
##  3  2 51 60  3