Matlab-退火算法

一、模拟退火概念

通俗来说模拟退火就是要走出舒适圈，多去试一试，万一成了,它是一个不断试错的过程。

舒适圈:当前处于局部最优状态。

试一试:随机试探新解，如果新解更好就选择新解，否则就以一定的概率选择新解。

万一成了:找到了更优解甚至最优解。

二、出现原因

由于贪心算法会陷入局部最优状态。

三、基本思想

爬山为例:

1. 在可见的范围内随机选择一点。
2. (1)如果该点比当前位置更高就会直接去该点。(2)如果比当前位置更低则就会以一定的概率来决定去不去该点。
3. 因此在当前的局部最优解的山峰上会有一定的概率继续往下走，从来发现更高的山峰，也就是更好的解。
   

关键点

1. 时间有限，需要设置终止条件，也就是退火的过程。
2. 算法需要做到，如果当前山峰越高，前往低点的概率越低。
   关键是概率怎么求

四、例题

在这个例题中起点和终点固定，因此可行解一共33!种情况，我们要做的就是找出最优解，也就是路径之和最少的一种。

##求解思路
(1) 设置一个初始问题T0和一个初始解x(0),对于初始解可以把原始数据排列直接作为初始解，也可以使用蒙特卡罗法求一个较好的初始解。
(2) 在当前温度下随机产生一个新解x`,其产生的规则不唯一，但是要随机性。

(3)求出新解和当前解的差值。
由这个新解 X和当前解的联系可以看出其差值就是:

(4)如果△f<0,说明新解比旧解的路程之和短，因此概率为1的接收此解。如果>=0,就以一定的概率接收，公式定义如下:


对于exp(-△f/Ti),当△f越小，意味着-△f/Ti越大，也就是exp(-△f/Ti)越大，它>=rand(0,1)的概率也就越大。
(5)不断重复上述过程，自行设置执行次数。
(6)退火过程:选定降温系数α,Ti+1=αTi；

算法实现

%读取exel
citys=readtable("citys.xlsx",'Range','b2:c35');
%转为double
city=table2array(citys);
%获取城市个数
n=size(city,1);  %1为矩阵行数 2为矩阵列数
%所有结点,终点看成第35的结点。
d=zeros(n+1,n+1);
%获取结点和结点之间的距离
for i=1:n
 for   j=1:n
     d(i,j)=distance(city(i,2),city(i,1),city(j,2),city(j,1),6371);%6371为地球半径
end 
end
%蒙特卡罗法求一个比较好的初始解。
path=[];%路径
length=inf; %inf表示正无穷
for i=1:5000
    %随机路径
    temp_path=[1,1+randperm(33),35];%起点和终点不变，2-34随机打乱顺序生成
    temp_length=0;
    %求随机路径总长度
    for j=1:34
        temp_length=temp_length+d(temp_path(j),temp_path(j+1));
    end
    if(temp_length<length)
        path=temp_path;
        length=temp_length;
    end

end
%模拟退火部分
%初始温度
T=1; %随便设置
%临界点
e=10^(-30);
%退火系数
alpha=0.999;
while T>e
    %随机生成uv
    %floor向下取整
    c=2+floor(33*rand(1,2)); %rand(1,2)表示随机生成2个数,33*rand(1,2)->[0,33)->[0,32]->[2,34]
    c=sort(c);
    u=c(1);
    v=c(2);
    %获取两个不同路径之间的差值
     df=d(path(u-1),path(v))+d(path(u),path(v+1))...
        -d(path(u-1),path(u))+d(path(v),path(v+1));
    if df<0
        %更新解,新解是uv之间的逆序
        path=[path(1:u-1),path(v:-1:u),path(v+1:35)];
        length=length+df;

   
   elseif exp(-df/T)>=rand
    path=[path(1:u-1),path(v:-1:u),path(v+1:35)];
        length=length+df;
    else
    end
    %更新温度
    T=T*alpha;
end
%画图
path(35)=1;     % 1号和35号都是北京，把35改成1，方便画图
%city(path ,1)，也就是将path中的结点对应的经度取出来
plot(city(path ,1), city(path ,2),'o-');
for i = 1:n
    %对每个城市进行标号
    text(city(i,1),city(i,2),['   ' num2str(i)]); %number to string
end
xlabel('东经')
ylabel('北纬')

路径图