模拟退火(simulated annealing,SA)算法解决TSP问题

标签：Metropolis end %% S2 S1 算法 annealing 模拟退火 simulated

模拟退火(simulated annealing,SA)算法

该算法的思想最早是由 Metropolis 等提出的。其出发点是基于物理中固体物质的退火过程与一般的组合优化问题之间的相似性。模拟退火法是一种通用的优化算法,其物理退火过程由以下三部分组成:

(1)加温过程：其目的是增强粒子的热运动,使其偏离平衡位置。当温度足够高时,固体将熔为液体,从而消除系统原先存在的非均匀状态。

(2)等温过程。对于与周围环境交换热量而温度不变的封闭系统,系统状态的自发变化总是朝自由能减少的方向进行的,当自由能达到最小时,系统达到平衡状态。

(3)冷却过程。使粒子热运动减弱,系统能量下降,得到晶体结构。

其中,加温过程对应算法的设定初温,等温过程对应算法的Metropolis抽样过程,冷却过程对应控制参数的下降。这里能量的变化就是目标函数,要得到的最优解就是能量最低态Metropolis准则是 SA算法收敛于全局最优解的关键所在,Metropolis准则以一定的概率接受恶化解,这样就使算法跳离局部最优的陷阱。
模拟退火算法为求解传统方法难处理的TSP问题提供了一个有效的途径和通用框架,并逐渐发展成一种迭代自适应启发式概率性搜索算法。模拟退火算法可以用以求解不同的非线性问题,对不可微甚至不连续的函数优化,能以较大概率求得全局优化解,该算法还具有较强的鲁棒性、全局收敛性、隐含并行性及广泛的适应性,并且能处理不同类型的优化设计变量(离散的、连续的和混合型的),不需要任何的辅助信息,对目标函数和约束函数没有任何要求。利用Metropolis算法并适当地控制温度下降过程,在优化问题中具有很强的竞争力,本文研究基于模拟退火算法的TSP算法。

SA 算法实现过程如下(以最小化问题为例)

(1) 初始化：取初始温度 $T_0$ 。足够大，令 $T = T_0$ ,任取初始解 $S_1$ ,确定每个 $T$ 时的迭代次数,即 Metropolis 链长 $L$

(2)对当前温度 $T$ 和 $k = 1,2,3...,L$ 重复步骤(3)~(6)。

(3)对当前解 $S_1$ ,随机扰动产生一个新解 $S_2$ 。
(4)计算 $S_2$ ,的增量 $df =f(S_2)-f(S_1)$ 其中f $f(S_1)$ 为 $S_1$ 的代价函数。

(5)若 $df<0$ ,则接受 $S_2$ 作为新的当前解,即 $S_1=S_2$ ;否则计算 $S_2$ 的接受概率 $exp(-df/T)$ ，即随机产生(0,1)区间上均匀分布的随机数rand,若 $exp(-df/T)$ >rand,也接受 $S_2$ 作为新的当前解, $S_1=S_2$ ;否则保留当前解 $S_1$ 。

(6)如果满足终止条件Stop,则输出当前解S,为最优解,结束程序。终止条件Stop通常为:在连续若干个Metropolis链中新解S,都没有被接受时终止算法,或是设定结束温度。否则按衰减函数衰减T后返回步骤(2)。

以上步骤称为 Metropolis 过程。逐渐降低控制温度,重复Metropolis过程,直至满足结束准则Stop,求出最优解。

TSP问题见利用混合粒子群算法的解决TSP问题一文

https://blog.csdn.net/pgpaojiao/article/details/137150051

详细步骤

1.算法流程：

2. 算法实现

(1)控制参数的设置

需要设置的主要控制参数有降温速率 $q$ 、初始温度 $T_0$ 、结束温度 $T_m$ 以及链长 $L$ 。

（2）初始解

对于n个城市的TSP问题,得到的解就是对1~n的一个排序,其中每个数字为对应城市的编号,如对10个城市的TSP问题{1,2,3,4,5,6,7,8,9,10),则|1|10|2|14|5|6|8|7|9|3 就是一个合法的解,采用产生随机排列的方法产生一个初始解S。

（3）解变换生成新的解

通过对当前解 $S_1$ ,进行变换,产生新的路径数组即新解,这里采用的变换是产生随机数的方法来产生将要交换的两个城市,用二邻域变换法产生新的路径,即新的可行解 $S_2$ 。

例如n=10时,产生两个[1,10]范围内的随机整数 $r_1$ 和 $r_2$ 确定两个位置,将其对换位置,如 $r_1$ =4, $r_2$ =7

得到的新的解

（4）Metropolis准则

若路径长度函数为 $f(S)$ 则当前解的路径为 $f(S_1)$ ,新解的路径为 $f(S_2)$ ,路径差为 $df = f(S_2)-f(S_1)$ ，Metropolis准则为

$\left.P=\left\{\begin{matrix}1,&\mathrm{d}f<0\\\exp\left(-\frac{\mathrm{d}f}{T}\right),&\mathrm{d}f\geqslant0\end{matrix}\right.\right.$

如果 $df<0$ ,则以概率1接受新的路径;否则以概率 $exp(df/T)$ 接受新的路径。

（5）降温

利用降温速率 $q$ 进行降温,即 $T=qT$ ,若 $T$ 小于结束温度,则停止迭代输出当前状态,否则继续迭代。

Matlab程序如下

主函数.m

clc;
clear;
close all;
%%
tic
T0=1000;   % 初始温度
Tend=1e-3;  % 终止温度
L=500;    % 各温度下的迭代次数（链长）
q=0.9;    %降温速率

%% 加载数据
load CityPosition1;
%%
D=Distanse(X);  %计算距离矩阵
N=size(D,1);    %城市的个数
%% 初始解
S1=randperm(N);  %随机产生一个初始路线

%% 画出随机解的路径图
DrawPath(S1,X)
pause(0.0001)
%% 输出随机解的路径和总距离
disp('初始种群中的一个随机值:')
OutputPath(S1);
Rlength=PathLength(D,S1);
disp(['总距离：',num2str(Rlength)]);

%% 计算迭代的次数Time
% Time=ceil(double(solve(['1000*(0.9)^x=',num2str(Tend)])));
syms x;
eq = 1000*(0.9)^x == num2str(Tend);
Time=ceil(double(solve(eq,x)));
count=0;        %迭代计数
Obj=zeros(Time,1);         %目标值矩阵初始化
track=zeros(Time,N);       %每代的最优路线矩阵初始化
%% 迭代
while T0>Tend
    count=count+1;     %更新迭代次数
    temp=zeros(L,N+1);
    for k=1:L
        %% 产生新解
        S2=NewAnswer(S1);
        %% Metropolis法则判断是否接受新解
        [S1,R]=Metropolis(S1,S2,D,T0);  %Metropolis 抽样算法
        temp(k,:)=[S1 R];          %记录下一路线的及其路程
    end
    %% 记录每次迭代过程的最优路线
    [d0,index]=min(temp(:,end)); %找出当前温度下最优路线
    if count==1 || d0<Obj(count-1)
        Obj(count)=d0;           %如果当前温度下最优路程小于上一路程则记录当前路程
    else
        Obj(count)=Obj(count-1);%如果当前温度下最优路程大于上一路程则记录上一路程
    end
    track(count,:)=temp(index,1:end-1);  %记录当前温度的最优路线
    T0=q*T0;     %降温
    fprintf(1,'%d\n',count)  %输出当前迭代次数
end
%% 优化过程迭代图
figure
plot(1:count,Obj)
xlabel('迭代次数')
ylabel('距离')
title('优化过程')

%% 最优解的路径图
DrawPath(track(end,:),X)

%% 输出最优解的路线和总距离
disp('最优解:')
S=track(end,:);
p=OutputPath(S);
disp(['总距离：',num2str(PathLength(D,S))]);
disp('-------------------------------------------------------------')
toc

计算距离矩阵.m

function D=Distanse(a)
%% 计算两两城市之间的距离
%输入 a  各城市的位置坐标
%输出 D  两两城市之间的距离
row=size(a,1);
D=zeros(row,row);
for i=1:row
    for j=i+1:row
        D(i,j)=((a(i,1)-a(j,1))^2+(a(i,2)-a(j,2))^2)^0.5;
        D(j,i)=D(i,j);
    end
end

画路径函数.m

function DrawPath(Chrom,X)
%% 画路径函数
%输入
% Chrom  待画路径   
% X      各城市坐标位置
R=[Chrom(1,:) Chrom(1,1)]; %一个随机解(个体)
figure;
hold on
plot(X(:,1),X(:,2),'o','color',[0.5,0.5,0.5])
plot(X(Chrom(1,1),1),X(Chrom(1,1),2),'rv','MarkerSize',20)
for i=1:size(X,1)
    text(X(i,1)+0.05,X(i,2)+0.05,num2str(i),'color',[1,0,0]);
end
A=X(R,:);
row=size(A,1);
for i=2:row
    [arrowx,arrowy] = dsxy2figxy(gca,A(i-1:i,1),A(i-1:i,2));%坐标转换
    annotation('textarrow',arrowx,arrowy,'HeadWidth',8,'color',[0,0,1]);
end
hold off
xlabel('横坐标')
ylabel('纵坐标')
title('轨迹图')
box on

生成新的解.m

function S2=NewAnswer(S1)
%% 输入
% S1:当前解
%% 输出
% S2：新解
N=length(S1);
S2=S1;                
a=round(rand(1,2)*(N-1)+1); %产生两个随机位置 用来交换
W=S2(a(1));
S2(a(1))=S2(a(2));
S2(a(2))=W;         %得到一个新路线

Metropolis准则函数.m

Mefunction [S,R]=Metropolis(S1,S2,D,T)
%% 输入
% S1：  当前解
% S2:   新解
% D:    距离矩阵（两两城市的之间的距离）
% T:    当前温度
%% 输出
% S：   下一个当前解
% R：   下一个当前解的路线距离
%%
R1=PathLength(D,S1);  %计算路线长度
N=length(S1);         %得到城市的个数

R2=PathLength(D,S2);  %计算路线长度
dC=R2-R1;   %计算能力之差
if dC<0       %如果能力降低 接受新路线
    S=S2;
    R=R2;
elseif exp(-dC/T)>=rand   %以exp(-dC/T)概率接受新路线
    S=S2;
    R=R2;
else        %不接受新路线
    S=S1;
    R=R1;
end

路线轨迹图.m

function DrawPath(Chrom,X)
%% 画路径函数
%输入
% Chrom  待画路径   
% X      各城市坐标位置
R=[Chrom(1,:) Chrom(1,1)]; %一个随机解(个体)
figure;
hold on
plot(X(:,1),X(:,2),'o','color',[0.5,0.5,0.5])
plot(X(Chrom(1,1),1),X(Chrom(1,1),2),'rv','MarkerSize',20)
for i=1:size(X,1)
    text(X(i,1)+0.05,X(i,2)+0.05,num2str(i),'color',[1,0,0]);
end
A=X(R,:);
row=size(A,1);
for i=2:row
    [arrowx,arrowy] = dsxy2figxy(gca,A(i-1:i,1),A(i-1:i,2));%坐标转换
    annotation('textarrow',arrowx,arrowy,'HeadWidth',8,'color',[0,0,1]);
end
hold off
xlabel('横坐标')
ylabel('纵坐标')
title('轨迹图')
box on

输出路径函数.m

function p=OutputPath(R)
%% 输出路径函数
%输入：R 路径
R=[R,R(1)];
N=length(R);
p=num2str(R(1));
for i=2:N
    p=[p,'—>',num2str(R(i))];
end
disp(p)

坐标转换.m

function varargout = dsxy2figxy(varargin)
if length(varargin{1}) == 1 && ishandle(varargin{1}) ...
                            && strcmp(get(varargin{1},'type'),'axes')   
    hAx = varargin{1};
    varargin = varargin(2:end);
else
    hAx = gca;
end;
if length(varargin) == 1
    pos = varargin{1};
else
    [x,y] = deal(varargin{:});
end
axun = get(hAx,'Units');
set(hAx,'Units','normalized'); 
axpos = get(hAx,'Position');
axlim = axis(hAx);
axwidth = diff(axlim(1:2));
axheight = diff(axlim(3:4));
if exist('x','var')
    varargout{1} = (x - axlim(1)) * axpos(3) / axwidth + axpos(1);
    varargout{2} = (y - axlim(3)) * axpos(4) / axheight + axpos(2);
else
    pos(1) = (pos(1) - axlim(1)) / axwidth * axpos(3) + axpos(1);
    pos(2) = (pos(2) - axlim(3)) / axheight * axpos(4) + axpos(2);
    pos(3) = pos(3) * axpos(3) / axwidth;
    pos(4) = pos(4) * axpos(4 )/ axheight;
    varargout{1} = pos;
end
set(hAx,'Units',axun)

结果

初始种群中的一个随机值:
7—>1—>8—>2—>12—>6—>14—>13—>3—>11—>5—>10—>4—>9—>7
总距离：72.9482

最优解:
11—>9—>10—>1—>2—>14—>3—>4—>5—>6—>12—>7—>13—>8—>11
总距离：29.3405
-------------------------------------------------------------
历时 1.392324 秒。

标签：Metropolis,end,%%,S2,S1,算法,annealing,模拟退火,simulated
From： https://blog.csdn.net/pgpaojiao/article/details/137202353

模拟退火(simulated annealing,SA)算法解决TSP问题