A classy approach to solving Traveling Salesman Problems effectively with Python | by Carlos J. Uribe

Here we won’t start from scratch. As stated earlier, we already developed the code that builds a Pyomo model of the TSP and solves it in sprint 3. And trust me, that was the hardest part. Now, we have the easier task of organizing what we did in a way that makes it general, hiding the details while keeping the essential elements visible. In a sense, we want the optimizer to look like a “magic box” that even users not familiar with math modeling are able to use to solve their TSP problems intuitively.

4.1. TravelingSalesmanOptimizer design

Our optimizer class will have “core” methods, doing the bulk of the work, and “superficial” methods, serving as the high-level interface of the class, which invoke the core methods underneath.

These are the steps that will lie at the core of the optimizer’s logic:

1. Create a Pyomo model out of a distance matrix. This is done by the _create_model method, which basically wraps the code of the proof-of-concept we already did. It accepts a dataframe of a distance matrix and builds a Pyomo model out of it. The only important difference between what we did and what we’re doing is that, now, the initial site is not hard-coded as simply "hotel", but is assumed to be the site of the first row in df_distances. In the general case, thus, the initial site is taken to be the first one in the coordinates dataframe⁴ df_sites. This generalization allows the optimizer to solve any instance.
2. (Attempt to) Solve the model. This is performed in the _optimize method inherited from BaseOptimizer, which returns True only if a solution is found.
3. Extract the solution from the model and parse it in a way that is easy to interpret and use. This happens inside _store_solution_from_model, which is a method that inspects the solved model and extracts the values of the decision variables, and the value of the objective function, to create the attributes tour_ and tour_distance_, respectively. This method gets invoked only if a solution exists, so if no solution is found, the “solution attributes” tour_ and tour_distance_ never get created. The benefit of this is that the presence of these two “solution attributes”, after fitting, will inform the user of the existence of a solution. As a plus, the optimal values of both the variables and objective can be conveniently retrieved at any point, not necessarily at the moment of fitting.

The last 2 steps — finding and extracting the solution — are wrapped inside the last “core” method, _fit_to_distances.

“But wait” — you might think — “As the name implies, _fit_to_distances requires distances as input; isn’t our goal to solve TSP problems using only coordinates, not distances?”. Yes, that’s where the fit method fits in. We pass coordinates to it, and we take advantage of GeoAnalyzer to construct the distance matrix, which is then processed normally by _fit_to_distances. In this way, if the user does not want to collect the distances himself, he can delegate the task by using fit. If, however, he prefers to use custom data, he can assemble it in a df_distances and pass it to _fit_to_distances instead.

4.2. TravelingSalesmanOptimizer implementation

Let’s follow the design outlined above to incrementally build the optimizer. First, a minimalist version that just builds a model and solves it — without any solution parsing yet. Notice how the __repr__ method allows us to know the name and number of sites the optimizer contains whenever we print it.

from typing import Tuple, List
class TravelingSalesmanOptimizer(BaseOptimizer):
"""Implements the Miller–Tucker–Zemlin formulation [1] of the 
Traveling Salesman Problem (TSP) as a linear integer program. 
The TSP can be stated like: "Given a set of locations (and usually 
their pair-wise distances), find the tour of minimal distance that 
traverses all of them exactly once and ends at the same location 
it started from. For a derivation of the mathematical model, see [2].
Parameters
----------
name : str
Optional name to give to a particular TSP instance
Attributes
----------
tour_ : list
List of locations sorted in visit order, obtained after fitting.
To avoid duplicity, the last site in the list is not the initial 
one, but the last one before closing the tour.
tour_distance_ : float
Total distance of the optimal tour, obtained after fitting.
Example
--------
>>> tsp = TravelingSalesmanOptimizer()
>>> tsp.fit(df_sites)  # fit to a dataframe of geo-coordinates
>>> tsp.tour_  # list ofsites sorted by visit order
References
----------
[1] https://en.wikipedia.org/wiki/Travelling_salesman_problem
[2] https://towardsdatascience.com/plan-optimal-trips-automatically-with-python-and-operations-research-models-part-2-fc7ee8198b6c
"""
def __init__(self, name=""):
super().__init__()
self.name = name
def _create_model(self, df_distances: pd.DataFrame) -> pyo.ConcreteModel:
""" Given a pandas dataframe of a distance matrix, create a Pyomo model 
of the TSP and populate it with that distance data """
model = pyo.ConcreteModel(self.name)
# a site has to be picked as the "initial" one, doesn't matter which 
# really; by lack of better criteria, take first site in dataframe 
# as the initial one
model.initial_site = df_distances.iloc[0].name
#===========  sets declaration  ===========#
list_of_sites = df_distances.index.tolist()
model.sites = pyo.Set(initialize=list_of_sites, 
domain=pyo.Any, 
doc="set of all sites to be visited (𝕊)")
def _rule_domain_arcs(model, i, j):
""" All possible arcs connecting the sites (𝔸) """
# only create pair (i, j) if site i and site j are different
return (i, j) if i != j else None 
rule = _rule_domain_arcs
model.valid_arcs = pyo.Set(
initialize=model.sites * model.sites,  # 𝕊 × 𝕊
filter=rule, doc=rule.__doc__)
model.sites_except_initial = pyo.Set(
initialize=model.sites - {model.initial_site}, 
domain=model.sites,
doc="All sites except the initial site"
)
#===========  parameters declaration  ===========#
def _rule_distance_between_sites(model, i, j):
""" Distance between site i and site j (𝐷𝑖𝑗) """
return df_distances.at[i, j]  # fetch the distance from dataframe
rule = _rule_distance_between_sites
model.distance_ij = pyo.Param(model.valid_arcs, 
initialize=rule, 
doc=rule.__doc__)
model.M = pyo.Param(initialize=1 - len(model.sites_except_initial),
doc="big M to make some constraints redundant")
#===========  variables declaration  ===========#
model.delta_ij = pyo.Var(model.valid_arcs, within=pyo.Binary, 
doc="Whether to go from site i to site j (𝛿𝑖𝑗)")
model.rank_i = pyo.Var(model.sites_except_initial, 
within=pyo.NonNegativeReals, 
bounds=(1, len(model.sites_except_initial)), 
doc=("Rank of each site to track visit order"))
#===========  objective declaration  ===========#
def _rule_total_distance_traveled(model):
""" total distance traveled """
return pyo.summation(model.distance_ij, model.delta_ij)
rule = _rule_total_distance_traveled
model.objective = pyo.Objective(rule=rule, 
sense=pyo.minimize, 
doc=rule.__doc__)
#===========  constraints declaration  ===========#
def _rule_site_is_entered_once(model, j):
""" each site j must be visited from exactly one other site """
return sum(model.delta_ij[i, j] 
for i in model.sites if i != j) == 1
rule = _rule_site_is_entered_once
model.constr_each_site_is_entered_once = pyo.Constraint(
model.sites, 
rule=rule, 
doc=rule.__doc__)
def _rule_site_is_exited_once(model, i):
""" each site i must departure to exactly one other site """
return sum(model.delta_ij[i, j] 
for j in model.sites if j != i) == 1
rule = _rule_site_is_exited_once
model.constr_each_site_is_exited_once = pyo.Constraint(
model.sites, 
rule=rule, 
doc=rule.__doc__)
def _rule_path_is_single_tour(model, i, j):
""" For each pair of non-initial sites (i, j), 
if site j is visited from site i, the rank of j must be 
strictly greater than the rank of i.
"""
if i == j:  # if sites coincide, skip creating a constraint
return pyo.Constraint.Skip
r_i = model.rank_i[i]
r_j = model.rank_i[j]
delta_ij = model.delta_ij[i, j]
return r_j >= r_i + delta_ij + (1 - delta_ij) * model.M
# cross product of non-initial sites, to index the constraint
non_initial_site_pairs = (
model.sites_except_initial * model.sites_except_initial)
rule = _rule_path_is_single_tour
model.constr_path_is_single_tour = pyo.Constraint(
non_initial_site_pairs,
rule=rule, 
doc=rule.__doc__)
self._store_model(model)  # method inherited from BaseOptimizer
return model
def _fit_to_distances(self, df_distances: pd.DataFrame) -> None:
self._name_index = df_distances.index.name
model = self._create_model(df_distances)
solution_exists = self._optimize(model)
return self
@property
def sites(self) -> Tuple[str]:
""" Return tuple of site names the optimizer considers """
return self.model.sites.data() if self.is_model_created else ()
@property
def num_sites(self) -> int:
""" Number of locations to visit """
return len(self.sites)
@property
def initial_site(self):
return self.model.initial_site if self.is_fitted else None
def __repr__(self) -> str:
name = f"{self.name}, " if self.name else ''
return f"{self.__class__.__name__}({name}n={self.num_sites})"

Let’s quickly check how the optimizer behaves. Upon instantiation, the optimizer does not contain any number of sites, as the representation string shows, or an internal model, and it’s of course not fitted:

tsp = TravelingSalesmanOptimizer("trial 1")
print(tsp)
#[Out]: TravelingSalesmanOptimizer(trial 1, n=0)
print(tsp.is_model_created, tsp.is_fitted)
#[Out]: (False, False)

We now fit it to the distance data, and if we don’t get a warning, it means that it all went well. We can see that now the representation string tells us we provided 9 sites, there’s an internal model, and that the optimizer was fitted to the distance data:

tsp._fit_to_distances(df_distances)
print(tsp)
#[Out]: TravelingSalesmanOptimizer(trial 1, n=9)
print(tsp.is_model_created, tsp.is_fitted)
#[Out]: (True, True)

That the optimal solution was found is corroborated by the presence of definite values in the rank decision variables of the model:

tsp.model.rank_i.get_values()

{'Sacre Coeur': 8.0,
'Louvre': 2.0,
'Montmartre': 7.0,
'Port de Suffren': 4.0,
'Arc de Triomphe': 5.0,
'Av. Champs Élysées': 6.0,
'Notre Dame': 1.0,
'Tour Eiffel': 3.0}

These rank variables represent the chronological order of the stops in the optimal tour. If you recall from their definition, they are defined over all sites except the initial one⁵, and that’s why the hotel does not appear in them. Easy, we could add the hotel with rank 0, and there we would have the answer to our problem. We don’t need to extract 𝛿ᵢⱼ, the decision variables for the individual arcs of the tour, to know in which order we should visit the sites. Although that’s true, we’re still going to use the arc variables 𝛿ᵢⱼ to extract the exact sequence of stops from the solved model.

💡 Agile doesn’t need to be fragile

If our only aim were to solve the TSP, without looking to extend the model to encompass more details of our real-life problem, it would be enough to use the rank variables to extract the optimal tour. However, as the TSP is just the initial prototype of what will become a more sophisticated model, we’re better off extracting the solution from the arc decision variables 𝛿ᵢⱼ, as they will be present in any model that involves routing decisions. All other decision variables are auxiliary, and, when needed, their job is to represent states or indicate conditions dependant on the true decision variables, 𝛿ᵢⱼ. As you’ll see in the next articles, choosing the rank variables to extract the tour works for a pure TSP model, but won’t work for extensions of it that make it optional to visit some sites. Hence, if we extract the solution from 𝛿ᵢⱼ, our approach will be general and re-usable, no matter how complex the model we’re using.

The benefits of this approach will become apparent in the following articles, where new requirements are added, and thus additional variables are needed inside the model. With the why covered, let’s jump into the how.

4.2.1 Extracting the optimal tour from the model

• We have the variable 𝛿ᵢⱼ, indexed by possible arcs (i, j), where 𝛿ᵢⱼ=0 means the arc is not selected and 𝛿ᵢⱼ=1 means the arc is selected.
• We want a dataframe where the sites are in the index (as in our input df_sites), and where the stop number is indicated in the column "visit_order".
• We write a method to extract such dataframe from the fitted optimizer. These are the steps we’ll follow, with each step encapsulated in its own helper method(s):

1. Extract the selected arcs from 𝛿ᵢⱼ, which represents the tour. Done in _get_selected_arcs_from_model.
2. Convert the list of arcs (edges) into a list of stops (nodes). Done in _get_stops_order_list.
3. Convert the list of stops into a dataframe with a consistent structure. Done in _get_tour_stops_dataframe.

As the selected arcs are mixed (i.e., not in “traversing order”), getting a list of ordered stops is not that straight-forward. To avoid convoluted code, we exploit the fact that the arcs represent a graph, and we use the graph object G_tour to traverse the tour nodes in order, arriving at the list easily.

import networkx as nx
# class TravelingSalesmanOptimizer(BaseOptimizer):
# def __init__()
# def _create_model()
# def _fit_to_distances()
# def sites()
# def num_sites()
# def initial_site()
_Arc = Tuple[str, str]
def _get_selected_arcs_from_model(self, model: pyo.ConcreteModel) -> List[_Arc]:
"""Return the optimal arcs from the decision variable delta_{ij}
as an unordered list of arcs. Assumes the model has been solved"""
selected_arcs = [arc 
for arc, selected in model.delta_ij.get_values().items()
if selected]
return selected_arcs
def _extract_solution_as_graph(self, model: pyo.ConcreteModel) -> nx.Graph:
"""Extracts the selected arcs from the decision variables of the model, stores 
them in a networkX graph and returns such a graph"""
selected_arcs = self._get_selected_arcs_from_model(model)
self._G_tour = nx.DiGraph(name=model.name)
self._G_tour.add_edges_from(selected_arcs)
return self._G_tour
def _get_stops_order_list(self) -> List[str]:
"""Return the stops of the tour in a list **ordered** by visit order"""
visit_order = []
next_stop = self.initial_site  # by convention...
visit_order.append(next_stop)  # ...tour starts at initial site
G_tour = self._extract_solution_as_graph(self.model)
# starting at first stop, traverse the directed graph one arc at a time
for _ in G_tour.nodes:
# get consecutive stop and store it
next_stop = list(G_tour[next_stop])[0]
visit_order.append(next_stop)
# discard last stop in list, as it's repeated (the initial site) 
return visit_order[:-1]
def get_tour_stops_dataframe(self) -> pd.DataFrame:
"""Return a dataframe of the stops along the optimal tour"""
if self.is_fitted:
ordered_stops = self._get_stops_order_list()
df_stops = (pd.DataFrame(ordered_stops, columns=[self._name_index])
.reset_index(names='visit_order')  # from 0 to N
.set_index(self._name_index)  # keep index consistent
)
return df_stops
print("No solution found. Fit me to some data and try again")

Let’s see what this new method gives us:

tsp = TravelingSalesmanOptimizer("trial 2")
tsp._fit_to_distances(df_distances)
tsp.get_tour_stops_dataframe()

The visit_order column indicates we should go from the hotel to Notre Dame, then to the Louvre, and so on, until the last stop before closing the tour, Sacre Coeur. After that, it’s trivial that one must return to the hotel. Good, now we have the solution in a format easy to interpret and work with. But the sequence of stops is not all we care about. The value of the objective function is also an important metric to keep track of, as it’s the criterion guiding our decisions. For our particular case of the TSP, this means getting the total distance of the optimal tour.

4.2.2 Extracting the optimal objective from the model

In the same manner that we didn’t use the rank variables to extract the sequence of stops because in more complex models their values wouldn’t coincide with the tour stops, we won’t use the objective function directly to obtain the total distance of the tour, even though, here too, both measures are equivalent. In more complex models, the objective function will also incorporate other targets, so this equivalence will no longer hold.

For now, we’ll keep it simple and create a non-public method, _get_tour_total_distance, which clearly indicates the intent. The details of where this distance comes from are hidden, and will depend on the particular targets that more advanced models care about. For now, the details are simple: get the objective value of the solved model.

# class TravelingSalesmanOptimizer(BaseOptimizer):
# def __init__()
# def _create_model()
# def _fit_to_distances()
# def sites()
# def num_sites()
# def initial_site()
# def _get_selected_arcs_from_model()
# def _extract_solution_as_graph()
# def _get_stops_order_list()
# def get_tour_stops_dataframe()
def _get_tour_total_distance(self) -> float:
"""Return the total distance of the optimal tour"""
if self.is_fitted:
# as the objective is an expression for the total distance, 
distance_tour = self.model.objective()  # we just get its value
return distance_tour
print("Optimizer is not fitted to any data, no optimal objective exists.")
return None

It may look superfluous now, but it’ll serve as a reminder to our future selves that there is a design for grabbing objective values we’d better follow. Let’s check it:

tsp = TravelingSalesmanOptimizer("trial 3")
tsp._fit_to_distances(df_distances)
print(f"Total distance: {tsp._get_tour_total_distance()} m")
# [Out]: Total distance: 14931.0 m

It’s around 14.9 km. As both the optimal tour and its distance are important, let’s make the optimizer store them together whenever the _fit_to_distances method gets called, and only when an optimal solution is found.

4.2.3 Storing the solution in attributes

In the implementation of _fit_to_distances above, we just created a model and solved it, we didn’t do any parsing of the solution stored inside the model. Now, we’ll modify _fit_to_distances so that when the model solution succeeds, two new attributes are created and made available with the two relevant parts of the solution: the tour_ and the tour_distance_. To make it simple, the tour_ attribute won’t return the dataframe we did earlier, it will return the list with ordered stops. The new method _store_solution_from_model takes care of this.

# class TravelingSalesmanOptimizer(BaseOptimizer):
# def __init__()
# def _create_model()
# def sites()
# def num_sites()
# def initial_site()
# def _get_selected_arcs_from_model()
# def _extract_solution_as_graph()
# def _get_stops_order_list()
# def get_tour_stops_dataframe()
# def _get_tour_total_distance()
def _fit_to_distances(self, df_distances: pd.DataFrame):
"""Creates a model of the TSP using the distance matrix 
provided in `df_distances`, and then optimizes it. 
If the model has an optimal solution, it is extracted, parsed and 
stored internally so it can be retrieved.
Parameters
----------
df_distances : pd.DataFrame
Pandas dataframe where the indices and columns are the "cities" 
(or any site of interest) of the problem, and the cells of the 
dataframe contain the pair-wise distances between the cities, i.e.,
df_distances.at[i, j] contains the distance between i and j.
Returns
-------
self : object
Instance of the optimizer.
"""
model = self._create_model(df_distances)
solution_exists = self._optimize(model)
if solution_exists:
# if a solution wasn't found, the attributes won't exist
self._store_solution_from_model()
return self
#====================  solution handling  ====================
def _store_solution_from_model(self) -> None:
"""Extract the optimal solution from the model and create the "fitted 
attributes" `tour_` and `tour_distance_`"""
self.tour_ = self._get_stops_order_list()
self.tour_distance_ = self._get_tour_total_distance()

Let’s fit the optimizer again to the distance data and see how easy it is to get the solution now:

tsp = TravelingSalesmanOptimizer("trial 4")._fit_to_distances(df_distances)
print(f"Total distance: {tsp.tour_distance_} m")
print(f"Best tour:\n", tsp.tour_)
# [Out]:
# Total distance: 14931.0 m
# Best tour:
# ['hotel', 'Notre Dame', 'Louvre', 'Tour Eiffel', 'Port de Suffren', 'Arc de Triomphe', 'Av. Champs Élysées', 'Montmartre', 'Sacre Coeur']

Nice. But we can do even better. To further increase the usability of this class, let’s allow the user to solve the problem by only providing the dataframe of sites coordinates. As not everyone will be able to collect a distance matrix for their sites of interest, the class can take care of it and provide an approximate distance matrix. This was done above in section 3.2 with the GeoAnalyzer, here we just put it under the new fit method:

# class TravelingSalesmanOptimizer(BaseOptimizer):
# def __init__()
# def _create_model()
# def _fit_to_distances()
# def sites()
# def num_sites()
# def initial_site()
# def _get_selected_arcs_from_model()
# def _extract_solution_as_graph()
# def _get_stops_order_list()
# def get_tour_stops_dataframe()
# def _get_tour_total_distance()
# def _store_solution_from_model()
def fit(self, df_sites: pd.DataFrame):
"""Creates a model instance of the TSP problem using a 
distance matrix derived (see notes) from the coordinates provided 
in `df_sites`.
Parameters
----------
df_sites : pd.DataFrame
Dataframe of locations "the salesman" wants to visit, having the 
names of the locations in the index and at least one column 
named 'latitude' and one column named 'longitude'.
Returns
-------
self : object
Instance of the optimizer.
Notes
-----
The distance matrix used is derived from the coordinates of `df_sites`
using the ellipsoidal distance between any pair of coordinates, as 
provided by `geopy.distance.distance`."""
self._validate_data(df_sites)
self._name_index = df_sites.index.name
self._geo_analyzer = GeoAnalyzer()
self._geo_analyzer.add_locations(df_sites)
df_distances = self._geo_analyzer.get_distance_matrix(precision=0)
self._fit_to_distances(df_distances)
return self
def _validate_data(self, df_sites):
"""Raises error if the input dataframe does not have the expected columns"""
if not ('latitude' in df_sites and 'longitude' in df_sites):
raise ValueError("dataframe must have columns 'latitude' and 'longitude'")

And now we have achieved our goal: find the optimal tour from just the sites locations (and not from the distances as before):

tsp = TravelingSalesmanOptimizer("trial 5")
tsp.fit(df_sites)
print(f"Total distance: {tsp.tour_distance_} m")
tsp.tour_
#[Out]:
# Total distance: 14931.0 m
# ['hotel',
# 'Notre Dame',
# 'Louvre',
# 'Tour Eiffel',
# 'Port de Suffren',
# 'Arc de Triomphe',
# 'Av. Champs Élysées',
# 'Montmartre',
# 'Sacre Coeur']

4.3. TravelingSalesmanOptimizer for dummies

Congratulations! We reached the point where the optimizer is very intuitive to use. For mere convenience, I’ll add another method that will be quite helpful later on when we do [sensitivity analysis] and compare the results of different models. The optimizer, as it is now, tells me the optimal visit order in a list, or in a separate dataframe returned by get_tour_stops_dataframe(), but I’d like it to tell me the visit order by transforming the locations dataframe that I give it directly—by returning the same dataframe with a new column having the optimal sequence of stops. The method fit_prescribe will be in charge of this:

# class TravelingSalesmanOptimizer(BaseOptimizer):
# def __init__()
# def _create_model()
# def sites()
# def num_sites()
# def initial_site()
# def _get_selected_arcs_from_model()
# def _extract_solution_as_graph()
# def _get_stops_order_list()
# def get_tour_stops_dataframe()
# def _get_tour_total_distance()
# def _fit_to_distances()
# def _store_solution_from_model()
# def fit()
# def _validate_data()
def fit_prescribe(self, df_sites: pd.DataFrame, sort=True) -> pd.DataFrame:
"""In one line, take in a dataframe of locations and return 
a copy of it with a new column specifying the optimal visit order
that minimizes total distance.
Parameters
----------
df_sites : pd.DataFrame
Dataframe with the sites in the index and the geolocation 
information in columns (first column latitude, second longitude).
sort : bool (default=True)
Whether to sort the locations by visit order.
Returns
-------
df_sites_ranked : pd.DataFrame 
Copy of input dataframe `df_sites` with a new column, 'visit_order', 
containing the stop sequence of the optimal tour.
See Also
--------
fit : Solve a TSP from just site locations.
Examples
--------
>>> tsp = TravelingSalesmanOptimizer()
>>> df_sites_tour = tsp.fit_prescribe(df_sites)  # solution appended
"""
self.fit(df_sites)  # find optimal tour for the sites
if not self.is_fitted:  # unlikely to happen, but still
raise Exception("A solution could not be found. "
"Review data or inspect attribute `_results` for details."
)
# join input dataframe with column of solution
df_sites_ranked = df_sites.copy().join(self.get_tour_stops_dataframe())    
if sort:
df_sites_ranked.sort_values(by="visit_order", inplace=True)
return df_sites_ranked

Now we can solve any TSP in just one line:

tsp = TravelingSalesmanOptimizer("Paris")
tsp.fit_prescribe(df_sites)

If we’d like to conserve the original order of locations as they were in df_sites, we can do it by specifying sort=False:

tsp.fit_prescribe(df_sites, sort=False)

And if we’re curious we can also check the number of variables and constraints the internal model needed to solve our particular instance of the TSP. This will be handy when doing debugging or performance analysis.

tsp.print_model_info()
#[Out]:
# Name: Paris
# - Num variables: 80
# - Num constraints: 74
# - Num objectives: 1