An algorithm for the open travelling salesman problem

Question

I have come up with a dynamic programming approach to the travelling salesman problem, but I don't have the background knowledge of whether this type of program is in use, or programming skills to figure out how good a program it is. These are the steps of my program:

1. For each city, connect it to the city that is closest to it (unless there is already a connection with that city). This step will result in one or more groups of cities, where a group is a set that is connected. (Marked in the diagram below as A, B, C, D.)

2. Repeat step 1, now at the level of groups. Each group has 2 active ends - the city at either end of the chain (Blue in the diagram above). To find the distance of a group to another group, take the minimum of the four possible combinations for joining the active ends to each other. Choose a random group. Compare the distance of that group to each other group, choose the one which is closest, and connect them to form a larger group (this means two cities will become inactive. Choose another random group, connect it to the one which is closest. Keep repeating until there is only one group.

This algorithm seems to be quite light on computational power. How good is it in generating close-to optimal solutions?

Answer 1

I wrote a small python script to do some simulations. On the top you can input the number of points as well as the number of iterations. I also changed the score to approximate / best, which lets us easily read of how far away your solution is to the reference. Note that this implementation is not necessarily very efficient, I just tried to implement what I understood from your description.

Some results: Up to N = 3 nodes, it unsurprisingly always finds the optimal solution. For N = 4 your solution seems to be always around 1.2% longer. For N = 6 it seems to be at arund 4.5% and for N = 7 around 7%.

Try it online!

#Number of points
N = 4
#Number of Iterations
iterations = 1000


import numpy as np
import scipy.spatial
import itertools


best_distances = []
approximate_distances = []
for k in range(iterations):
    #generate n uniformly random points in [0,1]x[0,1]
    p = np.random.rand(N,2) #use randn for normal distribution
    #p[:,0]*=0.1 #uncomment this for "squeezing" all the points into a long rectangle
    #compute adjecency atrix
    A = scipy.spatial.distance.squareform(scipy.spatial.distance.pdist(p, metric='euclidean'))
    
    #compute exact solution by brute force
    best_distance =  A.sum()#path cannot be longer than that
    best_permutation = None
    for permutation in itertools.permutations(list(range(N))):
        total_distance = sum([A[i][j] for (i,j) in zip(permutation[1:],permutation[:-1])])
        if total_distance < best_distance:
            best_distance = total_distance
            best_permutation = permutation
    
    #print("coordinates:\n", p)
    #print("shortest path: {}; length: {}".format(best_permutation, best_distance))
    best_distances.append(best_distance)
    
    
    #use suggested algorithm:
    class Group:
        def __init__(self, start, end=None, path=[], length=None):
            self.start = start
            self.end = start if end is None else end
            self.path = path
            self.length = length if length is not None else self.computePathDist()
    
        def getTotalPath(self):
            return [self.start] + self.path + [self.end]
                
        def computePathDist(self):
            tp = self.getTotalPath
            return sum([A[i][j] for (i,j) in zip(tp[1:], tp[:-1])])
        
        def reversed(self):
            rev = Group(start=self.end, end=self.start, path=list(reversed(self.path)), length=self.length)
            return rev
        
        def _append(self, g):#not inplace
            self_end = [] if self.start == self.end else [self.end] #edgecase
            g_start = [] if g.start == g.end else [g.start]#edgecase
            total_path = self.path + self_end + g_start + g.path
            total_length = self.length + A[self.end,][g.start] + g.length
            merged = Group(start=self.start, end=g.end, path=total_path, length=total_length)
            return merged
            
        def mergeWith(self, g, inplace=True):#inplace
            self_rev = self.reversed()
            g_rev = g.reversed()
            best_distance = A.sum()
            best_order = (None, None)
            for (a,b) in [(self,g), (self, g_rev), (self_rev, g), (self_rev, g_rev)]:
                if A[a.end][b.start] < best_distance:
                    best_distance = A[a.end][b.start] 
                    best_order = (a,b)
            (a,b) = best_order
            merged = a._append(b)
            if inplace:
                self.start = merged.start
                self.end = merged.end
                self.path = merged.path
                self.length = merged.length
                return self
            else:
                return merged
            
            
            
    groups = [Group(start=k, end=k, path=[], length=0) for k in range(N)]
    while len(groups) > 1:
        current = groups.pop()
        best_distance = A.sum()
        best_neighbor = None
        for g in groups:
            temp_merge = g.mergeWith(current, inplace=False)
            if temp_merge.length < best_distance:
                best_distance = temp_merge.length
                best_neighbor = g
        best_neighbor.mergeWith(current, inplace=True)
            
    #print("approximate shortest path: {}; length: {}".format(groups[0].getTotalPath(), groups[0].length))    
    approximate_distances.append(groups[0].length)
        
        
scores = np.array([a/b for b,a in zip(best_distances, approximate_distances)])  
print("average score: {} \t standard deviation: {}".format(scores.mean(), scores.std()))

Try it online!

Answer 2

I coded your heuristic in MATLAB (because I already happened to have a bunch of other TSP heuristics coded) and tested it on 50 random instances on the unit square with up to 200 nodes. I compared it to several other heuristics, as well as the Held-Karp lower bound. (In a handful of cases, the HK lower bound is the optimal objective function value, but in most cases it is strictly smaller.)

Here are the tour lengths plotted against the number of nodes. (I called your heuristic "Successive Grouping".)

For each heuristic, here is the average ratio between the tour length from the heuristic and the HK lower bound):

Successive Grouping            : 1.5890
Nearest Neighbor               : 1.2096
Nearest Insertion              : 1.1899
Farthest Insertion             : 1.0731
Cheapest Insertion             : 1.1558
Convex Hull                    : 1.0778
Minimum Spanning Tree          : 1.3078
Christofides                   : 1.1100

As you can see, your heuristic is not great. :) It's an interesting idea, though. I think part of the problem stems from the fact that after you connect all the pieces, you're left with an additional edge that must be added (from the first to the last node), and this edge is not considered at all during the heuristic, so it could be quite long. (In your example, it's quite long.)

By the way, your heuristic is not really a dynamic programming algorithm. It's just iterative.

Here's my MATLAB script, for the record:

% number of trials
T = 50;

% maximum number of nodes
maxN = 200;

% output arrays
N_arr = zeros(1, T);
SG_z = zeros(1, T);
NN_z = zeros(1, T);
NI_z = zeros(1, T);
FI_z = zeros(1, T);
CI_z = zeros(1, T);
CH_z = zeros(1, T);
MST_z = zeros(1, T);
C_z = zeros(1, T);
HK_z = zeros(1, T);
HK_code = zeros(1, T);

for t = 1:T

    if mod(t, 10) == 0
        display(t)
    end

    % generate instance
    N = randi(maxN);
    N_arr(t) = N;
    x = rand(N, 1);
    y = rand(N, 1);
    c = squareform(pdist([x y]));

    % solve with heuristics
    [~, SG_z(t)] = TSPSuccessiveGrouping(c);
    [~, NN_z(t)] = TSPNearestNeighbor(c, 1);
    [~, NI_z(t)] = TSPNearestFarthestInsertion(c, 1, true);
    [~, FI_z(t)] = TSPNearestFarthestInsertion(c, 1, false);
    [~, CI_z(t)] = TSPCheapestInsertion(c, 1);
    [~, CH_z(t)] = TSPConvexHull(c, x, y);
    [~, MST_z(t)] = TSPMST(c);
    [~, C_z(t)] = TSPChristofides(c);

    % solve with Held-Karp
    [HK_z(t), ~, HK_code(t), ~] = TSPHeldKarpSO(c, ones(1, N), NN_z(t));

end

% plot
close all
figure
[N_arr_sorted, I] = sort(N_arr);
plot(N_arr_sorted, SG_z(I), '-o');
hold on
plot(N_arr_sorted, NN_z(I), '-+');
plot(N_arr_sorted, NI_z(I), '-*');
plot(N_arr_sorted, FI_z(I), '-.');
plot(N_arr_sorted, CI_z(I), '-x');
plot(N_arr_sorted, CH_z(I), '-s');
plot(N_arr_sorted, MST_z(I), '-d');
plot(N_arr_sorted, C_z(I), '-^');
plot(N_arr_sorted, HK_z(I), '-v');
legend({'Successive Grouping', 'Nearest Neighbor', 'Nearest Insertion', 'Farthest Insertion', 'Cheapest Insertion', 'Convex Hull', 'MST', 'Christofides', 'Held-Karp Lower Bound'}, 'Location', 'southeast');
xlabel('N');
ylabel('Tour Length');

fprintf('Successive Grouping            : %.4f\n', mean(SG_z ./ HK_z));
fprintf('Nearest Neighbor               : %.4f\n', mean(NN_z ./ HK_z));
fprintf('Nearest Insertion              : %.4f\n', mean(NI_z ./ HK_z));
fprintf('Farthest Insertion             : %.4f\n', mean(FI_z ./ HK_z));
fprintf('Cheapest Insertion             : %.4f\n', mean(CI_z ./ HK_z));
fprintf('Convex Hull                    : %.4f\n', mean(CH_z ./ HK_z));
fprintf('Minimum Spanning Tree          : %.4f\n', mean(MST_z ./ HK_z));
fprintf('Christofides                   : %.4f\n', mean(C_z ./ HK_z));
fprintf('Held-Karp Lower Bound          : %.4f\n', mean(HK_z ./ HK_z));