数学建模之模拟退火算法(含matlab代码)

前言

本文以面向数学建模的同学为主，不适合想要深入研究此算法的同学，不涉及过多专业知识，尽量让大家在短时间内快速了解这个算法的大体内容，并能用matlab编程加以实现，有时间的话可以自己把代码敲一遍，代码关键部分都有注释，对自己的matlab编程能力还是有一定帮助的，没有时间的话，我也在文中介绍了哪些是可以直接套用的，哪些是需要依据具体问题修改的，有任何疑问，欢迎留言

模拟退火算法

简介

模拟退火算法来源于固体退火原理，是一种基于概率的算法，将固体加温至充分高，再让其徐徐冷却，加温时，固体内部粒子随温升变为无序状，内能增大，而徐徐冷却时粒子渐趋有序，在每个温度都达到平衡态，最后在常温时达到基态，内能减为最小

算法主体

概括

该算法有两个循环条件,一个是依据种群规模而定的L,一个是最低温度,达到L或者达到最低温度则结束循环，在循环内进行的操作是,我们有一个初始解，对初始解进行扰动产生新解，如果新解优于初始解，则用新解继续循环，在一定概率内接受较差解，防止局部最优，接受差解的概率随着温度的降低而降低，在整个算法过程中，温度由初始温度和温度衰减系数决定逐渐降低，整个解逐步稳定

算法难点

如何扰动产生新解是整个算法的难点,本文对此问题采取的方法是，我们对于该TSP问题的目标函数是总路程最短,比如有10个我们要经过的城市，路线为0->1->2->3->4->5->6->7->8->9->10->0,0表示的是初始点和结束点,我们对其产生的扰动是,在起始点和终点之间随机选取两个下标，在其中间的路线逆序,比如我们选择2和7,扰动过后路线为0->1->2->6->5->4->3->7->8->9->10->0,路线变化后整个总路程变化为d(2,3)+d(6,7)-d(2,6)-d(3,7).其中d(x,y)表示x到y的距离,如果该路程变化小于0,则接受新路线

各个变量设置

• T:初始温度，可直接同本文直接设置为1
• L:循环终止条件之一,一般为问题规模*100
• at:温度衰减系数,可直接同本文0.99
• e:最低温度,循环终止条件之一，可直接同本文10^(-30)

算法总结

该算法的最难点就是如何产生新解，其余几乎都一样，没啥太大变化，就是不断迭代，求最优的过程，还要稍微关注一下接受较差解的条件，本算法可以做的图也有两个，一个是最终结果图，一个是不断迭代，目标函数值趋近最优的变化图像，该算法也是个优化算法，每次运行的结果可能不同，进行该算法时，对初值稍微筛选一下，效果更好

问题和代码

问题

一个飞机从(70,40)出发依次经过100个点，然后返回起点，求最短路线(给出的数据坐标是经纬度)

数据

53.7121 15.3046 51.1758 0.0322 46.3253 28.2753 30.3313 6.9348
56.5432 21.4188 10.8198 16.2529 22.7891 23.1045 10.1584 12.4819
20.1050 15.4562 1.9451 0.2057 26.4951 22.1221 31.4847 8.9640
26.2418 18.1760 44.0356 13.5401 28.9836 25.9879 38.4722 20.1731
28.2694 29.0011 32.1910 5.8699 36.4863 29.7284 0.9718 28.1477
8.9586 24.6635 16.5618 23.6143 10.5597 15.1178 50.211 10.2944
8.1519 9.5325 22.1075 18.5569 0.1215 18.8726 48.2077 16.8889
31.9499 17.6309 0.7732 0.4656 47.4134 23.7783 41.8671 3.5667
43.5474 3.9061 53.3524 26.7256 30.8165 13.4595 27.7133 5.0706
23.9222 7.6306 51.9612 22.8511 12.7938 15.7307 4.9568 8.3669
21.5051 24.0909 15.2548 27.2111 6.2070 5.1442 49.2430 16.7044
17.1168 20.0354 34.1688 22.7571 9.4402 3.9200 11.5812 14.5677
52.1181 0.4088 9.5559 11.4219 24.4509 6.5634 26.7213 28.5667
37.5848 16.8474 35.6619 9.9333 24.4654 3.1644 0.7775 6.9576
14.4703 13.6368 19.8660 15.1224 3.1616 4.2428 18.5245 14.3598
58.6849 27.1485 39.5168 16.9371 56.5089 13.7090 52.5211 15.7957
38.4300 8.4648 51.8181 23.0159 8.9983 23.6440 50.1156 23.7816
13.7909 1.9510 34.0574 23.3960 23.0624 8.4319 19.9857 5.7902
40.8801 14.2978 58.8289 14.5229 18.6635 6.7436 52.8423 27.2880
39.9494 29.5114 47.5099 24.0664 10.1121 27.2662 28.7812 27.6659
8.0831 27.6705 9.1556 14.1304 53.7989 0.2199 33.6490 0.3980
1.3496 16.8359 49.9816 6.0828 19.3635 17.6622 36.9545 23.0265
15.7320 19.5697 11.5118 17.3884 44.0398 16.2635 39.7139 28.4203
6.9909 23.1804 38.3392 19.9950 24.6543 19.6057 36.9980 24.3992
4.1591 3.1853 40.1400 20.3030 23.9876 9.4040 41.1084 27.7149
第一列和第二列表示1到25个点的坐标，第三列和第四列表示26到50个点的坐标，以此类推,可直接复制到txt文档,文档命名为sj.txt，然后运行程序观察结果

matlab程序

sj0 = load('sj.txt'); % 载入目标数据
x = sj0(:,(1 : 2 : 8));
x = x(:); % 把25*4换成100*1,为了后面的合并
y = sj0(:,(2: 2: 8)); 
y = y(:);
sj = [x, y]; 
d1 = [70, 40];
sj = [d1; sj; d1]; % 合并初始点、中间点、终点
sj = sj * pi / 180; % 角度转化为弧度

% 以下是距离矩阵,存储任意两点之间距离
d = zeros(102);
for i = 1 : 101
    for j = i + 1 : 102
        d(i, j) = 6370 * acos(cos(sj(i, 1) - sj(j, 1)) * cos(sj(i, 2)) * cos(sj(j, 2)) ...
    + sin(sj(i, 2)) * sin(sj(j, 2)));
    end
end
d=d + d'; %d最后是个对称阵，为减少循环次数，所以加上其转置

% 以下确定较好的初始解
path = [];
long = inf;
rand('state', sum(clock)); %用时间设立随机数种子,防止随机数不变
for j = 1 : 1000
    path0 = [1, 1 + randperm(100), 102]; %randperm(100),1到100随机排列
    temp = 0;
    for i = 1 : 101
        temp = temp + d(path0(i), path0(i + 1)); %记录总行驶距离
    end
    if temp < long %更新较优解和新的路线
        path = path0;
        long = temp;
    end
end

% 以下是退火过程
e = 1e-30; %定的
L = 20000; %根据数据规模确定，一般为100*规模
at = 0.999; %定的
T = 1;
for k = 1 : L
    c = 2 + floor(100*rand(1,2)); % 产生新解，这两个解中的所有路线逆序
    c = sort(c); %升序，c为两个矩阵下标，必须升序，否则后面不能逆序改变
    c1 = c(1);
    c2 = c(2);
    % 路径长度改变量
    df = d(path(c1 - 1), path(c2)) + d(path(c1), path(c2 + 1)) - d(path(c1 - 1) ,path(c1)) ...
    -d(path(c2), path(c2 + 1));
    if df < 0 % 接受新路径的条件，总距离变小
        path = [path(1 : c1 - 1), path(c2 : -1 : c1), path(c2 + 1 : 102)];
        long = long + df;
    elseif exp(-df / T) >= rand %随机接收非最优解，防止陷入局部最优
        path = [path(1 : c1 - 1), path(c2 : -1 : c1), path(c2 + 1 : 102)];
        long = long + df;
    end
    T = T * at;
    if T < e %达到最低温度结束循环
        break;
    end
end
path, long % 路径及路径长度
xx = sj(path, 1);
yy = sj(path, 2);
plot(xx, yy, '-*') % 绘制路径

结果展示