应用GA和PSO算法求解10城市TSP问题.docx
《应用GA和PSO算法求解10城市TSP问题.docx》由会员分享，可在线阅读，更多相关《应用GA和PSO算法求解10城市TSP问题.docx（17页珍藏版）》请在冰豆网上搜索。
应用GA和PSO算法求解10城市TSP问题
应用GA和PSO算法求解10城市TSP问题
1问题描述
旅行团计划近期在城市A、B、C、D、E、F、G、H、I和J共10个城市间进行一次周游旅行，为了尽量节省旅行开支，希望能找到一条里程数最少或相对较少的旅行路线。
问题1，给定10个城市之间的公路里程如表1所示，并要求使用GA算法求解优化问题。
问题2，与问题1数据相同，要求使用PSO算法求解优化问题。
表1城市位置坐标（单位：
km）
横坐标
纵坐标
城市A
40
44.39
城市B
24.39
14.63
城市C
17.07
22.93
城市D
22.93
76.1
城市E
51.71
94.14
城市F
87.32
65.36
城市G
68.78
52.19
城市H
84.88
36.09
城市I
66.83
25.36
城市J
61.95
26.34
2使用GA算法求解
2.1算法描述
（1）编码和适应度函数
分别用1-10表示城市A-J，然后采用自然数编码方式为TSP问题编码，例如，旅程（16289105734）表示从城市A出发分别经过了F-B-H-I-J-E-G-C-D的一次旅行。
每一个问题的解及算法中的个体都可以计算相应的距离。
那么种群中的最小距离和最大距离也相应的可以确定。
选择种群个数为50。
根据种群中个体的距离并考虑使用自适应的标定方法，定义如下的适应度函数，
使用此适应度函数的后面的乘方次数可以调整来改变淘汰速度。
这里选择2，表示淘汰速度比较适中。
（2）定义算子
选择算子，根据适应度函数的大小进行选择。
计算每个个体被选中的概率
，以各个个体所分配到的概率值作为其遗传到下一代的概率，基于这些概率用赌盘选择法来产生下一代群体。
交叉算子，采用部分映射交叉（PartiallyMappedCrossover,PMX）方法实现算法交叉。
首先选取选需要交叉的区间段，然后确定区间段的映射关系，接下来交换区间段的遗传信息，此时未换部分的遗传信息会出现不合法的情况，因此根据将未换部分按映射关系进行交换。
交叉率为0.9。
变异算子，把一个染色体中的两个基因的交换作为变异算法。
在算法中随机的产生一个染色体中需要交换的两个基因的位置，将这两个基因进行交换来实现变异。
变异率为0.2。
（3）算法流程
根据上述的遗传算子的定义，并设置最大的迭代次数为1000，将GA算法流程叙述如下，
（i）随机生成初始种群。
（ii）按照本节
（1）中的公式计算各个个体适应度的值。
（iii）判断是否达到了最大的迭代次数。
若是，算法结束，输出计算结果；若否，转到（iv）。
（iv）按照本节
（2）中的选择算子进行选择操作。
（v）选择交叉宽度并随机确定交叉切点，按照本节
（2）中的交叉算子进行交叉操作。
（vi）按照本节
（2）中的变异算子进行变异操作。
（vii）将遗传算子产生的种群作为新的种群，转到（ii）。
2.2GA算法计算结果
使用Matlab编程实现2.1中的算法，计算结果如下。
图1GA算法过程的距离值变化
图2GA算法求解的10城市TSP问题的结果
最后的结果编码为（89102314567），解为H-I-J-B-C-A-D-E-F-G，距离为269.0671。
从计算结果可以看出，算法迭代到300步时已经收敛到了局部的最优值。
经过大量的计算发现至多迭代到300步，算法收敛到局部最优值。
经过进一步的验证发现，这个局部最优解也是全局最优解。
3使用PSO算法求解
3.1算法描述
（1）TSP问题的交换子和交换序
设n个节点的TSP问题的解序列为s=（ai），i=1…n。
定义交换子SO（i1,i2）为交换解S中的点ai1和ai2，则S’=S+SO（i1,i2）为解S经算子SO（i1,i2）操作后的新解。
一个或多个交换子的有序队列就是交换序，记作SS，SS=（SO1,SO2,…SON）。
其中，SO1,…,SON等是交换子，之间的顺序是有意义的,意味着所有的交换子依次作用于某解上。
若干个交换序可以合并成一个新的交换序，定义
为两个交换序的合并算子。
设两个交换序SS1和SS2按先后顺序作用于解S上，得到新解S’。
假设另外有一个交换序SS’作用于同一解S上，能够得到形同的解S’，可定义
和
属于同一等价集，在交换序等价集中，拥有最少交换子的交换序称为该等价集的基本交换序。
（2）TSP问题的速度和位置更新算式
根据
（1）中的交换子和交换序的定义，可以根据基本的PSO算法速度更新算式和位置更新算式，重新定义如下的速度和位置更新算式，
式中，
、
为[0,1]区间的随机数。
为粒子与个体极值的交换序，以概率
保留；
为粒子与全局极值的交换序，以概率
保留。
粒子的位置按照交换序
进行更新。
为惯性权重。
（3）算法流程
按照本节中的有关定义给出算法流程如下，
（i）初始化粒子群，给每一个粒子一个初始解
和随机的交换序
。
（ii）判断是否达到最大迭代次数1000。
若是，算法结束，输出结果；若不是，转到（iii）。
（iii）根据粒子当前位置计算下一个新解：
（a）计算A=
，A是一个基本交换序，表示A作用于
得到
；
（b）计算B=
，B是一个基本交换序；
（c）按照
（2）中的公式更新速度和位置。
（d）如果得到了更好的个体位置，更新
。
（iv）如果得到了更好的群体位置，更新
。
3.2PSO算法的计算结果
编程实现3.1中的算法，计算结果如下。
计算程序见附录。
图3PSO算法过程的距离值变化
图4PSO算法求解的10城市TSP问题的结果
最后的结果编码为（14567891023），解为A-D-E-F-G-H-I-J-B-C，距离为269.0671。
从计算结果可以看出，算法迭代到200步时已经收敛到了局部的最优值。
这个局部最优解也是全局最优解。
从收敛的速度的平均意义上而言，PSO算法与GA算法比收敛的更快。
附录
%GA算法的代码
%GA算法程序.
data=load（'pos10.txt'）;
a=[data（:
2）data（:
3）];
n=50;%种群数目
C=1000;%迭代次数
Pc=0.9;%交叉概率
Pm=0.2;%变异概率
%GA算法初始化.
[N,NN]=size（D）;
farm=zeros（n,N）;
fori=1:
n
farm（i,:
）=randperm（N）;
end
R=farm（1,:
）;
len=zeros（n,1）;
fitness=zeros（n,1）;
counter=0;
%GA算法迭代
whilecounterfori=1:
n
len（i,1）=myLength（D,farm（i,:
））;
end
maxlen=max（len）;
minlen=min（len）;
fitness=len;
fori=1:
length（len）
fitness（i,1）=（1-（（len（i,1）-minlen）/（maxlen-minlen+0.0001）））.^2;
end;
rr=find（len==minlen）;
R=farm（rr（1,1）,:
）;
FARM=farm;
%计算适应度函数
nn=0;
fori=1:
n
iffitness（i,1）>=rand
nn=nn+1;
FARM（nn,:
）=farm（i,:
）;
end
end
FARM=FARM（1:
nn,:
）;%FARM（nn*N）
[aa,bb]=size（FARM）;%（aa=nn）
%交叉
FARM2=FARM;
fori=1:
2:
aa
ifPc>rand&&iA=FARM（i,:
）;
B=FARM（i+1,:
）;
L=length（a）;
ifL<=10%确定交叉宽度
W=9;
elseif（（L/10）-floor（L/10））>=rand&&L>10
W=ceil（L/10）+8;
else
W=floor（L/10）+8;
end
p=unidrnd（L-W+1）;%随机选择交叉范围
fori=1:
W
x=find（a==b（1,p+i-1））;
y=find（b==a（1,p+i-1））;
[a（1,p+i-1）,b（1,p+i-1）]=exchange（a（1,p+i-1）,b（1,p+i-1））;
[a（1,x）,b（1,y）]=exchange（a（1,x）,b（1,y））;
end
FARM（i,:
）=A;
FARM（i+1,:
）=B;
end
end
%变异
FARM2=FARM;
fori=1:
aa
ifPm>=rand
FARM（i,:
）=mutate（FARM（i,:
））;
end
end
%群体更新
FARM2=zeros（n-aa+1,N）;
ifn-aa>=1
fori=1:
n-aa
FARM2（i,:
）=randperm（N）;
end
end
FARM=[R;FARM;FARM2];
[aa,bb]=size（FARM）;
ifaa>n
FARM=FARM（1:
n,:
）;
end
farm=FARM;
clearFARM
RR（counter+1）=myLength（D,R）
%统计迭代次数。
counter=counter+1;
end
%结果图形显示
figure
holdon
plot（[a（R
（1）,1）,a（R（10）,1）],[a（R
（1）,2）,a（R（10）,2）],'ms-','LineWidth',2,'MarkerEdgeColor','k','MarkerFaceColor','g'）;
scatter（a（:
1）,a（:
2）,'ms'）
gridon
holdon
fori=2:
length（R）
x0=a（R（i-1）,1）;
y0=a（R（i-1）,2）;
x1=a（R（i）,1）;
y1=a（R（i）,2）;
xx=[x0,x1];
yy=[y0,y1];
plot（xx,yy,'LineWidth',4,'MarkerEdgeColor','k','MarkerFaceColor','g'）
holdon
end
%结果输出
R
Rlength
%其他函数
functiona=mutate（a）
L=length（a）;
rray=randperm（L）;
[a（rray
（1））,a（rray
（2））]=exchange（a（rray
（1））,a（rray
（2）））;
function[x,y]=exchange（x,y）
temp=x;
x=y;
y=temp;
functionlen=myLength（D,p）
[N,NN]=size（D）;
len=D（p（1,N）,p（1,1））;
fori=1:
（N-1）
len=len+D（p（1,i）,p（1,i+1））;
end
%PSO算法代码
%%PSO算法代码
data=load（'pos10.txt'）
cityCoor=[data（:
2）data（:
3）];
n=size（cityCoor,1）;
cityDist=zeros（n,n）;
fori=1:
n
forj=1:
n
ifi~=j
cityDist（i,j）=（（cityCoor（i,1）-cityCoor（j,1））^2+...
（cityCoor（i,2）-cityCoor（j,2））^2）^0.5;
end
cityDist（j,i）=cityDist（i,j）;
end
end
individual=zeros（indiNumber,n）;
%初始化
nMax=1000;%迭代次数
indiNumber=50;%粒子个数
%%初始化个体和群体最优值
indiFit=fitness（individual,cityCoor,cityDist）;
[value,index]=min（indiFit）;
tourPbest=individual;
tourGbest=individual（index,:
）;
recordPbest=inf*ones（1,indiNumber）;
recordGbest=indiFit（index）;
xnew1=individual;
%%循环寻找最优路径
L_best=zeros（1,nMax）;
forN=1:
nMax
%更新最优值
indiFit=fitness（individual,cityCoor,cityDist）;
fori=1:
indiNumber
ifindiFit（i）recordPbest（i）=indiFit（i）;
tourPbest（i,:
）=individual（i,:
）;
end
ifindiFit（i）recordGbest=indiFit（i）;
tourGbest=individual（i,:
）;
end
end
[value,index]=min（recordPbest）;
recordGbest（N）=recordPbest（index）;
%%更新每一个粒子的位置
fori=1:
indiNumber
%个体最优值更新速度
c1=unidrnd（n-1）;
c2=unidrnd（n-1）;
whilec1==c2
c1=round（rand*（n-2））+1;
c2=round（rand*（n-2））+1;
end
chb1=min（c1,c2）;
chb2=max（c1,c2）;
cros=tourPbest（i,chb1:
chb2）;
ncros=size（cros,2）;
forj=1:
ncros
fork=1:
n
ifxnew1（i,k）==cros（j）
xnew1（i,k）=0;
fort=1:
n-k
temp=xnew1（i,k+t-1）;
xnew1（i,k+t-1）=xnew1（i,k+t）;
xnew1（i,k+t）=temp;
end
end
end
end
xnew1（i,n-ncros+1:
n）=cros;
dist=0;
forj=1:
n-1
dist=dist+cityDist（xnew1（i,j）,xnew1（i,j+1））;
end
dist=dist+cityDist（xnew1（i,1）,xnew1（i,n））;
ifindiFit（i）>dist
individual（i,:
）=xnew1（i,:
）;
end
%全体最优值更新速度
c1=round（rand*（n-2））+1;
c2=round（rand*（n-2））+1;
whilec1==c2
c1=round（rand*（n-2））+1;
c2=round（rand*（n-2））+1;
end
chb1=min（c1,c2）;
chb2=max（c1,c2）;
cros=tourGbest（chb1:
chb2）;
ncros=size（cros,2）;
forj=1:
ncros
fork=1:
n
ifxnew1（i,k）==cros（j）
xnew1（i,k）=0;
fort=1:
n-k
temp=xnew1（i,k+t-1）;
xnew1（i,k+t-1）=xnew1（i,k+t）;
xnew1（i,k+t）=temp;
end
end
end
end
xnew1（i,n-ncros+1:
n）=cros;
dist=0;
forj=1:
n-1
dist=dist+cityDist（xnew1（i,j）,xnew1（i,j+1））;
end
dist=dist+cityDist（xnew1（i,1）,xnew1（i,n））;
ifindiFit（i）>dist
individual（i,:
）=xnew1（i,:
）;
end
%%惯性项
c1=round（rand*（n-1））+1;
c2=round（rand*（n-1））+1;
whilec1==c2
c1=round（rand*（n-2））+1;
c2=round（rand*（n-2））+1;
end
temp=xnew1（i,c1）;
xnew1（i,c1）=xnew1（i,c2）;
xnew1（i,c2）=temp;
dist=0;
forj=1:
n-1
dist=dist+cityDist（xnew1（i,j）,xnew1（i,j+1））;
end
dist=dist+cityDist（xnew1（i,1）,xnew1（i,n））;
ifindiFit（i）>dist
individual（i,:
）=xnew1（i,:
）;
end
end
[value,index]=min（indiFit）;
L_best（N）=indiFit（index）;
tourGbest=individual（index,:
）;
end
%%其他函数
functiondist=dist（x,D）
n=size（x,2）;
dist=0;
fori=1:
n-1
dist=dist+D（x（i）,x（i+1））;
end
dist=dist+D（x
（1）,x（n））;
functionindiFit=fitness（x,cityCoor,cityDist）
m=size（x,1）;
n=size（cityCoor,1）;
indiFit=zeros（m,1）;
fori=1:
m
forj=1:
n-1
indiFit（i）=indiFit（i）+cityDist（x（i,j）,x（i,j+1））;
end
indiFit（i）=indiFit（i）+cityDist（x（i,1）,x（i,n））;
end
升级会员 