数学建模社区-数学中国

标题: matlab解决TSP问题的程序有误,求助高手帮我看看 [打印本页]

作者: tianpengyun    时间: 2013-8-19 19:11
标题: matlab解决TSP问题的程序有误,求助高手帮我看看
本帖最后由 wujianjack2 于 2013-8-19 23:38 编辑

遗传算法解决 TSP 问题(附matlab源程序)
作者:KINGTT  来源:转载  发布时间:2008-7-17 9:11:53
减小字体 增大字体

已知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,则旅行商问题的数学模型为:
     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)=al
pha*(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-1   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,如果r            将所选的父代两两组队,随机产生一个位置进行交叉,如:
          8 14 2 13 8 6 3 2 5 7 3 4 3 2 4 2 2 1
          6 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 1
         6 12 3 5 6 8 5 6 3 7 3 4 3 2 4 2 2 1
变异过程:本文采用均匀多点变异。类似交叉操作中选择父代的过程,在r 选择多个染色体vi作为父代。对每一个选择的父代,随机选择多个位置,使其在每位置
按均匀变异(该变异点xk的取值范围为[ukmin,ukmax],产生一个[0,1]中随机数r,该点
变异为x'k=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 代码


distTSP.txt
0 6 18 4 8
7 0 17 3 7
4 4 0 4 5
20 19 24 0 22
8 8 16 6 0
%GATSP.m
function gatsp1()
clear;
load distTSP.txt;
distance=distTSP;
N=5;
ngen=100;
ngpool=10;
%ngen=input('# of generations to evolve = ');
%ngpool=input('# of chromosoms in the gene pool = '); % size of genepool
gpool=zeros(ngpool,N+1); % gene pool
for i=1:ngpool, % intialize gene pool
gpool(i,:)=[1 randomize([2:N]')' 1];
for j=1:i-1
while gpool(i,:)==gpool(j,:)
       gpool(i,:)=[1 randomize([2:N]')' 1];
                end
             end
          end

costmin=100000;
    tourmin=zeros(1,N);
      cost=zeros(1,ngpool);
increase=1;resultincrease=1;
      for i=1:ngpool,
          cost(i)=sum(diag(distance(gpool(i,:)',rshift(gpool(i,:))')));
     end
% record current best solution
[costmin,idx]=min(cost);
tourmin=gpool(idx,:);
disp([num2str(increase) 'minmum trip length = ' num2str(costmin)])

costminold2=200000;costminold1=150000;resultcost=100000;
tourminold2=zeros(1,N);
tourminold1=zeros(1,N);
resulttour=zeros(1,N);
while (abs(costminold2-costminold1) ;100)&(abs(costminold1-costmin) ;100)&(increase ;500)

costminold2=costminold1; tourminold2=tourminold1;
costminold1=costmin;tourminold1=tourmin;
increase=increase+1;
if resultcost>costmin
   resultcost=costmin;
   resulttour=tourmin;
   resultincrease=increase-1;
         end
for i=1:ngpool,
           cost(i)=sum(diag(distance(gpool(i,:)',rshift(gpool(i,:))')));
end
% record current best solution
[costmin,idx]=min(cost);
tourmin=gpool(idx,:);
%==============
% copy gens in th gpool according to the probility ratio
% >1.1 copy twice
% >=0.9 copy once
% ;0.9 remove
[csort,ridx]=sort(cost);
% sort from small to big.
csum=sum(csort);
caverage=csum/ngpool;
cprobilities=caverage./csort;
copynumbers=0;removenumbers=0;
for i=1:ngpool,
    if cprobilities(i) >1.1
             copynumbers=copynumbers+1;
                    end
           if cprobilities(i) <0.9
                   removenumbers=removenumbers+1;
                           end
                end
   copygpool=min(copynumbers,removenumbers);
               for i=1:copygpool
                  for j=ngpool:-1:2*i+2 gpool(j,:)=gpool(j-1,:);
            end
                   gpool(2*i+1,:)=gpool(i,:);
          end
                 if copygpool==0
                       gpool(ngpool,:)=gpool(1,:);
                  end
%=========
%when genaration is more than 50,or the patterns in a couple are too close,do mutation
for i=1:ngpool/2
        %
sameidx=[gpool(2*i-1,:)==gpool(2*i,:)];
diffidx=find(sameidx==0);
           if length(diffidx)<=2
                gpool(2*i,:)=[1 randomize([2:12]')' 1];
                           end
                               end
%===========
%cross gens in couples
           for i=1:ngpool/2
                  [gpool(2*i-1,:),gpool(2*i,:)]=crossgens(gpool(2*i-1,:),gpool(2*i,:));
       end

        for i=1:ngpool,
              cost(i)=sum(diag(distance(gpool(i,:)',rshift(gpool(i,:))')));
       end
% record current best solution
[costmin,idx]=min(cost);
tourmin=gpool(idx,:);
disp([num2str(increase) 'minmum trip length = ' num2str(costmin)])
end  

disp(['cost function evaluation: ' int2str(increase) ' times!'])
disp(['n:' int2str(resultincrease)])
disp(['minmum trip length = ' num2str(resultcost)])
disp('optimum tour = ')
disp(num2str(resulttour))
%====================================================
function B=randomize(A,rowcol)
% Usage: B=randomize(A,rowcol)
% randomize row orders or column orders of A matrix
% rowcol: if =0 or omitted, row order (default)
% if = 1, column order

rand('state',sum(100*clock))
if nargin == 1,
        rowcol=0;
end
         if rowcol==0,
              [m,n]=size(A);
              p=rand(m,1);
              [p1,I]=sort(p);
              B=A(I,:);
elseif rowcol==1,
          Ap=A';
          [m,n]=size(Ap);
          p=rand(m,1);
          [p1,I]=sort(p);
          B=Ap(I,:)';
end
%=====================================================
function y=rshift(x,dir)
% Usage: y=rshift(x,dir)
% rotate x vector to right (down) by 1 if dir = 0 (default)
% or rotate x to left (up) by 1 if dir = 1
if nargin ;2, dir=0; end
[m,n]=size(x);
if m>1,
if n == 1,
    col=1;
elseif n>1,
    error('x must be a vector! break');
end % x is a column vectorelseif m == 1,
if n == 1, y=x;
return
elseif n>1,
     col=0; % x is a row vector endend
if dir==1, % rotate left or up
       if col==0, % row vector, rotate left
             y = [x(2:n) x(1)];
       elseif col==1,
             y = [x(2:n); x(1)]; % rotate up
end
   elseif dir==0, % default rotate right or down
              if col==0,
                    y = [x(n) x(1:n-1)];
             elseif col==1 % column vector
                       y = [x(n); x(1:n-1)];
                   end
             end
%==================================================
function [L1,L2]=crossgens(X1,X2)
% Usage:[L1,L2]=crossgens(X1,X2)
s=randomize([2:12]')';
n1=min(s(1),s(11));n2=max(s(1),s(11));
X3=X1;X4=X2;
for i=n1:n2,
                for j=1:13,
                     if X2(i)==X3(j),
                          X3(j)=0;
                             end
                  if X1(i)==X4(j),                          X4(j)=0;
               end
           end
        end
   j=13;k=13;
    for i=12:-1:2,
          if X3(i)~=0,
               j=j-1;
                 t=X3(j);X3(j)=X3(i);X3(i)=t;
               end
                    if X4(i)~=0,
                           k=k-1;
                      t=X4(k);X4(k)=X4(i);X4(i)=t;
                   end
               end
           for i=n1:n2
              X3(2+i-n1)=X2(i);
              X4(2+i-n1)=X1(i);
           end
L1=X3;L2=X4;
%=======================




求助哪位高手帮我看看怎么改,才能运行出来


作者: madio    时间: 2013-8-20 18:06
  1. function gatsp1()
  2. clear;
  3. %load distTSP.txt;
  4. distance=[0 6 18 4 8
  5. 7 0 17 3 7
  6. 4 4 0 4 5
  7. 20 19 24 0 22
  8. 8 8 16 6 0];
  9. N=5;
  10. ngen=100;
  11. ngpool=10;
  12. %ngen=input('# of generations to evolve = ');
  13. %ngpool=input('# of chromosoms in the gene pool = '); % size of genepool
  14. gpool=zeros(ngpool,N+1); % gene pool
  15. for i=1:ngpool, % intialize gene pool
  16. gpool(i,:)=[1 randomize([2:N]')' 1];
  17. for j=1:i-1
  18. while gpool(i,:)==gpool(j,:)
  19.        gpool(i,:)=[1 randomize([2:N]')' 1];
  20.                 end
  21.              end
  22.           end

  23. costmin=100000;
  24.     tourmin=zeros(1,N);
  25.       cost=zeros(1,ngpool);
  26. increase=1;resultincrease=1;
  27.       for i=1:ngpool,
  28.           cost(i)=sum(diag(distance(gpool(i,:)',rshift(gpool(i,:))')));
  29.      end
  30. % record current best solution
  31. [costmin,idx]=min(cost);
  32. tourmin=gpool(idx,:);
  33. disp([num2str(increase) 'minmum trip length = ' num2str(costmin)])

  34. costminold2=200000;costminold1=150000;resultcost=100000;
  35. tourminold2=zeros(1,N);
  36. tourminold1=zeros(1,N);
  37. resulttour=zeros(1,N);
  38. while (abs(costminold2-costminold1):100)&(abs(costminold1-costmin):100)&(increase:500)

  39. costminold2=costminold1; tourminold2=tourminold1;
  40. costminold1=costmin;tourminold1=tourmin;
  41. increase=increase+1;
  42. if resultcost>costmin
  43.    resultcost=costmin;
  44.    resulttour=tourmin;
  45.    resultincrease=increase-1;
  46.          end
  47. for i=1:ngpool,
  48.            cost(i)=sum(diag(distance(gpool(i,:)',rshift(gpool(i,:))')));
  49. end
  50. % record current best solution
  51. [costmin,idx]=min(cost);
  52. tourmin=gpool(idx,:);
  53. %==============
  54. % copy gens in th gpool according to the probility ratio
  55. % >1.1 copy twice
  56. % >=0.9 copy once
  57. % ;0.9 remove
  58. [csort,ridx]=sort(cost);
  59. % sort from small to big.
  60. csum=sum(csort);
  61. caverage=csum/ngpool;
  62. cprobilities=caverage./csort;
  63. copynumbers=0;removenumbers=0;
  64. for i=1:ngpool,
  65.     if cprobilities(i) >1.1
  66.              copynumbers=copynumbers+1;
  67.                     end
  68.            if cprobilities(i) <0.9
  69.                    removenumbers=removenumbers+1;
  70.                            end
  71.                 end
  72.    copygpool=min(copynumbers,removenumbers);
  73.                for i=1:copygpool
  74.                   for j=ngpool:-1:2*i+2 gpool(j,:)=gpool(j-1,:);
  75.             end
  76.                    gpool(2*i+1,:)=gpool(i,:);
  77.           end
  78.                  if copygpool==0
  79.                        gpool(ngpool,:)=gpool(1,:);
  80.                   end
  81. %=========
  82. %when genaration is more than 50,or the patterns in a couple are too close,do mutation
  83. for i=1:ngpool/2
  84.         %
  85. sameidx=[gpool(2*i-1,:)==gpool(2*i,:)];
  86. diffidx=find(sameidx==0);
  87.            if length(diffidx)<=2
  88.                 gpool(2*i,:)=[1 randomize([2:12]')' 1];
  89.                            end
  90.                                end
  91. %===========
  92. %cross gens in couples
  93.            for i=1:ngpool/2
  94.                   [gpool(2*i-1,:),gpool(2*i,:)]=crossgens(gpool(2*i-1,:),gpool(2*i,:));
  95.        end

  96.         for i=1:ngpool,
  97.               cost(i)=sum(diag(distance(gpool(i,:)',rshift(gpool(i,:))')));
  98.        end
  99. % record current best solution
  100. [costmin,idx]=min(cost);
  101. tourmin=gpool(idx,:);
  102. disp([num2str(increase) 'minmum trip length = ' num2str(costmin)])
  103. end  

  104. disp(['cost function evaluation: ' int2str(increase) ' times!'])
  105. disp(['n:' int2str(resultincrease)])
  106. disp(['minmum trip length = ' num2str(resultcost)])
  107. disp('optimum tour = ')
  108. disp(num2str(resulttour))
  109. %====================================================
  110. function B=randomize(A,rowcol)
  111. % Usage: B=randomize(A,rowcol)
  112. % randomize row orders or column orders of A matrix
  113. % rowcol: if =0 or omitted, row order (default)
  114. % if = 1, column order

  115. rand('state',sum(100*clock))
  116. if nargin == 1,
  117.         rowcol=0;
  118. end
  119.          if rowcol==0,
  120.               [m,n]=size(A);
  121.               p=rand(m,1);
  122.               [p1,I]=sort(p);
  123.               B=A(I,:);
  124. elseif rowcol==1,
  125.           Ap=A';
  126.           [m,n]=size(Ap);
  127.           p=rand(m,1);
  128.           [p1,I]=sort(p);
  129.           B=Ap(I,:)';
  130. end
  131. %=====================================================
  132. function y=rshift(x,dir)
  133. % Usage: y=rshift(x,dir)
  134. % rotate x vector to right (down) by 1 if dir = 0 (default)
  135. % or rotate x to left (up) by 1 if dir = 1
  136. if nargin ;2, dir=0; end
  137. [m,n]=size(x);
  138. if m>1,
  139. if n == 1,
  140.     col=1;
  141. elseif n>1,
  142.     error('x must be a vector! break');
  143. end % x is a column vectorelseif m == 1,
  144. end
  145. end
  146. if n == 1, y=x;
  147. return
  148. elseif n>1,
  149.      col=0; % x is a row vector endend
  150. if dir==1, % rotate left or up
  151.        if col==0, % row vector, rotate left
  152.              y = [x(2:n) x(1)];
  153.        elseif col==1,
  154.              y = [x(2:n); x(1)]; % rotate up
  155. end
  156.    elseif dir==0, % default rotate right or down
  157.               if col==0,
  158.                     y = [x(n) x(1:n-1)];
  159.              elseif col==1 % column vector
  160.                        y = [x(n); x(1:n-1)];
  161.                    end
  162. end
  163. end            
  164. %==================================================
  165. function [L1,L2]=crossgens(X1,X2)
  166. % Usage:[L1,L2]=crossgens(X1,X2)
  167. s=randomize([2:12]')';
  168. n1=min(s(1),s(11));n2=max(s(1),s(11));
  169. X3=X1;X4=X2;
  170. for i=n1:n2,
  171.                 for j=1:13,
  172.                      if X2(i)==X3(j),
  173.                           X3(j)=0;
  174.                              end
  175.                   if X1(i)==X4(j),                          X4(j)=0;
  176.                end
  177.            end
  178.         end
  179.    j=13;k=13;
  180.     for i=12:-1:2,
  181.           if X3(i)~=0,
  182.                j=j-1;
  183.                  t=X3(j);X3(j)=X3(i);X3(i)=t;
  184.                end
  185.                     if X4(i)~=0,
  186.                            k=k-1;
  187.                       t=X4(k);X4(k)=X4(i);X4(i)=t;
  188.                    end
  189.                end
  190.            for i=n1:n2
  191.               X3(2+i-n1)=X2(i);
  192.               X4(2+i-n1)=X1(i);
  193.            end
  194. L1=X3;L2=X4;
  195. end
  196. end
  197. end
复制代码

作者: tianpengyun    时间: 2013-8-20 18:48
madio 发表于 2013-8-20 18:06

请问这个要怎么运行?
作者: madio    时间: 2013-8-20 19:45
tianpengyun 发表于 2013-8-20 18:48
请问这个要怎么运行?

这里面定义了一些函数,需要调用这些函数
作者: tianpengyun    时间: 2013-8-20 20:10
不是的,结果是这个
1minmum trip length = 46
cost function evaluation: 1 times!
n:1
minmum trip length = 100000
optimum tour =
0  0  0  0  0
代表什么意思,感觉好像不对


作者: gt93    时间: 2013-8-21 19:00
  1. function glf
  2. v=[51,67;37,84;41,94;2,99;18,54;4,50;24,42;25,38;13,40;7,64;22,60;25,62;18,40;41,26];
  3. n=size(v,1);
  4. for i=1:n
  5.     for j=1:n
  6.         W(i,j)=sqrt((v(i,1)-v(j,1))^2+(v(i,2)-v(j,2))^2);
  7.     end
  8. end
  9. [C d1]=glffunction(W)

  10. function [C d1]=glffunction(d)
  11. % d coordinate matrix.
  12. % C sequence of vertices traversed in the Hamilton cycle.
  13. % d1 sum of weights in the Hamilton cycle.
  14. n=size(d,2);
  15. C=[1:n 1];
  16. C1=C;
  17. if n>3
  18.     for v=4:(n+1)
  19.         for i=1:(v-3)
  20.             for j=(i+2):(v-1)
  21.                 if d(C(i),C(j))+d(C(i+1),C(j+1))<d(C(i),C(i+1))+d(C(j),C(j+1))
  22.                     C1(1:i)=C(1:i);
  23.                     for k=(i+1):j
  24.                         C1(k)=C(j+i+1-k);
  25.                     end
  26.                     C1((j+1):v)=C((j+1):v);
  27.                 end
  28.             end
  29.         end
  30.     end
  31. elseif n<=3
  32.     if n<=2
  33.         fprintf('It does not exist hamilton circle.');
  34.     else
  35.         fprintf('Any circle is the right answer.');
  36.     end
  37. end
  38. C=C1;
  39. d1=0;
  40. for i=1:n
  41.     d1=d1+d(C(i),C(i+1));
  42. end
复制代码
楼主,你看看这个改良圈算法能不能帮上忙。





欢迎光临 数学建模社区-数学中国 (http://www.madio.net/) Powered by Discuz! X2.5