模拟退火算法的Python实现（非线性函数优化为例）

1、模拟退火算法（SA）介绍

模拟退火算法(Simulated Annealing，SA)是一种随机寻优算法，其出发点是基于物理中固体物质的退火过程与一般组合优化问题之间的相似性。

模拟退火算法从某一较高初温出发，伴随温度参数的不断下降，结合概率突跳特性在解空间中随机寻找目标函数的全局最优解，即在局部最优解能概率性地跳出并最终趋于全局最优。

（一）算法流程：

（二）Metropolis准则：

Metropolis准则用于判断是否接受新解。模拟退火中函数最优化的过程是f值不断降低的过程，但为了跳出局部最优解，当函数值变大时也以一定概率接受，且这个概率随着温度T的降低而变小，最终变化趋于稳定，这也是模拟退火的核心思想。

（三）非线性函数优化中扰动产生新解的过程

在非线性函数优化中，如何产生新的一组解是重要的问题，由于不能像TSP（旅行商）问题那样用多元置换、互换等方式简单产生新解，所以在这里我们将扰动与模拟退火的思想结合，扰动产生新解的幅度由系统温度决定，得到公式如下：

x1 = x + T ∗ R（R∈[−1,0)∪(0,1]）

其中，T为当前温度，R为[-1, 1]间不为零的一个随机数，以此公式来产生新解，并且随着温度的下降，新解的扰动幅度也会减小，最终趋于稳定。

2、Python实现（非线性函数优化为例）

（1）定义待优化函数

其中，x∈[0,10]，y∈[0,10]

import math

def func(x, y):
    num = 6.452 * (x + 0.125 * y) * (math.cos(x) - math.cos(2 * y)) ** 2
    den = math.sqrt(0.8 + (x - 4.2) ** 2 + 2 * (y - 7) ** 2)
    return - (num / den + 3.226 * y)

• 由于要求该函数的最大值，而模拟退火中求得的是目标函数的最小值，所以需要在函数前加上负号

（2）初始化SimulateAnnealing类

class SimulateAnnealing:
    def __init__(self, func, iter=10, T0=100, Tf=1e-8, alpha=0.99):
        """
        Annealing parameters
        :param iter: Number of internal cycles
        :param T0: Initial temperature
        :param Tf: Final temperature
        :param alpha: Cooling factor
        """
        self.func = func
        self.iter = iter
        self.alpha = alpha
        self.T0 = T0
        self.Tf = Tf
        self.T = T0
        self.x = [random() * 10 for i in range(iter)]
        self.y = [random() * 10 for i in range(iter)]
        self.history = {'f': [], 'T': []}

• iter：内循环迭代次数，目的是基于不同初值进行多次退火，避免陷入局部最优
• T0：初始温度，可以根据解的约束条件等适当进行调整
• Tf：温度终值，一般可以取1e-8、1e-13等较小的值，目的是降低收敛速度，提高精确度
• alpha：降温系数，一般取0.99，目的同上

（3）扰动产生新解的过程

    def generate_new(self, x, y):
        while True:
            x_new = x + self.T * (random() - random())
            y_new = y + self.T * (random() - random())
            if (0 <= x_new <= 10) & (0 <= y_new <= 10):
                break
        return x_new, y_new

• random() - random()可以得到[-1, 1]间的随机数， 并且降低了取到0的概率
• 重复这一过程，直到产生的新解满足约束条件

（4）Metropolis准则

    def Metrospolis(self, f, f_new):
        if f_new <= f:
            return 1
        else:
            p = math.exp((f - f_new) / self.T)
            if random() < p:
                return 1
            else:
                return 0

• 输入f_old与f_new，按Metropolis准则判断是否接受新解

（5）获取最优目标函数值

    def get_optimal(self):
        f_list = []
        for i in range(self.iter):
            f = self.func(self.x[i], self.y[i])
            f_list.append(f)
        f_best = min(f_list)
        idx = f_list.index(f_best)
        return - f_best, idx

• 用来返回当前解中最优的函数值

（6）主函数

    def run(self):
        count = 0
        # annealing
        while self.T > self.Tf:
            # iteration
            for i in range(self.iter):
                f = self.func(self.x[i], self.y[i])
                x_new, y_new = self.generate_new(self.x[i], self.y[i])
                f_new = self.func(x_new, y_new)
                if self.Metrospolis(f, f_new):
                    self.x[i] = x_new
                    self.y[i] = y_new
            # save to history
            ft, _ = self.get_optimal()
            self.history['f'].append(ft)
            self.history['T'].append(self.T)
            # cooling
            self.T = self.T * self.alpha
            count += 1
        # get optimal solution
        f_best, idx = self.get_optimal()
        print(f"F={f_best}, x={self.x[idx]}, y={self.y[idx]}")

• 模拟退火过程，当温度大于温度终值时重复迭代并降温
• 内循环迭代时，重复产生新解并代入Metropolis准则判断是否接受，接受则覆盖原来的解
• 将函数值和对应的温度存入历史记录，用以绘图
• 最后输出最优解及对应的函数值

（7）运行结果

if __name__ == '__main__':
    # run
    sa = SimulateAnnealing(func)
    sa.run()
    # plot
    sa.plot()
    sa.plot([0, 1], [99, 100])

实例化并运行，得到结果如下（为了更加直观展示结果，以下函数值均为退火过程中目标函数值的相反数）：

F=99.99522530386253, x=6.090887224877577, y=7.798623593311333

• 整个退火过程：

可以看到在温度较高时解和目标函数的变化幅度比较大，但随着温度的降低状态逐渐稳定，波动幅度减小，并且温度越低变化的频率越高，在这张图上很难看到温度小于1°时具体的函数变化。

• 温度在0-1之间，函数值在99-100之间的退火过程：

放大后可以看到温度小于1°时，函数值随着温度的下降不断波动上升，并且比先前更加稳定，最终趋近于最优值。

以上是模拟退火算法的Python实现，最后附上完整代码：

# !/usr/bin/env python
# -*- coding: utf-8 -*-
# @Author: gyy
import math
from random import random
import matplotlib.pyplot as plt

def func(x, y):
    num = 6.452 * (x + 0.125 * y) * (math.cos(x) - math.cos(2 * y)) ** 2
    den = math.sqrt(0.8 + (x - 4.2) ** 2 + 2 * (y - 7) ** 2)
    return - (num / den + 3.226 * y)

class SimulateAnnealing:
    def __init__(self, func, iter=10, T0=100, Tf=1e-8, alpha=0.99):
        """
        Annealing parameters
        :param iter: Number of internal cycles
        :param T0: Initial temperature
        :param Tf: Final temperature
        :param alpha: Cooling factor
        """
        self.func = func
        self.iter = iter
        self.alpha = alpha
        self.T0 = T0
        self.Tf = Tf
        self.T = T0
        self.x = [random() * 10 for i in range(iter)]
        self.y = [random() * 10 for i in range(iter)]
        self.history = {'f': [], 'T': []}

    def generate_new(self, x, y):
        while True:
            x_new = x + self.T * (random() - random())
            y_new = y + self.T * (random() - random())
            if (0 <= x_new <= 10) & (0 <= y_new <= 10):
                break
        return x_new, y_new

    def Metrospolis(self, f, f_new):
        if f_new <= f:
            return 1
        else:
            p = math.exp((f - f_new) / self.T)
            if random() < p:
                return 1
            else:
                return 0

    def get_optimal(self):
        f_list = []
        for i in range(self.iter):
            f = self.func(self.x[i], self.y[i])
            f_list.append(f)
        f_best = min(f_list)
        idx = f_list.index(f_best)
        return - f_best, idx

    def plot(self, xlim=None, ylim=None):
        plt.plot(sa.history['T'], sa.history['f'])
        plt.title('Simulate Annealing')
        plt.xlabel('Temperature')
        plt.ylabel('f value')
        if xlim:
            plt.xlim(xlim[0], xlim[-1])
        if ylim:
            plt.ylim(ylim[0], ylim[-1])
        plt.gca().invert_xaxis()
        plt.show()

    def run(self):
        count = 0
        # annealing
        while self.T > self.Tf:
            # iteration
            for i in range(self.iter):
                f = self.func(self.x[i], self.y[i])
                x_new, y_new = self.generate_new(self.x[i], self.y[i])
                f_new = self.func(x_new, y_new)
                if self.Metrospolis(f, f_new):
                    self.x[i] = x_new
                    self.y[i] = y_new
            # save to history
            ft, _ = self.get_optimal()
            self.history['f'].append(ft)
            self.history['T'].append(self.T)
            # cooling
            self.T = self.T * self.alpha
            count += 1
        # get optimal solution
        f_best, idx = self.get_optimal()
        print(f"F={f_best}, x={self.x[idx]}, y={self.y[idx]}")

if __name__ == '__main__':
    # run
    sa = SimulateAnnealing(func)
    sa.run()
    # plot
    sa.plot()
    sa.plot([0, 1], [99, 100])
    # result:
    # F=99.99522530386253, x=6.090887224877577, y=7.798623593311333