模拟退火算法(SA)

发布时间 2023-07-18 01:46:32作者: cxy8

求某个目标函数的最值

爬山法

首先我们通过爬山法来引出模拟退火算法
我们先看一个例子:求函数的最值
1.png
我们用爬山法解决这个问题的步骤
1、在解空间中随机生成一个初始解(图中小黄点就是我们生成的初始解)
2.png
2、根据初始解的位置,我们有两种决策,一种是向左邻域,一种是向右邻域走一步(走的步长越小越好,但是越小的话计算量会更大)
3、比较不同走法下目标函数的大小,决定下一步往哪个方向走,在上面的图上是向左走,目标函数的值更大,因此我们应该向左走一步
4、不断重复这个步骤,直到我们找到一个极大值点,此时我们结束搜索。

只要思考一下,我们可以明显的看出爬山法的缺点是特别容易找到局部最优解。爬山法的本质是贪心算法。这里我们的主角马上就要登场了,模拟退火算法能很好地解决上面的问题,模拟退火引入了一个概率帮助我们跳出局部最优解从而搜索到全局最优解。

模拟退火算法

3.png
模拟退火的核心思想是以一定的概率接收这个解
解决方法:
为了不直接拒绝xj,我们定义一个接收的概率p,p位于[0, 1]之间,且f(xj)和f(xi)越接近,p越大(这种想法很直观,新解和旧解的函数值越接近我们就越愿意接受新解)
p 正比于|1/(f(xj)- f(xi))|
4.png

这里当p=0时我们对应的不就是爬山法吗?,当p=1时,我们对应的就是蒙特卡洛(枚举)

5.png
6.png
7.png
模拟退火解决最值问题的步骤
8.png
9.png

1、在生成随机解的时候进行判断一下,当不满足的约束条件的时候之间再生成一个
第二种办法是加罚函数,当不满足约束条件时

10.png
这里就体现出来了模拟退火算法名字的来源,退火是温度随时间递减,但是这里我们Ct想要得到的是一个关于时间递增的函数所以我们这里就取了一下倒数。
这里的t怎么再编程中体现,很简单,这里的t我们可以看成迭代次数(循环)
11.png

12.png

另一个问题是什么时候停止搜索?
这里的准则并不唯一
1、我们可以设置达到指定的迭代次数,例如迭代200ci
2、达到指定的温度
3、我们找到最优解连续M(例如30次)次迭代都没有变化

如何生成新解呢?
13.png
模拟退火算法代码

%% 绘制函数图形
x = -3:0.1:3
y = 11*sin(x) + 7*cos(5*x);
figure(1);
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 = x_lb + (x_ub - x_lb)*rand(1, narvs);
% scatter是绘制二维散点图的函数(这里返回h是为了获得函数的句柄,未来我们对其位置进行更新)
h = scatter(x0, Obj_fun1(x0), "*r");

%% 初始化用来保存中间结果的x和y的取值
iter_x = x0;
iter_y = Obj_fun1(x0);

%% 模拟退火过程
for iter = 1 : maxgen % 外循环,我们这里采用的是指定最大迭代次数
   for i = 1 : LK % 内循环,在每个温度下开始迭代
       title(['当前温度:', num2str(T)])
       y0 = Obj_fun1(x0); % 计算当前解的函数值
       disp('**********************************分割线********************************')
       disp(['当前解的函数值:',num2str(y0)])
       y = randn(1, narvs); % 生成一行narvs列的N(0,1)随机数
       z = y / sqrt(sum(y .^ 2)) % 根据新解的产生规则计算x_new的值
       x_new = x0 + z * T; % 根据新解的产生挥着计算z
       % 如果这个新解的位置超出了定义域,就对其进行调整
       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); % 计算新解的函数值
       disp(['新解的函数值:', num2str(y1)])
       if y1 > y0 % 如果新解的函数值大于当前解的函数值
           disp('更新当前解')
           x0 = x1;
           iter_x = [iter_x; x1];
           iter_y = [iter_y; y1];
       else
           p = exp(-(y0 - y1) / T) % 根据MEtropolis准则计算一个概率
           disp(['新解的函数值小于当前解的函数值,接受新解的概率为:', num2str(p)])
           if rand(1) < p % 生成一个随机数和这个概率比较,如果该随机数小于这个概率
               disp('该随机数小于这个概率,那么我们接受新解')
               x0 = x1; % 更新当前解为新解
           end
       end
    h.XData = x0; % 更新散点图句柄的x轴数据
    h.YData = Obj_fun1(x0);
   end

   T = alfa*T; % 温度下降
   pause(0.01) % 暂停一段时间(单位:秒)后再接着画图
end

其中用到的函数代码

function y = Obj_fun1(x)
    y = 11*sin(x) + 7*cos(5*x);
end

TSP问题(旅行商问题)