An Elegant Approach to Effectively Solving Salesman Problems with Python | by Carlos J. Uribe

Here we will not start from scratch. As noted above, we have already developed code that builds a Pyomo model of the TSP and solves it in race 3. And believe me, that was the hardest part. Now we have the easier task of organizing what we did so that it is general, hiding the details and keeping the essential visible elements. In a sense, we want the optimizer to feel like a “magic box” that even users who are not familiar with mathematical modeling can use to solve their TSP problems intuitively.

4.1. TravelingSalesmanOptimizer design

Our optimizer class will have “core” methods, which will do most of the work, and “shallow” methods, which will serve as the class's high-level interface, which will invoke the underlying main methods.

These are the steps that will be at the center of the optimizer logic:

1. Create a Pyomo model outside of a distance matrix. This is done by _create_model method, which basically wraps the proof of concept code we already did. Accepts a data frame of a distance matrix and builds a Pyomo model from it. The only major difference between what we did and what we are doing is that, now, the initial site is not coded as simply "hotel"but it is fictional be the front row seat in df_distances. Thus, in the general case, the starting site is considered the first in the coordinate data frame⁴ df_sites. This generalization allows the optimizer to resolve any instance.
2. (Attempt to) Solve the model. This is done in the _optimize method inherited from BaseOptimizerthat returns True only if a solution is found.
3. Extract the solution from the model and analyze it. in a way that is easy to interpret and use. this happens inside _store_solution_from_modelwhich 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 is invoked only if there is a solutionso if no solution is found, the “solution attributes” tour_ and tour_distance_ never be created. The benefit of this is that the presence of these two “solution attributes”, after tuning, will inform the user of the existence of a solution. As an advantage, the optimal values ​​of both the variables and the objective can be conveniently recovered at any time, not necessarily at the time of adjustment.

The last 2 steps (find and extract the solution) are included within the last “basic” method. _fit_to_distances.

“But wait,” you might think, “as the name implies, _fit_to_distances requires distances as input; Isn't it our goal to solve TSP problems by using only coordinatesNo distances?”. Yes, that's where fit method fits in. we passed coordinates to it, and we take advantage 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 using fit. However, if you prefer to use custom data, you can assemble it into a df_distances and pass it to _fit_to_distances instead.

4.2. TravelingSalesmanOptimizer implementation

Let's follow the design described above to incrementally build the optimizer. First, a minimalist version that simply builds a model and solves it, without any solution analysis yet. Notice how the __repr__ The method allows us to know the name and number of sites the optimizer contains each time 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. When instantiated, the optimizer does not contain any number of sites, as shown by the rendering string, nor an internal model, and is of course not tuned:

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

Now we adjust it to the distance data and if we don't get a warning, it means everything went well. We can see that now the render chain tells us that we provided 9 sites, that there is an internal model, and that the optimizer adjusted 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)

The fact that the optimal solution was found is corroborated by the presence of defined values ​​in the range 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 ranking variables represent the chronological order of stops on the optimal route. If you remember from their definition, they are defined on all sites except the initial one⁵, and that is why the hotel does not appear on them. Easy, we could add the hotel with rank 0, and there we would have the answer to our problem. We do not need to extract 𝛿ᵢⱼ, the decision variables for the individual arcs of the route, to know in what order we should visit the sites. Although that is true, we will still use the arc variables 𝛿ᵢⱼ to extract the exact sequence of stops from the solved model.

Agile does not need to be fragile

If our only objective were to resolve the TSP, without looking at extend For the model to encompass more details of our real-life problem, it would be sufficient to use the range variables to extract the optimal path. However, since the TSP is only the initial prototype of what will become a more sophisticated model, we better extract the solution from the decision variables of the arc 𝛿ᵢⱼ, since they will be present in any model involving routing decisions. All other decision variables are auxiliary and, when necessary, their job is to represent states or indicate conditions that depend on the situation. TRUE decision variables, 𝛿ᵢⱼ. As you will see in upcoming articles, choosing the classification variables to extract the traversal works for a pure TSP model, but will not work for extensions of it that do optional to visit some places. Therefore, if we extract the solution of 𝛿ᵢⱼ, Our approach will be general and reusable, no matter how complex the model we are using..

The benefits of this approach will be evident in the following articles, where new requirements are added and therefore additional variables are needed within the model. With the because covered, let's jump to the as.

4.2.1 Extraction of the optimal model path

• Have the variable 𝛿ᵢⱼ, indexed by possible arcs (i, j), where 𝛿ᵢⱼ=0 means that the arc is not selected and 𝛿ᵢⱼ=1 means that the arc is selected.
• We want a data frame where the sites are in the index (like in our entry df_sites), and where the stop number is indicated in the column "visit_order".
• We wrote a method to extract said data frame from the installed optimizer. These are the steps we will follow, with each step encapsulated in its own helper method(s):

1. Extracts the selected arcs of 𝛿ᵢⱼ, which represents the path. Made in _get_selected_arcs_from_model.
2. Convert the list of arcs (edges) to a list of stops (nodes). Made in _get_stops_order_list.
3. Convert the list of stops into a data frame with a consistent structure. Made in _get_tour_stops_dataframe.

As the selected arcs are mixed (that is to say, not in “travel order”), obtaining an ordered list of stops is not so simple. To avoid complicated code, we take advantage of the fact that arcs represent a graphicand we use the graphic object G_tour to loop through the traversal nodes in order, reaching 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()

He visit_order The column indicates that we must 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 is trivial to have to return to the hotel. Well, now we have the solution in a format that is easy to interpret and work with. But the sequence of stops is not the only thing that matters to us. The objective function value is also an important metric to track, since it is the criterion that guides our decisions. For our particular case of the TSP, this means obtaining the total distance of the optimal route.

4.2.2 Extraction of the optimal model objective

In the same way that we do not use the range variables to extract the sequence of stops because in more complex models their values ​​would not coincide with the stops of the tour, we will not use the objective function. directly to obtain the total distance of the route, although here both measurements are equivalent. In more complex models, the objective function will also incorporate other goals.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 intention. The details of where this distance comes from are hidden and will depend on the particular objectives that most advanced models are interested in. For now, the details are simple: get the target value from 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 seem superfluous now, but it will serve as a reminder for our future that there is a design for capturing objective values ​​that we had 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 is about 14.9 km. Since both the optimal path and its distance are important, let's have the optimizer store them together whenever the _fit_to_distances the method is 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 simply created a model and solved it, we did not analyze the solution stored within the model. Now we will modify _fit_to_distances so that When the model solution is successful, 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_ The attribute will not return the data frame we made earlier, it will return the list with sorted 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 tune the optimizer to the distance data again 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. For more increase the usability of this class, let the user solve the problem by providing only the site coordinates data frame. Since not everyone will be able to compile a distance matrix for their sites of interest, the class can take care of this and provide an approximate distance matrix. This was done earlier in section 3.2 with the GeoAnalyzerhere we simply put it under the new one 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 optimal tour only from site locations (and not from distances like 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've gotten to the point where the optimizer is very intuitive to use. For the sake of convenience, I will add another method that will be very useful later when we do this (sensitivity analysis) and compare the results of different models. The optimizer as it stands now tells me the optimal visiting order in a list or in a separate data frame returned by get_tour_stops_dataframe()but I would like you to tell me the order of visit transforming the locations data frame which I give it directly, returning the same data frame with a new column that has the optimal sequence of stops. The method fit_prescribe will take 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()
# 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 A line:

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

If we wanted to preserve the original order of the locations as they were in df_siteswe can do it by specifying sort=False:

tsp.fit_prescribe(df_sites, sort=False)

And if we are curious, we can also check the number of variables and constraints that the internal model needs to solve our particular instance of the TSP. This will be useful when performing debugging or performance analysis.

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