模拟退火算法及Python实现

一、模拟退火算法

        模拟退火算法包含两个部分即Metropolis算法和退火过程,，分别对应内循环和外循环。外循环就是退火过程，将固体达到较高的温度（初始温度T（0）），然后按照降温系数alpha使温度按照一定的比例下降，当达到终止温度Tf时，冷却结束，即退火过程结束。
        Metropolis算法是内循环，即在每次温度下，迭代L次，寻找在该温度下能量的最小值（即最优解）。下图中所示即为在一次温度下，跌代L次，固体能量发生的变化。在该温度下，整个迭代过程中温度不发生变化，能量发生变化，当前一个状态x(n)的能量大于后一个状态x(n+1)的能量时，状态x(n)的解没有状态x(n+1)的解好，所以接受状态x(n+1)。但是如果下一状态的能量比前一个状态的能量高时，该不该接受下一状态呢？在这里设置一个接受概率P,即如果下一状态的能量比前一个状态的能量高，则接受下一状态的概率为P。

二、python实现

"""
模拟退火算法
问题描述，随机生成解，在解中得到函数y = (-10*x1) - (math.e ** ((-x2/10)-x3))
"""
import math
import numpy as np
from random import random
import matplotlib.pyplot as plt

def func(x1,x2,x3):      # 所要优化的函数
    res = (-10*x1) - (math.e ** ((-x2/10)-x3))
    return res

class SA:
    def __init__(self,func,iter=100,T0=100,Tf=0.01,alpha=0.9):
        self.func = func
        self.iter = iter                  # 内循环迭代次数，本代码中，迭代次数为100
        self.alpha = alpha                # 降温系数
        self.T0 = T0                      # 初始温度
        self.Tf = Tf                      # 温度终值
        self.T = T0                       # 当前温度
        self.x1 = np.random.randint(0,2,size=120)       # 随机120生成解的值
        self.x2 = np.random.randint(0,80,size=120)
        self.x3 = np.random.randint(0,120,size=120)
        self.most_best = []
        """
        x1取是否两种可能
        x2取0到80的整数
        x3取0到120的整数
        """
        self.history = {"f": [],"T":[]}

    def generate_new(self,x1,x2,x3):     # 扰动产生新解的过程
        while True:
            # pro = [(self.T)]
            x1_new = np.random.choice([0,1])
            x2_new = x2 + np.floor(self.T *(np.random.randint(0,5)-np.random.randint(0,5)))
            x3_new = x3 + np.floor(self.T *(np.random.randint(0,5)-np.random.randint(0,5)))
            if (x1_new in [0,1]) & (0 <= x2_new <= 80) & (0 <= x3_new <= 120):       #重复得到新解。知道产生的解满足约束条件
                break
        return x1_new,x2_new,x3_new
    
    def Metrospolis(self,f,f_new):   # Metrospolis法则
        if f_new <= f:
            return 1
        else:
            p = math.exp((f - f_new)/self.T)
            if random() < p:
                return 1
            else:
                return 0
            
    def best(self):       # 获取最优目标函数值
        f_list = []       # f_list数组保存每次迭代之后的值
        for i in range(self.iter):
            f = self.func(self.x1[i%2],self.x2[i%80],self.x3[i%120])
            f_list.append(f)
        f_best = min(f_list)
        idx = f_list.index(f_best)
        return f_best, idx       # f_best, idx分别表示在该温度下，迭代L次之后目标函数的最优解和最优解的下标
    
    def run(self):
        count = 0
        # 外循环迭代，当前温度小于终止温度的阈值
        while self.T > self.Tf:
            # 内循环迭代100次
            for i in range(self.iter):
                f = self.func(self.x1[i],self.x2[i],self.x3[i])       # f为迭代一次后的值
                x1_new, x2_new, x3_new = self.generate_new(self.x1[i], self.x2[i],self.x3[i])      #产生新解
                f_new = self.func(x1_new,x2_new,x3_new)   # 产生新解
                if self.Metrospolis(f,f_new):             # 判断是否接受新值
                    self.x1[i] = x1_new                 # 如果接受新值，则把新值存入x数组和y数组
                    self.x2[i] = x2_new
                    self.x3[i] = x3_new
            # 迭代L次记录在该温度下最优解
            ft,_ = self.best()
            self.history['f'].append(ft)
            self.history['T'].append(self.T)
            # 温度按照一定比例下降（冷却）
            self.T = self.T*self.alpha
            count += 1
        # 得到最优解
        f_best,idx = self.best()
        print(f"F={f_best}, x1={self.x1[idx]}, x2={self.x2[idx]},x3={self.x3[idx]}")
sa = SA(func)
sa.run()

plt.plot(sa.history['T'], sa.history['f'])
plt.title('SA')
plt.xlabel('T')
plt.ylabel('f')
plt.gca().invert_xaxis()
plt.show()

运行结果：

结果和预想中的一致。

鸣谢：模拟退火算法详细讲解（含实例python代码）_退火算法代码详解-CSDN博客