数学建模暑期集训23：模拟退火算法

Posted 2022-08-24 zstar_

模拟退火算法类属

和粒子群算法一样，模拟退火算法也属于启发式算法的一种。
启发式算法，可参照下面的定义。
启发式算法：在搜索最优解的过程中利用到了原来搜索过程中得到的信息，且这个信息会改进我们的搜索过程。

爬山法

模拟退火算法，可以算一种优化过的爬山法。

爬山法比较好理解，首先在函数图上随机选取一个点，之后再其左边或右边各选一点，若比该点大，则以大的点继续选择，整个过程类似于爬山。

问题在于，当爬到小山峰的时候，无法继续爬，这就导致陷入局部最优解。

模拟退火算法流程

模拟退火在爬山法的基础上，结合蒙特卡洛的思想，整个流程如下：

C t C_t Ct​的设置

在上面的算法流程中，并没有具体设置，控制是搜索范围，当越大时，接受B的概率越小，即目标点越不容易动。
因此，是关于t的递增函数，这样就使前期尽可能地在全局搜索，这样就会防止陷入局部最优解；后期减小搜索范围，提高搜索速度。（可以理解为渣男思想：前期广撒网，后期精确突击）

的设置，模拟的是退火原理。

引入 C t C_t Ct​后的模拟退火算法流程

在定义好的表达式之后，可以将模拟退火算法的流程进行补充。

算法循环可以设置为2层，即第一次在高温t的情况下进行遍历，之后逐渐降低温度t，也就是改变，再次遍历。

例题1：求给定函数的最大值或者最小值

目标函数Obj_fun1

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

模拟退火（含动画演示）

%% 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 = Obj_fun1(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的值);

例题2：旅行商问题

tic
rng(shuffle)  % 控制随机数的生成，否则每次打开matlab得到的结果都一样
% https://ww2.mathworks.cn/help/matlab/math/why-do-random-numbers-repeat-after-startup.html
% https://ww2.mathworks.cn/help/matlab/ref/rng.html
clear;clc
% 38个城市，TSP数据集网站(http://www.tsp.gatech.edu/world/djtour.html) 上公测的最优结果6656。
coord = [11003.611100,42102.500000;11108.611100,42373.888900;11133.333300,42885.833300;11155.833300,42712.500000;11183.333300,42933.333300;11297.500000,42853.333300;11310.277800,42929.444400;11416.666700,42983.333300;11423.888900,43000.277800;11438.333300,42057.222200;11461.111100,43252.777800;11485.555600,43187.222200;11503.055600,42855.277800;11511.388900,42106.388900;11522.222200,42841.944400;11569.444400,43136.666700;11583.333300,43150.000000;11595.000000,43148.055600;11600.000000,43150.000000;11690.555600,42686.666700;11715.833
以上是关于数学建模暑期集训23：模拟退火算法的主要内容，如果未能解决你的问题，请参考以下文章
 
 数学建模暑期集训18：经纬度转换为平面坐标
 数学建模暑期集训18：粒子群算法
 数学建模暑期集训26：遗传算法
 数学建模暑期集训28：元胞自动机
 数学建模暑期集训17：蒙特卡洛法
 数学建模暑期集训5：matlab求解常微分方程/偏微分方程