Multiple Traveling Salesman Problem (M-TSP) Genetic Algorithm (GA) Multi-Traveling Salesman Problem

Function [opt_rte, opt_brk, min_dist] = mtsp_ga (xy, dmat, salesmen ,...
Pop_size, num_iter, singles, show_plots, show_waitbar)
% Mtsp_ga Multiple Traveling Salesman Problem (M-TSP) Genetic Algorithm (GA)
% Finds a (near) Optimal Solution to the M-TSP by setting up a GA to search
% For the shortest route (least distance needed for the salesmen to travel
% To each city exactly once and return to their starting locations)
%
% Input:
% XY (float) is an nx2 matrix of city locations, where N is the number of cities
% Dmat (float) is an nxn matrix of city-to-city distances or costs
% Salesmen (scalar integer) is the number of salesmen to visit the cities
% Pop_size (scalar integer) is the size of the population (shocould be divisible by 8)
% Num_iter (scalar integer) is the number of desired iterations for the algorithm to run
% Singles (scalar logical) is a flag that, if zero, forces the salesmen
% Visit at least 2 cities
% Show_plots (scalar logical) is a flag that, if non-zero, generates a couple plots
% Show_waitbar (scalar logical) is a flag that, if non-zero, displays a waitbar
%
% Output:
% Opt_rte (integer array) is the best route found by the algorithm
% Opt_brk (integer array) is the list of Route break points (these specify the indices
% Into the route used to obtain the individual salesman routes)
% Min_dist (scalar float) is the total distance traveled by the salesmen
%
% Route/breakpoint details:
% If there are 10 cities and 3 salesmen, a possible route/break
% Combination might be: RTE = [5 6 9 1 4 2 8 10 3 7], brks = [3 7]
% Taken together, these represent the solution [5 6 9] [1 4 2 8] [10 3 7],
% Which designates the routes for the 3 salesmen as follows:
%. SALESMAN 1 travels from city 5 to 6 to 9 and back to 5
%. Salesman 2 travels from City 1 to 4 to 2 to 8 and back to 1
%. Salesman 3 travels from city 10 to 3 to 7 and back to 10
%
% Example:
% N = 35; dims = 2;
% Xy = 10 * rand (n, dims );
% Salesmen = 5;
% Pop_size = 80;
% Num_iter = 5e3;
% Singles = 0;
% Show_plots = 1;
% Show_waitbar = 1;
% A = reshape (xy, 1, n, dims );
% B = reshape (xy, N, 1, dims );
% Dmat = SQRT (sum (A (ones (n, 1), :,:)-B (:, ones (n, 1 ),:)). ^ 2, 3 ));
% [Opt_rte, opt_brk, min_dist] = mtsp_ga (xy, dmat, salesmen, pop_size ,...
% Num_iter, singles, show_plots, show_waitbar );
%
% Author: Joseph Kirk
% Email: jdkirk630 at gmail dot com
% Release: 1.0
% Release date: 3/4/08

If ~ Nargin % initialize example inputs
N = 50; dims = 2;
XY = 10 * rand (n, dims );
Salesmen = 7;
Pop_size = 80;
Num_iter = 1e4;
Singles = 0;
Show_plots = 1;
Show_waitbar = 1;
Try % compute the distance Matrix
Dmat = distmat (xy); % file exchange contribution #15145
Catch
A = reshape (xy, 1, n, dims );
B = reshape (xy, N, 1, dims );
Dmat = SQRT (sum (A (ones (n, 1),:,:)-B (:, ones (n, 1), :)). ^ 2, 3 ));
End
Else % process the given inputs
N = size (xy, 1 );
[NR, NC] = size (dmat );
If Nr ~ = N | NC ~ = N
Error ('invalid inputs XY and/or dmat ');
End
End

% Verify inputs
Salesmen = max (1, min (n, round (real (salesmen (1 )))));
If ~ Singles
Salesmen = min (floor (n/2), salesmen );
End
Pop_size = max (8 * Ceil (pop_size (1)/8 ));
Num_iter = max (2, round (real (num_iter (1 ))));
Singles = logical (singles (1 ));
Show_plots = logical (show_plots (1 ));
Show_waitbar = logical (show_waitbar (1 ));

% Plot the cities and the distance Matrix
If show_plots
Hfig = figure ('name', 'mtspga ', 'numbertitle', 'off ');
Subplot (2, 2, 1 );
Plot (XY (:, 1), XY (:, 2 ),'.');
For k = 1: N
Text (XY (k, 1), XY (K, 2), [''num2str (k)], 'fontweight ',' B ');
End
Title ('city location ')
Subplot (2, 2 );
Imagesc (dmat );
Colormap (flipud (Gray ));
Title ('distance matrix ');
End

% Initializations for Route break point selection
Num_brks = salesmen-1;
Addto = ones (1, N-2 * num_brks-1 );
For k = 2: num_brks
Addto = cumsum (addto );
End
Cum_prob = cumsum (addto)/sum (addto );

% Initialize the populations
Pop_rte = zeros (pop_size, n); % population of routes
Pop_brk = zeros (pop_size, num_brks); % population of breaks
For k = 1: pop_size
Pop_rte (K, :) = randperm (N );
Pop_brk (K, :) = randbreaks ();
End
Tmp_pop_rte = zeros (8, N );
Tmp_pop_brk = zeros (8, num_brks );
New_pop_rte = zeros (pop_size, N );
New_pop_brk = zeros (pop_size, num_brks );

% Select the colors for the plotted routes
CLR = [1 0 0; 0 0 1; 0.67 0 1; 0 1 0; 1 0.5 0];
If salesmen> 5
CLR = HSV (salesmen );
End

% Run the GA
Global_min = inf;
Total_dist = zeros (1, pop_size );
Dist_history = zeros (1, num_iter );
If show_plots
Pfig = figure ('name', 'current best solution', 'numbertitle', 'off ');
End
If show_waitbar
WB = waitbar (0, 'searching... Please wait ');
End
For iter = 1: num_iter
% Evaluate members of the population
For p = 1: pop_size
D = 0;
P_rte = pop_rte (p ,:);
P_brk = pop_brk (p ,:);
RNG = [[1 p_brk + 1]; [p_brk N] ';
For M = 1: salesmen
D = d + dmat (p_rte (RNG (m, 2), p_rte (RNG (M, 1 )));
For k = RNG (M, 1): RNG (m, 2)-1
D = d + dmat (p_rte (K), p_rte (k + 1 ));
End
End
Total_dist (p) = D;
End

% Find the best route in the population
[Min_dist, Index] = min (total_dist );
Dist_history (ITER) = min_dist;
If min_dist <global_min
Global_min = min_dist;
Opt_rte = pop_rte (index ,:);
Opt_brk = pop_brk (index ,:);
RNG = [[1 opt_brk + 1]; [opt_brk N] ';
% Plot the best route
If show_plots
Figure (pfig); hold off;
For M = 1: salesmen
RTE = [opt_rte (RNG (M, 1): RNG (m, 2) opt_rte (RNG (M, 1)];
Plot (XY (RTE, 1), XY (RTE, 2), '.-', 'color', CLR (M ,:));
Title (sprintf ('total distance = % 1.4f ', min_dist ));
Hold on;
End
End
End

% Genetic algorithm Operators
Rand_grouping = randperm (pop_size );
For P = 8: 8: pop_size
RTEs = pop_rte (rand_grouping (P-7: P ),:);
Brks = pop_brk (rand_grouping (P-7: P ),:);
Dists = total_dist (rand_grouping (P-7: p ));
[Ignore, idx] = min (dists );
Best_of_8_rte = RTEs (idx ,:);
Best_of_8_brks = brks (idx ,:);
Rte_ins_pts = sort (Ceil (N * rand (1, 2 )));
I = rte_ins_pts (1 );
J = rte_ins_pts (2 );
For k = 1:8% generate new solutions
Tmp_pop_rte (K, :) = best_of_8_rte;
Tmp_pop_brk (K, :) = best_of_8_brks;
Switch K
Case 2% flip
Tmp_pop_rte (K, I: J) = fliplr (tmp_pop_rte (K, I: J ));
Case 3% swap
Tmp_pop_rte (K, [I j]) = tmp_pop_rte (K, [J I]);
Case 4% slide
Tmp_pop_rte (K, I: J) = tmp_pop_rte (K, [I + 1: J I]);
Case 5% modify breaks
Tmp_pop_brk (K, :) = randbreaks ();
Case 6% flip, modify breaks
Tmp_pop_rte (K, I: J) = fliplr (tmp_pop_rte (K, I: J ));
Tmp_pop_brk (K, :) = randbreaks ();
Case 7% swap, modify breaks
Tmp_pop_rte (K, [I j]) = tmp_pop_rte (K, [J I]);
Tmp_pop_brk (K, :) = randbreaks ();
Case 8% slide, modify breaks
Tmp_pop_rte (K, I: J) = tmp_pop_rte (K, [I + 1: J I]);
Tmp_pop_brk (K, :) = randbreaks ();
Otherwise % do nothing
End
End
New_pop_rte (P-7: P, :) = tmp_pop_rte;
New_pop_brk (P-7: P, :) = tmp_pop_brk;
End
Pop_rte = new_pop_rte;
Pop_brk = new_pop_brk;
If show_waitbar, waitbar (ITER/num_iter, WB); End
End
If show_waitbar, close (WB); End

% Plot the results
If show_plots
RNG = [[1 opt_brk + 1]; [opt_brk N] ';
Figure (hfig); subplot (2, 2, 3); hold off;
For M = 1: salesmen
RTE = [opt_rte (RNG (M, 1): RNG (m, 2) opt_rte (RNG (M, 1)];
Plot (XY (RTE, 1), XY (RTE, 2), '.-', 'color', CLR (M ,:));
Title (sprintf ('total distance = % 1.4f ', min_dist ));
Hold on;
End
Subplot (2, 2, 4 );
Plot (dist_history, 'R', 'linewidth', 2 );
Title ('best solution history ');
Set (GCA, 'xlim', [1 num_iter], 'ylim', [0 1.1 * max (dist_history)]);
End

% Generate random set of break points
Function breaks = randbreaks ()
If singles % no constraints on breaks
Tmp_brks = randperm (n-1 );
Breaks = sort (tmp_brks (1: num_brks ));
Else % force breaks to be at least two apart
Num_adjust = find (RAND <cum_prob, 1)-1;
Spaces = Ceil (num_brks * rand (1, num_adjust ));
Adjust = zeros (1, num_brks );
For KK = 1: num_brks
Adjust (kk) = sum (spaces = KK );
End
Breaks = 2 * (1: num_brks) + cumsum (adjust );
End
End
End