模拟退火算法简介

2015-12-03

目录 [+]

1. 简介
2. 相关概念
3. 原理
4.算法步骤
5.源码

1. 简介

模拟退火算法 (Simulated Annealing，SA) 最早的思想是由 N. Metropolis 等人于1953年提出。1983年, S. Kirkpatrick 等成功地将退火思想引入到组合优化领域。它是基于 Monte-Carlo 迭代求解策略的一种随机寻优算法，其出发点是基于物理中固体物质的退火过程与一般组合优化问题之间的相似性。模拟退火算法从某一较高初温出发，伴随温度参数的不断下降,结合概率突跳特性在解空间中随机寻找目标函数的全局最优解，即在局部最优解能概率性地跳出并最终趋于全局最优。

2. 相关概念

退火

退火：将固体加热到足够高的温度，使分子呈随即排列状态，然后逐步降温使之冷却，最后分子以低能状态排列，固体达到某种稳定状态。

退火过程：

加温过程：增强粒子运动，消除系统原先可能存在的非均匀态。
等温过程：对于与环境换热而温度不变的封闭系统，系统状态的自发变化总是朝自由能减少的方向进行，当自由能达到最小时，系统达到平衡。
冷却过程：使粒子热运动减弱并渐趋有序，系统能量逐渐下降，从而得到低能的晶体结构。

Boltzman 概率分布

Boltzman 概率分布告诉我们：

同一个温度，分子停留在能量小状态的概率大于停留在能量大状态的概率
温度越高，不同能量状态对应的概率相差越小，温度足够高时，各状态对应概率基本相同。
随着温度的下降，能量最低状态对应概率越来越大，温度趋于0时，其状态趋于1

Metropolis 准则

固体在恒定温度下达到热平衡的过程可以用 Metropolis 方法加以模拟。

温度恒定时：若在温度T，当前状态i –> 新状态j，若$E_j < E_i$,则接受j为当前状态；否则，若概率 $p=e^{-(E_j-E_i)/(k_B*T)}$ 大于[0,1)区间的随机数，则仍接受状态j为当前状态；若不成立则保留状态i为当前状态。

温度变化时：$p=e^{-(E_j-E_i)/(k_B*T)}$ ,在高温下，可接受当前状态能量差较大的新状态；在低温下，只接受与当前状态能量差较小的新状态。

3. 原理

模拟退火算法来源于固体退火原理，将固体加温至充分高，再让其徐徐冷却，加温时，固体内部粒子随温升变为无序状，内能增大，而徐徐冷却时粒子渐趋有序，在每个温度都达到平衡态，最后在常温时达到基态，内能减为最小。根据 Metropolis 准则，粒子在温度T时趋于平衡的概率为e(-ΔE/(kT))，其中E为温度T时的内能，ΔE为其改变量，k为 Boltzmann 常数。用固体退火模拟组合优化问题，将内能E模拟为目标函数值f，温度T演化成控制参数t，即得到解组合优化问题的模拟退火算法：由初始解i和控制参数初值t开始，对当前解重复“产生新解→计算目标函数差→接受或舍弃”的迭代，并逐步衰减t值，算法终止时的当前解即为所得近似最优解，这是基于蒙特卡罗迭代求解法的一种启发式随机搜索过程。退火过程由冷却进度表(Cooling Schedule)控制，包括控制参数的初值t及其衰减因子Δt、每个t值时的迭代次数L和停止条件S。

4.算法步骤

step1

由一个产生函数从当前解产生一个位于解空间的新解；为便于后续的计算和接受，减少算法耗时，通常选择由当前新解经过简单地变换即可产生新解的方法，如对构成新解的全部或部分元素进行置换、互换等，注意到产生新解的变换方法决定了当前新解的邻域结构，因而对冷却进度表的选取有一定的影响。

step2

计算与新解所对应的目标函数差。因为目标函数差仅由变换部分产生，所以目标函数差的计算最好按增量计算。事实表明，对大多数应用而言，这是计算目标函数差的最快方法。

step3

判断新解是否被接受,判断的依据是一个接受准则，最常用的接受准则是 Metropolis 准则: 若Δt′<0则接受S′作为新的当前解S，否则以概率exp(-Δt′/T)接受S′作为新的当前解S。

step4

当新解被确定接受时，用新解代替当前解，这只需将当前解中对应于产生新解时的变换部分予以实现，同时修正目标函数值即可。此时，当前解实现了一次迭代。可在此基础上开始下一轮试验。而当新解被判定为舍弃时，则在原当前解的基础上继续下一轮试验。

5.源码

MATLAB

function [xo,fo] = Opt_Simu(f,x0,u,l,kmax,q,TolFun)
% 模拟退火算法求函数 f(x)的最小值点， 且 l <= x <= u
% f为待求函数，x0为初值点，l，u分别为搜索区间的上下限，kmax为最大迭代次数
% q为退火因子，TolFun为函数容许误差
%%%%算法第一步根据输入变量数，将某些量设为缺省值
if nargin < 7
    TolFun = 1e-8;
end
if nargin < 6
    q = 1;
end
if nargin < 5
    kmax = 100;
end
%%%%算法第二步，求解一些基本变量
N = length(x0); %自变量维数
x = x0;
fx = feval(f,x); %函数在初始点x0处的函数值
xo = x;
fo = fx;
%%%%%算法第三步，进行迭代计算，找出近似全局最小点
for k =0:kmax
    Ti = (k/kmax)^q;
    mu = 10^(Ti*100); % 计算mu
    dx = Mu_Inv(2*rand(size(x))-1,mu).*(u - l);%步长dx
    x1 = x + dx; %下一个估计点
    x1 = (x1 < l).*l +(l <= x1).*(x1 <= u).*x1 +(u < x1).*u; %将x1限定在区间[l,u]上
    fx1 = feval(f,x1);
    df = fx1- fx;
    if df < 0||rand < exp(-Ti*df/(abs(fx) + eps)/TolFun) %如果fx1<fx或者概率大于随机数z
        x = x1;
        fx = fx1;
    end
    if fx < fo
        xo = x;
        fo = fx1;
    end
end
function x = Mu_Inv(y,mu)
x = (((1+mu).^abs(y)- 1)/mu).*sign(y);