模拟退火法

模拟退火算法

简介

模拟退火算法是一种通用概率演算法，用于在一个大的空间内寻找命题的最优解。

原理

根据热力学原理，在温度为T时，出现能量差为dE的降温概率为P(dE),标识为
P(dE) = exp(dE/(kT))
k为常数，exp为自然数，且dE<0，公式指出：

1. 温度越高，出现一次能量差为dE的降温概率越大。
2. 温度越低，出现降温的概率最小
   由于dE<0 --> dE/kT<0,exp(dE/(kT))的取值为(0,1)，随着温度T的降低，P(dE)会逐渐降低。我们将一次向较差解的移动看做一次温度跳变过程，我们以概率P(dE)来接受这样的移动。
   演变：
   内能E模拟为目标函数值 f，温度T演化成控制参数 t，即得到解组合优化问题的模拟退火演算法
   总结
   若f( Y(i+1) ) <= f( Y(i) ) (即移动后得到更优解)，则总是接受该移动。
   若f( Y(i+1) ) > f( Y(i) ) (即移动后的解比当前解要差)，则以一定的概率接受移动，而且这个概率随着时间推移逐渐降低（逐渐降低才能趋向稳定

伪代码

while(T>T_min):
    de = f(Y(i+1))-f(Y(i))  # 
    if de>=0:
        Y(i+1)=Y(i)  # 接受状态转移
   else:
        if exp(de/T)>random(0,1):
            Y(i+1) = Y(i)   # 接受移动
   T=r*T   # 降温退火，r越大，降温越慢，r越小，降温越快
   i++

案例

求函数最小值

import matplotlib.pyplot as plt
import mathi
mport random
T_init = 100  # 初始最大温度
alpha = 0.95  # 降温系数
T_min = 1e-3  # 最小温度，即退出循环条件
def obj(x):   
    y = 10 * math.sin(5 * x) + 7 * math.cos(4 * x)    
    return -y
def SA(T_init, alpha, T_min):    
    T = T_init    
    x_new = random.random() * 10  # 初解    
    x_current = x_new    
    y_current = float('inf')    
    x_best = x_new    
    y_best = float('inf')    
    while T > T_min:        
        for i in range(100):            
        delta_x = random.random() - 0.5            # 自变量变化后仍要求在[0,10]之间            
        if 0 < (x_new + delta_x) < 10:                
            x_new = x_new + delta_x            
       else:                
            x_new = x_new - delta_x            
       y_new = obj(x_new)            
       
       if (y_new < y_current):                
            y_current = y_new                
            x_current = x_new                
            if (y_new < y_best):                    
                y_best = y_new                    
                x_best = x_new            
           else:                
                if random.random() < math.exp(-(y_new - y_current) / T): 
                    y_current = y_new                    
                    x_current = x_new                
               else:                    
                    x_new = x_current        
        T *= alpha    
   print('最优解', x_best, obj(x_best))SA(T_init, alpha, T_min)