【TSP】 Solving 76 city Travel agent Problem based on Matlab Ant Colony Algorithm

% % % % % % % % % % problem definition: the city location coordinates, calculating distance matrix InitOps % % % % % % % % % = []; [Location,DistMatrix,Ncities,Bestx,Lengx] = pr76init(InitOps); close all; figure (1);
hold on;
minx=min(DistMatrix(:,1));
maxx=max(DistMatrix(:,1));
miny=min(DistMatrix(:,2));
maxy=max(DistMatrix(:,2));
minm=min(minx,miny);
maxm=max(maxx,maxy);
l=(maxm-minm)/10;
for i=1:Ncities
    plot(Location(i,1),Location(i,2),'*b');
    text (Location(i,1)+l,Location(i,2)+l,num2str(i));
end
for i=1:Ncities- 1
    line([Location(Bestx(i),1),Location(Bestx(i+1),1)] , [Location(Bestx(i),2),Location(Bestx(i+1),2)]);end
line([Location(Bestx(1),1),Location(Bestx(Ncities),1)] , [Location(Bestx(1),2),Location(Bestx(Ncities),2)]) ;
grid on,title(['Initial roadmap -',num2str(Lengx)]),xlabel(X-coordinate),ylabel('ordinate');
legend('City location'); hold off ; % initialize random generator state rand('state',sum(100*clock)); % = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = % use the nearest neighbor method to construct an initial, calculated according to the information department of initial p = zeros (1,Ncities+1);
    p(1)=round(Ncities*rand+0.5); % p Store the numbers of all cities found so far1);
    count=2;
	whilecount <= Ncities NNdist= inf ; %NNdist Stores the distance between the found city and the city with the shortest distance pp= I; % I stores the current city id. Pp Stores the current found city idfor j= 1: Ncities
          	if (DistMatrix(i, j) < NNdist) & (j~=i) & ((j~=p) == ones(1,length(p))) % The target city is required to be -- short distance, and can not be the current city, also can not be the city that has been passed before NNdist= DistMatrix(I, j); pp= j ;end           
     	end
     	p(count)=pp; 
        i= pp ;
     	count= count + 1 ;
	end
    p=AcatspEval(p,DistMatrix,Ncities);
    len=p(1,Ncities+1);
	Q0=1/(Ncities*len); %%%%%%%%%% Set system parameters %%%%%%%%%% MaxNc=5000; % maximum algebra A=1; % pheromone factor B=2; % heuristic information factor P1=0.1; % initial value of local volatilization coefficient P2=0.1; % global volatility coefficient initial value R0=0.9; % selection probability M=10; % number of ants %%%%%%%%%% Initializing Pheromone, heuristic information matrix, determining initial position of ant and permit matrix %%%%%%%%%% Pheromone=Q0*ones(Ncities,Ncities); % Initial pheromone matrix; Heuristic=1./DistMatrix; % Temp=ones(1,Ncities);
Heuristic=1./ (1./Heuristic+diag(Temp));
RandL=round(rand(M,1)*Ncities+0.5); % Ant initial location Alocation0=zeros(M,Ncities+1); Deposit % M +1Alocation0(:,1)=RandL;
Allow0=repmat(1:Ncities,M,1); % Allows access to the city matrix initializationfor Ak=1:M
    Allow0(Ak,RandL(Ak))=0; End %%%%%%%%%% Run parameter initialization %%%%%%%%%% Nc=1; % first generation Lbestdis=inf; Cbestdis=inf; Fnewbest=0; Alocation=Alocation0; Allow=Allow0; Allow=Allow0; % allows the matrix to assign the initial value t1=clock; % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % ant colony algorithm initialization program end % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % of ant colony algorithm is the main loop start % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %while(Nc<=MaxNc) % M ants select Ncitiesfor Cityi=2:Ncities+1
          if Cityi<Ncities+1 
              for Ak=1:M
                  i=Alocation(Ak,Cityi- 1); J = % of the current city Select_for_aca (R0, Ak, I, to Allow A, B, Pheromone and Heuristic); % Select the next city based on Pijj
                  Alocation(Ak,Cityi)=j;
                  Allow(Ak,j)=0; % update permit matrix Pheromone(I,j)=(1-P1)*Pheromone(i,j)+P1*Q0; Pheromone(j, I)= Pheromone(I,j); endelse% returns to the starting cityfor Ak=1:M
                  i=Alocation(Ak,Cityi- 1); % Current city j=Alocation(Ak,1);
                  Pheromone(i,j)=(1-P1)*Pheromone(i,j)+P1*Q0; Pheromone(j, I)=Pheromone(I,j); end end end function D = dists(X1,X2,p,e) %DISTS Calculates distances between vectors of points. % D = dists(X1,X2,p,e) % X1 = n x d matrix of n d-dimensional points % X2 = m x d matrix of m d-dimensional points % D = n x m matrix of distances % = (n- 1) x 1 vector of distances between X1 points, if X2 = []
%     p = 2, Euclidean (default): D(i,j) = sqrt(sum((X1(i,:) - X2(j,:))^2% =))1, rectilinear: D(i,j) = sum(abs(X1(i,:) - X2(j,:))
%       = Inf, Chebychev dist: D(i,j) = max(abs(X1(i,:) - X2(j,:))
%       = (1 Inf), lp norm: D(i,j) = sum(abs(X1(i,:) - X2(j,:))^p)^(1/p)
%       = 'rad', great circle distance in radians of a sphere
%         (where X1 and X2 are decimal degree latitudes and longitudes)  
%       = 'mi' or 'sm', great circle distance in statute miles on the earth
%       = 'km', great circle distance in kilometers on the earth
%     e = epsilon for hyperboloid approximation gradient estimation
%       = 0 (default); no error checking if any non-empty 'e' input
%      ~= 0 => general lp used for rect., Cheb., and p outside [1.2] 
%
% Great circle distances are calculated using the Haversine Formula (from R.W.
% Sinnott, "Virtues of the Haversine", Sky and Telescope, vol. 68, no. 2.1984
% p. 159)

% Copyright (c) 1998 by Michael G. Kay
% Matlog Version 1.0 Apr- 3- 98.

% Input Error Checking *******************************************************
if nargin == 4 & ~isempty(e)		% No error checking is 'e' input
   [n,d] = size(X1);
   m = size(X2,1);
else					% Do error checking
   error(nargchk(1.4,nargin));
   e = 0;				% 'e' default
   [n,d] = size(X1);
   if n == 0 | ~isreal(X1)
      error('X1 must be non-empty real matrix');
   end
   if nargin < 2 | isempty(X2)		% Calc intra-seq dist b/w X1 points
      m = 0;				% X2 default
      if n < 2
	 error('X1 must have more than 1 point');
      end
   else					% Calc dist b/w X1 and X2 points
      [m,dX2] = size(X2);
      if m == 0 | ~isreal(X2)
	 error('X2 must be non-empty real matrix');
      end
      if d ~= dX2
	 error('Rows of X1 and X2 must have same dimensions');
      end
   end
   if nargin < 3 | isempty(p)
      p = 2;				% 'p' default
   elseif ~ischar(p)			% lp distances
      if length(p(:)) ~= 1 | ~isreal(p)
	 error(' ' 'p' ' must be a real scalar number');
      end
   elseif ischar(p)			% Great circle distances
      p = lower(p);
      if d ~= 2
	 error('Points must be 2-dimensional for great-circle distances');
      end
      if ~any(strcmp(p,{'rad'.'mi'.'sm'.'km'}))
	 error(' ' 'p' ' must be either ''rad,'' ''mi,'' ''sm,'' or ''km'"); end else error('''p'' not valid value');
   end
end
% End (Input Error Checking) ***********************************************

% Interchange if X2 is the only 1 point
intrchg = 0;
if n > 1 & m == 1tmp = X2; X2 = X1; X1 = tmp; m = n; n =1;
   intrchg = 1;
end

% 1-dimensional points
if d == 1
   if e == 0
      if m ~= 0
	 D = abs(X1(:,ones(1,m)) - X2(:,ones(1,n))'); else D = abs(X1(1:n-1) - X1(2:n))';	% X1 intra-seq. dist.
      end
   else
      if m ~= 0
	 D = sqrt((X1(:,ones(1,m)) - X2(:,ones(1,n))').^2 + e); else D = sqrt((X1(1:n-1) - X1(2:n)).^2 + e)';
      end
   end

% X1 only 1 point or intra-seq dist   
elseif n == 1 | m == 0
   if n == 1				% Expand X1 to match X2
      X1 = X1(ones(1,m),:);
      n = m;
   else					% X1 intra-seq. dist.
      X2 = X1(2:n,:);		               % X2 = ending points
      n = n - 1;
      X1 = X1(1:n,:); 	                       % X1 = beginning points
   end
   if p == 2				% Euclidean distance
      D = sqrt(sum(((X1 - X2).^2 + e)')); elseif ischar(p) % Great-circle distance X1 = pi*X1/180; X2 = pi*X2/180; D = 2*asin(min(1,sqrt(sin((X1(:,1) - X2(:,1))/2).^2 + ... cos(X1(:,1)).*cos(X2(:,1)).* ... sin((X1(:,2) - X2(:,2))/2).^2)))';
   elseif p == 1 & e == 0		% Rectilinear distance
      D = sum(abs(X1 - X2)'); elseif (p >= 1 & p <= 2) | (e ~= 0 & p > 0) % General lp distance D = sum((((X1 - X2).^2 + e).^(p/2))'.) ^ (1/p);
   elseif p == Inf & e == 0		% Chebychev distance
      D = max(abs(X1 - X2)'); else % Otherwise D = zeros(1,n); for j = 1:n D(j) = norm(X1(j,:) - X2(j,:),p); end endCopy the code

Version: 2014 a