python学智能算法（二）|模拟退火算法：进阶分析

【1】引言

前序学习进程中，已经对模拟退火算法进行了原理解释和最小值求解学习，相关文章链接为：

python学智能算法（一）|模拟退火算法：原理解释和最小值求解-CSDN博客

在此基础上，今天将进一步学习具体的数值处理退火过程，实现灵活运用该算法进行各式各样目标值求解。

【2】代码设计

在原理解释和最小值求解学习文章中，主要研究方向是把模拟退火原理用代码展示出来，这里面涉及到的操作步骤满足：

图1   模拟退火操作步骤

实际上，初始温度和最低温度都不参与计算，它们只是辅助决定要计算多少次。

真正实施的计算过程，是在初始温度时，从初始解开始，对初始解周围的区域进行搜索。当搜索次数达到预设数量时，温度降低。

在这个过程中，温度不参与计算，只是相当于一个条件，只要温度还没有到最低，搜索过程就要进行。

类似于温度就是一个跑道，计算过程就是粒子在奔跑。只要这个跑道还有直径，里面的粒子就要绕着跑道运动，且每个直径的跑道上，粒子运动的圈数相同。

图2   模拟退火说明

在此分析基础上，仅仅获得最有分析结果不再是唯一目标。

将代码搜寻过程中，选取的最优解和最优目标值同时展示出来，会更有助于更清晰地展示退火过程。

在前述文章中，算法运行的目标是获得最小值，使用的模拟退火算法核心代码为：

while temperature > final_temperature:
    for _ in range(num_iterations_per_temp):
        # 在当前解的邻域内生成新解
        new_solution = current_solution + np.random.uniform(-0.1, 0.1)
        # 计算新解的目标函数值
        new_energy = objective_function(new_solution)
        # 计算能量差
        delta_energy = new_energy - current_energy

        # 如果新解更优，直接接受
        if delta_energy < 0:
            current_solution = new_solution
            current_energy = new_energy
            # 更新最优解
            if new_energy < best_energy:
                best_solution = new_solution
                best_energy = new_energy
        # 如果新解更差，以一定概率接受
        else:
            acceptance_probability = np.exp(-delta_energy / temperature)
            if np.random.rand() < acceptance_probability:
                current_solution = new_solution
                current_energy = new_energy

这个计算过程，实际上可以取到的参数包括：

每次计算的解

每次计算的目标函数值/能量值

每次计算的最优解

每次计算的最优目标函数值/能量值

所以，如果将上述量尽可能展示出来，具体的退火过程会更为清晰。

【3】代码测试

要想展示各个具体的解和对应能量值，只需要定义相关的列表存储各个值即可，在这依然是先修改核心代码：

# 模拟退火算法
def simulated_annealing(initial_solution, initial_temperature, final_temperature, alpha, num_iterations_per_temp):
    # 当前解
    current_solution = initial_solution
    # 当前解的目标函数值
    current_energy = objective_function(current_solution)
    # 最优解
    best_solution = current_solution
    # 最优解的目标函数值
    best_energy = current_energy
    # 当前温度
    temperature = initial_temperature

    # 用于记录每一步的最优目标函数值
    energy_history = [best_energy]
    current_solution_history = [current_solution]
    current_energy_history = [current_energy]
    best_energy_history = [best_energy]
    best_solution_history =[best_solution]

    while temperature > final_temperature:
        for _ in range(num_iterations_per_temp):
            # 在当前解的邻域内生成新解
            new_solution = current_solution + np.random.uniform(-0.1, 0.1)
            # 计算新解的目标函数值
            new_energy = objective_function(new_solution)
            # 计算能量差
            delta_energy = new_energy - current_energy

            # 如果新解更优，直接接受
            if delta_energy < 0:
                current_solution = new_solution
                current_energy = new_energy
                # 更新最优解
                if new_energy < best_energy:
                    best_solution = new_solution
                    best_energy = new_energy
            # 如果新解更差，以一定概率接受
            else:
                acceptance_probability = np.exp(-delta_energy / temperature)
                if np.random.rand() < acceptance_probability:
                    current_solution = new_solution
                    current_energy = new_energy

            # 记录当前最优目标函数值
            energy_history.append(best_energy)
            # 更新当前解和能量的历史记录
            current_solution_history.append(current_solution)
            current_energy_history.append(current_energy)
            best_energy_history.append(best_energy)
            best_solution_history.append(best_solution)

        # 降低温度
        temperature *= alpha

    return best_solution, best_energy, energy_history, current_solution_history, current_energy_history, best_energy_history, best_solution_history

非常清晰，这里只是在两处位置分别新增了四行代码：

# 用于记录每一步的最优目标函数值
energy_history = [best_energy]
current_solution_history = [current_solution]
current_energy_history = [current_energy]
best_energy_history = [best_energy]
best_solution_history =[best_solution]

...

# 记录当前最优目标函数值
energy_history.append(best_energy)
# 更新当前解和能量的历史记录
current_solution_history.append(current_solution)
current_energy_history.append(current_energy)
best_energy_history.append(best_energy)
best_solution_history.append(best_solution)

新增的代码对运算没有任何影响，只是将参数值用列表的形式存储起来。

之后就比较简单，用动画的形式动态展示这些参数的关系就可以：

fig, ax = plt.subplots()  # 定义要画图
# 设置坐标轴范围
#ax.set_xlim(min(current_solution_history), max(current_solution_history))
#ax.set_ylim(min(current_energy_history), max(current_energy_history))
ax.set_xlim(-10,10)
ax.set_ylim(0,100)
# 定义散点图对象
scatter1 = ax.scatter([], [], c='b', label='当前解与当前能量关系')
scatter2 = ax.scatter([], [], c='g', marker="^",label='最优解与最优能量关系')
# 定义折线图对象来表示 best_energy 的变化
line, = ax.plot([], [], 'r-', label='当前解与最优能量关系')

num_frames = len(energy_history)
x = np.linspace(-10, 10, 400)
y = objective_function(x)
plt.plot(x, y, label='目标函数')
def update(frame):
    x = current_solution_history[:frame + 1]
    y = current_energy_history[:frame + 1]
    # 更新散点图的数据
    scatter1.set_offsets(np.c_[x, y])

    # 更新折线图的数据，这里使用迭代步数作为 x 轴，best_energy 作为 y 轴
    line_x = current_solution_history[:frame + 1]
    line_y = best_energy_history[:frame + 1]
    line.set_data(line_x, line_y)

    scatter_x=best_solution_history[:frame + 1]
    scatter_y =best_energy_history[:frame + 1]
    scatter2.set_offsets(np.c_[scatter_x, scatter_y])

    return scatter1, line,scatter2

ani = animation.FuncAnimation(fig, update, repeat=True,
                              frames=num_frames, interval=50)  # 调用animation.FuncAnimation()函数画图
plt.legend()

在这段代码中，分别定义了4个图像：

• 当前解和当前能量的关系，对应current_solution和current_energy；
• 最优解和最优能量的关系，对应best_solution和best_energy；
• 当前解和最优能量的关系，对应current_solution和best_energy；
• 目标函数。

代码运行后，获得的图像为：

图3   模拟退火计算过程参数变化

在图3所示的模拟退火计算过程参数变化图像中，可以清晰看到四个图像：

目标函数图像是不会变化的线条；

当前解和当前能量动画为蓝色散点图，蓝色散点会反复出现，甚至往横坐标减小的方向变化，这就是退火过程中不断搜索的证据；

最优解和最优能量动画为绿色散点图，绿色散点图有时候不会变化，这是因为当前解在搜索过程中，不能发现更好的最优能量，所以最优解和最优能量不会更新；

当前解和最优能量动画为红色线条，红色线条会沿着坐标轴减小到方向往水平方向延展，这就是因为当前解在搜索过程中，不能发现更好的最优能量，所以最优能量不会更新，只会保持之前取到的值。

图3清晰展现了模拟退火算法运行过程中数据的变化。

此时的完整代码为：

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

# 设置 matplotlib 支持中文
plt.rcParams['font.family'] = 'SimHei'
# 解决负号显示问题
plt.rcParams['axes.unicode_minus'] = False

# 定义目标函数，这里是 f(x) = x^2
def objective_function(x):
    return x ** 2

# 初始化参数
# 初始解
initial_solution = np.random.uniform(-10, 10)
# 初始温度
initial_temperature = 100
# 终止温度
final_temperature = 0.01
# 温度衰减系数
alpha = 0.95
# 每个温度下的迭代次数
num_iterations_per_temp = 5

# 模拟退火算法
def simulated_annealing(initial_solution, initial_temperature, final_temperature, alpha, num_iterations_per_temp):
    # 当前解
    current_solution = initial_solution
    # 当前解的目标函数值
    current_energy = objective_function(current_solution)
    # 最优解
    best_solution = current_solution
    # 最优解的目标函数值
    best_energy = current_energy
    # 当前温度
    temperature = initial_temperature

    # 用于记录每一步的最优目标函数值
    energy_history = [best_energy]
    current_solution_history = [current_solution]
    current_energy_history = [current_energy]
    best_energy_history = [best_energy]
    best_solution_history =[best_solution]

    while temperature > final_temperature:
        for _ in range(num_iterations_per_temp):
            # 在当前解的邻域内生成新解
            new_solution = current_solution + np.random.uniform(-0.1, 0.1)
            # 计算新解的目标函数值
            new_energy = objective_function(new_solution)
            # 计算能量差
            delta_energy = new_energy - current_energy

            # 如果新解更优，直接接受
            if delta_energy < 0:
                current_solution = new_solution
                current_energy = new_energy
                # 更新最优解
                if new_energy < best_energy:
                    best_solution = new_solution
                    best_energy = new_energy
            # 如果新解更差，以一定概率接受
            else:
                acceptance_probability = np.exp(-delta_energy / temperature)
                if np.random.rand() < acceptance_probability:
                    current_solution = new_solution
                    current_energy = new_energy

            # 记录当前最优目标函数值
            energy_history.append(best_energy)
            # 更新当前解和能量的历史记录
            current_solution_history.append(current_solution)
            current_energy_history.append(current_energy)
            best_energy_history.append(best_energy)
            best_solution_history.append(best_solution)

        # 降低温度
        temperature *= alpha

    return best_solution, best_energy, energy_history, current_solution_history, current_energy_history, best_energy_history, best_solution_history

# 运行模拟退火算法
best_solution, best_energy, energy_history, current_solution_history, current_energy_history, best_energy_history, best_solution_history= simulated_annealing(
    initial_solution, initial_temperature, final_temperature, alpha, num_iterations_per_temp)

print("最优解:", best_solution)
print("最优值:", best_energy)
print('energy_history=', energy_history)

fig, ax = plt.subplots()  # 定义要画图
# 设置坐标轴范围
#ax.set_xlim(min(current_solution_history), max(current_solution_history))
#ax.set_ylim(min(current_energy_history), max(current_energy_history))
ax.set_xlim(-10,10)
ax.set_ylim(0,100)
# 定义散点图对象
scatter1 = ax.scatter([], [], c='b', label='当前解与当前能量关系')
scatter2 = ax.scatter([], [], c='g', marker="^",label='最优解与最优能量关系')
# 定义折线图对象来表示 best_energy 的变化
line, = ax.plot([], [], 'r-', label='当前解与最优能量关系')

num_frames = len(energy_history)
x = np.linspace(-10, 10, 400)
y = objective_function(x)
plt.plot(x, y, label='目标函数')
def update(frame):
    x = current_solution_history[:frame + 1]
    y = current_energy_history[:frame + 1]
    # 更新散点图的数据
    scatter1.set_offsets(np.c_[x, y])

    # 更新折线图的数据，这里使用迭代步数作为 x 轴，best_energy 作为 y 轴
    line_x = current_solution_history[:frame + 1]
    line_y = best_energy_history[:frame + 1]
    line.set_data(line_x, line_y)

    scatter_x=best_solution_history[:frame + 1]
    scatter_y =best_energy_history[:frame + 1]
    scatter2.set_offsets(np.c_[scatter_x, scatter_y])

    return scatter1, line,scatter2

ani = animation.FuncAnimation(fig, update, repeat=True,
                              frames=num_frames, interval=50)  # 调用animation.FuncAnimation()函数画图
plt.legend()

ani.save('mnthsf.gif') #保存动画
plt.show()

【4】细节说明

上述代码中，为加快计算时间，所以减少了计算次数，每个温度下只重复计算5次。如果时间充足，可以扩大这个计算次数。

【5】总结

使用python进一步研究了模拟退火算法运行过程中数据的变化。