群智能算法-模拟退火

一.基本理念

模拟退火算法(Simulated Annealing,简称SA) 的思想最早是由Metropolis等提出的.其出发点就是基于物理中固体物质的退火过程与一般的组合优化问题之间的相似性。模拟退火法是一种通用的优化算法，其物理退火过程由以下三部分组成。

加温过程：我们这一步是为了将粒子变为活性粒子，在一开始温度最大的时候，恰恰是其活性最大的时候，其遍历的地方也就越多，这样便于我们将局部最优解的情况给剔除。

等温过程：对于每一个粒子我们都会向一个自由能变小的方向运行。

冷却过程：当自由能极小的时候我们就相对于说，就是求局部最优解。

加温过程对算法设定初温，等温过程对应算法的Metropolis抽样过程，冷却过程对应控制参数的下降。这里能量的变化

就是目标函数，我们要得到的最优解就是能量最低态。其中Metropolis准则是SA算法收敛于全局最优解的关键所在，

Metropolis准则以一定的概率接受恶化解，这样就使算法跳离局部最优的陷阱。

 

二.算法流程

 

1. 初始化：取初始温度T0足够大，令T = T0，任取初始解S1。

2. 对当前温度T，重复第（3）～（6）步。

3. 对当前解S1随机扰动产生一个新解S2。

4. 计算S2的增量df = f(S2) - f(S1)，其中f(S1)为S1的代价函数。

5. 若df < 0，则接受S2作为新的当前解，即S1 = S2；否则计算S2的接受概率exp(-df/T)， 即随机产生(0,1)区间上均

匀分布的随机数rand，若exp(-df/T) > rand，也接受S2作为新的当前解S1 = S2，否则保留当前解S1。

6. 如果满足终止条件Stop，则输出当前解S1为最优解，结束程序，终止条件Stop通常取为在连续若干个Metropolis

链中新解S2都没有被接受时终止算法或者是设定结束温度；否则按衰减函数衰减T后返回第（2）步。

三.实战(由于此算法相对于其他的算法的操作每那么多较简单直接实战) 

x = -3:0.1:3;
y = 11*sin(x) + 7*cos(5*x);//这个例子是跟一个大佬学的

模拟退火算法--matlab实现简单问题_模拟退火算法matlab_面壁十年忘何图的博客-CSDN博客

我们求在-3到3这个范围中，y函数的最大值。

%% SA 模拟退火: 求解函数y = 11*sin(x) + 7*cos(5*x)在[-3,3]内的最大值(动画演示）
tic
clear; clc
 
%% 绘制函数的图形
x = -3:0.1:3;
y = 11*sin(x) + 7*cos(5*x);
figure
plot(x,y,'b-')
hold on  % 不关闭图形，继续在上面画图

然后我们初始化

%% 参数初始化
narvs = 1; % 变量个数
T0 = 100;   % 初始温度
T = T0; % 迭代中温度会发生改变，第一次迭代时温度就是T0
maxgen = 200;  % 最大迭代次数
Lk = 100;  % 每个温度下的迭代次数
alfa = 0.95;  % 温度衰减系数
x_lb = -3; % x的下界
x_ub = 3; % x的上界
 

%%  随机生成一个初始解
x0 = zeros(1,narvs);
for i = 1: narvs
    x0(i) = x_lb(i) + (x_ub(i)-x_lb(i))*rand(1);    
end
y0 = 11*sin(x0) + 7*cos(5*x0); % 计算当前解的函数值
h = scatter(x0,y0,'*r');  % scatter是绘制二维散点图的函数（这里返回h是为了得到图形的句柄，未来我们对其位置进行更新）
 
%% 定义一些保存中间过程的量，方便输出结果和画图
max_y = y0;     % 初始化找到的最佳的解对应的函数值为y0
MAXY = zeros(maxgen,1); % 记录每一次外层循环结束后找到的max_y (方便画图）

之后就是我们的退火过程了

%% 模拟退火过程
for iter = 1 : maxgen  % 外循环, 我这里采用的是指定最大迭代次数
    for i = 1 : Lk  % 内循环，在每个温度下开始迭代
        y = randn(1,narvs);  % 生成1行narvs列的N(0,1)随机数
        z = y / sqrt(sum(y.^2)); % 根据新解的产生规则计算z
        x_new = x0 + z*T; % 根据新解的产生规则计算x_new的值
        % 如果这个新解的位置超出了定义域，就对其进行调整
        for j = 1: narvs
            if x_new(j) < x_lb(j)
                r = rand(1);
                x_new(j) = r*x_lb(j)+(1-r)*x0(j);
            elseif x_new(j) > x_ub(j)
                r = rand(1);
                x_new(j) = r*x_ub(j)+(1-r)*x0(j);
            end
        end
        x1 = x_new;    % 将调整后的x_new赋值给新解x1
        y1 = Obj_fun1(x1);  % 计算新解的函数值
        if y1 > y0    % 如果新解函数值大于当前解的函数值
            x0 = x1; % 更新当前解为新解
            y0 = y1;
        else
            p = exp(-(y0 - y1)/T); % 根据Metropolis准则计算一个概率
            if rand(1) < p   % 生成一个随机数和这个概率比较，如果该随机数小于这个概率
                x0 = x1; % 更新当前解为新解
                y0 = y1;
            end
        end
        % 判断是否要更新找到的最佳的解
        if y0 > max_y  % 如果当前解更好，则对其进行更新
            max_y = y0;  % 更新最大的y
            best_x = x0;  % 更新找到的最好的x
        end
    end
    MAXY(iter) = max_y; % 保存本轮外循环结束后找到的最大的y
    T = alfa*T;   % 温度下降
    pause(0.01)  % 暂停一段时间(单位：秒)后再接着画图
    h.XData = x0;  % 更新散点图句柄的x轴的数据（此时解的位置在图上发生了变化）
    h.YData = Obj_fun1(x0); % 更新散点图句柄的y轴的数据（此时解的位置在图上发生了变化）
end

然后就是我们要绘图了

 

disp('最佳的位置是：'); disp(best_x)
disp('此时最优值是：'); disp(max_y)
 
pause(0.5)
h.XData = [];  h.YData = [];  % 将原来的散点删除
scatter(best_x,max_y,'*r');  % 在最大值处重新标上散点
title(['模拟退火找到的最大值为', num2str(max_y)])  % 加上图的标题
 
%% 画出每次迭代后找到的最大y的图形
figure
plot(1:maxgen,MAXY,'b-');
xlabel('迭代次数');
ylabel('y的值');
toc

寻找过程

 寻找结果

迭代代数 

 折线图

 

 实战（2）我们依照这个方法我们不妨得一个三维山峦的最大值

初始化

clear gca
%%x=10;
% y=5;
% z=funs(x,y);
[x1,x2]=meshgrid(-5:0.1:5,-5:0.1:5);
y=x1 + x2 - 10*cos(pi*x1) - 10*sin(pi*x2);
%%y=(x1-5).^2+(x2-5).^2;
figure
mesh(x1,x2,y);
hold on
%%初始化参数
xs=zeros(100,2);
ys=zeros(100,1);
num=2;%%变量个数为2；
T0=100;%%初始温度
T=T0;%%每次迭代温度
maxgen=100;%%链长度
Lnum=100;%%每个温度状态下的链长度
alfa=0.998;%%温度衰减
x_lb=-5;
x_ub=5;
s=100;
x0=zeros(1,num);
for i=1:num
    x0(i)=x_lb+10*rand(1);
end
y0=funs(x0(1),x0(2));
h=scatter3(x0(1),x0(2),y0,s,"*r");
disp(h);
%%过程量的定义
max_y=y0;%%最优y值
MAXY=zeros(maxgen,1);%%每一代的

退火过程

%%退火过程
for iter=1:maxgen
    for i=1:Lnum
        y=randn(1,num);%%移动%%温度最大移动的位置范围就越大
        z=y/sqrt(sum(y.^2));%%移动
        x_new=x0+z*T;
        for j=1:num%%约束
            if x_new(j)<x_lb
            r=rand(1);
            x_new(j)=r*x_lb+(1-r)*x0(j);
            elseif x_new(j)>x_ub
                r=rand(1);
                x_new(j)=r*x_ub+(1-r)*x0(j);
            end
        end
        x1=x_new;
    y1=funs(x1(1),x1(2));
    if y1>y0%%退火法则
        x0=x1;
        y0=y1;
    else 
        p=exp(-(y0-y1)/T);
        if rand(1)<p
            x0=x1;
            y0=y1;
        end
    end
    if y0>max_y%%保留最大值
        max_y=y0;
        best_x=x0;
    end
    end
   xs(iter,:)=x0;
   ys(iter)=y0;
   MAXY(iter)=max_y;
   T=alfa*T;
%        h.XData=x0(1);
%     h.YData=x0(2);
%     h.ZData=funs(x0(1),x0(2));
   pause(0.1);
   %%disp(x0);
   %%disp(y0);
     h.XData=x0(1);%%换着输出值
    h.YData=x0(2);
     h.ZData=funs(x0(1),x0(2));
   %%scatter3(x0(1),x0(2),y0,s,"*r")

end

绘图

figure%%绘制其的折线图
plot3(xs(:,1),xs(:,2),ys);
h.XData=best_x(1);
    h.YData=best_x(2);
     h.ZData=max_y;

figure 
title("归");
plot(1:maxgen,MAXY);

寻找过程

 寻找结果

 迭代代数

 变化折线