启发式寻解
原理讲解
启发寻解:模拟退火
作者 : 老饼 日期 : 2022-09-14 19:00:46 更新 : 2022-09-14 19:15:51




   一.问题   

已知:
y=f(x)
求一x 使得f(x)值最小.
在无法求得精确解的情况下,则可用模拟退火算法对其寻找一个优秀解(不一定是最优的解,但至少局部最优)。
备注:这里的f(x) 不要求有表达式,可以是一个规则,更不要求可偏导



   二.算法背景   

借鉴物理退温过程设计的算法:

物理退温原理:


一高温物体,如果一下子把物体温度下降,这样的物体是比较脆的,也就是内部各粒子并未找到一个能量较小的位置。

但如果让它缓慢降温,粒子会在温度缓慢下降过程中,不断寻找能量更低的位置,最终,每个粒子都能找到一个能量较低的位置,从而使物体更加紧固(物体能量低,要改变它的结构,就需要更大能量,也即更紧固)。

为什么紧急降温不行,要缓慢降温呢,因为紧急降温每个粒子被逼马上收缩到一定位置。

而缓慢降温过程中,每个粒子由于温度还很高,它可以跃迁到不同位置(哪怕是能量更高的位置),而随着温度下降,它能跃迁的概率越来越小,随着缓慢降温,粒子刚开始的时候比较活跃,后慢慢的更倾向于收缩。

整个过程它有足够多的机会去找到更好的位置。


迁移到数学寻找极小值问题:


借鉴这一思想,迁移到我们的问题来,就是
x 对应 粒子
目标函数 对应 能量
我们的里的每个
就相当于物体的粒子,我们刚开始的时候,设置一个温度,在温度高的时候,x的调整比较自由,哪怕x的调整会让目标函数值更高,我们也用一定的概率接受它。
随着温度慢慢下降,接受更差的解的概率也降低。慢慢的,让每个x找到更好的位置。使整体目标函数值更小。




   三.算法流程   

1.初始化x,初始化温度。
2.外循环:
(2.1) 内循环(直接循环N次):
---------(2.11) 从x的邻域随机选择一个x作为新解x_{new},计算新解与旧解函数差 
---------(2.12) 按概率接受新解:如果  则接受新解,否则放弃新解,进入下一轮。
---------(2.13) 每次接受新解时,更新本轮内循环找到的最优最差解,更新历史找到的最优解。
(2.2) 降温(temp=temp*tempRate)
(2.3) 本轮找到的最好的解与最差的解差异不大,则跳出外循环
3.输出历史找到的最优解


   四.例子与代码   


  问题  


求 为何值时,
 取得最小值。
( 易知道,y的最小值在处取得,为0.)


下面展示用模拟退火法寻找解的具体过程,看结果是否与我们预期一致(python)。

# -*- coding: utf-8 -*-
"""
模拟退火算法
"""
import numpy as np

# ----定义能量计算函数(目标函数)----- 
def calE(x):
    return (x[0]-2)**2+(x[1]-3)**2

# -------初始化x---------- 
x=np.array([0,0])
E = calE(x)
teamp = 100

# ---------温度初始化可用此方法(取N次邻域能量差均值的10倍)------------- 
avgE= 0
try_cn = 20
for i in range(try_cn):
    new_x = x+ 0.1*(np.random.rand( 2)-0.5)
    new_E = calE(new_x)
    avgE = avgE+ abs(E-new_E)/try_cn
teamp =   avgE*10 
 
# -------------外循环------------------------- 
for t in range(10000):
    minE = E
    maxE = E
    
    # --------内循环----------------------- 
    for i in range(100):
        new_x = x+ 0.05*(np.random.rand( 2)-0.5)
        new_E = calE(new_x)
        
        # --------按概率接受新解------------------- 
        if(np.exp((E-new_E)/teamp)>np.random.rand()) :
            x = new_x
            E = new_E
            minE = min(E,minE)
            maxE = max(E,maxE)
            
    print("第"+str(t)+"轮:x="+str(x)+",E="+str(E)+",teamp:"+str(teamp))
    
    # -----------降温,并检查是否退出循环---------------------     
    teamp = teamp*0.5
    if( (maxE-minE)/maxE<0.01 ):
        break


结果:

第0轮:x=[-0.0369776  -0.16258671],E=14.151232416363044,teamp:2.1684776214061965
第1轮:x=[ 0.0299847  -0.07337761],E=13.326610236039077,teamp:1.0842388107030982
第2轮:x=[-0.05819167  0.26681517],E=11.706452281699834,teamp:0.5421194053515491
第3轮:x=[-0.00817767  0.4797082 ],E=10.384648303936515,teamp:0.27105970267577456
第4轮:x=[0.20577593 0.51937999],E=9.372715643129023,teamp:0.13552985133788728
第5轮:x=[0.50645925 0.86501793],E=6.788812407650605,teamp:0.06776492566894364
第6轮:x=[0.69372452 1.18815036],E=4.9891547382870245,teamp:0.03388246283447182
第7轮:x=[1.06131701 1.6479235 ],E=2.709236598898453,teamp:0.01694123141723591
第8轮:x=[1.47507248 2.17005339],E=0.9643602748259363,teamp:0.008470615708617955
第9轮:x=[1.88553456 2.6121867 ],E=0.16350149361940025,teamp:0.0042353078543089775
第10轮:x=[1.9784743  3.01239016],E=0.0006168717398273742,teamp:0.0021176539271544888
第11轮:x=[2.00887562 3.00266355],E=8.587110504778482e-05,teamp:0.0010588269635772444
第12轮:x=[2.0015052  2.98460822],E=0.00023917250869041312,teamp:0.0005294134817886222
第13轮:x=[2.01984423 2.98482052],E=0.0006242100682776004,teamp:0.0002647067408943111
第14轮:x=[2.00173833 2.98630884],E=0.00019046974378183272,teamp:0.00013235337044715555
第15轮:x=[1.99446963 3.00111664],E=3.183191654546613e-05,teamp:6.617668522357777e-05
第16轮:x=[1.99924394 3.00060683],E=9.39864724084946e-07,teamp:3.308834261178889e-05
第17轮:x=[1.99924394 3.00060683],E=9.39864724084946e-07,teamp:1.6544171305894443e-05
第18轮:x=[2.00036058 2.99942035],E=4.660124388700488e-07,teamp:8.272085652947222e-06

可见,最终的解为x = [2.00036058,2.99942035],与我们预期差异极小。




  五.对解无要求(求解TSP问题)


模拟退火并不受 f(x)和  x的表达形式所限制,看以下例子


(一) 问题:旅行商问题:


有N个城市 ,要求穿过每个城市,最终回到出发城市,问:寻找一条路径,使旅行总距离最小。


(二) 分析 


我们的x可以表示成以下形式: x = [0 1 5 3 2 4], 
它代表先到城市0,再到城市1,再到5.....最后由4回到0.
而总距离就是按这顺序行走的总距离。
这样的f(x) 和 x 并不具有非常严格的数值表达式。
而它们对模拟退火完全不受影响。只要定义好目标函数(能量值)的计算方法,和x的邻域解取值方法即可。


(三) 代码


下面是求解代码(python)

# -*- coding: utf-8 -*-
"""
模拟退火求解TSP旅行商问题
"""
import numpy as np
import matplotlib.pyplot as plt
import time
global cityLoc 
cityLoc = np.array([[3.64,2.68 ]  #beijing
,[4.18,1.76 ]  #shanghai
,[3.71,2.60 ]  #tianjin
,[1.33,3.30 ]  #wulumuqi
,[4.20,2.96 ]  #shenyang
,[3.92,1.82 ]  #nanjing
,[2.55,1.64 ]  #chengdu
,[2.37,1.02 ]  #kunming
,[3.43,2.09 ]  #zhengzhou
,[3.54,0.70 ]  #xianggang
,[3.51,1.62 ]  #wuhan
,[3.44,0.80 ]  #guangzhou
,[3.24,2.77 ]  #huhehaote
,[2.38,2.32 ]  #xiling
,[2.56,2.24 ]  #lanzhou
,[3.01,2.03 ]  #xian
,[2.79,2.51 ]  #yinchuan
,[4.03,1.16 ]  #fuzhou
,[3.47,0.70 ]  #aomen
,[1.30,1.69]]); #lasha
cityN= cityLoc.shape[0]

#------------定义能量计算方式(即目标函数值计算方法)------------------------------------------
def calE(x):
   loc = cityLoc[x] 
   E = np.sqrt(((loc[1:]-loc[0:-1])*(loc[1:]-loc[0:-1])).sum(1)).sum()
   E =E+np.sqrt(((loc[0]-loc[-1])*(loc[0]-loc[-1])).sum()) 
   return E

#------------定义取领域解的方法----------------------------
def getNeibX(x):
    new_x = x.copy()
    swap_idx = np.random.choice(np.arange(cityN), size=2, replace=False)
    new_x[swap_idx]=new_x[np.flip(swap_idx)]
    return new_x
       
x = np.random.choice(np.arange(cityN), size=cityN, replace=False)
E=  calE(x)

#---------温度初始化可用此方法(取N次邻域能量差均值的10倍)-------------
avgE= 0
try_cn = 20
for i in range(try_cn):
    new_x =  getNeibX(x)
    new_E = calE(new_x)
    avgE = avgE+ abs(E-new_E)/try_cn
teamp =   avgE*10 

#-------------外循环-------------------------
for t in range(1000):
    minE= E
    maxE= E
    
    #--------内循环-----------------------
    
    for i in range(3000):
        new_x = getNeibX(x)
        new_E = calE(new_x)
        
        #--------按概率接受新解-------------------
        
        if(np.exp((E-new_E)/teamp)>np.random.rand()):
            x = new_x
            E = new_E
            minE=min(E,minE)
            maxE=max(E,maxE)
    print("第"+str(t)+"轮:E="+str(E)+",teamp:"+str(teamp))
    print("minE="+str(minE)+",maxE:"+str(maxE))
    # print("x="+str(x))
    
    #--------降温,并检查是否退出循环--------   
    
    teamp= teamp*0.9
    
    if((maxE-minE)/maxE<0.001):
        break
    plt.clf()
    plt.plot(cityLoc[x,0],cityLoc[x,1],  color='r', marker='o', markerfacecolor='r', markersize=5)
    plt.scatter(cityLoc[:,0],cityLoc[:,1], marker='o', color='none',edgecolors='b',s=60)
    plt.draw()#注意此函数需要调用
    plt.pause(0.01)


运行结果:


结果分析:


20个城市,路径的组合有20!多种 ,是非常庞大的
但经模拟退火求解,却能自动找到一条不错的路径。
这展示了模拟退火的简单与强大!



联系老饼