Simulated Annealing in NetLogo

If you have played a little bit with the previous models in local search, you have realized that there are many cases where the hill climbing algorithms get stuck in not too good optimal states. For these cases there exist algorithms that sometimes can give relatively better solutions by allowing sometimes movements that goes against improving. We will see in this section the most common of them: Simulated Annealing.

This algorithm is inspired by a physical phenomenon observed in the metal annealing and crystallization processes. In this inspiration, an energy is associated with every state in our search space. This energy function reaches its minimum (usually \(0\), but not necessarily) in the solution we look for. Hence, the search again has been transformed into an optimization problem where a global minimum is searched. A searcher in this space will try to descend locally by comparing energies in order to reach a minimum... exactly in the same way as in hill climbing. The difference now is that we can provide an extra energy to the searcher (in the form of temperature) that allows it to move to worse states (with higher energy, so theoretically further from the minimum). With this little trick we can avoid the searcher to get stuck in a local minimum by jumping small hills to see if, behind them, there is a better optimal state.

If you would allow this behaviour indefinitely you would find the problem that the searcher will never stop in the optimal and will be jumping forever from one state to others, so our solution is to promote the exploration of the space allowing higher temperatures (extra energies) and, as the search time goes, reduce this temperature to avoid jumps and promote the exploitation of a good enough area of states (by hill-climbing similar ways).

As this method allows longer jumps than the usual hill-climbing, it can be applied to search spaces with more complex structures, bigger and roughness than the common one... but it still can find problems getting stuck in local optimum states.

You can find a nice example of implementation of this algorithm in the models library of NetLogo. Here we will present two more samples trying to make more clear how to use this algorithm for real problems.

Solving Sudoku (Download Model)

We will not introduce here what a Sudoku is... a simple query in any search engine will provide you with more information, examples, and interactive games than you can afford in your next years.

We will use Simulated Annealing to provide a model that generates valid Sudoku boards (a \(9\times 9\) board verifying all the restrictions). In this problem, one state will be any distribution of numbers on the grid that make use of 9 copies of every number between 1 to 9, not necessarily verifying the restrictions. Only to give an idea of the size of this state space we will say that the number of valid Sudoku grids for the standard 9×9 grid was calculated by Bertram Felgenhauer and Frazer Jarvis in 2005 to be 6,670,903,752,021,072,936,960. The result was derived through logic and brute force computation. Russell and Jarvis also showed that when symmetries were taken into account, there were 5,472,730,538 solutions. We leave to the reader the computation of the available states to browse in our search for one of the valid grid.

If you want to get more mathematical information about Sudoku, you can visit this wonderful page on Wikipedia.

This problem is very appropriated for local search methods, as we can transform easily states in similar states (by swapping 2 cells, for example), and measure the energy of the state by measuring the number of constraints present in the grid (in fact, we will compute in constraint variable some measure of the failed restrictions). Specifically, for every cell we compute three different values (that will be sumed up):

• We count the number of cells that has the same value in the same row.
• We count the number of cells that has the same value in the same column.
• We count the number of cells that has the same value in the same block.
• Finally, we sum up this three numbers, \(s\), and return \(e^s\). Later, we will use this constraints to decide probabilistically which cells will be swapped. By using the exponential we get to make exponentially more probable to take cells with higher error.

Taking into account that we will store the Sudoku grid on the patches (in a world of size \(9\times 9\)), and that we add a property to patches storing the block it belongs to, the code in NetLogo that implements this calculation is:

patches-own [
  value          ; Value in the cell
  block          ; Block 3x3 the cell belongs to
  constraints    ; Number of constraints the cell has in the current state:
                 ;    Number Cells in the same Row with the same value
                 ;    + Number Cells in the same Column with the same value
                 ;    + Number Cells in the same Block with the same value
                 ;   Finally, we take the exponential of this sum.
]
    
to compute-constraints
  ask patches [
    set constraints e ^ (same-in-row + same-in-column + same-in-block)
  ]
end
    
; Auxiliary reports to compute the different constraints
to-report same-in-row
  let x pxcor
  report -1 + length filter [z -> z = value] ([value] of patches with [pxcor = x])
end
    
to-report same-in-column
  let y pycor
  report -1 + length filter [z -> z = value] ([value] of patches with [pycor = y])
end
    
to-report same-in-block
  let b block
  report -1 + length filter [z -> z = value] ([value] of patches with [block = b])
end

We must consider that, when a cell has no constraints (it verifies all the restrictions) the energy it contributes to the system energy is \(1\) (\(e^0\)), consequently, the minimal energy we must reach to be sure that we have a valid Sudoku state is \(81\) (there are \(81\) cells).

Temperature and global-energy are global variables that store what their names show. The setup procedure simply creates an initial state (where several restrictions can fail) and sets the parameters to launch the search:

to setup
  clear-all
  ; Fill the initial data on patches
  ask patches [
    ; Block the cell belongs to (9 blocks 3x3)
    set block (3 * int (pycor / 3)) + 1 + int (pxcor / 3)
    ; Each block is filled with a slightly different color
    set pcolor scale-color black block -1 12
    ; Initial value in the cell: ordered in a row
    set value 1 + pxcor
    set plabel value
  ]
  ; Make a random swap of the cells
  repeat 100 [
    swap-values (one-of patches) (one-of patches)
  ]
  ; Compute the constraints of every cell
  compute-constraints
  ; The global energy of the system is the sum of local energies (constraints) of the cells
  set global-energy (sum [constraints] of patches )

  ; Update the view showing
  update-view
  ; Initial temperature
  set temperature 1
  ; Reset the tick counter
  reset-ticks
end
    
; Auxiliary swap procedure between values of cells
to swap-values [ p1 p2 ]
  let temp [ value ] of p1
  ask p1 [ set value [value] of p2 ]
  ask p2 [ set value temp ]
end
    
; Update View procedure. It gives a scale-red color to the label of the cell according to the
; relative energy of the cell
;    dark: correct, bright: higher energy
; It also plots the energy
to update-view
  let M max [constraints] of patches
  ask patches [
    set plabel-color scale-color red constraints 0 M
    set plabel value
  ]
  plot global-energy
end

The Go procedure, and auxiliary reports, that implement one step of the iterative process of the Simulated Annealing Algorithm are directly coded as follows (we have made use of the rnd extension, that allows to select one agent from an agentset proportionally to some internal property, in this case the constraints):

to go
  ; In every tick, with the same temperature, we make some swaps
  repeat 5
  [
    ; The swapping is made between 2 random cells. The cells are chosen proportionally to 
    ; their constraints, in order to change the worst with higher probability
    let p1 rnd:weighted-one-of patches [constraints]
    let p2 rnd:weighted-one-of patches [constraints]
    try-swap p1 p2
  ]
  ; Cool the system
  set temperature temperature * (1 - cooling-rate / 100)
  ; Compute the energy after the swaps
  set global-energy (sum [constraints] of patches )
  update-view
  ; The miminum energy of the system is 81
  if global-energy = 81 [stop]
  tick
end
    
; Swap procedure between cells p1 and p2
to try-swap [p1 p2]
  ; If they have the same value, there is no need to change them
  if ([value] of p1 = [value] of p2)
  [ stop ]
  ; Try the swap
  swap-values p1 p2
  ; Compute new energy after the swapping
  compute-constraints
  let new-energy (sum [constraints] of patches)
  ; If the swap is accepted
  ifelse (accept-swap? global-energy new-energy)
  [
    ; Set the new energy as global
    set global-energy new-energy
  ]
  [
    ; If not, swap them again and compute their constraints
    swap-values p1 p2
    compute-constraints
  ]
end
    
; A swap is accepted if it improves the energy, or randomly depending on the temperature of 
; the system (this is the central idea of the Algorithm)
to-report accept-swap? [ previous-energy new-energy]
  report (new-energy <= previous-energy) or ((random-float 1.0) < temperature)
end

With small changes in the code we can build from the previous model a general Sudoku Solver. You can find a simple proposal for this task in this link.

Solving Nonograms (Download Model)

(From Wikipedia) Nonograms, also known as Hanjie, Picross or Griddlers, are picture logic puzzles in which cells in a grid must be colored or left blank according to numbers at the side of the grid to reveal a hidden picture. In this puzzle type, the numbers are a form of discrete tomography that measures how many unbroken lines of filled-in squares there are in any given row or column. For example, a clue of [4 8 3] would mean there are sets of four, eight, and three filled squares, in that order, with at least one blank square between successive groups.

In this section we will give a model that will solve (in some cases) nonograms with a Simulated Annealing Algorithm. For that, we will use the patches as the cells of the grid where the nonogram must be solved, and some global variables to store information about it:

globals [
  rows-dist     ; Distribution of cells in rows, for example: [3] or [2 5 2]
  cols-dist     ; Distribution of cells in columns
  Global-Energy ; Global Energy of the system
  Temperature   ; Temperature of the system
]
    
patches-own [
  state           ; False/True = deactivated/active
  previous-state  ; Previous state when we try changes
]

One sample of nonogram can be the following (observe that we define the figure as a procedure that stores the rows and columns clues in the corresponding global variables):

to Figure1
  ;rows-dist: up to down
  set rows-dist [ [2] [1] [4] [5] [3 2] [2 1] [3] ]
  
  ;cols-dist: left to right
  set cols-dist [ [1] [3] [1 4] [4 1] [4] [3] [2] ]
end

The setup procedure prepares the model to run the algorithm. It makes use of auxiliary procedures show-states and update-plot that makes what we expect from them. Later we will analyze in detail the Compute-Global-Energy procedure that calculates the function to be minimized.

to setup
  ca
  ; Load Selected Figure (it simply fill rows-dist and cols-dist variables)
  ; Figure is a choose widget in the interface to select the figure to load
  run Figure
  ; Resize the world to fit the figure
  resize-world 0 (length cols-dist - 1) 0 (length rows-dist - 1)
  set-patch-size (25 * 16) / (max (list world-width world-height))
  ; Starts cell randomly (in the interface we choose the % of initially activated cells)
  ask patches [
    set state (random 100 < %-initial)
    set previous-state true
  ]
  ; Show the current states of the cells
  show-states
  ; Compute the current Energy of the system
  Compute-Global-Energy
  ; Update plot
  Update-Plot
  ; Reset the Temperature
  set Temperature 1
  reset-ticks
end

To compute the global energy of the state we sum the energies of all the rows and columns. The error of a row/column is obtained by comparing the clue of this line in the current state with the clue of the same line in the goal state. For that, the first step is to convert the line of states in a clue by grouping them (group and auxiliary-group reports are in charge of this process).

If we get \(v_1\) as the grouping of a line, and \(v_2\) is the goal grouping, the error of the line is computed as the euclidean distance between them. If they are not of the same length, the shortest is extended to equalize them by adding \(0\)'s, and a penalty is also considered proportional to the difference of lengths.

to Compute-Global-Energy
  let rows-error sum map [x -> Energy-Row x] (range world-height)
  let cols-error sum map [x -> Energy-Column x] (range world-width)
  set Global-Energy ( rows-error + cols-error )
end
    
; The energy of a row is computed from the error between the pattern of the goal row and 
; the current row
to-report Energy-Row [row]
  if row < 0 or row > max-pycor [report 0]
  let Row-states map [p -> [state] of p] (sort patches with [pycor = row])
  let Row-Pattern group Row-states
  let Row-Error Compute-Error Row-Pattern (item row rows-dist)
  report Row-Error
end
    
; The energy of a column is computed from the error between the pattern of the goal column 
; and the current column
to-report Energy-Column [col]
  if col < 0 or col > max-pxcor [report 0]
  let Col-States map [p -> [state] of p] (sort patches with [pxcor = col])
  let Col-Pattern group Col-States
  set Col-Pattern reverse Col-Pattern
  let Col-Error Compute-Error Col-Pattern (item col cols-dist)
  report Col-Error
end
    
; The error between 2 patterns
to-report Compute-Error [ v1 v2 ]
  ; Compute the differnce in lengths
  let dif abs (length v1 - length v2)
  ; Euqalize the patterns by adding 0`s ti the shortest
  if (length v1) < (length v2)
  [ set v1 sentence v1 (n-values dif [0]) ]
  if (length v2) < (length v1)
  [ set v2 sentence v2 (n-values dif [0]) ]
  ; Compute the euclidean distance between patterns
  let er sum (map [(?1 - ?2) ^ 2] v1 v2)
  ; Adding a penalty for the diference of lengths
  set er er + ( dif * Weight-Dif)
  report er
end
    
; Report to get the groups of a row/column of states:
;   [true true true false false true true false] -> [3 2]
; It works with reduce, leaving the false's and counting consecutive true's with the help of 
;   an auxiliary report. After that, we must remover the false's
to-report group [states]
  let res filter is-number? (reduce auxiliary-group (fput [false] states))
  ifelse res = []
  [ report [0] ]
  [ report res ]
end
    
to-report auxiliary-group [L x]
  ; false's are added directly
  ifelse x = false
  [report lput x L]
  [
    ; if x is a true, we must see if the last one was a true too or not. We recognize this 
    ; because L ends in a number that we must increment
    let a last L
    ifelse is-number? a
    [ report lput (a + 1) (bl L) ]
    [
      ; Otherwise, x is the first true of a serie... and we count it as 1
      report lput 1 L
    ]
  ]
end

As in the previous samples, the go procedure implements one step of the iterative process of the algorithm. As you can see in the code, we implement two different strategies to change a cell:

• In the first one we simply change the state of a cell (on / off).
• In the second one, according to the weights set in the interface, we apply one of the following operators:
  • Set off the cell.
  • Set on the cell.
  • Swap the value with one of its neighbours.

You can try different strategies and add them to the code easily.

Observe that, in every step, we randomly choose the cell to be changed. We leave to the reader the option to give a procedure where the cell is chosen in a different way (for example, higher probability for cells with higher errors), and check if the efficiency of the algorithm is improved.

to go
  ; Make a change on one cell
  ask one-of patches [
    ; We compute the part of the energy depending on this cell
    let old-energy Compute-Context-Energy pxcor pycor
    ; and store the previous states
    Store-Previous-States
    ; Apply the change following one of the strategies. It will affect only to the neighbors 
    ; of the cell
    if Strategy = "Change-me" [strategy1]
    if Strategy = "Probabilistic" [strategy2]
    ; Compute the new energy after the change
    let current-energy Compute-Context-Energy pxcor pycor
    ; If the change is not accepted, we reverse the changes made
    if not Accept-Change? current-energy old-energy [ Reverse-changes ]
  ]
  ; Update the view of the cells
  show-states
  ; Compute and plot the energy of the system
  Compute-Global-Energy
  Update-Plot
  ; Stop if the process has reached the goal
  if Global-Energy = 0 [stop]
  ; Cooling process
  set Temperature Temperature * .999
  tick
end
    
; Procedures to store (and reverse) states before any change
to Store-Previous-States
  ask neighbors4 [ set previous-state state ]
  set previous-state state
end
    
to Reverse-changes
  ask neighbors4 [ set state previous-state ]
  set state previous-state
end
    
; Report to accept or not the changes. It is the usual method
to-report Accept-Change? [ new old ]
  ifelse new < old
  [ report true ]
  [
    let prob exp ( ( old - new ) / Temperature )
    report random-float 1.0 < prob
  ]
end
    
; Report to evaluate the energy in the context (row and column) of the cell at (x,y)
to-report Compute-Context-Energy [ x y ]
  report sum map Energy-Row (Context y) + sum map Energy-Column (Context x)
end
    
; Auxiliary report to take the context of a location
to-report Context [x]
  report (list (x - 1) x (x + 1))
end