# /usr/bin/env python
import copy
from math import *
import numpy as np
from numpy.random import default_rng
import matplotlib.pyplot as plt
import time
[docs]def print_progress(step, total):
"""
Prints a progress bar.
Args:
step (int): progress counter
total (int): counter at completion
"""
message = "simulation progress ["
total_bar_length = 60
percentage = int(step / total * 100)
bar_fill = int(step / total * total_bar_length)
for i in range(total_bar_length):
if i < bar_fill:
message += "|"
else:
message += " "
message += "] "+str(percentage)+" %"
if step < total:
print(message, end="\r")
else:
print(message)
[docs]def read_route_from_file(filename):
"""
Read a route from file.
For a route of N cities, the file should contain
the numbers 0, ..., N-1 in some order, each number on its
own line.
The function returns an integer array containing the read numbers.
Args:
filename (str): name of the file to read
Returns:
array: sequence of integers representing the route
"""
f = open(filename)
lines = f.readlines()
f.close()
route = [ ]
for line in lines:
try:
route.append( int(line) )
except:
pass
return np.array(route)
[docs]def read_cities_from_file(filename):
"""
Read the coordinates of cities from a file.
Each line in the file should contain the x and y coordinates of one city.
Args:
filename (str): name of the file to read
Returns:
array: N x 2 array containing xy coordinates for all cities.
"""
f = open(filename)
lines = f.readlines()
f.close()
cities = [ ]
for line in lines:
parts = line.split()
if len(parts) > 0:
cities.append( [ float(parts[0]), float(parts[1]) ] )
return np.array(cities)
[docs]def show_route(cities, route):
"""
Plot the route.
Cities are plotted as a scatter plot in the given coordinates
and the route as city connecting lines. The route is drawn
as a closed loop that returns to the starting city in the end.
Args:
cities (array): city coordinates as N x 2 array
route (array): integer sequence telling the order in which the cities are visited
"""
# routeline will store the coordinates of the cities
# in the order specified by 'route'
routeline = []
for place in route:
routeline.append( cities[place] )
# to show a closed loop
routeline.append( cities[route[0]] )
routeline = np.array(routeline)
plt.plot(cities[:,0], cities[:,1], 'o')
plt.plot(routeline[:,0], routeline[:,1], '-')
plt.show()
[docs]def write_route_to_file(route, filename):
"""
Write the route to a file.
The route is written as a sequence of integers.
These integers are the indices of all the cities
and their order represents the order in which the
cities are visited.
The written file will have each integer on a separate line.
Args:
route (array): integer sequence telling the order in which the cities are visited
filename (str): name of the file to write
"""
writelines = ""
for i in route:
writelines += str(i)+"\n "
f = open(filename,'w')
f.write(writelines)
f.close()
[docs]def create_random_route(cities):
"""
Create a random route.
The route is represented as an array of integers, where
each integer is the index of a city from the list of cities and the order of
these indices represents the order in which the cities are visited.
Since each city is visited exactly once, the route must be a
vector containing the numbers 0, 1, ..., N-1 where N is the number
of cities. This function generates one such vector where the
numbers are randomly ordered.
Args:
cities (array): city coordinates as N x 2 array (only the length N matters)
"""
# an ordered list - this is a valid starting
# configuration as well
Ncities = len(cities)
route = np.linspace(0,Ncities-1,Ncities,dtype='int')
# You can randomly pick the cities one by one
# to form a random ordering, but there is a ready-made
# function to do this
random.shuffle(route)
return np.array( route )
[docs]def calculate_distances(cities):
"""
Calculate the distances between all pairs of cities.
Since the positions of all cities are given as xy coordinates on a plane,
the distances are calculated simply by Pythagorean theorem.
The distances are stored in an N x N array, where the element
(i,j) is the distance between cities i and j.
The function returns this array.
Args:
cities (array): city coordinates as N x 2 array
Returns:
array: city-city distances in a matrix
"""
Ncities = len(cities)
distances = np.array( [[0.0] * Ncities] * Ncities )
for i in range(len(cities)):
for j in range(len(cities)):
dx = cities[i,0]-cities[j,0]
dy = cities[i,1]-cities[j,1]
distances[i,j] = sqrt(dx*dx + dy*dy)
distances[j,i] = distances[i,j]
return distances
[docs]def get_trip_length(route, distances):
"""
Calculates the length of the route for a closed round-trip
that returns to the starting city in the end.
Here route is defined as an integer vector whose elements are
indices of all cities in visit order, and distances is a matrix
storing city-city distances for all city pairs.
More precisely, if route contains the elements
.. math ::
[ i_0, i_1, i_2, ... i_{N-1} ],
and we denote elements of the distances matrix as :math:`d_{i,j}`,
the distance between, say, the first and second city is :math:`d_{i_0, i_1}`
and the total trip length is calculated as
.. math ::
\sum_{n = 0}^{N-2} d_{i_n, i_{n+1}} + d_{i_{N-1}, i_0}.
The last term in the sum is the length of the final stretch of the trip
where the traveler returns to the starting location.
.. note ::
This function is incomplete!
Args:
route (array): integer sequence telling the order in which the cities are visited
distances (array): matrix storing the distances between all cities, as
calculated by :meth:`calculate_distances()`
Returns:
float: the total length of the trip
"""
trip_length = 0.0
# Loop over the pairs of cities on the route
# and add the city-city distance to the total
# trip length.
# Note! This is a closed loop. Remember to add the return trip
# back home from the last city on the route.
return trip_length
[docs]def change_route_randomly(route):
"""
Make random changes to the given route.
A new route is returned.
The original route is not modified.
.. note ::
This function is incomplete!
Args:
route (array): integer sequence telling the order in which the cities are visited
Returns:
array: the new route
"""
test_route = copy.copy(route)
# todo
return test_route
[docs]def Boltzmann_P(new, old, T):
"""
The Boltzmann probaility for accepting a change.
If the new value is better (smaller) than the old one,
the change is always accepted. In this case,
the function may return any number greater or equal to one.
If the new value is worse (larger) than the old one, the function
returns the Boltzmann probability, which is calculated as follows:
Denote the old and new values as :math:`x_\\text{old}` and :math:`x_\\text{new}`,
respectively. The change in the value is then
:math:`\\Delta x = x_\\text{new} - x_\\text{old}`, and the Boltzmann probability
for accepting such a change is
.. math ::
P_{\\text{old} \\to \\text{new}} = e^{-\\frac{\\Delta x}{T} }.
The result is physically analogous to the canonical ensemble in statistical
physics, where the probability that the system at temperature :math:`T`
is in a state :math:`i` with total energy :math:`E_i` is given by
.. math ::
P_i = \\frac{1}{Z} e^{-\\frac{E_i}{k_B T} },
where :math:`Z` is the partition function. Especially the ratio of the probabilites
between two states is
.. math ::
\\frac{P_j}{P_i} = \\frac{ \\frac{1}{Z}e^{-\\frac{E_j}{k_B T} } }{ \\frac{1}{Z}e^{-\\frac{E_i}{k_B T} } }
= e^{-\\frac{E_j - E_i}{k_B T} }
= e^{-\\frac{\Delta E}{k_B T} }.
The differences between the canonical distribution and the values
calculated by this function is that there is no real temperature in this function.
The argument T acts like temperature, but we are not simulating a physical system.
Therefore we set :math:`k_B = 1` since the temperature scale is arbitrary anyway.
.. note ::
This function is incomplete!
Args:
new (float): new value
old (float): old value
T (float): virtual temperature
Returns:
float: the Boltzmann probability
"""
return 0
[docs]def run_annealing_cycle(route, distances,
initial_temperature = 50.0,
trials = 10000,
final_temperature = 1.0):
"""
Perform simulated annealing optimization for the route length.
Annealing is a heat treatment where heating up a physical material,
say, a metal, increases the mobility of the particles that make up
the material. Especially defects become more mobile and may thus diffuse
out of the material or combine to eliminate each other. Thus,
annealing can remove defects from the material. Once the material is
cooled, it will contain very few defects.
Theoretically, annealing allows the system to sample the phase space
so that it can find a state that minimizes the total energy.
Simulated annealing uses the same idea for an abstract optimization
problem. If we treat the length of our route as if it was the energy of
a physical system, we can define a virtual temperature and allow the
system to evolve according to the probabilities given by statistical
physics. When our "temperature" is high, we will find new solutions.
Then, as "temperature" gradually decreases, the solution should
converge to some local optimum.
The annealing algorithm works by repeating the following steps:
* A copy of the current route is made to prevent accidental changes.
* This copy is modified using :meth:`change_route_randomly`.
* The length of the new route is measured using :meth:`get_trip_length`.
* The probability for accepting the new route is calculated using :meth:`Boltzmann_P`.
* The new route is randomly accepted or discarded according to this probability.
* If the new route was accepted, it becomes the current route. Otherwise
it is discarded and the old route remains current.
* If the current route is the shortest one yet found, it is saved as the best route.
.. note ::
This function is incomplete!
Args:
route (array): starting route
distances (array): matrix storing the distances between all cities, as
calculated by :meth:`calculate_distances()`
initial_temperature (float): virtual temperature in the beginning
final_temperature (float): virtual temperature in the end
trials (int): the number of times a new route is generated and tried
Returns:
array: the best route found during the simulation
"""
# the current temperature
temperature = initial_temperature
# store the best route found
best_route = route
shortest_length = get_trip_length(route, distances)
# the length of the current route
current_length = shortest_length
# change of temperature between trials
dt = (final_temperature - initial_temperature) / (trials-1)
successes = 0
# run the algorithm
for i in range(trials):
###################################################
# * make a new route by modifying the current one #
# * calculate the length of the new route #
# * calculate the change in route length, dL #
# * if the new route is shorter, accept it #
# * if the new route is longer, accept it with #
# the Boltzmann probability #
# * if the new route is accepted, make it the #
# current route #
# * test if the new route is the best route, and #
# and store it if it is #
# * if the new route is not accepted, discard it #
###################################################
test_route = copy.copy(route)
test_length = 0
# if the new route is the shortest yet, save it
if test_length < shortest_length:
shortest_length = test_length
best_route = test_route
# lower temperature
temperature += dt
return best_route
[docs]def optimize_route(initial_route, distances):
"""
Searches for as short a route as possible.
The optimization is carried out using several simulated annealing cycles
by repeatedly calling the function :meth:`run_annealing_cycle`.
It is beneficial to run several annealing cycles since there may be
several local optima and one annealing cycle can only converge to one of them.
So in order to search for the global optimum, one should anneal many times
to find different optima and finally choose the best one.
Args:
initial_route (array): starting route
distances (array): matrix storing the distances between all cities, as
calculated by :meth:`calculate_distances()`
.. note ::
Modify this function to adjust performance!
Returns:
array: the best route found during the simulation
"""
# Let's copy the initial route so we don't accidentally change it.
route = copy.copy(initial_route)
# Modify these parameters to find a reliable configuration.
cycles = 2
T_start = 10.0
T_end = 1.0
n_steps = 1000
print_progress(0, cycles)
for heatup in range(cycles):
# run optimization
route = run_annealing_cycle(route, distances, T_start, n_steps, T_end)
print_progress(heatup+1, cycles)
return route
[docs]def main(cityfile, routefile=None):
"""
The main program.
Reads the coordinates of all the cities from the given file and
also tries to read an initial guess for the route from a file.
The city file must exist, because it defines the problem to
be optimized.
The route file is optional and if no valid file is found,
a random route is created as a starting point.
Args:
cityfile (str): name of the file containing city coordinates
routefile (str): name of the optional file containing the initial
route as a sequence of integers (city indices)
"""
# read cities from a file
try:
cities = read_cities_from_file(cityfile)
except:
print("Could not read the cities from file "+str(cityfile))
quit()
# Calculate the distances between all cities and
# store them in a file.
# This way the distances need not be recalculated
# over and over again.
distances = calculate_distances(cities)
# Initialize the route.
#
# If the given file defines a valid route file, it is read.
# If this fails, create a random sequence.
try:
initial_route = read_route_from_file(routefile)
if len(initial_route) != len(cities):
print( "Incompatible route length and city count, ignoring route in "+str(routefile) )
initial_route = create_random_route(cities)
except:
if routefile is not None:
print("Could not read the route from file "+str(routefile) )
initial_route = create_random_route(cities)
# Run the optimization algorithm.
t_start = time.perf_counter()
best_route = optimize_route(initial_route, distances)
t_end = time.perf_counter()
# Print the results
print()
print( "shortest trip distance found: "+str(get_trip_length(best_route, distances)) )
print( "shortest trip found: "+str(best_route) )
print( "run time "+str( round( t_end - t_start, 2) )+" s" )
# Plot the best route and save it.
write_route_to_file(best_route, "best_route.txt")
show_route(cities, best_route)
if __name__ == "__main__":
random = default_rng()
main('A3_10_cities.txt', None)
else:
random = default_rng()