模拟退火（Simulated Annealing,SA）算法解决3维TSP问题

~容~若~

已于 2024-07-24 09:33:52 修改

阅读量2k

点赞数 26

文章标签：算法模拟退火算法机器学习

于 2024-07-24 09:14:08 首次发布

本文链接：https://blog.csdn.net/weixin_47616070/article/details/140652878

版权

模拟退火算法（Simulated Annealing，SA）是一种基于物理退火过程的随机优化算法，特别适用于全局优化问题。其核心思想来源于固体物理中的退火过程，即通过逐步降低温度来达到系统的最低能量状态。以下是模拟退火算法的主要特点和步骤：

主要特点

全局优化：模拟退火算法能够跳出局部最优，逐步寻找全局最优解。
随机性：算法通过随机搜索接受新解，可以避免陷入局部最优。
温度控制：温度逐步降低，控制搜索过程的随机性，最终收敛到一个稳定状态。

算法步骤

初始化：选择一个初始解 $(S)$ 和初始温度 $(T)$ 。
生成新解：在当前解 $(S)$ 的邻域内随机选择一个新解 $(S^{'})$ 。
评价新解：计算新解 $(S^{'})$ 和当前解 $(S)$ 的目标函数值 $(f (S))$ 和 $(f (S^{'}))$ 。
接受新解：根据以下规则决定是否接受新解 $(S^{'})$ ：
- 若 $(f (S^{'}) < f (S))$ ，则接受新解。
- 若 $\ge f(S) )$ ，则以概率 $\exp\left(\frac{-(f(S') - f(S))}{T}\right) )$ 接受新解。
更新温度：按某种方式降低温度 $(T)$ ，例如 $\alpha T )$ ，其中 $(\alpha )$ 是一个小于 1 的常数。
终止条件：当温度 $(T)$ 降到某个阈值或者满足其他终止条件时，停止算法。

应用领域

模拟退火算法广泛应用于组合优化问题，如旅行商问题（TSP）、背包问题、图像处理、生产调度等领域。

通过适当调整初始温度、降温策略和邻域搜索策略，模拟退火算法可以在不同的优化问题中表现出较好的性能。

代码通过逐步降低温度来优化路径，并使用随机扰动来探索可能的路径。以下是代码的结构和主要过程：

代码结构

数据加载和初始化
- 读取加工区域的坐标数据。
- 计算任意两个加工区域之间的距离矩阵。
- 设置初始温度、马尔可夫链长度和衰减参数。
- 初始化城市坐标结构体。
模拟退火算法主循环
- 在温度高于某个阈值时重复迭代：
  - 在每个温度下进行多次扰动实验。
  - 计算当前路径的总距离。
  - 生成随机扰动，即随机交换两个城市的位置。
  - 计算新路径的总距离。
  - 根据能量差值决定是否接受新路径。
  - 更新温度。
  - 绘制当前路径图。
  - 记录路径长度变化。
辅助函数
- 计算路径长度的函数 calculate_length。

主要过程

数据加载和距离矩阵计算

file_path = '..\chapter2\data.csv';
Location = csvread(file_path, 1, 1);
n = size(Location, 1);  % TSP问题的规模
D = zeros(n);  % 任意两个加工区域的距离间隔矩阵

for i = 1:n
    for j = 1:n
        D(i, j) = sqrt((Location(i, 1) - Location(j, 1))^2 + (Location(i, 2) - Location(j, 2))^2 + (Location(i, 3) - Location(j, 3)));
    end
end

初始参数设置

T = 100 * n;  % 初始温度
L = 100;  % 马尔可夫链长度
K = 0.99;  % 衰减参数

初始化城市坐标结构体

city = struct([]);
for i = 1:n
    city(i).x = Location(i, 1);
    city(i).y = Location(i, 2);
    city(i).z = Location(i, 3);
end

模拟退火算法主循环

l = 1;  % 统计迭代次数
figure(1);

while T > 0.001  % 停止迭代温度
    for i = 1:L
        len1 = calculate_length(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 = calculate_length(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;
    len(l) = calculate_length(city, n);
    T = T * K;
    
    for i = 1:n-1
        plot3([city(i).x, city(i+1).x], [city(i).y, city(i+1).y], [city(i).z, city(i+1).z], 'bo-');
        hold on;
    end
    plot3([city(n).x, city(1).x], [city(n).y, city(1).y], [city(n).z, city(1).z], 'ro-');
    title(['优化最短距离:', num2str(len(l))]);
    hold off;
    pause(0.005);
end

figure(2);
plot(len)
xlabel('迭代次数')
ylabel('目标函数值')
title('适应度进化曲线')

计算路径长度的函数

function len = calculate_length(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 + (city(i).z - city(i+1).z)^2);
    end
    len = len + sqrt((city(n).x - city(1).x)^2 + (city(n).y - city(1).y)^2 + (city(n).z - city(1).z)^2);
end

结果

在这里插入图片描述

完整代码

clc;
clear all;
%%%%%%%%%绘制加工区域中心图%%%%%%%%%%%%%%%%%
file_path = '..\chapter2\data.csv';
Location = csvread(file_path, 1, 1);
n=size(Location,1);                    %TSP问题的规模
D=zeros(n);                     %任意两个加工区域的距离间隔矩阵
%%%%%%%%%%%%%%%%%%%%%求任意两个城市距离间隔矩阵%%%%%%%%%%%%%%%%%%%%%
for i=1:n
    for j=1:n
        D(i,j)=((Location(i,1)-Location(j,1))^2+(Location(i,2)-Location(j,2))^2+(Location(i,3)-Location(j,3)))^0.5;
    end
end
T=100*n;                         %初始温度
L=100;                           %马可夫链长度
K=0.99;                          %衰减参数
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%城市坐标结构体%%%%%%%%%%%%%%%%%%%%%%%%%%
%S0=[22	78	122	101	49	131	80	21	100	97	103	99	96	37	12	93	31	102	98	115	41	44	19	20	109	107	111	112	64	35	46	42	75	113	16	71	74	72	69	70	73	68	67	63	29	76	120	50	9	105	62	130	126	124	128	129	125	127	123	84	34	33	66	60	116	59	57	53	54	58	56	52	13	94	25	121	117	26	11	132	85	87	91	88	89	92	86	90	28	15	36	38	17	48	134	39	2	114	61	135	82	95	10	65	77	1	5	3	8	7	4	6	55	79	24	51	104	133	14	23	18	27	106	81	110	108	30	118	45	47	40	136	83	32	119	43];%ACO结果
city=struct([]);
% for i=1:n
%     city(i).x=Location(S0(i),1);
%     city(i).y=Location(S0(i),2); 
%     city(i).z=Location(S0(i),3);
% end
for i=1:n
    city(i).x=Location(i,1);
    city(i).y=Location(i,2);
    city(i).z=Location(i,3);
end
l=1;                             %统计迭代次数
len(l)=calculate_length(city,n);
figure(1);

while T>0.001                    %停止迭代温度
    %%%%%%%%%%%%%%%%多次迭代扰动，温度降低之前多次实验%%%%%%%%%%%%%%%
    for i=1:L            
        %%%%%%%%%%%%%%%%%%%计算原路线总距离%%%%%%%%%%%%%%%%%%%%%%%%%
        len1=calculate_length(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=calculate_length (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;
    %%%%%%%%%%%%%%%%%%%%%%%%%计算新路线距离%%%%%%%%%%%%%%%%%%%%%%%%%%
    len(l)=calculate_length(city,n); 
    %%%%%%%%%%%%%%%%%%%%%%%%%%%温度不断下降%%%%%%%%%%%%%%%%%%%%%%%%%%
    T=T*K;   
    for i=1:n-1
        plot3([city(i).x,city(i+1).x],[city(i).y,city(i+1).y],[city(i).z,city(i+1).z],'bo-');
        hold on;
    end
    plot3([city(n).x,city(1).x],[city(n).y,city(1).y],[city(n).z,city(1).z],'ro-');
    title(['优化最短距离:',num2str(len(l))]);
    hold off;
    pause(0.005);
end
figure(2);
plot(len)
xlabel('迭代次数')
ylabel('目标函数值')
title('适应度进化曲线')










% 计算路径长度的函数
function len = calculate_length(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 + (city(i).z - city(i+1).z)^2);
    end
    len = len + sqrt((city(n).x - city(1).x)^2 + (city(n).y - city(1).y)^2 + (city(n).z - city(1).z)^2);
end