数学建模模拟退化算法TSP

模拟退火的计算过程分为

1. 对于模型求解最小值或者最大值，首先得建立一个函数关系式，然后在定义域范围内随机生成一个初始解x0
2. 在初始解附近生成另外一个解x_new,并保证新解也要在定义域内部
3. 确定什么情况下新解可以替代旧解
4. 重复上述过程N次，得到最优解

单元TSP问题描述

一个旅行商必须访问n个城市，这个城市是个完全图，旅行商需要恰好访问所有城市一次，并且回到最终的城市，城市与城市之间有一个旅行费用，旅行商希望总费用最小

基本概念 无向图：任意两个顶点构成的偶对是无序的 完全图：任意两个顶点之间都有方向互为相反的两条弧边相连接 图的表示：邻接矩阵adjacency matrix 邻接表adjacency list

下面利用算法计算来自数据共享网站的TSP求解过程，可视化显示未作出来

The Traveling Salesman Problemwww.math.uwaterloo.ca/tsp/index.html

%%准备1，输入数据
v=[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(v,1);%%计算需要走的城市个数

d=zeros(n);
for i=2:n
    for j=1:i
        x_i=v(i,1);
        y_i=v(i,2);
        x_j=v(j,1);
        y_j=v(j,2);
        d(i,j)=sqrt((x_i-x_j)^2+(y_i-y_j)^2);
    end
end
d=d+d';

%%准备2，设置基本参数
T0=1000;
T=T0;
N=1000;
Nn=500;
a=0.95;

• 第一步，生成随机解

%%第一步，随机生成可行解
path0=randperm(n);
iter_path=path0;
iter_result=tsp_result(path0,d);

• 由于路径距离的计算之后多次用到，于是可以另写成一个tsp_result.m的文件，用来专门计算距离

function  result = tsp_result(path,d)
% 输入：path:路径（1至n的一个序列），d：距离矩阵
    n=length(path);
    result = 0; % 初始化该路径走的距离为0
    for i=1:n-1  
        result = d(path(i),path(i+1)) + result;  % 按照这个序列不断的更新走过的路程这个值
    end   
    result = d(path(1),path(n)) + result;  % 别忘了加上从最后一个城市返回到最开始那个城市的距离
end

• 第二步，生成新解

%%第二步，生成新的解
for i=1:N
    for j=1:Nn
       result0=tsp_result(path0,d);%计算当前路径总长
       path1=tsp_newpath(path0);%利用方法生成新的路径
       result1=tsp_result(path1,d);

• 旅行商问题生成新解的方式和之前求解函数问题不一样，因为TSP的解可以看作是39个点的顺序，因此可以用“交换法”或者“移位法”或者“倒置法”，在实际计算过程中可以综合三种方法

function path1 = tsp_newpath(path0)
    % path0: 原来的路径
    n = length(path0);
    % 随机选择三种产生新路径的方法
    p1 = 0.33;  % 使用交换法产生新路径的概率
    p2 = 0.33;  % 使用移位法产生新路径的概率
    r = rand(1); % 随机生成一个[0 1]间均匀分布的随机数
    if  r< p1    % 使用交换法产生新路径 
        c1 = randi(n);   % 生成1-n中的一个随机整数
        c2 = randi(n);   % 生成1-n中的一个随机整数
        path1 = path0;  % 将path0的值赋给path1
        path1(c1) = path0(c2);  %改变path1第c1个位置的元素为path0第c2个位置的元素
        path1(c2) = path0(c1);  %改变path1第c2个位置的元素为path0第c1个位置的元素
    elseif r < p1+p2 % 使用移位法产生新路径
        c1 = randi(n);   % 生成1-n中的一个随机整数
        c2 = randi(n);   % 生成1-n中的一个随机整数
        c3 = randi(n);   % 生成1-n中的一个随机整数
        sort_c = sort([c1 c2 c3]);  % 对c1 c2 c3从小到大排序
        c1 = sort_c(1);  c2 = sort_c(2);  c3 = sort_c(3);  % c1 < = c2 <=  c3
        tem1 = path0(1:c1-1);
        tem2 = path0(c1:c2);
        tem3 = path0(c2+1:c3);
        tem4 = path0(c3+1:end);
        path1 = [tem1 tem3 tem2 tem4];
    else  % 使用倒置法产生新路径
        c1 = randi(n);   % 生成1-n中的一个随机整数
        c2 = randi(n);   % 生成1-n中的一个随机整数
        if c1>c2  % 如果c1比c2大，就交换c1和c2的值
            tem = c2;
            c2 = c1;
            c1 = tem;
        end
        tem1 = path0(1:c1-1);
        tem2 = path0(c1:c2);
        tem3 = path0(c2+1:end);
        path1 = [tem1 fliplr(tem2) tem3];   %矩阵的左右对称翻转 fliplr，上下对称翻转 flipud
    end
end

• 第三步，确定是否交换

%%第三步，是否交换
 if result1<result0
           path0=path1;
           iter_path=[iter_path;path1];
           iter_result=[iter_result;result1];
       else
           p=exp(-abs(result1-result0)/T);
           if rand(1)<p
               path0=path1;
               iter_path=[iter_path;path1];
                iter_result=[iter_result;result1];
           end
       end
    end
    T=a*T;
end

• 第四步，输出结果

%%第四步，显示结果
[best_result,ind]=min(iter_result);
min_path=iter_path(ind,:);
disp('最佳的方案是：'); disp(min_path)
disp('此时最优值是：'); disp(best_result)

添加图片注释，不超过 140 字（可选）

新建一个figure画出散点图可能效果会好一点