【翻译自 ： Simulated Annealing From Scratch in Python】

       【说明：Jason Brownlee PhD大神的文章个人很喜欢，所以闲暇时间里会做一点翻译和学习实践的工作，这里是相应工作的实践记录，希望能帮到有需要的人！】

       模拟退火是一种随机全局搜索优化算法。这意味着它将随机性作为搜索过程的一部分。这使得该算法适用于非线性目标函数，而其他局部搜索算法不能很好地运行。像随机爬山局部搜索算法一样，它修改单个解决方案并搜索搜索空间的相对局部区域，直到找到局部最优值为止。与爬山算法不同，它可能接受较差的解决方案作为当前的工作解决方案。接受更差解的可能性在搜索开始时就开始很高，并且随着搜索的进行而降低，从而使算法有机会首先为全局最优定位区域，避开局部最优，然后爬升到最优本身。

         在本教程中，您将发现用于功能优化的模拟退火优化算法。完成本教程后，您将知道：

模拟退火是用于函数优化的随机全局搜索算法。
如何在Python中从头开始实现模拟退火算法。
如何使用模拟退火算法并检查算法结果。

 

教程概述

       本教程分为三个部分： 他们是：

模拟退火
实施模拟退火
模拟退火工作示例

模拟退火

        模拟退火是一种随机全局搜索优化算法。该算法的灵感来自于冶金学中的退火工艺，该工艺将金属快速加热到高温，然后缓慢冷却，从而提高了强度并使其更易于使用。退火过程的工作原理是：首先在高温下激发材料中的原子，使原子大量移动，然后缓慢降低其兴奋性，使原子陷入新的，更稳定的结构。

        可以将模拟退火优化算法视为随机爬山的改进版本。随机爬山维护单个候选解，并在搜索空间中从候选中采取随机但受限制的大小的步骤。 如果新点比当前点好，则将当前点替换为新点。 此过程将继续进行固定数量的迭代。模拟退火以相同的方式执行搜索。 主要区别在于有时会接受不如当前点（差点）的新点。概率上接受一个更差的点，其中接受比当前解决方案更差的解决方案的可能性取决于搜索温度以及解决方案比当前解决方案差多少。

        搜索的初始温度作为超参数提供，并且随着搜索的进行而降低。尽管通常根据迭代次数来计算温度，但是可以使用许多不同的方案（退火计划）将搜索过程中的温度从初始值降低到非常低的值。计算温度的一个流行示例是所谓的“快速模拟退火”，其计算如下

                                         温度=初始温度/（迭代次数+ 1）
       在迭代次数从零开始的情况下，我们将迭代次数加1，以避免除以零误差。较差解的接受使用温度以及较差解和当前解的目标函数评估之间的差。使用该信息计算出一个介于0和1之间的值，表示接受较差解的可能性。然后使用随机数对该分布进行采样，如果该随机数小于该值，则表示可接受较差的解决方案。

       这被称为都市接受标准，并且为了最小化，其计算如下：

                                          criterion = exp( -(objective(new) – objective(current)) / temperature)          
          其中exp（）是e（数学常数），它被提高为提供的自变量的幂，而objective（new）和objective（current）是新（更差）和当前候选解决方案的目标函数评估。这样做的结果是，较差的解决方案更有可能在搜索的早期被接受，而在搜索的后期被接受的可能性较小。 目的是在搜索开始时的高温将帮助搜索找到全局最优值的盆地，而在搜索的后期，低温将帮助算法完善全局最优值。

       现在我们已经熟悉了模拟退火算法，让我们看看如何从头开始实现它。

实施模拟退火

         在本节中，我们将探索如何从头开始实现模拟退火优化算法。首先，我们必须定义目标函数和每个输入变量到目标函数的界限。 目标函数只是一个Python函数，我们将其命名为Objective（）。 边界将是一个2D数组，每个输入变量都具有一个维度，该变量定义了变量的最小值和最大值。例如，一维目标函数和界限将定义如下：

# objective function
def objective(x):
	return 0

# define range for input
bounds = asarray([[-5.0, 5.0]])

          接下来，我们可以将初始点生成为问题范围内的随机点，然后使用目标函数对其进行评估。

# generate an initial point
best = bounds[:, 0] + rand(len(bounds)) * (bounds[:, 1] - bounds[:, 0])
# evaluate the initial point
best_eval = objective(best)

         我们需要维护“当前”解决方案，这是搜索的重点，并且可以用更好的解决方案代替。

# current working solution
curr, curr_eval = best, best_eval

       现在我们可以遍历定义为“ n_iterations”的算法的预定义迭代次数，例如100或1,000。

# run the algorithm
for i in range(n_iterations):

       算法迭代的第一步是从当前的工作解决方案中生成新的候选解决方案，例如 迈出一步。这需要预定义的“ step_size”参数，该参数相对于搜索空间的边界。 我们将采用高斯分布的随机步骤，其中均值是我们的当前点，标准偏差由“ step_size”定义。 这意味着大约99％的步骤将在当前点的3 * step_size之内。

# take a step
candidate = solution + randn(len(bounds)) * step_size

        我们不必采取这种方式。 您可能希望使用0到步长之间的均匀分布。 例如：

# take a step
candidate = solution + rand(len(bounds)) * step_size

        接下来，我们需要对其进行评估。

# evaluate candidate point
candidte_eval = objective(candidate)

        然后，我们需要检查此新点的评估结果是否等于或优于当前最佳点，如果是，则用此新点替换当前最佳点。这与当前的工作解决方案不同，后者是搜索的重点。

# check for new best solution
if candidate_eval < best_eval:
	# store new best point
	best, best_eval = candidate, candidate_eval
	# report progress
	print('>%d f(%s) = %.5f' % (i, best, best_eval))

        接下来，我们需要准备替换当前的工作解决方案。第一步是计算当前解决方案和当前工作解决方案的目标函数评估之间的差异。

# difference between candidate and current point evaluation
diff = candidate_eval - curr_eval

        接下来，我们需要使用快速退火程序来计算当前温度，其中“ temp”是作为参数提供的初始温度。

# calculate temperature for current epoch
t = temp / float(i + 1)

         然后，我们可以计算接受性能比当前工作解决方案差的解决方案的可能性。

# calculate metropolis acceptance criterion
metropolis = exp(-diff / t)

         最后，如果目标函数评估更好（差异为负）或目标函数较差，我们可以接受新点作为当前工作解决方案，但是我们有可能决定接受它。

# check if we should keep the new point
if diff < 0 or rand() < metropolis:
	# store the new current point
	curr, curr_eval = candidate, candidate_eval

        就是这样。我们可以将该模拟退火算法实现为可重用函数，该函数以目标函数的名称，每个输入变量的界限，总迭代次数，步长和初始温度为参数，并返回找到的最佳解及其评估。

# simulated annealing algorithm
def simulated_annealing(objective, bounds, n_iterations, step_size, temp):
	# generate an initial point
	best = bounds[:, 0] + rand(len(bounds)) * (bounds[:, 1] - bounds[:, 0])
	# evaluate the initial point
	best_eval = objective(best)
	# current working solution
	curr, curr_eval = best, best_eval
	# run the algorithm
	for i in range(n_iterations):
		# take a step
		candidate = curr + randn(len(bounds)) * step_size
		# evaluate candidate point
		candidate_eval = objective(candidate)
		# check for new best solution
		if candidate_eval < best_eval:
			# store new best point
			best, best_eval = candidate, candidate_eval
			# report progress
			print('>%d f(%s) = %.5f' % (i, best, best_eval))
		# difference between candidate and current point evaluation
		diff = candidate_eval - curr_eval
		# calculate temperature for current epoch
		t = temp / float(i + 1)
		# calculate metropolis acceptance criterion
		metropolis = exp(-diff / t)
		# check if we should keep the new point
		if diff < 0 or rand() < metropolis:
			# store the new current point
			curr, curr_eval = candidate, candidate_eval
	return [best, best_eval]

               现在，我们知道了如何在Python中实现模拟退火算法，让我们看看如何使用它来优化目标函数。

模拟退火工作示例

        在本节中，我们将模拟退火优化算法应用于目标函数。首先，让我们定义目标函数。我们将使用一个简单的一维x ^ 2目标函数，其边界为[-5，5]。下面的示例定义了函数，然后为输入值的网格创建了函数响应面的折线图，并用红线标记了f（0.0）= 0.0的最优值

# convex unimodal optimization function
from numpy import arange
from matplotlib import pyplot

# objective function
def objective(x):
	return x[0]**2.0

# define range for input
r_min, r_max = -5.0, 5.0
# sample input range uniformly at 0.1 increments
inputs = arange(r_min, r_max, 0.1)
# compute targets
results = [objective([x]) for x in inputs]
# create a line plot of input vs result
pyplot.plot(inputs, results)
# define optimal input value
x_optima = 0.0
# draw a vertical line at the optimal input
pyplot.axvline(x=x_optima, ls='--', color='red')
# show the plot
pyplot.show()

          运行示例将创建目标函数的折线图，并清晰地标记函数的最优值。

         在将优化算法应用于问题之前，让我们花一点时间来更好地理解接受标准。首先，快速退火计划是迭代次数的指数函数。 通过为每次算法迭代创建温度图可以使这一点变得清晰。我们将使用10和100次算法迭代的初始温度，两者都是任意选择的。下面列出了完整的示例。

# explore temperature vs algorithm iteration for simulated annealing
from matplotlib import pyplot
# total iterations of algorithm
iterations = 100
# initial temperature
initial_temp = 10
# array of iterations from 0 to iterations - 1
iterations = [i for i in range(iterations)]
# temperatures for each iterations
temperatures = [initial_temp/float(i + 1) for i in iterations]
# plot iterations vs temperatures
pyplot.plot(iterations, temperatures)
pyplot.xlabel('Iteration')
pyplot.ylabel('Temperature')
pyplot.show()

        运行示例将计算每次算法迭代的温度，并创建算法迭代（x轴）与温度（y轴）的关系图。

        我们可以看到温度迅速地，指数地而不是线性地下降，因此在20次迭代之后，温度低于1，并且在其余的搜索过程中保持较低的温度。

        接下来，我们可以更好地了解大都市的接受标准如何随时间变化。回想一下，该标准是温度的函数，但也是新点的客观评估与当前工作解决方案相比有多大差异的函数。 这样，我们将为几个不同的“目标函数值差异”绘制标准，以查看其对接受概率的影响。下面列出了完整的示例。

# explore metropolis acceptance criterion for simulated annealing
from math import exp
from matplotlib import pyplot
# total iterations of algorithm
iterations = 100
# initial temperature
initial_temp = 10
# array of iterations from 0 to iterations - 1
iterations = [i for i in range(iterations)]
# temperatures for each iterations
temperatures = [initial_temp/float(i + 1) for i in iterations]
# metropolis acceptance criterion
differences = [0.01, 0.1, 1.0]
for d in differences:
	metropolis = [exp(-d/t) for t in temperatures]
	# plot iterations vs metropolis
	label = 'diff=%.2f' % d
	pyplot.plot(iterations, metropolis, label=label)
# inalize plot
pyplot.xlabel('Iteration')
pyplot.ylabel('Metropolis Criterion')
pyplot.legend()
pyplot.show()

         运行示例将使用每次迭代显示的温度（在上一节中显示）来计算每次算法迭代的都市接受标准。该图针对新的更差解与当前有效解之间的三个差异具有三条线。我们可以看到，解决方案越差（差异越大），则与我们预期的一样，无论算法迭代如何，模型接受较差解决方案的可能性越小。 我们还可以看到，在所有情况下，接受更差解的可能性都随算法迭代而降低。

         现在，我们更加了解温度和都会接受标准随时间的行为，让我们将模拟退火应用于我们的测试问题。首先，我们将播种伪随机数生成器。通常这不是必需的，但是在这种情况下，我想确保每次运行算法时都得到相同的结果（相同的随机数序列），以便稍后绘制结果。

# seed the pseudorandom number generator
seed(1)

          接下来，我们可以定义搜索的配置。在这种情况下，我们将搜索算法的1,000次迭代，并使用0.1的步长。 假设我们使用的是高斯函数来生成步长，这意味着大约99％的所有步长将位于给定点（0.1 * 3）的距离内，例如 三个标准差。我们还将使用10.0的初始温度。 搜索过程比初始温度对退火时间表更敏感，因此初始温度值几乎是任意的。

n_iterations = 1000
# define the maximum step size
step_size = 0.1
# initial temperature
temp = 10

        接下来，我们可以执行搜索并报告结果。

# perform the simulated annealing search
best, score = simulated_annealing(objective, bounds, n_iterations, step_size, temp)
print('Done!')
print('f(%s) = %f' % (best, score))

      完整实例如下：

# simulated annealing search of a one-dimensional objective function
from numpy import asarray
from numpy import exp
from numpy.random import randn
from numpy.random import rand
from numpy.random import seed

# objective function
def objective(x):
	return x[0]**2.0

# simulated annealing algorithm
def simulated_annealing(objective, bounds, n_iterations, step_size, temp):
	# generate an initial point
	best = bounds[:, 0] + rand(len(bounds)) * (bounds[:, 1] - bounds[:, 0])
	# evaluate the initial point
	best_eval = objective(best)
	# current working solution
	curr, curr_eval = best, best_eval
	# run the algorithm
	for i in range(n_iterations):
		# take a step
		candidate = curr + randn(len(bounds)) * step_size
		# evaluate candidate point
		candidate_eval = objective(candidate)
		# check for new best solution
		if candidate_eval < best_eval:
			# store new best point
			best, best_eval = candidate, candidate_eval
			# report progress
			print('>%d f(%s) = %.5f' % (i, best, best_eval))
		# difference between candidate and current point evaluation
		diff = candidate_eval - curr_eval
		# calculate temperature for current epoch
		t = temp / float(i + 1)
		# calculate metropolis acceptance criterion
		metropolis = exp(-diff / t)
		# check if we should keep the new point
		if diff < 0 or rand() < metropolis:
			# store the new current point
			curr, curr_eval = candidate, candidate_eval
	return [best, best_eval]

# seed the pseudorandom number generator
seed(1)
# define range for input
bounds = asarray([[-5.0, 5.0]])
# define the total iterations
n_iterations = 1000
# define the maximum step size
step_size = 0.1
# initial temperature
temp = 10
# perform the simulated annealing search
best, score = simulated_annealing(objective, bounds, n_iterations, step_size, temp)
print('Done!')
print('f(%s) = %f' % (best, score))

        运行示例将报告搜索进度，包括迭代次数，函数输入以及每次检测到改进时目标函数的响应。

        搜索结束时，找到最佳解决方案，并报告其评估结果。

        注意：由于算法或评估程序的随机性，或者数值精度的差异，您的结果可能会有所不同。 考虑运行该示例几次并比较平均结果。

        在这种情况下，我们可以看到在算法的1,000次迭代中约有20处改进，并且该解决方案非常接近于0.0的最佳输入，即f（0.0）= 0.0。

>34 f([-0.78753544]) = 0.62021
>35 f([-0.76914239]) = 0.59158
>37 f([-0.68574854]) = 0.47025
>39 f([-0.64797564]) = 0.41987
>40 f([-0.58914623]) = 0.34709
>41 f([-0.55446029]) = 0.30743
>42 f([-0.41775702]) = 0.17452
>43 f([-0.35038542]) = 0.12277
>50 f([-0.15799045]) = 0.02496
>66 f([-0.11089772]) = 0.01230
>67 f([-0.09238208]) = 0.00853
>72 f([-0.09145261]) = 0.00836
>75 f([-0.05129162]) = 0.00263
>93 f([-0.02854417]) = 0.00081
>144 f([0.00864136]) = 0.00007
>149 f([0.00753953]) = 0.00006
>167 f([-0.00640394]) = 0.00004
>225 f([-0.00044965]) = 0.00000
>503 f([-0.00036261]) = 0.00000
>512 f([0.00013605]) = 0.00000
Done!
f([0.00013605]) = 0.000000

          以线图的形式查看搜索的进度可能很有趣，该线图显示了每次改进后最佳解决方案的评估变化。

          每当有改进时，我们就可以更新simulated_annealing（）来跟踪目标函数评估，并返回此分数列表。

# simulated annealing algorithm
def simulated_annealing(objective, bounds, n_iterations, step_size, temp):
	# generate an initial point
	best = bounds[:, 0] + rand(len(bounds)) * (bounds[:, 1] - bounds[:, 0])
	# evaluate the initial point
	best_eval = objective(best)
	# current working solution
	curr, curr_eval = best, best_eval
	# run the algorithm
	for i in range(n_iterations):
		# take a step
		candidate = curr + randn(len(bounds)) * step_size
		# evaluate candidate point
		candidate_eval = objective(candidate)
		# check for new best solution
		if candidate_eval < best_eval:
			# store new best point
			best, best_eval = candidate, candidate_eval
			# keep track of scores
			scores.append(best_eval)
			# report progress
			print('>%d f(%s) = %.5f' % (i, best, best_eval))
		# difference between candidate and current point evaluation
		diff = candidate_eval - curr_eval
		# calculate temperature for current epoch
		t = temp / float(i + 1)
		# calculate metropolis acceptance criterion
		metropolis = exp(-diff / t)
		# check if we should keep the new point
		if diff < 0 or rand() < metropolis:
			# store the new current point
			curr, curr_eval = candidate, candidate_eval
	return [best, best_eval, scores]

         然后，我们可以创建这些分数的折线图，以查看搜索过程中发现的每个改进的目标函数的相对变化。

# line plot of best scores
pyplot.plot(scores, '.-')
pyplot.xlabel('Improvement Number')
pyplot.ylabel('Evaluation f(x)')
pyplot.show()

        结合在一起，下面列出了执行搜索并绘制搜索过程中改进解决方案的目标函数得分的完整示例。

# simulated annealing search of a one-dimensional objective function
from numpy import asarray
from numpy import exp
from numpy.random import randn
from numpy.random import rand
from numpy.random import seed
from matplotlib import pyplot

# objective function
def objective(x):
	return x[0]**2.0

# simulated annealing algorithm
def simulated_annealing(objective, bounds, n_iterations, step_size, temp):
	# generate an initial point
	best = bounds[:, 0] + rand(len(bounds)) * (bounds[:, 1] - bounds[:, 0])
	# evaluate the initial point
	best_eval = objective(best)
	# current working solution
	curr, curr_eval = best, best_eval
	scores = list()
	# run the algorithm
	for i in range(n_iterations):
		# take a step
		candidate = curr + randn(len(bounds)) * step_size
		# evaluate candidate point
		candidate_eval = objective(candidate)
		# check for new best solution
		if candidate_eval < best_eval:
			# store new best point
			best, best_eval = candidate, candidate_eval
			# keep track of scores
			scores.append(best_eval)
			# report progress
			print('>%d f(%s) = %.5f' % (i, best, best_eval))
		# difference between candidate and current point evaluation
		diff = candidate_eval - curr_eval
		# calculate temperature for current epoch
		t = temp / float(i + 1)
		# calculate metropolis acceptance criterion
		metropolis = exp(-diff / t)
		# check if we should keep the new point
		if diff < 0 or rand() < metropolis:
			# store the new current point
			curr, curr_eval = candidate, candidate_eval
	return [best, best_eval, scores]

# seed the pseudorandom number generator
seed(1)
# define range for input
bounds = asarray([[-5.0, 5.0]])
# define the total iterations
n_iterations = 1000
# define the maximum step size
step_size = 0.1
# initial temperature
temp = 10
# perform the simulated annealing search
best, score, scores = simulated_annealing(objective, bounds, n_iterations, step_size, temp)
print('Done!')
print('f(%s) = %f' % (best, score))
# line plot of best scores
pyplot.plot(scores, '.-')
pyplot.xlabel('Improvement Number')
pyplot.ylabel('Evaluation f(x)')
pyplot.show()

         运行示例将执行搜索，并像以前一样报告结果。创建一个线形图，显示在爬山搜索期间每个改进的目标函数评估。 在搜索过程中，我们可以看到目标函数评估大约有20个变化，其中初始变化较大，而随着搜索算法收敛于最优值，在搜索结束时变化很小至难以察觉的变化。