Tasked with optimising a complicated system, I went down the route of heuristic optimisation and ended up creating a simulated annealing algorithm in R.

I don't want to go into too much detail about how the model works, because that has been done well enough elsewhere, but essentially it is a smart way of searching through all the possible solutions to a problem to find the best one, while avoiding getting stuck in local minima. Taking inspiration from nature, it's based on the physical cooling of metals from a high temperature; when this is done very fast (e.g. by dropping the metal into water), the metal forms small grains (crystals), and still has a lot of energy trapped in the crystal structure. If this is done very slowly however, large grains form (more of molecules are aligned), and so the metal reaches a much lower (more optimal) energy state.

While it isn't guaranteed to find the absolute optimum solution, it will get you to a good solution fast, which is often what you really want.

undefined

Click permalink below to see the code!

Anyway, enough explanation. In this example, I'm going to optimise a wedding (or party, or ball, or bar mitzvah) seating plan. For the representation of the problem, I took some inspiration from this paper, although I am optimising with a different method.

So, I will add to this later, but for now, here is the main code:

source('./Projects/Seating Plan Problem/SeatingPlanFunctions.R')

# Example run:
N <- c('Bill','Sandra','Frank','Jane','Henry','Sarah','Jeremy','John','Rich','Geoff','Smithy','Neil','Arthur','Sophie','Emmet','Andy')

# Community matrix (which guests know which guests)
C <- matrix(c( 1,50,10,10,10,10,0,10,10,10,10,10,10,0,10,10,
               50,1,10,10,10,10,0,10,0,10,10,10,10,0,10,0,
               10,10,1,50,10,10,0,0,0,0,0,10,0,0,0,0,
               10,10,50,1,10,10,0,0,0,0,0,0,0,0,0,0,
               10,10,10,10,1,10,0,0,0,0,0,0,0,0,0,0,
               10,10,10,10,10,1,0,0,0,0,0,0,0,0,0,0,
               0,0,0,0,0,0,1,10,0,0,0,0,0,0,0,0,
               10,10,0,0,0,0,10,1,0,0,0,0,0,0,0,0,
               10,0,0,0,0,0,0,0,1,10,10,0,0,0,0,0,
               10,10,0,0,0,0,0,0,10,1,10,0,0,0,0,0,
               10,10,0,0,0,0,0,0,10,10,1,0,0,0,0,0,
               10,10,10,0,0,0,0,0,0,0,0,1,10,0,0,0,
               10,10,0,0,0,0,0,0,0,0,0,10,1,0,0,0,
               0,0,0,0,0,0,0,0,0,0,0,0,0,1,50,10,
               10,10,0,0,0,0,0,0,0,0,0,0,0,50,1,10,
               10,0,0,0,0,0,0,0,0,0,0,0,0,10,10,1), 16, 16, byrow=TRUE)

n <- nrow(C)

# Gender vector (not currently used)
S <- c(1,0,1,0,1,0,1,1,1,1,1,1,1,0,1,1)

# Table sizes
t <- c(3, 3, 4, 5, 4)

# Set seed for repeatability
set.seed(321)

initialSolution <- initialiseSolution(n, t)
initialSolution <- randomiseSolution(initialSolution)

finalSolution <- simulatedAnnealing(initialSolution, evaluateEnergy, getNeighbour, 2000, 0.0001, 100000)

colnames(finalSolution) <- N
finalSolution

...and here are the functions being called...

initialiseSolution <- function(n, t) {

  # check capacity is high enough!
  if(n>sum(t)) { stop('There is not enough table capacity for the number of guests') }
  
  # get number of tables
  nt <- length(t)
  
  # initialise a matrix with a row for each table, and a column for each guest
  initialSolution <- matrix(0, nt, n)
  
  j1 <- 0
  j2 <- 0
  for(i in 1:nt) {
    
    j1 <- j2 + 1
    j2 <- min(n,j1 + t[i] - 1)
    
    initialSolution[i,j1:j2] <- 1
    
    if (j2==n) break
    
    i <- i + 1
  }
  
  return(initialSolution)
  
}

randomiseSolution <- function(T) {
  return(T[,sample(c(1:ncol(T)))])
}

evaluateEnergy <- function(T) {
  CO <- t(T) %*% T  #co-seated matrix
  Energy <- -sum(CO * C)   #community Energy
  
  # Count of men per table
  #M<-T %*% G
  
  return(Energy)
}

getNeighbour <- function(T) {
  columnsToSwap <- ceiling(runif(2) * n)
  T[,columnsToSwap] <- T[,rev(columnsToSwap)]
  return(T)
}

simulatedAnnealing <- function(initialSolution, evaluateEnergyFunction, getNeighbourFunction, initialTemperature, coolingRate, maxIterations) {
  
  temperature <- initialTemperature
  currentSolution <- initialSolution
  bestSolution <- currentSolution
  bestEnergy <- evaluateEnergyFunction(bestSolution)
  iterations <- 0
  
  print(paste0("Starting annealing process. Initial energy: ", bestEnergy))
  
  while (temperature > 1) {
    if(iterations>=maxIterations) {
      print(paste0("Max iterations met. Exiting early at Temperature: ",temperature))
      break
    }
    
    currentEnergy<-evaluateEnergyFunction(currentSolution)
    neighbourSolution <- getNeighbourFunction(currentSolution)
    neighbourEnergy <- evaluateEnergyFunction(neighbourSolution)
    delta <- currentEnergy - neighbourEnergy
    
    if (delta>0) {
      currentSolution <- neighbourSolution
      currentEnergy <- neighbourEnergy
    } else if (exp(delta/temperature)>runif(1)) {
      currentSolution <- neighbourSolution
      currentEnergy <- neighbourEnergy
    }
    if (neighbourEnergy<bestEnergy) {
      bestSolution <- neighbourSolution
      bestEnergy <- neighbourEnergy
    }
    temperature <- temperature * (1-coolingRate)
    iterations <- iterations+1
  }
  
  print(paste0("Cooling finished. Number of iterations: ", iterations))
  print(paste0("Best energy solution found: ", bestEnergy))
  print(bestSolution)
  
  return(bestSolution)
  
}