数学建模暑期集训23：模拟退火算法
【摘要】
模拟退火算法类属
和粒子群算法一样，模拟退火算法也属于启发式算法的一种。 启发式算法，可参照下面的定义。 启发式算法：在搜索最优解的过程中利用到了原来搜索过程中得到的信息，且这个信息会改进我们的搜索过程...
模拟退火算法类属
和粒子群算法一样，模拟退火算法也属于启发式算法的一种。
启发式算法，可参照下面的定义。
启发式算法：在搜索最优解的过程中利用到了原来搜索过程中得到的信息，且这个信息会改进我们的搜索过程。
爬山法
模拟退火算法，可以算一种优化过的爬山法。
爬山法比较好理解，首先在函数图上随机选取一个点，之后再其左边或右边各选一点，若比该点大，则以大的点继续选择，整个过程类似于爬山。
问题在于，当爬到小山峰的时候，无法继续爬，这就导致陷入局部最优解。
模拟退火算法流程
模拟退火在爬山法的基础上，结合蒙特卡洛的思想，整个流程如下：
C t {C_t} Ct的设置
在上面的算法流程中， C t {C_t} Ct并没有具体设置， C t {C_t} Ct控制是搜索范围，当 C t {C_t} Ct越大时，接受B的概率越小，即目标点越不容易动。
因此， C t {C_t} Ct是关于t的递增函数，这样就使前期尽可能地在全局搜索，这样就会防止陷入局部最优解；后期减小搜索范围，提高搜索速度。（可以理解为渣男思想：前期广撒网，后期精确突击）
C t {C_t} Ct的设置，模拟的是退火原理。
引入 C t {C_t} Ct后的模拟退火算法流程
在定义好 C t {C_t} Ct的表达式之后，可以将模拟退火算法的流程进行补充。
算法循环可以设置为2层，即第一次在高温t的情况下进行遍历，之后逐渐降低温度t，也就是改变 C t {C_t} Ct，再次遍历。
例题1：求给定函数的最大值或者最小值
目标函数Obj_fun1
function y = Obj_fun1(x)
y = 11*sin(x) + 7*cos(5*x);
end
- 1
- 2
- 3
模拟退火（含动画演示）
%% 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的值');
toc
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
- 41
- 42
- 43
- 44
- 45
- 46
- 47
- 48
- 49
- 50
- 51
- 52
- 53
- 54
- 55
- 56
- 57
- 58
- 59
- 60
- 61
- 62
- 63
- 64
- 65
- 66
- 67
- 68
- 69
- 70
- 71
- 72
- 73
- 74
- 75
- 76
- 77
- 78
- 79
- 80
- 81
- 82
- 83
- 84
- 85
- 86
- 87
- 88
例题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.833300,41836.111100;11751.111100,42814.444400;11770.277800,42651.944400;11785.277800,42884.444400;11822.777800,42673.611100;11846.944400,42660.555600;11963.055600,43290.555600;11973.055600,43026.111100;12058.333300,42195.555600;12149.444400,42477.500000;12286.944400,43355.555600;12300.000000,42433.333300;12355.833300,43156.388900;12363.333300,43189.166700;12372.777800,42711.388900;12386.666700,43334.722200;12421.666700,42895.555600;12645.000000,42973.333300];
n = size(coord,1); % 城市的数目
figure % 新建一个图形窗口
plot(coord(:,1),coord(:,2),'o'); % 画出城市的分布散点图
hold on % 等一下要接着在这个图形上画图的
d = zeros(n); % 初始化两个城市的距离矩阵全为0
for i = 2:n
for j = 1:i
coord_i = coord(i,:); x_i = coord_i(1); y_i = coord_i(2); % 城市i的横坐标为x_i，纵坐标为y_i
coord_j = coord(j,:); x_j = coord_j(1); y_j = coord_j(2); % 城市j的横坐标为x_j，纵坐标为y_j
d(i,j) = sqrt((x_i-x_j)^2 + (y_i-y_j)^2); % 计算城市i和j的距离
end
end
d = d+d'; % 生成距离矩阵的对称的一面
%% 参数初始化
T0 = 1000; % 初始温度
T = T0; % 迭代中温度会发生改变，第一次迭代时温度就是T0
maxgen = 1000; % 最大迭代次数
Lk = 500; % 每个温度下的迭代次数
alpfa = 0.95; % 温度衰减系数
%% 随机生成一个初始解
path0 = randperm(n); % 生成一个1-n的随机打乱的序列作为初始的路径
result0 = calculate_tsp_d(path0,d); % 调用我们自己写的calculate_tsp_d函数计算当前路径的距离
%% 定义一些保存中间过程的量，方便输出结果和画图
min_result = result0; % 初始化找到的最佳的解对应的距离为result0
RESULT = zeros(maxgen,1); % 记录每一次外层循环结束后找到的min_result (方便画图）
%% 模拟退火过程
for iter = 1 : maxgen % 外循环, 我这里采用的是指定最大迭代次数
for i = 1 : Lk % 内循环，在每个温度下开始迭代
path1 = gen_new_path(path0); % 调用我们自己写的gen_new_path函数生成新的路径
result1 = calculate_tsp_d(path1,d); % 计算新路径的距离
%如果新解距离短，则直接把新解的值赋值给原来的解
if result1 < result0
path0 = path1; % 更新当前路径为新路径
result0 = result1;
else
p = exp(-(result1 - result0)/T); % 根据Metropolis准则计算一个概率
if rand(1) < p % 生成一个随机数和这个概率比较，如果该随机数小于这个概率
path0 = path1; % 更新当前路径为新路径
result0 = result1;
end
end
% 判断是否要更新找到的最佳的解
if result0 < min_result % 如果当前解更好，则对其进行更新
min_result = result0; % 更新最小的距离
best_path = path0; % 更新找到的最短路径
end
end
RESULT(iter) = min_result; % 保存本轮外循环结束后找到的最小距离
T = alpfa*T; % 温度下降
end
disp('最佳的方案是：'); disp(mat2str(best_path))
disp('此时最优值是：'); disp(min_result)
best_path = [best_path,best_path(1)]; % 在最短路径的最后面加上一个元素，即第一个点（我们要生成一个封闭的图形）
n = n+1; % 城市的个数加一个（紧随着上一步）
for i = 1:n-1
j = i+1;
coord_i = coord(best_path(i),:); x_i = coord_i(1); y_i = coord_i(2);
coord_j = coord(best_path(j),:); x_j = coord_j(1); y_j = coord_j(2);
plot([x_i,x_j],[y_i,y_j],'-b') % 每两个点就作出一条线段，直到所有的城市都走完
% pause(0.02) % 暂停0.02s再画下一条线段
hold on
end
%% 画出每次迭代后找到的最短路径的图形
figure
plot(1:maxgen,RESULT,'b-');
xlabel('迭代次数');
ylabel('最短路径');
toc
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
- 41
- 42
- 43
- 44
- 45
- 46
- 47
- 48
- 49
- 50
- 51
- 52
- 53
- 54
- 55
- 56
- 57
- 58
- 59
- 60
- 61
- 62
- 63
- 64
- 65
- 66
- 67
- 68
- 69
- 70
- 71
- 72
- 73
- 74
- 75
- 76
- 77
- 78
- 79
- 80
- 81
- 82
- 83
- 84
- 85
- 86
- 87
- 88
- 89
总结
模拟退火作为智能优化算法（启发式算法）中的一种，运行速度远高于蒙特卡洛法。
缺点：matlab无通用内置函数，应对不同问题需要根据需要自行编写代码。
文章来源: zstar.blog.csdn.net，作者：zstar-_，版权归原作者所有，如需转载，请联系作者。
原文链接：zstar.blog.csdn.net/article/details/119939080
【版权声明】本文为华为云社区用户转载文章，如果您发现本社区中有涉嫌抄袭的内容，欢迎发送邮件进行举报，并提供相关证据，一经查实，本社区将立刻删除涉嫌侵权内容，举报邮箱：
cloudbbs@huaweicloud.com
- 点赞
- 收藏
- 关注作者
评论（0）