I. Introduction to TSP

Traveling Salesman Problem, also known as TSP Problem, is one of the famous problems in mathematics. Suppose a traveling businessman wants to visit n cities. He must choose the route he wants to take. The limit of the route is that he can visit each city only once and return to the original city. The goal of path selection is to obtain the minimum path distance among all paths.Mathematical model of TSP

Introduction to genetic algorithm

1 the introduction Genetic algorithm theory2.1 Biological basis of genetic algorithm 2.2 Theoretical basis of genetic algorithm 2.3 Basic concepts of genetic algorithm 2.4 Standard genetic algorithm 2.5 Characteristics of genetic algorithm 2.6 Improvement direction of genetic algorithm 3. Genetic algorithm flow 4 Key Parameters

Three, some source code

function [ globalMin, optRoute, optCentroid] = bdtsp_ga_basic(popSize, numIter, xy, alpha, range  )
% FUNCTION: mclust_ga_basic calculates the minimum of the sum of the max
% distances for several clusters.  In essence, it is a clustering tool
% that minimizes the max distances between each cluster and the farthest
% member of that cluster. Genetic algorithm minimizes the sum of the max 
% distance of each of the centroid (cluster centers) then calculates the
% TSP distance from each centroid to each centroid.
% Inputs: 
%        population size and iterations
%        xy coordinates
%        speed of drone as a factor of blimp speed
%        range of drone
% Outputs:
%        number of centroids
%        assignement of xy coordinates to centroids
%        (note that one or more coordinates will become centroids)
% GA uses uses a tournament approach mutation type genetic algorithm.  
%  Initialize:
%  (1) Calculate the distance from each xy coordinate to every other xy
%  coordinate as distance matrix dmat.
%  Body:
%  (1) Randomly generate populations 
%  (2) Find min-cost of all pop (trials); keep best pop member and plot.
%  (3) Shuffle (reshuffle) pop for a new tournament
%  (4) Sub-group pop into groups of 4.
%      Find the best of the 4; Overwrite worst of 4 from sub-group pop
%  (5) Mutate the best of 4 (winner) in each sub-group
%  (6) Insert best of 4 (winner) and all mutations back into population
%  (7) If iteration budget remains, go to step 4, else terminate.
%  Termination: based on iteration budget.
% 
%  Example of inputs:
%  nStops  =   30;   % Number of delivery stops for blimp-drone
%  popSize =  500;   % Size of the population of trials.
%  numIter = 2500;   % Number of iterations of GA; iteration budget.
%  alpha   =    2;   % Speed of drone as a factor of blimp
%  range   =   10;   % Range of drone (i.e. 10 km)
%  xy      =   50*rand([nStops,2]);
%  bdtsp_ga_basic(popSize, numIter, xy, alpha, range )

% init variables
  minCost=inf; iter=0; 
  costHistory=[];costIteration=[];
% If null arguments to function call 
  showprogress=true;
if nargin < 5  
  showprogress=true;
  nStops=40;  popSize=800; numIter=4000; alpha=2; range=10;
  xy=45*rand([nStops,2]);
end
  
  % initialize distance matrix
    [nPoints , ~]=size(xy);
    meshg = meshgrid(1:nPoints); % calc. distance
    dmat = reshape(sqrt(sum((xy(meshg,:)-xy(meshg',:)).^2,2)),nPoints,nPoints);
  % Initialize the Population
    [n, ~]=  size(xy);
    pop  = zeros(popSize,n);
    popc = zeros(popSize,n);
    pop(1,:) = (1:n);
    for k = 2:popSize
        % random trials of of blimp-drone routing
        pop(k,:) = randperm(n);
    end
        
 % Run the GA
    globalMin     = Inf;
    totalCost     = zeros(1,popSize);
    costHistory   = zeros(1,numIter);
    costIteration = zeros(1,numIter);
    tmpPop        = zeros(4,n);
    newPop        = zeros(popSize,n);
  
 for iter = 1:numIter
        % Evaluate Each Population Member
        for p = 1:popSize
           % first point is always a centroid
             c = pop(p,1); popc(p,:)=zeros([1 n]);
              popc(p,1)=c; % first centroid
              sumDmax = dmat(c,pop(p,2)); % get dmax 
              dMax=sumDmax;
              for k=2:n  % Get max distances
                     % what is the distance of this from centroid
                     d = dmat(c,pop(p,k));
                   if  d>dMax || d>range%.4
                        % assign as a new centroid 
                         if k<n 
                          c  = pop(p,k);
                          dMax = dmat(c,pop(p,k+1));
                          popc(p,k)=c; % current centroid
                         else
                           c = pop(p,k-1);
                           dMax = dmat(c,pop(p,k));  
                           popc(p,k-1)=c; % current centroid
                         end
                         sumDmax = sumDmax + dMax;  % get next dmax 
                   end 
              end
              poptsp = popc(p,:).*(popc(p,:)>0);
              poptsp = poptsp(poptsp>0);
              ln=length(poptsp);
              d2= dmat(poptsp(ln),poptsp(1));
              for k=1:ln-1
                  d2 = d2 + dmat(poptsp(k),poptsp(k+1));
              end
              totalCost(p) = (2*sumDmax/alpha) + d2 ; 
        end
        
        % Find and keep the best layout in the population
        [minCost,index] = min(totalCost);
        costHistory(iter) = minCost;
        costIteration(iter) = iter;
        if minCost < globalMin
            globalMin = minCost;
            optRoute  = pop(index,:);   
            optCentroid = popc(index,:);
            layout      = optRoute(1:n);  % best layout
            centroids   = optCentroid(1:n); 
            if showprogress==true
              call_plot(xy,layout, centroids);
            end
        end
        
        % Genetic Algorithm Operators: tournament mutations
        % randomly reshuffle population for tournament - play different
        % teams on each iteration
        randomOrder = randperm(popSize);

        for p = 4:4:popSize
            % random reshuffle population, group by 4     
            laytes = pop(randomOrder(p-3:p),:); 
            csts   = totalCost(randomOrder(p-3:p));
                % what is the min layout?
            [~,idx] = min(csts); 
                % what is the best layout of 4
            bestOf4Layout = laytes(idx,:);
                % randomly select two layout insertion points and sort
            routeInsertionPoints = sort(ceil(n*rand(1,2)));
                I = routeInsertionPoints(1);
                J = routeInsertionPoints(2);
            for k = 1:4 % Mutate the best layout to get three new layouts; keep orig.
                % a small matrix of 4 rows of best layout
                tmpPop(k,:) = bestOf4Layout;
                switch k
                       % flip segment between two of the departments
                    case 2 % Flip
                        tmpPop(k,I:J) = tmpPop(k,J:-1:I);
                    case 3 % Swap departments
                        tmpPop(k,[I J]) = tmpPop(k,[J I]);
                    case 4 % Slide departments down
                       tmpPop(k,I:J) = tmpPop(k,[I+1:J I]);
                    otherwise % Do Nothing
                end
            end
             % using the original population, create a new population
            newPop(p-3:p,:) = tmpPop;
        end
        pop = newPop;
 end
 
     function call_plot( xy, ~,~)
      subplot(1,2,1)
      plot(xy(:,1), xy(:,2),'k.');
      hold on;
      idx2=(centroids>0).*centroids;
      idx2=idx2(idx2>0);
      cz = centroids; cnt=0;
      for kk= 1:length(cz)
          c=centroids(kk);
          if c>0
              cc=c;
              cnt=cnt+1;
              %RC=[ rand rand rand];
              RC=[  0  0 0];
          else
          ly=layout(kk);
          plot(xy([cc ly],1), xy([cc ly],2),'k');
          plot(xy(ly,1),xy(ly,2),...
           'MarkerSize',15,...
           'MarkerEdgeColor','black',...
           'MarkerFaceColor',RC);
          plot( xy(cc,1), xy(cc,2),'ks');
          end
      end
      idx2=[idx2 idx2(1)];
      plot(xy(idx2,1),xy(idx2,2),'k--+');
      hold off;
      title('Blimp Route/Drone Spokes');
      xlabel('x-coordinate'); ylabel('y-coordinate'); 
      drawnow;
      if iter>0
      subplot(1,2,2)
      dH=costHistory(costHistory>0);
      dI=costIteration(costIteration>0);
      plot(dI, dH,'k-'); 
      title(sprintf('Min-time: %1.1f',minCost));
      xlabel('iter'); ylabel('Cost');
      end
     end
end
    


Copy the code

4. Operation results

Matlab version and references

1 matlab version 2014A

[1] Yang Baoyang, YU Jizhou, Yang Shan. Intelligent Optimization Algorithm and Its MATLAB Example (2nd Edition) [M]. Publishing House of Electronics Industry, 2016. [2] ZHANG Yan, WU Shuigen. MATLAB Optimization Algorithm source code [M]. Tsinghua University Press, 2017.