一些解决TSP问题算法源代码.docx
《一些解决TSP问题算法源代码.docx》由会员分享，可在线阅读，更多相关《一些解决TSP问题算法源代码.docx（79页珍藏版）》请在冰点文库上搜索。
一些解决TSP问题算法源代码
模拟退火算法
模拟退火算法来源于固体退火原理，将固体加温至充分高，再让其徐徐冷却，加温时，固体内部粒子随温升变为无序状，
内能增大，而徐徐冷却时粒子渐趋有序，在每个温度都达到平衡态，最后在常温时达到基态，内能减为最小。
根据Metropolis
准则，粒子在温度T时趋于平衡的概率为e-ΔE/（kT>，其中E为温度T时的内能，ΔE为其改变量，k为Boltzmann常数。
用固体退
火模拟组合优化问题，将内能E模拟为目标函数值f，温度T演化成控制参数t，即得到解组合优化问题的模拟退火算法：
由初始
解i和控制参数初值t开始，对当前解重复“产生新解→计算目标函数差→接受或舍弃”的迭代，并逐步衰减t值，算法终止时的
当前解即为所得近似最优解，这是基于蒙特卡罗迭代求解法的一种启发式随机搜索过程。
退火过程由冷却进度表（Cooling
Schedule>控制，包括控制参数的初值t及其衰减因子Δt、每个t值时的迭代次数L和停止条件S。
3.5.1模拟退火算法的模型
模拟退火算法可以分解为解空间、目标函数和初始解三部分。
模拟退火的基本思想:
（1>初始化：
初始温度T（充分大>，初始解状态S（是算法迭代的起点>，每个T值的迭代次数L
（2>对k=1，……，L做第（3>至第6步：
（3>产生新解S′
（4>计算增量Δt′=C（S′>-C（S>，其中C（S>为评价函数
（5>若Δt′<0则接受S′作为新的当前解，否则以概率exp（-Δt′/T>接受S′作为新的当前解.
（6>如果满足终止条件则输出当前解作为最优解，结束程序。
终止条件通常取为连续若干个新解都没有被接受时终止算法。
（7>T逐渐减少，且T->0，然后转第2步。
算法对应动态演示图：
模拟退火算法新解的产生和接受可分为如下四个步骤：
第一步是由一个产生函数从当前解产生一个位于解空间的新解；为便于后续的计算和接受，减少算法耗时，通常选择由当
前新解经过简单地变换即可产生新解的方法，如对构成新解的全部或部分元素进行置换、互换等，注意到产生新解的变换方法
决定了当前新解的邻域结构，因而对冷却进度表的选取有一定的影响。
第二步是计算与新解所对应的目标函数差。
因为目标函数差仅由变换部分产生，所以目标函数差的计算最好按增量计算。
事实表明，对大多数应用而言，这是计算目标函数差的最快方法。
第三步是判断新解是否被接受,判断的依据是一个接受准则，最常用的接受准则是Metropo1is准则:
若Δt′<0则接受S′作
为新的当前解S，否则以概率exp（-Δt′/T>接受S′作为新的当前解S。
第四步是当新解被确定接受时，用新解代替当前解，这只需将当前解中对应于产生新解时的变换部分予以实现，同时修正
目标函数值即可。
此时，当前解实现了一次迭代。
可在此基础上开始下一轮实验。
而当新解被判定为舍弃时，则在原当前解的
基础上继续下一轮实验。
模拟退火算法与初始值无关，算法求得的解与初始解状态S（是算法迭代的起点>无关；模拟退火算法具有渐近收敛性，已在
理论上被证明是一种以概率l收敛于全局最优解的全局优化算法；模拟退火算法具有并行性。
3.5.2模拟退火算法的简单应用
作为模拟退火算法应用，讨论货郎担问题（TravellingSalesmanProblem，简记为TSP>：
设有n个城市，用数码1,…,n代表
。
城市i和城市j之间的距离为d（i，j>i,j=1,…,n．TSP问题是要找遍访每个域市恰好一次的一条回路，且其路径总长度为最
短.。
求解TSP的模拟退火算法模型可描述如下：
解空间解空间S是遍访每个城市恰好一次的所有回路，是{1，……，n}的所有循环排列的集合，S中的成员记为（w1,w2,…
…，wn>，并记wn+1=w1。
初始解可选为（1，……，n>
目标函数此时的目标函数即为访问所有城市的路径总长度或称为代价函数：
我们要求此代价函数的最小值。
新解的产生随机产生1和n之间的两相异数k和m，若k （w1,w2,…，wk,wk+1,…，wm,…，wn>
变为：
（w1,w2,…，wm,wm-1,…，wk+1,wk,…，wn>.
如果是k>m，则将
（w1,w2,…，wk,wk+1,…，wm,…，wn>
变为：
（wm,wm-1,…，w1,wm+1,…，wk-1,wn,wn-1,…，wk>.
上述变换方法可简单说成是“逆转中间或者逆转两端”。
也可以采用其他的变换方法，有些变换有独特的优越性，有时也将它们交替使用，得到一种更好方法。
代价函数差设将（w1,w2,……，wn>变换为（u1,u2,……，un>,则代价函数差为：
根据上述分析，可写出用模拟退火算法求解TSP问题的伪程序：
ProcedureTSPSA:
begin
init-of-T。
{T为初始温度}
S={1，……，n}。
{S为初始值}
termination=false。
whiletermination=false
begin
fori=1toLdo
begin
generate（S′formS>。
{从当前回路S产生新回路S′}
Δt:
=f（S′>>-f（S>。
{f（S>为路径总长}
IF（Δt<0>OR（EXP（-Δt/T>>Random-of-[0,1]>
S=S′。
IFthe-halt-condition-is-TRUETHEN
termination=true。
End。
T_lower。
End。
End
模拟退火算法的应用很广泛，可以较高的效率求解最大截问题（MaxCutProblem>、0-1背包问题（ZeroOneKnapsack
Problem>、图着色问题（GraphColouringProblem>、调度问题（SchedulingProblem>等等。
3.5.3模拟退火算法的参数控制问题
模拟退火算法的应用很广泛，可以求解NP完全问题，但其参数难以控制，其主要问题有以下三点：
（1>温度T的初始值设置问题。
温度T的初始值设置是影响模拟退火算法全局搜索性能的重要因素之一、初始温度高，则搜索到全局最优解的可能性大，但
因此要花费大量的计算时间；反之，则可节约计算时间，但全局搜索性能可能受到影响。
实际应用过程中，初始温度一般需要
依据实验结果进行若干次调整。
（2>退火速度问题。
模拟退火算法的全局搜索性能也与退火速度密切相关。
一般来说，同一温度下的“充分”搜索（退火>是相当必要的，但这
需要计算时间。
实际应用中，要针对具体问题的性质和特征设置合理的退火平衡条件。
（3>温度管理问题。
温度管理问题也是模拟退火算法难以处理的问题之一。
实际应用中，因为必须考虑计算复杂度的切实可行性等问题，常采
用如下所示的降温方式：
T（t+1>＝k×T（t>
式中k为正的略小于1.00的常数，t为降温的次数
使用SA解决TSP问题的Matlab程序：
functionout=tsp（loc>
%TSPTravelingsalesmanproblem（TSP>usingSA（simulatedannealing>.
%TSPbyitselfwillgenerate20citieswithinaunitcubeand
%thenuseSAtoslovethisproblem.
%
%TSP（LOC>solvethetravelingsalesmanproblemwithcities'
%coordinatesgivenbyLOC,whichisanMby2matrixandMis
%thenumberofcities.
%
%Forexample:
%
%loc=rand（50,2>。
%tsp（loc>。
ifnargin==0,
%ThefollowingdataisfromthepostbyJenniferMyers（jmyers@nwu.
edu>
edu>
%tocomp.ai.neural-nets.It'sobtainedfromthefigurein
%Hopfield&Tank's1985paperinBiologicalCybernetics
%（Vol52,pp.141-152>.
loc=[0.3663,0.9076。
0.7459,0.8713。
0.4521,0.8465。
0.7624,0.7459。
0.7096,0.7228。
0.0710,0.7426。
0.4224,0.7129。
0.5908,0.6931。
0.3201,0.6403。
0.5974,0.6436。
0.3630,0.5908。
0.6700,0.5908。
0.6172,0.5495。
0.6667,0.5446。
0.1980,0.4686。
0.3498,0.4488。
0.2673,0.4274。
0.9439,0.4208。
0.8218,0.3795。
0.3729,0.2690。
0.6073,0.2640。
0.4158,0.2475。
0.5990,0.2261。
0.3927,0.1947。
0.5347,0.1898。
0.3960,0.1320。
0.6287,0.0842。
0.5000,0.0396。
0.9802,0.0182。
0.6832,0.8515]。
end
NumCity=length（loc>。
%Numberofcities
distance=zeros（NumCity>。
%Initializeadistancematrix
%Fillthedistancematrix
fori=1:
NumCity,
forj=1:
NumCity,
distance（i,j>=norm（loc（i,-loc（j,>。
distance（i,j>=norm（loc（i,-loc（j,>。
end
end
%Togenerateenergy（objectivefunction>frompath
%path=randperm（NumCity>。
%energy=sum（distance（（path-1>*NumCity+[path（2:
NumCity>path（1>]>>。
%FindtypicalvaluesofdE
count=20。
all_dE=zeros（count,1>。
fori=1:
count
path=randperm（NumCity>。
energy=sum（distance（（path-1>*NumCity+[path（2:
NumCity>
path（1>]>>。
new_path=path。
index=round（rand（2,1>*NumCity+.5>。
inversion_index=（min（index>:
max（index>>。
new_path（inversion_index>=fliplr（path（inversion_index>>。
all_dE（i>=abs（energy-...
sum（sum（diff（loc（[new_pathnew_path（1>],>'.^2>>>。
end
dE=max（all_dE>。
dE=max（all_dE>。
temp=10*dE。
%Choosethetemperaturetobelargeenough
fprintf（'Initialenergy=%f\n\n',energy>。
%Initialplots
out=[pathpath（1>]。
plot（loc（out（,1>,loc（out（,2>,'r.','Markersize',20>。
axissquare。
holdon
h=plot（loc（out（,1>,loc（out（,2>>。
holdoff
MaxTrialN=NumCity*100。
%Max.#oftrialsata
temperature
MaxAcceptN=NumCity*10。
%Max.#ofacceptancesata
temperature
StopTolerance=0.005。
%Stoppingtolerance
TempRatio=0.5。
%Temperaturedecreaseratio
minE=inf。
%Initialvalueformin.energy
maxE=-1。
%Initialvalueformax.energy
%Majorannealingloop
while（maxE-minE>/maxE>StopTolerance,
minE=inf。
minE=inf。
maxE=0。
TrialN=0。
%Numberoftrialmoves
AcceptN=0。
%Numberofactualmoves
whileTrialNnew_path=path。
index=round（rand（2,1>*NumCity+.5>。
inversion_index=（min（index>:
max（index>>。
new_path（inversion_index>=
fliplr（path（inversion_index>>。
new_energy=sum（distance（...
（new_path-1>*NumCity+[new_path（2:
NumCity>
new_path（1>]>>。
ifrand/temp>,%
accept
it!
energy=new_energy。
path=new_path。
minE=min（minE,energy>。
maxE=max（maxE,energy>。
AcceptN=AcceptN+1。
end
TrialN=TrialN+1。
end
end
%Updateplot
out=[pathpath（1>]。
set（h,'xdata',loc（out（,1>,'ydata',loc（out（,2>>。
drawnow。
%Printinformationincommandwindow
fprintf（'temp.=%f\n',temp>。
tmp=sprintf（'%d',path>。
fprintf（'path=%s\n',tmp>。
fprintf（'energy=%f\n',energy>。
fprintf（'[minEmaxE]=[%f%f]\n',minE,maxE>。
fprintf（'[AcceptNTrialN]=[%d%d]\n\n',AcceptN,TrialN>。
%Lowerthetemperature
temp=temp*TempRatio。
end
%Printsequentialnumbersinthegraphicwindow
fori=1:
NumCity,
text（loc（path（i>,1>+0.01,loc（path（i>,2>+0.01,int2str（i>,...
'fontsize',8>。
end
又一个相关的Matlab程序
function[X,P]=zkp（w,c,M,t0,tf>
[m,n]=size（w>。
L=100*n。
t=t0。
clearm。
x=zeros（1,n>
xmax=x。
fmax=0。
whilet>tf
fork=1:
L
xr=change（x>
gx=g_0_1（w,x>。
gxr=g_0_1（w,xr>。
ifgxr<=M
fr=f_0_1（xr,c>。
f=f_0_1（x,c>。
df=fr-f。
ifdf>0
x=xr。
iffr>fmax
fmax=fr。
xmax=xr。
end
else
p=rand。
ifp
x=xr。
end
end
end
end
t=0.87*t
end
P=fmax。
X=xmax。
%下面的函数产生新解
function[d_f,pi_r]=exchange_2（pi0,d>
[m,n]=size（d>。
clearm。
u=rand。
u=u*（n-2>。
u=round（u>。
ifu<2
u=2。
end
ifu>n-2
u=n-2。
end
v=rand。
v=v*（n-u+1>。
v=round（v>。
ifv<1
v=1。
end
v=v+u。
ifv>n
v=n。
end
pi_1（u>=pi0（v>。
pi_1（u>=pi0（u>。
ifu>1
fork=1:
（u-1>
pi_1（k>=pi0（k>。
end
end
ifv>（u+1>
fork=1:
（v-u-1>
pi_1（u+k>=pi0（v-k>。
end
end
ifvfork=（v+1>:
n
pi_1（k>=pi0（k>。
end
end
d_f=0。
ifvd_f=d（pi0（u-1>,pi0（v>>+d（pi0（u>,pi0（v+1>>。
fork=（u+1>:
n
d_f=d_f+d（pi0（k>,pi0（k-1>>-d（pi0（v>,pi0（v+1>>。
end
d_f=d_f-d（pi0（u-1>,pi0（u>>。
fork=（u+1>:
n
d_f=d_f-d（pi0（k-1>,pi0（k>>。
end
else
d_f=d（pi0（u-1>,pi0（v>>+d（pi0（u>,pi0（1>>-d（pi0（u-1>,pi0（u>>-d（pi0（v>,pi0（1>>。
fork=（u+1>:
n
d_f=d_f-d（pi0（k>,pi0（k-1>>。
end
fork=（u+1>:
n
d_f=d_f-d（pi0（k-1>,pi0（k>>。
end
end
pi_r=pi_1。
遗传算法GA
遗传算法：
旅行商问题（travelingsalemanproblem,简称tsp>：
已知n个城市之间的相互距离，现有一个推销员必须遍访这n个城市，并且每个城市只能访问一次，最后又必须返回出发城市。
如何安排他对这些城市的访问次序，可使其旅行路线的总长度最短？
用图论的术语来说，假设有一个图g=（v,e>，其中v是顶点集，e是边集，设d=（dij>是由顶点i和顶点j之间的距离所组成的距离矩阵，旅行商问题就是求出一条通过所有顶点且每个顶点只通过一次的具有最短距离的回路。
这个问题可分为对称旅行商问题（dij=dji,,任意i,j=1,2,3，…,n>和非对称旅行商问题（dij≠dji,,任意i,j=1,2,3，…,n>。
若对于城市v={v1,v2,v3,…,vn}的一个访问顺序为t=（t1,t2,t3,…,ti,…,tn>,其中ti∈v（i=1,2,3,…,n>，且记tn+1=t1，则旅行商问题的数学模型为：
minl=σd（t（i>,t（i+1>>旅行商问题是一个典型的组合优化问题，并且是一个np难问题，其可能的路径数目与城市数目n是成指数型增长的，所以一般很难精确地求出其最优解，本文采用遗传算法求其近似解。
遗传算法：
初始化过程：
用v1,v2,v3,…,vn代表所选n个城市。
定义整数pop-size作为染色体的个数，并且随机产生pop-size个初始染色体，每个染色体为1到18的整数组成的随机序列。
适应度f的计算：
对种群中的每个染色体vi，计算其适应度，f=σd（t（i>,t（i+1>>.
评价函数eval（vi>：
用来对种群中的每个染色体vi设定一个概率，以使该染色体被选中的可能性与其种群中其它染色体的适应性成比例，既通过轮盘赌，适应性强的染色体被选择产生后台的机会要大，设alpha∈（0,1>，本文定义基于序的评价函数为eval（vi>=alpha*（1-alpha>.^（i-1>。
[随机规划与模糊规划]
选择过程：
选择过程是以旋转赌轮pop-size次为基础，每次旋转都为新的种群选择一个染色体。
赌轮是按每个染色体的适应度进行选择染色体的。
step1、对每个染色体vi,计算累计概率qi，q0=0。
qi=σeval（vj>j=1,…,i。
i=1,…pop-size.
step2、从区间（0,pop-size>中产生一个随机数r；
step3、若qi-1step4、重复step2和step3共pop-size次，这样可以得到pop-size个复制的染色体。
grefenstette编码：
因为常规的交叉运算和变异运算会使种群中产生一些无实际意义的染色体，本文采用grefenstette编码《遗传算法原理及应用》可以避免这种情况的出现。
所谓的grefenstette编码就是用所选队员在未选<不含淘汰）队员中的位置，如：
815216107431114612951813171
对应：
81421386325734324221。
交叉过程：
本文采用常规单点交叉。
为确定交叉操作的父代，从到pop-size重复以下过程：
从[0，1]中产生一个随机数r，如果r将所选的父代两两组队，随机产生一个位置进行交叉，如：
81421386325734324221
6123568563185633211
交叉后为：
81421386325185633211
612356856373
升级会员 