万字长文了解模拟退火算法原理及求解复杂约束问题（源码实现）

àì夳堔傛蜴生んèń 2022-10-08 05:47 270阅读 0赞

## 模拟退火算法原理 ##

退火这个词，其实是铁匠发明的。它的意思很简单，就是将铁匠炉烧热后，再把下边的火撤掉，让金属在炉子里边慢慢冷却。人们发现，这个缓慢的降温过程能消除金属内部的各种缺陷，使得其恢复能量最低的状态。后来，受到退火工艺的启发，研究人员把将退火的缓慢的降温思路推广开来，用于寻找复杂函数的全局最优值。举个简单的例子，氦原子在金属中空位附近运动时，其势能函数大概长这样：  
    举个简单的例子，氦原子在金属中空位附近运动时，其势能函数大概长这样：  
![20210620160202228.png_pic_center][]

把势能函数看成一个轨道，再把氦看成一个小球放到这个轨道上。那么，小球的高度就是它的重力势能。求这个函数的全局最小值，其实就是让小球滚到最深那个坑过程。  
![20210620160734766.png_pic_center][]

由于函数中存在非常多局域极小值（浅坑），常见的贪心算法（如梯度下降法，就是闭着眼睛沿坡往下滚）会陷入最近的局域极小值，很错过到全局最小值（深坑）。  
    那么，怎么才能让滚进深坑呢？  
    我们先晃一晃整个体系，让小球随机动起来。我们把晃动的剧烈程度称为温度T。在足够长的时间内。小球出现在i点的概率呈玻尔兹曼分布，与i点的高度(势能Ei)和温度T密切相关：

温度为0时，小球出现在全局最小值的概率为100%。  
    那么，把温度设为0不就行了？  
    当然不行。上面的概率成立的前提是提供无穷长的时间，使得小球能够遍历整个函数。  
    实际中，无穷长的时间当然是不现实的。小球从低概率区向高概率区转移，这个过程是需要消耗时间的。  
    高温下，转移速度快。但另一方面，高温下，小球落在各个坑里的概率差别比较小。  
    为了让小球能够尽快的转移到最低点，我们可以采取一个降温策略：

*  先设定一个较高的温度T，使得小球有足够的能量越过能垒，小球到处乱跑的概率比较大。
 *  随机移动小球的位置，检查其移动前后的能量差 ΔE。
 *  若 ΔE<0 ，小球能量降低，则接受此次移动。
 *  若 ΔE>0，小球能量升高，有一定概率接受此次移动，概率为：exp(-ΔE/(KBT)) 。这个判定方法称为Metropolis判据。
 *  逐渐降低温度T，使得小球滑向低处的概率逐渐增大，直到小球缓缓收敛到最低点，不再晃动。  
    ![20210621140708718.png_pic_center][]  
    ![2021062114083063.png_pic_center][]

通过这样一个逐渐降温的过程，小球最终有很大概率（并不是一定）能找到全局最低值。这个算法就是模拟退火算法。  
    模拟退火算法最大的优势在于其通用性：函数形式如何复杂我不管，我就一步步瞎跑，然后按Metropolis判据概率性的接受这一步就行。  
![20210621091335489.gif_pic_center][]

流程图见下：  
![2021062110581587.png_pic_center][]

--------------------

## 算例 ##

![20210621153412616.png_pic_center][]

算例1

计算函数f(x)=∑x^2(-20≤x≤20)的最小值，其中个体x的维数n= 10。 这是一个简单的平方和函数，只有一个极小点x=(0, 0, …0)， 理论最小值f(0，0, .，. 0)=0

MATLAB版求解

%%%%%%%%%%%%%%%%%%%%%%模拟退火算法解决函数极值%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%初始化%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    clear all;                      %清除所有变量
    close all;                      %清图
    clc;                            %清屏
    D=10;                           %变量维数 
    Xs=20;                          %上限                                
    Xx=-20;                         %下限
    %%%%%%%%%%%%%%%%%%%%%%%%%%%冷却表参数%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    L = 200;                        %马可夫链长度
    K = 0.998;                      %衰减参数
    S = 0.01;                       %步长因子
    T=100;                          %初始温度
    YZ = 1e-8;                      %容差
    P = 0;                          %Metropolis过程中总接受点
    %%%%%%%%%%%%%%%%%%%%%%%%%%随机选点 初值设定%%%%%%%%%%%%%%%%%%%%%%%%%
    PreX = rand(D,1)*(Xs-Xx)+Xx;
    PreBestX = PreX;
    PreX =  rand(D,1)*(Xs-Xx)+Xx;
    BestX = PreX;
    %%%%%%%%%%%每迭代一次退火一次(降温), 直到满足迭代条件为止%%%%%%%%%%%%
    deta=abs( func1( BestX)-func1(PreBestX));
    while (deta > YZ) && (T>0.001)
        T=K*T; 
        %%%%%%%%%%%%%%%%%%%%%在当前温度T下迭代次数%%%%%%%%%%%%%%%%%%%%%%
        for i=1:L  
            %%%%%%%%%%%%%%%%%在此点附近随机选下一点%%%%%%%%%%%%%%%%%%%%%
                NextX = PreX + S* (rand(D,1) *(Xs-Xx)+Xx);
                %%%%%%%%%%%%%%%%%边界条件处理%%%%%%%%%%%%%%%%%%%%%%%%%%
                for ii=1:D
                    if NextX(ii)>Xs | NextX(ii)<Xx
                        NextX(ii)=PreX(ii) + S* (rand *(Xs-Xx)+Xx);
                    end
                end            
            %%%%%%%%%%%%%%%%%%%%%%%是否全局最优解%%%%%%%%%%%%%%%%%%%%%%
            if (func1(BestX) > func1(NextX))
                %%%%%%%%%%%%%%%%%%保留上一个最优解%%%%%%%%%%%%%%%%%%%%%
                PreBestX = BestX;
                %%%%%%%%%%%%%%%%%%%此为新的最优解%%%%%%%%%%%%%%%%%%%%%%
                BestX=NextX;
            end
            %%%%%%%%%%%%%%%%%%%%%%%% Metropolis过程%%%%%%%%%%%%%%%%%%%
            if( func1(PreX) - func1(NextX) > 0 )
                %%%%%%%%%%%%%%%%%%%%%%%接受新解%%%%%%%%%%%%%%%%%%%%%%%%
                PreX=NextX;
                P=P+1;
            else
                changer = -1*(func1(NextX)-func1(PreX))/ T ;
                p1=exp(changer);
                %%%%%%%%%%%%%%%%%%%%%%%%接受较差的解%%%%%%%%%%%%%%%%%%%%
                if p1 > rand        
                    PreX=NextX;
                    P=P+1;         
                end
            end
        trace(P+1)=func1( BestX);    
        end
        deta=abs( func1( BestX)-func1 (PreBestX)); 
    end
    disp('最小值在点:');
    BestX
    disp( '最小值为:');
    func1(BestX)
    figure
    plot(trace(2:end))
    xlabel('迭代次数')
    ylabel('目标函数值')
    title('适应度进化曲线')
    
    %%%%%%%%%%%%%%%%%%%%%%%%%适应度函数%%%%%%%%%%%%%%%%%%%%%%%
    function result=func1(x)
    summ=sum(x.^2);
    result=summ;
    end

从图中可以看到迭代次数太多，是个体寻优，需要花费很多时间。相比其他群体寻优方法，模拟退火算法很差。  
![20210621091834545.png_pic_center][]  
![20210621091937982.png_pic_center][]

python版求解

#!/usr/bin/env python3
    # -*- coding: utf-8 -*-
    # @Author: yudengwu(余登武）
    # @Date : 2021/6/21
    #@email:1344732766@qq.com
    
    import numpy as np
    import matplotlib.pyplot as plt
    import matplotlib as mpl
    import matplotlib; matplotlib.use('TkAgg')
    mpl.rcParams['font.sans-serif'] = ['SimHei']  # 指定默认字体
    mpl.rcParams['axes.unicode_minus'] = False  # 解决保存图像是负号'-'显示为方块的问题
    
    
    #=================初始化参数============
    D=10                           #变量维数
    Xs=20                          #上限
    Xx=-20                         #下限
    #==冷却表参数==
    L = 200                       #马可夫链长度 #在温度为t情况下的迭代次数
    K = 0.998                     #衰减参数
    S = 0.01                      #步长因子
    T=100                         #初始温度
    YZ = 1e-7                     #容差
    P = 0                         #Metropolis过程中总接受点
    #====随机选点 初值设定====
    PreX = np.random.uniform(size=(D,1))*(Xs-Xx)+Xx
    PreBestX = PreX#t-1代的全局最优X
    PreX =  np.random.uniform(size=(D,1))*(Xs-Xx)+Xx
    BestX = PreX#t时刻的全局最优X
    
    #==============目标函数=============
    def func1(x):
        return np.sum([i**2 for i in x])
    
    
    #====每迭代一次退火一次(降温), 直到满足迭代条件为止===
    deta=np.abs(func1(BestX)-func1(PreBestX))#前后能量差
    
    trace=[]#记录
    while (deta > YZ) and (T>0.1):#如果能量差大于允许能量差 或者温度大于阈值
        T = K * T#降温
        print(T)
        #=在当前温度T下迭代次数=
        for i in range(L):#
            #=在此点附近随机选下一点==
            NextX = PreX + S * (np.random.uniform(size=(D, 1)) * (Xs - Xx) + Xx)
            #===边界条件处理
            for ii in range(D):#遍历每一个维度
                while NextX[ii]>Xs or NextX[ii]<Xx:
                    NextX[ii]=PreX[ii] + S* (np.random.random() *(Xs-Xx)+Xx)
    
            #===是否全局最优解 ===
            if (func1(BestX) > func1(NextX)):
                #保留上一个最优解
                PreBestX = BestX
                #此为新的最优解
                BestX = NextX
    
            #====Metropolis过程===
            if (func1(PreX) - func1(NextX) > 0):#后一个比前一个好
               #接受新解
               PreX = NextX
               P = P + 1
            else:
                changer = -1 * (func1(NextX) - func1(PreX)) / T
                p1 = np.exp(changer)
                #接受较差的解
                if p1 > np.random.random():
                    PreX = NextX
                    P = P + 1
            trace.append(func1(BestX))
        deta = np.abs(func1(BestX)- func1(PreBestX))#修改前后能量差
    
    
    
    print('最小值点\n',BestX)
    print('最小值\n',func1(BestX))
    plt.plot(trace,label='迭代曲线')
    plt.xlabel('迭代次数')
    plt.ylabel('目标函数值')
    plt.show()

![2021062110560460.png_pic_center][] ![20210621105515379.png_pic_center][]

--------------------

复杂约束求解：算例2![20210621141221508.png_pic_center][]

复杂约束求解：python版求解

在约束求解中，添加了惩罚项。

import numpy as np
    import matplotlib.pyplot as plt
    import matplotlib as mpl
    import matplotlib; matplotlib.use('TkAgg')
    mpl.rcParams['font.sans-serif'] = ['SimHei']  # 指定默认字体
    mpl.rcParams['axes.unicode_minus'] = False  # 解决保存图像是负号'-'显示为方块的问题
    
    
    #=================初始化参数============
    D=2                           #变量维数 ,对应x,y
    Xs=1                          #搜索区间上限
    Xx=-1                         #搜索区间下限
    #==冷却表参数==
    L = 300                       #马可夫链长度 #在温度为t情况下的迭代次数
    K = 0.98                     #衰减参数
    S = 0.01                      #步长因子
    T=100                         #初始温度
    YZ = 1e-8                     #容差(运行t-1,，t时刻的误差）
    P = 0                         #Metropolis过程中总接受点
    eloss=0.1                     #允许的惩罚项误差
    #====随机选点 初值设定====
    PreX = np.zeros(shape=(D,1))
    PreX[0]=np.random.uniform(1,2,1)#X范围【1,2】
    PreX[1]=np.random.uniform(-1,0,1)#y范围【-1,0】
    PreBestX = PreX#t-1代的全局最优X
    
    PreX = np.zeros(shape=(D,1))
    PreX[0]=np.random.uniform(1,2,1)#X范围【1,2】
    PreX[1]=np.random.uniform(-1,0,1)#y范围【-1,0】
    BestX = PreX#t时刻的全局最优X
    
    
    #==============目标函数=============
    def func1(X):
        A = 10
        pi = np.pi
        x = X[0]
        y = X[1]
        return 2 * A + x ** 2 - A * np.cos(2 * pi * x) + y ** 2 - A * np.cos(2 * pi * y)
    #==========惩罚项===========
    def calc_e(X):
        """计算个体的目惩罚项，X 的维度是 size * 2 """
        ee = 0
        """计算第一个约束的惩罚项"""
        e1 = X[0] + X[1] - 6
        ee += max(0, e1)
        """计算第二个约束的惩罚项"""
        e2 = 3 * X[0] - 2 * X[1] - 5
        ee += max(0, e2)
        return ee
    
    
    #====每迭代一次退火一次(降温), 直到满足迭代条件为止===
    deta=np.abs(func1(BestX)-func1(PreBestX))#前后能量差
    
    
    trace=[]#记录
    while (deta > YZ) and (T>0.1):#如果能量差大于允许能量差 或者温度大于阈值
        T = K * T#降温
        print(T)
        # =在当前温度T下迭代次数=
        for i in range(L):  #
            # =在此点附近随机选下一点==
            NextX = PreX + S * (np.random.uniform(size=(D, 1)) * (Xs - Xx) + Xx)
            # ===边界条件处理
            while NextX[0] > 2 or NextX[0] < 1:#x属于【1,2】
                NextX[0] =PreX[0] + S* (np.random.uniform() * (Xs - Xx) + Xx)
            while NextX[1] > 0 or NextX[0] < -1:  # y属于【-1,0】
                NextX[1] = PreX[1] + S* (np.random.uniform() * (Xs - Xx) + Xx)
    
            # ===是否全局最优解 ===
            if (func1(BestX) > func1(NextX)) and (calc_e(BestX)<eloss) and (calc_e(NextX)<eloss):#两者都没有违反约束，下一代适应度更好时选下一代
                # 保留上一个最优解
                PreBestX = BestX
                # 此为新的最优解
                BestX = NextX
            elif (calc_e(BestX)>eloss) and (calc_e(NextX)>eloss) and (func1(BestX) > func1(NextX)):#两者都违反约束，下一代适应度更好时（选下一代）
                # 保留上一个最优解
                PreBestX = BestX
                # 此为新的最优解
                BestX = NextX
            elif (calc_e(BestX)>eloss) and (calc_e(NextX)<eloss) :#上一代惩罚项不满足约束，下一代惩罚项满足约束（选下一代）
                # 保留上一个最优解
                PreBestX = BestX
                # 此为新的最优解
                BestX = NextX
    
    
            # ====Metropolis过程 当代温度T 前后个体判断===
            if (func1(PreX) - func1(NextX) > 0) and (calc_e(PreX)<eloss) and(calc_e(NextX)<eloss) :  # 两个惩罚项都满足约束，后一个的适应度更小（选后一个）
                # 接受新解
                PreX = NextX
                P = P + 1
            elif  (calc_e(PreX)>eloss) and(calc_e(NextX)>eloss) and (func1(PreX) - func1(NextX) > 0) :#两个惩罚项都不满足时，后一个适应度比前一个小（选后一个）
                # 接受新解
                PreX = NextX
                P = P + 1
            elif (calc_e(PreX)>eloss) and(calc_e(NextX)<eloss) :#前一个惩罚项不满足，后一个满足（选后一个）
                # 接受新解
                PreX = NextX
                P = P + 1
    
            else:
                changer = -1 * (func1(NextX) - func1(PreX)) / T
                p1 = np.exp(changer)
                # 接受较差的解
                if p1 > np.random.random():
                    PreX = NextX
                    P = P + 1
            trace.append(func1(BestX))
        deta = np.abs(func1(BestX) - func1(PreBestX))  # 修改前后能量差
    
    
    
    print('最优值X\n',BestX)
    print('最优目标函数值\n',func1(BestX))
    print('最优值对应的惩罚项',calc_e(BestX))
    plt.plot(trace,label='迭代曲线')
    plt.xlabel('迭代次数')
    plt.ylabel('目标函数值')
    plt.show()

![20210621153042307.png_pic_center][] ![20210621153131890.png_pic_center][]

![在这里插入图片描述][20210621153300483.jpg_pic_center]

作者：电气-余登武

[20210620160202228.png_pic_center]: /images/20221005/8b383f179b0d45719efbce7aa289ec12.png
[20210620160734766.png_pic_center]: /images/20221005/bf66c3980ff7426c98f72939f650deb0.png
[20210621140708718.png_pic_center]: /images/20221005/bcbe2745a21e4f6c9a34f15f02badbbf.png
[2021062114083063.png_pic_center]: /images/20221005/d70d4f6f2c984cc4bf975ee089e82856.png
[20210621091335489.gif_pic_center]: /images/20221005/841d470f9aff4782988edddb6e45fdf5.png
[2021062110581587.png_pic_center]: /images/20221005/aa6658de32324354b0cad5b1694200bc.png
[20210621153412616.png_pic_center]: /images/20221005/511038e47c5e444388ce15f8de56d37c.png
[20210621091834545.png_pic_center]: /images/20221005/b0592ce7c8ce4e50b754e92c1b4b016d.png
[20210621091937982.png_pic_center]: /images/20221005/509366203ce0400490673401fb79082c.png
[2021062110560460.png_pic_center]: https://img-blog.csdnimg.cn/2021062110560460.png#pic_center
[20210621105515379.png_pic_center]: /images/20221005/ea47abb03a864df4ada58c119ee912c9.png
[20210621141221508.png_pic_center]: /images/20221005/651c9bd1db8349f189b4883cc0d6eefd.png
[20210621153042307.png_pic_center]: https://img-blog.csdnimg.cn/20210621153042307.png#pic_center
[20210621153131890.png_pic_center]: /images/20221005/69a73584af374c0d96b3cbf933a5815e.png
[20210621153300483.jpg_pic_center]: /images/20221005/ed08dae4004f43deace2809fc8e2b7ae.png