# Simulated Annealing in NetLogo

Última modificación: 17 de Octubre de 2019, y ha tenido 618 vistas

Etiquetas utilizadas: || || || || ||

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.

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
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
; 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
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 SolverYou can find a simple proposal for this task in this link.

(From WikipediaNonograms, also known as HanjiePicross 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)
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 0s 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]
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
; 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`