利用模拟退火算法对TSP(旅行商问题)的实现，含matlab代码
前言：
为什么可以通过模拟退火算法去解决旅行商问题？
首先，模拟退火算法脱胎于物理退火过程，其通过设定初始温度，在温度不断地降低的过程中，通过某种手段去解空间中不断地寻找最优解（解空间：目标函数的离散化点集，最优解：极值）。
其次，在同等温度下，其也会进行多次改变并确定最优解。
实现的基本步骤：
1、控制参数的设置
需要设置的主要控制参数有降温速率q，最高温度Temper、最低温度1/Temper、以及链长。
2、编码
用整数排列编码的方式；对于n个城市的TSP（旅行商）问题，染色体分为n段，其中每一段对应城市编号，如对10个城市的TSP（旅行商）问题5-3-7-8-2-1-6-4-9-10就是一个合法的染色体。
3、种群初始化
在完成染色体编码后，需对种群进行初始化作为起始解，所以首先需要确定初始化种群的数目。种群中染色体的个数就是种群的数量。一般来说种群数量是一个常数，可以取50~200之间。
4、领域搜索算法
通过对当前解进行变换，产生的路径数组即为新解。此处含有四种领域搜索算子：
1、逆序搜索算子
s1:10-2-3-4-1-5-7-9-8-6,Ti=2,Tj=7
s1’:10-2-9-7-5-1-4-3-8-6
即将s1的染色体中的第Ti城市到第Tj城市之间进行倒序排列
2、1-0p交换搜索算子
s1:10-2-3-4-1-5-7-9-8-6,Ti=2,Tj=7
s1’:10-2-4-1-5-7-9-3-8-6
即将s1的染色体中的第3个数字放到第8个数字后进行排列
3、2-op交换搜索算子
s1:10-2-3-4-1-5-7-9-8-6,Ti=2,Tj=7
s1’:10-2-9-4-1-5-7-3-8-6
即将s1的染色体中的第3个数字与第8个数字进行交换
4、3-op交叉搜索算子
s1:10-2-3-4-1-5-7-9-8-6,Ti=2,Ti+1=3,Tj=5,Tj+1=7
s1’:10-3-2-4-7-5-1-9-8-6
即将Ti城市与Ti+1城市进行交换，将Tj城市与Tj+1城市进行交换
5、Metropolis准则
若路径长度函数为Z(s),则当前解的路径长度为Z(s1),新解路径长度为Z(s1’)，路径长度差为△Z=Z(s1’)-Z(s1),则Metropolis准则为
P=1,deltaZ<0
P=exp(-△Z/T),△>=0
之所以采用这个准则，那是因为在求最小值的时候。如果新解值小于当前解，我们全部采纳。如果新解大于当前解，我们按照一定的概率（exp(-△Z/T)<1）采纳，目的是为了避免陷入局部解！
6、降温
利用降温速率进行降温，即T=qT，若T小于结束温度，则停止迭代，输出当前状态，否则继续迭代。
实例：
以中国31个省会城市的最短旅行路径为例，其横纵坐标如表所示（左列为横坐标，右列为纵坐标）
2312 1304
1315 3639
2244 4177
1399 3712
1535 3488
1556 3326
1229 3238
1044 4196
790 4312
570 4386
1970 3007
1756 2562
1491 2788
1676 2381
695 1332
1678 3715
2179 3918
2370 4061
2212 3780
2578 3676
2838 4029
2931 4263
1908 3429
2376 3507
2643 3394
3201 3439
3240 2935
3550 3140
2357 2545
2826 2778
2975 2370
模拟退火参数设置见下表：
|降温速度| 初始温度 |结束温度 |链长
| 0.9 | 1000 |0.001 |200
1.主函数如下：
%SA_TSP(load('CityPosition.mat'),50,1000,100,100) 本句在运行代码前输入，以作为变量初始值
%计算最小值时，在同一温度下搜索，优先考虑改善情况，没有改善情况时，接受一定概率变差情况
%CityPosition为城市坐标，P为种群规模，Temper为最高温度
%StopN为恒温搜索次数，Stopn为连续未改变次数，类似禁忌长度
function [ShortestL,ShortestR]=SA_TSP(CityPosition,P,Temper,StopN,Stopn)
%%%%%%%%%%%%%%%%%%%%%种群及各参数初始化%%%%%%%%%%%%%%%%%%%%%%
load('CityPosition.mat'); %注意将上述城市的坐标输入到表格中并保存为CityPosition.mat文件
N=size(CityPosition,1); %计算城市个数，代表城市的表格CityPosition的行数，1代表行数，2代表列数
farm=zeros(P,N); %开辟空间farm，存放初始解，这里在初始化
Farm=zeros(P,N); %开辟空间Farm，存放变化后的新解，这里在初始化
len=zeros(1,P); %开辟空间len，存放初始解路径长度，这里在初始化
Len=zeros(1,P); %开辟空间len，存放变化后的路径长度新解，这里在初始化
Solver=zeros(1e3,N); %开辟Solver产生横1000、纵N的表空间，
Len_solver=inf*ones(1,1e3); %优秀路径长度
D=distance(CityPosition); %将城市的坐标矩阵转换为领接矩阵
for i=1:P
farm(i,:)=randperm(N); %构造初始可行解集合，产生最初旅游线路
end
T=Temper; %最高温度
step=1; %种群进化次数
%%%%%%%%%%%%%%%%%%%%降温的同时计算路径长度、搜索更新路线并比较改善情况%%%%%%%%%%%%%%%%%%%%
while T>=1/T %温度是否降到最低
Search=1; %固定温度下迭代计数器
Stop=1; %解没有改善的计数器
while Stop<Stopn && Search<StopN
for i=1:P %计算一次形变所有的城市的路程
len(i)=myLength(D,farm(i,:)); %计算路径长度
rr=randi(4); %产生一个1到4的随机数
switch (rr)
case 1 %搜索新解更新路线
Farm(i,:)=perform_reverse(farm(i,:));
case 2
Farm(i,:)=perform_1opt(farm(i,:));
case 3
Farm(i,:)=perform_2opt(farm(i,:));
case 4
Farm(i,:)=perform_3opt(farm(i,:));
end
Len(i)=myLength(D,Farm(i,:));
end
z=zeros(1,P); %比较解改善的情况
R=rand(1,P); %变差可接受的概率
if find((Len-len<z | exp((len-Len)/(T))>R)~=0) %保留改善和可以接受的变差解
farm(find((Len-len<z | exp((len-Len)/(T))>R)~=0),:)= Farm(find((Len-len<z | exp((len-Len)/(T))>R)~=0),:);
len(find((Len-len<z | exp((len-Len)/(T))>R)~=0))= Len(find((Len-len<z | exp((len-Len)/(T))>R)~=0));
[Temp_value,Temp_Index]=min(len);
Solver(step,:)=farm(Temp_Index(1),:);
Len_solver(step,:)=Temp_value(1);
step=step+1;
stop=1;
else
stop=stop+1;
end
Search=Search+1;
end
T=T*0.9;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%绘制降温过程中得到的最小路线图%%%%%%%%%%%%%%%%%%%%%%%
[Value,Index]=min(Len_solver);
ShortestL=Value(1);
ShortestR=Solver(Index(1),:);
subplot(1,2,1)
DrawRoute(CityPosition,ShortestR);
title('最优路径')
subplot(1,2,2)
plot(Len_solver)
xlabel('迭代次数')
ylabel('距离')
title('优化过程')
end
2.逆序搜索算子
%随机生成t1,t2，逆序其中间的代码
function route=perform_reverse(R)
N=length(R); %求解初始路径长度
route=R;
while (1)
t=unidrnd(N,1,2); %产生一组小于或等于N的，1行2列的随机整数
if abs(t(1)-t(2))>1
break;
end
end
t1=min(t);t2=max(t);
route(t1:t2)=R(t2:-1:t1); %采用逆序搜索算子，将t1到t2直降的数进行逆序排列
end
- 1-0p交换搜索算子
%输入截取一个代码插入另一个位置的新路径
function route=perform_1opt(R)
N=length(R);
while (1)
t=unidrnd(N,1,2);
if abs(t(1)-t(2))>1
break;
end
end
t1=min(t);t2=max(t);
if t1==1 && t2<N
route=[R(t1+1:t2),R(t1),R(t2+1:N)];
elseif t1==1 && t2==N
route=[R(t1+1:t2),R(t1)];
elseif t1>1 && t2==N
route=[R(1:t1-1),R(t1+1:t2),R(t1)];
else
route=[R(1:t1-1),R(t1+1:t2),R(t1),R(t2+1:N)];
end
end
4.2-op交换搜索算子
%交换选中的两个点的基因码
function route=perform_2opt(R)
N=length(R);
route=R;
while (1)
t=unidrnd(N,1,2);
if abs(t(1)-t(2))>1
break;
end
end
route(t(1))=R(t(2));
route(t(2))=R(t(1));
end
5.3-op搜索算子
%两两交换指定的基因码，如（i,i+1,j,j+1）->(i+1,i,j+1,j)
function route=perform_3opt(R)
N=length(R);
while (1)
t=unidrnd(N,1,2);
if abs(t(1)-t(2))>1
break;
end
end
t1=min(t);t2=max(t);
if t1==1 && t2==N
route=[R(t1),R(t1+2:t2),R(t1+1)];
elseif t1==1 && t2<N
route=[R(t1+1),R(t1),R(t1+2:t2-1),R(t2+1),R(t2),R(t2+2:N)];
elseif t1>1 && t2==N
route=[R(2:t1-1),R(t1+1),R(t1),R(t1+2:N-1),R(1),R(N)];
else
route=[R(1:t1-1),R(t1+1),R(t1),R(t1+2:t2-1),R(t2+1),R(t2),R(t2+2:N)];
end
end
6.计算距离矩阵
利用所给出的N个城市的坐标吗，计算N个城市之间的距离，得到矩阵NxN
%计算领接矩阵，输出无向图的临接矩阵
function D=distance(coordinate)
N=size(coordinate,1);
D=zeros(N,N);
for i=1:N
for j=1:N %采取欧式直角坐标法计算两点间的距离
%注：coordinate(i,1)代表行的第i个数，coordinate(j,1)代表行的第j个数
% coordinate(i,2)代表列的第i个数，coordinate(j,2)代表列的第j个数
bb=(coordinate(i,1)-coordinate(j,1)).^2+(coordinate(i,2)-coordinate(j,2)).^2;
D(i,j)=bb^(0.5);
D(j,i)=D(i,j);
end
end
end
%%%%所谓的邻接矩阵输出无向图，指的是两个地区相互之间的距离%%%%
%%%%求1-2-3的邻接矩阵%%%%%%%
%%%%则1-2-3
%%%%1
%%%%2
%%%%3
%%%%%%%%%%%%%%由上式构成的矩阵
7.绘制路线图
%绘制城市及旅行路线
function DrawRoute(Cityposition,R)
N=length(R);
p=[R(1:end),R(1)];
scatter(Cityposition(:,1),Cityposition(:,2),'o','LineWidth',2,'MarkerEdgeColor','k','MarkerEdgeColor','g');
hold on %描点
for i=1:N %画线
plot([Cityposition(p(i),1),Cityposition(p(i+1),1)],[Cityposition(p(i),2),Cityposition(p(i+1),2)],'ms-','LineWidth',2,'MarkerEdgeColor','k','MarkerEdgeColor','g')
hold on
end
end
8.计算路径长度
%计算路径长度
function len=myLength(D,farm)
len=0;
p=farm(1,:);
for i=1:(length(p)-1)
len=len+D(i,i+1);
end
end
扫一扫
5163
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
