模拟退火算法求解TSP问题

    %time:2016/8/31
    %author:zsou
    %模拟退火算法求解tsp问题
    clear;
    clc;
    [A,X,Y]=textread('berllin52.txt','%d%f%f',52,'headerlines',0);
    %距离矩阵52*52（任意两个节点可达）
    amount=52;
    x_temp1=X*ones(1,amount);
    x_temp2=x_temp1';
    y_temp1=Y*ones(1,amount);
    y_temp2=y_temp1';
    dist=sqrt((x_temp1-x_temp2).^2+(y_temp1-y_temp2).^2);
    
    %模拟退火参数
    t0=97;%初始温度
    tf=3;%终止温度
    markov_length=1000;%markov链长度
    a=0.99%温度衰减系数
    
    %生成初始值
    sol_new=1:amount;     %产生初始解
    sol_current=sol_new;   %当前解
    sol_best=sol_new;    
    E_current=inf;        %当前解对应的回路路径
    E_best=inf;
    
    t=t0;
    while t>tf
    for r=1:markov_length
        %产生随机扰动
        %等概率随机进行决定两次还是三次交换
        if rand<0.5
            %进行两交换
            id1=0;
            id2=0;
            while id1==id2
                id1=ceil(rand*amount);
                id2=ceil(rand*amount);
            end
            temp=sol_new(id1);
            sol_new(id1)=sol_new(id2);
            sol_new(id2)=temp;
        else
            %进行三交换
            id1=0;
            id2=0;
            id3=0;
            while id1==id2|id2==id3|id1==id3|abs(id1-id2)==1
                temp=[ceil(rand*amount) ceil(rand*amount) ceil(rand*amount)];
                temp=sort(temp);
                id1=temp(1);
                id2=temp(2);
                id3=temp(3);
            end
            templist=sol_new((id1+1):(id2-1));
            sol_new((id1+1):(id3-id2+id1+1))=sol_new((id2):(id3));
            sol_new((id3-id2+id1+2):(id3))=templist;
        end
    
        %检查是否满足约束
        %此处无与约束
    
        %计算目标函数值
        E_new=0;
        for i=1:amount-1
            E_new=E_new+dist(sol_new(i),sol_new(i+1));
        end
    
        if E_new<E_current
            E_current=E_new;
            sol_current=sol_new;
            if E_new<E_best
                E_best=E_new;
                sol_best=sol_new;
            end
        else
            %新解目标函数值大于当前解以一定概率接受
            if rand<exp(-(E_new-E_current)/t)
                E_current=E_new;
                sol_current=sol_new;
            else
                %未能接受新解时，则丢弃新解，将当前解作为下一次迭代的初始值
                sol_new=sol_current;
            end
        end
    end
    t=t*a;
    end
    disp('最优解为')
    disp(sol_best)
    disp('最短距离')
    disp(E_best)