模拟退火算法

模拟退火算法示例：最小化二元函数

本示例使用模拟退火算法寻找二元函数的全局最小值。我们选择一个经典的测试函数——Rosenbrock 函数，也被称为“香蕉函数”，其数学表达式为：

$$ f(x, y) = (a - x)^2 + b(y - x^2)^2 $$

其中 $ a = 1 $, $ b = 100 $。该函数在点 $ (1, 1) $ 处取得全局最小值 $ f(1,1) = 0 $。

这个函数具有以下特点：

• 全局最小值位于一个狭窄、平缓的抛物线形山谷中。
• 局部优化方法（如梯度下降）容易陷入局部最优或收敛缓慢。
• 非常适合作为全局优化算法的测试用例。

Python 实现代码

import random
import math
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np

def simulated_annealing_2d():
    """
    使用模拟退火算法求解二元函数 Rosenbrock 的最小值
    f(x, y) = (1 - x)^2 + 100*(y - x^2)^2
    全局最小值在 (1, 1)，f(1,1) = 0
    """

    # ======================
    # 1. 参数配置
    # ======================
    T0 = 1000.0        # 初始温度（高温利于探索）
    alpha = 0.995      # 降温因子（控制冷却速度，通常0.8~0.99）
    T_min = 0.1        # 最小温度（停止条件）
    max_iter = 10000   # 最大迭代次数（防止无限循环）
    search_range = (-5, 5)  # 搜索空间范围：x 和 y 都在 [-5, 5] 内

    # ======================
    # 2. 定义目标函数：Rosenbrock 函数
    # ======================
    def objective(x, y):
        """
        计算 Rosenbrock 函数的值
        这是一个经典的非凸测试函数，全局最小值在 (1,1)
        """
        a = 1
        b = 100
        return (a - x)**2 + b * (y - x**2)**2

    # ======================
    # 3. 初始化：随机生成初始解
    # ======================
    # 在搜索范围内随机初始化 (x, y)
    x = random.uniform(search_range[0], search_range[1])
    y = random.uniform(search_range[0], search_range[1])

    current_value = objective(x, y)  # 当前解的目标函数值
    best_x, best_y = x, y            # 全局最优解
    best_value = current_value       # 全局最优值

    # 记录迭代过程（用于可视化和分析）
    temperature_history = []
    best_value_history = []
    x_history = []
    y_history = []

    # ======================
    # 4. 主循环：模拟退火过程
    # ======================
    T = T0           # 当前温度
    iteration = 0    # 迭代计数器

    while T > T_min and iteration < max_iter:
        # 4.1 生成新解（对当前解进行扰动）
        # 使用高斯噪声扰动 x 和 y 坐标
        x_new = x + random.gauss(0, 1)  # 均值为0，标准差为1的高斯噪声
        y_new = y + random.gauss(0, 1)

        # 确保新解在搜索空间内（边界处理）
        x_new = max(min(x_new, search_range[1]), search_range[0])
        y_new = max(min(y_new, search_range[1]), search_range[0])

        # 4.2 计算新解的目标函数值
        new_value = objective(x_new, y_new)

        # 4.3 应用 Metropolis 准则决定是否接受新解
        if new_value < current_value:
            # 新解更优：直接接受
            x, y = x_new, y_new
            current_value = new_value
        else:
            # 新解较差：以概率 exp(-(ΔE)/T) 接受
            delta_E = new_value - current_value  # 能量差 ΔE > 0
            acceptance_prob = math.exp(-delta_E / T)

            if random.random() < acceptance_prob:
                # 以一定概率接受劣质解（避免陷入局部最优）
                x, y = x_new, y_new
                current_value = new_value

        # 4.4 更新全局最优解
        if current_value < best_value:
            best_x, best_y = x, y
            best_value = current_value

        # 4.5 降温（冷却策略）
        T *= alpha  # 温度按比例衰减
        iteration += 1

        # 4.6 记录状态（用于后续可视化）
        temperature_history.append(T)
        best_value_history.append(best_value)
        x_history.append(x)
        y_history.append(y)

    # ======================
    # 5. 输出结果
    # ======================
    print(f"=== 模拟退火算法结果 ===")
    print(f"最优解: x = {best_x:.6f}, y = {best_y:.6f}")
    print(f"理论最优解: (1.000000, 1.000000)")
    print(f"目标函数值: f(x,y) = {best_value:.6f}")
    print(f"与理论最小值误差: {abs(best_value):.2e}")
    print(f"迭代次数: {iteration}")
    print(f"最终温度: {T:.4f}")

    # ======================
    # 6. 可视化结果
    # ======================
    fig = plt.figure(figsize=(18, 6))

    # 子图1：温度变化曲线
    plt.subplot(1, 3, 1)
    plt.plot(temperature_history, 'b-', linewidth=2)
    plt.title('Temperature Decay')
    plt.xlabel('Iteration')
    plt.ylabel('Temperature')
    plt.grid(True, alpha=0.3)

    # 子图2：最优值收敛过程
    plt.subplot(1, 3, 2)
    plt.plot(best_value_history, 'g-', linewidth=2, label='Best Value')
    plt.axhline(y=0, color='r', linestyle='--', label='Global Minimum')
    plt.title('Best Value Convergence')
    plt.xlabel('Iteration')
    plt.ylabel('f(x,y)')
    plt.legend()
    plt.grid(True, alpha=0.3)

    # 子图3：解的轨迹（二维投影）
    plt.subplot(1, 3, 3)
    plt.plot(x_history, y_history, 'b-', alpha=0.7, linewidth=1, label='Solution Path')
    plt.scatter(x_history[0], y_history[0], color='green', s=100, marker='o', label='Start')
    plt.scatter(x_history[-1], y_history[-1], color='red', s=100, marker='x', label='End')
    plt.scatter(1, 1, color='black', s=150, marker='*', label='Global Optimum (1,1)')
    plt.xlim(search_range)
    plt.ylim(search_range)
    plt.xlabel('x')
    plt.ylabel('y')
    plt.title('Solution Trajectory in 2D Space')
    plt.legend()
    plt.grid(True, alpha=0.3)

    plt.tight_layout()
    plt.show()

    # 可选：绘制三维曲面图（显示函数形状）
    plot_3d_surface(objective, search_range)

    return (best_x, best_y), best_value

def plot_3d_surface(objective_func, search_range):
    """
    绘制目标函数的三维曲面图
    """
    fig = plt.figure(figsize=(10, 8))
    ax = fig.add_subplot(111, projection='3d')

    # 生成网格点
    x = np.linspace(search_range[0], search_range[1], 100)
    y = np.linspace(search_range[0], search_range[1], 100)
    X, Y = np.meshgrid(x, y)
    Z = np.zeros_like(X)

    # 计算每个点的函数值
    for i in range(X.shape[0]):
        for j in range(X.shape[1]):
            Z[i, j] = objective_func(X[i, j], Y[i, j])

    # 绘制曲面
    surf = ax.plot_surface(X, Y, Z, cmap='viridis', alpha=0.8, edgecolor='none')

    # 添加颜色条
    fig.colorbar(surf, ax=ax, shrink=0.5, aspect=5)

    # 设置标签和标题
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('f(x,y)')
    ax.set_title('Rosenbrock Function 3D Surface')

    plt.show()

# ======================
# 7. 运行程序
# ======================
if __name__ == "__main__":
    result = simulated_annealing_2d()

代码详细注释说明

1. 参数配置

• T0: 初始温度要足够高，确保早期能充分探索解空间。
• alpha: 降温速度不宜过快（一般0.9~0.99），否则可能错过全局最优。
• search_range: 定义了变量的取值范围，避免解超出合理区间。

2. 目标函数设计

• Rosenbrock 函数是典型的非凸函数，其等高线呈弯曲的“香蕉”形状。
• 全局最小值在 (1,1)，但存在许多局部极小值，对优化算法挑战性大。

3. 解的生成与更新

• 使用高斯噪声 random.gauss(0, 1) 进行扰动，模拟原子热运动。
• 边界检查确保新解仍在有效范围内。

4. Metropolis 准则

• 是算法的核心：允许以一定概率接受劣质解，从而跳出局部最优陷阱。
• 概率公式：P = exp(-ΔE/T)，其中 ΔE 是能量差，T 是当前温度。

5. 冷却策略

• 采用几何降温：T = T * alpha
• 也可以尝试其他策略如指数降温、线性降温等。

6. 可视化功能

• 温度衰减曲线验证冷却过程是否符合预期。
• 最优值收敛图显示算法是否趋于稳定。
• 解的轨迹图直观展示算法在二维空间中的搜索路径。
• 三维曲面图帮助理解函数的复杂地形。

预期输出示例

=== 模拟退火算法结果 ===
最优解: x = 0.999876, y = 0.999752
理论最优解: (1.000000, 1.000000)
目标函数值: f(x,y) = 0.000001
与理论最小值误差: 1.23e-06
迭代次数: 7289
最终温度: 0.1023

总结

通过这个二元函数示例，我们可以看到：

• 模拟退火算法能够有效处理复杂的非凸优化问题。
• 即使从随机初始点出发，也能逐步逼近全局最优解。
• 可视化工具帮助我们理解算法的行为和性能。
• 该算法特别适用于存在多个局部最优的复杂优化场景。

此实现可以轻松扩展到更高维的问题，只需调整变量数量和扰动方式即可。