Sheet metal texture

Simulated Annealing

Sometimes, a prob¬≠lem does¬≠n‚Äôt have a known, ef¬≠Ô¨Ā¬≠cient, de¬≠ter¬≠min¬≠is¬≠tic al¬≠go¬≠rithm. This hic¬≠cup is¬≠n‚Äôt nec¬≠es¬≠sar¬≠ily a big deal be¬≠cause of¬≠ten we only care about get¬≠ting an an¬≠swer that is close enough. One al¬≠go¬≠rithm that can give ap¬≠prox¬≠i¬≠mate val¬≠ues for hard to solve prob¬≠lems is sim¬≠u¬≠lated an¬≠neal¬≠ing. Before jump¬≠ing into the R code, let‚Äôs look at a mo¬≠ti¬≠vat¬≠ing ex¬≠am¬≠ple.

Motivation

The trav¬≠el¬≠ing sales¬≠man prob¬≠lem does¬≠n‚Äôt have a poly¬≠no¬≠mial-time de¬≠ter¬≠min¬≠is¬≠tic al¬≠go¬≠rithm and thus is a good can¬≠di¬≠date for us to ex¬≠plore. The prob¬≠lem name it¬≠self nearly gives away the de¬≠tails. Imagine a map with some Ô¨Āxed num¬≠ber of ran¬≠domly placed cities. We now want to hire a sales¬≠per¬≠son‚ÄĒlet‚Äôs name her Sally‚ÄĒwho needs to visit every city, with¬≠out re¬≠vis¬≠it¬≠ing any city, to sell some prod¬≠uct.

When the num¬≠ber of cities is small, we can use brute force tech¬≠niques to try every pos¬≠si¬≠ble per¬≠mu¬≠ta¬≠tion of city or¬≠der¬≠ings. Of course, for a rel¬≠a¬≠tively large amount of cities, try¬≠ing to Ô¨Ānd the route this way will be in¬≠tractable. We aren‚Äôt com¬≠pletely doomed though. Do we care if Sally takes the most op¬≠ti¬≠mal route? For this sit¬≠u¬≠a¬≠tion, not re¬≠ally. While we want Sally to take a fast route, it does¬≠n‚Äôt have to be the fastest route. Luckily, it turns out that sim¬≠u¬≠lated an¬≠neal¬≠ing can tractably gen¬≠er¬≠ate fast routes.

How does it work?

The trav¬≠el¬≠ing sales¬≠man prob¬≠lem is a lit¬≠tle com¬≠pli¬≠cated, so let‚Äôs set up a toy prob¬≠lem in¬≠stead. For this demon¬≠stra¬≠tion, we are go¬≠ing to Ô¨Ānd the global min¬≠i¬≠mum of the sine func¬≠tion for the do¬≠main [0, 2ŌÄ]. For the math¬≠e¬≠mat¬≠i¬≠cally in¬≠clined, you could use tech¬≠niques from cal¬≠cu¬≠lus class to Ô¨Ānd that the an¬≠swer is -1. Thus, if our al¬≠go¬≠rith¬≠m‚Äôs an¬≠swer is not in the same ball¬≠park as -1, then we know some¬≠thing went wrong.

State Spaces

Simulated an¬≠neal¬≠ing searches through a state space it¬≠er¬≠a¬≠tively. An it¬≠er¬≠a¬≠tive al¬≠go¬≠rithm be¬≠gins with an ini¬≠tial guess. This guess is then po¬≠ten¬≠tially up¬≠graded every it¬≠er¬≠a¬≠tion of the al¬≠go¬≠rithm to pro¬≠duce a bet¬≠ter ap¬≠prox¬≠i¬≠ma¬≠tion. For our demon¬≠stra¬≠tion, we‚Äôll end the al¬≠go¬≠rithm af¬≠ter a Ô¨Āxed num¬≠ber of it¬≠er¬≠a¬≠tions.

Taking care of the ini¬≠tial guess is easy‚Äďjust choose a ran¬≠dom valid state. From the ini¬≠tial guess, we‚Äôll need to Ô¨Ānd it‚Äôs neigh¬≠bor¬≠ing states and check to see if they‚Äôre bet¬≠ter. Disappointingly, there is no sin¬≠gle method to de¬≠rive neigh¬≠bor¬≠ing states. We‚Äôll have to come up with a method spe¬≠ciÔ¨Āc to our prob¬≠lem.

Here, I want to keep things sim¬≠ple. The im¬≠ple¬≠men¬≠ta¬≠tion be¬≠low adds a ran¬≠dom num¬≠ber from the range [-0.1, 0.1] to the cur¬≠rent state to gen¬≠er¬≠ate the suc¬≠ces¬≠sor. I also wrap the new state value around if it ex¬≠ceeds the range [0, 2ŌÄ].

gen_successor <- function(x) {
  TWO_PI <- 2 * pi
  successor <- x + runif(1, -.1, .1)
  
  if(successor < 0) {
    return(TWO_PI + successor)
  }
  else if(successor > TWO_PI) {
    return(0 + (successor - TWO_PI))
  }
    
  successor
}

Of course, you can ex­per­i­ment on your own. For ex­am­ple, try think­ing of new ways to gen­er­ate neigh­bor­ing states.

Temperature Schedule Go Brrr

This tem¬≠per¬≠a¬≠ture idea is what gives sim¬≠u¬≠lated an¬≠neal¬≠ing its name. The tem¬≠per¬≠a¬≠ture value of the al¬≠go¬≠rithm is anal¬≠o¬≠gous to the role of tem¬≠per¬≠a¬≠ture when an¬≠neal¬≠ing metal. To an¬≠neal metal, the metal is Ô¨Ārst heated and then slowly cooled over time. This cool¬≠ing process al¬≠lows the metal to have im¬≠proved duc¬≠til¬≠ity and re¬≠duced brit¬≠tle¬≠ness. By en¬≠dur¬≠ing a long cool¬≠ing pe¬≠riod, the metal is im¬≠proved. Similarly, the rate at which the ‚Äúheat‚ÄĚ of the al¬≠go¬≠rithm is de¬≠creased de¬≠ter¬≠mines the ac¬≠cu¬≠racy of the ap¬≠prox¬≠i¬≠ma¬≠tion.

The ‚Äúheat‚ÄĚ of the al¬≠go¬≠rithm af¬≠fects the prob¬≠a¬≠bil¬≠ity that it will make a bad move on pur¬≠pose. Surprisingly, the al¬≠go¬≠rithm some¬≠times needs to make bad moves to reach the global op¬≠ti¬≠mum. Why? There is a chance that it might get stuck on a lo¬≠cal op¬≠ti¬≠mum. By mak¬≠ing bad moves, the al¬≠go¬≠rithm can po¬≠ten¬≠tially stum¬≠ble its way out of lo¬≠cal op¬≠ti¬≠mum (but this is¬≠n‚Äôt guar¬≠an¬≠teed).

You can choose any num­ber to be the ini­tial tem­per­a­ture. The tem­per­a­ture sched­ule then de­ter­mines how the tem­per­a­ture value will de­crease every it­er­a­tion. Like choos­ing neigh­bor­ing states, there is no right way to de­crease the tem­per­a­ture. The only rule of thumb is that the slower the tem­per­a­ture is dropped, the bet­ter the ap­prox­i­ma­tion will be. The trade­off is that you’ll have to wait longer to get an an­swer. Thus, choose some­thing that short enough for your at­ten­tion span and long enough to get a good an­swer.

For our toy ex­am­ple, I choose to sub­tract an ar­bi­trary num­ber every it­er­a­tion.

schedule <- function(t) {
  t - .00009
}

Once again, feel free to get cre­ative and try dif­fer­ent ini­tial val­ues and sched­ules for your im­ple­men­ta­tion.

Probabilistically Making Moves

Every it¬≠er¬≠a¬≠tion a neigh¬≠bor¬≠ing state is gen¬≠er¬≠ated. The al¬≠go¬≠rithm will al¬≠ways re¬≠place the cur¬≠rent state with the new state if it‚Äôs bet¬≠ter. As we just dis¬≠cussed, the new state may be worse‚Äďdrag¬≠ging us fur¬≠ther away from an op¬≠ti¬≠mum. Yet, we want to ran¬≠domly make some bad moves with some prob¬≠a¬≠bil¬≠ity. The key to sim¬≠u¬≠lated an¬≠neal¬≠ing is that this prob¬≠a¬≠bil¬≠ity is not con¬≠stant. Rather, it on av¬≠er¬≠age de¬≠creases over time.

The prob¬≠a¬≠bil¬≠ity of choos¬≠ing the new state is de¬≠Ô¨Āned as the fol¬≠low¬≠ing for¬≠mula:

Two pa¬≠ra¬≠me¬≠ters af¬≠fect the prob¬≠a¬≠bil¬≠ity The ‚Äúchange in E‚ÄĚ pa¬≠ra¬≠me¬≠ter is the dif¬≠fer¬≠ence be¬≠tween the func¬≠tion ap¬≠plied with the po¬≠ten¬≠tial new state and the func¬≠tion ap¬≠plied with the cur¬≠rent state. The big¬≠ger this dif¬≠fer¬≠ence, the smaller the prob¬≠a¬≠bil¬≠ity. The next pa¬≠ra¬≠me¬≠ter is the tem¬≠per¬≠a¬≠ture value, which we‚Äôve al¬≠ready seen. Obviously, for every it¬≠er¬≠a¬≠tion of the al¬≠go¬≠rithm, the tem¬≠per¬≠a¬≠ture value will de¬≠crease‚Äďcaus¬≠ing the prob¬≠a¬≠bil¬≠ity to de¬≠crease slowly.

As the tem­per­a­ture gets closer to zero, the prob­a­bil­ity of choos­ing a bad move will also get closer to zero. Towards the end, the al­go­rithm will mostly take greedy moves to­wards bet­ter states. Once the tem­per­a­ture hits zero, then the al­go­rithm is done. We can get our an­swer and go home.

Put Everything Together

We’ve cov­ered all the im­por­tant parts of this al­go­rithm. Let’s put every­thing to­gether and see it in ac­tion. Here is the full code for my so­lu­tion.

library(ggplot2)

simulated_annealing <- function(gen_successor, get_value, schedule) {
  generated_states <- c()
  temp <- 10 # initial temperature
  current_state <- .5 # initial state
  current_state_value <- get_value(current_state)
  
  while(temp > 0) {
    generated_states <- append(generated_states, current_state)
    temp <- schedule(temp)
    next_state <- gen_successor(current_state)
    next_state_value <- get_value(next_state)
    change_in_e = next_state_value - current_state_value
    
    if(change_in_e < 0) {
      current_state <- next_state
      current_state_value <- next_state_value
    }
    else {
      probability <- exp(-change_in_e / temp)
      samp <- runif(1)
      
      if(samp <= probability) {
        current_state <- next_state
        current_state_value <- next_state_value
      }
    }
  }
  
  return(generated_states)
}

gen_successor <- function(x) {
  TWO_PI <- 2 * pi
  successor <- x + runif(1, -.1, .1)
  
  if(successor < 0) {
    return(TWO_PI + successor)
  }
  else if(successor > TWO_PI) {
    return(0 + (successor - TWO_PI))
  }
    
  successor
}

get_value <- function(x) {
  sin(x)
}

schedule <- function(t) {
  t - .00009
}

generated_states <- simulated_annealing(gen_successor, get_value, schedule)
index <- seq(length(generated_states))
sin_values <- sin(generated_states)

data <- data.frame(
  index = index, 
  x = generated_states, 
  sin = sin_values
  )

ggplot(data, aes(x = index, y = sin)) +
  geom_smooth() +
  ylab("Sine value") +
  xlab("Iteration")

I de­cided to cap­ture the state for every it­er­a­tion and graph the re­sult. As you can see, the ini­tial state is nowhere near -1. As the al­go­rithm pro­gresses, it be­gins to slowly con­verge. Towards the end, the prob­a­bil­ity of choos­ing bad moves is near 0, which is why the down­ward slope for later it­er­a­tions is steep.

Simulated annealing plot

For this run, I got a Ô¨Ā¬≠nal state of 4.703329. The cor¬≠re¬≠spond¬≠ing sine value is -0.9999590, which is close to -1. These two val¬≠ues match what we ex¬≠pected.