智能计算方法与实现2|模拟退火算法原理|工具箱及其应用

模拟退火算法原理

模拟退火算法

模拟退火算法（SimulatedAnnealing，SA）最早的思想是由N.Metropolis等人于1953年提出。
1983年，S.Kirkpatrick等成功地将退火思想引l入到组合优化领域
它是基于Monte-Carlo送代求解策略的一种随机寻优算法，其出发点是基于物理中固体物质的退火过程与一般组合优化问题之间的相似性

模拟退火算法从某一较高初温出发，伴随温度参数的不断下降，结合概率突跳特性在解空间中随机寻找目标函数的全局最优解，即在局部最优解能概率性地跳出并最终趋于全局最优。
模拟退火算法是一种通用的优化算法，理论上算法具有概率的全局优化性能，自前已在工程中得到了广泛应用，诸如VLSI、生产调度、控制工程、机器学习、神经网络、信号处理等领域

模拟退火算法是通过赋予搜索过程一种时变且最终趋于零的概率突跳性，从而可有效避免陷入局部极小并最终趋于全局最优的串行结构的优化算法

模拟退火核心思想

模拟退火算法从某一较高初温出发，伴随温度参数的不断下降，结合一定的概率突跳特性在解空间中随机寻找目标函数的全局最优解，即在局部最优解能概率性地跳出并最终趋于全局最优。
这里的”一定的概率”的计算参考了金属冶炼的退火过程，这也是模拟退火算法名称的由来。将温度T当作控制参数，目标函数值f视为内能E，而固体在某温度T时的一个状态对应一个解，然后算法试图随着控制参数T的降低，使目标函数f（内能E）也逐渐降低，直至趋于全局最小值（退火中低温时的最低能量状态），就像金属退火过程一样。

模拟退火数学原理

从上面我们知道，会结合概率突跳特性在解空间中随机寻找目标函数的全局最优解，那么具体的更新解的机制是什么呢？如果新解比当前解更优，则接受新解，否则基于Metropolis准则判断是否接受新解。
接受概率为：
$$P=\{\begin{array}{c}1,E_{t+1}<E_t\\ e^{\frac{-(E_{t+1}-E_t)}{kT}},E_{t+1}\ge E_t\end{array}$$
如上公式，假设当前时刻搜索的解为$x_t$，对应的系统能量（目标函数）为$E_t$，对搜索点施加随机扰动，产生新解$x_{t+1}$，相应地，系统能量为$E_{t+1}$，那么系统对搜索点从$x_t$到$x_{t+1}$转变的接受概率就为上公式

假设开始状态在A，随着送代次数更新到B局部最优解，这时发现更新到B时，能量比A要低，则说明接近最优解了，因此百分百转移，状态到达B后，发现下一步能量上升了，如果是梯度下降则是不充许继续向前的，而这里会以一定的概率跳出这个坑，这个概率和当前的状态、能量等都有关系，如果B最终跳出来了到达C，又会继续以一定的概率跳出来，直到到达D后，就会稳定下来。

模拟退火流程

算法实质分两层循环，在任一温度水平下，随机扰动产生新解，并计算目标函数值的变化，决定是否被接受。由于算法初始温度比较高，这样，使E增大的新解在初始时也可能被接受，因而能跳出局部极小值，然后通过缓慢地降低温度，算法就最终可能收敛到全局最优解，具体流程为：

1. 令$T=T_0$，表示开始退火的初始温度，随机产生一个初始解$x_0$，并计算对应的目标函数值$E(x_0)$
2. 令$T=kT$，其中$k$取值0到1之间，为温度下降速率；
3. 对当前解x施加随机扰动，在其邻域内产生一个新解$x_{t+1}$，并计算对应的目标函数值$E(x_{t+1})$，计算
   $$\Delta E=E(x_{t+1})-E(x_t)$$
4. 若$\Delta E<0$，接受新解作为当前解，否则按照概率$e^{-\Delta E/kT}$判断是否接受新解；
5. 在温度T下，重复L次扰动和接受过程，即执行步骤3和4；
6. 判断温度是否达到终止温度水平，若是则终止算法，否则返回步骤2.
   

其中有几个需要注意的点：

1. 初始点的选取对算法结果有一定的影响，最好是多次运行对结果进行综合判断
2. 在算法运行初期，温度下降快，避免接受过多的差结果。当运行时间增加，温度下降减缓，以便于更快稳定结果。
3. 当送代次数增加到一定次数时，结果可能已经达到稳定，但是距离算法结束还有一段时间。在设计程序时应该加入适当的输出条件，满足输出条件即可结束程序。

模拟退火算法求解TSP

clc; 
clear; 
close all;
%% 
tic
T0 = 1000; %初始温度 
Tend = 1e-3; %终止温度
L = 500; %各温度下的迭代次数（链长）
q = 0.9; %降温速率%%加载数据
load CityPosition;
%%
D = Distanse(X); %计算距离矩阵
N = size(D,1); %城市的个数
%%初始解
S1 = randperm(N); %随机产生一个初始路线%%画出随机解的路径图
DrawPath(S1,x) 
pause(0.0001)
%%输出随机解的路径和总距离 
disp('初始种群中的一个随机值：') 
OutputPath(S1);
Rlength = athLength(D,S1);
disp(['总距离：', num2str(Rlength)])%%计算迭代的次数Time
Time = ceil(exp(log(T0/Tend)/log(q))); 
count=0; %迭代计数
Obj = zeros(Time,1); %目标值矩阵初始化
track = zeros(Time，N); %每代的最优路线矩阵初始化%%迭代
while T0 > Tendcount = count + 1; %更新迭代次数temp = zeros(L, N+1); for k = 1:L%%产生新解S2 = NewAnswer(S1);%%Metropolis法则判断是否接受新解[S1, R] = Metropolis(S1, S2, D, T0);%Metropolis抽样算法temp(k,：) = [S1R]; %记录下一路线的及其路程 end%%记录每次迭代过程的最优路线[d0, index] = min(temp(：end)); %找出当前温度下最优路线 if count == 1 || d0 < Obj(count-1)obj(count) = d0；%如果当前温度下最优路程小于上一路程则记录当前路程 elseObj(count) = Obj(count-1); %如果当前温度下最优路程大于上一路程则记录上一路程 endtrack(count,：) = temp(index,1:end-1); %记录当前温度的最优路线 T0 = q * T0; %降温fprintf(1, '%d\n', count) %输出当前迭代次数
end

functionS2 = NewAnswer(S1)
%%输入
%S1:当前解
%%输出%S2：新解
N = length(S1); 
S2 = S1;
a = round(rand(1,2)*(N-1)+1);%产生两个随机位置用来交换 
W = S2(a(1));
S2(a(1)) = S2(a(2));
S2(a(2)) = W; %得到一个新路线

%%优化过程迭代图 
figure
plot(1: count, Obj) 
xlabel('迭代次数')
ylabel('距离') 
title('优化过程')%%最优解的路径图
DrawPath(track(end, :), X)%%输出最优解的路线和总距离
disp('最优解：')
S = track(end, :); 
p = OutputPath(S);
disp([总距离：num2str(PathLength(D,S)]);
disp('----------------------------------------------') 
toc

模拟退火算法工具箱及其应用

模拟退火算法工具箱

Matlab的Global OptimizationToolbox中集成了模拟退火算法，为了表示方便，我们称为模拟退火算法工具箱（SAT）。
它位于Matlab安装目录

\toolbox\globaloptim\globaloptim

下

下界，上界
[x, fval] = simulannealbnd(objectiveFunction, X0, lb, ub, options)

• x，最优解，
• fval，最优值
• simulannealbnd，主函数
• ObjectiveFunction，目标函数
• x0，初始解
• options，参数集
  参数设置
  options = saoptimset( …
  ‘Maxlter’，500，…%最大送代次数
  ‘StalllterLim’，500，…%算法终止条件’TolFun’,1e-100,…%适应度函数值偏差
  ‘AnnealingFcn’,@annealingfast, …%退火函数
  ‘Initial Temperature’, 100,…%初始温度
  ‘TemperatureFcn’,@temperatureexp,…%指数降温
  'Reanneallnterval，500，…%与Maxlter相同，表示不进行回火处理

'PlotFcns',{@saplotbestx, @saplotbestf, @saplotx, @saplotf,@saplottemperature);

PlotFcns，绘图函数
saplotbestx，最优解
saplotbestf，最优值
saplotx，解x
saplotf，x对应函数值
saplottemprature，温度

$$f(x_1,x_2)=20+x_1^2+x_2^2-10(\cos 2\pi x_1+\cos 2\pi x_2)$$
$$-1\le x_i\le 1,i=1,2$$

clear 
clcObjectiveFunction = @Rastrigin; % Function handle to the objective function 
X0 = [0.5 0.5]; %Starting point
lb = [-1 -1]; %Lowerbound 
ub = [1 1]; %Upperboundoptions = saoptimset(...
'Maxlter', 500,...
'StalllterLim', 500,...
'TolFun', 1e-100,...
'AnnealingFcn', @annealingfast,...
'InitialTemperature', 100,...
'TemperatureFcn', @temperatureexp,...
'Reanneallnterval', 5o0,...
'PlotFcns',{@saplotbestx, @saplotbestf, @saplotx, @saplotf, @saplottemperature});[x, fval] = simulannealbnd(ObjectiveFunction, x0, lb, ub, options)

x = 1.0e-04 * -0.1839    -0.2334
fval = 1.7517e-07

SAT求解压力容器设计问题

clear 
clcObjectiveFunction = @fitnessFunction; %Function handleto theobjective function
X0 = [50504545]; %Startingpoint 
lb = [001010]; %Lowerbound
ub = [100100100100]; %Upperbound
options = saoptimset(...
'Maxlter',500,...
'StalllterLim',500,...
'TolFun', 1e-100,...
'AnnealingFcn', @annealingfast,...
'InitialTemperature',100,...
'TemperatureFcn', @temperatureexp,...
'Reanneallnterval',5o0,...
'PlotFcns',{@saplotbestx, @saplotbestf, @saplotx, @saplotf, @saplottemperature})[x,fval]= simulannealbnd(ObjectiveFunction, X0, lb, ub, options)

x = 1.8048  0.8919  93.4926  33.6489
fval = 2.3785e+04