Theory
The Simulated Annealing method was introduced by Kirkpatrick et al. [1] in 1983. The annealing process of metals inspires this algorithm during the manufacturing process. The Simulated Annealing model is based on generating random neighbors from a starting point, similar to what occurs in the Hill Climbing Algorithm (see HC algorithm theory).
In the Simulated Annealing algorithm, the acceptance of the new solution (see new solution procedure in see HC algorithm theory) is given by a criterion that compares the system's energy given by equation (1). In this algorithm, the values of \(E_{i}\) are relative to the value of the objective function to \(i\) particle (current solution), i.e., \(E_{i}=of_{i}\).
| \[\Delta E = E_{new} - E_{cur}\] | (1) |
The value \(E_{new}\) is the value of the objective function for the newly generated neighbor, and \(E_{cur}\) is the value of the objective function for the current particle.
The solution will be accepted if \( E_{cur} > E_{new} \left( P(\Delta E,T)=1.00 \right) \). For solutions of type \(E_{cur} < E_{new}\) the acceptance will follow a certain probability given by equation (2) where \(T\) the annealing temperature.
| \[P(\Delta E,T) = e^{\frac{-\Delta E}{T}}\] | (2) |
At the end of the algorithm, the temperature is updated. \(\alpha\) is the cooling temperature adjustment factor. The methods for updating the temperature are (see Figure 1):
- Geometric [3]: See equation (3);
- Lundy [5]: See equation (4);
- Linear [4]: See equation (5);
- Exponential [2]: See equation (6);
| \[T_{i+1} = \alpha \cdot T_{i} \] | (3) |
| \[T_{i+1} = \frac{T_{i}}{1+ \alpha \cdot T_{i}}\] | (4) |
| \[T_{i+1} = T_{i} - \alpha \cdot T_{i} \] | (5) |
| \[T_{i+1} = T_{i} \cdot e^{-\alpha \cdot i}\] | (6) |
Figure 1. Colling schema.
Temperature is a crucial control parameter. Temperature is reduced systematically in SA procedures. Rapid cooling produces irregularities in the crystal structure that do not reach a minimum energy level. In contrast, a very slow cooling scheme provides ideal crystals with minimum attainable energy but may require prohibitive calculation time. Movements of molecules, particularly at high temperatures, are chaotic, and hence, there is a possibility of not reaching the minimum energy for these temperatures. In optimization SA, the mechanism mimics this by allowing acceptation of a worse solution. The probability of such acceptance depends on temperature and should approach zero for low temperatures. A chance of accepting temporarily a worse solution is an important feature of SA strategy. This mechanism should allow SA algorithm escaping from a local optimum [9].
Start temperature
The initial temperature \(T_0\) should be high enough large enough so that nearly all transitions are accepted at the first iterations. The initial probability of acceptance must not be close to one, neither must be close to zero. The probability of accepting a higher-cost solution was set to 0.80 [6,7,8]. The initial temperature is given by equation (7):
| \[T_0 = \frac{-\Delta E^+}{\ln{0.80}}\] | (7) |
\(\Delta E^+\) represents the energy of strictly positive transitions, i.e. \(\Delta E > 0 \).
Algorithm
1: Input initial parameters (temp, cov, n_population, x_lower, x_upper, obj_function, n_dimensions)
2: Input initial guess (x_pop)
3: Calculate of and fit (initial population)
4: for iter in range(n_iterations):
5: x_temp = neighbor solution equation hill climbing
6: delta_e = equation (1)
7: if delta_e > 0:
8: prob = 1
9: else:
10: prob = equation (2)
11: r = random number [0,1]
12: if r <= prob:
13: x_pop(iter+1) = x_temp
14: else:
15: x_pop(iter+1) = x_pop(iter)
16: temp_update = equation (3) - (6)
The Simulated Annealing algorithm can accept unfavorable moves. This characteristic depends on temperature and change in objective value ΔE.
See SA algorithm in METApy Framework.
Example 1
Use the Simulated Annealing optimization method to optimize the 2D sphere function. Use a total of 2 iterations to perform the optimization. Consider the limits \(\mathbf{x}_L = [-5.0, -5.0]\) and \(\mathbf{x}_U = [5.0, 5.0]\) for the problem design variables. Consider the initial guess (two agents) \(\mathbf{pop}_0 = [-0.74, 1.25]\) and \(\mathbf{pop}_1 = [3.58, -3.33]\). Use \(cov = 20%\), Gaussian random generator, \(T_0 = 15\) and geometric schedule (\(\alpha = 0.90\)).
Solution
Simulated Annealing 01 - report
Initial population
x0 = [-0.74, 1.25], of_pop 2.1101 - best solution
x1 = [3.58, -3.33], of_pop 23.9053
Iterations
Iteration: 1
Temperature: 15
Pop id: 0 - particle movement - mutation procedure
current x = [-0.74, 1.25]
Dimension 0: mean = -0.74, sigma = 0.14800000000000002, neighbor = -0.6602162206983541
Dimension 1: mean = 1.25, sigma = 0.25, neighbor = 1.2137927198636664
update x = [-0.6602162206983541, 1.2137927198636664], of = 1.9091782248672549, fit = 0.34373968272281763
energy = -0.20092177513274523, prob. state = 1
random number=0.154246375420666 <= prob. state=1 - accept this solution
Pop id: 1 - particle movement - mutation procedure
current x = [3.58, -3.33]
Dimension 0: mean = 3.58, sigma = 0.716, neighbor = 3.401325017838678
Dimension 1: mean = -3.33, sigma = 0.6659999999999999, neighbor = -3.4053363365858322
update x = [3.401325017838678, -3.4053363365858322], of = 23.165327442247097, fit = 0.04138160355533802
energy = -0.7399725577529033, prob. state = 1
random number=0.5668852469307849 <= prob. state=1 - accept this solution
update solutions
x0 = [-0.6602162206983541, 1.2137927198636664], of_pop 1.9091782248672549 - best solution
x1 = [3.401325017838678, -3.4053363365858322], of_pop 23.165327442247097
Iteration: 2
Temperature: 13.5
Pop id: 0 - particle movement - mutation procedure
current x = [-0.6602162206983541, 1.2137927198636664]
Dimension 0: mean = -0.6602162206983541, sigma = 0.1320432441396708, neighbor = -0.7783591115824063
Dimension 1: mean = 1.2137927198636664, sigma = 0.24275854397273328, neighbor = 1.4456134129098128
update x = [-0.7783591115824063, 1.4456134129098128], of = 2.69564104616811, fit = 0.27058905004772243
energy = 0.7864628213008551, prob. state = 0.9434079272824024
random number=0.21250149387436257 <= prob. state=0.9434079272824024 - accept this solution
Pop id: 1 - particle movement - mutation procedure
current x = [3.401325017838678, -3.4053363365858322]
Dimension 0: mean = 3.401325017838678, sigma = 0.6802650035677357, neighbor = 5.12575825256521
Dimension 1: mean = -3.4053363365858322, sigma = 0.6810672673171665, neighbor = -4.378437656968545
update x = [5.0, -4.378437656968545], of = 44.170716315960206, fit = 0.02213823648500941
energy = 21.00538887371311, prob. state = 0.21098784972914733
random number=0.5070777353026081 > prob. state=0.21098784972914733 - not accept this solution
update solutions
x0 = [-0.7783591115824063, 1.4456134129098128], of_pop 2.69564104616811 - best solution
x1 = [3.401325017838678, -3.4053363365858322], of_pop 23.165327442247097
Simulated Annealing - Guia de Cálculo
População Inicial
A população inicial é composta pelos seguintes vetores de solução:
- \( x_0 = [-0.74, 1.25], \ \text{of}_{\text{pop}} = 2.1101 \) - melhor solução
- \( x_1 = [3.58, -3.33], \ \text{of}_{\text{pop}} = 23.9053 \)
Temperatura Inicial Automática
O cálculo da temperatura inicial segue um procedimento automatizado, calculado conforme visto na equação (2). A média da temperatura calculada foi de \(T_0 = 36.061\) , com a soma dos valores de \(T_0\) igual a \(sum(T_0) = 24124.992855206798\). Por fim, o número de movimentos aceitos \(\Delta E > 0\) foi de \(669\).
Cálculos de Vizinhança
A cada iteração, o algoritmo gera novas soluções vizinhas a partir das soluções atuais da população, aplicando um procedimento de mutação. O cálculo da mudança de energia ΔE e a aceitação das novas soluções são feitos em cada iteração. Vamos analisar a Iteração 1.
Iteração 1:
Temperatura inicial: T0 = 36.061
Solução Pop id: 0
A solução atual é x0 = [-0.74, 1.25]. Para cada dimensão, uma mutação é aplicada, gerando uma nova solução vizinha:
- Dimensão 0: A solução vizinha para a dimensão 0 foi calculada com uma média de -0.74 e um desvio padrão de 0.74, gerando um valor de vizinhança de -0.8532.
- Dimensão 1: Para a dimensão 1, a média foi 1.25 e o desvio padrão foi 1.25, gerando a solução vizinha de 2.2558.
A nova solução gerada foi: xnew = [-0.8532, 2.2558]
O valor da função objetivo para essa nova solução foi of = 5.8166 e o valor de energia foi ΔE = 3.7065. A probabilidade de aceitação foi calculada usando a equação P(ΔE, T), que resulta em P = 0.9023. Um número aleatório foi gerado (0.4151), que foi menor que a probabilidade de aceitação, então a nova solução foi aceita.
Solução Pop id: 1
A solução atual é x1 = [3.58, -3.33]. Novamente, uma mutação foi aplicada:
- Dimensão 0: A solução vizinha para a dimensão 0 foi calculada com uma média de 3.58 e um desvio padrão de 3.58, gerando um valor de vizinhança de 3.4979.
- Dimensão 1: Para a dimensão 1, a média foi -3.33 e o desvio padrão foi 3.33, gerando a solução vizinha de -3.7698.
A nova solução gerada foi: xnew = [3.4979, -3.7698]
O valor da função objetivo foi of = 26.4471 e a energia foi ΔE = 2.5418. A probabilidade de aceitação foi P = 0.9319. O número aleatório gerado (0.479) foi menor que a probabilidade, então a nova solução também foi aceita.
Atualização das Soluções
Após a aceitação das novas soluções, a população é atualizada com as melhores soluções encontradas:
- \(x_0 = [-0.8532, 2.2558]\), \(of_{pop} = 5.8166\) - melhor solução
- \(x_1 = [3.4979, -3.7698]\), \(of_{pop} = 26.4471\)
O processo irá continuar iterativamente até que a temperatura seja suficientemente baixa ou até que o algoritmo tenha convergido para uma solução de energia mínima.
Reference list
| ID | Reference |
|---|---|
[1] | |
[2] | |
[3] | |
[4] | |
[5] | |
[6] | |
[7] | |
[8] | |
[9] |