模拟退火算法
一、模拟退火算法简介
模拟退火算法（Simulated Annealing, SA）是一种通用的随机搜索算法，是对局部搜索算法的扩展。与一般局部搜索算法不同，SA以一定的概率选择邻域中目标值相对较小的状态，是一种理论上的全局最优算法。
模拟退火算法是源于对热力学中退火过程的模拟，在某一给定初温下，通过缓慢下降温度参数，使算法能够在多项式时间内给出一个近似最优解。
二、模拟退火算法与金属退火过程对比
物理退火过程：
加温过程----固态溶解为液态的过程，分子的分布从有序的晶体态转化为无序的液态，消除原有的非均匀态。
等温过程----在退火中，需要保证系统在每一个温度下都达到充分的平衡状态。
冷却过程----液态凝固成固态晶体的过程。系统能量逐渐下降，从而得到低能有序的晶体结构。
模拟退火过程：
设定初始高温，相当于物理退火的加温过程。初始温度要足够高，在实际应用中，要根据以往的经验，通过反复实验来确定T0的值。
热平衡达到，相当于物理退火的等温过程。是指在一个给定温度下，SA用特殊的抽样策略进行随机搜索，最终达到平衡状态的过程。这是SA算法的内循环过程。
降温函数，相当于物理退火的冷却过程。用来控制温度的下降方式，这是SA算法的外循环过程。常用的降温函数有Tk+1=Tk-DT，Tk+1=Tk*r，其中r∈(0.95,0.99)。
1953年，Metropolis等提出了一种重要性采样法(Metropolis准则)，即以概率接受新状态。
Metropolis准则
在温度t，由当前状态i产生新状态j，两者的能量分别为Ei 和 Ej ,若Ei >Ej ，则接受新状态j为当前状态；否则，以一定概率：
来接受状态j，其中k为Boltzmann常数。当这种过程多次重复后，系统将趋于能量较低的平衡态。
三、模拟退火算法流程图
四、模拟退火算法求解TSP问题
参数的设定
•初始解:随机选定一条路径，如[A,B,C,D,E,A]
•初始温度T0=100
•终止温度Tf=20
•降温函数Tk+1=Tk-DT，DT=40
•热平衡达到:内循环次数n=3
代码实现如下:
%======================================================================================
%%%一个旅行商人要拜访全国31个省会城市，需要选择最短的路径%%%%
%%%模拟退火算法解决TSP问题%%%%%%%
clear all; %清除所有变量
close all; %清图
clc ; %清屏
C=[
1304 2312;
3639 1315;
4177 2244;
3712 1399;
3488 1535;
3326 1556;
3238 1229;
4196 1004;
4312 790;
4386 570;
3007 1970;
2562 1756;
2788 1491;
2381 1676;
1332 695;
3715 1678;
3918 2179;
4061 2370;
3780 2212;
3676 2578;
4029 2838;
4263 2931;
3429 1908;
3507 2367;
3394 2643;
3439 3201;
2935 3240;
3140 3550;
2545 2357;
2778 2826;
2370 2975
];%31个省会坐标
n=size(C,1); %TSP问题的规模，即城市数目
T=100*n; %初始温度
L=100; %马尔科夫链的长度
K=0.99; %衰减参数
%%%城市坐标结构体%%%%%%%
city=struct([]);
for i=1:n
city(i).x=C(i,1);
city(i).y=C(i,2);
end
l=1; %统计迭代次数
len(l)=func5(city,n); %每次迭代后路线的长度
figure(1);
while T>0.001
%%%%%%%%%%%%%%%多次迭代扰动，温度降低前多次试验%%%%%%%%
for i=1:L
%%%%%%%%%%%%%%%计算原路线总距离%%%%%%%%%
len1=func5(city,n);
%%%%%%%%%%%%%%%产生随机扰动%%%%%%%%%
%%%%%%%%%%%%%%%随机置换两个不同的城市的坐标%%%%%%%%%
p1=floor(1+n*rand);
p2=floor(1+n*rand);
while p1==p2
p1=floor(1+n*rand);
p2=floor(1+n*rand);
end
tmp_city=city;
%%交换元素
tmp=tmp_city(p1);
tmp_city(p1)=tmp_city(p2);
tmp_city(p2)=tmp;
%%%%%%%%%%%%%%%计算新路线总距离%%%%%%%%%
len2=func5(tmp_city,n);
%%%%%%%%%%%%%%%新老距离的差值，相当于能量%%%%%%%%%
delta_e=len2-len1;
%%%%%%%%%%%%%%%新路线好于旧路线，用新路线替代旧路线%%%%%%%%%
if delta_e<0
city=tmp_city;
else
%%%%%%%%%%%%%%%以一定概率选择是否接受%%%%%%%%%
if exp(-delta_e/T)>rand()
city=tmp_city;
end
end
end
l=l+1; %迭代次数加1
%%%%%%%%%%%%%%%计算新路线的距离%%%%%%%%%
len(l)=func5(city,n);
%%%%%%%%%%%%%%%温度不断下降%%%%%%%%%
T=T*K;
for i=1:n-1
plot([city(i).x,city(i+1).x],[city(i).y,city(i+1).y],'bo-');
hold on;
end
plot([city(n).x,city(1).x],[city(n).y,city(1).y],'ro-');
title(['优化最短距离:',num2str(len(l))]);
hold off;
pause(0.005);
end
figure(2)
plot(len);
xlabel('迭代次数');
ylabel('目标函数值');
title('适应度进化曲线');
%计算距离的函数
function len=func5(city,n)
len=0;
for i=1:n-1
len=len+sqrt((city(i).x-city(i+1).x)^2+(city(i).y-city(i+1).y)^2);
end
len=len+sqrt((city(n).x-city(1).x)^2+(city(n).y-city(1).y)^2);
end
%%%模拟退火算法解决TSP问题%%%%%%%
clear all; %清除所有变量
close all; %清图
clc ; %清屏
C=[
1304 2312;
3639 1315;
4177 2244;
3712 1399;
3488 1535;
3326 1556;
3238 1229;
4196 1004;
4312 790;
4386 570;
3007 1970;
2562 1756;
2788 1491;
2381 1676;
1332 695;
3715 1678;
3918 2179;
4061 2370;
3780 2212;
3676 2578;
4029 2838;
4263 2931;
3429 1908;
3507 2367;
3394 2643;
3439 3201;
2935 3240;
3140 3550;
2545 2357;
2778 2826;
2370 2975
];%31个省会坐标
n=size(C,1); %TSP问题的规模，即城市数目
T=100*n; %初始温度
L=100; %马尔科夫链的长度
K=0.99; %衰减参数
%%%城市坐标结构体%%%%%%%
city=struct([]);
for i=1:n
city(i).x=C(i,1);
city(i).y=C(i,2);
end
l=1; %统计迭代次数
len(l)=func5(city,n); %每次迭代后路线的长度
figure(1);
while T>0.001
%%%%%%%%%%%%%%%多次迭代扰动，温度降低前多次试验%%%%%%%%
for i=1:L
%%%%%%%%%%%%%%%计算原路线总距离%%%%%%%%%
len1=func5(city,n);
%%%%%%%%%%%%%%%产生随机扰动%%%%%%%%%
%%%%%%%%%%%%%%%随机置换两个不同的城市的坐标%%%%%%%%%
p1=floor(1+n*rand);
p2=floor(1+n*rand);
while p1==p2
p1=floor(1+n*rand);
p2=floor(1+n*rand);
end
tmp_city=city;
%%交换元素
tmp=tmp_city(p1);
tmp_city(p1)=tmp_city(p2);
tmp_city(p2)=tmp;
%%%%%%%%%%%%%%%计算新路线总距离%%%%%%%%%
len2=func5(tmp_city,n);
%%%%%%%%%%%%%%%新老距离的差值，相当于能量%%%%%%%%%
delta_e=len2-len1;
%%%%%%%%%%%%%%%新路线好于旧路线，用新路线替代旧路线%%%%%%%%%
if delta_e<0
city=tmp_city;
else
%%%%%%%%%%%%%%%以一定概率选择是否接受%%%%%%%%%
if exp(-delta_e/T)>rand()
city=tmp_city;
end
end
end
l=l+1; %迭代次数加1
%%%%%%%%%%%%%%%计算新路线的距离%%%%%%%%%
len(l)=func5(city,n);
%%%%%%%%%%%%%%%温度不断下降%%%%%%%%%
T=T*K;
for i=1:n-1
plot([city(i).x,city(i+1).x],[city(i).y,city(i+1).y],'bo-');
hold on;
end
plot([city(n).x,city(1).x],[city(n).y,city(1).y],'ro-');
title(['优化最短距离:',num2str(len(l))]);
hold off;
pause(0.005);
end
figure(2)
plot(len);
xlabel('迭代次数');
ylabel('目标函数值');
title('适应度进化曲线');
%计算距离的函数
function len=func5(city,n)
len=0;
for i=1:n-1
len=len+sqrt((city(i).x-city(i+1).x)^2+(city(i).y-city(i+1).y)^2);
end
len=len+sqrt((city(n).x-city(1).x)^2+(city(n).y-city(1).y)^2);
end
%=========================================================================================