NEIGHBOUR(T, best) - дава следваща точка в зависимост от температурата и настоящата най-добра. Например точка от нормалното разпределение около точка best с дисперсия Т.
ACCEPT(T, dE) - определя с каква вероятност ще приемем по-лоша точка. $$\displaystyle ACCEPT(T, \Delta E) = \mathrm{e}^{\frac{-\Delta E}{kT}}$$
Explore $\displaystyle T \to \infty \hspace{1pc} \lim_{T \to \infty} \frac{-\Delta E}{kT} = 0 \hspace{2pc} \mathrm{e}^{0} = 1$
Exploit $\displaystyle T \to 0 \hspace{1pc} \lim_{T \to 0} \frac{-\Delta E}{kT} = -\infty \hspace{2pc} \mathrm{e}^{-\infty} \to 0$
Обща схема на стратегията:
In [1]:
import operator
import numpy as np
import math
from numpy import random
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import matplotlib.animation as animation
%matplotlib nbagg
class SimulatedAnnealing():
def __init__(self, minimize, T_min, T_max, init,
neighbour, energy, accept, cooling):
self.compare = operator.lt if minimize else operator.gt
self.T_min = T_min
self.neighbour = neighbour
self.energy = energy
self.accept = accept
self.cooling = cooling
self.T = T_max
self.best = init()
self.step = 0
def iteration(self):
next = self.neighbour(self.T, self.best)
deltaE = self.energy(next) - self.energy(self.best)
if self.compare(deltaE, 0):
self.best = next
elif random.random() < self.accept(self.T, deltaE):
self.best = next
def cool(self):
self.T = self.cooling(self.T, self.best)
def execute(self):
while self.T > self.T_min:
self.iteration()
self.cool()
self.log()
self.step += 1
print('\nf(' + str(self.best) + ') = ' + str(self.energy(self.best)))
return self.best
def log(self):
if(self.step % 100 == 0):
print(str(self.step) + " T="+str(self.T) +
"\tf(" + str(self.best) + ") = " + str(self.energy(self.best)))
# standart configuration parameters and functions for Simulated Annealing
t_min = 1e-5
t_max = 1
def init():
return np.array([0.1, 1.4])
def neighbour(T, x):
#return x + random.uniform(-T, T, x.shape)
return random.normal(x, max(T, 0.1), x.shape)
def accept(T, deltaE, k=0.1):
return math.exp(-(deltaE)/k/T)
def cooling(T, best):
return T*0.99
In [2]:
# plot the Rosenbrock function
def rosenbrock(x, y, a=1, b=100):
return (a - x)**2 + b*(y - x**2)**2
def rosenbrock_energy(vector):
return rosenbrock(*vector)
# Design variables at mesh points
i1 = np.arange(-2, 2, 0.01)
i2 = np.arange(-1, 3, 0.01)
X, Y = np.meshgrid(i1, i2)
Z = rosenbrock(X, Y)
fig = plt.figure()
ax = plt.axes(projection='3d')
Gx, Gy = np.gradient(Z) # gradients with respect to x and y
G = (Gx**2+Gy**2)**.5 # gradient magnitude
N = G/G.max() # normalize 0..1
ax.plot_surface(X, Y, Z,facecolors=cm.jet(N))
plt.show()
In [3]:
# minimize the Rosenbrock function with Simulated Annealing
rosenbrock_min = SimulatedAnnealing(True, t_min, t_max, init,
neighbour, rosenbrock_energy, accept, cooling).execute()
In [12]:
# knapsack problem data set generator
DIM = 60
MAX_VALUE = 100000000
MAX_WEIGHT = 100000000
C = 1000000000
W = np.random.randint(MAX_WEIGHT, size=DIM)
V = np.random.randint(MAX_VALUE, size=DIM)
print(W)
print(V)
In [6]:
# maximization of the 0-1 Knapsack Problem with data taken from
# https://people.sc.fsu.edu/~jburkardt/datasets/knapsack_01/knapsack_01.html
# DIM = 15
# W = [70, 73, 77, 80, 82, 87, 90, 94, 98, 106, 110, 113, 115, 118, 120]
# V = [135, 139, 149, 150, 156, 163, 173, 184, 192, 201, 210, 214, 221, 229, 240]
# C = 750
# OPTIMUM = 1458
# OPTIMAL_KNAPSACK = [1, 0, 1, 0, 1, 0, 1, 1, 1, 0, 0, 0, 0, 1, 1]
DIM = 24
W = [382745, 799601, 909247, 729069, 467902, 44328, 34610, 698150,
823460, 903959, 853665, 551830, 610856, 670702, 488960, 951111,
323046, 446298, 931161, 31385, 496951, 264724, 224916, 169684]
V = [825594, 1677009, 1676628, 1523970, 943972, 97426, 69666, 1296457,
1679693, 1902996, 1844992, 1049289, 1252836, 1319836, 953277,
2067538, 675367, 853655, 1826027, 65731, 901489, 577243, 466257, 369261]
C = 6404180
OPTIMUM = 13549094
OPTIMAL_KNAPSACK = [1, 1, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1]
def kp_init():
return np.random.randint(2, size=DIM)
def kp_neighbour(T, knapsack):
new_knapsack = np.copy(knapsack)
new_knapsack[random.randint(DIM)] ^= 1
while kp_cost(new_knapsack) == 0:
new_knapsack[random.randint(DIM)] ^= 1
return new_knapsack
def kp_cost(knapsack):
total_weight = sum([W[i] for i, chosen in enumerate(knapsack) if chosen])
total_value = sum([V[i] for i, chosen in enumerate(knapsack) if chosen])
return 0 if total_weight > C else total_value
def kp_accept(T, deltaE, k=0.1):
return math.exp(deltaE / k / T)
kp_max = SimulatedAnnealing(False, t_min, t_max, kp_init,
kp_neighbour, kp_cost, kp_accept, cooling).execute()
mismatches = 0
for i, item in enumerate(kp_max):
mismatches += 1 if kp_max[i] != OPTIMAL_KNAPSACK[i] else 0
print(str(mismatches) + " items in the selection mismatch with the optimal solution")
print("{0:.3f}% less optimal then the most optimal solution".format((1 - kp_cost(kp_max)/OPTIMUM) * 100))
In [5]:
# animation of Rosenbrock minimization
sa = SimulatedAnnealing(True, t_min, t_max, init, neighbour, rosenbrock_energy, accept, cooling)
fig = plt.figure()
xmin, xmax = -2.0, 2.0
ymin, ymax = -1.0, 3.0
x1 = np.linspace(xmin, xmax, 100)
x2 = np.linspace(ymin, ymax, 100)
Xc, Yc = np.meshgrid(x1, x2)
Zc = rosenbrock(Xc, Yc)
ax = fig.add_subplot(111, aspect='equal', autoscale_on=False, xlim=(xmin, xmax), ylim=(ymin, ymax))
CS = ax.contour(Xc, Yc, Zc, levels=[0, 0.1, 1, 2, 10, 20, 100])
ax.clabel(CS, inline=1, fmt='%1.1f', fontsize=14)
ax.set_title('Non-Convex Function')
line, = ax.plot([], [], 'o-', lw=2)
time_text = ax.text(0.02, 0.95, '', transform=ax.transAxes)
energy_text = ax.text(0.02, 0.90, '', transform=ax.transAxes)
def init_animation():
line.set_data([], [])
time_text.set_text('')
energy_text.set_text('')
return line, time_text, energy_text
xs = []
ys = []
def animate(step):
global sa, xs, ys
sa.iteration()
sa.cool()
xs.append(sa.best[0])
ys.append(sa.best[1])
line.set_data(xs, ys)
time_text.set_text('time = %.1f' % step)
energy_text.set_text('temperature = %.3f J' % sa.T)
return line, time_text, energy_text
ani = animation.FuncAnimation(fig, animate, frames=25, interval=25, init_func=init_animation)
plt.show()
In [7]:
def spike_pit(x, y):
return 21.5 + x*np.sin(4*math.pi*x) + y*np.sin(20*math.pi*y)
xs = np.arange(-3, 12.1, 0.01)
ys = np.arange(4.1, 5.8, 0.01)
# xs = np.arange(10, 40, 0.01)
# ys = np.arange(-10, 10, 0.01)
X, Y = np.meshgrid(xs, ys)
Z = spike_pit(X, Y)
fig = plt.figure()
ax = plt.axes(projection='3d')
Gx, Gy = np.gradient(Z) # gradients with respect to x and y
G = (Gx**2+Gy**2)**.5 # gradient magnitude
N = G/G.max() # normalize 0..1
ax.plot_surface(X, Y, Z,facecolors=cm.jet(N))
# ax.plot_surface(X, Y, Z)
plt.show()
In [8]:
# minimize the spike pit
# https://www.wolframalpha.com/input/?i=minimize+21.5+%2B+x*sin(4*pi*x)+%2B+y*sin(20*pi*y)
def spike_init():
return np.array([30.0, 0.0])
def spike_energy(vector):
return spike_pit(*vector)
spike_pit_min = SimulatedAnnealing(False, t_min, t_max, spike_init,
neighbour, spike_energy, kp_accept, cooling).execute()
In [9]:
# animation of Spike Pit maximization
sa = SimulatedAnnealing(False, t_min, t_max, spike_init,
neighbour, spike_energy, kp_accept, cooling)
fig = plt.figure()
xmin, xmax = 10, 40
ymin, ymax = -10, 10
x1 = np.linspace(xmin, xmax, 100)
x2 = np.linspace(ymin, ymax, 100)
Xc, Yc = np.meshgrid(x1, x2)
Zc = spike_pit(Xc, Yc)
ax = fig.add_subplot(111, aspect='equal', autoscale_on=False, xlim=(xmin, xmax), ylim=(ymin, ymax))
CS = ax.contour(Xc, Yc, Zc, levels=[20, 30, 40, 50, 60])
ax.clabel(CS, inline=1, fmt='%1.1f', fontsize=14)
ax.set_title('Trigonometric Spikes')
line, = ax.plot([], [], 'o-', lw=2)
time_text = ax.text(0.02, 0.95, '', transform=ax.transAxes)
energy_text = ax.text(0.02, 0.90, '', transform=ax.transAxes)
def init_animation():
line.set_data([], [])
time_text.set_text('')
energy_text.set_text('')
return line, time_text, energy_text
xs = []
ys = []
def animate(step):
global sa, xs, ys
sa.iteration()
sa.cool()
sa.log()
sa.step += 1
xs.append(sa.best[0])
ys.append(sa.best[1])
line.set_data(xs, ys)
time_text.set_text('time = %.1f' % step)
energy_text.set_text('temperature = %.3f J' % sa.T)
return line, time_text, energy_text
ani = animation.FuncAnimation(fig, animate, frames=25, interval=25, init_func=init_animation)
plt.show()
In [ ]: