Search code examples
roptimizationgenetic-algorithmknapsack-problemparticle-swarm

Combinatorial optimization with discrete options in R


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,]

Solution

  • 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