一些解决TSP问题的算法及源代码.doc

模拟退火算法 模拟退火算法来源于固体退火原理，将固体加温至充分高，再让其徐徐冷却，加温时，固体内部粒子随温升变为无序状， 内能增大，而徐徐冷却时粒子渐趋有序，在每个温度都达到平衡态，最后在常温时达到基态，内能减为最小。根据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) 若t0，然后转第2步。 算法对应动态演示图： 模拟退火算法新解的产生和接受可分为如下四个步骤： 第一步是由一个产生函数从当前解产生一个位于解空间的新解；为便于后续的计算和接受，减少算法耗时，通常选择由当 前新解经过简单地变换即可产生新解的方法，如对构成新解的全部或部分元素进行置换、互换等，注意到产生新解的变换方法 决定了当前新解的邻域结构，因而对冷却进度表的选取有一定的影响。 第二步是计算与新解所对应的目标函数差。因为目标函数差仅由变换部分产生，所以目标函数差的计算最好按增量计算。 事实表明，对大多数应用而言，这是计算目标函数差的最快方法。 第三步是判断新解是否被接受,判断的依据是一个接受准则，最常用的接受准则是Metropo1is准则: 若t0则接受S作 为新的当前解S，否则以概率exp(-t/T)接受S作为新的当前解S。 第四步是当新解被确定接受时，用新解代替当前解，这只需将当前解中对应于产生新解时的变换部分予以实现，同时修正 目标函数值即可。此时，当前解实现了一次迭代。可在此基础上开始下一轮试验。而当新解被判定为舍弃时，则在原当前解的 基础上继续下一轮试验。 模拟退火算法与初始值无关，算法求得的解与初始解状态S(是算法迭代的起点)无关；模拟退火算法具有渐近收敛性，已在 理论上被证明是一种以概率l 收敛于全局最优解的全局优化算法；模拟退火算法具有并行性。 3.5.2 模拟退火算法的简单应用 作为模拟退火算法应用，讨论货郎担问题(Travelling Salesman Problem，简记为TSP)：设有n个城市，用数码1,n代表 。城市i和城市j之间的距离为d(i，j) i, j=1,nTSP问题是要找遍访每个域市恰好一次的一条回路，且其路径总长度为最 短.。 求解TSP的模拟退火算法模型可描述如下： 解空间 解空间S是遍访每个城市恰好一次的所有回路，是1，n的所有循环排列的集合，S中的成员记为(w1,w2 , ，wn)，并记wn+1= w1。初始解可选为(1，n) 目标函数 此时的目标函数即为访问所有城市的路径总长度或称为代价函数： 我们要求此代价函数的最小值。 新解的产生 随机产生1和n之间的两相异数k和m，若km，则将 (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问题的伪程序： Procedure TSPSA: begin init-of-T; T为初始温度 S=1，n; S为初始值 termination=false; while termination=false begin for i=1 to L do begin generate(Sform S); 从当前回路S产生新回路S t:=f(S)-f(S);f(S)为路径总长 IF(tRandom-of-0,1) S=S; IF the-halt-condition-is-TRUE THEN termination=true; End; T_lower; End; End 模拟退火算法的应用很广泛，可以较高的效率求解最大截问题(Max Cut Problem)、0-1背包问题(Zero One Knapsack Problem)、图着色问题(Graph Colouring Problem)、调度问题(Scheduling Problem)等等。 3.5.3 模拟退火算法的参数控制问题 模拟退火算法的应用很广泛，可以求解NP完全问题，但其参数难以控制，其主要问题有以下三点： (1) 温度T的初始值设置问题。 温度T的初始值设置是影响模拟退火算法全局搜索性能的重要因素之一、初始温度高，则搜索到全局最优解的可能性大，但 因此要花费大量的计算时间；反之，则可节约计算时间，但全局搜索性能可能受到影响。实际应用过程中，初始温度一般需要 依据实验结果进行若干次调整。 (2) 退火速度问题。 模拟退火算法的全局搜索性能也与退火速度密切相关。一般来说，同一温度下的“充分”搜索(退火)是相当必要的，但这 需要计算时间。实际应用中，要针对具体问题的性质和特征设置合理的退火平衡条件。 (3) 温度管理问题。 温度管理问题也是模拟退火算法难以处理的问题之一。实际应用中，由于必须考虑计算复杂度的切实可行性等问题，常采 用如下所示的降温方式： T(t+1)kT(t) 式中k为正的略小于1.00的常数，t为降温的次数 使用SA解决TSP问题的Matlab程序：function out = tsp(loc)% TSP Traveling salesman problem (TSP) using SA (simulated annealing).% TSP by itself will generate 20 cities within a unit cube and% then use SA to slove this problem.% TSP(LOC) solve the traveling salesman problem with cities% coordinates given by LOC, which is an M by 2 matrix and M is% the number of cities.% For example:% loc = rand(50, 2);% tsp(loc);if nargin = 0,% The following data is from the post by Jennifer Myers ()edu)% to comp.ai.neural-nets. Its obtained from the figure in% Hopfield & Tanks 1985 paper in Biological Cybernetics% (Vol 52, 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;endNumCity = length(loc); % Number of citiesdistance = zeros(NumCity); % Initialize a distance matrix% Fill the distance matrixfor i = 1:NumCity,for j = 1:NumCity,distance(i, j) = norm(loc(i, - loc(j, );distance(i, j) = norm(loc(i, - loc(j, );endend% To generate energy (objective function) from path%path = randperm(NumCity);%energy = sum(distance(path-1)*NumCity + path(2:NumCity) path(1);% Find typical values of dEcount = 20;all_dE = zeros(count, 1);for i = 1:countpath = 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_path new_path(1),).2);enddE = max(all_dE);dE = max(all_dE);temp = 10*dE; % Choose the temperature to be large enoughfprintf(Initial energy = %fnn,energy);% Initial plotsout = path path(1);plot(loc(out(, 1), loc(out(, 2),r., Markersize, 20);axis square; hold onh = plot(loc(out(, 1), loc(out(, 2); hold offMaxTrialN = NumCity*100; % Max. # of trials at atemperatureMaxAcceptN = NumCity*10; % Max. # of acceptances at atemperatureStopTolerance = 0.005; % Stopping toleranceTempRatio = 0.5; % Temperature decrease ratiominE = inf; % Initial value for min. energymaxE = -1; % Initial value for max. energy% Major annealing loopwhile (maxE - minE)/maxE StopTolerance,minE = inf;minE = inf;maxE = 0;TrialN = 0; % Number of trial movesAcceptN = 0; % Number of actual moveswhile TrialN MaxTrialN & AcceptN MaxAcceptN,new_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);if rand tffor k=1:Lxr=change(x)gx=g_0_1(w,x);gxr=g_0_1(w,xr);if gxr0x=xr;if frfmaxfmax=fr;xmax=xr;endelsep=rand;if pexp(df/t)x=xr;endendendendt=0.87*tendP=fmax;X=xmax;%下面的函数产生新解function d_f,pi_r=exchange_2(pi0,d)m,n=size(d);clear m;u=rand;u=u*(n-2);u=round(u);if un-2u=n-2;endv=rand;v=v*(n-u+1);v=round(v);if vnv=n;endpi_1(u)=pi0(v);pi_1(u)=pi0(u);if u1for k=1:(u-1)pi_1(k)=pi0(k);endendif v(u+1)for k=1:(v-u-1)pi_1(u+k)=pi0(v-k);endendif vnfor k=(v+1):npi_1(k)=pi0(k);endendd_f=0;if vnd_f=d(pi0(u-1),pi0(v)+d(pi0(u),pi0(v+1);for k=(u+1):nd_f=d_f+d(pi0(k),pi0(k-1)-d(pi0(v),pi0(v+1);endd_f=d_f-d(pi0(u-1),pi0(u);for k=(u+1):nd_f=d_f-d(pi0(k-1),pi0(k); endelsed_f=d(pi0(u-1),pi0(v)+d(pi0(u),pi0(1)-d(pi0(u-1),pi0(u)-d(pi0(v),pi0(1);for k=(u+1):nd_f=d_f-d(pi0(k),pi0(k-1); endfor k=(u+1):nd_f=d_f-d(pi0(k-1),pi0(k);endendpi_r=pi_1; 遗传算法GA 遗传算法： 旅行商问题(traveling saleman problem,简称tsp)：已知n个城市之间的相互距离，现有一个推销员必须遍访这n个城市，并且每个城市只能访问一次，最后又必须返回出发城市。如何安排他对这些城市的访问次序，可使其旅行路线的总长度最短？用图论的术语来说，假设有一个图 g=(v,e)，其中v是顶点集，e是边集，设d=(dij)是由顶点i和顶点j之间的距离所组成的距离矩阵，旅行商问题就是求出一条通过所有顶点且每个顶点只通过一次的具有最短距离的回路。这个问题可分为对称旅行商问题(dij=dji,任意i,j=1,2,3，,n)和非对称旅行商问题(dijdji,任意i,j=1,2,3，,n)。若对于城市v=v1,v2,v3,vn的一个访问顺序为t=(t1,t2,t3,ti,tn),其中tiv(i=1,2,3,n)，且记tn+1= t1，则旅行商问题的数学模型为：min l=d(t(i),t(i+1) （i=1,n）旅行商问题是一个典型的组合优化问题，并且是一个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-1rqi,则选择第i个染色体 ；step4、重复step2和step3共pop-size次，这样可以得到pop-size个复制的染色体。grefenstette编码：由于常规的交叉运算和变异运算会使种群中产生一些无实际意义的染色体，本文采用grefenstette编码遗传算法原理及应用可以避免这种情况的出现。所谓的grefenstette编码就是用所选队员在未选（不含淘汰）队员中的位置，如：8 15 2 16 10 7 4 3 11 14 6 12 9 5 18 13 17 1对应：8 14 2 13 8 6 3 2 5 7 3 4 3 2 4 2 2 1。交叉过程：本文采用常规单点交叉。为确定交叉操作的父代，从 到pop-size重复以下过程：从0，1中产生一个随机数r，如果rpc ，则选择vi作为一个父代。将所选的父代两两组队，随机产生一个位置进行交叉，如：8 14 2 13 8 6 3 2 5 7 3 4 3 2 4 2 2 16 12 3 5 6 8 5 6 3 1 8 5 6 3 3 2 1 1交叉后为：8 14 2 13 8 6 3 2 5 1 8 5 6 3 3 2 1 16 12 3 5 6 8 5 6 3 7 3 4 3 2 4 2 2 1变异过程：本文采用均匀多点变异。类似交叉操作中选择父代的过程，在rpm 的标准下选择多个染色体vi作为父代。对每一个选择的父代，随机选择多个位置，使其在每位置按均匀变异（该变异点xk的取值范围为ukmin,ukmax,产生一个0，1中随机数r，该点变异为xk=ukmin+r(ukmax-ukmin)）操作。如：8 14 2 13 8 6 3 2 5 7 3 4 3 2 4 2 2 1变异后：8 14 2 13 10 6 3 2 2 7 3 4 5 2 4 1 2 1反grefenstette编码：交叉和变异都是在grefenstette编码之后进行的，为了循环操作和返回最终结果，必须逆grefenstette编码过程，将编码恢复到自然编码。循环操作：判断是否满足设定的带数xzome，否，则跳入适应度f的计算；是，结束遗传操作，跳出。 Matlab程序：function bestpop,trace=ga(d,termops,num,pc,cxops,pm,alpha)%bestpop,trace=ga(d,termops,num,pc,cxops,pm,alpha)%d:距离矩阵%termops:种群代数%num:每代染色体的个数%pc:交叉概率%cxops:由于本程序采用单点交叉，交叉点的设置在本程序中没有很好的解决，所以本文了采用定点，即第cxops，可以随机产生。%pm:变异概率%alpha:评价函数eval(vi)=alpha*(1-alpha).(i-1).%bestpop:返回的最优种群%trace:进化轨迹%-%#版权所有！欢迎广大网友改正，改进！#%e-mail:%#%citynum=size(d,2);n=nargin;if n2disp(缺少变量！)disp(_开个玩笑_)endif n2termops=500;num=50;pc=0.25;cxops=3;pm=0.30;alpha=0.10;endif n3num=50;pc=0.25;cxops=3;pm=0.30;alpha=0.10;endif n4pc=0.25;cxops=3;pm=0.30;alpha=0.10;endif n5cxops=3;pm=0.30;alpha=0.10;endif n6pm=0.30;alpha=0.10;endif n7alpha=0.10;endif isempty(cxops)cxops=3;end t=initializega(num,citynum);for i=1:termopsl=f(d,t);x,y=find(l=max(l);trace(i)=-l(y(1);bestpop=t(y(1),:);t=select(t,l,alpha);g=grefenstette(t);g1=crossover(g,pc,cxops);g=mutation(g1,pm); %均匀变异t=congrefenstette(g);end -function t=initializega(num,citynum)for i=1:numt(i,:)=randperm(citynum);end-function l=f(d,t)m,n=size(t);for k=1:mfor i=1:n-1l(k,i)=d(t(k,i),t(k,i+1);endl(k,n)=d(t(k,n),t(k,1);l(k)=-sum(l(k,:);end-function t=select(t,l,alpha)m,n=size(l);t1=t;beforesort,aftersort1=sort(l,2);%fsort from l to ufor i=1:naftersort(i)=aftersort1(n+1-i); %change endfor k=1:n;t(k,:)=t1(aftersort(k),:);l1(k)=l(aftersort(k);endt1=t;l=l1;for i=1:size(aftersort,2)evalv(i)=alpha*(1-alpha).(i-1);endm=size(t,1);q=cumsum(evalv);qmax=max(q);for k=1:mr=qmax*rand(1);for j=1:mif j=1&rq(j-1)&r=q(j)t(k,:)=t1(j,:);endendend-function g=grefenstette(t)m,n=size(t);for k=1:mt0=1:n;for i=1:nfor j=1:length(t0)if t(k,i)=t0(j)g(k,i)=j;t0(j)=;breakendendendend-function g=crossover(g,pc,cxops)m,n=size(g);ran=rand(1,m);r=cxops;x,ru=find(ran=2for k=1:2:length(ru)-1g1(ru(k),:)=g(ru(k),1:r),g(ru(k+1),(r+1):n);g(ru(k+1),:)=g(ru(k+1),1:r),g(ru(k),(r+1):n);g(ru(k),:)=g1(ru(k),:);endend-function g=mutation(g,pm) %均匀变异m,n=size(g);ran=rand(1,m);r=rand(1,3); %dai gai jinrr=floor(n*rand(1,3)+1);x,mu=find(ranpm);for k=1:length(mu)for i=1:length(r)umax(i)=n+1-rr(i);umin(i)=1;g(mu(k),rr(i)=umin(i)+floor(umax(i)-umin(i)*r(i);endend-function t=congrefenstette(g)m,n=size(g);for k=1:mt0=1:n;for i=1:nt(k,i)=t0(g(k,i);t0(g(k,i)=;endend-又一个Matlab程序，其中交叉算法采用的是由Goldberg和Lingle于1985年提出的PMX(部分匹配交叉)，淘汰保护指数alpha是我自己设计的，起到了加速优胜劣汰的作用。%TSP问题（又名：旅行商问题，货郎担问题）遗传算法通用matlab程序%D是距离矩阵，n为种群个数,建议取为城市个数的12倍，%C为停止代数，遗传到第 C代时程序停止,C的具体取值视问题的规模和耗费的时间而定%m为适应值归一化淘汰加速指数 ,最好取为1,2,3,4 ,不宜太大%alpha为淘汰保护指数,可取为01之间任意小数,取1时关闭保护功能,最好取为0.81.0%R为最短路径,Rlength为路径长度function R,Rlength=geneticTSP(D,n,C,m,alpha) N,NN=size(D);farm=zeros(n,N);%用于存储种群for i=1:nfarm(i,:)=randperm(N);%随机生成初始种群endR=farm(1,:);%存储最优种群len=zeros(n,1);%存储路径长度fitness=zeros(n,1);%存储归一化适应值counter=0; while counter=alpha*randnn=nn+1;FARM(nn,:)=farm(i,:);endendFARM=FARM(1:nn,:); aa,bb=size(FARM);%交叉和变异while aanif nnnFARM=FARM(1:n,:);%保持种群规模为nend farm=FARM;clear FARMcounter=counter+1 end Rlength=myLength(D,R); function a,b=intercross(a,b)L=length(a);if L=rand&L10W=ceil(L/10);else W=floor(L/10);endp=unidrnd(L-W+1);%随机选择交叉范围，从p到p+Wfor i=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); endfunction x,y=exchange(x,y)temp=x;x=y;y=temp; % 计算路径的子程序function len=myLength(D,p)N,NN=size(D);len=D(p(1,N),p(1,1);for i=1:(N-1)len=len+D(p(1,i),p(1,i+1);end %计算归一化适应值子程序function fitness=fit(len,m,maxlen,minlen)fitness=len;for i=1:length(len)fitness(i,1)=(1-(len(i,1)-minlen)/(maxlen-minlen+0.000001).m;end 一个C+的程序：/c+的程序#include#includetemplateclass Graphpublic:Graph(int vertices=10)n=vertices;e=0;Graph()virtual bool Add(int u,int v,const T& w)=0;virtual bool Delete(int u,int v)=0;virtual bool Exist(int u,int v)const=0;int Vertices()constreturn n;int Edges()constreturn e;protected:int n;int e;templateclass MGraph:public Graphpublic:MGraph(int Vertices=10,T noEdge=0);MGraph();bool Add(int u,int v,const T& w);bool Delete(int u,int v);bool Exist(int u,int v)const;void Floyd(T*& d,int*& path);void print(int Vertices);private:T NoEdge;T* a;templateMGraph:MGraph(int Vertices,T noEdge)n=Vertices;NoEdge=noEdge;a=new T* n;for(int i=0;in;i+)ai=new Tn;aii=0;for(int j=0;jn;j+)if(i!=j)aij=NoEdge;templateMGraph:MGraph()for(int i=0;in;i+)deleteai;deletea;templatebool MGraph:Exist(int u,int v)constif(u0|vn-1|vn-1|u=v|auv=NoEdge)return false;return true;templatebool MGraph:Add(int u,int v,const T& w)if(u0|vn-1|vn-1|u=v|auv!=NoEdge)cerrBadInput!endl;return false;auv=w;e+;return true;templatebool MGraph:delete(int u,int v)if(u0|vn-1|vn-1|u=v|auv=NoEdge)cerrBadInput!endl;return false;auv=NoEdge;e-;return true;templatevoid MGraph:Floyd(T*& d,int*& path)d=new T* n;path=new int* n;for(int i=0;in;i+)di=new Tn;pathi=new intn;for(int j=0;jn;j+)dij=aij;if(i!=j&aijNoEdge)pathij=i;else pathij=-1;for(int k=0;kn;k+)for(i=0;in;i+)for(int j=0;jn;j+)if(dik+dkjdi