Why is this simulated annealing algorithm applied to the TSP not converging?

Question

For an exercise I have to find the optimal solution to this Travelling Salesman Problem, only difference is that the salesman may not visit the first coordinate twice. The optimum should lie around 1200 which is given. I don't see why this wouldn't converge. Also an important note is that the manhattan distance metric is used instead of the euclidian distance.

def shuffle_coords(coords):
     x0 = coords[0]
     coords = shuffle(coords[1:])
     return np.insert(coords,0,x0,axis = 0)

def distance(x,y):
    return abs(x[1] - y[1]) + abs(x[0] - y[0])

Here the coordinates are shuffled randomly.

def shuffle(array):
    df = pd.DataFrame(array)
    return df.sample(len(array)).to_numpy() 


def path_distance(path):
    dist     = []
    for i in range(1,len(path)):
        dist.append(distance(path[i],path[i-1]))
    return np.sum(dist)

Here the Simulated Annealing algorithm is initialised.

def SA_distance(path, T_0,T_min, alpha):
    T = T_0
    dist = path_distance(path)
    
    while T >  T_min:
        new_path = gen_subtour(path)
        diffF = path_distance(new_path) - dist
        if diffF < 0:
            path = new_path
            dist = path_distance(path)
 
        elif np.exp(-(diffF/T)) > random.uniform(0,1):
            path = new_path
            dist = path_distance(path)
      
        T = T * alpha
        print(dist,T)
    return dist,path

Here the random subtour is generated while keeping x0 in place.

def gen_subtour(path):
    subset = shuffle(np.delete(path,0,axis =0))
    subset = shuffle(path)
    if random.uniform(0,1) < 0.5:
        subset = np.flipud(subset)
    else:
        j = random.randint(1,(len(subset)-1))
        p = subset[j-1]
        q = subset[j]
        subset = np.delete(subset,[j-1,j],axis = 0)
        subset = np.insert(subset,0,p,axis = 0)
        subset = np.insert(subset,len(subset),q,axis = 0)
 
    return np.insert(subset,0,path[0],axis = 0)

def main():
    T_0     = 12
    T_min   = 10**-9
    alpha   =  0.999
    coords = np.array([[375, 375],[161, 190], [186, 169],[185, 124],
                       [122, 104],[109, 258], [55, 153] ,[120, 49],
                       [39, 85]  ,[59, 250] , [17, 310] ,[179, 265],
                       [184, 198]])
    path , distance = SA_distance(coords,T_0,T_min,alpha)

Answer 1

I have tested, and my first instinct is that the gen_subtour() is not working well, maybe it is changing the route with too many steps.

I would try different versions and check how the effect it. The SA scheme seem to be working well, it is the proposals where I think the error is.

Anyway, here some code that hopefully will help you test better.

I use the pdist to pre-calculate manhatten distance, i.e;

import numpy as np
from scipy.spatial.distance import pdist, cdist, squareform

coords = np.array([[375, 375],[161, 190], [186, 169],[185, 124],
                   [122, 104],[109, 258], [55, 153] ,[120, 49],
                   [39, 85]  ,[59, 250] , [17, 310] ,[179, 265],
                   [184, 198]])

Y = pdist(coords, 'cityblock')

distance_matrix = squareform(Y)
nodes_count = coords.shape[0]

And define start with

def random_start():
    """
        Random start, returns a state
    """
    a = np.arange(0,nodes_count)
    np.random.shuffle(a)
    return a

The objective function, where we have the version that we don't return to the origin;

def objective_function( route ):
    # uncomment when testing new/modify neighbors
    # assert check_all_nodes_visited(route)

    return np.sum( distance_matrix[route[1:],route[:-1]] )

I have here 3 kinds of proposals, based on a route;

def random_swap( route ):
    """
        Random Swap - a Naive neighbour function

        Will only work for small instances of the problem
    """
    route_copy = route.copy()

    random_indici = np.random.choice( route , 2, replace = False)
    route_copy[ random_indici[0] ] = route[ random_indici[1] ]
    route_copy[ random_indici[1] ] = route[ random_indici[0] ]

    return route_copy

def vertex_insert( route, nodes=1 ):
    """
        Vertex Insert Neighbour, inspired by

        http://www.sciencedirect.com/science/article/pii/S1568494611000573
    """
    route_copy = route.copy()
    random_indici = np.random.choice( route , 2, replace = False)
    index_of_point_to_reroute = random_indici[0]
    value_of_point_to_reroute = route[ random_indici[0] ]
    index_of_new_place = random_indici[1]
    route_copy = np.delete(route_copy, index_of_point_to_reroute)
    route_copy = np.insert(route_copy, index_of_new_place, values=value_of_point_to_reroute)
    return route_copy

def block_reverse( route, nodes=1 ):
    """
        Block Reverse Neighbour, inspired by

        http://www.sciencedirect.com/science/article/pii/S1568494611000573

        Note that this is a random 2-opt operation.
    """
    route_copy = route.copy()
    random_indici = np.random.choice( route , 2, replace = False)
    index_of_cut_left = np.min(random_indici)
    index_of_cut_right = np.max(random_indici)
    route_copy[ index_of_cut_left:index_of_cut_right ] = np.flip(route_copy[ index_of_cut_left:index_of_cut_right ])

    return route_copy

Optionally, you can to a 2-opt round after the SA, to make sure there are no crossing.

def swap_for_2opt( route, i, k):
    """
        Helper for 2-opt search
    """
    route_copy = route.copy()
    index_of_cut_left = i
    index_of_cut_right = k
    route_copy[ index_of_cut_left:index_of_cut_right ] = np.flip(route_copy[ index_of_cut_left:index_of_cut_right ])

    return route_copy

def local_search_2opt( route ):
    """
        Local Optimum with 2-opt

        https://en.wikipedia.org/wiki/2-opt

    """
    steps_since_improved = 0
    still_improving = True

    route = route.copy()

    while still_improving :
        for i in range( route.size - 1 ):
            for k in np.arange( i + 1, route.size ):
                alt_route = swap_for_2opt(route, i, k)

                if objective_function(alt_route) < objective_function(route):
                    route = alt_route.copy()
                    steps_since_improved = 0

            steps_since_improved += 1

            if steps_since_improved > route.size + 1:
                still_improving = False
                break

    return route

And use frigidum for SA

import frigidum


local_opt = frigidum.sa(random_start=random_start,
           objective_function=objective_function,
           neighbours=[random_swap, vertex_insert, block_reverse],
           copy_state=frigidum.annealing.naked,
           T_start=10**5,
           alpha=.95,
           T_stop=0.001,
           repeats=10**2,
           post_annealing = local_search_2opt)

The returned route is almost always 1145.

I have posted general hints and tips on frigidum homepage.