Hello everyone! In this post, I’m going to demonstrate the power of TSP solver (Travelling Salesman Problem solver) in Matlab for solving an instance with 95 cities. I think this is the easiest method to solve the travelling salesman problems because it takes only 2 minutes to set up and minimum programming skill is required. If you want to download this Matlab code, check the file at the end of this post.
To solve new instances, we just need to update the data on the locations of the cities in this file. You can remove or add cities to this file but don’t change this matrix format.
Let’s see the performance of this TSP solver in solving an instance with 95 cities.
For more videos like this, check my YouTube channel here.
Matlab code
Cities.m
% city x-coordinate y-coordinate
city_location =[1 0.8262 0.7081
2 0.4364 0.5108
3 1.3394 0.8963
4 0.0772 0.4408
5 0.9737 0.2785
6 1.0144 0.5909
7 0.0360 0.5589
8 0.3889 0.4151
9 0.4253 0.6931
10 0.6607 0.1569
11 0.8170 0.7803
12 0.4595 0.2220
13 1.3543 0.8458
14 0.9801 0.5578
15 0.5423 0.2251
16 0.6098 0.4689
17 0.4039 0.2918
18 0.6865 0.8605
19 0.8794 0.2835
20 0.4170 0.4546
21 0.7254 0.3622
22 1.0366 0.6892
23 0.5604 0.6681
24 0.5098 0.5728
25 0.4887 0.4451
26 1.0372 0.6505
27 1.0859 0.4751
28 0.1939 0.5533
29 1.0224 0.5410
30 1.0608 0.2639
31 1.0895 0.4802
32 1.2632 0.7448
33 0.9505 0.3659
34 0.8293 0.1964
35 0.2881 0.7257
36 1.2765 0.5436
37 0.1347 0.4889
38 0.7276 0.4553
39 0.1104 0.8924
40 0.3141 0.7307
41 0.9767 0.4790
42 0.4122 0.6522
43 0.1243 0.9597
44 0.8111 0.8375
45 1.0907 0.3591
46 0.4707 0.4196
47 1.3202 0.7990
48 1.0482 0.4825
49 0.4306 0.8337
50 1.3083 0.0921
51 0.3239 0.8318
52 0.4189 0.4308
53 1.2554 0.5347
54 1.1625 0.2308
55 1.1141 0.3047
56 0.7752 0.1563
57 0.2757 0.4364
58 0.7467 0.5833
59 0.9308 0.3728
60 0.9281 0.1572
61 0.4133 0.7987
62 0.3645 0.4795
63 0.6923 0.1625
64 0.8672 0.6948
65 0.6591 0.3720
66 0.3416 0.8148
67 0.3965 0.4103
68 0.2179 0.7020
69 0.8506 0.6116
70 1.1035 0.6327
71 0.6985 0.8571
72 0.2588 0.8284
73 1.0140 0.4459
74 1.0074 0.4604
75 0.3779 0.9040
76 1.0723 0.4144
77 1.1030 0.4771
78 0.7508 0.4176
79 1.3797 0.8626
80 0.7591 0.5679
81 0.1724 0.7621
82 0.2598 0.5338
83 0.9135 0.4912
84 0.9468 0.5770
85 0.9760 0.5432
86 0.6780 0.5350
87 0.6479 0.1021
88 0.7744 0.5702
89 1.2359 0.4867
90 1.1169 0.5528
91 0.3591 0.5679
92 0.1724 0.7021
93 0.2598 0.5338
94 0.3135 0.4912
95 0.9468 0.2770];
TSP_Solver.m
clc
close all
clear all
%Original source: https://www.mathworks.com
figure;
load('usborder.mat','x','y','xx','yy');
plot(x,y,'Color','red');
hold on
Cities % import the data stored in matrix: city_location
nStops = max(city_location(:,1));
stopsLon = city_location(:,2);
stopsLat = city_location(:,3);
plot(stopsLon,stopsLat,'*b')
idxs = nchoosek(1:nStops,2);
dist = hypot(stopsLat(idxs(:,1)) - stopsLat(idxs(:,2)), ...
stopsLon(idxs(:,1)) - stopsLon(idxs(:,2)));
lendist = length(dist);
Aeq = spones(1:length(idxs));
beq = nStops;
Aeq = [Aeq;spalloc(nStops,length(idxs),nStops*(nStops-1))];
for ii = 1:nStops
whichIdxs = (idxs == ii);
whichIdxs = sparse(sum(whichIdxs,2));
Aeq(ii+1,:) = whichIdxs';
end
beq = [beq; 2*ones(nStops,1)];
intcon = 1:lendist;
lb = zeros(lendist,1);
ub = ones(lendist,1);
opts = optimoptions('intlinprog','Display','off');
[x_tsp,costopt,exitflag,output]=intlinprog(dist,intcon,[],[],Aeq,beq,lb,ub,opts);
segments = find(x_tsp);
lh = zeros(nStops,1);
lh = updateSalesmanPlot(lh,x_tsp,idxs,stopsLon,stopsLat);
tours = detectSubtours(x_tsp,idxs);
numtours = length(tours);
A = spalloc(0,lendist,0);
b = [];
while numtours > 1
b = [b;zeros(numtours,1)];
A = [A;spalloc(numtours,lendist,nStops)];
for ii = 1:numtours
rowIdx = size(A,1)+1;
subTourIdx = tours{ii};
variations = nchoosek(1:length(subTourIdx),2);
for jj = 1:length(variations)
whichVar = and((sum(idxs==subTourIdx(variations(jj,1)),2)),...
(sum(idxs==subTourIdx(variations(jj,2)),2)));
A(rowIdx,whichVar) = 1;
end
b(rowIdx) = length(subTourIdx)-1;
end
[x_tsp,costopt,exitflag,output] = intlinprog(dist,intcon,A,b,Aeq,beq,lb,ub,opts);
lh = updateSalesmanPlot(lh,x_tsp,idxs,stopsLon,stopsLat);
tours = detectSubtours(x_tsp,idxs);
numtours = length(tours);
end
% Convert from binary solution to string-like solution
[x y]=size(x_tsp);
S=zeros(1,2);
k = 1;
for i = 1:x
if x_tsp(i,1)~= 0
S(k,:)=idxs(i,:);
k = k+1;
end
end
AA=S(1,:);
S(1,:) = [];
[x1 y1]=size(S);
[x2 y2]=size(AA);
for i = 1:x1
[x3 y3]=find(S == AA(1,y2));
if y3 == 1
AA(1,y2+1) = S(x3,2);
else
AA(1,y2+1) = S(x3,1);
end
S(x3,:)=[];
[x1 y1]=size(S);
[x2 y2]=size(AA);
end
Optimal_solution = AA
min_total_distance = costopt
title('Final solution');
hold off
fprintf('The solution gap = ')
disp(output.absolutegap)
fprintf('Note: If the solution gap = 0, the global optimal solution is found')
P/s: If you find the post useful, share it to remember and to help other people as well.
Hello, what is the function for updateSalesmanPlot?