模拟退火算法

模拟退火算法

我们知道，一个复杂的优化问题可能难以求出其最优解，或者是在短时间内无法求出。这时，有一个折中的办法就是利用启发式算法。其主要特点是在某些情况下可以提供一个快速的答案，但不一定是最优解。启发式算法通常基于经验和直觉，并利用特定的技巧或规则来找到问题的近似解。现阶段，启发式算法以仿自然体算法为主，主要有蚁群算法、模拟退火算法、遗传算法等。下面我们来介绍模拟退火算法。

模拟退火是一种全局优化算法，受到固体退火过程的启发。固体退火是一个物理过程，其中材料加热然后缓慢冷却，以达到低能量状态。

原理

1. 初始化阶段：选择一个初始解，并设置初始温度。
2. 迭代搜索：在当前解的邻域中随机选择一个新解。
3. 决策：基于某种准则接受或拒绝新解。
4. 退火过程：减少温度并返回第二步，直到满足某些停止准则。

在步骤3中，如果新解比当前解更好，就接受它。如果新解更差，仍然有可能接受它，但这取决于温度和解的差异程度。接受差解的概率通常由以下公式给出：

$$P(\text{接受})=\exp \left(-\frac{\text{差异}}{\text{温度}}\right)$$

案例求解

首先我们考虑一个简单的优化问题：在区间 $[-10,10]$ 中找到函数 $f(x)=x^2$的最小值。

Python 实现：

import math
import random

def objective_function(x):
    return x ** 2

def simulated_annealing():
    current_solution = random.uniform(-10, 10)
    current_value = objective_function(current_solution)

    T = 1000  # 初始温度
    T_min = 0.001
    alpha = 0.99  # 温度衰减率

    while T > T_min:
        # 随机选择一个新的解
        new_solution = current_solution + random.uniform(-1, 1)
        new_value = objective_function(new_solution)

        # 计算目标函数的差异
        delta = new_value - current_value

        if delta < 0 or random.random() < math.exp(-delta / T):
            current_solution = new_solution
            current_value = new_value

        T = T * alpha  # 温度衰减

    return current_solution, current_value

solution, value = simulated_annealing()
print(f"最优解: x = {solution}, f(x) = {value}")

所得结果如下：
$$\begin{gathered} 最优解:x=0.02350245227952774 \\ f(x)=0.0005523652631514788 \end{gathered}$$
可以看到，这个解接近于最优解。我们可以通过修改参数来提升精度，比如设置初始温度为10000度，得到如下结果：
$$最优解:x=-0.0009553046042742874,f(x)=9.126068869476529e-07$$
当然由于算法的随机性，每次运行得到的结果略有不同。

下面我们考虑一个更复杂的例子：
min ⁡    ∑ i = 1 50 k i s . t . { ∑ i = 1 50 ( 1 − β ) t i w i k i ≥ 2.82 × 1 0 4 ∑ i = 1 50 k i t i ≤ 8 × 6000 k i ∈ { 0 , 1 } \begin{equation*} \begin{split} &\min \,\, \sum\limits _{i=1}^{50} k_{i}\\ &s.t.\quad \left\{\begin{array}{lc} \sum\limits _ { i = 1 } ^ { 50 } \frac { ( 1 - { \beta } ) t _ { i } } { w _ { i } } k _ { i } \geq 2.82 \times 10 ^ { 4 } \\ \sum\limits _ { i = 1 } ^ { 50 } k _ { i } t _ { i } \leq 8 \times 6000 \\ k_{i}\in \left \{0,1\right \}\\ \end{array}\right. \end{split} \end{equation*} ​mini=1∑50​ki​s.t.⎩ ⎨ ⎧​i=1∑50​wi​(1−β)ti​​ki​≥2.82×104i=1∑50​ki​ti​≤8×6000ki​∈{0,1}​​​
其中，$t_i,w_i,\beta$均已知。

我们使用模拟退火算法求解：

import numpy as np
import random
import math
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation, PillowWriter


# 设置matplotlib支持中文
plt.rcParams['font.sans-serif'] = ['SimHei']  # 指定默认字体为黑体
plt.rcParams['axes.unicode_minus'] = False  # 解决保存图像时负号'-'显示为方块的问题


w = np.array([0.66, 0.6, 0.72, 0.66, 0.72, 0.6, 0.66, 0.6, 0.66, 0.66, 0.6, 0.6, 0.66, 0.66, 0.72, 0.72, 0.72, 0.6, 0.72, 0.6, 0.6, 0.6, 0.72, 0.6,
             0.72, 0.72, 0.72, 0.72, 0.72, 0.66, 0.66, 0.66, 0.66, 0.66, 0.66, 0.66, 0.6, 0.72, 0.72, 0.72, 0.72, 0.6, 0.72, 0.6, 0.6, 0.6, 0.72, 0.66, 0.72, 0.6])
beta = 0.01
t = np.array([1379.210046, 1478.695833, 1367, 1003.958333, 810.4083333, 2928.178571, 714.275, 705.5833333, 684.0630631, 570.825, 660.6375, 652.1583333, 569.3833333, 572.9666667, 542.9458333, 540.775, 525.4, 476.3969072, 422.3541667, 370.9625, 344.9458333, 457.2865497, 522.4175824, 1024.905405, 205.1,
             236.2416667, 194.1541667, 322.8407643, 173.4625, 171.6958333, 132.9375, 119.8458333, 147.5931373, 100.1708333, 109.7291667, 96.83333333, 28.95, 80.15416667, 78.50833333, 64.5125, 68.35833333, 27.1625, 87.13106796, 47.73362445, 7.525, 26.87083333, 7.106382979, 2.004854369, 1.908629442, 36.456621])


# 目标函数
def objective_function(k):
    return np.sum(k)

# 检查解是否满足约束条件


def constraints_satisfied(k):
    constraint1 = np.sum(((1 - beta) * t / w) * k) >= 2.82e4
    constraint2 = np.sum(k * t) <= 8 * 6000
    return constraint1 and constraint2


def simulated_annealing_animation():
    solutions = []
    values = []

    current_solution = np.random.choice([0, 1], 50)
    while not constraints_satisfied(current_solution):
        current_solution = np.random.choice([0, 1], 50)

    current_value = objective_function(current_solution)
    solutions.append(current_solution.copy())
    values.append(current_value)

    T = 500  # 初始温度
    T_min = 0.001
    alpha = 0.95  # 温度衰减率

    while T > T_min:
        new_solution = np.copy(current_solution)
        flip_index = random.randint(0, 49)
        new_solution[flip_index] = 1 - new_solution[flip_index]

        if constraints_satisfied(new_solution):
            new_value = objective_function(new_solution)
            delta = new_value - current_value

            if delta < 0 or random.random() < math.exp(-delta / T):
                current_solution = new_solution
                current_value = new_value

            solutions.append(current_solution.copy())
            values.append(current_value)

        T = T * alpha

    return solutions, values


solutions, values = simulated_annealing_animation()

# 设置动画
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8))

# 设定背景颜色和网格线
for ax in [ax1, ax2]:
    ax.set_facecolor("#f4f4f4")
    ax.grid(True, which='both', color='white', linewidth=0.5)
    ax.set_axisbelow(True)

# 辅助函数，用于根据解的状态确定颜色
def bar_colors(solution): 
    return ['#1f77b4' if k == 1 else '#ff7f0e' for k in solution]

bars = ax1.bar(range(len(solutions[0])), solutions[0], color=bar_colors(solutions[0]))
ax1.set_title("随时间变化的最佳解")
ax1.set_xlabel("索引")
ax1.set_ylabel("值 (0 或 1)")

line, = ax2.plot(values, color='#2ca02c', label="目标函数的值")
ax2.set_title("目标函数的值随时间变化")
ax2.set_xlabel("迭代次数")
ax2.set_ylabel("值")
ax2.legend()

# 更新动画的函数
def update(frame):
    for bar, k in zip(bars, solutions[frame]):
        bar.set_height(k)
        bar.set_color('#1f77b4' if k == 1 else '#ff7f0e')
    line.set_data(range(frame + 1), values[:frame + 1])
    ax2.set_xlim(0, frame + 1)
    ax2.set_ylim(min(values), max(values) + 1)
    return bars, line

ani = FuncAnimation(fig, update, frames=len(solutions), repeat=False)

plt.tight_layout()
plt.show()

# 保存动画为GIF
writer = PillowWriter(fps=20)
ani.save("动画.gif", writer=writer)

# 打印最后的结果
final_solution = solutions[-1]
final_value = values[-1]

print("最终的解：", final_solution)
print("目标函数的值：", final_value)

结果如下：
$$\begin{gathered} 最终的解:[11111111111111101011110100100001101000000000000000] \\ 目标函数的值:25 \end{gathered}$$

可见，模拟退火算法对复杂的规划问题是一个很好的解决方案。