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)
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.