在这里插入图片描述

作者:非妃是公主
专栏:《智能优化算法》
博客地址https://blog.csdn.net/myf_666
个性签:顺境不惰,逆境不馁,以心制境,万事可成。——曾国藩
在这里插入图片描述

专栏推荐

专栏名称专栏地址
软件工程专栏——软件工程
计算机图形学 专栏——计算机图形学
操作系统专栏——操作系统
软件测试专栏——软件测试
机器学习专栏——机器学习
数据库专栏——数据库
算法专栏——算法

对于一个问题,解决他的第一步是建立模型,利用数学的方法,将问题抽象成一个数学模型(可能是一个函数,一个算法,一个方程……)。第二步,就是对这个问题进行求解,也叫做模型求解。

但对于模型求解来讲,不同的模型有着不同大小的难度。传统的算法,比如梯度下降,很可能会陷入局部最优解等问题,对于更加全局性的最优解搜索,智能优化算法有着很好的效果。而模拟退火算法,就是智能优化算法中的一种,有着很好的效果。

模拟退火算法(Simulated Annealing,SA)最早由Metropolis等人于 1953 年提出。1983 年Kirkpatrick等人第一次使用模拟退火算法求解组合优化问题后,它就发表在了 Science 上。1直到今天,它依然被广泛使用,这篇文章将详细介绍模拟退火算法的基本原理,以及matlab的代码实现。


一、概论

模拟退火算法其实就是一个类似于仿生学的算法,模仿的就是物理退火的过程。我们炼钢的时候,如果我们急速冷凝,这时候的状态是不稳定的,原子间杂乱无章的排序,能量很高。而如果我们让钢水慢慢冷凝,很缓慢的降温,那么这个时候的状态就是很稳定的,各个分子都趋向于自己能量最低的位置。

而模拟退火算法,恰恰就是利用了物理退火这一过程的原理,求解一个优化目标(目标函数)的最小值。


二、物理退火

如果想说清楚模拟退火,必然绕不过物理退火!

1. 加温过程

增强粒子热运动,使其偏离平衡位置。当温度足够高时,固体将溶解为液体,从而消除系统原先可能存在的非均匀态,使随后进行的冷却过程以某一平衡态为起点。溶解过程于系统的能量增大过程相联系,系统能量也随温度的升高而增大。


2. 等温过程

通过物理学知识得知,对于与周围环境交换热量而温度不变的封闭系统,系统状态的自发变化总是朝着自由能降低的方向进行;当自由能达到最低时,系统达到平衡态。


3. 降温过程

随着温度降低,分子热运动减弱,趋向有序,系统能量逐渐降低,从而得到了低能量的晶体结构。


三、模拟退火原理

而模拟退火算法,就是要通过如上3个部分的操作,获得低能量的晶体结构(最优解)。

物理退火与模拟退火中的各个状态对应如下:

物理退火模拟退火
粒子状态每个状态对应一个(可行)解
能量最低态最优解
溶解过程设定初温(设定参数T的值)
等温过程一个温度下,多次采样(Metropolis采样过程)
冷却控制参数下降
能量目标函数

主要思想,在搜索区间进行随机游走(通过生成随机数实现),再利用Metropolis抽样准则,使随机游走逐渐收敛于局部最优解。在这里,温度是一个重要的控制参数,这个参数的大小控制了随机过程向局部或全局最优解移动的快慢。具体的过程如下。

首先,定义一个 p p p
p = e − E 2 − E 1 T p=e^{-\frac{E_2-E_1}{T}} p=eTE2E1

这个公式表示,系统从 E 1 E_1 E1变化到 E 2 E_2 E2,其概率为 p p p,就是上述公式左边的部分。

如果 E 2 < E 1 E_2<E_1 E2<E1,那么证明系统向更低的方向转移了,我们无条件接受此状态。

否则,以上述的概率,接受这个较坏的结果 E 2 E_2 E2。注意,这里的原则,不同于贪心,每次选择最好的,而是有一定的概率接受坏一些的。(这个概率受到能量差和温度两个因素的影响)

p ( 1 → 2 ) = { 1 E 2 < E 1 e − E 2 − E 1 T E 2 > E 1 p(1\rightarrow2)= \begin{cases} 1& E_2<E_1\\ e^{-\frac{E_2-E_1}{T}}& E_2>E_1 \end{cases} p(12)={1eTE2E1E2<E1E2>E1

这样通过一定的迭代次数,我们就会找到一个比较好的解了。


四、模拟退火的优点

总结一下模拟退火的优点如下:

  • 以一定的概率接受恶化解。看似在接受恶化解,其实是寻找了全局最优解。
  • 引进算法控制参数——温度。温度的引进,使得算法更加智能、灵活。在前期,进行跳跃搜索,更容易跳出局部最优解,提高搜索的全局性。在后期,进行幅度较小的搜索,更容易收敛。
  • 对目标函数要求少。模拟退火算法其实是一种搜索或者说枚举的方法,我们经过大量的实验,最终得到最优解(或者一定程度上的较优解)。由于这种做法,因此对目标函数的要求很低。不需要依赖什么,直接猜就好。

五、算法具体流程

学习的过程是广度优先,层层深入的,前面大致了解了模拟退火算法的原理及思想,下面我们看一看模拟退火算法的具体算法流程。


1. 整体流程

  1. 初始化温度 T 0 T_0 T0,初始解状态 X 0 X_0 X0(算法迭代起点)、每个 T T T值的迭代次数 L L L(也叫做Markov 链长度,就是在同一个温度下,我们猜测的次数)。
  2. k = 1 , … , L k=1, …, L k=1,,L做第(3)至第(6)步;
  3. 产生新解 X ′ X' X
  4. 计算增量 Δ E = E ( X ′ ) − E ( X ) ΔE=E(X')-E(X) ΔE=E(X)E(X),其中 E ( X ) E(X) E(X)为评价函数(越低越好);
  5. Δ E < 0 ΔE<0 ΔE0,则接受 X ′ X' X作为新的当前解,否则以概率 e − Δ E T e^{\frac{-ΔE}{T}} eTΔE接受 X ′ X' X作为新的当前解;
  6. 如果满足终止条件,则输出当前解为最优解,结束程序;
  7. T T T逐渐减小,且 T → 0 T\rightarrow 0 T0,然后转第 2 步。

流程图如下:

在这里插入图片描述


2. 细节处理

Ⅰ. 状态产生函数

主要就是通过在邻域内随机进行选择产生的。


Ⅱ. 初始温度

初始温度对于算法的影响较大,而且效果的好坏与求解问题的求解空间有关。初温越大,获得高质量解的可能性越高,但是相应的时间也就越长,具体可以在此处进行取舍。

一种很好的做法是,可以均匀抽样一组状态,以各状态目标值的方差作为初温。

  • 这样,如果方差越大,证明该问题越有可能陷入到局部最优解(多峰型优化),因此初温较高,防止陷入。

  • 如果方差很小,说明该问题不容易陷入局部最优(单峰型优化),因此,初温较低,可以快速收敛。


Ⅲ. 退温函数

退温函数即温度更新函数,用于在外循环中修改温度值。目前,最常用的温度更新函数为指数退温函数,即 T ( n + 1 ) = K × T ( n ) T(n+1)=K×T(n) T(n+1)=K×T(n),其中 0 < K < 1 0<K<1 0K1 K K K为一个非常接近于1的常数。


Ⅳ. Markov 链长度

Markov 链长度是在等温条件下进行迭代优化的次数,其选取原则是在衰减参数 T T T 的衰减函数已选定的前提下,还要产生随机数的次数,一般 L 取100~1000。


Ⅴ. 算法终止准则

算法停止的条件。常用的有,温度降低到一定的阈值结束,迭代一定的次数后结束,最优值连续保持不变(或者变化值 < δ <\delta <δ)时停止搜索。


六、仿真实例:模拟退火求解函数最小值

1. 题目

计算函数 f ( x ) = ∑ i = 1 n x i 2 ( − 20 < = x i < = 20 ) f(x)=\sum_{i=1}^n x_i^2 (-20<=x_i<=20) f(x)=i=1nxi2(20<=xi<=20)的最小值,其中个体 x x x的维数为 n = 10 n=10 n=10


2. 分析

这时一个平方和函数,只有一个极小值 x = ( 0 , 0 , . . . , 0 ) x=(0,0,...,0) x=(0,0,...,0),理论最小值 f ( 0 , 0 , . . . , 0 ) = 0 f(0,0,...,0)=0 f(0,0,...,0)=0

Markov链长度初始化为 L = 200 L=200 L=200,衰减参数设置为0.998,补偿因子为 S = 0.01 S=0.01 S=0.01,初始温度 T = 100 T=100 T=100,容差为 Y Z = 1 × 1 0 − 8 YZ=1\times 10^{-8} YZ=1×108(用于判断结束条件,如果最优值连续变化小于 Y Z YZ YZ,那么就终止搜索,输出这个较优解,其它参数上边都已经提到过了);随机产生初始解,并计算目标函数值。


3. 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 = (1 - S) * PreX + S * (rand(D,1) *(Xs-Xx)+Xx);
        %%%%%%%%%%%%%%%%%边界条件处理%%%%%%%%%%%%%%%%%%%%%%%%%%
        for ii=1:D
            if NextX(ii)>Xs || NextX(ii)<Xx
                NextX(ii)=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('最优解变化曲线')

其中,目标函数(适应度函数)fun1的定义如下:

%%%%%%%%%%%%%%%%%%%%%%%%%适应度函数%%%%%%%%%%%%%%%%%%%%%%%
function result=func1(x)
summ=sum(x.^2);
result=summ;

4. 求解结果及分析

求解的适应度变化曲线如下:

在这里插入图片描述

可以看到,随着迭代次数的增加,目标函数值下降速度还是很快的,说明算法收敛速度较好。

最终搜索结果如下:

在这里插入图片描述

可以看到,我们在理论最优解 0 附近找到了近似最优解 0.0056,求解效果还是不错的。

继续运行代码,求解结果如下:

在这里插入图片描述

可以看到求解过程以及最终的结果都不太相同,因此,对于模拟退火这种启发式智能优化算法,求解结果是有一定随机性的。但是,会随着迭代次数的增加,越发趋于稳定。


七、模拟退火的一些改进方向

  • 增加记忆功能
  • 增加升温或重升温过程。
  • 对每一当前状态,采用多次搜索策略,以概率接受区域内的最优状态,而不是标准模拟退火算法的单次比较方式。
  • 与其他搜索机制的算法(如遗传算法、免疫算法等)相结合。可以综合其他方法的优点,提高运行效率和求解质量。2

the end……

模拟退火算法到这里就要结束啦~~到此既是缘分,欢迎您的点赞评论收藏关注我,不迷路,我们下期再见!!

😘😘😘 我是Cherries,一位计算机科班在校大学生,写博客用来记录自己平时的所思所想!
💞💞💞 内容繁杂,又才疏学浅,难免存在错误,欢迎各位大佬的批评指正!
👋👋👋 我们相互交流,共同进步!

:本文由非妃是公主发布于https://blog.csdn.net/myf_666,转载请务必标明原文链接:https://blog.csdn.net/myf_666/article/details/129061318


  1. Kirkpatrick S,Gelatt C,Vecchi M.Optimization by simulated Anealing.Science,1983(220):671-680. ↩︎

  2. 包子阳,智能优化算法及其matlab实例.电子工业出版社. ↩︎

Logo

开放原子开发者工作坊旨在鼓励更多人参与开源活动,与志同道合的开发者们相互交流开发经验、分享开发心得、获取前沿技术趋势。工作坊有多种形式的开发者活动,如meetup、训练营等,主打技术交流,干货满满,真诚地邀请各位开发者共同参与!

更多推荐