模拟退火算法

模拟退火算法是一种基于概率的优化算法，其优点就是能够有效解决NP难问题，且能够避免陷入局部最优，且初值的选择对寻优结果没有太大的影响，接下来就为大家简单介绍一下模拟退火算法的基本原理。

此算法的提出最初是受到了固体由高温到低温的过程中，其内部分子状态及内部能量的变化规律的启发，因此得名叫做模拟退火算法。

在模拟退火算法中，解对应于固体的某一状态，函数值对应于固体处于某一状态时所具有的内能。在现实生活中固体随着温度的降低更趋向于处于一种低内能的状态。这便是模拟退火算法进行优化的原理。算法以固体所处的温度T为控制参数，随着T的下降使固体内能（目标函数值）也逐渐下降，直至趋于全局最小。故而可以看出，模拟退火算法是用来解决最小化问题的。

接下来介绍一些模拟退火算法中用到的符号：

：随机变量，代表内能

：固体所处的某一状态（优化问题的一个解）

：固体处于状态时所具有的内能（目标函数值）

：状态的集合

：中状态的个数

：固体所处的温度

，其中

上式表示固体所处温度为时，固体为状态（或固体内能为）的概率，其中为Boltzmann常数。

接下来先给出四个命题，并加以证明：

命题一：当固体所处温度一定时，系统处于低内能状态的概率大于系统处于高内能状态的概率。

证明：设，即固体处于状态2时的内能大于系统处于状态1时的内能

，问题得证。

命题二：当固体所处温度足够高时（很大时），固体处于每一种状态的概率基本相同，都等于。

证明：

得证。

命题三：若固体在某一温度下有两种或两种以上不同内能的状态，则处于最低内能状态的概率一定大于。

证明：设为最低内能状态。现假设，则根据命题一可知，，又因为存在且，则可以得知，然而这是不可能事件。因此假设不成立，原命题成立。

命题四：当固体温度趋于0时，固体处于最低内能状态的概率趋近于1。

证明：设为固体最低内能状态，为固体最低内能。

得证。

根据命题二可知，固体温度足够高时，由于每一种状态的概率都差不多，所以固体由低内能状态转为高内能状态也是很有可能发生的（算法中跳出局部最优解的前提条件）；由命题一中推导公式不难得知，随着固体所处温度的下降，处于高内能状态的概率与处于低内能状态的概率比值接近于0，即这两种状态的概率差距在拉大，因此，由低内能状态转向高内能状态的可能性也在下降。

算法步骤：

第一步：设定初始温度并令当前温度，以及设置终止温度并随机生成一个初始解（初始状态），并令当前状态

第二步：令，取值往往在0.5到0.99之间，根据具体问题而定

第三步：，计算当前状态的内能（目标函数值），根据当前状态（解）进行扰动，产生一个新状态，并计算新状态内能（目标函数值），接着计算

第四步：若，则新状态被接受，令；若，则以概率接收新状态，具体操作为随机产生一个0到1之间的随机数，如果，则接收新状态，令

第五步：在当前温度下，进行次扰动和接收过程，即循环执行第三步和第四步

第六步：判断是否小于，是则终止算法，否则跳到第二步继续执行

最后可以看出，算法实际上是两层循环嵌套，外层循环控制温度，内层循环来进行扰动产生新解。由于算法初始温度较高，这样，由一个低内能状态（目标函数值小的解或称之为较优解）转向一个高内能状态（目标函数值大的解或称之为较劣解）的可能性还是很大的，因此也很有可能调出局部最优解，之前的分析中对这一点也进行了相关说明。随着温度逐渐降低，算法最终由可能收敛到全局最优，这里说有可能的原因是因为，在温度很低时，虽然从地内能状态跳到高内能状态的可能性不大，但是也有可能发生。

以上便是对模拟退火算法的简要介绍，接下来附上本人用python编写的模拟退火的相关程序：

import numpy as np
from numpy import random as rd
import matplotlib.pyplot as plt
np.set_printoptions(suppress=True)


def objective_func(x, y):
    # 这是一个用来测试的目标函数
    return -2 * np.pi * np.exp(-0.2 * ((x ** 2 + y ** 2) / 2) ** 0.5) - np.exp(0.5 * (np.cos(2 * np.pi * x) + np.cos(2 * np.pi * y))) + 2 * np.pi


def objective_func2(x):
    # 这也是一个用来测试的目标函数
    return x ** 2


class SA(object):

    def __init__(self, alpha, T_0, T_f, obj_func, inner_loop_times, **kwargs):
        """

        :param alpha: 温度T的衰减系数，取值范围最好在0.8至0.99之间
        :param T_0: 初始温度
        :param T_f: 终止温度
        :param obj_func: 自定义的目标函数，要求求最小值
        :param kwargs: 决策变量取值范围，传入方式为'x1=(x1_min, x1_max), x2=(x2_min, x2_max), ...'，
                        此处决策变量的传入顺序需与目标函数中决策变量参数对应一致，例如def obj_func(x1, x2, ...)
        :param inner_loop_times: 算法中，内层循环的循环迭代次数
        """
        self.alpha = alpha
        self.T_f = T_f
        self.obj_func = obj_func
        self.kwargs = kwargs
        self.inner_loop_times = inner_loop_times
        self.T_0 = T_0
        # 确定初始状态
        self.x_i = [rd.uniform(values[0], values[1]) for values in self.kwargs.values()]
        # 以生成的初始状态（第一个解）作为目前为止最优状态self.best_x，其目标函数值作为目前为止最优函数值self.best_E
        self.best_x = self.x_i
        self.best_E = self.obj_func(*self.x_i)
        self.E_s = []

    def run(self):
        # 第一层循环控制温度
        while True:
            # 第二层循环
            for i in range(self.inner_loop_times):
                # 这个循环是为了生成与上一次状态不一样的状态，如果与上一次状态一样则重新生成
                while True:
                    self.x_j = [rd.uniform(values[0], values[1]) for values in self.kwargs.values()]
                    if self.x_j != self.x_i:
                        break
                E_i = self.obj_func(*self.x_i)
                E_j = self.obj_func(*self.x_j)
                delta_E = E_j - E_i
                if delta_E <= 0:
                    self.x_i = self.x_j
                else:
                    if rd.random() < np.exp(-delta_E / self.T_0):
                        self.x_i = self.x_j
                current_E = self.obj_func(*self.x_i)
                # 更新历史最优状态（解），及其对应的能量（函数值）
                if current_E < self.best_E:
                    self.best_E = current_E
                    self.best_x = self.x_i
            self.E_s.append(E_i)
            self.T_0 *= self.alpha
            if self.T_0 < self.T_f:
                break
        print("模拟退火过程中的历史最优解：\n", "解：", self.best_x)
        print("目标函数值：", self.best_E)
        print("模拟退火迭代最终确定的解：\n", "解：", self.x_i)
        print("目标函数值：", E_i)

    def draw(self):
        fig = plt.figure()
        ax = plt.subplot(1, 1, 1)
        ax.scatter(np.arange(1, len(self.E_s) + 1), self.E_s, marker="*", edgecolor="r", label="E trend", color="", alpha=0.5)
        ax.set_xlabel("temprature decline times")
        ax.set_ylabel("E")
        ax.legend()
        plt.show()


def main():
    T_f_0 = "0.01"
    T_f_1 = "0.000000001"
    T_f_2 = "0.000000000000000001"
    print("========第一个测试目标函数==============")
    print("---------T_f=%s----------" % T_f_0)
    sa = SA(0.99, 1000, float(T_f_0), objective_func, 200, x1=(-10, 5), x2=(-10, 5))
    sa.run()
    sa.draw()
    print("---------T_f=%s----------" % T_f_1)
    sa1 = SA(0.99, 1000, float(T_f_1), objective_func, 200, x1=(-30, 30), x2=(-30, 30))
    sa1.run()
    sa1.draw()
    print("========第二个测试目标函数==============")
    # 测试当温度足够低时以概率1趋近与能量最低（目标函数值最小）,即温度足够低时，模拟退火最终确定的解和历史最优解总是相同，
    # 如果模拟退火最终确定的解和历史最优不同，说明温度还不够低，还是有很大的概率使新的次好状态取代当前较好的状态
    print("---------T_f=%s---------" % T_f_2)
    sa2 = SA(0.99, 1000, float(T_f_2), objective_func2, 200, x1=(-30, 30))
    sa2.run()
    sa2.draw()


if __name__ == "__main__":
    main()